format: text

text:
# Strange Attractors in Multipath propagation: Detection and characterisation.

## I Introduction

Multipath propagation of radio waves in indoor or outdoor environments shows a highly irregular behavior as a function of time . The characterization of radio channels in mobile or in building propagation is important for addressing issues of design, coding, modulation and equalization techniques tailored specifically to combat time and frequency dispersion effects.
Irregular behavior of the received signal has prompted researchers in the past to model the channel with stochastic processes. One of the earliest linear models in this vein is the Turin et al. model  in which the impulse response of the channel is written as a superposition of replicas of the transmitted signal delayed and having altered amplitudes and phases.
A number of models exist differing in the numbers of replicas of the signal or in the type of probability distributions from which the amplitudes and phases are drawn. Also, different stochastic processes are used in the generation of delay times of received replicas. The popular choice for the amplitude probability density functions (PDF) are Rayleigh or Rice PDFs depending on whether a weak or strong line of sight propagation exists; nevertheless other PDFs have been used such the lognormal, Nakagami-m or uniform. The phases ought to be drawn from PDFs compatible with the ones selected for the corresponding amplitudes; nevertheless the most popular choice found in the literature is the uniform $`[02\pi ]`$ distribution. Delay times are usually extracted from either stationary of stationary Poisson processes although in some cases the Weibull PDF is used.

## II A HYPOTHESIS OF DETERMINISTIC CHAOS

Although an assumption of stochastic behavior in mobile or indoor propagation is ubiquitous, in the present work we set out to explore the hypothesis that the indoor communication channel displays deterministically chaotic behavior. This question is important in many respects and the tools to answer it readily exist. These tools are based on the determination of the correlation dimension of the strange attractor associated with the multipath profile considered as a real valued time-series x(t). Using delay coordinates  one forms the m-dimensional delay vector $`𝐗(t)=[x(t),x(tT)\mathrm{}x(t(m1)T]`$ with delay $`T`$ and computes the correlation sum C(r) which is the ratio of the number of pairs of delay vectors (the distance between which is less than r) to the total number of pairs. From this, the correlation dimension $`\nu `$ is defined as the logarithmic slope of C(r) versus r for small r. For a true stochastic process, $`\nu `$ increases with m without showing any saturation. In contrast, for deterministic chaos, $`\nu `$ saturates at a value, the next integer greater than which, represent the minimum number of non-linear recursion or differential equations from which it originates. If the profile turns out to be deterministic, the modulation, coding and detection/demodulation techniques ought to be adapted accordingly in order to account for this fact; otherwise one has to rely upon techniques capable of handling stochastic signals. Let us illustrate this by an analysis of multipath measurements we have made in an indoor environment.

## III EXPERIMENT AND CORRELATION DIMENSION ANALYSIS

The propagation environment from which the data are collected are hallways in the Engineering Building at The University of Calgary . The transmitted power was 10 dBm fed into a half wave dipole antenna with a matching balun. The receiving antenna was a cross-polarized dipole array. The co-polarized antennas (CPA) and crosspolarized antennas (XPA) profiles referred to in Figure 1 are from a single measurement run and points from both profiles were obtained in coincident pairs. The terms CPA and XPA simply refer to the relative state of polarisation between the transmit and receive antennas. The receiving hardware was specially developed to measure diversity characteristics and gives an accurate reference between distance from transmitting to receiving antennas and received signal strength. The same measurement procedure was employed as for the arbitrarily polarized data set. We have estimated the correlation dimension for the sets of data: Arbitrarily Polarized Antennas (APA, 6000 data points), co-polarized (CPA, 3800 points) and Cross-polarized antennas (XPA, 3800 points) by the box-counting method of Grassberger and Procaccia . There are a number of limitations and potential pitfalls with correlation dimension estimation that have been discussed by various authors . In addition, the number of operations it takes to estimate C(r) is $`O(N^2)`$ where $`N`$ is the number of collected experimental points. Recently Theiler  devised a powerful box-assisted correlation sum algorithm based on a lexicographical ordering of the boxes covering the attractor, reducing the number of operations to $`O(Nlog(N))`$ and incorporating several test procedures aimed at avoiding the previous pitfalls. Before we used Theiler’s algorithm, we made some preliminary tests against well known cases. We generated uniform random numbers, Gaussian random numbers, and numbers z(n) according to the logistic one-dimensional map at the onset of Chaos: z(n+1) = a z(n)(1-z(n)) with a = 3.5699456 and in the fully developed Chaotic regime at a=4. In the first two cases we found $`\nu `$ approximately equal to m as expected in purely stochastic series whereas $`\nu `$ saturated respectively at 0.48 and 0.98 (we used 2000 points only) for the logistic map indicating the presence of a low dimensional attractor and deterministic Chaos (the exact correlation dimensions for the logistic map is 0.5 at a=3.5699456 and 1 for a=4.). We tested as well the Hénon two dimensional map, the Lorenz three dimensional system of non-linear differential equations, the Rössler three and four dimensional systems as well as an infinite dimensional system, the delay-differential Mackey-Glass equation whose attractor dimension is tunable with the delay time. The results we found for the various correlation dimensions agreed with all the results known in the literature to within a few percents. Then we went ahead and examined the $`\nu `$ vs. m curves for the three sets of experimental data along with a set of 6000 Rayleigh and band-limited Rayleigh distributed numbers which constitute prototypic received envelopes. Our results, in Figure 1 show that, while $`\nu m`$ for the pure Rayleigh case (with a slope equal to one), and $`\nu m`$ for the band limited Rayleigh case (with a slope smaller than one) the $`\nu `$ vs. m curves for the three examined experimental sets of data start linearly with m then show saturation indicating the presence of a low dimensional attractor (whose dimension is about 4 for CPA and XPA data whereas it is slightly above 4 for the APA situation). This finding is in line with the Ruelle criterion  that sets an upper bound on the possible correlation dimension one can get from any algorithm of the Grassberger-Procaccia type. This upper bound is set by the available number of data points $`N`$ in the time series. The dimension that can be detected should be much smaller than $`2log_{10}(N)`$. Since we used 3800 and 6000 points respectively, the upper bound for the detectable correlation dimension in our case is about 7.16 to 7.56. We respect this bound since our correlation dimensions are around 4. Nevertheless the presence of Chaos is going to be confirmed through another route, the spectrum of the Liapunov exponents that will be discusssed next.

## IV SPECTRUM OF THE LIAPUNOV EXPONENTS

The spectrum of Liapunov exponents is very important in the study of dynamical systems. If the largest exponent is positive, this is a very strong indication for the presence of Chaos in the time series originating from the dynamical system. The reciprocal of this exponent is the average prediction time of the series and the sum of all the positive exponents (if more than one is detected like in hyper-chaotic systems such as the Rössler four dimensional system of non-linear differential equations or the large delay Mackey-Glass equation) is the Kolmogorov entropy rate of the system. The latter gives a quantitative idea about the information processes going on in the dynamical system. In addition, with the Kaplan-Yorke conjecture, the full spectrum gives the Hausdorff dimension of the strange attractor governing the long time evolution of the dynamical system. We have calculated the Liapunov exponents of the data with four different methods. Firstly, we determined the largest exponent $`\lambda _{max}`$ from the exponential separation of initially close points on the attractor and averaging over several thousand iterations. Second, we determined the largest Liapunov exponent from the correlation sum with the help of the relation $`C(r)r^\nu exp(mT\lambda _{max})`$ valid for large values of the embedding dimension $`m`$ and small values of r. Finally, we determined all exponents with two different methods: the Eckmann et al. method  and the Brown et al.’s . Our results for the spectrum of exponents is shown in figures 2 and 3. We tried several embedding delay times T and several approximation degrees for the tangent mapping polynomial (as allowed in the Brown et al.  algorithm ) without observing major changes in the spectrum. Several time series (Logistic map, Hénon, Lorenz, Rössler and Mackey-Glass) were tested for the sake of comparison to results obtained with the experimental data. In addition, the Liapunov exponents saturate smoothly as they should for large embedding dimension. Then we applied the Kaplan-Yorke conjecture to get the dimension of the attractor: Using the following typical numbers we obtained for the exponents $`\lambda _{max}=\lambda _1=18.06,\lambda _2=1.88,\lambda _3=8.85,\lambda _4=24.94,\lambda _5=68.80`$ and using the formula:

$$D=j+\frac{_i\lambda _i}{|\lambda _{j+1}|}$$
(1)
where the summation is over i=1,2…j. The $`\lambda _i`$’s are ordered in a way such that they decrease as i increases. We determine j from the conditions $`_i\lambda _i>0\text{and:}\lambda _{j+1}+_i\lambda _i<0`$. We get j=3 and D=3.44 for the strange attractor dimension (called its Liapunov dimension). The total sum of the Liapunov exponents is negative ( equal to -82.65) as it should be for dissipative systems with a strange attractor. The value of the attractor dimension will be confirmed from the spectrum of singularities or the mutifractal spectrum in the next section.

## V MULTIFRACTAL SPECTRUM

The generalized dimension may be used to characterize non-uniform fractals for which there are different scaling exponents for different aspects of the fractal, so-called multifractals. For these, there are two scaling exponents, one generally called $`\tau `$, for the support of the fractal, and one called q, for the measure of bulk of the fractal. In general, $`\tau (q)=(q1)D_q`$, where $`D_q`$ is the generalized dimension. Multifractals have been employed to characterize multiplicative random processes, turbulence, electrical discharge, diffusion-limited aggregation, and viscous fingering . Multifractals have this in common: there is a non-uniform measure (growth rate, probability, mass) on a fractal support. Besides the exponents, $`\tau `$ and q, and the generalized dimension $`D_q`$, there is another method for characterizing multifractals. This depends upon the use of the mass exponent $`\alpha `$, and the multifractal spectrum, $`f(\alpha )`$ . A graph of the multifractal spectrum explicitly shows the fractal dimension, f, of the points in the fractal with the mass exponent (or scaling index), $`\alpha `$. We have estimated $`D_q`$ by use of the generalised moments of the correlation sum with a window chosen carefully enough to avoid temporal correlation effects. We have developed a program that computes the generalized correlation sum using a box-assisted method. Our program is based on one written by Theiler . Several modifications had to be made to the straightforward box-assisted correlation sum method. In addition, our program allows for logarithmic scaling with the distance parameter r. From a log-log graph of the generalized correlation sums, appropriate scaling regions can be identified, for each order, q. Least-squares fits to these scaling regions yields a sequence of generalized correlation dimensions, $`D_q`$, for values of q between $`\pm \mathrm{}`$. We have found computation of $`D_q`$ for integers in the interval \[-10,10\] and $`D_{\mathrm{}}`$ and $`D_+\mathrm{}`$ to be sufficient. From the $`D_q`$, we calculate the $`\tau (q)=(q1)D_q`$. We then perform a Legendre transform to obtain the $`f(\alpha )`$ curve. We do this by first fitting a smooth curve (a hyperbola was considered to be adequate) to the $`\tau (q)`$ curve. With an analytic expression for the $`\tau (q)`$ curve, we can compute the Legendre transform in closed form. The domain of $`f(\alpha )`$ may be found from $`D_{\mathrm{}}`$ and $`D_+\mathrm{}`$; we assume that $`\alpha `$ is confined to this region, and that $`f(\alpha )`$ is 0 at these points. The values of $`f(\alpha )`$ for $`D_{\mathrm{}}\alpha D_+\mathrm{}`$ are calculated as min\[$`q\alpha \tau (q)`$\], the minimum being taken over q. We found that this procedure, although complex, corrects for the known numerical sensitivities of the Legendre transform. We checked that our method for obtaining $`f(\alpha )`$ gave the same results as those found in the literature for the Logistic map and the strongly dissipative circle map . The $`f(\alpha )`$ curve for the multipath data is shown in Figure 4. It may be seen that the peak value of $`f(\alpha )`$, corresponding to the box-counting dimension at q=0, is about 3.7. This is consistent with our above findings from the correlation dimension and the full spectrum of Liapunov exponents. Our further research in this area concerns the prediction of the received signal intensity, from our above hypothesis of the presence of deterministic chaos.

