text:
# Probing CMB Non-Gaussianity Using Local Curvature

## 1 Introduction

From the now well measured temperature fluctuations in the cosmological microwave background (CMB) we can gain some unvaluable constraints on the physics of the early universe. An important issue lies in determining whether these fluctuations are of Gaussian nature or not. Indeed, most of inflation scenarios predict a very low level of non-Gaussianity while models involving topological defects can give rise to significantly non Gaussian fluctuations. Even if a high degree of non-Gaussianity seems disfavored by recents measurements \[Netterfield et al. 2001, Lee et al. 2001, Halverson et al. 2001, Santos et al. 2001, Polenta et al. 2002\], the presence of topological defects at a significant level is not yet ruled out by observations, as advocated recently by e.g. Bouchet et al. (2000), who considered a mixed model where Cosmic Microwave Background (CMB) temperature fluctuations are seeded in part by cosmic strings.

There are numerous ways of probing non Gaussian features of a random field such as a CMB temperature fluctuations map. Since statistical properties of random Gaussian fields are entirely determined by their two-point correlation function, a natural approach consists in measuring higher order correlation functions or related statistics such as cumulants of the distribution. The measurements can be done in real and Fourier space \[Hinshaw et al. 1995, Ferreira et al. 1998, Magueijo 2000, Banday et al. 2000, Gangui & Martin 2000a, Gangui & Martin 2000b, Verde & Heavens 2001, Komatsu et al. 2001, Santos et al. 2001, Kunz et al. 2001\] or in the space of wavelet coefficients \[Aghanim & Forni 1999, Aghanim & Forni 2001\].

Alternatively, non Gaussian features of CMB maps can be probed by the topological analysis their excursions sets, defined for a random field $`T(\theta )`$ as $`A_u(T)=\{\theta |T(\theta )u\}`$. For example, for a sufficiently smooth and non degenerate random field it is possible to measure the Euler characteristic (or genus) of the excursion sets via a counting of critical points111The points where the gradient of the field cancels., classified according to their curvature, i.e. maxima, minima or saddle points in our 2D case (e.g. Adler 1981, p. 87). Critical points counting is a well known practical issue in the context of CMB analysis and specific predictions can be derived for smooth random Gaussian fields (e.g., Adler 1981), which allow one to constrain the degree of non Gaussianity of CMB maps by measuring e.g. the Euler characteristic \[Coles et al. 1989, Gott et al. 1990, Luo 1994, Smoot et al. 1994, Barreiro et al. 2001\] or peak statistics \[Bond and Efstathiou 1987, Coles et al. 1989, Novikov & Jorgensen 1996, Heavens & Sheth 1999, Heavens & Gupta 2001\].

Following an approach advocated recently in the large scale structure context (3D) by Colombi, Pogosyan & Souradeep (2000), we propose in this paper to extend the counting to ordinary points, i.e. to classify all the points according to their local curvature as belonging to “hills”, “lakes” or “saddles” (see Fig. 1 for an example). By measuring the relative abundances, $`𝒫_{\mathrm{hill}}`$, $`𝒫_{\mathrm{lake}}`$ and $`𝒫_{\mathrm{saddle}}`$, of these three types of points for various excursions sets and smoothing scales, i.e. by exploring the correlation between height and curvature222In this sense, our work is similar in spirit as in Takada (2001)., we are able to extract a mathematically well defined Gaussian signature, depending formally on only one specific ratio of spectral parameters. As a result, a comparison of the measured abundances with those predicted in the Gaussian case allows us to detect a certain level of non-Gaussianity.

This paper is organized as follows. In section § 2, after recalling some useful results about 2D Gaussian random fields, we derive analytic predictions for the abundances, by extending the work of Bond & Efstathiou (1987). In section § 3 we discuss technical issues involved in the practical measurements of $`𝒫_{\mathrm{hill}}`$, $`𝒫_{\mathrm{lake}}`$ and $`𝒫_{\mathrm{saddle}}`$. In section § 4 we build a $`\chi ^2`$ statistic that allows us to combine a set of measurements in order to quantify the likelihood of an image to be Gaussian, from the standpoint of this measure. We then apply the method to the case of noisy mixed models involving cosmic strings plus cold dark matter. Finally, results are discussed in section § 5.

## 2 Theoretical Considerations

In this section we recall some specific properties of 2D Gaussian random fields, following Bond & Efstathiou (1987) (BE87) (a 2D transcription of the 3D formalism of Bardeen et al. 1986). We extend the calculations of this work to obtain theoretical expressions for the abundances of interest, i.e. the fractions of space occupied respectively by hill, lake and saddle regions.

### 2.1 Facts about a Gaussian random field

Let us first consider a temperature fluctuation field, $`\delta _\mathrm{T}(𝐫)=T(𝐫)/\overline{T}1`$ whose $`2D`$ Fourier transform is $`\delta _\mathrm{T}(𝐤)d^2𝐤\delta _\mathrm{T}(𝐫)\mathrm{exp}(i𝐤.r)`$ and whose power spectrum $`P(k)`$ is

$$\delta _\mathrm{T}(𝐤)\delta _\mathrm{T}(𝐤^{})=(2\pi )^2P(𝐤)\delta _D(𝐤+𝐤^{}),$$
(1)
where $`\delta _D(𝐤)`$ is the Dirac distribution. We can define the moments of the power spectrum

$$\sigma _j^2\frac{1}{(2\pi )^2}d^2𝐤P(k)k^{2j}$$
(2)
and the ratio

$$\gamma \frac{\sigma _1^2}{\sigma _0\sigma _2}.$$
(3)
Note that $`\sigma _0`$ is the rms of the fluctuation field $`\delta _\mathrm{T}`$. Naturally, in the flat sky approximation, i.e. large multipole $`l`$ and small separation angles, the Fourier (flat) power spectrum is related to spherical harmonic power spectrum as follows

$$P(k)C_{\mathrm{}}\mathrm{and}k\mathrm{}.$$
(4)
Let us consider from now on only the normalized fluctuation field

