format: text

text:
# 1 Introduction

## 1 Introduction

Observations suggest that Ultra High Energy Cosmic Rays (UHECRs) with energies $`\mathrm{few}\times 10^{20}`$ eV may exist in Nature (for reviews, see e.g. ). Usually it is assumed that such particles are of astrophysical origin.***Alternatively, it has been suggested that they may originate from the nearby decay of very massive particles  or topological defects . While different candidate sources of UHECRs have been proposed, there is a consensus that they should be at distances farther than $`𝒪(100)`$ Mpc. Such distances are difficult to reconcile with the observed energies of UHECRs in simple scenarios. Indeed, if the UHECRs are protons, any flux above $`5\times 10^{19}`$ eV – the well-known Greisen-Zatsepin-Kuzmin (GZK) cutoff  – is strongly attenuated by the photo-pion production on the Cosmic Microwave Background (CMB) photons. If the UHECRs are photons, their flux is in turn strongly attenuated by the scattering against the extragalactic radio background, seemingly ruling them out as the candidates for the highest energy UHECRs .

A way around the GZK cutoff might be a mechanism where the cosmic rays mostly consist of some very weakly interacting particles, which could travel long distances without being significantly absorbed. The only experimentally observed particles which could do this are neutrinos. A few mechanisms with neutrinos as the mediators, traveling most of the distance between a source and us have been suggested. Such proposals, however, suffer from the “converse” problem: how does a sufficient fraction of the neutrinos convert into nucleons or photons near the Earth. It has been argued  that, for neutrino masses of order eV, the UHE neutrinos can interact significantly with the cosmological background neutrinos through the $`Z`$ boson resonance. This mechanism however may require very high incoming fluxes and a very strong local clustering of the background neutrinos to conform with the observations .

It is therefore natural to ask if other particles (from beyond the standard model (SM) of particle physics) could be such mediators. Recently, pseudoscalars arising from the spontaneous breaking of an axial symmetry, or axions for short (similar to the QCD-axion ), have been considered in an alternative explanation for the apparent dimming of distant type Ia supernovae (SNe) . The idea was that a fraction of the photons emitted by a supernova could convert into axions in the presence of an intergalactic magnetic field en route to the Earth, thereby dimming the supernova and making it look farther away than it is . Since these axions can travel very far without significant interactions, and they can convert to photons (and vice-versa) in background magnetic fields (see, e.g. ), they are a natural candidate for a mediator which could help evade the GZK cutoff.

Ref.  investigated whether the QCD-axion, which could arise from a solution of the strong CP problem , could act as a mediator of UHECRs. Because the photon-axion coupling is of the form

$$\frac{a}{4M}\stackrel{~}{F}F,$$
(1)
where $`a`$ is the axion field, and $`F`$ the electromagnetic field strength, the evolution in the presence of a background magnetic field $`B`$ induces “flavor” (photon-axion) conversions. The trilinear interaction (1) gives rise to a bilinear term which is not diagonal in “flavor” space. The probability for the photon-axion oscillation in a homogeneous magnetic field is

$`P_{a\gamma }`$ $`=`$ $`{\displaystyle \frac{1}{1+\frac{\mathrm{\Delta }m^4M^2}{4^2B^2}}}\mathrm{sin}^2\left[{\displaystyle \frac{y}{2}}\sqrt{{\displaystyle \frac{B^2}{M^2}}+{\displaystyle \frac{\mathrm{\Delta }m^4}{4^2}}}\right],`$ (2)
$`\mathrm{\Delta }m^2`$ $``$ $`m_\gamma ^2m_a^2`$ (3)
where $`m_\gamma `$ is the effective photon mass in the intergalactic plasma, $`m_a`$ is the axion mass, $``$ the energy of the particles, and $`y`$ the path length. For a QCD axion, $`10^5\mathrm{eV}\stackrel{<}{}m_{a,\mathrm{QCD}}\stackrel{<}{}10^1`$ eV, so the conversion $`a\gamma `$ is completely negligible . One way around this problem pursued in  was to assume that the axion does not convert into photons, but that it interacts with the Earth’s atmosphere due to its relatively large coupling to gluons. Here we will follow a different direction. In the photon-axion model for the dimming of the SNe a much lighter axion is invoked  (see also ). Then the axion mass will not suppress the mixing even for the optical photons contributing to the SNe luminosity, let alone for the UHECRs with energies larger by nearly 20 orders of magnitude. This observation is the starting point of this paper, where we will show that the coupled photon-axion system could indeed give rise to the highest energy cosmic rays. We will also compare the preferred region of the parameter space where the transmission of ultra-high energy photons is maximal to the parameters needed to explain the dimming of the SNe, and show that the UHECR flux due to the photon-axion mixing is especially enhanced for the choice of the parameters needed to account for the SNe dimming in . Since the mean free path $`l_{\mathrm{GZK}}`$ of the photons is very sensitive to their energy, and since the energy of the UHECRs is nearly 20 orders of magnitude higher than the ones of the SNe, this coincidence is very intriguing.

Our model for the transmission of UHECRs works as follows. Initially, a faraway astrophysical source releases a flux of very high energy unpolarized photons. A small fraction of these photons converts into axions before being depleted by the interactions with the background (radio and CMB) photonsThe photons which do interact with the background photons give rise to secondary photons which can also propagate along the beam with a lower energy than the primary ones. To simplify the present discussion we will neglect those and concentrate on particles which have never interacted with the background photons.. At distances $`dl_{\mathrm{GZK}}`$, defined as the mean free path of the photons in absence of the mixing with axions, the beam is mostly comprised of axionsAlternatively, one can start with only axions from the very beginning, as in ref. , assuming that they are produced directly at the source.. These axions will gradually convert back into photons, but most of these photons will be depleted by the interaction with background photons. However, the conversion rate per distance traveled is low, allowing most of the axions to travel unimpeded, and keep replenishing the photon beam. This will ensure that a fraction of the photons will survive over distances much larger than one would expect from the GZK length of UHE photons.

The paper is organized as follows: in Section 2 we describe the evolution of the photon-axion system in the presence of the GZK cutoff, making the (unrealistic) assumption of a homogeneous magnetic field of order $`10^9`$ Gauss. In this case the evolution of the system can be calculated exactly, and we show that there indeed exists a long-lived component in the photon beam. In Section 3 we examine a more realistic situation, where the magnetic field consists of domains of size $``$ Mpc, with a random orientation of the field inside each domain, again with a strength $`10^9`$ Gauss . We calculate the surviving intensity of the photon beam as a function of the distance traveled and the parameters of the photon-axion system, and show that the highest survival probabilities are expected for the same parameters that are preferred for the SNe dimming mechanism of . We also discuss possible peculiar signatures of this UHECR mediation mechanism arising from the fact that most of the observed photons would be generated from a predominantly axion beam within the last few magnetic domains around the Earth. We estimate the initial photon flux needed for this mechanism to work, and finally conclude in Section 4.

## 2 Evolution of the photon-axion system

The evolution of the photon-axion system propagating along the $`y`$ direction is governed by

$$\left\{\frac{d^2}{dy^2}+^2+\left(\begin{array}{cc}i\mathrm{\Gamma }()& i\frac{B}{M}\\ i\frac{B}{M}& 0\end{array}\right)\right\}\left(\begin{array}{c}|\gamma \\ |a\end{array}\right)=0$$
(4)
where we have Fourier-transformed the fields to the energy ($``$) picture. Here $`\stackrel{}{B}`$ is the extra-galactic magnetic field and $`B=\stackrel{}{e}\stackrel{}{B}|\stackrel{}{B}|`$ is its projection on the photon polarization $`\stackrel{}{e}`$. In a constant magnetic field $`|\gamma `$ denotes the photon polarization component along $`\stackrel{}{B}`$, $`|a`$ the axion, while the photon polarization $`|\gamma _{}`$ orthogonal to $`\stackrel{}{B}`$ decouples from $`|a`$ and evolves independently. In the next section we will discuss the evolution of the photon-axion beam through a sequence of domains with a different orientation of the magnetic field in each domain, where both photon polarizations must be retained. In both cases we take $`|\stackrel{}{B}|\mathrm{few}10^9`$ G. Throughout our analysis we can ignore the axion mass and the effective photon mass, since they are much smaller***We assume that $`m_a^2B/M`$; this ensures that the photon-axion oscillations are unsuppressed, but also excludes the QCD axion. than the beam energy $`10^{20}eV`$, and $`B/M`$. Instead we include the (energy dependent) decay rate $`\mathrm{\Gamma }`$, which parameterizes the decrease of the photon intensity due to their interaction with the background photons. If we rewrite $`\mathrm{\Gamma }/l_{\mathrm{GZK}}`$, in the absence of the coupling to the axions $`l_{\mathrm{GZK}}`$ gives the mean free path of the photons. For $`=310^{20}eV`$ (which we will consider in our numerical examples later on), $`l_{\mathrm{GZK}}6.4`$ Mpc .

For a constant $`B`$, the solution of Eqs. (4) is

$`\left(\begin{array}{c}|\gamma \\ |a\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}\frac{1+\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}& \frac{\delta }{\sqrt{14\delta ^2}}\\ \frac{\delta }{\sqrt{14\delta ^2}}& \frac{1\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}\end{array}\right)\mathrm{e}^{iy+\lambda _1y}\left(\begin{array}{c}|\gamma \\ |a\end{array}\right)_0+`$ (11)
$`+`$ $`\left(\begin{array}{cc}\frac{1\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}& \frac{\delta }{\sqrt{14\delta ^2}}\\ \frac{\delta }{\sqrt{14\delta ^2}}& \frac{1+\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}\end{array}\right)\mathrm{e}^{iy+\lambda _2y}\left(\begin{array}{c}|\gamma \\ |a\end{array}\right)_0,`$ (16)
where the subscript zero denotes the initial amplitudes at the source at $`y=0`$, and where we use

