text:
# Luminous infrared galaxies as possible sources of the UHE cosmic rays

## 1 Introduction

The existence of the UHE cosmic rays, after their discovery almost half a century ago, remains still puzzling. In advent of the new giant experiment named in honour of Pierre Auger, there are many proposals to explain their origin. UHE cosmic rays (UHECR) seem to be extragalactic because their arrival directions are not correlated with the Galactic disk. Thus, sites of their origin should be *different* from normal galaxies (like our Galaxy). In this work we check the hypothesis that the powerful luminous infrared galaxies (LIRGs) might be the UHECR sources. Large fraction of bright IR galaxies are found to be interacting systems, suggesting that collision and merging processes are mostly responsible in triggering the huge IR light emission . Fraction of interacting systems increases with IR luminosity and in the population of the most IR luminous objects in the Universe almost all appear as gravitationally interacting. Favourable environments for accelerating particles to UHE regime via the first order Fermi process are provided by amplified magnetic fields on the scale of tens kpc resulting from gravitational compression, as well as high relative velocities of galaxies and/or superwinds from multiple supernovae explosions . There have been several observational claims that colliding galaxies could be possible sites of the UHECR origin. Al-Dargazelli et al have proposed that some clustering of UHECR shower directions above 10 EeV111$`1EeV=1\times 10^{18}`$ eV are associated with nearby colliding galaxies. Uchihori et al have analysed combined world data from Northern hemisphere experiments and concluded that two triplets and a doublet lie in a vicinity of interacting systems. Takeda et al have noticed that the interacting galaxy VV 141 at z=0.02 is a possible candidate for the triplet of events above 50 EeV from the AGASA experiment. However, as we shall show, there is another favourable candidate for the origin of this UHECR cluster \- Arp 299 (Mrk171, VV 118a/b), a member of LIRG class of extragalactic objects. This system (RA$`=171.4^{}`$, $`\delta =58.8^{}`$), consisting of two interacting starburst galaxies, is the closest extragalactic object (distance 42 Mpc for $`H_0=75kms^1Mpc^1`$) with IR luminosity greater than $`5\times 10^{11}L`$ ($`2\times 10^{45}ergss^1`$) and it is the brightest IR source within 70 Mpc. Such high IR luminosity is related to young and violent star forming regions. There is also observational evidence of superwind outflows, large scale strong radio emissions and the estimated supernova rate is about 0.6 per year .
The aim of the present paper is to check whether the LIRGs could be the sites of origin of the UHECRs observed at the Earth. Using the PSCz catalogue  we construct the all-sky maps of UHE proton intensities originating in LIRGs, taking into account effects of particle propagation through the extragalactic medium and, as an example, possible influence of the regular galactic magnetic field (GMF). We check both hypotheses: origin in LIRGs and, on the other hand, the isotropic distribution of the experimental data above 40 EeV from AGASA , using various statistical tests.

## 2 Anisotropy calculations

### 2.1 The PSCz catalogue

The PSCz catalogue consists of almost 15000 IR galaxies with known redshifts, covering 84% of the sky. It should be noted that observational limitations cause that some of the extragalactic objects may not have measured redshifts or may be even unobserved in the dust obscured regions within the disk of the Galaxy (the so called zone of avoidance). Figure 1 (top) presents the distribution of all PSCz galaxies with known redshifts. We can see patches on the sky regions excluded from the catalogue. On inspection of the superimposed directions of the AGASA showers above 40 EeV we can find only a few cases where experimental events lie in the vicinity of excluded regions, which should not affect much our analysis. As the potential UHECR sources we have selected objects with luminosities in the far infrared (FIR) range exceeding LIRG limit $`L_{FIR}=10^{11}L_{}`$, that is over one order of magnitude greater than the estimated FIR luminosity of our Galaxy , and with distances up to 1 Gpc, giving the total number of 2811 sources (figure 1, bottom). A large fraction of such objects show to be in an apparent stage of collision and merging.

### 2.2 CR propagation

To predict the CR anisotropy expected in the model of LIRG origin we have assumed the following:

* Protons are injected at the sources with a spectrum $`dN/dEE^2`$ truncated at $`10^{21}`$ eV.
* CR luminosity at the source is proportional to its $`L_{FIR}`$ luminosity.
* CR propagate through the intergalactic medium, where they are scattered and suffer energy losses on the cosmic microwave background (CMB).

To calculate CR arrival directions we have derived an analytical formula for the distribution of the deflection angles, for the multiple small angle scattering in weak irregular magnetic fields with continuous energy losses taken into account (see Appendix A). Since strength and structure of the extragalactic magnetic field are largely unknown, we adopt here the upper limit for this field magnitude and coherence scale: $`B(l_c)^{1/2}=1nG(1Mpc)^{1/2}`$, as measured by Faraday rotation of radio signals from distant quasars . Because the analysis is strongly energy dependent, we consider two energy regions 40 to 80 EeV and above 80 EeV, where this limit is just about the Greisen-Zatsepin-Kuzmin (GZK) flux cut-off predicted from CR interactions with the 2.7 K CMB radiation .

### 2.3 Maps of the expected anisotropy

Equatorial maps of the expected intensity of the UHE protons originating in LIRGs for two energy regions, 40-80 EeV and above 80 EeV, with sky coverage and declination dependent exposure for the AGASA experiment  (see also Appendix B), are presented in figures 2 (top) and 3 with superimposed AGASA events.

In figure 2 (top) we can see that expected proton intensities for 40-80 EeV show good correlation with the distribution of the experimental events from AGASA. Especially, clearly visible is the region of the sky where high proton intensities from Arp 299 correlate with the AGASA triplet of events (RA$`170^{}`$, $`\delta 60^{}`$). At higher energies, because of the rapid proton energy losses on CMB, only sources located at distances not larger than 100 Mpc are visible on the map. No correlation of the data events above 80 EeV with the expected CR intensities can be seen (figure 3). There is an apparent group of UHE showers in the region with RA$`300^{}`$, where no possible sources exist. However some of the showers from this group lie close to the region heavily obscured by the Galactic disk.

### 2.4 Influence of Galactic magnetic field

The global regular GMF structure is not well known. Analyses of rotation measures of pulsars and extragalactic radio sources suggest that the GMF has a bisymmetric spiral (BSS) form with field directions reversed from arm to arm . This is also supported by the observed value of the pitch angle of the local field and number of field reversals within the Galactic disk . To examine the influence of the regular GMF on the extragalactic UHE fluxes we have chosen the bisymmetric spiral model with field reversals and odd parity (BSS-A) adopted from Stanev . Because the predictions are strongly model dependent, and the procedure applied gives a rather simplified picture, the analysis presented here is only an example of the possible extragalactic CR flux distortion. It is presented for the lower energy range 40-80 EeV, where the effect is stronger and easily visible.
Simulations of a large number of monoenergetic antiprotons ejected from the Earth and followed to the halo border through GMF give flux modification factors depending on directions in the ’extragalactic’ sky. Then, for an assumed flux distribution of protons at the Galactic borders we are able to calculate their fluxes at the Earth. Distortion effects i.e. reduction or magnification and shifting of the particle fluxes are visible in Figure 2 (bottom). The modification of the extragalactic UHECR intensities due to the influence of the GMF BSS-A model seems to worsen the correlation with the experimental data seen in the upper figure.

## 3 Statistical tests and results

### 3.1 Smirnov-Cramer-von Mises test

To check the hypotheses of isotropic and anisotropic distribution of the experimental data above 40 EeV from AGASA, we have used Smirnov-Cramer-von Mises (SCvM) free of binning test , modified for a 2-dimensional distribution analysis. SCvM test is based on comparing the cumulative distribution function $`F(X)`$ under hypothesis $`H_0`$ with the equivalent function $`S_N(X)`$ of the data. The considered statistics $`W^2`$ is defined as follows

$$W^2=_{\mathrm{}}^{\mathrm{}}\left(S_N(X)F(X)\right)^2f(X)𝑑X$$
(1)
where $`f(X)`$ is the probability density function corresponding to the hypothesis $`H_0`$. $`S_N(X)`$ is based on the experimental data (two coordinates of the UHECR events on the equatorial map) and always increases in steps of equal height, $`N^1`$, where $`N`$ is the total number of data. It is worth to note that this test is reliable even for small statistics ($`N3`$). Critical $`NW^2`$ values for the confidence level $`\alpha =`$0.1 and 0.05 are 0.347 and 0.461 respectively.

Each of the cumulative distribution functions, $`F(X)`$ and $`S_N(X)`$, is the integral of the probability density function over the rectangle area defined by the coordinate point and one of the corners of the map. In this way, from the 2-dimensional maps we have constructed a 1-dimensional $`F(X)`$ and $`S_N(X)`$, for four cases depending on the chosen corner of the map. Although the four obtained $`NW^2`$ values are not completely independent, they provide a valuable insight into the analysis. Results of tested hypotheses: isotropy and LIRGs anisotropy (with and without GMF), energies 40-80 EeV and above 80 EeV are shown in table 1. All the hypotheses in the energy range 40-80 EeV pass, assuming confidence level $`\alpha =`$0.1 ($`NW^2`$=0.347). However, with the presence of the GMF the agreement is worse. The hypotheses for energies above 80 EeV fail, especially in LIRG anisotropy scenario.
Irrespective of the above, table 1 shows that it is not easy to draw conclusions from this test as the $`NW^2`$ values calculated for different map corners scatter quite significantly. Thus, it is advisable to apply another, hopefully more powerful, statistical test

### 3.2 Eigenvector test

The idea is based on assigning unit directional vectors to the data points on the celestial sphere. By finding the normalized eigenvalues $`\tau _1`$, $`\tau _2`$, $`\tau _3`$ of the orientation matrix $`𝐓`$ constructed for the N unit data vectors it is possible to discriminate between the isotropic and some anisotropic distributions . Assuming $`0\tau _1\tau _2\tau _31`$ the empirical shape criterion $`\gamma =`$ $`\left[log_{10}(\tau _3/\tau _2)\right]/\left[log_{10}(\tau _2/\tau _1)\right]`$, and the strength parameter $`\zeta =log_{10}(\tau _3/\tau _1)`$ are used to discriminate girdle type from clustered distribution. Distributions of the girdle and cluster type plot with $`\gamma `$ below and above unity, respectively. Intermediate i.e. partly girdle, partly cluster distributions plot around the line $`\gamma =1`$. Isotropic distributions plot with strength $`\zeta `$ near zero.
In figure 4 (left) there are shown distributions of 500 samples, consisting of 47 events each, simulated under isotropic (dots) and LIRGs anisotropic (crosses) scenarios.

The open square denotes AGASA data, 47 events with energy 40-80 EeV. The two hypotheses, LIRG anisotropy and isotropy, correspond to the two regions in the $`\zeta \gamma `$ plane which overlap significantly for the 47 event samples. Mostly due to the small statistics considered, the eigenvector analysis is not sensitive enough to distinguish between isotropic and LIRG anisotropic origin of the experimental data. Hopefully with a sample of more than 500 events, to be collected in a few months of a full operational mode of the Auger observatory, the two distributions should separate (figure 4, right).
The eigenvector method has been already applied to the UHECR anisotropy analysis by Medina Tanco . The point corresponding to the ’old’ AGASA data (47 events above 40 EeV), calculated by this author, is shown on the $`\zeta \gamma `$ plane (thick cross in figure 4, left) lying well away from the simulated isotropic distribution. However, the point calculated by us for the same AGASA data (star) has a lower strength value and lies on the isotropy distribution. Thus, contrary to the result of Medina Tanco, the AGASA data do not show such a strong clustering (on the basis of the eigenvector analysis) to differ from the isotropy distribution. We have checked the correctness of our result by calculating the eigenvalues of the orientation matrix also analytically (which is possible for a 3x3 matrix).
We have not done the eigenvector analysis for the data above 80 EeV because of the sample of 11 events was too small.

### 3.3 Multiplet analysis

