Massive galaxy clusters like “El Gordo” hint at primordial quantum diffusion

Jose María Ezquiaga,<fnref n="1" />,²,* Juan García-Bellido,³,† and Vincent Vennin⁴,⁵,‡

¹Niels Bohr International Academy, Niels Bohr Institute,
Blegdamsvej 17, DK-2100 Copenhagen, Denmark

²NASA Einstein fellow; Kavli Institute for Cosmological Physics and Enrico Fermi Institute,
The University of Chicago, Chicago, IL 60637, USA

³Instituto de Física Teórica UAM-CSIC, Universidad Autónoma de Madrid, Cantoblanco, Madrid, 28049 Spain

⁴Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL,
CNRS, Sorbonne Université, Université Paris Cité, F-75005 Paris, France

⁵Laboratoire Astroparticule et Cosmologie, CNRS Université de Paris,
10 rue Alice Domon et Léonie Duquet, 75013 Paris, France

(Dated: February 21, 2023)

It is generally assumed within the standard cosmological model that initial density perturbations are Gaussian at all scales. However, primordial quantum diffusion unavoidably generates non-Gaussian, exponential tails in the distribution of inflationary perturbations. These exponential tails have direct consequences for the formation of collapsed structures in the universe, as has been studied in the context of primordial black holes. We show that these tails also affect the very-large-scale structures, making heavy clusters like “El Gordo”, or large voids like the one associated with the cosmic microwave background cold spot, more probable. We compute the halo mass function and cluster abundance as a function of redshift in the presence of exponential tails. We find that quantum diffusion generically enlarges the number of heavy clusters and depletes subhalos, an effect that cannot be captured by the famed $f_{NL}$ corrections. These late-universe signatures could thus be fingerprints of quantum dynamics during inflation that should be incorporated in N-body simulations and checked against astrophysical data.

**Introduction.** The standard cosmological model (ΛCDM), provides an excellent fit to the high-precision astrophysical and cosmological observations, in particular the Cosmic Microwave Background (CMB), the Large-Scale Structure (LSS) of the universe and the relative abundance of light elements. Its three main ingredients are: (1) General Relativity and the Cosmological Principle, (2) a universe made of baryonic matter, dark matter, radiation and dark energy, and (3) quasi scale-invariant, Gaussian initial density fluctuations.

However, hints for a few cracks start to emerge at different stages, in the form of moderate statistical tensions in parameter inference, e.g. the local expansion rate [1, 2], or via the existence of “extreme” objects or outliers, which are more frequently observed than what ΛCDM predicts. Those may be associated with either an extremely low value of the density field (such as the Eri-danus supervoid [3, 4], which seems to have a direct connection with the CMB cold spot [5]), or with extremely large values of the density field (such as massive galaxy clusters like *El Gordo* [6] – see however Ref. [7] for a recent smaller estimate of its mass – and the presence of galaxies and Quasi-Stellar Objects at extremely high redshifts, where according to standard ΛCDM there should not be any [8, 9]). In addition to the early structure formation issues, there are late time miss-matches at small scales such as the substructure problems [10], the too-

big-to-fail and the core-cusp problems [11], which could be alleviated by incorporating baryonic physics [12].

While most attempts to reconcile those potential issues focus on relaxing either the first or the second assumption mentioned above (i.e. modifying the laws of gravity, or in-voking the existence of additional components in the universe), a natural strategy to accommodate the existence of extreme objects within the ΛCDM paradigm would be to question the third assumption, namely the Gaussianity of the primordial density fluctuations. The reason is twofold: experimentally, there are more extreme objects than what Gaussian tails suggest, pointing toward the existence of heavier tails; theoretically, the typical mechanisms producing primordial cosmological perturbations anyway lead to non-Gaussian tails.

In the early universe indeed, vacuum quantum fluctuations are amplified by gravitational instability and stretched to large distances, giving rise to classical fluctuations in the density field, that later collapse into cosmological structures.<fnref n="1" /> At leading order in cosmological perturbation theory, it gives rise to Gaussian perturbations, in good agreement with CMB measurements [14]. However, the CMB gives access to large scales only and leaves small scales mostly unconstrained. Moreover, even at large scales, they restrict the statistics of the most likely fluctuations only, i.e. they reconstruct only the neighbourhood of the maximum of the underlying distribution functions, and say little about their tails.

<footnote>*Electronic address: jose.ezquiaga@nbi.ku.dk</footnote>

<footnote>†Electronic address: juan.garciabellido@uam.es</footnote>

<footnote>‡Electronic address: vincent.vennin@ens.fr</footnote>

<footnote n="1">¹ This mechanism is mostly studied in the context of inflation, but it also operates in most of its alternatives [13].</footnote>
Nonetheless, beyond linear order, those tails are ex-
pected to be non-Gaussian. The difficulty when charac-
terising the statistics of those tails is that they require
non-perturbative techniques. Perturbative approaches,
such as calculations of the bi- or tri-spectrum, and $f_{\text{NL}}$-
like parametrizations, are tailored to describe small de-
viations from Gaussianity around the maximum, not ac-
counting for the tails.

**Quantum diffusion and non-Gaussian tails.** Re-
cently, non-perturbative techniques have been developed
to study how quantum diffusion, the presence of which
is inevitable in scenarios where cosmological perturba-
tions have a quantum origin, modifies the expansion dy-
namics of the universe and thus affects the statistics
of density fluctuations. This can be done by combin-
ing three approaches to describe the dynamics of super-
Hubble degrees of freedom. First, the separate-universe
picture [15–18] (which is valid beyond slow roll [19, 20]),
according to which spatial gradients can be neglected on
super-Hubble scales, and each spatial point evolves inde-
pendently along the dynamics of an unperturbed uni-
verse. Second, stochastic inflation [21, 22], in which
quantum fluctuations act as a stochastic noise on the
classical, background evolution of each of these separate
universes. Third, the $\delta N$ formalism [18, 23–25], which
states that, in each of these separate universes, the local
fluctuation in the amount of expansion realized between
an initial flat hypersurface and a final hypersurface of
uniform energy density is nothing but the curvature per-
turbation. This gives rise to the stochastic-$\delta N$ formal-
ism [26–31], which provides a non-perturbative scheme
to compute the statistics of curvature perturbations on
super-Hubble scales. These methods now extend to the
calculation of the density contrast and the compaction
function [32].

