text:
# 𝐶𝑉𝐷 routes to 𝑀𝑔𝐵₂ conductors

## I Introduction

The discovery of superconductivity in $`MgB_2`$ at $`39.2K`$ opens some very interesting opportunities for practical applications of this material as a conductor in power applications. Typically, the operating temperature of these conductors would be at about half the superconducting transition temperature, $`T_c`$, or about $`20K`$ where closed cycle refrigerators are quite efficient. Early studies of the temperature and magnetic field dependence of the superconducting critical current density, $`J_c`$, indicated that grain boundaries in this material would carry very substantial supercurrents, $`5\times 10^5A/cm^2`$, in conductors prepared in a variety of ways. In addition, a very early attempt to explore the use of powder-in-tube methods by Jin and co-workers showed that practical conductors are feasible.

On the fundamental side, the measurement of an isotope effect with $`\alpha =0.26`$ where $`T_cM^\alpha `$ and $`M`$ is the isotopic mass of the boron atom indicated that phonons play a major role in the superconducting interaction in $`MgB_2`$. For clean samples, resistivity measurements at $`40K`$ indicated that $`\rho 1.0\mu ohmcm`$ which, combined with the calculated Fermi velocity,  $`v_F4.7\times 10^6m/s`$, gives a mean-free-path, $`\mathrm{}60nm`$. From the measured upper critical field, $`H_{c2}16T`$., the coherence distance of a clean sample can be determined from $`H_{c2}=\mathrm{\Phi }_o/[2\pi \xi _o^2]`$ to be $`\xi _o5nm`$. Here, $`\mathrm{\Phi }_o`$ is the flux quantum. With $`\mathrm{}`$ substantially greater than $`\xi _o`$, these materials are in the clean limit. Several groups found $`J_c`$ values in the range of $`500,000A/cm^2`$ in self field at $`5K`$, but these values fell below $`1,000A/cm^2`$ at $`1T`$ at $`5K`$. Broadly speaking, the irreversibility field, $`H_{irr}`$ was roughly half of $`H_{c2}`$.

The primary drawback for this material appeared to be the relatively low values of the upper critical field, $`H_{c2}=16T`$, and the resulting rapid drop in $`J_c`$ with field. The most straight forward way to raise $`H_{c2}`$ is to shorten the mean-free-path, $`\mathrm{}`$, which in turn shortens the coherence distance, $`\xi =\sqrt{\xi _o\mathrm{}}`$. Thus $`H_{c2}=\mathrm{\Phi }_o/[2\pi \xi ^2]`$ would rise with decreasing $`\mathrm{}`$. The introduction of point defects is an effective way to raise $`H_{c2}`$. In pure $`MgB_2`$, the volume pinning force seems to decrease as $`[1H/H_{irr}]^2`$ indicating that grain boundary pinning dominates. Hence, reducing the grain size should help. Alternatively, the strength of the pinning and $`J_c`$, can be increased by introducing precipitates. Typical dimensions might be an array of pinning sites with a characteristic dimension of the coherence distance, $`15nm`$, and a spacing between pinning centers about $`5`$ times larger.

Doping of $`MgB_2`$ has proven to be difficult because relatively few impurities go into the lattice in solid solution. Early reports by Slusky and co-workers show that $`Al`$ added by reacting powders in a solid state reaction substitute for $`Mg`$ in the $`MgB_2`$ lattice. There is a smooth change in lattice parameter and smooth decrease in $`T_c`$ for $`Mg_{1y}Al_yB_2`$ for $`y`$ out to $`0.10`$ and then a two phase region is entered. At $`y=0.1`$, $`T_c`$ is still about $`33K`$ and thus still relatively high.

For $`Ti`$ doping by solid state reaction of powders, it appears that relatively little $`Ti`$ is incorporated in the the $`MgB_2`$ grains. Work by Zhao and coworkersshowed that $`Ti`$ appears as very thin layers of $`TiB_2`$ precipitates on the $`MgB_2`$ grain boundaries. Samples were made via solid state reaction by grinding the $`Ti`$, $`Mg`$, and $`B`$ powders, pressing pellets, and heat treating. Microstructural analysis seemed to indicate that $`TiB_2`$ may act as a sintering agent leading to higher density, very fine $`MgB_2`$ grain-size of about $`10nm`$, and excellent grain-to-grain coupling. Some high resolution transmission electron microscopy evidence points to the thickness of some of these grain boundary layers being as low as one $`TiB_2`$ unit cell. Because the $`Ti`$ does not enter the unit cell, there is relatively little diminution in $`T_c`$ and very substantial increases in the critical current density, $`J_c`$. Optimum doping was in the range of a $`Mg/Ti`$ ratio of about $`0.1`$. One of the more important outcomes of this work was the observation that the $`Ti`$ seems to refine the grain structure to a grain size of about $`10nm`$.

For carbon doping, rather large amounts of carbon can replace boron in the structure. Takenobu and co-workers found that carbon substitutes for boron in $`MgB_{2x}C_x`$ out to $`x0.06`$ as indicated by a smooth decrease in a-axis lattice parameter from $`0.3085nm`$ at $`x=0`$ to $`0.3070nm`$ at $`x=0.06`$ while $`T_c`$ drops from $`38.2K`$ to about $`34.8K`$. There was relatively little change in the c-axis dimension. Here, samples were made by mixing powders of the elements and firing at $`900^{}C`$ for $`2h`$. At higher $`x`$-values, they found that the (110) X-ray peak of $`MgB_2`$ broadens indicating two phases. . Bharathi and co-workers, by contrast, found that reacting under high pressure $`Ar`$, permitted the extention of carbon content to $`x=0.5`$. They find a relatively smooth depression of $`T_c`$ from $`39K`$ at $`x=0`$ to about $`26K`$ at $`x=0.5`$. Using a somewhat different procedure, we have found that reacting $`B_4C`$ with $`Mg`$ vapor at $`950^{}C`$ for $`2h`$ gives a material with $`T_c17K`$. In other work, the change in lattice constant and $`T_c`$ with carbon content also has been reported to be somewhat less than these values.

For Si doping, a solid state reaction of $`Mg`$ and amorphous $`B`$ powder with nano-size $`Si`$ particles show substantial $`Mg_2Si`$ lines in the X-ray data and relatively little diminution in the transition temperature to about $`36K`$. This seems to imply that very little Si seems to enter the $`MgB_2`$ structure. In an alternate preparation method, Cooley and co-workers have shown $`Mg`$ can be reacted with $`SiB_6`$ at temperatures as low as $`500^{}C`$ to give $`Mg_2Si`$ precipitates in a $`MgB_2`$ matrix with $`T_c=37K`$. These preliminary experiments seem to indicate that the enhancements in flux pinning arise from $`Mg_2Si`$ precipitates rather than $`Si`$ substitution into the $`MgB_2`$ lattice.

The purpose of this review is to summarize progress in attempts to make $`MgB_2`$ conductors using a CVD method for depositing the starting boron material. There is a well established commercial scale method to make boron fiber in kilometer lengths for use in high strength composites. The primary goal of this work is to find ways to make doped boron fiber that will permit the preparation of long lengths of high performance $`MgB_2`$ conductor. Tasks include 1) the development of processing methods to convert pure boron fiber to $`MgB_2`$, 2) the development of methods to dope the starting boron fiber with point defects that will raise the upper critical field, $`H_{c2}`$, 3) the development of boron fibers that will produce nano-scale precipitates in the $`MgB_2`$, and 4) the development of methods to handle the brittle character of the superconductor.

Chemical vapor deposition of the starting boron is a particularly attractive route to materials preparation because the impurity atoms can possibly be introduced with atomic scale dispersion. Initial experiments have focussed on titanium-doping to produce small $`TiB_2`$ precipitate pinning centers with the next step to be the introduction of carbon point defects to shorten the electronic mean free path.

In a typical deposition process for boron fibers, a conducting filament such as tungsten or a carbon coated commercial silicon carbide fiber is introduced into a long glass chamber through a $`Hg`$ seal at a rate of a few $`cm/s`$. A flowing stream of $`BCl_3`$ and hydrogen gas move through the full length of the chamber. The fiber is heated electrically to temperatures in the $`11001300^{}C`$ range, and the boron is deposited in the hot zone. If a suitable carrier gas can be mixed with this gas stream, then impurity atoms such as $`Ti`$ or $`C`$ can be co-deposited with the boron to form filaments that can subsequently converted to impurity doped-$`MgB_2`$. Critical current densities were measured via magnetization as reported earlier.

## II Conversion of pure $`B`$-fiber to magnesium diboride

Soon after the discovery of superconductivity in $`MgB_2`$, small wire segments of commercial boron fiber were converted to $`MgB_2`$. These wire segments were only a few $`cm`$ long and rather fragile, but they showed exceptionally low resistivity of $`0.4\mu ohmcm`$ at $`40K`$ and a sharp superconducting transition at $`39.2K`$. In a subsequent series of upper critical field measurements, $`H_{c2}vsT`$ curves indicated an extrapolation to $`H_{c2}(T=0)`$ of about $`16T`$. Critical current densities, $`J_c`$ for a $`100\mu m`$ diameter pure boron fiber reacted in $`Mg`$ vapor at $`900^{}C`$ for $`2h`$ are shown for both $`5K`$ and $`25K`$ in Fig. 1. At $`5K`$, $`J_c5\times 10^5A/cm^2`$ in self field, and $`J_c`$ is about $`1000A/cm^2`$ at $`2.5T`$. At $`25K`$, $`J_c2\times 10^5A/cm^2`$ in self field, and $`J_c`$ is about $`1000A/cm^2`$ at $`1.4T`$. Also shown on Fig. 1 are data for a sample made by depositing $`25\mu m`$ of pure boron on a carbon substrate. The $`5K`$ data are the open circles and the $`25K`$ data are the solid triangles. Within the reproducibility of the experiments, the $`100\mu m`$ pure boron fibers give the same $`J_c`$ values as samples made from $`25\mu m`$ of pure B on a C substrate within the accuracy of the measurement. For all samples, the reaction is carried out at a single temperature, and the time can then be adjusted to optimize $`J_c`$ if that is desired.

A series of measurements were undertaken to determine the lowest temperatures where these boron fibers would transform to $`MgB_2`$ at a high rate. A sample reacted at $`800^{}C`$ for $`2h`$, as shown in Fig. 2, shows a $`Mg`$ region only part way across the boron fiber in the relatively light area near the edge of the filament. A more detailed analysis shows that this is a uniform wall of a high boron compound, $`MgB_x`$, progressing across the fiber. Energy dispersive spectra in the scanning electron microscope indicate that this phase has $`x`$ near $`6`$ or $`7`$. Because this phase is not the study of this work, further analysis was not undertaken. For this time and temperature, essentially no $`MgB_4`$ or $`MgB_2`$ forms. The white core of the filament is $`W_2B_5`$.

By contrast, if a similar pure boron fiber is heated in $`Mg`$ vapor at $`875^{}C`$ for $`2h`$, the fiber transforms to $`MgB_2`$ nearly everywhere as shown in Fig. 3. The darker gray areas are small regions of $`MgB_4`$ as shown by the EDS line scan along the line shown, where the circle with a cross is the starting point of the scan. The large dip in $`Mg`$ count rate is the $`W_2B_5`$ core and the $`30`$ percent dip in count rate to the left of the core illustrates the decrease in $`Mg`$ density going from $`MgB_2`$ to $`MgB_4`$. Also visible in the micrograph are several faults along the edge of the fiber and a large void near the right end of the portion of the $`W_2B_5`$ core. The filament is not completely flat, so that the $`W_2B_5`$ core does not show along the entire length of the filament for this particular depth of polish. These flaws are typical in $`MgB_2`$ filaments prepared from commercial pure $`B`$ fibers. In Fig. 4, an electron micrograph is shown in end view for the long time limit where a $`100\mu m`$ diameter fiber that has been reacted in $`Mg`$ vapor at $`950^{}C`$ for $`78h`$. The white core of $`W_2B_5`$ in the center is measured here to be $`16.7\mu m`$ in diameter and the starting $`100\mu m`$ boron fiber has grown to about $`135\mu m`$ diameter as $`Mg`$ is added. There is a small patch of $`MgB_4`$ in the upper left quadrant, and there the void pattern is very typical of all the samples grown.

