A Numerical Investigation on the High-Frequency Geometry of Spherical Random Eigenfunctions

TL;DR

This paper reviews and numerically investigates the high-frequency geometry of spherical random eigenfunctions, focusing on Lipschitz-Killing curvatures and critical points, confirming analytic predictions and exploring statistical applications like CMB analysis.

Contribution

It provides the first comprehensive numerical study of geometric functionals of spherical eigenfunctions, validating analytic formulas and examining cancellation phenomena.

Findings

Accurate analytic predictions for expected values and variances of geometric functionals.

Confirmation of Berry's cancellation phenomenon at specific thresholds.

Potential applications in cosmic microwave background data analysis.

Abstract

A lot of attention has been drawn over the last few years by the investigation of the geometry of spherical random eigenfunctions (random spherical harmonics) in the high frequency regime, i.e ., for diverging eigenvalues. In this paper, we present a review of these results and we collect for the first time a comprehensive numerical investigation, focussing on particular on the behaviour of Lipschitz-Killing curvatures/Minkowski functionals (i.e., the area, the boundary length and the Euler-Poincar\'e characteristic of excursion sets) and on critical points. We show in particular that very accurate analytic predictions exist for their expected values and variances, for the correlation among these functionals, and for the cancellation that occurs for some specific thresholds (the variances becoming an order of magnitude smaller - the so-called Berry's cancellation phenomenon). Most of these functionals can be used for important statistical applications, for instance in connection to the analysis of Cosmic Microwave Background data.

Tables

