Model-independent analyses of non-Gaussianity in Planck CMB maps using Minkowski Functionals

TL;DR

This paper develops a model-independent framework using Minkowski Functionals and Hermite expansions to analyze non-Gaussianity in Planck CMB maps, confirming weak non-Gaussian signals consistent with Gaussianity.

Contribution

It introduces a novel, model-independent method employing Minkowski Functionals and perturbative expansions to quantify non-Gaussianity in CMB data.

Findings

Confirmed weak non-Gaussianity in Planck 2015 maps

Validated the effectiveness of Hermite and perturbative expansions

Provided explicit formulas for cumulant-based coefficients

Abstract

Despite the wealth of $Planck$ results, there are difficulties in disentangling the primordial non-Gaussianity of the Cosmic Microwave Background (CMB) from the secondary and the foreground non-Gaussianity (NG). For each of these forms of NG the lack of complete data introduces model-dependencies. Aiming at detecting the NGs of the CMB temperature anisotropy $\delta T$, while paying particular attention to a model-independent quantification of NGs, our analysis is based upon statistical and morphological univariate descriptors, respectively: the probability density function $P(\delta T)$, related to ${\mathrm v}_{0}$, the first Minkowski Functional (MF), and the two other MFs, ${\mathrm v}_{1}$ and ${\mathrm v}_{2}$. From their analytical Gaussian predictions we build the discrepancy functions $\Delta_{k}$ ($k=P,0,1,2$) which are applied to an ensemble of $10^{5}$ CMB realization maps of the $\Lambda$CDM model and to the $Planck$ CMB maps. In our analysis we use general Hermite expansions of the $\Delta_{k}$ up to the $12^{th}$ order, where the coefficients are explicitly given in terms of cumulants. Assuming hierarchical ordering of the cumulants, we obtain the perturbative expansions generalizing the $2^{nd}$ order expansions of Matsubara to arbitrary order in the standard deviation $\sigma_0$ for $P(\delta T)$ and ${\mathrm v}_0$, where the perturbative expansion coefficients are explicitly given in terms of complete Bell polynomials. The comparison of the Hermite expansions and the perturbative expansions is performed for the $\Lambda$CDM map sample and the $Planck$ data. We confirm the weak level of non-Gaussianity ($1$-$2$)$\sigma$ of the foreground corrected masked $Planck$ $2015$ maps.

Tables

Table 1. Table 1: Table of the moments α n subscript 𝛼 𝑛 \alpha_{n} (main term, equation ( 111 )), cumulants ϰ n subscript italic-ϰ 𝑛 \varkappa_{n} , and dimensionless normalized cumulants C n subscript 𝐶 𝑛 C_{n} of the probability density function P ​ ( τ ) 𝑃 𝜏 P(\tau) .

(Units of $\mu K^n$)     $\Lambda$CDM sample  Full individual map range, no mask, 2∘fwhm, bin 13$\mu$K   c.f. A
$n$	$3$	$4$	$5$	$6$	$7$	$8$	$9$	$10$	$11$	$12$
${\alpha }_n$	$-110.103727$	$37985558.5$	$-4040069.36$	$6.8\times {10}^{11}$	$-1.5\times {10}^{11}$	$1.7\times {10}^{16}$	$-6.6\times {10}^{15}$	$5.4\times {10}^{20}$	$-3.1\times {10}^{20}$	$2.1\times {10}^{25}$
${\varkappa }_n$	$-110.103727$	$737.480513$	$-122231.715$	$-108751936.$	$785344215.$	$8.2\times {10}^{11}$	$1.1\times {10}^{14}$	$-1.1\times {10}^{16}$	$2.4\times {10}^{18}$	$8.4\times {10}^{20}$
$C_n$	$-5.187\times {10}^{-4}$	$5.82\times {10}^{-5}$	$-1.618\times {10}^{-4}$	$-2.4138\times {10}^{-3}$	$2.922\times {10}^{-4}$	$5.1255\times {10}^{-3}$	$1.14291\times {10}^{-2}$	$-1.96036\times {10}^{-2}$	$7.14927\times {10}^{-2}$	$0.4147370$
(Units of $\mu K^n$)     $\Lambda$CDM sample   Equal temperature range   (ETR $\pm$ 201$\mu K$), U73 mask, 2∘fwhm, bin 6$\mu$K
$n$	$0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$
${\alpha }_n$	$1.0001469$	$-1.4912\times {10}^{-3}$	$3515.73711$	$-113.475560$	$36119872.5$	$-2963938.41$	$5.9\times {10}^{11}$	$-7.8\times {10}^{10}$	$1.2\times {10}^{16}$	$-2.2\times {10}^{15}$
${\varkappa }_n$	$1.0$	$-1.4912\times {10}^{-3}$	$3515.73710$	$-97.7480482$	$-961350.317$	$741927.368$	$-1.3\times {10}^{10}$	$-3.4\times {10}^9$	$1.8\times {10}^{14}$	$-1.6\times {10}^{14}$
$C_n$	$1.0$	$-2.51\times {10}^{-5}$	$1.0$	$-4.689\times {10}^{-4}$	$-7.77766\times {10}^{-2}$	$1.0123\times {10}^{-3}$	$-0.2973678$	$-1.3108\times {10}^{-3}$	$1.1699687$	$-1.76138\times {10}^{-2}$

Table 2. Table 2: Table of coefficients a P ​ ( n ) subscript 𝑎 𝑃 𝑛 a_{P}(n) for the Λ Λ \Lambda CDM model, computed from table 1 using equation ( 2.2 ) in 1 s ​ t superscript 1 𝑠 𝑡 1^{st} line, then using equation ( 28 ) in 2 n ​ d superscript 2 𝑛 𝑑 2^{nd} line.

$\Lambda$CDM sample  Full individual map range, no mask, 2∘fwhm, bin 13$\mu$K   c.f. eqs. (2.2), (28) and model in A
$n$	$3$	$4$	$5$	$6$	$7$	$8$	$9$	$10$	$11$	$12$
$a_P\left(n\right)$	$-5.1872\times {10}^{-4}$	$5.825\times {10}^{-5}$	$-1.6184\times {10}^{-4}$	$-2.41112\times {10}^{-3}$	$2.9116\times {10}^{-4}$	$5.13035\times {10}^{-3}$	$1.153302\times {10}^{-2}$	$-1.964794\times {10}^{-2}$	$7.123993\times {10}^{-2}$	$0.41622908$
$a_P\left(n\right)$	$-5.1872\times {10}^{-4}$	$5.830\times {10}^{-5}$	$-1.6184\times {10}^{-4}$	$-2.40919\times {10}^{-3}$	$2.9116\times {10}^{-4}$	$5.19401\times {10}^{-3}$	$1.153302\times {10}^{-2}$	$-1.785385\times {10}^{-2}$	$7.123994\times {10}^{-2}$	$0.45721793$
$\Lambda$CDM sample   Equal temperature range   (ETR $\pm$ 201$\mu K$), U73 mask, 2∘fwhm, bin 6$\mu$K
$n$	$0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$
$a_P\left(n\right)$				$-5.4435\times {10}^{-4}$	$-7.777654\times {10}^{-2}$	$1.01233\times {10}^{-3}$	$-2.9736480\times {10}^{-1}$	$1.7101\times {10}^{-4}$	$1.38165948$
$a_P\left(n\right)$	$9.9534\times {10}^{-4}$	$-2.8\times {10}^{-7}$	$1.002019\times {10}^{-2}$	$-4.7175\times {10}^{-4}$	$5.20385\times {10}^{-3}$	$1.00469\times {10}^{-3}$	$-1.989114\times {10}^{-2}$	$9.19\times {10}^{-6}$	$-2.1129824\times {10}^{-1}$

Table 3. Table 3: Table of coefficients a 0 ​ ( n ) subscript 𝑎 0 𝑛 a_{0}(n) , computed from table 1 using equation ( 2.3 ).

$\Lambda$CDM   Full individual map range, no mask, 2∘fwhm, bin 13$\mu$K   c.f. equation (2.3) and model in A
$n$	$0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$
$a_0\left(n\right)$			$-1.7295\times {10}^{-4}$	$1.456\times {10}^{-5}$	$-3.237\times {10}^{-5}$	$-4.0185\times {10}^{-4}$	$4.159\times {10}^{-5}$	$6.4130\times {10}^{-4}$	$1.28145\times {10}^{-3}$
$\Lambda$CDM   Equal temperature range   (ETR $\pm$ 201$\mu K$), U73 mask, 2∘fwhm, bin 6$\mu$K
$n$	$0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$
$a_0\left(n\right)$			$-1.8145\times {10}^{-4}$	$-1.944413\times {10}^{-2}$	$2.0247\times {10}^{-4}$	$-4.95608\times {10}^{-3}$	$2.443\times {10}^{-5}$	$1.7270744\times {10}^{-1}$	$-1.54858\times {10}^{-3}$

Table 4. Table 4: Table of coefficients a 1 ​ ( n ) subscript 𝑎 1 𝑛 a_{1}(n) and a 2 ​ ( n ) subscript 𝑎 2 𝑛 a_{2}(n) .

$\Lambda$CDM sample  Full individual map range, no mask, 2∘fwhm, bin 13$\mu$K   c.f. A
$n$	$0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$
$a_1\left(n\right)$	$-2.21\times {10}^{-6}$	$-2.808\times {10}^{-5}$	$3.97274\times {10}^{-3}$	$-5.0873\times {10}^{-4}$	$-5.823\times {10}^{-5}$	$-6.643\times {10}^{-5}$	$-1.97922\times {10}^{-3}$	$5.059\times {10}^{-5}$	$6.37408\times {10}^{-3}$
$a_2\left(n\right)$	$1.927\times {10}^{-5}$	$-2.148\times {10}^{-5}$	$-1.7313\times {10}^{-4}$	$1.183631\times {10}^{-2}$	$-1.53488\times {10}^{-3}$	$-4.8941\times {10}^{-4}$	$-1.20822\times {10}^{-3}$	$-9.58844\times {10}^{-3}$	$-7.33397\times {10}^{-3}$

Table 5. Table 5: Table of μ 𝜇 \mu and σ 𝜎 \sigma values.

(Units of $\mu K$)     $\Lambda$CDM sample  Full individual map range, no mask, 2∘fwhm, bin 13$\mu$K
${\mu }_{px}$	${\sigma }_{0px}$	${\sigma }_{1px}$	${\sigma }_{1px}/{\sigma }_{0px}$
$-3.4\times {10}^{-7}$	$59.53348$	$2020.08398$	$33.93190$
(Units of $\mu K$)     $\Lambda$CDM sample   Equal temperature range   (ETR $\pm$ 201$\mu K$), U73 mask, 2∘fwhm, bin 6$\mu$K
${\mu }_{px}$	${\sigma }_{0px}$	${\sigma }_{1px}$	${\sigma }_{1px}/{\sigma }_{0px}$
$-1.4534\times {10}^{-3}$	$59.13275$	$2020.10372$	$34.16218$

Equations

$$F(\tau ):=\mathrm{prob}(\delta T<\tau )={\int }_{-\infty }^{\tau }P({\tau }^{\prime })\mathrm{d}{\tau }^{\prime },$$

$$F_C(\tau ):=\mathrm{prob}(\delta T\ge \tau )={\int }_{\tau }^{\infty }P({\tau }^{\prime })\mathrm{d}{\tau }^{\prime }=1-F(\tau ),$$

$$⟨f(\delta T)⟩:={\int }_{-\infty }^{\infty }f(\tau )P(\tau )\mathrm{d}\tau ,$$

$${\alpha }_n:=⟨(\delta T{)}^n⟩:={\int }_{-\infty }^{\infty }{\tau }^{n}P(\tau )\mathrm{d}\tau ;n=0,1,2,\cdots ,$$

$$m_n:=⟨(\delta T-\mu {)}^n⟩:={\int }_{-\infty }^{\infty }(\tau -\mu {)}^{n}P(\tau )\mathrm{d}\tau ,$$

$$M(x):=⟨{\text{rme}}^{x\delta T}⟩={\int }_{-\infty }^{\infty }{\text{rme}}^{x\tau }P(\tau )\mathrm{d}\tau =\sum _{n=0}^{\infty }{\alpha }_n\frac{x^n}{n!},$$

$$C(x):=\ln M(x)=\sum _{n=1}^{\infty }{\varkappa }_n\frac{x^n}{n!}.$$

$${\sigma }_0^2:=⟨(\delta T-\mu {)}^2⟩=m_2={\alpha }_2-{\mu }^2={\varkappa }_2,$$

$$m_{n}=\sum_{r=0}^{n}\left(\begin{array}[]{c}n\\ r\end{array}\right)\,(-1)^{n-r}\mu^{n-r}\alpha_{r}\quad;\quad n=0,1,2,\cdots\;\;,$$

$$\varkappa_{n}=m_{n}-\sum_{r=1}^{n-1}\left(\begin{array}[]{c}n-1\\ r-1\end{array}\right)\,\varkappa_{r}m_{n-r}\;\;.$$

$$m_3={\varkappa }_3={\alpha }_3-3\mu {\alpha }_2+2{\mu }^3,$$

$${\gamma }_1:=\frac{m_3}{{\sigma }_0^3}.$$

$${\gamma }_2:=\frac{m_4}{{\sigma }_0^4}-3.$$

$$⟨f(\delta T)⟩:={\int }_{-\infty }^{\infty }f(\tau )\mathrm{d}F(\tau ).$$

$$⟨f(\delta T)⟩:={\int }_{-\infty }^{\infty }f(\tau )F^{\prime }(\tau )\mathrm{d}\tau +\sum _{\ell }f({\tau }_{\ell })[F({\tau }_{\ell }+0)-F({\tau }_{\ell }-0)],$$