## III Titanium doping

In the early experiments with $`CVD`$ doping of $`B`$ with $`Ti`$, the material was deposited on either a $`W`$ substrate or a $`C`$ coated substrate, and a preliminary report of this work has been presented. Most of the experiments reported here, however, were done by simultaneously depositing $`B`$ and a few percent $`Ti`$ on a continuous $`75\mu m`$ diameter silicon carbide filament coated with several micrometers of glassy carbon (commercially available as SCS-9A). Deposition was accomplished in a reactor similar to that used for commercial boron filament production. The substrate filament was drawn through the reactor at $`100mm/s`$; the filament in the reactor was resistively heated to peak temperatures of approximately $`1100^{}C`$. The reactant gas was predominantly a stoichiometric mixture of hydrogen and boron trichloride at atmospheric pressure. Titanium dopant was added by bubbling the hydrogen component through titanium tetrachloride held at $`0^{}C`$. The resulting partial pressure of titanium tetrachloride was approximately $`2\times 10^3bar`$. The doped boron-coated substrate was spooled continuously upon emerging from the reactor. An X-ray powder diffraction study of the resulting fibers before $`Mg`$ reaction indicates comparable amounts of $`TiB`$ and $`TiB_2`$ as well as a small amount of $`B_{10}C`$. Subsequent reaction with $`Mg`$ vapor produced $`MgB_2`$ filaments with precipited $`TiB_2`$. The intensity of the X-ray lines indicates less $`TiB`$ in the reacted fiber than in the unreacted fiber.

A series of samples were prepared in which $`B`$ and $`Ti`$ were co-deposited on the carbon surface with thicknesses of $`4\mu m`$, $`6\mu m`$ and $`10\mu m`$. After reaction in $`Mg`$ vapor at both $`950^{}C`$ for $`2h`$ and at $`1100^{}C`$ for $`30min`$, the conversion of $`B`$ to $`MgB_2`$ was complete. An extensive study of the $`Ti`$ distribution was undertaken with $`EDS`$ in a $`SEM`$ where the beam size was approximately $`1`$ micrometer. Generally, the $`Ti/Mg`$ ratio was in the neighborhood of $`9\%`$ in the bulk of the material. This is illustrated in Fig. 5 where a sample with $`6\mu m`$ of $`B9\%Ti`$ is deposited on a carbon coated substrate. There are two regions of the sample, however, that are consistently high $`Ti`$. This is near the $`C`$ substrate and near the outer surface. Reaction kinetics indicate that the deposition of $`Ti`$ compared to $`B`$ is higher at lower temperatures. Hence, as the fiber moves through the deposition zone, the $`Ti`$ deposition will be faster than that of boron during the initial warm-up on entrance to the hottest regions and during the cool-down as the fiber exits the hottest regions. Another specimen illustrated in Fig. 6 shows a $`4\mu m`$ thick layer of $`B9\%Ti`$. Here, there are two rather distinct regions of $`Ti`$ content. In most of the boron layer, the $`Ti`$ content is about $`9\%`$, but sometimes there are bright lacy regions where the $`Ti`$ content is about $`30\%`$. By studying polished regions of the $`MgB_2`$ with polarized light, the grain structure is readily apparent. If a particular region of the sample is studied both with polarized light and with $`SEM`$ with $`EDS`$ area maps, the $`Ti`$ is found to be spread throughout the grains. The grain boundaries are not seen to be high in $`Ti`$ on these area maps. This result differs from the beautiful decoration of $`MgB_2`$ grain boundaries with $`TiB_2`$ seen for samples made from powders in solid state reactions. Titanium doping using $`CVD`$ methods gives a dispersion of $`TiB_2`$ in $`MgB_2`$ that is different from that seen for solid state reaction of powders.

## IV Critical current densities

The addition of $`Ti`$ greatly enhances the $`J_c`$ values if the level of $`Ti`$ doping is properly adjusted. As reported previously, anything approaching a $`Ti/Mg`$ ratio of 1 is far too much $`Ti`$. The volume fraction of superconductor is reduced by half, and the value of $`T_c`$ drops to about $`20K`$. In accord with earlier findings, a $`Ti/Mg`$ ratio in the neighborhood of $`0.09`$ greatly enhances both $`J_c`$ and $`H_{irr}`$ while depressing $`T_c`$ only a few degrees. At $`5K`$ and self field, $`J_c`$ increases from about $`4\times 10^5A/cm^2`$ for pure $`MgB_2`$ to about $`2\times 10^6A/cm^2`$ for the $`9\%Ti`$ sample. With increasing magnetic field, $`J_c`$ crosses $`10,000A/cm^2`$ at $`5T`$ with $`9\%Ti`$ as compared to $`1.4T`$ for the pure $`MgB_2`$ sample as shown by the solid triangles in Fig. 7. At $`25K`$, the $`9\%Ti`$ sample crosses $`10,000A/cm^2`$ at about $`1.3T`$ compared with the pure $`MgB_2`$ sample which crosses $`10,000A/cm^2`$ at $`0.9T`$.

In addition to simply adding $`Ti`$ to the boron fibers, there is a great deal of optimization that can yet be done. When the starting $`B_{0.91}Ti_{0.09}`$ fibers are prepared, there is a large amount of both $`TiB`$ and $`TiB_2`$ in the boron. During the reaction to form $`MgB_2`$, some of the $`TiB`$ is converted to $`TiB_2`$, and there may be some migration of $`C`$ from the substrate into the sample. Much more needs to be learned about both of these processes. An optimization of $`J_c`$ at $`25K`$ with $`B_{0.91}Ti_{0.09}`$ shows that a reaction time of $`30min`$ is best for a reaction temperature of $`1100^{}C`$ and a reaction time of $`2h`$ is best for a reaction temperature of $`950^{}C`$. An optimization for $`J_c`$ values at $`5K`$, shown in Fig. 8, indicates that $`30min`$ at $`1100^{}C`$ gives the highest $`J_c`$ values.

## V Thickness dependence

For pure $`B`$ fibers deposited on a $`15\mu m`$ $`W`$ core, a series of $`B`$ fibers having diameters ranging from $`100\mu m`$ to $`300\mu m`$ were fully reacted in $`Mg`$ vapor at $`950^{}C`$ for times ranging from $`2h`$ to $`48h`$. In all cases, the reaction to form $`MgB_2`$ increases the fiber diameter by about $`35\%`$. At $`5K`$ and self field, the $`J_c`$ values range from $`700,000A/cm^2`$ to $`300,000A/cm^2`$ and have a magnetic field dependence similar to those shown in Fig. 1. At $`25K`$, the data of Fig. 9 indicate that $`J_c(H)`$ nearly is independent of diameter. To the accuracy and reproducibility of the measurement, all of these fiber diameters have the same $`J_cvsH`$ plots.

For the $`B9\%Ti`$ deposited on a carbon coated substrate, samples having a thickness of $`4\mu m`$, $`6\mu m`$, and $`10\mu m`$ have been studied. We expected a difference in $`J_c`$ among these samples caused by different amounts of $`C`$ impurity or caused by a different distribution of $`TiB`$ or $`TiB_2`$ during the $`Ti`$ and $`B`$ deposition on the $`C`$ substrate. These are relatively small effects as shown on Fig. 10. These data were taken for samples reacted at $`1100^{}C`$ for $`15min`$. A reaction time less than the $`30min`$ for optimum $`J_c`$ was intended to limit possible diffusion of $`C`$ from the substrate. In broad terms, the $`5K`$ $`J_c`$ data are about $`3\times 10^6A/cm^2`$ at low field and the $`J_c`$ data cross the $`10,000A/cm^2`$ line at about $`4T`$. There is not a large thickness dependence of $`J_c`$, but at high fields, the thicker samples seem to show a higher $`J_c`$ values. At low fields, the thinner samples seem to show the higher $`J_c`$ values. Similar data have been reported elsewhere. There is lots of room for improvement by optimizing the $`Ti`$ content, the reaction temperatures, and the reaction times.

## VI Conclusions

Co-deposition of boron along with the desired impurities offers a distinctive method of doping. When these doped fibers are converted to $`MgB_2`$ by heating in $`Mg`$ vapor, the resulting $`TiB_2`$ is distributed relatively uniformly within the grains. This distribution is quite different from $`MgB_2`$ samples prepared by solid state diffusion of mixed powders. Critical current densities of $`2`$ to $`5`$ million $`A/cm^2`$ are obtained at $`5K`$ and self field. Critical current densities of $`10,000A/cm^2`$ are obtained at $`25K`$ and $`1.4T`$.

A number of important variables are still not well understood. Improvement is needed in the control of the distribution of $`TiB_2`$, and a more thorough study of carbon doping needs to be undertaken so that a quantitative contrast can be developed between the effects of $`C`$ substitution in the lattice with the effects of $`TiB_2`$ precipitates. Although there is a fairly wide sweet spot for the reaction of doped $`B`$ fibers to form doped $`MgB_2`$, there is lots of room for improvement. The performance of the conductor can yet be optimized with respect to $`Ti`$ content and carbon content as well as the time and temperature of $`Mg`$ reaction.

## VII Acknowledgments

Work is supported by the U.S. Department of Energy, Basic Energy Sciences, Office of Science, through the Ames Laboratory under Contract No. W-7405-Eng-82.


format: text

text:
# On Breaking Time Reversal in a Simple, Smooth, Chaotic System

## Abstract

Within random matrix theory, the statistics of the eigensolutions depend fundamentally on the presence (or absence) of time reversal symmetry. Accepting the Bohigas-Giannoni-Schmit conjecture, this statement extends to quantum systems with chaotic classical analogs. For practical reasons, much of the supporting numerical studies of symmetry breaking have been done with billiards or maps, and little with simple, smooth systems. There are two main difficulties in attempting to break time reversal invariance in a continuous time system with a smooth potential. The first is avoiding false time reversal breaking. The second is locating a parameter regime in which the symmetry breaking is strong enough to transition the fluctuation properties fully toward the broken symmetry case, and yet remain weak enough so as not to regularize the dynamics sufficiently that the system is no longer chaotic. We give an example of a system of two coupled quartic oscillators whose energy level statistics closely match those of the Gaussian unitary ensemble, and which possesses only a minor proportion of regular motion in its phase space.

Since its introduction into nuclear physics by Wigner wigner , random matrix theory (RMT) has grown to encompass a broad variety of applications, and can be viewed as a significant portion of the foundation of the statistical mechanics of finite systems. One of the central tenets of RMT emphasized by Wigner was that the ensembles carry no information other than that required by the symmetries of the system. A very important symmetry is time reversal which determines whether the ensemble is composed of real symmetric matrices, Gaussian orthogonal ensemble (GOE - good symmetry), or complex hermitian matrices, Gaussian unitary ensemble (GUE - broken symmetry). This distinction leads to quite different predictions for the statistical properties of both the energy levels and the eigenfunctions. In fact, Wigner proposed it as a test of time reversal invariance in the strong interaction using slow neutron resonance data tri . His suggestion was not fully realized until more than twenty years later fkpt .

Time reversal is well known to be an antiunitary symmetry in quantum mechanics, and Robnik and Berry generalized the criterium for expecting the statistical properties of the GOE to include any antiunitary symmetry rb . An elementary example of an antiunitary symmetry would be the product of some reflection symmetry and time reversal. Some of the early investigations of noninvariant systems were of Aharanov-Bohm chaotic billiards br , symmetry breaking quantum maps izrailev ), and combinations of magnetic and scalar forces selig .

