text:
# Predicting Multivariate Volatility

## Abstract

We suggest two classes of multivariate GARCH–models which are both easy to estimate and perform well in forecasting the covariance matrix of more than one hundred stocks. We apply methods from random matrix theory (RMT) to determine the number of principal components or the number of factors in the multivariate volatility models. In this way only statistically relevant information is used for the estimation of model parameters.

The basic input for the estimation of portfolio risk is a forecast of time varying covariance matrices. Although modern risk management is founded on the analysis of huge data bases and the use of sophisticated theoretical models, there has been only partial progress with respect to covariance forecasts. The reason for this lack of progress is the “curse of dimensionality”: in modern financial engineering, the dimension of the investment universe is comparable to the number of observations, and hence the number of parameters in the covariance or correlation matrix is of the same order as the number of data points.

Traditional multivariate statistics is not able to deal with this kind of problem. Both estimators and tests are consistent only in the limit of $`L/N`$ going to infinity, where $`L`$ is the number of observations and $`N`$ is the dimension of the observed vector Muirhead82 . In typical applications, one is dealing with ratios $`L/N`$ ranging from one to ten, and hence statistical fluctuations due to the lack of a sufficient number of observations dominate an empirically observed covariance matrix. This measurement noise induces randomness in the covariances and hence turns the covariance matrix at least partly in a random matrix.

A covariance matrix can be decomposed into the cross–correlation matrix C and the set of standard deviations $`\{\sigma _i^2\}`$. To identify the effects of randomness on the eigenvalue spectrum of an empirical cross–correlation matrix, one considers the null hypothesis of completely uncorrelated and identically distributed (i.i.d.) time series. The cross–correlation matrix calculated from these time series is a random matrix R from the so called Wishart ensemble. The eigenvalues of a random Wishart matrix fall in the interval $`[\lambda _{},\lambda _+]`$ with $`\lambda _\pm =1+N/L+2\sqrt{N/L}`$ Marcenko+67 . This eigenvalue distribution can be compared to that of the empirical cross–correlation matrix C. Agreement between R and C is a sign of randomness, whereas deviations indicate the presence of economically meaningful information. It was found that approximately 94 % of the empirically observed eigenvalues fall in the random matrix interval Laloux+99 . The random matrix nature of the main part of the spectrum was proven by comparing the eigenvalue statistics of empirical correlation matrices to the universal predictions of random matrix theory Plerou+99 . Only the largest eigenvalues and corresponding eigenvectors carry information about market correlations, business sectors and geographical correlations Gopikrishnan+01 ; Plerou+2002 . These economically meaningful correlations are stable in time and hence suitable for forecasting future correlations. Generally, correlations show a higher degree of time stability than covariances which are fluctuating due to time changing volatility.

The time evolution of the variance or volatility is often described with GARCH–models Engle95 , which relate the present variance to past variances and past returns. While GARCH–processes are well suited to describe the dynamics of a univariate volatility process Anderson97 , they are less successful in the multivariate setting. In addition to the “curse of dimensionality”, the application of multivariate GARCH–models is hampered by the need to make maximum likelihood estimates of a large number of parameters growing generically like the square of the number of time series. In this note, we suggest two classes of multivariate GARCH–models, which incorporate RMT based progress in the understanding of empirical cross–correlation matrices. These models are both easy to estimate and compare favorably to alternative models as the RiskMetrics covariance estimator Longerstaey+96 .

High dimensional correlation matrices

We estimate the sample covariance matrix

$$\mathrm{\Sigma }_{ij}=\frac{1}{(L1)}\underset{t=1}{\overset{L}{}}(ϵ_{i,t}ϵ_i)(ϵ_{j,t}ϵ_j)$$
(1)
and define the sample correlation matrix $`C_{ij}=\mathrm{\Sigma }_{ij}/(\sigma _i\sigma _j)`$ by normalizing with the standard deviations $`\sigma _i=\sqrt{\mathrm{\Sigma }_{ii}}`$. We calculate cross–correlation matrices for the 118 most frequently traded German companies from a) 1711 daily returns $`ϵ_{i,t}`$ starting 12/01/93 until 08/31/00, and b) 250 daily returns starting 09/01/99 until 08/31/00, diagonalize them, and rank-order their eigenvalues $`\lambda _i<\lambda _k`$ for $`i<k`$. The pdf of the eigenvalues of Ca is displayed in Fig. 1. For both Ca and Cb we find that three eigenvalues lie outside the random matrix interval, but only one of them clearly. The number of eigenvalues clearly outside the RMT–interval is the optimal number of principal components or factors that should be used to forecast the correlation matrix Rosenow2002 .

The largest eigenvalue $`\lambda _{\mathrm{max}}=14.5`$ for the seven–year interval is much larger than $`\lambda _{\mathrm{max}}=7.1`$ for the one–year interval. This confirms the result of Plerou+99 ; drozdz00 ; Plerou+2002  that the strength of market correlations changes in time. Although the use of time series extending over several years in theory increases the accuracy of correlation estimates, it does not capture the effect of time changing correlations and hence is not advisable.

In order to use only the statistically relevant information in the correlation matrix C, we diagonalize C by an orthogonal transformation $`𝖷`$ via $`𝖷^T\mathrm{𝖢𝖷}=\mathsf{\Lambda }`$. We construct a filtered correlation matrix $`𝖢_p^{}`$ by keeping only the $`p`$ largest eigenvalues in the diagonal matrix Gopikrishnan+01 ; Laloux+99 ; Rosenow2002

