Search for non-Gaussianities in the WMAP data with the Scaling Index Method

In the recent years, non-Gaussianity and statistical isotropy of the Cosmic Microwave Background (CMB) was investigated with various statistical measures, first and foremost by means of the measurements of the WMAP satellite. In this Review, we focus on the analyses that were accomplished with a measure of local type, the so-called Scaling Index Method (SIM). The SIM is able to detect structural characteristics of a given data set, and has proven to be highly valuable in CMB analysis. It was used for comparing the data set with simulations as well as surrogates, which are full sky maps generated by randomisation of previously selected features of the original map. During these investigations, strong evidence for non-Gaussianities as well as asymmetries and local features could be detected. In combination with the surrogates approach, the SIM detected the highest significances for non-Gaussianity to date.


I. INTRODUCTION
The inflationary phase, first proposed in 1981, is an important part of what is called the Standard Model of Cosmology. Since inflation occurred already a few moments after the Big Bang, where the Universe was extremely hot, dense and thus opaque, it is not possible to observe this short time period directly. The best way to achieve information about and to test theories for inflation, is to look at the temperature anisotropies of the Cosmic Microwave Background (CMB). While the simplest single field slow roll inflationary scenario predicts these fluctuations to be nearly Gaussian [1][2][3], a variety of more complex models could lead to a different result (e.g. [4][5][6][7]). By testing the Gaussianity of the CMB, it is possible to distinguish between different inflationary models, and therefore to shed some light on the physics of the very early Universe.
In this Review, we focus on another measure for non-Gaussianity analyses, namely the (weighted) scaling index method (SIM) [31,32]. The SIM is able to distinguish different structural behaviour of a data set in a local way. Scaling indices have already been used for texture discrimination [33] and feature extraction [34,35], time series analysis of stock exchanges [36] and active galactic nuclei [37,38], as well as structure analysis of bone images [39] and other different medical data, like biological specimens, skin cancer, computer tomographic images, and beat-to-beat sequences from electrocardiograms [40]. Investigations concerning the Gaussianity of the CMB by applying the SIM to simulated CMB maps, the WMAP 3-year, 5-year or 7-year data were performed in [32], [41], [42,43] and [44], respectively.
This Review is structured as follows: In section II, we present the scaling index technique for investigations on a sphere. The data sets that were used for investigations by means of scaling indices to date, including the technique of creating surrogate maps, are introduced in section III. The results of these analyses, in particular the detections of non-Gaussianities, asymmetries and local features, are outlined in section IV, followed by the conclusions in section V.

II. THE SCALING INDEX METHOD
Quite similar to wavelets, weighted scaling indices can be used to perform a local analysis of the data set, and can be calculated for different scales, which yields information about characteristic sizes of detected features. The measure has the ability of revealing the topological properties of an input map by detecting different structures in the data, as for example cluster-like or sheet-like structures, as well as filaments or walls. While wavelets are more sensitive to structures, which offer intensity variations of significant magnitude with respect to the existing noise, scaling indices also detect structural features which possess variations within the noise level, but not significantly higher or lower intensity values [41]. The In the first and the third plot, the full set of points is presented, while the second and fourth shows an x, z-projection of only the points with |y| < 0.05. Two different adjustment parameters were used, a = 0.075 (on the left) and a = 0.225 (on the right). The black circles represent the scaling ranges r = 0.075 and r = 0.225. Figure taken from [32].
Scaling indices investigate the spatial distribution of a previously prepared d-dimensional data set. In CMB investigations however, the fluctuations of the temperature maps are characterised by the values of the pixelised sky on a sphere. To be able to apply an analysis by means of scaling indices, one has to combine the temperature information with the two-dimensional spatial information of the map to create a three-dimensional point set, which includes all the information of the original map as spatial information only. This can be done by performing a preprocessing step, namely a transformation of the pixelised spherical sky S into three-dimensional space. Hereby, the pixels (θ i , φ i ), i = 1, ..., N pix , of S, where N pix denotes the number of pixels and (θ i , φ i ) latitude and longitude of the pixel i on the sphere, are converted into a three-dimensional jitter : Each temperature value T (θ i , φ i ) is assigned to one point p i , which is located in the radial direction through its pixel´s centre (θ i , φ i ), that is a straight line perpendicular to the surface of the sphere. Thus, the three-dimensional position vector of the new point p i reads as with R denoting the radius of the sphere and a describing an adjustment parameter. In addition, T and σ T characterise the mean and the standard deviation of the temperature fluctuations, respectively. The normalisation is performed to obtain for dR zero mean and a standard deviation of a. It is recommended to choose both R and a in a proper way to ensure a high sensitivity of the SIM with respect to the temperature fluctuations at a certain spatial scale. For CMB analysis, it turned out that this requirement is provided using R = 2 for the radius of the sphere and setting the adjustment parameter a to the value of the below introduced scaling range parameter r [41]. A CMB map transformed to a threedimensional point distribution is illustrated in Figure 1.
Hereby, two different adjustment parameters a were used in the embedding process. After this preprocessing step, the actual scaling index technique can be applied. In general, the SIM is a mapping that calculates for every point p i of the point set P a single value, which depends on the spatial position of p i in the group of the other points in which p i is embedded in. As already stated above, P is three-dimensional for the case of CMB analysis. For every point p i , the local weighted cumulative point distribution is defined as with r describing the scaling range, while s r (•) and d(•) denote a shaping function and a distance measure, respectively. To obtain our measure, the scaling index α( p i , r), we assume the following scaling law: The scaling index can therefore be computed as the logarithmic derivative of ρ( p i , r). Formally, this reads as In general, one is free to choose the shape of s r (•) and d(•). For the recent analyses that are discussed in this review, a set of quadratic gaussian shaping functions as well as the Euclidian norm were applied: Taking this into account, and using in addition the abbreviation d ij := d( p i , p j ), we obtain the final formula of FIG. 2: A simulated CMB map, in which the central regions were masked out and filled with nearly white noise, whereby the spatial noise patterns are preserved (see section III A) (upper left), and its scaling index responses α( pi, r) for three different scaling ranges: r = 0.05 (upper right), r = 0.15 (lower left) and r = 0.25 (lower right). Different values of α( pi, r) correspond to different types of structure in the underlying map. Small scaling ranges examine the behaviour of the small structures, while the characteristics of the larger structure is displayed by the higher scaling ranges. Note the different structure inside and outside the masked region of the simulated map, which is clearly identified by the scaling indices. These maps (and all following ones) are shown in a conventional scheme, namely the Mollweide projection in the Galactic reference frame with the Galactic Centre at the centre of the image and the longitude increasing from there to the left-hand side. the scaling indices: (1) In the resulting map α( p i , r), i = 1, ..., N pix , the structural behaviour of the underlying point set P becomes apparent, and different types of structure can be detected very easily. The values of α are related to structural characteristics in the following way: A point-or clusterlike structure leads to scaling indices α ≈ 0, filaments to α ≈ 1 and sheet-like structures to α ≈ 2. A uniform distribution of points would result in α ≈ 3. In between, curvy lines and curvy sheets produce 1 ≤ α ≤ 2 and 2 ≤ α ≤ 3, respectively. Underdense regions in the vicinity of point-like structures, filaments or walls feature α > 3. An example of a simulated CMB map and its scaling index response is shown in Figure 2. From equation (1), one can see that the scaling range parameter r can be chosen arbitrarily. This parameter weights the distances between our point of interest p i and the remaining points p j (see also definition of s r (x)). Therefore, we can make use of smaller or larger values for r to examine the different behaviour of the smallscale or large-scale structural configuration in the underlying map. For the analyses that were done by means of the scaling indices so far (see section IV), it was common to make use of the ten scaling range parameters r k = 0.025, 0.05, ..., 0.25, k = 1, 2, ...10. Figure 2 gives an example for the results of the SIM for three different values of r, applied onto a simulated CMB map.