There is no known mechanical-type system of a particle moving under the influence of a simple, closed, smooth potential whose dynamics is rigorously proven to be fully chaotic, independently of the question of anti-unitary symmetry; we exclude diffusive dynamics associated with random potentials. For some period of time, the $`x^2y^2`$ quartic potential was a prime candidate, but eventually stable trajectories were found swedes . Nevertheless, there exists a family of quartic potentials with $`x^2y^2`$ as a limiting case whose fluctuations closely approximate GOE statistics, and whose classical dynamics contain negligible phase space zones of stable trajectories. Therefore,to construct a close approximation of a Hamiltonian without an anti-unitary symmetry and GUE statistics, we begin by considering the symmetry preseverving two-degree-of-freedom coupled quartic oscillators whose Schrödinger equation is given by

$$\widehat{H}_0\mathrm{\Psi }(\stackrel{}{r})=\frac{\mathrm{}^2^2}{2m}\mathrm{\Psi }(\stackrel{}{r})+\widehat{V}(\stackrel{}{r})\mathrm{\Psi }(\stackrel{}{r})=E\mathrm{\Psi }(\stackrel{}{r})$$
(1)
The potential can be expressed as

$$\widehat{V}(\stackrel{}{r})=a(\lambda )\left[\frac{x^4}{b}+by^4+2\lambda x^2y^2\right]$$
(2)
where $`a(\lambda )`$ is a convenient constant, $`b1`$ lowers the symmetry, and $`\lambda `$ gives the strength of the coupling; $`(b=1,\lambda =1)`$ is equivalent to $`V(\stackrel{}{r})=x^2y^2`$ by a $`\pi /4`$ rotation of the $`(x,y)`$ coordinates. The corresponding classical Hamiltonian is

$$H_0(\stackrel{}{r},\stackrel{}{p})=\frac{\stackrel{}{p}^2}{2m}+V(\stackrel{}{r})$$
(3)
This system is symmetric under time reversal and reflections with respect to $`x`$ and $`y`$. For strong couplings, $`\lambda 0.6`$, the statistics have been shown to agree extremely well with the GOE predictions btu . As a first attempt to break time reversal symmetry, we add the following term to the quantum Hamiltonian ($`\widehat{H}=\widehat{H}_0+\widehat{H}_1`$),

$$\widehat{H}_1=\frac{i\mathrm{}ϵ}{2m}\left(r\frac{}{r}r+\frac{}{r}r^2\right)$$
(4)
where $`r`$ is the radial polar variable. $`\widehat{H}_1`$ breaks time reversal invariance in the Hamiltonian without altering the original reflection symmetries, and thus does not admit an antiunitary symmetry from any combination of reflection and time reversal. This term was chosen to maintain a scaling property of the eigensolutions of the quartic oscillators due to the homogeneity of the potential btu . However, $`\widehat{H}_1`$ gives us an excellent example of false symmetry breaking. In fact, it turns out that $`\widehat{H}_1`$ is given by the cross terms arising from a vector potential $`𝐀(\stackrel{}{r})=ϵr^2\widehat{r}`$ (which can also be expressed as the gradient of a scalar function); i.e.

$$\widehat{H}_1=\frac{i\mathrm{}}{2m}\left(𝐀(\stackrel{}{r})+𝐀(\stackrel{}{r})\right)$$
(5)
This means that, up to the addition to the potential of a term $`ϵ^2r^4`$, $`H_1`$ can be understood as deriving from a magnetic field $`𝐁=\times 𝐀=0`$, which obviously will not change any physical quantity. And indeed, it is possible to make a gauge transformation and rewrite the wave function using the Dirac substitution

$$\mathrm{\Psi }(\stackrel{}{r})=\mathrm{exp}\left(\frac{i}{\mathrm{}}^\stackrel{}{r}𝐀(\stackrel{}{r}^{})d\stackrel{}{r}^{}\right)\mathrm{\Psi }^{}(\stackrel{}{r})=\mathrm{e}^{\frac{ir^3}{3\mathrm{}}}\mathrm{\Psi }^{}(\stackrel{}{r})$$
(6)
in order to cancel $`\widehat{H}_1`$ plus the aforementioned $`ϵ^2r^4`$ term. Furthermore, the transformation is single valued since $`𝐁=0`$ implies

$$𝐀(\stackrel{}{r}^{})d\stackrel{}{r}^{}=0.$$
(7)
The expectation, following the considerations of rb , is therefore to find GOE statistics and not those of the GUE. See the upper panel of Fig. (1), which compares the number variance of the system (i.e. the variance of the number of levels found in an energy interval of width $`s`$ scaled locally to mean unit level density) with the predictions of the GOE and GUE.

The GOE statistics are closely matched to a mean level spacing and a bit beyond. The parameter $`ϵ`$ was chosen slightly greater than unity because this forced the eigenlevels through a couple of avoided crossings, which would have been sufficient to push the spectral fluctuations toward GUE if the symmetry was not being falsely broken.

From a classical perspective, the Hamiltonian, not incorporating the jauge transformation, is

$$H(\stackrel{}{r},\stackrel{}{p})=\frac{p_r^22ϵr^2p_r}{2m}+\frac{p_\theta ^2}{2mr^2}+V(r,\theta )$$
(8)
which itself appears to violate time reversal symmetry, as well as in Hamilton’s equations of motion:

$`\dot{r}`$ $`=`$ $`{\displaystyle \frac{1}{m}}(p_rϵr^2)\dot{p}_r={\displaystyle \frac{V(r,\theta )}{r}}+2{\displaystyle \frac{ϵ}{m}}rp_r`$
$`\dot{\theta }`$ $`=`$ $`{\displaystyle \frac{p_\theta }{mr^2}}\dot{p}_\theta ={\displaystyle \frac{V(r,\theta )}{\theta }}`$
However even without considering a canonical transformation, just by converting to the Lagrangian description of the dynamics using the left hand side equations, it turns out that

$$=\frac{1}{2}(m\dot{r}^2+mr^2\dot{\theta }^2)V(r,\theta )+\frac{ϵ^2}{2m}r^4+ϵr^2\dot{r}$$
(10)
The final term, which appears to break the symmetry, cannot enter the equations of motion. They are invariant under addition of any total time derivative. It is not always obvious, a priori, whether a symmetry breaking term leads to false symmetry breaking or not.

If we multiply the symmetry breaking term by any function of $`\theta `$, it can no longer be a total time derivative, nor can the vector potential be expressed as the gradient of a scalar function. Consider

$$\widehat{H}_1^{}=\frac{i\mathrm{}ϵ}{2m}\mathrm{cos}^2\theta \left(r\frac{}{r}r+\frac{}{r}r^2\right)$$
(11)
The vector potential becomes $`𝐀(\stackrel{}{r})=ϵr^2\mathrm{cos}^2\theta \widehat{r}`$, and therefore $`𝐁=2ϵxy/r0`$. The integrals $`^\stackrel{}{r}𝐀(\stackrel{}{r}^{})d\stackrel{}{r}^{}`$ and $`𝐀(\stackrel{}{r}^{})d\stackrel{}{r}^{}`$ are path dependent. The Dirac substitution is not useful, and false symmetry breaking is not an issue.

Driving the eigenvalues at numerically attainable energies through a couple of avoided crossings forces $`ϵ`$ to be chosen in the neighborhood of unity. This has a strong regularizing effect on the nature of the dynamics. See the surface of section in Fig. (2) for the case ($`\lambda =0.8,ϵ=1`$).

Only its spectrum gives number variance statistics close to the GUE results; see the lower panel of Fig. (1).

To summarize, we have given an example of a simple, continuous, dynamical system that comes close to generating GUE statistics. It is surprisingly difficult to find an essentially, fully chaotic system that does so. The pitfall of false time reversal breaking can lead to symmetries that are quite well hidden, and the addition of a vector potential to a dynamical system has a strong tendency to move the system away from fully developed chaos.

We gratefully acknowledge support from ONR grant N00014-98-1-0079 and NSF grant PHY-0098027.


text:
# Causality and the Skyrme Model

## 1 Introduction

The Skryme model continues to yield insights into pion and nucleon physics . Initially, the static properties of the model were studied but time-dependent processes are also important. Their study led to the realization that during the time evolution the classical equations of motion, which are highly non-linear, may cease to be hyperbolic . This is potentially worrying because although the equations of motion are relativistically covariant, their non-linear character makes it unclear whether Einstein causality is always maintained during the time evolution. The danger is particularly great during violent processes such as Skyrmion-anti-Skyrmion annihilation . Presumably if Einstein causality were to fail one would have to give up the model, on the grounds that the various approximations made to obtain it from microphysics, for example from QCD, can no longer be valid, since those models incorporate Einstein causality as a basic assumption. In this connection it is perhaps interesting to recall the suggestion that superluminal motion of pions may be possible in hadronic fluids  when one loop effects at finite temperature are taken into account for the linear sigma model.

The point of this note is to show that, regardless of whether hyperbolicity breaks down, the energy momentum tensor $`T_{\mu \nu }`$ of the $`SU(2)`$ Skyrme model satisfies the dominant energy condition  and hence, by a result of Hawking , causal behaviour is guaranteed.

The dominant energy condition is a pointwise condition on $`T_{\mu \nu }`$ and states that $`T_\nu ^\mu t^\nu `$ is future directed timelike or null for every future directed timelike vector $`t^\mu `$. Equivalently $`T_{\mu \nu }t^\nu s^\mu 0`$, for all pairs of future directed timelike vectors $`t^\mu `$ and $`s^\mu `$. Another equivalent formulation is that in all local Lorentz frames $`T_{00}|T_{\mu \nu }|`$ for all index pairs $`\mu \nu `$. Yet another formulation is that if we can diagonalize $`T_{\mu \nu }`$ relative to the spacetime metric $`\eta _{\mu \nu }`$, so that $`\eta _{\mu \nu }=\mathrm{diag}(1,1,1,1)`$ and $`T_{\mu \nu }=\mathrm{diag}(T_{00},T_{11},T_{22},T_{33})`$, then $`T_{00}|T_{11}|`$ etc. The quantities $`T_{11}`$ etc are called ‘principle pressures ’. Note that one cannot always diagonalize an arbitrary symmetric tensor relative to $`\eta _{\mu \nu }`$ because it is indefinite, but it is possible in the generic case.

The intuitive idea behind this condition is that if it holds, then relative to any local frame, the flow of energy and momentum is future directed timelike or null. In other words, locally matter flows no faster than light. Using this condition Hawking was able to show, for example that if $`T_{\mu \nu }`$ vanishes at one time outside a compact set $`K`$, then $`T_{\mu \nu }`$ vanishes outside the future of $`K`$. In other words the local condition guarantees that bulk motion is never superluminal. In fact the main aim of Hawking’s paper was rather broader. He showed that classical field theories satisfying the dominant energy condition do not permit particle creation ex nihilo even if the background metric is curved and time-dependent. In fact classical models of particle creation in external fields always seem to entail some loss of causality, for instance particle and anti-particle appearing at spacelike separation . In the quantum theory of course, there is no reason to expect that the dominant energy condition continues to hold and creation of particles can take place. In the case of the Skyrme model one still expects topological conservation laws to hold however. In the context of quantum gravity, this may lead to some constraints on the number of created Skyrmions and whether it is odd or even .

A remark that will will prove useful in the sequel is that if $`T_{\mu \nu }^1`$ and $`T_{\mu \nu }^2`$ are two energy momentum tensors satisfying the dominant energy condition, then so does $`aT_{\mu \nu }^1+bT_{\mu \nu }^2`$, for positive $`a`$ and $`b`$. Put more formally, the set of energy momentum tensors satisfying the dominant energy condition form a Lorentz-invariant convex cone in the space of all second rank symmetric tensors.

The Skyrme model is based on a map $`\varphi ^A(t,𝐱)`$ from spacetime to a target model $`M`$ with positive definite metric $`G_{AB}(\varphi )`$, with $`A=1,2,\mathrm{},\mathrm{dim}G`$. In the simplest case $`M`$ is $`SU(2)S^3`$ with the bi-invariant or round metric. A more ambitious model takes $`M=SU(3)`$. The Lagrangian $`L`$ and the energy-momentum tensor $`T_{\mu \nu }`$ are constructed from the pulled back metric $`L_{\mu \nu }=G_{AB}(\varphi )_\mu \varphi ^A_\nu \varphi ^B`$. Generically we can diagonalize $`L_{\mu \nu }`$ relative to $`\eta _{\mu \nu }`$. The eigenvalues are necessarily non-negative so we write them as $`T_{\mu \nu }=\mathrm{diag}(\lambda _0^2,\lambda _1^2,\lambda _1^2,\lambda _3^2)`$. The eigen-values $`\lambda _1^2,\lambda _2^2,\lambda _3^2`$ were introduced by Manton in the static case for which $`\lambda _0^2=0`$ .