Apart from analysing the large scale LIRG anisotropy (as in previous chapters) we will check this hypothesis by analysing probabilities of multiplets, i.e. groups of events with small angular separation.
In the available222AGASA claims recently  three doublets in the energy range (4-10)$`10^{19}`$eV but it is not clear whether the energy of one shower is below 8$`10^{19}`$ eV  AGASA data from 40-80 EeV there are two doublets (separation angles smaller than 2.5) and one triplet (with criteria taken from the AGASA publication , see also Appendix B). Figure 5 shows probability distribution of the number of doublets based on the $`10^4`$ samples of 47 events each, simulated from the LIRGs anisotropic map (crosses) for the cases of zero and one triplet. There are also similar probabilities for the isotropy scenario (circles), calculated analytically (Appendix B). From figure 5 we find that the probability of obtaining at least the observed number of multiplets from LIRGs anisotropy scenario is $`0.129`$ ($`7.1\times 10^2`$) for more than two (three) doublets, while from isotropy (with the same assumed response in declination of the AGASA experiment) it is only $`0.67\times 10^2(0.2\times 10^2`$). All AGASA events above 40 EeV have one triplet and six doublets, a collection even less probable to be obtained from isotropy.
The above analysis may indicate only that there is some anisotropy (which increases mainly the probability of triplet) but not necessarily that correlated with LIRGs.

A relationship between the set of the data and our hypothesis can be checked by looking for small angle correlation between the simulated and the real events. In figure 6 there are shown smoothed distributions of $`cos\theta `$ in the range $`0.91`$, where $`\theta `$’s are separation angles between the AGASA events, from one side (47 events, 40-80 EeV) and, from another side, 5000 events drawn from the isotropic (dotted line) or from the LIRG (solid line) distributions (giving 47$`\times `$5000 cosine values for each case). The obvious flat distribution of $`cos\theta `$ between the AGASA data and isotropic events (dotted line) with a small rise from left to right shows the effect of the AGASA declination efficiency. The dashed line represents $`cos\theta `$ distribution for $`\theta `$ taken between the events simulated from the LIRG map (47$`\times `$5000 events), averaged over $`10^3`$ such runs. Here, a significant rise in the range of a few degrees ($`\theta <8`$ deg, $`cos\theta >0.99`$) indicates a clustering of the simulated LIRG events. A similar rise is seen in the distribution of $`cos\theta `$ between AGASA data and events simulated from LIRG (solid line). Taking into consideration the sum of all cosines above 0.99, the resulting value lies within 1.67 $`\sigma `$ from the mean value obtained from the LIRG scenario, confirming quite a strong correlation between the real and simulated LIRG events.
In figure 7 there are presented distributions of 47 events from AGASA and equivalent samples of events simulated from the map in figure 2 (top). We have estimated that the probability of getting a doublet or a triplet in the vicinity of Arp 299 is about 0.5. We can see that samples of the 47 events simulated from LIRGs anisotropy scenario, show good resemblance to the experimental data distribution.

## 4 Astrophysical implications

The strongest evidence for UHECRs origin may come from the correlation between directions of the AGASA triplet of events with energies 53.5, 55.0 and 77.6 EeV and an energetic astrophysical source in the local Universe. So far Arp 299 is the best candidate as a member of LIRGs, the brightest infrared source within 70 Mpc and a system of colliding galaxies showing intense, violent starburst activity at only 42 Mpc away. It should be noted that Arp 299 appeared earlier as VV 118 in the list of candidates for the AGASA triplet presented by Takeda et al but was not recognized as a colliding, energetic system, and, as a result, was not given enough attention. These authors point to another object VV 141 being a colliding system at z=0.02. However, as Arp 299 is at a distance two times smaller and fulfils the necessary criteria for CR acceleration, we think that it is this object which could be the most probable CR source.
Recent studies on the nature of LIRGs suggest that compact strong radio emission may result from frequent multiple luminous radio supernovae . These objects are a poorly known class of supernovae with high nonthermal radio power indicating large kinetic energy input to accelerate particles . In systems with intense starburst activity in extremely dense molecular gas and strong magnetic field environment, ”nearly every supernova explosion results in a luminous radio SN with very high radio power” . It is also worth to note that there are some observations suggesting a relationship between gamma ray bursts (GRB) and supernovae explosions . Very recently Weiler et al have stated that GRB980425 and radio loud supernova SN1998bw are possibly related. Thus, it is not unreasonable to invoke here the intriguing hypothesis of the common origin of the two most energetic, mysterious phenomena in the Universe, UHECRs and GRBs.

## 5 Conclusions

We have shown that the available data from AGASA in the energy range 40 to 80 EeV are not in contradiction with the expected anisotropy of CR produced in LIRGs. After applying the GMF, the LIRG hypothesis passes as well, but an apparent worse agreement with the data, suggests that the GMF model used by us may not be appropriate. So, we may hope that, if point sources existed they would be useful for determining the global structure of the GMF. At energies above 80 EeV both the isotropic and the LIRGs anisotropic distribution hypotheses are rejected. This might be explained by an existence of a UHECR population of a different origin above GZK cut-off.
However, it seems that the SCvM test used here is not conclusive for distinguishing between the isotropy and the LIRGs distributions. Contrary to calculations done by other author, the eigenvector analysis has also appeared unable to do this for the existing sample sizes. Thus, we have considered a small scale clustering rather than a large scale anisotropy. Indeed, the analysis of multiplet probability gives good discrimination between isotropic and anisotropic scenarios (figure 5). The probability of the occurring of two doublets and one triplet (three doublets and one triplet) is 10 (20) times higher for the LIRG hypothesis than that for isotropy . The analysis of the distributions of the angular distances between the data and simulated LIRG events also supports an existence of a correlation (figure 6). The strongest argument for our hypothesis is the observation of the triplet from the direction of Arp299, a system of colliding galaxies. We have obtained, via simulations, a high probability of doublets and triplets from this source, estimated to be about 0.5.
Finally, with the prospect of new data to come from the Pierre Auger Observatory in a few years, the arrival direction distribution of UHECR should provide a much better clue to their origin.

The authors thank dr Jan Sroka for stimulating discussions. This work was supported by: KBN (Polish State Committee for Scientific Research) grant No 2PO3C 00618 and University of Lodz grant No 505/447.

## Appendix A. Small angle scattering

Let us consider the statistical process of the small angle scattering of a charged particle on $`N`$ randomly oriented magnetic cells. Here we derive formulae for deviation angles and time delays on simple assumptions that statistical variables $`\delta \stackrel{}{n}_k`$ (denoting change of unit vector along direction of flight occurring in the $`k`$-th cell) are: 1) independent and small, 2) cumulative change of direction is also small, 3) the number of cells $`N`$ is large.

As the particle energy changes along its trajectory, we have the relation

$$\stackrel{}{\vartheta }_k^2=\stackrel{}{\vartheta _N}^2\left(\frac{E_N}{E_k}\right)^2\stackrel{}{\vartheta _N}^2\left(\frac{E_{u=1}}{E_u}\right)^2$$
(2)
Here $`E_N`$ and $`E_k`$ denote particle energies in $`N^{th}`$ and $`k^{th}`$ cells and $`\stackrel{}{\vartheta }_N^2`$ is the mean square scattering angle in the last $`N_{th}`$ cell, which is at the detector. We also denote cells by a continuous variable $`u=k/N`$ instead of $`k`$ (i.e. $`u=1`$ indicates the detector). In this model the particle’s trajectory consists of $`N`$ segments of straight lines of equal length $`L=1/N`$. If the source of particles is placed in the centre of a sphere of radius $`R=1`$, then particle trajectories will end inside the sphere and the time delay is just the additional time the particle needs to reach the sphere. Looking back from the detector we miss the source by the deviation angle $`\eta `$ (figure A1). The formulae for the coordinates of the trajectory end: $`x`$,$`y`$, $`z`$ ($`z`$-axis is the particle initial direction), time delay $`T`$ and angle $`\eta `$ are shown below. The cumulative change of direction is described by the angle $`\mathrm{\Theta }_k`$ as: $`\stackrel{}{n_k}\stackrel{}{n_0}=_{i=1}^k\stackrel{}{\delta n_i}=_{i=1}^k\stackrel{}{\vartheta _i}=\stackrel{}{\mathrm{\Theta }_k}`$.

$`x={\displaystyle \frac{1}{N}}{\displaystyle \underset{k=1}{\overset{N}{}}}\mathrm{\Theta }_{xk};y={\displaystyle \frac{1}{N}}{\displaystyle \underset{k=1}{\overset{N}{}}}\mathrm{\Theta }_{yk};`$ (3)
$`z={\displaystyle \frac{1}{N}}{\displaystyle \underset{k=1}{\overset{N}{}}}(1{\displaystyle \frac{1}{2}}\mathrm{\Theta }_{xk}^2{\displaystyle \frac{1}{2}}\mathrm{\Theta }_{yk}^2);\stackrel{}{\eta }={\displaystyle \underset{k=1}{}}N{\displaystyle \frac{k}{N}}\stackrel{}{\vartheta _k}`$
$`T\mathrm{\Delta }z{\displaystyle \frac{1}{2}}(x^2+y^2)`$ (4)
$`where\mathrm{\Delta }z={\displaystyle \frac{1}{N}}{\displaystyle \underset{k=1}{\overset{N}{}}}({\displaystyle \frac{1}{2}}\mathrm{\Theta }_{xk}^2+{\displaystyle \frac{1}{2}}\mathrm{\Theta }_{yk}^2)`$
We can calculate statistical moments for the above random variables. Here are some of them:

$`\stackrel{}{\eta _{}}^2=N\stackrel{}{\vartheta _N}^2f_{12};\stackrel{}{\eta }^4\stackrel{}{\eta }^2^2=\stackrel{}{\eta }^2^2`$ (5)
$`T={\displaystyle \frac{N}{2}}\stackrel{}{\vartheta _N}^2(f_2f_{11})`$
The constants $`f_{11}`$, $`f_{12}`$ and $`f_2`$ are defined by the following integrals:

$`f_{11}={\displaystyle _0^1}(1u)^2A(u)u;f_{12}={\displaystyle _0^1}u^2A(u)u;`$ (6)
$`f_2={\displaystyle _0^1}(1u)A(u)u`$
where

$$A(u)=\left(\frac{E_{u=1}}{E_u}\right)^2$$

is the particle energy relation from source to detector assuming continuous energy losses. It should be noted that allowing for fluctuations in the energy losses would give larger deflection angles. The values of the derived moments have been positively checked in simulations with various functions $`A(u)`$.

## Appendix B. Probability of multiplets.

Here we derive analytical formulae for average numbers of doublet, triplet and quadruplet events on the assumption that the expected flux is known for any direction and single showers occur independently of each other. The last assumption means that a particular number of multiplets of a given kind (e.g. that of doublets) undergoes Poisson distribution with the expected value equal to the average of the expected (local) values over the entire map. Let us denote by $`\lambda (\delta ,\alpha )`$ ($`\delta `$ \- declination, $`\alpha `$ \- RA) the expected angular density of showers, i.e. the expected number of showers per unit solid angle within the measurement time ($`\lambda =const`$ for isotropy). The actual number of registered showers depends, of course, on the efficiency $`\eta (\delta ,\alpha )`$ of the particular air shower array. Assuming that only showers with zenith angles $`\mathrm{\Theta }`$ smaller than some $`\mathrm{\Theta }_z`$ are considered, we define here that $`\eta (\delta ,\alpha )=f(t)cos(\mathrm{\Theta }(t))`$ where $`f(t)=1`$ while $`\mathrm{\Theta }(\delta ,\alpha ,t)<\mathrm{\Theta }_z`$ and $`f(t)=0`$ otherwise. The brackets mean time average. It can be derived that:

$$\eta (\delta ,\alpha )=\left[\mathrm{sin}(\phi )\mathrm{sin}(\delta )\alpha _z+\mathrm{cos}(\phi )\mathrm{cos}(\delta )\mathrm{sin}(\alpha _z)\right]/\pi $$
(7)
where the angle $`\alpha _z`$ is defined by:

$$\mathrm{cos}(\alpha _z)=\frac{\mathrm{cos}(\mathrm{\Theta }_z)\mathrm{sin}(\phi )\mathrm{sin}(\delta )}{\mathrm{cos}(\phi )\mathrm{cos}(\delta )}$$
(8)
This formula agrees well with the experimental efficiency of the AGASA experiment . The expected angular density $`\rho (\delta ,\alpha )`$ of registered showers equals $`\eta (\delta ,\alpha )\lambda (\delta ,\alpha )`$. If the mean number of showers detected is N, then we have the normalizing relation $`N=\rho (\delta ,\alpha )\mathrm{\Omega }`$ and the average flux $`\lambda `$ is determined:

$$\lambda =\frac{N}{\eta (\delta ,\alpha )\mathrm{\Omega }}$$
(9)
We assume that the number of registered showers in an experiment is large enough to adopt it as the average number of showers which would have been registered by the detector in many such runs. The criterion for classifying two showers as a doublet is that the angle between them should be smaller than some given small value $`R`$. As the expected number of showers within the small angle $`R`$ around a given direction is $`\mu =\rho (\pi R^2)`$ (assuming constant $`\rho `$ within the circle), the local detection density of doublets can be written as $`\rho _2=\frac{1}{2}\rho \mu \mathrm{}^\mu `$ . The factor $`\frac{1}{2}`$ stands for the fact that both members of a doublet count as one doublet. Of course, the local detection density of singles, where there is no companion within the angle $`R`$ from the shower considered is $`\rho _1=\rho \mathrm{}^\mu `$.
The situation is more complex for the case of a triplet. Two criteria may be distinguished:

\- another two showers should be closer than $`R`$ to the shower considered (in the centre) - criterion applied by AGASA.

\- three showers should be within a circle with diameter $`R`$.

The 1st criterion is illustrated in the figure B1 (left). In the centre of the right circle there is the shower denoted as A and within the radius R there are another two showers, out of which the shower denoted as B is in a distance $`x`$ from the first one, in the centre of the left circle. The third shower may be either on the area where the two circles overlap (Ca) or outside it (Cb). In the first case around all three showers a triplet will be found (the same). In the second case the triplet is found only around shower A.
Let S(x) denote the ratio of overlapping area on the figure to the area of the circle. We have the following relations:

$$S(x)=\frac{R^2\left[2\alpha \mathrm{sin}(2\alpha )\right]}{\pi R^2}$$
(10)
where $`\mathrm{cos}(\alpha )=x/(2R)`$. The average value of the ratio $`S(x)`$ equals:

$$S(x)=2_0^RS(x)\frac{x}{R^2}x=1\frac{3\sqrt{3}}{4\pi }$$
(11)
Thus, the local detection density of triplets under 1st criterion can be written as :

$`\rho _3=\rho \left[{\displaystyle \frac{1}{3}}S(x)+\left(1S(x)\right)\right]{\displaystyle \frac{\mu ^2}{2!}}\mathrm{}^\mu `$ (12)
$`={\displaystyle \frac{1}{3}}\rho \left(1+{\displaystyle \frac{3\sqrt{3}}{2\pi }}\right){\displaystyle \frac{\mu ^2}{2!}}\mathrm{}^\mu `$
The mean numbers of doublets $`N_2`$ and triplets $`N_3`$ detected from all directions are found by integrating their detection densities over the whole solid angle.
We get:

$`N_2={\displaystyle \rho _2(\delta ,\alpha )\mathrm{\Omega }}={\displaystyle \frac{1}{2}}(\pi R^2){\displaystyle (\lambda \eta )^2\mathrm{}^\mu \mathrm{\Omega }}{\displaystyle \frac{1}{2}}(\pi R^2)\lambda ^2{\displaystyle \eta ^2\mathrm{\Omega }}`$ (13)
$`N_3={\displaystyle \rho _3(\delta ,\alpha )\mathrm{\Omega }}={\displaystyle \frac{1}{6}}(\pi R^2)^2\left(1+{\displaystyle \frac{3\sqrt{3}}{2\pi }}\right){\displaystyle (\lambda \eta )^3\mathrm{}^\mu \mathrm{\Omega }}`$
$`{\displaystyle \frac{1}{6}}(\pi R^2)^2\lambda ^3\left(1+{\displaystyle \frac{3\sqrt{3}}{2\pi }}\right){\displaystyle \eta ^3\mathrm{\Omega }}`$
where we have assumed that $`\lambda =const`$. Similar consideration as those for detection of triplets under the 1st criterion lead us to the corresponding formulae for quadruplets under the same criterion:

$`\rho _4=\rho \left[{\displaystyle \frac{1}{4}}S_{34}+{\displaystyle \frac{1}{2}}S_{14}+\left(1S_{34}S_{14}\right)\right]{\displaystyle \frac{\mu ^3}{3!}}\mathrm{}^\mu `$ (14)
$`N_4={\displaystyle \rho _4(\delta ,\alpha )\mathrm{\Omega }}={\displaystyle \frac{1}{6}}(\pi R^2)^3\lambda ^40.66{\displaystyle \eta ^4\mathrm{}^\mu \mathrm{\Omega }}`$
where $`S_{34}`$ is the probability that three showers randomly distributed within a circle with a forth shower in the centre will constitute the same quadruplet if *any* shower out of the three is chosen as the central one. The probability that only *one* shower out of the three, chosen as the central one, makes the same quadruplet possible is denoted $`S_{14}`$. It has been found by integration that $`S_{34}=0.274`$ and $`S_{14}=0.270`$.
The formulae for detection densities $`\rho _2,\rho _3,\rho _4`$ derived above allow for such shower configuration, where the shower B in the figure B1 (left) being a member of a doublet or a triplet is at the same time a member of a triplet or a quadruplet, respectively. After excluding such situations we get for doublets the following formula:

$`\rho _2^{^{}}=\rho _2\rho 3\left(1S(x)\right){\displaystyle \frac{\mu ^2}{2!}}\mathrm{}^\mu `$ (15)
Figure B1 (right) illustrates the situation for classifying three showers as a triplet for the 2nd criterion, i.e. three showers have to be within a circle with diameter $`R`$. The two circles there have diameter $`2r=R`$. Two showers from directions A and B are at the angular distance AB $`<R`$. To find the area where the third shower should be we draw the circle of radius $`R/2`$ through the point A and rotate it around point A until point B is inside the circle, so that points E goes to point E’. Similarly, we rotate the circle around point B until point A is inside the circle (point D goes to D’). To form a triplet the third shower must be located within area E’EDD’.
Let $`S_2(x)`$ be the ratio of the framed area to $`\pi R^2`$ for a given value AB$`=x`$. It follows that:

$$S_2(x)=\left[6\alpha +\pi \mathrm{sin}(2\alpha )\right]\frac{1}{4\pi }where\mathrm{cos}(\alpha )=\frac{x}{R}$$
(16)
Taking into account the distribution of distances AB we can calculate the average value of the $`S_2(x)`$:

$$S_2(x)=2_0^RS_2(x)\frac{x}{R^2}x=\frac{9}{16}$$
(17)
Further considerations for the local detection density of triplets under the 2nd criterion are only for the case $`\mu 1`$. So, in a similar way, we get that:

$`\rho _3=\rho {\displaystyle \frac{1}{3}}S_2(x){\displaystyle \frac{\mu ^2}{2!}}={\displaystyle \frac{3}{32}}\rho \mu ^2`$
$`N_3={\displaystyle \rho _3(\delta ,\alpha )\mathrm{\Omega }}={\displaystyle \frac{1}{6}}(\pi R^2)^2{\displaystyle \frac{9}{16}}{\displaystyle (\lambda \eta )^3\mathrm{\Omega }}={\displaystyle \frac{3}{32}}(\pi R^2)^2\lambda ^3{\displaystyle \eta ^3\mathrm{\Omega }}`$
Numbers of observed doublets, triplets, quadruplets etc undergo Poisson distributions with the corresponding mean (expected) values $`N_2,N_3,(N_2^{^{}},N_3^{^{}})`$, etc as calculated above. The formulae for the mean numbers of the multiplets have been positively checked in simulations. For the AGASA experiment the following values apply:
$`\mathrm{\Theta }_z=45^{};\phi =35^{}47^{};R=2.5^{}`$
We find that :
$`\eta \mathrm{\Omega }=\pi (1\mathrm{cos}^2\mathrm{\Theta }_z)=1.57sr;\eta ^2\mathrm{\Omega }=0.38sr;\eta ^3\mathrm{\Omega }=0.097sr;\eta ^4\mathrm{\Omega }=0.0255sr`$
For the expected value of showers N we use the actual number (registered in the energy range 40 to 80 EeV) - 47 showers and obtained that:
$`N_2=1.018;N_3=0.0283;N_2^{^{}}=0.954;N_3^{^{}}=0.0268`$ under 1st criterion; $`N_3=0.0087`$ under 2nd criterion; $`N_4=4.57\times 10^4`$ for quadruplets under 1st criterion.
Thus the probability of obtaining two doublets or more and one triplet or more, for the isotropic sky without point sources equals $`\left[1exp(N_2^{^{}})N_2^{^{}}exp(N_2^{^{}})\right]\left[1exp(N_3)\right]=0.2470.0279=6.9\times 10^3`$. We see that it is mainly the triplet event that makes the isotropy hypothesis very unlikely.

## References


format: text

text:
# 1 Exact resonance widths of a circle (radius 𝑅, refractive index 𝑛=3.29) at 35<𝑘𝑅<35.5 versus angular momentum 𝑚. Dashed line indicates 𝜒_𝑐 in Fresnel’s law, using Eq. (). Also shown are the TM widths. The Brewster angle is at 𝑚≈34 and causes a peak in the TE widths.

Martina Hentschel and Jens U. Nöckel111Current address: Department of Physics, University of Oregon, 1371 E 13th Avenue, Eugene, OR 97403

Published in Quantum Optics of Small Structures, edited by D.Lenstra, T.D.Visser and K.A.H.van Leeuwen (Edita KNAW, Amsterdam, 2000)

The sequential-reflection model in
deformed dielectric cavities

The stationary states of a microlaser are related to the decaying quasibound states of the corresponding passive cavity. These are interpreted classically as originating from sequential escape attempts of an ensemble of rays obeying a curvature-corrected Fresnel formula. Polarization-dependent predictions of this model, and its limitations for stable orbits in partially chaotic systems are discussed.

As a mechanism for achieving mode confinement, waveguiding by total internal reflection is ubiquitous in optics. However, in dielectric microresonators where three-dimensionally confined mode volumes are desired, there is always leakage because the ray picture, in which Fresnel’s formulas describe the outcoupling, acquires corrections. Leaky modes corresponding to classically confined rays can be found, e.g., in optical fibers as “spiral” modes (Poon), or in latterally structured cylindrical VCSELs (Ahn) as well as in microdisk lasers (McCall). The classically forbidden loss in such modes is analogous to tunneling through an effective potential barrier (Johnson).

The highest Q is achieved for modes which semiclassically correspond to rays almost at grazing incidence. Resonators with a circular cross section are a particularly simple realization of this requirement, because they exhibit whispering-gallery (WG) modes characterized by high intensity in an annular region near the surface. However, even Lord Rayleigh who first described the acoustic analog that gave the phenomenon its name, concluded (Strutt) that it requires only an everywhere positive curvature, not necessarily rotational symmetry. A rigorous proof of this is difficult because in the short-wavelength limit, this “clinging” of waves to the walls has to carry over to the ray picture, in which a generic oval cavity exhibits a transition to chaos (Lazutkin; Nöckel, 1996a). Notwithstanding, this problem is fundamental to microresonator design (Nöckel, 1997a), because the availability of high-Q modes is the foremost selection criterium in an otherwise unbounded space of potential resonator shapes (Angelow). Chaos can in fact make WG modes more useful, and moreover create other types of modes with desirable properties, such as the bowtie pattern whose confocal geometry points the way toward the strong-coupling regime in combination with focused emission (Gmachl, Morin).

