text:
# Young compact objects in the Solar vicinity

## 1 Introduction

Neutron stars (NSs) and black holes (BHs) are among the most interesting astrophysical sources. Usually NSs are observed as radio pulsars or as accreting objects in close binaries. Similarly, stellar mass BHs are observed when they accrete matter from a companion star. Here we focus on much more elusive sources, namely isolated NSs (which may show no radiopulsar activity) and isolated BHs.

An isolated NS can be relatively bright in soft X-rays due to its thermal emission during the first Myrs of its life, when it is still hot ($`T10^6`$ K) in the aftermath of the supernova (SN) explosion. Such objects are observed in the Solar proximity and in SN remnants . Older NSs (that is to say those which crossed the deathline in $`10^7`$ yr) are not expected to emit appreciable amounts of electromagnetic radiation in any energy band. However, accretion of the interstellar medium (ISM) may make them shine again as soft, faint X-ray sources (see e.g. Treves et al. ). Much in the same way, an isolated BH may be detected if it accretes from the ISM, or, possibly, revealed through microlensing .

In this paper we discuss the possible origin of close-by, isolated NSs and present some evidence that the seven radio-quiet, thermally emitting ROSAT sources (the “magnificent seven” ) may be characterized by different values of the stellar parameters (mass and magnetic field) with respect to ordinary radiopulsars. We first construct the Log N – Log S distribution for young close-by isolated NSs and compare it with present observations of close-by young isolated NSs of all types. Then we discuss how the alignment of the magnetic and rotation axes, together with the role of fall-back following the supernova event, may help in explaining the observed parameters of NSs. Finally in §3 we discuss how one can estimate an approximate positions of close-by young isolated BHs in the light of their possible detection at X/gamma-rays energies.

## 2 Isolated neutron stars

In this section we discuss isolated NSs. The material presented here is partly based on the results published in Popov et al. . In addition some recent results are also included .

### 2.1 Origin of close-by isolated NSs

To understand how the local population of isolated NSs originated, we construct their Log N – Log S distribution. The main components of our model are (see ): spatial distribution of NS progenitors, NS formation rate, NS cooling history, and a model of interstellar absorption (that is to say the spatial distribution of the ISM). In addition we calculate the dynamical evolution of NSs in the galactic potential. In brief our model can be described in the following way: NSs are born in the Galactic plane and in the Gould Belt (a local compound of stellar associations, see below); at birth they receive a kick velocity; we then follow the evolution of NSs in the Galactic potential; finally, we calculate the ROSAT count rate basing on cooling curves and an assumed model of interstellar absorption.

NSs are considered to be born with a constant rate: 20 NSs per Myr come from the Gould Belt, and 250 NSs per Myr from the Galactic plane (up to a limiting distance of 3 kpc from the Sun) with a uniform distribution. The Gould Belt is modeled as a disk of 500 pc radius with an inclination of 18 with respect to the Galactic plane. Its center is situated at 100 pc from the Sun in the Galactic anticenter direction. The central region ($`150`$ pc in radius) is devoided of newborn NSs (see Pöppel  and Torra et al. ).

To calculate the thermal evolution of NSs we use the data obtained by Sankt-Petersburg group (see Kaminker et al. , and the review by Yakovlev et al. ). The NS cooling depends on the star mass and we adopt a flat mass spectrum in the range $`1.1M_{}<M<1.8M_{}`$. A more standard spectrum with a sharp maximum around 1.35-1.4 $`M_{}`$ gives nearly the same result. Cooling curves take into account all neutrino processes but ignore neutron superfluidity in the crust and core since this is not expected to influence the final results significantly. Calculations for each NS are truncated when its temperature drops to $`10^5`$ K; this corresponds to a NS age of 4.25 Myrs for the lightest NSs ($`M=1.1M_{}`$) or less for more massive objects.

Since we expect the NS to emit most of its luminosity at UV/soft X-ray energies ($`20200`$ eV or $`T10^5`$$`10^6`$ K) interstellar absorption plays a crucial role as far as the observability of these sources is concerned. Any attempt to estimate the amount of observable cooling isolated NSs using unabsorbed flux greatly overestimates their number.

Our main results are presented in Fig. 1 where we compare the Log N – Log S for NSs born in the Gould Belt and the Galactic disk. All curves refer to the whole sky. As can be seen the contribution from NSs born in the disk fails to explain the observed distribution while the inclusion of the objects originating from the Gould Belt alone can match the observations. Absorption, the flat geometry of NS initial distribution and the finite extension of the Belt naturally explain the very flat (slope $`<1`$) Log N – Log S distribution.

Our calculations show that there may be at most a few dozens of unidentified close-by isolated NSs in the ROSAT All-Sky Survey (at count rate $`>`$ 0.015 cts s-1) depending on parameters of the model. Also there may be a few unidentified ROSAT isolated NSs (RINSs) with fluxes $`>0.1`$ cts s-1 at low Galactic latitudes (see also Schwope et al. ). Most objects should be observed at $`|b|<`$ 20 towards the directions of lower absorption. Some of them can have counterparts among unidentified gamma-ray sources (also possibly connected with the Gould Belt, see Grenier ). Identification of these objects can be important for choosing a correct cooling model and for determination of the mass spectrum of NSs.

### 2.2 Census of close-by young NSs

At present about 20 NSs are known which are younger than 4.25 Myrs and closer than 1 kpc to the Sun (see the Table). They include: the “magnificent seven” (radio-quiet, ROSAT isolated NSs with only thermal emission), Geminga and the Geminga-like object 3EG J1835 (pulsars the beams of which do not intersect the Earth), the “three musketeers” (Vela, PSR 0656+14, PSR 1055-52), PSR 1929+10 and seven young radio pulsars, which have not been detected in X-rays yet.

In addition to the observed sources, we expect about one hundred isolated NSs younger than $``$4 Myr inside 1 kpc. These NSs are not detected as radio pulsars, but tens of them can be identified in ROSAT data as dim sources (others are too old to be hot enough). Pulsar beaming can be responsible only for a fraction of these undetected (in the radio) young NSs (about 50-70% of young pulsars are not visible from Earth ), and most of RINSs should be really radio silent. This provides strong support to the arguments by Gotthelf and Vasisht , that “at least half of the observed young neutron stars follow an evolutionary path quite distinct from that of the Crab pulsar”.

### 2.3 Are RINSs of a different stock ?

An interesting feature of RINSs population is the detection of periods in the $``$ 10-20 s range (typical of SGRs/AXPs, i.e. of magnetars) for four objects. Present data allow to exclude any pulsation (down to a few % fraction) for at least one source, RX J1856.5-3754. V. Beskin (2001, private communication) suggested this could be due to the alignment of magnetic and spin axes (see for example Tauris & Manchester  for a discussion). Alignment is a process which leads to “period freezing” and a low pulsed fraction.

However in the case of coolers alignment should operate on short timescale, since the star cools down in $`1`$ Myr. For radio pulsars the timescale of alignment is about 10 Myrs or longer , so it seems unlikely that this mechanism is responsible for RINSs distribution of the pulsed fraction, unless one assumes that RINSs form a separate population from normal radio pulsars. To illustrate this let us assume that the alignment timescale is $`\tau _{align}(\mathrm{\Omega }_0^2\mathrm{cos}^2\alpha _0B_0^2)^1`$ (here $`\alpha _0`$ and $`B_0`$ are initial values of an angle between spin and magnetic axis and of magnetic field). The previous expression comes from magnetodipolar braking supplemented by the condition $`\mathrm{\Omega }_0\mathrm{cos}\alpha _0=\mathrm{\Omega }\mathrm{cos}\alpha `$. To explain the difference between $`\tau _{align}`$ in RINSs and radio pulsars RINSs need to have a different distribution in $`B_0`$ and/or $`\alpha _0`$. In this respect RINSs may come from the same population as radiopulsars but are characterized by different average properties, like e.g. higher values of the magnetic field and relatively higher surface temperatures. The latter would imply a lower mass for NSs of the same age.

RINSs are currently thought to be rather highly magnetized objects (in RX J0720.4-3125 the detected spind-down implies $`B2\times 10^{13}`$ G ). If this is indeed the case, then one has to explain why their $`B`$-field is a factor $`10`$ higher than the average value in radiopulsars ($`2\times 10^{12}`$ G). A possibility is that NSs with higher magnetic fields are hotter. It is known that less massive NSs cools more slowly because direct URCA processes are not effective (see e.g. ). This means that among NSs of the same age the lighter are the hotter, so to test our hypotesis we need to show that lighter NSs may support a stronger field. Such a correlation arises quite naturally if more massive NSs get their additional mass from fall-back. In this case their magnetic field can be significantly suppressed , so more massive NSs should have lower initial magnetic fields. Besides, strong initial magnetic field together with fast rotation can prevent strong fall-back (this is especially possible if the magneto-rotational mechanism of supernova explosion is valid, see ). Again this leads to the same correlation between mass and filed strength, i.e. NSs with stronger fields would have lower masses. Also, this picture makes room for long initial spin periods. The study of this (and possible other) correlation in isolated NSs may prove very useful in understanding the correct mechanism of the SN event which gave birth to these objects.

In summarizing §2, we stress that future determination of RINSs parallax and proper motion may help in tracing back their kinematical history and derive their age. It can give a clue to their mass determination basing on cooling curves (see Kaminker et al. ). Our results suggest that the fraction of low-mass NSs ($`M<\text{ }`$1.3 $`M_{}`$) may not be small. On the other hand there should be room for NSs with $`M>\text{ }1.4M_{}`$, because otherwise the number of bright objects would be too large.

## 3 Close-by young isolated black holes

In this section we base on the results published in Popov et al.  and Prokhorov and Popov .

SNae explosions produce not only NSs, but also BHs. It is commonly accepted that BHs are one order of magnitude less abundant than NSs. This estimate comes from the critical mass for BH formation and follows if one assumes that progenitors more massive than about 35 $`M_{}`$ ended as BHs. Having dozens of SNae in the close solar vicinity during the last 10 Myr we can expect several BHs to have formed during the same period in the solar neighborhood.

At present, 56 runaway stars are known within $``$ 700 pc from the Sun. Only a few of them result from star-star interactions, so the majority comes from SNae explosions in binary systems. If the above considerations are correct we can expect about 5 BHs formed in about 50 disrupted binaries.

Close-by massive runaway stars give us a chance to calculate the approximate position of close-by young isolated BHs. Among runaway stars the most massive are $`\lambda `$ Cep, $`\zeta `$ Pup, HIP 38518 and $`\xi `$ Per (see Hoogerwerf et al. ). Since their mass is $`>\text{ }33`$ $`M_{}`$, the companion (actually the primary in the original binary) was even more massive on the main sequence stage. So, the most likely product of the explosion of such a massive star should be a BH.

