text:
# Detection of non-Gaussianity in the WMAP 1-year data using spherical wavelets

## 1 Introduction

The study of the Gaussianity of the Cosmic Microwave Background (CMB) is one of the most powerful mechanisms for understanding the nature of the primordial density fluctuations: by estimating the probability distribution from the measured temperature fluctuations of the CMB, several models for the primordial density fluctuations can be rejected or accepted at a certain significance level. For instance, standard inflationary models predict that the temperature fluctuations of the CMB correspond to a Gaussian, homogeneous and isotropic random field; whereas non-standard inflation (Linde & Mukhanov 1997, Peebles 1997, Benardeau & Uzan 2002 and Acquaviva et al. 2002) and topological defects models (Turok & Spergel 1990 and Durrer 1999) usually predict non-Gaussian random fields. Moreover, Gaussianity analyses can also be used to study the effect of secondary anisotropies produced by the reionization of the universe (Ostriker & Vishniac 1986 and Aghanim et al. 1996), the Rees-Sciama effect due to the non-linear evolution of the large scale structure (Rees & Sciama 1968 and Martínez–González & Sanz 1990), gravitational lensing (Martínez–González et al. 1997, Hu 2000 and Goldberg & Spergel 1999), etc. One of the handicaps within the detection of non-Gaussian signatures is related to the contamination of the CMB signal due to foregrounds (like the Galactic emissions and the compact sources). These components produce *spurious* non-Gaussian signals that can contaminate any intrinsic signature.

In order to uncover the nature of the primordial density fluctuations (or to conclude that the anisotropies are secondary), it is essential to have high-resolution, low-noise and large-coverage CMB data. Recently, the Wilkinson Microwave Anisotropy Probe (WMAP) NASA satellite has reported results from the 1-year all-sky data (Bennet et al. 2003a). The WMAP science team (Komatsu et al. 2003) found the data to be consistent with Gaussianity. Two statistics were used: a measure of the phase correlations of temperature fluctuations (taking into account several combinations of the bispectrum) and the Minkowski functionals; moreover, improved limits for the non-linear coupling parameter were established. Other groups have presented different analyses of Gaussianity on WMAP data: Colley & Gott (2003) have done an independent genus topology (one of the three Minkowski functionals) study, by performing a stereographic projection of the WMAP data, finding also Gaussian compatibility. Chiang et al. (2003) presented a Gaussianity study (using a statistic based on the random-phase hypothesis) of the CMB map obtained from the WMAP data by Tegmark et al. (2003). They found a significant non-Gaussian contribution at high-order multipoles. The authors conclude that this non-Gaussianity is due to unresolved foreground contamination. Gaztañaga & Wagg (2003) have performed a 3-point angular correlation function analysis, finding good agreement with the Gaussian hypothesis. Higher-order moments (up to six) of the angular correlation function were calculated by Gaztañaga et al. (2003) which show the same Gaussian agreement. Another genus-statistic study has been done by Park (2003) using the foreground cleaned frequency maps of the WMAP data; a stereographic projection is also used finding a non-Gaussian detection by comparing the difference of the negative and positive genus —with respect to the Gaussian prediction— at different threshold levels. Several sources are considered in order to explain this detection, but higher signal-to-noise ratio data is needed to distinguish among them. Eriksen et al. (2003) have computed 2- and 3-point correlations finding north/south asymmetries in the WMAP data. These asymmetries appear at large and intermediate scales and are associated with a lack of structure in the northern hemisphere, while the southern one is compatible with the Gaussian model (additional asymmetries east/west are also found)111During the refereeing process new works have appeared in astro-ph like Eriksen et al. 2004..

In this work we present a wavelet analysis of the Gaussianity of the WMAP 1-year all-sky map performed in wavelet space. Wavelets have been extensively used in the analysis of the CMB data, not only concerning Gaussianity studies, but also in the component separation field (Tenorio et al. 1999, Cayón et al. 2001, Vielva et al. 2001a, Vielva et al. 2001b, Vielva et al. 2003 and Stolyarov et al. 2003) and in denoising techniques (Sanz et al. 1999a,b). Regarding the application of wavelets to Gaussianity studies, Hobson et al. (1999) have shown that the wavelet coefficients provide a significantly better detection of non-Gaussian features due to cosmic strings than the Minkowski functionals. Similarly, Aghanim et al. (2003) have shown that wavelets are more efficient to detect certain non-Gaussian features than other typical estimators like those based in Fourier analysis (bispectrum and trispectrum). Even more, due to the special nature of wavelets, the non-Gaussian sources can not only be detected but also identified in the maps, since the spatial information is kept in the wavelet space. Moreover, as it was shown in Barreiro & Hobson (2001) and in Martínez–González et al. (2002), the optimal combination of the information given by the wavelet coefficients at different scales, highly increases the power of wavelets to detect weak non-Gaussian signals.

The first application of wavelet techniques to detect non-Gaussian signatures was done by Pando et al. (1998), where the 2D flat Daubechies-4 wavelet was used to analyse the COBE-DMR data (see Mukherjee et al. 2000 for a critical review of this work). Lately, several works have used wavelets to study the Gaussianity of the COBE-DMR maps in the HEALPix scheme (Górski et al. 1999): Barreiro et al. (2000) have used the Spherical Haar Wavelet (SHW), whereas Cayón et al. (2001, 2003) have used the Spherical Mexican Hat Wavelet (SMHW). A critical comparison between the performances of both spherical wavelets for detecting non-Gaussian features (skewness and kurtosis) was done by Martínez–González et al. (2002).

The aim of this paper is to apply the technique developed in Cayón et al. (2001) and Martínez–González et al. (2002) to the 1-year all-sky WMAP data in order to detect non-Gaussian features. It is organized as follows. In Section 2 the process to generate the simulations and to reduce the data is explained. In Section 3 the SMHW is briefly described and the different statistics are defined. The SMHW analysis is presented in Section 4, and the possible sources of the detected non-Gaussianity are discussed in Section 5. Finally the conclusions are given in Section 6.

## 2 The WMAP data reduction and WMAP-like simulations

The NASA Wilkinson Microwave Anisotropy Probe (WMAP) satellite was launched in the summer of 2001 and in February 2003 the 1-year results were presented. The WMAP radiometers observe at 5 frequency bands: K-band (22.8 GHz, 1 receiver), Ka-Band (33.0 GHz, 1 receiver), Q-Band (40.7 GHz, 2 receivers), V-Band (60.8 GHz, 2 receivers) and W-Band (93.5 GHz, 4 receivers). All the papers, data and products generated by the WMAP team can be found in the *Legacy Archive for Microwave Background Data Analysis* (LAMBDA) web page 222http://cmbdata.gsfc.nasa.gov/. The noise and beam properties can also be found in the LAMBDA web page. The WMAP maps are presented in the Hierarchical, Equal Area and iso-Latitude Pixelization (HEALPix, Górski et al. 1999) at the $`N_{side}=512`$ resolution parameter. The number of pixels is given by $`12N_{side}^2`$.

In order to perform the Gaussianity study, the WMAP team (Komatsu et al. 2003) suggests combining all the maps produced by the receivers where the CMB is the dominant signal (Q-Band, V-Band and W-Band). We have followed the recommendation given by the WMAP team in the choice of the combined map used in the analysis. We would like to point out that none of the previous Gaussianity studies of the WMAP data have used this combined map (except the analysis performed by the WMAP team itself).

A very simple pipeline is given in Bennett et al. (2003b). At a given position in the sky ($`𝐱`$), the temperature is given by:

$$T(𝐱)=\underset{j=3}{\overset{10}{}}T_j(𝐱)w_j(𝐱),$$
(1)
where the indices $`j=3,4`$ refer to the Q-Band receivers, $`j=5,6`$ to the ones of the V-Band and, finally, $`j=7,8,9,10`$ correspond to the receivers of the W-Band (the indices $`j=1,2`$ are used for the K and Ka receivers, respectively). The noise weight $`w_j(𝐱)`$ is defined by:

$`w_j(𝐱)={\displaystyle \frac{\overline{w}_j(𝐱)}{_{j=3}^{10}\overline{w}_j(𝐱)}},`$ $`\overline{w}_j(𝐱)={\displaystyle \frac{N_j(𝐱)}{\sigma _{0}^{}{}_{j}{}^{}{}_{}{}^{2}}}`$ (2)
where $`\sigma _{0}^{}{}_{j}{}^{}`$ is the noise dispersion per observation and $`N_j(𝐱)`$ is the number of observations made by the receiver $`j`$ at position $`𝐱`$ (see Bennett et al. 2003a). Equation (1) provides a single map where the signal-to-noise has been increased. Although the CMB dominated frequencies have been chosen, we still have significant contribution due to Galactic foregrounds (thermal dust, free-free and synchrotron) as well as extragalactic point sources (negligible Sunyaev-Zel’dovich contribution due to galaxy clusters is expected, see Bennett et al. 2003b). In order to avoid the Galactic emissions, the WMAP team performed a foreground template fit described in Section 6 of Bennett et al. 2003b. The 94 GHz dust map of Finkbeiner et al. (1999) is used as the thermal dust template; the $`H_\alpha `$ map of Finkbeiner (2003) corrected for extinction through the $`E_{BV}`$ map of Schlegel et al. (1998) is used as the free-free template; finally, the synchrotron template is the 408 MHz Haslam et al. (1982) map. Hence, Equation (1) is modified by:

$$\widehat{T}(𝐱)=\underset{j=3}{\overset{10}{}}\widehat{T_j}(𝐱)w_j(𝐱),$$
(3)
where $`\widehat{T_j}(𝐱)`$ is the temperature at position $`𝐱`$ for the receiver $`j`$ after foreground correction. The parameters for the best-fit of the foreground templates are given in Bennett et al. (2003b). Even more, the foreground cleaned frequency maps can be found in the LAMBDA site.

The map is then degraded to resolution $`N_{side}=256`$ —since the very small scales are dominated by noise— and a mask is applied to avoid the contamination due to the strong emission at the Galactic plane and the contribution due to known radio point sources. This mask, that can also be found in the LAMBDA site, is called *Kp0* and keeps 76.8% of the sky. Finally, the residual monopole and dipole outside the mask are removed. The final map is shown in Figure 1.

In order to study the Gaussianity of this map, we have produced 10000 Gaussian simulations. Using CMBFAST (Seljak & Zaldarriaga 1996), we have calculated the $`C_{\mathrm{}}`$ given by the cosmological parameters estimated by the WMAP team (Table 1 of Spergel et al. 2003). Random Gaussian $`a_\mathrm{}m`$ of CMB realizations have been generated and convolved at each one of the WMAP receivers with the adequate beams. After the transformation from harmonic to real space, uncorrelated Gaussian noise realizations have been added following the number of observations per pixel ($`N_j(𝐱)`$) and the noise dispersion per observation ($`\sigma _{0}^{}{}_{j}{}^{}`$). We have combined all the maps following Equation (3). Finally, the 10000 simulations have been degraded to $`N_{side}=256`$, the *Kp0* mask has been applied, and the residual monopole and dipoles have been fitted and subtracted for each simulation independently.

## 3 The tool and the non-Gaussian estimators

Wavelet techniques have been shown to be very powerful for detecting non-Gaussianity in CMB data (Hobson et al. 1999, Aghanim et al. 2003). Due to the special nature of wavelets, a multi-scale study can be performed to amplify the signature of the non-Gaussian features dominating at a given scale. Moreover, the SMHW is ideal for the enhancement of non-Gaussian signatures with spherical symmetry. The SMHW has been already applied for non-Gaussian studies to the COBE-DMR data (Cayón et al. 2001, 2003) and to *Planck* simulations (Martínez-González et al. 2002). The SMHW can be obtained from the Euclidean Mexican Hat Wavelet (MHW) following the stereographic projection suggested by Antoine & Vanderheynst (1998). This projection ensures that the wavelet properties are kept and that the MHW is recovered in the small angle limit (see Martínez-González et al. 2002 for a graphical explanation of this extension). The SMHW satisfies the *compensation*, *admissibility* and *normalization* properties that define a wavelet and is given by:

$$\mathrm{\Psi }_S(y,R)=\frac{1}{\sqrt{2\pi }N(R)}\left[1+\left(\frac{y}{2}\right)^2\right]^2\left[2\left(\frac{y}{R}\right)^2\right]e^{y^2/2R^2},$$
(4)
where $`R`$ is the scale and $`N(R)`$ is a normalization constant:

$$N(R)R\left(1+\frac{R^2}{2}+\frac{R^4}{4}\right)^{1/2}.$$
(5)
The distance on the tangent plane is given by $`y`$ that is related to the polar angle ($`\theta `$) through:

$$y2\mathrm{tan}\frac{\theta }{2}.$$
(6)
At a given scale ($`R`$), several statistics can be defined. In this work we have used two simple non-Gaussian estimators: skewness ($`S(R)`$) and kurtosis ($`K(R)`$):

$`S(R)`$ $`=`$ $`{\displaystyle \frac{1}{N_R}}{\displaystyle \underset{i=1}{\overset{N_R}{}}}w_i(R)^3/\sigma (R)^3`$ (7)
$`K(R)`$ $`=`$ $`{\displaystyle \frac{1}{N_R}}{\displaystyle \underset{i=1}{\overset{N_R}{}}}w_i(R)^4/\sigma (R)^43,`$ (8)
where $`N_R`$ is the number of coefficients at scale R and $`\sigma (R)`$ is the dispersion of the wavelet coefficients at the scale R ($`w_i(R)`$):

$`\sigma ^2(R)`$ $`=`$ $`{\displaystyle \frac{1}{N_R}}{\displaystyle \underset{i=1}{\overset{N_R}{}}}w_i(R)^2.`$ (9)
In the previous $`S(R)`$, $`K(R)`$ and $`\sigma (R)`$ definitions, it is assumed that the $`N_R`$ wavelet coefficients ($`w_i(R)`$) have zero mean at each scale $`R`$.

## 4 The Spherical Wavelet analysis

In order to test the Gaussianity of the WMAP 1-year data, we have applied the following analysis using the SMHW.

We have performed 10000 simulations following the pipeline proposed by Bennett et al. (2003b) and already indicated in Section 2. Each one of these simulations have been convolved with the SMHW at different scales ($`R_1=13.7`$, $`R_2=25`$, $`R_3=50`$, $`R_4=75`$, $`R_5=100`$, $`R_6=150`$, $`R_7=200`$, $`R_8=250`$, $`R_9=300`$, $`R_{10}=400`$, $`R_{11}=500`$, $`R_{12}=600`$, $`R_{13}=750`$, $`R_{14}=900`$ and $`R_{15}=1050`$ arcmin). Acceptance intervals at certain significance levels $`\alpha `$ (e.g., 32%, 5% and 1%) have been established at each scale based on these simulations. These acceptance intervals are defined as the intervals which contain a probability $`1\alpha `$ and the remaining probability is the same above and below the interval, i.e. $`\alpha /2`$ at each side. The acceptance intervals have been determined by studying the distribution of the skewness and the kurtosis at each scale independently and calculating the corresponding percentiles. We have checked that the number of simulations performed (10000) is enough to establish those acceptance intervals with good precision. Finally, the same analysis has been applied to the WMAP map ($`\widehat{T}(𝐱)`$) plotted in Figure 1.

Since we are convolving a map with a large Galactic mask plus known point sources (the so-called *Kp0* mask) with the SMHW, we are introducing in the wavelet analysis a large number of pixels —those near the border of the mask— largely affected by the zero value of the mask. This effect is taken into account through the simulations, but, obviously, it is introducing an undesirable loss of efficiency. In order to avoid this effect, we will not consider pixels whose SMHW coefficients have a strong contamination from the mask. In other words, we *exclude* these pixels after the SMHW convolution. Since we are convolving the combined WMAP map ($`\widehat{T}(𝐱)`$) with different scales, the number of pixels to be excluded is growing with the scale. There could be several criteria to define this set of *exclusion masks*. We have tested several definitions and all of them lead to the same results. The results here presented were obtained with an *exclusion mask* defined as follows: at a given scale ($`R`$), the *exclusion mask* $`M(R)`$ is composed by the *Kp0* mask plus those pixels that are closer than $`2.5R`$ to any one of the pixels of the *Kp0* mask that are in the Galactic plane. In other words, at a given scale, the *exclusion mask* is an enlargement (proportional to the scale) of the *Kp0* mask from the Galactic plane (but note that the mask is not increased around the masked point sources). The *exclusion masks* corresponding to the 15 scales are plotted in Figure 2.

