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

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

Figure 1 shows the 10⁵ histograms envelope of the individual PDFs in the $\Lambda$CDM sample over a total temperature range  ±396.5 μ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 appendix A.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&F(\tau ):=\text{prob}(\delta T<\tau )={\int}_{-\infty}^{\tau}P\left({{\tau}^{\prime}}\right)\text{d}{{\tau}^{\prime}}\;,\end{eqnarray} \tag{ 1 }$$

respectively, the complementary cumulative distribution function ${{F}_{C}}(\tau )$ (CCDF)⁷,

$$\begin{eqnarray}&&{{F}_{C}}(\tau ):=\text{prob}(\delta T\geqslant \tau )={\int}_{\tau}^{\infty}P\left({{\tau}^{\prime}}\right)\text{d}{{\tau}^{\prime}}=1-F(\tau )\;,\end{eqnarray} \tag{ 2 }$$

satisfying $F(-\infty )=0$, $F(\infty )=1$, respectively, ${{F}_{C}}(-\infty )=1$, ${{F}_{C}}(\infty )=0$. Note that $F(\tau )$ is non-decreasing and $\text{prob}\left({{\tau}_{1}}\leqslant \tau <{{\tau}_{2}}\right)=F\left({{\tau}_{2}}\right)-F\left({{\tau}_{1}}\right)$.

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

$$\begin{eqnarray}&&\langle \,f(\delta T)\rangle :={\int}_{-\infty}^{\infty}f(\tau )P(\tau )\text{d}\tau \;,\end{eqnarray} \tag{ 3 }$$

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

$$\begin{eqnarray}&&{{\alpha}_{n}}:=\langle {{(\delta T)}^{n}}\rangle :={\int}_{-\infty}^{\infty}{{\tau}^{n}}P(\tau )\text{d}\tau \quad ;\quad n=0,1,2,\cdots \;,\end{eqnarray} \tag{ 4 }$$

the central moments m_n (i.e. the moments about the mean $\mu :={{\alpha}_{1}}$),

$$\begin{eqnarray}&&{{m}_{n}}:=\langle {{(\delta T-\mu )}^{n}}\rangle :={\int}_{-\infty}^{\infty}{{(\tau -\mu )}^{n}}P(\tau )\text{d}\tau \;,\end{eqnarray} \tag{ 5 }$$

and the cumulants ${{\varkappa}_{n}}:={{\langle {{(\delta T)}^{n}}\rangle}_{C}}$. The generating function of the moments ${{\alpha}_{n}}$ is given by:

$$\begin{eqnarray}&&M(x):=\langle {{\text{e}}^{x\delta T}}\rangle ={\int}_{-\infty}^{\infty}{{\text{e}}^{x\tau}}P(\tau )\text{d}\tau =\underset{n=0}{\overset{\infty}{\sum}}\,{{\alpha}_{n}}\frac{{{x}^{n}}}{n!}\;,\end{eqnarray} \tag{ 6 }$$

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

$$\begin{eqnarray}&&C(x):=\ln M(x)=\underset{n=1}{\overset{\infty}{\sum}}\,{{\varkappa}_{n}}\frac{{{x}^{n}}}{n!}\;.\end{eqnarray} \tag{ 7 }$$

The first few moments and cumulants are ${{\alpha}_{0}}={{m}_{0}}=1$, the mean values $\mu :={{\alpha}_{1}}={{\varkappa}_{1}}$, respectively, m₁  =  0, and the variance $\sigma _{0}^{2}$,

$$\begin{eqnarray}&&\sigma _{0}^{2}:=\langle {{(\delta T-\mu )}^{2}}\rangle ={{m}_{2}}={{\alpha}_{2}}-{{\mu}^{2}}={{\varkappa}_{2}}\;,\end{eqnarray} \tag{ 8 }$$

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

$$\begin{eqnarray}&&{{m}_{n}}=\underset{r=0}{\overset{n}{\sum}}\,\left(\begin{array}{c} n \\ r \end{array}\right)\,{{(-1)}^{n-r}}{{\mu}^{n-r}}{{\alpha}_{r}}\quad ;\quad n=0,1,2,\cdots \;,\end{eqnarray} \tag{ 9 }$$

and, for $n\geqslant 2$:

$$\begin{eqnarray}&&{{\varkappa}_{n}}={{m}_{n}}-\underset{r=1}{\overset{n-1}{\sum}}\,\left(\begin{array}{c} n-1 \\ r-1 \end{array}\right)\,{{\varkappa}_{r}}{{m}_{n-r}}\;.\end{eqnarray} \tag{ 10 }$$

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

$$\begin{eqnarray}&&{{m}_{3}}={{\varkappa}_{3}}={{\alpha}_{3}}-3\mu {{\alpha}_{2}}+2{{\mu}^{3}}\;,\end{eqnarray} \tag{ 11 }$$

respectively, the dimensionless skewness coefficient

$$\begin{eqnarray}&&{{\gamma}_{1}}:=\frac{{{m}_{3}}}{\sigma _{0}^{3}}\;.\end{eqnarray} \tag{ 12 }$$

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),

$$\begin{eqnarray}&&{{\gamma}_{2}}:=\frac{{{m}_{4}}}{\sigma _{0}^{4}}-3\;.\end{eqnarray} \tag{ 13 }$$

It should be stressed that the basic equations (1)–(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 )=\text{d}F/\text{d}\tau =-\text{d}{{F}_{C}}/\text{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:

$$\begin{eqnarray}&&\langle \,f(\delta T)\rangle :={\int}_{-\infty}^{\infty}f(\tau )\text{d}F(\tau )\;.\end{eqnarray} \tag{ 14 }$$

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):

$$\begin{eqnarray}&&\langle \,f(\delta T)\rangle :={\int}_{-\infty}^{\infty}f(\tau ){{F}^{\prime}}(\tau )\text{d}\tau +\underset{\ell}{\sum}\,f\left({{\tau}_{\ell}}\right)\left[F\left({{\tau}_{\ell}}+0\right)-F\left({{\tau}_{\ell}}-0\right)\right]\;,\end{eqnarray} \tag{ 15 }$$

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,

$$\begin{eqnarray}&&{{P}^{\text{G}}}(\tau )=\frac{1}{\sqrt{2\pi}{{\sigma}_{0}}}{{\text{e}}^{-\frac{{{(\tau -\mu )}^{2}}}{2\sigma _{0}^{2}}}}\;,\end{eqnarray} \tag{ 16 }$$

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

$$\begin{eqnarray}&&{{M}^{\text{G}}}(x)={{\text{e}}^{\;\mu x+\frac{\sigma _{0}^{2}}{2}\,{{x}^{2}}}}\;,\end{eqnarray} \tag{ 17 }$$

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

$$\begin{eqnarray}&&{{C}^{\text{G}}}(x)=\mu x+\frac{\sigma _{0}^{2}}{2}\,{{x}^{2}}\;.\end{eqnarray} \tag{ 18 }$$

Our numerical results for the CMB anisotropy are based on an ensemble of 10⁵ realizations (for details, see appendix A). It turns out that the mean value $\mu =\langle \delta T\rangle$ (calculated from pixels) is negligibly small, $\mu =\mathcal{O}\left({{10}^{-7}}~\mu \right.$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}^{\text{G}}}(\tau )$ vanish, $\alpha _{2n+1}^{\text{G}}=0$, and that the even moments are given by:

$$\begin{eqnarray}&&\alpha _{2n}^{\text{G}}=\frac{{{2}^{n}}}{\sqrt{\pi}}\; \Gamma \left(n+\frac{1}{2}\right)\sigma _{0}^{2n}=\frac{(2n)!}{{{2}^{n}}n!}\sigma _{0}^{2n}\;,\end{eqnarray} \tag{ 19 }$$