$$P^{\mathrm{G}}(\tau )=\frac{1}{\sqrt{2\pi }{\sigma }_0}{\text{rme}}^{-\frac{(\tau -\mu {)}^2}{2{\sigma }_0^2}},$$

$$M^{\mathrm{G}}(x)={\text{rme}}^{\mu x+\frac{{\sigma }_0^2}{2}{x}^2},$$

$$C^{\mathrm{G}}(x)=\mu x+\frac{{\sigma }_0^2}{2}{x}^2.$$

$${\alpha }_{2n}^{\mathrm{G}}=\frac{2^n}{\sqrt{\pi }}\Gamma \left(n+\frac{1}{2}\right){\sigma }_0^{2n}=\frac{(2n)!}{2^{n}n!}{\sigma }_0^{2n},$$

$${\Delta }_P(\tau ):=\frac{P(\tau )-P^{\mathrm{G}}(\tau )}{P^{\mathrm{G}}(0)},$$

$${\Delta }_P(\tau )={\text{rme}}^{-{\tau }^2/2{\sigma }_0^2}\sum _{n=3}^{\infty }\frac{a_P(n)}{n!}{\mathrm{He}}_n\left(\frac{\tau }{{\sigma }_0}\right),$$

$$C(x)=\frac{{\sigma }_0^2}{2}{x}^2+\Delta C(x),$$

$$\Delta C(x):=\ln \left[1+\sum _{n=3}^{\infty }\frac{a_P(n){\sigma }_0^n}{n!}{x}^n\right]=\sum _{n=3}^{\infty }\frac{{\varkappa }_n}{n!}{x}^n.$$

$$C_n:=\frac{{\varkappa }_n}{{\sigma }_0^n}=\frac{⟨(\delta T{)}^n{⟩}_C}{{\sigma }_0^n}:$$

$$a_P(3)={\gamma }_1;a_P(4)={\gamma }_2;a_P(5)=C_5;a_P(6)=10{\gamma }_1^2+C_6;$$

$$a_P(7)=35{\gamma }_1{\gamma }_2+C_7;a_P(8)=56{\gamma }_1{C}_5+35{\gamma }_2^2+C_8;$$

$$a_P(9)=280{\gamma }_1^3+126{\gamma }_2{C}_5+84{\gamma }_1{C}_6+C_9;$$

$$a_P(10)=2100{\gamma }_1^2{\gamma }_2+120{\gamma }_1{C}_7+210{\gamma }_2{C}_6+126C_5^2+C_{10};$$

$$a_P(11)=5775{\gamma }_1{\gamma }_2^2+4620{\gamma }_1^2{C}_5+165{\gamma }_1{C}_8+330{\gamma }_2{C}_7+462C_5{C}_6+C_{11};$$

$$a_P(12)=15400{\gamma }_1^4+5775{\gamma }_2^3+27720{\gamma }_1{\gamma }_2{C}_5+$$

$$9240{\gamma }_1^2{C}_6+220{\gamma }_1{C}_9+495{\gamma }_2{C}_8+792C_5{C}_7+462C_6^2+C_{12}.$$

$$a_P(n)=B_n(0,0,{\gamma }_1,{\gamma }_2,C_5,\cdots ,C_n),$$

$${\int }_{-\infty }^{\infty }{\text{rme}}^{-{\nu }^2/2}{\mathrm{He}}_m(\nu ){\mathrm{He}}_n(\nu )\mathrm{d}\nu =\sqrt{2\pi }n!{\delta }_{mn},$$

Full text

1 General context

In his March 26, 2010 talk (Observational Constraints on Primordial Non–Gaussianity) during the ‘Non–Gaussian Universe’ workshop at the Yukawa institute (YITP), Eiichiro Komatsu [1] concluded about the statistical analysis of the Cosmic Microwave Background (CMB): “So far, no detection of primordial non–Gaussianity of any kind by any method”111http://wwwmpa.mpa-garching.mpg.de/~komatsu/talks.html .. This conclusion, strengthened by the analysis of the last available CMB data at the time ($WMAP$ 7yr)—and the reader may judge, after reading this paper and in the future, whether we can say more after $Planck$—sounded like a challenge to the cosmological community. An ongoing challenge, because the high–precision $Planck$ CMB data are certainly not exhaustively analysed today. A challenge also as there are as many analytical definitions of ‘non–Gaussianity’ as there are different statistical descriptors calling for the application of unified statistical analysis tools. But mostly a challenge, because in the frame of the $\Lambda$CDM model and various inflationary models, a non–detectable up to a non–negligible primordial non–Gaussianity may be expected [2, 3, 4], as measured by bi– and tri–spectra [5]. In various other models of primordial physics we may expect different kinds of primordial non–Gaussianity for the CMB as, e.g., in string gas cosmology [6], or even large non–Gaussianity such as in some ekpyrotic phase models [7]. Furthermore, it is a challenge given that a sufficiently high tensor–to–scalar ratio $r$ should allow for a slight detection of primordial gravitational waves with the non–Gaussianity of the CMB polarization (B–mode) and temperature maps correlation function ($\langle BTT\rangle$ bispectrum)222However, small angular scales (l$>$200) temperature anisotropies have to be analysed to reach the convergence power spectrum of the B–modes; unfortunately, at such scales the Sunyaev–Zel’dovich effect and reionization scattering pollute the temperature lensing reconstruction. [9, 10] and see (inter alia) the chapter “$7$–Lensing and the CMB” in the book by Durrer [8].

The observable part of the CMB covers, typically, depending on the frequency, up to 70$\%$ of the celestial sphere, and one has to assume what the properties of the remaining invisible CMB are. For this problem there is no model–independent analysis of the CMB—a model is required to infer and build the CMB regions hidden beyond the galaxy mask and beyond each field source out of the main foreground mask. It is interesting to see that some of the CMB anomalies reviewed in [11], such as the cold spot or the hemispherical asymmetry, are very likely independent of the masked CMB reconstruction model. But, most of the various CMB anomalies, even of this magnitude, are consistent with a certain level of Gaussianity (p–value $>99\%$). Some of these anomalies are detected with the 2–point correlation function; some anomalies could be remedied by a cutoff between 60 and 160*∘* as in [12], statistically favouring multi–connected universe models with finite volume [13, 14]. Besides assumptions on the topology of the Universe, an assumption is needed for the geometry of the support manifold, commonly thought of as being an ideal sphere. This is related to the estimate of the relative motion vector of the observer to the CMB. Will peculiar–velocity analyses of larger and larger catalogues converge to this motion, or is there a global dipole of a non–idealized space form? What is the correct definition of “peculiar–velocities”? There can be significant differential expansion of space that is not allowed for in a Newtonian model of structure formation [15]. Another challenge to the CMB non–Gaussianity is the way to link the specific CMB intensity, $I(\nu)$, to the CMB temperature anisotropy, $\delta T$ [16], using the Taylor expansion at first order only; higher orders may eventually impair the evaluation of the primordial non–Gaussianity. We shall explicitly address near–Gaussian expansions in the present paper.

Given the variety of potential contaminations in the cosmic microwave background map, a good strategy to analyse it and discriminate primordial non–Gaussianity from secondary effects is not only to multiply the statistical methods, but to head for model–independent estimators. For that, integral geometry provides us with a general mathematical framework where a small set of descriptors allows for a complete morphological analysis over random fields such as the CMB temperature maps. The descriptive power of integral geometry relies on the polynomial of convex bodies in three dimensions, introduced by J. Steiner (1840) [17], its generalization to the mixed volume associated to a convex body by Minkowski [18], and the Blaschke problem and diagram. Then, the Bonnesen enhanced isoperimetric inequality (1921), the Aleksandrov (Fenchel) inequalities (1937), the Hadwiger works and theorem (1955, 1957) [19], the studies by Santaló (1976) [20], the statistical predictions for random fields by Adler (1981) [21] and the work by Tomita (1986) [22, 23] bring the key mathematical foundations to the Minkowski Functionals (henceforth ‘MFs’). Interesting theoretical and applied developments are found in the mathematical reviews by Groemer [24], Schneider [25] and Mecke [26].

The theory of $2D$ Gaussian random fields on the CMB 2–sphere was developed by Bond and Efstathiou (1987) [27]. They applied it to the number density of hot and cold spots and to the ellipticity of peaks. This has later been generalized for extrema counts and ellipticity contour lines by Aurich et al. [28] and also by Pogosyan et al. [29, 30, 31].

Minkowski Functionals comprise the by now well–known set of scalar functionals, being rotation and translation invariant, Minkowski additive333Minkowski additivity assigns a functional of the union of bodies to the functionals of the individual bodies minus the functionals of their intersection., and conditionally continuous. This set of MFs describes the morphology of any convex body444Even more generally, the morphology of non–convex bodies is made possible using the property of Minkowski additivity, extending the analysis to the convex ring. in a complete and unique way in the sense of Hadwiger’s theorem. The explicit introduction of the MFs into cosmology (describing the morphology of galaxy distributions) was made in statistics of large–scale structure using the Boolean grain description by Mecke et al. (1994) [32] with follow–up studies of galaxy catalogues [33, 34, 35]; Kerscher wrote a review on the MFs including applications to cosmology [36]. In 1997, Schmalzing and Buchert introduced the MFs for the excursion set approach, suitable for any density or temperature contour maps [37]. In 1998, Schmalzing and Górski are the first to apply the set of $3+1$ MFs of the CMB 2–sphere (curvature=$+1$) to COBE DMR data excursion sets on a quadriteralized spherical cube tesselation ($6~{}[2^{2(N-1)}]$ pixels) [38], implementing also the Gaussian premises predicted in 1990 by Tomita [39]. (See also [40].) Also, galaxy catalogues have been analysed with the excursion set approach [41, 42], and a generalization to vector–valued MFs has been proposed [43] and applied to the morphological evolution of galaxy clusters [44].

At present, many fundamental tools are available to take up the challenge in the broad sense of several different but unified statistical descriptors to ask: is the cosmic microwave background Gaussian? Looking at CMB non–Gaussianity within the MF approach was first undertaken by Winitzki and Kosowsky [45], and by Novikov et al [46]. Following work by Takada et al. [47] on the detectability of the weak lensing with the two–point correlation function, a further study predicts that the weak lensing effect could be detected directly with the MFs ${\mathrm{v}}_{1}$ and ${\mathrm{v}}_{2}$ [48]. Not only this capability of the MFs was confirmed in later work, but also the lensing–induced morphology changes in modified gravity theories could be detected. Furthermore, the lensing and the Sunyaev–Zel’dovich non–Gaussianity can be separately detected with MFs, all of this for CMB temperature maps at high resolution (up to $\ell=3000$) [49]. Specific, regional morphological features like the Cold Spot could be detected by local analyses with MFs [50].

Since then the majority of works rely on model assumptions, either by testing a given inflationary model as in [46], or explicitly replacing the analytic Gaussian premises for the MFs by the MFs of a grid of $\Lambda$CDM model maps: the so–evaluated deviations from non–Gaussianity yield what we below call the ‘difference of the normalized MFs’ (abridged by ${\rm Df}$) (see [51, 52, 53, 54, 55, 56], to mention only a few works in this context). These authors argue that the noise, the mask and the pixel effects are better taken into account by using a model ensemble rather than using analytical predictions for evaluating the data map non–Gaussianity.

Ade et al. [56] propose a rather exhaustive (claimed model–independent) investigation of the CMB isotropy and statistics and neither reach a clear rejection nor a confirmation of the standard FLRW cosmological model. It is then a natural next step to investigate perturbative models at a FLRW background (here the works by Matsubara [57, 58, 59] stand out as a sustained such attempt). In this paper we follow another route. We focus on and specify model–independent methods to quantify non–Gaussianity, and we compare with what is obtained when applying the standard model–dependent perturbative ansatz. General perturbative expansions are based on series of terms, some of which are solely specified in terms of the chosen model [60, 61, 62], leaving a certain degree of arbitrariness in the application of perturbation theory. Beyond perturbation theory taken in this broad sense, we can fundamentally explore the non–Gaussianity only when using non–perturbative expansions. We shall introduce these latter paying careful attention to some methodological details that may have a strong impact on the extremely weak level of CMB non–Gaussianity.

This article is structured as follows: we investigate in section 2 specific descriptors of the map ensemble, such as the probability density function, section 2.1, and derived descriptors of non–Gaussianity such as the discrepancy functions and Hermite expansions, section 2.2, putting the first Minkowski Functional into perspective in section 2.3. For the Hermite expansion coefficients we give closed expressions in terms of the cumulants using the complete Bell polynomials. In section 2.4 we demonstrate in general terms the accuracy of Hermite expansions of the discrepancy functions that apply to any non–Gaussianities of arbitrary magnitude, and we show that the Hermite expansions approach the discrepancy functions with a mean square error as small as needed. We then demonstrate that the Hermite expansions can be truncated to yield the perturbative expansions with hierarchical ordering of the cumulants, usually considered in the literature, in section 2.5. We discuss our descriptors for the remaining morphological Minkowski Functionals in section 2.6, and proceed by comparing the discrepancy function approach with other approaches in the literature in section 2.7. We apply our descriptors to the $Planck$ $2015$ data with and without masks in section 2.8, and we discuss the results by addressing the issue of the origin of non–Gaussianities in section 2.9. Some conceptual problems related to the range of analysis are raised in section 2.10, while a short situational analysis of the $f_{\rm NL}$ results is given in section 2.11. In section 3 we give a short conclusion. Basic definitions and general remarks on uncertainties, the construction of model maps, and the problem of discretization are addressed in A. Some important mathematical facts about cumulants, the Hermite and Edgeworth expansions, and the derivation of closed expressions for the Hermite expansion coefficients are given in B and C.

