text:
# Morphological Measures of nonGaussianity in CMB Maps

## 1 Introduction

The issue of Gaussianity of CMB maps plays a crucial role in testing assumptions about the early Universe. The simplest inflation models strongly favor Gaussianity of the primordial inhomogeneities (see Turner, 1997, and references therein), whereas other scenarios assuming cosmic strings or topological defects (see Vilenkin and Shellard, 1994, and references therein) predict non-Gaussian perturbations. Gaussianity is also a key underlying assumption of all experimental power spectrum analyses to date, entering into the computation of error bars (Tegmark, 1997; Bond and Jaffe, 1998), and therefore needs to be observationally tested. A third reason for studying Gaussianity of CMB maps is that it may reveal otherwise undetected foreground contamination.

Numerous tests of Gaussianity in the COBE maps (Colley, Gott, and Park, 1996; Kogut et al., 1996; Ferreira, Magueijo, and Górski, 1998; Pando, Valls-Gabaud, and Fang, 1998; Bromley and Tegmark, 1999; Novikov, Feldman, and Shandarin, 1999; Banday, Zaroubi, and Górski, 2000; Barreiro et al., 2000; Magueijo, 2000; Mukherjee, Hobson, and Lasenby, 2000; Aghanim, Forni, and Bouchet, 2001; Phillips and Kogut, 2001) have resulted in the general agreement that all non-Gaussian signals were of noncosmological origin. This was not unexpected because of COBE’s low ($`7^{}`$ deg) angular resolution.

Testing Gaussianity on smaller scales may bring much more interesting results. The first study of Gaussianity on degree scale showed the consistency of the QMASK map with the assumption of the Gaussianity (Park et al., 2001). However, this study tested only the Gaussianity of the topology of the map. The study of the MAXIMA-1 anisotropy data also showed the consistency with the Gaussianity in the range between 10 arc-minutes and 5 degrees (Wu et al., 2001). A total of 82 tests for Gaussianity were made in this study. Although performing many tests is always better than few it is not clear how independent these tests were.

The question of statistical independence of different non-Gaussianity tests is complicated. However, in some simple cases some of the tests can be statistically disentangled. For instance, consider the probability function, which is perhaps the simplest test for non-Gaussianity. It easy to transform the probability function to any given form simply by relabeling the contours using a monotonic function. This relabeling obviously does not have any effect on the morphology of the field since every contour line remains the same as well as the order of levels due to monotonicity of the transformation. However, this transformation may strongly affect the whole hierarchy of the $`n`$-point correlation functions. This is of course due to well known connection of the probability function to all $`n`$-point functions.

Genus represents only one statistic sensitive to non-Gaussianity. It may detect some types of non-Gaussianity and be not sensitive to the others. Consider, for instance, the following mapping of a Gaussian field $`\varphi (x,y)`$ into $`\mathrm{\Phi }(X,Y)`$

$$\mathrm{\Phi }(X,Y)=\varphi (x,y),$$
(1)
where

$`X`$ $`=`$ $`x+\tau u(x,y),`$
$`Y`$ $`=`$ $`y+\tau v(x,y).`$ (2)
Equation 1 guarantees that every level label remains unchanged, while the level lines are shifted and deformed according to eqs. 2. Functions $`u(x,y)`$ and $`v(x,y)`$ can be any smooth functions including random functions. At small $`\tau `$ before caustics have occurred, the mapping $`(x,y)(X,Y)`$ is single valued and therefore the contour map of $`\mathrm{\Phi }(X,Y)`$ is one-to-one image of $`\varphi (x,y)`$, although geometrically distorted. The genus of the $`\mathrm{\Phi }(X,Y)`$ field as a function of the level remains exactly the same (i.e. Gaussian) as that of the $`\varphi (x,y)`$ field because the mapping is continous and nonsingular and therfore preserves toplogy. On the other hand, the mapping $`(x,y)(X,Y)`$ can be neither isometric (conserving the lengths) nor area-preserving and therefore both the areas and the contours of the excursion sets may be distorted. As a result, the $`\mathrm{\Phi }(X,Y)`$ can be a strongly non-Gausssian field having exactly the Gaussian genus in the level parameterization. Perhaps, it is worth noting that a particular case of the above example ($`u(x,y)=\psi (x,y)/x`$ and $`v(x,y)=\psi (x,y)/y`$) corresponds to viewing the Gaussin contour map $`\varphi (x,y)`$ through a glass layer with varying thickness $`\psi (x,y)`$ (see e.g. Zel’dovich, Mamaev, and Shandarin (1983); Shandarin and Zel’dovich (1989)).

A one-dimensional illustration of the above field is shown in Fig. 1. A realization of the Gaussian field $`\varphi (x)`$ is shown by the solid line and of the non-Gaussian field $`\mathrm{\Phi }(x)`$ by the dashed line. The total area of the excursion set in 2D is analogous to the total length of the excursion set in 1D and the genus in 2D is analogous to the number of peaks in 1D. The level parameterization means the comparison of the number of peaks in the two fields at the same level. Both fields have the same number of peaks at every level as illustrated by the dotted and dash-dotted horizontal lines. However, the total lengths occupied by the peaks of the $`\varphi (x)`$ and $`\mathrm{\Phi }(x)`$ fields are different at a given level. The length parameterization in 1D means that the comparison of the numbers of peaks must be done at the same total length of the excursion sets which generally happens at different levels for the $`\varphi (x)`$ and $`\mathrm{\Phi }(x)`$ fields. For instant, the length of the excursion set of the Gaussian field at $`\varphi (x)=1.95`$ (the sum of two heavy line segments at level $`=1.95`$) equals the length of the excursion set of the non-Gaussian field at $`\mathrm{\Phi }(x)=1`$ (the sum of four heavy line segments at level $`=1`$) . One can see at this total length the number of peaks in the Gaussian field is two while in the non-Gaussian field it is four. Thus, the number of peaks parameterized by the total length does detect this type of non-Gaussianity while the number of peaks parameterized by the level does not. Similarly, the area parameterization of the genus in two dimensions does detect this type of non-Gaussianity while the level parameterization of the genus statistic does not.

If the mapping (eq. 2) is area-preserving than the genus of the $`\mathrm{\Phi }(X,Y)`$ field remains exactly Gaussian also in the area parameterization. Thus, a detection of such non-Gaussianity requires additional tests.

We use two different parameterization of the measured morphological characteristics and show that the parameterization by the total area of the excursion set (i.e. the cumulative probability function) has an advantage over parameterization by the temperature because it gives considerably smaller correlations between different measures. In particular, the total area parameterization detects non-Gaussianity by the genus statistic in the above example while the temperature parameterization does not. For additional discussion of this issue, see Shandarin (2002).

Along with the QMASK map (Xu et al. 2001), we analyze a thousand reference maps with the same sky coverage, noise properties and power spectrum as the QMASK map. For each map we compute seven morphological functions at a hundred temperature/area levels. These functions are:

1. the total area of the excursion set $`A(\mathrm{\Delta }T)`$,
2. the total contour length $`C(\mathrm{\Delta }T)`$ and $`C(A)`$ ,
3. the genus of the excursion set $`G(\mathrm{\Delta }T)`$ and $`G(A)`$ ,
4. the area of the largest (by area) region $`A_p(\mathrm{\Delta }T)`$ and $`A_p(A)`$ ,
5. its perimeter $`C_p(\mathrm{\Delta }T)`$ and $`C_p(A)`$,
6. genus $`G_p(\mathrm{\Delta }T)`$ and $`G_p(A)`$, and finally
7. the number of regions $`N_c(\mathrm{\Delta }T)`$ and $`N_c(A)`$ .

The first three are the global Minkowski functionals and the following three are the Minkowski functionals of the largest (by area) region. The largest region at some threshold spans the whole map, or “percolates”. Its Minkowski functionals are the best indicators of the percolation phase transition (Shandarin, 1983; Yess and Shandarin, 1996). Our notations reflect some of the jargon used in cluster analysis: $`A_p,C_p,G_p`$ for the Minkowski functionals of the percolating region and $`N_c`$ for the number of regions, often referred to as “clusters” in cluster analysis.

We then calculate the mean for each of these fourteen functions ($`\overline{M_i}(\mathrm{\Delta }T)`$ and $`\overline{M_i}(A),i=1,\mathrm{},7`$) computed from a thousand Gaussian maps and quantify the differences of every function for every map with respect to the mean functions in quadratic measure. Finally, we compute these differences for the QMASK map and test the Gaussianity hypothesis by comparing them to the thousand reference maps.

The rest of the paper is organized as follows. In Sec. 2 and 3 we briefly describe the QMASK map and the Gaussian simulations thereof. In Sec. 4 we define the morphological functions, their parameterization and the method of quantifying the departure from the Gaussianity. In Sec. 5 we summarize the results.

## 2 QMASK Map

The QMASK map, described in (Xu et al., 2001), combines all the information from the QMAP (Devlin et al., 1998; Herbig et al., 1998; de Oliveira-Costa et al., 1998) and Saskatoon (Netterfield et al., 1995, 1997; Tegmark et al., 1996) experiments into a single map at 30-40 GHz covering about 648 square degrees around the North Celestial Pole. The map consists of sky temperatures in 6495 sky pixels, conveniently grouped into a 6495-dimensional vector $`𝐱`$, with a FWHM angular resolution of $`0.68^{}`$. As detailed in Xu et al. (2001), all the complications of the map making and deconvolution process are encoded in the corresponding $`6495\times 6495`$ noise covariance matrix $`𝐍`$. The map has a vanishing expectation value $`𝐱=\mathrm{𝟎}`$ and a covariance matrix given by

$$𝐂\mathrm{𝐱𝐱}^t=𝐍+𝐒,$$
(3)
where the signal covariance matrix $`𝐒`$ is given by

$$𝐒_{ij}=\underset{\mathrm{}=2}{\overset{\mathrm{}}{}}\frac{\mathrm{}(\mathrm{}+1)}{4\pi }P_{\mathrm{}}(\widehat{𝐫}_i\widehat{𝐫}_j)e^{\theta ^2\mathrm{}(\mathrm{}+1)}C_{\mathrm{}},$$
(4)
the $`P_{\mathrm{}}`$ are Legendre polynomials, $`C_{\mathrm{}}`$ is the angular power spectrum, $`\widehat{𝐫}_i`$ is a unit vector pointing towards the $`i^{\mathrm{th}}`$ pixel, and $`\theta =`$FWHM$`/\sqrt{8\mathrm{ln}2}`$ is the rms beam width. When computing $`𝐒`$ in practice, we use the smooth power spectrum of (Xu, Tegmark, and de Oliveira-Costa, 2001) that fits the observed QMASK power spectrum measurements.

Since the raw map has a large and complicated noise component, we work with the Wiener-filtered version of the map in this paper, shown in Fig. 2 and defined as