$$\delta (𝐱)\frac{\delta _\mathrm{T}(𝐱)}{\sigma _0}.$$
(5)
At a given point $`𝐫`$, we note the gradient of the field $`\delta \eta `$ and the Hessian matrix, $`\zeta _{ij}\delta /x_ix_j`$. This symmetric real matrix can be diagonalised by applying a rotation of an angle $`\theta `$ and we note $`\lambda _1`$ and $`\lambda _2`$ the *opposite* of its eigenvalues, with $`\lambda _1\lambda _2`$. The normalized trace of the curvature matrix reads,

$$x\frac{\lambda _1+\lambda _2}{\sigma _2},$$
(6)
and the ellipticity $`e`$ is defined as

$$e=\frac{\lambda _1\lambda _2}{2\sigma _2x}.$$
(7)
Note that with this notation, *$`x`$ and $`e`$ should have the same sign*. If $`\lambda _20`$ then $`x0`$ and $`0e1/2`$, if $`\lambda _10`$ then $`x0`$ and $`1/2e0`$. If $`\lambda _1`$ and $`\lambda _2`$ have opposite signs, then neither $`x`$ nor $`e`$ are restricted.

For a Gaussian random field, the joint probability distribution function of $`\delta `$, $`\eta `$, $`x`$, $`e`$ and $`\theta `$ is given by \[Eq. (A1.6) in BE87\]

$`𝒫(\delta ,\eta ,x,e,\theta )d\delta dxded\theta d^2\eta =`$ (8)
$`e^{\delta ^2/2}{\displaystyle \frac{\mu d\delta }{\sqrt{2\pi }}}e^{\mu ^2(x\gamma \delta )^2/2}{\displaystyle \frac{dx}{\sqrt{2\pi }}}\times `$
$`e^{\eta ^2/\sigma _1^2}{\displaystyle \frac{d^2\eta }{\pi \sigma _1^2}}e^{4x^2e^2}8x^2ede{\displaystyle \frac{d\theta }{\pi }},`$
where we set

$$\mu (1\gamma ^2)^{1/2}.$$
(9)
Knowing this probability distribution function, it is easy to show that $`\delta ^2=x^2=\eta ^2=1`$, that $`\eta `$ is correlated with neither $`\delta `$ nor $`x`$ nor $`e`$, but that $`\delta x=\gamma `$, i.e. the height of a point is correlated with the curvature at this point. This lattest result is of particular interest to us since the so-induced non-trivial dependence on thresholding will generate a specific gaussian signature that we will exploit in the following.

### 2.2 Studying the local curvature

BE87 calculated the density probability of extrema for a Gaussian random field as a function of threshold. So using the properties of *first derivatives*, they open the way to further works \[Heavens & Sheth 1999, Heavens & Gupta 2001\] that illustrate the use of the number and spatial distribution of extrema to characterize Gaussian random fields.

We aim at extending this work to *second derivatives* by characterizing the local curvature *at any point* and then by comparing the abundance of $`3`$ defined types.

The local curvature is defined by the Hessian. By considering the sign of its eigenvalues ($`\lambda _1`$ and $`\lambda _2`$) we are led to distinguish three families of points :

* the “hill” points for which both eigenvalues are negative, i.e. $`x0`$ and $`0e1/2`$;
* the “lake” points for which both eigenvalues are positive, i.e. $`x0`$ and $`1/2e0`$;
* the “saddle” points for which the eigenvalues are of different signs.

The first, second and third family incorporate respectively maxima, minima and saddle points.

Extending the calculation of BE87 we now compute the probability that a point ($`\delta ,\eta ,x,e,\theta )`$ above a given threshold $`\delta _{\mathrm{th}}`$, i.e. such that $`\delta \delta _{\mathrm{th}}`$, belongs to any of these three classes. More precisely we calculate the quantities $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})𝒫(0x,0e1/2|\delta _{\mathrm{th}}\delta )`$, $`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}})𝒫(0x,0e1/2|\delta _{\mathrm{th}}\delta )`$ and $`𝒫_{\mathrm{saddle}}(\delta _{\mathrm{th}})=1𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}})`$.

To compute the quantity $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})`$ we first estimate the probability $`𝒫(0x,0e1/2,\delta _{\mathrm{th}}\delta )`$. Starting from the distribution, Eq. (8), for the cases of interest to us, we can perform straightforwardly both the integration $`_{\mathrm{}}^{\mathrm{}}d^2\eta /(\pi \sigma _1^2)`$ and $`_0^\pi 𝑑\theta /\pi `$. Doing so yields a probability function $`𝒫(\delta ,x,e)`$. Note that the fact that $`\theta [0,\pi ]`$ is due to the chosen ordering of the eigenvalues ($`\lambda _1\lambda _2`$) implicit in the distribution Eq. (8). The subsequent integration $`_0^{1/2}𝑑e`$ yields the differential density

$`𝒩_{\mathrm{hill}}(\delta ,x,0e1/2)d\delta dx=`$ (10)
$`e^{\delta ^2/2}{\displaystyle \frac{\mu d\delta }{(2\pi )^{1/2}}}e^{\mu ^2(x\gamma \delta )^2/2}(1e^{x^2}){\displaystyle \frac{dx}{(2\pi )^{1/2}}}.`$
Then we can still perform analytically the integration over $`x0`$ and we get

$`𝒩_{\mathrm{hill}}(\delta ,0x,0e1/2)d\delta =`$ (11)
$`e^{\delta ^2/2}{\displaystyle \frac{\mu d\delta }{2\sqrt{2\pi }}}[{\displaystyle \frac{1}{\mu }}(1+\mathrm{Erf}\left({\displaystyle \frac{\gamma \mu \delta }{\sqrt{2}}}\right))`$
$`{\displaystyle \frac{e^{\mu ^2\gamma ^2\delta ^2/(2+\mu ^2)}}{\sqrt{2+\mu ^2}}}(1+\mathrm{Erf}\left({\displaystyle \frac{\gamma \mu ^2\delta }{\sqrt{2}\sqrt{2+\mu ^2}}}\right))].`$
The integration over some threshold $`\delta _{\mathrm{th}}`$ cannot be performed analytically. However, knowing that