## VI NON LINEAR PREDICTION

We applied the above findings to the non linear prediction of multipath profiles considering that each point on the envelope of the measured signal is a function of some number of past points in the series, deviating from the traditional wave superposition approach. More precisely, we write:

$$y(n+1)=F[y(n),y(n1),y(n2),y(n3)\mathrm{}y(nm+1)]$$
(2)
with F an $`m`$ dimensional map and y(n) the value of the signal x(t) sampled at timestep n. Expressing F as a sum of sigmoidal and linear functions  we determine the unknown weights through the Marquardt least squares minimisation method  in order to achieve the best fit to the data. Our results comparing the onestep prediction to the multipath data analysed above are displayed in figure 5. The goodness of fit between the measured and predicted envelopes for a map dimension m=5 is another indication of the soundness of the approach. This is confirmed in Figure 6 where we display the normalised prediction error versus the embedding dimension. A minimum is observed in the prediction error for an embedding equal to 5 or 6. In Figure 6, we started always from the same initial weights and let the system run through 150 iterations searching for the least squares minimum for 200 data points and later on making 300 one step ahead predictions. For embedding dimension larger or equal to 7 the minimisation procedure stopped because of the presence of zero-pivot in the least-squares matrices. The presence of a minimum in the prediction error around an embedding equal to 5 or 6 complies again with the value of embedding dimension used previously in the correlation dimension analysis, the Liapunov spectrum and the Kaplan-Yorke conjecture.

## VII DISCUSSION

We stress that although the profiles we examined were found to be chaotic in all three experimental configurations with confirmations from the Liapunov spectrum and non linear prediction studies, indicating that we would be able to describe our data with a set of at most 5 non-linear differential or algebraic equations, investigation in other propagation situations is needed. Nevertheless, in our investigations we observed a significant amount of consistency between the various methods of detecting Chaos and characterising it using embedding dimensions $`m`$ beyond the minimum $`m_{min}`$ required by the Takens theorem ($`m_{min}>2d+1`$, where d is the dimension of the strange attractor). Assuming the hypothesis of the presence of Chaos in a given multipath profile is firmly established, many avenues become possible. For instance, one might consider devising ways for controlling the signal propagation by altering slightly some accessible system parameter and improving the performance characteristics of the channel . Shaw  has introduced the concept that a deterministically chaotic system can generate entropy. The consequences of this observation is important for the design of communication equipment when the channel is a chaotic system. For one, it implies that information at the receiver about the state of the channel is lost at a mean rate given by the Kolmogorov entropy. For another it implies that a channel estimator should be adapted to the mathematical nature of the set of non linear equations describing the channel as shown in the previous paragraph. Our studies up to this date have shown that this approach is valid in an indoor situation but not in an outdoor one. This might be due to the confined geometry one encounters inside buildings and the boundary conditions for the electromagnetic fields leading to a low dimensional system of non linear equations giving birth to the observed chaotic behaviour. Our studies in this direction are in progress.

Acknowledgements

We thank James Theiler, Jean-Pierre Eckmann and Reggie Brown for sending us their computer programs and correspondance, as well as Halbert White for some unpublished material. C.T. thanks Sunit Lohtia and Bin Tan for their friendly help with the manuscript.

Figure Captions

* Correlation dimensions vs embedding dimension: Full squares are for Rayleigh distributed points; full triangles are for handlimited Rayleigh distributed points. Full diamonds are for experimental results in the XPA case whereas open diamonds correspond to the CPA case and open squares to the APA case.
* Liapunov exponent spectrum from Eckmann et al.  method versus embedding dimension for the APA data (since APA data consist of the largest number of points, 6000).
* Liapunov exponent spectrum from Brown et al. method  for the same data as those of Fig.2. A linear tangent mapping is used to fit the dynamics. A very similar spectrum is obtained for a second order polynomial.
* Spectrum of generalised dimensions $`f(\alpha )`$ versus a for the APA data used in fig.2. The value at the maximum of $`f(\alpha )`$ corresponding to the Hausdorff dimension of the strange attractor agrees with the minimum bound obtained from fig.1 and with the Kaplan-Yorke conjecture (see text). The spectrum is obtained through embedding in 10 dimensions. This happens to be enough, given the values obtained for the various generalised dimensions.
* Measured envelope (APA data used in Fig.2 continuous curve) and its one step prediction (dashed curve) from the Neural Network fit to the five dimensional map F (eq.2). The training is over the first 200 first points.
* Normalised prediction error versus embedding. Starting from the same initial weights, we trained the neural network, for a given embedding, over the first 200 points with a Marquardt minimisation standard deviation parameter equal to 0.01 and total number of 150 iterations. Once, the parameters at the minimum error are found we made a one-step ahead prediction over the next 300 points and calculated the resulting squared error divided by the total of points. One sees a minimum for an embedding dimension around 5 or 6. For an embedding equal to 7 or larger, a large error or no convergence (null pivot encountered in the least square error matrices) were observed.


text:
# Extreme-Energy Cosmic Rays: Puzzles, Models, and Maybe Neutrinos

## I INTRODUCTION

An unsolved astrophysical mystery, now forty years old, is the origin and nature of the extreme energy cosmic ray primaries (EECRs) responsible for the observed events at highest energies, $`10^{20}`$ eV . About twenty events at $`10^{20}`$ eV have been observed by five different experiments . The origin of these events is a mystery, for there are no visible source candiates within 50 Mpc except possibly M87, a radio-loud AGN at $`20`$ Mpc, and Cen-A (NGC5128), a radio galaxy at 3.4 Mpc, and neither of these is in the direction of any of the observed events . Since the observed events display a large-scale isotropy, many sources rather than one source seem to be required. The nature of the primary particle is also mysterious, because interactions with the 2.73K cosmic microwave background (CMB) renders the Universe opaque to nucleons at $`E_{\mathrm{GZK}}5\times 10^{19}`$ eV, and double pair production on the cosmic radio background (CRB) renders the Universe opaque to photons at even lower energies. The theoretical prediction of the end of transparency for nucleons at $`E_{\mathrm{GZK}}0.5\times 10^{20}`$ eV is the famous “GZK cutoff” . Figure 1 shows a recent compilation of the AGASA data set, clearly extending beyond $`E_{\mathrm{GZK}}`$.

The main theory challenges in attempting to understand the super-GZK data are: (i) What cosmic source could have accelerated the primary particles to such extremely high energies? and (ii) If the sources are distant ($`\stackrel{>}{}100`$ Mpc), then how could their primaries have propagated through the cosmic background radiation without substantial energy loss? The acceleration mechanism either requires a Zevatron accelerator ($`1ZeV10^{21}`$ eV) , distant because such a source could not be missed if it were nearby; or speculative decaying super-massive particles (SMPs) or topological defects (TDs) with mass-scale $`\stackrel{>}{}10^{22}`$ eV, clustered nearby; or possibly magnetic monopoles accelerated by the cosmic magnetic fields. From distant Zevatrons, only neutrinos among the known particles can propagate unimpeded to earth. Exploiting this fact are the “Z-burst model,” and the strongly-interacting neutrino model. Adding to the drama and mystery at present is the observed large-scale isotropy and small-scale anisotropy. Surely these characteristics hint at a solution to the mystery of origin; protons, nuclei, and magnetic charges bend in cosmic magnetic fields, whereas photons and neutrinos do not. With more data above $`E_{\mathrm{GZK}}`$, several distinct telltale signatures including the isotropy/anisotropy will allow one to discriminate among the many models proposed for the origin and nature of the EECRs.

## II EVENT CLUSTERING

One revealing signature is already evident in the existing data sample. This is the pairing of events on the celestial sky. The AGASA experiment has already presented data strongly suggesting that directional pairing is occuring at higher than chance coincidence . Of the 47 published AGASA events above $`4\times 10^{19}`$ eV, 9 are contained in three doublets and one triplet with separation angle less than the angular resolution of $`2.5^{}`$. The chance probability of this clustering occuring in an isotropic distribution is less than 1%. The chance probability for the triplet alone is only 5%. Of the seven events above $`10^{20}`$ eV, three are counted among the doublet events. Most recently, AGASA has reported two more events above $`E_{\mathrm{GZK}}`$ . Each aligns in direction with a previous event, reducing further the probability for random clustering to 0.07% at 3$`\sigma `$. Comparisons of event directions in a combined data sample of four experiments further supports non-chance coincidences, especially in the direction of the SuperGlactic Plane .

Whether the pairing is random or dynamical, we shall also know in the near future. If the pairing turns out to be dynamical, I would argue that neutrinos are a favorite candidate for the primary particles. This is because photons have such a short ($`10`$ Mpc) abosrption length, and protons are bent by cosmic magnetic fields during their extragalactic journey.

If neutrinos are the primaries, they should point back to their sources, thereby enabling point-source astronomy for the most energetic sources of flux at and above $`10^{20}`$ eV. It was reported  that the first five events at $`10^{20}`$ eV did in fact point toward extragalactic compact radio-loud quasars, just the class of objects which could accelerate EECRs to ZeV energies via shock mechanisms. With the inclusion of subsequent data, this association is controversial . The jury awaits further evidence.

### A Coincidence or physics?

A random distribution is obtained by tossing $`n`$ events randomly into $`N(\mathrm{\Omega }/\pi \theta )^2=1045(\mathrm{\Omega }/1.0\mathrm{sr})(\theta /1.0^{})^2`$ angular bins, where $`\mathrm{\Omega }`$ is the solid angle on the celestial sphere covered by the experiment and $`\theta `$ is the bin half-angle. Each resulting event distribution is specified by the partition of the $`n`$ total events into a number $`m_0`$ of empty bins, a number $`m_1`$ of single hits, etc., among the $`N`$ angular bins. The probability to obtain a given event topology is

$$P(\{m_i\},n,N)=\frac{1}{N^n}\frac{N!n!}{_{j=0}m_j!(j!)^{m_j}}.$$
(1)
The variables in the probability are not all independent, as $`_{j=1}j\times m_j=n`$ and $`_{j=0}m_j=N`$. It is useful to use these constraints to rewrite this probability as

$$P(\{m_i\},n,N)=\frac{N!}{N^N}\frac{n!}{n^n}\underset{j=0}{}\frac{(\overline{m_j})^{m_j}}{m_j!},$$
(2)
where we have defined

$$\overline{m_j}N\left(\frac{n}{N}\right)^j\frac{1}{j!}.$$
(3)
When $`Nn1`$, a limit valid for the AGASA, HiRes, Auger and Telescope Array experiments, one finds

$$P(\{m_i\},n,N)\underset{j=2}{}\frac{(\overline{m_j})^{m_j}}{m_j!}e^{\overline{m_j}r^j(j2)!},$$
(4)
where $`r(Nm_0)/n1`$. The non-Poisson nature of eq. (4) is reflected in the factorials and powers of $`r`$ in the exponents.