If the present velocities of runaway stars are known, one can estimate their ages and places of birth. This has been done by Hoogerwerf et al. . To calculate the present position of a BH we have to know the binary parameters, i.e., the masses of stars before the explosion, the BH mass, the eccentricity of the orbit before the explosion, the orbit orientation, and finally the kick velocity of the BH. Some parameters can be inferred from the observation of the secondary star. We can assume a zero kick velocity for BHs and zero orbital eccentricity. Other parameters should be varied within assumed ranges (see details in Prokhorov and Popov .

We calculated approximate positions of isolated BHs for the four systems mentioned above and estimated error boxes where these BHs could be found. For $`\xi `$ Per and $`\zeta `$ Pup we obtained not very large error boxes inside each of which only one unidentified EGRET source is known. We suggest that these objects can be young isolated BHs. For the two other systems (HIP 38518 and $`\lambda `$ Cep) the present position of the BH is more uncertain and no definite conclusion on the possible detectability of the collapsed object can be drawn.

## 4 Conclusions

We have presented evidence that the seven radio-quiet ROSAT isolated NSs discovered so far can be connected with recent SNae explosions in the Gould Belt. These events produced nearby runaway stars and peculiar features in the local ISM including the Local Bubble. The relatively high local spatial density of young NSs is a natural consequence of the large number of massive progenitors in the Belt. The lack of a similar overabundance of active radiopulsars in the Solar vicinity lends further support to the claim that a large fraction ($`50\%`$ or more) of young NSs should be radio quiet. According to our results, the ROSAT All Sky Survey may contain about a few tens of unidentified RINSs. Moreover, it is possible that some unidentified RINSs with quite large flux ($`>0.1`$ cts s-1) are still hiding at low Galactic latitudes.

We also propose that massive runaway stars may be used to trace the present position of young close-by isolated BHs. Our calculations allowed to estimate with reasonable accuracy the positions of four such BHs. In two cases the error box is not too large and this may lead in the future to the positive identification of an isolated BH with X/gamma-rays observations.

## 5 Acknowledgments

We want to thank D.G. Yakovlev for the data on cooling curves and comments on them, and V.S. Beskin, M. Chieregato, A. Possenti and L. Zampieri for discussions. The work of S.P. was supported by the RFBR grant 02-02-06663 and by RSCI; that of M.P.by RFBR grant 01-15-99310. S.P. thanks Universities dell’Insubria and Milano-Bicocca for hospitality.


format: text

text:
# 1 Exotic scalar particles beyond the standard model and their couplings to combinations of two fermions.

UCRHEP-T345

October 2002

Baryon and Lepton Number Violation
with Scalar Bilinears
H. V. Klapdor-Kleingrothaus1, Ernest Ma2, and Utpal Sarkar3
1 Max-Planck-Institut für Kernphysik, D-69029 Heidelberg, Germany

2 Physics Department, University of California, Riverside, California 92521, USA
3 Physical Research Laboratory, Ahmedabad 380 009, India

## Abstract

We consider all possible scalar bilinears, which couple to two fermions of the standard model. The various baryon and lepton number violating couplings allowed by these exotic scalars are studied. We then discuss which ones are constrained by limits on proton decay (to a lepton and a meson as well as to three leptons), neutron-antineutron oscillations, and neutrinoless double beta decay.

Grand unified theories (GUTs) appear to be the most natural extensions of the standard model at the very high scale. In GUTs the gauge group is unified, so there is only one single gauge coupling constant to explain all the forces of the quarks and leptons, which are also treated at the same footing. However the scale of unification is much too high to be directly tested in the laboratory. Thus indirect searches become very important. Many predictions of GUTs were studied, but there is no supporting experimental evidence to date, with the possible exception of neutrinoless double beta decay . These GUT predictions also include new particles such as leptoquarks or diquarks and new interactions which violate baryon and lepton numbers. The latter are also required for an understanding of the present baryon asymmetry of the Universe. An observed Majorana neutrino mass also requires lepton number violation. In the context of these issues we present in this paper a general analysis of baryon and lepton number violation with scalar bilinears.

The scalar bilinears we consider are generalizations of the usual Higgs scalar doublet. They couple to two fermions of the standard model and are present in many of its extensions . The quantum numbers of these particles are thus fixed, but their masses and couplings could be arbitrary. Low-energy particle physics phenomenology (such as the non-observation of baryon number violation) could be used to constrain their masses and couplings . Baryogenesis with scalar bilinears has also been studied . Recently a general analysis of the neutrinoless double beta decay has been performed by constructing effective higher dimensional operators . Most of these operators may be realized in models where scalar bilinears are included, a few examples of which were demonstrated .

In Table 1 we list all the scalar bilinears allowed in any extension of the standard model. We then write down the different baryon and lepton number violating couplings of these scalars and what kind of interactions these different couplings may induce, leading to a simple classification of all the different possible baryon and lepton number violating interactions. We include the usual Higgs scalar doublet of the standard model in these interactions.

For the scalar bilinears, it is not possible to define definite baryon or lepton number in some cases. However, depending on their couplings to the combinations of two fermions, it is possible to assign baryon and lepton numbers separately. The different assignments are given in Table 2. For any particular interaction it will be possible to specify what processes it would mediate and what is the change in baryon and lepton numbers. This will allow us to classify them.

The quadratic terms with the scalars are trivial. For any quadratic term $`H_1H_2`$ to be invariant under the gauge interactions, either both of them are singlets (i.e. neutral under all gauge factors) or else $`H_1^{}=H_2`$. Since all the scalars we are discussing are non-singlets, none of the quadratic terms can have any $`BL`$ quantum number, although some of them may contribute to baryon or lepton number violation. Consider for example, $`Y_a^{}Y_a`$. Since $`Y_a`$ can couple to $`\overline{q}\overline{l}`$ as well as to $`qq`$, this term may induce a $`\mathrm{\Delta }B=1`$, $`\mathrm{\Delta }L=1`$ process, but $`BL`$ remains conserved.

The trilinear couplings are the most interesting ones. It was shown  that there are 38 possible trilinear couplings with the scalar bilinears. Although baryon and lepton numbers are not always uniquely specified, all of them carry a $`BL`$ quantum number $`\pm 2`$. As we shall see next, these trilinear operators can mediate four possible processes and hence can be classified in four categories

* $`B=1,L=1`$: In this case both baryon and lepton numbers are changed by 1, so all these couplings should give only proton decay, which violates $`(BL)`$ .
* $`B=0,L=2`$: This gives lepton number violating processes, without affecting the baryon number. So this contributes to the Majorana masses of the neutrinos. Some of them will contribute to the neutrinoless double beta decay , but others would involve neutrinos of different generations.
* $`B=1,L=3`$: These operators give a three-lepton decay mode of the proton . In some models these processes could dominate over the usual proton decay modes.
* $`B=2,L=0`$: In this case baryon number is changed by two unit with no change in lepton number. This corresponds to neutron-antineutron oscillations.

Some of the operators may contribute to two process types. There are no trilinear couplings, which conserve $`BL`$ (which includes of course the case $`B=0`$ and $`L=0`$).

There is only one trilinear coupling which can give proton decay into a lepton and a meson (or mesons) and is not involved in any other process. It is thus possible to assign $`B=1`$ and $`L=1`$ to this operator:

$$\left(\begin{array}{cc}B=1,& L=1\end{array}\right)$$

$$\begin{array}{c}𝒪_1=\mu _1X_a\varphi Y_d\hfill \end{array}$$
(1)
In the following $`BL`$ quantum numbers for the trilinear couplings are not mentioned explicitly since all of them have $`BL=\pm 2`$. All other couplings which contribute to the proton decay, can also mediate $`n\overline{n}`$ oscillations or the neutrinoless double beta decay.

There are only two operators which can give rise to the three-lepton decay mode of the proton:

$$\left(\begin{array}{cc}B=0,& L=2\end{array}\right)\left(\begin{array}{cc}B=1,& L=1\end{array}\right)\mathrm{and}\left(\begin{array}{cc}B=1,& L=3\end{array}\right)$$

$$\begin{array}{cc}𝒪_2=\mu _2Y_aY_c^{}\chi ^+\hfill & 𝒪_3=\mu _3Y_bY_c^{}\xi ^{}\hfill \end{array}$$
(2)
These operators will simultaneously allow the single-lepton decay mode of the proton ($`B=1,L=1`$) and neutrinoless double beta decay ($`B=0,L=2`$). There is no trilinear coupling which contributes to the three-lepton decay mode of the proton and no other process.

There are four trilinear couplings which carry $`B=2,L=0`$ and contribute only to $`n\overline{n}`$ oscillations and not to any other process:

$$\left(\begin{array}{cc}B=2,& L=0\end{array}\right)$$

$$\begin{array}{cc}𝒪_4=\mu _4\mathrm{\Delta }_a\mathrm{\Delta }_b\mathrm{\Delta }_b\hfill & 𝒪_5=\mu _5\mathrm{\Delta }_c\mathrm{\Delta }_a\mathrm{\Delta }_a\hfill \\ 𝒪_6=\mu _6\mathrm{\Delta }_cY_dY_d\hfill & 𝒪_7=\mu _7\mathrm{\Delta }_L\mathrm{\Delta }_L\mathrm{\Delta }_a\hfill \end{array}$$
(3)
Another four trilinear interactions contribute to both $`n\overline{n}`$ oscillations and proton decay:

$$\left(\begin{array}{cc}B=2,& L=0\end{array}\right)\mathrm{and}\left(\begin{array}{cc}B=1,& L=1\end{array}\right)$$

$$\begin{array}{cc}𝒪_8=\mu _8Y_cY_dY_d\hfill & 𝒪_9=\mu _9\mathrm{\Delta }_LY_bY_d\hfill \\ 𝒪_{10}=\mu _{10}\mathrm{\Delta }_bY_aY_d\hfill & 𝒪_{11}=\mu _{11}\mathrm{\Delta }_aY_dY_c\hfill \end{array}$$
(4)
There are several operators which contribute only to neutrinoless double beta decay:

$$\left(\begin{array}{cc}B=0,& L=2\end{array}\right)$$

$$\begin{array}{ccc}𝒪_{12}=\mu _{12}\varphi \varphi \chi ^{}\hfill & 𝒪_{13}=\mu _{13}\varphi \varphi \xi \hfill & 𝒪_{14}=\mu _{14}\chi ^{}\chi ^{}L^{++}\hfill \\ 𝒪_{15}=\mu _{15}\xi \xi L^{++}\hfill & 𝒪_{16}=\mu _{16}X_bX_a^{}\chi ^{}\hfill & 𝒪_{17}=\mu _{17}X_bX_a^{}\xi \hfill \\ 𝒪_{18}=\mu _{18}\mathrm{\Delta }_bX_a^{}X_a^{}\hfill & 𝒪_{19}=\mu _{19}\mathrm{\Delta }_cX_a^{}X_b^{}\hfill & 𝒪_{20}=\mu _{20}\mathrm{\Delta }_LX_a^{}X_a^{}\hfill \\ 𝒪_{21}=\mu _{21}\mathrm{\Delta }_a^{}\mathrm{\Delta }_b\chi ^{}\hfill & 𝒪_{22}=\mu _{22}\mathrm{\Delta }_a^{}\mathrm{\Delta }_L\xi \hfill & 𝒪_{23}=\mu _{23}\mathrm{\Delta }_b^{}\mathrm{\Delta }_c\chi ^{}\hfill \\ 𝒪_{24}=\mu _{24}\mathrm{\Delta }_L^{}\mathrm{\Delta }_c\xi \hfill & 𝒪_{25}=\mu _{25}\mathrm{\Delta }_a^{}\mathrm{\Delta }_cL^{}\hfill & \end{array}$$
(5)
There are few more operators contributing to neutrinoless double beta decay, which also allow proton decay:

$$\left(\begin{array}{cc}B=0,& L=2\end{array}\right)\mathrm{and}\left(\begin{array}{cc}B=1,& L=1\end{array}\right)$$

$$\begin{array}{ccc}𝒪_{26}=\mu _{26}Y_aY_d^{}\chi ^{}\hfill & 𝒪_{27}=\mu _{27}Y_bY_d^{}\xi \hfill & 𝒪_{28}=\mu _{28}Y_cY_d^{}L^{}\hfill \\ 𝒪_{29}=\mu _{29}X_aX_bY_c^{}\hfill & 𝒪_{30}=\mu _{30}X_a\varphi ^{}Y_a\hfill & 𝒪_{31}=\mu _{31}X_a\varphi ^{}Y_b\hfill \\ 𝒪_{32}=\mu _{32}X_aX_aY_a^{}\hfill & 𝒪_{33}=\mu _{33}X_aX_aY_b^{}\hfill & 𝒪_{34}=\mu _{34}X_bY_d\varphi ^{}\hfill \end{array}$$
(6)
The proton decay constraints would then restrict the couplings of these operators making them very much suppressed. However, if some of the couplings could be avoided with discrete symmetries, then these operators could also contribute to the neutrinoless double beta decay significantly. Finally there are four more operators which contribute to neutrinoless double beta decay, proton decay, and $`n\overline{n}`$ oscillations simultaneously:

$$\left(\begin{array}{cc}B=0,& L=2\end{array}\right),\left(\begin{array}{cc}B=1,& L=1\end{array}\right)\mathrm{and}\left(\begin{array}{cc}B=2,& L=0\end{array}\right)$$

$$\begin{array}{cc}𝒪_{35}=\mu _{35}Y_dY_aY_a\hfill & 𝒪_{36}=\mu _{36}Y_dY_bY_b\hfill \\ 𝒪_{37}=\mu _{37}\mathrm{\Delta }_aY_aY_a\hfill & 𝒪_{38}=\mu _{38}\mathrm{\Delta }_aY_bY_b\hfill \end{array}$$
(7)
The quartic terms belong to two broad classes, one with $`BL=0`$ and the other with $`BL=4`$. Again they can be categorized under certain subclasses. Those which are products of the usual quadratic couplings of the form $`\mathrm{\Phi }^{}\mathrm{\Phi }`$ belong to the trivial category with $`B=0`$, $`L=0`$ (and thus $`BL=0`$), which we shall not list here. All the quartic couplings involving only the dileptons ($`\chi ^{},L^{}`$, and $`\xi `$), the usual Higgs doublet $`\varphi `$, and the octet $`\mathrm{\Sigma }`$ also fall into this trivial category with $`B=0`$, $`L=0`$. There are also some other quartic couplings involving the diquarks and leptoquarks, which can be classified under this trivial category. These are

$$BL=0\left(\begin{array}{cc}B=0,& L=0\end{array}\right)$$

$$\begin{array}{ccc}\stackrel{~}{𝒪}_1=\chi ^+\chi ^+\xi \xi \hfill & \stackrel{~}{𝒪}_2=\varphi \varphi \chi ^+L^{}\hfill & \stackrel{~}{𝒪}_3=\mathrm{\Sigma }\mathrm{\Sigma }\chi ^+L^{}\hfill \\ \stackrel{~}{𝒪}_4=\varphi ^{}\varphi ^{}\xi L^{++}\hfill & \stackrel{~}{𝒪}_5=\mathrm{\Sigma }\mathrm{\Sigma }\xi ^{}L^{}\hfill & \stackrel{~}{𝒪}_6=\varphi ^{}\varphi ^{}\mathrm{\Sigma }\mathrm{\Sigma }\hfill \\ \stackrel{~}{𝒪}_7=\mathrm{\Delta }_a^{}X_aX_a\chi ^{}\hfill & \stackrel{~}{𝒪}_8=\mathrm{\Delta }_a^{}X_aX_a\xi \hfill & \stackrel{~}{𝒪}_9=\mathrm{\Delta }_b^{}X_bX_bL^{}\hfill \\ \stackrel{~}{𝒪}_{10}=\mathrm{\Delta }_c^{}X_bX_b\chi ^{}\hfill & \stackrel{~}{𝒪}_{11}=\mathrm{\Delta }_c^{}X_bX_b\xi \hfill & \stackrel{~}{𝒪}_{12}=\mathrm{\Delta }_L^{}X_bX_bL^{}\hfill \\ & \stackrel{~}{𝒪}_{13}=\mathrm{\Delta }_a^{}X_aX_bL^{}\hfill & \end{array}$$
(8)
There are several $`BL=0`$ quartic couplings which mediate baryon number violating processes. The simplest of them have $`B=2`$ and $`L=2`$. The hydrogen-antihydrogen ($`p+e^{}\overline{p}+e^+`$) oscillations and double proton decay into two positrons ($`p+pe^++e^+`$) are typical examples of such processes. There are 7 such interactions

$$BL=0\left(\begin{array}{cc}B=2,& L=2\end{array}\right)$$

$$\begin{array}{ccc}\stackrel{~}{𝒪}_{14}=\mathrm{\Delta }_a\mathrm{\Delta }_c\mathrm{\Delta }_cL^{}\hfill & \stackrel{~}{𝒪}_{15}=\mathrm{\Delta }_b\mathrm{\Delta }_b\mathrm{\Delta }_b\chi ^{}\hfill & \stackrel{~}{𝒪}_{16}=\mathrm{\Delta }_b\mathrm{\Delta }_b\mathrm{\Delta }_L\xi \hfill \\ \stackrel{~}{𝒪}_{17}=\mathrm{\Delta }_b\mathrm{\Delta }_L\mathrm{\Delta }_L\chi ^{}\hfill & \stackrel{~}{𝒪}_{18}=\mathrm{\Delta }_b\mathrm{\Delta }_L\mathrm{\Delta }_L\xi \hfill & \stackrel{~}{𝒪}_{19}=\mathrm{\Delta }_L\mathrm{\Delta }_L\mathrm{\Delta }_L\chi ^{}\hfill \\ & \stackrel{~}{𝒪}_{20}=\mathrm{\Delta }_L\mathrm{\Delta }_L\mathrm{\Delta }_L\xi \hfill & \end{array}$$
(9)
There are other quartic couplings also, which allow these processes with $`B=2`$ and $`L=2`$. They involve $`Y_a`$, $`Y_b`$, and $`Y_c`$, which couple to two quarks as well as to an antiquark and an antilepton. Hence some of these interactions can simultaneously allow $`B=1`$ and $`L=1`$ proton decay processes such as $`pe^+`$ or $`n\overline{\nu }`$ associated with a neutral meson and a lepton-antilepton pair. In all these interactions, the constraints from proton decay will not allow processes with $`B=2`$. These operators are

$$BL=0\left(\begin{array}{cc}B=2,& L=2\end{array}\right)\mathrm{and}\left(\begin{array}{cc}B=1,& L=1\end{array}\right)$$

$$\begin{array}{cc}\stackrel{~}{𝒪}_{21}=\mathrm{\Delta }_cY_aY_d\chi ^{}\hfill & \stackrel{~}{𝒪}_{22}=\mathrm{\Delta }_cY_bY_d\xi \hfill \end{array}$$
(10)
$$BL=0\left(\begin{array}{c}B=2\\ L=2\end{array}\right),\left(\begin{array}{c}B=1\\ L=1\end{array}\right)\mathrm{and}\left(\begin{array}{c}B=0\\ L=0\end{array}\right)$$

$$\begin{array}{ccc}\stackrel{~}{𝒪}_{23}=\mathrm{\Delta }_aY_aY_c\chi ^{}\hfill & \stackrel{~}{𝒪}_{24}=\mathrm{\Delta }_aY_bY_c\xi \hfill & \stackrel{~}{𝒪}_{25}=\mathrm{\Delta }_bY_aY_cL^{}\hfill \\ \stackrel{~}{𝒪}_{26}=\mathrm{\Delta }_bY_aY_a\chi ^{}\hfill & \stackrel{~}{𝒪}_{27}=\mathrm{\Delta }_bY_aY_b\xi \hfill & \stackrel{~}{𝒪}_{28}=\mathrm{\Delta }_bY_bY_b\chi ^{}\hfill \\ \stackrel{~}{𝒪}_{29}=\mathrm{\Delta }_bY_bY_b\xi \hfill & \stackrel{~}{𝒪}_{30}=\mathrm{\Delta }_cY_aY_aL^{}\hfill & \stackrel{~}{𝒪}_{31}=\mathrm{\Delta }_cY_aY_aL^{}\hfill \\ \stackrel{~}{𝒪}_{32}=\mathrm{\Delta }_LY_aY_a\xi \hfill & \stackrel{~}{𝒪}_{33}=\mathrm{\Delta }_LY_aY_b\chi ^{}\hfill & \stackrel{~}{𝒪}_{34}=\mathrm{\Delta }_LY_aY_b\xi \hfill \\ \stackrel{~}{𝒪}_{35}=\mathrm{\Delta }_LY_bY_b\chi ^{}\hfill & \stackrel{~}{𝒪}_{36}=\mathrm{\Delta }_LY_bY_b\xi \hfill & \stackrel{~}{𝒪}_{37}=\mathrm{\Delta }_LY_bY_cL^{}\hfill \\ \stackrel{~}{𝒪}_{38}=Y_aY_cY_d\chi ^{}\hfill & \stackrel{~}{𝒪}_{39}=Y_bY_cY_d\xi \hfill & \stackrel{~}{𝒪}_{40}=Y_cY_dY_dL^{}\hfill \end{array}$$
(11)
$$BL=0\left(\begin{array}{c}B=2\\ L=2\end{array}\right),\left(\begin{array}{c}B=1\\ L=1\end{array}\right),\left(\begin{array}{c}B=0\\ L=0\end{array}\right)\mathrm{and}\left(\begin{array}{c}B=1\\ L=1\end{array}\right)$$

$$\begin{array}{ccc}\stackrel{~}{𝒪}_{41}=Y_aY_aY_a\chi ^{}\hfill & \stackrel{~}{𝒪}_{42}=Y_aY_aY_b\xi \hfill & \stackrel{~}{𝒪}_{43}=Y_aY_aY_cL^{}\hfill \\ \stackrel{~}{𝒪}_{44}=Y_aY_bY_b\chi ^{}\hfill & \stackrel{~}{𝒪}_{45}=Y_aY_bY_b\xi \hfill & \stackrel{~}{𝒪}_{46}=Y_bY_bY_b\chi ^{}\hfill \\ \stackrel{~}{𝒪}_{47}=Y_bY_bY_b\xi \hfill & \stackrel{~}{𝒪}_{48}=Y_bY_bY_cL^{}\hfill & \end{array}$$
(12)
Another class of quartic couplings can give rise to $`BL=4`$ processes. Some of them can produce only $`B=2`$ and $`L=2`$ processes like $`n+ne^{}e^{}\pi ^+\pi ^+`$ or $`n\overline{n}`$ oscillation associated with charged leptons and charged pions. The couplings which can give rise to these processes are

$$BL=4\left(\begin{array}{cc}B=2,& L=2\end{array}\right)$$

$$\begin{array}{cc}\widehat{𝒪}_1=Y_dY_dY_dL^{++}\hfill & \widehat{𝒪}_2=\mathrm{\Delta }_aY_dY_dL^{++}\hfill \\ \widehat{𝒪}_3=\mathrm{\Delta }_a\mathrm{\Delta }_a\mathrm{\Delta }_aL^{++}\hfill & \widehat{𝒪}_4=\mathrm{\Delta }_a\mathrm{\Delta }_a\mathrm{\Delta }_b\chi ^+\hfill \end{array}$$
(13)
The other subclass of $`BL=4`$ processes have $`B=1`$ and $`L=3`$, which give rise to $`p\nu \nu \nu \pi ^+`$ for example. There are two such operators,

$$BL=4\left(\begin{array}{cc}B=1,& L=3\end{array}\right)$$

$$\begin{array}{cc}\widehat{𝒪}_5=X_aX_aX_a\varphi ^{}\hfill & \widehat{𝒪}_6=X_aX_aX_a\mathrm{\Sigma }^{}\hfill \end{array}$$
(14)
There are also some quartic couplings with $`BL=4`$, which allow both of the above two interactions. There are three such terms

$$BL=4\left(\begin{array}{cc}B=1,& L=3\end{array}\right)\mathrm{and}\left(\begin{array}{cc}B=2,& L=2\end{array}\right)$$

$$\begin{array}{ccc}\widehat{𝒪}_7=Y_aY_dY_d\chi ^+\hfill & \widehat{𝒪}_8=\mathrm{\Delta }_aY_aY_d\chi ^+\hfill & \widehat{𝒪}_9=\mathrm{\Delta }_aY_bY_d\xi ^{}\hfill \end{array}$$
(15)
In summary, we have made a general model-independent study of the possible baryon and lepton number violations beyond the standard model with scalar bilinears. Depending on the processes under consideration, they can be classified into only a few categories. We list all the trilinear and quartic couplings and discuss the processes to which these operators contribute.

Acknowledgement

Two of us (EM and US) thank the Max-Planck-Institut für Kernphysik for hospitality. The work of EM is supported in part by the U.S. Department of Energy under Grant No. DE-FG03-94ER40837.


text:
# 1 Introduction

## 1 Introduction

At present there are many possible models of $`\nu `$ masses and mixing . This variety is mostly due to the considerable experimental ambiguities that remain. In particular different solutions for solar neutrino oscillations are still possible. Although the Large Angle (LA) solution emerges as the most likely from present data, the other solutions LOW and Vacuum Oscillations (VO) are still not excluded (the small angle solution is very disfavoured by now and we will disregard it in most of the following discussion). Indeed no solution is actually leading to an imposingly good fit and, for example, the discrimination between LA and LOW is only based on a few hints which are far from compelling . Hopefully in a few months, when the first results from the KamLAND experiment  will be known, one will have decisive evidence on this matter. Here we tentatively disregard the possibility of a third neutrino oscillation frequency as indicated by the LSND experiment  but not confirmed by KARMEN  and to be finally checked by the MiniBOONE experiment  now about to start.

For model building there is an important quantitative difference between the LA solution on the one side and the LOW or VO solutions on the other side. While for all these solutions the relevant mixing angle $`\theta _{12}`$ is large, the value of the squared mass difference $`\mathrm{\Delta }m_{12}^2=m_2^2m_1^2`$ (with, by definition, $`m_2^2m_1^2`$) is very different for LA, LOW and VO: $`10^4eV^2`$, $`10^7eV^2`$ and $`10^{10}eV^2`$, respectively. Thus the gap with respect to the atmospheric neutrino oscillation frequency $`\mathrm{\Delta }m_{23}^2=m_3^2m_2^2`$, which is given by $`|\mathrm{\Delta }m_{23}^2|310^3eV^2`$, is moderate for LA and very pronounced for the other two solutions.

For the LOW and VO solutions the large frequency difference with respect to that of atmospheric neutrinos points to a hierarchical spectrum for the three light neutrinos. Possible hierarchical patterns are the normal hierarchy case $`m_3m_2>m_1`$ or the inverted hierarchy alternative $`m_2>m_1m_3`$ (in this case $`\mathrm{\Delta }m_{23}^2`$ is negative in our definition). Then a main problem is to explain the presence of large mixing angles between largely splitted mass states (in particular the atmospheric neutrino oscillation mixing angle $`\theta _{23}`$ is experimentally close to maximal). In hierarchical models the consistency of these usually opposed constraints is obtained by mechanisms that guarantee a vanishing or a strongly suppressed 23 sub-determinant. In the see-saw mechanism for neutrino masses this suppression can be naturally obtained, for example, in the so-called lopsided models  and/or by the dominance  of one eigenvalue in $`M^1`$, $`M`$ being the right-handed (RH) Majorana matrix. Models of this type have been studied and provide, as also quantitatively confirmed by our present analysis, an essentially unique framework for a successful description of both atmospheric and solar neutrino oscillations when the LOW or the VO solutions are adopted.

In the case of the LA solution the ratio of the solar and atmospheric $`\mathrm{\Delta }m^2`$ ranges is typically given by

$$r=\mathrm{\Delta }m_{12}^2/|\mathrm{\Delta }m_{23}^2|1/201/100.$$
(1)
For LA one can reproduce the data either in a nearly degenerate or in a hierarchical model. In a degenerate model, due to laboratory and cosmological bounds, the common value of $`m_i^2m^2`$ cannot exceed a few $`eV^2`$. But the actual value is probably well below this level because of the constraint imposed by neutrinoless double beta decay ($`0\nu \beta \beta `$)  that would otherwise require a strong cancellation, only possible for nearly maximal solar oscillation mixing . This fact, together with the general difficulty, in the absence of a specific mechanism, of obtaining too small values of $`\mathrm{\Delta }m_{12}^2/m^2`$, suggests that a moderate degeneracy is more likely. Typically we could have all $`m_i^2(\mathrm{few}10^310^2)eV^2`$ with one comparatively not-so-small splitting $`\mathrm{\Delta }m_{12}^210^4eV^2`$. Or, as a different example, we can have a (normal) hierarchical model with $`m_3^2310^3eV^2`$, $`m_2^2<10^4eV^2`$ and $`m_1^20`$ or an (inverse) hierarchical model with $`m_{1,2}^2310^3eV^2`$, and $`m_3^20`$. Actually a sufficient hierarchy (a factor of 5 in mass) can arise from no significant underlying structure at all. In particular, the see-saw mechanism, being quadratic in the Dirac neutrino masses, tends to enhance small fluctuations in the Dirac eigenvalue ratios. This is the point of view of anarchical models , where no structure is assumed to exist in the neutrino mass matrix and the smallness of $`r`$ is interpreted as a fluctuation. But one additional feature of the data plays an important role in this context and presents a clear difficulty for anarchical models. This is the experimental result that the third mixing angle $`\theta _{13}`$ is small, $`|U_{e3}|=|\mathrm{sin}\theta _{13}|<0.2`$ . So, for neutrinos two mixing angles are large and one is small. Instead in anarchical models all angles should apriori be comparable and not particularly small. Therefore this is a difficulty for anarchy and, for the survival of these models, it is necessary that $`\theta _{13}`$ is found very close to the present upper bound. Instead in hierarchical models the smallness of $`\theta _{13}`$ can be obtained as a reflection of the underlying structure in that some small parameter is present from the beginning in these models.

In this note we make a quantitative study of the ability of different models to reproduce the solar neutrino solutions. As a flexible testing ground we consider models based on SU(5)$`\times `$U(1)F. The SU(5) generators act “vertically” inside one generation, while the U(1)F charges are different “horizontally” from one generation to the other. If, for a given interaction vertex, the U(1)F charges do not add to zero, the vertex is forbidden in the symmetric limit. But the symmetry is spontaneously broken by the VEV $`v_f`$ of a number of “flavon” fields with non vanishing charge. Then a forbidden coupling is rescued but is suppressed by powers of the small parameters $`\lambda _f=v_f/\mathrm{\Lambda }`$ with the exponent larger for larger charge mismatch . We expect $`v_f>M_{GUT}`$ and, for the cut-off $`\mathrm{\Lambda }`$ of the theory, $`\mathrm{\Lambda }<M_{Pl}`$. In these models  the known generations of quarks and leptons are contained in triplets $`\mathrm{\Psi }_i^{10}`$ and $`\mathrm{\Psi }_i^{\overline{5}}`$, $`(i=1,2,3)`$ corresponding to the 3 generations, transforming as $`10`$ and $`\overline{5}`$ of SU(5), respectively. Three more SU(5) singlets $`\mathrm{\Psi }_i^1`$ describe the RH neutrinos. In SUSY models we have two Higgs multiplets, which transform as 5 and $`\overline{5}`$ in the minimal model. All mass matrix elements are of the form of a power of a suppression factor times a number of order unity, so that only their order of suppression is defined. We restrict for simplicity to integral charges: this is practically a forced choice for the LA case where the hierarchy parameter must be relatively large (so that $`\sqrt{\lambda }0(1)`$), while for the LOW and VO cases, where the hierarchy parameter is small, it is only motivated by the fact that enough flexibility is obtained for the present indicative purposes. There are many variants of these models : fermion charges can all be non negative with only negatively charged flavons, or there can be fermion charges of different signs with either flavons of both charges or only flavons of one charge. The Higgs charges can be equal, in particular both vanishing or can be different. We will make use of this flexibility in order to study the relative merits of anarchy versus various degrees and different patterns of hierarchy.

In this context we have studied in detail different classes of models: normal hierarchical models with and without automatic suppression of the 23 (sub)determinant, inverse hierarchy models and anarchical models. The normal hierarchical models without automatic suppression of the 23 determinant are clearly intermediate: in a sense in those cases anarchy is limited to the 23 sector. We denote them as partially hierarchical or semi-anarchical in the following. We also compare, when applicable, models with light neutrino masses dominated by the see-saw mechanism or by non renormalizable dim-5 operators. We construct our models by assigning suitable sets of charges for $`\mathrm{\Psi }_i^{10}`$, $`\mathrm{\Psi }_i^{\overline{5}}`$ and $`\mathrm{\Psi }_i^1`$. In all input mass matrices the coefficients multiplying the power of the hierarchy parameter are generated at random as real and complex numbers in a given range of values . We compare the case of real or complex parameters and we also discuss the delicate questions of the probability distribution for the coefficients and the stability of the results. We assign a merit factor to each model given by the percentage of success over a large sample of trials. For each model the value of the hierarchy parameter is adjusted by a coarse fitting procedure to get the best rate of success.

Our results can be summarized as follows. As expected, for the LOW and VO cases only hierarchical model provide a viable approach: in comparison the rate of success of anarchical and seminarchical models is negligible. But also for the LA solution we still find that hierarchical models are sizeably better in general. The most efficient ones are inverse hierarchy models with no see-saw dominance, which are more than 10 times better with respect to anarchy with see-saw (anarchy prefers the see-saw case by about a factor of 2). Among the see-saw dominance versions the most performant models remain the hierarchical ones (by a factor of about 4 with respect to anarchy with see-saw) with not much difference between inverse or normal hierarchy. Semi-Anarchical models are down by a factor of about 2 with respect to hierarchical models among the see-saw versions (but this value is less stable with respect to changes of the extraction procedure and, for example, it tends to be washed out going from complex to real coefficients). In all models the $`\mathrm{tan}^2\theta _{23}`$ distribution is in agreement with large mixing but it is not sharply peaked around 1 as for maximal mixing. Near maximal mixing is instead a prediction for solar neutrinos in inverse hierarchical models, so that their advantage with respect to other models would be rapidly destroyed if the data will eventually move in the direction away from maximality.

## 2 Framework

We consider a class of models with an abelian flavour symmetry compatible with SU(5) grand unification. Here we will not address the well-known problems of grand unified theories, such as the doublet-triplet splitting, the proton lifetime, the gauge coupling unification beyond leading order and the wrong mass relations for charged fermions of the first two generations. We adopt the SU(5)$`\times `$U(1)F framework simply as a convenient testing ground for different neutrino mass scenarios. In all the models that we study the large atmospheric mixing angle is described by assigning equal flavour charge to muon and tau neutrinos and their weak SU(2) partners (all belonging to the $`\overline{5}(l,d^c)`$ representation of SU(5)). Instead, the solar neutrino oscillations can be obtained with different, inequivalent charge assignments and both the LOW111The LOW and VO solutions can be fitted almost equally well by suitable models. Thus here we will focus mainly on the LOW case. and the LA solution can be reproduced.

A first class of models is characterized by all matter fields having flavour charges of one sign, for example all non negative. An important property of models in this class is that the light neutrino mass matrix $`L^Tm_\nu L`$ is independent from the charges of both the 10 $`(q,u^c,e^c)`$ and 1 $`\nu ^c`$ representations , even in the see-saw case when $`m_\nu =m_D^TM_R^1m_D`$. For $`m_\nu `$ entries the powers of the symmetry breaking parameter $`\lambda `$ are dictated by the charges F of the $`\overline{5}`$. Since in this case what really matters are charge differences, rather than absolute values, the equal charges for the second and third generations can be put to zero, without loosing generality:

$$\mathrm{F}(\overline{5})=(b,0,0)b0.$$
(2)
If $`b`$ also vanishes, then the light neutrino mass matrix will be structure-less and we will call anarchical (A) this sub-class of models. In a large sample of anarchical models, generated with random coefficients, the resulting neutrino mass spectrum can exhibit either normal or inverse hierarchy. Anarchical models clearly prefer the LA solution with a moderate separation between atmospheric and solar frequencies. They tend to favour large, not necessarily maximal, mixing angles, including $`U_{e3}`$, which represents a problem. Therefore, in anarchical models, $`U_{e3}`$ is expected to be close to the present experimental bound.

If $`b`$ is positive, then the light neutrino mass matrix will be structure-less only in the (2,3) sub-sector and we will call semi-anarchical (SA) the corresponding models. In this case, the neutrino mass spectrum has normal hierarchy. However, unless the (2,3) sub-determinant is accidentally suppressed, atmospheric and solar oscillation frequencies are expected to be of the same order and, in addition, the preferred solar mixing angle is small. Nevertheless, such a suppression can occur in a fraction of semi-anarchical models generated with random, order one coefficients. The real advantage over the fully anarchical scheme is represented by the suppression in $`U_{e3}`$.

In a second class of models matter fields have both positive and negative flavour charges. In these models, the light neutrino mass matrix will in general depend also on the charges of 10 and 1. A first sub-case arises when only the RH neutrino fields have charges of both signs. It has been shown that it is possible to exploit this feature to reproduce a neutrino mass spectrum with normal hierarchy and a natural gap between atmospheric and solar frequencies. Via the see-saw mechanism the (2,3) sub-determinant vanishes in the flavour symmetric limit. At the same time a large solar mixing angle can be obtained. Clearly this is particularly relevant for the LOW (or VO) solution. It is less clear to which extent the condition of vanishing determinant is needed for the LA solution and one of the purposes of the present paper is precisely to compare these models, which we will call hierarchical (H) with the anarchical and semi-anarchical models, that do not reproduce such a condition.

Finally, we can have fields with charges of both signs in both the 1 and the $`\overline{5}`$. In this context it is possible to reproduce an inverse hierarchical spectrum, with a large (actually, almost maximal) solar mixing angle and a large gap between atmospheric and solar frequencies . Also this sub-class of models, which we call inversely hierarchical (IH), are appropriate to describe both LOW and LA solutions.

The hierarchical and the inversely hierarchical models may come into several varieties depending on the number and the charge of the flavour symmetry breaking (FSB) parameters. Here we will consider both the case (I) of a single, negatively charged flavon, with symmetry breaking parameter $`\lambda `$ or that of two (II) oppositely charged flavons with symmetry breaking parameters $`\lambda `$ and $`\lambda ^{}`$. In case I, it is impossible to compensate negative F charges in the Yukawa couplings and the corresponding entries in the neutrino mass matrices vanish. Eventually these zeroes are filled by small contributions, arising, for instance, from the diagonalization of the charged lepton sector or from the transformations needed to make the kinetic terms canonical. In our analysis we will always include effects coming from the charged lepton sector, whereas we will neglect those coming from non-canonical kinetic terms.

Another important ingredient in our analysis is represented by the see-saw mechanism . Hierarchical models and semi-anarchical models have similar charges in the $`(10,\overline{5})`$ sectors and, in the absence of the see-saw mechanism, they would give rise to similar results. Even when the results are expected to be independent from the charges of the RH neutrinos, as it is the case for the anarchical and semi-anarchical models, the see-saw mechanism can induce some sizeable effect in a statistical analysis. For this reason, for each type of model, but the hierarchical ones (the mechanism for the 23 sub-determinant suppression is in fact based on the see-saw mechanism), we will separately study the case where RH neutrinos are present and the case where they are absent. When RH neutrinos are present, there are two independent contributions to the light neutrino mass matrix. One of them comes via the see-saw mechanism from the exchange of the heavy RH modes. The other one is provided by L-violating dimension five operators arising from physics beyond the cut-off. These contributions have the same transformation properties under the flavour group and, in general, add coherently. In our analysis we will analyze the case where the see-saw contribution is the dominant one (SS). The absence of RH neutrinos describes the opposite case, when the mass matrix is saturated by the non-renormalizable contribution (NOSS).

For each type of model we have selected what we consider to be a typical representative 222We made no real optimization effort to pick up the ‘most’ representative, but rather a model with a high success rate in its class. and we have collected in table 1 the corresponding charges. In the next section we will compare the performances of the following models: $`\mathrm{A}_{\mathrm{SS}}`$, $`\mathrm{A}_{\mathrm{NOSS}}`$, $`\mathrm{SA}_{\mathrm{SS}}`$, $`\mathrm{SA}_{\mathrm{NOSS}}`$, $`\mathrm{H}_{(\mathrm{SS},\mathrm{I})}`$, $`\mathrm{H}_{(\mathrm{SS},\mathrm{II})}`$, $`\mathrm{IH}_{(\mathrm{SS},\mathrm{I})}`$, $`\mathrm{IH}_{(\mathrm{SS},\mathrm{II})}`$, $`\mathrm{IH}_{(\mathrm{NOSS},\mathrm{I})}`$ and $`\mathrm{IH}_{(\mathrm{NOSS},\mathrm{II})}`$.

Anarchical, semi-anarchical and hierarchical models give rise to a mass matrix for light neutrinos of the type

$$m_\nu =\left(\begin{array}{ccc}ϵ^2& ϵ& ϵ\\ ϵ& 1& 1\\ ϵ& 1& 1\end{array}\right)(𝙰,\mathrm{𝚂𝙰},𝙷),$$
(3)
where all the entries are specified up to order one coefficients and the overall mass scale has been conventionally set to one. For anarchical models, $`ϵ=1`$. Then all the entries are uncorrelated numbers of order one and no particular pattern becomes manifest. For the semi-anarchical model of table 1, $`ϵ=\lambda `$. There is a clear distinction between the first row and column and the 23 block of the mass matrix, which is structureless as in the anarchical models. In particular, barring accidental cancellations, the sub-determinant in the 23 sector is of order one. Finally, the hierarchical model defined by the choice of charges in table 1, has $`ϵ=\lambda ^2`$. At variance with the anarchical or semi-anarchical models, the determinant of the 23 sector is suppressed by the see-saw mechanism and is of order $`\lambda \lambda ^{}`$.

The inversely hierarchical models are characterized by a neutrino mass matrix of the kind

$$m_\nu =\left(\begin{array}{ccc}ϵ^2& 1& 1\\ 1& \eta ^2& \eta ^2\\ 1& \eta ^2& \eta ^2\end{array}\right)(\mathrm{𝙸𝙷}),$$
(4)
where $`ϵ=\lambda `$ $`(\lambda ^2)`$ and $`\eta =\lambda ^{}`$ $`(\lambda ^2)`$ for the LA (LOW) solution. The ratio between the solar and atmospheric oscillation frequencies is not directly related to the sub-determinant of the block 23, in this case. The above mass matrices also receive an additional contribution from the diagonalization of the charged lepton sector, which, however, does not spoil the displayed structure. For completeness, we collect in table 2 the gross features of the models under consideration. Notice that the hierarchical models predict a ratio $`\mathrm{\Delta }m_{12}^2/|\mathrm{\Delta }m_{23}^2|`$ of order $`\lambda ^4`$ or $`\lambda ^6`$. In these cases it is possible to fit both the LOW and the LA solutions with an expansion parameter $`\lambda \lambda ^{}`$ that, within a factor of two, matches the Cabibbo angle. On the contrary, the inversely hierarchical model that reproduces the LA solution favours $`\mathrm{\Delta }m_{12}^2/|\mathrm{\Delta }m_{23}^2|O(\lambda ^2)`$. If we adopt this model to fit also the LOW solution, the corresponding values of $`\lambda \lambda ^{}`$ do not provide a decent description of the remaining fermion masses. For this reason, when analyzing the LOW solution, we have considered a separate set of charges for the inverted hierarchy.

## 3 Method and results

Abelian flavour symmetries predict each entry of fermion mass matrices up to an unknown dimensionless coefficient. These coefficients, that are free-parameters of the theory, are expected to have absolute values of order one. Aside from generalized kinetic terms, which we do not consider, the relevant mass matrices are specified by $`N=24`$ order one parameters: 9 from the charged lepton sector, 9 from Dirac neutrino mass matrix entries and 6 from Majorana RH neutrino mass matrices. When RH neutrinos are absent (that is, for the NOSS cases) the LH neutrino mass matrix contains 6 parameters and the relevant set reduces to $`N=15`$ parameters. We have analyzed the case of real or complex coefficients $`𝒫=(𝒫_1,\mathrm{}𝒫_N)`$, with absolute values generated as random numbers in an interval $``$ and random phases taken in $`[0,2\pi ]`$. To study the dependence of our results on $``$, we have considered several possibilities: $`[0.5,2]`$ (default), $`[0.8,1.2]`$, $`[0.95,1.05]`$ and $`[0,1]`$. In the case of real coefficients, which is studied for comparison, we allow both signs for the coefficients. For each model, only a portion $`V_S`$ of the volume $`V=(2\pi )^N`$ of the parameter space gives rise to predictions in agreement with the experimental data within the existing uncertainties. We may interpret $`V_S/V`$ as the success rate of the model in describing neutrino data. Clearly this portion shrinks to zero for infinitely good measurements. Therefore we are not interested in its absolute size, but rather in the relative sizes of $`V_S/V`$ in different models.

We evaluate the success rate of each model by considering, through a random generation, a large number of ‘points’ $`𝒫`$ and by checking whether the corresponding predictions do or do not fall in the experimentally allowed regions . To this purpose we perform a test based on four observable quantities: $`O_1=r\mathrm{\Delta }m_{12}^2/|\mathrm{\Delta }m_{23}^2|`$, $`O_2=\mathrm{tan}^2\theta _{12}`$, $`O_3=|U_{e3}||\mathrm{sin}\theta _{13}|`$ and $`O_4=\mathrm{tan}^2\theta _{23}`$. We take :

$$\begin{array}{c}410^6<r<410^5\hfill \\ |U_{e3}|<0.2\hfill \\ 0.52<\mathrm{tan}^2\theta _{12}<1.17\hfill \\ 0.33<\mathrm{tan}^2\theta _{23}<3.3\hfill \end{array}(\mathrm{LOW})$$
(5)
$$\begin{array}{c}0.01<r<0.2\hfill \\ |U_{e3}|<0.2\hfill \\ 0.24<\mathrm{tan}^2\theta _{12}<0.89\hfill \\ 0.33<\mathrm{tan}^2\theta _{23}<3.3\hfill \end{array}(\mathrm{LA})$$
(6)
The boundaries of these windows are close to the $`3\sigma `$ limits on the corresponding observable quantity. The test is successfully passed if $`O_i(𝒫)`$ are in the above windows, for all $`i=1,\mathrm{},4`$. We will study the dependence on the choice made in eqs. (5,6), by analyzing the distributions of the points generated for each observable. We then estimate $`V_S/V`$ from the ratio between the number of successful trials over the total number of attempts. Notice that we do not extend our test to the remaining fermion masses and mixing angles, though the models considered can also reproduce, at the level of order of magnitudes, charged lepton masses, quark masses and mixing angles. With the interval $``$ fixed at its reference value, $`[0.5,2]`$, the different models are compared at their best performance, after optimizing for each model the value of the symmetry breaking parameters $`\lambda `$ and $`\lambda ^{}`$.

### 3.1 LOW solution

In fig. 1 we compare the success rates for the LOW solution. In this figure the anarchical or semi-anarchical models do not appear simply because their rates of success are negligible on the scale of the figure.

It is clear that, if the future experimental results will indicate the LOW solution as the preferred one, then anarchical or semi-anarchical schemes will be completely inadequate to describe the data. It would be natural in that case to adopt a model where the large gap between the atmospheric and the solar oscillation frequencies is built in as a result of a symmetry. What is striking about the results displayed in fig. 1 is the ability of the IH models, both in the SS and in the NOSS versions, to reproduce the data. An appropriate choice for $`\lambda `$ ($`\lambda ^{}`$) is completely successful in reproducing $`r`$. $`U_{e3}`$ is numerically close to $`r`$ and easily respects the present bound. Moreover, $`\mathrm{tan}^2\theta _{12}`$ is very close to 1, due to the pseudo-Dirac structure of the 12 sector. Only the $`\mathrm{tan}^2\theta _{23}`$ distribution shows some flatness and contributes to deplete the success rate. It is worth stressing that $`\mathrm{tan}^2\theta _{12}`$ is so strongly peaked around 1, that any significant deviation from $`\mathrm{tan}^2\theta _{12}=1`$ in the data would provide a severe difficulty for the IH schemes. The H model has a smaller success rate (still much larger than that of the A and SA models), but has smoother distributions for the four observables and is less sensitive than IH to variations of the experimental data. The error bars in fig. 1 are dominated by the systematic effects, which have been estimated by varying the interval $``$. We considered four possibilities: $`=[0.5,2]`$ (default), $`=[0.8,1.2]`$, $`=[0.95,1.05]`$, and $`=[0,1]`$. To the highest (lowest) rate of each model we then add (subtract) linearly the statistical error. The latter is usually smaller than the systematic error. Further details are reported in appendix A.

### 3.2 LA solution

The success rates for the LA solutions are displayed in fig. 2 and 3, separately for the NOSS and SS cases. The two sets of models have been individually normalized to give a total rate 100. Before normalization the total success rates for NOSS and for SS were in the ratio 1.7:1 (see also appendix A). Although the gaps between the rates of different models are reduced compared to the LOW case, nevertheless a clear pattern emerges from these figures. The present data are most easily described by the IH schemes in their NOSS version. Their performances are better by a factor of 10-30 with respect to the last classified, the anarchical models. The ability of the IH schemes in describing the data can be appreciated from the distributions of the four observables, which, for $`\mathrm{IH}_{(\mathrm{NOSS},\mathrm{II})}`$ and $`\lambda =\lambda ^{}=0.25`$, are displayed in fig. 4.

The observables $`r`$ and $`U_{e3}`$ are strongly correlated. Actually, as discussed in ref.  and shown in table 2, in inversely hierarchical models $`U_{e3}`$ is typically of order $`r`$. Therefore, once $`\lambda `$ and $`\lambda ^{}`$ have been tuned to fit $`r`$, this choice automatically provides a good fit to $`U_{e3}`$. Moreover, similarly to the case of the LOW solution, $`\mathrm{tan}^2\theta _{12}`$ is peaked around 1. At present $`\mathrm{tan}^2\theta _{12}=1`$ is excluded for the LA solution, but, thanks to the width of the distribution, the experimentally allowed window is sufficiently populated. The width of the distribution is almost entirely dominated by the effect coming from the diagonalization of the charged lepton sector. Indeed, by turning off the small parameters $`\lambda `$ and $`\lambda ^{}`$ in the mass matrix for the charged leptons, we get a vanishing success rate for $`\mathrm{IH}_{(\mathrm{NOSS},\mathrm{I})}`$, whereas the rate for $`\mathrm{IH}_{(\mathrm{NOSS},\mathrm{II})}`$ decreases by more than one order of magnitude. It is worth stressing that even a moderate further departure of the window away from $`\mathrm{tan}^2\theta _{12}=1`$ could drastically reduce the success rates of the IH schemes. Finally, the $`\mathrm{tan}^2\theta _{23}`$ distribution is rather flat, with a moderate peak in the currently favoured interval. All the IH models, with or without see-saw, have distributions similar to those shown in fig. 4. In particular the $`\mathrm{tan}^2\theta _{23}`$ distribution of fig. 4 is qualitatively common to all U(1) models. This reflects the fact that the large angle $`\theta _{23}`$ is induced by the equal charges $`\mathrm{F}(\overline{5}_2)=\mathrm{F}(\overline{5}_3)`$, a feature shared by all the models we have considered. To appreciate the relevant correlations, the distributions for the model $`\mathrm{IH}_{(\mathrm{NOSS},\mathrm{II})}`$ are also displayed in fig. 5, as scatter plots in the planes $`(U_{e3},r)`$, $`(\mathrm{tan}^2\theta _{12},r)`$ and $`(\mathrm{tan}^2\theta _{23},r)`$.

Without see-saw mechanism, the next successful model is the semi-anarchical model SA, whose distributions are displayed in fig. 6 (in our notation there are no $`\mathrm{H}_{\mathrm{NOSS},\mathrm{I}(\mathrm{II})}`$ models). Compared to the IH case, the $`\mathrm{tan}^2\theta _{12}`$ distribution has no pronounced peak. Possible shifts in the central value of $`\mathrm{tan}^2\theta _{12}`$ would not drastically modify the results for the SA model. The $`U_{e3}`$ distribution is peaked around $`\lambda =0.2`$ with tails that exceed the present experimental bound. The $`r`$ distribution is centered near $`r=1`$. Finally, the anarchical scheme in its NOSS version is particularly disfavoured, due to its tendency to predict $`r`$ close to 1 and also due to $`U_{e3}`$, that presents a broad distribution with a preferred value of about 0.5.

The overall picture changes significantly if the LA solution is realized in the context of the see-saw mechanism, as illustrated in fig. 3. The IH models are still rather successful. Compared to the NOSS case, the $`\mathrm{IH}_{\mathrm{SS}}`$ models slightly prefer higher values of $`r`$ and, due to a smaller $`\lambda =\lambda ^{}`$, they have very narrow $`\mathrm{tan}^2\theta _{12}`$ distributions. This can be seen from the scatter plots of fig. 7. The other distributions are similar to those of figs. 4 and 5. We observe that while most of the points in fig. 7 are centered around $`U_{e3}O(\lambda ^2)`$, there is also a small region clustered at $`U_{e3}0.6`$.

Equally good or even better results are obtained by the $`\mathrm{H}_{(\mathrm{SS},\mathrm{II})}`$ model, with distributions shown in fig. 8 and 9. We see that, at variance with the IH models, $`\mathrm{tan}^2\theta _{12}`$ is not spiky, which results in a better stability of the model against variation of the experimental results. The preferred value of $`r`$ is close to the lower end of the experimental window. The $`U_{e3}`$ distribution is nicely peaked around $`\lambda ^2`$.

The $`\mathrm{SA}_{\mathrm{SS}}`$ model is significantly outdistanced from $`\mathrm{H}_{(\mathrm{SS},\mathrm{II})}`$, $`\mathrm{IH}_{(\mathrm{SS},\mathrm{II})}`$ and $`\mathrm{IH}_{(\mathrm{SS},\mathrm{I})}`$. It is particularly penalized by the $`U_{e3}`$ distribution, centered around $`\lambda =0.3`$. Finally, the least favoured models are $`\mathrm{H}_{(\mathrm{SS},\mathrm{I})}`$ and $`\mathrm{A}_{\mathrm{SS}}`$. The model $`\mathrm{H}_{(\mathrm{SS},\mathrm{I})}`$ fails both in $`U_{e3}`$ (see fig. 10), which tends to be too large for the preferred value of $`\lambda =\lambda ^{}=0.5`$ and in $`\mathrm{tan}^2\theta _{12}`$. The $`\mathrm{A}_{\mathrm{SS}}`$ model, as its NOSS version, suffers especially from the $`U_{e3}`$ distribution (see fig. 10) which is roughly centered at 0.5, with only few percent of the attempts falling within the present experimental bound. A large $`U_{e3}`$ can be regarded as a specific prediction of anarchy and any possible improvement of the bound on $`|U_{e3}|`$ will wear away the already limited success rate of the model. The distributions in $`\mathrm{tan}^2\theta _{12}`$ and $`\mathrm{tan}^2\theta _{23}`$ are equally broad and peaked around 1. Compared to the NOSS case, $`\mathrm{A}_{\mathrm{SS}}`$ has a better $`r`$ distribution, well located inside the allowed window.

As another criterion for evaluating the quality of a model, we address the issue of the stability of the observables $`O_i`$ with respect to small fluctuations of the set of coefficients $`𝒫`$. Notice that the random coefficients to be put in front of the powers of $`\lambda ,\lambda ^{}`$ stand for the combined result of a fundamental theory of flavour, present at a certain scale $`\mathrm{\Lambda }`$, and of an evolution from $`\mathrm{\Lambda }`$ down to $`m_Z`$. It would thus be natural if the physical observables $`O_i`$ were stable under small perturbations of the coefficients, $`\mathrm{\Delta }𝒫`$, around a given successful representative set $`𝒫_0`$.

This is illustrated for the seesaw case in fig. 11, where we compare the models $`\mathrm{IH}_{(\mathrm{SS},\mathrm{II})}`$, $`\mathrm{H}_{(\mathrm{SS},\mathrm{II})}`$ and $`\mathrm{A}_{\mathrm{SS}}`$. The blue dot refers to the observables $`O_i`$ which follow from a typical successful configuration for LA, $`𝒫_0`$. The 40 points in red, light blue and green are the result of adding to $`𝒫_0`$ random perturbations $`\mathrm{\Delta }𝒫`$ with $`|\mathrm{\Delta }𝒫|=|𝒫|/10`$ and random phases. The yellow dots correspond to $`|\mathrm{\Delta }𝒫|=|𝒫|/2`$. As appears from the scatter plots, $`\mathrm{IH}_{(\mathrm{SS},\mathrm{II})}`$ is very stable: as already argued in the above discussion, actually it is even dangerously stable with respect to the prediction for $`\mathrm{tan}^2\theta _{12}`$. If models with hierarchy display a sufficient degree of stability, anarchical ones are much less stable, in particular with respect to the predictions for $`U_{e3}`$ and $`\mathrm{tan}^2\theta _{12}`$. As shown in fig. 12, the non see-saw case enhances these features: $`\mathrm{IH}_{(\mathrm{NOSS},\mathrm{II})}`$ is extremely stable, while $`\mathrm{A}_{\mathrm{NOSS}}`$ is higly unstable. For the latter, in particular, even with the small fluctuation considered in the figures, $`\mathrm{tan}^2\theta _{12}`$ spans from .1 to 10. Thus, the criterion of stability supports the same ratings of models previously obtained by only considering total rates and distributions.

As a general comment we observe that our results are rather stable with respect to the choice of the interval $``$. With only one exception, namely the crossing between $`H_{(SS,I)}`$ and $`A_{SS}`$, the relative position of the different models, according to their ability in describing the data, does not vary when we shrink $``$, or when we extend it to cover the full circle of radius 1 333The results for $`=[0,n]`$ are independent on $`n`$, since changing $`n`$ amounts to perform a renormalization of all mass matrices by a common overall scale, which is not felt by the observables we have used in our analysis.. This stability would be partially upset by restricting to the case of real coefficients $`𝒫`$. In that case the relative rates are comparable to those obtained in the complex case only for sufficiently wide intervals $``$, typically $`[1/\sqrt{\lambda },\sqrt{\lambda }][\sqrt{\lambda },1/\sqrt{\lambda }]`$. If we further squeeze $``$ around $`\pm 1`$ the rates of all models tend to zero. In the real case, the distribution of $`\mathrm{tan}^2\theta _{23}`$ is very sensitive to the width of the extraction interval $``$. Indeed, in the flavour symmetry basis, charged leptons and light neutrinos mass matrices are both diagonalised by a 23 mixing angle which tends to $`\pi /4`$ when $``$ is squeezed to $`1`$. As an effect, in this limit the distribution of $`\theta _{23}`$ is almost empty around $`\pi /4`$ and presents two peaks at $`0`$ and $`\pi /2`$. On the contrary, in the complex case the distribution of $`\theta _{23}`$ is quite insensitive to the width of $``$ thanks to the smearing effect of the phases.

## 4 Conclusion

If a large gap between the solar and atmospheric frequencies for neutrino oscillations were finally to be established by experiment then this fact would immediately suggest that neutrino masses are hierarchical, similar to quark and charged lepton masses. In fact very small mass squared differences among nearly degenerate neutrino states are difficult to obtain and make stable under running in the absence of an ad hoc symmetry. Moreover, the presence of very small mixing angles would also indicate a hierarchical pattern. At present, the SA solution of solar neutrino oscillations is disfavoured so that 2 out of 3 mixing angles appear to be large and only one appears to be small, although the present limit is not terribly constraining. However, in the case of the LA solution for solar neutrino oscillations, the ratio $`r`$ of the solar to atmospheric frequencies is not so small and it has been suggested that possibly this solution does not require a symmetry to be generated. In the anarchical framework the smallness of $`U_{e3}`$, assumed not too pronounced, and that of $`r`$, are considered as accidental (the see-saw mechanism helps in this respect because the product of 3 matrices sizeably broadens the $`r`$ distribution). In the present study we have examined in quantitative terms the relative merits of anarchy and of different implementations of hierarchy in reproducing the observed features of the LA solution. For our analysis we have adopted the framework of SU(5)$`\times `$U(1)F which is flexible enough, by suitable choices of the flavour charges, to reproduce all interesting types of hierarchy and also of anarchy. This framework allows a statistical comparison of the different schemes under, as far as possible, homogeneous conditions. The rating of models in terms of their statistical rates of success is clearly a questionable procedure. After all Nature does not proceed at random and a particular mass pattern that looks odd could arise due to some deep dynamical reason. However, the basis for anarchy as a possible description of the LA solution can only be formulated in statistical terms. Therefore it is interesting to compare anarchy versus hierarchy on the same grounds. We have considered models both with normal and with inverse hierarchies, with and without see-saw, with one or two flavons, and compared them with the case of anarchy and of semi-anarchy (models where there is no structure in 2-3 but only in 2-3 vs 1). The stability of our results has been tested by considering different options for the statistical procedure and also by studying the effect on each type of model of small parameter changes.

Our conclusion is that, for the LA solution, the range of $`r`$ and the small upper limit on $`U_{e3}`$ are sufficiently constraining to make anarchy neatly disfavoured with respect to models with built-in hierarchy. If only neutrinos are considered, one might counterargue that hierarchical models have at least one more parameter than anarchy, in our case the parameter $`\lambda `$. However, if one looks at quarks and leptons together, as in the GUT models that we consider, then the same parameter that plays the role of an order parameter for the CKM matrix, for example, the Cabibbo angle, can be successfully used to reproduce also the LA hierarchy. On the one hand, it is interesting that the amount of hierarchy needed for the LA solution is just a small power of the Cabibbo angle. On the other hand, if all fermion masses and mixings are to be reproduced, a similar parameter is also needed in anarchical models, where all the quark and lepton mass structure arises from the charges of $`\mathrm{\Psi }_i^{10}`$. In comparison to anarchy even the limited amount of structure present in semi-anarchical models already improves the performance a lot. And the advantage is further increased when more structure is added as in inverse hierarchy models, or models with normal hierarchy and automatic suppression of the 23 determinant. In the see-saw case all these types of hierarchical models have comparable rates of success (except for $`\mathrm{H}_{(\mathrm{SS},\mathrm{I})}`$). In the non-see-saw versions of inverse hierarchy the performance is even better. An experimental criterion that could eventually decide between normal and inverse hierarchy models is the closeness of the solar angle $`\theta _{12}`$ to its maximal value. If the data moves away from $`\pi /4`$ the probability of inverse hierarchy will rapidly drop in comparison to hierarchical models.

## Acknowledgments

We would like to thank Francesco Vissani for useful discussions and for partecipating in the early stages of the present work. F.F. and I.M. thank the CERN TH divison for its hospitality during july and august 2002, when this project was developed. This work is partially supported by the European Programs HPRN-CT-2000-00148 and HPRN-CT-2000-00149.

## Appendix: raw data

In this section we list the results of our numerical analysis, for the case of complex random coefficients $`𝒫`$. In particular we detail the dependence of the success rates on the size of the window $``$ that specifies the absolute value of the coefficients $`𝒫`$. In most cases the ‘systematic’ error due to the choice of $``$ is larger than the statistical error. The latter is given by $`\delta P/P=1/\sqrt{N_s}`$, where $`P=100N_s/N`$ is the success rate and $`N_s`$ is the number of successes in $`N`$ trials. We chose $`N=10000`$, $`30000`$, $`50000`$ respectively, for the LOW solution, for LA without see-saw and for LA with see-saw. When fitting the LOW solution, the results for the anarchical and semi-anarchical models do not depend on $`\lambda `$ or $`\lambda ^{}`$.

## 5 Addendum

Since the completion of our analysis new important experimental results have been published, in particular the first data from KamLAND  and the new results from SNO . The main new information, of special relevance for the present analysis, is that the LOW solution is now discarded and the window corresponding to the LA solution has considerably narrowed. In particular the value of $`r`$ has settled around $`r1/35`$ and the solar angle $`\theta _{12}`$ is now more than $`5\sigma `$ away from maximal. For example, in ref.  from a comprehensive analysis of all the available data the following $`3\sigma `$ window was obtained :

$$\begin{array}{c}0.018<r<0.053\hfill \\ |U_{e3}|<0.23\hfill \\ 0.30<\mathrm{tan}^2\theta _{12}<0.64\hfill \\ 0.45<\mathrm{tan}^2\theta _{23}<2.57\hfill \end{array}$$
(7)
These values are to be compared with those in eq. (6), although in eq. (7) we have a $`3\sigma `$ window according to ref.  while the indicative ranges in eq. (6) were less precisely defined. These new experimental developments have a large enough impact on the present approach to make it certainly worthwhile and timely to update the results of our analysis taking the new information into account. As discussed in our paper, the deviation of the solar mixing angle from the maximal value is expected to produce a pronounced decrease of the success rate of the IH models. Also the considerable narrowing of the allowed $`r`$ range (the upper part of the LA solution is now disfavoured) should have important consequences on the relative effectiveness of A, SA, H and IH models. It is interesting to precisely quantify these issues and this addendum is devoted to this task.

Given the new values in eq. (7), the values of $`\lambda `$ and $`\lambda ^{}`$ have been reoptimized and in some cases their values are somewhat changed with respect to the previous numbers, as shown in Tables 7 and 8, which replace Tables 5 and 6. The values of the charges for the representation 10 have also been modified for those cases where a rather large value of $`\lambda `$ follows from the optimization (in order to maintain the correct amount of quark and lepton mass hierarchies), as shown in Table 9 that replaces Table 1. Since the neutrino mixing parameters are completely independent on the 10 charges, this change is only important for a better fit to quark and charged lepton masses and mixings.

The new versions of figs. 2 and 3 that describe the success rates separately for the NOSS and SS cases are presented in figs. 13 and 14. The main qualitative difference between the new and the old rates is that indeed, both in the NOSS and the SS cases, a large drop in the IH relative success rate is immediately apparent. This is due not only to the smaller upper limit on the solar mixing angle $`\theta _{12}`$ but also to the smaller value of $`r`$ (in fact, we recall from Table 2 that $`r`$ and $`\theta _{13}`$ are expected to be of the same order in IH models). From the updated histograms in figs. 13 and 14 we see that normal hierarchy models (with two oppositely charged flavons $`\mathrm{H}_{\mathrm{II}}`$) are neatly preferred over anarchy and inverse hierarchy in the context of these SU(5)$`\times `$U(1) models. In particular, in the SS case, the $`\mathrm{H}_{\mathrm{II}}`$ models with normal hierarchy and suppressed 23 sub-determinant are clearly preferred. We recall that for the chosen charge values the $`\mathrm{H}_{\mathrm{II}}`$ model is of the lopsided type.

In conclusion, with all the limitations of the present approach, it is interesting that the hierarchical models are preferred over anarchy and inverse hierarchy. In particular in the case of see-saw dominance the preferred models are those with natural suppression of the 23 subdeterminant (e.g those of the lopsided type or those with dominance of a light right-handed eigenvalue) which indeed provide a simple natural solution of all known constraints.


text:
# Expanded version of the contribution in “Field Theory and Statistical Mechanics, Rome 10-15 June 2002” Particle oscillations in external chaotic fields

## 1 INTRODUCTION

We review here the general formalism for the study of neutrino propagation in the presence of stochastic media . This formalism will allow the systematic derivation of evolution equations for averaged quantities as survival probabilities and higher order moments of them. New integrodifferential equations valid for finite correlation processes are obtained for the first time. For exponentially correlated processes a second order ordinary equation is obtained as a consequence. The Redfield equation valid for Gaussian delta-correlated noise is rederived in a simple way as a particular case: it can be obtained directly or as the zero-order term of an asymptotic expansion in the inverse of the correlation length. In the context of neutrino oscillations, it is shown that the presence of matter noise induces the appearance of a effective complex part in the matter density. Finally the formalism, together with the quantum correlation theorem, will be applied to the problem of the computation of distribution higher moments and non-equal-time correlation functions.

The formalism may be applied to any quantum system governed by similar Schrödinger equations: equations where a stochastic function appears multiplicatively in some part of the Hamiltonian. For simplicity of notation, in this work we will always deal with Hamiltonians which are finite dimensional operators in some Hilbert space, the equations will be obviously of the same type for general systems described by infinite-dimensional Hamiltonians.

## 2 THE GENERAL CASE

Let us consider a general system whose evolution is described by the following linear random Schrödinger equation:

$`i_tX`$ $`=`$ $`\rho (t)L(t)X,X(0)=X_0,`$ (1)
where $`X`$ is a vector (i.e., a multiflavor wave function) or a matrix (i.e., a density matrix) of arbitrary dimension. Equation (1) is defined in an interaction representation, any additive non-stochastic term has been solved for and absorbed in the definitions of $`X(t)`$ and $`L(t)`$. $`L(t)`$ is a general linear operator. It can be an ordinary matrix if $`X`$ represents a wave function. If $`X`$ is a density matrix, $`L(t)`$ is a commutator:

$$L_A(t)X[A(t),X].$$
(2)
Any other linear operator is admissible: obviously it is always possible to define an enlarged vector space where the action of $`L`$ is represented by a matrix.

We will assume that $`\rho (t)`$ in Eq.(1) is a scalar Gaussian process completely determined by its first two moments, which, without loss of generality can be taken as

$`\rho (t)=0,`$ $`\rho (t)\rho (t^{})=f(t,t^{}).`$ (3)
In other terms, $`\rho (t)`$ is characterized completely by the measure

$`[d\rho ]\mathrm{exp}{\displaystyle \frac{1}{2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\rho (t)f^1(t,t^{})\rho (t^{})𝑑t𝑑t^{}.`$ (4)
An important particular case is when the process is $`\delta `$-correlated, the correlation function is then of the form

$`f(t,t^{})=\mathrm{\Omega }^2\delta \left(tt^{}\right).`$ (5)
A convenient way of parametrizing a correlation function with a finite correlation length is to use an exponential function:

$`f(t,t^{})=\mathrm{\Omega }^2ϵ\mathrm{exp}\left(ϵtt^{}\right).`$ (6)
In this case the correlation length is defined as $`\tau =1/ϵ`$. The expression (5) is reobtained letting $`\tau 0`$ or $`ϵ\mathrm{}`$ in Eq.(6).

$`X(t)`$, the solution to Eq.(1), is a stochastic function. The objective of this work is to obtain equations for its ensemble average and higher moments in a systematic way . Let us remark that for us Eq.(1) is purely phenomenological, we assume that it is the result of a more complete microscopic analysis which can account for the randomness of $`\rho (t)`$ (see for example, Refs.).

To obtain a differential equation for the average of $`X`$ we make use of the following well known property: For any Gaussian process $`\rho (t)`$ characterized by a $`\delta `$-correlation function (6) and any functional $`F[\rho ]`$ we have the following relation :

$$F[\rho ]\rho (t)=\mathrm{\Omega }^2\frac{\delta F[\rho ]}{\delta \rho (t)}.$$
(7)
For a Gaussian process with an arbitrary correlation function \[Eq. (4)\] we have instead the general relation:

$$F[\rho ]\rho (t)=𝑑\tau \rho (t)\rho (\tau )\frac{\delta F[\rho ]}{\delta \rho (\tau )}.$$
(8)
Let us consider now the evolution operator $`U`$ for a particular realization of Eq.(1). By definition,

$$X(t)=U(t,t_0)X_0.$$
(9)
The operator $`U(t,t_0)`$ is a functional of $`\rho `$, it has the following formal expression in terms of a time ordered exponential:

$$U(t,t_0)=T\mathrm{exp}\left[i_{t0}^t\rho (\tau )L(\tau )𝑑\tau \right].$$
(10)
The functional derivative of $`U`$ with respect to $`\rho (t)`$ can be computed by direct methods. By differentiating term by term the series expansion for Eq.(10), the result is

$$\frac{\delta U(t,t_0)}{\delta \rho (\tau )}=iL(\tau )U(\tau ,t_0),t<\tau <t_0.$$
(11)
In order to obtain a differential equation for $`X`$ we observe that

$$i_tX=i_tX=L(t)\rho (t)X.$$
(12)
Combining together Eqs. (8),(9),(11) we easily obtain the result that the evolution equation for the ensemble average is in the general case an integrodifferential equation given by

$`i_tX(t)`$ $`=`$ $`i{\displaystyle _0^t}𝑑t^{}f(t,t^{})L(t)L(t^{})X(t^{}),X(0)=X_0.`$ (13)
This equation, which is exact and of very general validity, is the equation which we were looking for and one of the main results of the present work. Note that previously integrodifferential equations have been obtained, valid for particular cases or in particular limits, using heuristic ad hoc arguments (for example, in Ref. ). The derivation of Eq.(13) which has been done here is the rigorous justification for such approaches.

In some notable cases Eq.(13) can be reduced to an ordinary differential equation. In the next sections we will see how it can be reduced to, respectively, first and second order equations for the correlation functions given by Eqs. (5) and (6).

## 3 THE PARTICULAR $`\delta `$-CORRELATED CASE

For the particular case where the correlation function is of exponential type, a second order ordinary differential equation can be obtained as we will see below. On the other hand, for the simpler $`\delta `$-correlated case the evolution equation (13) becomes the ordinary differential equation

$`i_tX(t)`$ $`=`$ $`i\mathrm{\Omega }^2L^2(t)X(t).`$ (14)
Taking $`L`$ as a commutator, this last equation coincides with the Redfield equation derived by Ref. . Note that the effective ”Hamiltonian” appearing in the second part of Eq. (14) is not Hermitic anymore (an example of a fluctuation-dissipation effect).

In cases of interest for the neutrino oscillation problem, the original equation for the density matrix is of the slightly simpler form

$`i_tX`$ $`=`$ $`[H_0(t)+\rho (t)g(t)H_1,X],X(0)=X_0,`$ (15)
where we have written the commutator explicitly. $`\rho (t)`$ is a stochastic function as before and $`g(t)`$ an arbitrary scalar function. $`H_0,H_1`$ are Hamiltonian matrices, the former contains the average part of $`\rho (t)`$: $`H_0H_0^{}(t)+\rho _0(t)H_1`$, with $`\rho _0=\rho `$. The latter is assumed to be time independent.

For the problem described by Eq. (15) the corresponding Redfield equation for the averaged density matrix is of the form

$`i_tX`$ $`=`$ $`[H_0(t),X]i\mathrm{\Omega }^2g^2(t)[H_1,[H_1,X]],`$ (16)
$``$ $`H_0^{}XXH_0^++2i\mathrm{\Omega }^2g^2(t)H_1XH_1,`$ (17)
where in the last line the following effective Hamiltonians were defined:

$$H_0^\pm =H_0\pm ig^2(t)H_1^2.$$

The solution of what is called the “coherent” part of Eq. (17) (two first terms of the Hamiltonian) is accomplished by defining the average evolution operator:

$`U^\pm `$ $`=`$ $`T\mathrm{exp}{\displaystyle 𝑑\tau H_0^\pm (\tau )},U^{}=U^+^{}.`$ (18)
The coherent part of the density matrix is then

$`X_{coh}`$ $`=`$ $`U^{}X_0U^+^{}.`$ (19)
Defining a new “coherent” interaction representation by the relations

$`H_L`$ $`=`$ $`U^{}^1H_1U^{},H_R=U^{}H_1U^{}^1,X_I=U^{}X_{coh}U^{}^1,`$ (20)
we arrive at the equation we were looking for, the resolution of the original equation is equivalent to the resolution of the following one:

$`i_tX_I`$ $`=`$ $`2i\mathrm{\Omega }^2g^2(t)H_LX_IH_R.`$ (21)
There are some important particular cases where Eq. (14) can be solved or considerably simplified by taking into account the algebraic properties of a specific $`L`$ (in what follows $`k(t)`$ is always a scalar function).

(a) Let us assume that $`L`$ is such that $`L^2(t)=k(t)L(t)`$. This case appears in the computation of the average wave function with matter density noise. Equation (14) reduces to

$$i_tX=i\mathrm{\Omega }^2k(t)L(t)X.$$

The averaged equation is similar to the original one, the nonrandom part of the density is “renormalized” acquiring an imaginary term

$$\rho \rho _0i\mathrm{\Omega }^2k.$$

This is the density which will appear in the coherent effective Hamiltonians $`H_0^\pm `$.

(b) The case where $`L^2(t)=k(t)I`$ with $`I`$ the identity matrix appears in the computation of the averaged neutrino wave function under noisy magnetic spin-flavor precession. The resulting equation can trivially be integrated (to be compared with the previous case):

$$X(t)=\mathrm{exp}\left[\mathrm{\Omega }^2_0^t𝑑\tau k(\tau )\right]X_0.$$

The average wave oscillation is damped by a factor equivalent to the one first calculated by Nicolaidis . A similar damping also appears when computing the average density matrix from Eq. (21). From these differences of behavior with respect to case (a) it is expected that the presence of magnetic field noise can have some influence even if applied far from any resonance region.

(c) The case where $`L^4(t)=k(t)L^2(t)`$ appears in the computation of the average density matrix with matter or magnetic noise. We can obtain in this case the “conservation law”

$$\left[1\mathrm{\Omega }^2k(t)\right]L^2(t)_tX(t)=0.$$

$`L^2(t)`$ is not invertible because the operator $`L(t)`$ has a zero eigenvalue. The previous expression has proved to be of practical importance in some concrete numerical applications .

The presence of zero eigenvalues distinguishes cases (a) and (b) from (c); it can have consequences in the long term behavior of the respective ensemble averages as is shown elsewhere .

## 4 AN ASYMPTOTIC EXPANSION FOR EXPONENTIALLY CORRELATED SYSTEMS

We will see now how Eq. (14) can be obtained as a limiting case when the correlation length tends to zero. For this purpose we use an exponential correlation function as Eq. (6), the integrodifferential evolution equation becomes in this case

$`i_tX(t)`$ $`=`$ $`i\mathrm{\Omega }^2ϵ\mathrm{exp}(ϵt){\displaystyle _0^t}𝑑t^{}\mathrm{exp}(ϵ\tau )L(t)L(t^{})X(t^{}).`$ (22)
Let us compute the asymptotic expansion of the second term of Eq. (22) valid for large $`ϵ`$; the following expansion is valid for any function $`g(t)`$:

$`h(ϵ)`$ $``$ $`ϵ\mathrm{exp}(ϵt){\displaystyle _0^t}𝑑\tau \mathrm{exp}(ϵ\tau )g(\tau )g(t){\displaystyle \frac{g^{}(t)}{ϵ}}+{\displaystyle \frac{g^{\prime \prime }(t)}{ϵ^2}}+\mathrm{}`$ (23)
Inserting this expression into Eq. (22), we obtain the following expansion in powers of $`ϵ`$:

$`i_tX`$ $`=`$ $`i\mathrm{\Omega }^2L^2(t)X+i{\displaystyle \frac{\mathrm{\Omega }^2}{ϵ}}L(t)_t\left(L(t)X\right)+o\left({\displaystyle \frac{1}{ϵ^2}}\right).`$ (24)
To leading order in $`1/ϵ`$, we recover the expression corresponding to the $`\delta `$-correlated case. At next-to-leading order we get finite-correlation correction terms

$`i_tX`$ $`=`$ $`i\mathrm{\Omega }^2L^2(t)X+i{\displaystyle \frac{\mathrm{\Omega }^2}{ϵ}}L(t)L^{}(t)X+i{\displaystyle \frac{\mathrm{\Omega }^2}{ϵ}}L^2(t)_tX`$ (25)
or, equivalently,

$`\left[1{\displaystyle \frac{\mathrm{\Omega }^2}{ϵ}}L^2(t)\right]_tX=\left[\mathrm{\Omega }^2L^2(t)+{\displaystyle \frac{\mathrm{\Omega }^2}{ϵ}}L(t)L^{}(t)\right]X.`$ (26)
Finally, We get the following differential equation valid to order $`1/ϵ`$, making the assumption that the operator which multiplies the left term is invertible:

$`_tX`$ $`=`$ $`\left(\mathrm{\Omega }^2L^2(t)+{\displaystyle \frac{\mathrm{\Omega }^2}{ϵ}}L(t)L^{}(t)+{\displaystyle \frac{\mathrm{\Omega }^4}{ϵ}}L^4(t)\right)X.`$ (27)
This equation can be used for finite, but relatively large, correlation lengths. We see that, up to this degree of approximation, not only the ratio level of noise to correlation length ($`\mathrm{\Omega }^2/ϵ`$) is important. We have different regimes according to the value of $`\mathrm{\Omega }^2`$. The first term will be more important for a low noise amplitude ($`\mathrm{\Omega }^2<<1`$). For strong noise ($`\mathrm{\Omega }^2>>1`$) the second term, proportional in this case to $`L^4`$, will dominate.

## 5 AN EXACT DIFFERENTIAL EQUATION FOR EXPONENTIALLY CORRELATED SYSTEMS

In contrast with the approximate approach used in the previous section, we can actually derive a simple, ordinary second order differential equation for the case of exponential correlation. Let us assume that the original equation is of the same decomposable type as the one appearing in Eq. (15) but let us include other cases using the general notation

$`i_tX`$ $`=`$ $`\left(L_0(t)+\rho (t)g(t)L_1\right)X,X(0)=X_0,`$ (28)
with $`L_0,L_1`$ general linear operators as before, the latter time independent. In this case Eq. (13) is of the form

$`i_tX=L_0(t)Xi\mathrm{\Omega }^2ϵL_1^2e^{ϵt}g(t){\displaystyle _0^t}𝑑t^{}g(t^{})e^{ϵt^{}}X(t^{}).`$ (29)
Differentiating the equation once and performing some simple algebra we obtain the following ordinary second order differential equation:

$`\left[_t\lambda (t)\right]\left[i_tL_0(t)\right]X`$ $`=`$ $`i\mathrm{\Omega }^2ϵg^2(t)L_1^2X,`$
$`i(_tX)_0`$ $`=`$ $`L(0)X_0,X_0=X_0,`$ (30)
with

$$\lambda (t)ϵ+g^{}(t)/g(t).$$

Let us remark that this equation is exact and probably the most important result of this work from a practical point of view.

## 6 HIGHER ORDER MOMENTS AND THE QUANTUM REGRESSION THEOREM

Second order distribution moments, expressions of the type $`X_iX_j`$, or in general, moments of any order, can also be computed using equations similar to Eq. (13). The straightforward procedure is to define products $`X_{ij\mathrm{}k}=X_iX_j\mathrm{}X_k`$ and write differential equations for them using the constitutive equations for each of the $`X_i`$. The resulting equations are of the same type as Eq. (1). The similarity is obvious when one adopts a tensorial notation and defines products of the form $`X^{(n)}=X\mathrm{}X.`$ The task of obtaining evolution equations is especially simple within this notation.

Correlators of quantities at different times also appear in the computation of averages of quantities of physical interest, i.e., expected signal rates. In the most simple case expressions of the type $`X_i(t+\tau )X_j(t)`$. We will shortly show that the knowledge of equal-time moments and the application of the quantum correlation theorem is sufficient for the computation of this kind of correlators.

Our objective now is the computation of the correlation function appearing in two steps. First we define the generalized density matrix $`X^{(2)}`$ as the tensorial product of usual density matrices at two different times and energies:

$`X^{(2)}(E_1,E_2;t_1,t_2)`$ $``$ $`X(E_1,t_1)X(E_2,t_2).`$ (31)
The average of the element $`X_{1111}^{(2)}=X_{11}X_{11}`$ is evidently the probability correlation function we are looking for.

The differential equation for the equal time function $`X^{(2)}(E_1,E_2;t,t)`$ is obtained from the individual evolution equations for the matrices $`X_{1,2}X(E_{1,2})`$ \[indices (1,2) label expressions where $`E_1,E_2`$ appear respectively\]

$`_tX^{(2)}`$ $`=`$ $`H_1X^{(2)}+X^{(2)}H_2\rho LX^{(2)}.`$ (32)
Equation (32) is a random linear differential equation, linear in the stochastic variable $`\rho (t)`$. Applying the formalism developed in the previous section, we can immediately write the equation for the ensemble average $`X^{(2)}`$ \[Eqs. (1314)\]. Once we know the equal time correlator, we obtain the expression for any other pair $`(t,t^{})`$ using the quantum regression theorem  which reads as follows. For any vector Markov process $`Y`$, if the ensemble average of $`Y`$ fulfills a Schrödinger-like equation of the type

$`_tY(t)=G(t)Y(t),`$ (33)
with $`G(t)`$ an arbitrary matrix, then the second order correlations will obey the following equation:

$`_\tau Y_i(t+\tau )Y_l(t)={\displaystyle \underset{j}{}}G_{ij}(\tau )Y_j(t+\tau )Y_l(t).`$ (34)
Note that the quantum regression theorem is in principle not applicable to systems described by the integral equation Eq. (13). Nevertheless it is applicable to problems where the integral equation can be reduced to an ordinary differential equation such as Eq. (30). Such a second order equation can easily be expressed as a first order one by defining the auxiliary pair process $`Y=(X,_tX)`$.

For the simplest case where $`\rho `$ is a $`\delta `$ correlated process the equation for non-equal time correlators is explicitly given by

$`i_\tau X^{(2)}(t+\tau ,t)`$ $`=`$ $`i\mathrm{\Omega }^2L^2(\tau )X^{(2)}(t+\tau ,t).`$ (35)
The initial condition, $`X^{(2)}(t,t)`$, is the equal time correlator previously obtained. For exponentially correlated problems second order equations similar to Eq. (30) can immediately be obtained.

## 7 CONCLUSIONS AND FINAL REMARKS

In conclusion, in this work we have developed the general formalism for the study of neutrino propagation in stochastic media. This formalism has allowed the systematic derivation of evolution equations for averaged quantities as survival probabilities and higher order moments of them. New integrodifferential equations valid for finite correlation processes have been obtained for the first time. For exponentially correlated processes a second order ordinary equation is obtained as a consequence. The Redfield equation valid for Gaussian $`\delta `$-correlated noise is rederived in a simple way: it has been obtained directly and as the zero-order term of an asymptotic expansion in the inverse of the correlation length.

The formalism can be generalized in an obvious way to obtain the ensemble average of equations of a slightly more general type than Eq. (1), equations of the form

$`i_tX`$ $`=`$ $`{\displaystyle \underset{i}{}}\rho _i(t)L_i(t)X,X(0)=X_0.`$ (36)
These equations appear, for example, in the chaotic neutrino magnetic precession. The result for the case where the $`\rho _i`$ are $`\delta `$-correlated in time but mutually uncorrelated is simply

$`i_tX(t)`$ $`=`$ $`i{\displaystyle \underset{i}{}}\mathrm{\Omega }_i^2L_i^2(t)X(t).`$ (37)
The generalization to continuous Schrödinger equations in the presence of random potentials or random external forces is also obvious if we leave aside the mathematical differences coming from the appearance of infinite-dimensional Hilbert spaces. The formalism is of general application to any quantum system governed by similar Schrödinger equations: equations where a stochastic function appears multiplicatively in some term of the Hamiltonian.