$$\mathrm{\Lambda }_{ii}^\mathrm{p}=\{\begin{array}{c}\mathrm{\Lambda }_{ii}\mathrm{for}i>Np\hfill \\ 0\mathrm{for}iNp\hfill \end{array},$$
(2)
and by transforming back to the original basis $`𝖢^\mathrm{p}=𝖷\mathsf{\Lambda }^\mathrm{p}𝖷^T`$. For consistency, all diagonal elements of $`𝖢^\mathrm{p}`$ are set to one.

Model description

The basis of the Sliding Correlation Model (SCM) is the decomposition of the covariance matrix into the cross–correlation matrix C and standard deviations $`\{\sigma _i\}`$. Each individual time series is described by a univariate GARCH(1,1) Engle95  process defined by

$$\sigma _{i,t}^2=\alpha _{i,0}+\alpha _{i,1}\sigma _{i,t1}^2+\beta _{i,1}ϵ_{i,t1}^2.$$
(3)
Here, $`\sigma _{i,t}^2`$ and $`ϵ_{i,t}`$ are the volatility and the innovation of time series $`i`$ at time t, respectively. The innovations are modelled as a product $`ϵ_{i,t}=\sigma _{i,t}x_{i,t}`$ with normally distributed $`x𝒩(0,1)`$. The parameters $`\alpha _{i,0}`$, $`\alpha _{i,1}`$, and $`\beta _{i,1}`$ are obtained from maximum likelihood estimations. The cross–correlation matrix C is calculated as a moving average to accommodate the change of correlations in time. We use the filtering method described in the last section and keep only the $`p`$ largest principal components. The covariance matrix in the sliding correlation model SCM(p) of rank $`p`$ is described by

$$\mathrm{\Sigma }_{ij,t}=C_{ij,t}^\mathrm{p}\sigma _{i,t}\sigma _{j,t}.$$
(4)
This model is related to the constant correlation model Bollerslev90 , where a fixed and unfiltered correlation matrix is used in contrast to our $`C_{ij,t}^\mathrm{p}`$.

Factor Models are motivated by the observation that correlations are large in times of high volatility Plerou+2002 . In addition, general multivariate GARCH–models can be simplified to factor models, if the parameter matrices do not have full rank. As a special case with parameter matrices of rank one, a one–factor model was discussed in Engle1990 . We suggest a more general model with $`p`$ factors (abbreviated Factor(p))

$$\mathrm{\Sigma }_{ij,t}=\stackrel{~}{\sigma }_{i,t}^2\delta _{ij}+\underset{k=1}{\overset{p}{}}\lambda _i^{(k)}\lambda _j^{(k)}\mathrm{\Sigma }_t^{(k)}.$$
(5)
Here, $`\stackrel{~}{\sigma }_{i,t}^2`$ is the variance of the residual of time series $`i`$ after linear regression on the $`p`$ factors $`f_t^{(k)}`$ with $`k=1\mathrm{}p`$, and $`\lambda _i^{(k)}`$ is the regression coefficient (factor loading) of factor $`k`$ on time series $`i`$. The dynamics of the residuals is described by a GARCH(1,1)–process, and the dynamics of the factors is similarly given by

$$\mathrm{\Sigma }_t^{(k)}=\alpha _0^{(k)}+\beta _1^{(k)}\mathrm{\Sigma }_{t1}^{(k)}+\alpha _1^{(k)}(f_{t1}^{(k)})^2.$$
(6)
The factors $`f_t^{(k)}`$ are defined via an iterative procedure: we use the largest eigenvector $`𝐮^{(N)}`$ of $`𝖢`$ to define $`f_t^{(1)}=_{i=1}^Nu_i^{(N)}ϵ_{i,t}`$. Next, we perform a linear regression of $`f^{(1)}`$ on the $`ϵ_i`$ and recompute the correlation matrix from the residuals. From the new eigenvector corresponding to the largest eigenvalue, we define a factor $`f^{(2)}`$, perform another linear regression and so on. The number $`p`$ of factors is chosen according to the number of eigenvalues of the correlation matrix outside the RMT–interval, the time interval for the calculation of factors and factor loadings is not moving but fixed.

The data set for parameter estimation comprises seven years of daily data from 12/01/93 until 08/31/00, which are used for the maximum likelihood estimation of all GARCH parameters, for the calculation of factors and factor loadings of the models Factor(p)-seven, whereas only the last year of daily data is used for models Factor(p)-one. Similarly, the cross–correlation matrix of model SCM(p)-seven is calculated over a sliding time window of length seven year, and a window of length one year is used for model SCM(p)-one.

Our benchmark model is the commonly used Risk Metrics covariance estimator Longerstaey+96  defined by

$$\mathrm{\Sigma }_{ij,t}=0.94\mathrm{\Sigma }_{ij,t1}+0.06ϵ_{i,t}ϵ_{j,t}.$$
(7)
Test method

We judge the quality of a multivariate volatility model by its ability to forecast the daily variance of both equally distributed and minimum variance portfolios. The predicted portfolio volatility at time $`t`$ is given by

$$D_{t,\mathrm{predicted}}^2=\underset{i,j=1}{\overset{N}{}}m_{i,t}m_{j,t}\mathrm{\Sigma }_{ij,t}.$$
(8)
Here, $`m_{i,t}`$ is the fraction of the capital invested in stock $`i`$ at time $`t`$. For an equally distributed portfolio, we have $`m_{i,t}1/N`$, hence $`D^2`$ is just the average element of the covariance matrix. This test probes the ability of a model to correctly predict the average covariance between pairs of stock. On the other hand,

$$m_{i,t}=\underset{j=1}{\overset{N}{}}\mathrm{\Sigma }_{ij,t}^1/\underset{k,l=1}{\overset{N}{}}\mathrm{\Sigma }_{kl,t}^1$$
(9)
minimize the variance under the constraint that the total invested money is equal to one. The volatility prediction for minimum variance portfolios is a difficult test for volatility models, as the covariance matrix is first used to calculate the portfolio weights and then to estimate the variance of that portfolio. We compare the predicted volatility to the realized volatility, which is calculated from high frequency (hourly) data Anderson97 . The realized portfolio variance at time $`t`$ is calculated as

$$D_{t,\mathrm{realized}}^2=\underset{i,j=1}{\overset{N}{}}m_{i,t}m_{j,t}\left(g_{i,t}g_{j,t}g_{i,t}g_{j,t}\right)$$
(10)
where the expectation values are taken for 10 hourly returns $`g_i`$ per trading day. We perform an out of sample test from 09/01/00 until 08/31/01. We calculated the average predicted portfolio variance, the average realized variance, the mean square error (MSE) of the prediction and the mean average percentage error (MAPE).

Test results:

For equally distributed portfolios, the performance of the various factor models is very similar to each other. Model Factor(1)-seven, for example, overestimates the overall level of variance but is able to describe the effect of enhanced correlations in high volatility periods (Fig. 2).

The variance estimation of SCM(N)-seven and SCM(1)-seven is very similar to each other, both models severely overestimate the average level of correlations, see Fig. 2 for the results of SCM(1)-seven. The sliding correlation models with a one year time window for correlation estimates provide for the best estimates of realized variances, the RiskMetrics estimator shows comparable results. These results indicate that a good multivariate volatility model must be able to capture not only the change of volatility over time but also the change of the average correlation strength. For details see Table 1.

Minimum variance portfolio: Due to the measurement noise contained in its cross–correlation matrix, the model SCM(N)-one clearly underestimates the variance of minimum variance portfolios, the ratio of realized variance to predicted variance is 2.22. The filtered model SCM(1)-one performs much better in this respect, the ratio is given by 1.32 (Fig. 2d). In addition, the realized variance of the “filtered” model is 27% lower than that of the original model. Adding one more eigenvalue in the cross–correlation matrix does not change this result much. The performance of the various factor models with one or two factors and estimation periods of seven or one year is quite similar to each other, therefore we discuss the representative model Factor(1)-seven (detailed results can be found in Table 2). Although the ratio of realized to predicted variance is with 1.46 better than that of model SCM(N)-one, the realized variance itself is higher than that of the other models. Possibly the performance of the factor models could be increased if the factor loadings were estimated by the maximum likelihood method instead of using the regression coefficients. Details of the results can be found in Table 2. All these results compare favorably to the RiskMetrics covariance estimator, which underestimates the portfolio variance by a factor of 35 (Fig. 2c). Also, the results reported in Engle+01  for DCC-GARCH (overestimation of risk by several hundred percent) are much less precise than the estimates of SCM(1)-one, for example.

In summary, we have shown that the “curse of dimensionality” can be overcome by using appropriate statistical tools relying on RMT. Building on the distinction between statistically relevant information and noise in the correlation matrix, we have suggested two types of multivariate GARCH–models that are easy to estimate even for more than 100 time series. Both the sliding correlation model and the factor model studied have a much better prediction quality than the industry standard RiskMetrics covariance estimator.


format: text

text:
# Decaying neutron propagation in the Galaxy and the Cosmic Ray anisotropy at 1 EeV

## 1 Introduction

The distribution of the arrival directions of high energy Cosmic Rays (CRs) provides a useful tool to investigate their origin. This distribution is found to be very close to isotropic up to the highest energies measured. For energies below the knee ($`E_{knee}3\times 10^{15}`$ eV) the measured anisotropy is below $`0.1\%`$, showing a smooth increase for larger energies and approaching the percent level around $`10^{18}`$ eV . This isotropy is understood as due to the diffusive motion of charged particles in the galactic magnetic field for energies such that $`E/Z10^{17}`$ eV, with $`Z`$ the CR charge. Cosmic rays from galactic sources remain trapped in this field, leading to a flux building up at small energies and to the absence of correlations between the arrival direction of CRs and the source’s direction. At energies above $`10^{19}`$ eV, and up to the highest observed energies, the distribution of the events is compatible with isotropy. However, at these energies CRs are no more trapped in the galactic magnetic field, so that a correlation between the arrival and the source directions is to be expected. The absence of a significant correlation between ultra high energy events and the Galactic plane points then to extragalactic sources as the origin of these events.

In 1999 the AGASA group reported the measurement of the most significant anisotropy detected so far . They presented the results of a search for anisotropies in the arrival directions of cosmic rays with energy above $`10^{17}`$ eV, using data collected over 11 years. They found a strong anisotropy of amplitude $`4\%`$ in a small energy range: $`10^{17.9}10^{18.3}`$ eV. The phase of the first harmonic in right ascension was found to be $`\alpha _1300^{}`$. A two dimensional map showed that this anisotropy can be interpreted as due to excesses of events in two regions of size $`20^{}`$ near the galactic center (GC) and the Cygnus region with significance $`4\sigma `$ and $`3\sigma `$, respectively. A subsequent analysis by the same group, including an independent data set, confirmed the initial claim .

Taking into account that the excess reported by AGASA lies close to the galactic center, that is better seen from the south, the SUGAR group reanalyzed their data . The SUGAR experiment, located near Sidney, has a clear view of the GC, but has much poorer statistics than AGASA. They set a priori the energy range corresponding to the AGASA anisotropy and found a point like excess at $`7.5^{}`$ from the GC, and $`6^{}`$ away from the AGASA maximum. Due to the difference in direction and angular extension of the AGASA and SUGAR signals it is not yet clear whether they are consistent among each other and if they can result from a common origin.

A possible explanation of the anisotropy, already proposed by Hayashida et al. , is that it is due to neutron primary particles. Neutrons of $`10^{18}`$ eV have a Lorentz factor $`\gamma 10^9`$, so that their decay length is about 10 kpc, i.e. similar to the distance from the Earth to the galactic center. Hence, neutrons from a source near the galactic center can propagate up to the Earth before decaying, and they are undeflected by intervening magnetic fields, contrary to the case of charged primaries. The lack of significant anisotropies at smaller energies is explained within this scenario as due to the shorter decay length of the neutrons, and at larger energies by a cutoff in the acceleration energy (or due to a very steep spectrum) of the sources. They suggested that neutrons can be produced by accelerated heavy nuclei which disintegrated due to interactions with ambient matter or photons surrounding the source. Charged particles could stay confined by the magnetic fields present in the acceleration region, while neutrons would escape easily.

Medina-Tanco and Watson  proposed that a more likely way of generating a high energy neutron flux is through high energy protons via interactions with ambient protons or with infra red photons. They found that the neighborhood of the Galaxy’s central supermassive black hole could give the environment needed for these reactions to happen frequently enough. Nagataki and Takahashi showed that the proton-proton process was the most efficient one, and also studied the detectability of the $`\gamma `$ and $`\nu `$ emission associated to the production of neutrons in the GC .

The possibility that the GC is a source of high energy protons was studied by Levinson and Boldt . They showed that, under certain conditions, a compact black hole dynamo associated with the Sgr A\* source may accelerate protons up to energies $`10^{18}`$ eV, with a total luminosity of $`10^{38}`$ erg/s.

In this work, we study in detail the CR arrival direction distribution expected from a source of neutrons in the GC. The fact that the decay length of $`10^{18}`$ eV neutrons is similar to the distance to the proposed source implies that a sizable fraction of the neutrons will actually decay to protons before reaching the Earth, and the trajectories of these secondary protons will be bent by the galactic magnetic fields (GMF). These protons would arrive to the Earth from some preferred directions in the sky: those produced near the Earth would come from directions close to the GC one, while those produced in the inner Galaxy region would arrive preferentially from directions close to the spiral arm, as their trajectories wind around the regular magnetic field lines. Initial work in this direction was performed by Medina-Tanco and Watson  and here we extend it to take into account different allowed magnetic field strengths, we map the arrival direction distribution for different energies and take into account the angular exposure of experiments at different latitudes to study several features of the observations such as the energy dependence of the anisotropies, the compatibility of the observations of AGASA and SUGAR both regarding the angular extent of the observed features and the required source luminosity.

Other suggested explanations for the observed anisotropies include a diffusive flux of charged particles from the GC direction , the drift of heavy nuclei in the galactic magnetic field  and a source of protons in the galactic center .

## 2 Cosmic Ray propagation in the Galaxy

The trajectories of charged CRs will be determined by the galactic magnetic field, which has both a large scale regular component and a random component, both with strengths of a few $`\mu `$G. The magnetic lines of the regular field follow the spiral structure, reversing direction between consecutive arms. According to the preferred model for our galaxy, it is symmetric with respect to the galactic plane. The strength of the field can be modeled as

$`B_{sp}(\rho ,\theta ,z)=B_0{\displaystyle \frac{r_0}{2\rho }}\mathrm{tanh}^3\left({\displaystyle \frac{\rho }{\rho _1}}\right)\mathrm{cos}(\theta \beta \mathrm{ln}{\displaystyle \frac{\rho }{\xi _0}})`$
$`\left({\displaystyle \frac{1}{\mathrm{cosh}(z/z_1)}}+{\displaystyle \frac{1}{\mathrm{cosh}(z/z_2)}}\right),`$ (1)
where $`(\rho ,\theta ,z)`$ are cylindrical coordinates with origin in the galactic center, $`\xi _0=10.55`$ kpc is the value of $`\rho `$ for which the field is maximum in our spiral arm, $`\beta =1/\mathrm{tan}p=5.67`$, where $`p=10^{}`$ is the pitch angle and $`r_0=8.5`$ kpc is the distance from the Sun to the GC. The overall intensity is taken as $`B_0=3\mu `$G, corresponding to a local strength of the regular component close to $`2\mu `$G. The $`\mathrm{tanh}`$ term smoothes the $`1/\rho `$ behavior of the field near the center within a core radius taken as $`\rho _1=2`$ kpc. The vertical scale length associated to the disk field is taken as $`z_1=0.3`$ kpc, while the halo one as $`z_2=4`$ kpc.

The radial and azimuthal components are given by $`B_\rho =B_{sp}\mathrm{sin}p`$ and $`B_\theta =B_{sp}\mathrm{cos}p`$ respectively, while the strength of the field in the $`z`$ direction is assumed to vanish.

Besides the regular magnetic field structure, a random component due to turbulence in the interstellar plasma is known to exist, with the largest eddies having a scale of $`L_{max}100`$ pc. A usual assumption is that the spectrum of the field inhomogeneities is the same as that of the gas density, that is close to a Kolmogorov spectrum. Hence, the turbulent magnetic field can be described by a random field with a power law spectrum such that the magnetic energy density satisfies $`dE/dkk^{5/3}`$. It can then be written as a sum of Fourier modes

$$B_i(\stackrel{}{x})=\frac{d^3k}{(2\pi )^3}B_i(\stackrel{}{k})e^{i(\stackrel{}{k}.\stackrel{}{x}+\varphi _i(\stackrel{}{k}))},$$
(2)
where the phases $`\varphi _i(\stackrel{}{k})`$ are random and the amplitude of the modes is given by

$$B^2(k)=\frac{2}{3}B_{rms}^2k^{5/3}\frac{(2\pi /L_{max})^{2/3}}{1(L_{min}/L_{max})^{2/3}},$$
(3)
for $`2\pi /L_{max}k2\pi /L_{min}`$, and zero otherwise . $`B_{rms}`$ fixes the root mean square amplitude of the random component. In the numerical calculations the integral in Equation (2) was replaced by a sum over 200 modes. Each $`\stackrel{}{B}(\stackrel{}{k}_i)`$ had a random amplitude extracted from a Gaussian distribution with zero mean and root mean square value given by the square root of Equation (3) and the direction was chosen randomly in the plane perpendicular to $`\stackrel{}{k}_i`$, ensuring that $`\stackrel{}{}\stackrel{}{B}(\stackrel{}{x})=0`$. Since particles reaching the Earth at the energies under study never departed from the galactic plane by more than a few hundred pc, the detailed vertical profile of the random component is not relevant for the results, and hence a constant value for it was adopted throughout. Moreover, the random field was also taken as having a constant rms strength in the radial direction.

The effect of these field components on the motion of EeV protons can be estimated as follows. The gyroradius of a high energy CR ($`vc`$) with charge $`Ze`$ and energy $`E`$ traveling through an uniform magnetic field of strength $`B`$ is $`R_{\mathrm{kpc}}1.1E_{\mathrm{EeV}}/ZB_{\mu \mathrm{G}}`$. Hence, a proton with $`E=1`$ EeV in a $`B3\mu \mathrm{G}`$ field has a gyroradius $`R0.4`$ kpc. As the regular GMF is uniform over scales of the order of a few kpc and the thickness of the galactic disc is $``$ 1 kpc, the trajectory of an EeV proton is expected to wind around the magnetic field lines.

The turbulent component adds random deflections to the propagation. The rms deflection of a CR with energy $`E`$ propagating a distance $`L`$ through a turbulent magnetic field of strength $`B_{rms}`$ and maximum inhomogeneity scale $`L_{max}`$ (for $`LL_{max}`$) is given by

$$\delta _{rms}(L)27^{}\frac{\mathrm{EeV}}{E/Z}\frac{B_{rms}}{1\mu \mathrm{G}}\sqrt{\frac{L}{10\mathrm{kpc}}}\sqrt{\frac{L_{max}}{100\mathrm{pc}}}.$$
(4)
Thus, we expect that for EeV protons traveling in the galactic turbulent field, this will induce a significant dispersion in the arrival directions.

## 3 Results

### 3.1 Arrival direction distribution

Neutrons emitted from the GC travel along straight lines until the moment in which they decay into protons, and it is only the trajectories of these last which start to be bent by the GMF. Thus, in order to obtain the distribution of arrival directions at Earth we should follow the trajectories of protons injected in spheres of various radius around the GC and with velocities pointing out of the GC, and record the direction of those arriving to the Earth. These would have to be weighted by the probability of neutron decay at the corresponding radius. However, a more efficient way to compute this distribution is to back track antiprotons leaving the Earth in all different directions and record the points where their velocities point toward the GC. It is clear that if a neutron from the GC were to decay at that point, the proton produced would arrive to the Earth. We actually considered a source with finite radius, and obtained the distance from the GC ($`r`$) and the length of the segment ($`\mathrm{\Delta }r`$) of the proton trajectories for which the velocity pointed towards the source. The corresponding flux of protons from neutron decays that arrive from a given direction is then $`\mathrm{d}J/\mathrm{d}\mathrm{\Omega }e^{r/\gamma c\tau _n}\mathrm{\Delta }r/\gamma c\tau _n`$. Adding finally the contribution from neutrons which did not decay, the arrival distribution map can be constructed.

The left panels in Figure 1 show the arrival direction distribution obtained in this way

when only the regular magnetic field component is included and for different values of the energy between $`10^{17.9}`$ eV and $`10^{18.3}`$ eV. The most intense flux arrives from the GC direction ($`\alpha _{GC}=266.4^{}`$ and $`\delta _{GC}=28.9^{}`$) and corresponds to the direct neutrons. A thin (nearly horizontal) strip extending from the GC towards the north and south galactic poles (i.e. perpendicular to the galactic plane, shown as a dashed line) is due to neutrons decaying near the Earth in our spiral arm and in the nearest inner arm respectively. The more extended excess appearing for right ascensions between $`300^{}`$ and $`360^{}`$ correspond to neutrons that decayed in the inner galactic region, with the protons reaching the Earth through the spiral arm. We can see that the region of the sky where events are spread shrinks considerably with increasing energy, in particular, at $`10^{18.3}`$ eV almost all the events arrive outside the field of vision of AGASA ($`\delta >24^{}`$). Notice that the flux intensity of the plots with different energy should not be compared at this point, as the energy spectrum of the source has not yet been introduced.

The inclusion of the turbulent magnetic field component leads to a spreading in the arrival direction distribution as illustrated in the right panels of Figure 1, in which $`B_{rms}=3\mu `$G was adopted. The spread in the arrival directions increases with the strength of the random component, as indicated by Equation (4), and this is clearly seen in the left panels of Figure 2, which show this distribution for a fixed energy of $`10^{18}`$ eV and increasing intensity of the random field ($`B_{rms}=0,1,2`$ and $`3\mu `$G).

The right ascension density distribution for the $`B_{rms}=3\mu `$G random field maps of Figure 1 is shown in the left panels of Figure 3.

The sharp peak in $`\alpha 265^{}`$ is mainly due to the direct neutrons. This is relatively more important as the energy is increased, due to the increase in the neutron’s decay length. The broader excess associated to the secondary protons arriving along the spiral arm becomes strongly suppressed above $`10^{18}`$ eV due to the increasing difficulty for the regular field to confine the CRs along the spiral arms. This suppression is clearly expected to be more important as the value of $`B_{rms}`$ becomes stronger, while on the other hand a larger value of the regular component $`B_0`$ should lead to a signal more concentrated in the spiral arm direction and surviving up to larger energies.

In order to compare our results with AGASA’s measurements, it is necessary to take into account the *exposure* ($`\omega `$) of the observatory, which gives the time-integrated effective collecting area from each direction of the celestial sphere. The arrival direction distribution measured by an observatory should be proportional to the actual distribution multiplied by $`\omega `$.

The exposure is expected to be homogeneous in right ascension ($`\alpha `$) and to have a declination ($`\delta `$) dependence given by

$$\omega (\delta )\mathrm{cos}(a_0)\mathrm{cos}(\delta )\mathrm{sin}(\alpha _m)+\alpha _m\mathrm{sin}(a_0)\mathrm{sin}(\delta ),$$
(5)
where $`\alpha _m`$ is given by

$`\alpha _m=\{\begin{array}{cc}0\hfill & \text{if }\xi >1\hfill \\ \pi \hfill & \text{if }\xi <1\hfill \\ \mathrm{cos}^1(\xi )\hfill & \text{otherwise}\hfill \end{array}`$ (9)
and

$`\xi {\displaystyle \frac{\mathrm{cos}(\theta _m)\mathrm{sin}(a_0)\mathrm{sin}(\delta )}{\mathrm{cos}(a_0)\mathrm{cos}(\delta )}},`$
where $`a_0`$ is the detector latitude (in the case of AGASA $`a_0=35.78^{}`$) and $`\theta _m`$ is the maximum zenith angle observed. The analysis of AGASA data correspond to events with zenith angle $`\theta <60^{}`$.

Taking into account the exposure, the right ascension distribution expected to be measured at the AGASA latitude looks as shown in the right panels of Figure 3 for a random field amplitude $`B_{rms}=3\mu `$G. The main difference with the full sky results displayed in the left hand panels is the disappearance of the point-like density peak due to the direct neutrons. The right panels of Figure 2 also show this distribution for different values of $`B_{rms}`$ and for $`E=10^{18}`$eV, for the AGASA exposure. The increasing spread of the signal associated to the spiral arm, which is the dominating one in this case, is apparent.

Anisotropies in the CR distribution are often measured performing an harmonic analysis of the right ascension distribution of the events . This takes advantage of the uniform exposure in right ascension of most experiments. The right ascension distribution is related to the differential spectrum of CR, $`\mathrm{d}J/\mathrm{d}\mathrm{\Omega }(E)`$ through

$$F(\alpha )=d\delta \mathrm{cos}\delta \omega (\delta )\frac{\mathrm{d}J}{\mathrm{d}\mathrm{\Omega }}.$$

Expanding it as $`F(\alpha )=F_0+F_1\mathrm{cos}(\alpha \alpha _1)+\mathrm{},`$ the amplitude of the anisotropy is defined as the ratio of the dipole to the monopole intensities: $`\mathrm{\Delta }F_1/F_0`$. In Table 1, the phase of the first harmonic ($`\alpha _1`$)

from the analysis of the right ascension distribution are listed as a function of energy and turbulent magnetic field strength, taking into account AGASA’s exposure, and also for a hypotetical detector in the southern hemisphere (at a latitude $`35^{}`$, similar to the one of AUGER or Sidney).

### 3.2 Amplification

The total flux that we receive from a given source is modified by the presence of a magnetic field. This can be described by introducing the amplification, $`A(E)`$, which is the ratio of the flux of particles arriving to the Earth from any direction in the presence of the GMF, $`J(E)`$, and that in the absence of magnetic field, $`J_0(E)`$, i. e. $`A(E)=J(E)/J_0(E)`$. The amplification is a function of the energy and the mean effect turns out to be larger as the energy decreases, as expected since the CRs tend to be trapped more efficiently in the spiral arms of the GMF. At larger energies the amplification tends to unity, as protons are less affected by the magnetic field and a larger fraction of the neutrons arrive to the Earth before decaying. Figure 4 shows the amplification vs. energy for adopted values of the random magnetic field of $`B_{rms}=0`$ and $`3\mu `$G. The smaller amplification obtained in the presence of the random magnetic field is due to the enhancement of the escape of CRs from the spiral arms produced by this component.

### 3.3 Anisotropy

The anisotropy detected by AGASA arises in a narrow range of energies from $`E=10^{17.9}`$eV to $`E=10^{18.3}`$eV. It has been suggested that the lack of anisotropies at smaller energies could be explained in the model discussed here because of the shorter decay length of the neutrons and the deflection of the resulting protons in the galactic magnetic field. In order to test quantitatively this hypothesis we have to estimate the total anisotropy produced by the GC neutron source in the presence of the background of all the other sources of CR in that energy range.

The flux of CRs with energy around 1 EeV arriving to the Earth is thought to come from galactic sources, such as supernovae, that give the main contribution below the ankle ($`E_{ankle}5\times 10^{18}`$ eV), as well as extragalactic sources, which dominate above the ankle. This can be splitted into the flux from the GC source, $`\mathrm{d}J/\mathrm{d}\mathrm{\Omega }(E)`$, and the contribution from all the other sources, $`\mathrm{d}J_{bkg}/\mathrm{d}\mathrm{\Omega }(E)`$, i.e. $`\mathrm{d}J_{tot}/\mathrm{d}\mathrm{\Omega }=\mathrm{d}J_{bkg}/\mathrm{d}\mathrm{\Omega }+\mathrm{d}J/\mathrm{d}\mathrm{\Omega }`$.

The total (measured) differential energy spectrum of CRs can be approximated by

$$\frac{\mathrm{d}J_{tot}}{\mathrm{d}\mathrm{\Omega }}(E)=C\times \left(\frac{E}{6.3\times 10^{18}}\right)^{3.2\pm 0.05}$$
(10)
in the range of energies from $`4\times 10^{17}\mathrm{eV}<E<6.3\times 10^{18}\mathrm{eV}`$, and $`C=(9.23\pm 0.65)\times 10^{33}\mathrm{m}^2\mathrm{s}^1\mathrm{sr}^1\mathrm{eV}^1`$.

The intensity of the total flux received at the Earth from the source is $`J(E)=J_0(E)A(E)`$, where $`A(E)`$ is the amplification computed in the previous section and for the actual source flux we can consider a power law dependence with energy $`J_0(E)E^{2.2}`$, which is typical of Fermi acceleration in shocks.

If we assume that the main contribution to the anisotropy comes from the GC neutron source, and hence consider that the background CR flux is essentially isotropic, we can write the observed anisotropy as $`\mathrm{\Delta }(E)=\mathrm{\Delta }_s(E)f(E)`$, where $`f(E)F_{0(s)}/F_{0(bkg)}`$ is the ratio of the contributions to the right ascension monopole from the source flux and the background flux. For an anisotropic flux, as in this case, $`f(E)`$ depends on the exposure of the experiment. The anisotropy of the source flux, $`\mathrm{\Delta }_s(E)`$, can be computed from the harmonic analysis of the right ascension distribution of the GC events as discussed in Section 3.1. We show in Figure 5 the expected anisotropy as a function of the energy for different values of the turbulent magnetic field amplitude and exposures.

For a random field amplitude $`B_{rms}=3\mu `$G, the anisotropy grows sharply around $`E=10^{17.8}`$ as observed by AGASA. This is not the case for smaller values of the rms amplitude of the turbulent field as is apparent from the curve labeled $`B_{rms}=1\mu `$G in the Figure, as the CR arrival directions are not spread sufficiently at small energies to wash out the anisotropies. Notice that the effects associated to a harder GC source spectrum are partially compensated by the effects due to the increased amplification of the fluxes for decreasing energies, leading to a more moderate energy dependence. We have normalized the source flux in such a way that the anisotropy is 4$`\%`$ at $`E10^{18}`$eV for the AGASA exposure. This corresponds to a fraction $`f0.026`$ contributed by the GC source to the right ascension monopole at that energy.

The sharp fall in the anisotropy appearing at $`E10^{18.3}`$ eV for the AGASA curve arises because above that energy most CRs arrive as direct neutrons from the GC direction, which is out of the field of view of AGASA. The possible existence of a cutoff in the acceleration energy of the source around that range of values is actually not necessary to explain AGASA’s results.

### 3.4 Source luminosity

We can estimate the luminosity of the GC source which corresponds to the normalization of its flux at the Earth adopted in the previous Section to fit the 4$`\%`$ amplitude of the right ascension first harmonic measured by AGASA. This is given by

$$L_{GC}=4\pi r_0^2dEEJ_0(E),$$
(11)
where $`r_0`$ is the distance to the GC, and the source flux in the absence of magnetic fields is given by $`J_0(E)=J(E)/A(E)`$. The source flux $`J(E)`$ at $`10^{18}`$eV can be related to the background flux $`J_{bkg}`$ at that energy using the value of $`f(E)`$ obtained in the last Section. With the amplification computed in Section 3.2, the local flux at $`10^{18}`$eV results $`J_0(10^{18}\mathrm{eV})1.3\times 10^{30}m^2s^1\mathrm{eV}^1`$. For a power law energy spectrum of the source $`J_0(E)E^{2.2}`$, the total luminosity in the decade of energy from $`10^{17.5}`$ to $`10^{18.5}`$ then results $`L_{GC}(10^{17.5}10^{18.5})4\times 10^{36}`$ erg/s. This value is more than an order of magnitude lower than the maximum luminosity ($`10^{38}`$ erg/s) estimated for Sgr A\*  mentioned in the Introduction.

The expected flux of direct neutrons is a factor $`e^{r_0/\gamma c\tau _n}`$ of the source flux $`J_0(E)`$. The integral of this in the energy range from $`10^{17.9}`$ eV to $`10^{18.5}`$ eV corresponds then to a direct neutron flux of $`2\times 10^{13}\mathrm{m}^2s^1`$, which is not very different from the value $`(9\pm 3)\times 10^{14}\mathrm{m}^2\mathrm{s}^1`$ reported by SUGAR for their point-like excess .

In order to give an idea of the distribution of the excess of events over the ones expected for an isotropic distribution, we made a map (similar to those constructed by AGASA and SUGAR) in which we added to an isotropic background the CR density expected from a GC neutron source of the required intensity to fit the AGASA’s right ascension anisotropy. The map in Figure 6 displays the ratio of the intensity at each direction to the mean density at that declination (following AGASA’s analysis). We have smoothed the map with a window of $`3^{}`$ radius to take into account the typical angular resolution at these energies. Although the particular details of the map will depend on the GMF structure and strength, that are poorly known, we think that the main characteristic of the distribution are picked with the realistic model adopted. Notice however that the proton signal around the GC is sensitive to the actual location of the reversals in direction of the regular magnetic field. In the model adopted, in which the nearest reversal to the GC direction is at less than 1 kpc from the Earth, the excess observed at southern galactic latitudes is more significant than the one at northern latitudes, but this could change if the reversal was at a smaller galactocentric radius. Also, a vertical component of the magnetic field could displace these excesses in the direction parallel to the galactic plane.

## 4 Conclusions

We examined the possibility that a source of neutrons at the GC could be responsible for the anisotropy measured by AGASA and SUGAR. As the GC itself is out of the field of view of AGASA, this observatory would not detect direct neutrons from a source in the GC, but it would just see the protons produced by the neutron decays. As their trajectories are bent by the GMF, the distribution of protons arriving to the Earth depends strongly on the energy of the particles and on the structure and strength of the field, which are not well known. We assumed a particular model for the regular field, which is consistent with the main observational evidences, and considered different reasonable amplitudes for the turbulent magnetic field. We found that there is a sharp rise in the anisotropy at energies around 1 EeV, as observed by AGASA, for the largest amplitude of the turbulent magnetic field considered, $`B_{rms}=3\mu `$G. At energies larger than 2 EeV, essentially no particles are expected from the GC source in the field of view of AGASA, what explains the disappearance of the anisotropy independently of the source acceleration cutoff. The predictions are however very different for an observatory in the south. This one should detect the strong point-like signal from the direct neutrons in addition to the extended signal from the protons. At energies below 1 EeV the neutron signal is suppressed due to the reduced decay length of the neutrons and the expected anisotropy has a similar behavior as that of a northern observatory. At larger energies, and up to the maximum acceleration energy of the source, the point-like signal should be clearly seen. It is interesting to notice that the neutron source luminosity needed to fit the $`4\%`$ right ascension anisotropy detected by AGASA (due to the protons in this model) give rise to a point-like neutron flux consistent with the one measured by SUGAR. However, unless there is some systematic pointing error in their data, the offset of the signal from the GC direction cannot be explained. Moreover, for the GMF adopted, the phase of the first harmonic in right ascension obtained ($`330^{}`$) is somewhat larger than AGASA’s one. Another important point is that within this scenario the excess observed along the Cygnus region is also related to the central source (for this it is crucial that the regular magnetic field be along the spiral arm and not just azimuthal). This excess disappears for $`E>2\times 10^{18}`$ eV due to the impossibility for the GMF to confine the CRs along the spiral arms, independently of the source upper energy cutoff. The alternative explanation of a diffusive flux of charged particles from the GC direction  should give rise to a wide angle anisotropy, with a deficit of particles in the galactic anticenter direction and a cosine angular dependence with respect to the direction of maximum flux, quite different from the neutron source signal. More data from an observatory in the south would be very helpful to confirm or reject these hypothesis. Although the AUGER observatory  in construction in Argentina is focused to the study of CRs with energy above $`10^{19}`$ eV, it may be able to detect lower energy events (down to $`10^{18}`$ eV). A dedicated detector for the $`10^{17}`$ eV to $`10^{19}`$ eV energy range is also under study .

If the sharp peak in the GC direction is seen, confirming the neutron source hypothesis, then the distribution of the protons would give information about the structure of the GMF.

Work supported by CONICET and Fundación Antorchas.

## References


text:
# BOUNDS FOR THE MASS OF THE HEAVIEST RIGHT-HANDED NEUTRINO IN SO(10) THEORIES

DSF-12/2003

## Abstract

By relating the Dirac neutrino mass matrix to the mass of the charged fermions and assuming that the product of the masses of the two lightest left-handed neutrinos is of the order of $`\mathrm{\Delta }m_{sol}^2`$, we derive, within a leptogenesis scenario, a range of values for the mass of the heaviest right-handed neutrino, centered around the scale of $`BL`$ symmetry breaking in the $`SO(10)`$ theory with Pati-Salam intermediate symmetry.

Recently the neutrino oscillation p  solution to the solar neutrino problem has been strongly supported by two experimental results:

a) the rate for the neutral current reaction $`\nu +d\nu +p+n`$ measured at SNO sno , where the number of active neutrinos is in very good agreement with the flux of $`\nu _e`$ coming from the sun predicted by the standard solar model ssm ;