2 Probability Density Function and Minkowski Functionals—Descriptors for CMB non–Gaussianity and Application to Planck 2015 Data

The CMB anisotropy $\delta T$ is described by a scalar random field on the surface of last scattering. The CMB support manifold is idealized by the constant curvature $2-$sphere ${\cal S}^{2}$, but different supports of variable curvature may be envisaged for the CMB. Furthermore, the thickness of the surface of last scattering is not taken into account. The manifold ${\cal S}^{2}$ is very convenient and universally adopted and we use it in the present work too. Upon this chosen manifold, any statistical descriptor such as the Minkowski Functionals is otherwise model–independent.

2.1 The probability density function of the CMB anisotropy

An important statistical descriptor of the random variable $\delta T$ is its probability density function (PDF) or frequency function $P(\tau)\geq 0$, where we assume that $\tau$ can take on any real value, and $P(\tau)$ is continuous. (The more general case will be discussed below.) $P(\tau){\mathrm{d}}\tau$ gives the probability of finding $\delta T$ in between $\tau$ and $\tau+{\mathrm{d}}\tau$. Hence, $P(\tau)$ is normalized to unity.

Figure 1 shows the $10^{5}$ histograms envelope of the individual PDFs in the $\Lambda$CDM sample over a total temperature range $\pm$396.5$\mu K$, divided into $61$ bins of $13\mu K$. Detailed informations regarding the $\Lambda$CDM map sample generation, the numerical methodology, and the conventions adopted for the analysis all along this work are given in A.

From $P(\tau)$ one obtains the cumulative distribution function $F(\tau)$ (CDF),

[TABLE]

respectively, the complementary cumulative distribution function $F_{C}(\tau)$ (CCDF)555It will turn out that $F_{C}(\tau)$ is identical to the Minkowski Functional ${\mathrm{v}}_{0}$ (see below).,

[TABLE]

satisfying $F(-\infty)=0$, $F(\infty)=1$, respectively, $F_{C}(-\infty)=1$, $F_{C}(\infty)=0$. Note that $F(\tau)$ is non–decreasing and ${\rm prob}(\tau_{1}\leq\tau<\tau_{2})=F(\tau_{2})-F(\tau_{1})$.

For a continuous function $f:{\mathbb{R}}\rightarrow{\mathbb{R}}$, $f(\delta T)$ is again a random field whose expectation value is defined as:

[TABLE]

if the integral exists. An important role is played by the moments $\alpha_{n}$ of $\delta T$,

[TABLE]

the central moments $m_{n}$ (i.e., the moments about the mean $\mu:=\alpha_{1}$),

[TABLE]

and the cumulants $\varkappa_{n}:=\left\langle(\delta T)^{n}\right\rangle_{C}$. The generating function of the moments $\alpha_{n}$ is given by:

[TABLE]

from which one obtains the generating function of the cumulants $\varkappa_{n}$,

[TABLE]

The first few moments and cumulants are $\alpha_{0}=m_{0}=1$, the mean values $\mu:=\alpha_{1}=\varkappa_{1}$, respectively, $m_{1}=0$, and the variance $\sigma_{0}^{2}$,

[TABLE]

respectively, the standard deviation (uncertainty) $\sigma_{0}=\sqrt{\alpha_{2}-\mu^{2}}$. There are the following recurrence relations:

[TABLE]

and, for $n\geq 2$:

[TABLE]

Any odd non–vanishing central moment $m_{2n+1}$ of the CMB anisotropy $\delta T$ is a measure of the skewness of $P(\tau)$; the simplest of these is

[TABLE]

respectively, the dimensionless skewness coefficient

[TABLE]

This latter indicates a possible asymmetry of $P(\tau)$, i.e., whether the left tail ($\gamma_{1}<0$) or the right tail ($\gamma_{1}>0$) is more pronounced. Another important measure of the “tailedness” of $P(\tau)$ is the dimensionless excess kurtosis (sometimes simply called kurtosis or excess),

[TABLE]

It should be stressed that the basic equations (1), (2), and (3) will in general only hold for the theoretical continuous limit distributions of the CMB anisotropy in the sense of ensemble averages in the limit of an infinite ensemble (infinitely many realizations). For a single realization, as it is the case for the data obtained by $WMAP$ or $Planck$, or in computer simulations using a large but finite number of realizations, the distributions will in general not be continuous but rather discrete and, thus, the PDF will not satisfy $P(\tau)={\mathrm{d}}F/{\mathrm{d}}\tau=-{\mathrm{d}}F_{C}/{\mathrm{d}}\tau$, as is implied by (1), respectively (2). Instead, the above Riemann integrals have to be replaced by Riemann–Stieltjes integrals such that, for example, the expectation value (3) is given by:

[TABLE]

The last relation even holds in cases where $F(\tau)$ is ill–behaved, for example, if $F(\tau)$ has at most enumerably many jumps at the discrete points $\tau_{\ell}$, and as long as $F(\tau)$ is a CDF. In this case one has (if the integral and the sum over $\ell$ converge):

[TABLE]

where $F^{\prime}(\tau)$ is the almost everywhere existing derivative of $F(\tau)$. An example is the case where $F(\tau)$ is given as a histogram, i.e., it is a step function and, thus, $F^{\prime}(\tau)\equiv 0$ almost everywhere.

The non–Gaussianities, which are the subject of this paper, are defined as deviations from the analytically known Gaussian prediction for the PDF and for the Minkowski Functionals that we shall introduce below. In the case of the PDF, the Gaussian prediction (G) is given by the normal distribution,

[TABLE]

which has mean $\mu$ and variance $\sigma^{2}_{0}$. From (16) one derives with (6) the moment generating function,

[TABLE]

and with (7) the generating function for the Gaussian cumulants,

[TABLE]

Our numerical results for the CMB anisotropy are based on an ensemble of $10^{5}$ realizations (for details, see A). It turns out that the mean value $\mu=\left\langle\delta T\right\rangle$ (calculated from pixels) is negligibly small, $\mu=\mathcal{O}(10^{-7}\mu K$) for sky maps that cover the full sky. Since typical values for the standard deviation $\sigma_{0}$ are $\sigma_{0}\cong 59\mu K$, the ratio $\mu/\sigma_{0}$ is much smaller than $\tau/\sigma_{0}$ and, thus, we can put $\mu=0$ in (16) when comparing with the full PDF. We then obtain from equations (6) and (17) the well–known result that all odd moments of the Gaussian prediction $P^{\rm G}(\tau)$ vanish, $\alpha^{\rm G}_{2n+1}=0$, and that the even moments are given by:

[TABLE]

and increase with increasing $n$. (Note that $m_{n}^{\rm G}=\alpha_{n}^{\rm G}$.) Equation (19) gives $\alpha_{4}^{\rm G}=3\sigma_{0}^{4}$ and, thus, one obtains from (12) and (13) $\gamma_{1}^{\rm G}=\gamma_{2}^{\rm G}=0$, which shows that a non–vanishing value of the skewness coefficient $\gamma_{1}$ and/or of the excess kurtosis $\gamma_{2}$ are quantitative measures of non–Gaussianity. A comparison of equation (18) with equation (7) shows that all Gaussian cumulants vanish apart from $\varkappa^{\rm G}_{2}=\sigma_{0}^{2}$ and, therefore, any non–vanishing cumulant with $n\geq 3$ is a measure of non–Gaussianity.

2.2 Discrepancy functions and Hermite expansions

In the present paper we pursue a general model–independent approach to the PDF and to the Minkowski Functionals (MF) that depends, in general, on all higher–order poly–spectra. As a measure of non–Gaussianity using the PDF, we consider the dimensionless discrepancy function $\Delta_{P}(\tau)$, defined by:

[TABLE]

with $P^{\rm G}(0)=(\sqrt{2\pi}\sigma_{0})^{-1}={\rm max}\{P^{G}(\tau)\}$. Here, $P(\tau)$ denotes the ensemble average over a large number of realizations ($10^{5}$ in our case), compatible with $\mu=\langle\delta T\rangle=0$, and possessing the standard deviation $\sigma_{0}$. The Gaussian prediction is defined by equation (16) for $\mu=0$, and by identifying the standard deviation with the value $\sigma_{0}$ that is numerically obtained from $P(\tau)$, i.e., $\alpha_{n}^{\rm G}=\alpha_{n}$ for $n=0,1,2$, respectively $\varkappa_{1}=0$ and $\varkappa_{2}=\sigma_{0}^{2}$.

Although some models of inflation predict large non–Gaussianities for $P(\tau)$, there is clear evidence from $WMAP$ and $Planck$ data, [51, 52, 55, 56], that possible deviations from the Gaussian prediction $P^{\rm G}(\tau)$ are very small. Under very general conditions (for details see C), $\Delta_{P}(\tau)$ can be written as a product of a Gaussian and a function $h(\nu)\in L^{2}(\mathbb{R},w(\nu){\mathrm{d}}\nu)$, where $w(\nu)$ denotes the weight function, $w(\nu)=\exp(-\nu^{2}/2)$, expressed in terms of the dimensionless scaled temperature variable $\nu:=\tau/\sigma_{0}\in\mathbb{R}$.

The Hermite polynomials ${\rm He}_{n}(\nu)$ provide a complete orthogonal basis in the Hilbert space $L^{2}(\mathbb{R},w(\nu){\mathrm{d}}\nu)$. Therefore, $h(\nu)$ possesses a convergent Hermite expansion (see C) and we are led to

[TABLE]

which describes the non–Gaussian “modulations” of the PDF. Possible non–Gaussianities are parametrized by the real dimensionless coefficients $a_{P}(n)$, where the non–vanishing of any of them is a clear signature of non–Gaussianity. We would like to point out that the ${\rm He}_{n}(\nu)$ are the “probabilist’s Hermite polynomials” that are different from the Hermite polynomials $H_{n}(\nu)$ commonly used in physics and which are defined with respect to the weight function $\exp(-\nu^{2})$. In the cosmology literature the ${\rm He}_{n}(\nu)$ are used, but unfortunately denoted as $H_{n}(\nu)$! The relation between the two is ${\rm He}_{n}(\nu)=2^{-n/2}H_{n}(\nu/\sqrt{2}$), (see, e.g., [63]). Our condition on $h(\nu)$ is satisfied as long as $P(\tau)$ and, thus, also $\Delta_{P}(\tau)$ are piecewise continuous, and if $h(\nu)=o(\exp(\nu^{2}/4)/\sqrt{|\nu|})$ for $|\nu|\rightarrow\infty$ (see C). It is clear that the expansion (21) is particularly useful if only a few terms have to be taken into account such that the series can be cut off at a low value $n=N$, i.e., can be well–approximated by a polynomial of degree $N$.

In order to obtain a physical interpretation of the non–Gaussianity (NG) parameters $a_{P}(n)$, we insert $P(\tau)=P^{\rm G}(\tau)+[1/(\sqrt{2\pi}\sigma_{0})]\Delta_{P}(\tau)$ into the definition (6) of the moment generating function $M(x)$, which in turn gives with (7) the following generating function of the cumulants $\varkappa_{n}$ of $P(\tau)$ (see B):

[TABLE]

with

[TABLE]

Then, the cumulants are $\varkappa_{1}=0$, $\varkappa_{2}=\sigma_{0}^{2}$, and the higher cumulants, for $n\geq 3$, are uniquely determined by (23). Here are the first coefficients $a_{P}(n)$, expressed in terms of the skewness coefficient $\gamma_{1}$, equation (12), the excess kurtosis $\gamma_{2}$, equation (13), and the dimensionless normalized cumulants,

[TABLE]

[TABLE]

Note that there is the general closed expression in terms of the normalized cumulants ($n\geq 3$):

[TABLE]

where $B_{n}(x_{1},x_{2},\cdots,x_{n})$ denotes the $n^{\rm th}$ complete Bell polynomial (see B). An equivalent closed expression in terms of the moments is given in equation (129).

Table 1 shows the moments $\alpha_{n}$ (main term, equation (111) in A.5), the cumulants $\varkappa_{n}$ and the dimensionless normalized cumulants $C_{n}$ of $P(\tau)$. We here give the three first orders of each for the sample without mask, $\alpha_{0,1,2}=1,~{}0\mu\mathrm{K},~{}3558.31519\mu\mathrm{K}^{2}$; $\varkappa_{0,1,2}=1,~{}0\mu\mathrm{K},~{}3558.31519\mu\mathrm{K}^{2}$ and $C_{0,1,2}=1,0,1$, and we obtain $\sigma_{0}=59.65166\mu\mathrm{K}$, $\gamma_{1}=-5.1872\times 10^{-4}$ and $\gamma_{2}=5.825\times 10^{-5}$.

With the help of the orthogonality relation (see [63], p.775, and C),

[TABLE]

one derives from (21) the following integral representation for the $a_{P}(n)$’s ($n\geq 3$):

[TABLE]

If the discrepancy function $\Delta_{P}$ is known, one can compute from (28) the non–Gaussianity parameters $a_{P}(n)$ and then, from (2.2) and (26), the normalized cumulants $\gamma_{1},\gamma_{2},C_{n}$. Vice versa, one can compute the $a_{P}(n)$’s from (2.2) once the cumulants have been computed from the moments $\alpha_{n}=m_{n}$ of the PDF $P(\tau)$, using equations (4) and (10), or directly from a theory of the primordial CMB fluctuations as, e.g., given by the ansatz (94). With ${\rm He}_{2n+1}(0)=0$ and ${\rm He}_{2n}(0)=[(-1)^{n}(2n)!]/2^{n}n!$, we obtain from (21) and (2.2) for the discrepancy function at $\tau=0$:

[TABLE]

In table 2 we present the values for the first $a_{P}(n)$’s (see equations (2.2) and (28)), respectively for $\gamma_{1}$, $\gamma_{2}$, and $C_{5}$ for the $\Lambda$CDM model (model and computational details are given in A).

One observes that $a_{P}(6)=10\gamma_{1}^{2}+C_{6}<0$ implying $C_{6}<0$ in agreement with table 1. This is in contrast to hierarchical ordering (as assumed in perturbation theory [58]) where $a_{P}^{\rm HO}(J,6)=6a_{0}^{\rm HO}(J,5)=10\gamma_{1}^{2}>0$ in second– and third–order $(J=2,3)$ (see equations (2.5) and (2.5)). Only at fourth and higher order (i.e. $J\geq 4$) it holds that $a_{P}^{\rm HO}(J,6)=a_{P}(6)$ (see (2.5)). This will be discussed more in detail in subsection 2.5.

Figure 2 displays the discrepancy function $\Delta_{P}$ of the averaged PDF ($10^{5}$ map sample). $\Delta_{P}$ is calculated by equation (20) over the lattice defined by the mid–points of each segment in the histogram of $P(\tau)$ using $\sigma_{0}$ defined as the “main term” of $\alpha_{2}$ in equation (110). This figure also displays the expansion in Hermite polynomials (equation (21) limited to the order $8$).

Figure 3 shows the envelope of the $10^{5}$ discrepancy functions of the map sample.

2.3 The first Minkowski Functional ${\mathrm{v}}_{0}$

Let us discuss now the simplest morphological descriptor which is given by the first Minkowski Functional (MF). We consider the compact excursion set ${\mathcal{Q}}_{\nu}\in{\mathcal{S}}^{2}$ with boundary $\partial{\mathcal{Q}}_{\nu}(\nu:=\tau/\sigma_{0})$:

[TABLE]

We define the MF $V_{0}(\nu)$ as follows:

[TABLE]

with ${\mathrm{d}}a$ denoting the surface element on ${\mathcal{S}}^{2}$. In the following we shall be using the normalized MF ${\mathrm{v}}_{0}(\nu)$ that is normalized with respect to ${\rm area}({\mathcal{S}}^{2})=4\pi$, i.e.:

[TABLE]

where $F_{C}$ is the complementary cumulative distribution function (2). From the definition (32) follows also the relation

[TABLE]

Thus, the MF ${\mathrm{v}}_{0}(\nu)$ can be seen as a “smoothed” (cumulative) version of the PDF $P(\tau)$. Using the Gaussian prediction $P^{\rm G}(\tau)$ in equation (16), we obtain the Gaussian prediction for ${\mathrm{v}}_{0}(\nu)$ (setting $\mu=0$):

[TABLE]

in terms of the complementary error function [63], from which one derives ${\mathrm{v}}_{0}^{\rm G}(-\infty)=1$, ${\mathrm{v}}_{0}^{\rm G}(\infty)=0$, ${\mathrm{v}}_{0}^{\rm G}(0)=1/2$, and the symmetry relation ${\mathrm{v}}_{0}^{\rm G}(-\nu)=1-{\mathrm{v}}_{0}^{\rm G}(\nu)$. Furthermore, one has the asymptotic behaviour for $\nu\rightarrow\infty$,

[TABLE]

In analogy to equation (20), we define the dimensionless discrepancy function $\Delta_{0}(\nu)$ as follows:

[TABLE]

which in turn is expanded into Hermite polynomials,

[TABLE]

with the coefficients

[TABLE]

The above coefficients $a_{0}(n),\ n\geq 2$, measure a possible non–Gaussianity described by the MF ${\mathrm{v}}_{0}(\nu)$. With

[TABLE]

and the relation (see equations (21) and (33)),

[TABLE]

one obtains the following relation between the coefficients $a_{0}(n)$ and $a_{P}(n)$:

[TABLE]

and, thus, we arrive at the following exact coefficients (see equations (2.2) and (26)):

[TABLE]

In table 3 we present the values for the first $a_{0}(n)$’s of the $\Lambda$CDM model.

Figure 4 shows the first Minkowski Functional ${\mathrm{v}}_{0}(\nu)$ of the $\Lambda$CDM sample without mask, and figure 5 its discrepancy function $\Delta_{0}(\nu)$ together with the Hermite expansion (37) of $\Delta_{0}(\nu)$ (order $2$ to $5$) in $a_{0}(n)$.

2.4 The accuracy of the Hermite expansion of the discrepancy functions

The Hermite expansions (21) and (37) for the discrepancy functions $\Delta_{P}(\tau)$ and $\Delta_{0}(\nu)$, respectively, provide very convenient parametrizations and quantitative measures of the CMB non–Gaussianities. They hold under very general conditions (as discussed in C) and do not require that the NGs have to be small. The main assumption is the existence of a unique probability density function $P(\tau)$ (a very natural assumption, indeed, from a physical point of view), which is equivalent to demanding that the “Hamburger moment problem” is determinate (see e.g. [64, 65]). It is then guaranteed that the Hermite expansions are absolutely convergent. In subsections 2.2 and 2.3 we considered polynomial approximations of degree $N$ in $\tau$, respectively in $\nu$ of $\Delta_{P}(\tau)$ and $\Delta_{0}(\nu)$ by truncating the Hermite expansions at $n=N$. A nice property of this truncation is that the accuracy of the approximation is well under control, since there is, for example for $\Delta_{0}(\nu)$, the exact mean square error of equations (136) and (137),

[TABLE]

for $h_{0}(\nu):=e^{\nu^{2}/2}\Delta_{0}(\nu)$. The error ${\cal{E}}_{N}$ gets smaller and smaller if the degree $N$ increases and finally approaches zero in the limit $N\rightarrow\infty$ as a consequence of the completeness relation (Parseval’s equation),

[TABLE]

(here we use the weight function $w(\nu)=\frac{1}{\sqrt{2\pi}}e^{-\nu^{2}/2}$ in the definition of the inner product of the Hilbert space ${\rm I\!H}=L^{2}({\rm I\!R},w(\nu){\mathrm{d}}\nu)$, see C.) As an example, let us consider the approximation of degree $N=5$ applied to $\Delta_{0}(\nu)$ for which the mean square error is explicitly given in terms of the skewness coefficient $\gamma_{1}$, the excess kurtosis $\gamma_{2}$ and the higher cumulants $C_{5}$ and $C_{6}$ by

[TABLE]

using the coefficients $a_{0}(2,3,4,5)$ of (2.3). Figure 5 nicely illustrates that the polynomial approximation of degree 5 provides an excellent description of the discrepancy function $\Delta_{0}(\nu)$. Similar results exist for the discrepancy function $\Delta_{P}(\tau)$ as shown already in figure 2 at order $8$.

2.5 Hierarchical ordering and perturbation theory

Most of the commonly studied models of inflation predict weak primordial non–Gaussianities. Among these models there is a wide class where the normalized cumulants $C_{n}$ obey the additional property of hierarchical ordering (indexed below by HO), i.e., $C_{n}\sim\sigma_{0}^{n-2},n\geq 2$ (where we assume $C_{1}=0,C_{2}=1$). This suggests to expand the discrepancy function not into Hermite polynomials, but rather at fixed $\nu$ into a power series (perturbation theory) in $\sigma_{0}$. This approach has been pioneered by Matsubara [57, 58] who has derived second–order perturbative formulae in $\sigma_{0}$ for the discrepancy functions of the MFs. In the following, we show that the perturbation theory applied to $\Delta_{0}(\nu)$ is a direct consequence of the general Hermite expansion (37) that can be easily carried out to any order $J$ in $\sigma_{0}$, i.e., to the NG of order $\sigma_{0}^{J}$ (indexed below by ${\rm HO}(J),J\geq 1$). It turns out that the perturbation expansion ${\rm HO}(J)$ corresponds, at a given order $J$, to a double truncation of the Hermite expansion: i) a truncation of the Hermite expansion at $n=M_{0}(J):=3J-1$, ii) a truncation of the expansion coefficients $a_{0}(n)$ at order $J$ with respect to their expansion into a power series in $\sigma_{0}$. As an example, we give below the explicit formulae at second, third and fourth order, where the second–order formula is identical to Matsubara’s result. Analogous perturbative formulae hold for $\Delta_{P}(\tau)$, which have not been given before. As another new result we also present the exact mean square error of the perturbative formulae and compare it with the corresponding error (43) of the original Hermite expansion.

Let us assume that the hierarchical ordering holds and introduce the “renormalized cumulants” (or “HO–cumulants”),

[TABLE]

which are assumed to be of zeroth order in $\sigma_{0}$. Using the relation (41) and the explicit expression (26) for the coefficients $a_{P}(n)$ in terms of the complete Bell polynomials $B_{n}$, we obtain for $n\geq 2$,

[TABLE]

where $D_{n}(z)$ is a polynomial of degree $n$ in $z$ with $D_{0}(z):=1$ and

[TABLE]

Thus, $a_{0}(n)$ is a polynomial of degree $n-1$ in $\sigma_{0}$,

[TABLE]

Here, $j_{0}(n)$ denotes the smallest power of $\sigma_{0}$ which contributes to $a_{0}(n)$ and is given for $n\geq 2$ by

[TABLE]

where $n$ is parametrized as $n=n(k,l):=k+3l$ with $k=2,3,4$ and $l=0,1,2,\cdots$. (As a consequence of the symmetry relation (48), the coefficients $D_{n-1,j}$ vanish for $j$ even or odd depending on whether $n$ is even or odd.) The first few polynomials are explicitly given by:

[TABLE]

(Note that the highest power in $\sigma_{0}$ is for all $D_{n}$ given by $\frac{S_{n+2}}{n+2}z^{n}$ for $n\geq 1$.) Explicit expressions for the polynomials $D_{n}(z)$ for $n\geq 7$ are easily obtained either from the recurrence relation (122) or from the combinatorial expression (123).

In order to derive the perturbative expansion ${\rm HO}(J)$ for the discrepancy function $\Delta_{0}$ at any order $J\geq 1$, we insert in the Hermite expansion (37) the polynomial relation (47) for the expansion coefficients $a_{0}(n)$. At lowest order, $J=1$, one immediatly obtains from (47) and (49–2.5) the simple result:

[TABLE]

which is completely determined by the skewness coefficient $\gamma_{1}=S_{3}\sigma_{0}$. To obtain the ${\rm HO}(J)$–expansion for $J\geq 2$, we decompose the exact (general) Hermite expansion into three terms ($M_{0}(J)=3J-1$):

[TABLE]

Since $D_{n-1}(\sigma_{0})$ is a polynomial of degree $n-1$ in $\sigma_{0}$, the first sum in (53) represents (having $\nu$ fixed) a polynomial of degree $J$ in $\sigma_{0}$ and thus contributes to the ${\rm HO}(J)$–expansion. The last infinite series in (53), which is absolutely convergent, does not contribute at all to the ${\rm HO}(J)$–expansion, since the smallest power of $\sigma_{0}$ appearing in this series is already larger than $J$; (with $M_{0}(J)+1=n(3,J-1)$, equation (50) gives $j_{0}(M_{0}(J)+1)=J+1>J$.) We thus obtain the important result that the perturbative expansion at order $J$ necessarily implies a truncation of the Hermite expansion at $n=M_{0}(J)=3J-1$. It remains to discuss the second finite sum in (53) where the summation runs over $J+2\leq n\leq M_{0}(J)$. Inspection of the expression (49) for the polynomials $D_{n-1}(\sigma_{0})$ (see also the explicit expressions (2.5)) shows that they will contribute (for $n\geq J+2$ and $j_{0}(n)\leq J$) to this sum not only with powers $j\leq J$, but also with higher powers that are not admitted in the ${\rm HO}(J)$–expansion. Thus, the polynomials $D_{n-1}$ have to be replaced by truncated ones. In the case where $j_{0}(n)>J$, the polynomials $D_{n-1}$ do not contribute at all. This leads us to define, for $n\geq J+2$, the truncated polynomials $D_{n-1}^{{\rm HO}(J)}(z)$, which now depend also on $J$:

[TABLE]

We now explicitly give the truncated polynomials at order $J=2,3$ and $4$:

$J=2,M_{0}(2)=5$:

[TABLE]

$J=3,M_{0}(3)=8$:

[TABLE]

$J=4,M_{0}(4)=11$:

[TABLE]

We then obtain from (53) the general perturbative formula for the discrepancy function $\Delta_{0}(\nu)$ valid at any order $J\geq 1$:

[TABLE]

in terms of the $\rm HO$–Hermite expansion coefficients $a_{0}^{\rm HO}(J,n)$, which now also depend on $J$,

[TABLE]

As an example, we give the expansion coefficients of the hierarchical ordering at second, third and fourth order:

$J=2,M_{0}(2)=5$:

[TABLE]

$J=3,M_{0}(3)=8$:

[TABLE]

$J=4,M_{0}(4)=11$:

[TABLE]

(Here, we have replaced the $S_{n}^{\prime}s$ by the cumulants $C_{n}$ according to equation (46).)

Being a direct consequence of the general Hermite expansion (37), the perturbative formula (58), valid at arbitrary order $J$ in $\sigma_{0}$, has still the form of a Hermite expansion (truncated at $n=M_{0}(J)$). If we insert, however, for the Hermite coefficients $a_{0}^{\rm HO}(J,n)$ their definition (59) in terms of polynomials in $\sigma_{0}$, we obtain, by combining all terms of same power, the alternative version of the perturbative formula:

[TABLE]