$$𝐱_w\mathrm{𝐖𝐱},𝐖𝐒(𝐍+𝐒)^1.$$
(5)
This approach was also taken in Park et al (2001) and Wu et al (2001).

## 3 Mock Maps

In order to quantify the statistical properties of our Minkowski functionals, we need large numbers of simulated QMASK maps. We therefore generate one thousand mock Gaussian maps as follows. The covariance matrix $`𝐂_w`$ of the Wiener-filtered map is given by

$$𝐂_w𝐱_wx_w^t=(\mathrm{𝐖𝐱})(\mathrm{𝐖𝐱})^t=𝐖\mathrm{𝐱𝐱}^t𝐖^t=\mathrm{𝐖𝐂𝐖}^t,$$
(6)
and we Cholesky-decompose it as $`𝐂_w=\mathrm{𝐋𝐋}^t`$, where the matrix $`𝐋`$ lower triangular. It is straightforward to show that if $`𝐳`$ is a vector of independent Gaussian random variables with zero mean and unit variance (which is trivial to generate numerically), i.e., $`\mathrm{𝐳𝐳}^t=𝐈`$, the identity matrix, then the mock maps defined $`𝐱_m\mathrm{𝐋𝐳}`$ will have a multivariate Gaussian distribution with $`𝐱_m𝐱_m^t=𝐂_w`$.

## 4 Morphological Statistics

In this paper we use a set of morphological statistics based on Minkowski functionals. As suggested by (Novikov, Feldman, and Shandarin, 1999), one can use both global and partial Minkowski functionals in studies of CMB maps. In this study we use the global Minkowski functionals and those of the largest region by area (the so-called percolating region).

### 4.1 Global Minkowski functionals

Global Minkowski functionals were introduced into cosmology as quantitative measures of CMB anisotropy by Gott et al. (1990) although without reference to the Minkowski functionals and thus without stressing their unique role in differential and integral geometry. Mecke, Buchert, and Wagner (1994) were the first to place the studies of morphology of the large scale structure in the context of differential and integral geometry. In particular, they emphasized a powerful theorem by (Hadwiger, 1957) stating that with rather broad restrictions, the set of Minkowski functionals provides a complete description of the morphology (for further discussion see Kerscher, 1999).

Global Minkowski functionals describe the morphology of the excursion set at a chosen threshold: $`A(\mathrm{\Delta }T)`$, $`C(\mathrm{\Delta }T)`$ and $`G(\mathrm{\Delta }T)`$ are the total area, contour length and genus of the excursion set, respectively. The Minkowski functionals of Gaussian fields are known analytically as functions of the threshold. Assuming that the field is normalized, i.e. $`u=(\mathrm{\Delta }T<\mathrm{\Delta }T>)/\sigma _{\mathrm{\Delta }T}`$ so that $`<u>=0,<u^2>=1`$, the Minkowski functionals are

$`A(u)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[1\mathrm{erf}\left({\displaystyle \frac{u}{\sqrt{2}}}\right)\right],`$ (7)
$`C(u)`$ $`=`$ $`{\displaystyle \frac{1}{2R}}\mathrm{exp}\left({\displaystyle \frac{u^2}{2}}\right),`$ (8)
$`G(u)`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^{\frac{3}{2}}}}{\displaystyle \frac{1}{R^2}}u\mathrm{exp}\left({\displaystyle \frac{u^2}{2}}\right),`$ (9)
where $`R=\sqrt{2}/\sigma _1`$ is the characteristic scale of fluctuations in the field; $`\sigma _1`$ is the rms of the first derivatives (in statistically isotropic fields both derivatives $`u/x`$ and $`u/y`$ have equal rms). $`A(u)`$ is the fraction of the area in the excursion set and thus equals the cumulative probability function: $`A(u)=_u^{\mathrm{}}p(u^{})𝑑u^{}`$, where $`p(u)`$ is the probability function. $`C(u)`$ and $`G(u)`$ are the length of the contour and the genus of the excursion set per unit area respectively.

In this study we use somewhat different normalization and units. We measure the total contour length in the mesh units and we define $`G`$ as the number of isolated regions minus the number of holes in the excursion set.

### 4.2 Percolating region

The morphology of every isolated region in the excursion set can be characterized by three partial Minkowski functionals : the area, boundary contour length and genus of the region. In order to describe the morphology of the field in more detail, Novikov, Feldman, and Shandarin (1999) used the distribution functions of the partial Minkowski functionals at several level thresholds. Bharadwaj et al. (2000) used the averaged shape parameters in the study of the morphology in the LCRS slices. Here we choose to utilize only the Minkowski functionals of the largest region ($`A_p(\mathrm{\Delta }T),C_p(\mathrm{\Delta }T),G_p(\mathrm{\Delta }T)`$) in addition to the global Minkowski functionals . This is dictated by the relatively small size of the QMASK map: the maximum number of isolated regions is about 30 and therefore the distribution functions are very noisy. The Minkowski functionals of the largest, i.e. , percolating region of Gaussian fields are not known in an analytic form but can be easily measured in Monte Carlo simulations (Shandarin, 2002).

An additional comment may be useful. A naive idea of percolation invokes an image of the region connecting the oposite sides of the map. Thus, the question may arise: what does percolation mean in a map with complex boundaries and holes like the QMASK map (Fig. 2). In particular, which points must be connected at the percolation threshold in a map with irregular boundaries? The fact is that percolation is a phase transition which takes place regardless of the shape and topology of the map (see e.g. Stauffer and Aharony (1992)). It can be studied and relatively easily measured in maps of arbitrary shapes including holes. The Minkowski functionals of the largest region are particularly useful while the direct detection of the connection between the oposite sides is much more unstable even in simple square maps (Dominik and Shandarin, 1992). At percolation transition the area and perimeter of the largest region experience a rapid growth while the genus rapidly decreases (Shandarin, 2002). Similarly to many other statistics the percolation transition is affected mostly by the size of the map and its effective dimensions and much less by its shape or presence of holes. For example, if a three-dimensional map has the shape of a thin wedge then percolation becomes two-dimensional (Shandarin and Yess, 1998). Likewise, a two-dimensional map having the shape of a narrow strip effectively percolates as a one-dimensional map. More accurately, if the smallest size of the map is smaller than the scale of the field then the percolation transition effectively reduces the dimensions. A simple visual inspection of Fig. 2 shows that nothing of the kind of the problems discussed above is present in the QMASK map. The major drawback of the map is its size not the shape.

In addition, we wish to stress that we do not compare the percolation curves measured for the QMASK map with the standard Gaussian curves. Instead we compare them with the curves measured in the reference maps having exactly same irregular shape and holes as the QMASK map. Therefore, the question of accuracy of reproduction of the standard percolation properties in two-dimensional maps becomes largely irrelevant.

### 4.3 Numerical technique

The numerical technique used for measuring the Minkowski functionals is described in detail in (Shandarin, 2002). Here we outline the main features that may differ from numerical methods used by others (Coles, 1988; Gott et al., 1990; Schmalzing and Buchert, 1997; Winitzki and Kosowsky, 1997; Schmalzing and Górski, 1998; Park et al., 2001; Wu et al., 2001).

First, we use neither Koendernik invariants (Koenderink, 1984) nor Croftons’s formulae (Crofton, 1968) which are often used in other studies. A pixelized map is considered to be an approximation to a continuous field. This approach was used by Novikov, Feldman, and Shandarin (1999) but here we use a different and more refined version of the code computing Minkowski functionals . The major effect of this improvement is the numerical efficiency of the code.

The major features of our technique are illustrated by Fig. 3. Suppose we use a regular square lattice and the true contour corresponding to $`u=u_{th}`$ is a simple ellipse shown by the solid line in Fig. 3. The sites satisfying the threshold condition $`u>u_{th}`$ are shown by filled triangles, while the sites with $`u<u_{th}`$ are shown by empty triangles. The union of elementary squares (shown by heavy solid lines) placed on each site satisfying the threshold condition is widely used as an approximation to the region. One can easily count the number of elementary squares, number of edges of the elementary squares and number of vertices in the union. Then by using Crofton’s formulae one computes the area, perimeter and the Euler characteristic of the region. For instance, the area in this approximation is simply the sum of areas of the elementary squares. The perimeter is the number of external edges of the elementary squares. While the area converges to the right value as the lattice constant approaches zero, the perimeter does not (for a detailed discussion of this effect see e.g. Winitzki and Kosowsky (1997)). The perimeter of the set of the elementary squares obviously equals the perimeter of the rectangular region shown by the dotted line. As the lattice parameter approaches zero the perimeter of the region approches the perimeter of the large rectangle shown by the solid line which obviously is a wrong value. Assuming the isotropy of the map one can show that the mean length of a segment is $`4/\pi `$ times larger than its true value (Winitzki and Kosowsky, 1997) and can be corrected only statistically. Our method is completely free of this drawback and does not need such a correction.

In our approach we begin with the construction of the contour points for each isolated region of the excursion set. Suppose that $`u(i,k)<u_{th}`$ but $`u(i+1,k)u_{th}`$ then the contour point lies somewhere between ($`i,k`$) and ($`i+1,k`$) sites. Its coordinates ($`x,y_k`$) can be obtained by solving the linear interpolation equation

$$u_{th}=u(i,k)+\frac{u(i+1,k)u(i,k)}{x(i+1,k)x(i,k)}(xx(i,k))$$
(10)
for the coordinate $`x`$. We construct the contour points by linear interpolation of the field between the sites in both horizontal and vertical directions. These points are marked by the solid circles in Fig. 3. Then, computing the perimeter we assume that the contour points are connected by the segments of a straight line (dashed line). The contours made by this method converges to the right value of the perimeter as the lattice constant approaches zero.

The other major difference consists in allowing a direct linkage of diagonal grid sites in some cases depending on the values in the four neighboring sites. Consder a case when two diagonally adjacent pixels satisfy the threshold condition: $`u(i,k)u_{th}`$ and $`u(i+1,k+1)u_{th}`$ but the sites $`u(i+1,k)<u_{th}`$ and $`u(i,k+1)<u_{th}`$. One usually chooses either consider the $`(i,k)`$ and $`(i+1,k+1)`$ sites always directly linked or always not linked implying that the difference between the two alternatives disappears with vanishing the lattice parameter. Our definition of the direct linkage of two diagonally adjacent pixels is different from both. The diagonally adjacent pixels are directly linked if $`u_{mean}=(u(i,k)+u(i+1,k)+u(i,k+1)+u(i+1,k+1))/4u_{th}`$ and directly not linked otherwise. In the latter case they of course are allowed to be linked by the friend of friend procedure if they have proper friends. This method affects the number of regions, their sizes and the number of holes in the regions and therefore the values of all Minkowski functionals (including genus). The visual inspection of Fig. 4 and Fig. 4 in Park et al. (2001) reveals small but noticeable differences in the genus curves. This treatment of the direct diagonal linkage also considerably reduces the intrinsic anisotropy of a rectangular grid. Since the QMASK map occupies a small part of the sky, we ignore the intrinsic curvature of the map when computing the areas and perimeters.

If a part of the region boundary is the boundary of the map mask, then it does not contribute to the contour length of the region. Similarly, the holes in the map mask contribute neither into the genus nor the contour length.

The analysis of one realization of a Gaussian random field at a hundred levels takes about 0.7s, 2.9s, 12.3s, and 73s for $`64^2`$, $`128^2`$, $`256^2`$, and $`512^2`$ maps respectively on HP C240 workstation. It corresponds approximately to the $`tN_{pix}\mathrm{ln}(N_{pix})`$ dependence of the computational time on the size of the map.

### 4.4 Parameterization by the level of $`\mathrm{\Delta }T`$.

In addition to six Minkowski functionals , we computed the number of isolated regions $`N_c`$ (Novikov, Feldman, & Shandarin 1999). Thus, in total we compute seven morphological parameters

$$M_i(\mathrm{\Delta }T)=[A(\mathrm{\Delta }T),C(\mathrm{\Delta }T),G(\mathrm{\Delta }T),A_p(\mathrm{\Delta }T),C_p(\mathrm{\Delta }T),G_p(\mathrm{\Delta }T),N_c(\mathrm{\Delta }T)],i=1,\mathrm{},7$$
(11)
at a hundred threshold levels $`\mathrm{\Delta }T`$ for the QMASK map and the thousand reference maps. The six panels of Fig. 4 show six Minkowski functionals measured for the QMASK map (solid lines) and the mean values for a thousand Gaussian maps (heavy dashed lines). The thin dashed lines show the 68% and 95% intervals.

The genus curve (left bottom panel) is in good qualitative agreement with that obtained by Park et al. (2001). However, it is worth noting that we use a different normalization and a somewhat different variable on the horizontal axis: we use the level normalized to the rms $`u=\mathrm{\Delta }T/\sigma _{\mathrm{\Delta }T}`$, while Park et al. use a more complex quantity $`\nu _A`$ related to the excursion set of the Gaussian field. Park et al. describe their parameterization as follows: “the area fraction threshold level $`\nu _A`$ is defined to be the temperature threshold level at which the corresponding isotemperature contour encloses a fraction of the survey area equal to that at the temperature level $`\nu _A`$ for a Gaussian field” (the bottom of p.586 ). Some quantitative differences may be also ascribed to a different linkage scheme used by Park et al. (see the previous section). We wish to stress that the parameterization by $`\nu _A`$ has similar properties (including the noise correlation properties) with the parameterization by the area ($`A`$) discussed in the following section. The only difference between them is the geometrical interpretaion.

The number of isolated regions is shown in Fig. 5 for both the temperature and global area parameterizations. The latter is described in detail in Sec. 4.5.

In order to quantify the degree of non-Gaussianity of the QMASK map, we define and compute a functional that quantifies the “distance” between two functions in the functional space of measured statistics. The distance is defined by the integral

$`D_{ik}^{(\mathrm{\Delta }T)}`$ $`=`$ $`\left\{{\displaystyle _{\mathrm{}}^{\mathrm{}}}\left[M_{ik}(\mathrm{\Delta }T)\overline{M_i}(\mathrm{\Delta }T)\right]^2d(\mathrm{\Delta }T)\right\}^{1/2},`$ (12)
$`i`$ $`=`$ $`1,\mathrm{},7;k=1,\mathrm{},1000`$ (13)
where

$$\overline{M_i}(\mathrm{\Delta }T)=\frac{1}{1000}\underset{k=1}{\overset{1000}{}}M_{ik}(\mathrm{\Delta }T)$$
(14)
is the mean value of the statistic as a function of $`\mathrm{\Delta }T`$. The integral in eq. (13) is approximated by the sum over a hundred levels of $`\mathrm{\Delta }T`$. We also compute the distances of the each morphological curve from the Gaussian mean for the QMASK map. The third column of Table 1 shows the percentile of the QMASK map distances from the mean for each morphological statistic parameterized by $`\mathrm{\Delta }T`$. The highest departure from the mean is shown by the $`A(\mathrm{\Delta }T)`$ curve which is the cumulative probability function. It shows that only about 27% of Gaussian curves are closer to the Gaussian mean than the QMASK curve. Other parameters are even closer to the corresponding Gaussian mean.

### 4.5 Parameterization by the total area $`A`$.

The parameterization of the morphological parameters by the level of the field $`u`$ has serious drawbacks mentioned in Introduction. Although the parameterization by $`\nu _A`$ used in most studies of the genus is free of these drwbacks it seems to be more natural for morphological studies to use the total area $`A(\mathrm{\Delta }T)=A_{ES}(\mathrm{\Delta }T)/A_m`$ (where $`A_{ES}(\mathrm{\Delta }T)`$ is the area of the excursion set at $`\mathrm{\Delta }T`$ and $`A_mA_{ES}(\mathrm{})`$ is the total area of the map) as an independent parameter (Yess and Shandarin, 1996). It is one of the Minkowski functionals measured for the map, it represents the cumulative distribution functions of the field, and in addition it has a very simple geometrical interpretaaiton. It also has an additional advantage of decorrelating different statistics as discussed bellow. The $`\nu _A`$ parameterization must have similar properties because it is related to $`A`$ by a monotonic transform (the first eq. of 8).

The transformation from the level parameterization to the area parameterization is a highly nonlinear procedure illustrated in Fig. 6. The thick lines correspond to the QMASK map (solid) and a randomly chosen Gaussian map (dashed). The transformation from $`C(\mathrm{\Delta }T)`$ to $`C(A)`$ involves the $`A(\mathrm{\Delta }T)`$ function shown in the lower panel, which is different for each map. The thin solid lines show the transformation for the QMASK map and dashed lines for the Gaussian map. Note that in order to illustrate this transformation, we reversed the $`A`$ axis in the $`C=C(A)`$ plot.

As an illustration of the diference between the level and area parameterization consider a simple example of a strong nonlinear field. Let the field $`w(x,y)`$ is the exponent of a Gaussian field $`u(x,y)`$: $`w=\mathrm{exp}(\alpha u)`$. The global Minkowski functionals of such a field can be easily derived from eq. 8

$`A(w)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[1\mathrm{erf}\left({\displaystyle \frac{\mathrm{ln}w}{\alpha \sqrt{2}}}\right)\right],`$ (15)
$`C(w)`$ $`=`$ $`{\displaystyle \frac{1}{2R}}\mathrm{exp}\left({\displaystyle \frac{(\mathrm{ln}w)^2}{2}}\right),`$ (16)
$`G(w)`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^{\frac{3}{2}}}}{\displaystyle \frac{1}{R^2}}(\mathrm{ln}w)\mathrm{exp}\left({\displaystyle \frac{(\mathrm{ln}w)^2}{2}}\right).`$ (17)
Fig. 7 illustrates the diference between the Minkowski functionals of the Gaussian field $`u`$ and non-Gaussian field $`w`$ (with $`\alpha =1`$) if both are parameterized by the level. The percolation curves differ strongly as well. On the other hand both the perimeter and genus as the functions of the area remain identical for the both fields because the transformation simply relabels the levels but does not affect the geometry of the contours. This is an example of a trivial non-Gaussian field when all non-Gaussian information is stored in the one-point probability density function. The further discussion of such fields can be found in Shandarin (2002).