b) the disappearence of $`\overline{\nu }_e`$ found by KamLAND kamland , consistent only with the MSW msw  large mixing solution, which is also the preferred by the combined analysis of the previous experiments prev .

Together with the evidence for the oscillation of atmospheric neutrinos atm  and with the limit by CHOOZ chooz  for the disappearence of $`\overline{\nu }_e`$, it brings to a consistent description of the PMNS lepton mixing matrix p ; pmns

$$U\left(\begin{array}{ccc}\frac{\sqrt{3}}{2}& \frac{1}{2}& 0\\ \frac{1}{2\sqrt{2}}& \frac{\sqrt{3}}{2\sqrt{2}}& \frac{1}{\sqrt{2}}\\ \frac{1}{2\sqrt{2}}& \frac{\sqrt{3}}{2\sqrt{2}}& \frac{1}{\sqrt{2}}\end{array}\right)$$
(1)
and the square mass differences fogli

$$\mathrm{\Delta }m_{atm}^22.610^3\text{eV}^2,$$

$$\mathrm{\Delta }m_{sol}^2710^5\text{eV}^2.$$
(2)
More severe limitations on the sum of neutrino masses are coming from astrophysical observations wmap

$$|m_i|<0.7\text{eV}.$$
(3)
The smallness of neutrino masses with respect to charged fermions may be understood in the framework of the seesaw mechanism ss , which is a strong indication in favour of $`SO(10)`$ unification soio , where the existence of $`\nu _R`$ with very large Majorana masses is expected. These particles are also an essential ingredient of the baryogenesis from leptogenesis scenario fy , for which a lower limit of about $`510^8`$ GeV for the mass of the lightest $`\nu _R`$ has been found di . Indeed, even with a perfect knowledge of the PMNS matrix and of the square mass differences, there are several possibilities for the right handed neutrino mass matrix $`M_R`$, since the Dirac mass matrix $`M_D`$ is unknown, as well as the mass of the lightest $`\nu _L`$ and the relative phases of $`m_i`$.

