A non-perturbative test of consistency relations and their violation

TL;DR

This study uses N-body simulations to non-perturbatively verify large scale structure consistency relations derived from symmetry principles, and explores their violation with non-Gaussian initial conditions, aiding future observational constraints.

Contribution

It provides the first non-perturbative test of the consistency relations in the non-linear regime and demonstrates their violation with non-Gaussian initial conditions.

Findings

Consistency relations hold in the non-linear regime for Gaussian initial conditions.

Violations occur when initial conditions are non-Gaussian, specifically of the local fNL type.

The methodology enables constraining primordial non-Gaussianity using large scale structure data.

Abstract

In this paper, we verify the large scale structure consistency relations using N-body simulations, including modes in the highly non-linear regime. These relations (pointed out by Kehagias & Riotto and Peloso & Pietroni) follow from the symmetry of the dynamics under a shift of the Newtonian potential by a constant and a linear gradient, and predict the absence of certain poles in the ratio between the (equal time) squeezed bispectrum and power spectrum. The consistency relations, as symmetry statements, are exact, but have not been previously checked beyond the perturbative regime. Our test using N-body simulations not only offers a non-perturbative check, but also serves as a warm-up exercise for applications to observational data. A number of subtleties arise when taking the squeezed limit of the bispectrum--we show how to circumvent or address them. An interesting by-product of our investigation is an explicit demonstration that the linear-gradient symmetry is unaffected by the periodic boundary condition of the simulations. Lastly, we verify using simulations that the consistency relations are violated when the initial conditions are non-gaussian (of the local fNL type). The methodology developed here paves the way for constraining primordial non-gaussianity using large scale structure data, including (numerous) highly non-linear modes that are otherwise hard to interpret and utilize.

Tables

Table 1. Table 1: Detailed results of the likelihood fits with different sets of parameters ( a ¯ n subscript ¯ 𝑎 𝑛 \bar{a}_{n} ’s). The Bayesian information criterion clearly selects the ( a ¯ 0 , a ¯ 1 ) subscript ¯ 𝑎 0 subscript ¯ 𝑎 1 (\bar{a}_{0},\bar{a}_{1}) model.

${\overline{a}}_n$ included	${\overline{a}}_{-2}\left({10}^{-6}\text{Mpc}/h\right)$	${\overline{a}}_{-1}{\left({10}^{-2}\text{Mpc}/h\right)}^2$	BIC
${\overline{a}}_0$	$-$	$-$	99.63
${\overline{a}}_0,{\overline{a}}_1$	$-$	$-$	17.77
${\overline{a}}_0,{\overline{a}}_1,{\overline{a}}_2$	$-$	$-$	19.82
${\overline{a}}_{-2},{\overline{a}}_0$	$-30.4\pm 5.1$	$-$	65.50
${\overline{a}}_{-2},{\overline{a}}_0,{\overline{a}}_1$	$0.2\pm 6.7$	$-$	19.84
${\overline{a}}_{-1},{\overline{a}}_0$	$-$	$-42.3\pm 5.4$	39.57
${\overline{a}}_{-1},{\overline{a}}_0,{\overline{a}}_1$	$-$	$0.6\pm 10.3$	19.84
${\overline{a}}_{-2},{\overline{a}}_{-1},{\overline{a}}_0$	$69\pm 16$	$111\pm 17$	22.35
${\overline{a}}_{-2},{\overline{a}}_{-1},{\overline{a}}_0,{\overline{a}}_1$	$16\pm 55$	$26\pm 87$	21.83

Equations

$$\Phi \to \Phi +c$$

$$\vec{x}\to \vec{x}+\vec{n},\Phi \to \Phi -(\mathcal{H}{\vec{n}}^{\prime }+{\vec{n}}^{\prime \prime })\cdot \vec{x},\vec{v}\to \vec{v}+{\vec{n}}^{\prime }$$

$$\underset{\vec{q}\to 0}{\lim }{q}^2\frac{⟨{\delta }_{\vec{q}}{\delta }_{{\vec{k}}_1}\cdots {\delta }_{{\vec{k}}_N}{⟩}^{c^{\prime }}}{P_{\delta }(q)}=0,$$

$$\begin{array}{rl}& \underset{\vec{q}\to 0}{\lim }{\vec{\nabla }}_q\left[q^2\frac{⟨{\delta }_{\vec{q}}{\delta }_{{\vec{k}}_1}\cdots {\delta }_{{\vec{k}}_N}{⟩}^{c^{\prime }}}{P_{\delta }(q)}\right]\\ & =-\sum _{a=1}^N\frac{D({\eta }_a)}{D(\eta )}{\vec{k}}_a⟨{\delta }_{{\vec{k}}_1}\cdots {\delta }_{{\vec{k}}_N}{⟩}^{c^{\prime }},\end{array}$$

$$\frac{⟨{\delta }_{\vec{q}}{\delta }_{{\vec{k}}_1}\dots {\delta }_{{\vec{k}}_N}{⟩}^{c^{\prime }}}{P_{\delta }(q)}\mathrm{has}\mathrm{no}{q}^{-2}\mathrm{pole}\&\mathrm{no}{q}^{-1}\mathrm{pole}$$

$$\frac{⟨{\delta }^g{}_{\vec{q}}{\delta }^g{}_{{\vec{k}}_1}\dots {\delta }^g{}_{{\vec{k}}_N}{⟩}^{c^{\prime }}}{P_{{\delta }^g}(q)}\mathrm{has}\mathrm{no}{q}^{-2}\mathrm{pole}\&\mathrm{no}{q}^{-1}\mathrm{pole}$$