Using the area parameterization, we compute seven new morphological functions for each map:

$$M_i(A)=[\mathrm{\Delta }T(A),C(A),G(A),A_p(A),C_p(A),G_p(A),N_c(A)],i=1,\mathrm{},7.$$
(18)
Figure 8 shows the measured quantities as functions of $`A`$. Similarly to Fig. 4 six panels show six Minkowski functionals measured for the QMASK map (solid lines) and the mean values for a thousand Gaussian maps (heavy dashed lines). The thin dashed lines show the 68% and 95% intervals.

The difference of the QMASK curves from the Gaussian mean is quantified similarly to that described before by the distance between curves in the functional space:

$`D_{ik}^{(A)}`$ $`=`$ $`\left\{{\displaystyle _0^1}\left[M_{ik}(A)\overline{M_i}(A)\right]^2𝑑A\right\}^{1/2}`$ (19)
$`i`$ $`=`$ $`1,\mathrm{},7;k=1,\mathrm{},1000`$ (20)
where

$$\overline{M_i}(A)=\frac{1}{1000}\underset{k=1}{\overset{1000}{}}M_{ik}(A).$$
(21)
There is an important difference between the two distances defined by eq. 13 and 20. Using $`dA=p_k(\mathrm{\Delta }T)d(\mathrm{\Delta }T)`$ where $`p_k(\mathrm{\Delta }T)`$ is the estimate of the probability function derived from the $`k`$-th map one can rewrite eq. 20 as an integral over $`\mathrm{\Delta }T`$:

$`D_{ik}^{(A)}`$ $`=`$ $`\left\{{\displaystyle _0^1}\left[M_{ik}(A)\overline{M_i}(A)\right]^2𝑑A\right\}^{1/2}`$ (22)
$`=`$ $`\left\{{\displaystyle _{\mathrm{}}^{\mathrm{}}}\left[M_{ik}(\mathrm{\Delta }T)\overline{M_i}(\mathrm{\Delta }T)\right]^2p_k(\mathrm{\Delta }T)d(\mathrm{\Delta }T)\right\}^{1/2}`$ (23)
Comparing eq. (13) and (23), one notices that the distance between the curves computed in the area space (eq. 20) depends more on the bulk of the probability function and less on the tails of the probability function, since it is weighted by the probability function. This obviously makes the measure more robust because the measurements at extreme levels are subject to large fluctuations. It is worth stressing that the calculation of all morphological parameters $`M_{ik}`$ themselves (eq. 11 and 18) is independent of the type of parameterization because $`M_{ik}(A)M_{ik}(\mathrm{\Delta }T)`$ provided that $`A`$ is measured at corresponding $`\mathrm{\Delta }T`$: $`A=A(\mathrm{\Delta }T)`$. The only quantities that depend on the parameterization are the distances from the Gaussian expectation functions (eqs. 13 and 20) that are used for ranking the maps.

