Quijote-PNG: The Information Content of the Halo Mass Function

As we discussed in the Introduction, the problem of finding an optimal summary statistic that minimizes the error bars on a given cosmological or PNG parameter is unsolved. An alternative to using summary statistics is to perform field-level analysis. The goal with this kind of analysis is to maximize the amount of information that can be extracted without relying on summary statistics. While there are many types of methods to perform such analysis, in our case, we made use of GNNs (Battaglia et al. 2018). The advantages of GNNs over other methods are that (1) they do not impose a cut on scales, (2) symmetries (e.g., rotational and translational invariance) can be easily implemented, and (3) they can be more interpretable than other methods. Because of this, we decided to train the GNNs to perform field-level likelihood-free inference.

As a starting point, we run 1000 simulations, each containing 512³ particles in a periodic box of size 1 h⁻¹ Gpc. Each of those simulations has a different initial random seed but also a different value of ${f}_{\mathrm{NL}}^{\mathrm{local}}$ in the range −300, +300. The value of the cosmological parameters was the same in all simulations. We then trained a GNN to perform field-level likelihood-free inference on the value of ${f}_{\mathrm{NL}}^{\mathrm{local}}$. The architecture and training procedure are the same as those outlined in de Santi et al. (2023), Shao et al. (2022), and Villanueva-Domingo & Villaescusa-Navarro (2022).

From this exercise, we found that our model was able to infer the value of ${f}_{\mathrm{NL}}^{\mathrm{local}}$ with an error of $\sigma ({f}_{\mathrm{NL}}^{\mathrm{local}})\sim 35$ at z = 0. In an attempt to understand the behavior of the network, we trained a deep set model (Zaheer et al. 2017), where the only information we made use of for the halos was their masses, not their spatial positions. By training such a model, we found that its performance was almost identical to that of the GNNs. We thus concluded that the network was likely not using the clustering of the halos to perform the inference. Therefore, the network should be using the abundance of halos to infer ${f}_{\mathrm{NL}}^{\mathrm{local}}$.

To verify this, we trained a simple model consisting of fully connected layers on the HMF of the halo catalogs from the simulations. We found that this model performed almost as well as the GNN. From this exercise, we reached the conclusion that the HMF is a summary statistic that contains lots of information, likely more than clustering-based statistics, as the GNN did not use those to perform the inference. We emphasize that we trained the GNN using halo catalogs from simulations that only vary ${f}_{\mathrm{NL}}^{\mathrm{local}}$. Therefore, our results did not account for degeneracies with cosmological parameters that could degrade the constraints, as we shall see below.

This motivated a further analysis, illustrated in the following sections, in which we explicitly extract the power spectrum, bispectrum, and HMF from the Quijote data set, as well as their covariance and response to variations in both cosmological and PNG parameters, in order to perform a full Fisher matrix forecast on nonlinear scales.

In this section, we recall the main ingredients of our Fisher analysis pipeline, which was previously used in Jung et al. (2022b).

The Fisher information matrix, defined as

$$\begin{eqnarray}&&{F}_{\mathrm{ij}}={\left(\displaystyle \frac{\partial \bar{{\boldsymbol{s}}}}{\partial {\theta }_{i}}\right)}^{{\rm{T}}}{{\boldsymbol{C}}}^{-1}\left(\displaystyle \frac{\partial \bar{{\boldsymbol{s}}}}{\partial {\theta }_{j}}\right),\end{eqnarray} \tag{ 1 }$$

allows us to estimate the variance, ${\sigma }^{2}({\theta }_{i})=\sqrt{{\left({F}^{-1}\right)}_{{ii}}}$, of the optimal unbiased estimator of a given summary statistic s with covariance C assuming the statistic is Gaussian-distributed ¹⁸ and neglecting the dependence of C itself on the parameters (Carron 2013).

In this work, both the covariance and derivatives are computed from the simulations described in Section 2. The covariance matrix is evaluated using

$$\begin{eqnarray}&&\hat{{\boldsymbol{C}}}=\displaystyle \frac{1}{{n}_{r}-1}({\boldsymbol{s}}-\bar{{\boldsymbol{s}}}){\left({\boldsymbol{s}}-\bar{{\boldsymbol{s}}}\right)}^{{\rm{T}}},\end{eqnarray} \tag{ 2 }$$