$`\lambda _1`$ $``$ $`{\displaystyle \frac{1}{4l_{\mathrm{GZK}}}}\left[1+\sqrt{14\delta ^2}\right],`$
$`\lambda _2`$ $``$ $`{\displaystyle \frac{1}{4l_{\mathrm{GZK}}}}\left[1\sqrt{14\delta ^2}\right],`$
$`\delta `$ $``$ $`{\displaystyle \frac{Bl_{\mathrm{GZK}}}{M}}0.11\left({\displaystyle \frac{B}{1\mathrm{nG}}}\right)\left({\displaystyle \frac{10^{11}\mathrm{GeV}}{M}}\right)\left({\displaystyle \frac{l_{\mathrm{GZK}}}{\mathrm{Mpc}}}\right).`$ (17)
Notice that (16) is regular for all values of $`\delta `$. One can see this by rewriting (16) as

$`\left(\begin{array}{c}|\gamma \\ |a\end{array}\right)=\mathrm{e}^{iy}\mathrm{e}^{\frac{y}{4l_{\mathrm{GZK}}}}\left(\begin{array}{cc}𝒞𝒮& 2\delta 𝒮\\ 2\delta 𝒮& 𝒞+𝒮\end{array}\right)\left(\begin{array}{c}|\gamma \\ |a\end{array}\right)_0`$ (24)
$`𝒞\mathrm{cosh}\left[\sqrt{14\delta ^2}{\displaystyle \frac{y}{4l_{\mathrm{GZK}}}}\right],𝒮{\displaystyle \frac{\mathrm{sinh}\left[\sqrt{14\delta ^2}\frac{y}{4l_{\mathrm{GZK}}}\right]}{\sqrt{14\delta ^2}}},`$ (25)
with $`𝒞`$ and $`𝒮`$ real numbers for any choice of $`\delta `$.

We are interested in the mean free path of the photons in the presence of the coupling to the axions. The initial photon beam is taken to be composed of unpolarized photons, that is it contains an equal mixture of $`|\gamma `$ and $`|\gamma _{}`$, with the total intensity normalized to unity. Then the intensities of the surviving photons and axions after the distance $`y`$ away from the source, using eq. (16) $`I_\gamma `$ and $`I_a`$ are

$`I_\gamma `$ $`=`$ $`{\displaystyle \frac{1}{2}}\{\mathrm{e}^{\frac{y}{l_{\mathrm{GZK}}}}+{\displaystyle \frac{1}{4\left(14\delta ^2\right)}}[(1+\sqrt{14\delta ^2})\mathrm{e}^{\frac{1+\sqrt{14\delta ^2}}{4l_{\mathrm{GZK}}}y}`$
$`(1\sqrt{14\delta ^2})\mathrm{e}^{\frac{1\sqrt{14\delta ^2}}{4l_{\mathrm{GZK}}}y}]^2\},`$
$`I_a`$ $`=`$ $`{\displaystyle \frac{\delta ^2}{2\left(14\delta ^2\right)}}\left[\mathrm{e}^{\frac{1+\sqrt{14\delta ^2}}{4l_{\mathrm{GZK}}}y}\mathrm{e}^{\frac{1\sqrt{14\delta ^2}}{4l_{\mathrm{GZK}}}y}\right]^2.`$ (26)
The first term in (26) gives the intensity of the photons orthogonal to $`\stackrel{}{B}`$, which do not mix with the axion and thus simply decay away as is expected from the GZK mechanism, while the remaining terms give the intensities of the coupled axion-photon system. They have an easy interpretation both in the limit of large and small $`\delta `$. For a large mixing ($`\delta 1`$),

$`I_\gamma `$ $``$ $`{\displaystyle \frac{1}{2}}\mathrm{e}^{\frac{y}{l_{\mathrm{GZK}}}}+{\displaystyle \frac{1}{2}}\mathrm{e}^{\frac{y}{2l_{\mathrm{GZK}}}}\mathrm{cos}^2\left({\displaystyle \frac{\delta y}{2l_{\mathrm{GZK}}}}\right)`$
$`I_a`$ $``$ $`{\displaystyle \frac{1}{2}}\mathrm{e}^{\frac{y}{2l_{\mathrm{GZK}}}}\mathrm{sin}^2\left({\displaystyle \frac{\delta y}{2l_{\mathrm{GZK}}}}\right).`$ (27)
This displays the presence of an eigenstate which rapidly oscillates between the original $`|\gamma `$ and $`|a`$ states. In terms of the original parameters, the oscillation length is $`L_{\mathrm{osc}}M/B`$ and, not surprisingly, it is independent of $`l_{\mathrm{GZK}}`$. Thus, $`\delta l_{\mathrm{GZK}}/L_{\mathrm{osc}}`$ gives the number of the photon-axion oscillations within the mean free path $`l_{\mathrm{GZK}}`$. Thus we can view the coupled system as particles which are photon-like half of the time, and axion-like half of the time, so that their probability to interact with the background photons is one-half of the probability of the decoupled photons $`|\gamma _{}`$. This explains why the mean free path of this state is $`2l_{\mathrm{GZK}}`$.In the case of a magnetic field randomly oriented along different domains and of a strong mixing $`\delta `$, we find numerically a mean free path of $`3l_{\mathrm{GZK}}/2`$, showing that the quanta are equally ”shared” between the axion and the two photon polarizations states.

In the case of small mixing, $`\delta 1`$, eq. (26) reduces to

$`I_\gamma `$ $``$ $`{\displaystyle \frac{1}{2}}\mathrm{e}^{\frac{y}{l_{\mathrm{GZK}}}}+{\displaystyle \frac{1}{2\left(14\delta ^2\right)}}\left[\left(1\delta ^2\right)\mathrm{e}^{\frac{1\delta ^2}{2l_{\mathrm{GZK}}}y}\delta ^2\mathrm{e}^{\frac{\delta ^2y}{2l_{\mathrm{GZK}}}}\right]^2`$
$`I_a`$ $``$ $`{\displaystyle \frac{\delta ^2}{2}}\left[\mathrm{e}^{\frac{y}{2l_{\mathrm{GZK}}}}\mathrm{e}^{\frac{\delta ^2y}{2l_{\mathrm{GZK}}}}\right]^2.`$ (28)
This shows the presence of a long-lived mode, characterized by the mean free path $`l_{\mathrm{GZK}}/\delta ^2l_{\mathrm{GZK}}`$. To see how it arises, we consider a simple example depicted in Fig. 1, where we plot $`I_{\gamma ,a}`$ as functions of $`y`$, taking the parameters: $`l_{\mathrm{GZK}}=6.4\mathrm{Mpc},B=1\mathrm{n}\mathrm{G},M=410^{11}`$ GeV (this value of the photon-axion coupling is close to its experimental bound ) corresponding to $`\delta 0.18`$. These are also the values of $`B`$ and $`M`$ that have been used in .

The photon intensity $`I_\gamma `$ initially decreasesWe neglect the secondary photons formed by the scattering of the primaries against the background photons. as $`\mathrm{exp}\left(y/l_{\mathrm{GZK}}\right)`$. Meanwhile, some photons convert into axions, and the intensity $`I_a`$ (initially taken to be zero§§§This is a conservative assumption, since if initial axions are also produced at the source, then our mechanism will still work, the only difference being that the required sources can have lower intensities.) increases. For a small mixing, the quantity $`\delta l_{\mathrm{GZK}}/L_{\mathrm{osc}}`$ gives the amplitude for the process $`|\gamma |a`$ to occur within the mean free path of the photons. Thus, $`\delta ^2`$ is approximately the fraction of the original photons which converts into axions before being depleted by the interaction with the background photons. This is confirmed by Eq. (28) and by the numerical evolution shown, which indicate that the axion intensity peaks at $`\delta ^2`$. From this point on, $`I_a>I_\gamma `$ and the subsequent photon-axion oscillations lead to a decrease of $`I_a`$. The number of axions decreases by a fraction $`\delta ^2`$ per each distance $`\mathrm{\Delta }y=l_{\mathrm{GZK}}`$ traveled. To leading order in $`\delta ^2`$, this effect is not compensated by those photons which convert back into axions, since in the case of small mixing the photons in the beam are more likely to be dissipated by the background photons than oscillate back into axions. Thus, $`I_aI_{a,max}\left(1\delta ^2\right)^{y/l_{\mathrm{GZK}}}\delta ^2\mathrm{exp}\left(\delta ^2y/l_{\mathrm{GZK}}\right)`$.

This argument explains the presence of an eigenstate with the mean free path $`l_{\mathrm{GZK}}/\delta ^2l_{\mathrm{GZK}}`$ in Eq. (28). The axion intensity in this mode relative to the total initial photon intensity is small, of order $`\delta ^2`$, due to the small probability for $`|\gamma |a`$ oscillations. The photon-like component of this mode is suppressed by a further power of $`\delta ^2`$, because these photons arise from the double conversion $`|\gamma |a|\gamma `$. Most of these photons disappear quickly, within a distance $`yl_{\mathrm{GZK}}`$ because of the interaction with background photons. However, the incessant conversion of the surviving axions in the beam back into photons provides for a long-lasting source of photons, yielding a slowly decreasing photon intensity $`I_\gamma \delta ^4\mathrm{exp}\left(\delta ^2y/l_{\mathrm{GZK}}\right)`$. The rise of this “second population” of photons is clearly visible in Fig. 1. The distance $`y_{}`$ at which these photons reach the level of the surviving intensity of the primary photons in the beam which never underwent any interactions is roughly given by $`e^{y_{}/l_{GZK}}\delta ^4`$. In the present example, $`y_{}45`$ Mpc, in good agreement with the value shown in the figure.