Table 1. Table 1: Expected Values and Standard Deviation of the Euler-Poincaré characteristic. The theoretical expressions are given in Tables 6 and 7 in the Appendix, (the EPC corresponds to ℒ 0 ​ ( A u ​ ( f ℓ ; 𝕊 2 ) ) subscript ℒ 0 subscript 𝐴 𝑢 subscript 𝑓 ℓ superscript 𝕊 2 \mathcal{L}_{0}(A_{u}(f_{\ell};\mathbb{S}^{2})) ; the reported values are normalized dividing by 4 ​ π 4 𝜋 4\pi . Theoretical expressions at u = 0 𝑢 0 u=0 have yet to be determined.

Expected Values
	$\ell =100$		$\ell =300$		$\ell =500$		$\ell =700$		$\ell =900$
Threshold	Sim	Model	Sim	Model	Sim	Model	Sim	Model	Sim	Model
-3.0	-10.28	-10.53	-95.60	-95.38	-267.70	-264.88	-525.56	-519.01	-868.04	-857.79
-1.5	-155.50	-156.00	-1397.28	-1395.90	-3886.18	-3872.59	-7628.09	-7586.09	-12629.25	-12536.39
0.0	0.13	0.08	1.07	0.08	0.25	0.08	-1.13	0.08	1.16	0.08
1.5	155.62	156.16	1397.72	1396.05	3886.10	3872.75	7627.40	7586.25	12629.65	12536.55
3.0	10.44	10.69	95.96	95.54	268.15	265.04	525.77	519.17	868.19	857.95
Standard Deviation
	$\ell =100$		$\ell =300$		$\ell =500$		$\ell =700$		$\ell =900$
Threshold	Sim	Model	Sim	Model	Sim	Model	Sim	Model	Sim	Model
-3.0	4.23	4.26	21.63	22.05	46.80	47.39	77.76	78.46	119.51	114.36
-1.5	10.02	9.74	49.76	50.33	107.02	108.19	175.20	179.14	275.18	261.11
0.0	4.04	-	11.54	-	19.76	-	26.76	-	34.43	-
1.5	10.08	9.74	49.80	50.33	106.39	108.19	175.47	179.14	275.08	261.11
3.0	4.26	4.26	21.74	22.05	46.87	47.39	77.54	78.46	120.12	114.36

Table 2. Table 2: Expected Values and Standard Deviation of half of the boundary length functional, i.e. ℒ 1 ​ ( A u ​ ( f ℓ ; 𝕊 2 ) ) subscript ℒ 1 subscript 𝐴 𝑢 subscript 𝑓 ℓ superscript 𝕊 2 \mathcal{L}_{1}(A_{u}(f_{\ell};\mathbb{S}^{2})) ; the reported values are normalized dividing by 4 ​ π 4 𝜋 4\pi . Theoretical results are given in Tables 6 and 7 in the Appendix.

Expected Values
	$\ell =100$		$\ell =300$		$\ell =500$		$\ell =700$		$\ell =900$
Threshold	Sim	Model	Sim	Model	Sim	Model	Sim	Model	Sim	Model
-3.0	0.193	0.197	0.587	0.590	0.977	0.983	1.347	1.376	1.696	1.768
-1.5	5.736	5.768	17.225	17.246	28.667	28.724	39.965	40.202	51.077	51.681
0.0	17.764	17.766	53.101	53.121	88.417	88.477	123.727	123.832	159.052	159.187
1.5	5.732	5.768	17.232	17.246	28.667	28.724	39.961	40.202	51.079	51.681
3.0	0.193	0.197	0.589	0.590	0.977	0.983	1.347	1.376	1.695	1.768
Standard Deviation
	$\ell =100$		$\ell =300$		$\ell =500$		$\ell =700$		$\ell =900$
Threshold	Sim	Model	Sim	Model	Sim	Model	Sim	Model	Sim	Model
-3.0	0.088	0.089	0.150	0.153	0.193	0.198	0.226	0.234	0.266	0.265
-1.5	0.634	0.647	1.083	1.119	1.406	1.444	1.662	1.709	2.007	1.937
0.0	0.028	0.018	0.029	0.019	0.029	0.020	0.031	0.021	0.030	0.021
1.5	0.635	0.647	1.083	1.119	1.402	1.444	1.661	1.709	2.007	1.937
3.0	0.089	0.089	0.150	0.153	0.193	0.198	0.226	0.234	0.267	0.265

Table 3. Table 3: Expected Values and Standard Deviation of the Area functional. The Area is normalized to unity, i.e., divided by 4 ​ π 4 𝜋 4\pi . Theoretical expressions are summarized in Tables 6 and 7 in the Appendix.

Expected Values
	$\ell =100$		$\ell =300$		$\ell =500$		$\ell =700$		$\ell =900$
Threshold	Sim	Model	Sim	Model	Sim	Model	Sim	Model	Sim	Model
-3.0	0.9987	0.9987	0.9986	0.9987	0.9986	0.9987	0.9986	0.9987	0.9985	0.9987
-1.5	0.9336	0.9332	0.9330	0.9332	0.9325	0.9332	0.9320	0.9332	0.9315	0.9332
0.0	0.5000	0.5000	0.5000	0.5000	0.5000	0.5000	0.5000	0.5000	0.5000	0.5000
1.5	0.0663	0.0668	0.0671	0.0668	0.0675	0.0668	0.0680	0.0668	0.0685	0.0668
3.0	0.0013	0.0013	0.0014	0.0013	0.0014	0.0013	0.0014	0.0013	0.0015	0.0013
Standard Deviation
	$\ell =100$		$\ell =300$		$\ell =500$		$\ell =700$		$\ell =900$
Threshold	Sim	Model	Sim	Model	Sim	Model	Sim	Model	Sim	Model
-3.0	0.0007	0.0007	0.0004	0.0004	0.0003	0.0003	0.0003	0.0003	0.0002	0.0002
-1.5	0.0095	0.0097	0.0054	0.0056	0.0042	0.0043	0.0036	0.0037	0.0034	0.0032
0.0	0.0014	0.0013	0.0005	0.0004	0.0003	0.0003	0.0002	0.0002	0.0002	0.0001
1.5	0.0094	0.0097	0.0054	0.0056	0.0042	0.0043	0.0036	0.0037	0.0034	0.0032
3.0	0.0007	0.0007	0.0004	0.0004	0.0003	0.0003	0.0003	0.0003	0.0002	0.0002

Table 4. Table 4: Densities of the Expected Values of Critical Points, Extrema and Saddles. The simulation results are obtained for a bin width of Δ ​ u = 0.03 Δ 𝑢 0.03 \Delta u=0.03 . The theoretical results are given by the values of π 1 a ​ ( t ) superscript subscript 𝜋 1 𝑎 𝑡 \pi_{1}^{a}(t) defined in ( 32 )-(34).

	Critical		Extrema		Saddle
Threshold	Sim	Model	Sim	Model	Sim	Model
-3.0	0.0330	0.0318	0.0669	0.0635	0.0000	0.0000
-1.5	0.1648	0.1641	0.3101	0.3061	0.0232	0.0221
0.0	0.3471	0.3453	0.0000	0.0000	0.6855	0.6907
1.5	0.1642	0.1635	0.3057	0.3017	0.0264	0.0253
3.0	0.0332	0.0318	0.0672	0.0635	0.0000	0.0000

Table 5. Table 5: Total number ( u = − ∞ 𝑢 u=-\infty ) of critical points, extrema and saddles (Expected Values and their Standard Deviations). These mean values and standard deviations are computed from 1000 simulations. The theoretical expected value is 2 3 ​ ℓ ​ ( ℓ + 1 ) 2 3 ℓ ℓ 1 \frac{2}{\sqrt{3}}\ell(\ell+1) and the theoretical variance is 1 3 3 ​ π 2 ​ ℓ 2 ​ log ⁡ ℓ 1 superscript 3 3 superscript 𝜋 2 superscript ℓ 2 ℓ \frac{1}{3^{3}\pi^{2}}\ell^{2}\log\ell , see Table 8 in the Appendix.

	Critical		Extrema		Saddle
$\ell$	Mean	$\sigma$	Mean	$\sigma$	Mean	$\sigma$
100	11659.3300	52.7320	5830.6280	26.3011	5828.7020	26.4539
300	104306.4320	165.2952	52151.1920	82.6017	52155.2400	83.0788
500	289521.0780	275.0422	144729.9090	137.8649	144791.1690	139.0948
700	567436.9110	371.9949	283565.7720	187.7662	283871.1390	189.0990
900	937875.8670	479.1085	468449.3170	242.9301	469426.5500	247.7632

Table 6. Table 6: Theoretical expressions for the expected values and the variances of the three Lipschitz-Killing curvatures, Area ( ℒ 2 subscript ℒ 2 \mathcal{L}_{2} ), Half of the Boundary Length Area ( ℒ 1 subscript ℒ 1 \mathcal{L}_{1} ), Euler-Poincaré Characteristic ( ℒ 0 subscript ℒ 0 \mathcal{L}_{0} ). Threshold level u ≠ 0 𝑢 0 u\neq 0 .

LKC	Mean	Variance
${ℒ}_2\left(A_u(f_{\ell };{𝕊}^2)\right)$	$1-\Phi \left(u\right)$	$\frac{1}{8\pi }{u}^2{e}^{-u^2}\frac{1}{\ell }+O\left(\frac{\log \ell }{{\ell }^2}\right)$
${ℒ}_1\left(A_u(f_{\ell };{𝕊}^2)\right)$	$\frac{1}{4\sqrt{2}}{e}^{-u^2/2}\sqrt{\ell \left(\ell +1\right)}$	$\frac{1}{128}{u}^4{e}^{-u^2}\ell +O\left(\log \ell \right)$
${ℒ}_0\left(A_u(f_{\ell };{𝕊}^2)\right)$	$\frac{1}{4\pi }\sqrt{\frac{2}{\pi }}{e}^{-u^2/2}u\frac{\ell \left(\ell +1\right)}{2}+\frac{1}{2\pi }\left(1-\Phi \left(u\right)\right)$	$\frac{1}{128{\pi }^3}{u}^2{\left(u^2-1\right)}^2{e}^{-u^2}{\ell }^3+O\left({\ell }^2\log \ell \right)$

Table 7. Table 7: Theoretical expressions for the expected values and the variances of the three Lipschitz-Killing curvatures, Area ( ℒ 2 subscript ℒ 2 \mathcal{L}_{2} ), Half of the Boundary Length ( ℒ 1 subscript ℒ 1 \mathcal{L}_{1} ), Euler-Poincaré Characteristic ( ℒ 0 subscript ℒ 0 \mathcal{L}_{0} ). Threshold level u = 0 𝑢 0 u=0 .

LKC	Mean	Variance
${ℒ}_2\left(A_{u=0}(f_{\ell };{𝕊}^2)\right)$	$\frac{1}{2}$	$\frac{0.0188}{{\ell }^2}+o\left(\frac{1}{{\ell }^2}\right)$
${ℒ}_1\left(A_{u=0}(f_{\ell };{𝕊}^2)\right)$	$\frac{1}{4\sqrt{2}}\sqrt{\ell \left(\ell +1\right)}$	$\frac{1}{128}\frac{1}{16{\pi }^2}\left\{\log \ell +1.9542\right\}+O\left(1\right)$
${ℒ}_0\left(A_{u=0}(f_{\ell };{𝕊}^2)\right)$	$\frac{1}{4\pi }$	$O\left({\ell }^2\log \ell \right)$

Table 8. Table 8: Theoretical expressions for the expected values and the variances of the critical points. Threshold level u ≠ 0 𝑢 0 u\neq 0 and u = 0 𝑢 0 u=0 .

	Mean	Variance
${𝒩}_{-\infty }^c\left(f_{\ell }\right)$	$\frac{2}{\sqrt{3}}{\ell }^2+O\left(1\right)$	$\frac{1}{3^3{\pi }^2}{\ell }^2\log \ell +O\left({\ell }^2\right)$
${𝒩}_u^c\left(f_{\ell }\right)$	$\frac{2}{\sqrt{8\pi }}\ell \left(\ell +1\right)\left\{\frac{2}{\sqrt{6}}\Gamma (\frac{1}{2},\frac{3u^2}{2})+e^{-u^2/2}u\right\}$	${\ell }^3\frac{1}{8\pi }{e}^{-3u^2}{u}^2{\left(2+e^{u^2}\left(u^2-1\right)\right)}^2+O\left({\ell }^2\log \ell \right)$
${𝒩}_0^c\left(f_{\ell }\right)$	$\frac{1}{\sqrt{3}}\ell \left(\ell +1\right)+O\left(1\right)$	$\frac{1}{4{\pi }^{2}27}{\ell }^2\log \ell +O\left({\ell }^2\right)$

Equations

$$A_u(f):=\left\{x\in {\mathbb{S}}^2:f(x)\ge u\right\}\text{ .}$$

$${\Delta }_{{\mathbb{S}}^2}{f}_{\ell }+{\lambda }_{\ell }{f}_{\ell }=0,\text{mbox}\text{mbox}\text{mbox}{f}_{\ell }:{\mathbb{S}}^2\to \mathbb{R},$$

$$\dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}\bigg{\{}\sin\theta\dfrac{\partial}{\partial\theta}\bigg{\}}+\dfrac{1}{\sin^{2}\theta}\dfrac{\partial^{2}}{\partial\varphi^{2}},\mbox{ }0\leq\theta\leq\pi,\mbox{ }0\leq\varphi<2\pi,$$