The robustness of whispering-gallery type intensity patterns in the modes of convex resonators extends even to nonlinear media (Harayama). However, in that case the distinction to the widely studied phenomenon of vortex formation (Weiss) becomes washed out: a WG mode is also a vortex with a phase singularity at points of vanishing intensity; for a circular resonator where the field is proportional to a Bessel function $`J_m(kr)r^m`$ near the center $`r=0`$, the vorticity is simply the angular momentum quantum number $`m`$. Therefore, since our aim is to address the fundamental aspect of the shape dependence of high-Q modes in microresonators, we focus here on linear media where amplification is taken into account by a negative imaginary part of the refractive index $`\stackrel{~}{n}`$.

The model considered here can be derived from a homogenous cylinder by deforming its cross section and considering only propagation transverse to its axis. In this case, TE and TM polarization are decoupled and one has to consider only a scalar wave equation

$$^2\psi +\stackrel{~}{n}^2k^2\psi =0,$$
(1)
assuming a steady state time dependence so that $`k`$ is real. Here, $`\stackrel{~}{n}nin^{}`$ inside the resonator and $`\stackrel{~}{n}=1`$ outside, giving rise to an exterior and interior field, $`\psi _{\mathrm{ext}}`$ and $`\psi _{\mathrm{int}}`$, both of which are connected by the proper matching conditions at the dielectric interface, depending on polarization. For TM modes, $`\psi `$ denotes the electric field, which is parallel to the cylinder axis. In this case one finds that $`\psi `$ and its normal derivative are continuous at the interface, in analogy to quantum mechanics.

The system is open because it radiates energy into the environment via its modal losses. This openness increases as $`n1`$, and the closed-resonator limit is approached for $`n\mathrm{}`$. This can be understood from Fresnel’s formulas which imply total internal reflection for all angles of incidence $`\chi `$ satisfying $`\mathrm{sin}\chi >1/n\mathrm{sin}\chi _c`$ ($`\chi _c`$ is the critical angle). Equation (1) can be recast as

$$^2\psi +n^2\stackrel{~}{k}^2\psi =0,$$
(2)
where $`n`$ is the real part of $`\stackrel{~}{n}`$ as defined above, and $`\stackrel{~}{k}kikn^{}/n`$ is a complex wavenumber inside the cavity but reduces to $`\stackrel{~}{k}=k`$ outside. If instead of this we also had $`\stackrel{~}{k}=kikn^{}/n`$ outside, the solutions of Eq. (2) would be the quasibound states of the passive resonator, as they arise when one assumes a decaying time dependence $`\mathrm{exp}[ickt\gamma t]`$, where $`\gamma =ckn^{}/n`$.

For a quasibound (or metastable) state, the field at distances larger than $`c/(2\gamma )`$ from the cavity grows exponentially due to retardation, but within this physical range $`\psi _{\mathrm{ext}}`$ vanishes as $`\gamma 0`$, so that one can write $`\psi _{\mathrm{ext}}(𝐫)\gamma \zeta (𝐫)`$. If one expands the dependence of $`\psi _{\mathrm{int}}`$ and $`\zeta `$ on $`\gamma `$ in a Taylor series, then to linear order the $`\gamma `$-dependence of $`\zeta `$, but not that of $`\psi _{\mathrm{int}}`$, can be dropped in the full solution. Therefore, the stationary state of the active medium and the metastable decaying state are identical to first order in $`\gamma `$ within an area of order $`\gamma ^2`$.

This approximate equivalence establishes a connection to the study of S-matrix poles from which quasibound states arise, see (Schomerus), and to dissipation in quantum mechanics (Prigogine, Braun). The recent resurgence of interest in these problems is motivated to a significant extent by our lack of understanding of the quantum-to-classical transition, in particular in the presence of classical chaos. Precisely this constellation is also present in Eq. (2) when one considers its short-wavelength limit for the generic case of a deformed cavity.

In the context of laser resonators, there are three main differences to previous work on open quantum systems: firstly, we are interested in the properties of individual states of an open system, as opposed to a statistical ensemble; see also (Casati). Secondly, an important quantity that can be studied for such individual states is their emission directionality, which in other open systems of chemical or nuclear physics is averaged out. Finally, the classical limits of quantum mechanics with smooth potentials and optics with discontinuous refractive indices are qualitatively different (Kohler): the first yields deterministic Hamiltonian mechanics; the second leads to the probabilistic Fresnel formulas which moreover depend on polarization.

In principle, Eq. (2) can be solved numerically to find the discrete complex $`\stackrel{~}{k}`$ and the corresponding modes. One approach is based on the Rayleigh hypothesis (van den Berg, Barton) which in our implementation for quasibound states (Nöckel, 1996a) assumes that the fields can be expanded in cylinder functions as

$`\psi _{\mathrm{int}}(r,\varphi )`$ $`=`$ $`{\displaystyle \underset{m}{}}A_mJ_m(kr)e^{im\varphi },`$ (3)
$`\psi _{\mathrm{ext}}(r,\varphi )`$ $`=`$ $`{\displaystyle \underset{m}{}}B_mH_m^{(1)}(kr)e^{im\varphi }.`$
where a polar coordinate system with suitably chosen origin is used. These expansions always work inside some circle of convergence for $`\psi _{\mathrm{int}}`$ and outside some other circle for $`\psi _{\mathrm{ext}}`$, and for a large range of resonator shapes both convergence domains contain the dielectric interface where the matching conditions are imposed to obtain equations for the unknown coefficients $`A_m`$ and $`B_m`$.

Computational cost can be high here, especially at short wavelengths, and hence a semiclassical approximation can lead to simplifications while preserving physical insight. The ray picture is a cornerstone of classical optics, but its value in the study of open resonators only unfolds when the ray dynamics is studied in phase space (Nöckel, 1994, 1996a, 1996b; Mekis), because Fresnel’s formulas determine escape probabilities according to the angle of incidence $`\chi `$, not the position of impact. One can make use of the physical information contained in this picture in two ways: Either one starts from Eq. (2) and takes a short-wavelength limit (Narimanov); or alternatively, one starts from the classical dynamics and makes classical approximations that allow one to impose simple quantization conditions and thus make the connection to the resonator modes (Nöckel, 1997a). The question whether these different routes meet “in the middle” is not straightforward because the problem of semiclassical quantization in a generic deformed resonator is not completely solved as yet, owing to the coexistence of both regular and chaotic motion in their classical phase space.

Among the advantages of the ray-based approach (Nöckel, 1997a) are its flexibility and computational ease. However, in order for the prescription outlined in Ref. (Nöckel, 1997a) to correctly describe the limiting case of a circular cylinder, one must include the tunneling which in the circle is the only loss mechanism. This can be done in the ray picture with a curvature- and wavelength dependent “rounding” of Fresnel’s formulas which the simulation uses at each reflection along a ray path. The idea used in (Nöckel, 1997a) was to interpret the resonance widths of a circular cylinder in terms of a “sequential-tunneling” ansatz: if the intensity of a quasibound state decays as $`\mathrm{exp}[2\gamma t]`$, this can be interpreted in the ray picture as the result of $`\nu `$ sequential escape attempts with reflection probability $`p_0`$, where $`\nu `$ is the number of reflections the ray undergoes during the time $`t`$. In a circle of radius $`R`$, a trajectory characterized by the angle of incidence $`\chi `$ has $`\nu =ct/(nL)`$ reflections during $`t`$ ($`c/n`$ is the speed of light in the passive medium and $`L=2R\mathrm{cos}\chi `$ is the geometric path length between reflections). Therefore, one expects a decay law $`p_0^\nu =\mathrm{exp}[ct\mathrm{ln}p_0/(nL)]`$. Comparison with the wave result yields

$$p_0=\mathrm{exp}(2nL\gamma /c)$$
(4)
An analytic approximation for $`\gamma `$ in the circle with TM polarization has been derived in (Nöckel, 1997b),

$`\gamma `$ $``$ $`{\displaystyle \frac{c}{2nR}}\mathrm{ln}\left[{\displaystyle \frac{n1}{n+1}}\right]\times `$
$`{\displaystyle \frac{2}{\pi kR[J_m^2(kR)+Y_m^2(kR)]}}.`$
Using the semiclassical expression

$$m=nkR\mathrm{sin}\chi $$
(6)
for the angular momentum (Nöckel, 1996a), one then obtains the reflecticity in terms of purely classical variables, $`p_0(kR,\mathrm{sin}\chi )`$. It reduces to Fresnel’s formula in the limit of large radius of curvature $`R`$, and by construction reproduces the width of a mode in the circle if applied locally at each reflection in our classical ray model. The latter does not hold uniformly for a similar correction derived in (Snyder).

From the classical limit, it follows that there exists an upper bound on resonance widths for dielectric cavities with stepped index profiles and “quantum-mechanical” continuity conditions on $`\psi `$, because the reflectivity $`p_0`$, (for polarization perpendicular to the plane of incidence), is bounded away from zero. This minimum $`p_{0,min}`$ will limit the width of resonances in a cavity of characteristic size $`l`$ to $`\gamma _{max}=c\mathrm{ln}p_{0,min}/(2nl)`$. Smooth index profiles can also appear discontinuous on the scale of the wavelength but are eventually resolved as $`k\mathrm{}`$, allowing arbitrarily small reflectivities at perpendicular incidence.

However, extending these arguments to TE polarization where the electric field is in the plane of incidence, we furthermore conclude that a similar upper bound on the widths does not exist even for sharp interfaces. The reason is that Fresnel’s formula yields zero reflectivity at the Brewster angle $`\chi _B`$ at $`\mathrm{sin}\chi _B=(1+n^2)^{1/2}`$. The normal derivative of $`\psi `$ (which now represents the magnetic field) exhibits a jump proportional to $`n^2`$ at the dielectric interface – reminding us that this is a situation unique to optics. These general considerations have important implications for microresonator design especially at the large $`n`$ typical for semicondutors, because in that case $`\mathrm{sin}\chi _B1/n`$, i.e., the “hole” in the reflectivity for TE polarization approaches $`\chi _c`$ for total internal reflection from below. Taking tunneling due to finite curvature into account as in Eq. (4), the rounded Fresnel formula then exhibits reduced reflectivity even for incidence somewhat above $`\chi _c`$. One can approximately obtain the TE widths of the circle from $`\gamma =c\mathrm{Im}[\sigma ]/(nR)`$,

$$\sigma \mathrm{arctan}[\left(n\frac{H_{m1}^{(1)}(kR)}{H_m^{(1)}(kR)}\frac{m}{kR}(n1/n)\right)^1]$$
(7)
This is illustrated in Fig. 1, where refractive index and $`kR`$ are chosen close to those of Ref. (Gmachl). The reason is that the quantum-cascade material used there emits preferentially TM polarization, whereas the pioneering MQW microdisk lasers with sub-micron thickness permit guiding in the vertical direction only for TE modes \[cf. McCall (1992); there, TE/TM must be interchanged to get from the slab-waveguide to our cylinder convention\]. It is thus important to ascertain whether the identical oval lateral design of the quantum-cascade lasers in Ref. (Gmachl) would also permit a microdisk laser to operate in TE polarization.

The lasing mode in Fig. 3D of (Gmachl) was identified as a bowtie-shaped pattern corresponding to a periodic ray path with angle of incidence given by $`\mathrm{sin}\chi 1/n`$, i.e., directly at the critical angle. That this mode provides high Q can be seen by comparing to $`\gamma `$ in Fig. 1: assuming that the width $`\gamma _B`$ of a bowtie mode results from the sequential application of $`p_0`$ as determined for the circle, the argument leading to Eq. (4) implies that $`\gamma _B\gamma L/l`$ where $`L/l1.13`$ is the ratio of the classical path lengths between reflections in the WG orbit and bowtie, respectively. One sees that the TM line intersects the critical angle (corresponding to $`m=kR`$) at a much smaller width than the TE curve, and this Q-spoiling due to Brewster transmission is borne out by the actual TE resonances as well. This leads to the prediction that conventional microdisk lasers with a shape designed to yield a bowtie pattern just at $`\chi _c`$ as in Fig. 3D of (Gmachl) will not lase.