where n_r is the number of realizations at fiducial cosmology (15,000 here). Then, to obtain an unbiased estimate of the precision matrix, we apply the Hartlap correction factor (Hartlap et al. 2007),

$$\begin{eqnarray}&&{{\boldsymbol{C}}}^{-1}=\displaystyle \frac{{n}_{r}-{n}_{s}-2}{{n}_{r}-1}{\hat{{\boldsymbol{C}}}}^{-1},\end{eqnarray} \tag{ 3 }$$

where n_s is the length of the summary statistic vector s (note, however, that this correction is very small here, n_s ∼ 10², while n_r = 15,000).

The derivatives are calculated using finite difference,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \bar{{\boldsymbol{s}}}}{\partial {\theta }_{i}}=\displaystyle \frac{\bar{{\boldsymbol{s}}}({\theta }_{i}^{\mathrm{fid}}+\delta {\theta }_{i})-\bar{{\boldsymbol{s}}}({\theta }_{i}^{\mathrm{fid}}-\delta {\theta }_{i})}{2\delta {\theta }_{i}},\end{eqnarray} \tag{ 4 }$$

where we use the sets of 500 simulations where one parameter θ_i is displaced by ±δ θ_i with respect to its fiducial value. However, it was noticed in Coulton et al. (2023b) and Jung et al. (2022b) that this number of realizations was not sufficient to obtain fully converged derivatives of the halo power spectrum and bispectrum, leading to spuriously low predictions when analyzing jointly cosmological parameters and PNG amplitudes. To overcome this issue, a conservative approach to Fisher matrix computations was developed in Coulton & Wandelt (2023) and Coulton et al. (2023b) that is based on computing the Fisher matrix from maximally compressed statistics instead of working with the summary statistics directly.

As shown in Heavens et al. (2000) and Alsing & Wandelt (2018), the compressed quantity defined by

$$\begin{eqnarray}&&{\tilde{s}}_{i}=\displaystyle \frac{\partial \bar{{\boldsymbol{s}}}}{\partial {\theta }_{i}}{{\boldsymbol{C}}}^{-1}({\boldsymbol{s}}-\bar{{\boldsymbol{s}}})\end{eqnarray} \tag{ 5 }$$

conserves all of the statistical information about the parameter θ_i contained in the data vector s , if s follows a Gaussian likelihood (hence, the same assumption as for the Fisher matrix in Equation (1)). This compression uses the same ingredients as for the Fisher matrix computation (covariance and derivatives of s ), with the addition of the mean $\bar{{\boldsymbol{s}}}$ that is trivial to evaluate from the simulations at fiducial cosmology. Repeating the process for all parameters of interest in θ, one can then compute the Fisher matrix of the compressed statistics $\tilde{{\boldsymbol{s}}}$ by substituting it for s in Equation (1). In practice, one has to separate the initial data set into two subsets. The first is used to perform the compression (i.e., compute the derivatives in Equation (5)), and the second is compressed (i.e., s in Equation (5)) and then used to calculate the derivatives $\partial \tilde{{\boldsymbol{s}}}/\partial {\theta }_{i}$ and covariance $\hat{{\boldsymbol{C}}}$ of the compressed statistics to obtain a conservative estimation of the Fisher matrix. In this work, we use 80% and 20% of the simulations for the two steps, respectively, which have been verified to give optimal and numerically stable results. We repeat the procedure for many random splits of the data (between the two steps) and average the results to minimize the intrinsic variance of the method.

Finally, as shown in Coulton & Wandelt (2023), computing the following combination of the standard (overoptimistic) and compressed (conservative) Fisher matrices,

$$\begin{eqnarray}&&{F}^{\mathrm{combined}}=G({F}^{\mathrm{standard}},{F}^{\mathrm{compressed}}),\end{eqnarray} \tag{ 6 }$$

where G corresponds to the geometric mean defined by

$$\begin{eqnarray}&&G(A,B)={A}^{\tfrac{1}{2}}{\left({A}^{-\tfrac{1}{2}}{{BA}}^{-\tfrac{1}{2}}\right)}^{\tfrac{1}{2}}{A}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 7 }$$