$$B_{\delta }(q,k,\theta )=\sum _{n=-2}^{\infty }{a}_n(k,\theta )P_{\delta }(q)q^n.$$

$${\bar{B}}_{\delta }(q)$$

$$=\sum _{n=-2}^{\infty }\frac{{\sum }_i{a}_n(k_i,{\theta }_i)N_{i,q}}{{\sum }_i{N}_{i,q}}{P}_{\delta }(q)q^n$$

$$\equiv \sum _{n=-2}^{\infty }{\bar{a}}_n(q)P_{\delta }(q)q^n.$$

$${\bar{B}}_{\delta }(q)=\frac{{\sum }_i{B}_{\delta }(k_i,{\theta }_i,q)}{{\sum }_i}\equiv \sum _{n=-2}^{\infty }{\bar{a}}_n{P}_{\delta }(q)q^n.$$

$${\bar{a}}_n\equiv \frac{{\sum }_i{a}_n(k_i,{\theta }_i)}{{\sum }_i},$$

$${\bar{B}}_{\delta }(\bar{q})=\sum _{q\in \{\bar{q}\}}{\bar{B}}_{\delta }(q)=\sum _{n=-2}^{n_{\text{max}}}{\bar{a}}_n{M}_n(\bar{q}),$$

$${\mathcal{L}}_{(r)}\propto \frac{1}{\sqrt{\det \mathbf{C}}}\exp \left[-\frac{1}{2}{\mathbf{\Delta }}^{(r)}\cdot {\mathbf{C}}^{-1}\cdot {\mathbf{\Delta }}^{(r)}\right],$$

$${\bar{a}}_n$$

$${\sigma }_{{\bar{a}}_n}^2$$

$$\text{BIC}=-2\log {\mathcal{L}}_{\text{max}}+N_{\text{par}}\log N_q,$$

$$\begin{array}{rl}{a}_0^{\text{pt}}(k,\theta )& =\left(\frac{13}{7}+\frac{8}{7}{\cos }^2\theta \right)P_{\delta }(k)-{\cos }^2\theta k{P}_{\delta }^{\prime }(k).\end{array}$$

Full text

I Introduction

One of the key questions in modern cosmology concerns the initial condition of the universe. Are the primordial fluctuations consistent with what one would expect from single-field inflation? Or do they arise from a scenario in which additional light fields, besides the inflaton, play an important role? Or more radically, is some mechanism other than inflation at work?

The standard approach to answering these questions is to work with data in the linear or quasi-linear regime where perturbation theory can be relied upon to give reliable predictions. Modes in the non-linear regime (for instance, with momentum $k\,\gtrsim\,0.2$ $h$/Mpc in large scale structure data) are not utilized, even though they are abundant and measured with high precision.

The consistency relations offer an interesting alternative, where some of the information hidden in the nonlinear regime can be brought to light. First pointed out by Maldacena Maldacena (2003), consistency relations connect a squeezed $(N+1)$-point correlation function (squeezing means one of the momentum legs is soft) to an $N$-point function (see also Creminelli and Zaldarriaga (2004); Cheung et al. (2008)). More recent work pointed out additional consistency relations coming from new symmetries, clarified the assumptions behind consistency relations and emphasized their exact, non-perturbative nature, analogous to soft theorems in high energy physics Creminelli et al. (2012); Hinterbichler et al. (2012); Assassi et al. (2012); Kehagias and Riotto (2012); Hinterbichler et al. (2014); Goldberger et al. (2013); Hui et al. (2019). The non-perturbative nature of consistency relations is a mere curiosity for the microwave background since its fluctuations are small and linear, but becomes very interesting for large scale structure. Kehagias & Riotto and Peloso & Pietroni Kehagias and Riotto (2013); Peloso and Pietroni (2013) pointed out the relevant large scale structure consistency relations. It can be shown that of the infinite tower of general relativistic consistency relations Hinterbichler et al. (2014), two has non-trivial Newtonian, sub-Hubble limits Creminelli et al. (2013); Horn et al. (2014).

The study of large scale structure concerns, at a minimum, the following quantities: the mass fluctuation $\delta$, the peculiar velocity $\vec{v}$ and the gravitational potential $\Phi$. (One can further expand this list to include the galaxy count fluctuation $\delta^{g}$ and the galaxy peculiar velocity $\vec{v}^{\,g}$. The symmetries discussed below apply to them as well, where $\delta^{g}$ and $\vec{v}\,{}^{g}$ transform in the same way as $\delta$ and $\vec{v}$ do, see e.g. Kehagias et al. (2014a, 2015)) The dynamics of (sub-Hubble) fluctuations exhibits two non-linearly realized symmetries in a matter + cosmological constant universe.111The split into two separate symmetries here follows the discussion of Horn et al. (2015). One is a constant shift in the gravitational potential:

[TABLE]

where $c$ is independent of space but possibly a function of time. The other involves adding a linear gradient to the gravitational potential, together with a transformation of the spatial coordinates and the velocity field Horn et al. (2014):

[TABLE]

where $\vec{n}$ is independent of space but a function of time. Here ${}^{\prime}\equiv\partial/\partial\eta$ is the derivative with respect to the conformal time $\eta$, and $\mathcal{H}\equiv a^{\prime}/a$ is the comoving Hubble parameter, with $a$ being the scale factor. The above is a symmetry of the large scale structure dynamics for $\vec{n}$ having any time dependence, but the adiabatic mode condition Weinberg (2003) dictates that $\vec{n}$ must match the time-dependence of the linear growth factor, and likewise $c$ should match the corresponding time-dependence of the gravitational potential (see discussions in Horn et al. (2015); Hui et al. (2019) and point 2 below).