In $`SO(10)`$, where the left-handed fermions of each family are classified in a single representation (the spinorial 16), we may try to relate $`M_D`$ to the mass matrices of charged fermions. In the past years we have stressed ab  that in a class of $`SO(10)`$-inspired models where one assumes $`M_{R33}0`$, the largest matrix element of $`M_R`$, namely $`M_{R23}`$, takes a value of the order of the scale of $`BL`$ symmetry breaking found in the unified $`SO(10)`$ model b-l . Here we want to consider the limits on the value of the mass of the heaviest $`\nu _R`$ which follow from the hypothesis that, at the scale $`M_Z`$,

$$\text{Det}M_D=\frac{\text{Det}M_e\text{Det}M_u}{\text{Det}M_d},$$
(4)
where $`M_e`$, $`M_u`$ and $`M_d`$ are the mass matrices for charged leptons and quarks. Eqn.(4) is valid with the Higgs doublet belonging to 10 representations, but keeps its reliability for matrices like mm

$$M_u\left(\begin{array}{ccc}0& \sqrt{m_um_c}& 0\\ \sqrt{m_um_c}& m_c& \sqrt{m_um_t}\\ 0& \sqrt{m_um_t}& m_t\end{array}\right),M_d\left(\begin{array}{ccc}0& \sqrt{m_dm_s}& 0\\ \sqrt{m_dm_s}& m_s& \sqrt{m_dm_b}\\ 0& \sqrt{m_dm_b}& m_b\end{array}\right),$$

$$M_D\frac{m_\tau }{m_b}\left(\begin{array}{ccc}0& \sqrt{m_um_c}& 0\\ \sqrt{m_um_c}& 3m_c& \sqrt{m_um_t}\\ 0& \sqrt{m_um_t}& m_t\end{array}\right),M_e\frac{m_\tau }{m_b}\left(\begin{array}{ccc}0& \sqrt{m_dm_s}& 0\\ \sqrt{m_dm_s}& 3m_s& \sqrt{m_dm_b}\\ 0& \sqrt{m_dm_b}& m_b\end{array}\right),$$
(5)
where the factor $`m_\tau /m_b`$ comes from the renormalization from the scale of quark-lepton unification to the electroweak scale and the factor $`3`$ in the 2-2 position, typical of Higgs transforming as a $`(\mathrm{𝟏𝟓},\mathrm{𝟐},\mathrm{𝟐})`$ of $`SU(4)\times SU(2)\times SU(2)`$, realizes the Georgi-Jarlskog mechanism gj . For simplicity we omit the CP violating phase. We assume the normal hierarchy $`|m_3||m_2|,|m_1|`$. From Eqn.(4) and the inverse seesaw formula

$$M_RM_D^TM_L^1M_D,$$
(6)
where $`M_L`$ is the effective neutrino mass matrix, we get

$$|\text{Det}M_R|=\left(\frac{m_em_\mu m_\tau m_um_cm_t}{m_dm_sm_b}\right)^2\frac{1}{|m_1m_2m_3|}=1.610^{30}k\text{GeV}^3,$$
(7)
where $`k`$ is the ratio $`\mathrm{\Delta }m_{sol}^2/|m_1m_2|`$ and we have considered for the quark masses the values at the $`M_Z`$ scale fk , since the ratio $`\text{Det}M_u/\text{Det}M_d`$ is not modified by renormalization. It should be stressed that Eqn.(4) implies a larger value for $`\text{Det}M_D`$ than for $`\text{Det}M_e`$. Moreover, Eqn.(7) implies for the lightest and heaviest right-handed neutrino mass eigenvalues the inequalities

$$M_1<k^{1/3}1.210^{10}<M_3.$$
(8)
The first part is consistent with the lower limit found for $`M_1`$ in the baryogenesis from leptogenesis scenario di ,

$$M_1>510^8\text{GeV},$$
(9)
since we have