Following this process we obtain the results presented in Figure 3: an excess of kurtosis is detected at two consecutive scales: $`R_8=4.17^{}`$ and $`R_9=5^{}`$. The value of the kurtosis at $`R_8=4.17^{}`$ is a non-Gaussian detection outside the acceptance interval at the $`1\%`$ significance level. Even more, only 40 of the 10000 simulations present a kurtosis value larger or equal than the one detected in the combined WMAP map at this scale. This corresponds to a right tail probability of $`0.4\%`$. A similar result is found at scale $`R_9=5^{}`$ for which the right tail probability is also $`0.4\%`$.

We want to remark that this value is obtained for very different definitions of the *exclusion mask*. For instance, we have tested an enlargement of the *Kp0* mask from the Galactic plane scaling with the SMHW scale ($`R`$) as: $`0.7R`$, $`0.8R`$, $`0.9R`$, $`1.0R`$ and $`1.75R`$. Moreover, an enlargement of the whole *Kp0* mask not only from the Galactic plane but also from all the masked point sources has been also checked, leading to similar results.

In addition we have performed the SMHW analysis in each hemisphere to determine if the non-Gaussian detection is present in all the sky. Previous non-Gaussianity works on the WMAP data have reported some asymmetries between the northern and the southern hemisphere (Park 2003, Eriksen et al. 2003). As Figure 4 clearly shows, the kurtosis is located in the southern hemisphere. The non-Gaussian detection in the southern hemisphere occurs at the $`R_7=3.33^{}`$ and $`R_8=4.17^{}`$ SMHW scales. The value of the kurtosis at scale $`R_7=3.33^{}`$ is above the acceptance interval at the $`1\%`$ significance level. Only 11 of the 10000 simulations present values for the kurtosis equal or larger than the one detected at this scale, which corresponds to a right tail probability of $`0.1\%`$. A similar result is found for the scale $`R_8=4.17^{}`$ for which the right tail probability is $`0.2\%`$.

By looking at the scalogram of the wavelet coefficients in the northern and southern hemispheres (dispersions at the different scales $`R`$ in Figure 4), we can see that the northern hemisphere shows less structure than the southern one, at scales between $`4^{}`$ and $`12^{}`$. On the contrary, the northern hemisphere has more structure than the southern one at larger scales. However, the scalograms of both hemispheres are within the acceptance intervals of the Gaussian model.

Let us remark that, as mentioned before, previous works have reported asymmetries between the northern and southern hemispheres. For instance, Park (2003) has found different north/south genus behaviour which implies a non- Gaussian detection at the $`99\%`$. Eriksen et al. (2003) also found an asymmetric behaviour north/south regarding the structure level at large and intermediate angular scales: the northern hemisphere has a lack of structure whereas the southern one is compatible with the Gaussian model. However, our study reveals: 1) a direct measurement of non-Gaussian signatures related to the kurtosis of those structures with a typical size in the sky of around $`10^{}`$, with a probability $`0.4\%`$; 2) the southern hemisphere shows an excess of kurtosis with a probability $`0.2\%`$ (around the previous scale) whereas the northern one is compatible with the Gaussian model.

Finally, we would like to point out that we have also performed a wavelet analysis based on the Spherical Haar Wavelet (SHW) as it was proposed in Martínez- González et al. (2002). We did not find any non-Gaussian detection using the SHW. This is not so surprising, since (as it was shown in Martínez- González et al., 2002) the SHW is less efficient than the SMHW to detect non- Gaussian signals.

## 5 Discussion: Sources of non-Gaussianity

The detection of non-Gaussianity in the WMAP map ($`\widehat{T}(𝐱)`$) was presented in the previous Section. The kurtosis of the wavelet coefficients at scales $`R4^{}`$ was found to have a very low probability $`0.4\%`$. The non-Gaussianity detection is localized in the southern hemisphere.

The aim of this Section is to study, as far as we can, the source of this non- Gaussianity. Let us summarise the major hypotheses that have been assumed:

1. The CMB is a homogeneous and isotropic multivariate Gaussian random field on the sphere
2. The WMAP data are free from systematics
3. The WMAP data outside the Kp0 mask are not contaminated by foregrounds
4. The uncertainties in the cosmological parameters have a negligible effect on the results

Obviously, the most critical hypothesis is 1. If we find that hypotheses 2 to 4 have a negligible effect on the results, then we could conclude that hypothesis 1 is not verified. In fact, hypothesis 1 contains at least 2 hypotheses:

* The CMB is a homogeneous and isotropic random field on the sphere
* The CMB is a multivariate Gaussian random field (or equivalently the $`a_{lm}`$ are Gaussian)

In practice it is very difficult to distinguish between these two hypotheses. On the one hand, since we are using non-Gaussian estimators, we prefer to refer to a non-Gaussian detection which is produced by localised regions. More precisely, we have found a non-Gaussian detection that seems to be present in the southern hemisphere rather than in the northern one. On the other hand, homogeneity and isotropy is a more fundamental principle in cosmology than Gaussianity.

Systematic effects are a possible source of non-Gaussianity. In Subsection 5.1, we will study if there are any systematic effects: are any of the receivers producing the non-Gaussian feature?, or the noise?, or the beam?. In Subsection 5.2 we present tests to study the possible influence of the Galactic foregrounds on the detection. The influence of the accuracy of the estimated power spectrum of the WMAP data (or of the derived cosmological parameters) on our results is also studied (Subsection 5.3). Finally, in Subsection 5.4, we speculate about intrinsic fluctuations.

### 5.1 Systematics

Our first test to check the influence of systematics in the non-Gaussian detection consists on looking for any *rare* receiver. Instead of analysing the WMAP map ($`\widehat{T}(𝐱)`$) obtained by the weighted average of all the Q-Band, V-Band and W-Band receivers, we perform the SMHW analysis in each receiver map ($`\widehat{T}_{Q_1}(𝐱)`$, $`\widehat{T}_{Q_2}(𝐱)`$, $`\widehat{T}_{V_1}(𝐱)`$, $`\widehat{T}_{V_2}(𝐱)`$, $`\widehat{T}_{W_1}(𝐱)`$, $`\widehat{T}_{W_2}(𝐱)`$, $`\widehat{T}_{W_3}(𝐱)`$ and $`\widehat{T}_{W_4}(𝐱)`$) independently. If the results are the same for these maps, then we can conclude that the detection of non-Gaussianity is not produced by any particular receiver.

Results are presented in Figure 5. The $`K(R)`$ curve obtained for the WMAP map ($`\widehat{T}(𝐱)`$) is plotted together with the curves obtained for each of the receivers ($`K_{Q_1}(R)`$, $`K_{Q_2}(R)`$, $`K_{V_1}(R)`$, $`K_{V_2}(R)`$, $`K_{W_1}(R)`$, $`K_{W_2}(R)`$, $`K_{W_3}(R)`$, $`K_{W_4}(R)`$). We can see how the pattern for the $`K(R)`$ curve is almost perfectly followed by the other kurtosis curves.

We have also checked the influence of possible systematics related to instrumental features (noise and beams). We have performed the SMHW analysis on maps produced from the subtraction of receivers at the same frequency. These maps are almost free from CMB and foreground contribution (remaining residuals can be present due to the slightly different resolutions between the two subtracted channels). In particular, we have analysed the maps obtained by subtracting the two receivers in the Q-band ($`Q_1Q_2`$), the two receivers in V-Band ($`V_1V_2`$) and the map obtained from $`W_1W_2+W_3W_4`$. As seen in Figure 6, the patterns followed by the kurtosis in these three cases are compatible with simulations assuming a Gaussian CMB and performing the same operations as for the real data at each frequency (for each case, 1000 simulations were generated to establish the acceptance intervals at the 32%, 5% and 1% significance levels). In addition, they are completely different from the $`K(R)`$ curve obtained for the WMAP map ($`\widehat{T}(𝐱)`$).

Hence, these tests seem to indicate that systematic effects do not have a significant role in the non-Gaussian signature.

### 5.2 Foregrounds

Foregrounds have been considered as the next possible source of non-Gaussianity. We do not consider emission from point sources because the brightest radio sources have been previously masked and, moreover, the angular scale of the non-Gaussian detection is much larger than the WMAP angular resolution. As regards the expected Sunyaev-Zel’dovich effect contribution due to galaxy clusters, it has been shown to be almost negligible (Bennett et al. 2003b).

In order to study the frequency dependence of the non-Gaussian detection we have independently analysed the combined frequency maps: $`\widehat{T}_Q(𝐱)`$, $`\widehat{T}_V(𝐱)`$ and $`\widehat{T}_W(𝐱)`$. For the frequency range of Q, V and W (from 41 GHz to 94 GHz), the non-Gaussian pattern found for the kurtosis does not resemble any of the frequency dependence due to the Galactic foregrounds (synchrotron, free-free and thermal dust).