gives unbiased estimates of the Fisher error bars with a much smaller number of simulations. An illustration of the different convergences for the three methods is provided in Appendix C.

In addition to the halo power spectrum and bispectrum, we consider the HMF, defined as the number of dark matter halos per unit of comoving volume per unit of logarithmic mass bins.

We measure it in the Quijote simulations using 15 logarithmic bins corresponding to halo masses M between approximately 2.0 × 10¹³ and 4.6 × 10¹⁵ M_⊙ h^–1 (note, however, that we do not use the first two bins in the analyses presented in Section 4). To be exact, we use the same binning as in Bayer et al. (2021), where the counted halos each contain between 30 and 7000 dark matter particles. ¹⁹

In Figure 1, we show the impact of the three shapes of PNG on the HMF. Both the local and equilateral shapes increase the number of massive halos for a positive f_NL value (and decrease it for a negative f_NL) and have very degenerate signatures, while for orthogonal PNG, it is the opposite. For less massive halos, the effect of PNG changes sign (with the switch occurring for higher masses for orthogonal PNG, which is the only one that appears in the mass range of the plot at z = 1). This effect was already present in early works on the HMF with PNG simulations (see, e.g., LoVerde et al. 2008) and is due to the fact that, at fixed Ω_m , more massive halos can only appear at the expense of less massive halos and matter in smaller structures.

Zoom In Zoom Out Reset image size

As a preliminary exercise, in Figure 2, we show the constraining power of the HMF on the PNG amplitudes f_NL of the three shapes, assuming the exactly known cosmological parameters. As expected, the HMF is, in this case, extremely sensitive to the presence of PNG, leading to even tighter constraints than the power spectrum and bispectrum. For example, our Fisher forecast on PNG of the local type is $\sigma ({f}_{\mathrm{NL}}^{\mathrm{local}})\sim 30$ at z = 0, which is in very good agreement with the GNN and deep set results $\sigma ({f}_{\mathrm{NL}}^{\mathrm{local}})\simeq 35$ (see Section 3.1) and more than twice as small as the equivalent power spectrum + bispectrum forecast error bar.

Zoom In Zoom Out Reset image size

However, it is well known that there are large degeneracies between f_NL and several cosmological parameters, like σ₈ or Ω_m (Maturi et al. 2011), as can be verified in Figure 1.

When we jointly analyze all parameters, these degeneracies increase the errors significantly (by roughly 1 order of magnitude at z = 1 and slightly less at z = 0, where the change of sign of the f_NL derivative, seen in Figure 1, helps distinguish it from the response to variations in other cosmological parameters), making them larger than those achievable from the power spectrum and bispectrum combination.

While, as expected, the HMF alone does not produce competitive f_NL constraints in comparison with the power spectrum and bispectrum, it does remain interesting to investigate whether a combined analysis of all three statistics can produce significant improvements; this is the main point of the present work. Complementing our previous power spectrum + bispectrum analysis with the HMF can, in principle, benefit us in two ways. First of all, it directly adds extra information about the f_NL parameter; also, it could be useful to help break the important degeneracy between f_NL and the so-called b_ϕ bias parameter.

Before presenting our results, let us review and discuss the latter point in more detail. In the presence of local PNG, the halo density fluctuation field δ_h (z) can be written to leading order as follows (Dalal et al. 2008; Matarrese & Verde 2008; Slosar et al. 2008; McDonald 2008; Giannantonio & Porciani 2010; Desjacques & Seljak 2010a):

$$\begin{eqnarray}&&{\delta }_{h}(z)=\left[{b}_{1}(z)+\displaystyle \frac{3{{\rm{\Omega }}}_{m}{H}_{0}^{2}}{2D(z){k}^{2}}{b}_{\phi }{f}_{\mathrm{NL}}\right]{\delta }_{m}(z),\end{eqnarray} \tag{ 8 }$$

where δ_m is the matter density fluctuation, D(z) is the growth factor, and b₁ and b_ϕ are the bias parameters, defined respectively as the response of δ_h to mass density δ_m and primordial potential ϕ. It is evident in this relation that the scale-dependent signature depends on both b_ϕ and f_NL and that the two parameters are completely degenerate. This issue can be avoided if one assumes, as was generally done, the universality relation between b₁ and b_ϕ , that is,