The second column of Table 1 shows the percentile of the QMASK distances from the mean in the area parameterization. Comparing the values corresponding to different parameterizations, one can easily see that the parameterization by $`\mathrm{\Delta }T`$ suggests that the QMASK map is a more typical example of a Gaussian field than that using $`A`$ as an independent parameter. The parameterization by the level suggests that the most exceptional deviation of the QMASK map from Gaussianity is given by the $`A(\mathrm{\Delta }T)`$ curve, but more than 70% of Gaussian maps differ from the mean more than the QMASK map. Thus, the QMASK map looks like a very typical realization of the Gaussian field. Parameterizing by the total area, one may conclude that the Gaussian $`P(A)`$ and $`N_c(A)`$ curves differ more than the QMASK map in only about 21% and 18% of the cases respectively. Thus, this parameterization reveals that the QMASK map is not quite a typical realization of the Gaussian field, but is not very unusual either, quite consistent with Gaussianity.

It is perhaps worth noting that all statistics but the genus of the percolating region, $`G_p`$ show a greater difference of the QMASK map with respect to the corresponding expectation value of the set of Gaussian maps in the $`A`$ parameterization than that in $`\mathrm{\Delta }T`$ parameterization. Since the difference between two parameterizations for $`G_p`$ is not large (12% and 21%) it maybe a statistical fluctuation. One also may note that two columns of Table 1 appear anti-correlated. As we have discussed before the quantities shown in two columns of Table 1 have different sensitivity to the measurements at extreme $`\mathrm{\Delta }T`$ levels both highest and lowest. The numbers in the first column corresponding to the $`A`$ parameterization are more robust.

### 4.6 Cross-correlations

The parameterization of the morphological statistics by the total area makes the set of parameters $`M_i(A)`$ (eq. 18) less correlated than $`M_i(\mathrm{\Delta }T)`$ (eq. 11). Table 2 shows the correlation coefficients of the distances

$$r_{ij}=\frac{\underset{k=1}{\overset{1000}{}}(D_{ik}\overline{D_i})(D_{jk}\overline{D_j})}{\left[\left(\underset{k=1}{\overset{1000}{}}(D_{ik}\overline{D_i})^2\right)\left(\underset{k=1}{\overset{1000}{}}(D_{jk}\overline{D_j})^2\right)\right]^{1/2}}$$
(24)
for all 21 pairs of morphological statistics. The values above the diagonal show the correlations when all statistics were parameterized by the total area of the excursion set, $`A`$, while the values below the diagonal were computed using the parameterization by the level, $`\mathrm{\Delta }T`$.

The correlation coefficient shows for each pair of statistics how well the separation from the mean Gaussian curve in one statistic $`M_i`$ can be predicted from the other one $`M_j`$. The closer $`r_{ij}`$ to unity, the less independent are corresponding statistics. In this respect the $`\mathrm{\Delta }T`$ and $`A`$ parameterizations are very different. For instance, the lowest correlation in the class of the level parameterization is $`0.82`$ between $`G_p(\mathrm{\Delta }T)`$ and $`N_c(\mathrm{\Delta }T)`$ distances and 16 out of 21 cross-correlations are greater than $`0.9`$. In the class of the area parameterization, only one correlation coefficient is greater than $`0.82`$: that between $`G(A)`$ and $`N_c(A)`$ distances ($`0.85`$). Seventeen out of 21 pairs correlate at a level lower than $`0.5`$ thirteen of which at the level lower than $`0.3`$. Although absence of correlation does not prove the statistical independence of two sets of numbers, a correlation coefficient approaching unity indicates a strong statistical dependence.

## 5 Discussion

We present the results of a study of Gaussianity in the QMASK map. In agreement with the first study of topology of the degree-scale CMB anisotropy (Park et al., 2001), which was limited to the genus of QMASK, we conclude that the QMASK map is compatible with the assumption of Gaussian $`\mathrm{\Delta }T/T`$ fluctuations. In this study, we used six morphological statistics in addition to the genus, $`G`$: the areas $`A`$ and $`A_p`$ of the excursion set and the largest (i.e. percolating) region, the contour lengths or perimeters $`C`$ and $`C_p`$ of the excursion set and the largest region; the genus $`G_p`$ of the percolating region and the number of isolated regions $`N_c`$. According to its morphology, the QMASK map is not a very typical example of a Gaussian field, but not very exceptional either: about 20% of a thousand Gaussian maps have greater differences with at least one mean function (Table 1).

We show that the parameterization of the morphological statistic is very important. A naive parameterization by the temperature threshold $`\mathrm{\Delta }T`$ results in a highly correlated (almost degenerate) set of morphological characteristics: 16 out of 21 cross-correlation coefficients are greater than 0.9 and all the rest are greater than 0.82 (Table 2, bellow the diagonal). Representing the same morphological parameters as functions of the total area of the excursion set $`A`$ helps to decorrelate them: now 13 pairs correlate less than 0.3, four more correlate less than 0.5 and the rest less than 0.85 (Table 2, above the diagoanl).

Our measurements of the deviations of the morphological functions from the Gaussian mean strongly rely on the bulk of the probability function (see eq. 23) and therefore are robust to spurious effects.

In the $`A`$ parameterization the strongest deviations of the QMASK map from the Gaussian expectation values were shown by the number of peaks, $`N_c(A)`$ (82% i.e. the QMASK map deviates from the Gaussian expectation number of peaks more than 82% of Gaussian maps), total contour length of the excursion set, $`C(A)`$ (79%) and area of the percolating region, $`A_p(A)`$ (55%) while in the $`\mathrm{\Delta }T`$ parameterization by the total area of the excursion set i.e. the cumulative probability function, $`A(\mathrm{\Delta }T)`$ (27%), genus of the percolating region, $`G_p(\mathrm{\Delta }T)`$ (21%) and contour length of the percolating region, $`C_p(\mathrm{\Delta }T)`$ (20%).

The techniques that we have described here are computationally efficient ($`O(N_{pix}\mathrm{ln}N_{pix})`$), and should be useful also for much larger upcoming datasets such as that from the MAP satellite. Another promising area for future applications is computing Minkowski functionals from foreground maps and topological defect simulations, thereby enabling quantitative limits to be placed on the presence of such structures in CMB maps.

## 6 Acknowledgment

Support for this work was provided by NSF grant AST00-71213, NASA grant NAG5-9194, the University of Pennsylvania Research Foundation, and the University of Kansas General Research Fund. HAF wishes to acknowledge support from the National Science Foundation under grant number AST-0070702 and the National Center for Supercomputing Applications. We wish also to thank the referee for a constructive criticism and useful suggestions.


format: text

text:
# Heating of Gas in Galaxy Groups and Clusters

## 1. Introduction

In a simple universe governed only by gravity, the bremsstrahlung luminosity of galaxy clusters would scale as $`LT^2`$, where $`T`$ is the gas temperature (Kaiser 1986). However, the observed relation is steeper, $`LT^3`$ or $`LT^4`$ (e.g. Arnoud & Evrard 1999; Helsdon & Ponman 2000). Other relations, notably the one between entropy and temperature, also show deviation from self-similarity (Lloyd-Davies, Ponman, & Cannon 2000; Mohr, Mathiesen, & Evrard 1999).

These discrepancies can be explained if gas in clusters experienced some non-gravitational heating before (or when) it flowed into the cluster halo (Kaiser 1991; Evrard & Henry 1991). However, there is no general consensus about the details of the heating process.

Here we compare two pictures for the heating that have been proposed by a number of authors in recent years: the external heating scenario, in which the gas is heated at high redshifts, before it enters massive halos, and the internal heating scenario, in which the energy is injected when some or most of the gas is already inside a massive halo. In this latter case we assume that the sources of energy are SNII explosions. The motivation for considering SNII heating a viable form of heating is the success of our models in reproducing the evolution of ISM in NGC 4472 (Brighenti & Mathews 1999a). NGC 4472 is a giant elliptical at the center of a subcluster in Virgo, consistent with being the remnant of a galaxy group, probably stripped of the outer regions when it entered the Virgo cluster. We showed that SNII heating was able to reproduce the gas density and temperature profiles (and therefore also the entropy profile).

## 2. The Simulations

We use a modified version of the hydrocode ZEUS (Stone & Norman 1992) and we assume spherical symmetry. The code follows two fluids: a normal, collisional gas, and a collisionless fluid which represents the dark matter. Our groups and clusters evolve from a single top-hat perturbation in a $`\mathrm{\Lambda }`$CDM universe ($`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Lambda }=0.7`$, $`h=H_0/(100\mathrm{km}/\mathrm{s}\mathrm{Mpc})=0.725`$, $`\mathrm{\Omega }_\mathrm{b}=0.037`$). The dark matter accumulates in a Navarro, Frenk, & White (1996) halo by design. We consider a set of three objects with different virial masses (at the present time):

\- $`M_{\mathrm{vir}}=4.7\times 10^{13}`$ M$`_{}`$: “the group”;

\- $`M_{\mathrm{vir}}=2.2\times 10^{14}`$ M$`_{}`$: “the poor cluster”;

\- $`M_{\mathrm{vir}}=1.2\times 10^{15}`$ M$`_{}`$: “the rich cluster”;

We focus here mainly on groups, since lower mass system are more sensitive to heating, and our models, by ignoring the complex merging events, are anyway less appropriate to describe the formation and evolution of large clusters.

External pre-heating. We assume an extreme form of preheating: at very high redshift, $`z_\mathrm{h}=9`$, we reset everywhere the gas density to the mean baryon density, $`\rho =\overline{\rho }_\mathrm{b}(z=9)`$. At that epoch, the temperature is raised to some constant level $`T_\mathrm{h}`$. We consider 4 levels of heating: $`T_\mathrm{h}=10^4,5\times 10^6,10^7,3\times 10^7`$ K, corresponding to $`1.3\times 10^4,0.65,1.3,3.9`$ keV/particle. These amount of heating is characterized by the numbers $`1`$ to $`4`$ respectively, in our nomenclature. The entropy parameter $`S=T/n_\mathrm{e}^{2/3}`$ corresponding to these levels of heating, which depends on the heating epoch $`z_\mathrm{h}`$ through $`n_\mathrm{e}`$, is $`S=0.025,125,250,750`$ keV cm2.

Gas is allowed to cool and to dropout of the flow. This last process is modeled in the usual way (e.g. Sarazin & Ashe 1989) adding a sink term in the continuity equation $`\dot{\rho }_{do}=q\rho /t_{\mathrm{cool}}`$ with $`q=1`$.

Internal heating. In this series of models, heating is assumed to be the result of star formation occurring inside the group or cluster. Thus, we need to assume a schematic scenario for star formation in these systems. At $`z_{}=3`$ (2 Gyr after the big-bang) we form stars from cooled gas (conserving baryons) and release SNII energy inside the accretion shock radius $`r_{\mathrm{sh}}(z_{})`$. All the gas inside $`r_{\mathrm{sh}}`$ is assumed to get the same amount of energy per unit mass. Te total amount of energy released is $`E_{\mathrm{SN}}=\nu \eta _{\mathrm{Salpeter}}E_0(M_{}/M_{})`$, where $`\eta _{\mathrm{Salpeter}}0.007`$ is the number of SNII per unit solar mass predicted by a stellar population with a Salpeter IMF, $`E_0=10^{51}`$ erg is the kinetic energy released by a single SNII and $`M_{}`$ is the total mass of stars. The parameter $`\nu `$ controls the amount of heating (i.e. the number of SNII). We consider 4 models with $`\nu =0.5,1,2,4`$, labeled with numbers 1 to 4, in analogy with the external heating models.