While these techniques recover quasi-Gaussian distri-
butions close to their maximum, with $f_{\text{NL}}$-type correc-
tions, they also reveal the existence of systematic ex-
ponential tails [32–39], which strongly deviate from the
Gaussian profile (such heavy tails were also found in
Refs. [40–43] using different methods). More precisely,
the distribution function of the first-passage time $\mathcal{N}$ can
be expanded as $P(\mathcal{N}) = \sum_{n \ge 0} a_n(\mathbf{\Phi})e^{-\Lambda_n \mathcal{N}}$, in which
$\Lambda_n$ are the eigenvalues of the adjoint Fokker-Planck op-
erator associated with the stochastic problem under con-
sideration, and $a_n(\mathbf{\Phi})$ are coefficients that depend on
the initial configuration in field space (here denoted as
$\mathbf{\Phi}$). Far on the tail, the smallest eigen-value dominates,
$P(\mathcal{N}) \propto e^{-\Lambda_0 \mathcal{N}}$, which implies that large perturbations
are much more likely than what a Gaussian behavior,
$P_G \propto e^{-\Lambda \mathcal{N}^2}$, would suggest. In practice, these exponen-
tial tails are more important in models where quantum
diffusion dominates at some stages of the inflationary
dynamics (leading to smaller values of $\Lambda_0$, hence heav-
ier tails). Depending on the time at which this hap-
pens, they affect structures at different scales. If the
fluctuations are large enough, they may even collapse
into black holes upon horizon re-entry after inflation.

This is why non-Gaussian tails have been mostly studied
in the context of primordial-black-hole production (see
e.g. Refs. [32, 33, 36, 38, 44–49]).

Nonetheless, as we argue here, these heavy tails may
also play a key role in the formation of the LSS and
point toward potential solutions to some of the problems
of ΛCDM. Importantly, while the parameters $a_n$ and $\Lambda_n$
depend on the details of the model under consideration,
the existence of these exponential tails is ubiquitous and
arises in any model where quantum diffusion is at play. In
this sense, they are already embedded in the ΛCDM sce-
nario. Therefore, our approach does not rely on extend-
ing ΛCDM to solve the above-mentioned issues: our goal
is rather to point out that ΛCDM may already contain
the ingredients needed to explain those “anomalous” ob-
servations, provided we carefully compute the primordial
statistics beyond the perturbative level.

Heavy tails in the form of lognormal distributions are
already known to develop on sub-Hubble scales after in-
flation, due to gravitational collapse [50–52]. However,
the effect we are considering here is different: it leads
to primordial heavy tails, which are present even before
Hubble re-entry.

**Exponential tails in the primordial statistics of perturbations.** The details of the stochastic distribution associated with primordial perturbations depend on the specifics of the inflationary model (the number of fields, their potential, their kinetic coupling etc.). In order to describe the amplitude of fluctuations coarse-grained at a certain scale, one has to convolve the first-passage time distributions against backward distributions of the field value [32, 37]. Moreover, one must account for the non-linear mapping between the curvature perturbation and the density contrast [53], which further modifies distribution functions and can also introduce heavy tails [54, 55]. In this work, we do not aim at deriving predictions for specific models, but rather wish to explore generic consequences arising from the presence of heavy tails. This is why, in practice, we consider two normalized templates for the distribution function of the density contrast in comoving threading $\delta$,

$$
\begin{align}
P_2(\delta_k) &= -\frac{\pi}{2\mu^2} \vartheta'_2 \left( \frac{\pi\alpha_k}{2}, e^{-\frac{\pi^2}{\mu^2} \mathcal{D}_k} \right), \\
P_4(\delta_k) &= \frac{\pi}{2\mu^2 \alpha_k} \vartheta'_4 \left( \frac{\pi\alpha_k}{2}, e^{-\frac{\pi^2}{\mu^2} \mathcal{D}_k} \right). \tag{1}
\end{align}
$$

In these expressions, $\delta_k$ denotes the Fourier mode of
the density contrast, related to the positive variable $\mathcal{D}_k$
through the relation $\delta_k = \mathcal{D}_k - \langle \mathcal{D}_k \rangle$ where the mean
value is taken with respect to the distribution function
in question. These distributions depend on two param-
eters, $\alpha_k$ and $\mu$, the latter being scale independent to
reflect the fact that the eigenvalues $\Lambda_n$ do not depend on
the field configuration, hence on the scale [34]. Finally,
$\vartheta'_2$ and $\vartheta'_4$ are the derivatives of the elliptic theta func-
tions of the second and fourth kind respectively [56]. In
<figure><img src="image_1.png" /><figcaption>FIG. 1: Gaussian, elliptic, local-$f_{NL}$ and lognormal distributions, as a function of $\delta/\sigma$ where $\sigma$ is the standard deviation of the corresponding distribution. The free parameters of those distributions are set such that they all share the same value of $\sigma$ around the maximum and are given by $\alpha = 0.5$, $\mu = \pi$ and $\sigma = \sqrt{2}\alpha^2$, with the same $\alpha$ and $\mu$ for both elliptic functions. See Supplemental Material for how to match the shape around the maximum.</figcaption></figure>

what follows, they are refereed to as the “elliptic 2” and “elliptic 4” templates respectively. Such functions are often found in toy models of quantum diffusion [33, 34].