We have also studied the clearly foreground contaminated K and Ka channels in this test. The pattern of the kurtosis curve for these channels does not completely follow the one for the other WMAP channels. Whereas the peak around $`R_8`$ clearly appears, the behaviour is very different at scales below $`R_6`$. In addition, there is an offset along all the scales. We think that this effect corresponds to the very high foreground contamination that these channels show. In order to clarify this point, we have done the following exercise. We have generated an additional CMB map by subtracting the Ka map (at 33 GHz) from the K one (at 23 GHz), multiplying the first one by a factor of 2.65. This number corresponds to the expected increment of the synchrotron emission from 33 GHz to 23 GHz333A power law is assumed for the frequency dependence of the synchrotron emission: $`T_{syn}(\nu )T_{syn}(\nu _0)\left(\nu /\nu _0\right)^{2.7}`$, as it was proposed by Bennet et al. (2003b).. As seen in Figure 7, the pattern of the kurtosis curve for this map is the same as the one detected for the combined WMAP map ($`\widehat{T}(𝐱)`$). Therefore, the same pattern for the kurtosis is found, not only for the Q, V and W bands, but also for the whole WMAP frequency range (from 23 GHz to 94 GHz). This implies not only that the synchrotron does not seem to be the source of the non-Gaussian signature, but also that this detection does not show any significant frequency dependence.

An additional test to check the influence of the foregrounds was done by analysing a map where the CMB contribution is small: the two receivers of the Q- Band and the two ones of the V-Band are subtracted from the four receivers of the W-Band. In this map we have significant contributions from the foregrounds and the noise. The SMHW analysis of this map also shows (Figure 8) that, apparently, the foregrounds are not the source of the non-Gaussian detection since they are compatible with their own Gaussian bands. Even more, the pattern of the kurtosis curve is completely different from the one estimated from the combined WMAP map ($`\widehat{T}(𝐱)`$).

Finally, we have done the following exercise. We have estimated the foreground correction ($`F(𝐱)`$) subtracted from the data by the WMAP team just by doing: $`F(𝐱)=T(𝐱)\widehat{T}(𝐱)`$. This foreground correction is added to a Gaussian simulation of the combined WMAP map. The SMHW analysis was performed in this contaminated map. As shown in Figure 9, the curve for the kurtosis ($`K(R)`$) not only has a pattern completely different from the one detected in the combined WMAP map ($`\widehat{T}(𝐱)`$) but is also very similar to the uncontaminated Gaussian simulation. Even more, the foreground correction ($`F(𝐱)`$) was added twice to the Gaussian simulation, finding the same behaviour. This exercise —as the previous ones— also seems to discard the Galactic foreground emissions as the ones responsible for the non-Gaussian detection.

### 5.3 Influence of the uncertainties in the cosmological parameters

Since we use CMB Gaussian simulations following the power spectrum generated by CMBFast, we must also check how the uncertainties in the cosmological parameters (and, hence, in the power spectrum) could affect the results. In particular, how the 1% significance bands depend on these uncertainties? An exhaustive analysis would require a very sophisticated statistical framework that can accomodate uncertainties of cosmological parameter estimates. Instead we have adopted two simplified approaches that, nevertheless, we believe are good indicators of the effect on our results from the uncertainties in the $`C_{\mathrm{}}`$’s.

First, we have performed Gaussian simulations with three different power spectra, being all of them compatible with the $`1\sigma `$ errors of the WMAP estimated power spectrum. In Figure 10 the power spectrum given by CMBFast using the best-fit Cosmological Model of the WMAP team (solid red line) and the power spectrum estimation for the WMAP (cyan) have been plotted. The three additional power spectra are defined as follows: the so-called ’upper-limit’ and ’lower-limit’ power spectra are created from the best-fit power spectrum by adding (or subtracting) the $`1\sigma `$ error estimated by the WMAP team; the so- called ’zig-zag’ power spectrum is an intentionally exotic power spectrum that oscillates around the best-fit one, always within the ’upper-limit’ and ’lower- limit’ power spectra. By analysing Gaussian simulations following these power spectra, we can estimate the influence of the uncertainties of the power spectrum determination on the acceptance intervals at the $`1\%`$ significance level already established for the best-fit WMAP model. This is plotted in Figure 11. It can be seen that for small and intermediate scales (including the one at which the non-Gaussian detection is found), the uncertainties in the 1% significance bands are negligible. Only at large scales we can find some deviations (especially for the zig-zag power spectrum).

The second approach to test the influence of the uncertainties in the cosmological parameters has been done by whitening the data as well as the simulations. By doing this, most of the influence of the input power spectrum is eliminated. In other words, we also detect non-Gaussianity wherever the original power spectrum was. In this case, we also obtain a non-Gaussian feature with the same probability around the same scales.

### 5.4 Is it due to Intrinsic fluctuations?

The previous Subsections have shown that there are not strong evidences for the non-Gaussian detection due to neither systematic effects nor Galactic foreground emissions nor uncertainties in the cosmological parameters. Therefore, other sources for this detection may be considered.

In particular, intrinsic fluctuations can not be rejected. One of the powerful advantages of the wavelet analysis is the possible identification, on the sky, of the non-Gaussian source. In Figure 12 (left panel) the SMHW coefficients at the scale $`R_8`$ are plotted. The minimum value of this map is $`4.57\sigma (R_8)`$, where $`\sigma (R_8)`$ is the dispersion of the wavelet coefficients at scale $`R_8`$. Using the 10000 Gaussian simulations, we have checked that the probability for an extreme like this is $`1\%`$. The right panel of Figure 12 shows those SMHW coefficients that are above (in absolute value) $`3\sigma (R_8)`$. If these pixels are not taken into account in the kurtosis calculation, the estimation for this statistic is perfectly compatible with the Gaussian simulations. In fact, we have computed the acceptance intervals at the 32%, 5% and 1% significance levels for additional sets of 10000 Gaussian simulations, where those pixels with values above a given threshold are not taken into account in order to compute the kurtosis curve $`K(R)`$. The thresholds are: $`3.0\sigma _w(R)`$, $`3.5\sigma _w(R)`$, $`4.0\sigma _w(R)`$ and $`4.5\sigma _w(R)`$. These acceptance intervals have been used to check if the combined WMAP map satisfies the Gaussian hypothesis for each one of these different thresholds. The results are presented in Figure 13. It is clear that, as the threshold decreases, the combined WMAP data is more compatible with Gaussianity. In particular, for the $`3\sigma _w`$ threshold the WMAP data is compatible with the Gaussian hypothesis (within the acceptance interval at the 5% significance level).

A very interesting exercise was done in order to study the previous very cold spot (in SMHW space) in the southern hemisphere ($`b=57^{},l=209^{}`$). We have calculated the mean value of the cold spot (at $`R_8`$) at each one of the WMAP frequency channels. As shown in Figure 14, the frequency dependence of this spot follows the CMB one, and is not compatible with any of the Galactic foregrounds. Notice that if the 2.65 times Ka map is subtracted from the K one, we get an even better agreement with constant frequency dependence, given by the dot-dashed line. We have also studied the frequency dependence of this cold spot in real space (after filtering all the frequency maps with a Gaussian beam of FWHM = $`8^{}`$, the size corresponding to the SMHW scale $`R_8`$). A good agreement with the CMB behaviour is found for bands Q, V and W, whereas a clear contamination due to synchrotron appears between bands K and Ka. However, when the 2.65 times Ka map is subtracted from the K one, the CMB behaviour is recovered. We notice that whereas for the K and Ka bands in real space the synchrotron emission dominates over CMB, this is not the case for K and Ka bands in wavelet space, since the SMHW has diminished the synchrotron amplitude relative to the CMB one (due to the large scale variation of the former). Hence, the frequency analysis performed on the cold spot, also indicates that the foregrounds do not seem to be the source of the non-Gaussianity detection. This cold spot appears as a real structure in the sky, not only seen by WMAP, but also by COBE-DMR. If the COBE-DMR map is convolved with the SMHW at $`R_8`$, we find some pixels above $`3\sigma _w`$ around this point.