In more familiar units, SNII inject $`2.4\nu `$ keV/particle to the gas inside $`r_{\mathrm{sh}}`$ at $`z_\mathrm{h}`$. This is consistent with the global value $`0.22`$ keV/particle which is derived assuming a global star formation efficiency $`M_{}/M_{\mathrm{baryon}}=0.1`$ (Fukugita, Hogan, & Peebles 1998) and assuming that the SNII energy is shared among all the baryons. In our models, instead, SNII heat only the central part of clusters, similar to the models proposed by Loewenstein (2000). In these models the hydrodynamic equations are modified to take into account the mass and energy injected by stars and SNIa of the central, dominant galaxy. Full details about the simulations can be found in Brighenti & Mathews (2001).

## 3. Results

### 3.1. The $`LT`$ relation

After evolving to the current time, $`t=13`$ Gyrs, the location of our models can be compared with observed clusters in the $`LT`$ plot as shown in Fig. 1. In general, heating has a little effect on the emission weighted temperature (heated clusters are not hotter!). Paradoxically, models with maximum heating are always the coolest. Instead, the luminosity generally decreases as heating increases, a result of the lower mean gas density. Groups, having lower virial temperatures, are more affected by the heating, while rich clusters are quite insensitive to it.

The left panel of Fig. 1 shows the results for the external heating scenario. Models “1” without heating nicely follow the self-similar prediction $`LT^2`$. The group model with maximum heating (model “4”) lies among the observed groups, and the series of models “4” follows a relation $`LT^3`$, similar to the observed one. However, it requires that $`3.9`$ keV/particle are dumped in the gas, or that a entropy floor $`S750`$ keV cm2 is established at $`z_\mathrm{h}=9`$. The energy budget needed is a function of $`z_\mathrm{h}`$. If $`z_\mathrm{h}=5`$, $`1`$ keV/particle is sufficient to decrease the luminosity of groups to the observed level.

The internal heating models are less sensitive to the amount of energy injected (right panel of Fig. 1) and they tend to lie near the upper envelope of the observed $`LT`$ data, where strong cooling flows are dominated by a massive central galaxy as we have assumed. Groups models with $`\nu 1`$ have X-ray properties consistent with observed groups, provided the efficiency of SNII heating is high (we assumed efficiency $`=1`$). However, it is likely that a significant fraction of the SNII energy is lost by radiation, and a more realistic constraint may be $`\nu 23`$.

### 3.2. The $`ST`$ relation

The behavior of entropy for external heating models is illustrated in Fig. 2. In the left panel we compare the entropy evaluated at $`r=0.1r_{\mathrm{vir}}`$ with the data of Lloyd-Davies et al. (2000). The series with maximum heating (models “4”, indicated with the largest circles), which fits best in the $`LT`$ plot, has a much larger entropy than observed groups and poor clusters. Models “3” (which requires $`1.3`$ keV/particle) may be the best compromise.

The difficulty for external heating models to simultaneously fit the $`LT`$ and $`ST`$ relations is also illustrated by the radial entropy profiles (right panel of Fig. 2). The computed group entropy profiles are compared with the observed profile of NGC 2563 group (dot-dashed line; Trinchieri, Fabbiano, & Kim 1997); see also the profiles in Lloyd-Davies et al. (2000). Real groups have entropy profiles that increase monotonically with radius, while for model “4” $`S`$ is uniform and model “3” shows a flat entropy core for $`r<300`$ kpc. This behavior is due to the inefficiency of radiative cooling in the low density cores of these strongly heated groups.

Internal heating group models fall nicely among the observations in the $`ST`$ plane (Fig. 3), regardless of the heating parameter $`\nu `$. Strong radiative cooling regulates the entropy in the central regions to have similar values for all adopted $`\nu `$. Radial entropy profiles are also consistent with observations (Fig. 3, bottom panel), although our models are somewhat denser overall with slightly lower $`S(r)`$.

### 3.3. Baryon fractions and cooling times

Models experiencing external preheating at very early times differ from those heated internally by SNII in several other respects. Notably, it appears that internal heating removes baryons more efficiently than preheating at $`z_\mathrm{h}=9`$. We find for preheated groups that the baryon fraction at $`r_{\mathrm{vir}}900`$ kpc is $`f_\mathrm{b}0.11`$ for model “3”, almost equal to the cosmic baryon fraction assumed, $`0.123`$. Only for the maximally preheated model “4”, is the baryon fraction significantly lower than the cosmic one: $`f_\mathrm{b}0.07`$.

Groups heated internally by SNII have low baryon fraction even when a Salpeter IMF is assumed ($`\nu =1`$): $`f_\mathrm{b}0.075`$ ($`f_\mathrm{b}0.055`$ for $`\nu =2`$). The fraction of mass in gas at $`r_{\mathrm{vir}}`$ is $`0.04`$ and $`0.03`$ for $`\nu =1`$ and $`\nu =2`$, respectively.

A further distinction between external and internal heating models is the central cooling time which, if lower than the age of the system, may indicate the presence of a cooling flow. We find that models with preheating strong enough to fit observations in the $`LT`$ plot (models 3 and 4 in Fig. 1) never develop cooling flows (or, more precisely, never have $`t_{\mathrm{cool}}<\mathrm{age}`$), contrary to many observed groups and clusters. All internal heating models, instead, develop strong cooling flows. In particular, group models have $`\dot{M}110`$ M$`_{}`$ yr-1.

## 4. Heating-enrichment connection

Both the external and internal heating scenarios seem to require more energy than that provided, via SNII, by the observed stellar content with a Salpeter IMF. A reasonable requirement may be $`1`$ keV/particle (a value consistent with most models proposed in recent years). SNII also produce metals, so we should ask: is the number of SNII needed to heat the gas consistent with the number of SNII necessary to produce the metal content of the universe?

A Salpeter IMF (from 0.1 to 100 M$`_{}`$) produces $`0.007`$ SNII per M$`_{}`$ of stars formed. Assuming a global star formation efficiency $`\mathrm{\Omega }_{}/\mathrm{\Omega }_{\mathrm{baryon}}=0.1`$ (Fukugita et al. 1998), we get $`ϵ0.22\nu `$ keV/particle. Thus, to generate $`1`$ keV/particle we need $`\nu 5`$.

The present day universe has a global iron abundance $`<Z_{\mathrm{Fe}}>0.30.4`$ solar meteoritic units (Renzini 1997). Assuming again $`\mathrm{\Omega }_{}/\mathrm{\Omega }_{\mathrm{baryon}}=0.1`$ with an average Fe yield per SNII $`<y_{\mathrm{Fe}}>0.1`$ M$`_{}`$ (e.g. Gibson, Loewenstein, & Mushotzky 1997), the averaged metallicity produced by all SNII is $`<Z_{\mathrm{Fe}}>_{\mathrm{SNII}}0.053\nu `$ solar. To make the observed metallicity, a high SNII production efficiency is needed: $`\nu 56`$ (this value would be reduced somewhat if SNIa contribute a significant fraction of iron). Thus, it appear that both heating and cosmic metallicity may be produced by a stellar population with $`\nu 5`$. The agreement between these two estimates of $`\nu `$ supports SNII as the source of non gravitational heating. However, it should be noted that such a large production of SNII may be inconsistent with the chemical evolution of ISM in elliptical galaxies (Brighenti & Mathews 1999b).

## 5. Conclusions

External heating models fit the data in the $`LT`$ plot provided the preheating is sufficiently strong: $`14`$ keV/particle, depending on $`z_\mathrm{h}`$. However, successful models in the $`LT`$ plane have entropies that exceed observed values (a cautionary note: we are well aware that our models are approximate, and it’s possible that less extreme preheating scenarios may overcome the problems pointed out by the present work). The competing models with internal heating by SNII fit the whole set of X-ray observations better and more plausibly, but they likely require a production of energy per unit of stellar mass larger than that based on a Salpeter IMF. This may not be a severe demand since the cosmic metallicity itself requires such a higher number of SNII per unit of stellar mass.

## References

Allen, S.W., & Fabian, A.C. 1998, MNRAS, 297, 57

Arnoud, M., & Evrard, A.E. 1999, MNRAS, 305, 631

Brighenti, F., & Mathews, W.G. 1999a, ApJ, 512, 65

Brighenti, F., & Mathews, W.G. 1999b, ApJ, 515, 542

Brighenti, F., & Mathews, W.G. 2001, ApJ, 553, 103

Evrard, A.E., & Henry, J.P. 1991, ApJ, 383, 95

Fukugita, M., Hogan, C.J., & Peebles, P.J.E. 1998, ApJ, 503, 518

Gibson, B.K., Loewenstein, M., & Mushotzky, R.F. 1997, MNRAS, 290, 623

Helsdon, S.F., & Ponman, T.J. 2000, MNRAS, 315, 356

Kaiser, N. 1986, MNRAS, 222, 323

Kaiser, N. 1991, ApJ, 383,104

Lloyd-Davies, E.J., Ponman, T.J. 2000, & Cannon, D.B. MNRAS, 315, 356

Loewenstein, M. 2000, ApJ, 532, 17

Navarro, J.F, Frenk, C.S., & White, S.D. 1996, ApJ, 462, 563

Renzini, A. 1997, ApJ, 488, 35

Sarazin, C.L., & Ashe, G.A. 1989, ApJ, 345, 22

Stone, J.M., & Norman, M.L. 1992, ApJS, 80, 573

Trinchieri, G., Fabbiano, G., &Kim, D.-W. 1997, A&A, 318, 361


text:
# Giant Spin Relaxation Anisotropy in Zinc-Blende Heterostructures

## I Introduction

Spin relaxation processes have significant effect in optical and kinetic properties of semiconductors. They play important role in optical orientation of electrons and nuclei  and in anomalous magnetoresistance caused by weak localisation . Both theoretical calculations and experimental data analysis have been carried out assuming that one spin relaxation mechanism dominates only. Therewith in spite of the strong anisotropy of spin-orbit scattering, the relaxation times of spin lying in the plane of a heterostructure with zinc-blende lattice turn out to be independent on orientation with respect to crystallographic axes.