Typical values of effective area $`A`$ (km2 sr), celestial solid angle $`\mathrm{\Omega }`$, and angular resolution $`\theta _{\mathrm{min}}`$ for the existing and proposed EECR experiments  are shown in the following Table, where the incident flux $`F(E_{\mathrm{GZK}})=10^{19}\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1`$ has been used to estimate the number of events ($`n`$/yr) above $`E_{\mathrm{GZK}}=5\times 10^{19}`$ eV.

|  | AGASA | HiRes | Auger/TA | EUSO/OWL |
| --- | --- | --- | --- | --- |
| $`\mathrm{km}^2\mathrm{sr}`$ | 150 | 800 | $`6\times 10^3`$ | $`3\times 10^5`$ |
| $`n/\mathrm{yr}`$ | 5 | 30 | 200 | $`10^4`$ |
| $`\mathrm{\Omega }(\mathrm{sr})`$ | 4.8 | 7.3 | 4.8 | 4$`\pi `$ |
| $`\theta _{\mathrm{min}}`$ | $`3.0^{}`$ | $`0.5^{}`$ | $`1.0^{}`$ | $`1.0^{}`$ |

In Figure 2 are shown inclusive probabilities for the Auger experiment with 100 events, as determined by formulae eq. (4). Note that “inclusive probability” means the stated number of $`j`$-plets plus any other clusters. The 8-doublet probability is extremely sensitive to the angular binning; observation of a flatter dependence on angular bin-size could signal a non-chance origin for the clustering. The observation of two triplets with angular binning of less than $`2^{}`$ would constitute 3-sigma evidence for cosmic dynamics. The same conclusion would hold if a quadruplet within $`2.5^{}`$ is observed among the first 100 Auger events.

Figure 2: Inclusive probabilities for various clusters in a 100 event sample at Auger. The solid line is the exact result, the dashed line is the Poisson approximation

Without some modification, our analytic formula may not be directly relevant to HiRes data. Clear, moonless nights are required for detection of atmospheric fluorescence, and summer nights ($`18`$ hrs right ascension) are effectively 40% shorter than winter nights ($`6`$ hrs right ascension) for monocular HiRes . Accordingly, the HiRes efficiency versus Galactic longitude varies significantly, roughly as $`N(RA)=\overline{N}(1+ϵ\mathrm{sin}(2\pi RA/24\mathrm{h}\mathrm{r}\mathrm{s}))`$, with $`ϵ0.25`$. One must ask whether this sinusoidally-varying efficiency invalidates the analytic approach with its assumption of a constant efficiency. In the Monte Carlo approach there is an easy method to generate “background” data sets with the experimental efficiency properly incuded – one randomly permutes the RA coordinates of the real data to remove any dynamical correlation among the events. Fortunately, there is also a simple method to estimate the efficiency correction to an analytic Poisson distribution, if one assumes that although $`N`$ and $`n`$ vary with right ascension, $`\xi n/N`$ is constant. Then it is easy to show that the Poisson distribution

$$P(j,m)=\left(\frac{N\xi ^j}{j!}\right)^m\frac{1}{m!}e^{N\xi ^j/j!}$$
(5)
is corrected upwards by the factor

$$=\frac{1}{2\pi }_0^{2\pi }𝑑\varphi (1+ϵ\mathrm{sin}(\varphi ))^me^{\overline{N}ϵ\mathrm{sin}(\varphi )\xi ^j/j!}.$$
(6)
For $`\frac{ϵ\xi ^{j2}}{j!}\frac{n^2}{\overline{N}}1`$, the exponential is near unity and the correction has a closed form, $`_m(1ϵ^2)^{m/2}P_m(\frac{1}{\sqrt{1ϵ^2}})`$, with $`P_m`$ a Legendre polynomial. As a series, the correction is

$$_m=2^m\underset{k=0}{\overset{[m/2]}{}}\frac{()^k(2m2k)!}{k!(mk)!(m2k)!}(1ϵ^2)^k.$$
(7)
In particular, for $`n^2\overline{N}`$, we have $`_1=1`$ independent of $`ϵ`$, $`_2=1+ϵ^2/2,=1.03`$ for HiRes, $`_3=1+3ϵ^2/2,=1.09`$ for HiRes, and $`_4=1+3ϵ^2+3ϵ^4/8,=1.189`$ for HiRes. The lesson is that HiRes efficiency corrections appear negligible for small numbers of clusters, and we may proceed with our analytic analysis for the HiRes experiment, for which about 20 events at $`10^{20}`$ eV are expected when the first full year’s data is analyzed. We display in Figure (3) the inclusive probabilities for one or more, two or more, and three or more doublets; and one or more triplets, over a range of angular binning.

Figure 3: Inclusive probabilities for various clusters, given 20 events at HiRes. The solid line is the exact result, the dashed line is the Poisson approximation.

Note that for all except the 3 doublet configuration, the Poisson approximation using the mean values in Eq.(3) provides an estimate good to within 50% of the non-approximate form; for the (much suppressed) 3 doublet configuration, it overestimates the probability by about a factor of 3 in much of the angular region. For angular binning tighter than $`2^{}`$, an observation of two doublets among the first 20 events has a chance probability of less than 0.5%. Thus the observation of this topology could be construed as evidence (at the 3$`\sigma `$ level) for clustering beyond statistical. The observation of a triplet within $`3^{}`$ has a random probability of less than $`10^3`$, and hence observation of such a triplet would most likely signify clustered or repeating sources, or magnetic focusing effects. With the accumulation of 40 events (not shown in the figure), the appearance of two doublets has a probability of less than 0.5% for a correlation angle of $`1^{}`$ or less. This illustrates how the good angular resolution of HiRes may be used to detect non-statistical clustering with only a few observed clusters.

Projected event rates for the EUSO/OWL/AW experiments present a pleasant problem for our analytic formula. In the case where $`n>N1`$, relevant for the EUSO/OWL/AW experiments after a year or more of running, higher $`j`$-plets are common and the distribution of clusters can be rather broad in $`j`$. From eq. (3) we have $`(\overline{m_j}/\overline{m_{j1}})=(n/jN)(\pi n\theta ^2/j\mathrm{\Omega })`$. Already at $`j=1(2)`$, Stirling’s approximation to $`j!`$ is good to $`8\%(4\%)`$, and so for $`j1`$ we may approximate $`\overline{m_j}(N^3/2\pi en)^{\frac{1}{2}}(en/jN)^{j+\frac{1}{2}}`$. Extremizing this expression with respect to $`j`$, one learns that The most populated $`j`$-plet occurs near $`jn/N`$. Combining this result with the broad distribution expected for large $`n/N`$, one expects clusters with $`j`$ up to $`\mathrm{several}\times \frac{n}{N}`$ to be common in the EUSO/OWL/AW experiments. Probabilities and meaning for frequent, large clusters are somewhat difficult to assess numerically.

The random distributions displayed in figs. (2) and (3) are expected from some models, such as randomly-situated decaying super-massive particles (SMPs), or charged-particle or magnetic monopole primaries with directions randomized by incoherent cosmic magnetic fields. A complementary approach to the random probabilities shown here is to consider specific source models generating non-random angular distributions. Steps along this line of inquiry have recently been taken . Future progress in the field will involve comparisons of the random and non-random model predictions with the data.

## III MODELS FOR SUPER-GZK EVENTS

The conjectured origins of the super-GZK particles fall into four basic categories. These are (i) nearby accelerators, (ii) exotic primaries, (iii) exotic physical law, and (iv) neutrino primaries. Other reviews  have emphasized the first three categories (especially (i) and (ii)), so here we will be brief with those, and put more emphasis on the neutrino-primary option.

### A Nearby accelerators

Several types of sources have been proposed to exist within our Galactic halo. These include highly-ionized relativistic dust grains, Galactic supershocks, young neutron stars, magnetars (highly magnetized pulsars), decaying SMPs with GUT masses of order $`10^{15}`$ GeV or with inflation-motivated masses of order $`10^{13\pm 1}`$ GeV, topological defects  such as strings, Q-balls and vortons, and annihilating monopole-antimonopole bound states (monopolonium) . For rare sources emitting charged particles, such as magnetars, it is necessary to postulate that the primaries are iron nuclei to ease the acceleration requirement and to isotropize the flux in our galactic magnetic field. It is intriguing that the observed magnitude of CMB fluctuations fixes the reheat temperature following inflation to $`10^{13\pm 1}`$ GeV, which allows gravitational and thermal production of TDS and/or SMPs of just this mass (if they exist) ; this mass scale is just right for producing $`10^{20}`$ secondary particles via decay.

For some TD models, dimensional arguments and scaling laws seem to give an emission rate to short-lived SMPs consistent with the observed EECR rate without tuning exotic parameters. However, for most decaying sources such as SMPs, TDs, and monopolonium, it is necessary to tune the lifetime to be longer but not too much longer than the age of the Universe in order to maintain an appreciable secondary particle emission rate today. Discrete gauged symmetries  or hidden sectors  are introduced to stabilize the heavy particle. Then rather esoteric physics is needed to break the new symmetry super-softly to maintain the long lifetime. High-dimension operators, wormholes , and instantons  bibiteminstanton have been invoked for this purpose.

Models with sources mainly in the Galactic plane are disfavored by the lack of any observed planar asisotropy in the data. Models with sources clustered in the Galactic halo predict a dipole enhancement in the direction of the Galactic center , which will be tested by the Auger Observatory in a few years.

Possible sources outside our halo, but still relatively nearby, include the radio-loud quasar M87 at 18 Mpc , the similar Cen-A at 3.4 Mpc , rare nearby GRBs, and now-dormant rare AGNs (also called massive dark objects – MDOs). These sources are few at best, and strong magnetic fields must be postulated to isotropize and/or confine their emissions. In the case of GRBs, the identification of the red-shift of their host galaxies at typically $`\stackrel{>}{}1`$ renders a local occurence of a GRB highly improbable .

Finally, primordial black holes (PBHs) have also been suggested, but any sensible initial mass spectrum is unable to provide a sufficient number of PBHs in the final stage of decay today.

### B Exotic primaries

The GZK cutoff can be raised by simply postulating a primary hadron slightly heavier than the proton. The reason is kinematical – the cutoff energy varies as the square of the mass of the first excited resonant state. For this reason, Farrar has proposed that light supersymmetric baryons, made from a light gluino plus the usual quarks and gluons, may be the primary EECRs . Such a scenario renders the source-energetics issue even more challenging. In any event, terrestrial experiments seem to have recently closed the window on a possible light gluino.

Another interesting possibility for the primary EECR is the magnetic monopole . Any breaking of a semisimple gauge symmetry which occurs after inflation and which leaves unbroken a $`U(1)`$ symmetry group may produce an observable abundance of magnetic monopoles. The monopole mass is expected to be $`\alpha ^1`$ times the temperature $`T_c`$ of the symmetry breaking. At the time of the phase transition, roughly one monopole or antimonopole is produced per correlated volume. The resulting monopole number density today is

$$n_M10^{19}(T_c/10^{11}\mathrm{GeV})^3(l_H/\xi _c)^3\mathrm{cm}^3,$$
(8)
where $`\xi _c`$ is the phase transition correlation length, bounded from above by the horizon size $`l_H`$ at the time of the transition. In a second order or weakly first order phase transition, the correlation length is comparable to the horizon size. In a strongly first order transition, the correlation length is considerably smaller than the horizon size.

The kinetic energy imparted to a magnetic monopole on traversing a magnetic field $`B`$ is $`E_KgB\xi `$, where $`g=e/2\alpha `$ is the magnetic charge according to the Dirac quantization condition, and $`\xi `$ specifies the field’s coherence length. Given the magnitude and coherence length data for the cosmic magnetic fields, monopole kinetic energies in the range $`10^{20}`$ to $`10^{23}`$ eV are expected; the acceleration problem is naturally solved. Monopoles with $`M\stackrel{<}{}10^{14}`$ GeV should be relativistic, and carry the appropriate energy to qualify as candidates for the EECR primaries . Within field theory there exist many possibilities for an intermediate unification scale and intermediate-mass monopoles.