The two distributions are displayed in Fig. 1 as a function of $\delta/\sigma$ where hereafter $\sigma$ denotes the standard deviation of the distribution under consideration, and where they are compared with a Gaussian distribution, a local $f_{NL}$ distribution and a lognormal distribution. The free parameters of those distributions are set such that they are maximal at the same location and all share the same value of $\sigma$, see the Supplemental Material where this procedure is further detailed. Both elliptic profiles are endowed with a heavier upper tail, and with a lighter lower tail, than the Gaussian fit.

The local $f_{NL}$ parametrisation is defined as

$$\delta(x) = \delta_G(x) + \frac{3}{5} f_{NL} \left[ \delta_G^2(x) - \sigma_G^2 \right], \quad (2)$$

where $\delta_G$ has a Gaussian distribution function centered at zero and with dispersion $\sigma_G \equiv \langle \delta_G^2 \rangle^{1/2}$. From this expression, one can show that

$$P_{\text{NL}}(\delta) = \frac{1}{\sqrt{2\pi\sigma_G^2\Delta}} \left[ e^{-\frac{25(\sqrt{\Delta}-1)^2}{72f_{\text{NL}}\sigma_G^2}} + e^{-\frac{25(\sqrt{\Delta}+1)^2}{72f_{\text{NL}}\sigma_G^2}} \right], \quad (3)$$

where $\Delta(\delta) = 1 + 12/5 f_{NL}\delta + 36/25 f_{NL}^2\sigma_G^2$. As shown in Fig. 1, although $f_{NL}$ correctly describes the non-Gaussian corrections around the maximum, it fails to capture the highly non-Gaussian tails. Similarly, since $P_{NL}$ diverges when $\delta$ approaches $-3f_{NL}\sigma^2/5 - 5/(12f_{NL})$, it cannot

properly describe the small-$\delta$ statistics.<fnref n="2" /> Interestingly, the elliptic functions are more similar to a lognormal distribution than the perturbative $f_{NL}$ approximation.

**Implications for the Large-Scale Structure.** The simplest statistics to be extracted from the primordial density fluctuations is the one-point function, i.e. the number of collapsed objects. Following the Press-Schechter formalism [57], this is given by the probability that $\delta$ is above a given threshold $\delta_c$, $\beta = P(\delta > \delta_c) = 2 \int_{\delta_c}^{\infty} P(\delta) d\delta$. $\delta_c$ depends on the time of re-entry of the fluctuations and has been extensively explored in the literature [58–61]. For our purposes it will be enough to fix it to $\delta_c = 1.68$, as predicted by linear theory of spherical collapse [62]. From Fig. 1, it is clear that as $\nu \equiv \delta_c/\sigma$ increases, the number of collapsed objects is larger in heavy-tailed models than in the Gaussian case.

More precisely, let us study how structures distribute across different masses. This can be achieved with the Halo Mass Function (HMF), defined from the mass fraction $\beta$ as

$$\frac{dn}{d \ln M} = \frac{\rho_m}{M} \frac{d\beta}{d \ln M} = \frac{\rho_m}{M} \frac{d \ln \sigma^{-1}}{d \ln M} \nu \beta'(\nu), \quad (4)$$

where $M$ is the mass of the halo, $\rho_m$ is the energy density of matter and a prime denotes derivation with respect to $\nu$. Previous works have proposed to test $f_{NL}$ with the HMF, see e.g. [63–67]. Here we extend those results by exploring initial density perturbations with non-Gaussian tails.

Since present observations show a good agreement with the Gaussian hypothesis at galactic scales, we tune the free parameters of all considered distributions such that they peak at the same value and share the same standard deviation (see Supplemental Material for further details). As a consequence, the only difference in Eq. (4) comes from the term $\beta'(\nu)$. In a Gaussian distribution, one has $\beta'_G(\nu) = -2e^{-\nu^2/2}\sqrt{\pi}$, and similar expressions can be obtained for the other distributions.

Let us now study the redshift evolution of the number of halos. We can describe the HMF (4) as a function of redshift by writing $\rho_m(a) = \Omega_m(a)\rho_c$ with $\Omega_m(a) = \Omega_m/a^3$ and $\sigma(a) = \sigma(1)D(a)$, with the growth function $D(a) = \delta(a)/\delta(1)$ given by [68] $D(a) \propto a \times {}_2F_1[(w-1)/2w, -1/3w, (6w-5)/6w, 1-1/\Omega_m(a)]$, where $w$ is the equation-of-state of Dark Energy (set to) $w = -1$ and $a$ is the scale factor. By comparing the HMF at different redshifts with the abundance

<footnote n="2">² The local $f_{NL}$ in Eq. (2) is assumed to be constant, and thus cannot reproduce the scale-dependence of the amplitude of the tails that comes from the fact that quantum diffusion has a different effect at different points of the inflationary evolution. Other formulations of primordial non-Gaussianity take into account the “scale-dependence” of $f_{NL}$ through calculations of the full bispectrum, but they would still be perturbative.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 2: Halo mass function (i.e. differential number of halos per comoving volume) obtained from different distributions for the primordial density perturbations: Gaussian, elliptic 2 and 4, and local $f_{\text{NL}}$ (where $f_{\text{NL}}$ is fixed at the last scattering surface). Each column corresponds to the HMF at a different redshift. The bottom panels show the ratio between the HMF and the Gaussian result. Quantum diffusion affects both low- and high-mass ends of the HMF and become more significant at higher redshifts, making their signatures distinguishable from perturbative non-Gaussianities ($f_{\text{NL}}$). The normalization is fixed to match the Gaussian at $M = 10^{11}h^{-1}M_{\odot}$ and $z = 0$, where $h$ is the dimensionless Hubble constant, $h = H_0/(100 \text{ km/s/Mpc})$.</figcaption></figure>