Real heterostructures differ from investigated ideal objects in that several spin relaxation mechanisms exist . A spin relaxation mechanism due to only cubic in wave vector terms of the bulk Hamiltonian was investigated for rectangular quantum wells (QWs) in Ref. . It was noted that even in an asymmetrical GaAs QW, the efficiency of another mechanism due to linear in two dimensional wave vector terms is negligibly small. The authors of Ref.  analysing experimental data on anomalous magnetoresistance in InGaAs QWs demonstrated that the both mechanisms may be comparable in efficiency. But in Refs. , it was mentioned that the both mechanisms are additive in spin relaxation.

This communication is devoted to an investigation of spin relaxation processes when several mechanisms of spin-orbit scattering exist. We show that contributions of these mechanisms interfere and their simultaneous action leads to the strong anisotropy of spin relaxation even in the plane of a QW.

## II Theory

In zinc-blende semiconductors, spin relaxation of electrons is well known to be due to spin-orbit splitting of conduction band. In a bulk crystal, the splitting is cubic in wave vector. In a QW structure, the corresponding Hamiltonian has to be averaged over the motion along the growth axis. We consider the QW grown along $`z`$-direction parallel to  and choose $`x`$ and $`y`$ directions coinciding with crystallographic axes. At relatively small carrier concentrations, one can neglect cubic in 2D wave vector terms and the Hamiltonian has the form:

$$H_1=a_1(\sigma _xk_x\sigma _yk_y).$$
(1)
Here $`\sigma _i`$ ($`i=x,y`$) is the Pauli matrix, $`k_i`$ is the wave vector component in the plane of the QW and $`a_1`$ is a constant.

In asymmetrical heterostructures, there is a contribution to the Hamiltonian which is absent in the bulk :

$$H_2=a_2(\sigma _xk_y\sigma _yk_x),$$
(2)
where $`a_2`$ is the constant determined by heterointerface properties.

To calculate spin relaxation times, one has to consider the dynamics of spin density matrix, $`\rho `$:

$$\frac{\rho }{t}=\frac{i}{\mathrm{}}[H,\rho ],$$
(3)
where the total two dimensional Hamiltonian is:

$$H=\frac{\mathrm{}^2k^2}{2m}+V+H^{}.$$
(4)
Here $`m`$ is an effective electron mass, $`V(x,y)`$ is a scattering potential and

$`H^{}=H_1+H_2.`$

We assume that the scattering is elastic and independent on spin indices.

Since $`H^{}`$ is a small perturbation, the spin relaxation times turn out to be much longer than isotropisation times of momentum distribution of electrons. For this reason it is convenient to represent the density matrix as a sum :

$`\rho =\overline{\rho }+\rho ^{},\overline{\rho ^{}}=0,`$

where the bar means averaging over the directions of $`𝒌`$. Here $`\overline{\rho }`$ depends on $`\epsilon =\mathrm{}^2k^2/2m`$ and the anisotropic part of the density matrix is due to $`H^{}`$ only. Hence $`\rho ^{}`$ is proportional to $`H^{}`$, i.e. $`\rho ^{}`$ is small in comparison to $`\overline{\rho }`$.

Then in the first order in $`H^{}`$, the Eq. (3) has the form:

$$\frac{\rho ^{}}{t}=\frac{i}{\mathrm{}}[H^{},\overline{\rho }]\underset{𝐤^{}}{}W_{\mathrm{𝐤𝐤}^{}}\left[\rho ^{}(𝒌)\rho ^{}(𝒌^{})\right],$$
(5)
$$\frac{\overline{\rho }}{t}=\frac{i}{\mathrm{}}\overline{[H^{},\rho ^{}]}.$$
(6)
Here $`W_{\mathrm{𝐤𝐤}^{}}`$ is the probability for scattering from the potential $`V`$ from the state with $`𝒌`$ to the state with $`𝒌^{}`$.

One can see from the Eq. (5) that $`\rho ^{}`$ relaxes in the time which is of order of momentum relaxation time, but $`\overline{\rho }`$ relaxes in the longer time which is determined by $`H^{}`$ — see (6). At these long times, the Eq. (5) reduces to a quasistationary one:

$$\underset{𝐤^{}}{}W_{\mathrm{𝐤𝐤}^{}}\left[\rho ^{}(𝒌)\rho ^{}(𝒌^{})\right]=\frac{i}{\mathrm{}}[H^{},\overline{\rho }].$$
(7)
Finding $`\rho ^{}(𝒌)`$ from here and substituting it to the Eq. (6), one can obtain the closed equation for $`\overline{\rho }`$:

$$\frac{\overline{\rho }}{t}=\frac{1}{\mathrm{}^2}\underset{n}{}\tau _n[H_n^{},[H_n^{},\overline{\rho }]].$$
(8)
Here

$$H_n^{}=\frac{d\phi _𝐤}{2\pi }H^{}(𝒌)\mathrm{exp}(in\phi _𝐤),$$
(9)
where $`\phi _𝐤`$ is the angle between $`𝒌`$ and $`x`$ axis, and

$$\frac{1}{\tau _n}=𝑑\theta W_{\mathrm{𝐤𝐤}^{}}(1\mathrm{cos}n\theta ),$$
(10)
where $`\theta =\phi _𝐤\phi _𝐤^{}`$. The Eq. (8) clearly demonstrates that it is spin-orbit interaction which causes $`\overline{\rho }`$ relaxation.

After substituting $`\overline{\rho }`$ in a form

$`\overline{\rho }=f_0+{\displaystyle \frac{1}{2}}𝝈\mathit{\ae },`$

the Eq. (8) reduces to following equations:

$$\dot{f_0}(\epsilon ,t)=0,$$
(11)
$$\dot{\ae }_i(\epsilon ,t)=\mathrm{\Gamma }_{ij}(\epsilon )\ae _j(\epsilon ,t),$$
(12)
where

$$\mathrm{\Gamma }_{ij}=\frac{1}{\mathrm{}^2}\mathrm{Tr}\left\{\underset{n}{}\tau _n[H_n^{},[H_n^{},\sigma _j]]\sigma _i\right\}.$$
(13)
An initial condition may be derived considering the spin dynamics after the time $`\tau _n`$, but before the spin relaxation time. In the time $`\tau _n`$, the density matrix becomes isotropic but the spin relaxation processes do not start yet. Therefore :

$$f_0(\epsilon )=\frac{1}{2}[F_+(\epsilon )+F_{}(\epsilon )],$$
(14)
$$\mathit{\ae }(\epsilon )=𝒔[F_+(\epsilon )F_{}(\epsilon )],$$
(15)
where $`𝒔`$ is the unit vector along the spin and $`F_\pm (\epsilon )`$ are distribution functions of particles with the spin projection on $`𝒔`$ equal to $`\pm 1/2`$.

Taking into account that the spin density, $`𝑺(t)`$, is the average of $`\mathit{\ae }`$ over $`\epsilon `$, one can obtain the kinetic equation for it at the time longer than $`\tau _n`$:

$$\dot{S}_i=\frac{S_j}{\tau _{ij}},$$
(16)
where the tensor of reciprocal spin relaxation times is:

$$\frac{1}{\tau _{ij}}=\frac{𝑑\epsilon [F_+(\epsilon )F_{}(\epsilon )]\mathrm{\Gamma }_{ij}(\epsilon )}{𝑑\epsilon [F_+(\epsilon )F_{}(\epsilon )]}.$$
(17)
The Eq. (17) represents the extension of the results of Ref.  to the case of an arbitrary spin-orbit interaction $`H^{}`$ and takes into account the anisotropy of scattering.

Substituting $`H^{}=H_1+H_2`$ into (13) and then $`\mathrm{\Gamma }_{ij}`$ into (17) we have:

$$\frac{1}{\tau _{zz}}=C(a_1^2+a_2^2),\frac{1}{\tau _{zx}}=\frac{1}{\tau _{zy}}=0,$$
(18)
$$\frac{1}{\tau _{xx}}=\frac{1}{\tau _{yy}}=\frac{C}{2}(a_1^2+a_2^2),\frac{1}{\tau _{xy}}=Ca_1a_2,$$
(19)
where

$$C=\frac{8m}{\mathrm{}^4}\frac{𝑑\epsilon [F_+(\epsilon )F_{}(\epsilon )]\tau _1(\epsilon )\epsilon }{𝑑\epsilon [F_+(\epsilon )F_{}(\epsilon )]}.$$
(20)
Equations (1820) generalize the results  for the case of two spin relaxation mechanisms.

## III Discussion

It follows from the Eq. (19) that because of the two mechanisms ($`a_1a_20`$), the spin relaxation in the plane of the QW becomes anisotropic. It should be emphasized that if there is only one mechanism ($`a_1a_2=0`$), then the spin relaxation is isotropic in spite of the cubic symmetry of the Hamiltonian $`H_1`$ or $`H_2`$. Thus the cubic anisotropy of the real QW structure manifests itself due to the interference of two spin relaxation mechanisms.

The system (16) may be rewritten as follows:

$$\dot{S}_x\pm \dot{S}_y=\frac{S_x\pm S_y}{\tau _\pm },$$
(21)
where

$$\frac{1}{\tau _\pm }=\frac{C}{2}(a_1\pm a_2)^2.$$
(22)
The times $`\tau _+`$ and $`\tau _{}`$ describe the relaxation of the spin oriented along the directions  and \[1$`\overline{1}`$0\] respectively.

The most bright manifestation of spin relaxation anisotropy occurs if $`a_1=\pm a_2`$. In this case, one of the times $`\tau _+`$ or $`\tau _{}`$ becomes infinite. Therewith the other is equal to $`\tau _{zz}`$.

Efficiency of the mechanisms depends on both the material and the shape of the QW. It was shown that in a rectangular GaAs/AlGaAs QW, the mechanism (1) dominates , and in an asymmetrical InGaAs/AlAs QWs the mechanism (2) is the most important  or they are comparable .

It is clear that if both $`a_1`$ and $`a_2`$ are not equal to zero, then the spin sublevels split. In the work  it is shown that $`a_1`$ and $`a_2`$ may be comparable in magnitude and, hence, the spin splitting is strongly anisotropic.

The spin relaxation anisotropy results from the initial $`T_d`$ symmetry of the zinc-blende semiconductor. For this reason, the similar effect can take place in a strained bulk crystal. The corresponding Hamiltonian linear in 3D wave vector, $`𝒌`$, and components of an elastic strain tensor, $`u_{ij}`$, has the form:

$`H^{}(u)=A_1u_{ii}(\sigma _{i+1}k_{i+1}\sigma _{i+2}k_{i+2})`$ (23)
$`+A_2u_{ij}(\sigma _ik_j\sigma _jk_i).`$ (24)
Here $`i,j=x,y,z`$, $`i+3i`$, $`A_1`$ and $`A_2`$ are constants. Doing calculations for $`1/\tau _{ij}`$ in a way similar to (17), one can obtain three different spin relaxation times. It can be shown that the maximum anisotropy may be achieved if

