Production of ultra-energetic cosmic rays through the decay of super-heavy X particles

We present a new and complete numerical analysis of the decay of super-heavy X particles, assumed to be the origin of cosmic rays with energy beyond the GZK cut-off. The decay of X initiates a ``parton shower'', where we include all degrees of freedom contained in the Minimal Supersymmetric Standard Model (MSSM). Since at energies near $M_X$ all gauge couplings are of similar magnitude, we include all of them, as well as third generation Yukawa couplings. Technically the shower development is described through the DGLAP evolution of the relevant fragmentation functions (FFs). We also carefully treat the decay of the superparticles as well as heavy SM particles created in the shower. Nonperturbative physics is parameterized through the input values of the FFs, which we take from the literature. The final result is the the complete spectrum of all stable particles at the very end of the shower : protons, electrons, neutrinos, photons and neutralinos, for an energy range from $10^{-7} M_X$ to $M_X$. In particular, the flux of high-energy neutralinos is sizable ; it might serve as ``smoking gun'' signature for this kind of scenario.


Introduction
The origin of the Ultra Energetic Cosmic Rays (UHECR) is still an enigma for scientists. We neither know how and where they are produced, nor what they are. Photons with energy E γ ≥ 10 15 eV can be absorbed on the cosmic microwave background (CMB) through e + e − pair production; at E γ ∼ 10 20 eV, the main energy loss comes through interactions with the radio background. Electrons in addition loose energy through synchrotron radiation in intergalactic and galactic magnetic fields, the strengths of which are not known very well. Protons and heavy nuclei of energies above a few times 10 19 eV loose energy through interactions with the CMB; one would thus expect their spectrum to cut off at that energy (the so called GZK cut-off [1], [2]). However, a few events have been observed with energies beyond this cut-off [3]. This means that the UHECR should either have been created within a few interaction lengths from the Earth, i.e. at a distance < ∼ 100 Mpc, or they should have been produced at even higher energies than are observed now. The first possibility is in principle consistent with the "classical" or "bottom-up" explanation for the origin of UHECR, based on ideas already developed to explain the spectrum of CR at lower energies, i.e. through acceleration by electromagnetic fields; however, no object in our cosmological neighborhood is known which has a sufficiently strong magnetic field extending over a sufficiently large volume. The second possibility would require the existence of very massive and very long-lived X particles, with mass m X ≫ 10 20 eV and lifetime close to or larger than the age of the Universe. In these "topdown" scenarios the observed UHECR are the stable decay products of these ultra-massive particles (for reviews, see [4], [5]).
In this article we consider a generic top-down scenario. We provide an improved model of the decay of an X particle, which is independent of the origin and the nature of this particle. In leading order in perturbation theory X decays into a small number of very energetic particles. However, in reality this should be considered to be the starting point of a "parton shower", analogous to the formation of hadronic jets in decays of Z bosons studied at LEP. The existence of an energy scale m X much above the weak scale requires supersymmetry (SUSY) to stabilize the electroweak scale against radiative corrections. We therefore study this shower in the framework of the MSSM, which is the simplest potentially realistic supersymmetric extension of the Standard Model (SM). In contrast to previous works [6], [7], [8], we considered all gauge interactions as well as third generation Yukawa interactions, rather than only SUSY-QCD; note that at energies above 10 20 eV all gauge interactions are of comparable strength. The inclusion of electroweak gauge interactions in the shower gives rise to a significant flux of very energetic photons and leptons, which had not been identified in earlier studies. Moreover, we carefully modeled decays of all unstable particles. As a result, our treatment is the first that respects energy conservation.
In the next Section we describe our treatment of the parton shower in slightly more detail. In Sec. 3 we show some numerical results for the expected spectra of stable particles, and compare it with previous works. Finally, a brief summary and some conclusions are presented in Sec. 4.