## 3 Consequences for the UHECR

Until now we have considered the evolution of the photon-axion system in a constant magnetic field. However, the intergalactic magnetic field is unlikely to be completely uniform. Instead, it is more reasonable to take the magnetic fields with a strength close to the observed upper bound of few $`\times 10^9`$ G to be homogeneous within domains of a typical size $``$ Mpc , but with their orientation randomly varying from domain to domain (the same assumptions were made in ). We will now look at the evolution of the photon-axion beam in such a universe. Photons with polarization perpendicular to $`\stackrel{}{B}`$ do not couple to axions, so the photon-axion mixing involves photons whose helicity changes from domain to domain . If it were not for the interaction with the background photons, the quanta in the beam would be roughly equipartitioned between the axion component and the two photon polarizations after traversing several domains. This is what is expected to happen to the visible photons released in SNe explosions, as explained in . For the range of energies we are interested in, the interaction with the background photons cannot be neglected, and so the depletion effect described in the previous section must be taken into account.

It is easy to extend the formulae from the previous section to the general situation. If the beam is still taken to travel along the $`y`$ direction, it can be described by a vector in the basis $`(|\gamma _x,|\gamma _z,|a)`$, where $`|\gamma _{\{x,z\}}`$ denote photons polarized along the $`\{x,z\}`$ axis. Inside each domain, the transfer matrix is such that it can be rotated to the $`2\times 2`$ matrix (16) for the coupled axion-photon components, plus a third term $`\mathrm{exp}\left(y/2l_{\mathrm{GZK}}\right)`$ describing the evolution of the decoupled photon polarization orthogonal to $`\stackrel{}{B}`$.

The resulting transfer equation is:

$`\left(\begin{array}{c}\gamma _x\\ \gamma _z\\ a\end{array}\right)=\mathrm{e}^{iEy}\left[T_0\mathrm{e}^{\lambda _0y}+T_1\mathrm{e}^{\lambda _1y}+T_2\mathrm{e}^{\lambda _2y}\right]\left(\begin{array}{c}\gamma _x\\ \gamma _z\\ a\end{array}\right)_0`$ (35)
where $`\lambda _{1,2}`$ have been given above,

$$\lambda _0\frac{1}{2l_{\mathrm{GZK}}},$$
(36)
and the three transfer matrices are

$`T_0`$ $``$ $`\left(\begin{array}{ccc}\mathrm{sin}^2\theta & \mathrm{cos}\theta \mathrm{sin}\theta & 0\\ \mathrm{cos}\theta \mathrm{sin}\theta & \mathrm{cos}^2\theta & 0\\ 0& 0& 0\end{array}\right),`$ (40)
$`T_1`$ $``$ $`\left(\begin{array}{ccc}\frac{1+\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}\mathrm{cos}^2\theta & \frac{1+\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}\mathrm{cos}\theta \mathrm{sin}\theta & \frac{\delta }{\sqrt{14\delta ^2}}\mathrm{cos}\theta \\ \frac{1+\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}\mathrm{cos}\theta \mathrm{sin}\theta & \frac{1+\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}\mathrm{sin}^2\theta & \frac{\delta }{\sqrt{14\delta ^2}}\mathrm{sin}\theta \\ \frac{\delta }{\sqrt{14\delta ^2}}\mathrm{cos}\theta & \frac{\delta }{\sqrt{14\delta ^2}}\mathrm{sin}\theta & \frac{1\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}\end{array}\right)`$ (44)
$`T_2`$ $``$ $`\left(\begin{array}{ccc}\frac{1\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}\mathrm{cos}^2\theta & \frac{1\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}\mathrm{cos}\theta \mathrm{sin}\theta & \frac{\delta }{\sqrt{14\delta ^2}}\mathrm{cos}\theta \\ \frac{1\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}\mathrm{cos}\theta \mathrm{sin}\theta & \frac{1\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}\mathrm{sin}^2\theta & \frac{\delta }{\sqrt{14\delta ^2}}\mathrm{sin}\theta \\ \frac{\delta }{\sqrt{14\delta ^2}}\mathrm{cos}\theta & \frac{\delta }{\sqrt{14\delta ^2}}\mathrm{sin}\theta & \frac{1+\sqrt{14\delta ^2}}{2\sqrt{14\delta ^2}}\end{array}\right),`$ (48)
with $`\theta `$ being the angle between the magnetic field and the $`x`$ axis. The overall transfer matrix is then simply the product of the transfer matrices in each domain.

The evolution of the system is qualitatively similar to the one for a constant magnetic field illustrated in Fig. 1. However, the interference effects modify the quantitative results whenever the domain size is comparable to, or smaller than, the absorption length $`l_{\mathrm{GZK}}`$. For example, consider the evolution of the photon intensity illustrated in Fig. 2. The curve denoted “exact evolution” describes $`I_\gamma `$ as a function of the distance from the source. All parameters are as in Fig. 1 except that $`\stackrel{}{B}`$ is randomly oriented in each domain. The random orientation of $`\stackrel{}{B}`$ is responsible for the stochastic character of $`I_\gamma `$ seen in the figure. At these distances from the source the intensity in axions $`I_a`$ is much greater than the one in photons, decreasing smoothly as in Fig. 1. The overall intensity $`I_\gamma +I_a`$ is thus also smoothly decreasing. By comparing the value of $`I_\gamma `$ of the two figures (recalling that in both cases the photon intensity is normalized to unity at the source), one notices that $`I_\gamma `$ is more suppressed by a randomly pointing $`\stackrel{}{B}`$. This effect might be interpreted as a decrease of the effective mixing parameter $`\delta `$. As a consequence, fewer photons will convert into axions before being depleted by the interaction with the background photons. However, for the very same reason, the decrease of $`I_\gamma `$ at distances $`l_{\mathrm{GZK}}`$ appears to be much milder in the numerical simulations than the one computed for a constant magnetic field.

An interesting consequence for the transmission of UHECRs is that for sufficiently distant sources the only photons in the beam are the ones which have been generated by the axions within the previous $`45l_{\mathrm{GZK}}`$ distances. Otherwise they would have been depleted by the interactions with the background photons.***We remark that this concerns the primary photons, i.e. the ones which have not interacted with the background photons. Since within each domain only the polarization along the magnetic field is actually generated, the UHECRs reaching the Earth would be very sensitive to the orientation of $`\stackrel{}{B}`$ in domains within the sphere of radius $`5l_{\mathrm{GZK}}`$ from us. Consequently the incoming photons should have slightly different polarizations along different lines of sight. The magnitude of this effect depends on the relative size of the domains and on the value of $`l_{\mathrm{GZK}}`$. If the beam crosses several domains in the last few $`l_{\mathrm{GZK}}`$ distances, the differences in the polarization may average out.

We also note that the presence of domains of $`\stackrel{}{B}`$ with different orientation along different lines of sight can lead to differences in the intensities of UHECRs from different directions, even if the distribution of sources were homogeneous and isotropic. This is because only the magnetic field perpendicular to the line of travel contributes to the photon-axion mixing in each domain. Thus, we expect lower fluxes from regions in which the magnetic field is mostly oriented in the direction pointing toward the Earth. At present there is not enough experimental data about UHECRs and extragalactic magnetic fields to enable us to test these predictions. However, this is a very distinctive signature that may eventually be tested experimentally. We illustrate these effects in Fig. 2, where we have artificially removed the photons in the beam at a distance $`d=100`$ Mpc from the source, and compared the subsequent evolution to the exact one, without this artificial suppression. One can easily see that the two curves approach each other within a few $`l_{\mathrm{GZK}}`$ distances. Since in this case the evolution is stochastic, the curves themselves behave very differently according to different random choices for the orientations of $`\stackrel{}{B}`$. However, in all the different simulations the two curves approach each other very quickly, irrespectively of the particulars of the magnetic configuration.

In Fig. 2 the mild decrease of the photon intensity at large distances is hidden by its stochastic behavior. To see the decrease, one has to average the evolution over several configurations of $`\stackrel{}{B}`$. This averaging procedure is appropriate if we discuss UHECRs independently of their arrival directions, which have crossed different, uncorrelated configurations of the magnetic field. The remaining figures that we show have been obtained using this averaging procedure.

In Fig. 3 we display the evolution of the photon intensity as a function of the distance from the source (vertical axis) and of the mixing parameter $`\delta `$ (horizontal axis). As before $`l_{\mathrm{GZK}}=6.4`$ Mpc. The photon intensity drops rapidly both at low and high values of $`\delta `$. We note that the surviving photon intensity is maximized for $`\delta `$ between $`0.1`$ and order unity. This gives the region of preferred values of the photon-axion coupling parameter $`M`$ and the magnetic field $`B`$ where our mechanism is efficient (see Eq. (17)). Interestingly, the choice of parameters adopted in  to explain the SNe dimming in terms of the photon-axion conversion is precisely within this region.