of massive clusters, we can estimate whether e.g. “El Gordo” is a typical cluster or not, at a given redshift.

Our main results are presented in Fig. 2, where we display the HMF for the four distributions under consideration at three redshifts, $z = 0, 1, 7$. The bottom panels display the ratio of the HMF with respect to the Gaussian case, with the normalization fixed to the Gaussian at $M = 10^{11}h^{-1}M_{\odot}$. At $z = 0$ we observe two main effects from the exponential tails: an increase in the number of clusters ($M > 10^{13}M_{\odot}$) and a decrease in the number of substructures ($M < 10^{9}M_{\odot}$). Interestingly, the $f_{\text{NL}}$ reconstruction can partially mimic the increase in clusters, but not the decrease in substructures. This makes the predictions of quantum diffusion falsifiable. Equally important, the predictions of the elliptic HMF can potentially alleviate the shortcomings of ΛCDM. Moreover, we also observe that the redshift evolution of the HMF is a key discriminator of the nature of the primordial perturbations. An elliptic HMF predicts that more massive objects formed earlier, in agreement with the recent detection of massive, high-redshift objects (see e.g. [69] for a recent census of the age of young quasars).

In addition to the number of halos per unit mass, it is interesting to compute the number of clusters as a function of redshift. This can be probed directly for example with CMB data using the Sunyaev-Zeldovich (SZ) effect [70–72]. Our results are presented in Fig. 3, focusing on clusters with $M > 10^{15}M_{\odot}$ (see Supplemental Material for the detailed calculation). One can clearly

<figure><img src="image_3.png" /><figcaption>FIG. 3: Number of clusters with mass larger than $10^{15}M_{\odot}$ in redshift bins of $\Delta z = 0.1$. as a function of redshift for the Gaussian, elliptic 2 and 4, and $f_{\text{NL}}$ distributions. The normalization is fixed to the Gaussian case at $z = 0$. We highlight with a vertical dashed line the redshift of El Gordo ($z = 0.87$).</figcaption></figure>

see that, already beyond $z \sim 1$, the number of clusters is much enhanced when initial perturbations have heavy tails. This, again, shows the potential of this method to constrain the very early universe physics.

Observationally, in most cases we do not have direct access to the HMF, but rather to the amount of lumi-
nous matter. One thus needs to take into account the as-
trophysical systematics connecting these two. Recently,
constraints on $f_{\text{NL}}$ have been derived using UV galaxy
luminosity functions that marginalize over those system-
atics [73]. A natural extension of this work would thus be
to constrain the heavy tails from quantum diffusion with
these data. Moreover, the HMF at sub-galactic scales
could be probed analyzing the strong lensing rates and
magnifications [74].

**Future prospects.** ΛCDM relies on the assumption of Gaussian initial conditions. Although CMB observations tightly constrain the amount of non-Gaussianities at large scales, little is known about the primordial fluctuations at smaller scales. Several processes in the early universe could lead to non-Gaussian distributions. Notably, an inevitable exponential tail arises due to quantum diffusion during inflation. In this work we have studied the imprints these heavy tails leave in the number of halos and their mass function. We have found that they enhance the number of heavy clusters and deplete the number of sub-halos, and that this difference with respect to the standard Gaussian initial conditions becomes more important at high redshift, depending on the strength of quantum diffusion. This could be compared with current SZ catalogs (e.g. Fig. 18 in [72]) that did not find clusters of $M > 10^{16} M_{\odot}$ and has a redshift distribution peaking at $z < 1$. However, there are outstanding clusters like El Gordo with $M \sim 3 \times 10^{15} M_{\odot}$ at $z = 0.87$ [6], and deep voids like the Eridanus supervoid [5], and many more should soon be discovered with the James Webb Space Telescope (JWST) [75], Euclid [76] and the Vera Rubin Observatory (LSST) [77].

Let us note that the effect of quantum diffusion is
similar to having a lognormal initial density distribu-
tion. Such lognormal profiles are indeed typically ob-
tained from Gaussian initial conditions due to non-linear
gravitational collapse [50]. In fact, to speed up the com-
putation, many *N*-body codes start their evolution with
an already lognormal distribution [78, 79]. Here we find
that non-linear growth occurs at much earlier times, as
soon as non-Gaussian tails are present. This leaves spe-
cific features in the redshift-dependence of the statistics:
for high-redshift galaxies we expect highly non-Gaussian
statistics, much beyond what would be expected from
non-linear gravitational collapse in such a short time.

Using the HMF to probe primordial universe physics
requires further developments on various fronts. From
the observational side, we need to understand the sys-
tematics behind high-mass supergalactic structures at
low and high redshifts, which Euclid [76] and JWST [75]
observations may help to alleviate, as well as issues
with baryonic physics and Halo Occupation Distribu-
tions. The HMF is sensitive to the one-point statistics of
the density field only, but it would also be interesting to
consider observables probing non-Gaussianities in higher
correlators. From the theoretical side, we need to imple-
ment realistic physical models in our pipeline, and go be-

yond the simple phenomenological prescription adopted
here (although this is not expected to alter our qualita-
tive conclusions, it matters for quantitative details). Fi-
nally, on the numerical side, it would be necessary to run
*N*-body simulations with non-Gaussian initial conditions
of the type described above.

Altogether we have demonstrated the impact of quan-
tum diffusion on the LSS and how the HMF and cluster
abundances can probe early-universe physics. More im-
portantly, we have shown that, within the standard cos-
mological model itself, quantum diffusion is inevitable
during inflation, and some of the current tensions can be
alleviated thanks to the non-Gaussian nature of the tails
of primordial perturbations.

Acknowledgments