$$A_1u_{xx}=A_1u_{yy}=A_1u_{zz}/2=A_2u_{xy}/3$$
(25)
with the rest of $`u_{ij}=0`$. Therewith two spin relaxation times are equal to each other and the third is infinite. Note that the tensor $`u_{ij}`$ determined by (25) may be obtained by applying two uni-axial strains along the axes  and  and they are not restricted to uni-axial strain along any axes.

## IV Conclusion

The possibility for spin relaxation suppression was noted in Ref.  for a QW grown along  direction when the spin is oriented along the same axis. The present work shows that the spin relaxation rate also decreases for  direction, but in a QW grown in the symmetrical direction . Therefore this decrease takes place for the spin lying in the plane of the QW.

Analysing weak localisation effect, the authors of Refs.  showed that the mechanisms (1) and (2) suppress each other in anomalous magnetoresistance, but they are additive in spin relaxation. The present analysis shows that the suppression occurs in the spin relaxation also. Besides, we have found that spin relaxation is anisotropic even in the plane of the QW.

We thank U. Rössler, M.I. Dyakonov and A.N. Titkov for fruitful discussions.

This work was financially supported by Russian Foundation for Basic Research (grant 98-02-18424) and programm “Physics of Solid State Nanostructures” (grant 97-1035).


text:
# References

gr-qc/0107086

July 2001

Space-time quantum solves three experimental paradoxes

Giovanni AMELINO-CAMELIA

Dipartimento di Fisica, Universitá “La Sapienza”, P.le Moro 2, I-00185 Roma, Italy

ABSTRACT

I show that a Planck-scale deformation of the relativistic dispersion relation, which has been independently considered in the quantum-gravity literature, can explain the surprising results of three classes of experiments: (1) observations of cosmic rays above the expected GZK limit, (2) observations of multi-TeV photons from the BL Lac object Markarian 501, (3) studies of the longitudinal development of the air showers produced by ultra-high-energy hadronic particles. Experiments now in preparation, such as the ones planned for the GLAST space telescope, will provide an independent test of this solution of the three experimental paradoxes.

Theoretical physics has been puzzling over the structure of space-time at distance scales of the order of the Planck length111$`L_p\sqrt{\mathrm{}G/c^3}1.610^{33}cm`$, where $`\mathrm{}`$ is the reduced Planck constant, $`G`$ is the gravitational constant and $`c`$ is the speed-of-light constant $`L_p`$ for several decades . Unfortunately, there was no experimental counter-part for these sizeable theoretical effort. All effects predicted by Planck-scale theories are very small, since they are strongly suppressed  by the smallness of the Planck length, and this has kept the Planck-length structure of space-time beyond the reach of available experimental sensitivities. Only over the last 15 years some ideas for experimental investigations of this realm have started to emerge (see, e.g., Refs. \[4-16\]), relying on the remarkable sensitivities of advanced experiments now in preparation.

At a time when it appeared to be rather exciting  that some experiments could finally start exploring, in a few years, the structure of space-time at the Planck scale, it was recently argued (see, e.g., Refs. ) that we might be already witnessing the first manifestations of Planck-length physics, since quantum-gravity models can provide solutions for some experimental paradoxes that presently confront the astrophysics community: observed violations of the cosmic-ray GZK limit  and observed violations of the analogous 10-TeV limit  that applies to photons from Markarian 501 (a BL Lac object at a redshift of 0.034, i.e. $`157`$ Mpc).

I shall revisit these analyses of observations of cosmic rays and Markarian-501 photons and I shall then consider another independent experimental paradox which emerged from a very recent analysis  of data on the longitudinal development of the air showers produced by ultra-high-energy hadrons. Remarkably, I find that all three paradoxes are solved by the same phenomenological model of Planck-length physics, characterized by a deformed dispersion relation (without free parameters!). The preliminary evidence emerging from these three paradoxes could amount to establishing the first Planck-scale property of space-time. Because of the profound significance of such a discovery it is at present necessary to proceed very cautiously, and, accordingly, I shall also emphasize the aspects of these three paradoxes that still require further investigation. The situation will be fully clarified, as I discuss in the final part of this note, by forthcoming observations by the GLAST space telescope  which can provide an independent and robust test of the scenario that is being encouraged by the three experiments here analyzed.

Let me start briefly reviewing the three experimental paradoxes. They all involve the kinematic rules for particle production in a continuum classical space-time, but the relevant particle-production processes are different and the energy scales involved are also different:
Cosmic-ray paradox. Cosmic rays can interact with the Cosmic Microwave Background Radiation (CMBR), producing pions. Taking into account the typical energy of CMBR photons, and assuming the validity of the kinematic rules for the production of particles in our present, classical and continuous, description of space-time (conventional relativistic kinematics), one finds that these interactions should lead to an upper limit $`E<510^{19}`$eV, the GZK limit , on the energy of observed cosmic rays. Essentially, cosmic rays emitted with energies in excess of the GZK limit should loose energy on the way to Earth by producing pions, and, as a result, should still satisfy the GZK limit when detected in our observatories. Instead, several cosmic-rays above the GZK limit (with energies as high as $`310^{20}`$eV) have been observed .
Markarian-501 paradox. Just in the same way in which one obtains the GZK limit for cosmic rays, one also obtains a limit on the maximum energy of photons that can reach us from distant sources. The relevant process for establishing this limit is pair-production absorption of high-energy photons due to interactions with the Far Infrared Background Radiation (FIBR). For the high-energy photons emitted by Markarian 501, FIBR absorption should  become efficient around $`10`$TeV. Markarian-501 photons with higher energies should collide with FIBR photons, disappearing into an electron-positron pair, and should therefore not be able to reach our observatories. Instead, Markarian-501 photons with energies as high as $`24`$TeV have been observed .
Pion-stability paradox. The kinematic rules for the production of particles also govern the structure of the air showers produced by high-energy particles. In particular, they allow to predict some features of the longitudinal development of the showers, such as the probability distribution of the maximum depth of the showers. Experimental data on the longitudinal development of the air showers produced by ultra-high-energy hadronic primaries appear to be in disagreement  with these predictions. The analysis reported in Ref.  suggests that the observed longitudinal development of the air showers could be explained by assuming that ultra-high-energy neutral pions are much more stable than low-energy ones, as if, at ultra-high energies, the available phase space for decay in two photons was becoming smaller (perhaps, at some energy, even vanishingly small ) than the one predicted by conventional relativistic kinematics.

As announced, I intend to show that these paradoxes can be solved by adopting a deformation of the standard relativistic dispersion relation $`E^2=p^2+m^2`$. In the quantum-gravity literature there has been discussion of various mechanisms for the emergence of deformed dispersion relations. The most radically new of these scenarios is the one  in which a deformed dispersion relation is assumed to emerge as a reflection of the deformed symmetries of a quantum version222Of course, quantum versions of Minkowski space-time usually do not enjoy classical symmetries. In fact one of the schemes considered in Ref.  turned out to be connected with the $`\kappa `$-Minkowski noncommutative space-time . of (quasi-)flat space-time. Alternatively, deformed dispersion relations can also emerge as a property333Space-time foam could play the role of a dispersion-inducing environment . of the space-time foam background , as illustrated by the phenomenological model considered in Ref.  and by the analysis of loop-quantum-gravity “weave states”  reported in Ref. . In attempts to unify space-time physics with quantum mechanics one can also encounter deformed dispersion relations as a result of the presence of more ordinary (no foam) backgrounds; for example, string theory in certain magnetic-field-like backgrounds admits an effective-theory description in terms of a field theory in noncommutative geometry with associated emergence of deformed dispersion relations (see, e.g., Ref. .)

From the analysis reported here below, the careful reader will easily realize that the key ingredient for the solution of the mentioned three paradoxes is a deformation of the dispersion relation. It appears plausible that more than one of the quantum-gravity schemes that motivate the analysis of such deformations would provide solutions of the paradoxes. However, I shall here focus on the space-time-foam phenomenological scheme advocated in Ref. , based on the dispersion relation444The precise all-order function $`f`$ is not discussed in Ref. , but present and forthcoming experiments are anyway only sensitive to the leading-order correction to the presently-adopted dispersion relation $`E^2=p^2+m^2`$.

$$E^2=f(E,p;m;L_p)p^2+m^2L_pEp^2.$$
(1)
In fact, my analysis will be facilitated by the simplicity of this dispersion relation (for example, in other schemes  one should worry about a polarization dependence in the analysis of processes involving photons). Moreover, the fact that (1), unlike other proposed deformations of the dispersion relation, has no free parameters (but see closing remarks on alternatives with different sign and magnitude of the deformation term) renders particularly significant the fact that three independent experimental paradoxes find a common solution in this scheme.

Let me start with the cosmic-ray paradox. As first observed by Kifune , the deformed dispersion relation (1) would affect the kinematics of particle-production processes, such as photopion production ($`p+\gamma p+\pi `$, with the usual notation $`p`$, $`\gamma `$, $`\pi `$ to denote protons, photons and pions respectively) which, as mentioned, is relevant for the cosmic-ray paradox. Combining (1) with equations for the conservation of energy and momentum one finds that a collision between a proton of energy $`E`$ and a CMBR photon of (much smaller) energy $`ϵ`$ can produce a pion (and a proton) only if

$$E>\frac{2m_pm_\pi +m_\pi ^2}{4ϵ}+L_p\frac{(2m_p+m_\pi )^3m_\pi ^3}{256ϵ^4}\left(1\frac{m_p^2+m_\pi ^2}{(m_p+m_\pi )^2}\right),$$
(2)
where $`m_p`$ ($`m_\pi `$) is the proton (pion) mass. The $`L_p0`$ limit of this condition of course describes the conventional photopion-production threshold. In spite of the smallness of $`L_p`$ the correction term turns out to be significant for the cosmic-ray paradox (the magnitude of the correction term is suppressed by $`L_p`$ but is boosted by the large ratios $`m_p/ϵ`$, $`m_\pi /ϵ`$). In fact, one finds  that, according to (2), even at $`E310^{20}`$eV photopion production on CMBR photons is still not possible, providing an explanation for the fact that cosmic rays of such high energies are being observed.

As shown in Refs. , the Markarian-501 paradox can be explained in a completely analogous manner. Combining (1) with the relevant equations for the conservation of energy and momentum one finds that the process $`\gamma +\gamma e^{}+e^+`$ is only possible if

$$>\frac{m_e^2}{ϵ}+L_p\frac{m_e^6}{8ϵ^4},$$
(3)
where $`m_e`$ is the electron mass, $``$ is the energy of the (hard) photon emitted by Markarian 501, and $`ϵ`$ denotes again the energy of the (soft) background photon (here assumed to be a FIBR photon). The $`L_p0`$ limit of (3) of course describes the conventional pair-production threshold. The $`L_p`$-dependent term again represents a significant correction ($`m_e/ϵ1`$). Substituting for $`ϵ`$ some typical energies of FIBR photons ($`0.005`$ eV) one finds that the correction term is sufficient to forbid electron-positron pair production even well above $`E20`$TeV, consistently with the observations reported in Ref. .