$$k>\frac{710^5}{(2.3)^210^2}=1.310^3,$$
(10)
where the upper limit for $`|m_1m_2|`$ comes from Eqn.(3). Then, from Eqns.(8)-(10) we get

$$1.110^9\text{GeV}<k^{1/3}1.210^{10}<M_3<\frac{k1.610^{30}}{2510^{16}}=k6.410^{12}\text{GeV}.$$
(11)
In principle $`k`$ could be very large, since $`|m_1|`$ might be very small, but the fact that

$$\frac{\mathrm{\Delta }m_{sol}^2}{\mathrm{\Delta }m_{atm}^2}2.710^2,$$
(12)
and the large mixing angles in the PMNS matrix, are in favour of not so different neutrino masses. By taking $`k=1`$ the geometrical value for the lower and upper limit for $`M_3`$ is $`2.810^{11}`$ GeV, very near to the value $`2.710^{11}`$ GeV found for the scale of $`BL`$ symmetry breaking in the $`SO(10)`$ gauge theory with the intermediate $`SU(4)\times SU(2)\times SU(2)`$ symmetry b-l . The value for the scale of quark-lepton unification, smaller than the mass of the lepto-quarks mediating proton decay, improves the prediction for $`m_b/m_\tau `$ fk . Conversely, by taking that value for $`M_3`$ and the inequality (9) for $`M_{1,2}`$, one finds

$$|m_1m_2|<(310^2)^2\text{eV}^2$$
(13)
which implies

$$|m_i|<0.12\text{eV}.$$
(14)
In conclusion, the lower limit on $`M_1`$ in the baryogenesis via leptogenesis scenario and the assumption (7), which is rather reliable in the framework of $`SO(10)`$, and $`k=1`$, which is also reliable, since $`\sqrt{\mathrm{\Delta }m_{sol}^2}`$ is less than one order smaller than $`\sqrt{\mathrm{\Delta }m_{atm}^2}`$, estabilish a range for $`M_3`$ centered around the value of the $`BL`$ breaking scale in the $`SO(10)`$ model with Pati-Salam intermediate symmetry ps . Our conclusion depends also on the assumption (4), which is not necessary true, and from the hypothesis $`k1`$, which is reasonable. With smaller $`|m_1|`$ the lower limit found, which goes as $`|m_1|^{1/3}`$, has a slow dependance on $`|m_1|`$, while the upper limit, going as $`|m_1|^1`$, varies more rapidly.

Finally, the result we found, which supports the $`SO(10)`$ theory with Pati-Salam intermediate symmetry, does not exclude SUSY extensions, where lower values of $`M_R`$ than the scale of $`BL`$ breaking, which is higher in that case, may be obtained from non-renormalizable couplings nr  $`h\varphi _{\mathrm{𝟏𝟔}}\varphi _{\mathrm{𝟏𝟔}}\overline{\psi }_R^c\psi _L/M_P+\text{h.c.}`$, dumped by the Planck scale.


text:
# COMMENT ON “PHASE TRANSITION-LIKE BEHAVIOR IN A LOW-PASS FILTER”

Krivine and Lesne use an example taken from The Feynman Lectures in an attempt to illustrate that “many interesting physical properties can however be missed because of the improper use of mathematical techniques”. The supposedly incorrect mathematical procedure has to do with an ordering of limits. An infinite series that is convergent in the presence of a small parameter no longer converges when the parameter is set to zero before the series is summed.

The authors, correctly in my view, emphasize the physical importance of distinguishing between infinite systems and large finite systems. In their example discontinuities in certain physical quantities only exist (mathematically) for infinite systems.

I suggest, however, that the authors have demonstrated a different mathematical point than the one that they propose: infinite series live a life of their own and need not be constrained to be the limit of sequences of finite series. This point was made long ago by Borel and was probably known to Abel and Cauchy. I emphasize the point with an example of an infinite series of resistive elements that sum to a negative resistance. The infinite series represents different physics from any of the possible finite series.

Let $`\{R_i\}`$ be a set of resistors, each having resistance $`R_i=p^iRp>1`$ and $`R`$ an arbitrary resistance value. Then $`Z_n=_{i=0}^nR_i`$ is the resistance of a set of such resistors connected in series, and the value of $`Z_n`$ grows without bound as $`n`$ increases. Clearly, a quantity $`Z`$ defined by $`Z_{n=0}^{\mathrm{}}R_n`$ makes no sense as a limit of a convergent sequence of finite sums.

We may, however, emulate Feynman and define $`Z`$ from the recursive relation $`ZR=pZ`$ which follows from the definition of $`Z`$ and the fact that an infinite series less a finite set of its members is still an infinite series111The series, in fact, satisfies Hardy’s criteria of $`𝔛`$-summability. Solving the last equation for $`Z`$ leads to the result that

$$Z=R/(p1)$$
(1)
a negative resistance.

Feynman also shows us how to build such infinite-series resistors. One simply terminates a finite-series resistor having resistance $`Z_n`$ with a negative resistance having resistance $`p^{n+1}R/(p1)`$. Each such resistor will then have negative resistance $`Z`$.

When the quantity $`p`$ has values $`p<1`$, there is no difference between the limit of a sequence $`Z_n`$ of increasing $`n`$ and the value $`Z`$ obtained in Eq. 1. This does not mean that Eq. 1 is wrong, as the authors of Ref. 2 seem to imply. It does mean that the infinite sum involved represents two different physical situations when $`p<1`$ and when $`p>1`$, involving, respectively passive and active circuit elements.

This Comment is intended, however, to emphasize the mathematical fact that infinite (and infinitesimal) mathematical operations may be justified independently of arguments involving limits.

I am indebted to Cosmas Zachos for bringing the Borel reference to my attention. This work was supported by the U.S. Department of Energy, Division of High Energy Physics, Contract W-31-109-ENG-38.


text:
# 1 INTRODUCTION

## 1 INTRODUCTION

The quark level process $`bs\nu \overline{\nu }`$ represents a rare flavour-changing neutral-current (FCNC) decay which proceeds at the one-loop level in the Standard Model (SM) via “penguin” and “box” diagrams such as those shown in Fig. 1. The inclusive $`bs\nu \overline{\nu }`$ process is nearly free from theoretical uncertainties associated with strong interaction effects, permitting a fairly precise prediction of the SM branching fraction. The inclusive branching fraction, summed over the three neutrino flavours, is estimated to be $`(4.1_{1.0}^{+0.8})\times 10^5`$ . Since additional heavy particles would also contribute additional loop diagrams, various “New Physics” scenarios can potentially lead to significant enhancements to the SM branching fraction . Unfortunately, an experimental search for the inclusive $`bs\nu \overline{\nu }`$ process is extremely difficult in a $`B`$-factory environment due to the presence of two unobserved neutrinos which limit the available kinematic constraints that can be exploited in order to suppress other $`B`$ decay backgrounds.

Instead, we search for the exclusive $`B^{}K^{}\nu \overline{\nu }`$ decay mode, which proceeds via the $`bs\nu \overline{\nu }`$ process. The SM branching fraction for $`B^{}K^{}\nu \overline{\nu }`$ is estimated to be $`(B^{}K^{}\nu \overline{\nu })4\times 10^6`$ . The best published limit on the exclusive branching fraction is from CLEO  with a limit of $`(B^{}K^{}\nu \overline{\nu })2.4\times 10^4`$ at the $`90\%`$ confidence level. BABAR has already reported a preliminary upper limit $`(B^{}K^{}\nu \overline{\nu })9.4\times 10^5`$  based on $`50.7`$ fb-1 data. The two BABAR analyses use reconstruction methods which produce mutually exclusive data samples, permitting the two statistically independent results to be combined to obtain an improved limit.

## 2 THE BABAR DETECTOR AND DATASET

The data used in this analysis were collected with the BABAR detector at the PEP-II storage ring during the period from 2000 – 2002, corresponding to a total integrated luminosity of $`80.06`$ fb-1 collected at the $`\mathrm{{\rm Y}}(4S)`$ resonance. This “onpeak” sample is estimated to contain $`(86.9\pm 1.0)\times 10^6`$ $`B\overline{B}`$ pairs. This data set is supplemented by a sample of $`9.6`$ fb-1 of “offpeak” data collected approximately $`40`$ MeV below the $`\mathrm{{\rm Y}}(4S)`$ resonance, which is used to study continuum background sources due to $`e^+e^{}f\overline{f}`$ where $`f=u,d,s,c,e,\mu ,\tau `$.

The BABAR detector is an hermetic detector optimized to provide precision vertexing, charged and neutral particle reconstruction and particle identification in an asymmetric $`B`$-factory environment. Tracking is provided by a five-layer double-sided silicon vertex tracker (SVT), surrounded by a $`40`$-layer drift chamber (DCH) filled with a mixture of helium and isobutane. The SVT and DCH are situated within a $`1.5`$ T solenoidal field. $`K`$/$`\pi `$ separation is provided by a quartz ring-imaging Cherenkov detector (DIRC) located immediately outside of the DCH. The electromagnetic calorimeter (EMC) is used to measure energy and position of photons and electrons. Muon identification is achieved through segmentation and instrumentation of the iron of the magnetic flux return (IFR) using resistive plate chambers. A more detailed description of the BABAR detector can be found in .

A GEANT4  based Monte Carlo (MC) simulation is used to model the signal efficiency and physics backgrounds. MC samples equivalent to approximately three times the data luminosity were used to model $`B\overline{B}`$ events, and samples equivalent to approximately 1.5 times the data luminosity were used to model continuum events.

## 3 ANALYSIS METHOD

Due to the presence of two unobserved neutrinos, the $`B^{}K^{}\nu \overline{\nu }`$ decay mode 777Charge conjugate modes are implied throughout this paper, however the signal mode will always be denoted as a $`B^{}`$ decay, while the fully reconstructed tag $`B`$ will be denoted as a $`B^+`$ decay to avoid confusion. lacks the kinematic constraints which are usually exploited in $`B`$ decay searches in order to reject both continuum and $`B\overline{B}`$ backgrounds. Consequently, the approach which is used in this analysis is to first reconstruct the accompanying “tag” $`B^+`$, which is produced in association with the signal $`B^{}`$ through the process $`\mathrm{{\rm Y}}(4S)B^+B^{}`$, and then search for evidence of a $`B^{}K^{}\nu \overline{\nu }`$ decay among the tracks and clusters not associated with the reconstructed tag $`B`$. In order to avoid experimenter bias, the signal region in data is not examined (“blinded”) until the cuts are finalized.

We reconstruct the tag $`B^+`$ in a set of decay modes $`B^+\overline{D}^0X^+`$ where $`X^+`$ is a hadronic system composed of up to three charged mesons (either $`\pi `$ or $`K`$) and up to two $`\pi ^0`$ candidates. The $`\overline{D}^0`$ candidate is reconstructed in one of the three decay modes $`\overline{D}^0K^+\pi ^{}`$, $`\overline{D}^0K^+\pi ^{}\pi ^0`$ or $`\overline{D}^0K^+\pi ^{}\pi ^+\pi ^{}`$ . Candidate $`B^+`$ decays are identified by combining $`\overline{D}^0`$ candidates with sets of charged tracks and $`\pi ^0`$ candidates until the combination yields a $`B`$ candidate consistent with the kinematics expected for a true $`B`$ meson decay. We use the two kinematic variables $`m_{\mathrm{ES}}`$ and $`\mathrm{\Delta }E`$, defined by $`m_{ES}\sqrt{E_{\mathrm{beam}}^{}{}_{}{}^{2}\stackrel{}{p}_{\mathrm{B}}^{}{}_{}{}^{2}}`$ where $`\stackrel{}{p}_\mathrm{B}`$ is the momentum vector of the $`B`$ candidate and $`E_{\mathrm{beam}}`$ is the beam energy, and $`\mathrm{\Delta }EE_BE_{\mathrm{beam}}`$, where $`E_B`$ is the energy of the $`B`$ candidate. All these quantities are evaluated in the center of mass (CM) frame. If multiple $`B`$ candidates are identified within the kinematic region $`m_{\mathrm{ES}}>5.2`$ GeV$`/c^2`$ and $`1.8\mathrm{GeV}<\mathrm{\Delta }E<0.6`$ GeV, only the candidate for which $`\mathrm{\Delta }E`$ is closest to zero is retained.