The consistency relations corresponding to a shift of the gravitational potential by a constant and by a linear gradient are respectively:

[TABLE]

and

[TABLE]

where $P_{\delta}(q)$ is the mass power spectrum, $D$ is the linear growth factor and $c^{\prime}$ denotes the connected correlator with the overall $\delta$-function removed. The time dependence is as follows: the soft mode $\vec{q}$ is at time $\eta$ (likewise for $P_{\delta}(q)$) while the hard mode $\vec{k}_{a}$ is at time $\eta_{a}$. Several comments are called for on these two consistency relations.

1. The consistency relations are in general of an unequal time form. In this paper, we focus on the equal time limit, in which case the right hand side of the Eq. (4) vanishes. Thus, the content of the consistency relations is simple, that the equal time correlator

[TABLE]

in the $\vec{q}\rightarrow 0$ limit. The lack of a $1/q^{2}$ pole follows from the shift symmetry, and the lack of a $1/q$ pole follows from the linear gradient symmetry. That this statement is correct (for gaussian initial conditions) is easy to check in perturbation theory (see e.g. Peloso and Pietroni (2013); Kehagias and Riotto (2013); Horn et al. (2014)). But the consistency relations, as symmetry statements, are expected to be stronger than this. What we wish to accomplish in this paper is to test this statement in the non-perturbative regime using $N$-body simulations (i.e. with the hard momenta ${\vec{k}_{a}}$’s on nonlinear scales).222We focus on the equal time correlator largely for simplicity. There is also a practical reason for doing so: that the $1/q$ pole associated with the unequal time contributions (i.e. the right hand side of Eq. (4)) is naturally suppressed in observational data. Recall that the unequal times refer to the times of the hard modes ($\eta_{1}$ for $\vec{k}_{1}$, $\eta_{2}$ for $\vec{k}_{2}$ and so on); the hard modes are by definition short wavelength perturbations which also means their separation in time cannot be too big—keep mind that observational data are confined to the light cone. One can see from Eq. (4) that if the $\eta_{a}$’s are close to each other, one is almost summing the $\vec{k}_{a}$’s which yields zero. Nonetheless, it is worth asking how big of a $1/q$ pole one might inadvertently generate by measuring a correlator averaged over some survey volume, which inevitably spans a range of redshifts. Some care in defining the average might be useful to ensure it is negligible. It is also worth noting that the unequal-time contributions do not generate a $1/q^{2}$ pole. A $1/q^{2}$ pole can only appear with certain primordial non-gaussianities (see point 6 below).

2. It should be emphasized that the consistency relations are not statements merely about a strictly vanishing $\vec{q}$. Indeed, an exact $\vec{q}=0$ mode is not even observable. Rather, the consistency relations are statements about the absence of certain divergences as $\vec{q}$ is taken to be smaller and smaller, such as (5). This is why the so called adiabatic mode condition is crucial Weinberg (2003); Mirbabayi and Simonović (2016). This condition ensures that the symmetry in question, which in general originates as a gauge redundancy, generates a $\vec{q}=0$ mode that is smoothly connected to a physical mode of a small but finite $\vec{q}$. For more discussions on this point, see Hui et al. (2019).

3. The consistency relations Eqs. (3) and (4) take a particularly simple form in Lagrangian space where the corresponding “right hand side” vanishes even if the hard modes are at unequal times. See Horn et al. (2015) for a discussion.

4. The consistency relations take essentially the same form even in redshift space, as pointed out by Creminelli et al. (2014a). This means they can be profitably applied in galaxy surveys where the line-of-sight direction is almost always in redshift space.

5. There is the question of how galaxy biasing affects the consistency relations. As mentioned above, the relevant symmetries remain good symmetries even for the dynamics of galaxies (which can form, merge and so on).333Galaxy dynamics is of course different from mass dynamics: mass conservation is replaced by galaxy number density evolution that has a source (or sink) term; galaxies are subject to forces beyond gravity. The key observation is that as long as these new terms/forces depend only on mass/galaxy density, velocity gradients (or velocity difference between different species) and second derivatives of the gravitational potential (tidal forces), the symmetries espoused in Eqs. (1) and (2) hold. For instance, it is crucial the new forces on a galaxy do not depend on the absolute velocity, i.e. some form of equivalence principle (see point 6 below). There’s an additional requirement: that the squeezed momentum $\vec{q}$ must be sufficiently soft, that on that scale, gravity dominates (even though for the hard momenta $\vec{k}$’s, the dynamics can be complicated). See Peloso and Pietroni (2014); Horn et al. (2014) for further discussions.

Thus, the consistency relations Eqs. (3) and (4) remain valid even if the hard modes $\delta_{\vec{k}_{1}}\,,\dots,\,\delta_{\vec{k}_{N}}$ are replaced by galaxy density fluctuations $\delta^{g}{}_{\vec{k}_{1}}\,,\dots,\,\delta^{g}{}_{\vec{k}_{N}}$. The soft mode $\delta_{\vec{q}}$ can be replaced by $\delta^{g}{}_{\vec{q}}/b^{g}$ where $b^{g}$ is the galaxy bias; likewise $P_{\delta}(q)$ can be replaced by $P_{\delta^{g}}(q)/b^{g}{}^{2}$. In the soft limit, $b^{g}$ is expected to be a constant444This holds if the initial conditions were gaussian, an assumption that goes into the derivation of the consistency relations themselves. Or more precisely, this assumes single-field or single-clock initial conditions. See discussion in point 6 below., and thus the consistency relations Eqs. (3) and (4) are modified in a simple way. The equal time version (5) in fact takes the same form i.e. the equal time correlator:

[TABLE]

in the $\vec{q}\rightarrow 0$ limit.

6. Two important assumptions go into deriving the consistency relations. One is the equivalence principle, that on sufficiently large scales—i.e. $\vec{q}\rightarrow 0$—all objects fall at the same rate (whereas on small scales, different objects can be subject to different forces, such as pressure forces, etc). See Creminelli et al. (2014b); Horn et al. (2014) for a discussion. The other important assumption, which we focus on in this paper, is gaussian initial conditions. More precisely, it is the assumption that in the squeezed limit, the primordial connected $N$-point function vanishes for $N>2$, something that follows from single-clock inflation.555 The primordial consistency relations can be expressed as the vanishing of the squeezed $N$-point function if one accounts for the fact that the metric fluctuations enter into the definition of physical momenta. See Tanaka and Urakawa (2011); Pajer et al. (2013); Bravo et al. (2018a, b) for a discussion.

From the point of view of initiating cosmological $N$-body simulations, imposing gaussian initial conditions is sufficient to guarantee the validity of the consistency relations stated above, and this is what we adopt in this paper. It is not surprising that the consistency relations, or the precise form they take, are sensitive to initial conditions, since the symmetries underlying them are non-linearly realized or spontaneously broken—in other words, exactly how the initial conditions, or the “vacuum”, breaks the symmetries in question dictates the form of the consistency relations Hui et al. (2019). Examples that violate the stated consistency relations generally involve extra light fields during inflation, for instance the curvaton, a spectator scalar that dominates the curvature fluctuations Lyth (2006); Dvali et al. (2004); Kofman (2003).666Ultra-slow-roll inflation, while strictly a single field model, has essentially an extra clock due to the importance of what normally would be discarded as the decaying mode. See Namjoo et al. (2013); Martin et al. (2013); Finelli et al. (2018); Hui et al. (2019). The curvaton (or modulated-reheating) model motivates initial conditions of the local $f_{\rm NL}$ type (see §II.4), and we will examine how the consistency relations are violated in such a case. The ultimate goal would be to check consistency relations in observational data, and put a bound on local $f_{\rm NL}$ for instance. The robustness of the consistency relations means we can freely employ data in the highly nonlinear regime (the high momentum $\vec{k}$ modes), involving astrophysically realistic fluctuations, e.g. galaxies.

7. One might worry that the consistency relation could be violated by the finite size of the simulation box, especially for the symmetry transformation that involves shifting the gravitational potential by a term linear in $\vec{x}$ (Eq. (2)), which seems naïvely inconsistent with the periodic boundary condition of the simulations. However, from the point of view of the particles, all they see is the gradient of the potential, and the symmetry in question simply shifts this gradient by a constant, which does respect the periodic boundary condition. The fact that, as we will see, the consistency relations hold in the $N$-body data indeed confirms this expectation.

8. Lastly, it should be kept in mind that in the presence of features in the power spectrum (e.g. acoustic peaks), the bispectrum could present a $1/q$ behavior for mildly squeezed triangles, albeit recovering the behavior expected from Eq. (5) in the strict $q\to 0$ limit.777We are grateful to Marko Simonović for pointing this out. When such features are present a simple power series in $q$ will not suffice to describe the squeezed bispectrum, and the complete dependence should be taken into account Mirbabayi et al. (2014); Baldauf et al. (2015). For the range of $k$’s we are considering here, this is a negligible effect any way.

To summarize, the goal of this paper is twofold. First, we test the consistency relations (5) at equal time using the results of $N$-body simulations with gaussian initial conditions, focusing on the three-point function or bispectrum. To the best of our knowledge, this is the first time that the consistency relations have been verified for scales that are well within the non-perturbative regime.888 A different kind of consistency relation has been tested in Nishimichi and Valageas (2014) using $N$-body simulations as well. That interesting relation concerns the higher order coefficients of the low $q$ expansion of the bispectrum Kehagias et al. (2014b); Valageas (2014). More specifically, it concerns the $q^{0}$ behavior in the context of (5), and its derivation crucially rely on the hard observables being mass fluctuations as opposed to galaxy fluctuations. The consistency relations we focus on are instead more robust and valid even for galaxy observables. Secondly, we show that when the initial conditions for the primordial fields are non-gaussian of the local $f_{\rm NL}$ type, deviations from (5) are observed, as expected from theoretical arguments Peloso and Pietroni (2013); Horn et al. (2014); Valageas et al. (2017).

II Checking the consistency relations in $N$-body simulations

We describe in §II.1 our methodology, focusing in particular on how to obtain the bispectrum in the squeezed limit. This is followed by a discussion in §II.2 of how we fit the bispectrum with a power series in the squeezed momentum $q$. The results of the fit are presented in §II.3, for $N$-body simulations with gaussian initial conditions. We verify that the consistency relations are indeed satisfied, even though the high momentum modes are in the non-linear regime. We draw attention to, and comment on, the fact that the linear-gradient consistency relation (i.e. the lack of $1/q$ pole in (5)) is satisfied, despite the periodic boundary conditions of the simulations—which one might naïvely expect to invalidate the linear gradient symmetry of Eq. (2). We demonstrate in §II.4 that the consistency relations are violated for simulations with non-gaussian initial conditions of the local $f_{\rm NL}$ type.

II.1 Setup and details of the measurement

