text:
# Driven vortices in 3D layered superconductors: Dynamical ordering along the c-axis

\[

## Abstract

We study a 3D model of driven vortices in weakly coupled layered superconductors with strong pinning. Above the critical force $`F_c`$, we find a plastic flow regime in which pancakes in different layers are uncoupled, corresponding to a pancake gas. At a higher $`F`$, there is an “smectic flow” regime with short-range interlayer order, corresponding to an entangled line liquid. Later, the transverse displacements freeze and vortices become correlated along the c-axis, resulting in a transverse solid. Finally, at a force $`F_s`$ the longitudinal displacements freeze and we find a coherent solid of rigid lines.

\]

It is well-known that an external current can induce an ordering of the vortex structure in superconductors with pinning . For a long time, it was believed that the high-current phase would have crystalline order. Recently, it has been found that different kinds of order are possible at high currents, depending on pinning strength and dimensionality . This has led to numerous theoretical , experimental  and numerical studies . A crystal-like structure, which could be either a perfect crystal  or a Bragg glass , is only possible in $`d=3`$ at large drives. In $`d=2`$, or in $`d=3`$ for intermediate currents, a transverse glass is expected, with order only in the direction perpendicular to the driving force . In the equilibrium vortex phase diagram, the behavior of vortex line correlations along the direction of the magnetic field ($`c`$-axis) has been intensively discussed both experimentally  and theoretically . In the case of driven vortices, little is known on how the $`c`$-axis line correlations would behave in the different dynamical regimes. Here we will address this issue starting from the less favorable case: weakly coupled superconducting planes with strong pinning. We will show how the order along the $`c`$-axis and the in-plane structural order take place in a sequence of dynamical phases upon increasing current.

We study pancake vortices in a layered superconductor, considering the long-range magnetic interactions between all the pancakes and neglecting Josephson coupling . This model is adequate when the interlayer periodicity $`d`$ is much smaller than the in-plane penetration length $`\lambda _{}`$ . Previous simulations of driven vortices in 3D superconductors have been performed using Langevin dynamics of short-range interacting particles  or the driven isotropic 3D XY model .

The equation of motion for a pancake located in position $`𝐑_i=(𝐫_i,z_i)=(x_i,y_i,n_id)`$, ($`𝐳\widehat{c}`$), is:

$$\eta \frac{d𝐫_𝐢}{dt}=\underset{ji}{}𝐅_𝐯(\rho _{ij},z_{ij})+\underset{p}{}𝐅_𝐩(\rho _{ip})+𝐅,$$
(1)
where $`\rho _{ij}=|𝐫_i𝐫_j|`$ and $`z_{ij}=|z_iz_j|`$ are the in-plane and inter-plane distance between pancakes $`i,j`$, $`\rho _{ip}=|𝐫_i𝐫_p|`$ is the in-plane distance between the vortex $`i`$ and a pinning site at $`𝐑_𝐩=(𝐫_p,z_i)`$, $`\eta `$ is the Bardeen-Stephen friction, and $`𝐅=\frac{\mathrm{\Phi }_0}{c}𝐉\times 𝐳`$ is the driving force due to an in-plane current $`𝐉`$. We consider a random uniform distribution of attractive pinning centers in each layer with $`𝐅_𝐩=2A_pe^{(\rho /a_p)^2}𝐫/a_p^2`$, where $`a_p`$ is the pinning range. The magnetic interaction between pancakes $`𝐅_v(\rho ,z)=F_\rho (\rho ,z)\widehat{r}`$ is given by :

$`F_\rho (\rho ,0)`$ $`=`$ $`{\displaystyle \frac{A_v}{\rho }}\left[1{\displaystyle \frac{\lambda _{}}{\mathrm{\Lambda }}}\left(1e^{\rho /\lambda _{}}\right)\right]`$ (2)
$`F_\rho (\rho ,z_n)`$ $`=`$ $`{\displaystyle \frac{\lambda _{}}{\mathrm{\Lambda }}}{\displaystyle \frac{A_v}{\rho }}\left[e^{|z_n|/\lambda _{}}e^{R_n/\lambda _{}}\right].`$ (3)
Here, $`R=\sqrt{z^2+\rho ^2}`$ and $`\mathrm{\Lambda }=2\lambda _{}^2/d`$ is the 2D thin-film screening length. An analogous model to Eqs. (2-3) was used in . We normalize length scales by $`\lambda _{}`$, energy scales by $`A_v=\varphi _0^2/4\pi ^2\mathrm{\Lambda }`$, and time is normalized by $`\tau =\eta \lambda _{}^2/A_v`$. We consider $`N_v`$ pancake vortices and $`N_p`$ pinning centers per layer in $`N_l`$ rectangular layers of size $`L_x\times L_y`$, and the normalized vortex density is $`n_v=B\lambda _{}^2/\mathrm{\Phi }_0=(a_0/\lambda _{})^2`$. We consider $`n_v=0.29`$ with $`L_y=16\lambda _{}`$ and $`L_x=\sqrt{3}/2L_y`$, $`N_l=8`$ and $`N_v=64`$. We take a pinning range of $`a_p=0.2`$, a large pinning strengh of $`A_p/A_v=0.2`$, with a high density of pinning centers $`n_p=3.125n_v`$. The model of Eq.(2-3) is valid in the limit $`d\lambda _{}\mathrm{\Lambda }`$. We take $`d/\lambda _{}=0.01`$, which corresponds to BSCCO compounds . Moving pancake vortices induce a total electric field $`𝐄=\frac{B}{c}𝐯\times 𝐳`$, with $`𝐯=\frac{1}{N_vN_l}_i𝐯_i`$. We study the dynamical regimes in the velocity-force curve at $`T=0`$, solving Eq. (1) for increasing values of $`𝐅=F𝐲`$ . We use periodic boundary conditions both in the planes and in the $`z`$ direction and interactions between all pancakes in all layers are considered . The periodic long-range in-plane and inter-plane interaction is evaluated using Ref. . The equations are integrated with a time step of $`\mathrm{\Delta }t=0.01\tau `$ and averages are evaluated in $`16384`$ integration steps after $`2000`$ iterations for equilibration. Each simulation is

started at $`F=0`$ with a triangular vortex lattice and slowly increasing the force in steps of $`\mathrm{\Delta }F=0.1`$ up to values as high as $`F=8`$.

We start with a qualitative description of the different steady states that arise as a function of increasing force. In Figure 1(a-d) we show the vortex trajectories $`\{𝐑_i(t)\}`$ for typical values of $`F`$ by plotting the positions of the pancakes in five of the layers for all $`t`$. In Fig.1(e-h) we show the average in-plane structure factor $`S(𝐤)=\frac{1}{N_l}_n|\frac{1}{N_v}_i\mathrm{exp}[i𝐤𝐫_{ni}(t)]|^2`$, with $`𝐤=(k_x,k_y)`$. Above the depinning critical force $`F_c`$, we find the following dynamical phases. (i)Plastic flow ($`F_c<F<F_p`$): Pancakes flow in an intrincate network of “plastic” channels similar to the behavior found in 2D . The motion in different planes is completely uncorrelated, \[Fig.1(a)\] and there is no signature of order in the structure factor \[Fig.1(e)\]. (ii)Smectic flow ($`F_p<F<F_t`$): The motion organizes in “elastic” channels that are almost parallel and separated by a distance $`a_0`$, see Fig.1(b). Small and broad “smectic” peaks appear in $`S(𝐤)`$ for $`𝐤𝐅=0`$ \[Fig.1(f)\]. There are “activated” jumps of pancakes between channels. Along the $`c`$-direction the channels tend to align sitting on top of

each other between neighboring planes. (iii)Transverse solid ($`F_t<F<F_s`$): There are well defined channels in all the planes and the pancakes do not jump between channels \[Fig.1(c)\]. The structure factor has sharp smectic peaks and small “longitudinal” peaks ($`𝐤𝐅0`$) have appeared \[Fig.1(g)\]. The location of channels is correlated in the $`c`$-axis. (iv)Coherent solid ($`F>F_s`$): The channels become more straigth with small transverse wandering \[Fig.1(d)\]. The $`S(𝐤)`$ shows well defined peaks for all $`𝐤`$ in the reciprocal lattice \[Fig.1(h)\]. The dynamical phases (i)-(iii) are similar to the ones found previously in 2D thin films .

The characteristic forces $`F_c,F_p,F_t,F_s`$ separating the different dynamical phases are obtained from the analysis of the in-plane and out of plane structural and dynamical correlations. We show in Fig.2 the in-plane structure factor and temporal fluctuations, obtained in the same way as for 2D . In Fig.2(a) we plot the average velocity $`V=V_y(t)=\frac{1}{N_vN_l}_{n,i}\frac{dy_{ni}}{dt}`$, in the direction of the force as a function of $`F`$ and its corresponding derivative $`dV/dF`$ (differential resistance). The force $`F_p`$ corresponds to the peak in the differential resistance. We also see a small second maximum in $`dV/dF`$ for a force between $`F_t`$ and $`F_s`$ . In Fig.2(b) we plot the magnitude of the peaks in the in-plane structure factor. We show the peak height at $`G_1=2\pi /a_0\widehat{x}`$, corresponding to smectic ordering, and the average of the peaks corresponding to longitudinal ordering at $`𝐆_2=\pm 2\pi /a_0(1/2,\sqrt{3}/2)`$ and $`𝐆_3=\pm 2\pi /a_0(1/2,\sqrt{3}/2)`$. We see that at $`F_p`$ the smectic peak rises up from zero, then at $`F_t`$ it reaches an almost constant value and later at $`F_s`$ it has a small jump. The longitudinal peak has a small finite value for forces above $`F_p`$, and only at $`F_s`$ shows a significant increment towards a large value. Comparing with the previous 2D results , we can make the reasonable assumption that for $`F_p<F<F_t`$ there is only short-range smectic order (since there is activated transverse diffusion beteween elastic channels, see below), for $`F_t<F<F_s`$ there is probably quasi-long range smectic order but short range longitudinal order, and above $`F_s`$ there is both transversal and longitudinal order (quasilong-range or long-range). What is new, compared with the 2D thin film case , is that above a force $`F_s`$ there is a significant amount of longitudinal order. This may correspond either to a moving crystalline phase (if there is long-range order) or to a moving Bragg glass (if there is quasilong-range order) . We have verified that, for a given $`F>F_s`$, there is both longitudinal and transversal order for system sizes of $`N_l\times N_v=5\times 36,5\times 64,8\times 64,8\times 100,10\times 100`$. However, a detailed finite size analysis is not possible with these few small samples. We complement our discussion of the in-plane physics with the study of the temporal fluctuations, which are shown in Fig.2(c). We calculate the transverse diffusion coefficient $`D_x`$ from the average quadratic transverse displacements of vortices from their center of mass position $`(\overline{X}_n,\overline{Y}_n)`$, $`\frac{1}{N_vN_l}_i[x_i(t)\overline{X}_{n_i}(t)x_i(0)+\overline{X}_{n_i}(0)]^2D_xt`$. We find that $`D_x`$ is maximum at $`F_p`$ in coincidence with the peak in $`dV/dF`$. Below $`F_p`$ diffusion is through the intrincate network of plastic channels, above $`F_p`$ diffusion is through activated jumps between elastic channels. $`D_x`$ sharply drops to zero at $`F_t`$, indicating that transverse displacements are localized in the transverse solid phase . The drift from the center of mass of longitudinal displacements $`[\mathrm{\Delta }y(t)]^2=[y_i(t)\overline{Y}_{n_i}(t)y_i(0)+\overline{Y}_{n_i}(0)]^2`$ is superdiffusive for $`F<F_s`$, similar to the results observed in 2D films . For $`F>F_s`$ the longitudinal displacements become frozen in a constant value $`[\mathrm{\Delta }y(t)]^2<a_0/N_l`$, as it is shown in Fig.2(c). Since in-plane displacements are localized and there are large transversal and longitudinal Bragg peaks, we call this phase a coherent solid.