which has now (at fixed $\nu$) the form of a power series in $\sigma_{0}$ with coefficient functions ${\mathrm{w}}_{0}^{{\rm HO}(J)}(j,\nu)$. (Since the series in (58) is finite, the rearrangement as a power series in $\sigma_{0}$ is of course always possible.) The coefficient functions ${\mathrm{w}}_{0}^{{\rm HO}(J)}(j,\nu)$ are given as a linear combination of a finite number of Hermite polynomials. Up to the second order, they have been calculated by Matsubara [57, 58]. It is straightforward to obtain them at arbitrary order $J$ using the equations (49), (54) and (59). Here we give the coefficient functions up to fourth order:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The second–order formulae (65) agree with Matsubara’s result [57, 58] (who writes $S:=S_{3}$ and $K:=S_{4}$).

It is worthwhile to mention that the perturbative formula (63) closely resembles the so–called Edgeworth expansion (150) which gives a refinement of the classical central limit theorem (for details and references, see C). However, (150) is an asymptotic expansion in the sense of Poincaré where the role of $\sigma_{0}$ is played by the small dimensionless parameter $1/\sqrt{n}$ which can be made arbitrarily small and does not have a fixed finite value as in the case of the standard deviation $\sigma_{0}$ of the CMB anisotropy. In fact, the central limit theorem is precisely the statement that the limit $n\rightarrow\infty$ is asymptotically exactly Gaussian and thus the NGs in the Edgeworth expansions have no fundamental meaning, they just determine the rate of convergence to the Gaussian limit. In contrast, the NGs of the CMB—if they are non–zero and of primordial origin—contain genuine information on the underlying model of inflation.

Having shown that the two perturbative formulae (58) and (63) are identical, we discuss in the following only the Hermite expansion (58).

Let us compare the perturbative NG–coefficients, (2.5) respectively (2.5), with the complete NG–coefficients $a_{0}(n)$ given in equation (2.3). At second order ($J=2$) we see that the first two coefficients are identical to the complete coefficients $a_{0}(n)$, i.e., $a_{0}^{\rm HO}(2,n)=a_{0}(n)$ for $n=2,3$. The coefficient $a_{0}^{\rm HO}(2,4)$ vanishes because its complete value $C_{5}/5$ is of third order. Finally, the last coefficient $a_{0}^{\rm HO}(2,5)$ differs from $a_{0}(5)$ by the term $C_{6}/6$, which is of fourth order. At third order ($J=3$) one observes that now the first three coefficients are identical to the complete expressions, i.e., $a_{0}^{\rm HO}(3,n)=a_{0}(n)$ for $n=2,3,4$; $a_{0}^{\rm HO}(3,5)=a_{0}^{\rm HO}(2,5)$ still does not contain the term $C_{6}/6$; the coefficient $a_{0}^{\rm HO}(3,7)$ vanishes because the complete value for $a_{0}(7)$ is of order $\sigma_{0}^{4}$; $a_{0}^{\rm HO}(3,6)$ does not contain the term $C_{7}/7$ (of order $5$), and $a_{0}^{\rm HO}(3,8)$ does not contain the terms of order $5$, respectively $7$. Note, in particular, that a vanishing coefficient $a_{0}^{\rm HO}(J,n)$ implies that the associated contribution from the Hermite polynomial ${\rm He}_{n}(\nu)$ is absent in $\Delta_{0}^{\rm HO(J)}$ compared to $\Delta_{0}$. This is the case at second order with ${\rm He}_{4}(\nu)$, and at third order with ${\rm He}_{7}(\nu)$. Since the ${\rm He}_{n}(\nu)$’s have exactly $n$ distinct zeros, the omission of one or several of them can have an important influence on the shape of the discrepancy function.

Finally, we come to the important question about the accuracy of the perturbative expansion. To this purpose we consider in analogy to equation (43) the mean square error of the ${\rm HO}(J)$–expansion:

[TABLE]

where we have used the identity (137). Here, the first two terms are identical to the error ${\cal{E}}_{M_{0}(J)}$, i.e. the error of the complete Hermite expansion truncated at order $N:=M_{0}(J)$ (see equation (43)), which leads to the interesting formula:

[TABLE]

Here, we have used in the last sum $a_{0}(n)-a_{0}^{\rm HO}(J,n)\equiv 0$ for $2\leq n\leq J+1$ (see (59)). Since the last sum in (69) is strictly positive, one infers that we have at any order in perturbation theory:

[TABLE]

i.e., the mean square error of the ${\rm HO}(J)$–expansion is always larger than the mean square error of the complete Hermite expansion truncated at $n=M_{0}(J)$. For instance, at second–order of perturbation theory, $J=2$, we have ${\cal{E}}^{{\rm HO}(2)}>{\cal{E}}_{5}$, where ${\cal{E}}_{5}$ is explicitly given in (45). Precisely, we obtain from (69) through (2.3) and (2.5):

[TABLE]

One observes that both errors approach zero in the limit $J\rightarrow\infty$ which reflects the fact that the perturbative expansion becomes identical to the general untruncated Hermite expansion in this limit.

Test of the hierarchical ordering in perturbation theory

In figure 6 we compare the discrepancy function $\Delta_{P}$ with its Hermite expansion in terms of the coefficients $a_{P}(n)$ (computed from (2.2) and table 1) for $n=3$ to $6$, according to the following equation:

[TABLE]

In addition, in figure 6, we show the expansion of $\Delta_{P}$ with the assumption of hierarchical ordering (HO) in fourth–order perturbation theory, according to

[TABLE]

where the NG–coefficients $a_{P}^{\rm HO}(4,n)$ are obtained from equations (41) and (2.5) using $\gamma_{1}$, $\gamma_{2}$ and the cumulants $C_{n}$ in table 1.

We verified that, at second and third–order, the HO–expansion displays an almost purely odd function (central symmetry), while $\Delta_{P}$ (black dashed line) does not pass by the central point ($\tau=0$, $\Delta_{P}(0)$, see equation (29)), and we note that $\Delta_{Px}(\tau)$ derived from moments of pixels, not shown here, confirms this behaviour.

At fourth–order, however, the HO–expansion fits well the non–Gaussianity of the sample and reveals the same shift from perfect central symmetry. (For some further discussion on the perturbative model with hierarchical ordering, we refer the reader to C).

2.6 The Minkowski Functionals ${\mathrm{v}}_{1}(\nu)$ and ${\mathrm{v}}_{2}(\nu)$

According to Hadwiger’s theorem [19], there are three independent MFs on ${\mathcal{S}}^{2}$. The simplest case $V_{0}$, respectively ${\mathrm{v}}_{0}$ has already been discussed in subsection 2.3. In this section, we discuss the remaining two which are defined again with respect to the excursion set ${\mathcal{Q}}_{\nu}$, see equation (30), and are given in their normalized form by

[TABLE]

where ${\mathrm{d}}s$ denotes the line element along $\partial{\mathcal{Q}}_{\nu}$, and $\kappa(s)$ the geodesic curvature of $\partial{\mathcal{Q}}_{\nu}$. The Gaussian predictions for the MFs, ${\mathrm{v}}_{k}(\nu),k=1,2$, have been computed by Tomita [22, 23, 21] (we set $\mu=0$):

[TABLE]

Here, $\sigma_{0}$ denotes again the standard deviation of the CMB temperature anisotropy $\delta T({\boldsymbol{\hat{n}}})$ on ${\mathcal{S}}^{2}$, where ${\boldsymbol{\hat{n}}}={\boldsymbol{\hat{n}}}(\vartheta,\varphi)$ denotes a unit vector on ${\mathcal{S}}^{2}$ dependent on the coordinates $x^{1}=\vartheta\in[0,\pi]$ and $x^{2}=\varphi\in[0,2\pi]$. Then, the line element on ${\mathcal{S}}^{2}$ is given by ${\mathrm{d}}s^{2}=\gamma_{ij}{\mathrm{d}}x^{i}{\mathrm{d}}x^{j}$ with $\gamma_{11}=1$, $\gamma_{22}=\sin^{2}\vartheta$, $\gamma_{ij}=0$ otherwise, and $\gamma_{ik}\gamma^{kj}=\delta_{i}^{\ j}$. Furthermore, $\sigma_{1}^{2}$ is the variance of the gradient field $\boldsymbol{\nabla}\delta T=(\nabla^{1}\delta T,\nabla^{2}\delta T)$, i.e.,

[TABLE]

From the MFs ${\mathrm{v}}_{1}$ and ${\mathrm{v}}_{2}$ one can form, by a linear combination, two further interesting measures, the Euler characteristic $\chi(\nu)$, respectively the genus $g(\nu):=1-\frac{1}{2}\chi(\nu)$. On the excursion set ${\mathcal{Q}}_{\nu}$ with smooth boundary $\partial{\mathcal{Q}}_{\nu}$ the Gauss–Bonnet theorem holds ($K=1/R^{2}\equiv 1$ is the Gaussian curvature on the unit sphere with radius $R=1$):

[TABLE]

which gives

[TABLE]

As a measure of possible non–Gaussianities based on the MFs ${\mathrm{v}}_{k}(\nu),k=1,2$, we define, in analogy to the discrepancy functions $\Delta_{P}(\nu)$ (equation (21)) and $\Delta_{0}(\nu)$ (equation (37)), the following discrepancy functions (we here include also the case $k=0$):

[TABLE]

with $N_{0}=1/\sqrt{2\pi}$, $N_{1}={\rm max}\{{\mathrm{v}}_{1}^{\rm G}\}=[1/(8\sqrt{2})]\sigma_{1}/\sigma_{0}$, and $N_{2}=[1/(2(2\pi)^{3/2})]\sigma_{1}^{2}/\sigma_{0}^{2}$. Under the same assumptions made before (for $\Delta_{P}$ respectively $\Delta_{0}$, see C), we can expand the $\Delta_{k}$’s into a convergent Hermite expansion,

[TABLE]

with the dimensionless NG–coefficients $a_{k}(n)$ and $n_{0}=2$, $n_{1}=n_{2}=0$. With the help of (27) we obtain the integral representation:

[TABLE]

Again, it is assumed that the series (81) can be truncated at a low value $n=N$, where $N$ may depend on $k$. In the following, we shall again compare the general expansion (81) with the one derived under the assumption of hierarchical ordering ($k=0,1,2$),

[TABLE]

for which the second order NG–coefficients $a_{k}^{\rm HO}(2,n)$ have also been calculated by Matsubara [58, 59] for $k=1$ and $k=2$. The highest degree of the Hermite polynomials contributing for $J=2$ is given by $M_{1}(2)=6$ and $M_{2}(2)=7$.

The NG–coefficients are then given for $\Delta_{1}^{{\rm HO}(J)}(\nu)$ as follows:

$J=2,M_{1}(2)=6$:

[TABLE]

and for $\Delta_{2}^{{\rm HO}(J)}(\nu)$:

$J=2,M_{2}(2)=7$:

[TABLE]

Here, we introduced the three dimensionless skewness parameters $\gamma_{1}$, $S_{1}$, $S_{2}$ and the four dimensionless kurtosis parameters $\gamma_{2}$, $K_{1}$, $K_{2}$ and $K_{3}$ using the notation $\tau=\tau(\boldsymbol{\hat{n}}):=\delta T(\boldsymbol{\hat{n}})$ ($\gamma_{1}$ and $\gamma_{2}$ as in equations (12) and (13), respectively):

[TABLE]

(Note that the products of the field $\tau(\boldsymbol{\hat{n}})$ respectively of its derivatives are taken at the same point $\boldsymbol{\hat{n}}$ on ${\mathcal{S}}^{2}$.)

Table 4 shows the coefficients $a_{1}(n)$ and $a_{2}(n)$ calculated from (82).

Figure 7 shows the second Minkowski Functional, and figure 8 the discrepancy function $\Delta_{1}$ together with the Hermite expansion to order $8$. The third Minkowski Functional is shown in figure 9, and the discrepancy function $\Delta_{2}$ together with the Hermite expansion to order $8$ in figure 10.

2.7 Discrepancy functions and ${\rm Df}-$differences

Gaussian random fields have a specific signature (the Gaussian prediction) depending only on the choice of the descriptor. Non–Gaussian processes may generate strong departures from the Gaussian prediction as in the formation of large–scale structure; however, attempts to find general and specific analytic signatures of a statistical property sufficiently far away from Gaussianity are most of the time unsuccessful in the context of CMB analyses. It is clear that the values $\sigma_{0{\rm C\ell}}$ and $\sigma_{1{\rm C\ell}}$ are model–dependent and that their use in the formulae for the Gaussian prediction, equations (16), (35) and (76), biases the reference of Gaussianity in general, so that the $\sigma-$values from the moments of pixels (denoted by subscripts $\rm px$, $\rm x$) or from the moments of the PDF should rather be used in the discrepancy functions $\Delta_{k}(...)$ we defined above. As a reminder, we list here the whole set of discrepancy functions:

[TABLE]

These discrepancy functions are self–consistent and model–independent given that the terms used, $\sigma_{\rm S}$ and $\sigma_{\rm Spx}$ are model–independent; these four formulae ($\Delta_{k}(...)$) are straightforwardly applicable to a single data or sample map. In the case of determination for a sample $\rm S$ of maps, we define the sample discrepancy functions this way:

[TABLE]

As mentioned at the end of section 1, the $Planck$ collaboration applies a different method of calculation for the non–Gaussianity. While the analytic formulae are not given in the various papers [51, 52, 53, 54, 55, 56], we can rebuild them here from the explicit formulations in the text. In the case of a single data map $d$, the differences of the normalized MFs (denoted ${\rm Df}$) are given by

[TABLE]