and increase with increasing n. (Note that $m_{n}^{\text{G}}=\alpha _{n}^{\text{G}}$.) Equation (19) gives $\alpha _{4}^{\text{G}}=3\sigma _{0}^{4}$ and, thus, one obtains from (12) and (13) $\gamma _{1}^{\text{G}}=\gamma _{2}^{\text{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 _{2}^{\text{G}}=\sigma _{0}^{2}$ and, therefore, any non-vanishing cumulant with $n\geqslant 3$ is a measure of non-Gaussianity.

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:

$$\begin{eqnarray}&&{{ \Delta }_{P}}(\tau ):=\frac{P(\tau )-{{P}^{\text{G}}}(\tau )}{{{P}^{\text{G}}}(0)}\;,\end{eqnarray} \tag{ 20 }$$

with ${{P}^{\text{G}}}(0)={{\left(\sqrt{2\pi}{{\sigma}_{0}}\right)}^{-1}}=\text{max}\left\{{{P}^{\text{G}}}(\tau )\right\}$. Here, $P(\tau )$ denotes the ensemble average over a large number of realizations (10⁵ 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}^{\text{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}^{\text{G}}}(\tau )$ are very small. Under very general conditions (for details see appendix C), ${{ \Delta }_{P}}(\tau )$ can be written as a product of a Gaussian and a function $h(\nu )\in {{L}^{2}}\left(\mathbb{R},w(\nu )\text{d}\nu \right)$, where $w(\nu )$ denotes the weight function, $w(\nu )=\exp \left(-{{\nu}^{2}}/2\right)$, expressed in terms of the dimensionless scaled temperature variable $\nu :=\tau /{{\sigma}_{0}}\in \mathbb{R}$.

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

$$\begin{eqnarray}&&{{ \Delta }_{P}}(\tau )={{\text{e}}^{-{{\tau}^{2}}/2\sigma _{0}^{2}}}\;\underset{n=3}{\overset{\infty}{\sum}}\,\frac{{{a}_{P}}(n)}{n!}\text{H}{{\text{e}}_{n}}\left(\frac{\tau}{{{\sigma}_{0}}}\right)\;,\end{eqnarray} \tag{ 21 }$$

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 $\text{H}{{\text{e}}_{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 \left(-{{\nu}^{2}}\right)$. In the cosmology literature the $\text{H}{{\text{e}}_{n}}(\nu )$ are used, but unfortunately denoted as ${{H}_{n}}(\nu )$! The relation between the two is $\text{H}{{\text{e}}_{n}}(\nu )={{2}^{-n/2}}{{H}_{n}}\left(\nu /\sqrt{2}\right.$), (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\left(\exp \left({{\nu}^{2}}/4\right)/\sqrt{|\nu |}\right)$ for $|\nu |\to \infty$ (see appendix 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}^{\text{G}}}(\tau )+\left[1/\left(\sqrt{2\pi}{{\sigma}_{0}}\right)\right]{{ \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 appendix B):

$$\begin{eqnarray}&&C(x)=\frac{\sigma _{0}^{2}}{2}~{{x}^{2}}+ \Delta C(x)\;,\end{eqnarray} \tag{ 22 }$$

with

$$\begin{eqnarray} &&\Delta C(x):=\ln \left[1+\underset{n=3}{\overset{\infty}{\sum}}\,\frac{{{a}_{P}}(n)\sigma _{0}^{n}}{n!}\,{{x}^{n}}\right]\,=\,\underset{n=3}{\overset{\infty}{\sum}}\,\frac{{{\varkappa}_{n}}}{n!}\,{{x}^{n}}\;.\end{eqnarray} \tag{ 23 }$$

Then, the cumulants are ${{\varkappa}_{1}}=0$, ${{\varkappa}_{2}}=\sigma _{0}^{2}$, and the higher cumulants, for $n\geqslant 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,

$$\begin{eqnarray}&&{{C}_{n}}:=\frac{{{\varkappa}_{n}}}{\sigma _{0}^{n}}=\frac{{{\langle {{(\delta T)}^{n}}\rangle}_{C}}}{\sigma _{0}^{n}}\;\;:\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} &{{a}_{P}}(3)={{\gamma}_{1}}\quad ;\quad {{a}_{P}}(4)={{\gamma}_{2}}\quad ;\quad {{a}_{P}}(5)={{C}_{5}}\quad ;\quad {{a}_{P}}(6)=10\gamma _{1}^{2}+{{C}_{6}}\quad ; \\ {} & {{a}_{P}}(7)=35{{\gamma}_{1}}{{\gamma}_{2}}+{{C}_{7}}\quad ;\quad {{a}_{P}}(8)=56{{\gamma}_{1}}{{C}_{5}}+35\gamma _{2}^{2}+{{C}_{8}}\quad ; \\ {} & {{a}_{P}}(9)=280\gamma _{1}^{3}+126{{\gamma}_{2}}{{C}_{5}}+84{{\gamma}_{1}}{{C}_{6}}+{{C}_{9}}\quad ; \\ {} & {{a}_{P}}(10)=2100\gamma _{1}^{2}{{\gamma}_{2}}+120{{\gamma}_{1}}{{C}_{7}}+210{{\gamma}_{2}}{{C}_{6}}+126C_{5}^{2}+{{C}_{10}}\quad ; \\ {} & {{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}}+462{{C}_{5}}{{C}_{6}}+{{C}_{11}}\quad ; \\ {} & {{a}_{P}}(12)=15\,400\gamma _{1}^{4}+5775\gamma _{2}^{3}+27\,720{{\gamma}_{1}}{{\gamma}_{2}}{{C}_{5}}+9240\gamma _{1}^{2}{{C}_{6}}+220{{\gamma}_{1}}{{C}_{9}}+495{{\gamma}_{2}}{{C}_{8}} \\ &\qquad +792{{C}_{5}}{{C}_{7}}+462C_{6}^{2}+{{C}_{12}}\;. \end{array}\end{eqnarray} \tag{ 25 }$$

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

$$\begin{eqnarray}&&{{a}_{P}}(n)={{B}_{n}}\left(0,0,{{\gamma}_{1}},{{\gamma}_{2}},{{C}_{5}},\cdots,{{C}_{n}}\right)\;,\end{eqnarray} \tag{ 26 }$$

where ${{B}_{n}}\left({{x}_{1}},{{x}_{2}},\cdots,{{x}_{n}}\right)$ denotes the nth complete Bell polynomial (see appendix B). An equivalent closed expression in terms of the moments is given in equation (B.18).

Table 1 shows the moments ${{\alpha}_{n}}$ (main term, equation (A.17) in appendix 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 \text{K},~3558.315\,19~\mu {{\text{K}}^{2}}$; ${{\varkappa}_{0,1,2}}=1,~0~\mu \text{K},~3558.315\,19~\mu {{\text{K}}^{2}}$ and C_0,1,2  =  1, 0, 1, and we obtain ${{\sigma}_{0}}=59.651\,66~\mu \text{K}$, ${{\gamma}_{1}}=-5.1872\times {{10}^{-4}}$ and ${{\gamma}_{2}}=5.825\times {{10}^{-5}}$.

Table 1. Table of the moments ${{\alpha}_{n}}$ (main term, equation (A.17)), cumulants ${{\varkappa}_{n}}$, and dimensionless normalized cumulants C_n of the probability density function $P(\tau )$.

(Units of $\mu {{K}^{n}}$) $\Lambda$CDM sample Full individual map range, no mask, 2°fwhm, bin 13μK see appendix A
n	3	4	5	6	7	8	9	10	11	12
${{\alpha}_{n}}$	$-110.103\,727$	$37\,985\,558.5$	$-4040\,069.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.103\,727$	$737.480\,513$	$-122\,231.715$	$-108\,751\,936.$	$785\,344\,215.$	$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}}$
Cn	$-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.142\,91\times {{10}^{-2}}$	$-1.960\,36\times {{10}^{-2}}$	$7.149\,27\times {{10}^{-2}}$	$0.414\,7370$

(Units of $\mu {{K}^{n}}$) $\Lambda$CDM sample Equal temperature range (ETR  ±  201 μK), U73 mask, 2°fwhm, bin 6μK
n	0	1	2	3	4	5	6	7	8	9
${{\alpha}_{n}}$	$1.000\,1469$	$-1.4912\times {{10}^{-3}}$	$3515.737\,11$	$-113.475\,560$	$36\,119\,872.5$	$-2963\,938.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.737\,10$	$-97.748\,0482$	$-961\,350.317$	$741\,927.368$	$-1.3\times {{10}^{10}}$	$-3.4\times {{10}^{9}}$	$1.8\times {{10}^{14}}$	$-1.6\times {{10}^{14}}$
Cn	1.0	$-2.51\times {{10}^{-5}}$	1.0	$-4.689\times {{10}^{-4}}$	$-7.777\,66\times {{10}^{-2}}$	$1.0123\times {{10}^{-3}}$	$-0.297\,3678$	$-1.3108\times {{10}^{-3}}$	$1.169\,9687$	$-1.761\,38\times {{10}^{-2}}$

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

$$\begin{eqnarray}&&{\int}_{-\infty}^{\infty}{{\text{e}}^{-{{\nu}^{2}}/2}}\;\text{H}{{\text{e}}_{m}}(\nu )\text{H}{{\text{e}}_{n}}(\nu )\,\text{d}\nu =\sqrt{2\pi}\;n!\,{{\delta}_{mn}}\;,\end{eqnarray} \tag{ 27 }$$

one derives from (21) the following integral representation for the a_P(n)'s ($n\geqslant 3$):

$$\begin{eqnarray}&&{{a}_{P}}(n)=\frac{1}{\sqrt{2\pi}}{\int}_{-\infty}^{\infty}{{ \Delta }_{P}}\left({{\sigma}_{0}}\nu \right)\;\text{H}{{\text{e}}_{n}}(\nu )\,\text{d}\nu \;.\end{eqnarray} \tag{ 28 }$$

If the discrepancy function ${{ \Delta }_{P}}$ is known, one can compute from (28) the non-Gaussianity parameters a_P(n) and then, from (25) and (26), the normalized cumulants ${{\gamma}_{1}},{{\gamma}_{2}},{{C}_{n}}$. Vice versa, one can compute the a_P(n)'s from (25) 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 $\text{H}{{\text{e}}_{2n+1}}(0)=0$ and $\text{H}{{\text{e}}_{2n}}(0)=\left[{{(-1)}^{n}}(2n)!\right]/{{2}^{n}}n!$, we obtain from (21) and (25) for the discrepancy function at $\tau =0$:

$$\begin{eqnarray}\begin{array}{*{35}{l}} \quad \quad \quad {{ \Delta }_{P}}(0) & =\underset{n=2}{\overset{\infty}{\sum}}\,\frac{{{(-1)}^{n}}{{a}_{P}}(2n)}{{{2}^{n}}\,n!}=\frac{{{a}_{P}}(4)}{8}-\frac{{{a}_{P}}(6)}{48}+\frac{{{a}_{P}}(8)}{384}-\frac{{{a}_{P}}(10)}{3840}+\frac{{{a}_{P}}(12)}{46\,080}\mp \cdots \\ {} & =\frac{{{\gamma}_{2}}}{8}-\frac{5}{24}\gamma _{1}^{2}-\frac{{{C}_{6}}}{48}+\frac{35}{384}\gamma _{2}^{2}+\frac{7}{48}{{\gamma}_{1}}{{C}_{5}}+\frac{{{C}_{8}}}{384}\mp \cdots \;. \end{array}\end{eqnarray} \tag{ 29 }$$

In table 2 we present the values for the first a_P(n)'s (see equations (25) and (28)), respectively for ${{\gamma}_{1}}$, ${{\gamma}_{2}}$, and C₅ for the $\Lambda$CDM model (model and computational details are given in appendix A).

Table 2. Table of coefficients a_P(n) for the $\Lambda$CDM model, computed from table 1 using equation (25) in 1st line, then using equation (28) in 2nd line.

$\Lambda$CDM sample Full individual map range, no mask, 2°fwhm, bin 13μK see equations (25), (28) and model in appendix A
n	3	4	5	6	7	8	9	10	11	12
aP(n)	$-5.1872\times {{10}^{-4}}$	$5.825\times {{10}^{-5}}$	$-1.6184\times {{10}^{-4}}$	$-2.411\,12\times {{10}^{-3}}$	$2.9116\times {{10}^{-4}}$	$5.130\,35\times {{10}^{-3}}$	$1.153\,302\times {{10}^{-2}}$	$-1.964\,794\times {{10}^{-2}}$	$7.123\,993\times {{10}^{-2}}$	$0.416\,229\,08$
aP(n)	$-5.1872\times {{10}^{-4}}$	$5.830\times {{10}^{-5}}$	$-1.6184\times {{10}^{-4}}$	$-2.409\,19\times {{10}^{-3}}$	$2.9116\times {{10}^{-4}}$	$5.194\,01\times {{10}^{-3}}$	$1.153\,302\times {{10}^{-2}}$	$-1.785\,385\times {{10}^{-2}}$	$7.123\,994\times {{10}^{-2}}$	$0.457\,217\,93$

$\Lambda$CDM sample Equal temperature range (ETR  ±  201 μK), U73 mask, 2°fwhm, bin 6μK
n	0	1	2	3	4	5	6	7	8
aP(n)				$-5.4435\times {{10}^{-4}}$	$-7.777\,654\times {{10}^{-2}}$	$1.012\,33\times {{10}^{-3}}$	$-2.973\,6480\times {{10}^{-1}}$	$1.7101\times {{10}^{-4}}$	$1.381\,659\,48$
aP(n)	$9.9534\times {{10}^{-4}}$	$-2.8\times {{10}^{-7}}$	$1.002\,019\times {{10}^{-2}}$	$-4.7175\times {{10}^{-4}}$	$5.203\,85\times {{10}^{-3}}$	$1.004\,69\times {{10}^{-3}}$	$-1.989\,114\times {{10}^{-2}}$	$9.19\times {{10}^{-6}}$	$-2.112\,9824\times {{10}^{-1}}$

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

Figure 2 displays the discrepancy function ${{ \Delta }_{P}}$ of the averaged PDF (10⁵ 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 (A.16). This figure also displays the expansion in Hermite polynomials (equation (21) limited to the order 8).

Zoom In Zoom Out Reset image size

Figure 3 shows the envelope of the 10⁵ discrepancy functions of the map sample.

Zoom In Zoom Out Reset image size

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}}\left(\nu :=\tau /{{\sigma}_{0}}\right)$:

$$\begin{eqnarray}&&{{\mathcal{Q}}_{\nu}}:=\left\{{{\boldsymbol {\hat n}}}\,\in {{\mathcal{S}}^{2}}~|~\delta T\left({{\boldsymbol {\hat n}}}\,\right)\geqslant {{\sigma}_{0}}\nu \right\}\;.\end{eqnarray} \tag{ 30 }$$

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

$$\begin{eqnarray}&&{{V}_{0}}\left({{\mathcal{Q}}_{\nu}}\right):={{{\int}^{}}_{{{\mathcal{Q}}_{\nu}}}}\;\text{d}a=\text{area}\left({{\mathcal{Q}}_{\nu}}\right)\;,\end{eqnarray} \tag{ 31 }$$

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

$$\begin{eqnarray}&&{{\text{v}}_{0}}(\nu ):=\frac{1}{4\pi}{{V}_{0}}(\nu )=\text{prob}\left(\delta T\geqslant {{\sigma}_{0}}\nu \right)={{F}_{C}}\left({{\sigma}_{0}}\nu \right)={\int}_{-\infty}^{\infty} \Theta \left({{\tau}^{\prime}}-{{\sigma}_{0}}\nu \right)P\left({{\tau}^{\prime}}\right)\;\text{d}{{\tau}^{\prime}}\;,\end{eqnarray} \tag{ 32 }$$

where F_C is the complementary cumulative distribution function (2). From the definition (32) follows also the relation

$$\begin{eqnarray}&&P(\tau )=-\frac{1}{{{\sigma}_{0}}}\frac{\text{d}{{\text{v}}_{0}}(\nu )}{\text{d}\nu}{{|}_{\nu =\tau /{{\sigma}_{0}}}}\;.\end{eqnarray} \tag{ 33 }$$

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

$$\begin{eqnarray}&&\text{v}_{0}^{\text{G}}(\nu ):=\frac{1}{2}\text{erfc}\left(\frac{\nu}{\sqrt{2}}\right)\;,\end{eqnarray} \tag{ 34 }$$

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

$$\begin{eqnarray}&&\text{v}_{0}^{\text{G}}(\nu )=\frac{1}{\sqrt{2\pi}}\frac{{{\text{e}}^{-{{\nu}^{2}}/2}}}{\nu}\left[1-\frac{1}{{{\nu}^{2}}}+\mathcal{O}\left(\frac{1}{{{\nu}^{4}}}\right)\right].\end{eqnarray} \tag{ 35 }$$

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

$$\begin{eqnarray}&&{{ \Delta }_{0}}(\nu ):=\sqrt{2\pi}\left({{\text{v}}_{0}}(\nu )-\text{v}_{0}^{\text{G}}(\nu )\right)\;,\end{eqnarray} \tag{ 36 }$$

which in turn is expanded into Hermite polynomials,

$$\begin{eqnarray}&&{{ \Delta }_{0}}(\nu ):={{\text{e}}^{-{{\nu}^{2}}/2}}\,\underset{n=2}{\overset{\infty}{\sum}}\,\frac{{{a}_{0}}(n)}{n!}\,\text{H}{{\text{e}}_{n}}(\nu )\;,\end{eqnarray} \tag{ 37 }$$

with the coefficients

$$\begin{eqnarray}&&{{a}_{0}}(n)=\frac{1}{\sqrt{2\pi}}{\int}_{-\infty}^{\infty}{{ \Delta }_{0}}(\nu )\text{H}{{\text{e}}_{n}}(\nu )\,\text{d}\nu \;.\end{eqnarray} \tag{ 38 }$$

The above coefficients ${{a}_{0}}(n),~n\geqslant 2$, measure a possible non-Gaussianity described by the MF ${{\text{v}}_{0}}(\nu )$. With

$$\begin{eqnarray}&&\frac{\text{d}}{\text{d}\nu}\left[{{\text{e}}^{-{{\nu}^{2}}/2}}~\text{H}{{\text{e}}_{n}}(\nu )\right]=-{{\text{e}}^{-{{\nu}^{2}}/2}}~\text{H}{{\text{e}}_{n+1}}(\nu )\;,\end{eqnarray} \tag{ 39 }$$

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

$$\begin{eqnarray}&&\frac{\text{d}}{\text{d}\nu}\,{{ \Delta }_{0}}(\nu )=-{{ \Delta }_{P}}\left({{\sigma}_{0}}\nu \right)\;,\end{eqnarray} \tag{ 40 }$$

one obtains the following relation between the coefficients a₀(n) and a_P(n):

$$\begin{eqnarray}&&{{a}_{0}}(n)=\frac{{{a}_{P}}(n+1)}{n+1}\;\;\;(n\geqslant 2)\;,\end{eqnarray} \tag{ 41 }$$

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

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & \quad \quad {{a}_{0}}(2)=\frac{{{\gamma}_{1}}}{3}\quad ;\quad {{a}_{0}}(3)=\frac{{{\gamma}_{2}}}{4}\quad ;\quad {{a}_{0}}(4)=\frac{{{C}_{5}}}{5}\quad ;\quad {{a}_{0}}(5)=\frac{5}{3}\gamma _{1}^{2}+\frac{{{C}_{6}}}{6}\quad ; \\ {} & \quad \quad {{a}_{0}}(6)=5{{\gamma}_{1}}{{\gamma}_{2}}+\frac{{{C}_{7}}}{7}\quad ;\quad {{a}_{0}}(7)=7{{\gamma}_{1}}{{C}_{5}}+\frac{35}{8}\gamma _{2}^{2}+\frac{{{C}_{8}}}{8}\quad ; \\ {} & \quad \quad {{a}_{0}}(8)=\frac{280}{9}\gamma _{1}^{3}+14{{\gamma}_{2}}{{C}_{5}}+\frac{28}{3}{{\gamma}_{1}}{{C}_{6}}+\frac{{{C}_{9}}}{9}\quad ; \\ {} & \quad \quad {{a}_{0}}(9)=210\gamma _{1}^{2}{{\gamma}_{2}}+12{{\gamma}_{1}}{{C}_{7}}+21{{\gamma}_{2}}{{C}_{6}}+\frac{63}{5}C_{5}^{2}+\frac{{{C}_{10}}}{10}\quad ; \\ {} & \quad \quad {{a}_{0}}(10)=525{{\gamma}_{1}}\gamma _{2}^{2}+420\gamma _{1}^{2}{{C}_{5}}+15{{\gamma}_{1}}{{C}_{8}}+30{{\gamma}_{2}}{{C}_{7}}+42{{C}_{5}}{{C}_{6}}+\frac{{{C}_{11}}}{11}\quad ; \\ {} & \quad \quad {{a}_{0}}(11)=\frac{3850}{3}\gamma _{1}^{4}+\frac{1925}{4}\gamma _{2}^{3}+2310{{\gamma}_{1}}{{\gamma}_{2}}{{C}_{5}} \\ {} & \quad \quad \qquad \quad \,\,+\,770\gamma _{1}^{2}{{C}_{6}}+\frac{55}{3}{{\gamma}_{1}}{{C}_{9}}+\frac{165}{4}{{\gamma}_{2}}{{C}_{8}}+66{{C}_{5}}{{C}_{7}}+\frac{77}{2}C_{6}^{2}+\frac{{{C}_{12}}}{12}\;. \end{array}\end{eqnarray} \tag{ 42 }$$

In table 3 we present the values for the first a₀(n)'s of the $\Lambda$CDM model.

Table 3. Table of coefficients a₀(n), computed from table 1 using equation (42).

$\Lambda$CDM Full individual map range, no mask, 2°fwhm, bin 13μK see equation (42) and model in appendix A
n	0	1	2	3	4	5	6	7	8
a0(n)			$-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.281\,45\times {{10}^{-3}}$

$\Lambda$CDM Equal temperature range (ETR  ±  201 μK), U73 mask, 2°fwhm, bin 6μK
n	0	1	2	3	4	5	6	7	8
a0(n)			$-1.8145\times {{10}^{-4}}$	$-1.944\,413\times {{10}^{-2}}$	$2.0247\times {{10}^{-4}}$	$-4.956\,08\times {{10}^{-3}}$	$2.443\times {{10}^{-5}}$	$1.727\,0744\times {{10}^{-1}}$	$-1.548\,58\times {{10}^{-3}}$

Figure 4 shows the first Minkowski Functional ${{\text{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₀(n).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 appendix 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 sections 2.2 and 2.3 we considered polynomial approximations of degree N in τ, respectively in ν 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 (C.7) and (C.8),

$$\begin{eqnarray}&&{{\mathcal{E}}_{N}}:=\,\parallel {{h}_{0}}-\underset{n=2}{\overset{N}{\sum}}\,\frac{{{a}_{0}}(n)}{n!}\text{H}{{\text{e}}_{n}}{{\parallel}^{2}}\,=\,\parallel {{h}_{0}}{{\parallel}^{2}}-\underset{n=2}{\overset{N}{\sum}}\,\frac{{{\left({{a}_{0}}(n)\right)}^{2}}}{n!},\end{eqnarray} \tag{ 43 }$$

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

$$\begin{eqnarray}&&\underset{n=2}{\overset{\infty}{\sum}}\,\frac{{{\left({{a}_{0}}(n)\right)}^{2}}}{n!}=\,\parallel {{h}_{0}}{{\parallel}^{2}}:=\frac{1}{\sqrt{2\pi}}{\int}_{-\infty}^{+\infty}{{\text{e}}^{-{{\nu}^{2}}/2}}{{\left({{h}_{0}}(\nu )\right)}^{2}}\text{d}\nu,\end{eqnarray} \tag{ 44 }$$

(here we use the weight function $w(\nu )=\frac{1}{\sqrt{2\pi}}{{\text{e}}^{-{{\nu}^{2}}/2}}$ in the definition of the inner product of the Hilbert space $\text{I}\text{H}={{L}^{2}}\left(\text{I}\text{R},w(\nu )\text{d}\nu \right)$, see appendix 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₅ and C₆ by

$$\begin{eqnarray}&&{{\mathcal{E}}_{5}}=\,\parallel {{h}_{0}}{{\parallel}^{2}}-\left[\frac{\gamma _{1}^{2}}{18}+\frac{\gamma _{2}^{2}}{96}+\frac{C_{5}^{2}}{600}+\frac{{{\left(10\gamma _{1}^{2}+{{C}_{6}}\right)}^{2}}}{4320}\right],\end{eqnarray} \tag{ 45 }$$

using the coefficients a₀(2, 3, 4, 5) of (42). 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.

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\geqslant 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 ν 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 $\text{HO}(J),J\geqslant 1$). It turns out that the perturbation expansion $\text{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₀(J)  :=  3J  −  1, (ii) a truncation of the expansion coefficients a₀(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'),

$$\begin{eqnarray}&&{{S}_{n}}:=\frac{{{C}_{n}}}{\sigma _{0}^{n-2}},n\geqslant 3,\end{eqnarray} \tag{ 46 }$$

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\geqslant 2$,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{a}_{0}}(n)=\frac{{{a}_{P}}(n+1)}{n+1}=\frac{1}{n+1}{{B}_{n+1}}\left(0,0,{{\gamma}_{1}},{{\gamma}_{2}},{{C}_{5}},{{C}_{6}},\cdots,{{C}_{n+1}}\right) \\ \,\quad \\ \;=\frac{1}{n+1}{{B}_{n+1}}\left(0,0,{{S}_{3}}{{\sigma}_{0}},{{S}_{4}}\sigma _{0}^{2},\cdots,{{S}_{n+1}}\sigma _{0}^{n-1}\right),\quad i.e. \\ {{a}_{0}}(n)={{D}_{n-1}}\left({{\sigma}_{0}}\right),n\geqslant 2, \end{array}\end{eqnarray} \tag{ 47 }$$

where D_n(z) is a polynomial of degree n in z with D₀(z)  :=  1 and

$$\begin{eqnarray}&&{{D}_{n}}(-z)={{(-1)}^{n}}{{D}_{n}}(z).\end{eqnarray} \tag{ 48 }$$

Thus, a₀(n) is a polynomial of degree n  −  1 in ${{\sigma}_{0}}$,

$$\begin{eqnarray}&&{{a}_{0}}(n)=\underset{j={{j}_{0}}(n)}{\overset{n-1}{\sum}}\,{{D}_{n-1,j}}\left({{S}_{3}},{{S}_{4}},\cdots,{{S}_{n+1}}\right)\sigma _{0}^{j}\;,\;n\geqslant 2.\end{eqnarray} \tag{ 49 }$$

Here, j₀(n) denotes the smallest power of ${{\sigma}_{0}}$ which contributes to a₀(n) and is given for $n\geqslant 2$ by

$$\begin{eqnarray}&&{{j}_{0}}(n):=n(k,l)-(2l+1)=k+l-1\geqslant l+1\;,\end{eqnarray} \tag{ 50 }$$

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:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{D}_{n}}(z) & =\frac{{{S}_{n+2}}}{n+2}{{z}^{n}}~(n=1,2,3)\;\;;\;\;{{D}_{4}}(z)=\frac{5}{3}S_{3}^{2}{{z}^{2}}+\frac{{{S}_{6}}}{6}{{z}^{4}}\;\;; \\ {{D}_{5}}(z) & =5{{S}_{3}}{{S}_{4}}{{z}^{3}}+\frac{{{S}_{7}}}{7}{{z}^{5}}\;\;;\;\;{{D}_{6}}(z)=\left(7{{S}_{3}}{{S}_{5}}+\frac{35}{8}S_{4}^{2}\right){{z}^{4}}+\frac{{{S}_{8}}}{8}{{z}^{6}}\;. \end{array}\end{eqnarray} \tag{ 51 }$$

(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\geqslant 1$.) Explicit expressions for the polynomials D_n(z) for $n\geqslant 7$ are easily obtained either from the recurrence relation (B.11) or from the combinatorial expression (B.12).

In order to derive the perturbative expansion $\text{HO}(J)$ for the discrepancy function ${{ \Delta }_{0}}$ at any order $J\geqslant 1$, we insert in the Hermite expansion (37) the polynomial relation (47) for the expansion coefficients a₀(n). At lowest order, J  =  1, one immediatly obtains from (47) and (49)–(51) the simple result:

$$\begin{eqnarray} &&\Delta _{0}^{\text{HO}(1)}(\nu )={{\text{e}}^{-{{\nu}^{2}}/2}}\frac{{{a}_{0}}(2)}{2}\text{H}{{\text{e}}_{2}}(\nu )=\frac{{{\gamma}_{1}}}{6}\left({{\nu}^{2}}-1\right){{\text{e}}^{-{{\nu}^{2}}/2}}\;,\end{eqnarray} \tag{ 52 }$$

which is completely determined by the skewness coefficient ${{\gamma}_{1}}={{S}_{3}}{{\sigma}_{0}}$. To obtain the $\text{HO}(J)$-expansion for $J\geqslant 2$, we decompose the exact (general) Hermite expansion into three terms (M₀(J)  =  3J  −  1):

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{ \Delta }_{0}}(\nu )= & {{\text{e}}^{-{{\nu}^{2}}/2}}\left[\underset{n=2}{\overset{J+1}{\sum}}\,\frac{{{D}_{n-1}}\left({{\sigma}_{0}}\right)}{n!}\text{H}{{\text{e}}_{n}}(\nu )+\underset{n=J+2}{\overset{{{M}_{0}}(J)}{\sum}}\,\frac{{{D}_{n-1}}\left({{\sigma}_{0}}\right)}{n!}\text{H}{{\text{e}}_{n}}(\nu ) \right. \\ {} & \left. +\underset{n={{M}_{0}}(J)+1}{\overset{\infty}{\sum}}\,\frac{{{D}_{n-1}}\left({{\sigma}_{0}}\right)}{n!}\text{H}{{\text{e}}_{n}}(\nu )\right]. \end{array}\end{eqnarray} \tag{ 53 }$$

Since ${{D}_{n-1}}\left({{\sigma}_{0}}\right)$ is a polynomial of degree n  −  1 in ${{\sigma}_{0}}$, the first sum in (53) represents (having ν fixed) a polynomial of degree J in ${{\sigma}_{0}}$ and thus contributes to the $\text{HO}(J)$-expansion. The last infinite series in (53), which is absolutely convergent, does not contribute at all to the $\text{HO}(J)$-expansion, since the smallest power of ${{\sigma}_{0}}$ appearing in this series is already larger than J; (with M₀(J)  +  1  =  n(3, J  −  1), equation (50) gives ${{j}_{0}}\left({{M}_{0}}(J)+1\right)=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₀(J)  =  3J  −  1. It remains to discuss the second finite sum in (53) where the summation runs over $J+2\leqslant n\leqslant {{M}_{0}}(J)$. Inspection of the expression (49) for the polynomials ${{D}_{n-1}}\left({{\sigma}_{0}}\right)$ (see also the explicit expressions (51)) shows that they will contribute (for $n\geqslant J+2$ and ${{j}_{0}}(n)\leqslant J$) to this sum not only with powers $j\leqslant J$, but also with higher powers that are not admitted in the $\text{HO}(J)$-expansion. Thus, the polynomials D_n−1 have to be replaced by truncated ones. In the case where j₀(n)  >  J, the polynomials D_n−1 do not contribute at all. This leads us to define, for $n\geqslant J+2$, the truncated polynomials $D_{n-1}^{\text{HO}(J)}(z)$, which now also depend on J:

$$\begin{eqnarray}D_{n-1}^{\text{HO}(J)}(z):=\left\{\begin{array}{*{35}{l}} {} & {\sum}_{j={{j}_{0}}(n)}^{J}{{D}_{n-1,j}}\left({{S}_{3}},{{S}_{4}},\cdots,{{S}_{n+1}}\right){{z}^{j}} &,\;{{j}_{0}}(n)\leqslant J \\ {} & 0 &,\;{{j}_{0}}(n)>J. \end{array}\right. \end{eqnarray} \tag{ 54 }$$

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

J  =  2, M₀(2)  =  5:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & D_{3}^{\text{HO}(2)}(z)=0 \\ {} & D_{4}^{\text{HO}(2)}(z)=\frac{5}{3}~S_{3}^{2}{{z}^{2}}\;\;; \end{array}\end{eqnarray} \tag{ 55 }$$

J  =  3, M₀(3)  =  8:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & D_{4}^{\text{HO}(3)}(z)=D_{4}^{\text{HO}(2)}(z) \\ {} & D_{5}^{\text{HO}(3)}(z)=5~{{S}_{3}}{{S}_{4}}{{z}^{3}} \\ {} & D_{6}^{\text{HO}(3)}(z)=0 \\ {} & D_{7}^{\text{HO}(3)}(z)=\frac{280}{9}~S_{3}^{2}{{z}^{3}}\;\;; \end{array}\end{eqnarray} \tag{ 56 }$$

J  =  4, M₀(4)  =  11:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & D_{5}^{\text{HO}(4)}(z)=D_{5}^{\text{HO}(3)}(z) \\ {} & D_{6}^{\text{HO}(4)}(z)=\left(7~{{S}_{3}}{{S}_{5}}+\frac{35}{8}~S_{4}^{2}\right){{z}^{4}} \\ {} & D_{7}^{\text{HO}(4)}(z)=D_{7}^{\text{HO}(3)}(z) \\ {} & D_{8}^{\text{HO}(4)}(z)=210~S_{3}^{2}{{S}_{4}}{{z}^{4}} \\ {} & D_{9}^{\text{HO}(4)}(z)=0 \\ {} & D_{10}^{\text{HO}(4)}(z)=\frac{3850}{3}~S_{3}^{4}{{z}^{4}}\;. \end{array}\end{eqnarray} \tag{ 57 }$$

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

$$\begin{eqnarray} &&\Delta _{0}^{\text{HO}(J)}(\nu ):={{\text{e}}^{-{{\nu}^{2}}/2}}\underset{n=2}{\overset{{{M}_{0}}(J)}{\sum}}\,\frac{a_{0}^{\text{HO}}(J,n)}{n!}\text{H}{{\text{e}}_{n}}(\nu ),\end{eqnarray} \tag{ 58 }$$

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

$$\begin{eqnarray}a_{0}^{\text{HO}}(J,n):=\left\{\begin{array}{*{35}{l}} {} & {{a}_{0}}(n)={{D}_{n-1}}\left({{\sigma}_{0}}\right) & \text{, for}~2\leqslant n\leqslant J+1 \\ {} & D_{n-1}^{\text{HO}(J)}\left({{\sigma}_{0}}\right) & \text{, for}~J+2\leqslant n\leqslant {{M}_{0}}(J). \end{array}\right. \end{eqnarray} \tag{ 59 }$$

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

J  =  2, M₀(2)  =  5:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & a_{0}^{\text{HO}}(2,2)={{a}_{0}}(2)=\frac{{{\gamma}_{1}}}{3} \\ {} & a_{0}^{\text{HO}}(2,3)={{a}_{0}}(3)=\frac{{{\gamma}_{2}}}{4} \\ {} & a_{0}^{\text{HO}}(2,4)=0 \\ {} & a_{0}^{\text{HO}}(2,5)=\frac{5}{3}\gamma _{1}^{2}\;. \end{array}\end{eqnarray} \tag{ 60 }$$

J  =  3, M₀(3)  =  8:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & a_{0}^{\text{HO}}(3,2)={{a}_{0}}(2)=\frac{{{\gamma}_{1}}}{3} \\ {} & a_{0}^{\text{HO}}(3,3)={{a}_{0}}(3)=\frac{{{\gamma}_{2}}}{4} \\ {} & a_{0}^{\text{HO}}(3,4)={{a}_{0}}(4)=\frac{{{C}_{5}}}{5} \\ {} & a_{0}^{\text{HO}}(3,5)=\frac{5}{3}\gamma _{1}^{2} \\ {} & a_{0}^{\text{HO}}(3,6)=5{{\gamma}_{1}}{{\gamma}_{2}} \\ {} & a_{0}^{\text{HO}}(3,7)=0 \\ {} & a_{0}^{\text{HO}}(3,8)=\frac{280}{9}\gamma _{1}^{3}\;. \end{array}\end{eqnarray} \tag{ 61 }$$

J  =  4, M₀(4)  =  11:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & a_{0}^{\text{HO}}(4,n)={{a}_{0}}(n)\;\;\;(n=2,3,4,5) \\ {} & a_{0}^{\text{HO}}(4,6)=5{{\gamma}_{1}}{{\gamma}_{2}} \\ {} & a_{0}^{\text{HO}}(4,7)=7{{\gamma}_{1}}{{C}_{5}}+\frac{35}{8}\gamma _{2}^{2} \\ {} & a_{0}^{\text{HO}}(4,8)=\frac{280}{9}\gamma _{1}^{3} \\ {} & a_{0}^{\text{HO}}(4,9)=210\gamma _{1}^{2}{{\gamma}_{2}} \\ {} & a_{0}^{\text{HO}}(4,10)=0 \\ {} & a_{0}^{\text{HO}}(4,11)=\frac{3850}{3}\gamma _{1}^{4}\;. \end{array}\end{eqnarray} \tag{ 62 }$$

(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₀(J)). If we insert, however, for the Hermite coefficients $a_{0}^{\text{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:

$$\begin{eqnarray} &&\Delta _{0}^{\text{HO}(J)}(\nu )={{\text{e}}^{-{{\nu}^{2}}/2}}\underset{j=1}{\overset{J}{\sum}}\,\text{w}_{0}^{\text{HO}(J)}(\,j,\nu )\sigma _{0}^{j}\;,\end{eqnarray} \tag{ 63 }$$

which has now (at fixed ν) the form of a power series in ${{\sigma}_{0}}$ with coefficient functions $\text{w}_{0}^{\text{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 $\text{w}_{0}^{\text{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:

$$\begin{eqnarray}\begin{array}{*{35}{l}} J=1: & \text{w}_{0}^{\text{HO}(1)}(1,\nu )=\frac{{{S}_{3}}}{6}\text{H}{{\text{e}}_{2}}(\nu ), \end{array}\end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} J=2: & \text{w}_{0}^{\text{HO}(2)}(1,\nu )=\text{w}_{0}^{\text{HO}(1)}(1,\nu ) \\ {} & \text{w}_{0}^{\text{HO}(2)}(2,\nu )=\frac{{{S}_{4}}}{24}\text{H}{{\text{e}}_{3}}(\nu )+\frac{S_{3}^{2}}{72}\text{H}{{\text{e}}_{5}}(\nu ), \end{array}\end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} J=3: & \text{w}_{0}^{\text{HO}(3)}(1,\nu )=\text{w}_{0}^{\text{HO}(1)}(1,\nu ) \\ {} & \text{w}_{0}^{\text{HO}(3)}(2,\nu )=\text{w}_{0}^{\text{HO}(2)}(2,\nu ) \\ {} & \text{w}_{0}^{\text{HO}(3)}(3,\nu )=\frac{{{S}_{5}}}{120}\text{H}{{\text{e}}_{4}}(\nu )+\frac{{{S}_{3}}{{S}_{4}}}{144}\text{H}{{\text{e}}_{6}}(\nu )+\frac{S_{3}^{2}}{1296}\text{H}{{\text{e}}_{8}}(\nu ), \end{array}\end{eqnarray} \tag{ 66 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} J=4: & \text{w}_{0}^{\text{HO}(4)}(1,\nu )= & \text{w}_{0}^{\text{HO}(1)}(1,\nu ) \\ {} & \text{w}_{0}^{\text{HO}(4)}(2,\nu )= & \text{w}_{0}^{\text{HO}(2)}(2,\nu ) \\ {} & \text{w}_{0}^{\text{HO}(4)}(3,\nu )= & \text{w}_{0}^{\text{HO}(3)}(3,\nu ) \\ {} & \text{w}_{0}^{\text{HO}(4)}(4,\nu )= & \frac{{{S}_{6}}}{720}\text{H}{{\text{e}}_{5}}(\nu )+\frac{1}{720}\left({{S}_{3}}{{S}_{5}}+\frac{5}{8}S_{4}^{2}\right)\text{H}{{\text{e}}_{7}}(\nu ) \\ {} & {} & +\,\frac{S_{3}^{2}{{S}_{4}}}{1728}\text{H}{{\text{e}}_{9}}(\nu )+\frac{S_{3}^{4}}{31\,104}\text{H}{{\text{e}}_{11}}(\nu ). \end{array}\end{eqnarray} \tag{ 67 }$$

The second-order formulae (65) agree with Matsubara's result [57, 58] (who writes S  :=  S₃ and K  :=  S₄).

It is worthwhile to mention that the perturbative formula (63) closely resembles the so-called Edgeworth expansion (C.21) which gives a refinement of the classical central limit theorem (for details and references, see appendix C). However, (C.21) 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\to \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, (60) respectively (61), with the complete NG-coefficients a₀(n) given in equation (42). At second order (J  =  2) we see that the first two coefficients are identical to the complete coefficients a₀(n), i.e. $a_{0}^{\text{HO}}(2,n)={{a}_{0}}(n)$ for n  =  2, 3. The coefficient $a_{0}^{\text{HO}}(2,4)$ vanishes because its complete value C₅/5 is of third order. Finally, the last coefficient $a_{0}^{\text{HO}}(2,5)$ differs from a₀(5) by the term C₆/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}^{\text{HO}}(3,n)={{a}_{0}}(n)$ for n  =  2, 3, 4; $a_{0}^{\text{HO}}(3,5)=a_{0}^{\text{HO}}(2,5)$ still does not contain the term C₆/6; the coefficient $a_{0}^{\text{HO}}(3,7)$ vanishes because the complete value for a₀(7) is of order $\sigma _{0}^{4}$; $a_{0}^{\text{HO}}(3,6)$ does not contain the term C₇/7 (of order 5), and $a_{0}^{\text{HO}}(3,8)$ does not contain the terms of order 5, respectively 7. Note, in particular, that a vanishing coefficient $a_{0}^{\text{HO}}(J,n)$ implies that the associated contribution from the Hermite polynomial $\text{H}{{\text{e}}_{n}}(\nu )$ is absent in $\Delta _{0}^{\text{HO}\left(\text{J}\right)}$ compared to ${{ \Delta }_{0}}$. This is the case at second order with $\text{H}{{\text{e}}_{4}}(\nu )$, and at third order with $\text{H}{{\text{e}}_{7}}(\nu )$. Since the $\text{H}{{\text{e}}_{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 $\text{HO}(J)$-expansion:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\mathcal{E}}^{\text{HO}(J)}} & :=\,\parallel {{h}_{0}}-\underset{n=2}{\overset{{{M}_{0}}(J)}{\sum}}\,\frac{a_{0}^{\text{HO}}(J,n)}{n!}\text{H}{{\text{e}}_{n}}{{\parallel}^{2}} \\ {} & {} & =\,\parallel {{h}_{0}}{{\parallel}^{2}}-\underset{n=2}{\overset{{{M}_{0}}(J)}{\sum}}\,\frac{{{\left({{a}_{0}}(n)\right)}^{2}}}{n!}+\underset{n=2}{\overset{{{M}_{0}}(J)}{\sum}}\,\frac{{{\left({{a}_{0}}(n)-a_{0}^{\text{HO}}(J,n)\right)}^{2}}}{n!}, \end{array}\end{eqnarray} \tag{ 68 }$$

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

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{E}}^{\text{HO}(J)}}={{\mathcal{E}}_{{{M}_{0}}(J)}}+\underset{n=J+2}{\overset{{{M}_{0}}(J)}{\sum}}\,\frac{{{\left({{a}_{0}}(n)-a_{0}^{\text{HO}}(J,n)\right)}^{2}}}{n!}. \end{array}\end{eqnarray} \tag{ 69 }$$

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

$$\begin{eqnarray}&&{{\mathcal{E}}^{\text{HO}(J)}}>{{\mathcal{E}}_{{{M}_{0}}(J)}},\end{eqnarray} \tag{ 70 }$$

i.e. the mean square error of the $\text{HO}(J)$-expansion is always larger than the mean square error of the complete Hermite expansion truncated at n  =  M₀(J). For instance, at second-order of perturbation theory, J  =  2, we have ${{\mathcal{E}}^{\text{HO}(2)}}>{{\mathcal{E}}_{5}}$, where ${{\mathcal{E}}_{5}}$ is explicitly given in (45). Precisely, we obtain from (69) through (42) and (60):

$$\begin{eqnarray}&&{{\mathcal{E}}^{\text{HO}(2)}}={{\mathcal{E}}_{5}}+\frac{C_{5}^{2}}{600}+\frac{C_{6}^{2}}{4320}.\end{eqnarray} \tag{ 71 }$$

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

2.5.1. 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 (25) and table 1) for n  =  3 to 6, according to the following equation:

$$\begin{eqnarray}&&{{ \Delta }_{P}}(\tau )={{\text{e}}^{-{{\tau}^{2}}/2\sigma _{0}^{2}}}\;\underset{n=3}{\overset{6}{\sum}}\,\frac{{{a}_{P}}(n)}{n!}\;\text{H}{{\text{e}}_{n}}\left(\tau /{{\sigma}_{0}}\right)\;.\end{eqnarray} \tag{ 72 }$$

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray} &&\Delta _{P}^{\text{HO}(4)}(\tau )={{\text{e}}^{-{{\tau}^{2}}/2\sigma _{0}^{2}}}\,\underset{n=3}{\overset{12}{\sum}}\,\frac{a_{P}^{\text{HO}}(4,n)}{n!}\;\text{H}{{\text{e}}_{n}}\left(\tau /{{\sigma}_{0}}\right),\end{eqnarray} \tag{ 73 }$$

where the NG-coefficients $a_{P}^{\text{HO}}(4,n)$ are obtained from equations (41) and (62) 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 appendix C).

According to Hadwiger's theorem [19], there are three independent MFs on ${{\mathcal{S}}^{2}}$. The simplest case V₀, respectively v₀ 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

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\text{v}}_{1}}(\nu ):=\frac{1}{4\pi}{{V}_{1}}(\nu )=\frac{1}{4\pi}\frac{1}{4}{{{\int}^{}}_{\partial {{\mathcal{Q}}_{\nu}}}}\text{d}s=\frac{1}{16\pi}\,\text{length}\;\left(\partial {{\mathcal{Q}}_{\nu}}\right)\;, \end{array}\end{eqnarray} \tag{ 74 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\text{v}}_{2}}(\nu )=\frac{1}{4\pi}{{V}_{2}}(\nu )=\frac{1}{4\pi}\frac{1}{2\pi}{{{\int}^{}}_{\partial {{\mathcal{Q}}_{\nu}}}}\kappa (s)\,\text{d}s\;, \end{array}\end{eqnarray} \tag{ 75 }$$

where $\text{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, ${{\text{v}}_{k}}(\nu ),k=1,2$, have been computed by Tomita [21–23] (we set $\mu =0$):

$$\begin{eqnarray}&&\text{v}_{1}^{\text{G}}(\nu ):=\frac{1}{8\sqrt{2}}\frac{{{\sigma}_{1}}}{{{\sigma}_{0}}}\;{{\text{e}}^{-{{\nu}^{2}}/2}}\quad ;\quad \text{v}_{2}^{\text{G}}(\nu )=\frac{1}{2{{(2\pi )}^{3/2}}}\frac{\sigma _{1}^{2}}{\sigma _{0}^{2}}\;\nu \,{{\text{e}}^{-{{\nu}^{2}}/2}}\;.\end{eqnarray} \tag{ 76 }$$

Here, ${{\sigma}_{0}}$ denotes again the standard deviation of the CMB temperature anisotropy $\delta T\left({{\boldsymbol{\hat n}}}\,\right)$ 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 $\text{d}{{s}^{2}}={{\gamma}_{ij}}\text{d}{{x}^{i}}\text{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=\left({{\nabla}^{1}}\delta T,{{\nabla}^{2}}\delta T\right)$, i.e.

$$\begin{eqnarray}&&\langle {{\nabla}_{i}}\delta T\left({{\boldsymbol{\hat n}}}\,\right){{\nabla}_{j}}\delta T\left({{\boldsymbol{\hat n}}}\,\right)\rangle =\frac{\sigma _{1}^{2}}{2}\,{{\gamma}_{ij}}\;.\end{eqnarray} \tag{ 77 }$$

From the MFs v₁ and v₂ 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):

$$\begin{eqnarray}&&{{{\int}^{}}_{{{\mathcal{Q}}_{\nu}}}}K\,\text{d}a+{\int}_{\partial {{\mathcal{Q}}_{\nu}}}\kappa (s)\text{d}s={{V}_{0}}(\nu )+2\pi {{V}_{2}}(\nu )=2\pi \;\chi (\nu )\;,\end{eqnarray} \tag{ 78 }$$

which gives

$$\begin{eqnarray}&&\chi (\nu )=2{{\text{v}}_{0}}(\nu )+4\pi {{\text{v}}_{2}}(\nu )\quad \text{and}\quad g(\nu )=1-{{\text{v}}_{0}}(\nu )-2\pi {{\text{v}}_{2}}(\nu )\;.\end{eqnarray} \tag{ 79 }$$

As a measure of possible non-Gaussianities based on the MFs ${{\text{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):

$$\begin{eqnarray}&&{{ \Delta }_{k}}(\nu ):=\frac{{{\text{v}}_{k}}(\nu )-\text{v}_{k}^{\text{G}}(\nu )}{{{N}_{k}}}\;\;\;(k=0,1,2)\;,\end{eqnarray} \tag{ 80 }$$

with ${{N}_{0}}=1/\sqrt{2\pi}$, ${{N}_{1}}=\text{max}\left\{\text{v}_{1}^{\text{G}}\right\}=\left[1/\left(8\sqrt{2}\right)\right]{{\sigma}_{1}}/{{\sigma}_{0}}$, and ${{N}_{2}}=\left[1/\left(2{{(2\pi )}^{3/2}}\right)\right]\sigma _{1}^{2}/\sigma _{0}^{2}$. Under the same assumptions made before (for ${{ \Delta }_{P}}$ respectively ${{ \Delta }_{0}}$, see appendix C), we can expand the ${{ \Delta }_{k}}$'s into a convergent Hermite expansion,

$$\begin{eqnarray}&&{{ \Delta }_{k}}(\nu )={{\text{e}}^{-{{\nu}^{2}}/2}}\underset{n={{n}_{k}}}{\overset{\infty}{\sum}}\,\frac{{{a}_{k}}(n)}{n!}\,\text{H}{{\text{e}}_{n}}(\nu )\;\;;\;\;k=0,1,2\;,\end{eqnarray} \tag{ 81 }$$

with the dimensionless NG-coefficients a_k(n) and n₀  =  2, ${{n}_{1}}={{n}_{2}}=0$. With the help of (27) we obtain the integral representation:

$$\begin{eqnarray}&&{{a}_{k}}(n)=\frac{1}{\sqrt{2\pi}}{\int}_{-\infty}^{\infty}{{ \Delta }_{k}}(\nu )\,\text{H}{{\text{e}}_{n}}(\nu )\,\text{d}\nu \;.\end{eqnarray} \tag{ 82 }$$

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),

$$\begin{eqnarray} &&\Delta _{k}^{\text{HO}(J)}(\nu )={{\text{e}}^{-{{\nu}^{2}}/2}}\underset{n={{n}_{k}}}{\overset{{{M}_{k}}(J)}{\sum}}\,\frac{a_{k}^{\text{HO}}(J,n)}{n!}\text{H}{{\text{e}}_{n}}(\nu )\;,\end{eqnarray} \tag{ 83 }$$

for which the second order NG-coefficients $a_{k}^{\text{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₁(2)  =  6 and M₂(2)  =  7.

The NG-coefficients are then given for $\Delta _{1}^{\text{HO}(J)}(\nu )$ as follows:

J  =  2, M₁(2)  =  6:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & a_{1}^{\text{HO}}(2,0)=-\frac{{{K}_{3}}}{16} \\ {} & a_{1}^{\text{HO}}(2,1)=-\frac{{{S}_{1}}}{4} \\ {} & a_{1}^{\text{HO}}(2,2)=-\frac{1}{6}\left({{K}_{1}}+\frac{3}{8}S_{1}^{2}\right) \\ {} & a_{1}^{\text{HO}}(2,3)={{\gamma}_{1}} \\ {} & a_{1}^{\text{HO}}(2,4)={{\gamma}_{2}}-{{\gamma}_{1}}{{S}_{1}} \\ {} & a_{1}^{\text{HO}}(2,5)=0 \\ {} & a_{1}^{\text{HO}}(2,6)=10\gamma _{1}^{2}\;. \end{array}\end{eqnarray} \tag{ 84 }$$

and for $\Delta _{2}^{\text{HO}(J)}(\nu )$:

J  =  2, M₂(2)  =  7:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & a_{2}^{\text{HO}}(2,0)=-{{S}_{2}} \\ {} & a_{2}^{\text{HO}}(2,1)=-\frac{1}{2}\left({{K}_{2}}+{{S}_{1}}{{S}_{2}}\right) \\ {} & a_{2}^{\text{HO}}(2,2)=-{{S}_{1}} \\ {} & a_{2}^{\text{HO}}(2,3)=-{{K}_{1}}-{{\gamma}_{1}}{{S}_{2}} \\ {} & a_{2}^{\text{HO}}(2,4)=4{{\gamma}_{1}} \\ {} & a_{2}^{\text{HO}}(2,5)=5{{\gamma}_{2}}-10{{\gamma}_{1}}{{S}_{1}} \\ {} & a_{2}^{\text{HO}}(2,6)=0 \\ {} & a_{2}^{\text{HO}}(2,7)=70\gamma _{1}^{2}\;. \end{array}\end{eqnarray} \tag{ 85 }$$

Here, we introduced the three dimensionless skewness parameters ${{\gamma}_{1}}$, S₁, S₂ and the four dimensionless kurtosis parameters ${{\gamma}_{2}}$, K₁, K₂ and K₃ using the notation $\tau =\tau \left({{\boldsymbol{\hat n}}}\,\right):=\delta T\left({{\boldsymbol{\hat n}}}\,\right)$ (${{\gamma}_{1}}$ and ${{\gamma}_{2}}$ as in equations (12) and (13), respectively):

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & \quad \quad \quad \quad {{S}_{1}}:=\frac{{{\langle {{\tau}^{2}}\boldsymbol{}{{\nabla}^{2}}\tau \rangle}_{C}}}{{{\sigma}_{0}}\sigma _{1}^{2}}\;\;;\;\;{{S}_{2}}:=\frac{{{\langle |\boldsymbol{}\nabla \tau {{|}^{2}}\boldsymbol{}{{\nabla}^{2}}\tau \rangle}_{C}}{{\sigma}_{0}}}{\sigma _{1}^{4}}\;\;; \\ {} & {{K}_{1}}:=\frac{{{\langle {{\tau}^{3}}\boldsymbol{}{{\nabla}^{2}}\tau \rangle}_{C}}}{\sigma _{0}^{2}\sigma _{1}^{2}}\;\;;\;\;{{K}_{2}}:=\frac{2{{\langle \tau |\boldsymbol{}\nabla \tau {{|}^{2}}\boldsymbol{}{{\nabla}^{2}}\tau \rangle}_{C}}+{{\langle |\boldsymbol{}\nabla \tau {{|}^{4}}\rangle}_{C}}}{\sigma _{1}^{4}}\;\;;\;\;{{K}_{3}}:=\frac{{{\langle |\boldsymbol{}\nabla \tau {{|}^{4}}\rangle}_{C}}}{\sigma _{1}^{4}}. \end{array}\end{eqnarray} \tag{ 86 }$$

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

Table 4 shows the coefficients a₁(n) and a₂(n) calculated from (82).

Table 4. Table of coefficients a₁(n) and a₂(n).

$\Lambda$CDM sample Full individual map range, no mask, 2°fwhm, bin 13 μK see appendix A
n	0	1	2	3	4	5	6	7	8
a1(n)	$-2.21\times {{10}^{-6}}$	$-2.808\times {{10}^{-5}}$	$3.972\,74\times {{10}^{-3}}$	$-5.0873\times {{10}^{-4}}$	$-5.823\times {{10}^{-5}}$	$-6.643\times {{10}^{-5}}$	$-1.979\,22\times {{10}^{-3}}$	$5.059\times {{10}^{-5}}$	$6.374\,08\times {{10}^{-3}}$
a2(n)	$1.927\times {{10}^{-5}}$	$-2.148\times {{10}^{-5}}$	$-1.7313\times {{10}^{-4}}$	$1.183\,631\times {{10}^{-2}}$	$-1.534\,88\times {{10}^{-3}}$	$-4.8941\times {{10}^{-4}}$	$-1.208\,22\times {{10}^{-3}}$	$-9.588\,44\times {{10}^{-3}}$	$-7.333\,97\times {{10}^{-3}}$

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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\text{C}\ell}}$ and ${{\sigma}_{1\text{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 σ-values from the moments of pixels (denoted by subscripts $\text{px}$, $\text{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:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{ \Delta }_{P}}(\tau ):=\frac{P(\tau )-{{P}^{\text{G}}}(\tau )}{{{N}_{\text{P}}}},~ & \text{with}~~{{N}_{\text{P}}}=\frac{1}{{{\sigma}_{0}}\sqrt{2\pi}}\;; \end{array}\end{eqnarray} \tag{ 87 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{ \Delta }_{0}}(\nu ):=\frac{{{\text{v}}_{0}}(\nu )-\text{v}_{0}^{\text{G}}(\nu )}{{{N}_{0}}},\, & \text{with}\,\,{{N}_{0}}=\frac{1}{\sqrt{2\pi}},\, & \,\,\nu =\frac{\tau}{{{\sigma}_{0\text{px}}}}\;; \end{array}\end{eqnarray} \tag{ 88 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{ \Delta }_{1}}(\nu ):=\frac{{{\text{v}}_{1}}(\nu )-\text{v}_{1}^{\text{G}}(\nu )}{{{N}_{1}}},\, & \text{with}\,\,{{N}_{1}}=\frac{1}{8\sqrt{2}}\frac{{{\sigma}_{1\text{px}}}}{{{\sigma}_{0\text{px}}}},\, & \,\,\nu =\frac{\tau}{{{\sigma}_{0\text{px}}}}\;; \end{array}\end{eqnarray} \tag{ 89 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{ \Delta }_{2}}(\nu ):=\frac{{{\text{v}}_{2}}(\nu )-\text{v}_{2}^{\text{G}}(\nu )}{{{N}_{2}}},\, & \text{with}\,\,{{N}_{2}}=\frac{1}{2{{(2\pi )}^{3/2}}}\frac{\sigma _{1\text{px}}^{2}}{\sigma _{0\text{px}}^{2}},\, & \,\,\nu =\frac{\tau}{{{\sigma}_{0\text{px}}}}. \end{array}\end{eqnarray} \tag{ 90 }$$

These discrepancy functions are self-consistent and model-independent given that the terms used, ${{\sigma}_{\text{S}}}$ and ${{\sigma}_{\text{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 $\text{S}$ of maps, we define the sample discrepancy functions this way:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & P_{\text{S}}^{\text{G}}(\tau ):={{P}^{\text{G}}}\left({{\langle \mu \rangle}_{\text{S}}},\sqrt{{{\langle \sigma _{0}^{2}\rangle}_{\text{S}}}},~\tau \right), \end{array}\end{eqnarray} \tag{ 91 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & \text{v}_{i\text{S}}^{\text{G}}(\nu ):=\text{v}_{i}^{\text{G}}\left({{\langle {{\mu}_{\text{px}}}\rangle}_{\text{S}}},\sqrt{{{\langle \sigma _{1\text{px}}^{2}\rangle}_{\text{S}}}},\sqrt{{{\langle \sigma _{0\text{px}}^{2}\rangle}_{\text{S}}}},\,\nu \right). \end{array}\end{eqnarray} \tag{ 92 }$$

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–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 $\text{Df}$) are given by

$$\begin{eqnarray}&&\text{Df}_{i}^{d}(\nu ):=\frac{\text{v}_{i}^{d}(\nu )}{N_{i}^{d}}-{{\langle \frac{{{\text{v}}_{i}}(\nu )}{{{N}_{i}}}\rangle}_{\text{S}}},\end{eqnarray} \tag{ 93 }$$

where the reference map sample $\text{S}$ is Gaussian by construction. Compared with the discrepancy functions, the $\text{Df}$'s measure the distance to a Gaussian premise which may be different in amplitude and shape from our analytical one. The $\text{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 $\text{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.

Figures 11–14 show the statistical and morphological behaviour of each of the 10⁵ 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 \text{K}$. The average quantities μ, ${{\sigma}_{0}}$ and ${{\sigma}_{1}}$, as well as the discrepancy functions are calculated over the largest common and centered (on $\nu =0$) temperature range $\left[-201~\mu \text{K},+201~\mu \text{K}\right]$, 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σ cosmic variance limit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The benchmark study of Planck maps and the $\Lambda$CDM sample with U73 mask (figures 11–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 $100\,000$) peaking beyond the $6\sigma$-envelope. Peaking beyond $6\sigma$ we also count 354 maps for ${{ \Delta }_{0}}$.

Contrary to other works to unveil CMB non-Gaussianity, no specific model of non-Gaussianity has been put into the $100\,000$ maps of the $\Lambda$CDM sample probed all along this study (see appendix 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 $100\,000$ map sample without mask. Our previous analyses using different parameters N_side, l_max, $\text{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 7 yr) for the generation of the power spectrum to get a $100\,000$ map sample. The comparison with the NG of the observed CMB (figures 11–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 appendix 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.

When μ, ${{\sigma}_{0}}$, and ${{\sigma}_{1}}$ are evaluated over $\left[-{{\tau}_{m}},{{\tau}_{m}}\right]$, while the discrepancy functions are calculated over a smaller interval in τ (resp. in ν), this leads to an error on their shape and to an overestimate of the magnitude of NG, as we tested for the $100\,000$ map sample. Such a problem may arise when one works with the normalized temperature $\nu =\tau /{{\sigma}_{0}}$. While ${{\sigma}_{0}}\left({{\tau}_{\text{max}1}}\right)$ is known for the temperature range ${{\tau}_{\text{max}1}}$ giving ${{\nu}_{\text{max}1}}={{\tau}_{\text{max}1}}/{{\sigma}_{0}}\left({{\tau}_{\text{max}1}}\right)$, the need to impose a ν-range (${{\nu}_{\text{max}2}}$) different from ${{\nu}_{\text{max}1}}$ (e.g. for comparing with another map sample) is equivalent to solve for an implicit function f, ${{\sigma}_{0}}\left({{\nu}_{\text{max}2}}\right)=f\left({{\sigma}_{0}}\left({{\nu}_{\text{max}2}}\right)\right)$, 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 ν-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 τ and impose the same range to data and to the simulation maps.

For ${{ \Delta }_{0,1,2}}$ we first calculated μ's and σ's over the full τ-range (±201 $\mu \text{K}$) of the four different Planck maps. For each of these we obtain the following ${{\sigma}_{0}}$ values: NILC 51.532 $\mu \text{K}$, SEVEM 51.750 $\mu \text{K}$, SMICA 51.576 $\mu \text{K}$, Commander—Ruler 51.794 $\mu \text{K}$, while we obtain ${{\sigma}_{0}}=59.294$ $\mu \text{K}$ for the $\Lambda$CDM simulations over exactly the same (±201 $~\mu \text{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 [11, 66–71] and [56, 72]). When not treated correctly, the impact of this anomaly is significant for the calculation of the Gaussian premises ${{P}^{\text{G}}}$, equation (16), and $\text{v}_{k}^{\text{G}}$, equations (34) and (76), as they are functions of the standard deviation ${{\sigma}_{0}}$. Regarding v₁ and v₂ we made two interesting observations: firstly, the Gaussian premises $\text{v}_{1}^{\text{G}}$ and $\text{v}_{2}^{\text{G}}$ 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 \text{K}$. For the simulation sample, ${{\rho}_{ \Lambda CDM}}$ has a much smaller value: 34.162. Secondly, we observed that the ratio ρ is very stable for the sample with U73 mask and for Planck maps with U73 mask when passing from the largest common temperature range (±201 $\mu \text{K}$) to the smallest temperature range covering all the maps (±396.5 $\mu \text{K}$). On the other hand we noticed (see at the beginning of appendix A) that the ratio ρ 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}^{\text{G}}}$ and ${{ \Delta }_{P}}$, (ii) for v₀, $\text{v}_{0}^{\text{G}}$ and ${{ \Delta }_{0}}$ it required some iterations upon different temperature ranges to reach close to the ${{\nu}_{\text{max}}}$ values, and (iii) is simplified for v₁, $\text{v}_{1}^{\text{G}}$, v₂, $\text{v}_{2}^{\text{G}}$, ${{ \Delta }_{1}}$ and ${{ \Delta }_{2}}$ by assuming $\rho ={{\sigma}_{1}}/{{\sigma}_{0}}$ to be constant for the given change of temperature range.

Non-linear correction terms within the standard perturbation approach are commonly investigated by constraining the coefficients (bi-spectrum with f_NL and tri-spectrum with g_NL). Applied to Bardeen's curvature these constraints allow one 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 \left({{\boldsymbol{\hat n}}}\,\right)$ (${{\boldsymbol{\hat n}}}\,={{\boldsymbol{\hat n}}}\,(\vartheta,\varphi )$), is commonly expanded in real space around a Gaussian random field $\phi \left({{\boldsymbol{\hat n}}}\,\right)$ (of vanishing mean), using the expansion:

$$\begin{eqnarray} &&\Phi \left({{\boldsymbol{\hat n}}}\,\right)=\phi \left({{\boldsymbol{\hat n}}}\,\right)+f_{\text{NL}}^{\left(\text{local}\right)}\left({{\phi}^{2}}\left({{\boldsymbol{\hat n}}}\,\right)-\langle {{\phi}^{2}}\left({{\boldsymbol{\hat n}}}\,\right)\rangle \right)+g_{\text{NL}}^{\left(\text{local}\right)}\left({{\phi}^{3}}\left({{\boldsymbol{\hat n}}}\,\right)-\langle {{\phi}^{3}}\left({{\boldsymbol{\hat n}}}\,\right)\rangle \right)+\cdots.\end{eqnarray} \tag{ 94 }$$

In 2010, the limit derived from WMAP 7 yr was $f_{\text{NL}}^{\left(\text{local}\right)}\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_{\text{NL}}^{\left(\text{local}\right)}=2.7\pm 5.8$, $f_{\text{NL}}^{\left(\text{equilateral}\right)}=-42\pm 75$, and $f_{\text{NL}}^{\left(\text{orthogonal}\right)}=-25\pm 39$ [73]. Then, Planck 2015 combined temperature and polarization data provided $f_{\text{NL}}^{\left(\text{local}\right)}=0.8\pm 5.0$, $f_{\text{NL}}^{\left(\text{equilateral}\right)}=-4\pm 43$, and $f_{\text{NL}}^{\left(\text{orthogonal}\right)}=-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.

We consider a Gaussian window function with smoothing scale ${{\theta}_{\text{fwhm}}}\left({{}^{\circ}}\right)$,

$$\begin{eqnarray}&&{{\sigma}_{\text{G}}}:={{\theta}_{\text{fwhm}}}\frac{\pi}{180}\frac{1}{\sqrt{8\ln 2}}.\end{eqnarray} \tag{ A.1 }$$

For $\sigma _{\text{G}}^{2}\ll 1$ the Gaussian kernel is approximated by

$$\begin{eqnarray}&&{{\text{W}}_{\ell}}=\exp \left(-\frac{\sigma _{\text{G}}^{2}}{2}\ell (\ell +1)\right),\end{eqnarray} \tag{ A.2 }$$

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

$$\begin{eqnarray}&&{{\text{C}}_{\ell}}:=\langle \frac{1}{2\ell +1}\underset{m=-\ell}{\overset{+\ell}{\sum}}\,|{{a}_{\ell,m}}{{|}^{2}}\rangle,\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&\text{C}{{\text{w}}_{\ell}}={{\text{C}}_{\ell}}~\text{pxc}{{(\ell )}^{2}}\exp \left(-\ell (\ell +1)\sigma _{\text{G}}^{2}\right),\end{eqnarray} \tag{ A.4 }$$

and the angular power spectrum by

$$\begin{eqnarray}&&{{\left(\delta {{\text{T}}_{\ell}}\right)}^{2}}:=\left(\ell (\ell +1)/2\pi \right)\text{C}{{\text{w}}_{\ell}}.\end{eqnarray} \tag{ A.5 }$$

We consider the mean values,

$$\begin{eqnarray}&&\mu :={{\alpha}_{1}}~;~{{\mu}_{\text{px}}}:={{\langle \delta T\rangle}_{\text{px}}},\end{eqnarray} \tag{ A.6 }$$

and the variances,

$$\begin{eqnarray}&&\sigma _{0}^{2}:={{\alpha}_{2}}-{{\mu}^{2}}~;~\sigma _{0\text{px}}^{2}:={{\langle {{(\delta T)}^{2}}\rangle}_{\text{px}}}-\langle \delta T\rangle _{\text{px}}^{2}.\end{eqnarray} \tag{ A.7 }$$

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

$$\begin{eqnarray}&&\sigma _{0\text{C}\ell}^{2}:=\underset{\ell =2}{\overset{{{\ell}_{\text{max}}}}{\sum}}\,\text{C}{{\text{w}}_{\ell}}\frac{(2\ell +1)}{4\pi}.\end{eqnarray} \tag{ A.8 }$$

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

$$\begin{eqnarray}&&\sigma _{1\text{px}}^{2}:=\langle \left({{\nabla}_{1}}\delta T\left({{\boldsymbol{\hat n}}}\,\right)\right)\left({{\nabla}^{1}}\delta T\left({{\boldsymbol{\hat n}}}\,\right)\right)+\left({{\nabla}_{2}}\delta T\left({{\boldsymbol{\hat n}}}\,\right)\right)\left({{\nabla}^{2}}\delta T\left({{\boldsymbol{\hat n}}}\,\right)\right)\rangle,\end{eqnarray} \tag{ A.9 }$$

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

$$\begin{eqnarray}&&\sigma _{1\text{C}\ell}^{2}:=\underset{\ell =2}{\overset{{{\ell}_{\text{max}}}}{\sum}}\,\text{C}{{\text{w}}_{\ell}}\frac{(2\ell +1)\ell (\ell +1)}{4\pi}.\end{eqnarray} \tag{ A.10 }$$

The measures of the variance and variance of the local gradient, equations (A.7) and (A.9), contain an averaged statistical, respectively, morphological information of the random field. Since the model-dependent expressions, equations (A.8) and (A.10), 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 A1 shows the averaged values of μ, ${{\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⁵ map sample at ${{N}_{\text{side}}}=128$, ${{\ell}_{\text{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 μ, ${{\sigma}_{0}}$ and ${{\sigma}_{1}}$. And the variation of the ratio ${{\sigma}_{1}}$/${{\sigma}_{0}}$ is noticeable (e.g. $34.162\,18/33.931\,90=1.006\,787$) and affects the magnitude of the Gaussian premise of the Minkowski Functionals v₁ and v₂.

Table A1. Table of μ and σ values.

(Units of μK) $\Lambda$CDM sample Full individual map range, no mask, 2°fwhm, bin 13 μK
${{\mu}_{\text{px}}}$	${{\sigma}_{0\text{px}}}$	${{\sigma}_{1\text{px}}}$	${{\sigma}_{1\text{px}}}/{{\sigma}_{0\text{px}}}$
$-3.4\times {{10}^{-7}}$	$59.533\,48$	$2020.083\,98$	$33.931\,90$

(Units of μK) $\Lambda$CDM sample Equal temperature range (ETR  ±  201 μK), U73 mask, 2°fwhm, bin 6μK
${{\mu}_{\text{px}}}$	${{\sigma}_{0\text{px}}}$	${{\sigma}_{1\text{px}}}$	${{\sigma}_{1\text{px}}}/{{\sigma}_{0\text{px}}}$
$-1.4534\times {{10}^{-3}}$	$59.132\,75$	$2020.103\,72$	$34.162\,18$

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⁻¹ Mpc⁻¹, the Hubble constant today

h  =  H₀/(100 km s⁻¹ Mpc⁻¹)  =0.6774, the normalized Hubble constant today

${{T}_{0}}=2.7255\pm 0.0006$ K, present-day CMB temperature

${{ \Omega }_{b}}{{h}^{2}}=0.022\,30\pm 0.000\,14$, baryon density

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

${{ \Omega }_{\nu}}{{h}^{2}}=0.002\,09$, neutrino density

${{ \Omega }_{k}}=0$, constant-curvature density parameter

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

${{N}_{\text{eff}}}=3.04\pm 0.33$, number of massless neutrinos

${\sum}^{}{{m}_{\nu}}<0.194$ eV, neutrino mass eigenstates

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

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

n_S  =  0.9667  ±  0.0040, scalar spectral index

x_e  =  1, ionization fraction

z_*  =  1089.90  ±  0.23, redshift of decoupling

${{ \Omega }_{ \Lambda }}=0.6911$  ±  0.0062, cosmological constant density

${{ \Omega }_{m}}=0.3089$  ±  0.0062, total matter density

Age  =  13.799  ±  0.021 Gyr, age of Universe

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 ℓ ≤ 2500⁸ are no longer Gaussian distributed, but the post-Born corrections would be weakly detectable in low noise CMB temperature maps merging the contributions of hundreds of ℓ–values between ℓ  = 1000 and ℓ  =  2500. Our present study uses an ℓ–range [2,256], well below ℓ  =  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.

From the power spectrum we generate with synfast a 10⁵ map ensemble. Given the map resolution (${{N}_{\text{side}}}=128$, ${{\ell}_{\text{max}}}=256$, and ${{2}^{\circ}}$fwhm), the number of 10⁵ maps allows a statistical stability in the sense of the second averaged MF v₁, which becomes smooth and stable around 10⁵ maps, a second different sample of 10⁵ maps gives very similar results and doubling the sample brings no visible improvement in the shape of v₁. A number of 10⁵ 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⁵ maps must be considered as a minimum requirement. We shall not develop more on these studies in the present work.

One observes that in the case of 10⁵ realizations there is an interval in τ, $\left[{{\tau}_{-}},{{\tau}_{+}}\right]$ such that $P(\tau )=0$ for all realizations, if $\tau \notin \left[{{\tau}_{-}},{{\tau}_{+}}\right]$. In order to simplify the problem, we use a symmetric interval $\left[-{{\tau}_{m}},+{{\tau}_{m}}\right]$, where ${{\tau}_{m}}=\text{max}\left\{|{{\tau}_{-}}|,{{\tau}_{+}}\right\}$. The interval $\left[-{{\tau}_{m}},+{{\tau}_{m}}\right]$ 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 'τ-lattice', ${{\tau}_{l}}=-{{\tau}_{m}}+(2l-1)(\Delta \tau /2)$, $l=1,2,\cdots,(2L+1)$, in such a way that the (L  +  1)th bin is centered at ${{\tau}_{L+1}}=0$. (Example: ${{\tau}_{m}}=396.5\;~\mu \text{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)$):

$$\begin{eqnarray}&&P\left({{\tau}_{l}}\right):=\frac{ \#\text{pixels}\left\{\delta T\in \left[{{\tau}_{l}}-\frac{ \Delta \tau}{2},{{\tau}_{l}}+\frac{ \Delta \tau}{2}\right)\right\}}{{{N}_{\text{tot}}} \Delta \tau}\;,\end{eqnarray} \tag{ A.11 }$$

with N_tot the total number of pixels, and for τ within the lth 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⁵ realizations. $\langle P(\tau )\rangle$ is still constant in a given bin, denoted by $\langle P\left({{\tau}_{l}}\right)\rangle$ in the lth bin and, thus, the derivative $\text{d}\langle P(\tau )\rangle /\text{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 \left[-{{\tau}_{m}},{{\tau}_{1}}+ \Delta \tau /2\right)$, and a last jump at a value ${{\tau}_{2L+2}}\in \left({{\tau}_{2L+1}}- \Delta \tau /2,{{\tau}_{m}}\right]$, 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

$$\begin{eqnarray}&&\langle \,f(\delta T)\rangle =\underset{l=1}{\overset{2L+1}{\sum}}\,\langle P\left({{\tau}_{l}}\right)\rangle {\int}_{{{\tau}_{l}}- \Delta \tau /2}^{{{\tau}_{l}}+ \Delta \tau /2}f(\tau )\,\text{d}\tau \;.\end{eqnarray} \tag{ A.12 }$$

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$):

$$\begin{eqnarray}&&{{\alpha}_{n}}=\frac{1}{n+1}\underset{l=1}{\overset{2L+1}{\sum}}\,\langle P\left({{\tau}_{l}}\right)\rangle \left[{{\left({{\tau}_{l}}+ \Delta \tau /2\right)}^{n+1}}-{{\left({{\tau}_{l}}- \Delta \tau /2\right)}^{n+1}}\right]\;.\end{eqnarray} \tag{ A.13 }$$

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

$$\begin{eqnarray}&&{{\alpha}_{0}}=\underset{l=1}{\overset{2L+1}{\sum}}\,\langle P\left({{\tau}_{l}}\right)\rangle \, \Delta \tau =1\quad \left(\text{normalization}\right)\;,\end{eqnarray} \tag{ A.14 }$$

$$\begin{eqnarray}&&{{\alpha}_{1}}=\mu =\underset{l=1}{\overset{2L+1}{\sum}}\,\,{{\tau}_{l}}\langle P\left({{\tau}_{l}}\right)\rangle \, \Delta \tau =0\;,\end{eqnarray} \tag{ A.15 }$$

whereas, for $n\geqslant 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$:

$$\begin{eqnarray}&&{{\alpha}_{n}}=\underset{l=1}{\overset{2L+1}{\sum}}\,\,\tau _{l}^{n}\langle P\left({{\tau}_{l}}\right)\rangle \, \Delta \tau +\frac{2}{n+1}\underset{k=1}{\overset{[n/2]}{\sum}}\,\left(\begin{array}{c} n+1 \\ 2k+1 \end{array}\right)\underset{l=1}{\overset{2L+1}{\sum}}\,\tau _{l}^{n-2k}\langle P\left({{\tau}_{l}}\right)\rangle \,{{\left(\frac{ \Delta \tau}{2}\right)}^{2k+1}}\;.\end{eqnarray} \tag{ A.16 }$$

(Example for n  =  2: $\sigma _{0}^{2}:={{\alpha}_{2}}={\sum}_{l=1}^{2L+1}\tau _{l}^{2}\langle P\left({{\tau}_{l}}\right)\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 (A.16), if it is not taken into account. For a discussion of the correction term in the case of the MFs v₁ and v₂, 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,

$$\begin{eqnarray}&&\alpha _{n}^{\prime}=\underset{l=1}{\overset{2L+1}{\sum}}\,\,\tau _{l}^{n}\langle P\left({{\tau}_{l}}\right)\rangle \, \Delta \tau,\end{eqnarray} \tag{ A.17 }$$

we obtain that $\sigma _{0}^{\prime}=\sqrt{\alpha _{2}^{\prime}}=59.6516~\mu \text{K}$ to be compared with ${{\sigma}_{0\text{px}}}=59.5335~\mu \text{K}$ (second equation in (A.7) and see also table A1).