$$\mathbb{E}[f_{\ell }(x)f_{\ell }(y)]=P_{\ell }(\cos d(x,y))$$

$$d(x,y)=\arccos (⟨x,y⟩).$$

$$f(x)=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{a}_{\ell m}{Y}_{\ell m}(x)=\sum _{\ell =0}^{\infty }\sqrt{\frac{(2\ell +1)C_{\ell }}{4\pi }}{f}_{\ell }(x)\text{ ;}$$

$$a_{\ell m}={\int }_{{\mathbb{S}}^2}f(x){\bar{Y}}_{\ell m}(x)dx,\ell =1,2,....,m=-\ell ,...,\ell ,$$

$${\rho }_l(u)=(2\pi {)}^{-(l+1)/2}{H}_{l-1}(u)e^{-u^2/2}\text{ ,}$$

$$H_0(u)=1,H_1(u)=u,H_2(u)=u^2-1,\dots ;$$

$$H_{-1}(u)=\sqrt{2\pi }(1-\Phi (u))e^{u^2/2},$$

$${\rho }_0(u)$$

$${\rho }_1(u)$$

$$\left[\begin{array}[]{c}i+l\\ l\end{array}\right]=\left(\begin{array}[]{c}i+l\\ l\end{array}\right)\frac{\omega_{i+l}}{\omega_{i}\omega_{l}}\text{ , for }\omega_{i}=\frac{\pi^{i/2}}{\Gamma(\frac{i}{2}+1)}\text{ ,}$$

$$\lambda^{i/2}\mathbb{E}\mathcal{L}_{i}(A_{u}(f(x);D))=\sum_{l=0}^{\dim(D)-i}\left[\begin{array}[]{c}i+l\\ l\end{array}\right]\lambda^{(i+l)/2}\rho_{l}(u)\mathcal{L}_{i+l}(D)\text{ .}$$

$$\mathbb{E}{\mathcal{L}}_0(A_u(f_{\ell }(.);{\mathbb{S}}^2))=2\left\{1-\Phi (u)\right\}+\frac{{\lambda }_{\ell }}{2}\frac{u{e}^{-u^2/2}}{\sqrt{(2\pi {)}^3}}4\pi \text{ ;}$$

$$\mathbb{E}{\mathcal{L}}_1(A_u(f_{\ell }(.);{\mathbb{S}}^2))=\frac{\pi }{2}\frac{1}{\sqrt{2}}{\lambda }_{\ell }^{1/2}\frac{e^{-u^2/2}}{2\pi }4\pi =\frac{\pi }{\sqrt{2}}{\lambda }_{\ell }^{1/2}{e}^{-u^2/2}\text{ ;}$$

$$\mathbb{E}{\mathcal{L}}_2(A_u(f_{\ell }(.);{\mathbb{S}}^2))=4\pi \times \left\{1-\Phi (u)\right\}\text{ .}$$

$$\mathtt{Proj}[{\mathcal{L}}_k(A_u(f_{\ell };{\mathbb{S}}^2))|2]$$

$$=\frac{1}{2}\left[\begin{array}[]{c}2\\ k\end{array}\right]\left\{\frac{\lambda_{\ell}}{2}\right\}^{(2-k)/2}H_{1}(u)H_{2-k}(u)\phi(u)\frac{1}{(2\pi)^{(2-k)/2}}\int_{\mathbb{S}^{2}}H_{2}(f_{\ell}(x))dx+a_{k}(\ell),$$