where the reference map sample $\mathrm{S}$ is Gaussian by construction. Compared with the discrepancy functions, the ${\rm Df}$’s measure the distance to a Gaussian premise which may be different in amplitude and shape from our analytical one. The ${\rm Df}$’s may, under certain conditions, provide a direct statistical measure of the departure from a given cosmological model statistics, which is an interesting feature.

Let us make some remarks on this methodology addressing some important issues. The variance terms $\sigma_{S}$, entering the left and right denominators of equation (93), may differ in the sample map generation, if no special precaution is taken. We have compared the results of this ${\rm Df}-$method (using Gaussian $\Lambda$CDM map samples) with our discrepancy functions method, using the analytical Gaussian premises, and found no significant difference of non–Gaussianity. But, obviously, the dependence on the statistical properties of the map sample and the number of maps may yield different results. It is not clear how the non–Gaussianity of the reference map sample is calculated as the direct application of the formula (93) gives zero and in this case, the analytic Gaussian premise is probably used. If instead our discrepancy functions are applied to the sample, we do not employ different methods for comparing the non–Gaussianity of a CMB data map with a CMB map sample.

2.8 Discrepancy functions of the $\Lambda$CDM sample maps versus $Planck$ maps

The next figures, figure 11, figure 12, figure 13 and figure 14 show the statistical and morphological behaviour of each of the $10^{5}$ sample maps in the $\Lambda$CDM model compared with the four $Planck$ maps $NILC$, $SEVEM$, $SMICA$ and $Commander$—$Ruler$, all with the $U73$ mask and a bin width of $6\mu\mathrm{K}$. The average quantities $\mu$, $\sigma_{0}$ and $\sigma_{1}$, as well as the discrepancy functions are calculated over the largest common and centered (on $\nu=0$) temperature range $[-201\mu\mathrm{K},+201\mu\mathrm{K}]$, as given by the $Planck$ 2015 $SMICA$ map. In these figures, clearly, the non–Gaussianity of the average sample is weak. For each descriptor, many $\Lambda$CDM sample maps show a strong non–Gaussianity; several individual sample maps lying well beyond the 5$\sigma$ cosmic variance limit.

The benchmark study of $Planck$ maps and the $\Lambda$CDM sample with U73 mask (figures 11 to 14) compares all the different maps over the same temperature range ($\pm 201\mu K$) for a $6\mu K$ binning, and shows very populated (almost saturated) $4\sigma$–sample envelopes for each of the four discrepancy functions. In the case of $\Delta_{P}$ we count $240$ maps (among $100000$) peaking beyond the $6\sigma$–envelope. Peaking beyond $6\sigma$ we also count $354$ maps for $\Delta_{0}$.

2.9 Discussion: origin of non–Gaussianities

Contrary to other works to unveil CMB non–Gaussianity, no specific model of non–Gaussianity has been put into the $100000$ maps of the $\Lambda$CDM sample probed all along this study (see A). However, already without mask, as in the initial analysis of the sample explored in this paper, we find that the four discrepancy functions show small but clear departures from Gaussianity, and this for a map ensemble that is supposed to be highly if not completely Gaussian, since it is computed from the $\Lambda$CDM model using purely Gaussian initial conditions. What is the origin of these NGs? For a given ensemble of realizations one expects several extraneous mechanisms that generate supplementary NG such as: an increasing smoothing scale beyond the $1^{\circ}$ horizon angular scale on the surface of last scattering increases the NG amplitude (increasing the dependence between neighbouring pixels acts in the sense of enlarging the causal horizon radius); we also checked that smoothing leaves the shape of the discrepancy functions unaffected, while rising the amplitude of NG. Our present study shows that these numerical NGs are of non-negligible amplitude even for a $100000$ map sample without mask. Our previous analyses using different parameters $N_{\rm side}$, $l_{\rm max}$, $\rm fwhm$, without or with mask leave the shape of the descriptors $\Delta_{1}$ and $\Delta_{2}$ invariant, but with different amplitudes (we have no such invariance for $\Delta_{P}$ and $\Delta_{0}$, they detect the mask). We reached the same conclusion using different cosmological parameters (e.g. $WMAP$ 7yr) for the generation of the power spectrum to get a $100000$ map sample. The comparison with the NG of the observed CMB (figures 11 to 14) shows that numerous individual $\Lambda$CDM maps can possibly be of the same amplitude and shape of NG than the $Planck$ data map. These assessments would say that without any addition of supplementary NGs in each sample map, some sources of non–Gaussianity existing in the power spectrum used to generate the maps may well be detected by the discrepancy functions. The census of components is in A.3. Among them, probably some cannot be detected because we limit the *l–*range to $[2,256]$, but the Sachs–Wolfe effect, the Doppler effect and the integrated Sachs–Wolfe effect would contribute to the non–Gaussianity detected in the individual maps of the ensemble. A systematic study is needed to disentangle the NG effects of the various components.

2.10 Discussion: the case of $\sigma_{0}$

When $\mu$, $\sigma_{0}$, and $\sigma_{1}$ are evaluated over $[-\tau_{m},\tau_{m}]$, while the discrepancy functions are calculated over a smaller interval in $\tau$ (resp. in $\nu$), this leads to an error on their shape and to an overestimate of the magnitude of NG, as we tested for the $100000$ map sample. Such a problem may arise when one works with the normalized temperature $\nu=\tau/\sigma_{0}$. While $\sigma_{0}(\tau_{max1})$ is known for the temperature range $\tau_{max1}$ giving $\nu_{max1}=\tau_{max1}/\sigma_{0}(\tau_{max1})$, the need to impose a $\nu$–range ($\nu_{max2}$) different from $\nu_{max1}$ (e.g. for comparing with another map sample) is equivalent to solve for an implicit function $f$, $\sigma_{0}(\nu_{max2})=f(\sigma_{0}(\nu_{max2}))$, simply because two different map samples or two different maps will have in general different variances and variances of the gradient, even if they are explored over the exactly same temperature range. The estimation of the unknown $\sigma_{0}$ corresponding to the new $\nu$–range can be made using iterative methods, or by referring to a model that predicts the behaviour of $\sigma_{0}$ for given changes of the range. For the comparison (with U73 mask) of $Planck$ maps to the $\Lambda$CDM sample, $P(\tau)$ and $\Delta_{P}(\tau)$ show no difficulties as abscissas are in $\tau$ and impose the same range to data and to the simulation maps.

For $\Delta_{0,1,2}$ we first calculated $\mu$’s and $\sigma$’s over the full $\tau$–range ($\pm$ 201 $\mu\mathrm{K}$) of the four different $Planck$ maps. For each of these we obtain the following $\sigma_{0}$ values: $NILC$ 51.532$\mu\mathrm{K}$, $SEVEM$ 51.750$\mu\mathrm{K}$, $SMICA$ 51.576$\mu\mathrm{K}$, $Commander$—$Ruler$ 51.794$\mu\mathrm{K}$, while we obtain $\sigma_{0}=$59.294$\mu\mathrm{K}$ for the $\Lambda$CDM simulations over exactly the same ($\pm$ 201 $\mu\mathrm{K}$) range! This illustrates well the “anomalously” low variance of $Planck$ data already observed with $WMAP$ compared to a $\Lambda$CDM model map ensemble (see [66, 67, 68, 69, 70, 71, 11] and [56, 72]). When not treated correctly, the impact of this anomaly is significant for the calculation of the Gaussian premises $P^{\rm G}$, equation (16), and ${\mathrm{v}}^{\rm G}_{k}$, equations (34) and (76), as they are functions of the standard deviation $\sigma_{0}$. Regarding ${\mathrm{v}}_{1}$ and ${\mathrm{v}}_{2}$ we made two interesting observations: firstly, the Gaussian premises ${\mathrm{v}}^{\rm G}_{1}$ and ${\mathrm{v}}^{\rm G}_{2}$ are functions of the ratio $\rho=\sigma_{1}/\sigma_{0}$ in the prefactor of the exponential, and we verified that this ratio is very similar for the four $Planck$ maps, $\rho=36.657\pm 0.031$, even if the $\sigma_{0}$’s are not very similar, $\sigma_{0}=51.663\pm 0.131\mu\mathrm{K}$. For the simulation sample, $\rho_{\Lambda CDM}$ has a much smaller value: $34.162$. Secondly, we observed that the ratio $\rho$ is very stable for the sample with U73 mask and for $Planck$ maps with U73 mask when passing from the largest common temperature range ($\pm$ 201 $\mu\mathrm{K}$) to the smallest temperature range covering all the maps ($\pm$ 396.5 $\mu\mathrm{K}$). On the other hand we noticed (see at the beginning of A) that the ratio $\rho$ is no more stable but increasing when passing from no mask to U73 mask.

In summary, in the present work, the issue of the fair comparison of different CMB maps, (i) is treated by using exactly the same temperature range for $P$, $P^{\rm G}$ and $\Delta_{P}$, (ii) for ${\mathrm{v}}_{0}$, ${\mathrm{v}}^{\rm G}_{0}$ and $\Delta_{0}$ it required some iterations upon different temperature ranges to reach close to the $\nu_{\rm max}$ values, and (iii) is simplified for ${\mathrm{v}}_{1}$, ${\mathrm{v}}^{\rm G}_{1}$, ${\mathrm{v}}_{2}$, ${\mathrm{v}}^{\rm G}_{2}$, $\Delta_{1}$ and $\Delta_{2}$ by assuming $\rho=\sigma_{1}/\sigma_{0}$ to be constant for the given change of temperature range.

2.11 Discussion: $f_{NL}$ results

Non–linear correction terms within the standard perturbation approach are commonly investigated by constraining the coefficients (bi–spectrum with $f_{\rm NL}$ and tri–spectrum with $g_{\rm NL}$). Applied to Bardeen’s curvature these constraints allow to decide whether inflation models such as with single slow–roll scalar field are rejected or not. A scalar random field of the primordial period, the gravitational potential $\Phi({\boldsymbol{\hat{n}}})$ (${\boldsymbol{\hat{n}}}=\boldsymbol{\hat{n}}(\vartheta,\varphi)$), is commonly expanded in real space around a Gaussian random field $\phi({\boldsymbol{\hat{n}}})$ (of vanishing mean), using the expansion:

[TABLE]

In 2010, the limit derived from $WMAP$ 7yr was $f_{\rm NL}^{\rm(local)}\approx 30\pm 20$ at $1\sigma$, a very small non–Gaussianity almost compatible with a zero value allowing for single–scalar fields [1]. $Planck$ 2013 data, expected to give smaller error bars, have narrowed this to $f_{\rm NL}^{\rm(local)}=2.7\pm 5.8$, $f_{\rm NL}^{\rm(equilateral)}=-42\pm 75$, and $f_{\rm NL}^{\rm(orthogonal)}=-25\pm 39$ [73]. Then, $Planck$ 2015 combined temperature and polarization data provided $f_{\rm NL}^{\rm(local)}=0.8\pm 5.0$, $f_{\rm NL}^{\rm(equilateral)}=-4\pm 43$, and $f_{\rm NL}^{\rm(orthogonal)}=-26\pm 21$ ($68\%$ CL, statistical [4]). All these very similar outcomes are consistent with a vanishing $f_{NL}$, which is itself consistent with a very weak primordial non–Gaussianity. However, the interpretation of this characterization of CMB Gaussianity depends on the cosmological model and, in particular, on the type of inflation mechanism that is assumed in this model.

3 Conclusion

Detecting non–Gaussianity in the CMB temperature maps is still a great challenge. Whatever the descriptor, the discrepancy functions of the CMB by $Planck$ 2015 data are not zero but stay within $1\sigma$ or slightly leak into the $2\sigma$ cosmic variance band of current (model–dependent) ensembles. In our search for CMB non–Gaussianities we find systematic signatures of weak non–Gaussianity which would be of importance if the ensemble cosmic variance would be re-evaluated at smaller amplitudes. Compared to the $Planck$ maps we find that the $\Lambda$CDM simulations commonly used offer a systematically larger value for the variances $\sigma_{0}^{2}$ and $\sigma_{1}^{2}$ as observed in [66, 67, 68, 69, 70, 71, 11] and [56, 72]. This anomaly of the variances and its interpretation has to be explored in more detail in the prescriptions of model simulations. Furthermore, after the first main mask is applied, further different possible superimposed foreground and source masks have a big impact upon the magnitude of non–Gaussianities, showing that numerous tiny field sources contribute to residual foreground contamination and may imply a noticeable change in the values of the variances [74]. Several possible strategies should be explored to select ensembles on the basis of the actually observed statistical parameters within a constrained random field approach, in order to reach a deeper understanding of a statistical comparison with a single realization, or to accordingly constrain cosmological models.

Acknowledgements:

This work was conducted within the “Lyon Institute of Origins” under grant ANR–10–LABX–66. A part of this project was funded by the National Science Centre, Poland, under grant 2014/13/B/ST9/00845. The authors wish to thank David Wiltshire and François Bouchet for the invitation to this Focus Article, as well as Giuseppe Fanizza, Takahiko Matsubara, David Wiltshire, Wen Zhao, and an anonymous referee for valuable comments on the manuscript. We also wish to thank Pierre Mourier and Boud Roukema for interesting discussions, and we are grateful to Léo Michel–Dansac for his help concerning various computing facilities, and to Sven Lustig for a fruitful collaboration at an early stage of this work.