These ray arguments are known to yield large deviations from the true resonance widths when the modes under consideration are quantized on stable phase-space domains in a partially chaotic system, cf. Ref. (Nöckel, 1997a) where this was attributed to chaos-assisted tunneling. The latter yields enhanced outcoupling and hence the true widths are underestimated by the sequential ray picture. Therefore, the above Q-spoiling for TE modes is not counteracted by a correction of this nature. The prediction of an upper bound for TM widths is also not affected by chaos-assisted tunneling because it cannot be faster than the fastest classical process, which in turn is limited by $`p_0`$ at $`\mathrm{sin}\chi =0`$. Beyond this, however, quantitative widths for stable-orbit modes in mixed phase spaces are not provided by the ray model.

The disagreement is illustrated in Fig. 2 for a bowtie mode similar to the ones studied in (Narimanov), as a function of $`n`$, but at a deformation of $`ϵ=0.16`$ \[defined as in (Gmachl)\] and $`nkR119.8`$. Since $`nk`$ is the wavenumber inside the resonator, it should remain approximately independent of $`n`$ as long as the outcoupling can be taken into account in the form of a boundary phase shift intermediate between Dirichlet and Neumann. Indeed, for the state shown in Fig. 2, the change in $`nkR`$ in the plotted range of $`n`$ is only $`0.2`$. The length scale $`R`$ here is the radius of curvature at the points of reflection.

At small $`n`$ where $`\chi _c`$ is larger than the angle of incidence of the bowtie, escape is classically allowed in Fresnel’s formula and hence curvature corrections are unimportant. At $`n>3`$, the tunneling correction in Eq. (S0.Ex2) does improve on the classical Fresnel prediction ($`\gamma =0`$) but clearly still underestimates the true width. As tunneling in general is definable only with respect to a classical expectation, we could again label the discrepancy as chaos-assisted tunneling. However, a semiclassical theory starting from Eq. (2) which reproduces the exact behavior in Fig. 2 very well (Narimanov) can shed more light on the physics of the phenomenon.

References:

* et al., Phys. Rev. Lett. 82, 536
* , F. Laeri, T. Tschudi, Opt. Lett. 21, 1324
* , Appl. Opt. 36, 1312
* , P. L. Braun and F. Haake, Physica D 131, 265
* , G. Maspero and D. L. Shepelyanski, Phys. Rev. E 56, R6233
* , G. Maspero and D. L. Shepelyanski, Physica D 131, 311
* et al., Science 280, 1556
* , P. Davis and K. S. Ikeda, Phys. Rev. Lett. 82, 3803
* , J. Opt. Soc. Am. 10, 343
* , R. Blümel, Ann. Phys. 267, 249
* , KAM Theory and Semiclassical Approximations to Eigenfunctions, (Springer, New York)
* , A. F. Levi, R. E. Slusher, S. J. Pearton and R. A. Logan, Appl. Phys. Lett. 60, 289
* , J. U. Nöckel, G. Chen, A. D. Stone and R. K. Chang, Phys. Rev. Lett. 75, 2682
* , C. C. Yu and T. W. Mossberg, Phys. Rev. Lett. 73, 1489
* , G. Hackenbroich, P. Jaquod and A. Douglas Stone, cond-mat/9907109
* , A. D. Stone and R. K. Chang, Optics Letters 19, 1693
* , and A. D. Stone, in: Optical Processes in Microcavities, edited by R. K. Chang and A. J. Campillo (World Scientific, Singapore)
* , A.D.Stone, G.Chen, H. Grossman and R. K. Chang, Opt. Lett. 21, 1609
* , A. D. Stone, Nature 385, 45
* , Dissertation, Yale University
* , H. Schomerus, C. W. J. Beenakker, cond-mat/9905019
* , R. K. Chang and J. A. Lock, Opt. Lett. 23, 1105
* , Phys. Rep. 219, 93
* , J.D.Love, IEEE Trans. MTT-23, 134
* , and Baron Rayleigh, The Theory of Sound, Vol.II (Dover, New York)
* et al., Appl. Phys. B 68, 151 (1999)
* , J.T.Fokkema, IEEE Trans. Anten. Propag. AP-27, 577

Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Str. 38, 01187 Dresden, Germany


text:
# References

hep-ph/0203237

Minimal Schemes for Large Neutrino Mixings with Inverted Hierarchy

Hong-Jian Hea Electronic address: hjhe@physics.utexas.edu, Duane A. Dicusa Electronic address: phbd057@utxvms.cc.utexas.edu, John N. Ngb§§§ Electronic address: misery@triumf.ca

a Center for Particle Physics, University of Texas at Austin, Austin, Texas 78712, USA

b TRIUMF, 4004 Wesbrook Mall, Vancouver, BC V6T 2A3, Canada

## Abstract

Existing oscillation data point to nonzero neutrino masses with large mixings. We analyze the generic features of the neutrino Majorana mass matrix with inverted hierarchy and construct realistic minimal schemes for the neutrino mass matrix that can explain the large (but not maximal) $`\nu _e^{}\nu _\mu ^{}`$ mixing of MSW-LAM as well as the nearly maximal $`\nu _\mu ^{}\nu _\tau ^{}`$ mixing and the small (or negligible) $`\nu _e^{}\nu _\tau ^{}`$ transition. These minimal schemes are quite unique and turn out to be extremely predictive. Implications for neutrinoless double beta decay, tritium beta decay and cosmology are analyzed.
PACS numbers: 14.60.Pq, 12.15.Ff, 13.15.$`+`$g, 13.40.Em \[March, 2002\]

1. Introduction

The large, rather than small, neutrino mixings confirmed by atmospheric and solar oscillation experiments  over the recent years have brought neutrino physics to an exciting new era. It indicates that lepton flavor mixing is very different from quark flavor mixing, and neutrino mass generation may have a distinct origin from the traditional Dirac-type Yukawa interactions for the charged quarks and leptons in the standard model (SM). In fact, the neutrino masses can be naturally of Majorana nature, generated from either a seesaw mechanism  at high scales or a radiative mechanism around the weak scale .

The current global fit strongly favors Mikheyev-Smirnov-Wolfenstein Large Angle Mixing (MSW-LAM) in which the solar mixing angle $`\theta _{}`$ (giving the $`\nu _e^{}\nu _\mu ^{}`$ transition) is large but significantly deviates from the maximal value $`45^{}`$, i.e., $`25^{}\theta _{}39^{}`$ at 95% C.L.111The maximal value $`45^{}`$ is also excluded by the 99% C.L. limit of the MSW-LAM solution, $`24^{}\theta _{}43^{}`$ ., with a central value at $`\theta _{}32^{}`$ . On the other hand, the atmospheric data indicate a maximal mixing angle $`\theta _{\mathrm{atm}}`$ (representing the $`\nu _\mu ^{}\nu _\tau ^{}`$ transition), with the 95% C.L. limit $`33^{}\theta _{\mathrm{atm}}57^{}`$ and the central value $`\theta _{\mathrm{atm}}45^{}`$ . This is also supported by the K2K long baseline experiment . The Chooz and Palo Verde  long baseline reactor experiments (in combination with the mass range of atmospheric data ) bound $`\mathrm{sin}^2\theta _{\mathrm{chz}}0.04`$ at 95% C.L., where the angle $`\theta _{\mathrm{chz}}`$ measures the $`\nu _e^{}\nu _\tau ^{}`$ transition. Furthermore, the solar oscillations constrain the mass-square difference $`\mathrm{\Delta }_{}=|m_1^2m_2^2|`$ to be, $`1.8\times 10^5\mathrm{eV}^2\mathrm{\Delta }_{}4.1\times 10^4\mathrm{eV}^2`$, for MSW-LAM at 99% C.L., while the atmospheric oscillations confine the mass-square difference $`\mathrm{\Delta }_{\mathrm{atm}}=|m_{1,2}^2m_3^2|`$ as, $`1.3\times 10^3\mathrm{eV}^2\mathrm{\Delta }_{\mathrm{atm}}5.5\times 10^3\mathrm{eV}^2`$, at 99% C.L. This suggests two generic patterns for the light neutrino mass-eigenvalues, $`(m_1^{},m_2^{},m_3^{})0`$, namely, the “Normal Hierarchy” (called Type-A) and “Inverted Hierarchy” (called Type-B),

$$\mathrm{A}:m_1^{}<m_2^{}m_3^{};\mathrm{B}:m_1^{}m_2^{}m_3^{}.$$
(1)
Hereafter, we will focus on the Type-B scenario with inverted mass hierarchy. Our present goal is to construct a realistic scheme for the neutrino Majorana mass matrix, containing only a minimal set of parameters to describe the neutrino data, especially the non-maximal solar neutrino mixing à la MSW-LAM (which is hard  to realize in models with an approximate $`L_eL_\mu L_\tau `$ symmetry ). We show that such a Minimal Scheme can be quite uniquely derived and is highly predictive. Implications for neutrinoless double $`\beta `$ decay, tritium $`\beta `$ decay, and cosmology are analyzed.

2. Minimal Schemes for Neutrino Majorana Masses with Inverted Hierarchy

Consider the generic $`3\times 3`$ symmetric Majorana mass matrix $`M_\nu `$ for 3 light flavor-neutrinos $`(\nu _e,\nu _\mu ,\nu _\tau )`$, at the weak scale and with leptons in the mass-eigenbasis,

$$M_\nu =\left(\begin{array}{ccc}m_{ee}& m_{e\mu }& m_{e\tau }\\ m_{e\mu }& m_{\mu \mu }& m_{\mu \tau }\\ m_{e\tau }& m_{\mu \tau }& m_{\tau \tau }\end{array}\right).$$
(2)
With extra new heavy fields integrated out, Eq. (2) is the most general description of the Majorana masses of three active neutrinos based upon Weinberg’s unique dimension-5 effective operator , $`{\displaystyle \frac{𝒞_{ij}}{\mathrm{\Lambda }}}L_i^\alpha L_j^\beta H^\alpha ^{}H^\beta ^{}ϵ^{\alpha \alpha ^{}}ϵ^{\beta \beta ^{}}`$, which gives a mass term, $`\frac{1}{2}\nu ^TM_\nu \nu `$, with $`M_\nu ^{ij}=𝒞_{ij}v^2/\mathrm{\Lambda }`$, where $`H=v/\sqrt{2}`$ is the vacuum expectation value of the SM Higgs doublet. The mass matrix (2) contains nine independent real parameters, which can be equivalently chosen as three mass eigenvalues $`(m_1^{},m_2^{},m_3^{})0`$, three mixing angles $`(\theta _{12},\theta _{23},\theta _{13})`$, and three CP-violation phases $`(\varphi ,\varphi ^{},\varphi ^{\prime \prime })`$ with $`\varphi `$ the usual Dirac phase and $`(\varphi ^{},\varphi ^{\prime \prime })`$ the Majorana phases (which do not affect the neutrino oscillation). The neutrino mixing matrix $`VUU^{}`$ for diagonalizing $`M_\nu `$, via $`V^TM_\nu V=M_\nu ^{\mathrm{diag}}`$, contains six parameters (three rotation angles and three phases) and can be decomposed into a matrix $`U`$ (à la Cabibbo-Kobayashi-Maskawa) and a diagonal matrix $`U^{}`$ with only two Majorana phases,

$$U=\left(\begin{array}{ccc}c_1c_3& s_1c_3& s_3e^{i\varphi }\\ s_1c_2c_1s_2s_3e^{i\varphi }& c_1c_2+s_1s_2s_3e^{i\varphi }& s_2c_3\\ s_1s_2+c_1c_2s_3e^{i\varphi }& c_1s_2s_1c_2s_3e^{i\varphi }& c_2c_3\end{array}\right)$$
(3)
and $`U^{}=\mathrm{diag}(1,e^{i\varphi ^{}},e^{i\varphi ^{\prime \prime }})`$ . Here we use the notations $`(\theta _1,\theta _2,\theta _3)(\theta _{12},\theta _{23},\theta _{13})`$, for convenience. From the mass diagonalization, we can reconstruct the neutrino mass matrix $`M_\nu `$ via the relation,