The propagation problem is also naturally solved. The scattering cross-section for the monopole on the 3K and diffuse photon backgrounds is just classical Thomson, valid even for strong coupling: $`\sigma _\mathrm{T}=8\pi \alpha _M^2/3M^22\times 10^{43}(M/10^{10}\mathrm{GeV})^2\mathrm{cm}^2`$. The resulting mean free path for inverse Compton scattering is many orders of magnitude larger than the Hubble size of the Universe.

The relativistic monopole flux is simply

$$F_M=cn_M/4\pi 2\times 10^{19}(M/10^{10}\mathrm{GeV})^3(l_H/\xi _c)^3$$
(9)
per cm2-s-sr, which compares favorably with the integrated flux above $`10^{20}`$ eV, $`F_{\mathrm{data}}(>10^{20}\mathrm{eV})2\times 10^{20}\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1`$, and is comfortably below Parker’s upper bound $`F_{\mathrm{Parker}}=10^{15}\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1`$ for a cosmic monopole flux.

Signatures for EECR monopoles in our atmosphere and in ice are discussed at length in . Of particular interest as a model for the super-GZK primaries is the “baryonic monopole” . It is a bound state of three colored monopoles, confined by chromomagntic strings. Chromomagnetic strings can stretch when excited, but cannot break into quark-antiquark pairs. On first interaction in the atmosphere, the baryonic monopole of mass $`M`$ and energy $`E=\gamma M`$ stretches to create a huge geometrical cross-section of order $`\gamma /\mathrm{\Lambda }_{\mathrm{QCD}}^2(\gamma /10^6)\times 10^7`$ mb. Consequently, nearly all of the initial monopole energy is transferred to the atmospheric shower in a very short distance. A recent simulation of the baryonic-monopole showed good agreement with the lateral muon and hadron content of the highest-energy Yakutsk event, but less than good agreement with the longitudinal profile of the Fly’s Eye event .

Any confirmed directional pairing of events would appear difficult to achieve with the monopole model. Also, in the context of a model of the Galactic magnetic field, it has been shown that some memory of the local spiral arm direction, and an energy spectrum flatter than the observed one, are expected in the data if monopoles are the primaries . Although the directional criterion may not survive inclusion of extragalactic fields, the flatness of the spectrum probably does.

### C Exotic physical laws

The most remarkable proposals posit a breakdown of Lorentz invariance  or a breakdown of general relativity above some high scale. In the case of broken GR, it could be spacetime fluctuations (expected in a theory of quantum gravity) wiggling the on-shell dispersion relation  or appearing as $`1/M_P`$ operators which alter the physics. String theory provides motivation for GR-breaking at a possibly lower scale, the string scale $`M_S`$. Certainly the energy window of EECRs ($`10^{20}\mathrm{eV}/M_P10^8`$) is beyond that of terrestrial accelerators, and so ripe for speculation on new high-energy physics.

Signatures of models with no photo-pion production above $`E_{\mathrm{GZK}}`$ include the absence of a proton pile-up below $`E_{\mathrm{GZK}}`$, the absence of a cosmogenic neutrino flux, and possibly undeflected pointing of the primary back to its source.

### D Neutrino primaries

Turning to the possibility that the primaries may be neutrinos, one encounters an immediate obstacle: the SM neutrino cross-section is down from that of an electromagnetic or hadronic interaction by six orders of magnitude. This implies a low air-shower rate, and an accumulation of events at low altitudes (“penetrating” events) where the target density is highest. On the other hand, the neutrino-primary hypothesis is supported by the observed clustering discussed earlier. Two solutions to the small cross-section problem for primary neutrinos have been proposed.

#### 1 Neutrino annihilation to $`Z`$–bursts

Here it is proposed that the primary particles which propagate across cosmic distances above the GZK cutoff energy are neutrinos, which then annihilate with the cosmic neutrino background (CNB) within the GZK zone ($`D<D_{\mathrm{GZK}}50`$ Mpc) to create a “local” flux of nucleons and photons above $`E_{\mathrm{GZK}}`$, as shown in Fig. (4). It was noted many years ago that a cosmic ray neutrino arriving at earth from a cosmically distant source has an annihilation probability on the relic–neutrino background of roughly $`3.0h_{65}^1\%`$ (neglecting cosmic expansion) . The probability for a neutrino with resonant energy to annihilate to a $`Z`$-burst within distance $`D_{\mathrm{GZK}}`$ is then $`2.5\times 10^4`$ for a homogeneous CNB . The annihilation rate depends upon the CNB density, reliably predicted in the mean by Big Bang cosmology, and on the Standard Model (SM) of particle physics. The local annihilation rate is larger if our matter-rich portion of the Universe clusters neutrinos , or if there is an intrinsic CP-violating $`\nu \overline{\nu }`$ asymmetry .

Figure 4: Schematic diagram showing the production of a $`Z`$–burst resulting from the resonant annihilation of a cosmic–ray neutrino on a relic (anti)neutrino. If the $`Z`$–burst occurs within the GZK zone ($`50`$ to 100 Mpc) and is directed towards the earth, then photons and nucleons with energy above the GZK cutoff may arrive at earth and initiate super–GZK air–showers.

Each resonant neutrino annihilation produces a $`Z`$ boson with a 70% branching ratio into hadrons known to include on average about one baryon–antibaryon pair, seventeen charged pions, and ten neutral pions . The ten $`\pi ^0`$’s decay to produce twenty high–energy photons. For $`m_\nu `$ in the range $``$ 0.1 to 2 eV, the energy in this resonant “$`Z`$-burst” is fortuitously situated sufficiently above $`E_{\mathrm{GZK}}`$ at

$$E_\nu ^R=M_Z^2/2m_\nu =4(\mathrm{eV}/m_\nu )\times 10^{21}\mathrm{eV}$$
(10)
so as to produce photons and nucleons with energies exceeding $`E_{\mathrm{GZK}}`$.*** The resonant-energy width is narrow, reflecting the narrow width of the Z-boson: at FWHM $`\mathrm{\Delta }E_R/E_R\mathrm{\Gamma }_Z/M_Z=3\%`$.  The mean energies of the $`2`$ baryons and $`20`$ photons produced in the Z decay are easily estimated. Distributing the Z-burst energy among the mean multiplicity of 30 secondaries in Z-decay , one has

$$E_p\frac{E_R}{30}1.3\left(\frac{\mathrm{eV}}{m_j}\right)\times 10^{20}\mathrm{eV}.$$
(11)
The photon energy is further reduced by an additional factor of 2 to account for their origin in two-body $`\pi ^0`$ decay:

$$E_\gamma \frac{E_R}{60}0.7\left(\frac{\mathrm{eV}}{m_j}\right)\times 10^{20}\mathrm{eV}.$$
(12)
Even allowing for energy fluctuations about mean values, it is clear that in the Z-burst model the relevant neutrino mass cannot exceed $`2`$ eV. On the other hand, the neutrino mass cannot be too light or the predicted primary energies will exceed the observed event energies. In this way, one obtains the approximate 0.1 eV lower limit on the neutrino mass, when allowance is made for an order of magnitude energy-loss for those secondaries traversing 50 to 100 Mpc.

The challenging issue of how experiments might actually determine the absolute neutrino mass is discussed in .

If the $`Z`$–burst points in the direction of earth and occurs within the GZK distance, then one or more of the photons and nucleons in the burst may initiate a super–GZK air–shower at earth . For a sufficient cosmic neutrino flux, the hypothesis successfully explains the observed air–showers above $`E_{\mathrm{GZK}}`$. Comparisons of the model predictions to super-GZK data are available in .

The existence of neutrino mass in the desired range seems nearly guaranteed from the tritium decay upper bound  and the lower bounds inferred from the terrestrial neutrino oscillation experiments. The simplest explanation for the atmospheric neutrino results is neutrino oscillations driven by a mass–squared difference of $`\delta m_{\mathrm{atm}}^23\times 10^3\mathrm{eV}^2`$ , which implies a neutrino mass of at least 0.05 eV. Also, the recent LSND measurement appears to indicate a mass–squared difference $`\delta m_{\mathrm{LSND}}^2\stackrel{>}{}0.2\mathrm{eV}^2`$ , from which one deduces a neutrino mass of at least 0.5 eV. From these lower bounds on neutrino mass, one gets upper bounds on the $`Z`$–burst energy of $`10^{23}`$ and $`10^{22}`$ eV, respectively, just right for extending the air-shower spectrum an order of magnitude or two beyond the GZK cutoff!

A considerable cosmic neutrino flux above $`E_{\mathrm{GZK}}`$ is required for the $`Z`$–burst hypothesis to successfully explain the super–GZK events. The requirement is that the product of the resonant energy times the neutrino flux at the resonant energy per flavor, times the annihilation probability within the GZK zone (which may be as large as 0.025% to 1% due to neutrino clustering), times the photon and nucleon multiplicity per burst ($``$ 20), is comparable to the observed flux at $`10^{20}`$ eV. The resulting requirement on the neutrino flux is roughly $`E_RF_{\nu _j}(E_R)10^{18.5\pm 1}`$/cm2/s/sr. Such a neutrino flux at $`E10^{22}`$ eV is directly measureable in a teraton ($`10^{12}`$ ton) detector like EUSO/OWL/AirWatch, and possibly in a search for radio pulses produced by high energy neutrinos penetrating a small column–density of matter in the limb of the moon .

While certainly large, this required neutrino flux violates no existing limits. It has been pointed out  that this flux cannot extrapolate as $`E^2`$ to $`10^{17}`$ eV, for then it would violate the Fly’s Eye bound arising from nonobservation of penetrating horizonal (i.e., neutrino induced) air showers at that energy. It has also been pointed out  that local neutrino clustering is required to avoid generating a density of 30 MeV to 100 GeV photons in excess of the EGRET experimental bound , from distant Z-bursts undergoing electromagnetic cascading. The extreme-energy neutrino flux implied by the $`Z`$-burst model probably requires unusual source dynamics . Among the reasons to hope that Nature obliges is that resonant neutrino annihilation provides the best hope at present to actually measure the relic neutrino density.

#### 2 Strong $`\nu `$ Cross-section at $`\stackrel{>}{}E_{\mathrm{GZK}}`$

It is interesting and suggestive that the observed EECR flux beyond $`E_{\mathrm{GZK}}`$ is well matched by the flux predicted for cosmogenic neutrinos. This is not a complete coincidence. With the GZK cutoff, any continued nucleon flux beyond $`E_{\mathrm{GZK}}`$ is degraded in energy, photo-producing pions which in turn decay to produce cosmogenic neutrinos. The number of produced neutrinos compensates for their lesser energy, with the result that the neutrino flux matches well to the observed super-GZK flux. One may entertain the notion that the cosmogenic neutrinos are the super-GZK primaries, and that these neutrinos acquire a strong cross-section at $`10^{20}`$ eV.