Let us now discuss how the ordering along the $`c`$-axis takes place. We analyze the pair distribution function: $`g(\rho ,n)=\frac{L_xL_y}{N_v}_{ij}\delta (\rho \rho _{ij})\delta _{n,n_{ij}}`$. From $`g(\rho ,n)`$ we define a correlation function along c-axis $`C_z(n)=lim_{\rho 0}g(\rho ,n)`$. Short-range ordering will be given by a

finite $`C_z(n=1)`$, meaning that pancakes in neighboring planes are coupled and a “vortex line” can therefore be defined. In principle, an exponential decay $`C_z(n)\mathrm{exp}(n/\xi _z)`$ would define a correlation length for the vortex line . On the other hand, long-range ordering will be given by $`C_z(n\mathrm{})C_z^{\mathrm{}}>0`$. In Fig.3(a) we show $`C_z(n)`$ as a function of $`F`$ for $`n=1,2,3,4`$. We see that at $`F_p`$ there is an onset of short-range order along the c-axis with a finite $`C_z(n=1)`$. At higher forces between $`F_p`$ and $`F_t`$ the other $`C_z(n>1)`$ start to rise. The absence of correlations for $`F<F_p`$ means that pancake motion is completely random and uncorrelated between different planes. Therefore, we propose that the plastic flow regime corresponds to a pancake gas. Above $`F_p`$, in the smectic flow regime, it is possible to define a vortex line with short range correlations along the c-axis. Since there are in-plane jumps between elastic channels (i.e., cutting and reconnection of flux lines) we may consider this phase as an entangled line liquid. Above $`F_t`$, $`C_z(n)`$ is finite for all $`n`$ considered and tends to saturate upon increasing $`n`$. This indicates that vortex lines become more stiff above $`F_t`$. We also analyzed the $`c`$-axis correlation between averaged vortex densities. We first define $`\rho _v(𝐫,n,t)=\frac{1}{N_v}_i\delta (𝐫𝐫_{ni}(t))`$ taking a coarse-graining scale $`\mathrm{\Delta }r=a_0/2`$ (results do not vary much for $`\mathrm{\Delta }r=a_0/4`$). The regions where the average density $`\rho _v(𝐫,n)`$ is large define the paths of steady state vortex motion. We can thereby calculate the overlap function of vortex trajectories between different planes as $`O_n=C_\rho (n)/C_\rho (0)`$, with $`C_\rho (n)=\frac{L_xL_y}{N_l}[_m𝑑𝐫\rho _v(𝐫,m)\rho _v(𝐫,m+n)]1`$. This is shown in Fig.3(b). We see that $`O_n`$ also has an onset at $`F_p`$. For $`F_p<F<F_t`$, we have some overlap of the elastic channels that decreases with increasing $`n`$, consistent with the entangled line-liquid picture. More interestingly, at $`F_t`$ the overlap function $`O_n`$ becomes independent of $`n`$. This means that there is long-range $`c`$-axis coupling of the path of the elastic channels. When transverse displacements become localized in the $`x`$-direction, they also become locked in the $`c`$-direction. Thus, the freezing of in-plane transverse displacements occurs simultaneously with a transverse disentanglement of flux lines at $`F_t`$. A striking result is that we find $`O_n1`$ above $`F_s`$, i.e., a perfect c-axis coupling of elastic channels (within the scale $`a_0/4`$). Another interesting point to consider is the correlation of vortex velocities. If vortices in different planes move at different velocities, they will induce a Josephson voltage difference along the $`c`$-axis given by $`V_{n,n+1}(𝐫,t)=\frac{\mathrm{\Phi }_0}{2\pi c}\frac{d}{dt}\varphi _{n,n+1}(𝐫,t)`$, with $`\varphi _{n,n+1}`$ the superconducting phase difference between planes $`n`$ and $`n+1`$. A good approximation for pancakes at $`𝐫_{n,i}`$ is to write $`\varphi _{n,n+1}(𝐫,t)=_i[f(𝐫𝐫_{n,i})f(𝐫𝐫_{n+1,i})]`$ with $`f(𝐫)\mathrm{arctan}(x/y)`$. We can therefore estimate the c-axis voltage fluctuations as $`\delta ^2V_c=_n[V_{n,n+1}^2(𝐫,t)V_{n,n+1}(𝐫,t)^2]𝑑𝐫A_n[𝐕_n^2𝐕_n^2][𝐕_n𝐕_{n+1}𝐕_n𝐕_{n+1}]`$; with $`𝐕_n(t)=\frac{1}{N_v}_i𝐯_{n,i}(t)`$, and the constant $`A\mathrm{log}\mathrm{\Lambda }`$ if $`L>\mathrm{\Lambda }`$ or $`A\mathrm{log}(L)`$ otherwise. It is clear that $`\delta ^2V_c=0`$ for pancakes moving with the same velocity in all planes. We see in Fig.3(c) that the voltage fluctuations have a maximum at $`F_p`$. For $`F>F_p`$, $`\delta ^2V_c`$ decreases, and above $`F_s`$ it reaches an almost $`F`$-independent value. The fact that $`\delta ^2V_c`$ does not vanish above $`F_s`$ is consistent with the result that $`C_z(n)<1`$ for all values of $`F`$ in Fig.3(a). In other words, while transverse displacements are strongly correlated along the $`c`$ direction for large forces \[Fig.3(b)\], the longitudinal displacements in different planes are weakly correlated .

In conclusion, we have clearly distinguished different dynamical phases in 3D layered superconductors considering both in-plane and $`c`$-axis ordering . The onset of short-range $`c`$-axis correlations could be studied experimentally with plasma resonance measurements . The long-range ordering along the $`c`$-axis could be studied through simultaneous measurements of $`\rho _c`$ resistivity and in-plane current-voltage response .

We acknowledge discussions with L.N. Bulaevskii, P.S. Cornaglia, F. de la Cruz, Y. Fasano, M. Menghini. This work has been supported by ANPCYT (Proy. 03-00000-01034), by Fundación Antorchas (Proy. A-13532/1-96), Conicet, CNEA and FOMEC (Argentina); by CLC and CULAR (Los Alamos), and by the Director, Office of Adv. Sci. Comp. Res., Division of Mathematical, Information, and Computational Sciences of the U.S.D.O.E. (contract number DE-AC03-76SF00098).


format: text

text:
# Two-magnon Raman scattering in spin-ladder geometries and the ratio of rung and leg exchange constants

## I Introduction

Determining the exchange constants and establishing the microscopic spin Hamiltonian is an important step in understanding magnetic properties of exotic new materials. In materials which have predominantly a single exchange constant, the measurement of uniform susceptibility and the determination of Curie-Weiss parameter is sufficient to obtain the exchange constant. However, in recent years, many complex materials have been synthesized, which have more than one exchange constant. Furthermore, these exchange constants can be so large that the Curie-Weiss regime may not be experimentally accessible. In these cases alternative methods are needed to determine the exchange parameters. Examples of such materials are cuprate based spin-ladder materials, cousins of high temperature superconducting materials, for which there is a substantial literature for determining the exchange constants. The simplest methods for determining the exchange constants involve measuring the temperature dependence of uniform susceptibility, or Knight shift, or nuclear relaxation rates, which can then be compared with detailed theoretical calculations to obtain the exchange parameters. Such a fitting procedure is not very accurate, as reflected in the range of values that exist in the literature for these materials.

Recently, Sugai  noted that Raman scattering can be used to determine the ratio of rung and leg exchange constants in spin-ladder materials. He argued that by varying the polarization direction of incident and outgoing light, one might be able to shift the Raman spectra in ways that can be related to the different exchange constants. Our study is motivated by Sugai’s work. However, we find that the arguments used by Sugai to relate the position of the spectral peaks to the different exchange constants in the spin-ladder materials are incorrect. Within the Fleury-Loudon-Elliott  theory the spectral lineshape does not change with polarization at all. This gives a simple explanation why Sugai always finds a very small shift in the Raman spectra and sheds doubt on his inference that the exchange constants in all the ladder materials are close to unity.

It is well known that an accurate calculation of Raman spectra for low-dimensional spin-half antiferromagnets is very difficult due to quantum fluctuations.  And, short of a direct comparison of the spectra with theory, it would seem that Raman spectra cannot be used to determine the ratio of exchange constants. However, we show below that the Raman scattering intensities do depend on polarization geometries in a way that is easily calculated and related to the ratio of exchange constants. Unfortunately, they come with an unknown factor, about which certain assumptions need to be made before an estimate of the ratio of rung to leg exchange constants in the ladder material can be obtained.

## II Spin-Ladder Heisenberg Models

We begin with a system described by a Heisenberg model in ladder geometry, with the Hamiltonian,

$$H=J_r\underset{ij,r}{}\stackrel{}{S}_𝐢\stackrel{}{S}_𝐣+J_l\underset{ij,l}{}\stackrel{}{S}_𝐢\stackrel{}{S}_𝐣,$$
(1)
where $`J_r`$ and $`J_l`$ denote the coupling constants for the rung and leg bonds of the ladder, respectively, and the sums are only over the bonds in the indicated directions. Our primary goal is to examine the dependence of Raman scattering on the ratio of exchange constants $`J_r`$ and $`J_l`$.

Within the Fleury-Loudon-Elliott approach, magnetic Raman scattering is described by an effective Raman Hamiltonian or operator :