$$M_\nu =V^{}M_\nu ^{\mathrm{diag}}V^{},$$
(4)
which gives,

$$\begin{array}{ccc}m_{ee}\hfill & =& c_3^2\left[c_1^2m_1+s_1^2m_2^{}\right]+p^2s_3^2m_3^{},\hfill \\ m_{\mu \mu }\hfill & =& (s_1c_2\overline{p}c_1s_2s_3)^2m_1+(c_1c_2+\overline{p}s_1s_2s_3)^2m_2^{}+s_2^2c_3^2m_3^{},\hfill \\ m_{\tau \tau }\hfill & =& (s_1s_2+\overline{p}c_1c_2s_3)^2m_1+(c_1s_2\overline{p}s_1c_2s_3)^2m_2^{}+c_2^2c_3^2m_3^{},\hfill \\ m_{e\mu }\hfill & =& c_3\left[s_1c_1c_2(m_1m_2^{})\overline{p}s_2s_3(c_1^2m_1+s_1^2m_2^{})+ps_2s_3m_3^{}\right],\hfill \\ m_{e\tau }\hfill & =& c_3\left[s_1c_1s_2(m_1m_2^{})+\overline{p}c_2s_3(c_1^2m_1+s_1^2m_2^{})pc_2s_3m_3^{}\right],\hfill \\ m_{\mu \tau }\hfill & =& (s_1s_2+\overline{p}c_1c_2s_3)(s_1c_2\overline{p}c_1s_2s_3)m_1+(c_1s_2\overline{p}s_1c_2s_3)(c_1c_2+\overline{p}s_1s_2s_3)m_2^{}s_2c_2c_3^2m_3^{},\hfill \end{array}$$
(5)
where $`(m_2^{},m_3^{})(m_2^{}e^{i2\varphi ^{}},m_3^{}e^{i2\varphi ^{\prime \prime }})`$ and $`p=\overline{p}^{}e^{i\varphi }`$.

The precise form of the neutrino mass matrix $`M_\nu `$ in Eq. (2) should be predicted by an appropriate full theory where the mass-mechanism is known. On the other hand, Eq. (4) shows how $`M_\nu `$ can be fully reconstructed in terms of nine directly measurable quantities, the mass-eigenvalues, the mixing angles and the CP-phases. Before knowing the underlying full theory, this suggests an important and reliable bottom-up approach, namely, we ask: given the existing neutrino experiments, can we construct a simple, realistic $`M_\nu `$ with only a minimal set of input parameters which describes all the oscillation data? To be concrete, we will focus on the Type-B scenario with inverted mass hierarchy in Eq. (1).222For a very recent analysis of the normal hierarchy (Type-A) via a bottom-up approach, see Ref. . We will show that such a minimal scheme can be quite uniquely derived and is highly predictive. We can, of course, further extend or elaborate the Minimal Scheme with more fine structure and more input parameters if that is needed to match with an underlying theory (once specified). However, the essential structure of the Minimal Scheme and its capability for describing the existing oscillation data 333The result from Liquid Scintillation Neutrino Detector (LSND)  awaits confirmation by the Fermilab miniBooNE experiment  and will not be considered in the present study. will remain in any realistic extension.

2.1. Minimal Scheme of Type-B1

The neutrino mass matrices of inverted hierarchy (Type-B) can be classified into Type-B1 and -B2  which we will analyze in turn. We start from the simplest, naive mass matrix $`M_{\nu 0}`$ of Type-B1 ,

$$M_{\nu 0}[\mathrm{B1}]=\frac{\overline{m}_0^{}}{\sqrt{2}}\left(\begin{array}{ccc}\mathrm{\hspace{0.17em}0}& \mathrm{\hspace{0.17em}1}& \mathrm{\hspace{0.17em}1}\\ \mathrm{\hspace{0.17em}1}& \mathrm{\hspace{0.17em}0}& \mathrm{\hspace{0.17em}0}\\ \mathrm{\hspace{0.17em}1}& \mathrm{\hspace{0.17em}0}& \mathrm{\hspace{0.17em}0}\end{array}\right),$$
(6)
which generates a mass spectrum $`(m_1^{},m_2^{},m_3^{})=(m_1^{},m_2^{},m_3^{})=\overline{m}_0^{}(1,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}0})`$ and exact bi-maximal mixing, $`\theta _1=\theta _2=45^{}`$. This simple structure (6) is motivated by $`L_eL_\mu L_\tau `$ symmetry444For some recent non-minimal approaches with certain $`U(1)`$ flavor symmetry and additional fields, see Ref. .. (It was also shown to be generic for the minimal radiative Zee-model and its various extensions .) However, (6) is not realistic and is excluded by the solar oscillation data since it predicts $`\mathrm{\Delta }_{}=|m_1^2m_2^2|=0`$, and, more seriously, a maximal solar angle $`\theta _{}=\theta _1=45^{}`$ which is difficult to reconcile with the MSW-LAM . We observe that such a failure is due to the small but nonzero ratio of the two measured mass-square differences $`\mathrm{\Delta }_{}/\mathrm{\Delta }_{\mathrm{atm}}=𝒪(10^110^2)`$ and a moderate angular deviation $`(45^{}\theta _1)/\theta _1(45^{}32^{})/32^{}0.4`$, for MSW-LAM. Therefore, it is justified to take $`M_{\nu 0}`$ as our zeroth order mass matrix and build in the necessary Minimal Perturbations to make a realistic Type-B1 neutrino mass-matrix $`M_\nu =M_{\nu 0}+\mathrm{\Delta }M_\nu `$. Such minimal perturbations represent proper $`L_eL_\mu L_\tau `$ violation effects. What is the minimal set of extra parameters which we need for a realistic perturbation $`\mathrm{\Delta }M_\nu `$? First, we need a sizable parameter $`\kappa =𝒪(0.5)`$ to accommodate the solar angular deviation of $`(45^{}\theta _1)/\theta _10.4`$; second, we need a small parameter $`\delta ^{}|m_1^{}m_2^{}|/\overline{m}_0^{}𝒪(0.1)`$ to account for the minor mass ratio $`\mathrm{\Delta }_{}/\mathrm{\Delta }_{\mathrm{atm}}=𝒪(10^110^2)`$; finally, to ensure the Type-B mass spectrum (1) we should impose a condition $`|m_1^{}m_2^{}|/\overline{m}_0^{}m_3^{}/\overline{m}_0^{}=𝒪(\delta ^{})`$, which can be naturally realized only if we introduce an “interplay” parameter $`\delta `$ lying between $`\kappa `$ and $`\delta ^{}`$. In summary, to construct a realistic perturbation to $`M_{\nu 0}`$, we have to start with three dimensionless parameters $`(\kappa ,\delta ,\delta ^{})`$ satisfying the proper hierarchies,

$$\begin{array}{c}1>|\kappa |>|\delta |>|\delta ^{}|,|\kappa ||\delta ^{}|,\hfill \\ m_3^{}/\overline{m}_0^{}|m_1^{}m_2^{}|/\overline{m}_0^{}=𝒪(\delta ^{}).\hfill \end{array}$$
(7)
With these, we can almost uniquely determine the pattern of the perturbation $`\mathrm{\Delta }M_\nu `$, and derive the following Minimal Scheme-B1:

$$M_\nu [\mathrm{B1}]=\frac{\overline{m}_0^{}}{\sqrt{2}}\left(\begin{array}{ccc}\kappa & 1& 1\\ 1& \kappa & \delta \\ 1& \delta & \delta ^{}\end{array}\right).$$
(8)
The relative sign between 11- and 22-entry is uniquely fixed by the requirement $`|m_1^{}m_2^{}|/\overline{m}_0^{}=𝒪(\delta ^{})`$. Note that to affect $`\theta _1`$, $`\kappa `$ cannot be put in 12- and 21-entry as $`M_\nu `$ is symmetric. Another reason to arrange the 11-entry to be of $`𝒪(\kappa )`$ rather than $`𝒪(\delta ,\delta ^{})`$ comes from the generic observation about the nature of $`m_{ee}`$ by using the Type-B mass spectrum (1) and the general Eq. (5),

$$m_1^{}c_3^2|m_{ee}|m_1^{}c_3^2|\mathrm{cos}2\theta _1|m_1^{}|\mathrm{cos}2\theta _1|,$$
(9)
where the upper \[lower\] bound corresponds to the CP-conserving values of the Majorana phase $`\varphi ^{}=0,\mathrm{or},\pi [\pi /2,\mathrm{or},3\pi /2]`$ . For solar mixing within the 95% C.L. range, $`25^{}\theta _139^{}`$ (MSW-LAM), we deduce the lower limit,

$$|m_{ee}|/m_1^{}|m_{ee}|/\overline{m}_0^{}(0.210.64),$$
(10)
so that we can identify $`\sqrt{2}|m_{ee}^{}|/\overline{m}_0^{}=𝒪(\kappa )`$. It is important to note that, for general Type-B scenarios, the significant deviation of $`0.15(45^{}\theta _1)/\theta _10.8`$ for $`25^{}\theta _139^{}`$  à la MSW-LAM already requires a sizable $`m_{ee}^{}`$ which is potentially observable via $`0\nu \beta \beta `$-decay experiments , depending on the overall scale $`\overline{m}_0^{}`$ . As will be shown in Sec. 3.1, due to the condition (7) and the smallness of $`s_3^2`$, we can reduce the dimensionless inputs $`(\kappa ,\delta ,\delta ^{})`$ down to a single parameter $`\kappa `$ (or, equivalently, the solar angle $`\theta _{}\theta _1`$). This makes our minimal scheme-B1 highly predictive.

After a scan for all possible variations of the minimal Scheme-B1 under the condition (7), we find a few other acceptable minimal schemes with $`M_\nu (\overline{m}_0^{}/\sqrt{2})𝔐`$ and $`𝔐`$ given by

$$(\begin{array}{ccc}\kappa & 1& 1\\ 1& \delta ^{}& \delta \\ 1& \delta & \kappa \end{array}),(\begin{array}{ccc}\kappa & 1\delta ^{}& 1\\ 1\delta ^{}& \kappa & \delta \\ 1& \delta & 0\end{array}),\left(\begin{array}{ccc}\kappa & 1\delta ^{}& 1\\ 1\delta ^{}& 0& \delta \\ 1& \delta & \kappa \end{array}\right).$$
(11)
Here the first matrix has the same mass eigenvalues as Eq. (8); its rotation angles $`(\theta _3,\theta _2)`$ contain a sign flip for the small $`𝒪(s_3)`$ terms. The third matrix in Eq. (11) is a variation of the first matrix by relocating $`\delta ^{}`$; similarly, the second matrix above is constructed by relocating $`\delta ^{}`$ from Eq. (8). Hence, Eq. (11) differs from Eq. (8) only by small terms of $`𝒪(s_3,\delta ^{})`$, with no conceptual difference. We will focus on the minimal scheme (8) hereafter. We also note that all the realistic minimal schemes we find for the Type-B can have at most one independent texture zero \[cf. the above Eq. (11) and the following Eq. (13) with $`\xi ^{}=0`$ or $`\xi =0`$\]. A recent interesting analysis  classified viable schemes with two independent texture-zeros, which, as expected, do not contain Type-B schemes.

2.2. Minimal Scheme of Type-B2

The naive form of Type-B2 is defined as ,