III. DATA SETS
In the publications concerning scaling index analysis with CMB to date, different WMAP data sets as well as different techniques to handle foreground contaminated regions were applied. To test for the amount of non-Gaussianity in the WMAP data, simulations based on Gaussian random fields were constructed, which is the most common way in CMB analysis. In addition, socalled surrogate maps were generated. With the help of these surrogates, it is possible to test for the more specific hypothesis of uncorrelated phases. The surrogate method also offers the possibility to analyse the data in a scaledependent but model-independent way. In the following, we will give an overview of the used sets of maps.

FIG. 3:
The two plots on the left hand side illustrate the original 5-year WMAP data of the co-added VW-band (above) and the related colour-coded α-response (below) for a scaling range of r = 0.175. The Galactic Plane and several secondary spots are masked out by the KQ75-mask [45]. The equivalent plots for the mask-filling technique are arranged on the right hand side.

A. The WMAP data
In the current fourth data release of the WMAP team 1 [46], the recording of seven years of observation were made publicly available. One of the main challenges in all data releases as well as subsequent investigations is the handling of the heavily foreground contaminated regions, caused by different point sources and in particular the Galactic plane [45]. These approaches result in different maps, whereof the following ones were subject to an analysis by means of scaling indices: 1.) The ILC map The 7-years foreground-cleaned internal linear combination (ILC) map [45] is generated and provided by the WMAP team (in the following: ILC7). Measurements of all observed frequency bands are combined in terms of a linear combination. Different weights for the different bands as well as for different previously chosen fractions on the sky are determined to minimise the variance of the temperature fluctuations. The ILC7 map is downgraded to a resolution of 786432 pixels, which equals to N side = 256 in the employed HEALPix-software 2 [47].