There are several possibilities which can explain the non-Gaussian detection for the kurtosis of the SMHW coefficients. For instance, massive superstructures (like superclusters or Great Attractor-like structures) or large voids, can produce secondary anisotropies through the Rees-Sciama effect (Martínez-González & Sanz 1990). Even more, also primordial anisotropies can be considered. Non standard inflationary models could also explain the detected non-Gaussianity (see e.g. Acquaviva et al. 2002 and references therein). In addition, some kind of topological defects like monopoles or textures (Turok & Spergel 1990) could be present in the sky. On the contrary, non isotropic topological defects like cosmic strings may not produce this kind of non-Gaussianity, since their characteristic scale is around arcminutes and also the detection has been done with an isotropic wavelet.

## 6 Conclusions

We have presented the detection of non-Gaussianity in the WMAP 1-year data, in the kurtosis of the SMHW coefficients at scales around $`4^{}`$, which implies a size in the sky of around $`10^{}`$. At those scales, the kurtosis values corresponding to the WMAP combined map are outside the acceptance interval at the 1% significance level, with a right tail probability of $`0.4\%`$. The Gaussianity study was also done on both hemispheres independently, showing that the northern hemisphere is compatible with the Gaussian model, whereas the southern one presents a very clear non-Gaussian signal similar to the one detected in the all-sky analysis. This detection has a right tail probability of $`0.1\%`$ The number of performed simulations has been checked to be enough to establish accurately the acceptance interval at the 1% significance level for the kurtosis values at the different SMHW scales studied.

The independent analyses performed on each hemisphere also show that the northern hemisphere seems to have less structure than the southern one at scales between $`4^{}`$ and $`12^{}`$, whereas the southern one has less structure than the northern hemisphere at larger scales. However, both scalograms, the one for the north and the one for the south, seem to be compatible with the Gaussian model. Our non-Gaussian detection differs from other asymmetries reported up to date (Park 2003 and Eriksen et al. 2003), since a direct non-Gaussian signature has been found in the kurtosis of the wavelet coefficients at scale $`R10^{}`$.

Since those SMHW pixels clearly affected by the mask decrease the efficiency of the analysis, only those pixels with a small contribution from the Galactic mask (*Kp0*), at each SMHW scale, have been analysed. However, the non-Gaussian detection is quite insensitive to the particular choice of the *exclusion masks*, since we have tested that it appears outside the acceptance intervals at the same significance level for several choices.

We have performed several tests in order to identify the source of this non-Gaussian signature. Systematic effects and foregrounds have been carefully studied. The different channels have been checked, showing that all the CMB dominated ones (Q, V and W) have the same pattern for the kurtosis of the SMHW coefficients. Moreover, by analysing the maps produced by subtracting the data given by the receivers at each frequency —these maps are almost free of CMB and foreground contributions— we have shown that the non-Gaussian signal disappears. The non-Gaussian signal appears in the whole WMAP frequency range (from 23 GHz to 94 GHz) showing no frequency dependence. Even more, we added overestimated foreground contamination to Gaussian CMB maps, showing that the detected non-Gaussian signal does not appear in the SMHW analysis. This seems to indicate no correlation with Galactic foregrounds. Since the brightest point sources are masked and the non-Gaussian features appear at intermediate scales, we do not expect a significant contribution coming from this emission. The Sunyaev-Zel’dovich effect due to galaxy clusters is negligible for the WMAP 1-year data (only the most prominent local clusters like Coma can contribute at the $`2\sigma `$ level).

The uncertainties in the power spectrum estimation (and, hence, in the cosmological parameters) have also been taken into account, showing that the 1% significance bands for the kurtosis (at the scales of interest) are not affected by such uncertainties. Three additional power spectra covering the $`1\sigma `$ error band in the power spectrum determination were used. Even more, an additional test was done by performing a whitening of the maps. The non-Gaussian signal is also detected in this case. This gives robustness to the result, since possible errors in the estimation of the $`C_{\mathrm{}}`$ of the Gaussian model do not affect it.

Thanks to the properties of the SMHW, we can identify (in the sky) the possible source of the non-Gaussianity. We have found a cold spot in the southern hemisphere (in $`b=57^{},l=209^{}`$) with a probability $`1\%`$. The frequency dependence of this spot is compatible with the one of the CMB. A more detailed study of extrema will be presented in a future work. This study should be done in both real and wavelet spaces and could be very useful in order to clarify if the non-Gaussian signature is localized in certain positions on the sky or otherwise, it is produced by a non-Gaussian temperature distribution.

Finally, taking into account all the performed tests, intrinsic fluctuations (due to secondary anisotropies like the Rees-Sciama effect or due to primordial fluctuations produced by topological defects or non-standard inflation) cannot be rejected. The WMAP 2-year data will be, undoubtedly, very useful to confirm or not these results.