Combinatorial backgrounds from continuum processes are significantly reduced by requiring the thrust, computed using all tracks and clusters in the event, be less than 0.925 and that $`|\mathrm{cos}\theta _T|`$, the magnitude of the cosine of the angle between the thrust axis defined by tracks and clusters used to reconstruct the tag $`B`$ candidate and the thrust axis defined by all other tracks and clusters in the event, be less than 0.8. Correctly reconstructed $`B`$ meson candidates produce a peak in the $`m_{\mathrm{ES}}`$ distribution above a combinatorial background at the nominal $`B`$ mass as shown in Fig. 2. The tag $`B`$ yield is determined directly from data by determining the peaking component of the $`m_{ES}`$ distribution. Tag $`B`$ candidates which are reconstructed in MC events containing a true $`B^{}K^{}\nu \overline{\nu }`$ decay are found to possess very little combinatoric background, since there are few additional tracks and clusters in the event which can be randomly combined to produce combinatoric tag $`B`$ candidates. Events with tag $`B`$ candidates lying within a signal region defined by $`5.272\text{ GeV}/c^2\text{ }<m_{\mathrm{ES}}<5.288\text{ GeV}/c^2\text{ }`$ are retained for use in the search for $`B^{}K^{}\nu \overline{\nu }`$ decays. Events in the region $`5.225\text{ GeV}/c^2\text{ }<m_{\mathrm{ES}}<5.265\text{ GeV}/c^2\text{ }`$ are retained for use in background studies as discussed below. A discrepancy of approximately $`25\%`$ between the observed yield in data and the predicted yield in MC is corrected by scaling the peaking component of the MC simulation. Good agreement between data and MC is obtained once this correction is applied. The same scale factor, $`(75\pm 7)\%`$, is used to correct the tag $`B`$ yield in signal MC as well as MC estimates of peaking backgrounds. The quoted error of $`7\%`$ is due to the uncertainty in the estimation of the peaking component of the data and due to the discrepancy between MC and data in the shape of the non-peaking background component.

Once a reconstructed tag $`B`$ candidate has been identified, $`B^{}K^{}\nu \overline{\nu }`$ signal candidates are selected by considering all tracks and clusters in the event which are not used in the tag $`B`$ reconstruction. This set of tracks and clusters is referred to in the following as the “signal-side” of the event. Candidate events are required to possess exactly one signal-side reconstructed charged track with a charge which is opposite that of the tag $`B`$. The signal-side charged track multiplicity is plotted for signal MC and data in Fig. 3. The signal candidate track is required to satisfy particle identification criteria for charged kaons based on information from the tracking system and the DIRC. The kaon candidate is boosted into the CM frame assuming a kaon mass hypothesis, and the CM momentum, $`p_K^{}`$, is required to be greater than $`1.5`$ GeV$`/c`$. The average particle identification efficiency in the momentum range of interest is $`85\%`$ and the typical $`\pi K`$ misidentification rate is $`2\%`$. The $`p_K^{}`$ of signal candidate tracks is plotted in Fig. 4. We assume the kaon momentum spectrum is described by the decay model of reference  and correct the signal MC distribution (which is generated with a phase-space model) accordingly.

Constraints are also imposed on signal-side EMC clusters to reject events with significant neutral energy deposition in the calorimeter. $`B^{}K^{}\nu \overline{\nu }`$ events possess an average of two additional signal-side clusters. These are generally attributable to hadronic split-offs in the EMC, usually from pions or kaons associated with the tag $`B`$ side of the event, or to beam related backgrounds. Events possessing one or more $`\pi ^0`$ candidates, composed of combinations of two EMC clusters with lab frame energy greater than $`30`$ MeV which combine to produce an invariant mass in the range $`122\text{ MeV}/c^2<m_{\gamma \gamma }<145\text{ MeV}/c^2`$, are rejected. In addition to this $`\pi ^0`$ veto, a restriction is imposed on the total “extra” signal-side neutral energy, $`E_{\mathrm{extra}}`$, that is present in the event. $`E_{\mathrm{extra}}`$ is computed by summing the CM-frame energies of all signal-side EMC clusters with lab frame energy greater than $`30`$ MeV. Signal candidate events are required to possess $`E_{\mathrm{extra}}<300`$ MeV. The $`E_{\mathrm{extra}}`$ distribution obtained from signal MC is compared to data in Fig. 5.

The requirement that signal candidate events possess low signal-side track multiplicity and relatively little additional neutral energy in the EMC also tends to select $`B^+B^{}`$ events in which the tag $`B^+`$ has been correctly reconstructed, but the opposing $`B^{}`$ has one or more unreconstructed particles which have passed outside the detector acceptance in either the forward or backward directions. An additional cut is therefore imposed on the direction of the missing momentum vector in the CM-frame, $`p_{\mathrm{miss}}^{}`$, which is required to satisfy $`|\mathrm{cos}\theta _{p_{\mathrm{miss}}^{}}|<0.8`$.

Due to the low signal-side multiplicity, there is almost no tag $`B`$ reconstruction mode-dependence of the measured signal-side efficiency. The overall selection efficiency, $`ϵ_{\mathrm{tot}}`$, can therefore be expressed as the product of the tag $`B`$ selection efficiency, $`ϵ_{\mathrm{tag}}`$, and the signal-side efficiency $`ϵ_{\mathrm{sig}}`$. The overall efficiency is estimated to be $`ϵ_{\mathrm{tot}}=(0.046\pm 0.005)\%`$, while the $`ϵ_{\mathrm{sig}}`$ is estimated to be $`(35\pm 1)\%`$. The uncertainties on the efficiencies are due to both statistics and systematics.

Monte Carlo modeling of the signal efficiency and background estimates is validated by comparing the yields obtained in a number of data control samples. These samples include $`9.6`$ fb-1 of offpeak data, an $`m_{\mathrm{ES}}`$ sideband region defined by $`5.225\text{ GeV}/c^2\text{ }<m_{\mathrm{ES}}<5.265\text{ GeV}/c^2\text{ }`$, an $`E_{\mathrm{extra}}`$ sideband defined by $`0.5\text{ GeV }<E_{\mathrm{extra}}<1.5\text{ GeV}`$, and a “large $`m_{\mathrm{ES}}`$” sideband spanning the entire region defined by $`5.225\text{ GeV}/c^2\text{ }<m_{\mathrm{ES}}<5.265\text{ GeV}/c^2\text{ }`$ and $`E_{\mathrm{extra}}<1.5`$ GeV. In addition, we retain samples of events which pass all of the nominal signal selection requirements except that they are required to have a total of two or three charged tracks associated with the signal side of the event instead of only one. All data control samples were found to be in good agreement with the MC predictions, as shown in Table 1.

Backgrounds consist primarily of $`B^+B^{}`$ events in which the tag $`B^+`$ has been correctly reconstructed but in which the accompanying $`B^{}`$ decays to a high-momentum kaon and additional particles which are not reconstructed by the tracking detectors or calorimeter. Typically these events contain one or more $`K_L^0`$ and/or neutrinos, and frequently also one or more additional charged or neutral particles which pass outside of the tracking or calorimeter acceptance. This “peaking” background component is evaluated directly from $`B^+B^{}`$ MC and estimated to yield $`1.7\pm 0.6`$ events in $`80.06`$ fb-1 of data. This estimate is validated by comparison with the data in the $`E_{\mathrm{extra}}`$, $`N_{\mathrm{tracks}}=2`$ and $`N_{\mathrm{tracks}}=3`$ sideband, all of which have large contributions from peaking $`B^+B^{}`$ backgrounds. Due to the limited MC statistics, a smaller combinatorial component of the background is estimated by scaling the observed yield in the MC $`m_{\mathrm{ES}}`$ sideband into the signal region. This scaling assumes an “Argus function” shape for the $`m_{\mathrm{ES}}`$ distribution of the combinatorial component which is obtained from data using a “wrong-sign” tag $`B`$ sample in which the charge of the $`D^0`$ daughter kaon is inconsistent with the charge of the tag $`B`$ candidate, resulting in no significant peaking component in the $`m_{\mathrm{ES}}`$ distribution. Scaling the $`m_{\mathrm{ES}}`$ sideband into the signal region yields an additional background of $`1.0\pm 0.4`$ events, leading to an estimated total background of $`2.7\pm 0.7`$, where the quoted uncertainties are due to MC statistics.

## 4 SYSTEMATIC STUDIES

Estimates of systematic uncertainties are summarized in Table 2. Systematic uncertainties in the branching ratio determination are dominated by the statistical uncertainty on the background estimate, and by the uncertainty on the determination of the tag $`B`$ yield in data using $`(75\pm 7)\%`$ scaling correction determined from the comparison of $`B^+\overline{D}^0X^+`$ yields in data and MC. Since the tag $`B`$ yield is determined directly from data, no other systematic uncertainties associated with the tag $`B`$ need to be assigned. The signal track reconstruction and kaon identification procedure yield comparatively small uncertainties on the efficiency and background estimates. Uncertainties also arise from the MC modeling of the energy and multiplicity of low energy clusters in the EMC which would potentially produce a bias in the $`E_{\mathrm{extra}}`$ distribution. These uncertainties are estimated by evaluating the change in the MC efficiency and background estimates when low energy clusters are selectively removed from the MC until data and MC multiplicity distributions are in agreement. An additional systematic uncertainty is assigned to the efficiency estimate as a result of the re-weighting procedure used to correct the MC kaon momentum spectrum to be consistent with that predicted by reference . Note however that this uncertainty does not account for variations in the momentum spectrum which would result from the use of other theoretical models. A small uncertainty also enters the branching ratio limit calculation from the estimation of the number of $`B^+B^{}`$ events present in the data sample.

## 5 PHYSICS RESULTS

Unblinding the analysis revealed a total of $`N_{\mathrm{sel}}=\text{ 3}`$ events in the signal region, with an expected background, $`N_{\mathrm{bg}}`$, of $`2.7\pm 0.8`$ where the additional systematic uncertainties from table 2 have been combined with the MC statistical uncertainty. The $`m_{\mathrm{ES}}`$ distribution and signal kaon candidate momentum spectrum are shown in Fig. 6 and Fig. 7 respectively. The $`B^{}K^{}\nu \overline{\nu }`$ branching ratio is computed as follows:

$$(B^{}K^{}\nu \overline{\nu })=\frac{(N_{\mathrm{sel}}N_{\mathrm{bg}})}{2N_{B^+B^{}}ϵ_{\mathrm{tot}}}$$
(1)
where $`ϵ_{\mathrm{tot}}=(0.046\pm 0.005)\%`$ is the overall signal selection efficiency, and $`N_{B^+B^{}}`$ is the number of $`\mathrm{{\rm Y}}(4S)B^+B^{}`$ events in the data. $`N_{B^+B^{}}=(43.5\pm 0.5)\times 10^6`$ is obtained by assuming equal branching fractions for $`\mathrm{{\rm Y}}(4S)`$ decays into charged and neutral $`B`$ mesons. The central value of the branching ratio is determined to be  $`(B^{}K^{}\nu \overline{\nu })=(0.8\pm 2.0)\times 10^5`$ where the quoted uncertainty is from systematics only. The significance of the central value is somewhat less than $`1\sigma `$ and we therefore quote a $`90\%`$ confidence level limit. The branching ratio limit is computed using a frequentist approach based on reference . The confidence level for a given branching fraction limit “guess” is obtained by generating a large number of experiments in which the systematic uncertainties in the inputs to equation 1 are modeled by Gaussian distributions and the signal statistics are modeled by a Poisson distribution. The limit is set at the value of the branching fraction at which $`10\%`$ of the generated experiments produce a yield which is less than the observed data yield of three events. This procedure results in a limit of  $`(B^{}K^{}\nu \overline{\nu })<1.05\times 10^4`$ at the $`90\%`$ confidence level assuming the central model of reference . This limit is weakly model dependent, since the signal efficiency depends on the kaon momentum spectrum. Varying the model within the range specified in  results in a variation in the limit in the range $`(1.021.10)\times 10^4`$. Alternatively, the model of reference  yields a limit of $`9.5\times 10^5`$. Sensitivity to new physics will also be model dependent, however it should be noted that this analysis has an efficiency which is relatively uniform as a function of kaon momentum above the $`1.5`$ GeV$`/c`$ cut. The efficiency is zero for kaon momenta below $`1.5`$ GeV$`/c`$, which may effect the limit interpretation for exotic New Physics modes with significantly different kaon momentum spectra.