2.) The Needlet-based ILC map
For comparison we also included the map produced in [48], namely the 5-years needlet-based ILC map, which has been shown to be significantly less contaminated by foreground and noise than other existing maps obtained from WMAP data (in the following: NILC5). The NILC5 map is downgraded to N side = 256 as well.
3.) The masked band-wise maps Unlike to the two ILC maps from above, the single Q-, V-and W-bands of the WMAP-satellite as well as a co-added VW-map can shed light on the influence of the different wavelenght depending foregrounds onto the CMB signal. Although we work with the maps that are reduced by means of the Foreground Template Model proposed in [49] and [50], these maps still show strong foreground effects and differ from each other, making a band-wise analysis reasonable. To obtain the band-wise or combined maps, we accumulate the differencing assemblies Q1, Q2, V1, V2, W1, W2, W3, W4 via a noiseweighted sum [51]: In this equation, A characterises the set of required assemblies, e.g.
for the co-added VW-map A = {V 1, V 2, W 1, W 2, W 3, W 4}, while the noise per observation of the different assemblies, given in [52], is denoted by σ 0 . The co-added VW-map is created for the 3-year and the 5-year data, while the single Q-, V-and W-band maps are generated for the 5-year data only.
As for the ILC maps above, we decrease the resolution to N side = 256. In addition, the heavily foregroundaffected parts of the sky are cut out, using the Kp0mask for the WMAP 3-year, and the KQ75-mask for the WMAP 5-year observations [53]. Hereby, both of them have to be downgraded as well. We choose a conservative downgrading of the mask by taking only all pixels at N side = 256 that consist completely of non-mask-pixels at N side = 512. All downgraded pixels at N side = 256, for which one or more pixels at N side = 512 belonged to the Kp0-or KQ75-mask, respectively, are considered to be part of the downgraded mask as well. In doing so, 23, 5% of the sky is removed for the Kp0-mask, while the KQ75-mask is even more conservative with 28.4% (see upper left part of Figure 3). Finally, we remove the residual monopole and dipole by means of the appropriate HEALPix routine applied to the unmasked pixels only.
4.) The mask-filled band-wise maps Just cutting out the masked regions like above spoils the results of the scaling index method: Instead of a more or less uniform distribution, the α-values in the regions around the mask now detect a sharp boundary with no points in the masked area, into which the scaling regions extend (see Figure 1). This results in lower values of α. The effect can clearly be seen in the α-response of the masked VW-band 5-year WMAP-data in the lower left corner of Figure 3. A solution to this problem is to fill the masked areas with suitable values, that prevent the low outcome at the edges of the mask. This can be accomplished by applying the following two steps: At first, we fill the masked regions with Gaussian noise, whose standard deviation for each pixel corresponds to the pixel noise made available by the WMAP-team: Here, σ (θ,φ) denotes the pixel noise of the pixel which is located in the direction (θ, φ). Then, we scale the expectation value and the variance as a whole to the empirical mean µ rem and variance σ 2 rem of the remaining regions of the original temperature map: where R and M stand for the non-masked and masked region of the map respectively, and N R as well as N M denote their number of pixels. Thus, we filled the mask with (nearly) white Gaussian noise whose mean and standard deviation equal the respective terms of the remaining map, whereby the spatial noise patterns are preserved. With this filling technique, we obtain a complemented data set instead of just excluding the masked regions. Boundary effects caused by the mask can be eliminated. The right column of Figure 3 shows the filling method as well as the corresponding α-response.

B. Simulations
A simple approach to evaluate the amount of non-Gaussianity in the WMAP data is to compare the measured data with maps that fulfil the Gaussian hypothesis. For the band-wise analysis, it is important to create simulations for each respective band. The proceeding is hereby as follows: We take the best fit ΛCDM power spectrum C l , derived from the respective WMAP 3-year or 5-year data only, and the according window function for each differencing assembly (Q1-Q2, V1-V2, W1-W4), as again made available on the LAMBDA-website. With these requisites, one can create Gaussian random fields mimicking the Gaussian properties of the best fit ΛCDM -model and including the WMAP-specific beam properties by convolving the C l´s with the window function. For every assembly, we add Gaussian noise to these maps with a particular variance for every pixel of the sphere. This variance depends on the number of observations N i (θ, φ) in the respective direction and the noise dispersion per observation, σ 0,i . After this procedure, the co-added VW-band (for the 3-year and 5year analyses) as well as the Q-, V-and W-bands (for the 5-year investigations only) can be summarised using equation (2) from above and decreased to the resolution of N side = 256. The respective Kp0-or KQ75-mask is cut out and the residual monopole and dipole removed, just as for the WMAP-data above. For comparison with the mask-filled data maps, the filling method described above is additionally applied onto the 5-year simulations as well.