Ultimately we would like to find how the photon-axion mixing affects the spectrum of UHECRs observed on Earth. A precise computation, which requires the inclusion of the secondary particles produced when a photon in the beam interacts with a background photon, is beyond the scope of this paper. However, a rough but illustrative estimate can be obtained relatively simply. We begin by considering an explicit example of the evolution of a primary photon spectrum from a single source at $`d=500`$ Mpc from the Earth, shown in Fig 4. The intensity at the source is assumed to scale as $`E^3`$ (the overall normalization is irrelevant for the present discussion). Even without the mixing with axions the spectrum of primary photons reaching the Earth is strongly modified due to their interaction with the background photons. The reason is that the mean free path $`l_{\mathrm{GZK}}`$ is strongly dependent on their energy. For the energy range shown in Fig. 4, it ranges from $`0.3`$ Mpc for $`E=10^{18}`$ eV up to $`100`$ Mpc for $`E=10^{22}`$ eV (we have used the results of  summarized in , see these references for details). The strong increase of $`l_{\mathrm{GZK}}`$ in this interval explains why the second spectrum in the figure is peaked at high energies, although the original spectrumThe local maximum at $`E210^{19}`$ eV corresponds to a local maximum of $`l_{\mathrm{GZK}}`$, due to the fact that at smaller energies the incoming photons are mostly affected by CMB photons, while at higher energies they are most affected by radio photons. See  for details. scaled as $`E^3`$.The curves shown in Fig. 4 confirm a much lower depletion of primary photons when the mixing with axions is taken into account (as in all the previous cases, the axion-photon mass parameter is taken to be $`M=410^{11}`$ GeV).

We finally relate this decreased depletion of $`I_\gamma `$ to an apparent absence of a GZK cutoff in the UHECR spectrum. In Fig. 5 we combined a spectrum of protons of extragalactic origin with the spectrum of primary photons. Both spectra have been computed for an isotropic and homogeneous distribution of sources at distances $`100\mathrm{Mpc}<d<1000\mathrm{Mpc}`$ from the Earth. For both populations, the intensity at the source is again taken to have a $`IE^3`$ dependence. Notice that the proton spectrum $`I_p`$ exhibits the GZK cut-off,For the proton evolution, we have followed the results of  summarized in . while – due to their coupling with the axions – this is not the case for the photon spectrum $`I_\gamma `$. By an appropriate relative choice of $`I_p`$ and $`I_\gamma `$ at the source, it is thus possible to match the two spectra such that the photon flux “catches up” precisely where the proton flux exhibits the GZK cutoff. The overall spectrum can thus be extended at energies well above $`E_{\mathrm{GZK}}`$.

Before concluding, we stress the limitations of the computations summarized in the last two figures. The main one is related to the fact that $`I_\gamma `$ does not include the secondary photons arising when the photons in the beam scatter against the background photons. This will lead to a significant modification of $`I_\gamma `$. Fig. 5 has been obtained for a ratio $`I_\gamma /I_p=25,000`$ at the source. This would amount to a flux of about $`10^{23}\left(m^2srsGeV\right)^1`$ on Earth, if the photons were travelling unimpeded. We expect that a smaller ratio will be required once secondary photons are taken into account. Secondly, in all of the above computations we have neglected the decrease of the energy (with a consequent change of $`l_{\mathrm{GZK}}`$) of the quanta in the beam due to the cosmological redshift (for sources at $`d=1000`$ Mpc, this correspond to about a $`30\%`$ decrease of the energy). Thus, the computation summarized in Fig. 5 should be understood only as an order of magnitude estimate and as an illustration of the effect which can be expected. The significance of this estimate will be highly dependent on which sources for the photons one is assuming. The problem of accelerating particles at such energies in different astrophysical environments is completely independent from the problem of their propagation across intergalactic media. Thus it is sensible to separate these problems from each other. Therefore in the present paper we have focused only on the latter, the propagation problem.

We close this section with a note on the nature of the UHECRs. Once UHECRs hit the atmosphere of the Earth, they give rise to extensive air showers. Different incident particles generate distinct longitudinal profiles, that is the amount of particles as a function of the amount of atmosphere penetrated by the cascade. Several analyses, for example the one reported by the Fly’s Eye Collaboration in , suggest a change from an iron dominated composition of the incoming particles at $`10^{17}`$ eV to a proton dominated composition near $`10^{19}`$ eV. At higher energies, the current statistics are too poor for a conclusive claim, and one still cannot exclude photons as a relevant fraction of UHECRs above the GZK cut-off . The most accurate method to determine the nature of UHECR of such energies is the study of the muon content of inclined showers (that is, the ones with incident zenith angle $`\theta >60^0`$). The strongest bounds quoted in the literature  have been obtained in , where inclined air showers recorded by the Haverah Park  detector have been analyzed. The results of  are affected by uncertainties in the parameterization of the flux of cosmic rays. Assuming the parameterization given in , the authors of  obtain that less than $`48\%`$ of the observed events above $`10^{19}`$eV can be photons with a $`95\%`$ confidence level. For energies higher than $`4\times 10^{19}`$eV, they find instead that the content of photons has to be less than $`50\%`$ of the overall flux. Different results are obtained based on the parameterization given in the more recent work . The two upper bounds change to $`25\%`$ for energies above $`10^{19}`$eV, while only to $`70\%`$ for energies greater than $`4\times 10^{19}`$eV .

Hence, we see that a mixed composition of proton and photons in the UHECR above the GZK cutoff is compatible with the present experimental limits. Thus the mechanism we propose here is in agreement with the data at the moment: it relies precisely on the assumption that both protons and photons give a sizeable contribution to the UHECR at energies around $`10^{19}`$eV. The present experimental results do not yield clear-cut constraints at higher energies, as manifested by the discrepancy between the analysis based on the different parameterizations of  and . On the contrary, the higher statistics and better energy determination expected in the forthcoming Auger Observatory  have the power of greatly improving the present bounds. Therefore it may either discover photon UHECRs or yield significantly stronger bounds on the dynamics of photon-axion mixing we employ here. It will thus be of interest to subject the mechanism we propose to more scrutiny.

One more concern about photons being UHECRs could be that most of the photons that do interact with the background (the secondary photons) would cascade down to GeV energies and contribute to the diffuse gamma ray background. If the survival probability of the UHE photons were too low, then the model would have predicted too many secondary photons and thus a gamma ray background that is too large. The average survival probability from Fig. 3 for the optimal values of $`\delta `$ are of order $`10^4`$. The ratio of total energies contained in the 1 GeV gamma ray background to the energies contained in UHECRs at $`10^{20}`$ eV can be obtained from comparing the measured EGRET flux  to the various measurements of UHECR fluxes. This ratio turns out to be of order $`10^310^4`$, depending on whether a high or low number for the UHECR flux is adopted. The lower value seems to fit in our model very well, while for the higher value there could seemingly be a marginal conflict (of at most an order of magnitude) with the gamma ray background bounds. However, the ratio of primary to secondary photons can be higher than the survival probabilities in Fig. 3 for various different reasons. First, there could be an initial axion component in the beam which would dramatically change the prediction for the ratio of primary to secondary photons. Note that if the initial beam contained an order unity fraction of the energy in axions, the primary photon to secondary photon ratio would also be of order unity. Thus it is clear that even a relatively small admixture of axions in the initial beam can drastically improve the primary to secondary photon ratio. Second, the conversion rate within our galaxy is much larger since the galactic magnetic field is $`10^3`$ times bigger than the assumed value of the intergalactic magnetic field. If a somewhat larger value of magnetic field persists for a few kiloparsecs, then the survival probabilities will be increased by a factor of few, and this on its own could be sufficient to bring down the number of secondary photons. Lastly, as we have stressed above, only a fraction of UHECRs may be photons. Then the total number of secondary photons would fall below the diffuse gamma ray bounds even if the primary to secondary photon ratio stays the same. In this way the signatures of our mechanism would still persist without any violations of the diffuse gamma ray background bounds.

Finally, a photon component of UHECRs may be related to the modification of their spectrum by the interaction with the background photons, once the secondary photons are also taken into account in the propagation codes. After scattering with the background photons, an incoming photon of energy $`_0`$ produces secondary photons which accumulate at energies $`<_0`$. Starting with some flux at the source and normalizing the final spectrum to the observed amount of events above the GZK cutoff, one would typically find an excess of events at lower energies . However here the mixing with the axions can significantly modify the final spectrum precisely in the direction of “balancing out” this excess at lower energies. Indeed, the energy of the photons continuously supplied by the axions is given by the initial energy at the source decreased only by the expansion of the Universe, and it is greater than that of the secondary photons present in the beam. A detailed computation of the final spectrum, including both primary and secondary photons will be an important test of the viability of the mechanism proposed here.

## 4 Conclusions

An interesting property of axions, which arise in many extensions of the SM, is that they mix with photons in an external magnetic field. Photon-axion conversions in extragalactic magnetic fields can have relevant consequences for astrophysics and cosmology. We have shown that the same parameters (axion mass, coupling scale, magnetic field configurations) which provide an alternative explanation for the dimming of distant SNe  could also imply an enhancement of the flux of super-GZK photons arriving to the Earth from faraway sources. Whether this enhancement is observable clearly depends on the distribution of the sources and on the intensity of the high-energy photon spectrum they emit. We have seen here that, for realistic assumptions, super-GZK photons can be detectable, and may provide an explanation for the origin of super-GZK events in the UHECR spectrum, if the required sources of very high energy photons exist. The viability of this mechanism can be tested by improved simulations of the photon propagation, and ultimately by the increase of UHECR data expected over the next few years.

## Acknowledgments

We thank J. Arons, J. Feng, H.P. Nilles, and G. Raffelt for useful discussions, and to Dimitry Semikoz for useful comments on the manuscript. C.C. is supported in part by the DOE OJI grant DE-FG02-01ER41206 and in part by the NSF grant PHY-0139738. N.K. is supported in part by a Research Innovation Award from the Research Corporation. J.T. is supported by the U.S. Department of Energy under contract W-7405-ENG-36.


text:
# Tunnelling Amplitudes from Perturbation Expansions

## I Introduction