Authors kindly thank Dr. Eiichiro Komatsu for very useful comments regarding the data reduction and suggestions and Julio Edgar Gallegos for useful discussions. PV thanks Universidad de Cantabria (UC) for a postdoc grant. RBB thanks UC and Spanish Ministerio de Ciencia y Tecnología (MCYT) for a Ramón y Cajal contract. We acknowledge partial financial support from the Spanish MCYT project ESP2002- 04141-C03-01. We kindly thank IFCA for providing us with its Grid Wall cluster to generate and analyse the WMAP simulations. We also thank Centro de Supercomputación de Galicia (CESGA) for providing the Compaq HPC320 supercomputer to run part of the analysis. We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA). Support for LAMBDA is provided by the NASA Office of Space Science. This work has used the software package HEALPix (Hierarchical, Equal Area and iso-latitude pixelization of the sphere, http://www.eso.org/science/healpix), developed by K.M. Gorski, E. F. Hivon, B. D. Wandelt, J. Banday, F. K. Hansen and M. Barthelmann. We acknowledge the use of the software package CMBFAST (http://www.cmbfast.org) developed by U. Seljak and M. Zaldarriaga.


format: text

text:
# Optimizing omnidirectional reflection by multilayer mirrors

## 1 Introduction

At optical frequencies, metallic mirrors reflect strongly for any angle of incidence and any polarization. However, they display dissipative losses, which constitutes a drawback in practical applications. Photonic crystals  were originally proposed by Yablonovitch  to solve this problem: they are periodically microstructured materials that reflect in stop (forbidden) bands within which light propagation is not possible in an infinite structure. Since they are made from transparent materials, photonic crystals can be almost free of losses at any prescribed frequency.

In the one-dimensional case, a photonic crystal is nothing other than a periodic multilayer. Although much attention has been payed to dielectric Bragg mirrors consisting of alternating low- and high-index layers, certain aspects of the reflection by periodic layered media are universal: for $`N`$ periods the reflectance goes to unity as $`N^2`$ at the band edges, while tends to unity exponentially at the band gaps .

This means that, in practice, not very many periods are needed to have a stop band. One is then led to consider stacks of $`N`$ periods (which are often called finite periodic structures) and apply to them conditions that are valid only when the system is strictly infinite. Recently , we have put forward this problem and provided an alternative framework for dealing with these finite periodic structures: the trace of the basic period allows us to classify them into three classes with quite different properties.

Concerning the performance of these structures, it is indisputable that Bragg quarter-wave stacks (designed for normal incidence) are the most thoroughly studied in connection with omnidirectional reflection . In spite of this, the current interest in extreme ultraviolet  and soft x-ray  optics is driving a great deal of work on new methods for optimizing the design of multilayer mirrors. In addition to the simple but cumbersome optimization by eye, only recently more sophisticated techniques have been started to be used : relevant examples include the downhill simplex algorithm , the systematic search in the parameters space , the simulating annealing , the needle variation technique , the Levenberg-Marquardt algorithm , or genetic algorithms . The aim of this paper is to provide an alternative optimization criterion that has the virtue of using a simple analytical figure of merit with a very clear physical meaning.

## 2 Notations and general relations

We start by examining the properties of the basic period of our structure, which is assumed to be lossless. The field amplitudes at each side (called ambient and substrate) of the unit cell are related by the $`2\times 2`$ complex transfer matrix $`𝖬_{as}`$ that can be expressed as

$$𝖬_{as}=\left(\begin{array}{cc}1/T_{as}& R_{as}^{}/T_{as}^{}\\ R_{as}/T_{as}& 1/T_{as}^{}\end{array}\right)\left(\begin{array}{cc}\alpha & \beta \\ \beta ^{}& \alpha ^{}\end{array}\right),$$
(1)
where $`R_{as}`$ and $`T_{as}`$ are, respectively, the overall reflection and transmission coefficients for a wave incident from the ambient. Note that for identical ambient and substrate media we have $`det𝖬_{as}=+1`$, which is equivalent to $`|R_{as}|^2+|T_{as}|^2=1`$.

We take as known the theory of reflection from multilayers  and its main result for our purposes, namely that strong reflection will occur when $`[\mathrm{Tr}(𝖬_{as})]^2>4`$ (these conditions, one for the $`s`$ polarization and one for the $`p`$ polarization, locate the band stops for each basic polarization).

When we consider a finite periodic system that consists of $`N`$ basic periods, it is possible to show that in the stop bands the reflectance takes the general form

$$^{(N)}=\frac{|\beta |^2}{|\beta |^2+[\mathrm{sinh}(\chi )/\mathrm{sinh}(N\chi )]^2},$$
(2)
where $`\mathrm{cosh}(\chi )=\mathrm{Re}(\alpha ).`$ We are considering only positive values of $`\mathrm{Re}(\alpha )`$ since negative values can be treated much in the same way.

In practice, it is usual that the basic period of the structure consists of two thin homogeneous dielectric slabs with low, $`n_{\mathrm{}}`$, and high, $`n_h`$, indices of refraction and corresponding thicknesses $`d_{\mathrm{}}`$ and $`d_h`$. These Bragg structures are also appropriately denoted as \[LH\]N, where $`N`$ is the total number of periods. In such a case, the condition $`[\mathrm{Tr}(𝖬_{as})]^2>4`$, required to have a stop band, can be written as

$$|\mathrm{Re}(\alpha )|=|\mathrm{cos}\delta _{\mathrm{}}\mathrm{cos}\delta _h\mathrm{\Lambda }_\mathrm{}h\mathrm{sin}\delta _{\mathrm{}}\mathrm{sin}\delta _h|>1,$$
(3)
where $`\delta _i=(2\pi /\lambda )\mathrm{\Delta }_i`$ is the phase shift of a wave of wavelength in vacuum $`\lambda `$ in traversing the layer $`i`$th and $`\mathrm{\Delta }_i`$ is the corresponding optical path, of value

$$\mathrm{\Delta }_i=n_id_i\mathrm{cos}\theta _i=d_i\sqrt{n_i^2\mathrm{sin}^2\theta _0},$$
(4)
$`\theta _0`$ being the angle of incidence. For simplicity, we have assumed that the system is imbedded in air.

The function $`\mathrm{\Lambda }_\mathrm{}h`$ is

$$\mathrm{\Lambda }_\mathrm{}h=\frac{1+r_\mathrm{}h^2}{1r_\mathrm{}h^2},$$
(5)
where $`r_\mathrm{}h`$ is the Fresnel reflection coefficient for the interface $`\mathrm{}`$-$`h`$. This function $`\mathrm{\Lambda }_\mathrm{}h`$ is frequency independent but takes different forms for $`s`$ and $`p`$ polarizations. However, one can check that, irrespective of the angle of incidence, the following relation for both basic polarizations holds:

$$\frac{\mathrm{\Lambda }_\mathrm{}h(p)}{\mathrm{\Lambda }_\mathrm{}h(s)}=\left(\frac{n_{\mathrm{}}}{n_h}\right)^2<1.$$
(6)
Due to the restriction (6), whenever Eq. (3) is fulfilled for $`p`$ polarization, it is always true also for $`s`$ polarization. In consequence, the $`p`$-polarization stop bands are more stringent than the corresponding $`s`$-polarization ones.

Because of historical reasons , typical use of dielectric mirrors has been evaluated at normal incidence, with layers at a quarter-wavelength thick (at the design frequency):

$$d_{\mathrm{}}=\frac{\lambda _{\mathrm{}}}{4}=\frac{\lambda }{4n_{\mathrm{}}},d_h=\frac{\lambda _h}{4}=\frac{\lambda }{4n_h}.$$
(7)
The optical paths are equal and maximum reflection occurs at the frequency

$$\omega _0=\frac{\pi }{2}\frac{c}{n_{\mathrm{}}d_{\mathrm{}}}=\frac{\pi }{2}\frac{c}{n_hd_h},$$
(8)
which is the center of the stop band.

## 3 Optimization strategy

In order to explain the optimization criterion we wish to propose, in Fig. 1 we have plotted the reflectance $`^{(N)}`$ of a Bragg \[LH\]N structure as a function of the angle of incidence for several values of the number of periods $`N`$ and for $`p`$ polarization. The layer thicknesses are chosen as in Eq. (7).

It is clear from this figure the well-known fact that the reflectance tends to unity as $`N`$ increases. It is then hardly surprising that the usual designs found in the literature use the quarter-wave thicknesses at normal incidence without raising any doubt about such assumption. In practice not very many periods (say $`N10`$) are used in the visible, and one may be tempted to ask whether other thicknesses could improve the performance of the structure. Such a problem could be attacked by a straightforward computation of the reflectance as a function of layer thicknesses. With the a priori information that the optimum condition is close to a quarter wave stack, it will not take too much computational effort to find a reasonable solution . However, our goal here is to provide a more systematic way of dealing with this question.

A reasonable option for optimizing the system is that, when $`\theta _0`$ varies from 0 to $`\pi /2`$, the area under the curve $`1^{(N)}`$ (which is the transmittance of the system) would be as small as possible. Therefore, once the materials and the wavelength are fixed, we propose

$$𝒜(d_{\mathrm{}},d_h)=_0^{\pi /2}[1^{(N)}(d_{\mathrm{}},d_h,\theta _0)]𝑑\theta _0,$$
(9)
as figure of merit for the periodic structure. Alternatively, given the characteristic dip appearing in the reflectance for $`p`$ polarization, one could also impose that this dip would be as smaller as possible. We have numerically checked that both criteria give essentially the same results.

For $`s`$ polarization there is no dip. The different behaviour at oblique incidences of $`s`$ and $`p`$ reflectances has been analyzed previously . Since for a given stack the reflectance $`^{(N)}`$ is greater for $`s`$ polarization than for $`p`$ polarization for every incidence angle (they are equal only at normal and grazing incidences), we argue that once the area $`𝒜(d_{\mathrm{}},d_h)`$ is optimized for $`p`$ polarization, it is also improved for $`s`$ polarization, which seems quite plausible.

We have employed an easy-to-use quasi-Newton algorithm in order to find the minimum of the function (9), subject only to the physical conditions $`0<d_{\mathrm{}}/\lambda _{\mathrm{}}1`$ and $`0<d_h/\lambda _h1`$, because the periodic character of the solutions. The numerical results of this optimization are shown in Table 1. Note that the use of adimensional thicknesses $`d_{\mathrm{}}/\lambda _{\mathrm{}}`$ and $`d_h/\lambda _h`$ simplifies the presentation of the results, although by no means the results are universal: they apply only to this particular system because we are not considering dispersion in refractive indices (in all the paper we take the wavelength in vacuum $`\lambda =10\mu `$m). In any case, the proposed optimization strategy is independent of the particular example under study.

For a better understanding of the performance of our optimized structure, in Fig. 2 we have plotted the reflectance $`^{(N)}`$ (with $`N=4`$ and $`N=7`$) and $`p`$ polarization for both the quarter-wave thicknesses as in Eq.(7) and the optimum thicknesses in Table 1. The improvement is remarkable: in terms of areas, we have 20 % and 100 %, respectively.

To test the tolerance of the optimum thicknesses against small deviations, in Fig. 3 we have a contour plot of $`𝒜(d_{\mathrm{}},d_h)`$ as a function of the adimensional thicknesses $`d_{\mathrm{}}/\lambda _{\mathrm{}}`$ and $`d_h/\lambda _h`$ for the same structure as in Fig. 2 with $`N=7`$ periods. The elliptical contours delimit the range of thicknesses giving a definite value of the area. Moreover, if we take the projection of the major axis of the ellipse on the coordinate axes as a qualitative measure of the maximum tolerance for a given area, we conclude from Fig. 3 that one must be more careful in controlling the thickness $`d_{\mathrm{}}`$ than $`d_h`$.

In the same spirit, one could ask about how critical is the behaviour of the area under deviations in the values of refractive indices. In Fig. 4 we have plotted a contour plot of the area as a function of $`n_{\mathrm{}}`$ and $`n_h`$ for the optimum thicknesses $`d_{\mathrm{}}`$ and $`d_h`$ given in Table 1 for $`N=7`$. We find again the same kind of elliptical contours as in the previous figure, but now the variation is much smoother, indicating that, roughly speaking, the role of refractive indices is not so crucial as the role of thicknesses.

Obviously, as $`N`$ grows the improvement in the area is larger. In fact, table 1 suggests a considerable improvement that, to some extent, may be illusory: it only concerns a decrease in the integrated transmission, while for the reflection, the improvement is not so impressive. However, the important point is that the improvement is just in the dip of the reflectance. We conclude finally that the method is especially appropriate for moderate values of $`N`$ ($`N10`$), which constitute a typical experimental situation.

The results presented so far hold only for the given ratio of refractive indices. To show that the method can be employed for arbitrary parameters of the multilayers, in Fig. 5 we have plotted the optimum adimensional thickness $`(d_{\mathrm{}}/\lambda _{\mathrm{}})_{\mathrm{opt}}`$ as a function of the refractive indices for $`N=7`$, while in Fig. 6 we have represented the optimum values of $`(d_h/\lambda _h)_{\mathrm{opt}}`$. We have assigned a zero thickness whenever the condition $`[\mathrm{Tr}(𝖬_{as})]^24`$, required to have a stop band, is not fulfilled for some values of the incidence angle $`\theta _0`$. The abrupt step is the same in both figures and gives the boundary of omnidirectional reflection for the stack.

As a general feature, we note that $`(d_{\mathrm{}}/\lambda _{\mathrm{}})_{\mathrm{opt}}<(d_h/\lambda _h)_{\mathrm{opt}}`$ for every pair of allowed refractive indices. In Fig. 5, the optimum thickness varies pronouncedly with $`n_{\mathrm{}}`$ until $`n_{\mathrm{}}1.5`$ and then it is almost constant, while decreases slowly with $`n_h`$ in the region of interest. In Fig. 6 we observe that the optimum thickness decreases for a fixed value of $`n_{\mathrm{}}`$ (resp. $`n_h`$) when $`n_h`$ (resp. $`n_{\mathrm{}}`$) increases, and has a flat variation. The same qualitative behaviour is also observed for other values of $`N`$.

In summary, we have proposed a simple figure of merit that allows for an improvement in the performance of omnidirectional reflectors, especially when the number of periods is not too large. The method can be easily extended to other materials and wavelengths and predicts optimum thicknesses that can depart from the usual quarter-wave design.


text:
# On the Real Spectra of Calogero Model with Complex Coupling

## Abstract

We study the eigenvalue problem of the rational Calogero model with the coupling of the inverse-square interaction as a complex number. We show that although this model is manifestly non-invariant under the combined parity and time-reversal symmetry $`𝒫𝒯`$, the eigenstates corresponding to the zero value of the generalized angular momentum have real energies.

preprint: SINP/TNP/03-37

The standard practice in quantum mechanics is to consider self-adjoint operators so that the corresponding spectrum is real and consequently, the time-evolution of the states is unitary. However, it has been found recently that a non-Hermitian Hamiltonian invariant under the combined $`𝒫𝒯`$ symmetry can have real spectrum pt1 ; pt2 ; pt3 ; ali ; ddt ; ex1 ; ex2 ; ex3 ; bm ; zn . The spectrum is entirely real if the $`𝒫𝒯`$ symmetry is unbroken. On the other hand, part of the spectrum is real if the $`𝒫𝒯`$ symmetry is broken spontaneouslypt1 .

In this Letter we study the rational Calogero model calo3 ; pr  with complex coupling for the many-body interaction. Several extensions of the Calogero model with non-Hermitian Hamiltonians have been studied in the literature bm ; zn . However, all these models respect $`𝒫𝒯`$ symmetry and lead to real spectrum. The Hamiltonian we study is identical to that of the original $`N`$-body rational Calogero model except that the coefficient of the inverse square interaction is complex. This model is non-Hermitian, manifestly non-invariant under $`𝒫𝒯`$, but nevertheless is shown to admit a partly real spectrum, with the ground state energy always real.

Following the analysis of Calogero calo3 , the eigenvalue equation of the rational Calogero model with complex coupling for the many-body inverse-square interaction can be reduced to the eigenvalue equation of an effective single particle Hamiltonian containing harmonic and inverse-square interaction. We show that for certain values of the complex coupling of the inverse square interaction, this effective Hamiltonian can indeed be Hermitian when the generalized angular momentum calo3  is zero. This however happens only for the strongly attractive values of the inverse square coupling term in the effective Hamiltonian. Such a strongly attractive single particle system in absence of the harmonic term has already been analyzed by Case case , who obtained the corresponding spectrum in terms of an undetermined parameter which has the physical interpretation of a cutoff in the coordinate space. In our analysis, we explicitly introduce a cutoff which prevents the particles to occupy the same position simultaneously, i.e. it serves as a cutoff in the coordinate space. In the analysis of ref. case , the spectrum of the strongly attractive model in absence of the harmonic term is unbounded from below. In the model under consideration here, which includes the harmonic term, the spectrum is however bounded from below for finite nonzero values of the cutoff.

The Hamiltonian of the rational Calogero model is given by

$$H=\underset{i=1}{\overset{N}{}}\frac{^2}{x_i^2}+\underset{ij}{}\left[\frac{a^2\frac{1}{4}}{(x_ix_j)^2}+\frac{\mathrm{\Omega }^2}{16}(x_ix_j)^2\right]$$
(1)
where $`a`$ and $`\mathrm{\Omega }`$ are constants, $`x_i`$ is the coordinate of the $`i^{\mathrm{th}}`$ particle and units have been chosen such that $`2m\mathrm{}^2=1`$. The constant $`a`$ is conventionally taken to be real, which is required for the operator $`H`$ to be Hermitian. In this Letter we shall however show that for certain complex values of $`a`$, the eigenvalue problem

$$H\mathrm{\Psi }=E\mathrm{\Psi }$$
(2)
still admits a real spectrum with normalizable solutions. In the discussion below, we therefore take $`a=a_R+ia_I`$ where $`a_R`$ and $`a_I`$ are the real and imaginary parts of $`a`$. The Hamiltonian $`H`$ is invariant under parity $`x_ix_i`$, while non-invariant under the time-reversal symmetry $`ii`$. Thus, the combined $`𝒫𝒯`$ symmetry is not respected by $`H`$.

Following calo3 , we consider the above eigenvalue equation in a sector of configuration space corresponding to a definite ordering of particles given by $`x_1x_2\mathrm{}x_N`$. The translation-invariant eigenfunctions of the Hamiltonian $`H`$ can be written as

$$\mathrm{\Psi }=\underset{i<j}{}\left(x_ix_j\right)^{a+\frac{1}{2}}\varphi (r)P_k(x),$$
(3)
where $`x(x_1,x_2,\mathrm{},x_N)`$,

$$r^2=\frac{1}{N}\underset{i<j}{}(x_ix_j)^2$$
(4)
and $`P_k(x)`$ is a translation-invariant as well as homogeneous polynomial of degree $`k(0)`$ which satisfies the equation

$$\left[\underset{i=1}{\overset{N}{}}\frac{^2}{x_i^2}+\underset{ij}{}\frac{2(a+\frac{1}{2})}{(x_ix_j)}\frac{}{x_i}\right]P_k(x)=0.$$
(5)
The existence of complete solutions of (5) for real $`a`$ has been discussed by Calogero calo3 . For both real and complex $`a`$, $`P_0(x)`$ is a constant and is a solution of Eqn. (5). For $`k0`$, the solutions of Eqn. (5) for complex $`a`$ are obtained by analytic continuation in the parameter space $`a`$.

Substituting Eqn. (3) in Eqn. (2), using Eqns. (4\- 5) and making a further substitution $`\varphi =r^{(\mu +\frac{1}{2})}\mathrm{\Phi }`$ we get

$$\stackrel{~}{H}\mathrm{\Phi }=E\mathrm{\Phi },$$
(6)
where

$$\stackrel{~}{H}=\frac{d^2}{dr^2}+\frac{(\mu ^2\frac{1}{4})}{r^2}+\omega ^2r^2,$$
(7)
with $`\omega ^2=\frac{1}{8}\mathrm{\Omega }^2N`$ and

$$\mu =k+\frac{1}{2}(N3)+\frac{1}{2}N(N1)(a_R+ia_I+\frac{1}{2}).$$
(8)
The operator $`\stackrel{~}{H}`$ can be interpreted as the effective Hamiltonian of a particle in the combined harmonic plus inverse-square interaction. Let us now choose the real part of $`a`$ as

$$a_R=\frac{N3}{N(N1)}\frac{1}{2},$$
(9)
and keep $`a_I`$ as an arbitrary real number. It may be noted that the wavefunction $`\mathrm{\Psi }`$ is square-integrable for $`a_R+\frac{1}{2}>\frac{1}{2}`$ us1 ; we , which is indeed satisfied for the above choice of $`a_R`$ for $`N3`$. With the above choice of the complex coupling $`a`$, for $`k=0`$ we see that

$$\mu =\frac{i}{2}N(N1)a_Ii\nu $$
(10)
is purely imaginary and the operator $`\stackrel{~}{H}`$ takes the form

$$\stackrel{~}{H}=\frac{d^2}{dr^2}\frac{(\nu ^2+\frac{1}{4})}{r^2}+\omega ^2r^2,$$
(11)
which is Hermitian and is expected to have real eigenvalues. However, for the same choice of $`a_R`$, when $`k0`$, $`\mu `$ is in general a complex quantity. Thus the eigenvalues of $`\stackrel{~}{H}`$ in the $`k0`$ sector in general would not be real. Below we shall first find the eigenvalues of the operator $`\stackrel{~}{H}`$ for $`k=0`$ and then provide a suitable interpretation of the states in the $`k0`$ sector.

We now proceed to solve the eigenvalue problem given in Eqn. (6) with $`\stackrel{~}{H}`$ given by Eqn. (11). In this case, the inverse-square interaction is necessarily in the strongly attractive regime for $`\nu 0`$. As mentioned before, we put a cutoff $`r_0`$ in the coordinate space which leads to the boundary condition $`\mathrm{\Phi }(r)=0`$ for $`r=r_0`$. We also demand that $`\mathrm{\Phi }(r)\mathrm{L}^2(R^+)`$. Defining $`q=\omega r^2`$ and $`\mathrm{\Phi }=q^{\frac{1}{4}}\chi (q)`$, the eigenvalue equation $`\stackrel{~}{H}\mathrm{\Phi }=E\mathrm{\Phi }`$ can be written as,

$$\frac{d^2\chi }{dq^2}+\left[\frac{1}{4}+\frac{1}{4q^2}(1+\nu ^2)+\frac{E}{4wq}\right]\chi =0,\chi (q_0=wr_0^2)=0.$$
(12)
The above equation can be identified as the Whittaker’s equation, the two linearly independent solutions of which are given by

$$W_{\frac{E}{4\omega },\pm \frac{i\nu }{2}}(q)=e^{\frac{q}{2}}q^{\frac{1\pm i\nu }{2}}M(\frac{1\pm i\nu }{2}\frac{E}{4\omega },1\pm i\nu ,q),$$
(13)
where $`M(\frac{1\pm i\nu }{2}\frac{E}{4\omega },1\pm i\nu ,q)`$ is Kummer’s function abr . The general solution of Eqn. (12) satisfying the boundary condition at $`q_0`$ is given by

$$\chi (q)=A\left[W_{\frac{E}{4\omega },\frac{i\nu }{2}}(q)W_{\frac{E}{4\omega },\frac{i\nu }{2}}(q_0)W_{\frac{E}{4\omega },\frac{i\nu }{2}}(q)W_{\frac{E}{4\omega },\frac{i\nu }{2}}(q_0)\right],$$
(14)
where $`A`$ is the normalization constant. In the limit $`q\mathrm{}`$, we have abr

$$\frac{M(a,b,q)}{\mathrm{\Gamma }(b)}\frac{e^{i\pi a}q^a}{\mathrm{\Gamma }(ba)}+\frac{e^qq^{ab}}{\mathrm{\Gamma }(a)}+O\left(\frac{1}{q}\right).$$
(15)
Thus, as $`q\mathrm{}`$,

$$\chi (q)Ae^{\frac{q}{2}}\left[\frac{\mathrm{\Gamma }(1+i\nu )}{\mathrm{\Gamma }(\frac{1+i\nu }{2}\frac{E}{4\omega })}W_{\frac{E}{4\omega },\frac{i\nu }{2}}(q_0)\frac{\mathrm{\Gamma }(1i\nu )}{\mathrm{\Gamma }(\frac{1i\nu }{2}\frac{E}{4\omega })}W_{\frac{E}{4\omega },\frac{i\nu }{2}}(q_0)\right].$$
(16)
The square integrability of $`\chi `$ can thus be ensured if the quantity in the parenthesis in Eqn. (16) is identically zero, i.e.

$$\frac{W_{\frac{E}{4\omega },\frac{i\nu }{2}}(q_0)}{W_{\frac{E}{4\omega },\frac{i\nu }{2}}(q_0)}=\frac{\mathrm{\Gamma }(1i\nu )}{\mathrm{\Gamma }(1+i\nu )}\frac{\mathrm{\Gamma }(\frac{1+i\nu }{2}\frac{E}{4\omega })}{\mathrm{\Gamma }(\frac{1i\nu }{2}\frac{E}{4\omega })}.$$
(17)
We are interested in the solution of the eigenvalue equation in the limit of a small cutoff, i.e. when $`q_00`$. In this limit, we have $`W_{\frac{E}{4\omega },\frac{i\nu }{2}}(q_0)q_0^{\frac{1i\nu }{2}}`$ abr . The energy eigenvalues are thus determined from the equation,

$$q_0^{i\nu }=\frac{\mathrm{\Gamma }(1i\nu )}{\mathrm{\Gamma }(1+i\nu )}\frac{\mathrm{\Gamma }(\frac{1+i\nu }{2}\frac{E}{4\omega })}{\mathrm{\Gamma }(\frac{1i\nu }{2}\frac{E}{4\omega })},$$
(18)
or equivalently from

$$e^{i(\nu \mathrm{ln}q_0+2\theta )}=e^{2i\alpha },$$
(19)
where $`\theta `$ and $`\alpha `$ are the arguments in the polar representation of $`\mathrm{\Gamma }(1i\nu )`$ and $`\mathrm{\Gamma }(\frac{1+i\nu }{2}\frac{E}{4\omega })`$ respectively. The spectrum is obtained by solving Eqn. (19) graphically. We have studied Eqn. (19) using Mathematica for different values of $`\nu `$ and $`q_0`$. For a fixed $`\nu `$ and $`q_0`$, there is one negative energy bound state and infinitely many positive energy bound states. Moreover, the energy levels are not equispaced. The absolute value of the negative energy depends on the choice of $`q_0`$. For smaller values for $`q_0`$, $`|E|`$ for the negative energy eigenstate increases. We plot our results in Fig.1 and Fig. 2.

The following points about the spectrum in the $`k=0`$ sector may be noted:
1) As shown above, even for certain complex values of the inverse-square coupling, the rational Calogero model admits a real spectrum in the $`k=0`$ sector. The real part of the relevant coupling is $`N`$ dependent and is given by Eqn. (9) whereas the imaginary part of the coupling is arbitrary. Thus there is a range of values of the parameter $`a`$ for which the spectrum in the $`k=0`$ sector is real.