Calculation of the fragmentation functions
Consider the two-body decay of an ultra-massive X particle of mass m X ≫ 10 20 eV into a qq pair, in the framework of the MSSM. This triggers a parton shower, which can be understood as follows. The q andq are created with initial virtuality Q X ∼ M X 2 . Note that Q 2 X > 0. The initial particles can thus reduce their virtuality, i.e. move closer to being on-shell (real particles), by radiating additional particles, which have initial virtualities < Q X .These secondaries then in turn initiate their own showers. The average virtuality and energy of particles in the shower decreases with time, while their number increases (keeping the total energy fixed, of course). As long as the virtuality Q is larger than the electroweak or SUSY mass scale M SUSY , all MSSM particles can be considered to be massless, i.e. they are all active in the shower. However, once the virtuality reaches the weak energy scale, heavy particles can no longer be produced in the shower; the ones that have already been produced will decay into SM particles plus the lightest superparticle (LSP), which we assume to be a stable neutralino. Moreover, unlike at high virtualities, at scales below M SUSY the electroweak interactions are much weaker than the strong ones; hence we switch to a pure QCD parton shower at this scale. At virtuality around 1 GeV, nonperturbative processes cut off the shower evolution, and all partons hadronize. Most of the resulting hadrons, as well as the heavy τ and µ leptons, will eventually decay. The end product of X decay is thus a very large number of stable particles: protons, electrons, photons, the three types of neutrinos and LSPs; we define the FF into a given particle to include the FF into the antiparticle as well, i.e. we do not distinguish between protons and antiprotons etc.
Note that at most one out of these many particles will be observed on Earth in a given experiment. This means that we cannot possibly measure any correlations between different particles in the shower; the only measurable quantities are the energy spectra of the final stable particles, dΓ X /dE P , where P labels the stable particle we are interested in. † This is a wellknown problem in QCD, where parton showers were first studied. The resulting spectrum can be written in the form [9] where I labels the MSSM particles into which X can decay, and we have introduced the scaled energy variable x = 2E/m X ; for a two-body decay, dΓ( z . All the nontrivial physics is now contained in the fragmentation functions (FFs) D P I (z, Q 2 ). Roughly speaking, they encode the probability for a particle P to originate from the shower initiated by another particle I, where the latter has been produced with initial virtuality Q. This implies the "boundary condition" which simply says that an on-shell particle cannot participate in the shower any more. For reasons that will become clear shortly, at this stage we have to include all MSSM particles J in the list of "fragmentation products". The evolution of the FF with increasing virtuality is described by the well-known DGLAP equations [9]: where α KI is the coupling between particles I and K, and the splitting functions P KI describes the probability for particle K to have been radiated from particle I. As noted earlier, for Q > M SUSY ∼ 1 TeV we allow all MSSM particles to participate in the shower. Since we ignore first and second generation Yukawa couplings, we treat the first and second generations symmetrically. I, J, K in eq.(3) thus run over 30 particles: 6 quarks q L , u R , d R , t L , t R , b R , 4 leptons l L , e R , τ L , τ R , 3 gauge bosons B, W, g, the two Higgs fields of the MSSM, and all their superpartners. Here we sum over all color and SU(2) indices (i.e., we assume unbroken SU(2) symmetry), and we ignore violation of the CP symmetry, so that we can treat the antiparticles exactly as the particles. All splitting functions can be derived from those listed for SUSY-QCD in [10]; explicit expressions will be given in a later paper. The starting point of this part of the calculation is eq.(2) at Q = M SUSY . This leads to 30 × 30 FFs D J I , which describe the shower evolution from Q X to M SUSY .
At scales Q < M SUSY all interactions except those from QCD can be ignored. Indeed, at scales Q < Q 0 ≃ 1 GeV QCD interactions become too strong to be treated perturbatively. leading to confinement of partons into hadrons. This nonperturbative physics cannot be computed yet from first principles; it is simply parameterized, by imposing boundary conditions on the D h i (z, Q 2 0 ), where h stands for a long-lived hadron and i for a light parton (quark or gluon). Here we used the results of [11], where the FFs of partons into protons, neutrons, pions and kaons are parameterized in the form Nx α (1 − x) β , using fits to LEP results. Note that these functions are only valid down to x = 0.1; for smaller x, mass effects become relevant at LEP energies. However, at our energy scale these effects are still completely irrelevant, even at x = 10 −7 . We chose a rather simple extrapolation at small x of the functions given in [11], of the form N ′ x α ′ . We computed the coefficients N ′ and α ′ by imposing energy conservation. We had to use the set of NLO FF given in [11], although our evolution equations are only written at the leading order; the LO set violates energy conservation badly, especially for the gluon FF. This choice renders difficult the comparison with previous results [8] [6], where the LO set was used. Apart from this difference, our results for the pure QCD case or for the SUSY QCD evolution seem to agree rather well. Starting from these modified input distributions, and evolving up to Q = M SUSY using the pure QCD version of eq.(3), leads to FFs D h i (x, M 2 SUSY ) which describe the QCD evolution (both perturbative and non-perturbative) at Q < M SUSY .
Finally, the two calculations have to be matched together. First we note that "switching on" SU(2) and SUSY breaking implies that we have to switch from weak interaction eigenstates to mass eigenstates. This is described by unitary transformations of the form D S I = J |c SJ | 2 D J I , with S |c SJ | 2 = J |c SJ | 2 = 1; here S stands for a physical particle. We used the ISASUSY code [12] to compute the SUSY mass spectrum and the |c SJ | corresponding to a given set of SUSY parameters. The decay of all massive particles S into light particles i is then described by adding S D S I ⊗ P iS to the FFs D i I of light particles i. We assumed that the x−dependence of the functions P iS originates entirely from phase space. In this fashion each massive particle S is distributed over massless particles i, with weight given by the appropriate S → i decay branching ratio. Note that this step often needs to be iterated, since heavy superparticles often do not decay directly into the LSP, so that the LSP is only produced in the second, third or even fourth step. All information required to model these cascade decays have again been taken from ISASUSY. The effects of the pure QCD shower evolution can now be included by one more convolution, Finally, the decay of long-lived but unstable particles µ, τ, n, π, K has to be treated; this is done in complete analogy with the decays of particles with mass near M SUSY .
Note that the linearity of the evolution equation (3) allowed us to factorize the problem in the fashion described above. We integrate these equations numerically using the Runge-Kutta method. The various FFs are modeled as cubic splines, with about 100 x values distributed logaritmically between x = 10 −7 and 0.5, and again logarithmically in 1 − x for x between 0.5 and 1. This allows us to compute the CR spectrum down to energies of the order of 10 18 eV even if M X = M GUT ≃ 10 16 GeV. We should stress that a good test of our algorithm was the possibility of checking energy conservation at all steps. We found that the "loss" of energy due to numerical approximations (the loss to particles with x < 10 −7 being negligible) never exceeds a few percent.

Results and discussion
We are now ready to present results for the FFs of any particle I of the unbroken MSSM produced at high virtuality through the decay of X into one of the 7 stable particles P . This requires 30 × 7 = 210 FFs for any set of SUSY parameters, which enable us to predict any CR spectrum near the source for any decay mode of X into MSSM particles. We saw earlier that for 2-body decays of X, the spectrum of P is directly given by the relevant FFs; eq.(1) shows that any N−body decay [8] can be treated with just one more convolution. We computed all 210 FFs for several sets of soft SUSY breaking parameters, but we found that most features depend very little on this choice. Here we only give results for the "low tan β gaugino region" (ratio of Higgs vevs tan β = 3.6, gluino mass mg ∼ 400 GeV, supersymmetric Higgs mass parameter µ ∼ 500 GeV, slepton masses around 140 GeV, squark masses around 300 GeV, CP-odd Higgs boson mass m A = 180 GeV and trilinear soft breaking parameter A t ∼ 1 TeV). As usual, we show results for x 3 × D P I (x, M X ), appropriate for comparing with the flattening of the observed spectrum (beyond the region of the knee characterized by a power law spectrum with power n ≃ −2.7). We take M X = 10 16 GeV, as appropriate for a GUT interpretation of the X particle. According to ref. [13], such a large value of M X is compatible with existing data if most UHECR originate at cosmological distances. Since the FFs evolve only logarithmically with Q, the final results for M X = 10 12 − 10 13 GeV [14,8] would not differ too much from the ones shown here, if they are expressed in terms of the scaling variable x.  Results for the fragmentation functions of first generation SU(2) doublet (s)quarks are shown in fig. 1. * For small x < ∼ 0.01, the FFs turn out to be very similar for these two cases. However, at large x > ∼ 0.1 some differences appear. In particular, an initial squark produces many more hard LSPs than a quark does; in the former case, LSPs carry ∼ 25% of the original squark energy, while they only carry ∼ 7% of the energy of an initial quark. Similarly, an initial squark produces a significantly larger flux of very energetic neutrinos, since neutrinos are frequently produced in sparticle decays; for tan β ≫ 1 the majority of these very energetic neutrinos would be ν τ , but for the given small value of tan β all three neutrino species are produced with equal abundance. In contrast, in case of an initial quark the three neutrino fluxes are of similar size only at very large x; these neutrinos come from the radiation of SU(2) × U(1) Y gauge bosons early in the shower, and their subsequent decay. At smaller values of x the ν τ component drops * Recall that we assume exact SU (2) invariance for Q > M SUSY and treat first and second generation (s)quarks symmetrically, so that q L is actually composed of u L , d L , c L , s L and their antiparticles, and similarly forq L . quickly relative to the ν µ and ν e components, since only the latter are produced in the decays of light mesons, in the ratio of roughly 2:1. Neither ultra-energetic LSPs nor neutrinos have as yet been observed. Among the particles with large cross sections on air, at large x we observe a significantly smaller FF into protons for an initial squark than for an initial quark. The reason is that a quark can directly fragment into a proton, but a squark will decay first, thereby distributing its energy over a larger number of softer particles. Note also that the electromagnetic component (sum of FFs into photons and electrons † ) is always bigger than the FF into protons. This is partly due to direct photon emission during the early stages of the shower; this effect, which was not included in earlier analyses, leads to a very hard component in the photon flux. This component is about two † In scenarios where an electromagnetic cascade can take place during the propagation in the extragalactic medium, the spectrum of primary electrons contributes to the spectrum of UHE photons (see [4]).
times larger for an initial quark than for a squark, due to the different splitting functions. Since experiments like Haverah Park [15] presently favor the hypothesis of hadronic primaries [16], this result could lead to problems for top-down models. However, interactions with the intergalactic medium will affect the electromagnetic component even more than the proton flux. Hence top-down models could still be compatible with a primary CR spectrum composed mostly of protons if the flux is dominated by X decays at cosmological distances.
Analogous results for SU(2) doublet (s)leptons of the first or second families are shown in fig. 2. In case of an initial lepton, the resulting flux for x > ∼ 0.1 remains dominated by leptons; due to our assumption of exact SU(2) symmetry at Q > M SUSY , electrons, ν e and ν µ contribute equally here. The flux of both LSPs and photons is much higher than that of protons in this case. This remains true for an initial slepton, where the LSP component is the dominant one for x > ∼ 0.2, followed by (charged or neutral) leptons and photons. In both cases the peak of the proton flux (after multiplication with x 3 ) is shifted down in energy by about a factor of 2 compared to the results of fig.1, and is reduced in size by about one order of magnitude for large x. The dominance of the electromagnetic component is therefore much more prominent for primary X decays into (s)leptons than for decays into (s)quarks.
In fig. 3 we compare the spectra of protons and photons as predicted by QCD, SUSY QCD and the full MSSM, for hadronic 2-body X decays; in the latter two cases we show results for X → q LqL as well as X →q LqL . We see that introducing superparticles into the parton shower reduces the FF into protons by about a factor of 2 at large x. Including superparticles not only opens up new channels, i.e. introduces new splitting functions in eq.(3), it also implies a larger SU(3) gauge coupling at high Q, which speeds up the shower evolution through standard splitting functions. Including in addition the full set of gauge and third generation Yukawa interactions has relatively little impact on the proton spectrum. Note, however, the tail of very energetic photons produced by the electroweak gauge interactions. As already seen in fig. 1, at large x the flux of both protons and photons is somewhat smaller for X → q decays than for X →q. Unfortunately it is difficult to directly compare fig. 3 with earlier analyses. As noted in the previous Section, we employ somewhat different input FFs at Q 0 = 1 GeV; moreover, most earlier papers did not include the contribution from sparticle decays.
We also compared the spectra of LSPs and leptons for the same set of assumptions as in fig. 3. The LSP flux is very insensitive to the inclusion of electroweak gauge and Yukawa interactions, but changes by a factor of ∼ 5 at x = 0.5 between X → qq and X →qq decays. The FFs into leptons gain a very hard component (which, however, falls approximately linearly as x → 1) once electroweak gauge interactions are included; moreover, if X decays produce a primary squark, subsequent sparticle decays give rise to a large flux of hard leptons at x ∼ 0.25.
So far we have focused on the large−x region of the FFs. All FFs increase rapidly as x decreases; this is not apparent in the figures, since the increase is over-compensated by the x 3 scaling factor used there. This increase towards small x can have quite dramatic effects. For example, we found that the total multiplicity of all stable particles with x ≥ x min can approximately be described by N tot (x min ) ∼ 2 · x −0.68 min , if X decays into quarks and/or squarks. This means that a single X decay will give rise to ∼ 10 5 stable particles with energy exceeding 5 · 10 −8 M X ! The result for pure QCD is only about a factor of 2 lower. The reason is that the  shower evolution between M SUSY and 1 GeV is not affected by the presence of superparticles. In this region QCD interactions are strongest; in the pure QCD case this part of the shower evolution therefore contributes more to the final multiplicity than the evolution at Q > 1 TeV. The large multiplicity also implies that a full Monte Carlo simulation of this shower [14] would require a very large numerical effort; this illustrates the advantage of using fragmentation functions.
We saw earlier that at large x, LSPs play an important role even if the primary X decay does not involve superparticles. We found that the integrated LSP multiplicity can be parameterized as N LSP (x min ) ∼ 0.65x −0.32 min for x min < ∼ 0.01, if X decays into two (s)quarks. This is somewhat smaller than the result of ref. [17], N LSP (x min ) ∼ 0.6x −0.4 min . It is also much smaller than the multiplicity of neutrinos, which is about half of the total multiplicity discussed above. The best hope for detecting ultra-energetic LSPs, which would be a "smoking gun" signature for any top-down scenario [17], might therefore lie in the LSPs with energies not far below M X , where the neutrino background is smallest. This conclusion is strengthened if most UHECR originate at cosmological distances, since then the interaction of the original protons with the intergalactic medium will produce [4] pions, and hence additional neutrinos; of course, these will also have energies well below M X . LSPs essentially do not interact with the intergalactic medium. Nevertheless medium effects will be crucial for turning our FFs into LSPs into predictions for LSP fluxes on Earth, since the flux of protons (possibly plus photons), which is affected by the intergalactic medium, always has to be normalized to the measured value. Moreover, bottom-up models originally produce almost no tau neutrinos. However, atmospheric neutrino data indicate large mixing between ν µ and ν τ , in which case the original ν µ flux will be distributed equally between ν µ s and ν τ s after at most a few kpc of propagation. Therefore a very energetic ν τ flux can unfortunately by itself not discriminate between these kinds of models.
We also saw in figs. 1 and 2 that a flux of neutrinos with energies above those of the most energetic protons is a generic feature of top-down models. This signal is especially pronounced if primary X decays produce squarks and/or (s)leptons; it will be strengthened by interactions of the protons with the medium, which move the protons to lower energies. In contrast, in bottom-up models one expects the neutrino flux to cut off at significantly lower energies than the original proton flux. However, even in this kind of model the neutrino flux might extend to higher energies than the observed proton flux once the protons have traveled for a few tens of Mpc. A very hard component in the neutrino flux can therefore only be considered to be a signature for top-down models if it can be shown that most UHECR are produced "locally", or if the neutrino spectrum is as shown in the upper fig. 2, with a sharp peak at the highest energy.

Summary and Conclusions
In this paper we presented a new and quite complete algorithm for treating the decay of a superheavy X particle in the framework of the MSSM. We improved previous analyses by including the full set of MSSM gauge and third generation Yukawa interactions, and by carefully modeling the decays of heavy sparticles and particles. We applied this algorithm to the "topdown" solution of the problem of UHECR. We gave the results in terms of fragmentation functions for any possible decay product of X into stable particles, and studied the consequences for a few primary decay modes of the X particle itself. An important result of this complete study, compared to a simplified SUSY QCD treatment, was the prediction of sizable fluxes of leptons and photons at energies above the peak of the proton spectrum. We also discussed possible ways to distinguish between bottom-up and top-down scenarios, and concluded that the cleanest signal would probably come from the detection of LSPs with energies above the peak of the proton spectrum. Of course, detailed analyses are required before we can conclude that these LSPs are indeed detectable on the background of neutrinos.
So far the only UHECR that have been detected are hadrons, or perhaps photons. If propagation effects can be neglected, our results can be used directly to fit the mass M X of the primary. If X decays primarily into (s)quarks, the result of such a fit would presumably resemble that of [8], if we assume that all events with energy above a few times 10 19 eV originate from X decay. This would favor a relatively "low" M X ∼ 10 12 GeV, well below the GUT scale. On the other hand, ref. [13] found a much higher value of M X , near the GUT scale, if X decays at cosmological distances. However, ref. [13] used a relatively crude model for X decay. ‡ It might be worthwhile to re-do this analysis, which requires the careful treatment of propagation effects. Finally, there is some evidence that there are two different sources for post-GZK events. For example, it has been claimed [18] that the number of "doublet" and "triplet" events (which originate from a small patch in the sky) is too high to be compatible with an essentially isotropic distribution of sources, which is the most plausible assumption for a top-down model § . Current statistics is poor, but there is some indication that this clustering of events occurs only at energy below 10 20 eV. Moreover, most of the experiments (Agasa, Yakutsk, and the Fly's eye collaboration) indicate that there might be a cut-off in the spectrum, at an energy ∼ 4 · 10 19 eV (see [21] for a review). This might also hint towards a two-component explanation for the events before and above the cut-off [22]. If so, it might be best to only use events with energy well above the GZK cutoff, say with E > 10 20 eV, in the fit of M X ; in that case the small number of events would presumably allow large values of M X even if most X decays are "local". In any case, when expressed in terms of the dimensionless scaling variable x = 2E/M X , the fragmentation functions only depend logarithmically on M X .
If M X ≫ 10 21 eV the events that have been observed so far would all have x ≪ 1. The fluxes at small x are much higher than in the large x region. Even if the observed UHE events come from the large x region, i.e. if M X ∼ 10 21 eV, care has to be taken not to over-produce particles at lower energies. In particular, stringent limits exist on the fluxes of photons in the TeV energy region, and on neutrinos in the multi-TeV region. The interpretation of these limits in the framework of a given top-down model again depends on where X decays occur, i.e. if propagation effects are important or not. We note here that in our case the evolution equation (3) can still be applied at x as small as 10 −7 . Coherence effects, which give rise to the "MLLA" description of the parton shower, become large [23] at yet smaller values of x ∼ Q h /M X , where Q h characterizes the hadron mass scale where strong interactions become non-perturbative. We also checked that the precise form of our small−x extrapolation of the nonperturbative FF's is not important, since the small−x behavior of the final FF's is essentially determined by the perturbative evolution. However, perturbative higher order corrections might be sizable at small x; this effect is currently being investigated.
To summarize, we presented a first treatment of superheavy particle decay that includes all relevant particle physics effects at the leading-log order of perturbation theory; in particular, for the first time we were able to account for the entire energy released in the decay. This provides a tool which can be used for detailed tests of the hypothesis that the most energetic cosmic rays originate from the decay of ultra-massive particles. If true, this would give us experimental access to energies well beyond the range that can ever be tested by experiments working at Earth-based colliders. ‡ For example, their input FFs violate energy conservation badly. § However, clustering of events can be explained in top-down models with "clumpy halo" [19,20]