1. Tunnelling processes govern many important physical phenomena. Their theoretical description requires the calculation of the contribution of critical bubbles to the partition function, including their fluctuation entropy Langer . The latter can be found at most to one-loop order. Higher loop effects are prohibitively complicated Hagen . It would be of great advantage tunnelling amplitudes could be derived from ordinary perturbation expansions around the free theory since these can be performed to many loops Verena . The difficulty arising in such a program is that tunnelling amplitudes are described by the analytic continuation of divergent Borel-summable power series expansions in the coupling constant $`g`$ to negative $`g`$ where they become non-Borel-summable. None of the currently known resummation schemes Verena ; Resumm  is able to deal with such expansions. Some time ago it was suggested that a resummation is possible by variational perturbation theory Karr . However, the imaginary parts calculated there gave accurate imaginary parts only in the sliding regime of larger negative $`g`$ and did not invade into the proper tunnelling regime of small $`g`$ dominated by critical bubbles. A separate variational treatment of an instanton calculation was set up to cover this region Tunn .

In this note we show that the non-Borel-summable series can be evaluated with any desired precision by an appropriate continuation of variational perturbation theory Hagen ; Verena  in to the complex $`g`$-plane.

Variational perturbation theory has a long history refs ; finiteg ; KleinertJanke ; Guida . It is based on the introduction of a dummy variational parameter $`\mathrm{\Omega }`$ on which the full perturbation expansion does not depend, while the truncated expansion does. An optimal $`\mathrm{\Omega }`$ is then selected by the principle of minimal sensitivity Stevenson , requiring the quantity of interest to be stationary as function of the variational parameter. The optimal $`\mathrm{\Omega }`$ is usually taken from a zero of the derivative with respect to $`\mathrm{\Omega }`$. If the first derivative has no zero, a zero of the second derivative is chosen. For Borel-summable series, these zeros are always real, in contrast to statements in the literature Neveu ; RULES ; Braaten ; Ramos  advocating the use of complex zeros. Complex zeros, however, produce in general wrong results for Borel-summable series, as was recently demonstrated in Ref. AntiBraaten .

The purpose of this note is to show that there does exist a wide range of applications of complex zeros in the previously untreatable field of non-Borel-summable series. These arise typically in tunnelling problems, and we shall see that variational perturbation theory provides us with an efficient method for evaluating these series and rendering their real and imaginary parts with any desirable accuracy if only enough perturbation coefficients are available. The choice of the complex zeros is dictated by the requirement to achieve at each order the least oscillating imaginary part when approaching the tip of the cut. We shall call this selection rule the principle of minimal sensitivity and oscillations.

2. For an introduction to the method consider the exactly known partition function of an anharmonic oscillator in zero spacetime dimensions, which is a simple integral representation of a modified Bessel function $`K_\nu (z)`$:

$`Z(g)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x\mathrm{exp}(x^2/2gx^4/4)`$ (1)
$`=`$ $`e^{1/8g}(4\pi g)^{1/2}K_{1/4}(1/8g),`$
For small $`g<0`$, the function $`Z(g)`$ and its inverse $`D(g)Z^1(g)`$ have a divergent non-Borel-summable power series:

$$D(g)=\underset{l=0}{\overset{\mathrm{}}{}}a_lg^l,$$
(2)
In the strong-coupling regime, there exists a convergent expansion

$`D(g)=g^\alpha {\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}b_lg^{\omega l},`$ (3)
with $`\alpha =1/4`$ and $`\omega =1/2`$. In the context of critical phenomena, the exponent $`\omega `$ coincides with the the Wegner exponent Wegner  of approach to scaling strong . The $`L`$th variational approximation depending on the variational parameter $`\mathrm{\Omega }`$ is given by the truncated series Hagen ; Verena

$`D_{\mathrm{var}}^{(L)}(g,\mathrm{\Omega })=\mathrm{\Omega }^p{\displaystyle \underset{j=0}{\overset{L}{}}}\left({\displaystyle \frac{g}{\mathrm{\Omega }^q}}\right)^jϵ_j(\sigma ),`$ (4)
where $`q=2/\omega =4`$, $`p=\alpha q=1`$. Introducing the parameter $`\sigma =\mathrm{\Omega }^{q2}(\mathrm{\Omega }^21)/g`$, the re-expansion coefficients are

$`ϵ_j(\sigma )={\displaystyle \underset{l=0}{\overset{j}{}}}a_l\left({\displaystyle \genfrac{}{}{0pt}{}{(plq)/2}{jl}}\right)(\sigma )^{jl}.`$ (5)
Following the principle of minimal sensitivity, we have to find the zeros of the derivative of $`_\mathrm{\Omega }D_{\mathrm{var}}^{(L)}(g,\mathrm{\Omega })`$, which happen to coincide with the zeros of a function of the variable $`\sigma `$ only, to be denoted by $`\zeta ^{(L)}(\sigma )`$:

$`\zeta ^{(L)}(\sigma )={\displaystyle \underset{l=0}{\overset{L}{}}}a_l\left({\displaystyle \genfrac{}{}{0pt}{}{(plq)/2}{Ll}}\right)(plq+2l2L)(\sigma )^{Ll}.`$ (6)
For a proof of this remarkable property see JKsig , and the textbook Hagen  (p. 291).

As an example, we take the weak-coupling expansion to $`L=16`$th order and calculate real and imaginary parts for the non-Borel region $`2<g<0.008`$ selecting the zero of $`\zeta ^{(16)}(\sigma )`$ according to the principle of minimal sensitivity and oscillations. The result is shown in Fig. 1. In order to point out how well the variational result approximates the essential singularity of the imaginary part $`\mathrm{exp}(1/4g)`$ for small $`g<0`$, we have removed this factor. The agreement is excellent down to very small $`g`$.

3. Let us now turn to the nontrivial problem of summing the instanton region of the anharmonic oscillator for $`g<0`$. The divergent weak-coupling perturbation expansion for the ground state energy of the anharmonic oscillator with the potential $`V(x)=x^2/2+gx^4`$ is, to order $`L`$;

$`E_{0,\mathrm{weak}}^{(L)}(g)={\displaystyle \underset{l=0}{\overset{L}{}}}a_lg^l,`$ (7)
where $`a_l=(1/2,3/4,21/8,333/16,30885/128,\mathrm{})`$. The expansion is obviously not Borel-summable for $`g<0`$, but will now be evaluated with our new technique, proceeding in the same way as for the above test function $`D(g)`$ via Eqs. (4) through (6). The known strong-coupling growth parameters are $`\alpha =1/3`$ and $`\omega =2/3`$, so that $`p=1`$ and $`q=3`$ in Eq. (6) which will guarantee the correct scaling properties for $`g\mathrm{}`$. To order $`L=64`$ we obtain from the optimal zero of $`\zeta ^{(64)}(\sigma )`$ the logarithm of the scaled imaginary part

$`l(g):=\mathrm{log}\left[\sqrt{\pi g/2}E_{0,\mathrm{var}}^{(64)}(g)\right]1/3g,`$ (8)
shown in Fig. 2 for $`0.2<g<0.006`$. The point $`g=0.006`$ is the closest approach to the singularity at $`g=0`$ for $`L=64`$ before the onset of oscillations.

Let us compare our curve with the expansion derived from instanton calculations ZINNJ :

$`f(g)`$ $`=b_1gb_2g^2+b_3g^3b_4g^4+\mathrm{},`$ (9)
with coefficients $`b_1=3.95833,b_2=19.344,b_3=174.21,b_4=2177`$. This expansion is divergent and non-Borel-summable for $`g<0`$. Remarkably, we are able to extract this expansion from our data points. Since our result is given by a convergent expansion, the fitting procedure will depend somewhat on the interval chosen over which we fit our curve by a power series. A compromise between a sufficiently long interval and the runaway of the divergent instanton expansion is obtained for a lower limit $`g>.0229\pm .0003`$ and an upper limit $`g=0.006`$. Fitting a polynomial to the data, we extract the following first three coefficients:

$`b_1=3.9586\pm .0003,b_2=19.4\pm .12,b_3=135\pm 18.`$ (10)
4. The agreement of our curves with the exact results in Figs. 1 and of our expansion coefficients in (10) with the exact ones in (9) demonstrates that our method is capable of probing deeply into the instanton region of the coupling constant.

Let us end by remarking that another procedure of summing non-Borel series can be deduced from a development in the first of Refs. strong  (see also Chapter 20 of the textbook Verena ) One may derive a strong-coupling expansion of the type (3) from the divergent weak-coupling expansion, which can then be continued analytically to negative $`g`$ by a simple rotation of the power $`g^{\omega l}`$ to $`e^{i\pi \omega l}(g)^{\omega l}`$. This method is applicable only in the sliding regime. In Fig. 2, we have plotted the resulting curve to order $`L=9`$. The present work fills the gap between the sliding regime and the tunnelling regime by extending variational perturbation theory to all $`g`$ arbitrarily close to zero, without the need for a separate treatment of the tunnelling regime.

There exists, of course, a wealth of possible applications of this simple theory, in particular to quantum field theory where variational perturbation theory has so far yielded the most accurate critical exponents from Borel-summable series strong ; Verena ; Lipa .

Acknowledgment
This work was partially supported by ESF COSLAB Program and by the Deutsche Forschungsgemeinschaft under Grant Kl-256.


text:
# Substructure in the Andromeda Galaxy Globular Cluster System

## 1 Introduction

In hierarchical structure formation scenarios, galaxies form at least in part from the merger and accretion of smaller clumps of matter (e.g., Searle & Zinn, 1978; Cole et al., 2000). Clumping in the early universe is believed to have resulted from small perturbations in the primordial density field, instabilities that triggered the growth of structure from the “bottom up”. Such primordial clumps would then evolve into dwarf galaxies like those observed in the present epoch, or may merge to form more massive galaxies.