The cosmic-ray paradox and the Markarian-501 paradox involve different energy scales and different collision processes, but admit the same type of description (they require an increase in the theory estimate of the threshold energy) and, as shown, admit a common solution based on (1). The pion-stability paradox emerging from the analysis reported in Ref.  is of a different type, since it involves a particle-decay process rather than a collision process. While for the cosmic-ray and the Markarian-501 paradoxes solutions based on the deformed dispersion relation (1) had already been discussed in the literature , previous studies of schemes leading to (1) had not noticed the associated emergence of neutral-pion increased stability, and actually had not noticed any implications for particle-decay processes. This is the main technical/theory result reported in the present note, and, remarkably, its phenomenology implications for pion decay into photons are in agreement with the indication that has emerged from the analysis reported in Ref. . Let me focus the analysis of the implications of (1) for particle decay directly on the example relevant for the pion-stability paradox: the process $`\pi \gamma +\gamma `$ (for other particle-decay processes one can of course proceed in strict analogy). My observation is based on the kinematical condition that establishes a relation between the energy $`E_\pi `$ of the incoming pion, the opening angle $`\theta `$ between the outgoing photons and the energy $`E_\gamma `$ of one of the photons (the energy $`E_\gamma ^{}`$ of the second photon is of course not independent; it is given by the difference between the energy of the pion and the energy of the first photon). This relation is found, as usual, by combining the dispersion relation, here assumed to be described by (1), with the equations for the conservation of energy and momentum. One finds

$`m_\pi ^2`$ $`=`$ $`[2E_\gamma E_\gamma ^{}+L_pE_\pi E_\gamma E_\gamma ^{}][1\mathrm{cos}(\theta )]+2L_pE_\pi E_\gamma E_\gamma ^{}`$ (4)
$`=`$ $`[2E_\gamma (E_\pi E_\gamma )+L_pE_\pi E_\gamma (E_\pi E_\gamma )][1\mathrm{cos}(\theta )]+2L_pE_\pi E_\gamma (E_\pi E_\gamma ).`$
In the $`L_p0`$ limit (the limit that corresponds to our present classical picture of space-time) this kinematical condition of course reproduces the corresponding result for conventional relativistic kinematics. The $`L_p`$-dependent correction term starts to be significant at pion energies of order $`(m_\pi ^2/L_p)^{1/3}`$. When $`E_\pi ^3>2m_\pi ^2/L_p`$ one finds that some of the values of $`E_\gamma `$ which correspond to viable decay processes according to conventional relativistic kinematics are no longer available to the decay process (in particular, since $`1\mathrm{cos}(\theta )`$ is always positive, one must exclude all values of $`E_\gamma `$ such that $`m_\pi ^22L_pE_\pi E_\gamma (E_\pi E_\gamma )<0`$). As one easily sees from (4), this reduction of the available phase space starts rather quietly (only a very small reduction) at $`E_\pi (2m_\pi ^2/L_p)^{1/3}10^{15}`$eV, but gets stronger and stronger as the pion energy increases. This picture of pion decay, with the associated depletion of the number of photons produced by ultra-high-energy neutral pions, would explain the puzzling experimental data discussed in Ref. .

Having added the pion-stability paradox to the cosmic-ray and Markaria-501 paradoxes, we now have three paradoxes that are solved by adopting the Planck-scale-deformed dispersion relation (1). In the literature one does not find any other indication of departures from a classical space-time picture and it is easy to check  that in the (huge number of) experiments that are consistent with the classical dispersion relation $`E^2=p^2+m^2`$ the correction term introduced in (1) is completely negligible (in order to compensate for the smallness of $`L_p`$ the physical context must involve unusually large hierarchies between some relevant energy/length scales, such as the ratios $`m_p/ϵ`$, $`m_e/ϵ`$ and $`E_\pi /m_\pi `$ encountered in the analysis of the paradoxes). Therefore, we seem to be confronted with the exciting perspective of having to replace the classical space-time picture with a new (quantum) picture involving the Planck length. I should stress however that, while ordinarily the combined indications of three independent experiments are considered to be sufficient for drawing definitive conclusions, the three experiments on which my analysis is based are still subject to some residual elements of doubt, and some prudence may be appropriate.

The cosmic-ray paradox is well established, in the sense that there can be no residual doubt concerning the fact that cosmic rays with energies beyond the GZK limit are being detected. There are however some alternative (not less speculative ) possible explanations of the cosmic-ray paradox, which exploit the fact that we are unable to identify the astrophysical sources of these cosmic rays and that the identification of ultra-high-energy cosmic ray as protons is still subject to (however small ) margins of uncertainty.

With respect to the Markarian-501 paradox the residual elements of doubt are complementary to the case of cosmic rays. In fact, we have a clear identification of the observed particles as photons and equally clear is the identification of Markarian 501 as the source. However, while measurements of the CMBR have become more and more accurate over the years, measurements of the FIBR have reached a satisfactory level of accuracy only very recently (see, e.g., Ref. ) and the robustness and interpretation of these recent experimental results may still be subject to further scrutiny. This is of course significant for establishing the Markarian-501 paradox, since the likelyhood that a Markarian-501 photon above threshold would reach our detectors also depends on the (density of the) FIBR.

Concerning the robustness of the pion-stability paradox, besides the need for more accurate data (some of the graphs that support the analysis reported in Ref.  show data points with significant error bars), a key reason of residual concern resides, in this author’s opinion, in the role that quantum chromodynamics (QCD) playes in the structure of the air showers produced by hadronic particles. The analysis reported in Ref.  appears to provide rather convincing evidence of the fact that within the presently-favoured phenomenological model of the relevant QCD processes (a model which has been found to be reliable in other contexts) the assumption of increased pion stability at ultra-high energies leads to improved agreement with the data on the longitudinal development of the air showers. However, certain quantitative estimates based on QCD are still rather challenging (in spite of the fact that QCD has been well understood conceptually and at the level of the formalism for several years) and it appears to be reasonable to wonder whether one should also explore the possibility of modifying the presently-favoured phenomenological model of QCD processes, without introducing an increase in pion stability.

In summary, the status of the three paradoxes is rather robust, but each of them (to different degrees) is still not completely immune from potential weaknesses. Perhaps a greater level of confidence should be attributed to the analysis here being reported by considering the consistency of the combined indications of the three experiments. In particular, from a strictly phenomenological viewpoint one could also contemplate deformed dispersion relations with the same structure of (1) but with the opposite sign choice for the correction term and/or with a deformation scale which is significantly different from the Planck scale. But for all three paradoxes the solutions require the same sign choice555With the opposite sign choice the two threshold conditions here considered would go in the opposite direction, predicting that the process is allowed at even lower energies than in the conventional theory (in clear contrast with the experimental information). As one can easily see by changing the sign in front of $`L_p`$ in (4), the opposite sign choice would also not predict the increase of pion stability here discussed., so in order to assume that the evidence emerging from these paradoxes is the result of the preliminary nature of the experimental data one should assume that the independent inaccuracies of these data have somehow conspired to point all in the same direction. Similarly, the evidence emerging from the three paradoxes also has a significant level of internal consistency for what concerns the deformation length scale. It is easy to see that the requirement of explaining the three paradoxes imposes that this length scale cannot be much smaller than $`L_p`$ (the solution of the Markarian-501 paradox is lost already by decreasing the deformation length scale by a factor of 100 or so, while the solutions of the cosmic-ray paradox and of the pion-stability paradox have a few more orders of magnitude margin). On the other hand the deformation length scale certainly cannot be much larger than $`L_p`$ because otherwise a disagreement would emerge with data at lower energies . Therefore the requirement of explaining the three paradoxes, besides imposing a consistent sign choice, also imposes that the deformation length scale be within a few orders of magnitude of the Planck length, just as one would expect in light of the quantum-space-time arguments that support (1).

While some prudence is certainly appropriate, we are clearly confronted with growing experimental evidence in favour of the exciting perspective of having to modify our present classical description of the short-distance structure of space-time. The issue will be completely settled within a few years by experiments such as the ones planned for the GLAST space telescope. As discussed in detail in Refs.  these experiments are sensitive to the implications of (1) for the structure of bursts of gamma rays that we detect from distant astrophysical sources, an effect that (since it does not involve particle-production processes) is completely independent from the effects here analyzed in association with (1). The expected sensitivity levels of GLAST (which even extend several orders of magnitude beyond the Planck-length choice of the deformation length scale) are such that (1) will be either fully confirmed or completely rejected .

If the GLAST verdict does confirm the growing evidence that is emerging from the experimental paradoxes here considered, theoretical physics will find itself in a situation that is amusingly analogous to the one that was created, a century ago, by the Michelson-Morley experiments (which can be described as providing evidence in favour of the dispersion relation $`E^2=p^2`$ for photons, and forced a revolution in the description of space-time, abandoning the Galileo-Newton picture). This will include the need to establish which of the space-time pictures that support deformed dispersion relations is actually realized in Nature. At present the fact that (1), without any free parameter, explains all observations appears to be significant. But even if (1) is indeed realized in Nature we will still have to consider two main alternative scenarios that can lead to (1). The scenario on which I focused here, which, as mentioned, is based  on some conjectured properties of space-time foam, is strongly characterized by the fact that space-time foam could have a preferred frame , just like the classical aspects of the geometry of the space-time of our Universe are such that it is possible to identify a preferred frame666At the fundamental level the theory does not have a preferred frame, but of course the field distributions that correspond to a given solution of the equations of dynamics can be used to identify a preferred frame (a frame in which these field distributions acquire a certain characteristic form.) (a convenient frame for most applications is the one in which the CMBR has the simplest properties ). The second scenario that supports (1) is based  on the opposite assumption: the quantum features of space-time would not have a preferred frame, in which case the role of $`L_p`$ in (1) should be enforced as an observer-independent condition. As mentioned, this second scenario leads  to a fundamental picture of space-time that is noncommutative , and accordingly Lorentz transformations between different frames would be governed by a quantum-algebra version of the Lorentz algebra that had emerged  as part of the mathematical-physics programme of systematic studies of quantum deformations of classical groups and algebras.

It is rather compelling that the simple first scenario, here considered, manages to explain the three paradoxes, in spite of its highly constrained structure (no free parameters). Still, it will of course be interesting to compare the paradoxes with other schemes leading to deformed dispersion relations, and particularly the mentioned scenario  that supports (1) as an observer-independent property of (noncommutative) space-time. This more delicate (both technically, because of the complexity of noncommutative geometry, and phenomenologically) analysis is postponed to future studies .