$$_R=\underset{ij}{}J_{ij}^{}(\widehat{ϵ}_{in}\widehat{r}_{ij})(\widehat{ϵ}_{out}\widehat{r}_{ij})\stackrel{}{S}_𝐢\stackrel{}{S}_𝐣,$$
(2)
where the $`\widehat{r}_{ij}`$ are unit vectors along the bond directions, and $`\widehat{ϵ}_{in}`$ and $`\widehat{ϵ}_{out}`$ are unit vectors indicating the direction of polarization of the incident and scattered light, respectively. The $`J_{ij}^{}`$ are constants representing the strength of the Raman scattering interaction between spins $`i`$ and $`j`$. In previous studies of the spin-ladder geometry , it has been assumed that $`J_{ij}`$ is nearest-neighbor and is a constant for all nearest-neighbor bonds. Thus it can be taken out of the summation. Such a constant simply sets the overall scale for the scattering and does not influence any other result. However, we believe that if the rung and leg exchange constants are not equal, a priori, we cannot assume that the ratio of $`J^{}`$ along the rung and leg directions to be equal. Thus we proceed here with a more general $`J_{ij}^{}`$, and later consider possible scenarios for their values, which would play an important role.

Following Sugai, it is most useful to consider the case where the incident and scattered light have parallel polarization directions, both lying in the plane of our 2D system. Thus, $`\widehat{ϵ}_{in}=\widehat{ϵ}_{out}`$. Most generally, we can denote the polarization with an angle $`\theta `$ with respect to the vertical bonds, which makes the effective Raman Hamiltonian

$$_R(\theta )=\mathrm{cos}^2\theta \underset{ij,r}{}J_{ij}^{}\stackrel{}{S}_𝐢\stackrel{}{S}_𝐣+\mathrm{sin}^2\theta \underset{ij,l}{}J_{ij}^{}\stackrel{}{S}_𝐢\stackrel{}{S}_𝐣.$$
(3)
As a notational simplification, we define $`_r=_R(0)`$ and $`_l=_R(\frac{\pi }{2})`$. Also, for reasons of symmetry, we assume that all the $`J_{ij}^{}`$ in each summation are the same (we call them $`J_r^{}`$ and $`J_l^{}`$, respectively) so

$$_R(\theta )=J_r^{}(\mathrm{cos}^2\theta )_r+J_l^{}(\mathrm{sin}^2\theta )_l.$$
(4)
The two-magnon Raman scattering intensity as a function of frequency can be expressed using Fermi’s golden rule,

$$I(\omega ,\theta )=\underset{n}{}^{}|\psi _n|_R(\theta )|\psi _0|^2\delta (\omega (E_nE_0)),$$
(5)
where $`|\psi _n`$ and $`E_n`$ are eigenvectors and eigenvalues of the Hamiltonian, and the prime indicates that the ground state is excluded from the sum.

Our key result follows from the following simple consideration: using the terminology defined above, we can express Eq. (1) as

$$_l=\frac{1}{J_l}(HJ_r_r),$$
(6)
which can be substituted into Eq. (4) to give

$$_R(\theta )=J_r^{}(\mathrm{cos}^2\theta \frac{J_rJ_l^{}}{J_lJ_r^{}}\mathrm{sin}^2\theta )_r+\frac{J_l^{}\mathrm{sin}^2\theta }{J_l}H.$$
(7)
The second term of the sum is a multiple of the Hamiltonian, and thus cannot contribute to the scattering. The first term is proportional to $`_r`$ for all angles. Thus, within this theory, the observed two-magnon scattering spectrum will have the same line shape and peak position for all angles. This result is in direct contradiction with the arguments of Sugai.

The intensity of the spectrum is given by the expression

$$I(\omega ,\theta )=(\mathrm{cos}^2\theta \frac{J_rJ_l^{}}{J_lJ_r^{}}\mathrm{sin}^2\theta )^2I(\omega ,0).$$
(8)
This variation of angle can be used to determine the ratio of exchange constants, provided one either knows $`J_r^{}/J_l^{}`$ or can relate it to $`J_r/J_l`$. This point is discussed a little later. In principle, one can perform the experiments by varying $`\theta `$ continuously, to obtain the above variation. One simply needs to keep the polarization direction of incoming and outgoing light fixed parallel to each other in the plane and rotate the sample. To find the ratio of exchange constants, it is sufficient to consider two angles:

$$\frac{J_r}{J_l}=\frac{J_r^{}}{J_l^{}}\sqrt{\frac{I(\omega ,\pi /2)}{I(\omega ,0)}}.$$
(9)
A simple procedure for determining the ratio of coupling constants could be to measure the maximum two-magnon Raman scattering intensity, then rotate the material through an angle of 90 and measure the intensity again. The ratio of intensities would give the ratio of exchange constants, provided that one can reasonably estimate the value of the ratio of $`J_r^{}/J_l^{}`$.

## III Analysis of Experimental Results

Experimental evidence does support the idea that the spectra do not vary much in shape as a function of angle, but do vary in intensity. For example, Sugai and Suzuki  measured the Raman spectra for the spin-1/2 two-leg ladder materials La6Ca8Cu24O41 and Sr14Cu24O41 in two different configurations separated by an angle of 90. One configuration had the incident and scattered polarizations parallel to the legs of the ladder, and the other had them parallel to the rungs. The peak of the observed spectra for the two configurations of La6Ca8Cu24O41 were found at 3004 cm-1 and 2948 cm-1, a relative difference of less than 2% which can be neglected. The spectral shapes are largely identical as well, but the intensities are vastly different. The ratio of the intensities of the peaks of these two spectra is approximately 0.52, which is quite significant. The same measurements performed on Sr14Cu24O41 yielded peaks at 3006 cm-1 and 3004 cm-1, and an intensity ratio of approximately 0.39. The peaks are even closer together than the others studied, and the intensity ratio is smaller. Popović and collaborators  present scattering data for Sr14Cu24O41 that shows a similar lack of peak location shift, as do measurements by Sugai and collaborators  on LaCuO2.5. The small observed shift in the spectral peak could be due to impurities, phonons, inter-ladder couplings as well as due to other intra-ladder interactions not included here.  Note that because our result derives from an operator relation, it is not sensitive to long range ordering in the system, and should be valid for small interladder couplings.

There are several possiblities for the value of $`J_r^{}/J_l^{}`$ that should be considered. We shall begin by examining the possibility that $`J_r^{}/J_l^{}=J_r/J_l`$, as follows from a large-U perturbative treatment and suggested by the work of Moriya  and Shastry and Shraiman.  If we assume that $`J_r^{}/J_l^{}=J_r/J_l`$, Eq. (8) becomes $`I(\omega ,\theta )=(\mathrm{cos}2\theta )^2I(\omega ,0)`$. Thus, all of the maxima should be identical in magnitude. This conclusion is directly contradicted by all of the experimental results discussed previously, where the ratio was found to vary by more than a factor of 2. Clearly, direct proportionality does not exist.

A second possibility is that all of the nearest-neighbor Raman operator coupling constants are equal to one another. This has been assumed universally by almost all previous studies of these ladder materials. The bonds along the rungs and legs of the ladder are nearly identical, which is perhaps what led to this assumption. However, it leaves open the question of why this ratio will remain unity when the ratio $`J_r/J_l`$ deviates significantly from unity. Assuming $`J_r^{}/J_l^{}1`$, Eq. (9) becomes

$$\frac{J_r}{J_l}=\sqrt{\frac{I(\omega ,\pi /2)}{I(\omega ,0)}}.$$
(10)
A third possibility can be motivated by perturbation theory. If the two-magnon Raman process involves a direct exchange it is a second order process, whereas the usual superexchange could be mediated by non-magnetic intermediate ions and thus could be a fourth order process in a large U expansion. In this case, a more natural relationship between the Raman and Heisenberg coupling constants is $`J^{}\sqrt{J}`$. So,

$$\frac{J_r}{J_l}=\frac{I(\omega ,\pi /2)}{I(\omega ,0)}.$$
(11)
Using these two possibilities, we computed the ratio of the Heisenberg coupling constants for the materials studied by Sugai and his collaborators, using his published data. The results for these materials are shown in Table I.

Let us discuss these results in light of previous studies . Some authors have adopted the point of view that this ratio is close to unity, and have used that as the starting point of their analysis . On the other hand, some local density approximation calculations  find the ratio to be closer to $`0.5`$. A number of other studies which allow the ratio to vary, also find values close to $`0.5`$ , although a range of experimental values from 0.5 to 1.13 have been quoted . Brehmer et al.  make an interesting attempt to reconcile the conflicting viewpoints, keeping the ratio to be unity but allowing an additional biquadratic ring-interaction in the Hamiltonian.

In the Sugai and Suzuki paper, the authors used the incorrect argument for the energy shift to estimate the ratio of exchange couplings along the rungs and the legs of the ladder. Their ratio was determined to be 0.95 for La6Ca8Cu24O41 and 1 for Sr14Cu24O41. It is now clear that Sugai and collaborators obtained values close to unity because, to a good approximation, the spectra do not shift at all with change in polarization direction. Using their data and Eq. (10) and Eq. (11), we calculate the ratio to be less than one in all cases. It is interesting to note that using Eq. (11), our calculated ratio agrees well with the conclusion that $`J_r/J_l0.5`$.

## IV Conclusion

We have seen how Fleury-Loudon-Elliott theory predicts that the shape of the Raman spectra in the spin-ladder geometry, with nearest-neighbor Heisenberg exchange constants, does not change with the polarization directions. Instead, there exists a relationship between the intensity of two-magnon Raman scattering in different polarizations and the ratio of Heisenberg exchange constants. With some suitable assumptions about the Raman hamiltonian, we have used it to estimate the ratio of rung to leg exchange constant in several cuprate materials. The full dependence of the intensity on the polarization direction can be experimentally verified and should serve as a test for the Fleury-Loudon-Eliott theory.

###### Acknowledgements.

The authors would like to thank T. Imai and A. Slepoy for their helpful comments. This work was supported in part by the National Science Foundation under grant number DMR-9986948.


text:
# 1 Introduction

## 1 Introduction

Magnetic monopoles (MMs), nuclearites (nuggets of strange quark matter + electrons) and Q-balls (supersymmetric states of squarks, sleptons and Higgs fields) could be part of the cold dark matter (DM) located in the halo of galaxies. They would have typical velocities of $`\beta c10^3c`$; galactic magnetic fields may accelerate magnetic monopoles to larger velocities.

In this note we show some qualitative pictures of the structures of (i) superheavy GUT magnetic monopoles, (ii) intermediate mass magnetic monopoles, (iii) nuclearites and (iv) Q-balls. We discuss in particular the relation between the mass and the size of these particles.

Based on astrophysical considerations, we discuss phenomenological limits on the fluxes in the cosmic radiation, of magnetic monopoles, nuclearites and Q-balls.

## 2 Magnetic monopoles

The concept of magnetic monopole may be traced back to the origin of magnetism. In 1931 Dirac introduced the MMs in order to explain the quantization of the electric charge, obtaining the formula $`eg=n\mathrm{}c/2`$, from which $`g=ng_D=n\mathrm{}c/2e=68.5e=3.2910^8`$ c.g.s. symmetric system . A MM possessing also an electric charge is called a dyon. A MM and an atomic nucleus may form a bound system with both magnetic and electric charges: also this system is called a dyon.