There are indications that such mergers are ongoing, as is evident from observations of the Sagittarius dwarf galaxy’s current interaction with the Milky Way (Ibata, Gilmore & Irwin, 1995; Dohm-Palmer et al., 2001; Newberg et al., 2002). There are also indications of tidal distortions in M31’s dwarf satellite NGC 205 (Walterbos & Kennicutt, 1987; Bender, Paquet & Nieto, 1991). Furthermore, a wide-field imaging survey of M31 by Ibata et al. (2001) has revealed the presence of a giant stream of metal-rich stars within the halo of M31. This feature may have originated from tidal interactions with M31’s dwarf companions, NGC 205 and M32. Either the process of satellite accretion is not uncommon, or we are observing our Galaxy and M31 at particularly unique stages of their evolution. Recent studies by Reitzel & Guhathakurta (2002) and Ferguson et al. (2002) have also found signs of substructure and tidal trails in the M31 halo which support the possibility of past accretion events.

Globular clusters (GCs) have been identified within some of the dwarf satellites of the Milky Way (Da Costa & Armandroff, 1995; Hodge et al., 1999; Oh, Lin & Richer, 2000) and M31 (Grebel, Dolphin & Guhathakurta, 2000), as well as in the dwarf members of distant galaxy clusters (e.g. Durrell et al., 1996; Miller et al., 1998). Based on their positions, distances and radial velocities, there are four Milky Way globular clusters which are believed to be associated with the Sagittarius dwarf (Ibata, Gilmore & Irwin, 1994; Ibata et al., 1997; Bellazzini, Ferraro, & Ibata, 2003). Ultimately, these clusters — M54, Arp 2, Terzan 7 and Terzan 8 — will have their orbits randomized by the strong potential of the Milky Way and they will presumably take their place as bona fide members of the Galactic globular cluster system.

Therefore, in the context of hierarchical galaxy formation, one would expect that the globular cluster systems (GCSs) of massive galaxies could include the cannibalized globular clusters of past-accreted dwarfs. For instance, this seems to be a plausible mechanism for building up the GCSs of giant elliptical galaxies (e.g. Côté, Marzke, & West, 1998). Such a scenario raises the question of whether the globular cluster systems of merged dwarfs would remain distinct dynamical entities or instead become hopelessly mixed with the larger galaxy’s GCS. If they remain recognizable over significant time-scales, satellite GCS remnants could provide us with a powerful dynamical probe of the formation and evolution of galaxies.

The first significant effort at distinguishing sub-groups within the GCS of an external galaxy was performed for the M31 system by Ashman & Bird (1993). These authors suggested that an apparent sub-clustering in the positions and velocities of M31 globular clusters represents a surviving signature of progenitor gas clouds from which the galactic halo itself may have formed. Ashman & Bird employed a technique in which globular clusters were grouped based on deviations from global mean velocities and the velocity dispersions between each cluster and its $`N`$ nearest neighbors. Their sample included 144 M31 GCs with spectroscopic data from Huchra, Brodie & Kent (1991), many of which had large velocity uncertainties ($`\pm 50`$ km s-1).

In many respects, the M31 globular cluster system provides an ideal target in the search for substructure and grouping. More than 435 confirmed globular clusters have been identified in M31 (Barmby et al., 2000). A recent spectroscopic study of the M31 GCS by Perrett et al. (2002) has now pushed the number of clusters with radial velocities and metallicities beyond 300. In addition, Perrett et al. (2002) quote typical velocity uncertainties for the majority of the clusters of $`\pm 12`$ km s-1, a significant improvement over many previous studies. The time is right to take a closer look for signatures of hierarchical formation within the M31 globular cluster system.

There are several important things to consider in this type of study. One must define what is to be regarded a significant “grouping” in the available parameter space. In other words, appropriate grouping criteria must be established a priori. Furthermore, one must investigate the likelihood that any observed groups may be mere coincidental associations rather than actual substructure. A greater difficulty is the interpretation of the significance of the grouping results: do the groups represent true merger remnants that have survived to the present epoch?

In this paper, we develop a technique to identify potential dynamical remnants within a galaxy’s globular cluster system. This group-finding algorithm is described in Section 2. In Section 3, we apply the group-finding technique to the sample of M31 GCS positions, velocities and metallicities from Perrett et al. (2002). We identify several potential groups of clusters in this sample. Monte Carlo simulations are performed on randomized globular cluster populations to characterize the significance of the identified groups. Section 4 describes our use of N-body simulations to investigate whether these groups could represent actual merger remnants. Finally, the results of this study are summarized in Section 5.

## 2 The Technique of Group Identification

In most applications, simulated data has the advantage over observational data in that it provides full three-dimensional positional and dynamical information. Observational data is, of course, less ideal: we must cope with projected positions and line-of-sight velocities and we can, at best, make generalized assumptions about the true dynamical behavior of any system. Additionally, every observed parameter carries with it an uncertainty, the value for which may vary from object to object even within the same dataset. These are some of the considerations that must be made when devising a group-finding algorithm for observational data.

### 2.1 Group-Finding Algorithm

The basic technique used in our group identification is a friends-of-friends algorithm. This algorithm links together particles that fall within a specified distance in parameter space. Each distinct set of joined particles constitutes a group. In the case of the M31 globular cluster system observations, the parameter space is defined by $`X`$ and $`Y`$ (projected) positions, radial velocities $`v_r`$, and metallicity \[Fe/H\].

The appropriate weighting of these parameters must be determined in advance in order to set the group-linking criteria. Since each parameter has units, each must be converted to an equivalent unitless quantity on the same scale before calculating distances in parameter space. For example, if $`X`$ and $`Y`$ each span a range of $`100`$ arcmin, the velocities range over 600 km s-1, and metallicity spans only 3 dex, an appropriate metric must be devised such that the group is not defined principally based on \[Fe/H\] owing to its comparatively small numerical range. In essence, the metric puts all of the parameters on to the same effective scale. In addition to the metric, the group-finding algorithm requires a linking length which defines the maximum allowable distance between successive group members in the specified parameter space.

Therefore, for each GC, the algorithm links a nearby GC to it if

$$(\mathrm{\Delta }s)^2\underset{i}{}\gamma _i^2(\mathrm{\Delta }x_i)^2<^2,$$
(1)
where $`\mathrm{\Delta }s`$ is the total separation (in parameter space) of the neighbor to the GC in question, $`\gamma _i`$ is the scaling or weighting factor defined for parameter $`i`$, $`\mathrm{\Delta }x_i`$ is the separation between the GC and its neighbor in parameter $`i`$, and $``$ is the linking length. Traditionally, $``$ is specified as a fraction, $`\mathrm{\Lambda }`$, of the mean inter-particle spacing for the entire system:

$$=\mathrm{\Lambda }\overline{\mathrm{\Delta }s}$$
(2)
The advantage of using a friends-of-friends technique over average density methods is that it can follow elongated structures, features that may be expected in tidal disruptions.

The error bars associated with each input parameter are accommodated by the addition of an “uncertainty ellipsoid” surrounding each GC in the available parameter space. A GC’s parameter-space neighbor is linked with it if the minimum distance between the uncertainty ellipsoids is less than the specified linking length.

### 2.2 Metric and Linking-Length Determinations

In order to define a suitable metric to be used in the definition of a group, a method was devised to test the group-finding algorithm on simulated associations using a range of scalings and linking lengths. Groups of objects with various population sizes, spatial distributions (e.g. spherical, linear and curved streams), galactocentric positions and orientations were created with general parameters similar to those of Local Group dwarf galaxies. The simulated groups were given internal velocity dispersions of $`\sigma _v=6`$ to $`12`$ km s-1, mean metallicities of $``$\[Fe/H\]$`=1.2`$ to $`1.6`$ dex, and \[Fe/H\] spreads of $`\sigma _{[\mathrm{Fe}/\mathrm{H}]}=0.3`$ to $`0.5`$ dex (Mateo, 1998; Da Costa et al., 2000).

The GCs in these simulated groups were then superimposed over a background designed to mimic a realistic distribution of M31 globular clusters. This background of several hundred GCs consisted of the following components:

* a random, metal-poor, spherical halo distribution,
* a rotating, metal-rich thick disk distribution inclined at $`12.3^{}`$ to the line-of-sight, and
* a rotating, intermediate-metallicity ellipsoidal bulge distribution.

Projected positions, radial velocities and metallicities for the fabricated group members and additional background objects were analysed using the group-finding algorithm. The procedure was cycled through a range of linking lengths and parameter scalings to determine the optimal combination for group retrieval. A minimum group membership of $`n_{GC}=4`$ clusters was specified to define a group.

Setting the $`X`$ and $`Y`$ scaling to $`\gamma _X=\gamma _Y=1.0`$ arcmin-1 to define the baseline, it was found that the best recovery rate occurred when $`\mathrm{\Lambda }=0.1`$ and the velocities were scaled by $`\gamma _v=0.13`$ km$`^1`$s for cases when metallicities were not considered in the group recovery. With the addition of metallicity information, the optimal velocity and \[Fe/H\] scalings were found to be $`\gamma _v=0.14`$ km$`^1`$s and $`\gamma _{[\mathrm{Fe}/\mathrm{H}]}=9.7`$ dex-1, respectively, with a linking length given by $`\mathrm{\Lambda }=0.09`$. These optimal linking length and parameter scaling values were used to search for grouping in the globular cluster system of M31.

## 3 Substructure in the M31 GCS

### 3.1 The GCS Sample

The best available spectroscopic sample of $`300`$ GC positions, radial velocities and metallicities as defined in Perrett et al. (2002) was examined for signs of grouping. Quoted errors on individual velocities and metallicities were used to define the uncertainty ellipsoids about each data point in parameter space. Uncertainties on the target positions were typically around 0.2 – 1.0 arcseconds; these were deemed to be negligible and thus position uncertainties were not incorporated in the analysis.