This work is based on observations obtained with $Planck$ (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. Some of the results in this paper have been derived using the HEALPix package [75], available at http://healpix.sourceforge.net. The CMB power spectra are calculated from the cosmological parameters using $CAMB$ software written by Lewis and Challinor (http://lambda.gsfc.nasa.gov/toolbox/tb camb form.cfm) from the original Boltzmann codes by Bertschinger, Ma and Bode resumed by Seljak and Zaldarriaga. The $CAMB$ ReadMe (2016) is available here (http://camb.info/readme.html) and $CAMB$ Notes by A Lewis (2014) here (http://cosmologist.info/notes/CAMB.pdf).

We acknowledge the use of the gfortran compiler and the gnuplot graphics utility software under various Linux operating systems.

Appendix A Definitions and notations for the CMB analysis and CMB map generation

Statistical and morphological descriptors of a random field are a priori based on model–independent definitions for a given support manifold. But these descriptors are bias– and error–dependent according to the impact they have on the amplitude and shape of the non–Gaussianity, problems we shall touch upon in this appendix. Numerical biases and errors originate from the pixelation and finite–temperature resolution of the CMB data maps and sample maps, but also from the way infinitesimal calculations are translated into discrete algorithms. These numerical issues affect the descriptors, the fundamental quantities of the field and also the Gaussian predictions.

In what follows we recall some fundamental quantities of the CMB random field as they are used in the numerical analysis. The temperature anisotropy $\delta T$ we use in this paper is dimensionful, $\delta T:=T-T_{0}$ in units of $\mu K$ with $T$ the absolute local temperature, and $T_{0}=2.7255\pm 0.0006$K, is the mean CMB temperature, which can be considered as the monopole component of the CMB maps. The dimensionless CMB temperature anisotropy $(T-T_{0})/T_{0}$ is used in linear perturbation theory, in the Sachs–Wolfe formula (see for instance [76, 77]). The temperature anisotropy $\delta T$ derives from the measure of the spectral radiance $I(\nu({\rm GHz}))$ by the instruments $HFI$ and $LFI$ of the $Planck$ probe [16], and depend on the modelling of several effects: relativistic Doppler–Fizeau, and Sunyaev–Zel’dovich. The magnitude of these effects is a function of the cosmological model, therefore, the CMB spectral radiance and consequently the CMB temperature rest–frame estimations are model–dependent. The suffixes, px or Cℓ apply to the quantities calculated: from the moments of the probability distribution function (no suffix), $P(\tau)$, directly from the pixels, or from the angular power spectrum, respectively.

A.1 Basic quantities of the CMB random field

We consider a Gaussian window function with smoothing scale $\theta_{\rm fwhm}(^{\circ})$,

[TABLE]

For $\sigma_{\mathrm{G}}^{2}\ll 1$ the Gaussian kernel is approximated by

[TABLE]

and once the $a_{\ell,m}$ coefficients of the expansion of $\delta T$ on ${\cal S}^{2}$ in spherical harmonics are calculated from the map, given the pixelation correction ${\rm pxc}(\ell)$, the corrected multipole moments are given by

[TABLE]

[TABLE]

and the angular power spectrum by

[TABLE]

We consider the mean values,

[TABLE]

and the variances,

[TABLE]

For a homogeneous and isotropic Gaussian random field, $\sigma_{\rm 0C\ell}^{2}$ is independent of the absolute direction of light provenance; the variance ensemble average is

[TABLE]

The variance of the local gradient is (compare the definition in equation (77))

[TABLE]

and, only in a homogeneous and isotropic model, for a Gaussian random field, the ensemble average of $\sigma_{1{\rm C\ell}}^{2}$ is

[TABLE]

The measures of the variance and variance of the local gradient, equations (101) and (103), contain an averaged statistical, respectively, morphological information of the random field. Since the model–dependent expressions, equations (102) and (104), are valid only for homogeneous–isotropic Gaussian fields, any significant difference between these two methods indicates a preliminary detection of non–Gaussianity and the degree of conformity with the model for the CMB. For the $\Lambda$CDM sample we notice such a difference between the variances calculated from the pixels and the variances calculated from the angular power spectrum. Model–dependent predictions for $\sigma_{0C\ell}$ and $\sigma_{1C\ell}$ and further discussions can be found in [28].

Table 5 shows the averaged values of $\mu$, $\sigma_{0}$ and $\sigma_{1}$ over the CMB map sample in the $\Lambda$CDM model. Two cases are shown in this table and all the further tables for our $\Lambda$CDM $10^{5}$ map sample at $N_{\rm side}=128$, $\ell_{\rm max}=256$, $2^{\circ}$fwhm: without mask at bin width $13\mu$K and over the full temperature range of the sample; with $U73$ mask subtraction at bin width $6\mu$K and over the limited temperature range fixed by the $Planck$ $SMICA$ map. The effect of passing from no mask to the mask $U73$ is upon $\mu$, $\sigma_{0}$ and $\sigma_{1}$. And the variation of the ratio $\sigma_{1}$/$\sigma_{0}$ is noticeable (e.g. $34.16218/33.93190=1.006787$) and affects the magnitude of the Gaussian premise of the Minkowski Functionals ${\mathrm{v}}_{1}$ and ${\mathrm{v}}_{2}$.

A.2 The cosmic microwave background in the $\Lambda$CDM model.

In order to construct the model map sample, the CMB angular power spectrum is generated with $CAMB$ from the cosmological parameters of the $\Lambda$CDM model according to $Planck$ 2015 [78], p.31, table 4, last column, and Review of particle physics [79]:

$H_{0}=67.74\pm 0.46$ km$s^{-1}$Mpc*-1*, the Hubble constant today

$h=H_{0}$/($100$kms*-1Mpc-1*) $=0.6774$, the normalized Hubble constant today

$T_{0}=2.7255\pm 0.0006$ K, present–day CMB temperature

$\Omega_{b}h^{2}=0.02230\pm 0.00014$, baryon density

$\Omega_{c}h^{2}=0.1188\pm 0.0010$, Cold Dark Matter density

$\Omega_{\nu}h^{2}=0.00209$, neutrino density

$\Omega_{k}=0$, constant–curvature density parameter

$Y_{P}=0.249_{-0.026}^{+0.025}$, helium fraction

$N_{\rm eff}=3.04\pm 0.33$, number of massless neutrinos

$\sum m_{\nu}<$ 0.194 eV, neutrino mass eigenstates

$\tau$=0.066 $\pm$ 0.012, reionization optical depth

$z_{re}$=8.8 ${}_{-1.1}^{+1.2}$, redshift of the reionization

$n_{S}$=0.9667 $\pm$ 0.0040, scalar spectral index

$x_{e}$=1, ionization fraction

$z_{*}$=1089.90 $\pm$ 0.23, redshift of decoupling

$\Omega_{\Lambda}$=0.6911 $\pm$ 0.0062, cosmological constant density

$\Omega_{m}$=0.3089 $\pm$ 0.0062, total matter density

$Age$= 13.799 $\pm$ 0.021 Gyr, age of Universe

A.3 The CMB radiation components

The full Boltzmann physics is implemented in the $CAMB$ software; thus, the maps of the CMB sample in the $\Lambda$CDM model include the following effects: ordinary Sachs–Wolfe effect / Doppler effect / Silk damping / Reionization / Polarization of photons / Neutrinos / Integrated Sachs–Wolfe effect ($ISW$) / Lensing. The weak lensing effect upon the CMB is treated in [8] and in [80] beyond the Born approximation, which is usually applied at first–order of perturbation theory. For a second–order treatment the deflection angles, still assumed small if $\ell\leqslant 2500$666Two times the typical angular deviation equals $4$ arc minutes (CMB lensing due to structures at $z<20$ in a flat $FLRW$ model); this corresponds to $\ell_{max}=180^{\circ}/(4/60)=2700$., are no more Gaussianly distributed, but the post–Born corrections would be weakly detectable in low noise CMB temperature maps merging the contributions of hundreds of $\ell$–values between $\ell=1000$ and $\ell=2500$. Our present study uses an $\ell$–range [2,256], well below $\ell=1000$. One notices that an evaluation of the Sunyaev–Zeldovich effect ($SZ$) on the CMB radiation is not implemented in $CAMB$, the sources responsible of the $SZ$ effect being more or less excluded when using an appropriate foreground mask. It is clear that, in this list, some of the effects are not primordial but depend on model–dependent knowledge of the matter distribution and the physical assumptions in the concordance model, i.e., mainly the $ISW$ effect, the lensing and the neutrinos.

A.4 Map ensemble and statistical stability

From the power spectrum we generate with $synfast$ a $10^{5}$ map ensemble. Given the map resolution ($N_{\rm side}=128$, $\ell_{\rm max}=256$, and $2^{\circ}$fwhm), the number of $10^{5}$ maps allows a statistical stability in the sense of the second averaged MF ${\mathrm{v}}_{1}$, which becomes smooth and stable around $10^{5}$ maps, a second different sample of $10^{5}$ maps gives very similar results and doubling the sample brings no visible improvement in the shape of ${\mathrm{v}}_{1}$. A number of $10^{5}$ maps is satisfactory in the sense of the statistical stability for the MFs themselves, but less for the PDF. However, once we come to the discrepancy functions, the stability is definitely impaired when using a smaller number of maps, and $10^{5}$ maps must be considered as a minimum requirement. We shall not develop more on these studies in the present work.

A.5 Discretization: definition of the $\tau$–lattice

One observes that in the case of $10^{5}$ realizations there is an interval in $\tau$, $[\tau_{-},\tau_{+}]$ such that $P(\tau)=0$ for all realizations, if $\tau\notin[\tau_{-},\tau_{+}]$. In order to simplify the problem, we use a symmetric interval $[-\tau_{m},+\tau_{m}]$, where $\tau_{m}={\rm max}\{|\tau_{-}|,\tau_{+}\}$. The interval $[-\tau_{m},+\tau_{m}]$ is divided into $2L+1$ bins of equal bin width $\Delta\tau:=2\tau_{m}/(2L+1)$, where the mid–points of the bins are given by the “$\tau-$lattice”, $\tau_{l}=-\tau_{m}+(2l-1)(\Delta\tau/2)$, $l=1,2,\cdots,(2L+1)$, in such a way that the $(L+1)^{\rm th}$ bin is centered at $\tau_{L+1}=0$. (Example: $\tau_{m}=396.5\;\mu\mathrm{K}$, $L=30$, $\Delta\tau=13\mu K$ $\Rightarrow(2L+1)=61$ bins.) Then, the discretized PDF for a given realization, which is now a step–function, i.e. piecewise continuous, can be represented as a histogram that is defined by ($l=1,2,\cdots,(2L+1)$):

[TABLE]

with $N_{\rm tot}$ the total number of pixels, and for $\tau$ within the $l^{\rm th}$ bin, i.e., $\tau\in\left[\tau_{l}-\frac{\Delta\tau}{2},\tau_{l}+\frac{\Delta\tau}{2}\right)$ (half open interval!).

The ensemble average (mean) $\langle P(\tau)\rangle$ is given by the arithmetic mean of all the PDFs of the $10^{5}$ realizations. $\langle P(\tau)\rangle$ is still constant in a given bin, denoted by $\langle P(\tau_{l})\rangle$ in the $l^{\rm th}$ bin and, thus, the derivative ${\mathrm{d}}\langle P(\tau)\rangle/{\mathrm{d}}\tau$ is zero almost everywhere, but $\langle P(\tau)\rangle$ will have, in general, $(2L+2)$ jumps, namely at the points $\tau_{l}+\Delta\tau/2$, $l=1,2,\cdots 2L$ and at the points $\mp\tau_{m}$ (iff $\langle P(\tau)\rangle$ possesses this special property!) In general, there will be a first jump at a value $\tau_{0}\in[-\tau_{m},\tau_{1}+\Delta\tau/2)$, and a last jump at a value $\tau_{2L+2}\in(\tau_{2L+1}-\Delta\tau/2,\tau_{m}]$, and these contributions have to be treated separately.

If we ignore the last subtlety, the expectation value (3) of a given random field $f(\delta T)$ is exactly given by the finite sum

[TABLE]

As an important example, this yields with $f(\tau)=\tau^{n}$ the exact formula for the moments $\alpha_{n}$ of $\delta T$ (see equation (4) for $n=0,1,2,\cdots$):

[TABLE]

For $n=0$ and $n=1$ one obtains the expected values

[TABLE]

[TABLE]

whereas, for $n\geq 2$, one gets a decomposition into a “main term” and an “exact correction” in the form of a finite series in the bin width $\Delta\tau$:

[TABLE]

(Example for $n=2$: $\sigma_{0}^{2}:=\alpha_{2}=\sum_{l=1}^{2L+1}\tau_{l}^{2}\langle P(\tau_{l})\rangle\Delta\tau+(\Delta\tau)^{2}/12$.)

In many papers on MFs the integral in equation (3) and similar integrals are approximated by the “main term” and, thus, the results suffer from an error, which is exactly given (in the case of (3)) by the “correction” in (110), if it is not taken into account. For a discussion of the correction term in the case of the MFs ${\mathrm{v}}_{1}$ and ${\mathrm{v}}_{2}$, see [81]. This error can be made small, iff $\Delta\tau$ or $\Delta\nu$ is chosen small enough in principle, which is not the case, i.e., for $\Delta\tau=13\mu K$. Looking at the “main term” (mentioned above) without correction,

[TABLE]

we obtain that $\sigma_{0}^{\prime}=\sqrt{\alpha^{\prime}_{2}}=59.6516\mu\mathrm{K}$ to be compared with $\sigma_{0\rm px}=59.5335\mu\mathrm{K}$ (second equation in (101) and see also table 5).

Appendix B The generating functions of the moments and cumulants, and a closed expression for the NG–parameters ${\bf a}_{\bf P}{\bf(n)}$

From the definition (20) of the discrepancy function $\Delta_{P}(\tau)$ and the Hermite expansion (21), one obtains with (6) for the generating function $M(x)$ of the moments $\alpha_{n}$ (with $\alpha_{0}=1$, $\alpha_{1}=\mu=0$, $\alpha_{2}=\sigma_{0}^{2}$):

[TABLE]

with

[TABLE]

and

[TABLE]

Using the generating function of the Hermite polynomials [82],

[TABLE]

equation (B) can be rewritten as

[TABLE]

which leads with the orthogonality relation (27) to

[TABLE]

and, thus, with (112), (113), to the factorization

[TABLE]

with

[TABLE]

The cumulant generating function, defined in (7), is then additive,

[TABLE]

as given in equations (22) and (23) of the main text. It is convenient to write $\Delta C(x):=\ln K(x)$ in terms of the dimensionless variable $z:=\sigma_{0}x$,

[TABLE]

where $C_{n}:=\varkappa_{n}/\sigma_{0}^{n}$ are the normalized (dimensionless) cumulants. To obtain closed expressions for the $\rm NG$–coefficients $a_{P}(n)$ in terms of the cumulants $C_{n}$, we recall the generating function of the complete Bell polynomials $B_{n}$ [83] ($B_{0}=1$):

[TABLE]

By expanding the exponential and comparing in (120) the terms of the same power in $z$, it is not difficult to obtain, e.g., the first four Bell polynomials: $B_{1}(x_{1})=x_{1}$, $B_{2}(x_{1},x_{2})=x_{1}^{2}+x_{2}$, $B_{3}(x_{1},x_{2},x_{3})=x_{1}^{3}+3x_{1}x_{2}+x_{3}$, $B_{4}(x_{1},x_{2},x_{3},x_{4})=x_{1}^{4}+6x_{1}^{2}x_{2}+4x_{1}x_{3}+3x_{2}^{2}+x_{4}$. A comparison between (119) and (120) yields the closed expression ($x_{1}=x_{2}=0,x_{3}=\gamma_{1},x_{4}=\gamma_{2},x_{n}=C_{n}$ for $n\geq 5$):

[TABLE]

as given in the main text in equation (26). The explicit expressions for $B_{3}$ and $B_{4}$ give $a_{P}(3)=B_{3}(0,0,\gamma_{1})=\gamma_{1}$, $a_{P}(4)=B_{4}(0,0,\gamma_{1},\gamma_{2})=\gamma_{2}$, in agreement with (2.2). In order to obtain the higher Bell polynomials, one can either use the recurrence relations,

[TABLE]

or the combinatorial expression

[TABLE]

Here, the sum is over all partitions $\pi(n)$ of $n$, i.e., over all positive integers $a_{m}$ such that $\sum_{m=1}^{n}m\;a_{m}=n$. The multinomial coefficients

[TABLE]

are given, for $n=1,2,\cdots,10$, in table 24.2 in [63]. The Bell polynomials have the nice property that their coefficients are integers and, therefore, the NG coefficients $a_{P}(n)$ of the PDF discrepancy function are linear combinations of the normalized cumulants with integer coefficients, as seen in equation (2.2). Since the NG–coefficients $a_{0}(n)$ of the discrepancy function $\Delta_{0}$ are related to the $a_{P}(n)$’s by $a_{0}(n)=\frac{a_{P}(n+1)}{n+1},n\geq 2$, it follows that the coefficients of the $a_{0}(n)$’s are, in general, rational numbers (see equation (2.3)).

Finally, we give an alternative closed formula for the expansion coefficients $a_{P}(n)$, which does not express them in terms of the cumulants $C_{n}$ as in equation (121), but rather in terms of the normalized moments ${\hat{\alpha}}_{n}$ (see equation (4)),

[TABLE]

Replacing $\Delta_{P}$ in equation (28) by its definition (20), we obtain:

[TABLE]

which gives, with the definition (3) and the orthogonality relation (27), ($n=0,1,2,\dots$):

[TABLE]

From (127) follows immediately $a_{P}(0)=0$, $a_{P}(1)={\hat{\alpha}}_{1}=\mu/\sigma_{0}=0$, and $a_{P}(2)=\alpha_{2}/\sigma_{0}^{2}-{\hat{\alpha}}_{0}=0$, and for the required $a_{P}(n)$’s with $n\geq 3$ using the expansion [63],

[TABLE]

the closed expression

[TABLE]

By virtue of this formula, the integral in the original definition (28) is replaced by a finite sum and, thus, one obtains in combination with the formula (4) for the moments simple closed expressions for the Hermite expansion coefficients.

Appendix C General Hermite and Edgeworth expansions

Here, we summarize some well–known mathematical facts (see, e.g., [84, 85, 86, 87, 88, 89]), which are used in section 2 for expanding the various discrepancy functions $\Delta_{P}(x)$, $\Delta_{0}(x)$, etc., in terms of Hermite polynomials.

Let $f(x)$, $\mathbb{R}\rightarrow\mathbb{R}$, be square integrable with respect to a positive weight function $w(x)$, i.e., $f(x)\in\mathbb{H}:=L^{2}(\mathbb{R},w(x){\mathrm{d}}x)$. For $f,g\in\mathbb{H}$, we define the inner product

[TABLE]

which satisfies the Schwarz inequality $|(f,g)|\leq||f||\cdot||g||$, where $||f||$ denotes the norm of $f$, $||f||:=\sqrt{(f,f)}$. Let $\psi_{0}(x),\psi_{1}(x),\cdots\in\mathbb{H}$ be an orthonormal system satisfying the orthogonality relation $(m,n\in\mathbb{N}_{0})$,

[TABLE]

and let $f\in\mathbb{H}$ be any function. Then, the numbers

[TABLE]

are called the expansion coefficients (“Fourier coefficients”) of $f$ with respect to the $\psi_{n}$’s. From the relation

[TABLE]

one obtains $\sum_{n=0}^{N}b_{n}^{2}\leq||f||^{2}$, and since the right–hand–side of the last inequalities is independent of $N$, we obtain Bessel’s inequality

[TABLE]

This proves that the sum of the squares of the expansion coefficients $b_{n}$ always converges.

Consider now, for a given function $f(x)\in\mathbb{H}$ the following linear combination:

[TABLE]

with constant coefficients $\gamma_{n}$ and fixed $N$. Then, there arises the question under which conditions the approximation (135) can be considered as an approximation “in the mean” such that the mean square error

[TABLE]

is as small as possible. From the identity

[TABLE]

it follows immediately that ${\cal E}_{N}$ takes on its least value for $\gamma_{n}=b_{n}\ (n=0,1,\cdots,N)$. If the error ${\cal E}_{N}$ converges for every piecewise continuous function $f\in\mathbb{H}$ to zero as $N$ goes to infinity, then the orthonormal system $\{\psi_{n}\}$ is said to be complete, i.e., it provides a complete basis of $\;\mathbb{H}$, and Bessel’s inequality (134) becomes an equality for every piecewise continuous function $f$,

[TABLE]

which is known as completeness relation (also called Parseval’s equation).

It is important to note that the completeness of the system $\{\psi_{n}\}$, expressed by the equation

[TABLE]

does not necessarily imply that $f(x)$ can be expanded in a series in the functions $\psi_{n}(x)$. The expansion

[TABLE]

is, however, valid if the series in (140) converges uniformly and, thus, the limit in (139) can be carried out under the integral.

In this paper, our main concern is not the theoretical problem to find an expansion (140) of $f(x)$, but rather the practical problem to obtain a representation (135) of $f(x)$ with $\gamma_{n}=b_{n}$ with a small number of terms, $N<10$, say, which provides a fairly good approximation by minimizing the error (136). To this end, it is convenient to choose an orthogonal basis $\{\phi_{n}(x)\}$ in $\mathbb{H}$, where each $\phi_{n}(x)$ is a polynomial of degree $n$. For the Gaussian weight function $w^{\rm G}(x):=\exp{(-x^{2}/2)}=\sqrt{2\pi}P^{\rm G}(x)$, it turns out that the polynomials are uniquely determined (up to a multiplicative constant in each polynomial) by the Hermite polynomials ${\rm He}_{n}(x)$ (c.f. [63]),

[TABLE]

which provide a complete orthogonal (not orthonormal) system in $\mathbb{H}=L^{2}(\mathbb{R},\exp{(-x^{2}/2)}\ {\mathrm{d}}x)$, satisfying the orthogonality relation (see (27))

[TABLE]

with $h_{n}:=\sqrt{2\pi}n!\,$. In order that the piecewise continuous function $f(x)$ satisfies $||f||<\infty$ in $\mathbb{H}$, it must obey the asymptotic condition

[TABLE]

The approximation (135) becomes then the polynomial approximation

[TABLE]

where the expansion coefficients $a(n)=\frac{1}{\sqrt{2\pi}}(f,\phi_{n})$ are explicitly given by

[TABLE]

The completeness relation (138) then reads:

[TABLE]

which implies the asymptotic behaviour

[TABLE]

It is important to bear in mind that, even in the case when the series (140) is uniformly convergent, it by no means follows that the $N^{th}$ partial sum (144) is the best selection of $N$ terms for representing the function $f(x)$. Even though (144) gives the best fit in the sense of least squares by minimizing the error (136), it may be that some other measure of approximation is better suited for a given problem [85]. All the more this may be the case if the series (140) is divergent, in which case one may ask whether there exists an asymptotic expansion in the sense of Poincaré (c.f. [86]).

In the particular case of the Hermite expansion (144), many authors have worked on asymptotic expansions since quite a long time, mainly in the context of probability theory, statistics, number theory, and mathematical aspects of insurance risk (c.f. [85, 87, 88, 89], and references therein). The results relevant to us in this paper are connected with the attempts to give a refinement of the classical central limit theorem in probability theory, and are often, historically incorrect, referred to as Charlier or Gram–Charlier A–series and Edgeworth expansion, although they had been introduced by Tchebychev already before [89]. Since Matsubara’s expansion of the MFs in [57, 58] (see also [60] and [61]), based on the assumption of hierarchical ordering (HO), is formally closely related to the Edgeworth expansion, we shall summarize the main properties of the latter. It turns out that the Edgeworth expansion furnishes, in the case of the central limit theorem, a genuine asymptotic expansion with a well–defined remainder term.

Let $x_{1},x_{2},\cdots,x_{n}$ be independent and identically distributed random variables with a common continuous distribution function, such that every $x_{k}$ has zero mean, standard deviation $\sigma$, and a third absolute moment $\beta_{3}=\langle|\delta T|^{3}\rangle$. Consider the standardized sum variable

[TABLE]

and let $P_{n}(x)$ be the PDF of $y_{n}$. The central limit theorem then asserts that, as $n\rightarrow\infty$, under appropriate conditions, $P_{n}(x)$ tends to the Gaussian (normal) PDF, $P^{\rm G}(x)=\frac{1}{\sqrt{2\pi}}\exp{(-x^{2}/2})$. The Edgeworth expansion is a refinement of this by giving also the rate of convergence to the Gaussian limit. Define the discrepancy function

[TABLE]

Then, the Edgeworth (E) expansion [85, 87, 88, 89] is the following asymptotic expansion:

[TABLE]

where $q_{k}(x)$ is a polynomial of degree $3k$, which only depends on the normalized cumulants $C_{n}\;(n\geq 3;C_{3}=\gamma_{1},C_{4}=\gamma_{2})$, and which can be expressed as a linear combination of the Hermite polynomials ${\rm He}_{n}(x)$. For $k=1,\cdots,4$ they are explicitly given by:

[TABLE]

The expansion (150), $\delta_{n}^{E}(x)$, taking the first $K$ terms into account, is an asymptotic expansion of $\delta_{n}(x)$ in powers of $1/\sqrt{n}$ with a remainder of the same order as the first term neglected (i.e., it is an “asymptotic expansion to $K$ terms” as defined by Poincaré, c.f. [86]). Thus, it gives a correction to the central limit theorem, with a well–defined error of order $n^{-(K+1)/2}$, in the case where $n$ is finite but large $(1\ll n<\infty)$. Due to the Hermite expansions (C) of the polynomials $q_{k}(x)$, the Edgeworth expansion (150) can be interpreted as a particular rearrangement of the following polynomial approximation:

[TABLE]

where the ordering is not with respect to powers of $1/\sqrt{n}$, but rather according to the order of the Hermite polynomials. (The $n-$dependence of the coefficients $a^{E}(m)$ is not explicitly noted.) By a comparison of (152) with (150,C) one immediately reads off the first $12$ expansion coefficients (here we choose $K=4$):

[TABLE]

(Note that the coefficients $a^{E}(m)$, $m\geq 7$, change if one goes to higher order in $1/\sqrt{n}$, e.g., for $K=5$ they receive additional contributions of order $n^{-5/2}$ and $n^{-3}$.) If the coefficients (C) are compared with the coefficients $a_{P}(m)$ in equation (2.2), derived for a general Hermite expansion for the non–Gaussianities as in (21), it is not difficult to derive the law that governs the size, with respect to powers of $1/\sqrt{n}$, of the coefficients $a^{E}(m)$ associated with the asymptotics (150) of the central limit theorem. It simply says that every cumulant $C_{r}$ in (2.2) has to be replaced by $C_{r}/n^{(r-2)/2}$ (for fixed $n\gg 1,r\geq 3$).

This is precisely the law, rigorously proven for the central limit theorem, which in some models of inflation is built in as hierarchical ordering (HO). But then, the role of the small dimensionless expansion parameter $1/\sqrt{n}$ is played by the standard deviation $\sigma_{0}$ which, for the CMB anisotropy, has a fixed value and cannot be made arbitrarily small as in the case of the central limit theorem. It is, thus, important to perform numerical checks, as done in this paper. Furthermore, one should keep in mind that the Edgeworth expansion, applied to a given PDF, can suffer from the fact that the PDF is not correctly normalized and it can exhibit undesirable properties such as negative probabilities. In particular, the approximation deteriorates in the tails.