The energy losses of MMs and of dyons in matter, in the Earth and in different detectors were calculated in ref. . An extensive bibliography of MMs is given in ref. .

### 2.1 GUT magnetic monopoles

In the context of the Grand Unified Theory (GUT) of electroweak and strong interactions, MMs with magnetic charges $`g=ng_D`$ are predicted to appear at the cosmic time of $`10^{34}`$ s, during the phase transition corresponding to the spontaneous breaking of the Grand Unified group . They should have been produced as point defects or in extremely high energy collisions of the type $`e^+e^{}M\overline{M}`$ . One of the main problems with GUT MMs is the too large abundance predicted by the standard cosmology. Models with inflation at the end of the GUT epoch reduce dramatically their number and we would be left mainly with MMs produced in very high energy collisions.

The GUT MM mass is related to the mass of the X vector bosons, carriers of the unified interaction, by the relation $`M_MM_X/\alpha `$, where $`\alpha 0.025`$ is the dimensionless unified coupling constant and $`M_X10^{14}÷10^{16}`$ GeV/c2. Thus GUT magnetic monopoles should have $`M_M10^{16}`$ GeV/c2. Due to their large masses, these MMs cannot be produced with existing and foreseen accelerators. They must be searched for in the penetrating cosmic radiation, using large area detectors . GUT MMs should be characterized by relatively low velocities and relatively large energy losses. Direct searches for GUT MMs gave flux upper limits of few $`10^{16}`$ cm-2s-1sr-1 ; several indirect limits were obtained by different experiments .

The predicted spatial structure of a GUT MM is illustrated in Fig. 1 . The various regions are described in the figure caption. A GUT magnetic monopole may catalyse proton decay ; the cross section for this process could be relatively large because of the size of the fermion-antifermion condensate, which may contain terms violating baryon number conservation (see Fig. 1).

### 2.2 Intermediate mass monopoles

Magnetic monopoles with masses of $`10^{10}÷10^{12}`$ GeV/c2 are predicted by theories with an intermediate mass scale  and would appear in the early universe at a time considerably later than the GUT time. Also these MMs are topological point defects; an undesirable large number of relatively light monopoles may be gotten rid of by means of higher dimensional topological defects (strings, walls, textures) .

A possible structure of an intermediate mass MM is illustrated in Fig. 2. Notice that it has a larger core compared to the structure of a GUT monopole. The various regions correspond to: (i) intermediate mass scale of $`R10^{25}`$ cm; (ii) the electroweak scale; (iii) the condensate of fermion-antifermion pairs which could be the same as for a GUT monopole, but without any term violating baryon number conservation (thus these MMs cannot catalyse proton decay); (iv) the confinement region; (v) the outside region.

The number of intermediate mass monopoles could be considerably higher than that of GUT monopoles and galactic magnetic fields could accelerate them to high velocities. It has even be assumed that few of them could reach extremely high energies, interact in the upper earth atmosphere and lead to the highest energy cosmic ray showers .

## 3 Nuclearites

Strange Quark Matter (SQM) should consist of aggregates of $`u,d`$ and $`s`$ quarks in approximately equal proportions . The SQM is a colour singlet, thus it may have only integer electric charges. The overall neutrality of SQM is ensured by an electron cloud which surrounds it, forming a sort of atom. SQM may be the ground state of QCD.

The aggregates of $`u,d,s`$ quarks will be denoted with the terms strangelet, quark bag, SQM and nuclearite core; we shall use the word nuclearite for the core+electrons system.

Strangelets could have been produced shortly after the Big Bang and may have survived as remnants; they could also appear in violent astrophysical processes, such as neutron star collisions. Nuclearites should have a constant density , $`\rho _N=M_N/V_N3.510^{14}`$ g cm-3, somewhat larger than that of atomic nuclei.

They should be stable for all baryon numbers in the range between ordinary heavy nuclei and neutron stars ($`A10^{57}`$) . Nuclearites could contribute to the cold dark matter.

The structure of SQM can be described in terms of a bag model . In order to equilibrate the chemical potential of the quark species, SQM should have a number of $`s`$ quarks slightly lower than the number of $`u`$ or $`d`$ quarks . Thus the nuclearite core should have a positive electric charge which would be balanced by a number of electrons $`N_e\left(N_dN_s\right)/3`$, where $`N_d,N_s`$ and $`N_e`$ are the numbers of quarks $`d,s`$ and electrons, respectively, assuming $`N_d=N_u`$ .

In the following, the radii $`R_N`$ refer only to the nuclearite core (strangelet); the core+electronic cloud (nuclearite) system should have the constant radius of $`1`$ Å for $`R_N<1`$ Å $`=10^5`$ fm.

For $`R_N10^5`$ fm, all electrons must be inside the quark bag.

For $`10^4<R_N<10^5`$ fm, a fraction of the electrons are inside the quark bag, another fraction is external and gives the global dimension of $`10^8`$ cm $`=10^5`$ fm to the nuclearite core + electrons system. In this condition a nuclearite is similar to a Bohr atom.

Fig. 3 illustrates qualitatively the space distributions for the system of nuclearite core+electrons. Notice that for $`R_N>10^5`$ fm, we picture the nuclearite as a sort of Thomson atom.

For nuclearite masses larger than $`1.510^9`$ g $`10^{15}`$ GeV/c2, the relation between mass and radius should be (mass $``$ volume, $`M_NV_NR_N^3`$)

$$R_N=\left(\frac{3M_N}{4\pi \rho _N}\right)^{1/3}$$
(1)
For $`M_N=1.510^9`$ g, the nuclearite should have a radius

$$R_N=\left(\frac{3}{4\pi }\frac{M_N}{3.510^{14}}\right)^{1/3}=8.810^6M_N^{1/3}=8.810^6(1.510^9)^{1/3}10^8\text{cm}=10^5\text{fm}=1\text{Å}$$
(2)
For $`M_N=10^{18}`$ GeV/c2 one has $`R_N10^5`$ fm $`\left(1000\right)^{1/3}=10^6`$ fm.

Assuming that formula (1) is valid also for $`M_N<10^{15}`$ GeV/c2, one has, for the nuclearite core, the following values

$$M_N=10^{18}\text{GeV/c}^2R_N10^6\text{fm}$$

$$M_N=10^{15}\text{GeV/c}^2R_N10^5\text{fm}=1\text{Å}$$

$$M_N=10^9\text{GeV/c}^2R_N10^3\text{fm}$$

$$M_N=10^3\text{GeV/c}^2R_N10\text{fm}$$

Fig. 4 gives the dependence of the core radius $`R_N`$ and of the radius of the core+electrons system on the nuclearite mass $`M_N`$.

The expected relation between mass $`M_N`$ and charge $`Z_N`$ for nuggets of strange quark matter is shown in Fig. 5. The solid line is the curve computed in ref. , the dashed line is a rough interpolation taking into account that for $`M_N10^{15}`$ GeV/c2 all the electrons are inside the quark bag. Nuclearites with $`\beta 10^3`$ could be detected by scintillators and nuclear track detectors, independently of the charge of the nuclearite core.

A “curious” problem arose in the discussion of the nuclearite production in heavy ion high energy colliders: there was a fear that possibly produced SQM would grow in size and destroy the earth. This was proven to be inaccurate in ref. .

The present flux upper limits on nuclearites in the cosmic radiation are discussed in Section 5 and given in ref. .

## 4 Q-balls

Q-balls  are aggregates of squarks $`\stackrel{~}{q}`$, sleptons $`\stackrel{~}{l}`$ and Higgs fields . The scalar condensate inside a Q-ball core has a global baryon number $`Q`$ (may be also lepton number) and a specific energy much smaller than 1 GeV per baryon. We assume that the $`Q`$ numbers of quarks and squarks are equal to 1/3 ($`Q_q=Q_{\stackrel{~}{q}}=1/3`$) or 2/3 ($`Q_q=Q_{\stackrel{~}{q}}=2/3`$). Protons, neutrons and may be electrons can be absorbed in the condensate.

The vacuum expectation value inside a Q-ball core develops along “flat directions” of the potential . These flat directions are parametrized by combinations of squarks and sleptons that are electrically neutral (otherwise they would not be flat directions). Supposing that the different flavour squarks are not mass degenerate, their numbers inside the Q-ball core would not be equal for the same reason as in the nuclearite core case.

By assuming that the baryon number is packed inside a Q-ball core, one can get upper limits for the Q-ball quantum number $`Q`$ and for the Q-ball mass $`M_Q`$: $`Q10^{30}`$ and $`M_Q10^{25}`$ GeV/c2, respectively; Q-balls with $`M_Q<10^8`$ GeV/c2 are unstable .

In the early universe only neutral Q-balls were produced: SENS (Supersymmetric Electrically Neutral Solitons), which do not have a net electric charge, are generally massive and may catalyse proton decay. SENS may obtain an integer positive electric charge absorbing a proton in their interactions with matter; thus we may have SECS (Supersymmetric Electrically Charged Solitons), which have a core electric charge, have generally lower masses and the Coulomb barrier prevents the capture of nuclei. SECS have only integer charges because they are colour singlets. Some Q-balls which have sleptons in the condensate can also absorb electrons. The squarks $`\stackrel{~}{q}`$ inside the scalar potential bag have essentially zero mass.

The possible structures of SECS and SENS are shown in Fig. 6 a,b.

SENS may also interact with a proton of the interstellar medium, catalyse the proton decay leading to the emission of $`23`$ pions (or kaons) and transform quarks into squarks via the reaction $`qq\stackrel{~}{q}\stackrel{~}{q}`$. Thus some SENS may become SECS with a charge $`+1`$ emitting $`2\pi ^0`$ with total energy of about 1 GeV.

When a SENS enters the earth atmosphere, it could absorb, for example, a nucleus of nitrogen which gives it the positive charge of $`+7`$ (SECS with $`Z=+7`$). The next nucleus absorption is prevented by Coulomb repulsion. If the Q-ball can absorb electrons at the same rate as protons, the positive charge of the absorbed nucleus may be neutralized by the charge of absorbed electrons. The incoming SENS remain neutral most of the time. Electrons may be absorbed via the reaction $`u+ed+\nu _e`$. If, instead, the absorption of electrons is slow or impossible, the Q-ball carries a positive electric charge after the capture of the first nucleus in the atmosphere (SECS). Other SENS could “swallow” entire atoms (remaining SENS). If a SENS could absorb an electron, it would acquire a negative charge (SECS with $`Z=1`$). In the following we shall neglect the possibility that a neutral Q-ball (SENS) becomes charged (SECS) and viceversa.

The Q-ball mass $`M_Q`$, core size $`R_Q`$ and global quantum number $`Q`$ are related by the following relations , in theories with a bag potential $`V=\varphi \left(M_S^4\right)`$,