The group-finding algorithm was used to analyze the data based on the parameter scalings and linking lengths determined in Section 2.2. Again, a minimum of four clusters was required in order to define a group. A search for groups was conducted in position/velocity/metallicity parameter space for the sample of 301 globular clusters with published spectroscopic metallicities. Roughly one dozen groups were initially flagged within the GC population. The robustness of these groups was then tested against small ($`10\%`$) changes in the specified parameter scalings and linking length. Groups that were not found to be stable against these perturbations were rejected.

At this point it should be noted that the inclusion of metallicity information may not necessarily improve the effectiveness of the group-finding results. The \[Fe/H\] uncertainties for many of the M31 globular clusters remain quite large ($`\pm 0.30.9`$ dex). This reduces the significance of metallicity as a linking parameter due to the increased limits of each uncertainty ellipsoid in this dimension. Furthermore, it has been demonstrated that the dwarf galaxies of the Local Group exhibit a wide variety of star formation histories and metallicities (see Mateo, 1998, and references therein). For example, Grebel & Guhathakurta (1999) report a large ($`1`$ dex) spread in stellar metallicities within the M31 satellites And VI and And VII. Similar spreads have been found in the globular cluster systems of some dwarf elliptical galaxies (J. Lotz, private communication). As a result, proximity in metallicity is not necessarily a very suitable criterion for group identification within a hierarchical formation scenario.

For these reasons, the complete available sample of 321 M31 globular cluster positions and velocities was also examined for grouping without the inclusion of \[Fe/H\] information. Again, roughly a dozen potential groups were identified and then tested for stability against small changes in the grouping criteria. The majority of the resulting groups correspond to those found when metallicity was included as a linking parameter.

Perrett et al. (2002) identified significant bimodality in their sample of M31 GC metallicities, in confirmation of earlier reports made by Huchra, Brodie & Kent (1991) and Ashman & Bird (1993). The centrally concentrated, rapidly rotating, metal-rich cluster system is consistent with a bulge population (Perrett et al., 2002). The group-finding results in the full GCS data sample both with and without metallicities revealed a couple of relatively large groups that correspond to a sizable fraction of this metal-rich bulge population. A third grouping analysis was therefore also performed on the sample of 231 GCs in the metal-poor population as defined by Perrett et al. (2002). Metallicity was not used as a linking parameter in this test. Several tentative groups persisted within the central region even within this metal-poor population, yet the memberships of all but one of these groups were not robust to small variations in the linking criteria. Ultimately, groups with average galactocentric radii within the central $`10\mathrm{}`$ of the galaxy were rejected due to the higher probability of false detections, an effect which is confirmed by the Monte-Carlo simulations presented later in Section 3.3. Additionally, any cannibalized dwarf GCSs would be extremely difficult to recover within the dense inner regions of the galaxy; this effect will also be noted in the numerical simulation results discussed in Section 4.

### 3.2 Grouping Results

In the M31 GCS, there is evidence for the presence of ten unique globular cluster groups of $`n_{GC}4`$ members with mean galactocentric radii beyond $`10\mathrm{}`$. The positions of these groups are shown in Figure 1, and the data for the group members are provided in Table 1. Five of these groups were identified in all three search cases described above. Three additional groups (Groups 7, 8 and 10) resulted from a search of the full sample of 321 globular clusters without the inclusion of \[Fe/H\] information. Two further groups can be considered as tentative identifications: Group 3 was found only when metallicity was included as a linking parameter, and Group 4 resulted only from the search of the metal-poor population. All of the groups listed in Table 1 are stable to small perturbations in the linking criteria.

An eleventh group is also listed in Table 1; despite the fact that this group lies close to the galaxy center, the mean radial velocity of its group members ($`155`$ km s-1) is considerably higher than the systemic velocity of M31 ($`300`$ km s-1). This group deviates significantly from the GCS rotation curve determined by Perrett et al. (2002), and in fact is moving counter to the general rotation at that major axis position.

The association of a group with NGC 205 is clearly evident in Figure 1 (Group 5 of Table 1). The mean velocity of the NGC 205 group is $`v_r=259`$ km s-1 with $`\sigma _v=78`$ km s-1, a result that is consistent with NGC 205’s systemic velocity of $`242\pm 3`$ km s-1 (de Vaucouleurs et al., 1991).

It is interesting to note that none of the groups recovered in this analysis overlap with the associations found by Ashman & Bird (1993). This may be partly attributable to the larger sample of higher precision velocities used here. Ashman & Bird (1993) also indicated that their method may produce false hits in the presence of GCS rotation. In an effort to reduce this affect, they applied a KMM mixture-modelling technique to remove the rotating metal-rich (disk/bulge) population from their grouping analysis, thus decreasing their sample from 144 GCs to 95 “halo” objects (Huchra, Brodie & Kent, 1991). However, Perrett et al. (2002) found that even the metal-poor cluster population demonstrates significant rotation, a factor which may also have biased the Ashman & Bird (1993) grouping results.

The majority of the clusters in the inner and outer subsystems found by Saito & Iye (2000) are included in the metal-rich (bulge) population inferred from the KMM results of (Perrett et al., 2002). The bulk of inner subsystem found by Saito & Iye (2000) was recovered here as two central groups; the memberships of these inner groups were not stable and these are not listed in Table 1 for reasons discussed previously.

### 3.3 Monte-Carlo Simulations

A series of Monte-Carlo simulations were performed to investigate the likelihood that the identified M31 globular cluster groups are mere chance associations on the sky. Randomized populations of 321 globular clusters were generated based on the globular cluster positions and radial velocities from Perrett et al. (2002). Metallicities were not included in the randomized tests due to the large \[Fe/H\] uncertainties in the data. Furthermore, the eight non-tentative outer groups in our sample were found without incorporating any metallicity information in the group finding.

First, the cumulative distribution observed in projected position along the major axis of the galaxy was used to select a random cluster X position with roughly the same underlying distribution as that seen in the M31 GCS. Next, the cumulative distribution in projected minor axis location of the 70 nearest objects in X was used to generate a randomized (but realistic) value of Y position for each target. Similarly, the cumulative distribution in radial velocity for the 70 nearest objects in X was used to generate a randomized $`v_r`$ for each target. This procedure generated a population of globular clusters which shares similar observed properties with the observed data — including spatial distribution and rotation profile — while eliminating any small-scale correlations that may exist in the dataset.

One thousand realizations of the randomized GC populations were performed and each was tested for grouping in position and velocity parameter space. The results of these Monte Carlo tests are shown in Figure 2. The top panel of Figure 2 shows the radial distribution of the groups found in the randomized data. Roughly half of these “false hits” were recovered within the central $`10\mathrm{}`$ of the galaxy, a region which appears to be too dense to yield reliable group detections. The fact that the spatially compact (metal-rich) GC population demonstrates rapid rotation also adds to the likelihood of detecting chance conglomerations in position-velocity parameter space. Therefore, neglecting all of the groups with mean radii within $`10\mathrm{}`$ of the galaxy center, we obtain the number distribution of random-data groups shown in the lower panel of Figure 2.

If we exclude the two tentative groups and the central group, the remaining eight groups observed in the actual M31 dataset (Table 1) lie above $`97\%`$ of the randomized runs. This suggests that the substructure found in the M31 GCS is significant overall. If we further exclude the group associated with NGC 205, this result still lies above $`91\%`$ of the randomized runs. Therefore, it is unlikely that all of the groups are chance associations, although the majority may indeed be. For example, there is nearly a $`60\%`$ chance that three or more of the groups are not real. With the currently available data, it is not possible to determine if the individual groups are dynamically bound structures, nor can we gauge the specific degree of internal contamination by background clusters.

## 4 Modelling the Fate of the Globular Cluster Systems

Based on the discovery of apparent substructure in the M31 GCS, it is worthwhile to investigate the plausibility that the proposed groups are a result of merger events. To accomplish this, it is instructive to consider the fate of a dwarf satellite’s small globular cluster population during a merger with a massive galaxy. N-body simulations were performed to follow the evolution of a dwarf galaxy GCS over timescales of $`\tau 1`$ Gyr as it interacts with an external potential representing a massive spiral galaxy. The form of the galactic potential and the modelling of the dwarfs are discussed in Sections 4.1 and 4.2, respectively. The interpretation of the simulation results is presented in Section 4.3.

### 4.1 The Galactic Potential

In this study, the galactic potential was assumed to be rigid and time-independent, and was built from three components: a spheroidal bulge potential (Hernquist, 1990), a Miyamoto & Nagai (1975) thick disk, and a logarithmic halo. The forms of these potentials were based on the modelling of the Sagittarius dwarf by Johnston, Spergel & Hernquist (1995) and Helmi & White (2001). The potentials are:

$`\varphi _{\mathrm{disk}}`$ $`=`$ $`{\displaystyle \frac{GM_{\mathrm{disk}}}{\sqrt{R^2+(a+\sqrt{z^2+b^2})^2}}}`$ (3)
$`\varphi _{\mathrm{bulge}}`$ $`=`$ $`{\displaystyle \frac{GM_{\mathrm{bulge}}}{(r+c)}}`$ (4)
$`\varphi _{\mathrm{halo}}`$ $`=`$ $`v_{\mathrm{halo}}^2\mathrm{ln}(r^2+d^2)`$ (5)
where $`G`$ is the gravitational constant, $`r`$ and $`R`$ represent radii in spherical and cylindrical coordinates, respectively, and $`z`$ is the height above the galactic plane. The masses ($`M`$), scale lengths ($`a`$, $`b`$, $`c`$, $`d`$) and velocity ($`v_{\mathrm{halo}}`$) in the above potentials were selected to reproduce the characteristics of the Milky Way. The components of our Galaxy are better understood than those of M31 and we can expect the overall effects of satellite disruption within M31 to be quite similar. The adopted component parameters are provided in Table 2 and have been selected following the work of Johnston (1998).