$$\begin{eqnarray}&&{b}_{\phi }=2{\delta }_{c}({b}_{1}-1),\end{eqnarray} \tag{ 9 }$$

where δ_c is the critical density for collapse. However, it was recently pointed out in Barreira (2020, 2022) that such a relation does not accurately describe the bias of either galaxies, selected by stellar mass, or halos, selected by concentration. Therefore, b_ϕ is not exactly determined anymore, and this reintroduces the b_ϕ –f_NL degeneracy problem. To overcome the issue, different studies have focused on using simulations to produce accurate priors on b_ϕ (Lazeyras et al. 2023) and exploiting the multitracer technique (Barreira & Krause 2023; Karagiannis et al. 2023; Sullivan et al. 2023). In the present context, the idea is instead to try and break the degeneracy by exploiting the information in the HMF—which selects all halos in each given mass bin—and its direct dependence on f_NL and not on b_ϕ .

For clarity, we split the discussion of our results into two parts. Initially, we assume universality in the b_ϕ (b₁) relation using Equation (9), and we measure the sheer extra information content in the HMF, in the absence of the b_ϕ –f_NL degeneracy. ²⁰ Later on, we instead treat b_ϕ as a free parameter.

The outcome of the first part of the analysis (assuming universality in b_ϕ (b₁)) is illustrated in Figures 3 and 4 (see also Table 2). We see that by adding the HMF, the error bars on σ₈ and ${f}_{\mathrm{NL}}^{\mathrm{equil}}$ become roughly twice as small as the power spectrum + bispectrum result. Moreover, there is also a noticeable improvement for Ω_m and ${f}_{\mathrm{NL}}^{\mathrm{ortho}}$. For ${f}_{\mathrm{NL}}^{\mathrm{local}}$, there is instead no clear improvement; this seems to be due to the fact that in this case, the information content is totally dominated by the power spectrum contribution via scale-dependent bias. Such a contribution is instead smaller for the orthogonal shape and absent for the equilateral case, making the HMF inclusion more important for these scenarios, especially the equilateral one.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. The 1σ Constraints on the Cosmological Parameters and PNG Amplitudes at z = 1 Obtained by Combining the Information of the Halo Power Spectrum, Bispectrum, and Mass Function, Each Measured from the Quijote and Quijote-png Simulations

	bϕ Fixed	No Prior on bϕ	Planck Priors
σ8	0.012	0.013	0.005
Ωm	0.018	0.017	0.002
ns	0.075	0.075	0.003
h	0.072	0.071	0.017
${f}_{\mathrm{NL}}^{\mathrm{local}}$	40	89	34
${f}_{\mathrm{NL}}^{\mathrm{equil}}$	203		136
${f}_{\mathrm{NL}}^{\mathrm{ortho}}$	85		79
${M}_{\min }/{10}^{13}$	0.019	0.045	0.009

Download table as:  ASCIITypeset image

Note that we consider only halos with masses above ∼4 × 10¹³ M_⊙ h^–1 in the HMF, which is larger than the fiducial ${M}_{\min }=3.2\,\times \,{10}^{13}{M}_{\odot }/h$ used to study the power spectrum and bispectrum. This means that the HMF is not sensitive at all to small variations of ${M}_{\min }$ around the fiducial value. However, through cross-correlated terms with the other summary statistics, the error bars on ${M}_{\min }$ are almost 2 orders of magnitude smaller. ²¹ In Appendix C, we verify the numerical stability of our results by varying the number of simulations used.

It is interesting to check which halo mass range gives the largest contribution to the observed improvements. To this purpose, we repeat the analysis by varying ${M}_{\min }^{\mathrm{HMF}}$, the lowest mass bin of the HMF used to evaluate Fisher matrices. Our results are displayed in Figure 5, which highlights different behaviors for the different parameters considered. Most importantly, for the two PNG parameters ${f}_{\mathrm{NL}}^{\mathrm{equil}}$ and ${f}_{\mathrm{NL}}^{\mathrm{ortho}}$, halos of intermediate masses (∼2–6 × 10¹⁴ M_⊙ h^–1 at z = 0 and slightly smaller at z = 1) play a significant role in the observed improvement of constraints, while less massive halos, despite being more numerous, have a much smaller effect. However, the situation is different for cosmological parameters like σ₈ and Ω_m , where those same less massive halos contain most of the information.