C. Surrogates
A comparison with simulated CMB maps represents the most obvious and common approach to search for non-Gaussianities in the data set. However, it is also possible to create maps, so-called surrogates, that are similar to the original map except for one (or more) previously selected feature(s) which is (are) randomised. By comparing the data with this set of maps, one focuses on the deviations caused by the randomisation of these feature(s). The whole proceeding is therefore model-independent. One way of applying this method in terms of a scale-dependent search for non-Gaussianities has been proposed and discussed in [42] and [44]. In the following, we describe the various steps for generating surrogate maps in more detail: Consider a CMB map T (θ, φ), where T (θ, φ) is Gaussian distributed, and its Fourier transform. The complex valued Fourier coefficients a m can be written as a m = |a m |e iφ m with φ m = arctan (Im(a m )/Re(a m )). The linear or Gaussian properties of the underlying random field are contained in the absolute values |a m |, whereas all higher-order correlations (HOCs) -if present -are encoded in the phases φ m and the correlations among them. Having this in mind, a versatile approach for testing for scale dependent non-Gaussianities relies on a scale-dependent shuffling procedure of the phase correlations followed by a statistical comparison of the so-generated surrogate maps.
However, the Gaussian shape of the histogram of the temperature distribution and the randomness of the set of Fourier phases in the sense that they are uniformly distributed in the interval [−π, π], are a necessary prerequisite for the application of the surrogate-generating algorithm, which we propose in the following. To fulfil these two conditions, we perform the following preprocessing steps. First, the maps are remapped onto a Gaussian distribution in a rank-ordered way. This means that the amplitude distribution of the original temperature map in real space is replaced by a Gaussian distribution in a way that the rank-ordering is preserved, i.e. the lowest value of the original distribution is replaced with the lowest value of the Gaussian distribution etc. By applying this remapping we automatically focus on HOCs induced by the spatial correlations in the data while excluding any effects coming from deviations of the temperature distribution from a Gaussian one.
To ensure the randomness of the set of Fourier phases we performed a rank-ordered remapping of the phases onto a set of uniformly distributed ones followed by an inverse Fourier transformation. These two preprocessing steps only have marginal influence to the maps (see Figure 4). The main effect is that the outliers in the temperature distribution are removed. Due to the large number of temperature values (and phases) we did not find any significant dependence of the specific Gaussian (uniform) realisation used for remapping of the temperatures (phases). The resulting maps may already be considered as a surrogate map and we named it zeroth order surrogate map. The first and second order surrogate maps are obtained as follows: At first, one generates a first order surrogate map, in which any phase correlations for the scales, which are not of interest, are randomised. This is achieved by a random shuffle of the phases φ m for / ∈ ∆ = [ min , max ], 0 < m ≤ and by performing an inverse Fourier transformation.
In a second step, a chosen number of realisations of second order surrogate maps are generated for the first order surrogate map, in which the remaining phases φ m with ∈ ∆ , 0 < m ≤ are shuffled, while the already randomised phases for the scales, which are not under consideration, are preserved. Note that the Gaussian properties of the maps, which are given by |a m |, are exactly preserved in all surrogate maps.
In [42], the surrogates method was applied only to the -range ∆ = [2,20], while in [44], the analysis were extended to smaller scales as well: Three more -intervals, namely ∆ = [20,60], ∆ = [60, 120] and ∆ = [120, 300] were considered. The choice of 60 as min and max is somewhat arbitrary, whereas the min = 120 and max = 300 for the last -interval was selected in such a way that the first peak in the power spectrum is covered. Going to even higher does not make much sense, because the ILC7 map is smoothed to one degree FWHM. In principle, one could include higher , since some other maps -especially NILC5 -are not smoothed. But to allow for a consistent comparison of the results obtained with the different observed and simulated input maps, a restriction to only investigate -intervals up to max = 300 is applied. Besides this two-step procedure aiming at a dedicated scale-dependent search for non-Gaussianity, one can also test for non-Gaussianity using surrogate maps without specifying certain scales. In this case there are no scales, which are not of interest, and the first step in the surrogate map making procedure becomes dispensable. The zeroth order surrogate map is to be considered here as first order surrogate and the second order surrogates are generated by shuffling all phases with 0 < m ≤ for all available , i.e. in our case ∆ = [2,1024].
Finally, for calculating scaling indices to test for higher order correlations, the surrogate maps were degraded to N side = 256 and residual monopole and dipole contributions were subtracted. In contrast to a comparison between the data and simulated maps, which reveals all kinds of deviations from Gaussian random fields, the statistical comparison of the two classes of surrogates focusses on possible HOCs on certain scales, and the question if these have left traces in the first order surrogate maps, which were then deleted in the second order surrogates.

The WMAP three-year data
Finding differences between observed and simulated CMB maps which fulfil the Gaussian hypothesis of the best fitting ΛCDM model is a strong indication of the existence of non-Gaussianities in the CMB. The WMAP 3-year data was compared to N = 1000 simulations. The probability densities P (α) of the scaling indices for one selected scale (r = 0.175) are displayed in Figure 5 for the WMAP 3-year data and a subset of 20 simulations. The probability density for the WMAP data is shifted towards higher values, which indicates that the underlying temperature fluctuations for the observed data resemble more 'unstructured', that is, random and uniform fluctuations in comparison to the simulations. This effect is more pronounced in the northern hemisphere of the galactic coordinate system than in the southern. Furthermore, the histograms of the simulations are slightly broader, indicating that the simulations exhibit a larger structural variability than the observed data. These effects can more exactly be quantified by calculating the mean and standard deviation for the distribution of scaling indices as calculated for different scaling ranges. For scales larger than r = 0.1 the mean of the scaling indices is always systematically higher for WMAP than for the simulations. The effect is much more pronounced in the northern hemisphere. For the standard deviation we observe for the same scales significantly lower values for WMAP in the northern hemisphere and slightly higher ones for the southern sky. For the full sky these two effects cancel each other so that the observed deviations to lower values are no longer so significant. Besides the mean and standard deviation we additionally considered a combination of these two test statistics, namely a diagonal χ 2 -statistic at a given scale r k , where M 1 (r k ) = α(r k ) , M 2 (r k ) = σ α(r k ) . These statistics are computed for both the simulations and the observed moments. The σ-normalised deviations S of the WMAP data from the simulations (with Y = α(r k ) , σ α(r k ) , χ 2 ) are shown in Figure 6. The mean and standard deviation σ are obtained by summing over N = 1000 simulations. The percentages p of the simulations with higher (lower, respectively) results of the scale-independent diagonal χ 2 -statistics are calculated as well. The levels for the detection of non-Gaussianities (NGLs) are very high and do not fall below 99% for any scale. Even higher values for both the significances and the confidence levels are found, if one only considers the northern hemisphere. For scales larger than r = 0.15 none of the simulations was found to have a higher values for α than the observation. For the southern hemisphere, however, both the significances and the confidence levels for the smaller radii are slightly higher than for the northern sky but continuously decrease for higher radii r. For the standard deviation we find slightly different results. In a transition regime r ≈ 0.075 the width of the distribution of α is practically the same for the observation and the Monte Carlo sample. On smaller scales σ α is higher for WMAP, on larger scales r we observe higher standard deviations for the simulations. This effect is more pronounced in the northern hemisphere. For the largest scales the differences for σ between simulations and observation diminishes. Especially, for the southern hemisphere no signatures for deviations from Gaussianity are identified at larger scales using σ α . The behaviour of the χ 2 -statistics as a function of the scale parameters r can -as expected -be regarded as a superposition of the two underlying statistics < α > and σ α .
Some readers might argue that the selection of certain moments (mean, standard deviation, χ 2 ) and scales r k for highest significance, represents an a posteriori choice analysing the data. Although a choice might be well motivated by the results obtained with simulations, we are also using statistics that are a priori. In order to test for non-Gaussianity, we calculated scale-independent diagonal χ 2 statistics, where we considered only one (mean or standard deviation) or both measures, and summed over N r considered length scales r k , k = 1, ..., 10.
There is an ongoing discussion, whether a diagonal χ 2statistic or the ordinary χ 2 -statistic, which takes into account correlations among the different random variables through the covariance matrix is the better suited measure. On the one hand it is important to take into account correlations among the test statistics, on the other hand it has been argued by [11] that the calculation of the inverse covariance matrix may become numerically unstable when the correlations among the variables are strong making the ordinary χ 2 statistic sensitive to fluctuations rather than to absolute deviations. For the WMAP 3year data we follow the reasoning of [11] and choose a diagonal χ 2 -statistics, because also in our case the moments are highly correlated leading to high values in the off-diagonal elements of the cross-correlation matrix. However, if the chosen model is a proper description of the data, any combination of measures should yield statistically the same values for the observations and the simulations. Also for the a priori scale-independent test statistics, where some unimportant scales contribute to the final value of χ 2 , we find significant signatures for non-Gaussianities in the northern sky. We detect non-Gaussianity for the full sky at a level of 96.9% regarding the mean, 96.5% for the standard deviation and 97.3% for a combination of mean and standard deviation. For the northern hemisphere, the signatures of non-Gaussianities are more pronounced and we obtain 97.7% (mean), 99.5% (standard deviation) and 98.9% (combination), whereas the Southern hemisphere is more consistent with Gaussianity [94.2% (mean), 70.0% (standard deviation) and 91.6% (combination)]. These differences between Northern and Southern hemispheres induce pronounced asymmetries, which can be interpreted as a global lack of structure in the Northern hemisphere, which is consistent with previous findings (see below).
If we select only those pixels for P (α) which have |b| > 30 (b galactic latitude), well outside the galactic plane, we get higher values for P (α). The disturbing edge effects of the Kp0-mask are almost totally removed. Only now we detect a localised anomaly in the Southern hemisphere when analysing the α spectra. We identify this signature as the Cold Spot, which was also already detected in the first-year WMAP data by [18]. More local features are identified with the 5-year data and will be discussed in more detail in the following section. The probability densities P (α) of the selected pixels well outside the galactic plane are very similar to the former ones (see Figure 5). The same holds for the significances for non-Gaussianity.