Limits on the strength of the neutrino cross-section at $`10^{20}`$ eV can be inferred from existing data. Heuristically, one argues as follows. The GZK process ensures that there is a flux of cosmogenic neutrinos at $`10^{20}`$ eV, with an easily calculated flux. If the neutrino cross-section were weak, an experiment looking for penetrating air-showers initiated by the cosmogenic neutrinos would see nothing. If the cross-section were strong enough, the neutrino could not penetrate the atmosphere at all. So the fact that the Fly’s Eye experiment saw no penetrating showers tells us that the neutrino cross-section is either strong or weak; the mid-range is excluded. The vertical depth of our atmosphere is $`x_\mathrm{v}=1033\mathrm{g}/\mathrm{cm}^2`$, and the horizontal depth $`x_\mathrm{h}`$ is about 36 times greater. In terms of the mean free path $`\lambda `$ of a particle with cross-section $`\sigma `$, one has $`x_\mathrm{v}/\lambda =\sigma /1.6\mathrm{mb}`$, and $`x_\mathrm{h}/\lambda =\sigma /44\mu \mathrm{b}`$. Thus, an estimate of the excluded cross-section is $`40\mu \mathrm{b}`$ to $`1\mathrm{m}\mathrm{b}`$. A more careful calculation has been performed, with the result that $`20\mu \mathrm{b}`$ to $`1\mathrm{m}\mathrm{b}`$ is excluded . Hypothetical high-energy neutrino cross-sections in excess of a mb remain viable.

The idea that neutrinos, indeed, all particles, may have a strong interaction at a high but observable energy scale is not new . However, some recent ideas concerning new interactions relate well to the $`10^{20}`$ eV scale. One idea is that leptons are bound states of dual QCD gluons, which reveal themselves just above the electroweak (EW) scale at parton-parton $`\sqrt{s}`$ TeV . Another idea is that grand unification occurs precociously at $`\sqrt{s}`$ TeV, because of extra dimensions or other reasons, and a neutrino above this threshold becomes strongly-interacting via leptoquark resonances . A third idea is that the exchange of a towers of Kaluza-Klein (KK) modes from extra compactified dimensions lead to a strong neutrino cross-section above $`\sqrt{s}`$ TeV . In all three cases, it is the combination of a low $``$ TeV scale for radically new physics and a quickly rising spectrum of new states (possibly increasing exponentially, $`\rho e^{\sqrt{s/s_0}}`$) that provides a rapid turn-on of a strong cross-section for the neutrino. Through unitarity, the new threshold at $``$ TeV has consequences for cross-sections at lower energies , but they are not dramatic.

The KK exchange model may fail  in that the KK modes couple to neutral currents, and the scattered neutrino carries away 90% of the incident energy per interaction, thereby elongating the shower profile. But if the neutrino cross-section can be made large enough, $`\stackrel{>}{}20`$ mb, then multiple scattering within a nucleus will effect a sufficiently large energy transfer and save the model . Independent of the neutral current issue, the dual QCD and TeV-scale unification models seem to provide viable explanations of the super-GZK data. However, a recent calculation for the rate of rise of the low-scale unification cross-section in a string context is not encouraging .

Signatures for these models include directional pointing back to the EECR source, longitudinal shower profiles differing somewhat from those of a proton or a gamma, and a strong correlation between observed energies and zenith angle. The latter signature should show an inverse proportionality between the neutrino-air interaction length and the rising (with energy) neutrino cross-section.

## IV MODEL SIGNATURES

There are several telltale discriminators to be sought in higher statistics data. These will eventually eliminate most (perhaps all!) of the models so far proposed for the super-GZK events. We list some discriminating signatures and discuss them.

### A Small-scale anisotropies and pointing

The discriminatory power of small-scale clustering was discussed already in §II A. Here we add some detail to the discussion. In traversing a distance $`D`$, a charged particle interacting with magnetic domains having coherence length $`\lambda `$ will bend through an energy-dependent angleOn average, half of the interactions of a super-GZK nucleon with the CMB change the isospin. At energies for which $`c\tau `$ of the neutron is small compared to the interaction mfp of $`6`$ Mpc, the neutron decays back to a proton with negligible energy loss and the bending-angle formula is unchanged. However, at the energy $`6\times 10^{20}`$, $`c\tau `$ for the neutron is comparable to the interaction mfp, so at higher energies the nucleon bending-angle is reduced by $`\stackrel{<}{}2`$.

$$\delta \theta 0.5^{}\times \frac{ZB_{nG}}{E_{20}}\sqrt{D_{\mathrm{Mpc}}\lambda _{\mathrm{Mpc}}}.$$
(13)
Here $`B_{nG}`$ is the magnetic field in units of nanogauss, $`E_{20}`$ and $`Z`$ are the particle energy in units of $`10^{20}`$ eV and charge, and the lengths $`D`$ and $`\lambda `$ are given in units of Mpc. It is thought likely that coherent extragalactic fields are nanogauss in magnitude , in which case super-GZK primaries from $`\stackrel{<}{}50`$ Mpc will typically bend only a few degrees (but note that protons at $`10^{19}`$ eV will bend through $`30^{}`$). Thus, local models either postulate many invisible sources isotropically-distributed with respect to the Galaxy to provide the roughly isotropic flux observed above $`E_{\mathrm{GZK}}`$, or postulate a large extragalactic magnetic field to isotropize over our Northern Hemisphere the highest-energy particles from a small number of sources . Among the latter category, some models postulate helium or iron nuclei as the primaries, to increase the bending by the charge factors 2 and 26, respectively. For those models invoking randomly distributed, decaying super-massive particles (SMPs)  or topological defects (TDs)  as sources, and models invoking a large magnetic field with considerable incoherent component, one may expect a nearly chance distribution of observed events on the sky. However, there may be some clustering even in these models, due to possible small-scale density fluctuations in the local SMP or TD distributions , or due to possible caustics in the projection of large-scale extragalactic magnetic fields on our sky (assuming the incoherent magnetic fields are sufficiently small) . In the SMP and TD models, a high photon fraction in the primary composition further enhances clustering possibilities.

From the point of view of opening a new window to astronomy, those models in which the primaries do point back to their distant, active sources are the most interesting. These models are few in number. They are the Z-burst model, the strongly-interacting neutrino model, and the quantum-gravity/LI-violating models.

### B Large-scale anisotropies

On large scales, one seeks associations of the CR directions with the Galactic halo (to be revealed by a dipole anisotropy favoring the direction of the Galactic center) or the local galactic magnetic field, with matter distributions in nearby galactic or super-galactic clusters such Virgo, or with possible large coherent galactic or extragalactic magnetic fields . For large-scale studies, the Southern hemisphere Auger experiment will prove invaluable for several reasons. It offers coverage of potential sources and matter distributions, and galactic and extragalactic magnetic fields, not available from the North. Moreover, it offers a view of our Galactic center which will provide a North-South dipole discriminator for or against a halo-centered population of sources such as magnetars or halo-bound SMPs. Southern Auger will also discriminate the M87-source model  wherein EECRs are channeled by a hypothesized galactic magnetic wind into the Northern hemisphere, and the Cen-A source model  which also yields a dipole anisotropy. Of course, an orbiting experiment with $`4\pi `$ vision like EUSO/OWL will be an even better instrument for multipole analysis.

### C Energy-direction-time correlations

Because bending of charged particle trajectories by intervening magnetic fields increases as particle energy decreases, one may learn about the strength and geometry of extragalactic fields from relative time delays and angular correlations of particles from a common source. One may also learn about the source. Quantitatively, the increase in path length due to bending leads to a relative increase in travel time of $`\delta t/t(\delta \theta )^2`$, for small bending angle. Adding the contributions from the coherent magnetic domains then yields

$$\delta t300D_{\mathrm{Mpc}}\left(\frac{ZB_{nG}\lambda }{E_{20}}\right)^2\mathrm{yrs}$$
(14)
for the time delay. The time separation at earth is obtained by taking differences in eq. (14); to first order in $`\delta E`$ it is already large:

$$t_1t_2600D_{\mathrm{Mpc}}\left(\frac{ZB_{nG}\lambda }{E_{20}}\right)^2\left(\frac{\delta E}{E_{20}}\right)\mathrm{yrs}.$$
(15)
The correlation in energy and time becomes even more significant when it is remembered that the higher energy primary has an even higher mean energy in transit, before losses on the 2.7K background.

Surprisingly, one AGASA event-pair has the higher-energy $`1.06\times 10^{20}`$ primary arriving about 3 years after a $`0.44\times 10^{20}`$ primary. Assuming these primaries originate form a common source, a possible explanation is that their source has a duration of at least $`3/(1+z)`$ years (the red-shift factor is due to time-dilation). Such a long-lived source does not occur in one-time burst models (e.g. GRBs). It may occur in decay models (e.g. SMPs, TDs, monopolonium) if sub-clustering exists on small-angular scales withing halo clusters . A recent paper  notes that such sub-clusters may actually be observable as micro-lenses for stars and background galaxies. Counter-intuitive reverse pairing with the earlier arrival time for the lower energy charged-primary can also occur in certain magnetic field configurations, as shown in .

### D Composition of the primaries

Another signature to be sought is a statistical identification of the nature of the primaries as a function of their energies. Three methods have been identified to distinguish photon-initiated showers from hadron-initiated showers. One method relies on the longitudinal profile of the event, particularly the depth at shower maximum $`x_{\mathrm{max}}`$. The longitudinal profile of the Fly’s Eye event at $`3\times 10^{20}`$ eV (well-measured by its nitrogen fluorescence trail) is ill-fit by a photon primary, well-fit by an iron primary, and somewhat fit by a proton primary . The second method relies on measurement of muon number, with a high muon number purportedly favoring a hadron primary over a photon, and vice versa. A recent study of the muon content of showers above $`10^{19}`$ eV seems to favor nucleons over photons . However, caution is warranted with this method, in that some simulations show little difference in the muon-content of showers from nucleon primaries versus photon primaries ; and photons have themselves a significant partonic component at high energy. Nevertheless, it is safe to say that photon primaries appear disfavored at present, but more data is needed before conclusions are drawn.

The third method of gamma identification will will rely on a predicted characteristic N-S vs. E-W gamma asymmetry. This quadrupole asymmetry results from the polarization-dependent interaction of the gammas with the earth’s magnetic field.

It is interesting to mention that the highest-energy Yakutsk event has an unusually high muon number; the only model so far which successfully explains this invokes a magnetic baryonic-monopole as the primary particle . Basically, a relativistic baryonic-monopole of mass $`M`$ and incident energy $`E`$ showers like a giant nucleus of atomic number $`A_{\mathrm{eff}}M/\mathrm{\Lambda }_{\mathrm{QCD}}10^6(M/\mathrm{PeV})`$ and the same $`\gamma =E/M`$. Intermediate mass monopoles therefore generate many, many charged pions which decay to muons.

All models wherein the primaries arise from QCD jets produce many more pions than baryons. The neutral pions in turn produce gammas. In these models, the ratio of gammas to baryons is typically of order ten at the source. Even allowing for the shorter attenuation length of photons relative to nucleons, the measurement of the primary compositon for super-GZK events becomes an excellent discriminator for models with jet-producing sources. These models include $`Z`$-bursts, and decaying SMPs and TDs.

### E Possible $`E_{\mathrm{max}}`$ energy cutoff

The predicted of a cutoff at $`E_{\mathrm{GZK}}`$ is wrong. Does Nature provide an alternative cutoff within our reach? Or do the data continue beyond our reach? The shock-jock experts claim that it is difficult for conventional shock-acceleration mechanisms to produce ZeV proton energies; for this class of model, an $`E_{\mathrm{max}}`$ below $`10^{22}`$ eV is certainly expected . Decaying SMP models also have a natural cutoff, at half of the SMP mass. This could be as high as $`E_{\mathrm{max}}10^{24}`$ eV for a long-lived GUT-mass particle, but would be lower for other postulated SMPs. In the $`Z`$-burst model there is a natural cutoff related to the tiny mass of neutrinos: $`E_{\mathrm{max}}=4(\mathrm{eV}/m_\nu )`$ ZeV. Implications from the atmospheric SK data are that this cutoff is at most $`7\times 10^{22}`$ eV.