We use a suite of $N$-body simulations consisting of $N_{r}=40$ realizations with gaussian initial conditions. The box size is $L=2.4$ Gpc/$h$ comoving, with $1280^{3}$ particles. The cosmological parameters are $\Omega_{\Lambda}=0.75$, $\Omega_{m}=0.25$ (of which $\Omega_{b}=0.04$), $h=0.7$, $n_{s}=1$ and $\sigma_{8}=0.8$. We analyze the simulation outputs at redshift $z=0$. For further details on the simulations, see Scoccimarro et al. (2012).

A prime observable of focus is the bispectrum, in the so-called squeezed limit, i.e. when one of the legs (in momentum space) is soft. We are particularly interested in what happens when that leg, labeled by the momentum $q$, becomes softer and softer as other quantities that label the relevant momentum-space triangle are kept fixed. A convenient parametrization is to take them to be the highest momentum leg, labeled by $k$, and its angular separation from the soft leg, labeled by $\theta$. With this choice, $\theta$ is between $\pi/2$ and $\pi$ (see Fig. 1).999By restricting ourselves to $k$ being the highest momentum and $\theta$ between $\pi/2$ and $\pi$, we are implicitly assuming parity invariance: that two triangles related to each other by a reflection have the same bispectrum.

We will have more to say about the choice of parametrization below.

The bispectrum in the squeezed limit can then be expressed as a power series in the soft mode,

[TABLE]

We will truncate this power series at some finite $n$, with the understanding that this is a good approximation for small values of $q$—the precise $n$ at which we truncate will be determined by the goodness-of-fit to the data. The consistency relations (5) tell us that, for equal time correlators, one has $a_{-2}(k,\theta)=a_{-1}(k,\theta)=0$. The goal of this paper is to check this prediction. We wish to do it in a way that does not assume any knowledge of the coefficients $a_{n}$. They are known robustly only within perturbation theory, that is, if $k$ is not too large. For large $k$’s, non-linearity, or baryonic physics in the case of galaxy observables (in anticipation of applications to observational data), makes it difficult to robustly predict $a_{n}$. Thus we carry out the analysis without prior assumptions on them.

The key feature we exploit is that Eq. (7) takes a factorized form: for each $n$, the dependence on the soft momentum $q$ is factorized from the dependence on the hard momentum $k$ (and $\theta$). The coefficients that contain the $k$ and $\theta$ dependence, $a_{n}(k,\theta)$, can be treated as free parameters when fitting the bispectrum. As a simplifying procedure, since we are not ultimately interested in the $k$ and $\theta$ dependence of the bispectrum or $a_{n}$, we average over all possible values of $k$ and $\theta$ when we measure the bispectrum for a given $q$.101010We will later check this procedure by varying the range of $\theta$ over which we average.

At this point, a subtlety occurs because of the discrete nature of the Fourier modes in a finite volume. Let us focus on the coefficient $a_{n}$ for a particular $n$. Our procedure is effectively to compute some averaged version of $a_{n}$ by summing over all possible $k$’s and $\theta$’s at a fixed $q$, i.e. summing over all triangles which has one momentum leg of magnitude $q$. The issue is this: within our set of discrete Fourier or momentum modes, for a given $q$, not all possible $k$’s and $\theta$’s are actually allowed—in fact, the span of possible $k$’s and $\theta$’s would depend on the value of $q$ in a subtle way; this means the averaged $a_{n}$ would end up inheriting a subtle $q$ dependence. This $q$ dependence cannot be predicted without prior knowledge or assumption of how $a_{n}$ depends on $k$ and $\theta$. It is useful to concretely see how this comes about by dividing the $k$’s and $\theta$’s into bins, labeled by $i$. For instance, a bin centered at $(k_{i},\theta_{i})$ might have contributions from $N_{i,q}$ triangles. Note how the $q$ dependence is “sneaked” in through the fact that $N_{i,q}$ depends on $q$. In this language, averaging over all possible triangles for a given $q$ amounts to computing the following:

[TABLE]

We are thus left with an averaged $a_{n}$, which we call $\bar{a}_{n}$, that has an unwanted $q$ dependence which cannot be predicted without making assumptions about how $a_{n}$ behaves for high $k$’s. Thus, imagine we fit the $N$-body data with $n$ up to, for example, $1$. Even if one puts aside the possible $q$ dependence of $\bar{a}_{-2}$ and $\bar{a}_{-1}$ (which for Gaussian initial conditions are expected to vanish), the unknown $q$ dependence of $\bar{a}_{0}$ and $\bar{a}_{1}$ is problematic.

This way of spelling out the problem also suggests its cure. The above averaging weighs each $i$-th bin by the number of triangles in it, $N_{i,q}$. We can instead weigh each bin equally (or for that matter, use any other weights as long as they do not depend on $q$):111111In practice, this means when we loop through the triangles for a given $q$, we weigh them by $1/N_{i,q}$.

[TABLE]

The coefficients $\bar{a}_{n}$ are now given by

[TABLE]

and are truly independent of the soft momentum. They are treated as free parameters in our fit of the data.

To simplify the analysis, we also bin in $q$. In particular, if $\{\bar{q}\}$ is a bin with average soft momentum $\bar{q}$, the binned version of the bispectrum (9) is

[TABLE]

where $M_{n}(\bar{q})=\sum_{q\in\{\bar{q}\}}P_{\delta}(q)q^{n}$ is also measured from the data to avoid any theoretical bias. The value of $n_{\text{max}}$ is something we have to experiment with: qualitatively, the more squeezed our triangles are (smaller $q$’s), the lower is the $n_{\text{max}}$ we need. In the next section we will rely on data to find how many $\bar{a}_{n}$’s we need to account for higher order corrections to the small $q$ expansion.

