The small scale clustering properties of ultra-high energy cosmic rays as a constraint on the particle charge and intervening magnetic fields
Patrick Younk
Fort Collins
CO 80523
United States
The observations of anisotropy above 57 EeV and of a suppression of flux above 40 EeV are supportive of two long held hypotheses: (1) the ultra-high energy cosmic ray (UHECR) sources are extragalactic, and (2) the Greisen-Zatsepin-Kuz'min effect is operational. If these hypotheses are true, we should expect that the flux of UHECR from the few brightest sources comprises a large fraction of the total flux. In fact, it is likely that groups of two or more events that are associated with a single source have already been observed. Such groups of events will cluster at an angular scale $\bar{s}$, the apparent angular size of the brightest UHECR sources. This implies that the value of $\bar{s}$ is embedded within the autocorrelation spectrum of event directions. An indication of $\bar{s}$ in a certain energy range can be used to infer limits on the particle charge and intervening magnetic fields, both of which are not well constrained. This is possible since $\bar{s}$ is similar in scale to the average magnetic deflection. In this work, we describe a method for obtaining an indication of $\bar{s}$ in the energy range $E > 57~\mbox{EeV}$. A metric $m$, which indicates whether structure in the arrival directions extends to small angular scales, is defined. Then, Monte Carlo simulations are used to estimate the expected range of $m$ for two different scenarios, $\bar{s} = 2.5^\circ$ and $\bar{s} = 10^\circ$. The results of the simulations are not strongly dependent on input parameters that are poorly constrained, such as the spatial distribution of sources. We show that $m$ differentiates between the two scenarios for a data set that contains 92 events above energy $E > 57~\mbox{EeV}$. Such a data set will be available from the Pierre Auger Observatory in late 2010 or early 2011. We discuss the implications of different values of $m$ being observed in the near future and the corresponding limits on the value of $\bar{s}$, the particle charge, and intervening magnetic fields.
SourceSize15.tex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EXAMPLE FILE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% FOR JHEPcls 3.1.0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\documentclass[published, 12pt]{JHEP3} % 10pt is ignored!

\JHEP{00(2007)000}