### F CR flux above vs. below $`E_{\mathrm{GZK}}`$

A “smoothness” variable such as $`R_jF_j(E>E_{\mathrm{GZK}})/F_j(E<E_{\mathrm{GZK}})`$ for each primary species $`j`$=nucleon, photon, iron nuclei, neutrino, etc. may be revealing . For primaries with a GZK cutoff, $`F(E<E_{\mathrm{GZK}})`$ samples sources from the whole volume of the Universe, and may even include cascade products from $`F(E>E_{\mathrm{GZK}})`$, whereas $`F(E>E_{\mathrm{GZK}})`$ samples just the GZK volume; for primaries without a GZK cutoff, $`F(E<E_{\mathrm{GZK}})`$ and $`F(E>E_{\mathrm{GZK}})`$ sample sources from the whole volume of the Universe. Lumps, bumps, and gaps in the spectrum near $`E_{\mathrm{GZK}}`$ are a consequence of some models. For example, hadron and neutrino pile-ups just below $`E_{\mathrm{GZK}}`$ are expected from the photo-pion production process which occurs above $`E_{\mathrm{GZK}}`$ . Other sources of neutrino pile-ups have also been suggested . Present data show continuity. Smoothness studies of various models require simulation, and are just beginning.

### G Spectral index above $`E_{\mathrm{GZK}}`$

One means of achieving more events above $`E_{\mathrm{GZK}}`$ is to postulate a flattening of the primary proton spectrum at highest energies. With more data, the extreme-energy spectrum will be measured. A flattening of the $`E^{2.7}`$ power law inferred from just below $`E_{\mathrm{GZK}}`$ would indicate and constrain new sources.

### H Measureable neutrino flux above $`E_{\mathrm{GZK}}`$

If there is a new source of primary neutrinos above $`E_{\mathrm{GZK}}`$, or if the neutrino cross-section becomes strong at super-GZK energies, then there is the possibility that the primary neutrino flux can be measured in an EUSO/OWL/AW-sized detector. This possiblity was discussed in section III D.

### I Diffuse $``$GeV gamma-ray flux

An upper limit on the diffuse gamma-ray flux between 30 MeV and 100 GeV has been published by the EGRET experiment of the now defunct Gamma Ray Observatory . This limit has serious discriminatory power for models where SMPs or TDs or extremely boosted massive particles decay to quark-antiquark jets which then hadronization to produce the EECRs . This is because QCD jets via $`\pi ^0`$ production and decay produce very high energy gammas which initiate an electromagnetic cascade on the cosmic radio, microwave, infrared, and magnetic field backgrounds. For models with jet-production distributed over cosmic distances, such as some TD models and the $`Z`$-burst model with a homogeneous distribution of relic neutrino targets, the cascade has the distance to reach completion; the end result is gamma power (energy/time) in the EGRET range roughly an order of magnitude below the total power of the original sources, and comparable to the power in EECR neutrinos. Such models are disfavored. More local models, such as SMPs bound to our halo, and the $`Z`$-burst model with a local over-density of relic neutrinos, are not impacted at present, but may be tested with GLAST, the next generation gamma-ray observatory.

## V CONCLUSIONS

The ultimate explanation for the puzzles in EECRs will provide a surprise at a minimum, and possibly radically new physics at a maximum. The extreme energies of events already observed cannot be approached by terrestrial accelerators. Thus, there is ample motivation to build the next generation of CR detectors, Auger and Telescope Array, and to plan even farther beyond for a teraton detector like EUSO, OWL and AirWatch. At $`10^{20}`$ eV, AGASA provides about an event per year, and HiRes about an event per month. Auger and TA will see two such events per week, while EUSO/OWL/AW has the potential to collect such an event every two hours. As the data sample grows, statistical studies will reveal signatures that discriminate among the many galactic and extragalactic sources so far proposed to resolve the EECR puzzles.

It is quite possible that neutrino primaries are responsible for the EECRs. If so, it appears that the weakly-interacting neutrino either grows a very strong cross-section at $`10^{20}`$ eV, or it annhilates on the relic neutrinos left-over from the hot phase of the Big Bang. Another possibility is that free magnetic charges exist and are the EECR primaries; for magnetic monopoles, the Universe is transparent and cosmic magnetic fields provide a natural acceleration mechanism.

Clearly, we live in exciting EECR times, and we possess the technology to prove this is so. The resolution of our puzzles is forthcoming, as on-going and future experiments will provide us with the statistics to discriminate among the many interesting models.

## Acknowledgements

We acknowledge fruitful collaboration with P. Biermann, H. Goldberg, T. Kephart, H. Paes, and Dipthe Wick on some of the material presented here; and discussions with G. Farrar, G. Gelmini, A. Kusenko, D. McKay, A. Olinto, J. Ralston, G. Sigl, F. Stecker, and E. Zas. This work was supported in part by the U.S. Department of Energy grant no. DE-FG05-85ER40226.


text:
# Rheology of a confined granular material

## Abstract

We study the rheology of a granular material slowly driven in a confined geometry. The motion is characterized by a steady sliding with a resistance force increasing with the driving velocity and the surrounding relative humidity. For lower driving velocities a transition to stick-slip motion occurs, exhibiting a blocking enhancement whith decreasing velocity. We propose a model to explain this behavior pointing out the leading role of friction properties between the grains and the container’s boundary.
PACS numbers: 46.55.+d, 45.70.-n, 81.05.Rm

Dynamics of granular materials lacks of an established unified picture. A great diversity of mechanical and rheological behaviors were reported depending on whether they are vibrated, slowly sheared, avalanching on a surface, flowing in a hopper or falling in a chute. For dense granular assemblies, experiments reveal strain localization , , non local rheological properties and aging phenomena,. This complex phenomenology could possibly be due to the presence, at a mesoscopic level, of a disordered contact force network with unraveled mechanical properties ,. Moreover, dramatic effects on the mechanical strength induced by slight changes in compaction were reported, ,. Another source of difficulty is the dynamical contribution of contact forces and there is, so far, only few studies on the macroscopic emergence of these aging properties due to the slow plastic deformation of contacts and the influence of surrounding humidity ,,,.

Dynamical behavior of slowly driven granular materials was investigated recently by different groups both in compression and/or in shear experiments ,,,,. Here, we investigate the rheology of a granular assembly confined in a cylindrical column and pushed vertically from the bottom. The resistance to vertical motion as well as the blocking/unblocking transitions, reveal a phenomenology possibly shared by many confined granular assemblies, as for example, gouge sheared between two faults, pipe flows, compaction under stress or dense granular paste extrusion. A previous investigation of the same display but in 2D has already shown a rich phenomenology. For 3D granular assemblies, we get similar behavior but in this letter, we choose to report only on the most simple dynamical situation involving low friction grains confined in a column with rather frictional walls. Here, we have a weak coupling with dilatancy effects due to shearing at the walls and therefore, relations between the granular nature of the bulk (i.e. the stress redirection) and the solid friction properties at the walls are the most clearly revealed.

The grains are dry, non cohesive and monodisperse steel beads of diameter $`d=1.58`$ $`mm`$ piled into a vertical cylinder in duralumin of diameter $`D=36`$ $`mm`$. The column is closed at the bottom by a movable piston avoiding contact with the column (diameter mismatch is $`0.5`$ $`mm`$). A force probe of stiffness $`k=40000`$ $`N.m^1`$ is located under the piston and is pushed at a constant driving velocity $`V`$ (between $`5`$ $`nm.s^1`$ and $`100`$ $`\mu m.s^1`$) via a stepping motor (see inset of fig.1a). The resistance force $`F`$ encountered by the piston is measured as a function of time. We monitor also the relative humidity ($`RH`$) and the surrounding temperature. The central question here is an attempt to estimate the relative influence of the bulk mechanical properties with respect to frictional properties of the walls. To address directly this issue we built a special device (the slider) designed to apply a constant normal load ($`F_N=2`$ $`N`$) on three steel beads sliding vertically on the cylinder’s wall (see inset of fig.3a). Then, the dynamical evolution of the resistance force encountered by the piston pushing the grains is compared with the slider’s friction resistance driven in the same conditions. We observe two distinct regimes(fig.1a): for high driving velocities, the motion is characterized by a steady-sliding and a constant pushing force; for low velocities, the system undergoes a dynamic instability and then a stick-slip motion occurs. The transition between these behaviors is similar to the inertial regime of Heslot et al.  and details will be reported elsewhere.

For a vertically pushed granular assembly, the driving force exerted by the piston is screened. To evaluate this effect, the mean resistance force is measured as a function of the packing height (see fig.1b). For this dataset the driving velocity $`V`$ corresponds to a steady and continuous sliding of the grains. The resistance force $`F`$ increases very rapidly with the column’s height $`H`$. This strong resistance to motion is due to the horizontal redirection of stresses in association with solid friction at the side walls. Following the standard Janssen’s screening picture, , the force $`F`$ exerted by the grains on the piston can be modelled by the relation:

$$F_ϵ=\varrho g\lambda \pi R^2\times ϵ(\mathrm{exp}(ϵ\frac{H}{\lambda })1)$$
(1)
where $`\varrho `$ is the mass density of the granular material, $`R`$ is the cylinder radius and $`g`$ the acceleration of gravity. The length $`\lambda =R/2K\mu `$ is the so-called effective screening length, where $`K`$ is the Janssen’s parameter rendering the average horizontal redirection of vertical stresses and $`\mu `$ can either be the dynamic or the static coefficient of friction of beads at the cylinder’s wall. When $`ϵ=+1,`$ friction is fully mobilized downwards (our pushing experiment) and when $`ϵ=1,`$ friction is fully mobilized upwards. It is easily seen from (1) that when $`ϵ=+1`$, any slight change in $`\mu `$ or $`K`$ is exponentially amplified with a drastic influence on the pushing force $`F`$. In the case of steel beads, we found that, starting from a dense or a loose packing, the final average steady state compacity $`\stackrel{\mathrm{\_}}{\nu }`$ does not change ; we have $`\stackrel{\mathrm{\_}}{\nu }62.5\%`$ for all velocities and relative humidities $`RH`$ tested. In the steady state regime, the experimental data obtained for a given pushing velocity $`V`$ can be fitted with relation (1) by adjusting only one parameter i.e. $`p_{+1}=K\times \mu _d`$ where $`\mu _d`$ is the dynamic coefficient of friction at velocity $`V`$. For the relative humidity $`RH=42\%`$, we obtain $`p_{+1}=0.140\pm 0.001`$ at $`16`$ $`\mu m.s^1`$ and $`p_{+1}=0.146\pm 0.001`$ at $`V_{up}=100`$ $`\mu m.s^1`$. As a check of consistency, we performed the following dynamical experiment. First, the granular column is pushed upwards in order to mobilize the friction forces downwards and to reach the steady state compacity. Starting from this situation, the friction forces are reversed at the walls by moving the piston downwards at a constant velocity $`V_{down}=16\mu m.s^1`$. Following relation (1), this procedure should imply a change of $`ϵ`$from $`1`$ to $`1`$, and consequently, the dynamical force on the piston should decrease from $`F_{+1}`$ to $`F_1`$. In fig.1b the pushing force $`F_1`$ is measured for different packing heights $`H`$. Injecting the preceding value of $`p_{+1}(16`$ $`\mu m.s^1)`$ into (1) with $`ϵ=1,`$ we check on fig1b, that the theoretical expectation agrees quite well with the experimental data of $`F_1`$ versus $`H`$. Note that in a previous study it was found that the Janssen’s picture has a general tendency to underestimate the stress below a granular column . But in the present situation, with low friction steel beads, this model though elementary, seems a good base for analysis. Nevertheless, a question is still that the fitting parameters $`p=K\times \mu `$ extracted from the model are unable to sort between what comes out from wall-bead interactions ($`\mu `$) and what comes out from bulk properties ($`K`$). Actually, from series of static Janssen experiments we extracted $`K\mu _s`$ . Independently, the static coefficient of friction $`\mu _s`$ of our steel beads on duralumin was measured (in the short time limit) using the sliding angle of a three beads tripod. We get $`\mu _s=0.170\pm 0.005`$ and $`K=1.08\pm 0.05`$ is extracted from this procedure. This $`K`$ value is consistent with previous measurements done on a granular column at this compacity. Note that, if we tentatively assume a constant value for $`K`$ in static and dynamic experiments, the dynamical coefficient $`\mu _d`$ can be extracted from the measured values of $`p_{+1}`$. For instance, for $`RH=42\%`$ we extract the values, $`\mu _d(16\mu m.s^1)=0.130\pm 0.005`$ and $`\mu _d(100\mu m.s^1)=0.135\pm 0.005`$. This would imply a slight increase of the bead/wall friction with the driving velocity. This result is going to be directly tested in the following, using the three beads slider device.