$$𝒫(\delta _{\mathrm{th}}\delta )=\frac{1}{2}\mathrm{Erfc}\left(\frac{\delta _{\mathrm{th}}}{\sqrt{2}}\right),$$
(12)
it is easy to evaluate numerically the quantity

$$𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})=\frac{_{\delta _{\mathrm{th}}}^{\mathrm{}}𝒩_{\mathrm{hill}}(\delta ,0x,0e1/2)d\delta }{𝒫(\delta _{\mathrm{th}}\delta )}.$$
(13)
We can still obtain analytically the limiting value $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}}=\mathrm{})`$, i.e. the fraction of “hill” points in the absence of threshold. We find

$$𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}}=\mathrm{})=\frac{1}{2}\left(1\frac{1}{\sqrt{3}}\right)0.2113.$$
(14)
An analogous calculation can be performed for “lake” points. It leads to the differential density

$`𝒩_{\mathrm{lake}}(\delta ,x0,1/2e0)d\delta =`$ (15)
$`e^{\delta ^2/2}{\displaystyle \frac{\mu d\delta }{2\sqrt{2\pi }}}[{\displaystyle \frac{1}{\mu }}\mathrm{Erfc}({\displaystyle \frac{\gamma \mu \delta }{\sqrt{2}}})`$
$`{\displaystyle \frac{e^{\mu ^2\gamma ^2\delta ^2/(2+\mu ^2)}}{\sqrt{2+\mu ^2}}}\mathrm{Erfc}\left({\displaystyle \frac{\gamma \mu ^2\delta }{\sqrt{2}\sqrt{2+\mu ^2}}}\right)]`$
from which we can deduce numerically $`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}})`$ as in Eq. (13). The asymptotic value $`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}}=\mathrm{})`$ can also be obtained analytically, and from parity argument one finds

$`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}}=\mathrm{})`$ $`=`$ $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}}=\mathrm{})`$
$`=`$ $`{\displaystyle \frac{1}{2}}(1{\displaystyle \frac{1}{\sqrt{3}}})0.2113.`$
Two remarkable properties of Gaussian random fields emerge from the above analytical calculations:

* the evolution of the “hill”, “lake” and “saddle” point fractions as functions of threshold $`\delta _{\mathrm{th}}`$ depends only on the spectral parameter $`\gamma `$;
* the asymptotic value for $`\delta _{\mathrm{th}}\mathrm{}`$ is independent of any spectral parameter, i.e. is the same for any Gaussian random field. This would not be true if we were considering only maxima, minima or saddle points (see BE87).

To illustrate these results, we draw on Fig. 2 the functions $`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}})`$ and $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})`$ for a set of $`\gamma `$ values. Except in the limit $`\delta _{\mathrm{th}}\mathrm{}`$, the evolution of these functions with $`\delta _{\mathrm{th}}`$ is rather sensitive to $`\gamma `$.

At this point, it is important to be aware that till now, we have supposed an idealistic, infinite resolution experiment, i.e. we did not take into account any *beam smearing* effect that would affect any real measurement. Consequently, in order to compare the predictions to any true measurements, it would be necessary to incorporate this effect in our calculations. To do so, we should consider the convolved field, $`\delta B`$, where $`B`$ stands for the instrumental beam response, instead of the idealistic fluctuation field, $`\delta `$. Since linear transformations such as the convolution by an arbitrary beam function do not change the Gaussian nature of a random field, the beam smearing should appear only as a dependence of $`\gamma `$ on the smoothing scale $`\sigma _\mathrm{b}`$, that we will note $`\gamma (\sigma _\mathrm{b})`$. Furthermore, although the problem of beam convolution in its full generality is very intricate, it turns out to be analytically tractable for a Gaussian, symmetric beam, a reasonable approximation in practice. We will illustrate this point in more details below.

## 3 Confronting Predictions and Measurements on Simulations

In order to test our ability to measure the functions $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})`$ and $`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}})`$,333Again, we restrict here our analysis to these two functions, since function $`𝒫_{\mathrm{saddle}}(\delta _{\mathrm{th}})`$, which is equal (by construction) to $`1𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}})`$, does not yield any further information. and to measure the incertainties for maps of limited extent, we performed several measurements on simulated maps. We henceforth limit ourselves to small simulated square patches of the sky of width $`12.5\mathrm{deg}`$. These patches are considered as being flat and pixelized with a $`512\times 512`$ Cartesian mesh of resolution $`\theta _{\mathrm{pix}}=1.5^{}`$.