Lastly, let us comment on our bispectrum triangle parametrization, described in Fig. 1. One alternative Lewis (2011) is to parametrize the triangle in terms of $\vec{q}$, $\vec{k}+\vec{q}/2$ (where $\vec{k}$ is one of the two high momentum legs), and the angle between them, say $\beta$. Assuming invariance under parity (a reflection of the triangle), the bispectrum should be unchanged under ${\,\rm cos\,}\beta\rightarrow{\,\rm cos\,}(\pi-\beta)$. Thus, if $q$ enters into the bispectrum only through ${\,\rm cos\,}\beta$, the squeezed bispectrum should contain only even powers of $q$, as suggested by Lewis (2011). This appears to be true in some cases but not in others—for instance, it can be checked in perturbation theory that if the initial conditions were of the local $f_{\rm NL}$ type, the squeezed mass bispectrum depends on the transfer function at the soft-momentum $q$ (which equals $1$ when $q=0$, but has corrections with both even and odd powers of $q$) Peloso and Pietroni (2013).121212Even in those cases where the squeezed bispectrum parametrized according to Lewis (2011) appears as an even-power series in $q$, the information does not allow us to for instance infer $\bar{a}_{-1}$ from $\bar{a}_{-2}$ (these are parameters we are ultimately interested in, using our parametrization). In those cases, $\bar{a}_{-1}$ is not directly related to $\bar{a}_{-2}$ but is instead related to the average over $k$ and $\theta$ of the derivatives of $a_{-2}$ with respect to $k$ and $\theta$. We thank Antony Lewis for discussions on different parametrizations of the bispectrum triangle.

We thus do not find a significant advantage for using the alternative parametrization, although the bispectrum can be analyzed that way if one wishes.

II.2 Details of the fit

Our analysis is performed averaging the bispectrum as in Eq. (9) over hard modes ranging from $k=0.52$ $h$/Mpc to $k=0.65$ $h$/Mpc (corresponding to $k=199k_{f}$ to $k=251k_{f}$, with $k_{f}=2\pi/L\simeq 2.6\times 10^{-3}$ $h/$Mpc being the fundamental mode), and over all the relative angles, $\theta\in[\pi/2,\pi]$131313As a simple check, we also repeated our analysis averaging over $\theta\in[\pi/2,3\pi/4]$. The conclusions are largely unchanged, suggesting that adjusting the angular weighting does not have a significant impact on the outcome, and that the vanishing of $\bar{a}_{-2}$ and $\bar{a}_{-1}$ for gaussian initial conditions is not the result of accidental cancellation when averaging over angles.. As one can see from Fig. 2 the hard modes we are considering are well within the non-linear regime.

We also measure the bispectrum for soft momenta ranging from $q=3k_{f}$ to $q=19k_{f}$, with a bin size $\Delta k=2k_{f}$. The choice of this window for $q$’s requires some explanations. The high end $q=19k_{f}$ is chosen to include as many modes as possible (thus minimizing error bars on $\bar{a}_{-2}$ and $\bar{a}_{-1}$), while still staying within the squeezed limit such that the expansion in Eq. (9), truncated at $n$ of a few, is a good approximation.

The low end $q=3k_{f}$ is chosen because the procedure of eliminating the unwanted $q$ dependence in $\bar{a}_{n}$, described in Sec. II.1, is actually not perfect. Recall that we form bins, labeled as $(k_{i},\theta_{i})$, and compute $\bar{B}_{\delta}$ and $\bar{a}_{n}$ (Eqs. (9) and (10)) by giving these bins equal weights. In doing so, it is important that each bin is actually not empty, that there are triangles that fall into them. Thus, the bins have to be sufficiently wide. But there is some tension between using wide bins and using the bin-averaged $a_{n}$ as a fair representation of how $a_{n}$ varies with $k$ and $\theta$.141414In other words, within a wide bin, the precise set of $k$’s and $\theta$’s that fall into that bin would depend on $q$, and thus we are not achieving the goal of eliminating the unwanted $q$ dependence. In Appendix  A, we show a test of our procedure for a particular model of $a_{n}(k,\theta)$ (one motivated by perturbation theory), and check to what extent our procedure yields $\bar{a}_{n}$ that is truly $q$ independent. We find that this works well as long as $q\geq 3k_{f}$. Hence we restrict our analysis only to soft modes such that $q\geq 3k_{f}$.

To determine the best-fit values for the parameters $\bar{a}_{n}$ we maximize the following likelihood for each realization, $r=1,\dots,N_{r}$:

[TABLE]

We define the vector $\bm{\Delta}^{(r)}=\bar{\bm{B}}_{\delta}^{(r)}-\sum_{n}\bar{a}^{(r)}_{n}\bar{\bm{M}}^{(r)}_{n}$, and the covariance matrix $C_{ij}=\langle\Delta_{i}\,\Delta_{j}\rangle$. All vectors run over the soft momenta, $f_{i}\equiv f(\bar{q}_{i})$, and the angular brackets stand for an average over the available realizations, i.e. $\langle f_{i}\rangle=\frac{1}{N_{r}-1}\sum_{r}f_{i}^{(r)}$ (for instance, the covariance matrix is obtained by averaging over realizations).