$$a_{k}(\ell)=\left\{\begin{array}[]{cc}O_{p}(\ell)&\text{for }k=0,\\ 0&\text{for }k=1,2\end{array}\right.\text{ .}$$

$$\frac{{\lambda }_{\ell }}{2}{\int }_{{\mathbb{S}}^2}{H}_2(f_{\ell }(x))dx$$

$${\mathcal{L}}_2({\mathbb{S}}^2)=\frac{{\lambda }_{\ell }}{2}\times 4\pi =\frac{{\lambda }_{\ell }}{2}{\int }_{{\mathbb{S}}^2}{H}_0(f_{\ell }(x))dx\text{ .}$$

$$\mathtt{Proj}[{\mathcal{L}}_k(A_u(f_{\ell };{\mathbb{S}}^2))|0]$$

$$=\left[\begin{array}[]{c}2\\ k\end{array}\right]\left\{\frac{\lambda_{\ell}}{2}\right\}^{(2-k)/2}H_{1-k}(u)\phi(u)\frac{1}{(2\pi)^{(2-k)/2}}\int_{\mathbb{S}^{2}}H_{0}(f_{\ell}(x))dx+b_{k}(\ell)\text{ ,}$$

$$b_{k}(\ell)=\left\{\begin{array}[]{cc}2(1-\Phi(u))=O(1)&\text{for }k=0,\\ 0&\text{for }k=1,2\end{array}\right..$$

$$\mathtt{Proj}[{\mathcal{L}}_2(A_u(f_{\ell };{\mathbb{S}}^2))|0]$$

$$\mathtt{Proj}[{\mathcal{L}}_2(A_u(f_{\ell };{\mathbb{S}}^2))|2]=\frac{1}{2}{\left\{\frac{{\lambda }_{\ell }}{2}\right\}}^0\left[H_0(u)H_1(u)\varphi (u)\right]{\int }_{{\mathbb{S}}^2}{H}_2(f_{\ell }(x))dx\text{ ;}$$

$$\text{Var}\left({\mathcal{L}}_2(A_u(f_{\ell };{\mathbb{S}}^2))\right)$$

$$=$$

$$\text{Var}\left({\int }_{{\mathbb{S}}^2}{H}_2(f_{\ell }(x))dx\right)=16{\pi }^2\frac{2}{2\ell +1}.$$

$$D_{\ell }=2{\mathcal{L}}_2(A_{u=0}(f_{\ell };{\mathbb{S}}^2))-4\pi ;$$

Full text

**Keywords and Phrases: ***CMB, Data Analysis, Minkowski Functionals, Gaussian Kinematic Formula, Spherical Harmonics

***AMS Classification: ***60G60, 62M15, 53C65, 42C10, 33C55

*PACS Numbers: 98.80.Es, 95.75.Mn, 95.75.Pq, 02.50.-r

1 Introduction

The geometry of the excursion sets for random fields on the sphere (to be defined below) has been the object of rather intense research over the last decade or so. It is well-known that these geometrical properties can be characterized in terms of the so-called Lipschitz-Killing curvatures (or equivalently, Minkowski functionals), which in the two dimensional case correspond to the excursion area, (half) the boundary length and the Euler-Poincaré characteristic (connected regions minus number of holes) of the excursion set. A comprehensive description of Lipschitz-Killing curvatures for excursion sets of random fields is given in the excellent monograph by [1]; these functionals can be computed on real data by means of accurate and numerically efficient algorithms [16, 15, 13]. They have been widely applied in the analysis of experimental data, especially in a Cosmological framework, see i.e., [29, 27, 11, 38, 28, 35, 49] and the references therein.

A lot of mathematical efforts have been spent since the ’80s on the characterization of expected values of these functionals under Gaussianity, culminating in the discovery of the beautiful Gaussian Kinematic Formula ([1]); comparing these expected values with realizations allows the implementation of a number of tests for Gaussianity and Isotropy (see again [37, 49]).

While the behaviour of expected values is now fully understood, it is clear that the implementation of more sophisticated, hence more sensitive, testing procedures requires further knowledge, in particular the variance of these functionals and therefore the possibility to establish Central Limit Theorems with correct normalization factors. Establishing a Central Limit result requires of course the exploitation of a suitable notion of asymptotic behaviour; in the framework of spherical fields, the only relevant notion seems to be the one of High-Frequency asymptotics. In particular, it is well-known that isotropic random fields on the sphere can be decomposed by means of the Spectral Representation Theorem into the sum of orthogonal components, each of them corresponding to a different multipole $\ell$. The behaviour of geometric functionals in the high-frequency/high energy limit for these components has been studied by several authors in recent years, starting from [30, 31] for the number of connected components, [48] for the nodal length, and then including, among others, [25] and [21] for the excursion area, [26], [24] and [40] for the Defect, [4] for the Euler-Poincaré characteristic, [48], [22], [20], [41] for the distribution of the nodal length, [7] for the critical values and [8] for the total number of critical points (see also [17], [43], [10, 19, 42, 34] for related works covering also the 2-dimensional torus, [18] for the 3-dimensional torus and [33] for planar random waves).

Our aim in this paper is to present a unified overview of the literature, and especially to perform a detailed numerical investigation to verify the practical relevance of these results when investigating spherical Gaussian maps. We shall address several issues concerning not only the expected value and variances of Minkowski functionals, but also their cross-correlation across different level sets. The theoretical predictions which have so far been produced are validated for the first time from a numerical point of view, and moreover their domain of applicability is clarified. Indeed, in terms of the variances the theoretical expressions which are obtained should be viewed as leading terms in series expansions of the variances over different “chaos" components; as such, the approximation depends on the rate of convergence to zeroes of the terms which are dropped. These rates are known to be polynomial in some cases (namely, those corresponding to non-zero levels) and logarithmic in others (those corresponding to zero levels); this duality is mirrored in the numerics that we shall present below.

Minkowski functionals are not the only objects of interest in this paper. Indeed, some other recent contributions have derived neat analytic formulae for the expected number and the variance of critical points on the same spherical harmonics components as considered earlier for Minkowski functionals. We are hence providing for the first time numerical evidence also on these statistics.

The plan of the paper is as follows: in Section 2, we review the results on the expected values and variances for the Lipschitz-Killing curvatures; in Section 3 we discuss the behaviour of critical points, again reviewing the analytic results that are currently available, while Section 4 is devoted to the analysis of the correlation among these different functionals. Section 5 describes our implementation algorithms and presents the numerical results; we then draw some conclusions and present directions for future work.

2 Characterization of Excursion Sets for Random Spherical Harmonics

In the case of the two-dimensional sphere, the excursion sets $A_{u}(f)$ of a given (possibly random) function $f$ are defined for any real number $u$ as

[TABLE]

Of course, in the limit where we take $u=-\infty,$ we have that $A_{u}(f)=\mathbb{S}^{2}$. In this paper, we shall be concerned with random eigenfunctions $f_{\ell}$ which satisfy the Helmhotz equation

[TABLE]

where $\Delta_{\mathbb{S}^{2}}$ is the Laplace-Beltrami operator on $\mathbb{S}^{2}$, defined as usual as

[TABLE]

and $\lambda_{\ell}=\ell(\ell+1),\mbox{ }\ell=0,1,\dots$. For a given eigenvalue $\lambda_{\ell},$ the corresponding eigenspace is the $(2\ell+1)-$dimensional space of spherical harmonics of degree $\ell.$ The random fields $\{f_{\ell}(x),x\in\mathbb{S}^{2}\}$ are Gaussian and isotropic with zero mean $\mathbb{E}[f_{\ell}(x)]=0$ and variance $\mathbb{E}[f_{\ell}(x)^{2}]=1.$ The covariance function is given by

[TABLE]

where $P_{\ell}$ are the Legendre polynomials and $d(x,y)$ is the spherical geodesic distance between $x$ and $y,$ i.e.

[TABLE]

Spherical random eigenfunctions are of interest because they can also be interpreted as the harmonics/Fourier components of data observed on the sphere. Indeed, let us first recall the well-known Spectral Representation Theorem for spherical random fields, which states that the following identity holds, in the $L^{2}$ sense:

[TABLE]

here, the sequence $\{C_{\ell}\}$ denotes the so-called angular power spectrum of the field. The spherical harmonics coefficients may be computed from the field $f(.)$ by means of the inverse transform

[TABLE]

with $\mathbb{E}[a_{\ell m}]=0$ and $\mathbb{E}|a_{\ell m}|^{2}=C_{\ell}$. The inverse transform (3) is only feasible for unmasked (full-sky) data, a condition which is usually considered rather difficult to meet for astrophysical experiments such as those concerning Cosmic Microwave Background radiation. Rather recently, however, full-sky maps were produced for instance by [3] and by the Planck collaboration in its 2018 release (see [36]).

Let us now recall again the definitions of the Lipschitz-Killing Curvatures (LKCs), which correspond to Minkowski functionals up to a different indexing and normalization factors; in two dimension, they are given by (a) the Euler-Poincaré characteristic (written $\mathcal{L}_{0}(A_{u}(f))$), e.g. the number of connected regions minus the number of holes; (b) half the boundary length of the excursion regions (written $\mathcal{L}_{1}(A_{u}(f))$); the area of the excursion regions (written $\mathcal{L}_{2}(A_{u}(f))$), which corresponds to the first Minkowski functional. The expected values of these functionals when evaluated on the excursion sets of Gaussian fields have been fully characterized by the Gaussian Kinematic Formula (GKF), see [1].

We now need the family of functions $\rho_{l}(u)$, for $l\in\mathbb{N}$, defined as

[TABLE]

where $H_{k}(u)$, $k\in\mathbb{N}$, denotes as usual the family of Hermite polynomials, that is,

[TABLE]

it is convenient to define also

[TABLE]

where $\Phi(u)$ is the Gaussian cumulative distribution function, whence

[TABLE]

[1] write these components as $\mathcal{M}_{k}([u,\infty))=\frac{1}{\sqrt{2\pi}}H_{k}(u)e^{-u^{2}/2}$ and denote them Gaussian Minkowski functionals. The so-called “flag" coefficients are instead given by

[TABLE]

that is, $\omega_{i}$ represents the area of the $i-$dimensional unit ball, $\omega_{1}=2,$ $\omega_{2}=\pi,$ $\omega_{3}=\frac{4}{3}\pi$ and $\Gamma(\cdot)$ being the Gamma function $\Gamma(n+1)=n\Gamma(n)$. As a last ingredient, we write $\lambda$ for the parameter which represents the second derivative of the covariance function at the origin.

We are now ready to present the general expression for the expected value of Lipschitz-Killing curvatures of a process $f$ on a manifold $D$, i.e., the Gaussian Kinematic Formula which reads (Theorem 13.2.1 in [1]):

[TABLE]

As an application of the previous result, let us consider the Fourier components $\{f_{\ell}(\cdot)\}_{\ell=1,2,\dots}$ normalized to have variance one; the GKF yields immediately (compare [23], Corollary 5, see also [9])

[TABLE]

[TABLE]

and

[TABLE]

Of course, in order to exploit Lipschitz-Killing curvatures/Minkowski functionals to implement data analysis tools the expected value by itself is not sufficient, but we need also analytic expression for the variance. The latter was derived in some recent results by [7, 4]; see [45] for a review.

For our purposes, the results in these papers can be summarized as follows; the asymptotic behaviour of each of the three Lipschitz-Killing curvatures, evaluated on the excursion sets of random spherical harmonics, is dominated by a single, fully degenerate component, which can be written as:

[TABLE]

[TABLE]

where

[TABLE]

Here, and in the sequel, we use $\mathtt{Proj}[.|q]$ for the projection of random quantities on the so-called Wiener chaoses of order $q$; the latter are spaces generated by linear combinations of Hermite polynomials of order $q$, computed in $f_{\ell}$ and its derivatives (we refer to [32, 20], [4] and the references therein for more discussions and details). It is also important to notice that $\frac{\lambda_{\ell}}{2}=P_{\ell}^{\prime}(1)$ represents the derivative of the covariance function of random spherical harmonics at the origin, so that the term

[TABLE]

can be viewed as a (random) measure of the sphere induced by the Riemannian metric, somewhat in analogy with the interpretation given for the Gaussian Kinematic Formula on the expected value in the book by [1]; recall indeed that for eigenfunctions $f_{\ell}$ on the sphere $\mathbb{S}^{2}$ the term ${\mathcal{L}}_{2}(\mathbb{S}^{2})$ which appears in (10) is exactly given by the area of the sphere with radius $\left\{\frac{\lambda_{\ell}}{2}\right\}^{1/2},$ i.e.,

[TABLE]

As was noted in [4], the Gaussian Kinematic Formula can be rewritten with a very similar expression to (14), i.e.:

[TABLE]

[TABLE]

where

[TABLE]

More explicitly (see also [25, 26, 21], [39]), we have the following analytic expressions for the leading term components of the LKCs (expected values and dominant stochastic term):

a) Excursion Area ($k=2$) As explained above, the expected value for the excursion area can be obtained (as for the other Lipschitz-Killing curvatures) by a simple application of the Gaussian Kinematic Formula, which yields:

[TABLE]

the leading term in the fluctuations is provided by (see [25, 21])

[TABLE]

with an asymptotic variance which is given by

[TABLE]

where we have used the fact that (see [25, 26, 39])

[TABLE]

Analogous results, although with different constants, can be established on subdomains of the sphere (see [47]). For $u=0$, we obtain a quantity equivalent to the so-called defect (see [26]) i.e.,

[TABLE]

the expected value is immediately seen to be zero, while it can be shown that the variance is given by

[TABLE]

where the constant $C$ can be computed as

[TABLE]

(see equation (25), [26]), and

[TABLE]

with

[TABLE]

being the $J_{0}$ Bessel function. In the Appendix, we perform a numerical investigation on the value of the constant $C$; more precisely, to obtain a precision of $1.0\times 10^{-4}$, it is sufficient to sum the terms in (20) until $q=20$, obtaining the value

[TABLE]

The constants $C_{q}$ are obtained by numerical integration, whereas for $C_{3}$ the exact value is computed in [25] and it is given by

[TABLE]

Remark 2.1**.**

It is easily seen that $50\%$ of the contribution of the sum in (23) comes from the first term, which is $0.0613$. Moreover, the sum of the first and second term is $0.0860$ (see Appendix), and thus, $80\%$ of the variance for the defect is explained by the third and fifth chaoses alone.

Summing up, for $u=0$ the leading term in equation (16) disappears and we have to use the higher order approximation to find that

[TABLE]

In all the above equations, normalizing the area by $4\pi$ divides out the $16\pi^{2}$ term; this is the normalization that we shall adopt in the tables to follow in Section 5 (Tables 1, 2, 3 and Figures 1, 2).

b) (Half) The Boundary Length ($k=1$) Let us now consider the boundary length of excursion regions. To compute the expected value, it is enough to exploit the Gaussian Kinematic Formula; as before, note that we shall normalize by $4\pi$ in the simulations (see Table 2) so that we obtain