2) The spectrum here for a given value of $`q_0`$ and $`\nu `$ is bounded from below. In order to see this, consider $`E=`$ with $`\mathrm{}`$. In this limit, the argument of the $`\mathrm{\Gamma }(\frac{1+i\nu }{2}\frac{E}{4\omega })`$ can be approximated by $`\frac{\nu }{2}\psi (\frac{1}{2}+\frac{}{4w})`$ and Eqn. (18) can be written as

$$\nu \mathrm{ln}q_02\theta =\nu \psi \left(\frac{1}{2}+\frac{}{4w}\right),$$
(20)
where $`\psi `$ denotes the digamma function abr . For positive values of its argument, the digamma function is monotonically increasing. Thus, for any given finite but nonzero value of $`q_0`$ and $`\nu `$, $`E\mathrm{}`$ is not a part of the spectrum.

3) This system generically admits a single negative energy eigenstate. Existence of negative energy eigenstates in the Calogero model is associated with the self-adjoint extension of the corresponding Hamiltonian us1 ; we . We have however not studied the self-adjointness of the effective Hamiltonian $`\stackrel{~}{H}`$ here.

As mentioned before, the parameter $`\mu `$ is complex when $`k0`$. In that case the effective Hamiltonian would not be Hermitian and the states with $`k0`$ are expected to have complex eigenvalues. These states with complex values of the energy would decay under time evolution. Thus the spectrum obtained here has a part which is stable under time evolution and a decaying part corresponding to complex eigenvalues. This situation is similar to what happens in the case of many nuclei, where only a few energy states are stable and the others are resonances with characteristic decay widths. It would be interesting to find such physical applications of the model described here.

###### Acknowledgements.

The work of PKG is supported(DO No. SR/FTP/PS-06/2001) by the SERC, DST, Govt. of India, under the Fast Track Scheme for Young Scientists:2001-2002.