This section is organized as follows: we first describe the measurement principles (§ 3.1) and apply them to noise free realization of Cold Dark Matter (CDM) ($`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`h=0.5`$, $`\mathrm{\Omega }_b=0.05`$ and $`n=1`$) Gaussian maps generated from their power spectrum (taking into account both the uniform distribution of the phases of the Fourier modes and the Rayleigh distribution of their modulus) 444The power spectrum has been computed making use of the publicly available CMBFAST code \[Seljak & Zaldarriaga 1996\] available at http://physics.nyu.edu/matiasz/ CMBFAST/cmbfast.html (§ 3.2); in particular, we examine the Gaussianity of the distribution function of measurements in order to be able to perform later $`\chi ^2`$ tests; finally we discuss practical issues related to finite resolution and finite volume effects (§ 3.3).

### 3.1 Principles of measurements

The measurements consist in determining the fractions $`𝒫_{\mathrm{hill}}`$ and $`𝒫_{\mathrm{lake}}`$ for an ensemble of excursion sets (subsets for which $`\delta \delta _{\mathrm{th}}`$) of an image smoothed at different scales $`\sigma _\mathrm{b}`$. To insure sufficient differentiability, i.e. at least continuity of second order derivatives, we choose a Gaussian smoothing window.

Similarly as in Colombi et al. (2000) for the 3D case, to measure the local curvature at a given point, we find the quadratic form which fits this point and its $`8`$ closest neighbors. We deduce from the quadratic form coefficients the local values of the second derivatives, i.e. the Hessian matrix, then diagonalize this matrix and determine the sign of its eigenvalues $`\lambda _1`$ and $`\lambda _2`$. Both easy to implement and fast, this method turns out to be quite robust in the presence of noise.

### 3.2 Example: application to noise free realizations

As a first example, we apply this measurement technique to $`200`$ simulated noise free realizations of CDM like temperature maps smoothed by a Gaussian beam of $`\sigma _\mathrm{b}=5\theta _{\mathrm{pix}}=7.5^{}`$. In each realization, we measure both $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})`$ and $`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}})`$ for the ensemble of excursion sets defined by $`\delta _{\mathrm{th}}=3.0`$, $`2.5`$, $`2.0\mathrm{}3.0`$, as illustrated by Fig. 3. For comparison purpose, we plotted on this figure the theoretical curves corresponding to the expected value of the parameter $`\gamma (\sigma _\mathrm{b}=5^{})=0.48`$. This effective $`\gamma `$ can be easily calculated for a Gaussian beam, since in this context we just have to replace $`P(k)`$ (an exact theoretical input here) by $`P(k)\mathrm{exp}(k^2\sigma _\mathrm{b}^2)`$ in Eqs. (2) and (3), (see Fig. 6 and section § 3.3 for a mode detailled illustration of this calculation). With this choice of $`\gamma `$, the theoretical curves fit very well the data. The dispersion over the measurements is very tight except at high threshold where we enter a rare events regime, as indicated by the larger $`1\sigma `$ error bars.

To illustrate more visually the protocol, Fig. 5 displays one CDM realization (noise free again) smoothed at $`\sigma _\mathrm{b}=3\theta _{\mathrm{pix}}`$, the corresponding excursion set for $`\delta 1`$ as well as the local curvature map, where $`hill`$, $`lake`$ and $`saddle`$ points are identified by respectively black, white and grey points. A cross examination of these images provides some insight on the processe involved, i.e. the correlation of extrema and local curvature.

To further investigate the dispersion over the measured values of $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})`$ at a given $`\delta _{\mathrm{th}}`$, we drawn on Fig. 4 the probability distribution function of the hill point fraction measured from the 200 realizations \[similar results can be obtained for $`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}})`$\]. Two thresholds are considered, $`\delta _{\mathrm{th}}=2`$ (left pannel) and $`\delta _{\mathrm{th}}=0`$ (right pannel). The dashed line on each pannel corresponds to the Gaussian limit with same standard deviation as the one actually measured. It fits qualitatively well the measurements, so the previously drawn $`1\sigma `$ errors are meaningful estimates of the cosmic variance errors.555However, one must keep in mind that error bars for different values of $`\delta _{\mathrm{th}}`$ are correlated with each other, but this will be taken into account in the analyses conducted later. We checked wether the distribution function of the measurements is Gaussian for other values of $`\delta _{\mathrm{th}}`$ and found that it is indeed the case except at high threshold, $`\delta _{\mathrm{th}}2.5`$, where one enters the rare events regime.

Thus, the measurements in this *idealistic* noise free case are in very good agreement with theoretical expectations. Furthermore, the cosmic errors associated with these measurements are very tight for a Gaussian field. This last point will be of particular interest to us since it makes the evolution of functions $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})`$ and $`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}})`$ with $`\delta _{\mathrm{th}}`$ a sharp signature of random Gaussian fields.

### 3.3 Measuring $`\gamma `$: spurious effects and available dynamic range

In the above discussion, we used a *predicted* value of $`\gamma `$ to check if the measured “hill”, “saddle” and “lake” point fractions agreed with theoretical predictions for a Gaussian random field. Therefore, the only thing we proved so far is that for a finite realization of a Gaussian random field smoothed with a Gaussian of width $`\sigma _\mathrm{b}=3\theta _{\mathrm{pix}}`$, whose (unsmoothed) power spectrum is perfectly known, the predicted value $`\gamma `$ agrees very well with the measured points. In practice however, the knowledge of the power spectrum might be less accurate and we want to develop a non-Gaussianity test which can be performed both independently or in combination with the power spectrum measurement. Thus it is important to determine from the measurements of these fractions themselves an effective value, $`\gamma _{\mathrm{eff}}`$, that fits well the data. We will see that for a finite realization of a Gaussian random field, there always exists a value $`\gamma _{\mathrm{eff}}`$ which is such that measurements agree very well with analytic predictions.

As we shall see below, the determination of the $`\gamma _{\mathrm{eff}}`$ value is easy and this fact, by itself, is highly significant and is the key point of our paper: for a map in which temperature fluctuations are partly seeded by cosmic strings, we shall see that it is actually not possible to find a value of $`\gamma _{\mathrm{eff}}`$ leading to point fractions matching the measurements in all the available dynamic range – namely, for all possible values of $`\delta _{\mathrm{th}}`$ and $`\sigma _\mathrm{b}`$. However, as we shall see later, it is not needed to accurately determine the real value of $`\gamma `$ to efficiently constrain the level of non Gaussianity of a map.

However, studying the difference between $`\gamma `$ and $`\gamma _{\mathrm{eff}}`$ can help to determine the available dynamic range, i.e. the set of values of $`\delta _{\mathrm{th}}`$ and $`\sigma _\mathrm{b}`$ in which one can trust the measurements. We already noticed in previous section that the density threshold should be small enough, $`\delta _{\mathrm{th}}2.5`$, to avoid entering the rare event regime, where the cosmic distribution function of measurements becomes non Gaussian. Here, we are concerned by two effects that we ignored previously:

* Pixelization effects, or equivalently, effects of finite resolution: as discussed above, it is important to smooth the data to insure sufficient differentiability. The smoothing scale should be large enough compared to the pixel size to avoid anisotropies or discreteness effects brought by the pixelization.
* Finite volume and edge effects: our experimental maps cover a rather small fraction of the sky, due to the limited dynamical range in the cosmic-string simulations from \[Bouchet et al. 1988, Bennett & Bouchet 1990\]. Therefore, to avoid reducing too much the number of statistically independent regions of the map, the smoothing scale should not be too large. In a more realistic experiment covering a large fraction of the sky, finite volume and edge effects should not be as much as of a concern.