[TABLE]

Likewise, using results in [39], [48], [22], we have for the leading stochastic term

[TABLE]

and using again (18) the variance can easily be seen to be

[TABLE]

Again, in the simulations below (see Table 2), normalizing the boundary length by $4\pi$ divides by a factor $16\pi^{2}$, leading (up to negligible terms) to a variance of order $\frac{\ell}{128}u^{4}e^{-u^{2}}$.

For $u=0$ the leading term in the previous expression disappears (the so-called Berry’s cancellation phenomenon, see [2], [48]) and the variance is of smaller order; more precisely, we have that ([48])

[TABLE]

(the same happens for shrinking subdomains of the sphere, see [46]). It is important to notice that the difference between the leading and remainder terms is here only of logarithmic order, and we hence expect a less precise approximation (in relative terms) in the simulations. On the other hand, it should also be noted that the variances at stake are much smaller than for $u\neq 0$, and then the absolute error in the simulations will turn out to be particularly small.

Hence, when $u=0$, the leading term in (25) disappears and the nodal length is asymptotic to the sample trispectrum, namely

[TABLE]

where $h_{\ell,4}=\int_{\mathbb{S}^{2}}H_{4}(T_{\ell}(x))\,dx$, which is logarithmic and hence we derive (27). To be clear, as given in [22],