The WMAP five-year data
The empirical probability densities P (α) of the scaling indices (calculated here with r = 0.2) for WMAP 5-year data and respective simulations in Figure 7 show again a shift of the WMAP data to higher values, which becomes particularly apparent in the northern hemisphere of the galactic coordinate system. Comparing the non-filling and the mask filling method (see section III A), the histograms of the latter feature a higher maximum as well as higher values for large α, but lower probabilities for α ∈ [2.0, 2.5]. The obvious reason for this shift is the fact that the filled mask does not reduce the α-values of its surroundings as it was the case with the former method. Now, the outcome of these regions is influenced by the white noise and is therefore allocated at higher values.
We also calculated the σ-normalised deviations S and the percentages p of the simulations with higher (lower, respectively) results of the scale-independent diagonal χ 2 -statistics, again for N = 1000 simulations. High deviations are found, particularly in the northern hemisphere. We derive evidence for non-Gaussianity with a probability of up to 97.3% for the mean when regarding the KQ75-masked full sky and summing up over all considered length scales by means of a diagonal χ 2 -statistics. Looking at only the northern or southern hemisphere of the galactic coordinate system, we obtain up to 98.5% or 96.6%, respectively. For the standard deviation, these results appear as 95.6% for the full sky (99.7% north, 89.4% south) and for a χ 2 -combination of both measurements as 97.4% (99.1% north, 95.5% south). We obtain larger deviations from Gaussianity when looking at separate scale lengths. In general, all occurring characteristics match the findings of the analysis of the WMAP 3year data. This indicates that the results are not based on some time-dependent effects. Since the 5-year data features lower error bars than the 3-year data, it is also improbable that both results are induced by noise effects only.
Evidence for north-south asymmetry in the WMAP data was already detected using the angular power spectrum [54,55] and higher order correlation functions [29] spherical wavelets [18], local curvature analysis [56], twodimensional genus measurements [57] as well as all three Minkowski functionals [11], correlated component analysis [58], spherical needlets [21], frequentist analysis of the bispectrum [59], two-point correlation functions [60,61] and Bayesian analysis of the dipole modulated signal model [62]. To take a closer look at asymmetries in the WMAP 5-year data in our investigations, we perform an analysis of rotated hemispheres and detect an obvious asymmetry in the data: For each scale we calculate the mean α(r k ) and standard deviation σ α(r k ) of the map of scaling indices α( p i , r) (or α(θ, φ; r k )) of the full sky and a set of 3072 rotated hemispheres. Every pixel centre of the full sky consisting of 3072 pixels (N side = 16) in the HEALpix [47] pixelisation scheme marks a new northern pole of one of the different hemispheres. The differences between the results of the original as well as the simulated maps are again quantified by the σ-normalised deviation S (see equation 4). Every hemisphere of the set of 3072 hemispheres delivers one deviation value S, which is then plotted on a sky map at the northern pole of the respective hemisphere.
Thus, the colour of each pixel in the corresponding Figure 8 expresses the positive or negative σ-normalised deviation S(r) of the hemisphere around that pixel in the WMAP-data compared to the hemispheres around that pixel in the simulations. We apply this analysis for the co-added VW-band as well as for the single bands Q, V and W, whereas for the VW-band we use both the original and the mask filling method, but for the single bands the filling method only. In all charts of Figure 8 we can detect an obvious asymmetry in the data: The largest deviations between the data and the simulations are exclusively obtained for rotations pointing to northern directions relative to the galactic coordinate system. The maximum value for S(r) of the χ 2 analysis (right column of Figure 8) using the mask-filling method on the co-added VW-band is obtained in the reference frame pointing to (θ, φ) = (27 • , 35 • ), which is close to the galactic north pole. This proximity to the pole is consistent to the results of [56] and [41], as well as to those findings of [54] and [29] that consider large angular scales. For the standard deviation (central column of Figure 8), the northern and southern hemispheres offer different algebraic signs. The negative S(r) of the north implies a lower variability than the simulations in this region, while the south shows a converse behaviour. The fact that the plots using the new method show slightly lower values for S(r) than the ones using the old method may be explained by the fraction of pure noise values within every rotated hemisphere, that diminish the degree of difference between the data and the simulations.
While the Q-band is heavily foreground-affected, first of all by synchrotron radiation as well as radiation from electron-ion scattering ("free-free emission"), the Wband is mainly distorted by dust emission. The V-band is affected by all three of these foregrounds, even though less than the other bands. Despite the different influences on the different bands, we obtain the same signatures of non-Gaussianity in all single bands as well as in the co-added VW-band. The correlations c of the different bands are high (c ≥ 0.95). Therefore, we conclude that the measured asymmetry is not the result of a foreground influence but has to be concluded of thermal origin.
An interesting anomaly in the CMB data is that there are small regions which show very high or very low values in some local structure analysis. One of the first of these local features, the well-known Cold Spot at (θ, φ) = (147 • , 209 • ), was first detected by [18] in 2004 by using a wavelet analysis. Scaling indices were able to redetect the Cold Spot in the WMAP 3-year data (see section IV A 1). Furthermore, it was identified using using amongst others wavelet analysis [63][64][65][66] or the Kolmogorov stochasticity parameter [23]. Furthermore, there have been some investigations which, in addition to the re-detection of the first spot, detected secondary spots via directional [67][68][69] or steerable wavelets [70], needlets [21] and again the Kolmogorov stochasticity parameter [24]. These spots could be the result of some yet not fully understood physical process. For the Cold Spot lots of theories already exist which try to explain its origin by second-order gravitational effects [71,72], a finite universe model [73], large dust-filled voids [74][75][76][77], cosmic textures [66], non-Gaussian modulation [78], topological defects [79], textures in a brane world model [80] or an asymptotically flat Lemaître-Tolman-Bondi model [81,82].
For our investigations concerning spots in the WMAP data we use the mask-filling method of section III A. Boundary effects caused by the mask are eliminated, which allows hidden effects to emerge. We extend the analysis of scaling indices by applying two different approaches to detect anomalies: The first one is to calculate the σ-normalised deviation of every pixel on the α-response of the CMB map. For a given scale parameter r, this is achieved by comparing the scaling index α( p i , r) of each vector p i , i = 1, ..., N pix , of the original data with the mean of the corresponding values α ( p i , r), = 1, ..., N sim , of the simulations depending on their standard deviation, where N sim denotes the number of the simulations. Formally, this reads as: with The results are illustrated in the upper left part of Figure  9.
The second approach smoothes the α-maps of the original and simulated data by computing for every pixel the mean value of its surroundings given by some specified maximum distance, which equals 3 • in our analysis. We apply the pixel-wise deviations S i,r again on the resulting maps. The outcome of this procedure is shown in the upper right part of Figure 9. In the lower left plot of the same Figure only the deviations S i,r ≤ −3.0 are illustrated to gain yet another clearer view on the interesting areas. We identify several local features on the map.
The first approach clearly shows the Cold Spot and indicates some secondary spots in the southern as well as in the northern hemisphere. These get confirmed in the plot of the smoothing method, where we obtain a deviation of up to −7σ for several clearly visible areas: In the southern hemisphere we detect a cold spot at (θ, φ) = (124 • , 320 • ) and another one at (θ, φ) = (124 • , 78 • ). Both were already detected with the above mentioned directional and steerable wavelet as well as with a needlet analysis. The former one is a hot spot in these investigations. In our analysis, the latter spot actually appears as two spots close to each other, which is in agreement with [21]. We discover another southern cold spot at (θ, φ) = (120 • , 155 • ) which is very close to the mask. This spot represents a good example for the use of the mask filling method since it is located at the edge of the non-masked region: The influence of the mask is diminishing the results of the calculation of the scaling indices in the area of this spot. This becomes obvious if one recalls the lower left plot of Figure 3, in which the coordinates of the spot would be completely located in a "blue" region with low α-values. Since the results of the scaling indices of local features show a similar, namely lower-valued, behaviour, an overlapping like that could prevent the detection of such spots close to the mask. By using the mask filling method, the detection of this cold spot on the edge of the mask is equivalent to a detection in an unmasked region, and therefore reliable. The spot at (θ, φ) = (136 • , 173 • ), described by [67] and [21], is not traced in our analysis. In the northern hemisphere, our investigation shows two other cold spots at (θ, φ) = (49 • , 245 • ) and (θ, φ) = (68 • , 204 • ), which do not correspond with the so-called northern cold spot of [24], but with the results of [67], where again one of them is a hot spot. Also [21] locates one of these two spots. All these results were achieved with an analysis of the VW-band, but we find similar results in a single band analysis.
If the considered spots really depend on some yet not completely understood, maybe secondary, physical effect, they should not be implemented in a testing for intrinsic non-Gaussianity. For this reason, we modify the 5-year KQ75-mask by additionally excluding all above mentioned spots. A small peculiarity at the edge of the mask next to the Cold Spot as well as three very small blurs in the right half of the lower left Mollweide projection in Figure 9 are not considered, since we regard their appearance as insufficient for being a distinctive feature. The modification of the KQ75-mask is illustrated in the lower right part of Figure 9.
We now apply this new mask to the α-response of both the WMAP data as well as the simulations and repeat the analysis from above. When excluding all these spots from the analysis, the deviation from Gaussianity increases, which shows that the discovered local anomalies are not the reason of the global detection of non-Gaussianity, but actually were damping the deviations on average. The results of the σ-normalised deviations S are illustrated in Figure 10. An increase of S(r) in comparison to the former analysis where the usual KQ75-mask was used (results not shown, please refer to [41] ) is in particular present in the southern hemisphere, where we detected more local features than in the north. The largest increase takes place in the co-added VW-band, where we now reach deviations of up to 4.0 for the χ 2 -combination in a full-sky analysis (former maximum: 2.9) and to the extend of 6.0 in an analysis of the northern hemisphere (former maximum: 5.5). Also the single bands as well as all scale-independent diagonal χ 2 -statistics show without exception a greater evidence for non-Gaussianity.