The practical measurement of $`\gamma _{\mathrm{eff}}`$ is made straightfoward by noticing that the dependence of the function $`\stackrel{~}{𝒫}_{\mathrm{hill}}(\delta _{\mathrm{th}}=1)`$ on $`\gamma `$ is very well approximated by an exponential law666Note that such a fit is also appropriate for values of $`\delta _{\mathrm{th}}`$ different from 1. in the domain of interest to us, i.e. $`0.4\gamma 0.95`$: the fit

$$\mathrm{ln}\stackrel{~}{𝒫}_{\mathrm{hill}}(\delta _{\mathrm{th}}=1)a+b\gamma $$
(17)
with $`a=1.4099`$ and $`b=1.2395`$ is accurate to $`0.6\%`$. The choice of the particular value $`\delta _{\mathrm{th}}=1`$ is an ad hoc compromise coming from the competition between two effects: (i) the dependence on $`\gamma `$ of the quantity $`\stackrel{~}{𝒫}_{\mathrm{hill}}(\delta _{\mathrm{th}})`$ increases with $`\delta _{\mathrm{th}}`$ (Fig. 2), but (ii) the uncertainty on the determination of $`\stackrel{~}{𝒫}_{\mathrm{hill}}(\delta _{\mathrm{th}})`$ increases with $`\delta _{\mathrm{th}}`$ (Fig. 3).

To test the effects of finite coverage and finite resolution, we examine a CDM case, where we generated once again some Gaussian field realisations from their $`C_{\mathrm{}}`$. Assuming as before that the field has been smoothed by a Gaussian of width $`\sigma _\mathrm{b}`$, we can easily compute the prediction for the function $`\gamma (\sigma _\mathrm{b})`$.

In Fig. 6, the measured values of $`\gamma _{\mathrm{eff}}`$ is displayed in each case as a function of smoothing scale. An average over 50 realizations of $`512\times 512`$ pixels maps is performed. The error bars on the figure represent the corresponding scatter, which increases with $`\sigma _\mathrm{b}`$ as expected, due to finite volume effects. The solid line corresponds to the theoretical function $`\gamma (\sigma _\mathrm{b})`$. Agreement between $`\gamma _{\mathrm{eff}}`$ and the analytic prediction is good, except at small scales where pixelisation effects contaminate this measurements. We however see that pixelisation effects become negligible when $`\sigma _\mathrm{b}/\theta _{\mathrm{pix}}>3`$. Note that there should be an upper bound as well for $`\sigma _\mathrm{b}`$ as discussed above, since cosmic errors increase with scale. Our $`\chi ^2`$ analysis below will anyway naturally take that into account by giving a lesser weight to scales with larger errors. This reasonable agreement between the measured $`\gamma _{\mathrm{eff}}`$ at a given scale $`\sigma _\mathrm{b}`$ and the expected $`\gamma (\sigma _\mathrm{b})`$ validates the approach we will use in the next section which consists in finding from the fraction themselves an ad hoc $`\gamma _{\mathrm{eff}}`$, and deducing from it an agreement with Gaussianity.

## 4 Testing Non-Gaussianity : Principle and Application to Mixed Models

¿From the previous section, we conclude that in the case of a smoothed Gaussian random field, the functions $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})`$ and $`𝒫_{\mathrm{void}}(\delta _{\mathrm{th}})`$ can be measured accurately and fit very well the analytical predictions provided $`\gamma `$ be considered as an adjustable parameter. Furthermore, their probability distribution function is well approximated by a Gaussian if $`\delta _{\mathrm{th}}2.5`$, which now allows us to define a rigorous measurement protocol based on $`\chi ^2`$ analysis that we shall apply to simulated data in the trustable scale range determined above, namely $`\sigma _\mathrm{b}3`$.

In this section, we first detail how we use the functions $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})`$ and $`𝒫_{\mathrm{void}}(\delta _{\mathrm{th}})`$ to build a $`\chi ^2`$ statistic testing Gaussianity (§ 4.1). Then we apply the method to simple simulated mixed models involving cosmic strings plus CDM and see how it can be used to bound quantitatively the relative contribution of cosmic strings (§ 4.2).

### 4.1 The method