$$M_{\nu 0}[\mathrm{B2}]=m_0^{}\left(\begin{array}{ccc}1& 0& 0\\ 0& 1/2& 1/2\\ 0& 1/2& 1/2\end{array}\right),$$
(12)
which has a Type-B mass-spectrum $`(m_1^{},m_2^{},m_3^{})=(m_1^{},m_2^{},m_3^{})=m_0^{}(1,1,0)`$, a maximal mixing angle $`\theta _2=45^{}`$, and vanishing $`(\theta _1,\theta _3)`$. To be realistic, the zeroth order matrix (12) has to be properly perturbed for generating the observed small but nonzero $`\mathrm{\Delta }_{}=|m_1^2m_2^2|`$ and the large (rather than maximal) mixing $`\theta _1`$ for MSW-LAM. Using the fact of $`m_1^{}m_2^{}m_3^{}`$ and the general relation (5), we find all the perturbations in $`\mathrm{\Delta }M_\nu =M_\nu M_{\nu 0}`$ to be of $`𝒪(m_1^{}m_2^{},m_3^{},s_3^2m_0^{})`$ or smaller. This suggests a completely different perturbation structure from Type-B1, namely, we use only small perturbation parameters of $`𝒪(|m_1^{}m_2^{}|/m_0^{},m_3^{}/m_0^{})`$ and the large solar mixing angle $`25^{}\theta _139^{}`$ (95% C.L.) can be naturally generated from an $`𝒪(1)`$ ratio of two small perturbation parameters. Inspecting the structure of $`\mathrm{\Delta }M_\nu =M_\nu M_{\nu 0}`$ for Type-B2 and using Eq. (5), we are quite uniquely led to the following Scheme-B2,

$$M_\nu [\mathrm{B2}]=m_0^{}\left(\begin{array}{ccc}1+\delta & \xi & \xi ^{}\\ \xi & 1/2& 1/2\\ \xi ^{}& 1/2& 1/2\end{array}\right),$$
(13)
where we impose

$$\frac{m_3^{}}{m_0^{}}\frac{|m_1^{}m_2^{}|}{m_0^{}}=𝒪(\delta ,\xi ,\xi ^{})1.$$
(14)
We can define a truly Minimal Scheme-B2 by setting $`\xi ^{}=0`$ . Choosing $`\xi ^{}=\xi `$ will give, $`\theta _3=𝒪(\delta ^2,\xi ^2,s_3^2,\delta s_3^{},\xi s_3^{})0`$, and is consistent with the Chooz bound . But Type-B2 generally has negligible $`s_3^{}`$ even for nonzero $`\xi ^{}\xi `$. In Eq. (13), there is no need to perturb the $`2\times 2`$ block of $`\nu _\mu \nu _\tau `$ as the maximal mixing is favored by the atmospheric data ; also the $`2\times 2`$ block of $`\nu _e\nu _\mu `$ invokes two small perturbations so that an $`𝒪(1)`$ ratio is generated to explain the non-maximal solar mixing (MSW-LAM). Unlike the Type-B1, the scheme-B2 has larger $`m_{ee}^{}=m_0^{}(1+\delta )m_0^{}`$ at the zeroth order and is more sensitive to the $`0\nu \beta \beta `$ experiments.

3. Analysis of the Minimal Schemes and Predictions for Neutrino Oscillation

In this section, we systematically solve the diagonalization equations in (5) for the Minimal Scheme-B1 (8) and -B2 (13) \[$`\xi ^{}=0`$\] with CP-conservation. We then study their predictions for the neutrino oscillations.

3.1. Analyzing the Minimal Scheme of Type-B1

The parameters $`(\kappa ,\delta )`$ will be retained up to all orders without approximation. But, from the solar and Chooz oscillation data, it is justified to treat the small parameters $`(\delta ^{},s_3^{})`$ as perturbations to first power and ignore terms of $`𝒪(\delta _{}^{}{}_{}{}^{2},s_3^2)𝒪(10^2)`$ or smaller. As will be shown below, the expansion of $`s_3^{}`$ also plays a key role for eliminating $`\delta ^{}`$ from inputs.

From Eq. (8), we deduce the mass-eigenvalues of $`M_\nu `$, up to $`𝒪(\overline{\delta ^{}})`$,

$$m_{1,2}^{}=m_0^{}\left[1\frac{1}{2}\left(1\frac{x}{\mathrm{\hspace{0.17em}2}+\omega }\right)\overline{\delta ^{}}\right],m_3^{}=m_0^{}\frac{x}{\mathrm{\hspace{0.17em}2}+\omega }\overline{\delta ^{}},$$
(15)
where we expand $`\overline{\delta ^{}}`$ to first order and define,

$$\begin{array}{c}\overline{\kappa }\frac{\kappa }{\overline{\omega }},\overline{\delta }\frac{\delta }{\overline{\omega }},\overline{\delta ^{}}\frac{\delta ^{}}{\overline{\omega }},\overline{\omega }\sqrt{2+\omega }=\sqrt{\frac{2}{1\overline{\kappa }^2\overline{\delta }^2}},\hfill \\ \\ \omega \kappa ^2+\delta ^2=\frac{\mathrm{\hspace{0.17em}2}(\overline{\kappa }^2+\overline{\delta }^2)}{\mathrm{\hspace{0.17em}1}(\overline{\kappa }^2+\overline{\delta }^2)},m_0^{}\overline{m}_0^{}\sqrt{1+\frac{\omega }{2}}.\hfill \end{array}$$
(16)
The parameter $`x=𝒪(1)`$ in Eq. (15) will be determined by the consistency condition,

$$\left(1\delta ^2\right)\kappa 2\delta +\left(1+\kappa ^2\right)\delta ^{}=x\delta ^{},$$
(17)
due to the requirement $`m_3^{}/\overline{m}_0^{}=𝒪(\delta ^{})`$ in Eq. (7). With the definition of Eq. (16), we can rewrite the neutrino mass matrix (8) scaled by $`m_0^{}`$,

$$M_\nu [\mathrm{B1}]=m_0^{}\left(\begin{array}{ccc}\overline{\kappa }& \overline{\omega }^1& \overline{\omega }^1\\ \overline{\omega }^1& \overline{\kappa }& \overline{\delta }\\ \overline{\omega }^1& \overline{\delta }& \overline{\delta ^{}}\end{array}\right).$$
(18)
From Eq. (15), we deduce

$$\frac{\mathrm{\Delta }_{}}{\mathrm{\Delta }_{\mathrm{atm}}}=\frac{|m_1^2m_2^2|}{|m_{1,2}^2m_3^2|}\mathrm{\hspace{0.17em}2}\left[1\frac{x}{2+\omega }\right]\overline{\delta ^{}}.$$
(19)
With the mass-eigenvalues given, we can then solve for the mixing angles by substituting $`M_\nu `$ \[cf. Eq. (8) or Eq. (18)\] into the six diagonalization equations in (5) and expanding $`(\delta ^{},s_3^{})`$ systematically to first order. At this order, we find that only five of the six equations are independent. Note that we have three dimensionless parameters $`(\kappa ,\delta ,\delta ^{})`$ in $`M_\nu `$ (in which the overall scale $`\overline{m}_0^{}`$ is irrelevant to the diagonalization) and three mixing angles $`(\theta _1,\theta _2,\theta _3)`$. Hence, from the five equations, we can solve five out of the six parameters as functions of a single dimensionless input parameter which will be chosen as the angle $`\theta _1`$ (measured in the solar oscillation). To explicitly understand this nontrivial reduction of input parameters, we first note that even though we have three dimensionless inputs $`(\kappa ,\delta ,\delta ^{})`$ in (8), the condition $`m_3^{}/\overline{m}_0^{}=𝒪(\delta ^{})`$ in (7) \[or, (17)\] relates $`\kappa `$ and $`\delta `$ at zeroth order of $`\delta ^{}`$ so that only two inputs among $`(\kappa ,\delta ,\delta ^{})`$ are independent under the expansion of $`\delta ^{}`$. Then, we summarize two relevant relations derived from Eq. (5) \[and Eq. (18)\],

$$\begin{array}{ccc}\overline{\delta }^2\hfill & =& (1r)\overline{\kappa }\overline{\delta ^{}}+𝒪(s_3^2,\overline{\delta ^{}}^2),\hfill \\ \overline{\kappa }\hfill & =& 2\overline{\delta }(12r)\overline{\delta ^{}}+𝒪(s_3^2),\hfill \end{array}$$
(20)
with $`r{\displaystyle \frac{x}{2+\omega }}`$. Now we see that the absence of $`𝒪(s_3^{})`$ term in Eq. (20) and the smallness of $`s_3^2(0.04`$ ) lead us to have three constraints \[two in Eq. (20) and one in Eq. (17)\] among the four parameters $`(\overline{\kappa },\overline{\delta },\overline{\delta ^{}},x)`$ . This feature remains if we include higher order terms via iteration. This makes our scheme end up with a single input for all mixings and thus extremely predictive. After a lengthy and careful derivation, we arrive at the following complete set of solutions of our Minimal Scheme-B1, up to $`𝒪(\overline{\delta ^{}},s_3^{})`$ ,

$$\begin{array}{c}\overline{\kappa }=\frac{8}{9}\mathrm{cos}2\theta _1,\overline{\delta }=\frac{1}{2}\overline{\kappa },\overline{\delta ^{}}=\frac{\overline{\kappa }}{4(1r)},\hfill \\ \theta _2=\frac{\pi }{4}\frac{34r}{8(1r)}\overline{\kappa }^2,\theta _3s_3^{}=\frac{34r}{8(1r)}\overline{\kappa }\sqrt{1\overline{\kappa }^2},\hfill \end{array}$$
(21)
where we have,

$$2+\omega =\frac{2}{1\frac{5}{4}\overline{\kappa }^2},r=\frac{1\frac{5}{4}\overline{\kappa }^2}{\mathrm{\hspace{0.17em}2}(1\overline{\kappa }^2)},x=\frac{1}{\mathrm{\hspace{0.17em}1}\overline{\kappa }^2}.$$
(22)
Finally, inputting the solar angle $`\theta _{}(\theta _1)`$, we deduce the following numerical predictions,

$$\begin{array}{c}25^{}\theta _139^{},[\mathrm{Input}\mathrm{of}\mathrm{MSW}\mathrm{LAM},95\%\mathrm{C}.\mathrm{L}.];\hfill \\ 39.8^{}\theta _244.5^{},0.13s_3^{}\theta _30.046,\hfill \\ 0.57\overline{\kappa }0.19,0.29\overline{\delta }0.092,0.25\overline{\delta ^{}}0.091;\hfill \\ 0.28\frac{\mathrm{\Delta }_{}}{\mathrm{\Delta }_{\mathrm{atm}}}0.095,\hfill \end{array}$$
(23)
and also $`0.44r0.50,\mathrm{\hspace{0.17em}1.48}x1.01.`$ The results in Eq. (23) agree well with the oscillation data , i.e., $`33^{}\theta _2\theta _{\mathrm{atm}}57^{}`$ and $`s_3^{}=\mathrm{sin}\theta _{\mathrm{chz}}0.2`$ at 95% C.L., and $`4.2[3.3]\times 10^3\mathrm{\Delta }_{}/\mathrm{\Delta }_{\mathrm{atm}}\mathrm{\hspace{0.17em}0.17}[0.32]`$ at 95%\[99%\] C.L. The on-going KamLAND experiment  will more precisely test the MSW-LAM parameter space (though it will not be sensitive to $`s_3^{}`$ ). We have checked the numerical accuracy of the above solutions by substituting Eq. (23) back into Eqs. (5) and (8) and evaluating the difference of the two sides in each equation. We find that the difference (uncertainty) is always less than $`0.015(0.0009)`$ for $`\theta _1=25^{}(39^{})`$. Thus, our systematical expansion works well up to $`𝒪(\overline{\delta ^{}},s_3^{})`$, as expected, since the ignored terms are of $`𝒪(\overline{\delta ^{}}^2,s_3^2)`$ and become smaller for larger $`\theta _1`$ as shown in Eq. (23). For Type-B schemes, the mass scale $`m_0^{}`$ is generally bounded by

$$0.036\mathrm{eV}m_0^{}m_{1,2}^{}\mathrm{\Delta }_{\mathrm{atm}}^{1/2}\mathrm{\hspace{0.17em}0.074}\mathrm{eV},[99\%\mathrm{C}.\mathrm{L}.].$$
(24)
3.2. Analyzing the Minimal Scheme of Type-B2

We now turn to the minimal scheme-B2 in (13) with $`\xi ^{}=0`$, which has the mass-eigenvalues, up to $`𝒪(\delta ,\xi )`$,

$$m_{1,2}^{}=m_0^{}\left[1+\frac{\delta }{2}\pm \frac{1}{2}\sqrt{\delta ^2+2\xi ^2}\right],m_3^{}=0.$$
(25)
Substituting (13) into (5), expanding up to $`𝒪(\delta ,\xi ,s_3^{})`$ and using (25), we derive the solutions,