JGB acknowledges funding from the Research Project PGC2018-094773-B-C32 (MINECO-FEDER) and the Centro de Excelencia Severo Ochoa Program CEX2020-001007-S. JME is supported by the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 847523 INTERACTIONS, and by VILLUM FONDEN (grant no. 53101 and 37766). He was also supported by NASA through the NASA Hubble Fellowship grant HST-HF2-51435.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555; and by the Kavli Institute for Cosmological Physics through an endowment from the Kavli Foundation and its founder Fred Kavli.

Supplemental Material

In this supplemental material, we discuss how the el-
liptic templates can be approximated by other profiles,
and where such approximations are expected to be valid.

I. LOCAL $f_{\text{NL}}$ DISTRIBUTION

Let us consider the elliptic profiles given in Eq. (1) of
the main text, and see how an effective $f_{\text{NL}}$ distribution,
see Eq. (3) of the main text, can be associated to them.

A. Matching the first moments

A first approach consists in arranging the parameters
of the $f_{\text{NL}}$ distribution such that the three first moments
of both profiles coincide. The first moments of the $f_{\text{NL}}$
distribution can be calculated either from evaluating the
expectation value of Eq. (2) taken to some integer power,
or from integrating Eq. (3) of the main text directly. For
the three first moments, one obtains

$$ \begin{aligned} \langle \delta \rangle_{\text{NL}} &= 0, \\ \langle \delta^2 \rangle_{\text{NL}} &= \sigma_G^2 \left( 1 + \frac{18}{25} f_{\text{NL}}^2 \sigma_G^2 \right), \\ \langle \delta^3 \rangle_{\text{NL}} &= \frac{18}{5} f_{\text{NL}} \sigma_G^4 \left( 1 + \frac{12}{25} f_{\text{NL}}^2 \sigma_G^2 \right). \end{aligned} \quad (5) $$

These relations can be inverted to obtain $\sigma_G^2$ and $f_{\text{NL}}$ in terms of $\langle \delta^2 \rangle$ and $\langle \delta^3 \rangle$, leading to

$$ \begin{aligned} \sigma_G^2 &= \langle \delta^2 \rangle_{\text{NL}} \left( X + \frac{1}{X} - 1 \right) \\ f_{\text{NL}} &= \frac{5}{3\sqrt{2\langle \delta^2 \rangle_{\text{NL}}}} \frac{\sqrt{2 - X - \frac{1}{X}}}{|X + \frac{1}{X} - 1|} \text{sign}(\langle \delta^3 \rangle_{\text{NL}}) \end{aligned} \quad (6) $$

where

$$ X = \left[ 1 + \frac{\langle \delta^3 \rangle_{\text{NL}}^2}{4 \langle \delta^2 \rangle_{\text{NL}}^3} \left( \sqrt{1 - 8 \frac{\langle \delta^2 \rangle_{\text{NL}}^3}{\langle \delta^3 \rangle_{\text{NL}}^2}} - 1 \right) \right]^{1/3} . \quad (7) $$

One may note that, when $8 \langle \delta^2 \rangle_{\text{NL}}^3 > \langle \delta^3 \rangle_{\text{NL}}^2$, X is complex, but in that case one can readily check that $|X| = 1$, hence $X + 1/X$ is real.

From the three first moments of a given distribution, one can thus extract the corresponding $\sigma_G^2$ and $f_{\text{NL}}$ parameters, hence the $f_{\text{NL}}$ distribution that reproduces those moments.

## B. Matching the behaviour around the maximum

The above procedure provides a good fit around the mean value of the reference distribution. However, elliptic profiles are such that the location of the mean and the location of the maximum are substantially different, and in practice it is more efficient to use the $f_{\text{NL}}$ distribution that best describes the behaviour around the maximum of the PDF.

In the regime where $\Delta \gg (\sigma_G f_{\text{NL}})^4$, the second branch in Eq. (3) of the main text can be neglected, and one can approximate

$$ P_{\text{NL}}(\xi) \approx \frac{1}{\sqrt{2\pi\Delta}} \exp \left[ - \frac{(\sqrt{\Delta} - 1)^2}{2\bar{\sigma}^2} \right], \quad (8) $$

where we have defined $\xi = \delta/\sigma_G$ and $\bar{\sigma} = 6f_{\text{NL}}\sigma_G/5$. The maximum of this distribution is at $2\sqrt{\Delta_{\text{max}}} = 1 + \sqrt{1 - 4\bar{\sigma}^2}$, hence

$$ \xi_{\text{max}} = \frac{\Delta_{\text{max}} - 1 - \bar{\sigma}^2}{2\bar{\sigma}}. \quad (9) $$

This allows one to fix the $f_{\text{NL}}$ parameter such that the $f_{\text{NL}}$ distribution peaks at a given value of $\xi$. In Fig. 1 of the main text, we apply this procedure and display

several $f_{\text{NL}}$ distributions (corresponding to several values of $\sigma_G$) that share the same maximum location with the elliptic profiles. In order to accommodate the heavy tail, one can see that a large value of $f_{\text{NL}}\sigma$ needs to be used. However, when increasing $f_{\text{NL}}\sigma_G$, the agreement at $\xi < \xi_{\text{max}}$ becomes worse, and the point where the $f_{\text{NL}}$ distribution diverges gets dangerously too close to the maximum of the distribution. From this, one concludes that there is no $f_{\text{NL}}$ distribution that provides a reliable approximation of the elliptic profile both close to its maximum and along its tail.

In order to fix $\sigma_G$, one can further expand Eq. (8) around $\xi = \xi_{\text{max}}$, and one finds

$$ P_{\text{NL}}(\xi) = P_{\text{NL}}(\xi_{\text{max}}) \left[ 1 - \frac{(\xi - \xi_{\text{max}})^2}{2\sigma_{\text{max}}^2} + \dots \right] \quad (10) $$

where