Note that the vector $\bm{\Delta}$ depends on the fit parameters $\bar{a}_{n}$, and so the covariance matrix itself depends on the parameters. We use an iterative procedure (akin to the Newton-Raphson algorithm) to determine the optimal $\bar{a}_{n}$’s that maximize the likelihood. First, we determine $C_{ij}$ with the $\bar{a}_{n}$’s set to zero. The maximization of the likelihood can thus be done analytically, because the remaining dependence on $\bar{a}_{n}$ shows up only in the exponent of the likelihood in the standard $\chi^{2}$ fashion (essentially equivalent to fitting the slope of a straight line). The resulting best-fit $\bar{a}_{n}$’s are plugged back into the definition of $C_{ij}$, and the whole procedure is repeated again to obtain a new set of best-fit $\bar{a}_{n}$’s. So on and so forth until convergence is achieved.

Once this is done, the final value of the ML estimators and their uncertainties is computed from the average and variance over realizations, i.e.

[TABLE]

Note that the likelihood analysis itself, applied to each realization, does yield an error estimate, but we deem $\sigma_{\bar{a}_{n}}^{2}$ estimated from the spread between independent realizations as more reliable. For one thing, the likelihood analysis treats the data vector as Gaussian distributed, which is an approximation. The desire to have an accurate error estimate is why we analyze the realizations one at a time, as opposed to using all of them in one go.151515Analyzing the realizations all at once would give us essentially the same final best-fit $\bar{a}_{n}$, but would not let us reliably estimate the associated errorbar.

To determine the goodness of the fit, we rely on the Bayesian information criterion (BIC) Schwarz et al. (1978); Liddle (2004):

[TABLE]

which has been shown to be dimensionally consistent, i.e. not to favor overfitted models Liddle (2004). Here $\mathcal{L}_{\text{max}}$ is the maximum likelihood combining all realizations, $N_{\text{par}}$ the number of parameters of the model and $N_{q}$ the number of data points used. A model with the lowest BIC represents the best compromise between maximizing likelihood and minimizing the number of parameters.

II.3 Results

Let us now present the results of our analysis. In Fig.  3 we report some of the fits including different sets of parameters as well as the corresponding residues. In Table 1 we compare all the models we have tested.

Focus first on the first three models in the table, which do not involve $\bar{a}_{-2}$ nor $\bar{a}_{-1}$. We see that the first model, involving $\bar{a}_{0}$ alone, is not a good fit to the $N$-body data (from both the BIC value in the table and from Fig. 3). Adding $\bar{a}_{1}$ greatly improves the fit, while further adding $\bar{a}_{2}$ does not lower the BIC score. Recall that in our power-series fit of the squeezed bispectrum as a function of a range of soft momenta (Eq. (11)), we do not know a priori how many higher order terms we need. This exercise tells us it is sufficient to stop at $\bar{a}_{1}$ (but necessary to include it), with the kind of precision and the range of soft momenta we have.

The rest of the models in the table involve $\bar{a}_{-1}$ and/or $\bar{a}_{-2}$. In all cases, the BIC score worsens. The inferred values for $\bar{a}_{-1}$ and $\bar{a}_{-2}$ are consistent with zero, except for the $(\bar{a}_{-2},\bar{a}_{-1},\bar{a}_{0})$ model. For this model, the fit prefers non-zero values for $\bar{a}_{-2}$ and $\bar{a}_{-1}$ to compensate for the lack of a $\bar{a}_{1}$ term. Note however this model has a worse BIC score compared to the $(\bar{a}_{0},\bar{a}_{1})$ model. It is also reassuring that the $(\bar{a}_{0},\bar{a}_{1})$ fits the data well, with residues that have no clear trend with momenta (see the right panel of Fig. 3).

We conclude from this exercise that the $N$-body data, with gaussian initial conditions, are consistent with a vanishing value for $\bar{a}_{-2}$ and for $\bar{a}_{-1}$, confirming expectations from the consistency relations.

II.4 Violation of the consistency relations from non-gaussian

initial conditions

In this section we show that, when the initial conditions for the cosmological fields are non-gaussian (of the local $f_{\rm NL}$ type), statistically significant deviations from the consistency relations in Eq. (5) are observed.161616Specifically, the local $f_{\rm NL}$ model is this: the primordial Bardeen potential $\Phi_{p}(\vec{x})=\phi(\vec{x})+f_{\rm NL}(\phi(\vec{x})^{2}-\langle\phi^{2}\rangle)$, where $\phi$ is a gaussian random field. The Bardeen potential (after multiplication by the transfer function) gives the gravitational potential used in initializing $N$-body simulations (see Scoccimarro et al. (2012) for details). A primordial non-gaussianity of this type is motivated by the curvaton and modulated reheating models Lyth (2006); Dvali et al. (2004); Kofman (2003).

We employ a smaller set of $N_{r}=12$ realizations with the same cosmological parameters as before, but with an initial matter distribution characterized by a local non-gaussian parameter $f_{\text{NL}}=100$. The details of the measurement and the analysis are the same as in Sections II.1 and II.2. For the sake of checking whether or not deviations from consistency relations occur we limit our analysis to models that only include $\bar{a}_{0}$ and $\bar{a}_{1}$ (in addition to possibly $\bar{a}_{-2}$ and $\bar{a}_{-1}$). We leave a more detailed study of the realizations with non-gaussian initial condition for future work Esposito et al. .

