text:
# Ruling out four-neutrino oscillation interpretations of the LSND anomaly?

## 1 Introduction

The atmospheric neutrino data atm-exp , skatm , macro , including the most recent 1489 Super-K data sample provide strong evidence for $`\nu _\mu `$ oscillations into an active neutrino (mainly $`\nu _\tau `$) with $`\mathrm{\Delta }m^2atm2\times 10^3\text{eV}^2`$ solat02 . On the other hand, apart from confirming, once again, the long-standing solar neutrino problem sksol , chlorine , sage , gallex\_gno , sno01 , the recent results from the Sudbury Neutrino Observatory (SNO) sno02  have given strong evidence that solar neutrinos convert mainly to an active neutrino flavor. This suggests that an extension of the Standard Model of particle physics is necessary in the lepton sector, capable of incorporating the required $`\nu _e`$ conversion. Although certainly not yet unique Miranda:2000bi , Guzzo:2001mi , the most popular explanation of the solar neutrino data is provided by the active large mixing angle (LMA) neutrino oscillation hypothesis, characterized by a neutrino mass-squared difference $`\mathrm{\Delta }m^2sol10^4\text{eV}^2`$ solat02 .

In contrast, reactor and accelerator neutrino data KARMEN , bugey , CHOOZ , PaloV , CDHS  lead to negative results. However, the LSND experiment LSND , LSND2001  has provided positive results, which may or may not be confirmed by the forth-coming MiniBooNE experiment MiniBooNE . The required neutrino mass-squared difference $`\mathrm{\Delta }m^2lsnd0.2`$ eV2 is in conflict with the neutrino mass-squared differences indicated by solar and atmospheric data in a minimal three-neutrino picture. Four-neutrino models ptv92 , pv93 , cm93 , 4nuModels , 4nuextra , Ioannisian:2001mu , Hirsch:2000xe  potentially account for all current oscillation data. The status of four-neutrino descriptions has been presented in Ref. Maltoni:2001bc . An exhaustive list of four-neutrino references can be found in Ref. giuntiwebp .

The six possible four-neutrino mass spectra are shown in Fig. 1. For the case $`\mathrm{\Delta }m^2lsnd\mathrm{\Delta }m^2atm,\mathrm{\Delta }m^2sol`$, which we tacitly assume in order to reconcile the LSND anomaly with solar and atmospheric data, these schemes can naturally be divided into two very different classes, usually called (3+1) and (2+2) mass schemes barger00 . It is important to note that (3+1) mass spectra include the three-active neutrino scenario as limiting case. In this case solar and atmospheric neutrino oscillations are explained by active neutrino oscillations, with mass-squared differences $`\mathrm{\Delta }m^2sol`$ and $`\mathrm{\Delta }m^2atm`$, and the fourth neutrino state gets completely decoupled. We will refer to such limiting scenario as (3+0). In contrast, the (2+2) spectrum is intrinsically different, as there must be a significant contribution of the sterile neutrino either in solar or in atmospheric neutrino oscillations or in both.

In recent studies Maltoni:2001bc  it has been realized that there is considerably tension in four-neutrino fits to the global data. In such global four-neutrino analyses one is faced with the problem that there are different data sets, which all give very good fits if analyzed separately. Problems arise due to the combination of the different data sets in a four-neutrino framework. In this letter we re-evaluate the status of four-neutrino interpretations of the LSND anomaly in the light of the recent solar sksol , sage , sno02  as well as atmospheric skatm  neutrino results. To evaluate the quality of the fit we will apply appropriate statistical methods, which are especially sensitive to the combination of different data sets in a global analysis. We find that both the SNO NC result as well as the 1489-day Super-K atmospheric neutrino data strongly reject against sterile neutrino conversions. This essentially rules out (2+2) mass schemes. We also re-evaluate the status of (3+1) schemes by considering two different analyses of LSND data LSND2001 , Church:2002tc . We will further elaborate the result of previous studies barger00 , 3+1early , BGGS , peres , carlo , GS , Maltoni:2001mt  that in (3+1) schemes the LSND signal is in serious disagreement with bounds from short-baseline (SBL) experiments reporting no evidence for oscillations KARMEN , bugey , CDHS . The net result is that all four-neutrino descriptions of the LSND anomaly, both in (2+2) as well as (3+1) realizations, are strongly disfavoured by the data. This brings the LSND anomaly to a theoretically more puzzling status. We note that also cosmology put strong constraints on four-neutrino schemes (for recent analyses see Ref. cosmology ).

The plan of the paper is as follows. In Sec. 2 we briefly describe our parameterization and approximations used in the global four-neutrino analysis. In Secs. 3 and 4 we summarize the solar and atmospheric neutrino data and their analysis solat02 . In Sec. 5 we describe the SBL data we are using, and we compare the two different analyses of the LSND data LSND2001 , Church:2002tc . In Sec. 6 we show that (2+2) oscillation schemes are ruled out because of the tension between solar and atmospheric data, whereas in Sec. 7 we present our re-analysis of the disagreement in SBL data in (3+1) oscillation schemes in light of the two LSND samples. In Sec. 8 we show that the goodness of fit is very bad in all four-neutrino cases. Furthermore, we compare the relative quality of the fit for the schemes (3+1) and (2+2), as well as the three active neutrino case (3+0). The quantitative statistical criteria to analyze the conflict between different data sets, which are used in Secs. 6, 7 and 8, are formulated in the appendix A. In summary, we find that all four-neutrino descriptions of the LSND anomaly both in (2+2) as well as (3+1) realizations are highly disfavoured due to recent data, as mentioned in our conclusions, Sec. 9.

## 2 Four-neutrino oscillations

In this letter we interpret current neutrino oscillation data in a four-neutrino framework. This requires the existence of a fourth light neutrino, which must be sterile to account for the LEP invisible Z width data ptv92 , pv93 , cm93 , 4nuModels , giuntiwebp . The simplest theoretical models add just one such light $`SU(2)U(1)`$ singlet neutrino and ascribe its lightness either to an underlying protecting symmetry, such as lepton number ptv92 , pv93 , or to the existence of extra dimensions 4nuextra , Ioannisian:2001mu . Symmetry breaking is required in order to generate solar and atmospheric oscillations, characterized by splittings which arise due to small tree–level Hirsch:2000xe  or radiative effects ptv92 , pv93 . Many other four-neutrino models have been considered 4nuModels , giuntiwebp .

Quite generally, the four-neutrino charged current leptonic weak interaction is specified as by a rectangular $`3\times 4`$ lepton mixing matrix $`K=\mathrm{\Omega }U`$, with $`\mathrm{\Omega }`$ diagonalizing the $`3\times 3`$ charged lepton mass matrix and $`U`$ diagonalizing the $`4\times 4`$ Majorana neutrino mass matrix. Altogether the elements of the mixing matrix $`K`$ are characterized by 6 mixing angles and 3 physical phases Schechter:1980gr , which could lead to CP violation in oscillation phenomena Schechter:1981gk . For simplicity we neglect CP-violating phases, whose effects are expected to be small in the experiments we consider, due to the strong hierarchy of the mass-squared differences required by the experimental data. This leaves us with nine parameters altogether relevant for the description of CP conserving four-neutrino oscillations: 6 mixing angles contained in the mixing matrix and 3 mass-squared differences.

As shown in Ref. Maltoni:2001bc  it is very important to adopt a convenient choice for these parameters in order to factorize the analysis. Here we adopt all the notation and conventions introduced in Ref. Maltoni:2001bc , which are based on the use of physically relevant quantities. The 6 parameters $`\mathrm{\Delta }m^2sol`$, $`\theta sol`$, $`\mathrm{\Delta }m^2atm`$, $`\theta atm`$, $`\mathrm{\Delta }m^2lsnd`$, $`\theta lsnd`$ are similar to the 2-neutrino mass-squared differences and mixing angles and are directly related to the oscillations in solar, atmospheric and the LSND experiments. For the remaining 3 parameters we use $`\eta _s,\eta _e`$ and $`d_\mu `$. Here, $`\eta _s(\eta _e)`$ is the fraction of $`\nu _s(\nu _e)`$ participating in solar oscillations, and $`1d_\mu `$ is the fraction of $`\nu _\mu `$ participating in atmospheric oscillations (for exact definitions see Ref. Maltoni:2001bc ). For the analysis we adopt the following approximations:

* We make use of the hierarchy $`\mathrm{\Delta }m^2sol\mathrm{\Delta }m^2atm\mathrm{\Delta }m^2lsnd`$. This means that for each data set we consider only one mass-squared difference, the other two are set either to zero or to infinity.
* In the analyses of solar and atmospheric data (but not for SBL data) we set $`\eta _e=1`$, which is justified because of strong constraints from reactor experiments bugey , CHOOZ , PaloV .

Due to these approximations (well justified by reactor data and in oder to fit in LSND with solar + atmospheric data) the parameter structure of the problem gets rather simple Maltoni:2001bc . The parameter dependence of the different data sets solar, atmospheric and SBL is illustrated in Fig. 2. We see that only $`\eta _s`$ links solar and atmospheric data and $`d_\mu `$ links atmospheric and SBL data. All the other parameters are “private” to one data set.

## 3  Solar neutrinos

Here we borrow from Ref. solat02  a recent re-analysis of present solar data including a study of degree of rejection against the sterile neutrino oscillation hypothesis. The experimental data used are the solar neutrino rate of the chlorine experiment Homestake chlorine , the rates of the gallium experiments SAGE sage , GALLEX and GNO gallex\_gno , as well as the 1496-day Super-Kamiokande data sample sksol  in form of 44 recoil electron zenith-energy bins. Moreover, we include the data from the recent SNO charged current and neutral current event measurements, as well as day-night information in the form of 34 data points of the observed energy spectrum sno02 . For the solar neutrino fluxes we use the Standard Solar Model fluxes ssm , including its prediction of the 8B flux. The first hints in favour of the LMA solution Gonzalez-Garcia:1999aj , which followed mainly from the flatness of the Super-K spectrum, have now become a robust result, thanks to the new SNO data and additional Super-K spectrum data. In contrast with the pre-SNO-NC situation, presently all non-LMA solutions are strongly rejected, so that in this paper we will restrict the solar oscillation parameters to the LMA region.

Another important consequence of the recent SNO-NC data is specially relevant for us, namely, that the degree of rejection against solar neutrino oscillations to sterile neutrinos has now become much more significant. This results can be summarized in Fig. 3. The line labeled “after” shows $`\mathrm{\Delta }\chi ^2sol`$ as a function of $`0\eta _s1`$, after recent solar neutrino up-dates, most important SNO-NC sno02 . We minimize with respect to the other two solar neutrino oscillation parameters $`\mathrm{\Delta }m^2sol`$ and $`\theta sol`$, required only to lie in the LMA region. For comparison we show also the line labeled “before”, corresponding to the data used in Ref. Maltoni:2001bc . From this figure one can place a mass scheme-independent limit on the parameter $`\eta _s`$ from solar data:

$$\text{solar data (LMA):}\eta _s0.45\text{(99\% C.L.)},$$
(1)
compared to $`\eta _s0.52`$ obtained in Ref. Maltoni:2001bc . Moreover, the pure sterile solution ($`\eta _s=1`$) and large values of $`\eta _s`$ are now much more disfavoured than previously.