Zoom In Zoom Out Reset image size

Accounting for the effects of b_ϕ in our methodology is not straightforward, since b_ϕ cannot be explicitly included as an input parameter in our simulations, and this does not allow us to directly compute the numerical derivative ∂ s /∂b_ϕ . To circumvent this issue in a simple way and be able to perform a first test of the ability of the HMF to remove degeneracies between b_ϕ and ${f}_{\mathrm{NL}}^{\mathrm{local}}$, we decide here to work under the conservative assumption that these two parameters are fully degenerate at the level of the halo power spectrum and bispectrum. In other words, we assume that $\partial {\boldsymbol{s}}/\partial {b}_{\phi }\propto \partial {\boldsymbol{s}}/\partial {f}_{\mathrm{NL}}^{\mathrm{local}}$, where s is either the power spectrum or the bispectrum. For the HMF, we instead set the derivative with respect to b_ϕ equal to zero, as it does not depend on this parameter, and compute the ${f}_{\mathrm{NL}}^{\mathrm{local}}$ derivative as usual.

In Figure 6, we show the 1σ Fisher constraints obtained in this assumption and compare them with the "ideal" (b_ϕ fixed) constraints derived in the previous section for different ${k}_{\max }$ (see also Table 2).

Zoom In Zoom Out Reset image size

The most important result here is that the inclusion of the HMF makes it possible to break the b_ϕ –${f}_{\mathrm{NL}}^{\mathrm{local}}$ degeneracy to a level that allows us to produce meaningful ${f}_{\mathrm{NL}}^{\mathrm{local}}$ constraints without resorting to any prior information on b_ϕ . The final ${f}_{\mathrm{NL}}^{\mathrm{local}}$ forecast is, however, degraded by a factor of ∼2.5 with respect to the idealized, b_ϕ fixed case that was shown in Figure 11 in Appendix C. In order to achieve this constraining level, it is also crucial to include the information from the power spectrum and bispectrum at nonlinear scales (k between 0.2 and 0.5 h Mpc⁻¹), as it helps break degeneracies with several cosmological parameters (Ω_m in particular).

We corroborate our findings with a simulation-independent analysis based on the halo model (for a review, see Cooray & Sheth 2002; Asgari et al. 2023). Within this framework, we describe the HMF and halo power spectrum following Takada & Spergel (2014) up to ${k}_{\max }=0.2\,h\,{\mathrm{Mpc}}^{-1}$. We use the HMF and bias from Tinker et al. (2010) using M_200,m directly as the mass definition in the mass integration. In the power spectrum analysis of the simulations, the halos are considered pointlike; thus, we use a Dirac delta as the halo profile. Thanks to the low ${k}_{\max }$ we use, the two-halo term dominates the signal, and this approximation is appropriate. The effect of PNG—here we only consider the local model—is included as a correction to the HMF parameterized according to LoVerde & Smith (2011) and through the scale-dependent halo bias shown in Equation (8). While we are aware that the M_200,m mass does not match the FOF mass used in the rest of the paper, we still consider the HMF divided in 10 bins logarithmically spaced between 3.2 × 10¹³ and 3.2 × 10¹⁵ M_⊙ h^–1 as observable. We bin the halo power spectrum in 30 bins logarithmically spaced between 6.3 × 10⁻³ and 0.2 h Mpc⁻¹. We choose a relatively low ${k}_{\max }$ to ensure that nonlinearities are negligible at this stage. In the HMF–halo power spectrum covariance, for which we again follow Takada & Spergel (2014), only the Gaussian terms are included at present. A more refined analysis, including a wider range of scales and masses, the complete covariance, the uncertainties on the parameterization of the HMF, and, crucially, the bispectrum, will be presented in a future work (Ravenni et al. 2023, in preparation).

The results are shown in Figure 7, which highlights a very good agreement between our preliminary theoretical computations and the purely simulation-based forecast. This result confirms that a joint analysis including the HMF is an interesting approach that deserves further investigation and could be adopted as a complementary strategy to those already implemented in the literature to address the b_ϕ –${f}_{\mathrm{NL}}^{\mathrm{local}}$ degeneracy issue.