B. Comparison with Surrogates
We compare the first and second order surrogate maps by calculating the σ-normalised deviations S (similar to equation 4) between the two classes of surrogates for a set of (now) 768 hemispheres to test for NGs and asymmetries in the ILC7 map and the NILC5 map. Figure 11 shows the deviations S per hemisphere for the mean value S( α(r k ) ), k = 2, 6, 10 for the ILC7 map as derived from the comparison of the different classes of surrogates for the scale-independent surrogate test and for the four selected -ranges. The following striking features become immediately obvious: First, various deviations representing features of non-Gaussianity and asymmetries can be found in the S-maps for the ILC7 map. These features can nearly exactly be reproduced when the NILC5 map is taken as input map, whose results are illustrated in Figure 12. Second, we find for the scale-independent surrogate test (top rows in Figures 11 and 12) large isotropic deviations for the scaling indices calculated for the small scale r 2 . The negative values for S indicate that the mean of the scaling indices for the first order surrogate is smaller than for the second order surrogates. This systematic trend can be interpreted such that there is more structure detected in the first order surrogate than in the second order surrogate maps. Obviously, the random shuffle of all phases has destroyed a significant amount of structural information at small scales in the maps. Third, for the scale-dependent analysis we obtain for the largest scales (∆ = [2,20]) (second lines in Figures 11 and 12) highly significant signatures for non-Gaussianities and ecliptic hemispherical asymmetries at the largest r−values. These results are perfectly consistent with those obtained for the WMAP 5-year ILC map and the foreground removed maps generated by [25] on the basis of the WMAP 3-year data (see [42]). The only difference between this study and our previous one is that we now obtain higher absolute values for S ranging now from −4.00 < S < 3.72 for the ILC7 map and −4.36 < S < 4.50 for the NILC5 map as compared to −3.87 < S < 3.51 for the WMAP 5-year ILC map. Thus, the cleaner the map becomes due to better signal-to-noise ratio and/or improved map making techniques the higher the significances of the detected anomalies, which suggests that the signal is of intrinsic CMB origin. Fourth, we also find for the smallest considered scales (∆ = [150, 300]) large isotropic deviations for the scaling indices calculated for a small scaling range r very similar to those observed for the scale-independent test. Fifth, we do not observe very significant anomalies for the two other bands (∆ = [20,60] and ∆ = [60,120]) being considered in this study. Thus, the results obtained for the scale independent surrogate test can be interpreted as a superposition of the signals identified in the twobands covering the largest (∆ = [2,20]) and smallest ∆ = [120, 300]) scales. Figure 13 shows the probability densities derived for the full sky and for (rotated) hemispheres for the scaling indices at the largest scaling range r 10 for the first and second order surrogates for the -interval ∆ = [2,20]. We recognise the systematic shift of the whole density distribution towards higher values for the upper hemisphere and to lower values for the lower hemisphere. As these two effects cancel each other for the full sky, we do no longer see significant differences in the probability densities in this case. Since the densities as a whole are shifted, the significant differences between first and second order surrogates found for the moments cannot be attributed to some salient localizable features leading to an excess (e.g. second peak) at very low or high values in otherwise very similar P (α)-densities. Rather, the shift to higher (lower) values for the upper (lower) hemisphere must be interpreted as a global trend indicating that the first order surrogate map has less (more) structure than the respective set of second order surrogates. The seemingly counterintuitive result for the upper hemisphere is on the other hand consistent with a linear hemispherical structure analysis by means of a power spectrum analysis, where also a lack of power in the northern hemisphere and thus a pronounced hemispherical asymmetry was detected [55,56]. However, it has to be emphasised that the effects contained in the power spectrum are -by construction -exactly preserved in both classes of surrogates, so that the scaling indices measure effects that can solely be induced by HOCs thus being of a new, namely non-Gaussian, nature. Interestingly though, the linear and nonlinear hemispherical asymmetries seem to be correlated with each other. The density distributions derived from the ILC7 and NILC5 map are clearly shifted against each other. The differences between these two maps can be attributed to e.g. the smoothing of the ILC7 map. However, the systematic differences between first and second order surrogates induced by the phase manipulations prevail in all cases -irrespective of the input map. The results for the deviations |S(r)| for the full sky and rotated upper and lower hemisphere are shown for all considered -ranges and all scales r in Figures 14 and 15. Using scale-independent χ 2 α ,σα -statistics combining the mean and the standard deviation and summing up over all scales r, the largest values for S are found for the largest ∆ = [2,20] and smallest scales ∆ = [120, 300] and for the scale-independent NGs. For the full sky non-Gaussianity is detected with a probability of up to 94.2% (χ 2 α ,σα ) and 99.8% for the northern and southern hemispheres.
To test whether all these signatures are of intrinsic cosmic origin or more likely due to foregrounds or systematics induced by e.g. asymmetric beams or map making, we performed the same surrogate and scaling indices analysis for five additional maps described in [44]. This set of tests to investigate whether and to what extend the detected anomalies can be explained by systematics cannot convincingly rule out the intrinsic nature of the anomalies for the low case, while the ILC map making procedure and/or residual noise in the maps can also lead to NGs at small scales.