$$ \sigma_{\text{max}} = \frac{\left(1 + \sqrt{1 - 4\bar{\sigma}^2}\right)^2}{2\sqrt{2\left(1 - 4\bar{\sigma}^2 + \sqrt{1 - 4\bar{\sigma}^2}\right)}}. \quad (11) $$

This can be used to match the curvature of the PDF around its maximum with a given reference distribution. It leads to the $f_{\text{NL}}$ profile with $f_{\text{NL}}\sigma_G = 0.1$ as plotted in Fig. 1 of the main text.

## II. GAUSSIAN AND LOGNORMAL DISTRIBUTIONS

Similarly to what was done above, the elliptic profiles [see Eq. (1) of the main text] can be approximated by a lognormal (LN) or a Gaussian (G) distribution,

$$ P(x) = A \text{LN}(x, \rho, \sigma) = A e^{-\frac{\sigma^2}{2}} G(x, \rho, \sigma_G), \quad (12) $$

where $x = \pi^2 D / \mu^2$,

$$ \text{LN}(x, \rho, \sigma) = \frac{1}{\rho \sigma \sqrt{2\pi}} \exp \left[ - \frac{\ln(x/\rho)^2}{2\sigma^2} - \frac{\sigma^2}{2} \right] $$

$$ \text{G}(x, \rho, \sigma_G) = \frac{1}{\sigma_G \sqrt{2\pi}} \exp \left[ - \frac{(x - \rho)^2}{2\sigma_G^2} \right] $$

and $\sigma_G = \rho\sigma$. In these parametrizations, $\rho$ stands for the value of $x$ where the PDF is maximal. It is the solution of the equations

$$ \sum_{n=0}^{\infty} (2n + 1)^3 e^{-\rho n(n+1)} \sin \left[ (2n + 1) \frac{\pi}{2} \alpha \right] = 0, \quad (13) $$

for the elliptic 2 distribution, and

$$ \sum_{n=0}^{\infty} (-1)^n n^3 e^{-\rho^2} \sin(n\pi\alpha) = 0 \quad (14) $$
for the elliptic 4 distribution. When $\alpha$ is small, those equations have approximate solutions $\rho_2(\alpha) = \pi^2\alpha^2/6$ for the elliptic 2 distribution and $\rho_4(\alpha) = \pi^2(1-\alpha)^2/6$ for the elliptic 4 distribution (these approximations turn out to be reliable up until $\alpha \simeq 0.6$), otherwise those equations have to be solved numerically. From here, the height of the elliptic distributions at their maxima can be inferred, $P_{2/4}^{\max}(\alpha) = P_{2/4}[\rho_{2/4}(\alpha), \alpha]$, and equated with $Ae^{-\sigma^2/2}/(\sigma_G\sqrt{2\pi})$.

The curvature of the elliptic distributions around their maximum can also be computed according to

$$ \sigma_G^{(2)}(\alpha) = \left[ - \frac{\partial^2 \ln P_2}{\partial x^2} \bigg|_{x=\rho_2(\alpha)} \right]^{-1/2}, \quad (15) $$

$$ \sigma_G^{(4)}(\alpha) = \left[ - \frac{\partial^2 \ln P_4}{\partial x^2} \bigg|_{x=\rho_4(\alpha)} \right]^{-1/2}. $$

In the same small-$\alpha$ limit as above, they boil down to $\sigma_G^{(2)}(\alpha) \simeq \sqrt{2}\alpha^2$ and $\sigma_G^{(2)}(\alpha) \simeq \sqrt{2}(1-\alpha)^2$. This allows one to set the value of $\rho$, $\sigma_G$ (hence $\sigma$) and $A$.

This procedure is performed in Fig. 1 of the main text, where one can check that the agreement around the maximum is indeed excellent, but that the agreement between the lognormal and the elliptic profiles on the tail is also reasonable.

### III. CONNECTING THE PRIMORDIAL SPECTRUM TO THE HALO MASS FUNCTION

The halo mass function determines the number of collapsed halos of a given mass. It is computed from the fraction of collapsed objects $\beta$ as described in Eq. (5) in the main text. We can connect the primordial density perturbations to the halo mass with the following procedure. From a given dimensionless power spectrum $\Delta(k)$ we can compute the real-space variance

$$ \sigma^2(R) = \int \Delta^2(k) W(kR) d\ln k, \quad (16) $$

where $W(r) = 3j_1(r)/r$ is a window function with $j_n$ being the spherical Bessel functions of the first kind. For

spherical collapse the radius can be linked directly to the enclosed mass via $R = (M/M_*)^{1/3}$, with the reference mass scale $M_* = 4\pi\Omega_m\rho_c/3$ and critical energy density $\rho_c = 2.77 \times 10^{11} h^{-1} M_{\odot}/(h^{-1} \text{Mpc})^3$. This allows one to compute $\sigma(M)$ at $z=0$. As described in the main text we can compute the variance at any other redshift including the growth factor, i.e. $\sigma(M, z) = \sigma(M, z=0) \cdot D(z)$.

For Gaussian primordial fluctuations the dimensionless power spectrum today is given by

$$ \Delta_{\text{Gauss}}^2(k) = \frac{k^3}{2\pi^2} A_s \cdot k^n \cdot T^2(k), \quad (17) $$

where $T(k)$ is the transfer function [80] approximated by

$$ T(k) \approx \frac{1}{1 + (k/k_*)^{1.6}} \quad (18) $$

and $A_s = 2 \cdot 10^{-9}$, $n_s = 0.96$ and $k_* = 0.0426/(h^{-1}\text{Mpc})$. The above equations fully determine $\sigma_{\text{Gauss}}(M, z)$.

### IV. CLUSTER COUNTS

The number of halos in a given mass range as a function of redshift is an interesting cosmological observable. From the halo mass function we can compute the number of clusters with masses above a certain threshold, $M_{\text{thr}}$, in a given redshift range, $[z, z + \Delta z]$, as