Let us assume as before that we have at our disposal an observed temperature map $`\delta _\mathrm{T}`$, on which we measured the fractions $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}_\mathrm{i}},\sigma _{\mathrm{b}_\mathrm{i}})`$ and $`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}_\mathrm{i}},\sigma _{\mathrm{b}_\mathrm{i}})`$ for a set of threshold values $`\{\delta _{\mathrm{th}_\mathrm{i}}\}`$ and smoothing scales $`\{\sigma _{\mathrm{b}_\mathrm{i}}\}`$. If the power spectrum of this random field is known, we can in principle define analytically a set of values of $`\gamma `$ for each smoothing scale, $`\{\gamma _{\sigma _{\mathrm{b}_\mathrm{i}}}\}`$. Assuming furthermore that the field is Gaussian and that our estimators of $`𝒫_{\mathrm{hill}}`$ and $`𝒫_{\mathrm{lake}}`$ are unbiased, it is easy to deduce from this $`\gamma `$ set some *theoretical expectations* for the set, $`\stackrel{~}{𝒫}_\mathrm{X}(\delta _{\mathrm{th}_\mathrm{i}},\sigma _{\mathrm{b}_\mathrm{i}})𝒫_\mathrm{X}(\delta _{\mathrm{th}_\mathrm{i}},\sigma _{\mathrm{b}_\mathrm{i}})`$, where $`\mathrm{X}`$ stands from now on for “lake” or “hill”. To quantify the distance between these theoretical predictions and the measurements, we introduce the standard $`\chi ^2`$ statistic. Since these measurements can obviously not be considered as independent, we have to introduce the full theoretical variance-covariance matrix

$$C_{\mathrm{II}^{}}(𝒫_\mathrm{I}\stackrel{~}{𝒫}_\mathrm{I})(𝒫_\mathrm{I}^{}\stackrel{~}{𝒫}_\mathrm{I}^{}),$$
(18)
where we use the short hand notation $`𝒫_\mathrm{I}=𝒫_\mathrm{X}(\delta _{\mathrm{th}_\mathrm{i}},\sigma _{\mathrm{b}_\mathrm{i}})`$ and $`𝒫_\mathrm{I}^{}=𝒫_\mathrm{X}^{}(\delta _{\mathrm{th}_\mathrm{i}^{}},\sigma _{\mathrm{B}_\mathrm{i}^{}})`$. Then, the statistic

$$\chi ^2\underset{\mathrm{II}^{}}{}(𝒫_\mathrm{I}\stackrel{~}{𝒫}_\mathrm{I})(C_{\mathrm{II}^{}})^1(𝒫_\mathrm{I}^{}\stackrel{~}{𝒫}_\mathrm{I}^{}).$$
(19)
is expected to follow a $`\chi ^2`$-distribution, as illustrated by a practical example below. Indeed results of § 3.2 suggest that the distribution function of the measured $`𝒫_\mathrm{I}`$ is nearly Gaussian if $`\delta _{\mathrm{th}}2.5`$.

In practice, we have to introduce two more subtleties in this protocol.

First, as discussed extensively in § 3.3, the practical realization of a Gaussian random field yields a measured function $`𝒫_\mathrm{I}`$ matching the theory, $`\stackrel{~}{𝒫}`$, but with an effective value of $`\gamma `$. Taking into account some additional Gaussian noise, as discussed in more details below, would make this effect even stronger. Rigorously, we should modelize this $`\gamma _{\mathrm{eff}}`$ effect in terms of statistical bias on the estimators of the function $`𝒫_\mathrm{X}(\delta _{\mathrm{th}},\sigma _\mathrm{b})`$, but we proceed here differently, for simplicity: we extract the $`\gamma _{\mathrm{eff}}(\sigma _\mathrm{b})`$ from the data by fitting the analytic prediction for function $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}}=1,\sigma _\mathrm{b})`$ to the measured one. Thus, our “theory” depends itself on the data through $`\gamma _{\mathrm{eff}}`$. As already explained, the critical test for non Gaussianity will be on the evolution of the functions $`𝒫_\mathrm{X}(\delta _{\mathrm{th}},\sigma _\mathrm{b})`$ with $`\delta _{\mathrm{th}}`$, which, once $`\gamma _{\mathrm{eff}}`$ is determined, is entirely fixed. The practical measurement of $`\gamma _{\mathrm{eff}}`$ is performed as explained in the previous section. In principle we could include the value of $`\delta _{\mathrm{th}}`$ used to measure $`\gamma _{\mathrm{eff}}`$ as a varying parameter in our $`\chi ^2`$ test, to optimize the analysis. Given the level of accuracy of our numerical experiments and since we just aim here to illustrate the method in a simple and convincing way, we did not feel necessary to do so.

Second, even if in principle we could compute analytically or numerically the variance-covariance matrix, for the sake of simplicity we evaluated it using a Monte-Carlo method based on a few hundreds of Gaussian realizations having the same power-spectrum and noise properties as the input map. So for each model we will consider below (pure CDM, mixed model with CDM + cosmic strings, and pure cosmic strings), the covariance matrix $`C`$ should be different in each case. However, in practice we checked that this matrix coefficients depend very slightly on $`\gamma _{\mathrm{eff}}`$ in the range of interest to us, i.e. $`0.4\gamma 0.95`$ so that we consider for simplicity only the same $`C`$ matrix in each case. Another issue concerning this matrix is that it might be singular. Indeed, since fractions measured at various thresholds and different scales can be highly correlated, e.g. the values at very low threshold where the dependence on $`\gamma `$ is very weak, the studied $`\{\delta _{\mathrm{th}_\mathrm{i}}\}`$ and $`\{\sigma _{\mathrm{b}_\mathrm{i}}\}`$ sets have to be restricted so that $`C`$ is invertible.

### 4.2 Application on mixed models for MAP and Planck

As an illustration of the accuracy of our statistics, we apply it in the framework of MAP777http://map.gsfc.nasa.gov and Planck Surveyor888http://astro.estec.esa.nl/SA-general/Projects/Planck/ experiments to a mixed model where temperature fluctations are seeded in part by cosmic strings and otherwise by adiabatic inflationary perturbations in standard CDM model (e.g. Bouchet et al. 2001 and references therein). Our approach is rather simple, since we add together a CDM map and a temperature fluctuation map obained from ray-tracing in cosmic string simulations \[Bouchet et al. 1988, Bennett & Bouchet 1990\] and neglect possible cross-correlations between such maps. This implicitely supposes the existence of a standard CDM inflationary epoch followed by a phase transition during which cosmic strings appear and then imprint supplementary fluctuations on the dark matter distribution. Neglecting cross-correlations mentionned just above is equivalent to assume that the distribution of cosmic string as well as their dynamical evolution is not coupled with that of cold dark matter, which should be a good assumption at the level of approximation considered in this paper \[Linde & Riotto 1997, Contaldi et al. 2001\]. We neglect also any contribution to the fluctuations prior to the last scattering surface (lss) so that the comparaison with results from Bouchet et al. 2001, who considered this contribution, is not immediate at small scales.

To take into account the noise expected in the best channels of MAP and Planck Surveyor experiments (respectively $`12.8\mu \mathrm{K}/\mathrm{K}`$ per $`0.3\times 0.3\mathrm{deg}^2`$ pixels and $`2\mu \mathrm{K}/\mathrm{K}`$ per $`8^{}\times 8^{}\mathrm{deg}^2`$), we add an extra Gaussian white noise $`n`$ to our square maps. Note again, that as long as the noise is Gaussian, its presence can be simply taken into account by a change in the value of $`\gamma _{\mathrm{eff}}`$. Thus, our simulated temperature maps read