This analysis can be combined with the result of a preliminary BABAR search for $`B^{}K^{}\nu \overline{\nu }`$ reported previously , based on a sample of $`50.7`$ fb-1 of data. This analysis yielded a limit of $`(B^{}K^{}\nu \overline{\nu })9.4\times 10^5`$ at the $`90\%`$ confidence level assuming the model of reference . In this analysis, the tag $`B`$ was reconstructed in a set of semileptonic $`B`$ decay modes of the form $`B^+\overline{D}^0\mathrm{}^+\nu X^0`$ were $`X^0`$ can be either nothing, or one or more photons which are consistent with the decay products of higher mass open charm states such as $`\overline{D}^0\overline{D}^0\gamma /\pi ^0`$. Two events were observed, and were treated as signal for the limit determination (i.e. no background subtraction was performed). Because this previous analysis required an identified high momentum lepton for the tag $`B`$ reconstruction, while in the present work the tag $`B`$ modes are purely hadronic, the two analyses are by construction statistically independent and can be readily combined. Using the same frequentist method as described above and ignoring correlated systematic uncertainties which are expected to be small, we obtain a combined limit of  $`(B^{}K^{}\nu \overline{\nu })<7.0\times 10^5`$ at the $`90\%`$ confidence level.

## 6 SUMMARY

We have performed a search for the rare FCNC decay $`B^{}K^{}\nu \overline{\nu }`$ using a method in which the accompanying tag $`B`$ meson is reconstructed into a set of hadronic final states. Using a data sample corresponding to an integrated luminosity of $`80.06`$ fb-1, we observe a total of three signal candidate events, consistent with the background expectation of $`2.7\pm 0.8`$. We determine a preliminary limit of the branching fraction  $`(B^{}K^{}\nu \overline{\nu })<1.05\times 10^4`$ at the $`90\%`$ confidence level. Combining our result with the result from a previous and independent BABAR search for $`B^{}K^{}\nu \overline{\nu }`$ yields a combined preliminary limit of  $`(B^{}K^{}\nu \overline{\nu })<7.0\times 10^5`$. The result is consistent with the SM expectation of $`(B^{}K^{}\nu \overline{\nu })4\times 10^6`$.

## 7 ACKNOWLEDGMENTS

We are grateful for the extraordinary contributions of our PEP-II colleagues in achieving the excellent luminosity and machine conditions that have made this work possible. The success of this project also relies critically on the expertise and dedication of the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and the kind hospitality extended to them. This work is supported by the US Department of Energy and National Science Foundation, the Natural Sciences and Engineering Research Council (Canada), Institute of High Energy Physics (China), the Commissariat à l’Energie Atomique and Institut National de Physique Nucléaire et de Physique des Particules (France), the Bundesministerium für Bildung und Forschung and Deutsche Forschungsgemeinschaft (Germany), the Istituto Nazionale di Fisica Nucleare (Italy), the Foundation for Fundamental Research on Matter (The Netherlands), the Research Council of Norway, the Ministry of Science and Technology of the Russian Federation, and the Particle Physics and Astronomy Research Council (United Kingdom). Individuals have received support from the A. P. Sloan Foundation, the Research Corporation, and the Alexander von Humboldt Foundation.


text:
# Tuning Heavy Fermion Systems into Quantum Criticality by Magnetic Field

## I Introduction

Quantum critical points (QCPs) are of extensive current interest to the physics of correlated electrons, as the proximity to a QCP provides a route towards non-Fermi liquid (NFL) behavior. While a broad range of correlated electron materials are being studied in this context, heavy fermions (HF) metals have been playing an especially important role since a growing list of them have been found to explicitly display a magnetic QCP Stewart . HF metals contain a dense lattice of certain lanthanide ($`4f`$) or actinide ($`5f`$) ions which are, at sufficiently low temperatures ($`TT_K`$, $`T_K`$: Kondo temperature), strongly coupled by the conduction electrons, yielding the formation of heavy quasiparticles (QP). They are highly suited to study quantum-critical behavior since they can be tuned continuously from an antiferromagnetic (AF) to a paramagnetic metallic state by the variation of a single parameter, i.e. the strength of the $`4f`$-conduction electron hybridization $`g`$, which can be modified by the application of either external pressure or chemical substitution.
The origin of non-Fermi liquid (NFL) behavior in HF systems has been studied intensively in the past decade but is still unclear up to now Stewart . In particular two different scenarios are discussed for the QCP, where long-range AF order emerges from the HF state; one in which NFL behavior arises from Bragg diffraction of the electrons off a critical spin-density wave (SDW) Hertz ; Millis ; Moriya , the other in which the bound-state structure of the composite heavy fermions breaks up at the QCP resulting in a collapse of the effective Fermi temperature Coleman ; Si . In the SDW scenario, assuming threedimensional (3D) spinfluctuations, singular scattering occurs only along certain ”hot lines” connected by the vector q of the near AF order while the remaining Fermi surface still behaves as a Landau Fermi liquid (LFL). Therefore, the low-temperature specific heat coefficient $`C(T)/T`$ that measures the QP mass is expected to show an anomalous temperature dependence $`C(T)/T=\gamma _0\alpha \sqrt{T}`$, but remains finite at the QCP Moriya . A diverging QP mass, as evident from the $`C(T)/T\mathrm{log}(T_0/T)`$ behavior found, e.g. in the prototypical system CeCu6-xAux for $`x_c=0.1`$ (Ref. Loehneysen ), would arise only if truly 2D critical spinfluctuations render the entire Fermi surface ”hot” Rosch . On the other hand, measurements of the inelastic neutron scattering on CeCu5.9Au0.1 Schroeder  showed that the critical component of the spin fluctuations is almost momentum independent leading to the proposal of the locally critical scenario Si ; Schroeder . Since $`T`$-dependent measurements at the QCP alone provide no information on how the heavy quasiparticles decay into the quantum critical state it is necessary to tune the system away from the magnetic instability in the LFL state and to follow the QP properties upon approaching the QCP. We will show that magnetic fields can be used for this purpose.
The aim of this article is to review a number of experiments on the HF systems CeNi2Ge2 Gegenwart CNG  (section II) and YbRh2Si2 Trovarelli Letter  (section III) that both crystallize in the tetragonal ThCr2Si2 structure. They are ideally suited to study AF QCPs since they are located very near to the magnetic instability, and since the effect of disorder is minimized in these high quality single crystals with very low residual resistivities. We focus on low temperature specific heat, $`C`$, and electrical resistivity, $`\rho `$, measurements at various magnetic fields that allow to follow the QP mass and QP-QP scattering cross section upon approaching the QCP from the field-induced FL side. Furthermore, the high field dc-magnetization $`M(B)`$ has been measured for both systems in order to investigate the suppression of the Kondo effect by the polarization of the magnetic moments.

## II $`\text{CeNi}_2\text{Ge}_2`$

As first reported by Knopp et al. Knopp , CeNi2Ge2 crystallizing in the tetragonal ThCr2Si2 structure is one of the very few HF systems exhibiting a non-magnetically ordered non-superconducting ground state. From both the low temperature specific heat coefficient of about 0.35 JK-2mol-1 and the flat minimum at $`T30`$ K in the temperature dependence of the quasielastic line width $`\mathrm{\Gamma }(T)`$, measured by inelastic neutron scattering, a Kondo temperature of 30 K has been deduced Knopp . Furthermore, a broad maximum is observed at 30 K for the magnetic susceptibility $`\chi (T)`$ measured along the magnetic easy direction parallel to the crystallographic $`c`$-axis, resembling the prototypical Kondo-lattice system CeRu2Si2 ($`T_K=10`$ K) Haen . In contrast to CeRu2Si2, for CeNi2Ge2 the closer inspection of the low temperature thermodynamic and transport properties, as will be discussed below, revealed pronounced NFL behavior at temperatures well below $`T_K`$ Gegenwart CNG . The NFL effects are clearly related to the very near magnetic instability: increasing the lattice parameter by isoelectronic substitution of the Ni by the larger Pd atomes in Ce(Ni1-xPdx)2Ge2 induces long-range AF order below $`T_N=2`$ K for $`x=0.2`$ (Ref Knebel ). The extrapolation of $`T_N(x)`$ towards zero temperature suggests a very small critical concentration $`x_c`$.

The low temperature specific heat coefficient measured with a high-quality polycrystalline CeNi2Ge2 sample Gegenwart CNG  is plotted in Fig. 1 on a logarithmic temperature scale. Instead of being constant, as expected for a FL, $`C(T)/T`$ strongly increases upon cooling from 6 K to 0.4 K. Similar behavior has been reported in Refs. Knopp ; Aoki ; Koerner ; Steglich2000 ; Cichorek . Below 2.5 K, this increase can be well described by $`C(T)/T=\gamma _0c\sqrt{T}`$ with $`\gamma _0=0.43`$ JK-2mol-1 (see dotted line in Fig. 1). Such a temperature dependence has been predicted by the 3D-SDW theory for systems at an AF QCP Hertz ; Millis ; Moriya .

The application of magnetic fields to CeNi2Ge2 is found to gradually reduce the low-$`T`$ specific heat coefficient (see Fig. 1) and induces LFL behavior. For $`B2`$ T a constant $`\gamma _0(B)=C(T,B)/T`$ value is observed at lowest temperatures. Upon warming, $`C(T,B)/T`$ passes through a broad maximum at a characteristic temperature $`T^{}(B)`$ that increases linearly with increasing magnetic field Gegenwart CNG ; Aoki ; Koerner . In the inset of Fig. 1, we analyze the magnetic field dependence of the low-$`T`$ specific heat coefficient $`\gamma _0(B)`$ that measures the heavy quasiparticle mass in the LFL regime. Upon reducing the magnetic field towards zero, we are able to follow the evolution of this true $`T=0`$ property upon tuning the system towards the QCP. We observe (see inset Fig. 1) that the field dependence of the specific heat coefficient $`\gamma _0(B)`$ (at $`T=0`$) is very similar to the temperature dependence of $`C/T`$ at $`B=0`$.

Next we turn to the low temperature resisitivity of CeNi2Ge2 . In Fig. 2a three polycrystalline samples with different residual resistivities are compared Gegenwart CNG . The ”standard-quality” sample with $`\rho _0=2.7\mu \mathrm{\Omega }`$cm follows $`\mathrm{\Delta }\rho (T)=\rho \rho _0=\beta T^{1.5}`$ for more than two decades in temperature between 20 mK and 2.5 K, in perfect agreement with the 3D-SDW prediction Hertz ; Millis ; Moriya . For high-quality samples with a ten times lower residual resisitivty, the exponent is slightly smaller, i.e. deviates even stronger from FL behavior. A similar observation has been made on slightly-off stoichiometric Ce1+xNi2+yGe2+z polycrystals where the Ni-Ge ratio was varied by a few at-%. Exponents of $`ϵ=1.5`$ for $`\rho _03\mu \mathrm{\Omega }`$cm and $`1.3ϵ<1.5`$ for 0.17 $`\mu \mathrm{\Omega }`$cm $`\rho _03\mu \mathrm{\Omega }`$cm have been found Gegenwart Physica . The onset of superconductivity at very low tempeatures (Fig. 2b) observed by several different investigations Gegenwart CNG ; Grosche ; Gegenwart Physica ; Koerner ; Braithwaite  has not been detected by any bulk probe so far.