[TABLE]

and, setting $L:=\ell+\frac{1}{2}$, since

[TABLE]

we compute the last integral numerically, exploiting Matlab. We report some values in the table below.

[TABLE]

More explicitly, it was shown in [26] that

[TABLE]

to find a better approximation, we evaluate numerically the constant

[TABLE]

see the Appendix for some analytic results. Hence, we conclude that, up to smaller order terms

[TABLE]

Then, the variance of the scaled sample trispectrum $\mathcal{M}_{\ell}$ is asymptotically given by

[TABLE]

Finally, let us recall that these results, as in [48], [39], [22] and [46] refer to the boundary length, not to the first Lipschitz-Killing curvature; there is hence a difference of a factor 2 in the expected value, and a factor 4 in the variance. The values in the Table 2 refer to the Lipschitz-Killing curvature, hence they have been normalized accordingly.

c) Euler-Poincaré characteristic ($k=0$) The Euler-Poincaré characteristic (EPC) for random spherical harmonics was investigated by [6], [4] among others, where the following expressions are given for the expected value and the second chaotic component:

[TABLE]

[TABLE]

All the EPC equations are normalized by $4\pi$ in the simulations, hence the $16\pi^{2}$ term is divided out.

Given these results, [4] showed that the variances of LKCs are dominated by the variance of the second order Wiener chaos; indeed, for the Euler-Poincaré characteristic the expected value and variance are given, for an interval $I\subset\mathbb{R}$, by

[TABLE]

[TABLE]

and in particular for semi-intervals of the form $I=[u,\infty)$ one obtains

[TABLE]

[TABLE]

Note that, after normalizing the Lipschitz-Killing curvatures by their expected value, their relative variances converge to zero as the frequency increases, so that relative fluctuations become negligible on small scales (Tables 1, 2, 3).

3 Characterization of Critical Points for Random Spherical Harmonics

As a further tool of investigation, we shall consider in this paper also the behaviour of critical points for random spherical harmonics, which has recently been fully characterized by [7, 8, 5], among others.

More precisely, by definition critical points, extrema and saddles are, respectively, given by:

[TABLE]

[TABLE]

[TABLE]

where we used $a=c,e,s$ to label critical points, extrema and saddles respectively.

We now recall the following results on the expectations and variances:

For every interval $u\in\mathbb{R}$, we have, as $\ell\rightarrow\infty$,

[TABLE]

where $a=c,e,s$ and for the density functions

[TABLE]

Similarly, for every $u\in\mathbb{R}$, as $\ell\rightarrow\infty$,

[TABLE]

where,

[TABLE]

The leading constants for the variances can be written more explicitly as

[TABLE]

Note that also in this case, the second component is the leading term of the expansion and it is important to stress how the leading terms in the variances cancel in all cases at the threshold $u=-\infty$; in other words, the variance is smaller when we focus on the total number of critical points (see [8]). This is again a form of the so-called “Berry’s cancellation phenomenon", which we have also discussed earlier for the Lipschitz-Killing curvatures. Indeed, the behaviour of critical points and saddles can be shown to be dominated by the second order chaotic component, which takes the form (see [4])