$$M_Q=\frac{4\pi \sqrt{2}}{3}M_SQ^{3/4}5924M_S\left(\text{TeV}\right)Q^{3/4}\left(\text{GeV}\right)$$
(3)
$$R_Q=\frac{1}{\sqrt{2}}M_S^1Q^{1/4}1.410^{17}M_S^1\left(\text{GeV}^1\right)Q^{1/4}\text{(cm)}$$
(4)
where the parameter $`M_S`$ is the energy scale of the SUSY breaking symmetry. In the following we shall use $`M_S=100`$ GeV/c2 and $`M_S=1000`$ GeV/c2.

From Eq. (3) we have

$$Q^{1/4}=\left(\frac{3}{4\pi \sqrt{2}}\frac{M_Q}{M_S}\right)^{1/3}$$
(5)
Placing Eq. (5) into Eq. (4), we have

$$R_Q=\frac{1}{\sqrt{2}}M_S^{4/3}\left(\frac{3M_Q}{4\pi \sqrt{2}}\right)^{1/3}$$
(6)
Note that, as for nuclearites, we have $`R_QM^{1/3}`$, once $`M_S`$ is fixed.

For $`M_S=100`$ GeV/c2, we have from Eq. (6):

$$R_Q=\frac{1}{\sqrt{2}}\left(\frac{3M_Q}{4\pi \sqrt{2}}\right)^{1/3}$$
$`\left(7\text{a}\right)`$
$$R_Q\left(\text{fm}\right)0.39M_Q\left(\text{GeV/c}^2\right)^{1/3}=0.39\left(M_Q0.197\right)^{1/3}\text{fm}=0.227M_Q^{1/3}\text{fm}$$
$`\left(7\text{b}\right)`$
The dependence of $`R_Q`$ from $`M_Q`$ is shown in Fig. 7 for $`M_S=100`$ GeV/c2 and 1000 GeV/c2.

The Q-balls have been considered as possible cold dark matter candidates; their core sizes should be only one order of magnitude larger than a typical atomic nucleus .

Flux limits on Q-balls come mainly from the astrophysical dark matter limit given in Fig. 8. SECS with $`\beta 10^3`$ and $`M_Q<10^{13}`$ GeV/c2 could reach an underground detector from above, SENS also from below . SENS may be detected by their almost continuous emission of charged pions (energy loss of about 100 GeV g-1cm2), while SECS may be detected by their large energy losses yielding light in scintillators, and possibly ionization. The energy losses of Q-balls in matter were computed in ref. . Flux upper limits on SECS could be deduced from the limits for dyons with the same electric charge; flux limits on SENS could be obtained from limits on MMs which catalyse proton decay.

## 5 Astrophysical limits

Magnetic monopoles, nuclearites and Q-balls could be components of the galactic cold dark matter, required by the rotation curves of the stars in the outskirts of our galaxy and of other galaxies. Assuming a local DM energy density of $`\rho =0.3`$ GeV/cm3 and that MMs, nuclearites and/or Q-balls could be part of it and have typical velocities of $`\beta 10^3`$, we can obtain upper limits on their flux in the cosmic radiation.

Fig. 8 shows the flux upper limits for nuclearites and for Q-balls versus their mass, assuming that they saturate the local DM density. Clearly, if their percentages are only $`1\%`$ of the local DM, then the limits have to be lowered by a factor of 100.

An experiment like MACRO has placed upper limits on the fluxes of MMs and nuclearites and could place similar upper limits for SECS at the level of few $`10^{16}`$ cm-2s-1sr-1 for masses larger than approximately $`10^{14}`$ GeV/c2 , see Fig. 7. The experiment should also be capable to place limits on the MM catalysis of proton decay.

The SLIM experiment could reach a level of sensitivity of few $`10^{15}`$ cm-2s-1sr-1 for masses larger than approximately $`10^8`$ GeV/c2 .

For superheavy MMs the DM limits are considerably larger than the present experimental limits. But monopoles can gain energy in the galactic field; from the survival of the galactic field one has more stringent limits at the level of $`10^{15}`$ cm-2s-1sr-1 (Parker limit ), or even orders of magnitude lower (Extended Parker limits ). The limits, both experimental and expected, are less well known for intermediate mass monopoles.

## 6 Conclusions

We have discussed qualitative pictures of the structures of superheavy and intermediate mass magnetic monopoles, of nuclearites and of Q-balls. In particular we have given the possible structure of these objects as function of their mass. We concluded with astrophysical considerations on their flux limits in the cosmic radiation. For nuclearites and Q-balls we assumed that they could be part of the cold DM and have $`\beta =10^3`$. Magnetic monopoles may be accelerated by galactic and intergalactic magnetic fields to higher velocities; for MMs, DM limits are considerably higher than limits based on the influence of magnetic fields.

Acknowledgments

We thank many members of the MACRO Collaboration for cooperation and for providing information and criticisism. We thank A. Kusenko for many stimulating discussions.


text:
# 1 Introduction

## 1 Introduction

The origin of the highest energy cosmic rays is a problem that has persisted for 4 decades since the pioneering measurements at Volcano Ranch . There is some consensus that, above the spectrum’s ankle at about $`10^{18.5}`$ eV, they originate outside the disk of the Galaxy. For particles of such high magnetic rigidity, sources in the Galaxy’s disk would presumably cause an obvious anisotropy in arrival directions that is not observed. Evidence for a composition changing to lighter particles at the ankle  strengthens this argument (particles having even greater rigidity because of lesser charge) and supports the view that cosmic rays with energies above the ankle are of extragalactic origin. The sources of those particles remain to be identified.

The observations of cosmic rays with energies above the expected GZK cutoff  should be a powerful clue to the nature of the sources. The Fly’s Eye  and AGASA  measured air showers with energies well above the GZK threshold. Recent reports  suggest that the spectrum might continue without a strong GZK effect. These super-GZK results have posed several related, but distinguishable, puzzles:

1. How are particles produced with such prodigious energy?
2. Why do the arrival directions of those particles not point back to recognizable sources in our local part of the universe?
3. Why is the intensity of particles above $`6\times 10^{19}`$ eV not more strongly suppressed?

New attempts to improve the observational data include the recently commissioned High Resolution Fly’s Eye (HiRes) , the Pierre Auger Observatory , the proposed Telescope Array Project , and Airwatch/OWL  plans for future space-based detectors of atmospheric air showers. The Auger Observatory and the potential space-based detectors will have exposure to the entire sky, which will open new possibilities for anisotropy analysis. These methods will be explored here in the context of the better controlled exposure of the Auger surface arrays.

Hillas  pointed out years ago that there are few astrophysical sites that can produce a large enough “electrical potential” $`\beta BL`$ to accelerate even highly charged nuclei to $`10^{20}`$ eV. (Here $`\beta =v/c`$ is the relative velocity of moving media with magnetic field strength $`B`$ and size $`L`$. In the case of statistical acceleration, including Fermi shock acceleration, the same product $`\beta BL`$ governs the maximum particle rigidity even though particles do not pass through a monotonic change in electrical potential.) Puzzle #1 is exacerbated in many contexts by synchrotron radiation and/or pion photoproduction.

Puzzle #2 can be resolved by invoking stronger-than-expected extragalactic magnetic fields, but that does not readily simplify puzzle #3. The expected suppression of particles above the GZK cutoff is based on travel time (or total distance traveled) rather than the straight-line distance to the sources. Particles below the GZK threshold have been accumulating over billions of years, whereas the mean age of particles well above the threshold cannot be greater than tens of millions of years .

One possible inference from the lack of an observable GZK spectral break is that sub-GZK particles may not be much older than 30 million years either, in which case the GZK effect would not significantly suppress the cosmic rays above the threshold relative to those below. It is not reasonable to suppose that the sources of high energy cosmic rays turned on so recently, so the young age of sub-GZK cosmic rays would require an intergalactic mechanism for dissipating their energy. Known mechanisms (e.g. nuclear interactions, synchrotron radiation, $`e^\pm `$ production, pion photoproduction via infrared and visible background photons, etc.) do not rob energy from sub-GZK particles rapidly enough. Perhaps high energy cosmic rays are attenuated through interactions with the unknown dark matter of the universe .

More conventional approaches to puzzle #3 are to defeat the GZK cutoff with a very hard extragalactic spectrum (e.g. from topological defect annihilation ) or to evade it by invoking neutrinos  or non-standard particles  that are immune to the microwave background radiation. Others  conclude that the sources must be localized to the Galaxy, but distributed in a halo large enough that galactic anisotropy has not become obvious.

The cumulative cosmic ray observations at this time are not sufficient to sort out the possibilities. AGASA and HiRes are currently building up the world’s total exposure at the highest energies. With better statistics and better measurements, the observations could soon lead to a breakthrough that identifies the sources of the highest energy cosmic rays. This same hope has been expressed for decades through the course of numerous experiments, however, and the puzzles have only become deeper mysteries. The answers may not come easily, and we should prepare the best possible analyses of the energy spectrum, particle mass distribution, and arrival directions.

Careful determinations of the energy spectrum and mass composition can be used to weed out classes of theories, but these tools are not likely to yield a clear signature for picking out a unique theory. A positive identification of the cosmic ray sources requires seeing their fingerprint in the sky. This may come in the form of arrival direction clusters  that identify discrete sources, or it may come as a large-scale celestial pattern that characterizes a particular class of potential sources. In the worst case, we might discover that the arrival directions are isotropic and the sources still elude positive identification. In that case, observers must strive for the best possible upper limit on anisotropy.

The Auger Project’s surface arrays will provide the best search for anisotropy fingerprints. Their combined exposure function on the celestial sphere will be unambiguous because they operate continuously and are not sensitive to atmospheric variability. Continuous operation means the celestial exposure function is uniform in right ascension. By having observatories appropriately located in both the southern and northern hemispheres, the exposure does not vary strongly with declination either.

The methods described in this paper are applicable to any observatory with full-sky coverage. They are not limited to the Auger Project, although the specific Auger site locations are used in the example simulations that are reported here.

Coverage of the full sky could be achieved piecemeal by combining results from different experiments. There is serious risk of spurious results from such meta-analyses, however, unless the exposures, energy resolutions, and detector systematics are perfectly understood and correctly incorporated in the analysis. The reliable approach is to use identical detectors in both hemispheres or the same (orbiting) detector for both hemispheres.

This paper seeks to evaluate the sensitivity of a full-sky observatory to large-scale anisotropy patterns and how that sensitivity depends on the number of arrival directions in a data set. Large-angle fingerprinting will be needed if there are many contributing sources or if the flux from each single source is diffused over a large solid angle due to magnetic deflection of the charged cosmic rays. If, instead, there are point sources to be detected, then the advantage of a full-sky observatory is in mapping the entire celestial sphere with comparable sensitivity in all regions.