As demonstrated for a high-quality sample in Fig. 3, the application of magnetic fields forces the low-temperature resistivity to turn into a $`\mathrm{\Delta }\rho =A(B)T^2`$ behavior consistent with the specific heat results discussed above. The cross-over temperature, detected for the three different polycrystals, increases proportionally to $`B^{0.65}`$ (Fig. 4a), whereas the coefficient $`A(B)`$ in the LFL regime diverges as $`A(B)B^{0.6}`$ (Fig. 4b). Again, this field dependence of the $`T=0`$ scattering cross-section for QP-QP scattering is very similar to the temperature dependence of the corresponding $`A(T)=\mathrm{\Delta }\rho (T)/T^2T^{0.5}`$ at $`B=0`$. Thus, in the approach to the QCP, both as a function of $`T`$ and $`B`$, the QP mass and QP-QP scattering cross section grow consistent with the predictions of the SDW scenario for 3D critical AF spinfluctuations.

The specific heat and electrical resistivity measurements indicate that the quantum-critical fluctuations, responsible for the NFL behavior in CeNi2Ge2, are extremely sensitive to the application of magnetic fields. Fields of the order of 1 T only, are sufficient to induced a LFL state with a strongly field-dependent QP mass and QP-QP scattering cross section. On the other hand, the magnetic fields necessary to suppress the Kondo effect by a polarization of the magnetic moments are much larger. We studied the high-field magnetization $`M(B)`$ of CeNi2Ge2 at the Dresden pulsed-field facility up to 50 T (Fig. 5). In order to minimize the effect of eddy currents, induced by the large $`dB/dt`$ rate in this high-purity material with very low residual resistivity, the sample has been powdered and the grains have been mixed with electrical insulating paraffin. A steady field of 5 T, applied above the paraffin’s melting temperature, has been used to orient the powder. After subsequent cooling, the susceptibility of the oriented powder agreed well with that for a single crystal measured along the easy direction along the $`c`$-axis. As shown in Fig. 5, $`M(B)`$ shows a step-like anomaly around 43 T, corresponding to a peak in the susceptibility that grows in magnitude and sharpens upon reducing the temperature. A similar behavior has also been reported by Fukuhara et al. Fukuhara  who measured an (unoriented) free powder sample. They proposed that the origin of this metamagnetic anomaly is related to that observed in CeRu2Si2 Haen  and other HF systems. Apparently, this metamagnetic transition is clearly separated from the quantum-critical behavior observed in very small magnetic fields $`B<1`$ T.

## III $`\text{YbRh}_2\text{Si}_2`$

We next consider the HF metal YbRh2Si2 for which pronounced NFL phenomena, i.e., a logarithmic divergence of $`C(T)/T`$ and a quasi-linear $`T`$-dependence of the electrical resistivity below 10 K, have been observed above a low-lying AF phase transition Trovarelli Letter . We use magnetic fields to tune the system through the QCP. This investigation was motivated by previous studies on field-induced NFL behavior in the doped AF systems CeCu6-xAgx Heuser98a ; Heuser98b  and YbCu5-xAlx Seuring . YbRh2Si2 is best suited to tune into quantum criticality by magnetic field, because (i) the ordering temperature $`T_N=70`$ mK is the lowest among all undoped HF systems at ambient pressure, (ii) already a very small critical magnetic field $`B_c=0.06`$ T is sufficient to push $`T_N`$ towards zero temperature, and (iii) the influence of disorder is minimized in very clean single crystals with a residual resistivity of about 1 $`\mu \mathrm{\Omega }`$cm.

We first discuss the magnetic properties of YbRh2Si2 (Fig. 6). At high temperatures ($`T200\text{K}`$) the magnetic susceptibility of YbRh2Si2 measured along both major crystallographic directions ($`a`$ and $`c`$) follows a Curie-Weiss law with effective magnetic moments very close to the value of free Yb3+ ($`\mu _{\text{eff}}=4.5\mu _\text{B}`$); but due to the strong magnetocrystalline anisotropy there is a marked difference in the respective extrapolated values for the Weiss temperatures $`\mathrm{\Theta }^a9\text{K}`$ and $`\mathrm{\Theta }^c180\text{K}`$ Trovarelli1 . At $`T=2\text{K}`$ the magnetic susceptibility measured along the basal plane is about 20 times larger compared to the value measured with applied field parallel to the $`c`$-axis. Fig. 6a is indicating that in YbRh2Si2 the Yb3+ moments form an ”easy-plane” square lattice with a strongly anisotropic response. This anisotropy is also reflected in the isothermal magnetization (Fig. 6c): Whereas for fields applied in the easy plane a strongly nonlinear response is observed with a polarized moment of 1 $`\mu _B`$/Yb at 10 T, the magnetization along the hard ($`c`$-axis) is much smaller and almost perfectly linear up to 56 T. No indication of any anomaly is observed for this field orientation Custers . On the other hand, a clear kink is visible in the easy-plane magnetization at 10 T. Corresponding anomalies have been observed in the isothermal magnetoresistance and magnetostriction as well unpublished . According to specific heat and electrical resistivity measurements as a function of temperature, the HF state is completely suppressed for $`B10`$ T and the system behaves as a metal with polarized local magnetic moments unpublished . In the following, we will focus on the quantum-critical behavior in small fields, well below 10 T.

Temperature-dependent measurements of the ac-susceptibility Trovarelli Letter ; Gegenwart YRS  and electrical resistivity Gegenwart YRS  have revealed that the AF order at $`T_N=70`$ mK (Fig. 6b) is completely suppressed by a critical magnetic field $`B_c`$ of 0.06 T, applied in the easy plane, only. The isothermal magnetization, measured at a temperature $`TT_N`$, shows a kink at $`B_c`$, indicative for a continuous second-order phase transition (Fig. 6d). The value of the ordered magnetic moment derived from the magnetic entropy Gegenwart YRS  and $`\mu `$SR experiments Ishida usr  is well below 0.1 $`\mu _B`$/Yb, in accordance with the small $`M(B_c)`$. Thus a large fraction of the local moments appears to remain fluctuating within the easy plane in the AF ordered state. Their continuous polarization gives rise to a strong curvature in $`M(B)`$ for $`B>B_c`$. The low value of $`B_c`$ indicates that the AF and field-aligned states are nearly degenerate in YbRh2Si2. Thus, low-lying ferromagnetic (q $`=0`$) fluctuations are expected as indeed found by recent 29Si-NMR experiments Ishida NMR .

In Fig. 7 we show the evolution of the low-temperature resistivity upon applying magnetic fields along and perpendicular to the easy magnetic plane. At small magnetic fields the Néel temperature, determined from the maximum value of $`d\rho /dT`$, shifts to lower temperatures and vanishes at $`B_c=0.06`$ T (in the easy magnetic plane) and 0.66 T ($`Bc`$). At $`B=B_c`$, the resistivity follows a linear $`T`$-dependence down to the lowest accessible temperature of about 20 mK. This observation provides striking evidence for field-induced NFL behavior at the critical magnetic fields applied along both crystallographic directions Gegenwart YRS . Fields larger than $`B_c`$ induce LFL behavior in the specific heat coefficient $`C(T)/T=\gamma _0(B)`$, electrical resistivity $`\mathrm{\Delta }\rho (T)=A(B)T^2`$ and magnetic susceptibility $`\chi (T)=\chi _0(B)`$ below a characteristic temperature $`T^{}(B)`$ that increases linearly with increasing $`B`$ (Ref. Gegenwart YRS ).

In the following, we analyze the evolution of these coefficients that characterize the field-induced LFL state upon reducing $`B`$ towards the QCP at $`B_c`$. As shown in Fig. 8, the coefficient $`A(B)`$ diverges as $`(BB_c)^1`$. Since $`A(B)`$ measures the QP-QP scattering cross section, this divergence indicates that the whole Fermi surface undergoes singular scattering at the field-induced QCP Gegenwart YRS . Furthermore, as shown in the inset of Fig. 8, we observe a constant ”Kadowaki-Woods” ratio $`K=A/\gamma _0^2`$ with a value typical for other HF systems KW  in $`B0.5`$ T. Since $`\chi _0`$ scales with $`A`$ and $`\gamma _0`$, also the QP mass and Pauli susceptibility diverge upon approaching the QCP. The very large Sommerfeld-Wilson ratio $`R=(\chi _0/\gamma _0)(\pi ^2k_B^2/\mu _0\mu _{eff}^2)`$ of about 14 $`(\mu _{eff}=1.4\mu _B)`$ indicates the strongly enhanced (q $`=0`$) susceptibility in the field-aligned state of YbRh2Si2.

## IV Conclusion

We have performed a comparative study of thermodynamic, magnetic and transport measurements on the two clean, stoichiometric and isostructural HF compounds CeNi2Ge2 and YbRh2Si2 which have been tuned into quantum criticality by magnetic fields. Whereas YbRh2Si2 shows very weak AF order at $`T_N=70`$ mK, suppressed by a tiny critical field $`B_c=0.06`$ T, CeNi2Ge2 is located almost directly at the QCP ($`B_c0`$). We have used magnetic fields $`b>0`$ ($`b=BB_c`$) to induce a FL state in both materials and followed the field dependence of the QP mass being proportional to $`\gamma _0(b)`$ and QP-QP scattering cross section (proportional to $`A(b)`$) upon approaching the QCP by reducing $`b0`$. For CeNi2Ge2 the evolution of both properties is in perfect agreement with the predictions of the 3D-SDW scenario, whereas for YbRh2Si2 the observed linear $`T`$-dependence of the electrical resistivity at $`b=0`$ as well as the $`1/b`$ divergence of $`A(b)`$ would fit the SDW prediction only, by assuming truly 2D critical spinfluctuations that render the entire Fermi surface ”hot”. The same scenario, however, predicts a logarithmic divergence of the QP mass at the QCP Kotliar , i.e. $`C(T)/T\mathrm{ln}(T)`$ at $`b=0`$ and $`\gamma _0(b)\mathrm{ln}(b)`$ in the field-induced FL state. Very recently, careful specific heat experiments on YbRh2(Si1-xGex)2 with a (nominal) Ge-concentration of $`x=0.05`$ (Ref. Trovarelli Physica ) for which $`T_N`$ and $`B_c`$ ($`c`$) are reduced by a slight volume expansion to 20 mK and 0.025 T, respectively, have been performed Custers Nature . At $`b=0`$ a power-law divergence $`C(T)/TT^{1/3}`$ is observed below about 0.5 K. Additionally, $`\gamma _0(b)`$ in the field-induced LFL state for $`b>0`$ diverges as $`b^{1/3}`$ for fields lower than 0.5 T. The stronger than logarithmic mass divergence is incompatible with the SDW scenario. The latter predicts the ”Kadowaki-Woods” ratio, $`K`$, to divergence as $`K_{SDW}(b\mathrm{ln}^2(b))^1`$. A truly constant $`K`$ as found above $`B=0.5`$ T in YbRh2Si2 (see inset Fig. 8) would indicate that the scattering amplitude remains local in the approach to the QCP. The weak field dependence $`Kb^{1/3}`$ observed for YbRh2(Si0.95Ge0.05)2 for $`b0`$ implies that the characteristic length scale of the scattering amplitude renormalizes more slowly than expected by the SDW scenario and thus favors a locally-critical scenario for YbRh2Si2 (Ref. Custers Nature ). According to the model by Si et al. (Ref. Si ), the local QCP occurs for systems with strongly anisotropic, i.e. quasi-2D critical spinfluctuations and results in the emergence of spatially local critical excitations that co-exist with the spatially extended critical spin fluctuations. Thus, to understand the different behaviors observed for CeNi2Ge2 and YbRh2Si2 one would assume the magnetic fluctuations in the relevant energy range to be 3D for the former and 2D for the latter system. For CeNi2Ge2 recent inelastic neutron scattering experiments determined the critical fluctuations and found them to behave 3D, with no indications for a 2D-like anisotropy Kadowaki . Similar experiments on YbRh2Si2 have not been performed yet.

To conclude, tuning the quantum critical points in CeNi2Ge2 and YbRh2Si2 by magnetic fields revealed that for the former system can be well described within the itinerant 3D-SDW model, whereas a locally critical model seems to be appropriate for the latter system. For both systems, the quantum-critical behavior observed at low magnetic fields is well separated from a high-field scale at which the HF state is suppressed by a full polarization of the magnetic moments.

We are grateful to C. Pépin, P. Coleman, and Q. Si for stimulating discussions. Work in Dresden has been supported in part by the Fonds der Chemischen Industrie.