[TABLE]

and similarly for saddles.

Because this second-order chaos component (and hence the leading term in the variance) vanishes at $u=-\infty,0$, the next component becomes of interest; it can be shown that this term is proportional to the fourth-order chaos, and indeed, for the total number of critical points, it holds that (see [8])

[TABLE]

moreover, it is also possible to consider separately extrema (minima and maxima) and saddles, yielding

[TABLE]

and

[TABLE]

4 On Correlations

The results presented in the previous sections can be summarized as follows:

1. For general threshold $u\neq 0$, the fluctuations around the proper expected values for the area, the boundary length and the Euler-Poincaré characteristic of excursion regions is dominated by a single stochastic term, which is proportional to the so-called second order Wiener chaos; namely $h_{\ell;2}=\int_{\mathbb{S}^{2}}H_{2}(f_{\ell}(x))\,dx$.

2. At $u=0$, this term is disappearing; the boundary length is then dominated by the fourth-order chaos, i.e., a single term which is proportional to $h_{\ell;4}=\int_{S^{2}}H_{4}(f_{\ell}(x))dx$. For the excursion area and the Euler-Poincaré, this term is disappearing as well and lower order terms are dominant.

3. Likewise, the critical points above general threshold levels $u$ are dominated by a single term, proportional to $h_{\ell;2}$; this term disappears for $u=-\infty$, where the total number of critical points is dominated by a single term proportional to $h_{\ell;4}$.

Note also that the variance of $h_{\ell;2}$ is of order $O(\frac{1}{\ell})$, the variance of $h_{\ell;4}$ is of order $O(\frac{\log\ell}{\ell^{2}})$, and the variance of all other chaoses (for $q=3,5,6,7...$) is of order $O(\frac{1}{\ell^{2}}$). As a consequence, we expect almost perfect correlation for all statistics which are dominated by $h_{\ell;2}$; some correlation (but not too strong, given the logarithmic rate) for statistics dominated by $h_{\ell;4}$; no correlation for statistics which are dominated by chaoses of different order. These conjectures are indeed very well confirmed by the numerical evidence that we shall present in the Section below (Fig. 5).

5 Numerical results

In this section we describe the comparison of the analytical results outlined in the previous sections to the corresponding results from simulations. In order to implement this comparison, we generated 1000 Gaussian realizations of random spherical eigenfunctions/spherical harmonics for different values of the multipoles $\ell$, ranging from $\ell=100$ to $\ell=900$. These values for the multipoles $\ell$ are representative of the resolution which can be currently achieved by satellite experiments such as Planck (see [36]); for instance, these eigenfunctions could be taken to be the spherical Fourier component of a simulated CMB map, according to a standard routine provided by the HEALpix [14] package. The simulations algorithms are described more fully in the subsection to follow.

Simulations and Algorithm

We used the HEALpix synfast routine to simulate a Gaussian realization map starting from a given power-spectrum. In practice, we used the so-called best-fit Planck power spectrum to generate the maps, and then we extracted the multipoles to focus on, normalizing their variance to unity. Of course, our results are independent from the choice of the input power spectrum, and indeed it would be possible to generate directly the single eigenfunctions at a given multipole.

A single multipole map $f_{\ell}(x)$ is obtained by using the HEALpix alm2map routine, after having extracted the proper subset of coefficients $a_{\ell m}$. In all cases the map resolution parameter $N_{side}$ is set to twice the value of the corresponding multipole. As mentioned earlier each map is normalized to have unit variance.

It is very important to notice that each functional is normalized “per unit area", i.e., all the reported values have been standardized dividing by $4\pi$. Both the expected values and the variances are affected in the obvious way.

We compute the three Minkowski functionals, which are equivalent to the LKCs up to constant factors, and critical point counts from these normalized multipole maps. In short, the area, i.e. the first Minkowski functional, is simply computed by evaluating the number of pixels above a certain threshold. The perimeter length, the second Minkowski functional, is computed by tracing isocontour lines in pixel space. For a sufficiently high-resolution map, pixels around isocontour lines have different signs relative to the contour line, after normalizing the lines to zero. To measure the length of these lines, sets of four pixels are compared; when at least two of them have different signs, the locations where the contour line enters and exits these sets of pixels are determined and the length is iteratively calculated by standard dot product. For the Euler-Poincaré characteristic, the third Minkowski functional, we used the Fortran implementation of the algorithms described in Appendix G of [13] (see also [12]). This algorithm is based on the Gauss-Bonnet theorem - where the Euler characteristic of a region is obtained by integrating the curvature over the boundary surface. Given we are working on a pixelized surface, the surface curvature of an excursion region can be thought of as concentrated in the corners of the pixels that are at the boundary between the pixels above and below the threshold. This is true as any continuous deformation of the region conserve the topology. What is needed is, therefore, to devise a strategy that assigns appropriate curvature weights for each boundary grid vertex - Appendix G of [13] explains in more detail the strategy used in the Fortran code. Once the weights are assigned, the sum of the weights over all the vertices gives us the Euler characteristic of the excursion set.

Our detailed investigation using different algorithms to compute the Euler-Poincaré characteristic showed that for a map defined at a given $N_{\rm side}$, the maximum multipole for which a percent numerical accuracy can be obtained is $\ell_{max}\sim N_{\rm side}/3$. While it would be possible to cover larger values, we do not believe this is essential for our purpose in this paper.

Results

We first proceed to report in Tables 1, 2, 3, a numerical comparison of the expected values computed on simulated maps and their analytical predictions, and likewise, Monte Carlo estimates of root mean squared errors and their analytical predictions. We stress that the fit is truly remarkable: the percentage errors are smaller than 1% for most non-zero values of the threshold parameter. It should be recalled here that the analytical predictions for expected values are exact, while for the variances we are only giving the leading term in a series of positive addends; for non-zero values of $u$, the neglected terms in the variance (as mentioned earlier) are a factor $\ell$ smaller than the leading one, as mirrored in the simulations.