Full-sky coverage is crucial for large-scale anisotropy analysis. It makes it possible to do integrals over the sky, so the powerful tools of multipole moments and angular power spectra are available. With full sky coverage, cosmic ray anisotropy analysis will be similar to gamma ray burst anisotropy analysis. The numbers of events will be comparable, the direction error boxes will be comparable, the exposure non-uniformities will be comparable, and in both cases events come from all parts of the sky. All of the techniques that were employed to search for anisotropy in the BATSE data  can be applied to a full-sky cosmic ray data set.

With a cosmic ray detector in only one hemisphere, there is a solid angle hole in the sky where the detector has zero exposure despite the Earth’s rotation. A zero-exposure hole makes it impossible to do integrals over the whole celestial sphere. No matter how many events the detector collects overall, it will never determine any multipole moment. A single-hemisphere detector can test hypotheses like, “Does the observed distribution match better what would be accepted from the clustering of radio galaxies toward the supergalactic plane or what would be accepted from an isotropic distribution?” It can also make qualified measurements like, “Assuming the anisotropy is a perfect quadrupole with axial symmetry, fit for the axis orientation that best explains the observed celestial distribution.”

The role of an observatory, however, should be to map the sky and make results available in a form which is readily usable without knowledge of the detector properties and which is independent of any theoretical hypothesis. Low-order multipole tensors (or spherical harmonic coefficients) can summarize the large-scale information. The angular power spectrum reveals if there is clumpiness on smaller scales. These results can be tabulated so that theorists can test arbitrary models quantitatively without privileged access to the data. With approximately uniform exposure, even eyeball inspection of arrival direction scatter plots can show large-scale patterns that are hidden when steep exposure gradients dominate the scatter plots.

While the role of an observatory should be to map the sky and determine the patterns without preconceived expectations, it is nevertheless worthwhile to consider what might be learned by measuring the low order multipoles or the angular power spectrum.

Monopole. There is no information about anisotropy patterns in the monopole scalar by itself. It is simply the sky integral of the cosmic ray intensity. That is information already present in the energy spectrum. A pure monopole intensity distribution is equivalent to isotropy. The strength of other multipoles relative to the monopole is a measure of anisotropy.

Dipole. A pure dipole distribution is not possible because the cosmic ray intensity cannot be negative in half of the sky. A “pure dipole deviation from isotropy” means a superposition of monopole and dipole, with the intensity everywhere $`0`$.

A predominantly dipole deviation from isotropy might be expected if the sources are distributed in a halo around our Galaxy, as has been suggested . In this case, there is a definite prediction that the dipole vector should point toward the galactic center.

An approximate dipole deviation from isotropy could be caused by a single strong source if magnetic diffusion or dispersion distributes those arrival directions over much of the sky. In general, a single source would produce higher-order moments as well.

A dipole moment is measurable in the microwave radiation due to Earth’s motion relative to the universal rest frame . If we are moving relative to the cosmic ray rest frame, a dipole moment should exist also in the cosmic ray intensity (the Compton-Getting effect). At lower energies, this may occur if the sun and Earth are moving relative to the galactic magnetic field or if the cosmic rays are not at rest with respect to the galactic field. For extragalactic cosmic rays, a Compton-Getting dipole is expected if the Galaxy is moving relative to the intergalactic field or if the cosmic rays themselves are streaming in intergalactic space. In any case, the expected velocities would be small ($`v/c\stackrel{<}{}10^3`$), and the Compton Getting anisotropy  ($`\frac{I_{max}I_{min}}{I_{max}+I_{min}}=(\gamma +2)(v/c)`$) should be $`\stackrel{<}{}0.005`$. (Here $`\gamma `$ is the differential spectral index, which is roughly 3.) An anisotropy of one-half percent would require high statistics for detection (see section 3).

A larger dipole anisotropy might be produced by a cosmic ray density gradient. If the magnetic field is disorganized, the gradient produces streaming by diffusion and the Compton-Getting dipole vector is parallel to the density gradient. If there is a regular magnetic field, however, the expected dipole vector $`\stackrel{}{D}`$ can be perpendicular to both the gradient and the field direction, $`\stackrel{}{D}\stackrel{}{}\rho \times \stackrel{}{B}`$. The direction of strongest intensity corresponds to the arrival direction of particles whose orbit centers are located in the direction of increased density.

Quadrupole. An equatorial excess in galactic coordinates or supergalactic coordinates would show up as a prominent quadrupole moment. A measurable quadrupole is expected in many scenarios of cosmic ray origins, and is perhaps to be regarded as the most likely result of a sensitive anisotropy search.

In general, a quadrupole tensor is characterized by 3 relative eigenvalues with associated orthogonal eigenvectors. In the case of axial symmetry, there is a single non-degenerate eigenvector that gives the symmetry axis. An axisymmetric “prolate” distribution would be hot spots at antipodal points of the sky, whereas an “oblate” distribution has the excess concentrated toward the equator that is perpendicular to the symmetry axis.

The axis of an oblate quadrupole distribution might differ from the galactic axis or the supergalactic axis if we are embedded in a magnetic field that systematically rotates the arrival directions.

The angular power spectrum. Spherical harmonic coefficients for a function on a sphere are the analogue of Fourier coefficients for a function on a plane. Variations on an angular scale of $`\theta `$ radians contribute amplitude in the $`\mathrm{}=1/\theta `$ modes just as variations of a plane function on a distance scale of $`\lambda `$ contribute amplitude to the Fourier coefficients with $`k=2\pi /\lambda `$.

For cosmic ray anisotropy, we might look for power in modes from $`\mathrm{}=1`$ (dipole) out to $`\mathrm{}60`$, higher order modes being irrelevant because the detector will smear out any true variations on scales that are smaller than its angular resolution. For charged cosmic rays, magnetic dispersion will presumably smear out any point source more than the detector’s resolution function. Even at the highest observed particle energies, there is unlikely to be any structure in the pattern of arrival directions over angles smaller than $`3^{}`$. The interesting angular power spectrum is therefore probably limited to $`\mathrm{}\stackrel{<}{}20`$.

## 2 Exposure

For a cosmic ray observatory, exposure is a function on the celestial sphere. Measured in units $`km^2yr`$, it gives the observatory’s time-integrated effective collecting area for a flux from each sky position. In this paper, the relative exposure $`\omega `$ is usually the function of interest. That will be a dimensionless function on the sphere whose maximum value is 1. In other words, $`\omega `$ at any point of the sky is a fraction between 0 and 1 given by the exposure at that point divided by the largest exposure on the sky.

In other contexts, the term “exposure” refers to the total exposure integrated over the celestial sphere. It then has units $`km^2sryr`$. For example, in determining the cosmic ray energy spectrum, one divides the number of cosmic rays observed in each energy bin by the total exposure for that energy. (In general, an observatory’s exposure is energy dependent.) If there were evidence that the energy spectrum were not uniform over the sky, then we would need to use the exposure’s dependence on celestial position to map the spectrum over the sky.

Since the spectrum is defined by the number of observed events divided by total exposure, one can use the measured spectrum to get the expected number of cosmic rays for any given total exposure. In the case of the Auger surface arrays, the continuous acceptance is approximately $`14,000km^2sr`$, independent of energy above $`10^{19}`$ eV. After operating for 5 years, they will have a total exposure of $`70,000km^2sryr`$. The integral cosmic ray intensity above $`10^{19}`$ eV is approximately $`0.5/(km^2sryr)`$, and it falls roughly like $`E^2`$ (perhaps less rapidly, but the energy dependence is not well determined above $`6\times 10^{19}`$ eV.) Using this simple $`E^2`$ dependence gives the following estimates for Auger cosmic ray counts after 5 years:

* 35,000 above $`10^{19}`$ eV (believed to be mostly extragalactic)
* 2200 above $`4\times 10^{19}`$ eV (compared to 47 in the AGASA cluster analysis)
* 350 above $`10^{20}`$ eV (above the GZK threshold region)
* 35 above $`3.2\times 10^{20}`$ eV (highest energy measured so far)

How any number of detected cosmic rays are distributed on the sky depends on both the true celestial anisotropy and the observatory’s relative exposure $`\omega `$.

The relative exposure can be calculated as follows for a detector at a single site with continuous operation. Full-time operation means that there is no exposure variation in sidereal time and therefore constant exposure in right ascension. Suppose the detector is at latitude $`a_0`$ and that it is fully efficient for particles arriving with zenith angles $`\theta `$ less than some maximum value $`\theta _m`$. (Full efficiency means the zenith angle acceptance depends on zenith angle only due to the reduction in the perpendicular area given by $`cos(\theta )`$.) This results in the following dependence on declination $`\delta `$:

$$\omega (\delta )cos(a_0)cos(\delta )sin(\alpha _m)+\alpha _msin(a_0)sin(\delta ),$$

where $`\alpha _m`$ is given by

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

and

$$\xi \frac{cos(\theta _m)sin(a_0)sin(\delta )}{cos(a_0)cos(\delta )}.$$

The upper left plot of figure 1 shows the resulting declination dependence for a site at $`a_0=35^{}`$ and another site at $`a_0=+39^{}`$, which are the latitudes for the two Auger observatories. The detectors are assumed to be fully efficient out to $`\theta _m=60^{}`$, and no arrival directions are counted from larger zenith angles. The combined exposure is also shown. The maximum is at the north pole direction, which is always detectable at the northern site, although the effective detector area is reduced by $`cos(51^{})`$ for flux arriving from the north pole direction.

The lower plots in figure 1 are scatter plots of the accepted cosmic rays for each site, where directions have been sampled from an isotropic distribution but accepted according to each detector’s exposure. (A sampled cosmic ray direction is accepted if a randomly sampled number between 0 and 1 is less than the relative exposure $`\omega `$ for that direction.) There are 10,000 accepted arrival directions in each of those two plots. Shown in the upper right plot is the superposition of all 20,000 events from the combined observatory. The combined distribution is not uniform, but has the modest declination dependence indicated in the upper left plot.

## 3 Dipole sensitivity

The objective here is to study the sensitivity of a full-sky observatory to a dipole deviation from isotropy. How well can the dipole be measured? How does that accuracy depend on the number of arrival directions in the data set? How does it depend on the amplitude of the dipole anisotropy?

For a dipole deviation from isotropy, the cosmic ray intensity varies over the sky as

$$I(\stackrel{}{u})=\frac{N}{4\pi }(1+\alpha \stackrel{}{D}\stackrel{}{u}).$$

Here $`\stackrel{}{u}`$ is a unit vector defining the celestial direction, $`\frac{N}{4\pi }`$ is the average intensity, $`\stackrel{}{D}`$ is the dipole direction unit vector, and $`\alpha `$ is its (non-negative) amplitude. In order for the cosmic ray intensity to be nowhere negative, $`\alpha `$ must lie in the range $`0\alpha 1`$. The amplitude $`\alpha `$ gives the customary measure of anisotropy amplitude: $`\alpha =(I_{max}I_{min})/(I_{max}+I_{min})`$.

The dipole can be recovered from the celestial intensity function by

$$\alpha \stackrel{}{D}=\frac{3}{N}I(\stackrel{}{u})\stackrel{}{u}𝑑\mathrm{\Omega }.$$