Consider, to begin with, a non-linear sigma model without Skyrme term. This has

$$L^{\mathrm{sigma}}=\frac{1}{2}L_\mu ^\mu ,T_{\mu \nu }^{\mathrm{sigma}}=L_{\mu \nu }+\frac{1}{2}L_\tau ^\tau \eta _{\mu \nu }.$$
(1)
In our special frame

$$T_{\mu \nu }^{\mathrm{sigma}}=\mathrm{diag}\frac{1}{2}(\lambda _0^2+\lambda _1^2+\lambda _2^2+\lambda _3^2,\lambda _0^2+\lambda _1^2\lambda _2^2\lambda _3^2,\lambda _0^2\lambda _1^2+\lambda _2^2\lambda _3^2,\lambda _0^2\lambda _1^2\lambda _2^2+\lambda _3^2).$$
(2)
By inspection, the sigma model energy momentum tensor satisfies the dominant energy condition. Although this guarantees causality, it certainly does not guarantee non-singular evolution. Indeed it is known that finite energy initial data can collapse to a singularity in finite time . One even has an explicit solution illustrating this collapse. If $`\chi ,\alpha ,\beta `$ are polar coordinates on $`S^3`$, and $`t,r,\theta ,\varphi `$ polar coordinates on spacetime, this solution is given, for $`t<0`$, by $`\chi =2\mathrm{tan}^1(r/t)`$ for $`rt`$ and $`\chi =\frac{\pi }{2}`$ for $`r>t`$. This represents a one parameter family of inverse stereographic projections.

In order to prevent collapse, one adds an additional term, $`L^{\mathrm{Skyrme}}`$ to the Lagrangian. For $`M=SU(2)`$, the Skyrme term may be written in suitable units as

$$L^{\mathrm{Skyrme}}=L_{\mu \nu }L^{\mu \nu }L_\mu ^\mu L_\nu ^\nu .$$
(3)
The contribution to the energy momentum tensor is

$$T_{\mu \nu }^{\mathrm{Skyrme}}=4\left(L_{\mu \nu }L^\sigma L_\sigma L_{\mu \lambda }L_\nu ^\lambda \right)+\eta _{\mu \nu }L^{\mathrm{Skyrme}}.$$
(4)
In the adapted frame

$$T_{00}^{\mathrm{Skryme}}=2\left(\lambda _1^2\lambda _2^2+\lambda _2^2\lambda _3^2+\lambda _3^2\lambda _1^2+\lambda _0^2\lambda _1^2+\lambda _0^2\lambda _2^2+\lambda _0^2\lambda _3^2\right),$$
(5)
$$T_{11}^{\mathrm{Skyrme}}=2\left(\lambda _1^2\lambda _2^2+\lambda _1^2\lambda _3^2\lambda _2^2\lambda _3^2\right)+2\lambda _0^2\left(\lambda _2^2+\lambda _3^2\lambda _1^2\right).$$
(6)
Clearly $`T_{00}^{\mathrm{Skyrme}}|T_{11}^{\mathrm{Skyrme}}|`$, which is sufficient to show that $`T_{\mu \nu }^{\mathrm{Skyrme}}`$ satisfies the dominant energy condition, and hence by the convexity property, so does the complete energy momentum tensor $`T_{\mu \nu }=T_{\mu \nu }^{\mathrm{sigma}}+T_{\mu \nu }^{\mathrm{Skyrme}}`$.

Since,

$$L_\mu ^\mu =T_{}^{\mathrm{sigma}}{}_{\mu }{}^{\mu },L_{\mu \nu }=T_{}^{\mathrm{sigma}}{}_{\mu \nu }{}^{}\frac{1}{2}\eta _{\mu \nu }T_{}^{\mathrm{sigma}}{}_{\lambda }{}^{\lambda },$$
(7)
one has the nice formula

$$T_{\mu \nu }^{\mathrm{Skyrme}}=4T_{\mu \lambda }^{\mathrm{sigma}}T_{}^{\mathrm{sigma}}{}_{\nu }{}^{\lambda }+\eta _{\mu \nu }T_{\lambda \rho }^{\mathrm{sigma}}T_{}^{\mathrm{sigma}}{}_{}{}^{\lambda \rho }.$$
(8)
However this does not immediately suggest a less basis dependent proof that the dominant energy condition holds. An interesting question is whether the breakdown of hyperbolicity is associated with the boundary of the cone. More practically, it might be illuminating to calculate the eigen-values $`\lambda _0^2,\lambda _1^2,\lambda _2^2,\lambda _3^2`$ during the numerical evolution as a diagnostic tool.

From (5, 6) one

$$T_{00}^{\mathrm{Skyrme}}+T_{11}^{\mathrm{Skyrme}}+T_{22}^{\mathrm{Skyrme}}+T_{33}^{\mathrm{Skyrme}}0.$$
(9)
One easily checks from (2) that (9) also holds for the sigma model energy momentum tensor. It follows that $`T_{\mu \nu }`$ belongs to the relativistically invariant convex cone of symmetric second rank tensors for which

$$\left(T_{\mu \nu }+\frac{1}{2}\eta _{\mu \nu }T_\lambda ^\lambda \right)t^\mu t^\nu 0,$$
(10)
for all future directed timelike or null vectors $`t^\mu `$. In other words, the strong energy condition  is also satisfied. Physically this means that the gravitational field of Skyrmion matter is always attractive.

A different extension of the non-linear sigma model Lagrangian to include higher powers of first derivatives is the Born-Infeld-Skyrme model. This has

$$L^{\mathrm{Born}\mathrm{Infeld}}=\sqrt{det\eta _{\mu \nu }}\sqrt{det(\eta _{\mu \nu }+L_{\mu \nu })}.$$
(11)
One finds that

$$T_{00}^{\mathrm{Born}\mathrm{Infeld}}=\frac{1}{\sqrt{1\lambda _0^2}}\left(\sqrt{(1+\lambda _1^2)(1+\lambda _2^2)(1+\lambda _3^2)}1\right),$$
(12)
and

$$T_{11}^{\mathrm{Born}\mathrm{Infeld}}=1\frac{\sqrt{1\lambda _0^2}}{1+\lambda _1^2}\sqrt{(1+\lambda _1^2)(1+\lambda _2^2)(1+\lambda _3^2)}.$$
(13)
Note that the arbitrary constant in the Lagrangian has been chosen to make it and the energy momentum tensor vanish for a constant field configuration. It follows from (12,13) that $`T_{00}^{\mathrm{Born}\mathrm{Infeld}}|T_{11}^{\mathrm{Born}\mathrm{Infeld}}|`$, and so the Born-Infeld-Skyrme energy momentum tensor also satisfies the dominant energy condition. However, unlike the usual Skyrme model, a scaling argument reveals there are no smooth static solutions with positive energy.

Finally we note that it follows easily from (12) and (13) that the Born-Infeld-Skryme energy momentum tensor also satisfies the strong energy condition. Thus neither Skyrme matter nor Born-Infeld-Skyrme matter is a suitable candidate for the dark energy currently expanding the universe. To violate the strong energy condition one may shift the arbitrary constant mutiple of $`g_{\mu \nu }`$ in $`T_{\mu \nu }`$. This looks less arbitrary if one also introduces a potential as one does in the case of tachyons and considers

$$L=V(\varphi )\sqrt{det(g_{\mu \nu }+L_{\mu \nu })}.$$
(14)
Now one finds that this can lead to a period of acceleration  but as it stands this is probably not a viable model either for the inflaton or for the dark energy.

## 2 Acknowledgements

Some of the ideas in this paper emerged during unpublished work with Michael Green. I would like to thank him and Nick Manton for many helpful and enjoyable discussions.


text:
# Arrival Distribution of Ultra-High Energy Cosmic Rays: Prospects for the Future

## 1 INTRODUCTION

The AGASA observation of ultra high energy cosmic rays (UHECRs) above $`10^{19}`$ eV reveals at least two features. The cosmic-ray energy spectrum does not show the GZK cutoff (Greisen, 1966; Zatsepin & Kuz’min, 1966) because of photopion production with the photons of the cosmic microwave background (CMB) and extends above $`10^{20}`$ eV (Takeda et al., 1998). On the other hand, their arrival distribution seems to be isotropic on a large scale with a statistically significant small scale clustering (Takeda et al., 1999). The current AGASA data set of 57 events above 4 $`\times 10^{19}`$ eV contains four doublets and one triplet within a separation angle of 2.5. Chance probability to observe such clusters under an isotropic distribution is only about 1 $`\%`$ (Hayashida et al., 2000).

Recently, the High Resolution Fly’s Eye (HiRes; Wilkinson et al., 1999) reports the cosmic ray flux with the GZK cut-off around $`10^{20}`$ eV (Abu-Zayyad et al., 2002). At present, it is very difficult to draw a conclusion about the existence or non-existence of the GZK cutoff, because both the two experiments have detected only a handful of events above this energy. On the other hand, there are new large-aperture detectors under development, such as South and North Auger project (Capelle et al., 1998), the EUSO (Benson & Linsley, 1982), and the OWL (Cline & Stecker, 2000) experiments. The detection or non-detection of the GZK cutoff in the cosmic-ray spectrum remains open to investigation by these future generation experiments.

Potential models of UHECR origin are constrained by their ability to reproduce the measured energy spectrum and the arrival distribution observed by the AGASA. In our recent work (Yoshiguchi et al., 2003a, hereafter Paper I), we perform numerical simulations for propagation of UHE protons in intergalactic space, and examine whether the present AGASA observation can be explained by a bottom-up scenario in which the source distribution of UHECRs is proportional to that of galaxies. We use the Optical Redshift Survey (ORS; Santiago et al., 1995) to construct realistic source models of UHECRs. We can construct realistic source models of UHECRs by using the galaxy sample, because the astrophysical candidates of UHECR sources, such as active galactic nuclei (AGN, Halzen & Zas, 1997), gamma-ray bursts (GRB, Waxman, 1995, 2000), and colliding galaxies (Cesarsky, 1992; Smialkowski, Giller, & Michalak, 2002) are in connection with the galaxies. For example, AGNs are considered as supermassive black holes in the centers of the galaxies. A GRB occurs in a galaxy, and so on.

In Paper I, we calculate both the energy spectrum and arrival directions of UHE protons, and compare the results with the AGASA observation. We find that the arrival distribution of UHECRs become to be most isotropic as restricting sources to luminous galaxies $`(M_{\mathrm{lim}}=20.5)`$. This is because luminous galaxies in the Local Super Cluster (LSC) distribute outward than faint galaxies, contrary to general clusters of galaxies (Yoshiguchi et al., 2003b). However, it is not isotropic enough to be consistent with the AGASA observation, even for $`M_{\mathrm{lim}}=20.5`$. In order to obtain sufficiently isotropic arrival distribution, we randomly select sources, which contribute to the observed cosmic ray flux, from the ORS sample more luminous than $`20.5`$ mag, and find that the isotropic arrival distribution of UHECRs can be reproduced in the case that the number fraction of $`1/50`$ of the sample is selected as UHECR sources. In terms of the source number density, this constraint corresponds to $`10^6`$ Mpc-3.

We further find that the small scale anisotropy can not be well reproduced in the case of strong extragalactic magnetic field ($`B10`$ nG). This is because the correlation at small scale between events which originate from a single source is eliminated, or the correlation continues to larger angle scale, due to large deflection when propagating in the EGMF from sources to the earth.