\JHEPspecialurl{http://jhep.sissa.it/JOURNAL/JHEP3.tar.gz}

\usepackage{epsfig,multicol,bbm}

%Begin special definitions for Instructions file
\newcommand{\ttbs}{\char'134}%\backslash for \tt
\newcommand\fverb{\setbox\fverbbox=\hbox\bgroup\verb}
\newcommand\fverbdo{\egroup\medskip\noindent%
			\fbox{\unhbox\fverbbox}\ }
\newcommand\fverbit{\egroup\item[\fbox{\unhbox\fverbbox}]}
\newbox\fverbbox
\newcommand{\jhepname}{JHEP}
%end

\title{The small scale clustering properties of ultra-high energy cosmic rays as a constraint on the particle charge and intervening magnetic fields}

\author{Patrick Younk\\
	Colorado State University, Fort Collins, CO 80523, United States\\
	E-mail: \email{Patrick.Younk@colostate.edu}}

\received{} 		%%
\revised{}
\accepted{}		%% These are for published papers.


\preprint{\hepth{9912999}}	% OR: \preprint{Aaaa/Mm/Yy\\Aaa-aa/Nnnnnn}
			  	% Use \hepth etc. also in bibliography.



\abstract{The observations of anisotropy above 57 EeV and of a suppression of flux above 40 EeV are supportive of two long held hypotheses: (1) the ultra-high energy cosmic ray (UHECR) sources are extragalactic, and (2) the Greisen-Zatsepin-Kuz'min effect is operational. 
If these hypotheses are true, we should expect that the flux of UHECR from the few brightest sources comprises a large fraction of the total flux.
In fact, it is likely that groups of two or more events that are associated with a single source have already been observed.
Such groups of events will cluster at an angular scale $\bar{s}$, the apparent angular size of the brightest UHECR sources. 
This implies that the value of $\bar{s}$ is embedded within the autocorrelation spectrum of event directions.
An indication of $\bar{s}$ in a certain energy range can be used to infer limits on the particle charge and intervening magnetic fields, both of which are not well constrained. 
This is possible since $\bar{s}$ is similar in scale to the average magnetic deflection.
In this work, we describe a method for obtaining an indication of $\bar{s}$ in the energy range $E > 57~\mbox{EeV}$. A metric $m$, which indicates whether structure in the arrival directions extends to small angular scales, is defined. 
Then, Monte Carlo simulations are used to estimate the expected range of $m$ for two different scenarios, $\bar{s} = 2.5^\circ$ and $\bar{s} = 10^\circ$. 
The results of the simulations are not strongly dependent on input parameters that are poorly constrained, such as the spatial distribution of sources.
We show that $m$ differentiates between the two scenarios for a data set that contains 92 events above energy $E > 57~\mbox{EeV}$.
Such a data set will be available from the Pierre Auger Observatory in late 2010 or early 2011.
We discuss the implications of different values of $m$ being observed in the near future and the corresponding limits on the value of $\bar{s}$, the particle charge, and intervening magnetic fields.}

\keywords{ultra-high energy cosmic rays, magnetic fields}

\dedicated{For Mary,\\ you still fill up my senses.}




\begin{document} 

%\maketitle  IS IGNORED %%%%%%%%%%%


\section{Introduction}

The latest results from the Pierre Auger Collaboration offer some important clues regarding the sources of 
cosmic rays (CR) with ultra-high energy. The observed flux suppression above approximately 40 EeV \cite{Pierre3, HiResSpec} and anisotropy above 57 EeV \cite{Pierre1, Pierre2} are consistent with two long held hypotheses: (1) the ultra-high energy CR sources are extragalactic, and (2) the Greisen-Zatsepin-Kuz'min (GZK) energy loss processes are occurring. If these hypotheses are true, the single brightest CR source in the sky is expected to contribute a large fraction of the total flux. In the energy range $E > E_{th} = 57~\mbox{EeV}$, the relative flux of the brightest CR source is expected to be near the 10\% level~\cite{Younk, YounkICRC}.

The Pierre Auger Observatory has measured approximately 60 CR over 57 EeV. Therefore, it is likely that the data set from this observatory already contains groups of CR events that were accelerated in a common source. The events in each of these groups will cluster on an angular scale $\bar{s}$, the apparent angular size of the brightest CR sources.
This implies that the value of $\bar{s}$ is embedded within the autocorrelation spectrum of the arrival directions, and that an indication of its value is currently possible.
An indication of $\bar{s}$ is an important step toward charge particle astronomy, i.e., the measurement of astronomical objects with CR messenger particles.

In this work, we define the angular size of a CR source as the angular extend over which CR from the source arrive.
The angular size $s$ is the radius that contains 95\% of the source flux. This definition assumes that the aspect ratio of the source is close to unity.  If the measured flux density (particles per solid angle) is described by a Gaussian distribution $dn/d\Omega = A \exp{(-x^2/2a^2)}$, where $x$ is the angular distance from the center of the distribution, then $\int_0^{2.0a} dn = 0.95$, and the angular size of the source is $s = 2.0 a$.
%
%\begin{equation}
%$$s = 2.0 a.$$
%\label{eq:SourceSize1}
%\end{equation}
%

A source that appears point-like, $s \approx 0$, when imaged with neutral particles (i.e., conventional photonic astronomy) will generally appear as an extended object when imaged with charged particles, $s >> 0$. 
The difference has nothing to do with the size of the emitting region. Rather, it is due to how CR are deflected by magnetic fields en route to earth. 
For a given energy range the value of $\bar{s}$ contains information on the CR charge and intervening magnetic fields.
In general, $\bar{s}$ is similar in scale to the average magnetic deflection of CR from the brightest sources.
This suggests the following analysis chain, which we detail in this paper. The observed small scale clustering properties of ultra-high energy CR are used to constrain the physical parameter $\bar{s}$. Then $\bar{s}$ is used to constrain the parameter space that describes the allowed distribution of particle charge and the allowed structure of the intervening magnetic fields.

This paper is organized as follows.
In section 2, we discuss the expected range of $\bar{s}$ and how its value depends on the particle charge and the intervening magnetic fields. Then in section 3, we 
predict, via Monte Carlo simulation, the small scale clustering properties of ultra-high energy CR under two different plausible assumptions for the value of $\bar{s}$. The small scale clustering properties are quantified with a metric $m$.
%We show that $m$ is a strong predictor of the value of $\bar{s}$ even for data sets that contain less than 100 events. 
In section 4, we discuss how different observed values of $m$ constrain the value of $\bar{s}$ and, in turn, the particle charge and intervening magnetic fields.


\section{Source size as a function of particle charge and intervening magnetic fields}

Throughout this work, we only consider CR observed in the energy range $E > E_{th} = 57~\mbox{EeV}$.  In this energy range, the spectral slope of the diffuse background is $\gamma \approx 5$ \cite{Pierre3, HiResSpec}. If the observed steepening of the energy spectrum beyond approximately 40 EeV is caused by the GZK effect, then we should expect relatively nearby sources to have a flatter spectrum than the diffuse background. For this reason, we assume the energy spectrum of the brightest sources is $dn/dE \propto E^{-\gamma}$ with $\gamma = 3$. We treat all CR as having the same charge $Z$.

Many details of the intervening magnetic fields are not well known. 
For the treatment here, we use a model that contains 
only the most well established structural details. (The model is necessarily simplistic.)
From this, we obtain an expected value of $\bar{s}$.
Then, we discuss the uncertainties of our simple model and how they affect the expected value of $\bar{s}$.


\subsection{A simple model for the magnetic deflection of CR}

We start by assuming that deflections in the galactic magnetic field are dominant over deflections in the extragalactic magnetic fields. The galactic magnetic field is commonly described as the superposition of a regular component and a turbulent component. 
The regular component is known to follow the spiral structure of the galaxy and to have a magnitude of $B \approx 2~\mu\mbox{G}$ near the solar system \cite{Han, Men, Sun, Beck}. For ultra-high energy CR arriving at earth, the average path length $L$ through the field is on the order of a few kiloparsecs. The value of $L$ depends on the galactic latitude of the source and extent of the magnetic field out of the galactic plane. A CR in the regular field is deflected by an angle
%
$$ \delta = 3.2^\circ \frac{10^{20}~\mbox{EeV}}{E/Z} \left| \int_0^L \left( \frac{\mbox{d}\textbf{x}}{3~\mbox{kpc}} \times \frac{\textbf{B}}{2~\mu\mbox{G}}  \right) \right|. $$


For the regular field, the correlation length $G_r \approx L$. Therefore, it is reasonable to approximate the regular field as perfectly uniform.
For the case that there is no turbulent component, the image of a point source is a line segment on the celestial sphere with one end at the source location. The distribution of arrival directions along the line segment as a function of distance from the source is given by
%
%\begin{equation}
$$\delta = 3.2^\circ \frac{10^{20}~\mbox{eV}}{E/Z} \frac{L}{3~\mbox{kpc}} \frac{B}{2~\mu\mbox{G}}.$$
%\label{eq:RegularField}
%\end{equation}
%


The turbulent field is thought to originate from supernova and their remnants, which transfer energy to the interstellar medium on scales of several tens of parsecs. The characteristic correlation length of the turbulent field is $L_c \approx 50~\mbox{pc}$, and its rms magnitude $B_{rms}$ is approximately twice that of the regular field \cite{Beck, Ohno}. The rms deflection for a CR in a turbulent field is
%
%\begin{equation}
$$\delta_{rms} = 0.6^\circ \frac{10^{20}~\mbox{eV}}{E/Z} \frac{B_{rms}}{4~\mu\mbox{G}} \sqrt{\frac{L}{3~\mbox{kpc}}} \sqrt{\frac{G_t}{50~\mbox{pc}}},$$
%\label{eq:TurbField}
%\end{equation}
%
where $G_t$ is the characteristic correlation length of the turbulent component.

For CR with low enough magnetic rigidity $E/Z$ such that $G/L < \delta ~\mbox{radians}$, many independent paths through the magnetized region will connect a source and an observer. This is known as magnetic lensing \cite{Harari2}. For our model, 
the lensing threshold for the regular component is $E/Z \approx 10~\mbox{EeV}$, while for the turbulent component it is $E/Z \approx 100~\mbox{EeV}$. 
From our starting assumptions, the minimum magnetic rigidity of the CR is $E/Z = 57~\mbox{EeV}$. 
This magnetic rigidity is within the lensing region for the turbulent field, but not the regular field. Because the assumed particle energy distribution is a steeply falling power-law with nearly all the CR with $E/Z \approx 70~\mbox{EeV}$, we treat all CR as being lensed by the turbulent component.

In our model $G_t << L$, which implies that lensing by the turbulent component can be characterized as a diffusive process.
We approximate the effects of the turbulent component with a Gaussian point spread function $\exp(x^2/2\delta_{rms}^2)$, where $x$ is the angular distance from the source direction.
For the case that there is no regular component, the image of a point source is a circular cloud centered on the source location.


If we consider the superposition of the regular and turbulent components, we expect that a CR source will appear like 
an oblong cloud. In Fig. 1, we show the expected arrival direction distribution for our model with $B = 2 \mu\mbox{G}$, $B_{rms} = 4 \mu\mbox{G}$, $L = 3~\mbox{kpc}$, $G_t = 50~\mbox{pc}$, and $Z=1$.

\DOUBLEFIGURE[t]{PSize1.eps, width=.5\textwidth}
{PSize2.eps, width=.5\textwidth}{The expected flux density distribution for an ultra-high energy CR source observed in the energy range $E>57~\mbox{EeV}$. $Z=1$ and $dn/dE \propto E^{-3}$ is assumed.}
{Relative flux as a function of distance from the center of the flux density distribution shown in Fig. 1.}


The shading represents relative intensity in five equally spaced logarithmic intervals.  The mean of the distribution is shown with a cross. The source location is shown with an x. The angular separation between the source location and the mean of the distribution is $3.8^\circ$. 
In Fig. 2, we show the distribution of arrival directions as a function of distance from the mean.
95\% of the events are within $2.5^\circ$ of the mean of the distribution.
The dashed circle in Fig. 1 is concentric with the location of the mean and has a radius of $2.5^\circ$.  

Based on our simple model and assuming $Z=1$, it is plausible that $\bar{s} = 2.5^\circ$. This is a factor of 3 or 4 greater than the typical angular resolution of an ultra-high energy CR observatory \cite{AugerAngleRes, HiResAngleRes}. In the range that our magnetic lensing assumptions are valid, a good approximation is $\bar{s} \propto BB_{rms}LZ$


\subsection{The plausible range of $\bar{s}$}

The estimate of $\bar{s}$ given in the last section is relatively uncertain. The uncertainty comes from several areas.
First, uncertainties in the structure of the galactic magnetic field (e.g., Ref. \cite{Men}) and the latitude of the brightest sources amount to uncertainties on $\bar{s}$ of about a factor of 2. 

Second, the extragalactic magnetic fields are not well constrained. It has been shown \cite{Ryu} that filaments associated with the large scale structure may harbor nG magnetic fields with correlation lengths on the order of 1 Mpc. These structures have characteristic sizes measured in Mpc and would effectively scatter 60 EeV protons. CR from a source within or behind a filament may be diffused over several degrees on the sky.

Third, the charge of the particle is unknown. Presently, this consideration is particularly interesting.
The latest composition results from the Pierre Auger Observatory \cite{Bellido}
may indicate that heavy nuclei such as iron ($Z=26$) are common at the highest energies.
This indication is not certain because it is not currently possible to decouple ultra-high energy composition measurements from the phenomenology of high energy particle physics. The compatibility of iron with the observed anisotropy \cite{Pierre1, Pierre2} seems to be only possible if there are few bright sources at high galactic latitude, the galactic magnetic field is confined to the thin disk, and the extragalactic magnetic fields are negligible. Even if this is the case, we must expect $\bar{s} > 10^\circ$. If the particles are not heavy nuclei, it could be that extensive air showers initiated by ultra-high energy protons develop more like what it expected of  showers initiated by iron. Either result is intriguing.

The above considerations indicate that the true value of $\bar{s}$ may be near $1^\circ$ or many tens of degrees.
The value of $\bar{s}$ is very much an open question. 



\section{Small scale clustering as a function of source size}

Our method uses a small scale clustering metric $m$ to differentiate between different values of $\bar{s}$.
We use a Monte Carlo algorithm to predict the range of $m$ for two source size scenarios: (1) Small source, $\bar{s} = 2.5^\circ$, and (2) Large source, $\bar{s} = 10^\circ$. The small source scenario represents a situation where $Z=1$, the galactic magnetic fields are mainly confined to the thin disk, and the extragalactic fields are negligible. The large source scenario represents a situation with $Z>1$ and/or the magnetic deflections in the galactic halo and the intergalactic medium are significant. 

The Monte Carlo algorithm generates sets of events with an assumed value for $\bar{s}$ and calculates an expected range of $m$.
By comparing the value of $m$ from an observed data set to the expected ranges of $m$ obtained from the Monte Carlo calculations, we can determine which source size scenario is favored.
The small scale clustering metric $m$, the models used to generate the simulated event sets, the Monte Carlo algorithm, and the results of the algorithm are described in the following sections.

\subsection{The Small scale clustering metric}

To indicate the level of clustering at an angular scale $\chi$, we calculate the number of pairs with angular separation less than $\chi$ weighted by the inverse of their angular separation and by $1/\chi$
%
$$ M(\chi) = \frac{1}{\chi} \sum^{n}_{i=2} \sum^{i-1}_{j=1} \Theta(\chi - \beta_{ij})/\beta_{ij}, $$
%
where $\beta_{ij}$ is the angular separation between event $i$ and $j$, $\Theta$ is the step function, and $n$ is the number of events. The motivation to weight each pair by the inverse of its angular separation comes from the fact that for an isotropic distribution of events, the expected number of event pairs with a angular separation $\beta$ is $\left\langle dn_p/d\beta \right\rangle \propto \beta$. This is valid for small $\beta$. We weight $M$ by $1/\chi$ so that $\left\langle M(x) \right\rangle \approx \left\langle M(y) \right\rangle$ for an isotropic distribution of events, where $x$ and $y$ are any  small angle.


The metric we use to differentiate between the small and large source size scenarios is then
%
\begin{equation}
m = M(2.5^\circ)/M(10^\circ).
\label{eq:ClusteringMetric}
\end{equation}
%
If $\bar{s} \geq 10^\circ$, we expect $m \approx 1$. If $\bar{s} < 10^\circ$, we expect $m>1$. 
If $\bar{s} < 2.5^\circ$, we expect $m \approx 10^\circ/2.5^\circ = 4$. In this way, $m$ indicates whether structure in the arrival directions of cosmic rays extends to small angular scales.

If $n$ is too small, there is a possibility of $m = 0/0$, which is undefined. We work with events sets where $n$ is large enough for this not to be a concern.
 

\subsection{The model for the distribution of sources}

We simulate $N = \rho/V$ sources, where $\rho$ is the number density of sources and $V$ is the volume of space considered in the simulation. We assume the sources are associated with large galaxies, galaxies near the knee of the luminosity function, yet not all large galaxies need contain an active source. This choice is motivated by the fact that most stars and super massive black holes (i.e., potential energy sources for particle acceleration) are contained within large galaxies.  

The distribution of nearby sources is constructed from the PSCz Catalog \cite{Sanders}.
The PSCz catalog contains over 15,000 galaxies measured by the Infrared Astronomical Satellite (IRAS) across 84\% of the sky.
To translate redshift $z$ (radial velocity) into distance, we use Hubbles's law $d = cz/H_0$ where $c$ is the speed of light and $H_0 = 71~\mbox{(km/s)/Mpc}$. The PSCz catalog is flux limited, not volume or luminosity limited. Therefore, the number density of sources decreases with distance. The expected number density of a source catalog in the absence of clustering is termed the selection function $\phi(d)$. For the PSCz catalog, the selection function  is well determined \cite{Sanders}. To form a distribution with an average source density $\rho$, we weight each galaxy in the catalog by $\rho / \phi(d)$. 

In Fig.~\ref{fig:PSC1}, we show with a dotted line the number density relative to $\rho$ after this weighting is applied.
There are regions of over-density and under-density, signifying the presence of clusters and voids. 
In particular, the local density of galaxies is greater than the average, signifying the Milky Way's association with the Virgo super cluster and the local group. At large distances, a large volume of space is averaged over so the number density of galaxies is necessarily near the average.


\FIGURE{\epsfig{file=fig200.eps,width=.8\textwidth} 
        \caption{Source number density relative to the average source density $\rho$ vs. distance. The dotted line is the PSCz catalog weighted by $\rho/\phi$, where $\phi$ is the selection function. The solid line is the PSCz weighted by $\psi \rho/\phi$, where $\psi$ is a function designed to reproduce the observed local over-abundance of large galaxies.}
	      \label{fig:PSC1}}

The local over-density of potential sources is an important input parameter to our simulations if the actual number density of sources is relatively high (i.e., $\rho \geq 10^{-3}~\mbox{Mpc}^{-3}$). 
For example, a high local over-density means a greater probability for a few nearby (and thereby bright) sources.
These nearby sources may contribute a significant fraction of the flux even though the total number of sources within the GZK sphere may be relatively large. 

The distribution of nearby large galaxies is known quite well. Accurate distance estimates are available that do not depend on recessional velocity, and so are not affected by local fluctuations from the average Hubble flow.
From Ref. \cite{Karachentsev2} there are 30 large galaxies (i.e., with an absolute blue magnitude $M_B < -19.5$) that are within 8 Mpc. The density of galaxies with $M_B < -19.5$ averaged over a large volume is approximately $6 \times 10^{-3}~\mbox{Mpc}^{-3}$ (e.g., see Ref. \cite{Marzke}). This implies the number of large galaxies within 8 Mpc is 2 times greater than average. This is in good agreement with the PSCz catalog weighted by $\rho / \phi(d)$.

However, the number density of large galaxies at smaller radii is underestimated by this model.
For example, there are nine large galaxies that are within 4 Mpc \cite{Karachentsev2}. This is approximately 6 times greater than the global average. Furthermore, the number of large galaxies within 2 Mpc is 10 times greater than average.
We list the details of large galaxies that are within 4 Mpc in Table~\ref{table1}. 

\TABULAR{c c c c c c c}{
	\hline	 Name	   & Type & $M_B$  & Declination & $D$   & Distance  & Nuclear \\
	                 &      & (mag)  &             & (Mpc) & Method    & Activity \\ \hline
	\hline	 Milky way & Sbc     & -20.80 & $-29^\circ$  & 0.01 & cep & - \\
		       M31      & Sb      & -21.58 & $+41^\circ$  & 0.77 & cep & Seyfert2 \\
		       M81      & Sb      & -21.06 & $+69^\circ$  & 3.63 & cep & LINER \\
		       M82      & Sdm     & -19.63 & $+69^\circ$  & 3.53 & rgb & Starburst \\
		       NGC5128  & -S$0^0$ & -20.77 & $-43^\circ$  & 3.66 & rgb & BL Lac?\\
		       NGC4945  & Scd     & -20.51 & $-49^\circ$  & 3.6  & mem & Seyfert \\
		       IC342    & Sc      & -20.69 & $+68^\circ$  & 3.28 & cep & Starburst \\
		       Maffei2  & Sbc     & -20.15 & $+60^\circ$  & 2.8  & tf  & Starburst \\
		       NGC253   & Sc      & -21.37 & $-25^\circ$  & 3.94 & rgb & Seyfert \\
		      
	\hline}
	{\label{table1}List of large galaxies within 4 Mpc. Absolute blue magnitude and distance are taken from Ref.~\cite{Karachentsev2}. Here cep = distance from the luminosity of Cepheids, rgb = distance from the luminosity of the tip of the red giant branch, men = distance from group membership, tf = distance from the Tully-Fisher relation.}


If we weight the PSCz only by $\rho / \phi(d)$, the local over-density of large galaxies is under estimated by a factor of about 2.5. This is because the PSCz samples mainly small (low luminosity) galaxies at small distances and the local region has a higher ratio of large galaxies to small galaxies than the average. To include this observation into our model in a simple manner, we weight the PSCz catalog by an addition factor
$$ \psi =
\cases{ 
2.5       & : $d < 4~\mbox{Mpc}$ \cr
4 -0.375d & : $4 < d < 8~\mbox{Mpc.}$ \cr
1.0       & : $d > 8~\mbox{Mpc.}$ \cr}
\label{eq:psi} $$
%
In Fig.~\ref{fig:PSC1}, we show with a solid line the number density for the PSCz catalog weighted by $\psi \rho/\phi$. 

Beyond 60 Mpc, the number density of the PSCz catalog rapidly drops below that of large galaxies. Our model for the distribution of sources, therefore, follows the PSCz catalog weighted by $\psi \rho/\phi$ out to a distance of 60 Mpc. For distances greater than 60 Mpc, we assume the sources are evenly distributed (i.e., every location has equal probability of containing a source). This is a valid approximation since the characteristic size of super clusters is a few tens of Mpc.

The validity of this model rests on the assumption that the sources are distributed in proportion to the number density of large galaxies. Although this assumption is well motivated, it is not hard to imagine other scenarios. For example, if the sources are hyper-nova, we should expect the sources to be distributed in proportion to the massive star formation rate. It is not clear how closely correlated this is with the number density of large galaxies. To test the sensitivity of our results to the source distribution model, we also test a model where the sources are distributed evenly throughout the whole test volume.


\subsection{Other models used in the Monte Carlo algorithm}

To calculate energy losses associated with CR propagation, we assume that the particles are protons. We take into account energy losses due to inelastic interactions with background radiation fields (the GZK effect). We do not consider 
energy losses due to the expansion of the universe. For the source distances in our simulation, these redshift losses are negligible. 
 
We model the flux from a source as $q = A L/4 \pi d^2$, where $A$ is the GZK attenuation factor, $L$ is the source luminosity, and $d$ is the source distance. The GZK attenuation factor is calculated, as in Ref.~\cite{Harari}, as the fraction of particles originally above $E_{th}$ that still have an energy above $E_{th}$ after traversing a distance $x$. For a source with a proton injection spectrum $dn/dE \propto E^{-\alpha}$, the attenuation factor is $A(E_{th},x) = ((E_{th} + \Delta E)/E_{th})^{1-\alpha}$, where $\Delta E$ is the average energy loss of a proton traversing a distance $x$. The energy loss $\Delta E$ is calculated by solving $dE/dx = -E/ \lambda$, where $\lambda$ is the proton energy loss length in the extragalactic medium (e.g., see \cite{Protheroe}). This is the so-called continuous energy loss approximation. As suggested by Ref.~\cite{Allard}, we assume $\alpha = 2.6$.

We describe the source luminosity distribution with a differential luminosity function $\phi = dN/dL$. 
The luminosity function is totally unconstrained. Therefore, we consider the simplest case, that all the sources have the same luminosity; i.e., $\phi = \delta(L_0 - L)$. This is not a realistic luminosity function, but it represents an important limit. It is the luminosity function that proportions the flux the most evenly among the sources.
In our simulations, a very similar effect occurs when the source number density is increased.
Since we test a wide range of source number densities, there is no need to test a wide range of luminosity function models.

The angular size of a source is modeled as a 2-D Gaussian.  
That is, the expected event distribution from a source is $dn/d\Omega \propto \exp{(-2x^2/\bar{s}^2)}$, where $x$ is the angular distance from the source direction.
We assume that every source has the same characteristic size.
We also assume that the centroid of the flux density distribution is centered at the source location.
In reality, the centroid would be offset from the source location because of the regular component of the galactic magnetic field. This approximation does not introduce a bias since we only require the PSCz catalog to model the general clustering behavior of the sources, not the absolute locations of the sources.

The expected event distribution is weighted by the sky exposure for a particular CR observatory.
The details of the sky exposure have little impact on our results. Here we suppose the most general case, an observatory with equal exposure to all parts of the sky.



\subsection{The Monte Carlo algorithm}

The number density of ultra-high energy sources is not known. The lack of an apparent source in our galaxy implies $\rho < 6 \times 10^{-3}~\mbox{Mpc}^{-3}$, the approximate number density of large galaxies. The existence of CR events with energy greater than $10^{20}~\mbox{EeV}$ and the GZK effect imply $\rho \geq 10^{-6}~\mbox{Mpc}^{-3}$, the number density where we expect at least a few sources within 100 Mpc. We test four values of $\rho$, logarithmically spaced within this range of density.

Our test volume is a sphere centered at earth with radius $D$, the distance that nearly all CR with energy above $E_{th}$
are expected to originate within. The value of $D$ is calculated as in Ref.~\cite{Younk, Harari}. For 
$E_{th} = 57~\mbox{EeV}$, we use $D = 250~\mbox{Mpc}$. We have checked that increasing $D$ does not change the results of the Monte Carlo algorithm. 

To generate a single event data set, we randomly disperse $N = \rho 4 \pi D^{3}/3$ sources throughout the test volume according to one of the distribution models described in section 3.2. We calculate the expected distribution of events from each source using either the small source size scenario or the large source size scenario. These distributions are weighted by the source flux and sky exposure.
The expected distribution of events from all sources is then the sum of the individual distributions from each source.
From this distribution, we randomly generate $n = 92$ events. This value of $n$ corresponds to the number of events with $E>57~\mbox{EeV}$ expected to be measured by the Pierre Auger Observatory at its fully deployed southern site over a four year time span. This number of events will probably be reached in late 2010 or early 2011.
In Fig.~\ref{fig:SampleDist1}, we show a Monte Carlo realization for the small source scenario, and 
in Fig.~\ref{fig:SampleDist2}, we show a Monte Carlo realization for the large source scenario.


\FIGURE{\epsfig{file=plotAAA1.eps,width=0.8\textwidth} 
        \caption{The expected arrival direction distribution for a single Monte Carlo realization assuming that the CR sources are associated with large galaxies and that $\bar{s} = 2.5^\circ$ and $\rho = 10^{-4}~\mbox{Mpc}^{-3}$. The figure is in galactic coordinates. The blue shading represents the relative probability of a CR arriving from a certain sky direction, darker corresponding to a higher probability. There are 5 levels of shading, equally spaced on a logarithmic scale. The dotted line marks the location of the super galactic plane. The red crosses are 92 event directions, randomly pulled from this distribution. For this set of events directions, $m = 2.9$.}
	      \label{fig:SampleDist1}}

\FIGURE{\epsfig{file=plotAAA2.eps,width=0.8\textwidth} 
        \caption{The expected arrival direction distribution for a single Monte Carlo realization with $\bar{s} = 10^\circ$. All other conditions are the same as Fig.~\ref{fig:SampleDist1}. For this set of events directions, $m = 1.6$.}
      	\label{fig:SampleDist2}}
 
 
After a simulated event set is generated, the value of $m$ is calculated with Eq.~\ref{eq:ClusteringMetric}. 
We calculate the expected range of $m$ by running 1000 Monte Carlo simulations.
This process is repeated for different values of $\rho$ and $\bar{s}$.


\subsection{Results}

In Fig.~\ref{fig:Results1}, we show the expected values of $m$ as a function of source density.
We show the results for the two source size scenarios and the two source distribution models.
The error bars represent the 10-90\% quantiles for the distribution of $m$. 

The expected value of $m$ for the small source scenario is near 3, while the expected value of $m$ for the large source scenario is near 1. 
The only situation where the ranges of $m$ for the two source size scenarios
substantially overlap is where the sources are modeled as being evenly distributed with a number density of $10^{-3}~\mbox{Mpc}^{-3}$.


\FIGURE{\epsfig{file=fig800.eps,width=.8\textwidth} 
        \caption{The expected values of $m$ as a function of source number density. The solid markers are the results using the large galaxy source model. The open markers are the results using the evenly distributed source model. The markers are slightly offset from each other on the x-axis for clarity.
The error ranges represent the 10-90\% quantile range.}
	      \label{fig:Results1}}


\section{Discussion and Conclusions}

The expected value of $m$ is not a strong function of source number density. This is an important benefit since the source number density is not well constrained.
In contrast, the numerator and denominator of $m$ (i.e., $M(2.5^\circ)$ and $M(10^\circ)$, respectively) are both 
strong functions of the source number density.

The differences between the large galaxy source model and the evenly distributed source model are not significant for source number densities $\rho \leq 10^{-4}~\mbox{Mpc}^{3}$. At higher source number densities, the nearest sources are expected to be within 10 Mpc. Thus, the assumed local over-density of sources becomes important in this density range.
The large galaxy source model predicts a larger value of $m$ since there is a greater probability in this model for a few nearby (and thereby bright) sources to contribute a large fraction of the flux, even though the total number of sources may be relatively high.

With regard to the observed value of $m$, there are two regions to consider.
Here, we describe the regions for a data set with $n=92$.
The first region is where $m > 1.6$. If the observed value of $m$ falls within this region, $\bar{s} < 10^\circ$ is favored. This signifies that the particles are not heavy nuclei like iron ($Z = 26$) and that magnetic bending in the galactic halo and in the intergalactic medium is not substantially greater than in the thin disk. 

The second region is where $m < 1.6$ ($m < 0.90$ for the evenly distributed source model). If the observed value of $m$ falls within this region, $\bar{s} > 2.5^\circ$ is favored. If $m << 1.6$, then $\bar{s} >> 2.5^\circ$ is plausible.
There are a large number of physical scenarios that are compatible with this situation. For example, perhaps $Z>1$ or perhaps magnetic bending in the galactic halo and intergalactic medium are substantially greater than in the thin disk. 

It is interesting that limits on $\bar{s}$ (and, in turn, $Z$ and the intervening magnetic fields) can be obtained even before the first CR source becomes apparent. Consider the following example.
For a data set of 92 events drawn from a isotropic distribution, $\left\langle  M(2.5^\circ) \right\rangle \approx 1.0~\mbox{degree}^{-2}$.
This value of $M(2.5^\circ)$ is equivalent to having 2 pairs with an angular separation of $0.8^\circ$. 
For a data set of 92 events where we observe 6 pairs (e.g., either six isolated pairs, two triples, or one quadruple) with an angular separation of $0.8^\circ$ and the clustering at higher angular scales is equivalent to the isotropic expectation, $m=2$.
This value of $m$ is firmly in the small source size region.
We should be attentive to what the data is telling us regarding the value of $\bar{s}$.


To summarize, we have described a method that gives an indication of the apparent angular size of CR sources.
The method can differentiate between a small source size ($\bar{s} = 2.5^\circ$) and a large source size ($\bar{s} = 10^\circ$) with a data set containing 92 events with energy $E > 57~\mbox{EeV}$. 
Such an analysis will be possible in the near future.
The method is not sensitive to input parameters that are not well constrained, e.g. the spatial distribution of sources and the source luminosity function.
An indication of the apparent angular size of CR sources can be used to set limits on the particle charge $Z$ and the intervening magnetic fields, both of which are not well constrained by other methods.





\acknowledgments{
I wish to thank my fellow Pierre Auger collaborators. I regret that there are too many of you to list by name. You have given me great support in this foray into cosmic ray physics. 

Some of the results in this paper have been derived using the HEALPix \cite{Gorski} software package. I also acknowledge the support of the National Science Foundation (award number PHY-083808) and the Michigan Space Grant Consortium.}


\begin{thebibliography}{999}

\bibitem{Pierre3} Pierre Auger Collaboration [J. Abraham, et al.], \textit{Evidence for suppression of the flux of cosmic rays above $4 x 10^{19}$ eV}, \prl{101}{2008}{061101} [astro-ph/08064302].

\bibitem{HiResSpec} R. U. Abbasi et al., \textit{First Observation of the Greisen-Zatsepin-Kuzmin Suppression}, \prl{100}{2008}{101101} [astro-ph/0703099].

\bibitem{Pierre1} Pierre Auger Collaboration [J. Abraham, et al.] 2007, \textit{Correlation of the highest energy cosmic rays with nearby extragalactic objects}, Science, 318, 939 [astro-ph/07112256].

\bibitem{Pierre2} Pierre Auger Collaboration [J. Abraham, et al.], \textit{Correlation of the highest-energy cosmic rays with the positions of nearby active galactic nuclei}, \app{29}{2008}{188} [astro-ph/07122843].

\bibitem{Younk} P. Younk, \textit{Estimating the Flux of the Brightest Cosmic Ray Source above $57 x 10^{18}$ eV}, Astrophys. J. Lett., 696, L40 [astro-ph/09033981].

\bibitem{YounkICRC} P. Younk, \textit{Prospects for Charge Particle Astronomy above $57 \times 10^{18}$ eV}, Proc. 31st Int. Cosmic Ray Conf. (Lodz, 2009).

\bibitem{Han} J. L. Han, R. N. Manchester, A.G. Lyne, G.J. Qiao, \& V.W. Straten, \textit{Pulsar Rotation Measures and the Large-Scale Structure of the Galactic Magnetic Field}, Astrophys. J., 642 (2006), 868 [astro-ph/0601357].

\bibitem{Men} H. Men, K. Ferriere, \& J.L. Han, \textit{Observational constraints on models for the interstellar magnetic field in the Galactic disk}, Astron. and Astrophys., 486 (2008), 819 [astro-ph/08053454].

\bibitem{Sun} X. H. Sun, W. Reich, A. Waelkens, \& T.A. En\ss lin, \textit{Radio observational constraints on Galactic 3D-emission models}, Astron. and Astrophys., 477 (2008), 573 [astro-ph/07111572].

\bibitem{Beck} R. Beck, \textit{Galactic and Extragalactic Magnetic Fields}, Proceedings of the 4th International Meeting on High Energy Gamma-Ray Astronomy, pp. 83-96 [astro-ph/08102923].

\bibitem{Ohno} H. Ohno, S. Shibata, \textit{The Random Magnetic-field in the Galaxy},  Mon. Not. R. Astron. Soc., 262, 953.

\bibitem{Harari2} D. Harari, S. Mollerach, \& Roulet, E., \textit{Lensing of ultra-high energy cosmic rays in turbulent magnetic fields}, J. High Energy Phys., 03 (2002), 045 [astro-ph/0202362].

\bibitem{AugerAngleRes} J. Abraham et al. [The Pierre Auger Collaboration], \textit{Angular Resolution of the Pierre Auger Observatory}, Proc. 29th Int. Cosmic Ray Conf. (Pune, 2005), V7, p. 17.

\bibitem{HiResAngleRes} R. U. Abbasi et al. [The HiRes Collaboration], \textit{Study of Small-Scale Anisotropy of Ultra-High-Energy Cosmic Rays Observed in Stereo by the High Resolution Fly's Eye Detector}, Astrophys. J. Lett., 610, L73.

\bibitem{Ryu} D. Ryu, H. Kang, J. Cho, \& S. Das, \textit{Turbulence and Magnetic Fields in the Large-Scale Structure of the Universe}, Science, 320(5878):909�912, May 2008 [astro-ph/08052466].

\bibitem{Bellido} J. Bellido [for the Pierre Auger Collaboration], \textit{Measurement of the average depth of shower maximum and its fluctuations with the Pierre Auger Observatory}, Proc. 31st Int. Cosmic Ray Conf. (Lodz, 2009).

\bibitem{Sanders} W. Saunders et al., \textit{The PSCz catalogue}, Mon. Not. R. Astron. Soc., 317 (2000), 55.

\bibitem{Karachentsev2} I. D. Karachentsev, V. E. Karachentseva, W. K. Huchtmeier, D. I. Makarov, \textit{A Catalog of Neighboring Galaxies}, Astron. J., 127, 2031.

\bibitem{Marzke} R. Marzke, L. Costa, M. Geller, \textit{The Galaxy Luminosity Function at $z < 0.05$: Dependence on Morphology}, Astrophys. J., 503, 617 [astro-ph/9805218].

\bibitem{Harari} D. Harari, S. Mollerach, \& Roulet, E., \textit{On the ultrahigh energy cosmic ray horizon}, J. Cosmology Astropart. Phys., 11 (2006), 012 [astro-ph/0609294].

\bibitem{Protheroe} R. J. Protheroe \& P.A. Johnson, \textit{Propagation of ultra high energy protons and gamma rays over cosmological distances and implications for topological defect models}, \app{4}{1996}{253} [astro-ph/9506119].

\bibitem{Allard} D.Allard, E. Parizot, \& A. V. Olinto, \textit{On the transition from galactic to extragalactic cosmic-rays: Spectral and composition features from two opposite scenarios}, \app{27}{2007}{61}.

\bibitem{Gorski} K.M. Gorski, E. Hivon, A.J. Banday, B.D. Wandelt, F.K. Hansen, M. Reinecke, \& M. Bartelmann, \textit{HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere}, Astrophys. J., 622, 759 [astro-ph/0409513].

\end{thebibliography}


\end{document}

0909.3483.extracted_bib
\begin{thebibliography}{999}

\bibitem{Pierre3} Pierre Auger Collaboration [J. Abraham, et al.], \textit{Evidence for suppression of the flux of cosmic rays above $4 x 10^{19}$ eV}, \prl{101}{2008}{061101} [astro-ph/08064302].

\bibitem{HiResSpec} R. U. Abbasi et al., \textit{First Observation of the Greisen-Zatsepin-Kuzmin Suppression}, \prl{100}{2008}{101101} [astro-ph/0703099].

\bibitem{Pierre1} Pierre Auger Collaboration [J. Abraham, et al.] 2007, \textit{Correlation of the highest energy cosmic rays with nearby extragalactic objects}, Science, 318, 939 [astro-ph/07112256].

\bibitem{Pierre2} Pierre Auger Collaboration [J. Abraham, et al.], \textit{Correlation of the highest-energy cosmic rays with the positions of nearby active galactic nuclei}, \app{29}{2008}{188} [astro-ph/07122843].

\bibitem{Younk} P. Younk, \textit{Estimating the Flux of the Brightest Cosmic Ray Source above $57 x 10^{18}$ eV}, Astrophys. J. Lett., 696, L40 [astro-ph/09033981].

\bibitem{YounkICRC} P. Younk, \textit{Prospects for Charge Particle Astronomy above $57 \times 10^{18}$ eV}, Proc. 31st Int. Cosmic Ray Conf. (Lodz, 2009).

\bibitem{Han} J. L. Han, R. N. Manchester, A.G. Lyne, G.J. Qiao, \& V.W. Straten, \textit{Pulsar Rotation Measures and the Large-Scale Structure of the Galactic Magnetic Field}, Astrophys. J., 642 (2006), 868 [astro-ph/0601357].

\bibitem{Men} H. Men, K. Ferriere, \& J.L. Han, \textit{Observational constraints on models for the interstellar magnetic field in the Galactic disk}, Astron. and Astrophys., 486 (2008), 819 [astro-ph/08053454].

\bibitem{Sun} X. H. Sun, W. Reich, A. Waelkens, \& T.A. En\ss lin, \textit{Radio observational constraints on Galactic 3D-emission models}, Astron. and Astrophys., 477 (2008), 573 [astro-ph/07111572].

\bibitem{Beck} R. Beck, \textit{Galactic and Extragalactic Magnetic Fields}, Proceedings of the 4th International Meeting on High Energy Gamma-Ray Astronomy, pp. 83-96 [astro-ph/08102923].

\bibitem{Ohno} H. Ohno, S. Shibata, \textit{The Random Magnetic-field in the Galaxy},  Mon. Not. R. Astron. Soc., 262, 953.

\bibitem{Harari2} D. Harari, S. Mollerach, \& Roulet, E., \textit{Lensing of ultra-high energy cosmic rays in turbulent magnetic fields}, J. High Energy Phys., 03 (2002), 045 [astro-ph/0202362].

\bibitem{AugerAngleRes} J. Abraham et al. [The Pierre Auger Collaboration], \textit{Angular Resolution of the Pierre Auger Observatory}, Proc. 29th Int. Cosmic Ray Conf. (Pune, 2005), V7, p. 17.

\bibitem{HiResAngleRes} R. U. Abbasi et al. [The HiRes Collaboration], \textit{Study of Small-Scale Anisotropy of Ultra-High-Energy Cosmic Rays Observed in Stereo by the High Resolution Fly's Eye Detector}, Astrophys. J. Lett., 610, L73.

\bibitem{Ryu} D. Ryu, H. Kang, J. Cho, \& S. Das, \textit{Turbulence and Magnetic Fields in the Large-Scale Structure of the Universe}, Science, 320(5878):909�912, May 2008 [astro-ph/08052466].

\bibitem{Bellido} J. Bellido [for the Pierre Auger Collaboration], \textit{Measurement of the average depth of shower maximum and its fluctuations with the Pierre Auger Observatory}, Proc. 31st Int. Cosmic Ray Conf. (Lodz, 2009).

\bibitem{Sanders} W. Saunders et al., \textit{The PSCz catalogue}, Mon. Not. R. Astron. Soc., 317 (2000), 55.

\bibitem{Karachentsev2} I. D. Karachentsev, V. E. Karachentseva, W. K. Huchtmeier, D. I. Makarov, \textit{A Catalog of Neighboring Galaxies}, Astron. J., 127, 2031.

\bibitem{Marzke} R. Marzke, L. Costa, M. Geller, \textit{The Galaxy Luminosity Function at $z < 0.05$: Dependence on Morphology}, Astrophys. J., 503, 617 [astro-ph/9805218].

\bibitem{Harari} D. Harari, S. Mollerach, \& Roulet, E., \textit{On the ultrahigh energy cosmic ray horizon}, J. Cosmology Astropart. Phys., 11 (2006), 012 [astro-ph/0609294].

\bibitem{Protheroe} R. J. Protheroe \& P.A. Johnson, \textit{Propagation of ultra high energy protons and gamma rays over cosmological distances and implications for topological defect models}, \app{4}{1996}{253} [astro-ph/9506119].

\bibitem{Allard} D.Allard, E. Parizot, \& A. V. Olinto, \textit{On the transition from galactic to extragalactic cosmic-rays: Spectral and composition features from two opposite scenarios}, \app{27}{2007}{61}.

\bibitem{Gorski} K.M. Gorski, E. Hivon, A.J. Banday, B.D. Wandelt, F.K. Hansen, M. Reinecke, \& M. Bartelmann, \textit{HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere}, Astrophys. J., 622, 759 [astro-ph/0409513].

\end{thebibliography}