$$\delta _\mathrm{T}=(1\beta )\delta _\mathrm{T}^{\mathrm{string}}+\beta \delta _\mathrm{T}^{\mathrm{CDM}}+n,$$
(20)
where $`1\beta `$ represents the fraction of the signal seeded by cosmic strings, given the fact that we impose $`\mathrm{rms}(\delta _\mathrm{T}^{string})=\mathrm{rms}(\delta _\mathrm{T}^{CDM})=\mathrm{rms}(\delta _\mathrm{T}^{COBE})`$ for the fields smoothed at a scale corresponding to the COBE-DMR beam, i.e. $`\sigma _\mathrm{b}=7^{}/\sqrt{8\mathrm{ln}2}`$. Note that we differ from Bouchet et al. 2001 in the relative normalisation of the two contributions since they define their relative normalisation with $`C_{\mathrm{}}`$, i.e. they considered a family of models parametrized by $`C_{\mathrm{}}=(1\alpha )C_{\mathrm{}}^{\mathrm{string}}+\alpha C_{\mathrm{}}^{\mathrm{CDM}}`$, where both $`C_{\mathrm{}}^{\mathrm{string}}`$ and $`C_{\mathrm{}}^{\mathrm{CDM}}`$ are COBE normalised and where $`C_{\mathrm{}}^{\mathrm{string}}`$ includes contributions from before and after the lss. If we wanted however to express our results, in terms of $`\beta `$, we should consider $`1\alpha =(1\beta )^2/((1\beta )^2+\beta ^2)`$.

We now construct the previously defined $`\chi ^2`$ statistics for various values of $`\beta `$ and try to figure out what is the smallest “string like contribution” that our method should be able to detect with a high statistical significance.

#### 4.2.1 Constructing the $`\chi ^2`$ and adequacy of the $`\chi ^2`$ probability distribution