V. CONCLUSIONS
In this Review, we gave an overview of the application of scaling indices in CMB analysis to date. The SIM is a measure that detects different forms of topological be-FIG. 13: The probability densities P (α) of the scaling indices for the first (black) and second order surrogates (coloured) of WMAP 7-year data, calculated for the largest scaling range r10 and for the −interval ∆ = [2,20]. Yellow (green) curves denote the densities for 20 realizations of second order surrogates derived from the ILC7 (NILC5) map. The reference frame for defining the upper and lower hemispheres is chosen such that the difference ∆S = Sup − S low becomes maximal for α of the respective map and respective scale r. haviour in the data, which turned out to be very useful for identifying deviations from Gaussianity and statistical isotropy in the spherical data set of the microwave background radiation. In the following, the large number of deviations from Gaussianity and statistical isotropy detected by means of the SIM is summarised and the resulting conclusions are drawn.
By comparing the 3-year and 5-year measurements of the WMAP satellite with simulated CMB maps, several clear non-Gaussianities as well as asymmetries were detected: The spectrum of scaling indices of the data is systematically broader and shifted towards higher values than the one of the simulations, yielding highly significant deviations of the mean, the standard deviation and a χ 2combination. These effects can be interpreted as too few structure and structural variations in the temperature anisotropies as measured by WMAP compared to the predicted ones within the concordance model, which is in agreement with previous results (e.g. [11,29,55,62,83]). By performing an analysis of rotated hemispheres, the rotations pointing to northern directions show by far higher FIG. 14: Deviations |S(r)| for the ILC7 map and the considered −intervals as a function of the scale parameter r for the full sky (black) and the upper (red) and lower (blue) hemisphere. The '+' symbols denote the results for the mean α(r k ) , the' * ' symbols for the standard deviation σ α(r k ) and the boxes for the χ 2 -combination of α(r k ) and σ α(r k ) .
FIG. 15: Same as Figure 14, but for the NILC5 map. deviations from Gaussianity for the mean and the χ 2 analysis than rotations pointing to the south. For the standard deviation, the rotated hemispheres show a negative outcome in the north and a positive in the south. This implies that the north possesses a more consistent pattern than the simulations, while the south shows a converse behaviour.
All these results are consistent in different ways: Since the detected effects are the same for the 3-year as well as for the 5-year WMAP data, they can be concluded to be time-independent. In addition, the findings are nearly the same for the different bands that were analysed for the 5-year data, which leads to the conclusion that the foreground influence only plays a minor role. Furthermore, the usage of the mask-filling method, again applied to the 5-year data only, reduces the distorting influence of the mask. Since this leads to similar results as well, the detected deviations from Gaussianity and statistical isotropy must therefore be taken to be of cosmological origin so far.
In addition to these findings, several local features including the Cold Spot could be detected with the scaling indices, which turns out to be another advantage of this method. The fact that most of them are located in the southern hemisphere confirms the conclusions concerning the asymmetries from above. Nearly all detected spots are in agreement with former analyses (e.g. [18,21,67]), which confirms the existence of these local anomalies.
By accomplishing a comparison of the CMB data with surrogate maps, one focuses on the more specific assumption of random and uncorrelated phases, which is part of the Gaussian hypothesis. In addition, this method offers the possibility of a scale-dependent analysis. The scaling indices are the first measure which is used in combination with this surrogates approach. For an analysis of the 5-and 7-year observations of the WMAP satellite, the results are as follows: Highly significant non-Gaussianities could be detected, again by performing an analysis of rotated hemispheres, for the very large scales and for the -interval covering the first peak in the power spectrum. The results show the most significant evidence of non-Gaussianity in the CMB to date, and disagree strongly with predictions of isotropic cosmologies in single field slow roll inflation. Several checks on systematics were performed, which lead to the conclusion that the findings are of cosmological origin.
For smaller scales (i.e. higher -ranges), it turns out that phase correlations can be easily induced by the ILC map making procedure, so that it is difficult to disentangle possible intrinsic anomalies from effects induced by the preprocessing of the data. In this case, more tests are required to further determine of these high-anomalies.
The SIM is the only measure in CMB analysis that was used in combination with the surrogates technique so far. Further studies that combine the surrogates method with different measures, as for example Minkowski functionals, could support these investigations and produce even more reliable results. In addition, the upcoming data of the PLANCK satellite offers an independent measurement of the CMB, and will allow investigations concerning higher -bands.