In our case, the observed intensity function consists of $`N`$ discrete arrival directions, each associated with a relative exposure $`\omega _i`$. The components of the dipole vector are then estimated by

$$\alpha D_a=\frac{3}{𝒩}\underset{i=1}{\overset{N}{}}\frac{1}{\omega _i}u_a^{(i)},$$

where $`u_a^{(i)}`$ denotes a component of the $`i`$th vector, and $`𝒩`$ is the simple sum of the $`N`$ weights $`\frac{1}{\omega _i}`$. (These dipole components are linear combinations of the three spherical harmonic coefficients with $`\mathrm{}=1`$.)

To test this method’s sensitivity to a dipole of amplitude $`\alpha `$ when there are N directions in the data set, one can produce an ensemble of artificial data sets of this type (with random dipole directions $`\stackrel{}{D}`$). For each data set, use the above formula to estimate the dipole vector, and record the difference of the estimated $`\alpha `$ from the input alpha and also the angle between the estimated direction and the input dipole direction. These error distributions describe the measurement accuracy. The RMS deviation from the true $`\alpha `$ is a single number to characterize the amplitude measurement accuracy, and the average space angle error summarizes the accuracy of determining the dipole direction.

This procedure can be repeated for different values of $`N`$ and different values of $`\alpha `$. For any pair ($`N,\alpha `$) the ensemble of simulation data sets yields the amplitude resolution and direction resolution as above.

To generate an individual simulation data set, one samples $`N`$ arrival directions on the celestial sphere. First, a direction is sampled from the assumed celestial distribution with dipole deviation from isotropy. Then the detector acceptance is applied by rejecting the sampled direction if a random number is greater than the relative exposure $`\omega `$ for that direction. This continues until the data set has $`N`$ arrival directions. The data set then reflects both the presumed celestial anisotropy and the detector’s non-uniform exposure.

These methods yield the results summarized in figure 2. The number of arrival directions was increased by factors of two: $`N`$ = 250, 500, 1000, 2000, 4000, 8000, 16000, 32000. For each $`N`$, amplitudes were studied at $`\alpha `$ = 0.1, 0.2,…, 1.0.

The upper left plot of figure 2 shows that the amplitude is determined to an accuracy of about 0.1 with 250 directions, improving to approximately 0.01 with 32,000 directions.

The upper right plot shows the dipole direction resolution as a function of the number of arrival directions. The mean error is less than 10 degrees for all cases (250 directions or more) if the amplitude is nearly 1, and it is less than 10 degrees regardless of the amplitude if the number of directions is 16,000 or more.

The lower left plot shows that the mean dipole direction error decreases as the dipole amplitude increases. A strong amplitude yields a good direction determination even for a small number of directions in the data set. For example, with only 250 cosmic ray arrival directions, the dipole direction is determined to better than 20 degrees if the amplitude exceeds 0.5. To get that same resolution with $`\alpha `$=0.1, you need a data set with more than 4000 directions.

For a fixed number of arrival directions, the RMS error in the amplitude has little dependence on the amplitude. That is to say, you can distinguish amplitudes 0.85 and 0.90 as well as you can distinguish 0.10 and 0.15. For the purpose of detecting an anisotropy (as opposed to measuring it), the relevant quantity is the amplitude divided by the RMS error, which is the number of sigmas deviation from isotropy. That quantity increases with the amplitude $`\alpha `$. It can also be expected to increase in proportion to $`\sqrt{N}`$ as the number of arrival directions increases. The lower right plot shows that

$$\frac{\alpha }{\mathrm{\Delta }\alpha }0.65\alpha \sqrt{N}.$$

For $`N`$=250, for example, the deviation from isotropy increases from 1-sigma to 10-sigma as $`\alpha `$ increases from 0.1 to 1.0. For $`N`$=8000, the range is from 6-sigma to 60-sigma. Etc.

To achieve a 5-sigma detection of a Compton-Getting anisotropy amplitude of 0.005 would require 2.4 million arrival directions. An anisotropy amplitude of 0.2 from a galactic halo distribution of sources, however, could be detected at the 5-sigma level with 1500 arrival directions.

## 4 Quadrupole sensitivity

A quadruople deviation from isotropy is characterized by an intensity function on the celestial sphere given by

$$I(\stackrel{}{u})=𝐐(\stackrel{}{u},\stackrel{}{u}),$$

where $`\stackrel{}{u}`$ is an arbitrary direction unit vector and $`𝐐`$ is a symmetric 2nd order tensor. Its trace gives the monopole moment. Its other 5 independent components in any coordinate basis are determined from the $`\mathrm{}=2`$ spherical harmonic coefficients $`a_{2m}`$.

Denoting by $`\xi _i`$ the three eigenvalues of $`𝐐`$ and the three (unit) eigenvectors by $`\stackrel{}{q_i}`$, the intensity function has the form

$$I(\stackrel{}{u})=\xi _1(\stackrel{}{q_1}\stackrel{}{u})^2+\xi _2(\stackrel{}{q_2}\stackrel{}{u})^2+\xi _3(\stackrel{}{q_3}\stackrel{}{u})^2.$$

To keep the number of studied variables manageable, consideration will be limited to axisymmetric oblate intensity functions, as might be expected from sources in the galactic disk or near the supergalactic plane. Let $`\stackrel{}{q}`$ denote the eigenvector that defines the symmetry axis, and let $`\xi `$ be the ratio of its eigenvalue to those in the symmetry plane, so the intensity function on the sphere is of the form

$$I(\stackrel{}{u})\xi (\stackrel{}{q}\stackrel{}{u})^2+(\stackrel{}{u}_{}\stackrel{}{u}_{}),$$

where $`\stackrel{}{u}_{}\stackrel{}{u}(\stackrel{}{u}\stackrel{}{q})\stackrel{}{q}`$ is the part of $`\stackrel{}{u}`$ perpendicular to the symmetry axis, and $`0\xi 1`$. The anisotropy amplitude $`\alpha `$ is related to $`\xi `$ by

$$\alpha \frac{I_{max}I_{min}}{I_{max}+I_{min}}=\frac{1\xi }{1+\xi }\xi =\frac{1\alpha }{1+\alpha }.$$

The objective here is to test how accurately the anisotropy amplitude $`\alpha `$ and the symmetry axis direction $`\stackrel{}{q}`$ can be determined from a data set of arrival directions. How does the accuracy depend on the number of directions $`N`$ and on the amplitude $`\alpha `$?

The method of investigation is the same as for the dipole sensitivity study. For each pair ($`N`$,$`\alpha `$), an ensemble of simulation data sets are produced, each with a randomly chosen direction for its symmetry axis. For each data set, the arrival directions are sampled from the relative intensity function (with quadrupole anisotropy), and each direction is accepted with probability equal to the relative exposure $`\omega `$ evaluated at that direction.

The anisotropy amplitude $`\alpha `$ and symmetry axis direction $`\stackrel{}{q}`$ are estimated for each simulation data set. The tensor with components

$$S_{ab}Iu_au_b𝑑\mathrm{\Omega }$$

($`u_a`$ denoting a component of $`\stackrel{}{u}`$) has the same eigenvectors as $`Q_{ab}`$ and the components $`S_{ab}`$ can be estimated by

$$S_{ab}=\frac{1}{𝒩}\underset{i=1}{\overset{N}{}}\frac{1}{\omega _i}u_a^{(i)}u_b^{(i)}.$$

The eigenvectors and eigenvalues of this symmetric matrix are then found. The symmetry axis is taken to be defined by the eigenvector with the smallest eigenvalue. Let $`\mathrm{\Delta }`$ be that smallest eigenvalue subtracted from the average of the other two (which should be equal, corresponding to directions in the symmetry plane). The eigenvalue $`\xi `$ of the intensity tensor $`𝐐`$ is given in terms of $`\mathrm{\Delta }`$ by

$$\xi =\frac{210\mathrm{\Delta }}{2+5\mathrm{\Delta }}.$$

Then the anisotropy amplitude $`\alpha `$ is gotten by $`\alpha =\frac{1\xi }{1+\xi }`$.

Results for ensembles with different ($`N,\alpha `$) values are presented in figure 3 in complete analogy with the dipole results presented in figure 2. The RMS error in the amplitude and the average space-angle error in the direction of the symmetry axis both decrease as $`N`$ increases. The symmetry axis is also seen to be better determined as the anisotropy amplitude increases for any fixed number of arrival directions $`N`$. The sensitivity for detecting anisotropy, as shown in the lower right plot, is given by

$$\frac{\alpha }{\mathrm{\Delta }\alpha }0.45\alpha \sqrt{N}.$$

With the definition of anisotropy amplitude $`\alpha \frac{I_{max}Imin}{I_{max}+I_{min}}`$, detecting a quadrupole anisotropy requires more data than for a dipole anisotropy of the same amplitude. Twice as many cosmic ray arrival directions ($`(.65/.45)^22`$) are needed for the same resolution.

## 5 Spherical harmonics

For any data set of arrival directions (with full-sky exposure), the anisotropy patterns can be fully characterized by the set of spherical harmonic coefficients $`a_\mathrm{}m`$, in terms of which the intensity function over the sphere is given by

$$I(\theta ,\varphi )=\underset{l=1}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}a_\mathrm{}mY_\mathrm{}m(\theta ,\varphi ).$$

The coefficients $`a_\mathrm{}m`$ are given by

$$a_\mathrm{}m=I(\theta ,\varphi )Y_\mathrm{}m(\theta ,\varphi )𝑑\mathrm{\Omega }.$$

Real-valued spherical harmonics are used in this paper, so the coefficients are real. The real-valued $`Y_\mathrm{}m`$ functions are obtained from the complex ones by substituting