Even if one would like to include as many threshold values and smoothing scales as possible in the analysis, this is in practice not necessary. Indeed, some measurements are, at least in the Gaussian case, highly correlated with each other, due to e.g. some strong constraints (in the very low threshold regime, $`\delta _{\mathrm{th}}\mathrm{}`$, both $`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}})`$ and $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})`$ tend to an asymptotic value independent of the spectral parameters). Therefore, one can restrict the set of values of $`(\sigma _B,\delta _{\mathrm{th}})`$ by requiring that the $`C`$ matrix (see Eq.18) be non-singular.

By investigating the numerical properties of $`C`$, we find for all the maps 999As already stated, we checked that this choice is adequate for various $`C`$ matrices, i.e. corresponding to different $`\gamma `$s, altough we use only one in practice. that $`\delta _{\mathrm{th}}=2`$,$`1.5\mathrm{}`$, $`+1.5`$ and $`\sigma _\mathrm{b}/\theta _{\mathrm{pix}}=6`$, $`10`$ and $`14`$ (corresponding to repectively $`21^{}`$, $`35^{}`$ and $`49^{}`$ on the sky) are satisfying values for the Planck case, whereas $`\delta _{\mathrm{th}}=1.5`$,$`1.0\mathrm{}`$, $`+1.5`$ and $`\sigma _\mathrm{b}/\theta _{\mathrm{pix}}=8`$ and $`12`$ (corresponding to respectively $`28^{}`$ and $`42^{}`$) are reasonable for MAP. Of course, this choice of the set of values of $`(\sigma _B,\delta _{\mathrm{th}})`$ is strongly influenced by the fact that our map is rather small: for a full sky survey, the sampling would probably be different.

Given these restrictions, it is now important to explore the properties of the probability distribution function of the $`\chi ^2`$ statistics defined in Eq. (19). As stated previously, we expect it to follow a $`\chi ^2`$-distribution, at least in the Gaussian case, since the distribution of measured values of $`𝒫_I`$s is well described by a Gaussian (see § 3.1). We tested this by generating 300 realizations of noisy (Planck level) CDM like only temperature fluctuations for which we measure $`\chi ^2`$. Its probability distribution is drawn on Fig. 7, where we superimposed a theoretical $`\chi ^2`$-distribution with $`47`$ degrees of freedom (dof) ($`48`$ measurements: $`6\times 4`$ measured hill points fraction and $`6\times 4`$ measured void point fractions; $`1`$ parameter, $`\gamma `$). The agreement is obviously very good, comforting the analysis procedure and making this $`\chi ^2`$ statistic adequate for assessing Non-Gaussianity.101010Here, we do not test if the $`\chi ^2`$ estimator always follows a $`\chi ^2`$-distribution, i.e. also for non-Gaussian maps. Such a test is not necessary since the prior we use to analyse the data assumes underlying Gaussianity (it would be impossible to do otherwise in practice).

#### 4.2.2 Ability to distinguish a “Non-Gaussianity Level” with a realistic noise level

Making use of this $`\chi ^2`$ statistics, we are now in position to determine how well we can distinguish non-Gaussian signatures.

In Figs. 8, 10 and 10, we display some measurements of functions $`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}})`$ and $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})`$ in the Planck case. Three values of $`\beta `$ are considered: pure CDM with $`\beta =1`$, hybrid model with $`\beta =0.3`$ and pure string model with $`\beta =0`$. Again, On each figure, the Gaussian limit obtained from the procedure explained in detail in § 4.1 is compared to the measurements. The error bars on the theoretical predictions correspond to then the $`1\sigma `$ dispersion over 300 realizations of pure random Gaussian fields with same power-spectrum and noise properties as the data, as explained in end of § 4.1. Very good agreement is found for CDM, as expected, while the hydrid model exhibits some discrepancies, especially for $`𝒫_{\mathrm{lake}}`$. These discrepancies are overwhelming in the string only model. Note that since we fit the measured value of $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}}=1)`$ to estimate the Gaussian prediction, one expects non Gaussian features to show up more in the measured shape of function $`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}})`$ than that of function $`𝒫_{\mathrm{hill}}(\delta _{\mathrm{th}})`$, as indeed seen at least for the mixed model. Note also that if for the mixed model, the measured $`𝒫_{\mathrm{lake}}(\delta _{\mathrm{th}})`$ tends to lie below the Gaussian prediction, the reverse happens for the pure string model (this seemingly counterintuitive result stems from our comparing the data to a “theoretical” curve which is different in the 2 cases) . Finally, it is important as well to recall that data points plotted on Figs. 810 represent only 2/3 of all the points used for the $`\chi ^2`$ computation since we consider one more smoothing scale, $`\sigma _\mathrm{b}/\theta _{\mathrm{pix}}=6`$.

To illustrate more qualitatively these measurements, Figures 11 and 13 display, similarly as in Fig. 5, an example of initial smoothed noise free map, its thresholded counterpart (with $`\delta _{\mathrm{th}}=1`$) and the corresponding local curvature maps with lake, saddle and hill points. Figure 11 corresponds to the hybrid model case ($`\beta =0.3`$) while Figure 13 corresponds to the string only model ($`\beta =0.0`$). A strong similarity can be seen between the mere CDM and the hybrid models despite the differences in measured $`𝒫_{\mathrm{lake}}`$ functions. Note also the peculiar patterns in the pure string curvature map: there seem to be rather extended saddle point regions, which calls for other specific pattern detection statistics.

Table 1 summarizes the results obtained for $`\chi ^2`$ using the noise level of the best channel of either Planck or MAP. It shows, for various values of $`\beta `$, the measured $`\chi ^2`$, the $`\chi ^2`$ per dof, as well as the probability (noted Prob.) that a random variable, in a $`\chi ^2`$-distribution with the relevant (number of) dof (shown above to be adequate) is greater than the obtained $`\chi ^2`$. This quantity gives therefore the significance of the detection.

Roughly speaking, the smaller the $`\chi ^2`$, the better is the agreement with a Gaussian distribution. Reversely, the higher the obtained $`\chi ^2`$, the more unlikely does the analyzed map follow a Gaussian distribution.

Note that the numbers displayed in Table 1 represent in no way the best Planck and MAP can provide with our method, since we only considered 0.3% of the sky for one single channel. Nevertheless, these figures give some trends of the ability of our method, the purpose of this paper.

With these limitations in mind, we see that MAP can get a meaningful detection (1% level) of the non-Gaussianity induced by string like background with $`\beta =0.7`$, corresponding to a number $`1\alpha =0.15`$, a number slightly lower than the constrains from Bouchet et al. (2000), $`1\alpha =0.18`$. Any smaller non-Gaussian contribution is not detected with a good significance. On the other hand, we see that Planck could obtain, using this method, a highly significant detection level, i.e. smaller than 0.1%, for $`\beta `$ between 0.8 and 0.9, i.e. $`0.051\alpha 0.12`$.

## 5 Discussion and Conclusion

In this paper, we have tested a statistic based on local curvature measurements by counting three well defined types of points according to the sign of the eigenvalues of the Hessian of a temperature map. In the first section, we computed the relevant theoretical expectations. We validated them and demonstrated our ability to properly measure them in the next section and eventually introduced this statistic as an accurate non-Gaussianity test (by performing a $`\chi ^2`$ like analysis) in the last sections, where we applied it to the case of mixed models with realistic noise levels. Even if our physical simulation might be somewhat simplistic (e.g. no pre-lss string contribution, no foregrounds residuals, etc.), they showed clearly that we are able by using this technique in the Planck context, to distinguish a non-Gaussian background amplitude around 10% in temperature rms (5% in $`C_{\mathrm{}}`$) with a good significance using only 0.3% of the sky and one single channel. Consequently, these results are definitely positive and, as it is, this method seems to be ready to be applied on true data.

Naturally, in practice this test should be used in conjunction with other more standard statistical measurements, like 2-, 3-, 4-,…points correlation functions or their harmonic transform, power spectrum, bispectrum, trispectrum …or even already presented more sophisticated ones like the distribution of wavelet coefficients, etc. This array of methods might actually help in identifying the source of non-gaussianity detected in a realistic context. Indeed, in practical situations, important non-Gaussian effects might be induced by CMB contaminants, in particular galactic foregrounds, residual map stripping or point sources, or any instrumental systematic error.

### 5.1 Possible advantages

One important advantage of these statistic is its local nature, i.e. the fact that we are only interested by the 8 neighboring pixels. The sphericity of the sky is thus not an issue for this method and so is the pixelisation, provided that the smoothing scale is about 3 times the pixel size111111Which is convenient, since pixelisation of low noise experiments are usually done at a size similar to one-third of the FWHM of the beam. (see §3.3).

A second advantage of the method lies in the fact that it is mathematically well-defined and simple since it depends only on one single spectral parameter, $`\gamma `$ that can be either measured from the fractions themselves or from the $`C_{\mathrm{}}`$ measurement (see §3.3). Thus measurements are well understood and well controlled from a theoretical point of view.

### 5.2 Outlooks

As we mentioned earlier, by moving to second order derivatives, this work opens the way to many more characterization of random Gaussian field. We expect that an analogous work could be done by considering critical points only, which would lead, e.g. , to a measurement of the Euler characteristic as a function of threshold. However, we expect this restriction to be less efficient in probing non-Gaussianity since it appears that by keeping all types of points we benefit more from the full structure of the probability distribution function. This point will be demonstrated in a future work. In addition, before applying the method to real data, it will also be interesting to compare it’s ability to that of other estimators on the same test protocole, but extended to include other families of models, as e.g. \[Lyth & Wands 2001\].

## Acknowledgments

We are grateful to Roland Triay for organizing a prolific “Cosmological school” at Porquerolles where this work has been initiated, to Francis Bernardeau for useful discussions at an early stage of this work and to Christophe Pichon for valuable help on the use of Mathematica and keyboards. S.C. thanks Jean-Philippe Belial for quite useful suggestions on the methodological approach. We are grateful to Jérôme Martin and Patrick Peter for a careful reading of the manuscript.