However, we should comment on the studies by other workers. Isola & Sigl (2002) and Sigl, Miniati, & Ensslin (2003) study the propagation of UHE protons in strong EGMF ($`1\mu `$ G) in the LSC, assuming local enhancement of UHECR sources in the LSC. The strong EGMF of $`1\mu `$ G leads to substantial deflections of UHECRs, which are better for explaining the observed isotropic distribution of UHECRs. However, the consistency of small-scale anisotropy and also large-scale isotropy predicted by their scenarios with the AGASA observation is marginal and somewhat worse than that predicted by our scenario in Paper I. There are also suggestions that the observed energy spectra have the imprints of UHE proton interaction with the CMB photons as beginning of the GZK cutoff at $`E4\times 10^{19}`$ eV (Berezinsky, Gazizov, & Grigorieva, 2002b, 2003; Marco, Blasi, & Olinto, 2003). In the presence of local enhancement of UHECR sources in the LSC, this feature of the observed spectra seems to be difficult to be reproduced. The AGASA observation may imply that UHECRs propagate along nearly straight lines in intergalactic space. However we can not draw any firm conclusion because of the limited amount of data. The next generation experiments will clarifies the situation with a large amount of data.

If local enhancement of UHECR sources in the LSC (Sigl, Lemoine, & Biermann, 1999; Lemoine, Sigl, & Biermann, 1999; Isola & Sigl, 2002; Sigl, 2002; Sigl, Miniati, & Ensslin, 2003) is disfavored from the observations, there is no way that explains the observed extension of the cosmic-ray spectrum beyond the GZK cutoff. Our conclusion in Paper I is that a large fraction of cosmic rays above $`10^{20}`$ eV observed by the AGASA experiment might originate in the top-down scenarios, or that the energy spectrum measured by the Hires experiment might be better.

As mentioned above, many future new experiments are under development. These experiments will provide us a large number of data, at least below $`10^{20}`$ eV, and allow us to discuss the features of the arrival distribution of UHECRs and determine which model of the UHECR origin can explain the observations with better statistical significance. In this paper, we predict the arrival distribution of UHECRs above $`4\times 10^{19}`$ eV with the event number expected by future experiments in the next few years. We perform event simulations with the source model which can explain the current AGASA observation. It is noted that our prediction is not the exact arrival directions of each UHECR but the statistical features of the arrival distribution, because there are degrees of freedom of randomly selecting the UHECR sources from the ORS sample. At first, we examine how much the future experiments decrease the statistical uncertainty of the cosmic-ray spectrum at the highest ($`10^{20}`$ eV) energies (Marco, Blasi, & Olinto, 2003). Next, we calculate the harmonic amplitude and the two point correlation function, and demonstrate that observational constraints on the model of UHECR origin become severer by new experiments. We further investigate the relation between the number of events clustered at a direction and the distance of their source. Such analysis has never performed before. Implications of the results are discussed in detail. Finally, we also predict the arrival distribution of UHECRs above $`10^{20}`$ eV, which is expected to be observed if the current HiRes spectrum is correct, and discuss their statistical features and implications.

In section 2, we describe our method of calculation. Results are shown in section 3. In section 4, we summarize the main results and discuss their implications.

## 2 NUMERICAL METHOD

### 2.1 Numerical Simulation

This subsection provides the method of Monte Carlo simulations for propagation of UHE protons in intergalactic space. We use the same numerical approach used in Paper I. Detailed explanations are presented in Paper I.

UHE protons below $`8\times 10^{19}`$ eV lose their energies mainly by pair creations and above it by photopion production (Berezinsky & Grigorieva, 1988; Yoshida & Teshima, 1993) in collisions with photons of the CMB. The pair production can be treated as a continuous loss process considering its small inelasticity ($`10^3`$). We adopt the analytical fit functions given by Chodorowski, Zdziarske, & Sikora (1992) to calculate the energy loss rate for the pair production on isotropic photons. On the other hand, protons lose a large fraction of their energy in the photopion production. For this reason, its treatment is very important. We use the interaction length and the energy distribution of final protons as a function of initial proton energy which is calculated by simulating the photopion production with the event generator SOPHIA (Mucke et al., 2000).

The EGMF are little known theoretically and observationally. There is the upper limit for the strength and correlation length of the universal EGMF, $`Bl_\mathrm{c}^{1/2}<1\mathrm{n}\mathrm{G}(1\mathrm{M}\mathrm{p}\mathrm{c})^{1/2}`$, as measured by Faraday rotation of radio signals from distant quasars (Kronberg, 1994). However, simple analytical arguments based on magnetic flux freezing, and large scale structure simulations passively including the magnetic field (Kulsrud et al., 1997) demonstrate that the magnetic field is most likely as structured as are the baryons. The local EGMF as strong as $`1\mu `$G in sheets and filaments of large scale galaxy distribution, such as in the LSC, are compatible with existing upper limits on Faraday rotation (Ryu, Hang, & Biermann, 1998; Blasi, Burles, & Olinto, 1999). It is suspected that the arrival distribution of UHECRs depends on the fields in the immediate environment of the observer.

However, we may expect that the effects of such strong EGMF on UHECR arrival directions may be small in our source scenario. In the source model by Isola & Sigl (2002) and Sigl, Miniati, & Ensslin (2003), all sources are assumed to be in the LSC. Thus, in their model the effect of the strong EGMF on the deflection angles of UHECRs is reported to be large. On the contrary, we show in Paper I that the source number density $`10^6`$ Mpc-3 is favored in order to explain the arrival distribution of UHECRs observed by the AGASA. In this case, there is no source in the LSC. Accordingly, almost all of the paths of the observed UHECRs from sources to the earth are not in the LSC but outside it in our model, because the structure of the strong EGMF in the LSC would be like two-dimensional sheet of large scale galaxy distribution, as shown in Sigl, Miniati, & Ensslin (2003). We also note that Faraday rotation measurement gives not the absolute value of the strength of the EGMF, but only its upper limit.

We further find that small scale clustering can not be well reproduced in the case of strong EGMF $`(B>10\mathrm{n}\mathrm{G})`$. If the local strong EGMF affects the arrival directions of UHECRs, small scale clustering observed by the AGASA may not be obtained. Thus, we assume that the local strong EGMF in the LSC, even if it exists with the structure like two-dimensional sheet, does not affect the arrival directions of UHECRs, and adopt a homogeneous random turbulent magnetic field with $`(B,l_\mathrm{c})=(1\mathrm{n}\mathrm{G},1\mathrm{M}\mathrm{p}\mathrm{c})`$.

Of course, we can not draw conclusion on the effects of the strong EGMF, considering the limited amount of observed data. However, since the model of local enhancement of UHECR sources (Isola & Sigl, 2002; Sigl, Miniati, & Ensslin, 2003) predicts the emergence of large-scale anisotropy which reflects the spatial structure of the LSC, we will be able to draw firm conclusion by comparing the prediction of our source scenario presented in this paper with the results of the future experiments. This is one of the purposes of the present work.

We also note that numerical simulations of UHECR propagation in inhomogeneous EGMF over cosmological distances are highly time-consuming. (On the other hand, Isola & Sigl (2002) and Sigl, Miniati, & Ensslin (2003) perform numerical simulations of UHECR propagation over $``$ 50 Mpc.) With the assumption of homogeneous EGMF, we can perform numerical simulation of UHECR propagation under spherical symmetry. Provided that the EGMF is inhomogeneous, however, the propagation of UHECRs from a single point source has no longer spherical symmetry. As a result, we have to specify the earth position in the universe. We also have to choose the detector (earth) size so small enough for us to accurately calculate the arrival directions. In this case, the number fraction of UHECRs arriving at the earth to injected ones is extremely small. This requires the number of particle to be propagated several orders of magnitudes higher than that used in this study, which takes enormous CPU time. Detailed study on the effects of the strong EGMF in the LSC are beyond the scope of this paper and left for future study.

We assume turbulent magnetic field with power-law spectrum $`B^2(k)k^{n_H}`$ for $`2\pi /l_\mathrm{c}k2\pi /l_{\text{cut}}`$ and $`B^2(k)=0`$ otherwise, where $`l_{\text{cut}}`$ characterizes the numerical cut-off scale. We use $`n_H=11/3`$ corresponding to the Kolmogorov spectrum. Physically one expect $`l_{\text{cut}}l_\mathrm{c}`$, but we set $`l_{\text{cut}}=1/8\times l_\mathrm{c}`$ in order to save the CPU time. The universe is covered with cubes of side $`l_\mathrm{c}`$. For each of the cubes, Fourier components of the EGMF are dialed on a cubic cell in wave number space, whose side is $`2\pi /l_\mathrm{c}`$, with random phases according to the Kolmogorov spectrum, and then Fourier transformed onto the corresponding cubic cell in real space. We create the EGMF of 20 $`\times `$ 20 $`\times `$ 20 cubes of side $`l_\mathrm{c}`$, and outside it, adopt the periodic boundary condition in order to reduce storage data for magnetic field components. Similar methods for the turbulent magnetic fields have been adopted (Sigl, Lemoine, & Biermann, 1999; Lemoine, Sigl, & Biermann, 1999; Isola & Sigl, 2002). In this study, we neglect the effects of the galactic magnetic field. We will conduct studies on its effects in forthcoming paper.

Finally, we explain how the energy spectrum and the arrival directions of UHECRs are calculated. At first, protons with a flat energy spectrum are injected isotropically at a given point within the range of ($`10^{19.5}`$ \- $`10^{22}`$)eV. 5000 protons are injected in each of 26 energy bins, that is, 10 bins per decade of energy. Then, UHE protons are propagated in the EGMF over $`1`$ Gpc for $`15`$ Gyr. Weighted with a factor corresponding to a $`E^2`$ power law spectrum, this provides distribution of energy, deflection angle, and time delay of UHECRs as a function of the distance from the initial point. In this paper, we use the distribution of energies and deflection angles integrated over the time delay, assuming that the cosmic ray flux at the earth is stationary. With this distribution, we can calculate the energy spectrum and the arrival directions of UHECRs injected at a single UHECR source. Then, summing contributions from all the sources (see the section 2.2), we obtain the angular probability distributions of UHECRs as a function of their energies. According to this angular probability distributions, we simulate the energy spectrum and the arrival directions of UHECRs.

### 2.2 Source Distribution

In this study, we assume that the source distribution of UHECRs is proportional to that of the galaxies. We use the realistic data from the ORS (Santiago et al., 1995) galaxy catalog. As mentioned in section 1, we show in Paper I that the arrival distribution of UHECRs observed by the AGASA can be reproduced in the case that the number fraction of $`1/50`$ of the ORS galaxies more luminous than $`M_{\mathrm{lim}}=20.5`$ is selected as UHECR sources. We consider only the prediction of this source model throughout the paper. It is unknown how much an ultimate UHECR source contribute to the observed cosmic ray flux. In paper I, we thus consider the two cases in which all galaxies are the same, and they inject cosmic rays proportional to their absolute luminosity. However, we find that the results in the two cases do not differ from each other, as far as we focus on the luminous galaxies as UHECR sources. Accordingly, we restrict ourselves to the case that all galaxies inject cosmic rays with the same amount.

In order to calculate the energy spectrum and the distribution of arrival directions of UHECRs realistically, there are two key elements of the galaxy sample to be corrected. First, galaxies in a given magnitude-limited sample are biased tracers of matter distribution because of the flux limit. Although the sample of galaxies more luminous than $`20.5`$ mag is complete within 80 $`h^1`$ Mpc, it does not contain galaxies outside it for the reason of the selection effect, where $`h`$ is the Hubble constant divided by 100 km s-1 and we use $`h=0.75`$. We distribute sources of UHECRs outside 80 $`h^1`$ Mpc homogeneously, and calculate their amount from the number of galaxies inside it. Second, our ORS sample does not include galaxies in the zone of avoidance ($`|b|<20^{}`$). In the same way, we distribute UHECR sources in this region homogeneously, and calculate its number density from the number of galaxies in the observed region.

### 2.3 Statistical Methods

In this subsection, we explain the three statistical quantities, the harmonics analysis for large scale anisotropy (Hayashida et al., 1999), the two point correlation function for small scale anisotropy, and the correlation value for investigation of the correlation between the events and their sources defined in our previous study (Ide et al., 2001).