## 4 Atmospheric neutrinos

For the atmospheric data analysis solat02  we use the following data from the Super-Kamiokande experiment skatm : $`e`$-like and $`\mu `$-like samples of sub- and multi-GeV data, as well as the up-going stopping and through-going muon data. The Super-K data correspond to the 1489-day sample presented at Neutrino 2002 skatm . In contrast to previous analyses Maltoni:2001bc , GMPV , Concha2+2  here we adopt a more precise ten-bin presentation of the corresponding contained zenith-angle distributions. Further, we use the data from MACRO macro , including the recent update of their up-going muon sample (10 angular bins). Since $`\mathrm{\Delta }m^2sol\mathrm{\Delta }m^2atm`$ we adopt the standard approximation of neglecting $`\mathrm{\Delta }m^2sol`$ when performing the analysis of the atmospheric data111For detailed studies of this issue see Ref. Gonzalez-Garcia:2002mu .. For further details of the atmospheric neutrino analysis see Refs. GMPV , Concha2+2  and references therein.

Let us consider the impact of atmospheric data alone on the parameter $`\eta _s`$. As shown in Ref. Maltoni:2001bc  the behaviour of the atmospheric data as a function of $`\eta _s`$ depends on the mass schemes considered, with three inequivalent cases222Note that $`\eta _s`$ is defined as the fraction of sterile neutrinos in solar oscillations and is only indirectly related to the fraction of sterile neutrinos in atmospheric oscillations. This relation is different for the various mass schemes Maltoni:2001bc . (3+1)a,d, (3+1)b,c and (2+2)A,B. The qualitative behaviour of (3+1)a,d and (3+1)b,c schemes (see Fig. 1) is similar, however, small differences appear due to matter effects in atmospheric oscillations. In the following we will always minimize between these two cases and refer simply to the (3+1) scheme.

We present in Fig. 4 the $`\mathrm{\Delta }\chi ^2atm`$ profile displayed as a function of $`0\eta _s1`$ minimized with respect to all other parameters, before and after the latest data of Ref. skatm  have been included. We see that the new data significantly improve the sensitivity of atmospheric data with respect to a sterile component. Especially for (2+2) schemes small values of $`\eta _s`$ are highly disfavoured, which will be crucial for the exclusion of (2+2) spectra. From the figure we obtain the following bounds on $`\eta _s`$ at 99% C.L.:

$$\text{atmospheric data:}\{\begin{array}{cc}\eta _s0.32\hfill & \text{for (3+1) schemes},\hfill \\ \eta _s0.65\hfill & \text{for (2+2) schemes}.\hfill \end{array}$$
(2)
## 5 Reactor and accelerator data

We divide the data from reactor and accelerator experiments in data from experiments reporting no evidence (NEV) for oscillations and in the LSND experiment. In this paper we compare the implications of two different analyses of the LSND data. First, we use the likelihood function obtained in the final LSND analysis LSND2001  from their global data with an energy range of $`20<E_e<200`$ MeV and no constraint on the likelihood ratio $`R_\gamma `$ (see Ref. LSND2001  for details). This sample contains 5697 events including decay-at-rest (DAR) $`\overline{\nu }_\mu \overline{\nu }_e`$, and decay-in-flight (DIF) $`\nu _\mu \nu _e`$ data. We refer to this analysis as LSND global. Second, we use the LSND analysis performed in Ref. Church:2002tc  based on 1032 events obtained from the energy range $`20<E_e<60`$ MeV and applying a cut of $`R_\gamma >10^5`$. These cuts eliminate most of the DIF events from the sample, leaving mainly the DAR data, which are more sensitive to the oscillation signal. We refer to this analysis as LSND DAR.

In both cases the collaboration has provided us with the likelihood function obtained in their analyses LSNDpriv , eitel  which we convert into a $`\chi ^2`$ according to $`\chi ^22\mathrm{ln}`$ (see Ref. Maltoni:2001bc  for details). In Fig. 5 we compare the 99% C.L. regions obtained from the two LSND analyses. The LSND DAR data prefers somewhat larger mixing angles, which will lead to a stronger disagreement of the data in (3+1) oscillation schemes. Furthermore, the differences in $`\chi ^2`$ between the best fit point and no oscillations for the two analyses are given by $`\mathrm{\Delta }\chi _{\mathrm{no}\mathrm{osc}}^2=29`$ (global) and $`\mathrm{\Delta }\chi _{\mathrm{no}\mathrm{osc}}^2=47`$ (DAR). This shows that the information leading to the positive oscillation signal seems to be more condensed in the DAR data.

Concerning the NEV data, we take into account the detailed information from the SBL disappearance experiments Bugey bugey  and CDHS CDHS , as described in Ref. GS . The constraints from KARMEN KARMEN  are included by means of the likelihood function provided by the collaboration eitel . Further, we include the results on the $`\overline{\nu }_e`$ survival probability from the long-baseline experiments CHOOZ CHOOZ  ($`P=1.01\pm 0.028\pm 0.027`$) and Palo Verde PaloV  ($`P=1.01\pm 0.024\pm 0.053`$) as explained in Ref. Maltoni:2001bc .

## 6 (2+2) oscillations: ruled out by solar and atmospheric data

The strong preference of oscillations into active neutrinos in solar and atmospheric oscillations leads to a direct conflict in (2+2) oscillation schemes Maltoni:2001bc , peres , Concha2+2 . This is evident from the comparison of Eqs. (1) and (2). We will now show that thanks to the new SNO solar neutrino data sno02  and the improved SK statistic on atmospheric neutrinos skatm  the tension in the data has become so strong that (2+2) oscillation schemes are, in contrast to Ref. Maltoni:2001bc , essentially ruled out.

In Fig. 6 we show the $`\chi ^2`$ as a function of $`\eta _s`$ for solar data and for atmospheric combined with SBL data333For completeness we use here (atmospheric+NEV+LSND global) data. However, the impact of SBL data on the $`\chi ^2`$-dependence on $`\eta _s`$ is very small, as can be seen by comparing the corresponding lines in Figs. 6 and 4. Also the differences between LSND global and LSND DAR are negligible in this case.. Furthermore, we show the $`\chi ^2`$ of the global data defined by

$$\overline{\chi }^2(\eta _s)\mathrm{\Delta }\chi ^2sol(\eta _s)+\mathrm{\Delta }\chi _{\mathrm{atm}+\mathrm{sbl}}^2(\eta _s).$$
(3)
In this section we investigate quantitatively the degree of disagreement of the two data sets solar and atmospheric+SBL in (2+2) oscillation schemes using the statistical methods described in detail in the appendix.

First, in the method which we are calling parameter consistency test (PC) we calculate the C.L. at which still some value of the common parameter $`\eta _s`$ is allowed by both data sets. From Fig. 6 we find

$$\chi _{\mathrm{PC}}^2=12.2.$$
(4)
This implies that only if we take both data sets at the 99.95% C.L. a value of $`\eta _s`$ exists, which is contained in the allowed regions of both sets.

Second, the parameter goodness of fit (PG) makes use of the $`\overline{\chi }^2`$ defined in Eq. (3). This criterion evaluates the g.o.f. of the combination of data sets, without being diluted by a large number of data points, as it happens for the usual g.o.f. criterion. For the case of two data sets constraining one single common parameter this method is equivalent to subtracting the best fit values for the parameter obtained from the two data sets and calculating the C.L. at which the difference is consistent with zero. According to Eq. (12) we find

$$\chi _{\mathrm{PG}}^2\overline{\chi }_{\mathrm{min}}^2=23.5.$$
(5)
Evaluated for 1 d.o.f.444Applying Eq. (12) we find for the number of d.o.f. $`3sol+7_{\mathrm{atm}+\mathrm{sbl}}9=1`$ corresponding to the single common parameter $`\eta _s`$. such a high $`\chi ^2`$-value leads to the marginal parameter g.o.f. of $`1.3\times 10^6`$.

From Eqs. (4) and (5) we conclude that (2+2) oscillation schemes are highly disfavoured, if not ruled out: solar data and atmospheric+SBL data are in disagreement at the 3.5$`\sigma `$ level according to the PC test and at the 4.8$`\sigma `$ level according to the PG. This is a very robust result, independent of whether LSND is confirmed or disproved. We note, however, that in the unlikely event that the KamLAND reactor experiment kamland  will see no signal and rule out the LMA solution, the fit of (2+2) schemes may improve. In that case, although the solar fit quality will substantially deteriorate since we are left with non-LMA solar solutions, the agreement with atmospheric data will improve, especially for the quasi-vacuum case solat02 .

## 7 (3+1) oscillations: strongly disfavoured by SBL data

It is known for a long time that (3+1) mass schemes are disfavoured by the comparison of SBL disappearance data with the LSND result 3+1early , BGGS . Here we extend the first quantitative analyses of Refs. GS , Maltoni:2001mt  by evaluating whether the two LSND data sets (global LSND2001  and DAR Church:2002tc ) are consistent with all the rest of the global sample of oscillation data555Here the NEV data plays a key role, but also atmospheric data is important in the lower regions of $`\mathrm{\Delta }m^2lsnd`$ BGGS , Maltoni:2001mt . For completeness we include also solar data, although the impact on the LSND amplitude is very small. (SOL+ATM+NEV) by using the statistical methods explained in detail in the appendix.

Following previous works GS , Maltoni:2001mt  and the original analyses of Bugey bugey  and CDHS CDHS  we perform the analysis for fixed values of $`\mathrm{\Delta }m^2lsnd`$. Therefore, the two data sets we are comparing depend on $`1_{\mathrm{LSND}}`$ and $`8_{\mathrm{SOL}+\mathrm{ATM}+\mathrm{NEV}}`$ parameters and have only the single parameter $`\theta lsnd`$ in common. For each value of $`\mathrm{\Delta }m^2lsnd`$ the analysis is completely analogous to the one of the previous section. In Fig. 7 we show the results of the PC test and the PG test for both LSND data sets as a function of $`\mathrm{\Delta }m^2lsnd`$. The horizontal lines indicate the C.L. of 95.45%, 99.73% and 99.9937%, corresponding to 2, 3 and 4 standard deviations, respectively. From this figure we observe that the disagreement is roughly one standard deviation worse with the LSND DAR analysis compared to the LSND global data analysis. The larger values of $`\mathrm{sin}^22\theta lsnd`$ preferred by the DAR data are in stronger disagreement with the bounds of the NEV experiments.

The upper panel of Fig. 7 shows that for the global LSND data a value of $`\theta lsnd`$ consistent with the rest of the oscillation data exists at best at the 2$`\sigma `$ level. For LSND DAR parameter consistency is possible only at a C.L. greater than 99% (for most $`\mathrm{\Delta }m^2lsnd`$ values, only at more than 3$`\sigma `$). The parameter g.o.f. is shown in the lower panel. For the global LSND data it is always worse than 3$`\sigma `$, except for $`\mathrm{\Delta }m^2lsnd6`$ eV2 it reaches 1%. For LSND DAR it is for most values of $`\mathrm{\Delta }m^2lsnd`$ below 0.01%. The weakest disagreement in all cases occurs for the “LSND island” around $`\mathrm{\Delta }m^2lsnd6`$ eV2. Note that this does not correspond to the best fit value in a global fit (see next section). For the LSND global data the best fit point is at $`\mathrm{\Delta }m^2lsnd=27`$ eV2 with a PC at 97% C.L. and a PG 0.2%. For LSND DAR the best fit point is at $`\mathrm{\Delta }m^2lsnd=0.91`$ eV2 with a PC at 3$`\sigma `$ and an PG of $`3\times 10^5`$.

## 8 Comparing (3+1), (2+2) and (3+0) hypotheses

With the methods developed in Ref. Maltoni:2001bc  we are able to perform a global fit to the oscillation data in the four-neutrino framework. This approach allows to statistically compare the different hypotheses. Let us first evaluate the g.o.f. of (3+1) and (2+2) spectra with the help of the PG method described in the appendix. We divide the global oscillation data in the 4 data sets SOL, ATM, NEV and LSND. According to Eq. (12) we consider the $`\chi ^2`$-function

$`\overline{\chi }^2`$ $`=`$ $`\mathrm{\Delta }\chi ^2sol(\theta sol,\mathrm{\Delta }m^2sol,\eta _s)+\mathrm{\Delta }\chi ^2atm(\theta atm,\mathrm{\Delta }m^2atm,\eta _s,d_\mu )`$ (6)
$`+`$ $`\mathrm{\Delta }\chi ^2nev(\theta lsnd,\mathrm{\Delta }m^2lsnd,d_\mu ,\eta _e)+\mathrm{\Delta }\chi ^2lsnd(\theta lsnd,\mathrm{\Delta }m^2lsnd),`$
and the relevant number of d.o.f. is given by

$$3_{\mathrm{SOL}}+4_{\mathrm{ATM}}+4_{\mathrm{NEV}}+2_{\mathrm{LSND}}9_{\mathrm{parameters}}=4.$$
(7)
This corresponds to the 4 parameters $`\eta _s,d_\mu ,\theta lsnd,\mathrm{\Delta }m^2lsnd`$ describing the coupling of the different data sets.

The results of this analysis is shown in Tab. 1. We show the contributions of the 4 data sets to $`\chi _{\mathrm{PG}}^2\overline{\chi }_{\mathrm{min}}^2`$ for (3+1) and (2+2) oscillation schemes and for both LSND analyses. As expected we observe that in (3+1) schemes the main contribution comes from SBL data due to the tension between LSND and NEV data in these schemes. For (2+2) oscillation schemes a large part of $`\chi _{\mathrm{PG}}^2`$ comes from solar and atmospheric data, however, also SBL data contributes significantly. This comes mainly from the tension between LSND and KARMEN Church:2002tc , which does not depend on the mass scheme and, hence, also contributes in the case of (2+2). Therefore, the values of $`\chi _{\mathrm{PG}}^2`$ given in Tab. 1 for (2+2) schemes are higher than the one of Eq. (5), where the tension in SBL data is not included.

In the last row of Tab. 1 we show the PG for the various cases. The best g.o.f. is obtained with the global LSND data in the (3+1) scheme. However, even in this best case the g.o.f. is only 0.56%. It is fair to say that all other cases are essentially ruled out by the PG criterion. In general the fit is worse with the LSND DAR analysis. Due to the larger values of $`\mathrm{sin}^22\theta lsnd`$ preferred by this data the disagreement between LSND and NEV experiments is significantly stronger than with the global data. We observe that the fit gets worse, if one selects LSND data with higher sensitivity to the oscillation signal. This seems rather discouraging for the four-neutrino hypothesis as an interpretation of LSND.

Although we have seen that none of the four-neutrino mass schemes can provide a reasonable good fit to the global oscillation data including LSND, it might be interesting to consider the relative status of the three hypotheses (3+1), (2+2) and the three-active neutrino scenario (3+0). This can be done by comparing the $`\chi ^2`$ values of the best fit point (which is in the (3+1) scheme) with the one corresponding to (2+2) and (3+0). In Tab. 2 we show the $`\mathrm{\Delta }\chi ^2`$ of (2+2) and (3+0) with respect to (3+1) for LSND global and DAR data.

First we observe that for LSND global (2+2) is strongly disfavoured with respect to (3+1) with a $`\mathrm{\Delta }\chi ^2=17.8`$. The reason for the big change with respect to the value of $`\mathrm{\Delta }\chi ^2=3.7`$ found in Ref. Maltoni:2001bc  is the improved sensitivity of solar (SNO NC) and atmospheric (SK 1489-days) data to a sterile component. With this new data now (3+1) schemes are clearly preferred over (2+2): for 4 d.o.f. a $`\mathrm{\Delta }\chi ^2=17.8`$ implies that (2+2) is ruled out at 99.87% C.L.. Further we see from Tab. 2 that now (2+2) is only slightly better than (3+0) for LSND global. Using LSND DAR data (2+2) is somewhat less disfavoured with respect to (3+1). The reason for this is that for LSND DAR (3+1) gets much worse compared to LSND global, whereas the fit of (2+2) deteriorates only slightly (see Tab. 1). Hence, the difference becomes smaller. Still, (3+1) is clearly preferred: for 4 d.o.f. a $`\mathrm{\Delta }\chi ^2=11.5`$ implies that (2+2) is ruled out at 97.9% C.L. with respect to (3+1). Of course, for LSND DAR (3+0) is very strongly disfavoured with a $`\mathrm{\Delta }\chi ^2=28.5`$ because of the high significance of the oscillation signal in this data set, which cannot be accounted for in the (3+0) case.

## 9 Conclusions

Prompted by the recent SNO neutral current result and recent atmospheric data, we have re-analyzed the global status of all current neutrino oscillation data in the framework of four-neutrino schemes, in a way similar to that of Ref. Maltoni:2001bc . Our present update includes the global solar neutrino data, the most recent sample of atmospheric neutrino experiments, and a detailed treatment of experiments reporting no evidence (NEV) for oscillations (KARMEN, Bugey, CDHS as well as information from CHOOZ and Palo Verde). Besides the analysis of the global LSND data, we also perform an analysis based mainly on the decay-at-rest (DAR) data sample, which provides a higher sensitivity to the oscillation signal in LSND. We have identified the parameter consistency (PC) test and the parameter goodness of fit (PG) as useful statistical methods to evaluate the compatibility of different data sets in a global analysis.

We have found that in (2+2) mass schemes recent solar and atmospheric data are compatible only at the 3.5$`\sigma `$ according to the PC test, and at 4.8$`\sigma `$ according to the PG method. In (3+1) mass schemes the disagreement of LSND with the rest of the oscillation data has been evaluated for both LSND analyses as a function of $`\mathrm{\Delta }m^2lsnd`$. For the global LSND data PC is achieved only at 2$`\sigma `$ and the PG is below 1% for all values of $`\mathrm{\Delta }m^2lsnd`$. Using the LSND DAR sample the disagreement is always stronger: we find a PC only at 99% C.L. and for most values of $`\mathrm{\Delta }m^2lsnd`$ the PG is worse than 4$`\sigma `$.

We have evaluated the g.o.f. of a global fit in terms of the PG by dividing the data in SOL, ATM, LSND and NEV samples. In the best case we find a PG of only 0.56%. This value occurs for the (3+1) mass scheme and the global LSND data. Using the LSND DAR analysis we get a PG of $`7.6\times 10^5`$ for (3+1). For (2+2) oscillation schemes the situation is worse: we find a very bad PG of $`1.6\times 10^6`$ for LSND global and $`3.5\times 10^7`$ for LSND DAR. Concerning the relative status of the hypotheses (3+1), (2+2) and the three-active neutrino case (3+0) we find that for the LSND global data (2+2) and also (3+0) are strongly disfavoured with respect to (3+1) with a $`\mathrm{\Delta }\chi ^2=17.8`$ and $`\mathrm{\Delta }\chi ^2=20.0`$, respectively. For LSND DAR we find for (2+2) a $`\mathrm{\Delta }\chi ^2=11.5`$ relative to (3+1), and the high significance of the oscillation signal condensed in the DAR sample leads to the huge value of $`\mathrm{\Delta }\chi ^2=28.5`$ of (3+0) relative to (3+1).

The exclusion of four-neutrino oscillation schemes of the (2+2)-type is based on the improved sensitivity of solar and atmospheric neutrino experiments to oscillations into a sterile neutrino, thanks to recent experimental data. This is a very robust result, independent of whether the LSND experiment is confirmed or disproved. The exclusion of (3+1) schemes depends somehow on the used LSND data sample666We want to stress that, although (3+1) models themselves are not ruled out, they do not offer an acceptable framework for a combined description of current oscillation data, if LSND is included. In contrast, any model of the (2+2)-type is not viable if the gap separating the two pairs of neutrino states is large compared to $`\mathrm{\Delta }m^2atm`$.. Furthermore, it heavily relies on the results of negative SBL experiments, especially on the Bugey and CDHS disappearance experiments. Therefore, if LSND should be confirmed by MiniBooNE, it will be crucial to improve the experimental data on SBL $`\stackrel{()}{\nu }_e`$ and/or $`\stackrel{()}{\nu }_\mu `$ disappearance. If the signal in LSND should indeed stem from oscillations in a (3+1) mass scheme the SBL $`\stackrel{()}{\nu }_e`$ and $`\stackrel{()}{\nu }_\mu `$ disappearance amplitude must be just on the edge of the sensitivity of current experiments.

In summary, the interpretation of the global neutrino oscillation data including LSND in terms of four-neutrino mass schemes – in either (3+1) or (2+2) realizations – is strongly disfavoured. In the best case we obtain a parameter g.o.f. of only 0.56%. The analysis we have presented brings the LSND anomaly to a theoretically more puzzling status indeed. We want to note that also introducing more sterile neutrinos participating in the oscillations is unlikely to substantially improve the situation peres . If LSND should be confirmed by the results of the MiniBooNE experiment the situation will become even more puzzling. Examples of more far-fetched alternative explanations are the possibility of lepton number violating muon decay Babu:2002ic  or large CPT violation in the neutrino sector CPT . Such scenarios will be crucially tested at the upcoming experiments MiniBooNE MiniBooNE  and KamLAND kamland .

## Appendix A Statistical compatibility of different data sets

In the global analysis of oscillation data in the framework of four-neutrino schemes we are faced with the problem that for a given mass scheme some data are in conflict with each other. In this appendix we present two statistical methods to evaluate quantitatively whether such data sets are consistent with each other.

Suppose a theoretical model characterized by $`P`$ parameters $`x=(x_1,\mathrm{},x_P)`$ and $`K`$ statistically independent experiments (or data sets) described by the theory. Each of these experiments may consist of many data points; experiment $`k`$ observes $`N_k`$ data points. Furthermore, experiment $`k`$ may in general depend only on $`P_kP`$ parameters $`x^k`$, which is a sub-set of all parameters $`x`$ of the full theory, leading to the $`\chi ^2`$-function

$$\chi _k^2(x^k)=(\chi _k^2)_{\mathrm{min}}+\mathrm{\Delta }\chi _k^2(x^k).$$
(8)
The three terms are distributed with $`N_k,N_kP_k,P_k`$ degrees of freedom (d.o.f.). The global $`\chi ^2`$ is given by

$$\chi _{\mathrm{tot}}^2(x)=\underset{k=1}{\overset{K}{}}\chi _k^2(x^k)=\chi _{\mathrm{min}}^2+\mathrm{\Delta }\chi ^2(x).$$
(9)
Usually the goodness-of-fit (g.o.f.) is evaluated by considering the value of $`\chi _{\mathrm{min}}^2`$, which is distributed with $`NP`$ d.o.f., where $`N=_kN_k`$.

Let us assume, that all the experiments considered separately lead to a good fit: $`(\chi _k^2)_{\mathrm{min}}N_kP_k`$. However, there may be some disagreement between different experiments. If the number of data points is very large the disagreement may be completely washed out and the usual g.o.f. would be very good, i.e. $`\chi _{\mathrm{min}}^2NP`$ despite the disagreement.777For the global four-neutrino analysis this effect was discussed explicitly in Ref. Maltoni:2001bc . Here we propose two methods which are more useful in global analyses with many data points, like in the four-neutrino case.

### A.1 Parameter consistency test (PC)

Let us consider the case of two data sets ($`K=2`$), which have $`\overline{P}P`$ parameters $`\overline{x}`$ in common. Then at a given C.L. the allowed regions for the parameters $`\overline{x}`$ according to each data set are determined by two $`\chi ^2`$ functions $`\mathrm{\Delta }\chi _1^2(\overline{x})`$ and $`\mathrm{\Delta }\chi _2^2(\overline{x})`$, where we minimize with respect to the other parameters. Now we calculate $`\chi _{\mathrm{PC}}^2`$, defined as the minimal $`\mathrm{\Delta }\chi ^2`$ value where there is still some overlap between the allowed parameter regions of the two data sets. This $`\chi _{\mathrm{PC}}^2`$ has to be evaluated for $`\overline{P}`$ d.o.f. to obtain the C.L. at which there exists still some values of the parameters such that both data sets are consistent.

We call this criterion parameter consistency test. It can be applied if the global data can be divided into two potentially conflicting sets. Similar methods have been used in Refs. GS , Maltoni:2001mt  to evaluate the consistency of LSND and NEV experiments in (3+1) schemes and in Ref. Church:2002tc  to check the compatibility of LSND and KARMEN.

### A.2 Parameter goodness of fit (PG)

Let us write the global $`\chi ^2`$ Eq. (9) in the following way:

$$\chi _{\mathrm{tot}}^2(x)=\underset{k=1}{\overset{K}{}}(\chi _k^2)_{\mathrm{min}}+\underset{k=1}{\overset{K}{}}\mathrm{\Delta }\chi _k^2(x^k).$$
(10)
To evaluate the g.o.f. we now use only the second term:

$$\overline{\chi }^2(x)\underset{k=1}{\overset{K}{}}\mathrm{\Delta }\chi _k^2(x^k)=\overline{\chi }_{\mathrm{min}}^2+\mathrm{\Delta }\overline{\chi }^2(x).$$
(11)
The first term in Eq. (10) describes the g.o.f. of all the experiments taken alone, whereas $`\overline{\chi }_{\mathrm{min}}^2`$ in Eq. (11) corresponds to the contribution to the total minimum $`\chi _{\mathrm{min}}^2`$ which stems from the combination of the data sets. The three terms in Eq. (11) are distributed with $`_kP_k,_kP_kP,P`$ d.o.f., respectively. Hence, to check the g.o.f. we just have to evaluate

$$\chi _{\mathrm{PG}}^2\overline{\chi }_{\mathrm{min}}^2\text{for}\underset{k=1}{\overset{K}{}}P_kP\text{d.o.f..}$$
(12)
We call this criterion parameter goodness of fit and in general it is more stringent than the parameter consistency test (for an illustration see Fig. 6).

If $`M_pK`$ experiments depend on the parameter $`x_p`$ we have the obvious relation $`_{k=1}^KP_k=_{p=1}^PM_p`$, and the relevant number of d.o.f. given in Eq. (12) can be written as $`_{p=1}^PM_pP=_{p=1}^P(M_p1)`$. From this relation we see that parameters relevant for only one experiment ($`M_p=1`$) do not contribute because they are minimized and contribute only to the individual $`(\chi _k^2)_{\mathrm{min}}`$ which are not considered in $`\chi _{\mathrm{PG}}^2`$. What is tested is only the parameter dependence of the data sets relevant for the combination.

A simple example: Let us consider the case of only one parameter $`y`$, which is determined by 2 experiments. Hence, we have $`K=2`$ and $`P=P_1=P_2=1`$ and $`\chi _{\mathrm{PG}}^2=\text{min}\left[\mathrm{\Delta }\chi _1^2(y)+\mathrm{\Delta }\chi _2^2(y)\right]`$, which has to be evaluated for 1 d.o.f.. The following statement is easy to prove: Let $`y_1`$ and $`y_2`$ be the best fit values of the two experiments ($`\mathrm{\Delta }\chi _k^2(y_k)=0`$). Then the C.L. at which $`y_1y_2`$ is consistent with zero is identical to the parameter goodness of fit.

## Acknowledgments

We are grateful to W.C. Louis and G. Mills for providing us with the table of the global LSND likelihood function. Further, we thank K. Eitel for providing us with the likelihood functions of the LSND DAR and the KARMEN analyses and for very useful discussions. This work was supported by Spanish grant PB98-0693, by the European Commission RTN network HPRN-CT-2000-00148 and by the European Science Foundation network grant N. 86. T.S. has been supported by the DOC fellowship of the Austrian Academy of Science and, in the early stages of this work, by a fellowship of the European Commission Research Training Site contract HPMT-2000-00124 of the host group. M.M. was supported by contract HPMF-CT-2000-01008 and M.A.T. was supported by the M.E.C.D. fellowship AP2000-1953.


format: text

text:
# On the applicability of jellium model to the description of alkali clusters

## I Introduction

During the last decade, investigation of the detailed structure and properties of small sodium clusters attracted a lot of attention (see, e.g., MetCl99 ; LesHouches ; StructNa  and references therein), because namely the sodium clusters were used in such important experimental work as the discovery of metal cluster electron shell structure Knight84  and the observation of plasmon resonances Br89 ; Selby89 ; Selby91 . These experiments were definitely among those, which clearly demonstrated that atomic clusters and small nanoparticles are in fact new physical objects possessing their own properties. The novelty of cluster physics is also greatly connected with the fact that cluster properties explain the transition from single atoms or molecules to solid state. Comprehensive survey of the field can be found in review papers and books, see, e.g., deHeer93 ; Brack93 ; BrCo94 ; Haberland94 ; Guet97 ; MetCl99 ; LesHouches .

With the discovery of electronic shell structure in free alkali clusters Knight84 ; Ekky84  the essential role of the quantized motion of delocalized valence electrons in the mean field created by ions in a cluster has been understood. Under different experimental conditions, the detailed ionic structure has been found not to affect the properties of alkali and other simple metal clusters very much (see, e.g., deHeer93 ; Brack93 ; BrCo94  for review). This behavior suggests the validity of a jellium model, defined by a Hamiltonian which treats the electrons in the usual quantum mechanical way, but approximates the field of the ionic cores by treating them as a uniform positively charged background. This naturally leads to a description of the electron density in terms of single particle wave functions that extend over the entire cluster.

Initially, jellium calculations for metal clusters were based on the density functional formalism with the use of pseudopotentials for the description of electron relaxation effects and lattice structure Martins81 . Fully self-consistent calculations for spherical jellium metal clusters have been performed within the framework of the spin-density-functional method Hintermann83  and the Kohn-Sham formalism for the self-consistent determination of electron wave functions Ekky84 ; Ekky85 . The Hartree-Fock (HF) scheme for the self-consistent determination of the electron wave functions of spherical jellium metal clusters was also introduced later in GJ92 ; IIKZ93 .

Shortly after the discovery of electronic shell structure in free alkali clusters Knight84 ; Ekky84  it was realized that the detailed size dependence of ionization potentials and other characteristics of small metallic clusters can be understood as a consequence of non-spherical cluster shapes Clemenger85  by analogy with the nuclear shell model. Direct evidence for cluster deformation was achieved in experiments on photoabsorption, where splitting of plasmon resonances caused by cluster deformation was observed (see deHeer93 ; Brack93 ; BrCo94  and references therein).

Kohn-Sham calculations for spherical Ekky84 ; Ekky85  metal clusters have been generalized for spheroidal Ekky88 ; Ekky91  and more general axial shapes Montag95  of light clusters. Light clusters of arbitrary shapes have been studied by means of the ”ultimate” jellium model Koskinen95 ; Koskinen95a . The Hartree-Fock approach was generalized for axially deformed cluster systems in LSSCG00 ; Lyalin01a . The existence of the different shape isomers were discussed in Montag95 ; Koskinen95a ; Hirschmann93 . It was shown that the shape of the magic clusters can deviate from the sphere if higher multipole deformations are taken into account. For example, the magic cluster $`Na_{40}`$ is not spherical if octupole deformations are allowed. This possibility have been pointed out in Hamamoto91  and confirmed by the ‘ultimate’ jellium Reimann97  and the Born-Oppenheimer local-spin density molecular dynamics method Rytkonen98 . It has been shown that alkali-metal clusters have similar shapes with small atomic nuclei Hakkinen97 . This similarity is a universal result of the density-functional theory.

The shell-correction method known from nuclear physics Strutinsky67  has been used in a number of papers in studying spheroidal ground state deformations, energetics, stability towards various fragmentation channels of metal clusters Yannouleas93a ; Yannouleas93b ; Reimann93 ; Bulgac93 . The role of higher multipole and tri-axial deformations has been elucidated in Lauritsch91 ; Frauendorf93 ; Yannouleas95 ; Frauendorf95 ; Frauendorf96 . Using the shell-correction method, the investigation of metal cluster electronic properties such as binding energies, ionization potentials, electron affinities, energetics of fission channels and systematic comparison of the theoretical results with the available experimental data have been done in Yannouleas95  (see also Yannouleas\_in\_MetCl  and references therein for a review).

Dynamical jellium model for metal clusters, which treats simultaneously collective vibrational modes (volume vibrations, i.e. breathing, and shape vibrations) of the ionic jellium background in a cluster, quantized electron motion and interaction between the electronic and ionic subsystems was developed in GSG99 ; GISG00 . This model allowed the widths of electron excitations in metal clusters beyond the adiabatic approximation to be described.

The jellium model provides a very useful basis for studying various collision processes, such as photabsorption Ivanov96 , photoionization Bertsch91 ; C20-C60 , elastic GCSG97 ; GEMS98  and inelastic scattering GEMS98 ; GIS97 ; GISG98 ; GIPS00 , electron attachment CGIS98 ; CGIS99 , photon emission GS97 ; GIS98 , atomic cluster fission process LSGS02 ; LSG02  and others, involving metal clusters. On the basis of the jellium model one can develop ab initio many-body theories, such as the random phase approximation with exchange or the Dyson equation method and effectively solve many-electron correlation problem even for relatively large cluster systems containing up to 100 atoms or even more. Review of these methods in their application to the electron scattering of metal clusters one can find in AVSol . As elucidated in the papers cited above, many-electron correlations are quite essential for the correct description of various characteristics of the cluster systems.

Structural properties of small metal clusters have been widely investigated using quantum chemistry ab initio methods. Here we refer to the papers Boustani87 ; Bonacic88 ; Boustani88 ; Guest88 ; Martins85 ; Spieg98 ; Nagueira99 ; Gutierrez01 ; StructNa , in which optimized geometries, binding energies, ionization potentials, electron structure and electron transport properties of small lithium and sodium clusters have been calculated.

In spite of the fact that both jellium model results and results of ab initio frameworks do exist in literature there have been performed no systematic comparison of the results of the two different theoretical schemes so far. We fill this gap in our present paper and demonstrate that such a comparison is rather illustrative and explains essential physical aspects of the formation of various cluster characteristics and properties. Also, we compare the results of the Hartree-Fock and local-density approximation (LDA) deformed jellium models and on this basis elucidate the role of many-electron correlation effects in the formation of cluster deformations.

On the basis of comparison of the jellium model results with those derived within the ab initio theoretical framework and experiment we elucidate the level of applicability of the jellium model to the description of alkali cluster properties. For neutral and singly charged sodium clusters with $`N20`$, we have calculated on the basis of Hartree-Fock and LDA deformed jellium model the binding energies per atom, ionization potentials, parameters of deformation and optimized values of the Wigner-Seitz radii. These characteristics are compared with the results derived from the ab initio all-electron theoretical framework for the calculation of the cluster ionic and electronic structure based on the density functional theory as well as on the post Hartree-Fock perturbation theory on many-electron correlation interaction. Comparison performed in our work demonstrates the great role of cluster deformations in the formation cluster properties and the quite reasonable level of applicability of the deformed jellium model.

Our paper is organized as follows. In section II, we provide a brief review of theoretical approaches and methods used in the calculation. In section III, we present and discuss jellium and ab initio results and make their comparison. In section IV, we draw a conclusion to this paper.

We use the atomic system of units $`\mathrm{}=|e|=m_e=1`$ in this paper.

## II Theoretical methods

In this work we calculate the binding energies per atom, ionization potentials, deformation parameters and the optimized values of the Wigner-Seitz radii for neutral and singly charged sodium clusters with $`N20`$ using the Hartree-Fock and LDA deformed jellium model. The jellium model results are compared with those derived in StructNa  on the basis of ab initio all-electron simulations of the cluster electronic and ionic structure. Below we present a brief review of the theoretical methods used in our work. Since the main part of calculations have been done within the framework of the deformed jellium model, we focus in our brief review on the jellium model approach rather than on the details of the more sophisticated ab initio methods for which we refer to StructNa .

### II.1 Two-centered jellium model and cluster shape parameterization

According to the main postulate of the jellium model, the electron motion in a metallic cluster takes place in the field of the uniform positive charge distribution of the ionic background. Originally, the Hartree-Fock model for metal-cluster electron structure has been worked out in the framework of spherically symmetric jellium approximation in GJ92 ; IIKZ93 . It is valid for clusters with closed electronic shells that correspond to magic numbers (8, 20, 34, 40, 58,…). For metal clusters with arbitrary number of valence electrons an open-shell two-center jellium Hartree-Fock approximation has been developed (see LSSCG00 ; Lyalin01a ; LSGS02 ). The two-centered jellium method treats the quantized electron motion in the field of the spheroidal ionic jellium background, whose principle diameters $`a`$ and $`b`$ can be expressed as follows:

$$a=\left(\frac{2+\delta }{2\delta }\right)^{2/3}R,b=\left(\frac{2\delta }{2+\delta }\right)^{1/3}R.$$
(1)
$`R=r_sN^{1/3}`$ is the radius of an undeformed spherical cluster with $`N`$ atoms, $`r_s`$ is the Wigner-Seitz radius, which for the bulk sodium is equal to 4.0. The deformation parameter $`\delta `$ characterizes the families of the prolate ($`\delta >0`$), and the oblate ($`\delta <0`$) spheroids of equal volume $`V_c=4\pi ab^2/3=4\pi R^3/3`$.

The electrostatic potential $`U`$ of the ionic background can be determined from the solution of the corresponding Poisson’s equation:

$$\mathrm{\Delta }U(𝐫)=4\pi \rho (𝐫),$$
(2)
where

$$\rho =\{\begin{array}{cc}\rho _c,\hfill & (x^2+y^2)/b^2+z^2/a^21\hfill \\ 0,\hfill & (x^2+y^2)/b^2+z^2/a^2>1\hfill \end{array}$$
(3)
describes a uniform distribution of the ions charge density in the volume of the cluster. Here $`\rho _c=Z_c/V_c`$ is the ionic charge density inside the cluster, and $`Z_c`$ is the total charge of the ionic core.

### II.2 Hartree-Fock and LDA formalism

The Hartree-Fock equations can be written out explicitly in the form (see, e.g., Lindgren )

$$\left(\mathrm{\Delta }/2+U+U_{HF}\right)a>=\epsilon _aa>.$$
(4)
The first term here represents the kinetic energy of electron $`a`$, and $`U`$ its attraction to the cluster core. The Hartree-Fock potential $`U_{HF}`$ represents the average Coulomb interaction of electron $`a`$ with the other electrons in the cluster, including the non-local exchange interaction, and $`\epsilon _a`$ describes the single electron energy.

According to the density-functional theory, the ground state energy reaches its minimum as a function of the density of the system at the exact density Hohenberg-Kohn . A self-consistent method for calculation of the electronic states of many-electron systems was proposed by Kohn and Sham Kohn-Sham . This method leads to the Kohn-Sham LDA self-consistent equations:

$$\left(\mathrm{\Delta }/2+U+U_H+V_{xc}\right)a>=\epsilon _aa>.$$
(5)
Here $`U_H`$ is the Hartree potential, which represents the direct Coulomb interaction of electron $`a`$ with other electrons in the cluster, but does not take into account the non-local exchange effects, while $`V_{xc}`$ is the phenomenological density dependent local exchange-correlation potential. The important feature of the LDA method consist in the fact that it takes into account many-electron correlations (see, e.g., Parr-book ; Fulde-book  for review). In the present work we use the Gunnarsson and Lundqvist model Gunnarsson  for the LDA electron exchange-correlation energy density $`ϵ_{xc}`$, which reads as

$$ϵ_{xc}(\rho _{el}(𝐫))=\frac{3}{4}\left(\frac{9}{4\pi ^2}\right)^{1/3}\frac{1}{r_s(𝐫)}0.0333G\left(r_s(𝐫)/11.4\right).$$
(6)
Here $`r_s(𝐫)=(3/4\pi \rho _{el}(𝐫))^{1/3}`$ is a local Wigner-Seitz radius, while $`\rho _{el}(𝐫)`$ is the electron density in the cluster, and the function $`G(x)`$ is defined by following relation:

$$G(x)=(1+x^3)\mathrm{ln}\left(1+\frac{1}{x}\right)x^2+\frac{x}{2}\frac{1}{3}.$$
(7)
The exchange-correlation energy density $`ϵ_{xc}(\rho _{el}(𝐫))`$, defines the LDA exchange-correlation potential $`V_{xc}(\rho _{el}(𝐫))`$ as

$`V_{xc}(\rho _{el}(𝐫))`$ $`=`$ $`{\displaystyle \frac{\delta \left[\rho _{el}(𝐫)ϵ_{xc}(\rho _{el}(𝐫))\right]}{\delta \rho _{el}(𝐫)}}=`$
$`\left({\displaystyle \frac{9}{4\pi ^2}}\right)^{1/3}{\displaystyle \frac{1}{r_s(𝐫)}}0.0333\mathrm{ln}\left(1+{\displaystyle \frac{11.4}{r_s(𝐫)}}\right).`$
The Hartree-Fock (4) and LDA (5) equations have been solved in the system of the prolate spheroidal coordinates KPS  as a system of coupled two-dimensional second order partial differential equations. The partial differential equations have been discretized on a two-dimensional grid and the resulting system of linear equations has been solved numerically by the successive overrelaxation method SOR . The third dimension, the azimuthal angle has been treated analytically.

An important characteristic of the cluster, which defines its stability is the total energy $`E_{tot}(N,\delta )`$. The total energy depends on the size $`N`$ of the cluster and its core deformation $`\delta `$. The total energy $`E_{tot}(N,\delta )`$ is equal to the sum of the electrostatic energy of the ionic core $`E_{core}(N,\delta )`$ and the energy of the valence electrons $`E_{el}(N,\delta )`$:

$$E_{tot}(N,\delta )=E_{core}(N,\delta )+E_{el}(N,\delta ).$$
(9)
The electrostatic energy of the cluster ionic core is equal to

$$E_{core}(N,\delta )=\frac{1}{2}_V\rho (𝐫)U(𝐫)d𝐫.$$
(10)
In the HF approximation, the electronic energy $`E_{el}(N,\delta )`$ is given by the general expression Lindgren :

$`E_{el}^{HF}(N,\delta )`$ $`=`$ $`{\displaystyle \underset{a}{}}<a\mathrm{\Delta }/2+Ua>+`$ (11)
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{abk}{}}q_aq_b\left[c(abk)F^k(a,b)+d(abk)G^k(a,b)\right],`$
where $`a`$ and $`b`$ run over all shells. The values $`F^k(a,b)`$ and $`G^k(a,b)`$ in the Eq. (11) are the Coulomb and exchange Slater integrals, $`q_a`$ and $`q_b`$ are the occupation numbers for orbitals $`a`$ and $`b`$, respectively. The Hatree-Fock coefficients $`c(abk)`$ and $`d(abk)`$ for the Coulomb and exchange energy contributions depend on the occupation numbers (see for details Lindgren ).

In the framework of LDA the electronic energy of the system is given by Hohenberg-Kohn ; Kohn-Sham :

$`E_{el}^{LDA}(N,\delta )`$ $`=`$ $`{\displaystyle \underset{a}{}}<a\mathrm{\Delta }/2+Ua>+`$ (12)
$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\rho _{el}(𝐫)\rho _{el}(𝐫^{})}{|𝐫𝐫^{}|}d𝐫d𝐫^{}}+{\displaystyle \rho _{el}(𝐫)ϵ_{xc}(\rho _{el}(𝐫))d𝐫},`$
where the latter term represents the exchange-correlation energy.

### II.3 All-electron ab initio Hartree-Fock and LDA calculations

When performing all-electron ab initio Hartree-Fock or LDA calculations for metal clusters, one has to solve equations which are symbolically the same to the Hartree-Fock (4) and the Kohn-Sham (5) equations initially written for the jellium model case. In ab initio calculations, instead of using the spheroidal jellium parameterization (3) for the distribution of the ionic charge density and solving the equations only for the valence electrons, one has to deal with the exact Coulomb forces of all the ions and to solve equations (4) or (5) assuming that all electrons are present in the system. The explicit form of the exchange-correlation density functional in equations (5) can be chosen differently. There are many different functional forms, although there is no unique one. For example in StructNa , the ab initio results have been obtained with the use of the gradient-corrected Becke-type three-parameter exchange functional Becke88  paired with the gradient-corrected Lee, Yang and Parr correlation functional (B3LYP) LYP ; Parr-book . The post Hartree-Fock theories accounting for many-electron correlations, such as for example the Møller-Plesset perturbation theory of the fourth order (MP4) MP  have also been used in StructNa  for metal cluster simulations.

In this paper we do not present the explicit forms of these functionals and omit the discussion of the methods used for the solution of the Hartree-Fock and LDA equations in the ab initio framework. Instead, we refer to the papers, where these functionals and methods are presented (see, e.g., StructNa ; Becke88 ; LYP ; Gaussian98\_man ; VWN ; PerWan ; MP ; LL3  and references therein).

## III Numerical results and discussion

In this section we present the results of systematic numerical calculations performed for the small sodium clusters in the size range $`N20`$ on the basis of deformed Hartree-Fock and LDA jellium model. We determine the binding energies per atom, ionization potentials and cluster deformations. These Hartree-Fock and LDA results are compared with each other, with the available experimental data, and with the results derived from the ab initio theoretical framework. This comparison elucidates the important role of many-electron correlations in metal clusters and establishes the level of applicability of the deformed jellium model.

### III.1 Total cluster energy minimization

In the Hartree-Fock and LDA jellium models the total cluster energy (11) and (12) depends on the deformation parameter $`\delta `$ introduced in (1).

Varying the total cluster energy on $`\delta `$ one can find its minimum at certain $`\delta `$ for each electron configuration considered. This value of $`\delta `$ corresponds to the stable geometrical configuration of the cluster.

In figures 1 and 2 we present the total cluster energy per atom $`E_{tot}(N,\delta )/N`$ calculated as a function of deformation parameter $`\delta `$ for the neutral and singly charged sodium clusters in the size range $`N20`$.

Solid lines in figures 1 and 2 present $`E_{tot}(N,\delta )/N`$ for those electronic configurations of clusters, that result in the absolute minimum of the total cluster energy. The $`\delta `$-value corresponding to the minimum of $`E_{tot}(N,\delta )/N`$ characterizes the cluster shape at the equilibrium point. Dashed curves show $`E_{tot}(N,\delta )/N`$ for the electronic configurations providing the energy minimum the most closely located to the absolute minimum of the total energy. These figures demonstrate that for many clusters with open electron shells both oblate and prolate isomers are possible and have close energies. The type of deformation that develops in the cluster is determined by the type of the corresponding electronic configuration. If the electronic orbitals characterizing the chosen electronic configuration are alongated with respect to the cluster axis of symmetry then the prolate deformation of the ionic jellium background is preferable. In the opposite case, when electronic orbitals lay mostly at the plane perpendicular to the cluster axis of symmetry the oblate deformation of the cluster becomes energetically more favorable.

It is interesting that the deformed jellium model predicts the existence of a non-spherical isomer for the $`Na_{20}`$ magic cluster. So for $`Na_{20}`$, the second minimum in $`E_{tot}(N,\delta )/N`$ arises at $`\delta 0.37`$, while the electronic configuration keeps the same. In spite of the fact that this second local minimum is energetically unfavorable, this prediction is qualitatively correct, because it corresponds to the result of ab initio calculations StructNa , which prove the existence of the two cluster isomers (pyramid and deformed structure based on the two linked icosahedrons) with the close energies. The pyramid cluster isomer possesses the tetrahedral group of symmetry, which is rather high. Thus, this isomer is analogous to the spherical cluster in the jellium picture. The deformed icosahedral structure corresponds to the prolate jellium cluster configuration. This correspondence would probably be even better if one allows the tri-axial deformations in the jellium approach.

The physical reason for the non-monotonous behaviour of the total energy $`E_{tot}(N,\delta )/N`$ with increasing the deformation parameter $`\delta `$ is the strong mixing of the highest occupied $`4\sigma `$ state with the unoccupied $`5\sigma `$ state. These levels exhibit a rather peculiar behaviour with the deformation parameter $`\delta `$ avoiding each other at $`\delta 0.37`$. Such an avoided crossing-point is linked to the Wigner’s no-crossing theorem wigner , which states that two energy levels of the same symmetry cannot cross. Thus, the strong interaction of the $`4\sigma `$ state with the more prolate type $`5\sigma `$ state results in the appearance of the second minima at $`\delta 0.37`$ for $`Na_{20}`$ cluster.

Figure 2 demonstrates that similar behaviour of the total energy upon deformation parameter observe for singly charged cluster ions. This figure shows that in the jellium approach the closed shell cluster ions, $`Na_9^+`$ and $`Na_{21}^+`$, have the spherical shape, $`\delta =0`$, at the equilibrium point, similarly to the neutral magic clusters $`Na_8`$ and $`Na_{20}`$.

### III.2 Cluster deformations

In the axially deformed jellium model clusters can either be spherical or have a shape of ellipsoid of revolution (spheroid). The spheroidal shape can be of the two types prolate or oblate, depending on the sign of the deformation parameter $`\delta `$ introduced in (1). In the ab intio approach, the cluster shape is determined by the optimized coordinates of all the ions and it can be characterized by the tensor $`R_{ij}`$

$$R_{ij}=x_ix_j$$
(13)
Here, the summation is performed over all ions in the system. The principle values of this tensor $`R_{xx}`$, $`R_{yy}`$ and $`R_{zz}`$ define the dimensions $`R_x`$, $`R_y`$ and $`R_z`$ of the ionic charge distribution in the cluster along the principle axes $`x`$, $`y`$ and $`z`$. Note that tensor $`R_{ij}`$ is closely connected with the cluster moment of inertia tensor and the quadrupole moment tensor of the ionic distribution.

The tensor $`R_{ij}`$ can also be defined for the jellium model. In this case, sum in (13) should be replaced by the integral and the integration to be performed over the homogeneous spheroidal distribution of the ionic density in the cluster. Then, the principal values of the tensor $`R_{ij}`$ can easily be determined. The result of this calculation reads as

$$R_{xx}=R_{yy}=\frac{b^2}{5}N,R_{zz}=\frac{a^2}{5}N.$$
(14)
Here, $`a`$ and $`b`$ are the principle diameters of the spheroid defined in (1).

In figure 3, we present the principle values $`R_{xx}`$, $`R_{yy}`$ and $`R_{zz}`$ calculated for the neutral sodium clusters with $`N<20`$ in the framework of the deformed jellium model according to (14). The diameters $`a`$ and $`b`$ have been determined by minimizing the total cluster energy in the LDA approximation (12) as explained in the previous subsection. This calculation has been performed with $`r_s=4.0`$ (Fig. 3a), which corresponds to the density of the bulk sodium, and with optimized value of the Wigner-Seitz radius (Fig. 3b). The cluster energy minimization on Wiger-Seitz radius will be discussed in section III.6 in more detail. The LDA deformed jellium model results are shown in figure 3 by the filled triangles. The filled triangles pointing up correspond to $`R_{xx}`$=$`R_{yy}`$, while those pointing down to $`R_{zz}`$. The opened triangles are the results of the all-electron ab initio framework derived in StructNa  with the use of the B3LYP density functional. The opened triangles pointing up and down show $`R_{xx}`$ and $`R_{zz}`$ respectively, while the opened triangles pointing right represent $`R_{yy}`$. In figure 4 we present the results of similar calculations performed for singly charged ions. The notations used in figure 4 are the same as in figure 3.

Comparison of the results presented in figures 3 and 4, demonstrate that the optimization of the cluster energy on the Wigner-Seitz radius does not change much the neutral cluster geometries. For cluster ions, the alteration of the cluster shape with the variation of the Wigner-Seitz radius is more noticeable, although it does not improve the agreement of the jellium model and ab initio results. This comparison demonstrates that the Wigner-Seitz radius variation in the jellium model does not actually improve the quality of the model.

Figures 3 and 4 demonstrate rather good agreement of the jellium model and ab initio results. In most of the cases the jellium model predicts correctly the type of the dominant cluster deformation, prolate or oblate one. Of course, ab initio calculations include tri-axial deformations of the cluster, which turned out to be noticeable for the clusters with the open subshells and play important role for clusters like $`Na_{12}Na_{14}`$, $`Na_{17}`$, $`Na_{13}^+`$, $`Na_{14}^+`$, $`Na_{20}^+`$. The axially symmetric deformed jellium model does not take into account tri-axial deformations and thus in this case always $`R_{xx}=R_{yy}`$.

The axially symmetric jellium model gives the wrong type of deformation in the open shell clusters, like $`Na_5`$, $`Na_{16}Na_{19}`$, $`Na_6^+`$, $`Na_{17}^+Na_{20}^+`$. However, it is necessary to note that for all these clusters there are almost degenerate oblate and prolate isomers within the axially symmetric jellium model, as it is shown in figures 1 and 2 and in tables 1 and 2 (see Appendix). Thus, accounting for tri-axial deformations in these clusters plays the crucial role as it becomes clear from the comparison of the jellium and ab initio results.

For the magic clusters $`Na_8`$ and $`Na_{20}`$, the principle values $`R_{xx}=R_{yy}=R_{zz}`$ are almost identical in both approaches, which demonstrates the closeness to the sphericity of the ab initio magic cluster shapes. Note that for the magic cluster ion, $`Na_9^+`$, some small deformation remains in the ab initio approach, while in the jellium approach it turns out to be spherical. This demonstrates that the ionic framework of the cluster is not that deformable as it follows from the jellium model.

In Appendix, we compiled in tables the optimized Wigner-Seitz radii $`r_s`$, total energies per atom $`E_{tot}(N,\delta )/N`$, deformation parameters $`\delta `$ and the second derivatives of the total energy on cluster deformation at the $`\delta `$-point corresponding to minimum of the total energy. The later characteristic of the cluster, $`^2E_{tot}(N,\delta )/\delta ^2`$, is directly connected to the frequency of cluster surface vibrations (see, e.g., GISG00  for details). All values which are presented in the tables 1 and 2 have been calculated in the LDA deformed jellium model for neutral and singly charged clusters.

### III.3 Binding energies per atom

In this paper we calculate the dependence of binding energy per atom in the deformed jellium model and compare it with ab initio results from StructNa . The binding energies per atom for the neutral and singly charged clusters are defined as follows:

$`E_b/N`$ $`=`$ $`E_1E_N/N`$ (15)
$`E_b^+/N`$ $`=`$ $`\left((N1)E_1+E_1^+E_N^+\right)/N,`$ (16)
where $`E_N`$ and $`E_N^+`$ are the total energies of a neutral and singly-charged N-atomic jellium cluster respectively.

Figure 5 shows the dependence of the binding energy per atom for neutral (Fig. 5a) and singly charged (Fig. 5b) clusters as a function of cluster size calculated in the deformed jellium model. We compare the calculated dependences with the ab initio results from StructNa  obtained by the B3LYP and MP4 methods. In figure 5 we show the jellium model results obtained with bulk and optimized values of the Wigner-Seitz radius. It is seen that the cluster optimization on the Wigner-Seitz radius brings the cluster energies down and makes them closer to the ab initio results.

Figure 5 demonstrates that the general trend of the curves calculated within the jellium framework turns out to be very close to the one obtained from the ab initio calculation. The similarity of the the jellium and ab initio curves is higher for $`N10`$. In the region $`10N20`$ small discrepancy in the behaviour of the curves can be attributed to the tri-axial cluster deformations taken into account in the ab initio approach and omitted in the axially symmetric jellium model.

Note that the jellium model results for both neutral and singly charged sodium clusters are somewhat closer to the predictions of the $`MP4`$ method. This method is based on the accounting of the many-electron correlations up to the fourth order of the perturbation theory and is free of any adjustable parameters.

Figure 5 demonstrates that in spite of the simplicity, the jellium model turns out to be rather reliable approximation able to reproduce reasonably well the dependence of binding energy per atom for both neutral and singly charged sodium clusters.

### III.4 On the role of exchange and many electron correlation interaction.

To illustrate the importance of the exchange and many-electron correlation interaction in the cluster, we plot on figure 6 the total energy per atom calculated both in the HF and LDA spherical jellium approximations for neutral (Fig. 6a) and singly charged (Fig. 6b) spherical sodium clusters as a function of cluster size. The spherical cluster shape was assumed in this calculation for the sake of simplicity. Figure 6 shows the significant difference between the HF and LDA results for both neutral and singly charged clusters, which is the result of the accounting for the many-electron correlation interaction within LDA.

Figure 6 demonstrates that the LDA total energy dependence possesses the local minima at the shell closings $`N=2,8,20`$ for neutral clusters and $`N=3,9,21`$ for singly charged ones. The HF total energy curve has the extra minima at the half-shell closings, i.e. $`N=5,13,19`$ for neutral clusters and $`N=6,14,20`$ for singly charged cluster ions. The appearance of these extra minima is the result of the more accurate accounting for the exchange interaction within the HF approximation. We found that with increasing cluster size within the given shell the favorable electronic configuration of the cluster changes resulting in the formation of the extra minima on the HF total energy dependence. Qualitatively, this situation can be understood on the basis of the Hund’s rule which states that the lowest energy level in the system at a fixed electronic configuration is characterized by the maximum value of total spin and the maximum possible (at this spin) angular momentum. As an illustration, figure 7 shows the $`\delta `$-dependencies of the total energy per atom for $`Na_5`$ cluster obtained in the Hartree-Fock (Fig. 7a) and LDA (Fig. 7b) deformed jellium model for different electronic configurations. Figure 7a clearly demonstrates that the extra minima on HF total energy curve correspond to the electronic states in which spins of all the electrons from the open shell are parallel. The LDA framework, with the exchange-correlation potential (II.2), which we have used in our work does not take into account the spin polarization effects and thus we found no minima on the LDA total energy curve at closing of the half filled shells (see Fig. 7b).

### III.5 Ionization potentials

Another important characteristic of the cluster system is its ionization potential. The ionization potential is determined by the energy needed to take an electron out of the cluster. It is equal to:

$$V_i=E_N^+E_N$$
(17)
In figure 8, we present the ionization potential of neutral sodium clusters calculated within the jellium model as a function of cluster size. We compare the jellium model results with those obtained in StructNa  using ab initio theoretical framework and with the available experimental data deHeer93 . This comparison demonstrates that the jellium model reproduces correctly most of the essential features of the ionization potential dependence on $`N`$. Some discrepancy, like in the region $`11N14`$, can be attributed to the neglection of the tri-axial deformation in the axially symmetric jellium model.

In spite of the fact that ab initio results are closer to the experimental points, one can state quite satisfactory agreement of the jellium model results with the experimental data, which illustrates correctness of the jellium model assumptions and its applicability to the description of sodium clusters.

Figure 8 also demonstrates the role of cluster deformations on the formation of the odd-even oscillations in the dependence of the cluster ionization potential on $`N`$. Indeed, for spherically symmetric clusters this dependence turns out to be monotonous within the range of the given shell Ivanov96  contrary to the experimental observations. Allowing for the cluster deformation and introducing a single deformation parameter $`\delta `$, we have achieved much better agreement of theoretical results with the experimental data as it is clear from figure 8.

### III.6 Wigner-Seitz radius variation

Calculations of the cluster total energy are usually performed at the certain value of the Wigner-Seitz radius $`r_s`$. The bulk value of the Wigner-Seitz radius for sodium is equal to 4.0. However, one can also perform the calculation minimizing the total cluster energy by variation of the Wigner-Seitz radius.

Figure 9 demonstrates the dependence of the optimized Wigner-Seitz radii on cluster size calculated for neutral and singly charged sodium clusters within the HF and LDA approximations. This figure shows that the alteration of the optimized $`r_s`$ values is much larger for the cluster ions as compared to the neutral clusters. For neutral clusters, the optimized values are somewhat larger than the bulk value $`r_s=4.0`$ in both LDA and HF approximations. The LDA dependence goes closer to the bulk limit. With increasing $`N`$ this dependence approaches the bulk limit, being very close to it also for the magic numbers $`N=8`$ and $`N=20`$, which is another manifestation of the shell effect.

The difference between the optimized values of the Wigner-Seitz radii obtained in the HF and LDA approximations can be atributed to the manifestation of the many-electron correlation interaction in the system. To illustrate this fact, we plot in figure 10 the difference between the optimized Wigner-Seitz radii calculated in the LDA and HF approximations for neutral and singly charged sodium clusters. These dependences have the prominent peculiarities. The origin of these peculiarities is the same as for those in figure 6 (see section III.4 and discussion therein), although now we consider deformed cluster systems.

## IV Conclusion

In this paper we performed systematic calculation of various characteristics of neutral and singly charged sodium clusters with $`N20`$ on the basis of the deformed HF and LDA jellium models. We compared the results of our calculations with the ab initio results obtained in StructNa  and with the available experimental data. From these comparisons, we have established the level of applicability of the jellium model to the description of various cluster characteristics.

Our consideration shows that the deformed jellium model provides qualitatively correct description of the sodium clusters and their ions. The quantitatively reliable results with the accuracy below than 10 per cent one can expect from the jellium model description providing one allows for the tri-axial cluster deformations.

We have performed our calculations using HF and LDA approximations. Comparison of the results of the two approaches allowed us to demonstrate the importance of the many-electron correlations in the formation of cluster characteristics and properties.

We have performed our calculations for sodium clusters. However, most of the conclusions should be applicable to other alkali clusters, potassium for example. The level of applicability of the jellium approach to other metals, like alkali-earth, requires a separate careful consideration.

###### Acknowledgements.

The authors acknowledge support from the Alexander von Humboldt Foundation, DAAD and the Russian Academy of Sciences (Grant 44).

*

## Appendix A Tables

In Appendix, we present tables of the optimized Wigner-Seitz radii, total energies per atom, deformation parameters and the second derivatives of the total energy on cluster deformation calculated in the LDA deformed jellium model for neutral and singly charged sodium clusters.

## References


text:
# EXPERIMENTAL STATUS OF 𝑲→𝝅𝝂𝝂̄

## 1 Introduction

The primary focus in kaon physics today is the two golden modes: $`K^+\pi ^+\nu \overline{\nu }`$  and $`K_L^{}\pi ^{}\nu \overline{\nu }`$ . These modes are interesting as there is essentially no theoretical ambiguity in extracting fundamental CKM parameters from measurements of the branching ratios. The intrinsic theoretical uncertainty in B($`K^+\pi ^+\nu \overline{\nu }`$ ) is $``$7% and is even smaller in B($`K_L^{}\pi ^{}\nu \overline{\nu }`$ ), only $``$2%; in both cases the hadronic matrix element can be extracted from the $`K^+\pi ^{}e^+\nu _e`$  ($`K_{e3}`$) decay rate.

The unitarity of the CKM matrix can be expressed as

$$V_{us}^{}V_{ud}+V_{cs}^{}V_{cd}+V_{ts}^{}V_{td}=\lambda _u+\lambda _c+\lambda _t=0$$

with the three vectors $`\lambda _iV_{is}^{}V_{id}`$ converging to form a very elongated triangle in the complex plane. The length of first vector $`\lambda _u=V_{us}^{}V_{ud}`$ is precisely determined from $`K_{e3}`$ decay. The length of the third side $`\lambda _t=V_{ts}^{}V_{td}`$ is measured by $`K^+\pi ^+\nu \overline{\nu }`$ and the height of the triangle, $`Im\lambda _t`$, can be measured by $`K_L^{}\pi ^{}\nu \overline{\nu }`$ . Branching ratio measurements of the two $`K\pi \nu \overline{\nu }`$ modes, along with the well known $`K_{e3}`$, will completely determine the unitarity triangle.

Comparison of CKM parameters as measured from the golden $`K\pi \nu \overline{\nu }`$ , $`B_d^{}\psi K_S^{}`$ and $`\mathrm{\Delta }M_{B_d}/\mathrm{\Delta }M_{B_s}`$ modes, provide the best opportunity to over-constrain the unitary triangle and to search for new physics. In particular, comparisons of

* $`|V_{td}\text{ }|`$ from $`K^+\pi ^+\nu \overline{\nu }`$ and from the ratio of the mixing frequencies of $`B_d`$ and $`B_s`$ mesons $`\mathrm{\Delta }M_{B_d}/\mathrm{\Delta }M_{B_s}`$  ,
* $`\beta `$ from B($`K_L^{}\pi ^{}\nu \overline{\nu }`$ )/B($`K^+\pi ^+\nu \overline{\nu }`$ ) and from the time dependent asymmetry in the decay $`B_d^{}\psi K_S^{}`$

offer outstanding opportunities to explore the Standard Model (SM) picture of CP–violation.

The SM prediction for the $`K\pi \nu \overline{\nu }`$ branching ratios are B($`K^+\pi ^+\nu \overline{\nu }`$ ) = $`(0.72\pm 0.21)\times 10^{10}`$ and B($`K_L^{}\pi ^{}\nu \overline{\nu }`$ ) = $`(0.26\pm 0.12)\times 10^{10}`$ . In addition, an upper limit on B($`K^+\pi ^+\nu \overline{\nu }`$ ) can be derived with small uncertainty from the current limit on $`B_s`$ mixing; this limit is B($`K^+\pi ^+\nu \overline{\nu }`$ ) $`<1.32\times 10^{10}`$ .

## 2 $`𝑲^\mathbf{+}\mathbf{}𝝅^\mathbf{+}𝝂\overline{𝝂}`$

The Alternating Gradient Synchrotron (AGS) as a high-energy physics (HEP) research facility has had a broad and rich program in kaon physics, culminating in the observation of two $`K^+\pi ^+\nu \overline{\nu }`$ events by E787. With the recent successful commissioning of the Relativistic Heavy Ion Collider (RHIC), the primary role of the AGS has shifted to become an injector of heavy ions for RHIC. Nevertheless, the AGS remains the highest intensity proton synchrotron in the world and is designed to be available for $``$20 hours/day when not filling RHIC, and as such retains an important role in the US HEP program. DOE and BNL have approved and agreed to fund one new HEP experiment to run at the AGS between RHIC fills: the E949 experiment seeks to make a precise measurement of the branching ratio B($`K^+\pi ^+\nu \overline{\nu }`$ ). While DOE approved E949 to run for 60 weeks, the FY03 budget of the President of the United States does not include running for E949. At this point the E949 experiment will be terminated after only 12 weeks of running, unless the US Congress restores funding. The fate of the E949 experiment will not be known until Congress initiates and then completes the appropriations process between June and October of 2002.

The completion of the Main Injector (MI) at FNAL allows for the simultaneous operation of a fixed target program along with the Tevatron collider. The next step in the pursuit of $`K^+\pi ^+\nu \overline{\nu }`$ will be the CKM experiment, which plans to use a modest fraction of the MI protons ($`5\times 10^{12}`$), extracted over a 1-second spill with minimal bunching of the proton beam, and will push the $`K^+\pi ^+\nu \overline{\nu }`$ sensitivity to the current limits of theoretical precision.

### 2.1 E787

The first $`K^+\pi ^+\nu \overline{\nu }`$ signal was observed in the 1995 data sample of the E787 experiment . No new events were seen in the data sample from 1996–97 , and with a background of $`0.08\pm 0.02`$ events and a signal of one event a branching ratio of B($`K^+\pi ^+\nu \overline{\nu }`$ ) = $`(1.5_{1.2}^{+3.4})\times 10^{10}`$ was measured. That event was in fact in a very clean region of the predefined $`K^+\pi ^+\nu \overline{\nu }`$ signal region with a SM-signal to background ratio of 35. An analysis of the final E787 data sample from the 1998-99 run has recently been reported . With a measured background of $`0.066_{0.025}^{+0.044}`$, one new event was observed. The final plot of range vs. energy from the combined E787 1995–99 data sample for events passing all other cuts is shown in Figure 1.

The branching ratio, as determined from these two events, is

$$B(K^+\pi ^+\nu \overline{\nu }\text{ })=(1.57_{0.82}^{+1.75})\times 10^{10}.$$
(1)
This branching ratio is a factor of two higher than expected in the SM and is higher than allowed by the current limit on $`B_s`$ mixing. Of course, the uncertainty on the BR measurement is large due to limited statistics and new data from the E949 experiment are eagerly awaited.

The new event found in the 1998 data sample is in a relatively clean region of the accepted signal region: the SM signal to background ratio for this event is 3.6. An event display for this event, as well as the previous event, is shown in Figure 2.

Limits on $`|V_{td}|`$ and $`\lambda _t`$ can be obtained (these are 1-$`\sigma `$ limits except for $`Im\lambda _t`$ which is 90% CL),

$`0.007<`$ $`|V_{td}|`$ $`<0.030,`$ (2)
$`2.9\times 10^4<`$ $`|\lambda _t|`$ $`<1.2\times 10^3,`$
$`0.88\times 10^3<`$ $`Re\lambda _t`$ $`<1.2\times 10^3,`$
$`Im\lambda _t`$ $`<1.1\times 10^3.`$
Even with the large statistical error, this new measurement provides a non-trivial contribution to global fits of the CKM parameters . The constraints on $`\lambda _t`$ from this result are shown in Figure 3. The constraints from the other golden B modes, $`\mathrm{\Delta }M_{B_d}/\mathrm{\Delta }M_{B_s}`$ and $`B_d^{}\psi K_S^{}`$ are shown on the same plot. One can immediately see that $`|\lambda _t|`$ is tightly constrained to a narrow crescent by $`K^+\pi ^+\nu \overline{\nu }`$ and $`\mathrm{\Delta }M_{B_d}/\mathrm{\Delta }M_{B_s}`$ .

New data from the successor to E787, E949, will make a significant contribution to our knowledge of the CKM parameters.

In addition, E787 has searched for the decay $`K^+\pi ^+\nu \overline{\nu }`$ in the pion kinematic region below the $`K^+\pi ^+\pi ^{}`$ ($`K_{\pi 2}`$) peak . This region contains more of the $`K^+\pi ^+\nu \overline{\nu }`$ phase space, but is complicated by a significant background from $`K^+\pi ^+\pi ^{}`$ decays with the $`\pi ^+`$ scattering in the scintillating fiber target and shifting its kinematics into the search region. The data from the 1996 run of E787 has been analyzed ($``$20% of the entire E787 data sample). One event was observed in the search region, consistent with the background estimate of $`0.73\pm 0.18`$. This implies an upper limit on B($`K^+\pi ^+\nu \overline{\nu }`$ ) $`<4.2\times 10^9`$ (90% C.L.), and is consistent with the 2 events observed above the $`K_{\pi 2}`$ peak and the SM spectrum. Some additional reduction of the background levels in the remaining E787 data may be possible, but the major focus will shift to the new E949 experiment, which has significantly enhanced photon veto capabilities that will further suppress this background. In addition, the next experiment after E949, CKM at FNAL, will be essentially free of this background since there is no stopping target.

### 2.2 E949

E949 is an upgraded version of the E787 experiment, planning to capitalize on the full AGS beam to collect $`K^+\pi ^+\nu \overline{\nu }`$ data at 14 times the rate of the E787 run in 1995. The new detector has substantially upgraded photon veto capabilities, enhanced tracking, triggering, monitoring, and DAQ capability, and will run at a higher AGS duty factor and a lower kaon momentum (with an increased fraction of stopped kaons). It has been designed to reach a sensitivity of at least 5 times beyond E787 and observe 5–10 SM events. The background level for E949 measurement of B($`K^+\pi ^+\nu \overline{\nu }`$ ) above the $`K_{\pi 2}`$ peak is reliably projected from E787 data to be $``$10% of the Standard Model signal.

E949 should see up to 10 SM events (or 20 events at the branching ratio measured by E787) within the next couple of years. This is an exciting opportunity to make a significant contribution to quark mixing and CP–violation that should be fully exploited. A history of the search for $`K^+\pi ^+\nu \overline{\nu }`$ is shown in Figure 4.

### 2.3 CKM

The next step towards a precision measurement of B($`K^+\pi ^+\nu \overline{\nu }`$ ) will be the CKM experiment at FNAL. CKM has been given scientific (Stage–1) approval by FNAL and could be running by 2007. CKM plans to use a novel technique for $`K^+\pi ^+\nu \overline{\nu }`$ : a decay in flight experiment, with redundant kinematic constraints from a conventional momentum spectrometer and a novel velocity spectrometer based on RICH counters. CKM expects to observe 100 SM signal events in a two year run, using the Main Injector simultaneously with the Tevatron. The background is expected to be $``$10% of the SM signal, predominantly from $`K^+\pi ^+\pi ^{}`$ . A vacuum of $`10^6`$ Torr is required to minimize backgrounds from kaon interactions in the decay volume. CKM will require less than 20% of the flux from the Main Injector, but will require a slow extracted spill of $``$1 second duration. CKM will run with a 50 MHz 22 MeV/c beam with RF separators and about 70% K<sup>+</sup> purity. Figure 4 shows a projected measurement of B($`K^+\pi ^+\nu \overline{\nu }`$ ) from CKM, assuming the current central value of the branching ratio.

## 3 $`𝑲_𝑳^{\mathbf{}}\mathbf{}𝝅^{\mathbf{}}𝝂\overline{𝝂}`$

The decay $`K_L^{}\pi ^{}\nu \overline{\nu }`$ is even cleaner theoretically and is CP-violating. However, it is even more challenging experimentally as all of the particles involved are neutral.

Presently, the best limit on $`K_L^{}\pi ^{}\nu \overline{\nu }`$ is derived in a model-independent way  from the E787 measurement of $`K^+\pi ^+\nu \overline{\nu }`$ :

$`B(K_L^{}\pi ^{}\nu \overline{\nu }\text{ })`$ $`<`$ $`4.4\times B(K^+\pi ^+\nu \overline{\nu }\text{ })`$
$`<`$ $`1.7\times 10^9(90\%\mathrm{CL}).`$
Of course, it is desirable to observe this mode directly in order to extract a second constraint on the CKM matrix parameters. The current best direct limit is derived from a KTeV search for high transverse momentum $`\pi ^{}`$’s decaying via $`\pi ^{}e^+e^{}\gamma `$. From the full 1997 data set, KTeV observed no events with an expected background of $`0.12_{0.04}^{+0.05}`$ and set a 90%-CL limit  of

$$B(K_L^{}\pi ^{}\nu \overline{\nu }\text{ })<5.9\times 10^7.$$
(4)
Due to the small $`\pi ^{}e^+e^{}\gamma `$ branching ratio all future experiments plan to use the more copious $`\pi ^{}\gamma \gamma `$ mode. KTeV also made a search in this mode in a special one day test, with a highly collimated ‘pencil’ beam and observed one background event, most likely from a neutron interaction, and set a 90%-CL limit  of B($`K_L^{}\pi ^{}\nu \overline{\nu }`$ ) $`<1.6\times 10^6`$.

### 3.1 E391

The next generation of $`K_L^{}\pi ^{}\nu \overline{\nu }`$ experiments will start with E391a at the High Energy Accelerator Research Organization (KEK) which hopes to reach a sensitivity of $`3\times 10^{10}`$. This experiment will use a technique similar to KTeV, with a pencil beam, high quality calorimetry and very efficient photon vetos. This is the first experiment dedicated to searching for $`K_L^{}\pi ^{}\nu \overline{\nu }`$ and aims to close the window for non-SM contributions to the decay. It will also serve as a test bed for the experimental techniques necessary to observe $`K_L^{}\pi ^{}\nu \overline{\nu }`$ in future experiments. Beam tests were started in 2001 and the first data-taking run is scheduled for 2003. This experiment plans to eventually move to the Japanese Hadron Facility (JHF) in $``$2007, and attempt to push to a sensitivity of $`𝒪`$(1000) events.

### 3.2 KOPIO

The National Science Board of the National Science Foundation (NSF) has approved the construction of the two new large experiments at the BNL AGS: KOPIO and MECO, as components of the Rare Symmetry Violation Proposal (RSVP). RSVP is planned to be one of the next Major Research Equipment construction projects at the NSF.

The KOPIO experiment is designed to discover the $`K_L^{}\pi ^{}\nu \overline{\nu }`$ decay and measure its branching ratio to $``$20%. KOPIO will make use of a time-of-flight technique to measure the momentum of the $`K_L`$ and will operate at a large targeting angle to improve the $`p_K`$ resolution and soften the neutron spectrum to reduce $`\pi ^{}`$ hadroproduction. All possible aspects of the decay will be measured: the photon directions will be measured in a pre-radiator, the times and energy will be precisely measured in a Shashlyk calorimeter, the time of the kaon’s creation is determined by the proton bunch width from the AGS ($``$250 psec). KOPIO will have a substantial photon veto system and make the same sort of background measurements as E787, directly from the data, with two independent tools for attacking the major background, $`K_L^{}\pi ^{}\pi ^{}`$ , through both kinematics and photon veto. KOPIO expects to observe 50 SM events, with a background of 50%. This will allow a determination of $`Im\lambda _t`$ to 10%. KOPIO is expected to start data collection in $``$2006.

## 4 Conclusion

The next decade will be an exciting time for improved understanding of CP-violation and quark mixing. It is quite likely that precise measurements of all four golden modes will be made: $`B_d^{}\psi K_S^{}`$ at the B-factories; $`\mathrm{\Delta }M_{B_d}/\mathrm{\Delta }M_{B_s}`$ , most likely at the Tevatron; and $`K\pi \nu \overline{\nu }`$ at BNL, KEK and FNAL. These measurements will allow a precise determination of CKM parameters and provide a critical test of the SM picture of CP-violation.

## 5 Acknowledgements

I would like to thank Takao Inagaki, Bob Tschirhart, Taku Yamanaka, Takeshi Komatsubara, Shojiro Sugimoto, David Jaffe and Laurie Littenberg for discussions and/or comments on this paper. This work was supported in part under US Department of Energy contract #DE-AC02-98CH10886.