$$ n_{\text{clusters}} = \int_{z}^{z+\Delta z} \int_{M_{\text{thr}}}^{\infty} \frac{\mathrm{d}n}{\mathrm{d}\ln M} \mathrm{d}\ln M \mathrm{d}z. \quad (19) $$

Then, the number of clusters per comoving volume is

$$ N_{\text{clusters}} = n_{\text{clusters}} \cdot \Delta V_c \quad (20) $$

where

$$ \Delta V_c = \int_z^{z+\Delta z} \frac{4\pi d_L^2 d_H}{(1+z)^2 E(z)} dz \quad (21) $$

with $d_L = (1+z)d_H \int dz/E(z)$ being the luminosity distance, $d_H = c/H_0$ the Hubble distance and $E(z) = H(z)/H_0$ the dimensionless rate of expansion.

[1] W. L. Freedman, Nature Astron. 1, 0121 (2017), 1706.02739.

[2] L. Verde, T. Treu, and A. G. Riess, Nature Astron. 3, 891 (2019), 1907.10625.

[3] F. Finelli, J. García-Bellido, A. Kovács, F. Paci, and I. Szapudi, Mon. Not. Roy. Astron. Soc. 455, 1246 (2016), 1405.1555.

[4] A. Kovács and J. García-Bellido, Mon. Not. Roy. Astron. Soc. 462, 1882 (2016), 1511.09008.

[5] A. Kovács et al., Mon. Not. Roy. Astron. Soc. 510, 216 (2022), 2112.07699.

[6] E. Asencio, I. Banik, and P. Kroupa, Mon. Not. Roy. Astron. Soc. 500, 5249 (2020), 2012.03950.

[7] J. Kim, M. J. Jee, J. P. Hughes, M. Yoon, K. HyeongHan, F. Menanteau, C. Sifon, L. Hovey, and P. Arunachalam, Astrophys. J. 923, 101 (2021), 2106.00031.

[8] A. H. Gonzalez, S. A. Stanford, M. Brodwin, C. Fedeli, A. Dey, P. R. M. Eisenhardt, C. Mancone, D. Stern, and
G. Zeimann, Astrophys. J. **753**, 163 (2012), 1205.3788.

[9] S. L. Finkelstein et al., Nature **502**, 524 (2013), 1310.6031.

[10] A. Del Popolo and M. Le Delliou, Galaxies **5**, 17 (2017), 1606.07790.

[11] O. Newton, M. Cautun, A. Jenkins, C. S. Frenk, and J. Helly, Mon. Not. Roy. Astron. Soc. **479**, 2853 (2018), 1708.04247.

[12] J. Zavala and C. S. Frenk, Galaxies **7**, 81 (2019), 1907.11775.

[13] R. H. Brandenberger, Int. J. Mod. Phys. Conf. Ser. **01**, 67 (2011), 0902.4731.

[14] Y. Akrami et al. (Planck), Astron. Astrophys. **641**, A9 (2020), 1905.05697.

[15] D. Wands, K. A. Malik, D. H. Lyth, and A. R. Liddle, Phys. Rev. D **62**, 043527 (2000), astro-ph/0003278.

[16] I. M. Khalatnikov and A. Y. Kamenshchik, Class. Quant. Grav. **19**, 3845 (2002), gr-qc/0204045.

[17] D. H. Lyth and D. Wands, Phys. Rev. D **68**, 103515 (2003), astro-ph/0306498.

[18] D. H. Lyth, K. A. Malik, and M. Sasaki, JCAP **05**, 004 (2005), astro-ph/0411220.

[19] C. Pattison, V. Vennin, H. Assadullahi, and D. Wands, JCAP **1907**, 031 (2019), 1905.06300.

[20] D. Artigas, J. Grain, and V. Vennin, JCAP **02**, 001 (2022), 2110.11720.

[21] A. A. Starobinsky, Phys. Lett. **117B**, 175 (1982).

[22] A. A. Starobinsky, Lect. Notes Phys. **246**, 107 (1986).

[23] M. Sasaki and E. D. Stewart, Prog. Theor. Phys. **95**, 71 (1996), astro-ph/9507001.

[24] M. Sasaki and T. Tanaka, Prog. Theor. Phys. **99**, 763 (1998), gr-qc/9801017.

[25] D. H. Lyth and Y. Rodriguez, Phys. Rev. Lett. **95**, 121302 (2005), astro-ph/0504045.

[26] K. Enqvist, S. Nurmi, D. Podolsky, and G. I. Rigopoulos, JCAP **0804**, 025 (2008), 0802.0395.

[27] T. Fujita, M. Kawasaki, Y. Tada, and T. Takesako, JCAP **1312**, 036 (2013), 1308.4754.

[28] T. Fujita, M. Kawasaki, and Y. Tada, JCAP **1410**, 030 (2014), 1405.2187.

[29] V. Vennin and A. A. Starobinsky, Eur. Phys. J. **C75**, 413 (2015), 1506.04732.

[30] V. Vennin, H. Assadullahi, H. Firouzjahi, M. Noorbala, and D. Wands, Phys. Rev. Lett. **118**, 031301 (2017), 1604.06017.

[31] H. Firouzjahi, A. Nassiri-Rad, and M. Noorbala, JCAP **01**, 040 (2019), 1811.02175.

[32] Y. Tada and V. Vennin, JCAP **02**, 021 (2022), 2111.15280.

[33] C. Pattison, V. Vennin, H. Assadullahi, and D. Wands, JCAP **1710**, 046 (2017), 1707.00537.

[34] J. M. Ezquiaga, J. García-Bellido, and V. Vennin, JCAP **03**, 029 (2020), 1912.05399.

[35] V. Vennin, Habilitation thesis (2020), 2009.08715.