Zoom In Zoom Out Reset image size

As highlighted in Section 4.2, removing the degeneracies of the HMF using the information from the halo power spectrum and bispectrum significantly improves the constraints on PNG of the equilateral type. In this section, we push the idea further by assuming strong but realistic priors on cosmological parameters based on CMB measurements from Planck.

We use the same Gaussian likelihood based on the Planck CMB data (Aghanim et al. 2020) as in Uhlemann et al. (2020) in Figure 8 in addition to our HMF, power spectrum, and bispectrum measurements to derive 1σ Fisher constraints (see also Table 2). For both ${f}_{\mathrm{NL}}^{\mathrm{local}}$ and ${f}_{\mathrm{NL}}^{\mathrm{equil}}$, it improves these constraints, while the effect is smaller for ${f}_{\mathrm{NL}}^{\mathrm{ortho}}$. Note also that the effect is strongest when the HMF is also considered in the analysis, meaning it removes degeneracies between the PNG and cosmological parameters at the level of the HMF. Concerning numerical convergence with the number of simulations used to compute the derivatives, including these Planck priors also improves it significantly, where only ${f}_{\mathrm{NL}}^{\mathrm{equil}}$ is not optimally constrained for the power spectrum + bispectrum case, and all parameters have converged when we add the HMF information.

Zoom In Zoom Out Reset image size

The trispectrum is defined as

$$\begin{eqnarray}&&\langle \delta ({{\boldsymbol{k}}}_{1})\delta ({{\boldsymbol{k}}}_{2})\delta ({{\boldsymbol{k}}}_{3})\delta ({{\boldsymbol{k}}}_{4})\rangle =T({k}_{1},{k}_{2},{k}_{3},{k}_{4},{K}_{a},{K}_{b}),\end{eqnarray} \tag{ A1 }$$

where k_i = ∣ k _i ∣, K_a = ∣ k ₁ + k ₂∣, and K_b = ∣ k ₁ + k ₃∣. Estimating the full trispectrum is computationally highly challenging; so, in this work, we measure trispectra averaged over K_b , i.e.,

$$\begin{eqnarray}&&{ \mathcal T }({k}_{1},{k}_{2},{k}_{3},{k}_{4},{K}_{a})\propto \displaystyle \sum _{{K}_{b}}T({k}_{1},{k}_{2},{k}_{3},{k}_{4},{K}_{a},{K}_{b}).\end{eqnarray} \tag{ A2 }$$

A binned version of this can be estimated as

$$\begin{eqnarray}&&\begin{array}{l}\hat{{ \mathcal T }}({k}_{a},{k}_{b},{k}_{c},{k}_{d},{K}_{E})=\displaystyle \frac{1}{{N}_{a,b,c,d,E}}\displaystyle \int \displaystyle \prod _{i\,=\,1,4}\displaystyle \frac{{{\rm{d}}}^{3}{k}_{i}}{{\left(2\pi \right)}^{3}}\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{3}{K}_{a}}{{\left(2\pi \right)}^{3}}\\ \times {\left(2\pi \right)}^{3}{\delta }^{(3)}({{\boldsymbol{k}}}_{1}+{{\boldsymbol{k}}}_{2}-{{\boldsymbol{K}}}_{a}){\left(2\pi \right)}^{3}{\delta }^{(3)}({{\boldsymbol{K}}}_{a}-{{\boldsymbol{k}}}_{3}-{{\boldsymbol{k}}}_{4}){W}_{a}({{\boldsymbol{k}}}_{1})\\ \times {W}_{b}({{\boldsymbol{k}}}_{2}){W}_{c}({{\boldsymbol{k}}}_{3}){W}_{d}({{\boldsymbol{k}}}_{4}){W}_{E}({{\boldsymbol{K}}}_{a})\delta ({{\boldsymbol{k}}}_{1})\delta ({{\boldsymbol{k}}}_{2})\delta ({{\boldsymbol{k}}}_{3})\delta ({{\boldsymbol{k}}}_{4}),\end{array}\end{eqnarray} \tag{ A3 }$$

where W_a (k) selects modes that lie within binned a, and N_a,b,c,d,E is the normalization. In this work, we use 14 equally spaced bins between k = 0.0102 and 0.193 h Mpc^–1. By utilizing