The harmonic analysis to the right ascension distribution of events is the conventional method to search for global anisotropy of cosmic ray arrival distribution. For a ground-based detector like the AGASA and the Auger, the almost uniform observation in right ascension is expected. The $`m`$-th harmonic amplitude $`r`$ is determined by fitting the distribution to a sine wave with period $`2\pi /m`$. For a sample of $`n`$ measurements of phase, $`\varphi _1`$, $`\varphi _2`$, $`\mathrm{}`$, $`\varphi _n`$ (0 $`\varphi _i2\pi `$), it is expressed as

$$r=(a^2+b^2)^{1/2}$$
(1)
where, $`a=\frac{2}{n}\mathrm{\Sigma }_{i=1}^n\mathrm{cos}m\varphi _i`$, $`b=\frac{2}{n}\mathrm{\Sigma }_{i=1}^n\mathrm{sin}m\varphi _i`$. We calculate the harmonic amplitude for $`m=14`$ from a set of events generated according to predicted probability density distribution of arrival directions of UHECRs.

If events with total number $`n`$ are uniformly distributed in right ascension, the chance probability of observing the amplitude $`r`$ is given by,

$$P=\mathrm{exp}(k),$$
(2)
where

$$k=nr^2/4.$$
(3)
The current AGASA 57 events is consistent with isotropic source distribution within 90 $`\%`$ confidence level (Takeda et al., 1999; Hayashida et al., 2000). We therefore compare the harmonic amplitude for $`P=0.1`$ with the model prediction, and estimate the event number at which large scale anisotropy of the UHECR arrival distribution become significant.

The two point correlation function $`N(\theta )`$ contains information on the small scale anisotropy. We start from a set of generated events or the actual AGASA data. For each event, we divide the sphere into concentric bins of angular size $`\mathrm{\Delta }\theta `$, and count the number of events falling into each bin. We then divide it by the solid angle of the corresponding bin, that is,

$`N(\theta )={\displaystyle \frac{1}{2\pi |\mathrm{cos}\theta \mathrm{cos}(\theta +\mathrm{\Delta }\theta )|}}{\displaystyle \underset{\theta \varphi \theta +\mathrm{\Delta }\theta }{}}1[\mathrm{sr}^1],`$ (4)
where $`\varphi `$ denotes the separation angle of the two events. $`\mathrm{\Delta }\theta `$ is taken to be $`1^{}`$ in this analysis. The AGASA data shows strong correlation at small angle $`(2^{})`$ with 5 $`\sigma `$ significance of deviation from an isotropic distribution (Takeda et al., 1999; Hayashida et al., 2000).

We use the correlation value, defined in our previous study (Ide et al., 2001), in order to investigate statistically the similarity between the arrival distribution of UHECRs and the source distribution. The correlation value, $`\mathrm{\Xi }`$, between two distributions $`f_\mathrm{e}`$ and $`f_\mathrm{s}`$, is defined as

$$\mathrm{\Xi }(\mathrm{f}_\mathrm{e},\mathrm{f}_\mathrm{s})\frac{\rho (\mathrm{f}_\mathrm{e},\mathrm{f}_\mathrm{s})}{\sqrt{\rho (\mathrm{f}_\mathrm{e},\mathrm{f}_\mathrm{e})\rho (\mathrm{f}_\mathrm{s},\mathrm{f}_\mathrm{s})}},$$
(5)
where

$$\rho (f_a,f_b)\underset{j,k}{}\left(\frac{f_a(j,k)\overline{f}_a}{\overline{f}_a}\right)\left(\frac{f_b(j,k)\overline{f}_b}{\overline{f}_b}\right)\frac{\mathrm{\Delta }\mathrm{\Omega }(j,k)}{4\pi }.$$
(6)
Here subscripts $`j`$ and $`k`$ discriminate each cell of the sky, $`\mathrm{\Delta }\mathrm{\Omega }(j,k)`$ denotes the solid angle of the ($`j,k`$) cell, and $`\overline{f}`$ means the average of $`f`$. In equation 5, $`f_\mathrm{e}`$ and $`f_\mathrm{s}`$ represent the distribution of the simulated events and the sources, respectively. In this study, the size of the cell is chosen to be $`1^{}\times 1^{}`$. The meaning of $`\mathrm{\Xi }`$ is as follows. By definition, $`\mathrm{\Xi }`$ ranges from $`1`$ to $`+1`$. When $`\mathrm{\Xi }=`$ +1(-1), two distributions are exactly same (opposite). When $`\mathrm{\Xi }=`$ 0, we can not find any resemblance between two distributions.

## 3 RESULTS

### 3.1 Statistical Significance of the Energy Spectrum at $`10^{20}`$ eV

Before we discuss the future prospect of the UHECR arrival distribution, we examine how much the future experiments decrease the statistical uncertainty of the cosmic-ray spectrum at the highest energy range $`(10^{20}\mathrm{eV})`$. In figure 3.1, we show the energy spectrum predicted by a specific source scenario in the case that the number fraction of $`1/50`$ of the ORS galaxies more luminous than $`M_{\mathrm{lim}}=20.5`$ is selected as UHECR sources. The injection spectrum is taken to be $`E^2`$.

Berezinsky, Gazizov, & Grigorieva (2002a) show that predicted flux of UHECRs fall short of the observed flux below $`10^{19.5}`$ eV in the case of injection spectrum $`E^2`$. However, there may be UHECR production sites in which maximum energy of cosmic ray achieved is lower than $`10^{19.5}`$ eV. These components may substantially contribute to the cosmic ray flux below $`10^{19.5}`$ eV. Throughout the paper, we assume that these components bridge the gap between the observed flux and the predicted one with the injection spectrum $`E^2`$, and restrict ourselves to cosmic rays only above $`10^{19.5}`$ eV.

In figure 3.1, event simulations are performed with the fixed event numbers above $`4\times 10^{19}`$ eV. Typically, we perform 10000 such simulations. The shaded region indicates the statistical error due to the finite number of simulated events. The spectrum measured by the HiRes is consistent with our model prediction, while that of the AGASA is not. However, the statistical significance of deviation from the prediction of our source scenario is about only $`2\sigma `$. The region of the energy spectrum dominated by statistical fluctuation is moved to higher energies with increasing the data, as shown in figure 3.1. It is noted that the future experiment, such as the Auger and the EUSO/OWL, would detect $`500`$ and $`5000`$ events above $`4\times 10^{19}`$ eV per year, respectively. This high statistics will allow us to conclude the presence or absence of the GZK cutoff in the cosmic-ray spectrum in the next few years. Similar conclusion is obtained in Marco, Blasi, & Olinto (2003).

### 3.2 Arrival Distribution of UHECRs above 4$`\times 10^{19}`$eV

In this subsection, we present the results of the numerical calculations for the arrival distributions of UHECRs. Figure 1 shows realizations of the UHECR arrival direction above $`4\times 10^{19}`$ eV predicted by a specific source scenario in the case that the number fraction of $`1/50`$ of the ORS galaxies more luminous than $`M_{\mathrm{lim}}=20.5`$ is randomly selected as UHECR sources. Distribution of selected sources within $`200`$ Mpc is also shown as circles of radius inversely proportional to their distances. Only the sources within $`100`$ Mpc are shown with bold lines. Throughout the paper, we show the results only for this specific source scenario. We have checked that the results do not depend very much on the realization of the source selection, unless extremely nearby sources ($`<30\mathrm{M}\mathrm{p}\mathrm{c}`$) are accidentally selected. We note that the current AGASA data set include $`49`$ events with energies of $`4\times 10^{19}10^{20}`$ eV in the range of $`10^{}\delta 80^{}`$. This event number corresponds to $`100`$ events in figure 1, where we do not restrict the arrival directions of UHECRs in any range of $`\delta `$.

A visual inspection of figure 1 reveals no significant large scale anisotropy. We show the harmonic amplitude as a function of the event number for m$`=14`$ in figure 2. We plot the average over all trial of the event realizations from the calculated probability distribution with the statistical error. In order to obtain the average and the variance, we dial the simulated sets of events $`201000`$ times depending on the total event number. The region below the solid line in this figure is expected from the statistical fluctuation of isotropic source distribution with the chance probability larger than 10$`\%`$. Of course, this region become smaller with increasing the event number.

At present, our source model predicts the harmonic amplitude consistent with the isotropic source distribution. However, future experiments will separate our model prediction from the isotropic source at the confidence level of $`90\%`$ with the event number of the order $`10^3`$. It is also found that the amplitude for $`m=1`$ is smaller than that for another values of $`m`$. The amplitude for larger $`m`$ quantifies anisotropy on smaller scale. Therefore, this dependence on $`m`$ reflects the fact that the arrival distribution shown in figure 1 becomes to reveal the distribution of their sources with increasing the event number, in the manner that the event clusterings occur at the directions of the nearby ($`<100\mathrm{Mpc}`$) sources. This feature of the UHECR arrival distribution is discussed below in detail.

In this source scenario, the small scale anisotropy observed by the AGASA is also reproduced, as is evident from figure 3, where the two point correlation functions of the simulated events are shown. The errorbars represent the statistical fluctuations due to the finite number of the simulated events. For the event number of $`100`$, we also show the two point correlation function on the AGASA 49 events in the energy range of $`4\times 10^{19}10^{20}`$ eV, multiplied by factor 2 in order to compensate the difference of the range of $`\delta `$ between the observation and the numerical calculation. Slight shrinkage of $`N(\theta )`$ at the smallest angle bin is due to manner of dividing the sphere into concentric bins when taking the data of numerical simulations of UHECR propagation.

Since the small scale anisotropy is due to the point-like nature of UHECR sources, the angle scale at which there is strong correlation between the events corresponds to the deflection angle of UHECRs in propagation in the EGMF from sources to the earth. Figure 3 demonstrates that, for the amount of data expected with next generation experiments, the statistical uncertainty will considerably decrease. This clarify the angle scale at which the events have strong correlation with each other. Accordingly, we will be able to determine the strength of the universal EGMF from the two point correlation function of observed UHECRs with the sufficient amount of data.

As mentioned above, our source model predicts the statistically significant small scale anisotropy, which correlate with the sources located within $`100`$ Mpc, with a large number of data (see, figure 1). What is deduced from this feature on the origin of UHECRs? To begin with, we quantitatively examine the relation between the event distribution and the source distribution. The upper panel of figure 3.2 shows the correlation value defined in section 2.3 between the distribution of events above $`4\times 10^{19}`$ eV and the distribution of sources within 100, 200, and 300 Mpc as a function of the event number. The results for the events above $`10^{20}`$ eV are discussed in the next subsection. The errorbars are same ones as that in Figure 2.

Clearly visible in this figure is that the correlation of the event distribution with the sources are strongest for the sources within $`100`$ Mpc. This strong correlation is due to the event sets which cluster in the direction of the sources within 100 Mpc (See Figure 1). The number of clustered events fluctuates every realizations, and this number is a critical factor for the correlation with the sources within 100 Mpc. On the other hand, this number do not affect the correlation with the sources within larger distances very much, because there are a number of sources in this case. Therefore, the statistical error is smaller for the correlation with the sources at larger distances. When the event number is close to the order of several $`\times 10^3`$, the correlation values begin to converge, and the final values can be estimated. We emphasize that the expected event rate by the Auger observation is $`500`$ per year. After several years from the operation of the Auger, we would be able to know the source distribution within $``$ 100 Mpc.

Here we note that the UHECR sources outside $`107(=80h^1)`$ Mpc in our model are randomly distributed, as mentioned in section 2.2. We should compare the results of the correlation function between the simulated events and the actual galaxy distribution with that between the observed events and the actual galaxy distribution. For the reason of the flux limit, we can not do so using the ORS galaxy sample. However, there is the galaxy survey with a limiting magnitude much deeper than that of the ORS, Sloan Digital Sky Survey (SDSS, Stoughton et al., 2002). We conducted a study of the galaxy number density based on a comparison of the observed number counts between the ORS and the SDSS Early Data Release (Yoshiguchi et al., 2003b). As data obtained by the SDSS accumulate, we will be able to know the actual galaxy distribution of much larger volume, and make precise comparisons between numerical calculations and the observations.