For a given height $`H`$ of beads ($`M=380g`$ i.e. $`H=4.3R`$), we study extensively how the pushing force depends on the driving velocity $`V`$ and on the surrounding relative humidity ($`RH`$). We worked in the range $`35\%<RH<75\%`$, and also in dry air ($`RH<3\%`$). Note that except for the dry situation, we did not have a mean of regulation of this last parameter ($`RH`$) but we record its values around the experimental set-up. All data are shown on a same series for similar humidity values (within $`3\%`$). As it was already mentioned, the motion is characterized by a steady sliding above a critical velocity (fig.1a). The mean force level in this regime increases slowly with velocity but surprisingly strongly with $`RH`$ (fig.2a). For example at $`V=100`$ $`\mu m.s^1`$, the resistance force is raised by $`35\%`$ for a change of $`RH`$ from $`53\%`$ to $`72\%`$. Now we perform the same series of experiment but with the three beads slider, in order to test directly the wall/bead friction properties. At a given $`RH`$, we indeed observe velocity strengthening for the sliding of individual steel beads, corroborating qualitatively the general trend observed on the granular column. But now, we may go one step further by testing directly the possibility of a quantitative agreement within the Janssen’s model. If we compare these data to the values of $`\mu `$ extracted from (1) (see fig.3a), assuming the static value of $`K=1.08\pm 0.05`$ for all velocities, we observe that the increase of $`\mu `$ with $`V`$ is significantly less important in the case of the granular column than what is directly measured using the slider device. Actually, if we suppose a logarithmic increase of $`\mu `$ with velocity $`V`$, i.e. $`\mu b\times \mathrm{log}(V)`$, we find, at $`RH=40\%,`$ $`b=(2.7\pm 0.2)\times 10^2`$ for the slider and $`b=(1.2\pm 0.2)\times 10^2`$ in the case of the granular column; in dry conditions, we find $`b=(2.4\pm 0.2)\times 10^2`$ for the slider and $`b=(1.1\pm 0.2)\times 10^2`$ in the case of the granular column. It means that the increase of $`F`$ with $`V`$ cannot be entirely attributed to friction effects at the walls, and that the dynamics may have also an effect on force transmission (i.e. on $`K`$). In the framework of a Janssen’s analysis it would mean that $`K(V)`$ would slightly decrease when velocity increases. Using a simple Hertz law to estimate contact interactions, we find the depths of penetration of steel beads in duralumin to be around $`\delta 30nm`$, whereas in the slider case, we estimated $`\delta 1\mu m`$ which is the order of duralumin roughness. Therefore, it is also possible that contacts are not both in the same loading regime and then friction laws could be slightly different. Importantly, we have also found that an increase of humidity has quite a strong influence on the friction properties (fig.2a). Using the inverted Janssen’s model (eq. (1)), when assuming the redirection parameter $`K`$ unchanged by humidity, we recover that the dependence on $`RH`$ in the slider experiment is consistent with the enhancement of the friction forces measured in the granular column. In a future series of experiments, we will try to bridge the gap to controlled values of humidity close to $`100\%`$.

Now let us consider slow driving velocities where the system undergoes a dynamical instability. A stick-slip motion occurs (see fig.1a) with a narrow gaussian distribution of slip force amplitudes. In fig.2b, we display the mean maximum and mean minimum resistance forces (resp. $`F_{\mathrm{max}}`$ and $`F_{\mathrm{min}}`$) as a function of the driving velocity, for $`m=380g`$ of beads (height $`H=7.7cm`$), and relative humidity $`RH=45\pm 3\%`$. The mean amplitude of the slip events $`\mathrm{\Delta }F=F_{\mathrm{max}}F_{\mathrm{min}}`$ increases strongly when velocity reaches values as small as $`5`$ $`nm.s^1`$. We propose a model where this enhanced blocking effect can be simply interpreted by an aging effect of the contacts at the side walls. Friction coefficients of solid on solid contacts are known to evolve logarithmically with waiting time $`t`$ : $`\mu _s(t)=\mu _s^0+\beta _s\mathrm{log}(t)`$. According to fig.2a, we observe no noticeable variation of $`F_{\mathrm{min}}`$ with velocity, for given height and $`RH`$. Therefore, we will consider in the following $`F_{\mathrm{min}}`$ to be a constant. Starting at the onset of blocking $`t=0`$, the force exerted by the force probe during a stick event is $`F(t)=F_{\mathrm{min}}+kVt`$. The time elapsed during a sticking event is:

$$t_{stick}=\frac{F_{\mathrm{max}}F_{\mathrm{min}}}{k\times V}.$$
(2)
The slip occurs when $`F(t)`$ reaches the maximum force sustainable by the granular material at time $`t`$, given by (1) with $`ϵ=+1`$. The aging properties of the friction at the wall are included in the time evolution of the static coefficient of friction $`\mu _s(t).`$ Then we write $`F(t_{stick})=F_{\mathrm{max}}`$, i.e.:

$`F_{\mathrm{max}}`$ $`=`$ $`{\displaystyle \frac{\varrho g\pi R^3}{2K(\mu _s^0+\beta _s\mathrm{log}(t_{stick}))}}`$ (3)
$`\times (\mathrm{exp}(2K(\mu _s^0+\beta _s\mathrm{log}(t_{stick})){\displaystyle \frac{H}{R}})1).`$
This exponential amplification of the logarithmic aging, due to stress redirection at the walls, gives an effective power-law : $`F_{\mathrm{max}}t_{stick}^\alpha `$, with $`\alpha ={\displaystyle \frac{2\mathrm{log}(\text{e})HK\beta _s}{R}}`$. On fig.3a , for $`RH=45\pm 3\%`$, we display $`\mu _s`$ extracted from (3) as a function of the time of stick. We assume $`K=1.08`$ independent both of the waiting time and the driving velocity. We actually observe a logarithmic aging for waiting times $`3000s`$, with a coefficient $`\beta _s=1.8\times 10^2\pm 2\times 10^3`$, value consistent with many previous reports; in the last decade, aging is strongly increased and we have $`\beta _s6\times 10^2`$. Note that our experiment is not a ”clean” aging experiment since the applied loads and the shear forces are not constant in time and along the vertical direction. Furthermore and consistently with the finding of refs  and , we clearly observe that the aging properties are strongly affected by a variation of the relative humidity $`RH`$ (see inset of fig.2a).

In conclusion, we investigated the dynamical behavior of a granular column pushed vertically from the bottom. This model experiment is suited to understand the rheology of slowly driven granular assemblies in confined geometries. Overall, the pushing force data are analyzed consistently using an inverted Janssen’s law. At such slow driving velocities we show that, all the non trivial dynamical properties exhibited by the granular rheology (including a strong dependence on relative humidity) can be dominantly attributed to the dynamical properties of solid on solid friction. In addition, the model seems to indicate the presence of a dynamical structural effect induced in the bulk at higher driving velocities.

We acknowledge many discussions with Profs R.P Behringer, C.Caroli, T. Baumberger and the financial support of the CNRS-PICS grant $`\mathrm{\#}563`$.


text:
# Introduction

## Introduction

The phenomena of Glass transition has been studied extensively through theoretical  models, computer simulations  as well as scattering experiments . If a liquid is cooled beyond the freezing point avoiding the first order transition to an ordered crystalline state, it continues in the amorphous state with solid like properties. In this supercooled state the liquid develops very long relaxation times and the transport coefficients (viscosity) abruptly increase. A microscopic mechanism for this slow dynamics using a statistical mechanical approach has been proposed through the Mode coupling theory (MCT). This is based upon the idea  that the transport coefficients in the supercooled state gets a feedback from the nonlinear coupling of slow modes in the fluid. The effect of the structure on the dynamics is studied from the wave vector dependence of the model. MCT predicts a two step relaxation process for glassy dynamics : power law ($`\beta `$-relaxation) at intermediate time scales which crosses over to stretched exponential ($`\alpha `$-relaxation) decay of density correlation function over long time scales. In the ideal model  there is a sharp transition from ergodic to a nonergodic phase at a critical density. In the later version of more careful analysis  the sharp transition is finally removed from the ergodicity restoring mechanisms coming from the coupling of the density fluctuations with currents in the compressible fluid. The final ergodicity restoring mechanisms in a supercooled liquid behaves like a collective slow mode. There has been phenomenological extension of the theory to consider extra slow modes that develop in the liquid in the supercooled region. In the present work we consider an extension  of the simple feed back mechanism of ideal MCT through the inclusion of the extra slow mode of defect density in the hydrodynamic description for the isotropic liquid. Here we use this model to study the dynamics in the simple One component Lennard-Jones fluid in the supercooled state. In this respect we like to note that previous MCT studies on LJ liquid  show that one component system freezes into a nonergodic phase at a relatively low density.

The defect density considered in the model can be associated with the vacancies in the solids or void or free volumes in the fluid. The movement of these free volumes were considered crucial for the transport in the liquid . Following Ref.  we consider the defects to be moving in a metastable potential well. Here the extent of coupling between defect density and the particle density governs the location of the ideal glass transition point that occurs as a result of the feedback mechanism proposed in MCT. Thus at a constant pressure we observe a transition line whereas a transition point is predicted in the simple MCT. Phase diagram of coupling parameters with temperature is studied through nonergodicity parameter. We compute the relaxation dynamics of fluid in the power law regime and calculate the power law exponents. From the present model the temperature dependence of the relaxation exponents arises naturally as a consequence of the critical line. The paper is organized as follows. In next section we briefly describe the model considered followed by section III where we discuss briefly the static strucutre factor used in the calculation. In section IV the results for the power law relaxation as predicted by the mode coupling theory is presented. In last section a brief discussion of the work is presented.

## The Model Studied

In the mode coupling theory we focus on the normalized density-density correlation function $`(\psi _q(t))`$. The dynamical evolution of $`\psi _q(t)`$ is written in terms of an integro-differential equation,