$$e^{im\varphi }\{\begin{array}{cc}\sqrt{2}sin(m\varphi )\hfill & m<0\hfill \\ 1\hfill & m=0\hfill \\ \sqrt{2}cos(m\varphi )\hfill & m>0\hfill \end{array}$$

For a set of $`N`$ discrete arrival directions with non-uniform relative expousre $`\omega (\stackrel{}{u})`$, the estimate for $`a_\mathrm{}m`$ is given by

$$a_\mathrm{}m=\frac{1}{𝒩}\underset{i=1}{\overset{N}{}}\frac{1}{\omega _i}Y_\mathrm{}m(\stackrel{}{u}^{(i)})$$

where $`\omega _i`$ is the relative exposure at arrival direction $`\stackrel{}{u}^{(i)}`$ and $`𝒩`$ is the simple sum of the weights $`\frac{1}{\omega _i}`$.

The upper left plot in figure 4 shows a scatter plot of 2921 directions to extragalactic infrared sources with $`z<0.01`$, obtained from the NASA/IPAC Extragalactic Database (NED) . There are certainly selection effects in these directions, but they are used here only as an example of an anisotropic celestial distribution. The spherical harmonic coefficients $`a_\mathrm{}m`$ for this distribution are plotted to the right of that scatter plot. There are 440 coefficients plotted for $`1\mathrm{}20`$. Each set of $`2\mathrm{}+1`$ coefficients is plotted for each $`\mathrm{}`$ over an interval of 0.4 units on the abscissa. These $`a_\mathrm{}m`$ constitute a “fingerprint” of the anisotropy. They define a celestial intensity function that is a smoothed version of the scatter plot. The prominent coefficients $`a_{11}`$ and $`a_{22}`$ in this example result from the strong excess of the Virgo Cluster seen in the left central part of the scatter plot. Virgo is at declination 12.7 degrees and right ascension 187 degrees, and the $`a_\mathrm{}m`$ coefficients are derived here using that equatorial coordinate system (not the plotted supergalactic coordinate system).

To illustrate how well the $`a_\mathrm{}m`$ coefficients characterize the anisotropy, arrival directions can be sampled from the intensity function that they define. The lower left plot is a scatter plot with the same number of directions (2921) based on the relative intensity function

$$I=\underset{l=1}{\overset{20}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}a_\mathrm{}mY_\mathrm{}m.$$

The Auger exposure function has also been imposed. The lower left plot should not be identical to the upper left plot because of the Auger exposure simulation as well as the random sampling from the smoothed celestial anisotropy function. It clearly does have the same primary features, however.

The lower right plot in figure 4 summarizes how well the anisotropy is determined this way as a function of the number of arrival directions. The filled circles in that plot correspond to the displayed scatter plot with 2921 arrival directions. The $`a_\mathrm{}m`$ coefficients were derived from that simulation data set (using the relative exposure weights) and compared with those derived from the infrared source distribution. For each $`\mathrm{}`$, the RMS difference in the $`a_\mathrm{}m`$ values is plotted. One can see that the typical error in $`a_\mathrm{}m`$ is small compared to the significant coefficients in the upper right plot that characterize the anisotropy. The asterisks in the lower right plot are derived in the same way using a simulation with 1000 arrival directions. The open circles are the RMS coefficient differences resulting from a simulation with just 250 sampled arrival directions. The anisotropy fingerprint in this example is still measurable with 250 directions, although the RMS uncertainty in the $`a_\mathrm{}m`$ coefficients grows as $`1/\sqrt{N}`$ as $`N`$ decreases.

## 6 The angular power spectrum

The angular power spectrum is the average $`a_\mathrm{}m^2`$ as a function of $`\mathrm{}`$:

$$C(\mathrm{})=\frac{1}{2l+1}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}a_\mathrm{}m^2.$$

The power in mode $`\mathrm{}`$ is sensitive to variations over angular scales near $`1/\mathrm{}`$ radians. The angular power spectrum provides a quick and sensitive method to test for anisotropy and to determine its magnitude and characteristic angular scale(s).

As an example, consider the distribution of galaxies shown in the upper left plot of figure 5. These are all galaxies with redshift $`z<0.01`$ (also obtained from NED), and they are plotted in the supergalactic coordinate system. The Virgo cluster is the highest density region toward the left in this plot. There are 7321 galaxy positions plotted. The angular power spectrum (with no exposure correction) is shown out to $`\mathrm{}=20`$ in the upper right plot of that figure. There is excess power at all $`\mathrm{}`$-values, but especially for the dipole and quadrupole moments ($`\mathrm{}`$ = 1, 2) due to the high intensity from Virgo and other parts of the supergalactic plane.

The lower left plot in figure 5 indicates how sensitivity to the power spectrum is affected by the number of arrival directions (with non-uniform exposure as is expected for the Auger Observatory). Open circles in that plot are the power spectrum derived using a data set of 250 arrival directions. Those directions were obtained by randomly sampling from the 7321 galaxy directions (without replacement) and rejecting a sampled direction if a random number between 0 and 1 fell above the relative exposure $`\omega `$ evaluated at that direction. The 250 directions therefore represent a simulation Auger data set if each galaxy direction were to have equal probability of being a cosmic ray arrival direction. The open circles in the lower left plot are a decent approximation to the power spectrum (filled circles) of the “true” power spectrum defined by all 7321 directions with uniform exposure. The triangles in that plot represent the power spectrum for 1000 arrival directions sampled in the same way from the 7321, and the squares are obtained from 4000 sampled with the Auger exposure. It is clear that the approximation to the true power spectrum improves as the data set gets richer, but the gross information is already present with 250 arrival directions.

The lower right plot of figure 5 indicates how much power is expected due to fluctuations when directions are sampled from an isotropic intensity (and biased for the non-uniform expousre). The power is the same for all $`\mathrm{}`$-values and decreases like $`1/N`$ as the number of arrival directions increases. The power spectrum of the galaxy distribution is well above this noise level for all $`\mathrm{}`$-values for the cases $`N`$ = 4000 or 16,000. For 250 directions, only the first prominent harmonics in the lower left plot are clearly above the noise level indicated in the lower right plot.

Any class of candidate objects (e.g. active galaxies, or active galaxies with giant radio hot spots) has a celestial distribution that can be compared with a cosmic ray map when the whole sky has been surveyed with adequate sensitivity. Full information about the celestial distribution is provided by the set of coefficients $`a_\mathrm{}m`$. They can be tabulated out to $`\mathrm{}=20`$ in a list of 441 numbers (including the monopole). The angular power spectrum is a coordinate-independent gross summary of the features present in the celestial distribution. For example, you may learn from it that there is a large quadrupole moment, but you do not learn if the quadrupole has axial symmetry or the orientations of its principal axes. Full anisotropy information is given by the 441 $`a_\mathrm{}m`$ coefficients, not the 20 $`C(l)`$ powers.

The magnitude of the angular power $`C(l)`$ for larger $`\mathrm{}`$-values may contain useful information in the case that cosmic rays come from a limited number of discrete sources. The solid angle extent of the typical source affects the power at large values of $`\mathrm{}`$. Figure 6 displays an example in which there are 50 sources of equal flux with positions sampled randomly on the sky. Three different sky plots are shown, corresponding to different hypotheses about how much the arrival directions are dispersed from the source direction. In the upper left plot, sampled arrival directions are accepted only if they lie within $`10^{}`$ of one of the sources. In the upper right plot, they are accepted only within 5 degrees of a source, and only within 1.5 degrees in the lower left plot. The graph in the lower right shows the power spectra for the three different simulations. The power at low $`\mathrm{}`$-values is governed by the chance pattern in the distribution of source positions. For $`\mathrm{}\stackrel{>}{}10`$, however, the power clearly increases as the amount of source smearing decreases. The high end of the measurable angular power spectrum is sensitive to anisotropy structure on that finer scale.

## 7 Discussion

There is great advantage in a cosmic ray observatory having exposure to the entire celestial sphere, especially if the relative exposure is nearly uniform. In that case, scatter plots of arrival directions are immediately interpretable, and eyeball evaluations can readily identify discrete sources or large-scale patterns. Discrete sources will be identified with equal sensitivity anywhere in the sky. If no such sources are found, the flux upper limits will be uniform over the sky.

At the highest energies there is no proven anisotropy. Unlike the COBE anisotropy analysis, it is not necessary to subtract a large known dipole pattern and a myriad of uninteresting foreground sources. Any cosmic ray deviations from isotropy will be of immediate interest. The search for cosmic ray anisotropy is more similar to the case of gamma-ray bursts, where expectations and early claims of anisotropy were not supported by additional data.

The role of an observatory is to map the sky and make the results available to the scientific community. This is highly challenging for an observatory without full-sky coverage. Measurements in that case are made with different sensitivity in different parts of the sky, and nothing at all can be said about a large hole where the exposure is zero. Certainly it is not possible to perform the full-sky integrations that are required to measure the multipoles of the celestial cosmic ray intensity. In this paper, frequent use is made of the inverse of the relative exposure, $`1/\omega `$. Such methods obviously fail if the relative exposure anywhere becomes infinitesimal or zero.

To underscore the difficulty of anisotropy analysis without full-sky coverage, one can cite the work by Wdowczyk and Wolfendale . In that paper, the authors argue that the cosmic ray intensity measurements support a model of excess arrivals from equatorial galactic latitudes. The argument is based on the same data that two experimental groups had previously used in support of a gradient in galactic latitude that suggested an excess from southern latitudes relative to northern latitudes. In effect, because those northern detectors had poor exposure for southern galactic latitudes, Wdowczyk and Wolfendale were able to argue that the data supported a quadrupole distribution rather than a dipole distribution. Neither a dipole moment nor a quadrupole moment can be measured without full-sky coverage.

The techniques outlined in this paper pertain to any full-sky detector. Non-uniformity in celestial exposure is not hard to handle, provided it is well determined and there is adequate exposure to all parts of the sky. The true cosmic ray intensity is mapped with a sensitivity that depends primarily on the total number of detected arrival directions. This number $`N`$ is related to particle energy and observing time for any detector of known acceptance. This relationship for the Auger Observatory was given in the Introduction.

While complete information about anisotropy is encoded in the $`a_\mathrm{}m`$ coefficients (tied to some specified coordinate system), important gross properties of the anisotropy are characterized by the (coordinate independent) angular power spectrum $`C(l)`$. One can tell, for example, if there is a strong dipole or quadrupole moment. Such large-scale patterns are expected in many theories. It should be noted, however, that $`C(1)`$ gives the dipole moment but not its direction. It is obviously important whether the dipole points toward the galactic center, toward Virgo, toward Cen A, or in some unexpected direction. Similarly, all components of the quadrupole tensor are of interest, not just the average of their squares, $`C(2)`$. Nevertheless, the angular power spectrum provides a powerful tool for discriminating between viable and non-viable theories without detailed investigation. Also, the higher-order moments of the angular power spectrum can quantitatively characterize whatever clumpiness may exist in a map of arrival directions.

The techniques and examples mentioned in this paper are only representative of the powerful analysis methods that become possible with full-sky observatories. Data sets from such observatories will open a rich field of anisotropy study. The primary goal seen from the present time is the discovery of the highest energy cosmic ray origins. If that objective is accomplished (perhaps even before full-sky data sets are available), then the observed patterns from a known source distribution will be analyzed to infer properties of magnetic fields in the galactic halo and in intergalactic space.

A full-sky observatory has the ability to summarize completely the anisotropy information with the use of the spherical harmonic expansion coefficients $`a_\mathrm{}m`$. A table of coefficients (perhaps with multiple columns for different energy cuts) will provide the whole story. As has been done in this paper, those coefficients will be reliably corrected for the observatory’s non-uniform exposure. Detailed anisotropy analysis will no longer require privileged access to detector data. The published anisotropy fingerprint encoded in the $`a_\mathrm{}m`$ spherical harmonic coefficients can be matched against any theoretical suspect by any interested investigator.

Acknowledgement. It is a pleasure to thank Brian Fick for many helpful discussions about these matters.