In order to further investigate the relation between the number of the clustered events and the source distance, we calculate the number of events above $`4\times 10^{19}`$ eV within $`2.5^{}`$ from the direction of each source of Figure 1. The angle $`2.5^{}`$ roughly corresponds to both the observational error of arrival directions and the deflection angle of UHECR when propagating in the EGMF ($`B=1\mathrm{n}\mathrm{G}`$) over $`100`$ Mpc. The result is shown in Figure 4. The solid angle of the observer viewed from distant sources are inversely proportional to the square of the distance. Thus, contribution to cosmic-ray flux from a source closer to our galaxy is larger than that from a distant source. This is reflected in Figure 4, where the event number at the direction of nearby sources are larger than that of distant sources. A source at $``$ 280 Mpc which have larger number of event in its direction happens to be located at the direction of a closer source.

We again note that the event number $`100`$ corresponds to the observed one by the AGASA experiment. From figure 4, there must be a source at the direction of triplet event sets within 100 Mpc from us at $`2\sigma `$ confidence level. This implies that the triplet observed by the AGASA would originate from sources within $`100`$ Mpc. Indeed, Smialkowski, Giller, & Michalak (2002) show that there are merger galaxies Arp 299 (NGC 3690 $`+`$ IC 694) at the direction of the AGASA triplet at $`42`$ Mpc. Considering the analysis presented here, Arp 299 is the best candidate of the UHECR source.

Increasing the event number decrease the statistical uncertainty as is evident from figure 4, and thus the relation between the source distance and the number of clustered events in its direction becomes clear. As data accumulate by future experiment, we can know the distance of the source which contributes to a clustered event set by using this relation. Performing this procedure for all the event clusterings, we would be able to determine the distribution of UHECR sources within about 100 Mpc.

### 3.3 Arrival Distribution of UHECRs above $`10^{20}`$ eV

In this subsection, we present the results of the arrival distribution of UHECRs above $`10^{20}`$ eV. It is noted that our source model predicts the cosmic-ray spectrum with the GZK cutoff, unless other components are introduced at this energy range. Accordingly, the features of arrival distributions that we present here are expected to be observed by future experiments if the current HiRes spectrum is correct. If the AGASA spectrum is correct, UHECRs of top-down origin may dominate the cosmic-ray flux at this energy range. However, we may be able to extract UHECRs of bottom-up origin from ones of top-down origin by using the feature of arrival distribution of bottom-up UHECRs, as discussed in the final section.

In figure 5, we show the arrival distribution of UHECRs above $`10^{20}`$ eV predicted by the source model of figure 1. This figure is same as figure 1, but only the events with energies above $`10^{20}`$ eV are shown. Compared to figure 1, we find that the arrival directions are concentrated in the few directions. This can be also seen from figure 6, which is the same figure as figure 3, but for only UHECRs above $`10^{20}`$ eV. At $`\theta >4^{}`$, the values of $`N(\theta )`$ are almost equal to $`0`$. This is because the sources of UHECRs above $`10^{20}`$ eV must be located in a limited volume of radius at most $`100`$ Mpc because of the pion production. Sources at larger distances can not contribute to the cosmic-ray spectrum above this energy.

This feature of the correlation between the events and the sources can be also understood from the lower panel of figure 3.2. Note that the horizontal axis of this figure is the event number above $`4\times 10^{19}`$ eV. Since the number of events above $`10^{20}`$ eV is very small, the statistical error is larger for $`E_{\mathrm{th}}=10^{20}`$ eV than $`E_{\mathrm{th}}=4\times 10^{19}`$ eV. From this figure, it is clear that the correlation with the sources within $`100`$ Mpc is stronger for $`E_{\mathrm{th}}=10^{20}`$ eV than that for $`E_{\mathrm{th}}=4\times 10^{19}`$ eV. On the other hand, this dependence becomes opposite in the case of the sources within $`200`$ and $`300`$ Mpc, because UHECRs below $`10^{20}`$ eV can come from a much larger volume ($`1`$ Gpc3).

We also calculated the harmonic amplitude and the number of events within $`2.5^{}`$ from the directions of the sources for UHECRs above $`10^{20}`$ eV. The volume from which UHECRs at this energies can originate is much smaller than that for $`E<10^{20}`$ eV. Although the dependence of the harmonic amplitude on random selection of UHECR sources from the ORS sample is relatively large for this reason, significant anisotropy on large angle scale may be observed when future experiments detect about $`\mathrm{several}\times 10`$ cosmic rays above $`10^{20}`$ eV. The dependence of the number of events above $`10^{20}`$ eV within $`2.5^{}`$ from the directions of the sources on their distances is almost same as that for above $`4\times 10^{19}`$ eV.

## 4 SUMMARY AND DISCUSSION

In this paper, we predicted the arrival distribution of UHECRs above $`4\times 10^{19}`$ eV with the total number of events expected by next generation experiments in the next few years. We performed event simulations using the ORS galaxy sample to construct a source model of UHECRs, which can explain the current AGASA observation below $`10^{20}`$ eV (Paper I). It is noted that our prediction is not the exact arrival directions of each UHECR but the statistical features of the arrival distribution, because there are degrees of freedom of randomly selecting the UHECR sources from the ORS sample and this sample does not contain galaxies outside $`107(=80h^1)`$ Mpc. However, we will be able to know the actual galaxy distribution of much larger volume by using the SDSS galaxy sample (Stoughton et al., 2002) with increasing amount of data.

At first, we calculated the harmonic amplitude and the two point correlation function for the simulated event sets. We found that significant anisotropy on large angle scale will be observed when $`10^3`$ cosmic rays above $`4\times 10^{19}`$ eV are detected by future experiments. The Auger array will detect cosmic rays with this event number in a few years after its operation. The local enhancement model of UHECR sources (Isola & Sigl, 2002; Sigl, Miniati, & Ensslin, 2003) predicts large-scale anisotropy of UHECR arrival distribution that reflects the spatial structure of the LSC, which can not be seen from our model prediction (see figure 1). Thus we will be able to determine which model is favored from the observations in a few years. The statistics of the two point correlation function will also increase, and the angle scale at which there is strong correlation between the events corresponds to deflection angle of UHECR in propagating in the EGMF. Thus, it is expected that we will be able to know the strength of the EGMF, using the two point correlation function for observed arrival distribution of UHECRs, if our source and magnetic field model is supported by the future experiments.

Next, we investigated the relation between the number of clustered events and the distance of the source at their direction. We found that the C2 triplet observed by the AGASA (Hayashida et al., 2000) may originate from the source within 100 Mpc from us. Indeed, Smialkowski, Giller, & Michalak (2002) show that there are merger galaxies Arp 299 (NGC 3690 $`+`$ IC 694) at the direction of the AGASA triplet at $`42`$ Mpc. Considering the analysis presented here, Arp 299 is the best candidate of the UHECR source.

When the event number increases, the statistical uncertainty decreases as is evident from figure 4, and thus the relation between the source distance and the number of clustered events in its direction becomes clear. Using this relation, we will be able to determine the distribution of UHECR sources within about 100 Mpc.

Identification of the sources of UHECRs is extremely important. At first, this will provide some kinds of information about poorly known parameters which influence the propagation of UHECRs, such as extragalactic and galactic magnetic field, chemical composition of observed cosmic rays. Furthermore, this will give an invaluable information on mechanisms and physical conditions which lead to acceleration of cosmic rays to energies of order $`10^{20}`$ eV. In particular, we showed that there was the strong correlation between the arrival distribution of UHECR above $`10^{20}`$ eV and the source distribution within $`100`$ Mpc. If the cosmic-ray spectrum measured by the HiRes experiment is correct, the UHECR arrival distribution similar to figure 5 will be observed by future experiments. We will be able to know the maximum energies achieved by cosmic rays in each identified object. If the AGASA spectrum is correct, and cosmic-ray flux is dominated by the component of top-down origin. However, the number density of the supermassive particles, whose decay product can be observed UHECRs above $`10^{20}`$ eV, is estimated as $`10^{36}`$ Mpc-3 in our galactic halo in order to explain the observed flux, as we discussed in Paper I. In this case, there would be no small scale anisotropy of arrival distribution of UHECRs above $`10^{20}`$ eV which are generated by the top-down mechanisms. On the other hand, the arrival directions of bottom-up UHECRs are strongly concentrated in the few directions, as presented in the previous section. We may be able to extract UHECRs of bottom-up origin from ones of top-down origin, and obtain an information on the maximum energies of cosmic rays.

In Paper I, we showed that the number density of UHECR sources may be $`10^6`$ Mpc-3 in order to explain the UHECR arrival distribution observed by the AGASA experiments. Nevertheless, we could not know which of the astrophysical objects mainly contributes to the observed cosmic-ray flux. However, with the method to identify the sources of UHECR developed in this paper, we will reveal their origin and obtain an useful information on acceleration mechanism to the highest energy.

This research was supported in part by Giants-in-Aid for Scientific Research provided by the Ministry of Education, Science and Culture of Japan through Research Grant No.S14102004 and No.14079202.


text:
# Absolute Drell-Yan Dimuon Cross Sections in 800 GeV/𝑐 𝑝𝑝 and 𝑝𝑑 Collisions

Present address: \]Indiana University, Bloomington, IN 47405

Present address: \]University of Illinois, Urbana, IL 61801

Present address: \]University of Cincinnati, Cincinnati, OH 45221

Present address: \]Florida State University, Tallahassee, FL 32306

Present address: \]University of Illinois, Urbana, IL 61801

Present address: \]University of Illinois, Urbana, IL 61801

On leave from \]Kurchatov Institute, Moscow, Russia

FNAL E866/NuSea Collaboration

## Abstract

The Fermilab E866/NuSea Collaboration has measured the Drell-Yan dimuon cross sections in 800 GeV/$`c`$ $`pp`$ and $`pd`$ collisions. This represents the first measurement of the Drell-Yan cross section in $`pp`$ collisions over a broad kinematic region and the most extensive study to date of the Drell-Yan cross section in $`pd`$ collisions. The results indicate that recent global parton distribution fits provide a good description of the light antiquark sea in the nucleon over the Bjorken-$`x`$ range $`0.03x<0.15`$, but overestimate the valence quark distributions as $`x1`$.

Drell-Yan dimuon production provides valuable information about the partonic structure of hadrons that is complementary to deep-inelastic scattering (DIS) studies because it distinguishes between quarks and antiquarks. The Drell-Yan cross section for $`p`$A$`\mu ^+\mu ^{}X`$ may be written in terms of the parton distribution functions (PDFs) of the colliding hadrons as

$`M^3{\displaystyle \frac{d^2\sigma }{dMdx_F}}`$ $`=`$ $`{\displaystyle \frac{8\pi \alpha ^2}{9}}{\displaystyle \frac{x_1x_2}{x_1+x_2}}`$
$`\times `$ $`{\displaystyle \underset{q}{}}e_q^2\left[q_1(x_1)\overline{q}_2(x_2)+\overline{q}_1(x_1)q_2(x_2)\right]`$
where the subscript 1 (2) denotes the beam (target) hadron. Measurement of the invariant mass, $`M^2=x_1x_2s`$, and Feynman-$`x`$, $`x_F=2p_L/\sqrt{s}=x_1x_2`$, of the muon pair (where $`p_L`$ is the longitudinal momentum of the muon pair, and $`\sqrt{s}`$ the center-of-mass energy of the hadrons) determines the momentum fraction $`x_{1(2)}`$ of the beam (target) parton. Equation Absolute Drell-Yan Dimuon Cross Sections in 800 GeV/$`c`$ $`pp`$ and $`pd`$ Collisions is exact to leading order, but its general features are also preserved by next-to-leading-order (NLO) calculations.