$$\mathrm{\Delta }_q\ddot{\psi }_q(t)+\mathrm{\Gamma }_q^0\dot{\psi }_q(t)+\psi _q(t)+_0^t𝑑s\mathrm{\Gamma }_q^{mc}(ts)\dot{\psi }_q(s)=0,$$
(1)
where the constant $`\mathrm{\Delta }_q=\beta mS(q)/q^2`$ with $`m`$ as the mass of the liquid particle and $`\beta `$ is the Boltzman factor. S(q) is the static structure factor for the liquid. The bare transport coefficient $`\mathrm{\Gamma }_q^0`$  is defined in a dimensionless form as

$$\mathrm{\Gamma }_q^0=\frac{2\mathrm{\Delta }}{3t_E}\left[1j_0(q^{})+2j_2(q^{})\right],$$
(2)
where $`q^{}=q\sigma `$ and $`j_n(x)`$ is the $`n^{th}`$ order spherical Bessel function . The Enskog time $`t_E`$  for the LJ fluid is approximated using the equivalent hard sphere diameter. The generalized transport coefficient $`\mathrm{\Gamma }_q(t)`$ has the bare ($`\mathrm{\Gamma }_q^0`$) as well as the mode coupling ($`\mathrm{\Gamma }_q^{mc}(t)`$) contributions from the nonlinear coupling of density fluctuations. Following Ref. , the mode coupling part of the generalized transport coefficient is given by,

$$\mathrm{\Gamma }_q^{mc}(t)=\frac{d\stackrel{}{k}}{(2\pi )^3}[\stackrel{~}{V}^{(1)}(q,k)\psi _k(t)\varphi _{k_1}(t)+\stackrel{~}{V}^{(2)}(q,k)\psi _k(t)\psi _{k_1}(t)]$$
(3)
with $`\stackrel{}{k_1}=\stackrel{}{q}\stackrel{}{k}`$. Here $`\varphi _{k_1}(t)`$ is the normalized defect auto correlation function and gives the contribution to the transport coefficient coming from the coupling of the defect density with the density fluctuations. The mode coupling vertices with full wave vector dependence are given by

$$\stackrel{~}{V}^{(1)}(q,k)=\frac{S(q)}{n_0}\left[2y\stackrel{~}{U}(q,k)S(k)S(k_1)+\kappa n_0c(k)S(k)\right]$$
(4)
and

$$\stackrel{~}{V}^{(2)}(q,k)=\frac{S(q)}{2n_0}\left[\stackrel{~}{U}^2(q,k)2y\stackrel{~}{U}(q,k)\right]S(k)S(k_1)$$
(5)
These equations are obtained  from the nonlinear equations for the extended set of hydrodynamic variables. The vertex function $`\stackrel{~}{U}`$ in (5) is given by,

$$\stackrel{~}{U}(q,k)=\frac{n_0}{q}[\widehat{q}.kc(k)+\widehat{q}.(k_1)c(k_1)]$$
(6)
where $`n_0`$ is the particle number density. It is assumed that the defects are moving in a metastable double well potential with $`\kappa `$ and $`y`$ as two dimensionless parameters characterizing the coupling of the defect densities with the particle density and the depth of the potential well. Defect density is assumed to be weakly interacting with the mass density and is treated as a variable similar to mass density. For $`\kappa =y=0`$, (3) reduces to the standard result, . The presence of the linear term in the mode coupling vertex comes through the coupling with the very slowly varying defect correlation . In the next section we study the implications of this feed back mechanism that arises from the coupling of the density fluctuations as well as the defect density coupling. For this the static structure factor is taken as an input - measuring the extent of structural effects on relaxation. We use the standard model  for the Lennard-Jones Structure factor. We give a brief description of computation of the static structure factor next.

## Static Structure Factor

In mode coupling equations interparticle interaction potential enters as an input via static structure factor $`S(q)`$. For a Lennard-Jones (LJ) fluid the perturbation approach described in  is being used for the computation of the structure factor. The full Lennard-Jones potential,

$$U(r)=4ϵ[(\sigma /r)^{12}(\sigma /r)^6]$$
(7)
is divided into a repulsive part $`U_0(r)`$ and an attractive part $`U_A(r)`$ as,

$$U_0(r)=\{\begin{array}{cc}U(r)+ϵ,\hfill & r<r_0\hfill \\ 0,\hfill & rr_0\hfill \end{array}$$

$$U_A(r)=\{\begin{array}{cc}ϵ\hfill & r<r_0\hfill \\ U(r),\hfill & rr_0\hfill \end{array}$$

with $`r_0=2^{1/6}\sigma `$. A trial system with hard sphere diameter $`dr_0`$ and an attractive potential energy part $`U_A(r)`$ is taken. One then finds the appropriate diameter $`d`$ through perturbation approach until $`e^{\beta U_T(r)}g_T(r)=e^{\beta U(r)}g(r)`$. Where, subscript $`T`$ represents the trial system. Once the hard sphere diameter $`d`$ is known the static structure factor for the LJ liquid is approximated by,

$`S(q)=S_T(q)+\rho {\displaystyle 𝑑𝐫[e^{\beta U_T(r)}g_T(r)][e^{\beta U(r)}e^{\beta U_T(r)}]e^{i𝐪.𝐫}},`$ (8)
For the trial hard sphere system the Percus-Yevick  solution with Verlet-Wiess correction  is used.

## Power law exponents

The ideal transition from an ergodic phase (fluid) to a non-ergodic (glass) phase is characterized in terms of non-ergodicity parameter defined as, $`\psi _q(t\mathrm{})=f(q)`$. In the long time limit or $`z0`$ limit of the Laplace-transformed form of the equation (1) reduces to

$$\frac{f(q)}{1f(q)}=\frac{d^3k}{(2\pi )^3}[\stackrel{~}{V}^{(1)}(q,k)f(k)+\stackrel{~}{V}^{(2)}(q,k)f(k)f(k_1)]$$
(9)
where vertices $`\stackrel{~}{V}^{(1)}`$ and $`\stackrel{~}{V}^{(2)}`$ are given by equation (4) and equation (5). These set of coupled equations are solved self-consistently for $`f(q)`$. In solving the integral equations a maximum cutoff for the wave vector integration $`\mathrm{\Lambda }\sigma =30.0`$ with 200 grid points is used. A set of critical parameters ($`\kappa _c,y_c,\rho _c,T_c`$) are identified below which all $`f(q)^{}s`$ converge to zero values (ergodic phase) and above that critical point all $`f(q)^{}s`$ converge to nonzero values (nonergodic phase). The critical points ($`y_c`$, $`\kappa _c`$, $`\rho _c`$, $`T_c`$) in $`T\kappa `$ phase space are located by fixing $`y_c=0.15`$ and varying $`\kappa `$ and $`T`$. Here we study the equations with the thermodynamic parameters changing in a way so that the volume change with the temperature keeping the Pressure constant, $`P^{}=P\sigma ^3/ϵ=0.5`$ as shown in Fig. 1. Then, at a fixed pressure, we observe a line of transition giving different critical temperatures with different values of the metastability parameters. In the ideal model under similar conditions a single transition point is obtained. Along this transition line at a fixed pressure we calculate, from the solutions of the MCT equations, the power law exponent for relaxation in the LJ fluid as a function of temperature and density. We study the relaxation over time scales longer than the microscopic time scales $`t>(\mathrm{\Omega }_q\tau )^1`$, where $`\tau `$ is the Lennard-Jones time $`(m\sigma ^2/ϵ)^{\frac{1}{2}}`$, but shorter compared to the $`\alpha `$-relaxation time scales. The dynamical equation (1), for density correlation function $`\psi _q(t)`$ is solved self consistently along the transition line shown in Fig. 1. The metastable potential well for defects used for computing vertex functions is shown in Fig. 2 at critical parameters ($`T_c=0.544,y_c=0.15,\kappa _c=0.038`$). The observed relaxation behavior is given by,

$$\psi _q(t)=f_c(q)+Ch(q)(t/t_o)^a$$
(10)
where C is a constant and $`t_o`$ is the time scale over which the power law relaxation persists. In Fig. 3 we show the decay of density correlation function with time ($`t/\tau )`$. The fit to the power law exponent is shown by dotted lines and the solid line represents the result from MCT equations for temperature $`T_c=.544`$ in the power law regime. Temperature dependence of the power law exponent $`a`$ is studied at a constant pressure $`P^{}=0.5`$. The exponent $`a`$ is calculated at the transition temperatures for fixed $`y_c=0.15`$ and the corresponding $`\kappa _c^{}s`$. Fig. 4 shows the variation of the power law exponent $`a`$ with temperature. We also calculate the wave vector dependence of the power law exponent $`a`$. Fig. 5 shows the exponent $`a`$ along the wave vector $`q\sigma `$ at $`T=0.544`$ and $`\rho =0.958`$. In the inset the static structure factor is shown to depict the alignment of the maxima of $`S(q\sigma )`$ with minima of exponent $`a`$ and vice-versa. The arrows are shown at the maxima and minima. The power law exponent is minimum at the peak of the static structure factor, which reflects the characteristic slow relaxation at the diffraction maximum for liquids.

## Discussion

In the simple mode coupling theory there is a single transition point with a fixed $`T_c`$. As a result of which the power law exponent determined in its vicinity is temperature independent. In the present case through the coupling of the density fluctuations with defect motions a line of dynamic instability is predicted and this results in a temperature dependent exponent. Here we consider the simple one component LJ system and we obtain an oscillatory behavior of the power law exponent a - as shown in Fig. 4. For hard sphere system a similar oscillatory behavior is found  from the solutions of the MCT equations. This was also shown to agree with the power law exponent that was calculated from fitting of the light scattering data on Colloids. A similar behavior was observed in the wave vector dependence of the exponent as well. Computer simulations of one component systems at metastable densities are particualrly rare due to problem of crystallization and at this time we are unable to find simulation data to look for such behavior. However simulation results for binary mixtures exist and we studied the data reported in ref.  without restricting to a fixed non-ergodicity parameter (NEP) as in the ideal model that gives rise to a single exponent. We do observe in this case a qualitatively same oscillatory behavior for the exponent $`a`$ of the simulation data as shown in Figure 6. The NEP that is obtained in our fitting also shows qualitatively similar behavior as predicted from the extended model studied here. The present model can be extended to the $`\alpha `$-relaxation regime as well where the combined role of the density fluctuation and the defect fluctuation need to be considered.

## Acknowledgement

SS acknowledges the support of University Grants Commission (UGC) of India for the research.

## Figure Captions

Fig 1 : Phase diagram in $`T\kappa `$ space at $`y_c=0.15`$ along a constant pressure line. Temperature is expressed in units of $`ϵ/k_B`$.

Fig 2 : Variation of depth of the potential with $`n/n^{}`$ for $`T_c^{}=0.544`$. In the figure $`h^{}(n)`$ represent the dimensionless quantity $`h(n)\beta ϵn\sigma ^3`$.

Fig 3 : Decay of normalized density-density correlation function $`\psi _q(t/\tau )`$ with time ($`t/\tau `$) at $`T=0.544`$. Where $`\tau `$ is the Lennard-Jones time $`(m\sigma ^2/ϵ)^{\frac{1}{2}}`$. Dotted lines show the power law fit (10) (see text).

Fig 4 Variation of power law exponent $`a`$ as a function of temperature $`T`$. $`T`$ is expressed in units of $`ϵ/k_B`$.

Fig 5 : Variation of exponent $`a`$ as a function of wave number $`q\sigma `$ for a packing fraction $`T=0.544`$. In inset static structure factor is shown along wave number $`q\sigma `$.

Fig 6 : Variation of power law exponent $`a`$ with temperature as obtained by fitting 10(see text) to the self correlation function data from Ref.  at wave vector $`q=7.22`$.