$$\begin{array}{c}\theta _1=\frac{\pi }{4}\frac{1}{2}\mathrm{arcsin}\frac{\delta }{\sqrt{\delta ^2+2\xi ^2}},\theta _2=\frac{\pi }{4},\theta _3=\frac{\xi }{\sqrt{2}},(\delta ,\xi )>0;\hfill \\ \mathrm{and}\frac{\mathrm{\Delta }_{}}{\mathrm{\Delta }_{\mathrm{atm}}}=\frac{|m_1^2m_2^2|}{|m_{1,2}^2m_3^2|}=2\sqrt{\delta ^2+2\xi ^2}.\hfill \end{array}$$
(26)
We see that $`\theta _2`$ is maximal at this order. The sizable deviation of $`\theta _1{\displaystyle \frac{\pi }{4}}`$ is indeed naturally generated by an $`𝒪(1)`$ ratio of two small parameters $`(\delta ,\xi )1`$. \[Allowing $`\xi ^{}0`$, the corresponding formulas for Eqs. (25)-(26) can be directly obtained by the simple replacements, $`\xi \xi \xi ^{}`$ for $`\theta _3`$ and $`\xi \xi +\xi ^{}`$ for all other quantities.\] Using the inputs for LAM , $`25^{}\theta _{}\theta _139^{}`$ and $`4.2\times 10^3\mathrm{\Delta }_{}/\mathrm{\Delta }_{\mathrm{atm}}0.17`$ at 95% C.L., we deduce,

$$1.2\delta /\xi \mathrm{\hspace{0.17em}0.3},1.1\times 10^3\xi \mathrm{\hspace{0.17em}0.06}.$$
(27)
Thus we have

$$8\times 10^4\theta _3\mathrm{\hspace{0.17em}0.04}.$$
(28)
The mass scale $`m_0^{}\mathrm{\Delta }_{\mathrm{atm}}^{1/2}`$ is bounded as in Eq. (24). As mentioned above, allowing nonzero $`\xi ^{}=𝒪(\xi )`$, we can derive,

$$s_3^{}\theta _3=\frac{\xi ^{}\xi }{\sqrt{2}},$$
(29)
which remains of the same order. Hence, $`|\theta _3|𝒪(10^2)`$ generally holds for Type-B2, implying negligible CP-violation from the Dirac phase $`\varphi `$ . It has been shown  that combining the data from two near-future long baseline accelerator experiments , the Main Injector Neutrino Oscillation Search (MINOS) and the Imaging Cosmic And Rare Underground Signals (ICARUS), may place a 95% C.L. lower bound, $`\theta _3𝒪(0.05)`$ (when $`\theta _3`$ lies within their combined sensitivity), which could possibly discriminate Type-B2 from Type-B1.

4. Implications for Neutrinoless Double $`\beta `$-Decay, Tritium $`\beta `$-Decay
and Cosmology

The oscillation data may already give a strong hint on the neutrino mass scale \[cf. Eq. (24)\] so long as the neutrino masses exhibit the hierarchy structure \[cf. Eq. (1)\], but the possibility of three nearly degenerate neutrinos ($`m_1^{}m_2^{}m_3^{}`$) could allow a higher scale. Hence, the laboratory experiment on neutrinoless double $`\beta `$-decay ($`0\nu \beta \beta `$)  is indispensable to pin down the absolute mass scale, as well as the Majorana nature of active neutrinos. For the Minimal Scheme-B1 and -B2, we have,

$$m_{ee}^{}[\mathrm{B1}]=\overline{\kappa }m_0^{},m_{ee}^{}[\mathrm{B2}]m_0^{}.$$
(30)
Thus, using Eqs. (30) and (24), we derive, at 99% C.L.,

$$0.014\mathrm{eV}m_{ee}^{}[\mathrm{B1}]0.029\mathrm{eV},0.036\mathrm{eV}m_{ee}^{}[\mathrm{B2}]0.074\mathrm{eV},$$
(31)
where we input the central value of solar fit (LAM), $`\theta _132^{}`$, for Type-B1. The $`m_{ee}^{}[\mathrm{B2}]`$ is already sensitive to the current $`0\nu \beta \beta `$ measurement . The experiments of $`0\nu \beta \beta `$ decay , such as the on-going Nemo3 and the upcoming Cuore, can probe $`|m_{ee}^{}|0.1`$ eV, while the near-future measurements at Genius, Exo, Majorana and Moon aim at a sensitivity of $`|m_{ee}^{}|0.01`$ eV, which is decisive for testing the whole mass range (31) of Type-B1 and -B2 schemes.

Tritium $`\beta `$-decay requires , $`m_{\nu _e}^{}<2.2`$ eV, at 95% C.L., where $`m_{\nu _e}^{}(M_\nu ^{}M_\nu )_{ee}^{1/2}m_{1,2}^{}m_0^{}`$ for Type-B. This is well above the range given in Eq. (24). The sensitivity of $`H^3`$ $`\beta `$-decay could eventually reach $`m_{\nu _e}^{}0.5`$ eV .

The latest cosmology measurements of the power spectrum for the Cosmic Microwave Background (CMB), Galaxy Clustering and Lyman Alpha Forest  put a 95% C.L. upper bound on the neutrino masses, $`_jm_j^{}4.2`$ eV. This gives, for our Type-B schemes, $`m_1^{}m_2^{}m_0^{}2.1`$ eV, which is about the same as the tritium $`\beta `$-decay bound. The newest analysis  from the 2dF Galaxy Redshift Survey arrives at an upper bound, $`_jm_j^{}2.2`$ eV, which results in, $`m_1^{}m_2^{}m_0^{}1.1`$ eV, for Type-B schemes. Stronger constraints of $`_jm_j0.4`$ eV are expected from the forthcoming Microwave Anisotropy Probe (MAP) and Planck satellite experiments .

It is interesting to note that the neutrino mass scale may also be determined from the so-called $`Z`$-bursts  due to the resonant annihilation of ultra high energy neutrinos with cosmological relic (anti-)neutrinos into $`Z`$ bosons (whose decay produces protons and photons). In the most plausible case where the ordinary cosmic rays are protons of extragalactic origin, the required neutrino mass range is ,

$$0.01[0.02]\mathrm{eV}m_\nu ^{}(\mathrm{heaviest})\mathrm{\hspace{0.17em}3.0}[2.1]\mathrm{eV},\mathrm{at}99\%\mathrm{C}.\mathrm{L}.[95\%\mathrm{C}.\mathrm{L}],$$
(32)
which is compatible with the neutrino oscillation bound (24) for the inverted mass hierarchy.

5. Conclusions

In this study, we have considered two essential and distinct scenarios for the neutrino Majorana mass matrix with inverted hierarchy, called Type-B1 \[cf. Eq. (6)\] and -B2 \[cf. Eq. (12)\]. For Type-B1, we start with the form of Eq. (6) at zeroth order and perturb it into the realistic form of Eq. (8) with three parmeters $`(\kappa ,\delta ,\delta ^{})`$ under the hierarchy (7) that is necessary to correctly predict the oscillation data, especially, the non-maximal solar neutrino mixing of MSW-LAM. The sizes of $`\delta ^{}|m_1^{}m_2^{}|/m_0^{}\mathrm{\Delta }_{}/\mathrm{\Delta }_{\mathrm{atm}}1`$  and $`s_3^{}=\mathrm{sin}\theta _{\mathrm{chz}}1`$  justify the expansion of $`(\delta ^{},s_3^{})`$, which enables us to reduce the number of inputs down to a single parameter $`\kappa `$, or, equivalently, $`\theta _{}(\theta _1)`$. Thus, using only the measured solar angle $`\theta _{}`$ as input, we predict the atmospheric mixing angle, $`\theta _{\mathrm{atm}}(\theta _2)`$, the value of $`\theta _{\mathrm{chz}}(=\theta _3)`$, and the mass ratio $`\mathrm{\Delta }_{}/\mathrm{\Delta }_{\mathrm{atm}}`$, in complete agreement with the existing data. We also note that the Minimal Scheme-B1 in Eq. (8) points to a generic way for naturally extending the minimal Zee-model  in which $`M_\nu `$ exhibits the following structure,

$$M_\nu ^{\mathrm{Zee}}=\left(\begin{array}{ccc}0& m_{e\mu }^{}& m_{e\tau }^{}\\ m_{e\mu }^{}& 0& m_{\mu \tau }^{}\\ m_{e\tau }^{}& m_{\mu \tau }^{}& 0\end{array}\right),$$
(33)
where the pattern $`m_{e\mu }^{}m_{e\tau }^{}m_{\mu \tau }^{}`$ can be realized , which ensures an approximate $`L_eL_\mu L_\tau `$ symmetry. The necessity of modifying Eq. (33) in the minimal Zee-model for accommodating the MSW-LAM was noted recently . Our minimal construction of Scheme-B1 in Eq. (8) demonstrates a generic way to extend $`M_\nu ^{\mathrm{Zee}}`$ under an appropriate perturbation \[cf. Eq. (7)\]. The choice of $`m_{ee}^{}=(\overline{m}_0^{}/\sqrt{2})\kappa =m_0^{}\overline{\kappa }`$ is due to the general observation in Eq. (9) for Type-B1 and a sizable $`\mathrm{cos}2\theta _1(0.210.64)`$ for $`25^{}\theta _139^{}`$ (95% C.L., MSW-LAM); while $`m_{\mu \mu }^{}m_{ee}^{}`$ is enforced by the Type-B mass-spectrum, $`m_{1,2}^{}|m_1^{}m_2^{}|m_3^{}`$. The $`\kappa `$ terms in $`M_\nu [\mathrm{B1}]`$ represent the generic leading modification to the minimal Zee-model (33).

For Type-B2, we start with the leading order mass matrix (12) and find the perturbation structure in Eq. (13) based on the general relations in Eq. (5) and the smalless of $`\theta _{\mathrm{chz}}(=\theta _3)`$ . The realistic Minimal Scheme-B2 contains only two small parameters $`(\delta ,\xi )1`$, defined as in Eq. (13), when $`\xi ^{}=0`$ . In contrast to Typy-B1, the non-maximal solar mixing angle $`\theta _{}`$ is naturally accommodated by a ratio of two small parameters, $`\delta /\xi =𝒪(1)`$, while the atmospheric mixing angle remains maximal. Using the measured values of solar angle $`\theta _{}`$ and mass ratio $`\mathrm{\Delta }_{}/\mathrm{\Delta }_{\mathrm{atm}}`$, we derive the ranges for $`\theta _3(=\theta _{\mathrm{chz}})`$ and the perturbation parameters $`(\delta ,\xi )`$. The angle $`\theta _3`$ is found to be of $`𝒪(10^2)`$ or smaller. Combining the data from both MINOS and ICARUS experiments  may reach the sensitivity  to discriminate between the minimal Type-B1 and Type-B2 schemes.

The overall neutrino mass scale $`m_0`$ for the inverted mass hierarchy is quite uniquely fixed by the atmospheric neutrino data on the mass-squared difference $`\mathrm{\Delta }_{\mathrm{atm}}=|m_{1,2}^2m_3^2|m_{1,2}^2m_0^2`$ \[cf. Eq. (24)\]. Thus, the mass matrix of our minimal scheme-B1 or -B2 is known and highly predictive. Some implications of the Type-B1 and -B2 minimal schemes for the neutrinoless double $`\beta `$-decay, tritium $`\beta `$-decay and cosmology are given above.

Note Added: After the submission of this work, a new announcement  from the Sudbury Neutrino Observatory (SNO) collaboration appeared on April 20, 2002, which further confirms the MSW-LAM as the best solution to the solar neutrino oscillations.

Acknowledgments
We are happy to thank Vernon Barger, Danny Marfatia and Alexei Yu. Smirnov for valuable discussions. This work was supported in part by U.S. Department of Energy under grant DE-FG03-93ER40757 and by Natural Science and Engineering Research Council of Canada.