The result of our analysis is unambiguous. The model with just $\bar{a}_{0}$ and $\bar{a}_{1}$ (i.e. no poles in the soft limit), which fits very well the bispectrum in the case of gaussian initial conditions, is not a good description of the data obtained from $f_{\text{NL}}=100$. From the fit we obtain $\text{BIC}=74.12$, much larger than the one reported in Table 1 for the gaussian case. Moreover, Fig. 4 shows that the likelihood fit for this model is not a good description of the data, which is also confirmed by the fact that the residues exhibit a parabolic pattern around zero. Indeed, introducing either $\bar{a}_{-2}$ or $\bar{a}_{-1}$ to the fit one obtains values that are statistically different from zero: $\bar{a}_{-2}=(11.3\pm 1.8)\times 10^{-5}$ Mpc/$h$ with BIC$=30.39$, or $\bar{a}_{-1}=(16.5\pm 2.8)\times 10^{-3}$ (Mpc/$h$)2 with BIC$=36.79$.171717 For completeness, let us mention several additional models we investigated: the $(\bar{a}_{0},\bar{a}_{1},\bar{a}_{2})$ model has a BIC of $56.5$, and the $(\bar{a}_{-2},\bar{a}_{0},\bar{a}_{1},\bar{a}_{2})$/ $(\bar{a}_{-1},\bar{a}_{0},\bar{a}_{1},\bar{a}_{2})$/ $(\bar{a}_{-2},\bar{a}_{-1},\bar{a}_{0},\bar{a}_{1},\bar{a}_{2})$/ $(\bar{a}_{-2},\bar{a}_{-1},\bar{a}_{0},\bar{a}_{1})$ models have respectively a BIC score of $31.8,34.3,32.2,31.2$.

The inferred values for $\bar{a}_{-2}$ and $\bar{a}_{-1}$ can in principle be turned into an estimate of $f_{\rm NL}$, which we leave for future work.

This shows that, in presence of a non-gaussian distribution (of the local $f_{\rm NL}$ type) for the initial cosmological fields, the consistency relations in Eq. (5) are violated as expected Peloso and Pietroni (2013); Horn et al. (2014); Valageas et al. (2017).

III Discussion

The search for primordial non-gaussianities has been so far a challenging task. This is partly due to our lack of theoretical control over observables that are outside the linear regime. Consistency relations are non-perturbative statements that follow solely from symmetry arguments and, as such, might provide a key tool to overcome these difficulties.

In this paper we successfully test them, for the first time, in a regime well outside the domain of perturbation theory. In doing so, we highlight and solve a number of technical and conceptual subtleties associated with the analysis of the bispectrum in the squeezed regime from $N$-body simulations, whose systematic study has been lacking from the literature (see Nishimichi and Valageas (2014) for an exception, though see footnote 8).

Moreover, we show how in presence of non-gaussian initial conditions of the local $f_{\rm NL}$ type, significant deviations from the standard consistency relations are observed. This is the first step towards extracting constraints on $f_{\rm NL}$ from observational data, using the consistency relations (or violations thereof). The appeal is that with this method, (non-linear) modes that are normally discarded can now be used. Several issues need to be investigated before this goal can be realized. They include: checking the consistency relations (1) for biased observables such as halos in $N$-body simulations or galaxies in hydrodynamic simulations, and (2) including redshift space distortions. As explained in Sec. I, the consistency relations are expected to be robust against these complications, but it would be useful to test the expectations against simulations—our simple exercise presented in this paper suggests there could well be subtleties that need to be understood and addressed.

Acknowledgements.

The authors are grateful to A. Joyce, A. Lewis, A. Nicolis, R. Penco, M. Pietroni, M. Simonović and S. Wong for interesting and useful discussions, and to M. Abitbol for illuminating insights on the statistical analysis. The work done by A.E. is supported by the Swiss National Science Foundation under contract 200020-169696 and through the National Center of Competence in Research SwissMAP. The work done by L.H. is supported in part by the NASA grant NXX16AB27G and the DOE grant DE-SC011941.

Appendix A Checking soft momentum factorization

In this Appendix we show that the procedure outlined around Eqs. (9) and (10) might not eliminate the unwanted $q$ dependence if $q$ is extremely small.

As an explicit check let us consider the result obtained in perturbation theory. When the hard mode $k$ is within the linear regime one can easily show that the squeezed limit of the bispectrum (see e.g. Bernardeau et al. (2002); Horn et al. (2014)) gives

[TABLE]

Let us then consider the toy model $f(q,k,\theta)=a_{0}^{\text{pt}}(k,\theta)q$ (i.e. $f$ is our toy bispectrum for which $a_{0}^{\text{pt}}$ as a function of $k$ and $\theta$ is exactly known), and bin it as in Eqs. (9) and (11) over $k\in[79,91]k_{f}$ and all relative angles. Let us call the result $\bar{f}$. Recall that the worry was that, with a discrete set of triangles, the average over $k$ and $\theta$ of $a_{0}^{\text{pt}}$ (we call this $\bar{a}_{0}^{\text{pt}}$) would secretly depend on $q$. In our toy example, since the $k$ and $\theta$ dependence of $a_{0}^{\text{pt}}$ are precisely known, we can compute this average exactly without reference to the particular triangles we happen to have in our discrete grid. If our procedure to remove the unwanted $q$ dependence works, then it should be $\bar{f}(\bar{q})=\bar{a}_{0}^{\text{pt}}\bar{q}$.

In Figture 5 we report the result of our measurement. As one can see, the procedure works reasonably well only if the soft momentum is $q\gtrsim 3k_{f}$. In the analysis reported in Section II.3 we therefore exclude the lowest momentum bin. See Sec. II.2 for a discussion of why our procedure does not perfectly remove the unwanted $q$ dependence.

Finally, it can be checked that, in presence of local $f_{\text{NL}}$ non-gaussianities, the factorization in Eq. (11) holds better, even for low momenta. Indeed, the dominant $a_{-2}(k)$ term is expected to have no $\theta$-dependence and only a mild $k$-dependence Peloso and Pietroni (2013).