$$\begin{eqnarray}&&{\delta }^{(3)}({{\boldsymbol{k}}}_{1}+{{\boldsymbol{k}}}_{2}+{{\boldsymbol{k}}}_{3})=\int {{\rm{d}}}^{3}{{xe}}^{i{\boldsymbol{x}}\cdot {\boldsymbol{k}})}\delta ({\boldsymbol{k}}),\end{eqnarray} \tag{ A4 }$$

we efficiently implement the estimator by first computing

$$\begin{eqnarray}&&{\delta }_{{W}_{a}}({\boldsymbol{x}})=\int \displaystyle \frac{{{\rm{d}}}^{3}{k}_{a}}{{\left(2\pi \right)}^{3}}{e}^{i{\boldsymbol{x}}\cdot k}\delta ({\boldsymbol{k}}){W}_{a}({\boldsymbol{k}}),\end{eqnarray} \tag{ A5 }$$

then computing

$$\begin{eqnarray}&&\begin{array}{l}{D}_{{ab}}({\boldsymbol{K}})=\displaystyle \int {{\rm{d}}}^{3}{{xe}}^{-i{\boldsymbol{x}}\cdot {\boldsymbol{K}}}{\delta }_{{W}_{a}}({\boldsymbol{x}}){\delta }_{{W}_{b}}({\boldsymbol{x}}),\end{array}\end{eqnarray} \tag{ A6 }$$

and then the estimate is given by

$$\begin{eqnarray}&&\begin{array}{l}\hat{{ \mathcal T }}({k}_{a},{k}_{b},{k}_{c},{k}_{d},{K}_{E})\\ =\ \displaystyle \frac{1}{{N}_{a,b,c,d,E}}\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{3}K}{{\left(2\pi \right)}^{3}}{D}_{{ab}}({\boldsymbol{K}}){D}_{{cd}}(-{\boldsymbol{K}}){W}_{E}({\boldsymbol{K}}).\end{array}\end{eqnarray} \tag{ A7 }$$

The normalization is obtained by evaluating this estimator (without the N_a,b,c,d,E term) on maps with δ( k ) = 1.

To perform a stringent test of the ICs, we study the difference between the trispectrum of the ICs with f_NL ≠ 0 and f_NL = 0, i.e.,

$$\begin{eqnarray}&&\begin{array}{l}{\hat{{ \mathcal T }}}^{\mathrm{diff}}({k}_{a},{k}_{b},{k}_{c},{k}_{d},{K}_{E})={\hat{{ \mathcal T }}}^{{f}_{\mathrm{NL}}\ne 0}({k}_{a},{k}_{b},{k}_{c},{k}_{d},{K}_{E})\\ \,\,-\,{\hat{{ \mathcal T }}}^{{f}_{\mathrm{NL}}=0}({k}_{a},{k}_{b},{k}_{c},{k}_{d},{K}_{E}).\end{array}\end{eqnarray} \tag{ A8 }$$

This cancels the leading noise contribution to the trispectrum measurement.

The results are shown in Figure 9. For equilateral, there is no detectable trispectrum. For orthogonal NG, there are small hints of a trispectrum signal. As this measurement uses 200 simulations and a method to cancel the cosmic variance, it is likely that this small trispectrum is negligible. However, the local case shows significant evidence of a trispectrum. This is not unexpected. Local PNG is generated in these simulations by

$$\begin{eqnarray}&&{\rm{\Phi }}({\boldsymbol{x}})={\phi }^{G}({\boldsymbol{x}})+{f}_{\mathrm{NL}}\left({\phi }^{G}{\left({\boldsymbol{x}}\right)}^{2}-\langle ({\phi }^{G}({\boldsymbol{x}})\rangle \right),\end{eqnarray} \tag{ A9 }$$

where Φ^G ( x ) is the Gaussian primordial potential. This generates a primordial trispectrum known as τ_NL (Kogo & Komatsu 2006). In many inflationary models, τ_NL is generated with local NG; thus, the trispectrum seen here is physical.

Zoom In Zoom Out Reset image size

These trispectra measurements suggest that unphysical higher-order N-point functions are not significant in our simulations.