[36] D. G. Figueroa, S. Raatikainen, S. Rasanen, and E. Tomberg, Phys. Rev. Lett. **127**, 101302 (2021), 2012.06551.

[37] K. Ando and V. Vennin, JCAP **04**, 057 (2021), 2012.02031.

[38] C. Pattison, V. Vennin, D. Wands, and H. Assadullahi (2021), 2101.05741.

[39] G. Rigopoulos and A. Wilkins, JCAP **12**, 027 (2021), 2107.05317.

[40] G. Panagopoulos and E. Silverstein (2019), 1906.02827.

[41] A. Achucarro, S. Cespedes, A.-C. Davis, and G. A. Palma, JHEP **05**, 052 (2022), 2112.14712.

[42] F. Kuhnel and D. J. Schwarz (2021), 2101.10340.

[43] Y.-F. Cai, X.-H. Ma, M. Sasaki, D.-G. Wang, and Z. Zhou (2022), 2207.11910.

[44] S. Clesse and J. García-Bellido, Phys. Rev. **D92**, 023524 (2015), 1501.07565.

[45] M. Kawasaki and Y. Tada, JCAP **1608**, 041 (2016), 1512.03515.

[46] J. M. Ezquiaga and J. García-Bellido, JCAP **1808**, 018 (2018), 1805.06731.

[47] M. Biagetti, G. Franciolini, A. Kehagias, and A. Riotto, JCAP **07**, 032 (2018), 1804.07124.

[48] D. G. Figueroa, S. Raatikainen, S. Rasanen, and E. Tomberg (2021), 2111.07437.

[49] N. Kitajima, Y. Tada, S. Yokoyama, and C.-M. Yoo, JCAP **10**, 053 (2021), 2109.00791.

[50] P. Coles and B. Jones, Mon. Not. Roy. Astron. Soc. **248**, 1 (1991).

[51] A. Taruya, M. Takada, T. Hamana, I. Kayo, and T. Futamase, Astrophys. J. **571**, 638 (2002), astro-ph/0202090.

[52] S. Hilbert, J. Hartlap, and P. Schneider, Astron. Astrophys. **536**, A85 (2011), 1105.3980.

[53] I. Musco, Phys. Rev. D100, 123524 (2019), 1809.02127.

[54] M. Biagetti, V. De Luca, G. Franciolini, A. Kehagias, and A. Riotto (2021), 2105.07810.

[55] V. De Luca and A. Riotto (2022), 2201.09008.

[56] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, *The NIST Handbook of Mathematical Functions* (Cambridge Univ. Press, 2010), ISBN 9780521140638 0521140633 9780521192255 0521192250.

[57] W. H. Press and P. Schechter, Astrophys. J. **187**, 425 (1974).

[58] T. Harada, C.-M. Yoo, and K. Kohri, Phys. Rev. D **88**, 084051 (2013), [Erratum: Phys.Rev.D 89, 029903 (2014)], 1309.4201.

[59] S. Young, C. T. Byrnes, and M. Sasaki, JCAP **07**, 045 (2014), 1405.7023.

[60] C. Germani and I. Musco, Phys. Rev. Lett. **122**, 141302 (2019), 1805.04087.

[61] C.-M. Yoo, T. Harada, J. Garriga, and K. Kohri, PTEP **2018**, 123E01 (2018), 1805.03946.

[62] P. J. E. Peebles, *The large-scale structure of the universe* (1980).

[63] S. Matarrese, L. Verde, and R. Jimenez, Astrophys. J. **541**, 10 (2000), astro-ph/0001366.

[64] N. Dalal, O. Dore, D. Huterer, and A. Shirokov, Phys. Rev. D **77**, 123514 (2008), 0710.4560.

[65] M. LoVerde and K. M. Smith, JCAP **08**, 003 (2011), 1102.1439.

[66] S. Yokoyama, N. Sugiyama, S. Zaroubi, and J. Silk, Mon. Not. Roy. Astron. Soc. **417**, 1074 (2011), 1103.2586.

[67] A. Moradinezhad Dizgah, M. Biagetti, E. Sefusatti, V. Desjacques, and J. Noreña, JCAP **05**, 015 (2021), 2010.14523.

[68] A. Bueno belloso, J. Garcia-Bellido, and D. Sapone, JCAP **10**, 010 (2011), 1105.4825.

[69] A.-C. Eilers, J. F. Hennawi, F. B. Davies, and R. A. Simcoe (2021), 2106.04586.

[70] P. A. R. Ade et al. (Planck), Astron. Astrophys. **594**, A27 (2016), 1502.01598.

[71] L. E. Bleem et al. (SPT, DES), Astrophys. J. Suppl. **247**, 25 (2020), 1910.04121.

[72] M. Hilton et al. (ACT, DES), Astrophys. J. Suppl. **253**,
3 (2021), 2009.11043.

[73] N. Sabti, J. B. Muñoz, and D. Blas, JCAP 01, 010 (2021), 2009.01245.

[74] D. Gilman, A. Benson, J. Bovy, S. Birrer, T. Treu, and A. Nierenberg, Mon. Not. Roy. Astron. Soc. 512, 3163 (2022), 2112.03293.

[75] J. P. Gardner et al., Space Sci. Rev. 123, 485 (2006), astro-ph/0606175.

[76] L. Amendola et al., Living Rev. Rel. 21, 2 (2018), 1606.00180.

[77] v. Ivezić et al. (LSST), Astrophys. J. 873, 111 (2019), 0805.2366.

[78] A. Agrawal, R. Makiya, C.-T. Chiang, D. Jeong, S. Saito, and E. Komatsu, JCAP 10, 003 (2017), 1706.09195.

[79] C. Ramírez-Pérez, J. Sanchez, D. Alonso, and A. Font-Ribera (2021), 2111.05069.

[80] J. Lesgourgues (CLASS) (2011), arXiv:1104.2932.