The Fermilab E866/NuSea Collaboration has measured the Drell-Yan cross sections in 800 GeV/$`c`$ $`pp`$ and $`pd`$ collisions. This represents the first study of the Drell-Yan cross section in $`pp`$ collisions over a broad kinematic range and the most comprehensive measurement in $`pd`$ collisions. Previous publications Hawk98 ; Tow01  have described the E866 measurement of the $`\overline{d}/\overline{u}`$ ratio in the proton as a function of $`x`$, based on the $`x_2`$ dependence of the Drell-Yan cross section ratio $`\sigma ^{pd}/2\sigma ^{pp}`$, and those results have been included in recent global PDF fits cteq5 ; cteq6 ; mrst98 ; mrst01 ; grv98 .

We present measurements of the $`pp`$ and $`pd`$ Drell-Yan dimuon cross sections for pair mass in the ranges $`4.2<M<8.7`$ GeV or $`10.85<M<16.85`$ GeV and $`0.05<x_F<0.8`$. At the large values of $`x_F`$ appropriate for much of the E866 data, the first term in the sum in Eq. Absolute Drell-Yan Dimuon Cross Sections in 800 GeV/$`c`$ $`pp`$ and $`pd`$ Collisions dominates and the cross sections are primarily sensitive to the valence distributions in the proton beam and the antiquarks at small $`x`$ in the proton and deuteron targets. Thus, these measurements probe the partonic structure of the nucleon in two important kinematic regions. Recent global PDF fits use the precise HERA DIS cross sections to fix the magnitude of the light quark sea at small $`x`$ and the E605 $`p`$Cu cross sections e605  to fix the antiquark distributions for $`0.14x0.3`$. In contrast, Ref. mrst02,  noted that the magnitude of the antiquark distributions near $`x0.04`$ is relatively unconstrained by current data, and this limits the precision of calculations of the $`W`$ production cross section at the Tevatron. There has also been considerable recent interest in understanding the valence quark distributions in the nucleon as $`x1`$, which play a key role in searches for physics beyond the Standard Model, e.g., production of additional vector bosons as predicted by left-right symmetric models. To date, this region is constrained only by DIS cross sections, some of which suffer from limited statistics, while many others involve substantial ambiguities associated with corrections for the nuclear targets used Yang99 ; Kuhl00 . (See Ref. Zhang02,  for a recent review.) The present results provide new constraints on the magnitudes of the antiquark sea for $`0.03x<0.15`$ and of the valence quarks for $`x1`$.

E866 used a 3-dipole magnetic spectrometer employed previously in E605 e605 , E772 and E789, modified by the addition of new detectors at the first tracking station. An 800 GeV/$`c`$ proton beam bombarded identical target flasks containing liquid hydrogen, liquid deuterium and vacuum that were alternated every few minutes. After passing through the target, the beam was intercepted by a copper beam dump that was followed by a thick hadron absorber, ensuring that only muons traversed the spectrometer’s detectors. Drell-Yan events were recorded using three different spectrometer magnet settings, chosen to focus low-, intermediate- and high-mass muon pairs and provide acceptance from below the $`J/\psi `$ to above 15 GeV. A detailed description may be found in Ref. Tow01, .

The present analysis used 55,000 $`pp`$ and 121,000 $`pd`$ Drell-Yan events, approximately half the statistics of the E866 $`\overline{d}/\overline{u}`$ study Tow01 . The $`\overline{d}/\overline{u}`$ analysis was optimized to achieve minimum relative uncertainties in $`\sigma ^{pd}/2\sigma ^{pp}`$, whereas the present analysis was optimized for absolute measurements of the cross sections. Therefore, more stringent fiducial cuts were adopted in the present analysis, eliminating events for which the absolute acceptance of the spectrometer could not be reliably determined. This minimized systematic uncertainties at the cost of statistical precision, especially for $`x_20.15`$. In contrast, 7% of the events accepted in this analysis came from data sets that were excluded previously Tow01 ; jason . Several improvements were made to the event reconstruction, notably involving the treatments of energy loss and multiple scattering in the absorber, that led to a more precise reconstruction of the dimuon kinematics and a better match between the real and simulated events.

The spectrometer acceptance was calculated separately for each data set as a function of mass, $`x_F`$, and transverse momentum ($`p_T`$) using a detailed Monte Carlo simulation of the spectrometer. The virtual photon azimuthal production and decay angles were assumed to be distributed isotropically, and the polar decay angles were assumed to be distributed according to $`1+\mathrm{cos}^2\theta _d`$, consistent with theoretical expectations and previous Drell-Yan angular distribution studies McGau99 ; Chang99 ; Bro01 . Drell-Yan Monte Carlo events were thrown with realistic kinematic distributions in mass, $`x_F`$, and $`p_T`$, then reweighted to provide a precise match to the observed experimental distributions. Triply-differential cross sections, $`d^3\sigma /dMdx_Fdp_T`$, were calculated for each kinematic bin, then integrated over $`p_T`$ to obtain the invariant cross sections, $`M^3d^2\sigma /dMdx_F`$. The transverse momentum acceptance of the spectrometer extended to $`p_T57`$ GeV/$`c`$, but it was nonetheless necessary to extrapolate the integration to large $`p_T`$ where the Monte Carlo showed zero acceptance. The extrapolation contributed well under 1% in most cases, and was never more than 5%.

A detailed study was performed of the point-to-point systematic uncertainties in the measured cross sections. The two dominant contributions were the statistical uncertainty in the Monte Carlo event samples that were used for the acceptance calculations and the absolute field strength of the spectrometer magnets. Smaller contributions came from the uncertainties in the hodoscope, drift chamber, and trigger efficiencies, the composition and density of the targets, and the extrapolations to large $`p_T`$. Finally, some large-$`x_F`$ events contained a muon that passed very close to the edge of the beam dump, increasing the uncertainty in the energy loss and multiple scattering corrections. An additional $`\pm `$ 5% systematic uncertainty was added to the final results for the affected kinematic bins. The total point-to-point systematic uncertainties within any ($`M`$,$`x_F`$) bin were strongly correlated between the $`pp`$ and $`pd`$ cross sections. In addition to the point-to-point systematic uncertainties, there was a $`\pm `$ 6.5% overall normalization uncertainty, associated with the calibration of the beam intensity. See Ref. jason,  for additional details.

Figure 1 shows the Drell-Yan $`pp`$ and $`pd`$ invariant cross sections per nucleon for selected $`x_F`$ bins. The results agree with previous 800 GeV/$`c`$ Drell-Yan cross section measurements in $`p`$Cu collisions by E605 e605 . They also agree with previous measurements in $`pd`$ collisions by E772 e772  for $`x_F<0.3`$. At larger $`x_F`$ and small $`M`$, the E772 cross sections are systematically larger than the present results. E772 quoted a larger point-to-point systematic uncertainty in this kinematic region, but the results of the two experiments differ by more than the combined systematic uncertainties would predict. The largest differences are in the region $`0.4<x_F<0.6`$ and $`4.2<M<7`$ GeV, where inconsistencies between the E772 results and expectations from global PDF fits were noted previously mrst98 . This region was studied during E866 with very different acceptances by the low- and high-mass spectrometer settings, and the results are consistent.

Figure 1 also shows the results of next-to-leading-order calculations of the Drell-Yan cross sections based on the CTEQ6 cteq6  and MRST2001 mrst01  global PDF fits. The agreement with the global fits is very good over the entire kinematic region. This agreement may be quantified by computing a $`K^{}`$-factor, which we define to be the ratio of the experimental cross section to a NLO prediction. Table 1 shows $`K^{}`$-factors for several recent global PDF fits. With the exception of GRV98, all of the recent PDF fits predict the absolute magnitude of the Drell-Yan cross sections to within the $`\pm `$ 6.5% normalization uncertainty.

While the overall normalization is well reproduced, there are systematic deviations between the measurements and the predictions that are reflected in the large $`\chi ^2`$ values. To elucidate these deviations, it is useful to examine the experimental cross sections separately as functions of $`x_1`$ and $`x_2`$. Most of the events have $`x_1x_2`$, which implies that the $`x_1`$ dependence is primarily sensitive to the valence quarks in the proton beam and the $`x_2`$ dependence to the antiquarks in the target. Figure 2 shows the ratios of the experimental cross sections to NLO calculations using the MRST2001 global PDF fit, separately as a function of $`x_2`$ and $`x_1`$, averaged over the other momentum fraction. The uncertainties in the NLO calculations from the PDF fit are also shown. After accounting for the normalization uncertainty, the MRST2001 partons provide a very good description of the $`x_2`$ dependence of the $`pp`$ cross sections over the full range, and a good description of the $`pd`$ cross sections for $`x_2<0.15`$. The CTEQ6 fits describe the $`x_2`$ dependence equally well. This indicates that the current PDFs interpolate the $`\overline{u}`$ and $`\overline{d}`$ distributions successfully between the HERA measurements at small $`x`$ and the E605 measurements for $`x>0.14`$, but these new results provide tighter constraints on the magnitude of the antiquark distributions for $`0.03x<0.15`$ than have existed to date mrst02 .

Figure 2b shows the ratios of the cross sections to NLO calculations using the MRST2001 PDF fits as a function of $`x_1`$. The qualitative behavior of the CTEQ6 distributions is quite similar: both PDF fits overestimate the valence quarks by $``$ 15–20% at large $`x_1`$, and the discrepancy between the data and current PDFs appears to be larger for the $`pd`$ cross sections than for the $`pp`$ cross sections. The acceptance of these data includes an anti-correlation between $`x_1`$ and $`x_2`$ Tow01 , so most of the events at large $`x_1`$ also have $`x_2<0.1`$. The NLO calculations describe the average $`x_2`$ dependence of the cross sections quite well for $`x_2<0.1`$, so the discrepancy in Fig. 2b is isolated to the subset of the low-$`x_2`$ events that also have large $`x_1`$. Further comparisons between the measured cross sections and the NLO calculations, in which the $`x_1`$ and $`x_2`$ dependences are examined for limited ranges of the other variable jason , reinforce this conclusion. Table 2 shows the results of fits to the $`x_1`$ dependence of the $`K^{}`$-factors with a phenomenological form, $`K^{}(x_1)=\alpha \left(1x_1\right)^\beta `$, for both PDF sets.

Unlike large-$`x`$ DIS off deuterium targets, these data are in a kinematic region where nuclear-dependence corrections are known to be small Alde90 . The $`pp`$ and $`pd`$ cross sections at large $`x_1`$ constrain slightly different linear combinations of $`u(x)`$ and $`d(x)`$, with greater sensitivity to $`u`$ quarks. The results imply that the $`u`$ quark distributions in CTEQ6 and MRST2001 are overestimated as $`x1`$. Recent fits to H1 charged- and neutral-current data have also indicated that the global PDF fits appear to overestimate the $`u`$ quark distribution near $`x0.65`$, albeit with limited statistics Zhang02 . The Drell-Yan cross sections may also point to problems with the $`d/u`$ ratio as $`x1`$ Yang99 ; Kuhl00 . However, determination of the $`d/u`$ ratio at large $`x`$ from the present data will require a fit to the full two-dimensional invariant cross sections together with the rest of the current world data.

Figure 2 shows that the present results fall well outside the experimental uncertainty mrst02  for the valence quarks in the MRST2001 PDF fit. A similar comparison shows that the present results track the lower limit on the experimental uncertainty cteq6  in the CTEQ6 PDF fit. These discrepancies reflect additional uncertainties in the PDF fits cteq6 ; mrst02 , and they imply that future PDF fits will see a substantial correction to the $`u`$ and $`d`$ quark distributions at large $`x`$. The reduced valence quark distributions at large $`x`$ implied by these results will also lead to predictions of smaller cross sections for some signatures for new physics that are being sought at the Tevatron.

We thank the Fermilab Particle Physics, Beams, and Computing Divisions for their assistance in performing this experiment. We also thank W.K. Tung of the CTEQ Collaboration for providing the code to calculate the next-to-leading-order Drell-Yan cross sections. This work was supported in part by the U.S. Department of Energy.


text:
# Nothing but Relativity

## Acknowledgements:

I thank P. Bhattacharjee, B. P. Das, J. Samuel, D. Sen and S. Sinha for patiently listening to my arguments, commenting on them, and suggesting some references which I could have missed otherwise.