The representation of the massive galaxy by a fixed potential neglects any exchange of energy between the satellite and its parent galaxy. However, typical masses of dwarf satellites are several orders of magnitude smaller than that of either the Milky Way or M31. As a result, the effects of dynamical friction and energy exchange between the satellite and the massive galaxy can be neglected over the timescales of interest here (Helmi & White, 2001; Binney & Tremaine, 1987).

### 4.2 The Model Dwarf

The dwarf satellite galaxy in each simulation was modelled as a Hernquist sphere with a density distribution given by

$$\rho (r)=\frac{M_{\mathrm{sat}}}{2\pi }\frac{s}{r}\frac{1}{(r+s)^3}$$
(6)
(Hernquist, 1990), which corresponds to a cumulative mass distribution of

$$M(r)=M_{\mathrm{sat}}\frac{r^2}{(r+s)^2}.$$
(7)
$`M_{\mathrm{sat}}`$ is the total satellite mass and $`s`$ is a scale length that defines the half-mass radius of the distribution: setting $`M(r_{1/2})=M_{\mathrm{sat}}/2`$ in Equation 7 gives $`r_{1/2}=(1+\sqrt{2})s`$. This analytical mass model is an approximation of the de Vaucouleurs $`R^{1/4}`$ profile applicable to galaxies and bulges (de Vaucouleurs, 1953).

For each satellite, roughly $`\mathrm{10\hspace{0.17em}000}`$ particles were distributed according to this density profile. Dwarf galaxies of several masses were considered within the range from $`10^7M_{}`$ to $`10^9M_{}`$. Ten of the particles were given a somewhat more extended distribution about the dwarf and were designated to be its globular clusters; these GCs were each assigned masses of $`5\times 10^5M_{}`$. The remaining mass of the simulated dwarfs was equally divided among the rest of the particles. The number of particles used to represent the dwarf was sufficient to model its disruption, while being small enough to efficiently run the required number of simulations.

A series of $`85`$ simulations were run in order to explore a reasonable range of parameter space by altering the initial conditions of the incident satellite. Table 3 lists the ranges of initial conditions that were tested. The total dwarf masses ($`M_{\mathrm{sat}}`$) were selected to span a range between a typical dwarf spheroidal and a dwarf elliptical and include the mass of its 10 globular clusters. The range in satellite half-mass radii ($`R_{1/2}`$) was chosen to correspond to typical values of dwarf galaxy core radii — a reasonable assumption if light follows mass. $`R_{\mathrm{GC}}`$ is the scale radius of the dwarf’s globular cluster system. The satellites were assigned a starting galactocentric distance ($`D`$) and a position angle ($`\theta `$) above the plane of the disk. The azimuthal angle was taken to be $`\varphi =0`$, since the potential is axisymmetric. The parameters $`\alpha `$ and $`\beta `$ dictate the initial partition of velocity in the $`XZ`$ and $`XY`$ plane, respectively, in units of circular velocity such that $`\alpha ^2+\beta ^2=1`$ for circular orbits.

These N-body simulations were performed using the Barnes-Hut treecode111Treecode version 1.4 distributed by J. Barnes. and modified to include an external potential. The treecode algorithm is based on the commonly-used hierarchical force method described by Barnes & Hut (1986).

It should be noted that the simulations presented herein are intended to be illustrative, and not to reproduce the detailed characteristics of this complex interaction. The addition of evolving potentials for the bulge, disk and halo components and the incorporation of effects such as dynamical friction, disk-shocking, and gas dynamics would permit a more comprehensive investigation of the merger event itself.

### 4.3 Grouping in the Simulation Results

Using the scaling and linking lengths determined in the previous step ($`\mathrm{\Lambda }=0.1`$, $`\gamma _v=0.13`$ km$`^1`$s), the N-body simulation results were analysed to investigate the likelihood of identifying the cannibalized dwarf’s globular clusters as a group in parameter space. The three-dimensional positions and velocities were projected into an orientation similar to M31 ($`i=77.7^{}`$), and were superimposed on the background components described in Section 2.2. The group-finding algorithm was then applied to the simulation results at various timesteps.

Two parameters, $`Q`$ and $`Q_{max}`$, were defined in order to quantify the retrieval efficiency of the satellite’s GCS:

$$Q=\frac{1}{N}\underset{k}{}\left(\frac{n_{GC}}{n_g}\right)_k$$
(8)
$$Q_{max}=\mathrm{max}\left(\frac{n_{GC}}{n_g}\right)_k$$
(9)
where $`N`$ is the total number of groups containing at least one of the dwarf’s original globular clusters, $`n_{GC}`$ is the number of the globular clusters in group $`k`$, and $`n_g`$ is the total number of objects (including background) that were identified in group $`k`$. The sum in $`Q`$ is calculated only over those groups that contain one or more globular clusters and reflects the general contamination of the GC groups by background objects. $`Q_{max}`$ is the maximum ratio of the number of satellite GCs in a group to the total number of objects allocated to that group, and therefore indicates whether a bona fide group has been found. For perfect group recovery with no contamination, $`Q=Q_{max}=1`$.

To demonstrate, let us examine a particular simulation with an initial dwarf position angle $`90^{}`$ ($`\alpha =0.9`$, $`\beta =0`$) with $`M_{\mathrm{sat}}=10^9M_{}`$ at a timestep of 1.2 Gyr. The group-finding algorithm returned globular cluster group memberships as follows: 6 of the 10 dSph GCs were in the same group and had 1 additional background object included, 1 GC was allocated to a background group with 3 additional members, 1 GC was allocated to a background group with 4 additional members, and 2 GCs were not identified as members of any group. For this case with $`N=3`$ groups containing one or more globular clusters:

$$Q=\frac{1}{3}\left[\frac{6}{7}+\frac{1}{4}+\frac{1}{5}\right]=0.44$$

$$Q_{max}=\mathrm{max}(\frac{6}{7},\frac{1}{4},\frac{1}{5})=0.86$$

The high $`Q_{max}`$ indicates that it is very likely that we can detect a group that contains a significant fraction of the disrupted satellite’s GCS, here with only $`14\%`$ contamination by background objects. The moderate to low $`Q`$ ratio reflects the fact that some of the satellite’s GCs were allocated to other groups to which they do not truly belong.

The orbital periods of the satellite debris within the potential of the parent galaxy are of order 1 Gyr. The probability of detection clearly depends to a large degree on the relic’s location in its orbit, as there is a decreased likelihood of finding a group if the GCS debris is fore or aft of the bulge, or obscured in the rotating disk particles. Therefore, the likelihood that the detection of a group in a galaxy’s GCS reveals the remnant of a real accretion event is comparatively low if the group is superimposed over the dense part of the disk or the bulge of the galaxy. Not surprisingly, satellite GCS remnants are also most easily observable at large galactocentric radii due to their significantly lower degree of disruption over the relevant timescales.

In reality, we rarely know the orbital parameters of past (or even present) satellite galaxies. We can therefore examine the average $`Q`$ values of a large sample of the simulations having a variety of initial position angles and incident velocities to get some idea of the likelihood of group detection for a satellite on a random orbit. The mean and median $`Q`$ and $`Q_{max}`$ values for satellites with an initial galactocentric radius of $`40`$ kpc are shown in Figure 3. The average and median $`Q`$ values decrease with time, as expected. The $`Q`$ results indicate that it is not likely that we can recover the complete GCS of the satellite dwarf in a single uncontaminated group. However, the $`Q_{max}`$ plot shows that one significant grouping of clusters associated with the cannibalized dwarf can generally be found up to several Gyr after the initial encounter.

## 5 Conclusions

The numerical simulation results presented in this study reveal that it is possible to identify some fraction of a disrupted satellite’s initial GCS based on grouping in projected position and radial velocity. This can be accomplished for timescales of $`\tau \mathrm{\hspace{0.17em}1}`$ Gyr after the initial encounter, although generally not without some significant degree of background contamination in the group membership.

In the Milky Way, there are three-dimensional data available that trace a distinct stellar stream associated with the Sagittarius dwarf merger. It is then relatively easy to pick out any globular clusters that coincide with this stream in both distance and velocity, although even here the results are somewhat controversial (for example, see Minniti, Meylan & Kissler-Patig, 1996). At the distance of M31, we have projected data and are only now beginning to uncover signs of stellar streams and substructure (Ibata et al., 2001). These observed streams extend well beyond the spatial limits of the current cluster sample (on a line connecting NGC 205 and M32), and thus cannot be directly matched with any of the GC groups identified in this study. Large-area surveys (e.g., Massey, 2002; Lee et al., 2002) will help to provide the necessary additional stellar astrometry and to reveal any new GC candidates in M31. Nonetheless, Lynden-Bell & Lynden-Bell (1995) demonstrated that it is possible to use streams of globular clusters to trace the orbits of potential merger remnants within the Milky Way. Streams or groups of globular clusters in nearby massive galaxies may also provide a viable means of identifying the orbits of any recently merged satellites.

We have found significant signs of grouping in projected position and radial velocity within the globular cluster system of M31. Do these M31 globular cluster groups represent the signatures of past accretion events? Without accurate distances or dynamics, the correlation of M31 GCS groups with past accretion events remains tenuous. In followup studies, the identification of a stellar stream coincident with a grouping of globular clusters would lend the necessary credibility to the notion that a galaxy’s GCS — and by extension the galaxy itself — formed at least partly in a hierarchical manner.

## 6 Acknowledgments

The authors wish to thank Larry Widrow and Geraint Lewis for useful comments and related discussions. The authors are also grateful to Josh Barnes for making his treecode software available, and to the anonymous referee for the helpful suggestions. This work was supported in part by funding from the Natural Sciences and Engineering Research Council of Canada.