The analytic approximation for the variances in the case $u=0$ is slightly worse, in relative terms, but actually even better, in absolute terms. This was explained earlier in Section 2; in short, the variances at $u=0$ are an order of magnitude smaller than at other thresholds, because the leading term cancels, and new elements become dominant (the fourth-order chaos, in the case of the nodal length). Thus, focussing for instance on the boundary length, here the dominant term is larger than the neglected ones only by a logarithmic factor; as a consequence, variances tend to be underestimated (a similar phenomenon occurs for the total number of critical points, see below). In absolute terms, the discrepancy between simulations and analytic results for the nodal length is in the order of $10^{-3}/10^{-1}$, to be compared with expected values in the order of $10/10^{2}$, so that the relative error is in the order of $10^{-3}$.

The results for critical points (Tables 4, 5 and Figures 3, 4) are, in our view, equally impressive, with relative errors in the order of $10^{-3}/10^{-4}$, and absolute ones in the order of $10/10^{2}$, to be compared with expected values that run in the hundreds of thousands ($10^{4},10^{5}$).

To help visualization, we produced some plots that compare the analytic predictions with the realizations; more precisely, in Figure 1 we compare the multipole space analytical results (red curve) given in Section 2 with that of the simulations (black curve - mean of the simulations). The $68\%,95\%$ and $99\%$ Confidence Intervals are shown from dark to light grey bounds. From the top to the bottom rows, the figures show the plots of the results corresponding to multipoles $\ell=500,700,900$. We stress that our fit is extremely accurate even at low multipole values; we also note the improved concentration around the expected values at higher-multipoles.

In figure 5 we present our evidence on cross-correlations. As expected, correlations are very close to one (in absolute value) for any pair of random statistics evaluated at non-zero thresholds, including area, boundary length, Euler-Poincaré characteristic and the number of critical points; considering extrema (maxima and minima) and saddles separately would yield the same outcome. The simulations also confirm uncorrelation when expected, for instance between the nodal length (which is dominated by the fourth-order chaos, see [22]) and the defect, which is dominated by odd order chaoses (see [26]).

All these results have potential for applications in the statistical analysis of random fields, for instance when testing for nonGaussianity and isotropy or to search point-like sources/impurities in Cosmic Microwave Background radiation data. We do not address these issues in the present work, but we leave them as avenues for further research.

6 Appendix

In this Appendix, we report some numerical computations on constants needed for higher-order approximations on the behaviour of the limiting variances.

Recall first that, denoting, as usual, $h_{\ell,q}:=\int_{\mathbb{S}^{2}}H_{q}(T_{\ell}(x))\,dx$, it is known (see for example [26]), that for $q=3$ and $q\geq 5$, one has

[TABLE]

with

[TABLE]

being the $J_{0}$ Bessel function. Moreover for $q=2,4,$ the order of magnitude of the corresponding variance is larger, namely:

[TABLE]

[TABLE]

Defining, as in (21),

[TABLE]

for the values $L=50$, 100, 200, we find, exploiting Matlab, the following numerical evaluations.

[TABLE]

It can be seen from Figure 6 and Figure 7 (realized for $L=100$) that the behavior of $C_{q}$, for $q\geq 5$ (odd or even), is well approximated by

[TABLE]

Remark 6.1**.**

In [24], the asymptotic behavior of the coefficients $a_{k}$ is proved to be (using Stirling approximation) $\frac{a}{k^{3/2}}$, where the constant $a$ can be computed to be $\frac{1}{2\sqrt{\pi}}$; in view of (39) the product $a_{k}C_{2k+1}$ behaves as $\frac{1}{2\sqrt{\pi}k^{3/2}}\times\frac{2}{2k+1}$ and therefore as $\frac{1}{2\sqrt{\pi}k^{5/2}}$. Indeed, figure 8 compares the points $a_{k}C_{2k+1}$ with the function $\frac{1}{2\sqrt{\pi}q^{5/2}}$; the fit appears very good, for reasonably large values of $q$.

Let us now try to improve the numerical approximation for the variance of the fourth-order chaos. We have shown numerically that

[TABLE]

Actually, to check the validity of this result we can estimate the difference semi-analytically. Indeed, splitting the domain of the integral in $[0,10]$ and $[10,L]$, we obtain

[TABLE]

Now, exploiting the expansion of $J_{0}(x)$, as $x\rightarrow+\infty$ (see [44]), namely:

[TABLE]

and substituting it on the second term of (41), we get that the left-hand side in (41) is given by

[TABLE]

Expanding the fourth power we obtain

[TABLE]

which is equal to

[TABLE]

Solving the integral of the square of the cosine, we get

[TABLE]

and then the logarithm terms cancel, leading to the expression

[TABLE]

Now, using the fact that

[TABLE]

and

[TABLE]

where $C_{i}(\cdot)$ and $S_{i}(\cdot)$ are the cosine and sine integral functions, respectively, we can approximate (42) with

[TABLE]

and computing the value of $\int_{0}^{10}J_{0}(\psi)^{4}\psi\,d\psi$ numerically, we find that (46) is equal to 0.298, in very good agreement with the value in (40).

Finally, we summarize, in the following tables, the analytic formulas used, for the Lipschitz-Killing curvatures, in the simulations.

We stress again that here the area of the sphere has been normalized to $|\mathbb{S}^{2}|=1$, hence, to obtain these statistics when $|\mathbb{S}^{2}|=4\pi$ we need to multiply the mean for $4\pi$ and the variance for $16\pi^{2}$, for the area and the Euler-Poincaré characteristic. For the boundary length we recall that there is a further factor 2 to take into account (boundary length $=2\mathcal{L}_{1}(A_{u}(f_{\ell};\mathbb{S}^{2}))$), so that we need to multiply for $4\pi\times 2$ to obtain the expected value and for $16\pi^{2}\times 4$ for the variance. The asymptotic behavior at $u=0$ of the Euler Poincaré characteristic is easily seen to be $O(\ell^{2}\log\ell)$ (exploiting results on extrema and saddles) but a rigorous evaluation of the leading constant is still missing.

We conclude reporting in the following table the formulae exploited for expected values and variances for the critical points, recalling that

[TABLE]

7 Acknowledgments

We acknowledge the use of the National Energy Research Scientific Computing Center (NERSC) super-computing facilities. Maps and results have been derived using the HEALpix (http://healpix.jpl.nasa.gov) software package developed by [14]. Some of the theoretical results exploited in this paper were derived in collaboration with Igor Wigman, to whom we are grateful for many insightful discussions.