text:
# X-ray iron line variability: constraints on the inner accretion disk

## 1. Introduction

AGN are observed to be copious X-ray emitters. These X-rays are thought to originate from the innermost regions of an accretion disk around a central supermassive black hole. Since the accretion disk itself is expected to be an optical/UV emitter, the most likely mechanism producing the X-rays is inverse Compton scattering of these soft photons in a hot and tenuous corona that sandwiches the accretion disk. Thus, in principle, the study of these X-rays should allow the immediate environment of the accreting black hole as well as the exotic physics, including strong-field general relativity, that operates in this environment to be probed.

In the past decade, X-ray astronomy has begun to fulfill that promise. Guided by observations with the Ginga, ASCA, RXTE and BeppoSAX satellites, there is a broad consensus that X-ray irradiation of the surface layers of the accretion disk in a class of AGN known as Seyfert 1 galaxies gives rise to fluorescent K$`\alpha `$ emission line of cold iron via the process of “X-ray reflection”. Accompanying the iron line is a reflected continuum possessing a characteristic shape that is determined by the competing effects of photoelectric absorption and Compton scattering.

In this contribution, I will review the basic physical processes leading to the observed X-ray reflection signatures, and the evidence leading us to believe that they do indeed originate from the innermost regions of the accretion disk. We will then discuss three recent observational campaigns on NGC 4258 (M106), NGC 5548 and MCG–6-30-15. The purpose of presenting results from these campaigns is to describe the successes, and mysteries, resulting from recent studies of X-ray reflection signatures.

## 2. X-ray reflection spectra and broad iron lines

The basic physics of X-ray ‘reflection’ can be understood by considering a hard X-ray (power-law) continuum illuminating a semi-infinite slab of cold gas. In this context, ‘cold’ is taken to mean that metal atoms are essentially neutral, but H and He are mostly ionized. When a hard X-ray photon enters the slab, it is subject to a number of possible interactions: Compton scattering by free or bound electrons, photoelectric absorption followed by fluorescent line emission, or photoelectric absorption followed by Auger de-excitation. A given incident photon is either destroyed by Auger de-excitation, scattered out of the slab, or reprocessed into a fluorescent line photon which escapes the slab.

Figure 1 shows the results of a Monte Carlo calculation which includes all of the above processes (Reynolds 1996; based on similar calculations by George & Fabian 1991). Due to the energy dependence of photoelectric absorption, incident soft X-rays are mostly absorbed, whereas hard photons are rarely absorbed and tend to Compton scatter back out of the slab. This gives the reflection spectrum a broad hump-like shape. In addition, there is an emission line spectrum resulting primarily from fluorescent K$`\alpha `$ lines of the most abundant metals. The iron K$`\alpha `$ line at 6.4 keV is the strongest of these lines.

For most geometries relevant to this discussion, the observer will see this reflection component superposed on the direct (power-law) primary continuum. Under such circumstances, the main observables of the reflection are a flattening of the spectrum above approximately 10 keV (as the reflection hump starts to emerge) and an iron line at 6.4 keV. For solar or cosmic abundances and a plane-parallel slab geometry, the expected equivalent width of the iron line is 150–200 eV (George & Fabian 1991; Reynolds, Fabian & Inoue 1995).

A major advance came with the launch of ASCA and its medium-resolution CCD spectrometers. A long (4.5 day) observation of the bright Seyfert galaxy MCG$``$6-30-15 in July-1994 allowed the iron line profile to be determined with some accuracy. The resulting line profile is shown in Fig. 2 (published by Tanaka et al. 1995). The line was found to be extremely broad (almost $`10^5\mathrm{km}\mathrm{s}^1`$ FWZI) and asymmetric in so far that it possesses an extensive red-wing. Such a broad and asymmetric line is expected if the X-ray reflection is occurring in the inner regions of an accretion disk — strong line-of-sight Doppler shifts, transverse Doppler shifts and gravitational redshifts combine to produce an extensive low-energy wing and a sharp truncation of the line at high-energies (Fabian et al. 1989). Tanaka et al. (1995) showed that the MCG$``$6-30-15 data is in good agreement with the disk model provided the inclination of the disk is $`\theta 27^{}`$ and line fluorescence occurs down to $`6GM/c^2`$, the innermost stable orbit around a Schwarzschild black hole. This result is the best evidence to date for a radiatively-efficient accretion disk around a black hole in any object.

Subsequent studies of large samples of objects by Nandra et al. (1997) confirmed the presence of these features in many other Seyfert 1 galaxies, and show that there is a tendency for the iron lines to indicate face-on accretion disks (as expected for Seyfert 1 nuclei from the unified Seyfert scheme).

## 3. Alternative models for broad iron lines

The claim that iron line studies are probing the region within a few gravitational radii of the black hole is a bold one, and should be tested against other models at every opportunity. Given the quality of data, the July-1994 MCG$``$6-30-15 line profile has become a test bed for such comparisons.

Fabian et al. (1995) examined many alternative models including lines from mildly relativistic outflows, the effect of absorption edges on the observed spectrum, and broadening of the line via comptonization. Fabian et al. found that none of these models were viable alternatives for the MCG$``$6-30-15 line profile. The idea of producing the broad line via Comptonization has been revived recently by Misra & Kembhavi (1997) and Misra & Sutaria (1999). They argue that the spectrum initially consists of a narrow iron line superposed on a power-law continuum and that Comptonization in a surrounding cloud with optical depth $`\tau 4`$ produces the broad line. The Comptonizing cloud must be both cold ($`kT<0.5`$ keV in order to predominately downscatter rather than upscatter the line photons), and fully-ionized (since no strong iron absorption edges are seen in the continuum spectrum). The cloud is kept fully ionized and yet cool by postulating that the continuum source has a very luminous optical/UV component.

There are strong arguments against such a model. Since the power-law continuum emission also passes through any such Comptonizing cloud, one would observe a break in the continuum spectrum at $`E_{\mathrm{br}}m_\mathrm{e}c^2/\tau ^230`$ keV. Such a break is not observed in the BeppoSAX (Guainazzi et al. 1999) or RXTE data (Lee et al. 1999) for MCG–6-30-15 (see Misra 1999). Also, both continuum variability (which is seen on timescales as short as 100 s) and ionization arguments limit the size of the Comptonizing cloud in MCG–6-30-15 to $`R<10^{12}`$ cm. The essence of this ionization argument is that the ionization parameter at the outer edge of the cloud (which, for a fixed cloud optical depth, scales with cloud size as $`\xi 1/R`$) must be sufficiently high that all abundant metals, including iron, are fully ionized (Fabian et al. 1995; Reynolds & Wilms 2000). In the case of MCG–6-30-15, these constraints on the cloud size turn out to so tight that the postulated optical/UV component required to Compton cool the cloud would violate the black body limit (Reynolds & Wilms 2000). Comptonization moreover provides a poor fit (Ruszkowski & Fabian 2000). Hence, we consider the Comptonization model for the broad iron line not to be viable.

In another alternative model, Skibo (1997) has proposed that energetic protons transform iron in the surface of the disk into chromium and lower $`Z`$ metals via spallation which then enhances their fluorescent emission. With limited spectral resolution, such a line blend might appear as a broad skewed iron line. This model suffers both theoretical and observational difficulties. On the theoretical side, high-energy protons have to be produced and slam into the inner accretion disk with a very high efficiency (Skibo assumes $`\eta =0.1`$ for this process alone). On the observational side, it should be noted that the broad line in MCG–6-30-15 (Tanaka et al 1995) is well resolved by the ASCA SIS (the instrumental resolution is about 150 eV at these energies) and it would be obvious if it were due to several separate and well-spaced lines spread over 2 keV.

## 4. Case study I — the low luminosity AGN NGC 4258

Low luminosity AGN (LLAGN) are almost certainly systems which are accreting at very low Eddington fractions. This places them in a regime of accretion which has been the subject of active theoretical work for many years. It was realized by several authors that when the accretion rate is low (relative to the Eddington rate), an accretion disk may switch into a hot, radiatively-inefficient mode (Ichimaru 1977; Rees 1982; Narayan & Yi 1994; Narayan & Yi 1995). In essence, the plasma becomes so tenuous that the timescale for energy transfer from the protons to the electrons (via Coulomb interactions) becomes longer than the inflow timescale. The energy remains as thermal energy in the protons (which are very poor radiators) and gets advected through the event horizon of the black hole. These are the so-called Advection Dominated Accretion Flows (ADAFs111In recent years there has been a large amount of work on advective accretion flows, with a corresponding expansion in terminology and acronims (e.g. ADIOS, BDAF, and CDAF). See Prof. Narayan’s contribution in these proceedings for a fully discussion of the various types of advective accretion flow. Here, we use the term ADAF to encompass all such accretion flows which are radiatively inefficient.). ADAFs are to be contrasted with ‘standard’ radiatively-efficient accretion disks in which the disk remains cool and geometrically thin all of the way down to the black hole (Shakura & Sunyaev 1973; Novikov & Thorne 1974).

LLAGN provide an ideal laboratory in which to study ADAFs. If the central regions of LLAGN accretion disks are operating in an ADAF mode, we would not expect to observe broad iron lines since the ADAF is far too hot to contain combined iron ions. On the other hand, if we see broad iron lines from a LLAGN, this would be strong evidence for radiatively efficient disks even at these low accretion rates. However, the fact that LLAGN are faint has hampered such detailed X-ray studies of these sources.

In an attempt to study one LLAGN in some detail, we performed a 200 ksec joint ASCA/RXTE observation of the nearby LLAGN NGC 4258. A full discussion of this observation is presented in Reynolds, Nowak & Maloney (2000). Here, we summarize the iron line properties of this object.

As can be seen in Fig. 3, we clearly detect a iron line with an equivalent width of $`W_{K\alpha }=107_{37}^{+42}`$ eV and a velocity width of $`<22000\mathrm{km}\mathrm{s}^1`$ (FWHM) which is much narrower than those seen in higher luminosity Seyfert galaxies. If it is assumed that the observed iron line originates from the accretion disk, then we can use relativistic “disk-line” models to place constraints on the location of the line emission. Using such models it is found that the narrowness of the line implies that the bulk of the line emission originates from a distance $`r>100GM/c^2`$ from the black hole. Such a large line emitting radius is consistent with the ADAF model for NGC 4258 (Lasota et al. 1996) in which the thin fluorescing accretion disk goes through an transition into an ADAF at about $`r100GM/c^2`$. However, the data are also consistent with a radiatively efficient fluorescing disk extending down to the radius of marginal stability and a large ($`r100GM/c^2`$) X-ray corona. In this latter case, there would be faint broad wings to the line profile which could be found with a high signal-to-noise XMM–Newton spectrum. One way or another, the iron line properties of this LLAGN are different to those of many higher-luminosity Seyfert galaxies which display a broad iron line indicating radiative disks and small ($`10GM/c^2`$) X-ray emitting regions.

We cannot rule out the possibility that most (or all) of the observed iron line originates from matter that is not associated with the accretion disk. However, simple arguments lead us to disfavor such a scenario. Consider iron line emission in a geometrically-thick torus surrounding the accretion disk of NGC 4258. An upper limit to the column density of this structure along our line of sight to the AGN is given by the observed column density of $`N_\mathrm{H}10^{23}\mathrm{cm}^2`$. If we suppose that this torus is in the same plane as the accretion disk (so that we are also viewing it edge-on), it is plausible to assume that we are looking through the optically-thickest part of the torus. We can then show that the maximum equivalent width produced by fluorescence of this material is $`W_{\mathrm{Fe},\mathrm{max}}65`$ eV (Reynolds, Nowak & Maloney 2000). However, there are several reasons why this upper limit would almost certainly not be achieved. Firstly, modeling of the accretion disk warp strongly suggests that our line of sight intercepts the disk and that the bulk of the column density which obscures the AGN originates in the disk (Herrnstein, priv. communication). Therefore, we might expect significantly smaller column densities along lines of sight that have smaller inclinations angles relative to the accretion disk. Secondly, the accretion disk may well occult half of this fluorescing cloud, thereby reducing this prediction further. Thus, the true iron line from surrounding non-disk material may well be reduced from our naive prediction by a factor of several.

Despite these arguments, let us now suppose that the observed iron line is not associated with the accretion disk. Our data then permit the presence of a “Seyfert-like” broad iron line which is buried in the noisy continuum data. If we suppose NGC 4258 to possess a normal “Seyfert-like” broad iron line (with an intrinsic equivalent width in the range 200–400 eV, modified by the light bending and limb-darkening effects that are important in edge-on systems such as NGC 4258), we can set some interesting constraints (see Fig. 4). In particular, we can rule out models in which the limb-darkening is very weak (as would be the case if, for example, the $`\tau _\mathrm{e}=1`$ surface of the accretion disk is very filamentary rather than planar). Of course, higher signal-to-noise data from XMM–Newton are required if one is to constrain the nature of the disk in this object any further.

## 5. Case study II — an EUV/X-ray campaign on NGC 5548

In the rest of this contribution, we turn to the properties of higher luminosity systems. Firstly, we shall discuss the classical Seyfert galaxy NGC 5548.

We222The team of investigators consisted of Omer Blaes (UCSB), James Chiang (Colorado), Greg Madjeski (NASA-GSFC), Pawel Magdziarz (deceased), Hermann Marshall (MIT-CSR), Mike Nowak (Colorado) and myself. observed NGC 5548 with EUVE, ASCA, and RXTE in June/July 1998. The main part of the campaign consisted of four observations (with all of these satellites simultaneously) separated by approximately one week. One of these observations (18-June-1998) was substantially longer than the others ($`110`$ ksec of good data, as compared with an average of $`30`$ ksec for the other observations). This long observation provided us with intensive monitoring over a period of 3 days from the EUV to hard X-rays. The various aspects of this campaign are discussed in detail in Chiang et al. (2000).

Although not strictly related to X-ray reflection studies, it is worth mentioning the continuum variability properties found during this campaign. Figure 5 shows the EUVE, ASCA and RXTE light curves for our long observation. The most distinct feature is a step at $`2\times 10^4`$ s. By computing the z-transformed discrete correlation function (ZDCF; Alexander 1997) between these various data, we find clear evidence that the EUV emission leads the RXTE-band emission by $`30`$ ksec, and the ASCA-band emission by $`15`$ ksec. If confirmed, this results has at least two important implications. Firstly, the EUV emission cannot simply be reprocessed hard X-ray emission since these scenario would predict the EUV to lag the hard X-rays. Secondly, if we make the assumption that the observed lags are due to the Compton upscattering of seed photons, we can constrain the size of the Comptonizing region to be about $`10GM/c^2`$ for a $`10^8M_{}`$ black hole.

Returning to the issue of X-ray reflection, this campaign produced a major surprise. Before discussing the data, it is worth recapping our theoretical expectations from simple X-ray reflection models. For a fixed disk/corona geometry and ionization state, we expect the iron line equivalent width to remain constant provided one is probing timescales longer than the reverberation timescale of the system (which is always the case for spectral studies with ASCA and RXTE). If the iron line equivalent width does change due to, for example, a change in the geometry of the system, then the relative strength of the reflection continuum would be expected to vary in step. In other words, we expect a strict proportionality between the iron line equivalent width and relative strength of the reflection continuum.

Figure 6 shows the various spectral parameters, measured from the RXTE data, for NGC 5548 as a function of the 2–10 keV continuum flux of the source. We found that the equivalent width of the iron line declined as the continuum source flux increased — in fact, the iron line flux was consistent with being constant. Simultaneous spectroscopy with ASCA showed the iron line to be broadened. Hence, the line does appear to be originating from the central regions ($`r<50GM/c^2`$) of the accretion disk and, for any reasonable size black hole mass, the constancy of the line flux cannot be due to light travel time effects. Moreover, the relative reflection fraction ($`R`$) fails to track the changes in the iron line equivalent width. Indeed, these quantities appear to be weakly anti-correlated (see Fig. 7), contrary to the simple reflection models.

## 6. Case study III — the X-ray campaign on MGC–6-30-15

A very long, simultaneous ASCA/RXTE observation of the Seyfert 1 galaxy MCG–6-30-15 by an (almost) independent team of researchers333The team of investigators consisted of Niel Brandt (PSU), Andy Fabian (Cambridge), Kazushi Iwasawa (Cambridge), Julia Lee (MIT) and myself. revealed very similar spectral results to those obtained for NGC 5548. Figure 8 shows the clear anti-correlation between iron line equivalent width and the relative strength of the reflection continuum found by Lee et al. (2000). To produce this figure, almost nine days of RXTE data have been binned into four flux-sorted spectra. Unlike with NGC 5548, these data strongly require a change in the relative strength of the reflection continuum as well as the iron line equivalent width. The details of this observation are discussed by Lee et al. (1999, 2000).

This particular RXTE campaign yielded data of an almost unprecedented nature — the observation lasted over 100 consecutive orbits of RXTE, giving us dense coverage ($`>50\%`$) for almost 9 days. Reynolds (2000) took advantage of this unique dataset and performed a detailed search for delays between a continuum band (2–4 keV) and the iron line band (5–7 keV) in an attempt to search for iron line reverberation effects. In essence, the method of Press, Rybicki & Hewitt (1992) was used to interpolate across the gaps in the continuum light curve (caused primarily by Earth occultation) and estimate the errors in the interpolation. Figure 9a shows a small portion of the light curve, together with the Press et al. (1992) reconstruction. The continuous (interpolated) continuum light curve can then be folded through a trial transfer function and compared with the line-band light curve in a $`\chi ^2`$ sense (e.g., Fig. 9b shows the $`\chi ^2`$ surface that results from such a procedure). By varying the parameters of the trial transfer function, we can place constraints on any time delays between bands. No time delays were found on timescales down to 500–1000 s. Instead, it was found that the line flux remained constant, thereby extending the result of Lee et al. (2000) down to these timescales.

## 7. Evidence for flux-correlated changes in the disk’s ionization

To summarize the mysteries resulting from the NGC 5548 and MCG–6-30-15 campaigns, it is found that:

1. The iron line flux, rather than the equivalent width, appears to remain approximately constant as the continuum X-ray source undergoes rapid variability. In NGC 5548, this result is derived by direct X-ray spectral fitting and so applies on timescales of $`50`$ ksec (which is approximately timescale that can be probed via such methods). In MCG–6-30-15, the spectral fitting of Lee et al. (2000) combined with the analysis of Reynolds (2000) leads us to believe that such behavior occurs on all timescale from 500 ksec down to 500 s. This cannot be due to light travel delays.
2. The iron line equivalent width and relative strength of the reflection continuum are anti-correlated, in contradiction to simple X-ray reflection models.

Changing the ionization state of the surface of the accretion disk is one of the few ways of breaking the expected proportionality between the iron line equivalent and the (inferred) relative strength of the reflection continuum.

Reynolds (2000) proposed a toy model to explain these RXTE results (Fig. 10). In this model, there is a radial ionization gradient on the disk surface due to the radial dependence of both the ionizing flux and the density of the surface layers. The surface of the inner disk, which is strongly irradiated by the primary X-ray flux, is highly ionized and acts as a pure Compton mirror. It thereby fails to imprint distinct (atomic) features in the observable spectrum. By contrast, the outer disk surface (starting at $`2050GM/c^2`$) possesses a relatively low ionization parameter and produces ‘neutral-like’ X-ray reflection signatures such as those calculated by George & Fabian (1991). The interesting and uncertain processes occur in the partially ionized intermediate zone of the accretion disk. In this region, there are two possible physical mechanisms by which the iron line can be strongly suppressed relative to the reflected continuum. These mechanisms are described in the next two paragraphs. Reynolds (2000) noted that, if the boundaries between these zones are defined by some given set of ionization parameters, the zones will move as the primary continuum flux varies. If the primary continuum flux increases, the fluorescing outer zone will move outwards and so subtends a smaller solid angle as seen by the primary X-ray source. The observed iron line equivalent width will then be seen to decrease. However, the measured relative strength of the reflection continuum may remain roughly constant if the intermediate zone still produces bound-free absorption features in the reflection continuum.

As mentioned above, there are at least two mechanisms which may operate in the intermediately ionized zone of the accretion disk to suppress iron line emission relative to the reflection continuum. Nayakshin, Kazanas & Kallman (2000) showed that the well-known thermal instability of photoionized plasma could lead to the formation of a highly ionized skin on the surface of the accretion disk. Compton scattering in this skin would smear out sharp spectral features such as the iron line, possibly to the level where line photons from this region would no longer be identified with the iron line. However, very broad spectral features such as the iron edge and broad hump that observationally define the reflection continuum, would be much less affected by Compton scattering in the skin. Thus, from the observational point of view, the presence of such a skin might reduce the iron line strength while leaving the strength of the reflection continuum unaffected.

The second line destruction mechanism is that of resonant trapping followed by Auger destruction (Ross & Fabian 1993). A fluorescent iron line photon is typically created at an optical depth of $`\tau _\mathrm{e}1`$. However, when the bulk of the ions are at least as ionized as Fe xvii, the presence of a vacancy in the L-shell of the iron ions leads to resonant scattering of the line photons with a rather small path length. If, furthermore, the bulk of the ions possess at least two L-shell electrons (i.e. Fe xxiii or less), then the Auger effect acts as a photon destruction mechanism. The combination of resonant ‘trapping’ and Auger ‘destruction’ might reduce the observed iron line to a very small strength.

Clearly, more work is needed on the physics and modeling of ionization accretion disks in order to assess whether such a toy model can explain the mysterious X-ray variability that we see in current RXTE data.

## 8. Down to the reverberation timescale

The rapid X-ray variability of many Seyfert galaxies leads us to believe that the primary X-rays are emitted during dramatic flare-like events in the accretion disk corona. When a new flare becomes active, the hard X-rays from the flare will propagate down to the cold disk and excite iron fluorescence. Due to the finite speed of light, the illumination from the flare sweeps across the disk, and the reflected X-rays act as an ‘echo’ of this flare. Such flaring will cause temporal changes in the iron line profile and strength due to the changing illumination pattern of the disk and, more interestingly, time delays between the observed flare and the its fluorescent echo. This latter effect is directly analogous to the optical broad line reverberation that is a major theme of this meeting. The observational study of iron line reverberation is beyond the capabilities of current instruments (even, probably, XMM–Newton) but should be within reach of future high-throughput X-ray observatories such as Constellation-X and XEUS.

In recent work, it has been shown that interesting diagnostics of the black hole spin and X-ray source geometry are encoded within the iron line reverberation pattern (Reynolds et al. 1999; Young & Reynolds 2000; Ruzskowski 2000). However, there are many complications and issues that will affect our ability to extract these reverberation signatures from real data. Firstly, there are likely to be several X-ray flares active at any one time which will give overlapping reverberation signatures. Moreover, due to the fact that the X-ray emitting corona is probably extended on the same spatial scales as the inner disk, each X-ray emitting flare will have its own distinct transfer function. This breaks the linearity of the transfer problem. Young & Reynolds (2000) showed that if there are only a small number of powerful flares, the individual transfer functions can be identified and separated in Constellation-X data. However, we must examine more powerful statistical methods for extracting useful information from this system.

Possible reverberation signatures will be affected by short timescale ionization changes. The models of Reynolds et al. (1999) and Young & Reynolds (2000) assume that the disk has a fixed, fairly cold, ionization structure. However, as discussed extensively above, RXTE argues for ionization changes in the disk surface that are correlated with the continuum flux. On a local scale, this effect may be driven by ionization fronts that sweep across the disk when an X-ray flare erupts. In future work, these effects must be assessed so that instruments and observing strategies can be planned to study reverberation effects. It must be noted that, even though RXTE data argue for almost constant iron line flux due to the ionization changes, we would still expect to see variability of the line strength and/or profile on either the recombination timescale in the disk surface or the light-crossing timescale from the flare to the disk, whichever is longest.

Given all of these complications and difficulties, why bother to pursue the idea of iron line reverberation? The answer to this question is straightforward — until we have the technology to directly image the immediate environment of a black hole (and this is at least 20 years away), X-ray iron line reverberation is simply our best hope of mapping mapping out the detailed astrophysical environment of supermassive black holes, including the effects of the black hole’s spin.

## 9. Conclusions

In this contribution, we present three case studies in order to describe the current state of X-ray reflection studies in AGN. For the low-luminosity AGN NGC 4258, we find that the iron line is much narrower than is typically found in higher luminosity AGN. We argue that this is evidence for either a truncated cold accretion disk (possibly due to a transition to an ADAF at $`r100GM/c^2`$) or a large ($`r100GM/c^2`$) X-ray emitting corona surrounding the accretion disk. We also present results for the higher luminosity Seyfert nuclei in NGC 5548 and MCG–6-30-15. In both of these sources, RXTE shows that the iron line equivalent width decreases with increasing luminosity. Furthermore, the iron line equivalent width is found to be anticorrelated with the relative strength of the reflection continuum, contrary to all simple reflection models. It is proposed that continuum-flux correlated changes in the ionization of the accretion disk surface can explain this spectral variability. Finally, we address the issue of X-ray iron line reverberation in the light of these complicating factors.

## 10. Acknowledgements

The author appreciates support from Hubble Fellowship grant HF-01113.01-98A. This grant was awarded by the Space Telescope Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA under contract NAS 5-26555. He also appreciates support from NASA under LTSA grant NAG5-6337, and the National Science Foundation under grants AST-9529170 and AST-9876887.

## References

Alexander T., 1997, in Astronomical Time Series, ed. D. Maoz et al., (Dordrecht: Kluwer), 163

Chiang J. et al., 2000, ApJ, 528, 292

Fabian A. C., Rees M. J., Stellar L., White N. E., 1989, MNRAS, 238, 729

Fabian A. C. et al. 1995, MNRAS, 277, L11

George I. M., Fabian A. C., 1991, MNRAS, 249, 352

Guainazzi M. et al., 1999, A&A, 341, L27

Ichimaru S., 1977, ApJ, 214, 840

Lasota J.-P., Abramowicz M. A., Chen X., Krolik J., Narayan R. 1996, ApJ, 462, 142

Lee J., Fabian A. C., Reynolds C. S., Brandt W. N., Iwasawa K., 2000, MNRAS, in press

Lee J., Fabian A. C., Brandt W. N., Reynolds C. S., Iwasawa K., 1999, MNRAS, 310, 973

Misra R., 1999, IUCAA preprint 32/99

Misra R., Kembhavi A. K., 1998, ApJ, 499, 205

Misra R., Sutaria F. K., 1999, ApJ, 517, 661

Nandra K., George I. M., Mushotzky R. F., Turner T. J., Yaqoob T., 1997, ApJ, 477, 602

Nayakshin S., Kazanas D., Kallman T. R., 2000, ApJ, submitted (astro-ph/9909359)

Narayan R., Yi I., 1995, ApJ, 452, 710

Narayan R., Yi I., 1994, ApJ, 428, L13

Novikov I. D., Thorne K. S., 1973, in Black Holes, eds C. DeWitte & B. S. DeWitte (Gordon and Breach Science Publicaters, New York), P.344

Press W. H., Rybicki G. B., Hewitt J. N., 1992, ApJ, 385, 404

Rees M. J., 1982, in Riegler G., Blandford R. D., eds, The Galactic Center, Am. Inst. Phys., New York, p.166.

Reynolds C. S., 1996, PhD thesis, University of Cambridge

Reynolds C. S., 2000, ApJ, 533, 811

Reynolds C. S., Wilms J., 2000, ApJ, 533, 821

Reynolds C. S., Fabian A. C., Inoue H., 1995, MNRAS, 276, 1311

Reynolds C. S., Nowak M. A., Maloney P. R., 2000, ApJ, in press (astro-ph/0004068)

Reynolds C. S., Young A. J., Begelman M. C., Fabian A. C., 1999, ApJ, 514, 164

Ross R. R., Fabian A. C., 1993, MNRAS, 261, 74

Ruszkowski M., 2000, MNRAS, 315, 1

Ruszkowski M., Fabian A. C., 2000, MNRAS, submitted

Shakura N. L., Sunyaev R. A., 1973, A&A, 24, 337

Skibo, J.G., 1997, ApJ, 478, 522

Tanaka Y. et al., 1995, Nat, 375, 659

Young A. J., Reynolds C. S., 2000, ApJ, 529, 101


format: text

text:
# REFERENCES

Experimental test of local realism using non-maximally entangled states

Marco Genovese, G. Brida and C.Novero Istituto Elettrotecnico Nazionale Galileo Ferraris, Str. delle Cacce 91, I-10135 Torino

E. Predazzi Dip. Fisica Teorica Univ. Torino and INFN, via P. Giuria 1, I- 10125 Torino

Abstract

In this paper we describe a test of Bell inequalities using a non- maximally entangled state, which represents an important step in the direction of eliminating the detection loophole. The experiment is based on the creation of a polarisation entangled state via the superposition, by use of an appropriate optics, of the spontaneous fluorescence emitted by two non-linear crystals driven by the same pumping laser.

Entanglement is the main resource of Quantum Information Processing and thus deserves to be widely studied. One of the most characteristic properties of entanglement are the correlations between space-like separated entangled systems studied by Einstein-Podolsky-Rosen, for which Bell showed that no Local Hidden Variable Theory can reproduce all the results of Quantum Mechanics and, in particular, that some inequality, valid for every LHV theory, is violated by QM.

Many interesting experiments have been devoted to a test of Bell inequalities, the most interesting of them using photon pairs , leading to a substantial agreement with quantum mechanics and disfavouring LHV theories. But, up to now, no experiment has yet been able to exclude definitively such theories , because of the low total detection efficiency, which requires the hypothesys that the observed sample of particle pairs is a faithful subsample of the initial set of pairs. In the 90’s a big progress in the direction of eliminating this loophole has been obtained by using parametric down conversion (PDC) processes, which overcomes some former limitations, as the poor angular correlation of atomic cascade photons . The first experiments were performed with type I PDC, which gives phase and momentum entanglement and can be used for a test of Bell inequalities using two spatially separated interferometers . The use of beam splitters, however, strongly reduces the total quantum efficiency. In alternative, a polarisation entangled state can be generated , but in most of the used configurations, half of the initial photon flux is lost.

Recently, an experiment where a polarisation entangled state is directly generated, has been realised using Type II PDC . This scheme has permitted, at the price of delicate compensations for having identical arrival time of the ordinary and extraordinary photon, a much higher total efficiency than the previous ones, which is, however, still far from the value of $`0.81`$ required for eliminating the detection loophole for a maximally entangled state. A large interest remains therefore for new experiments increasing total quantum efficiency in order to reduce and finally overcome the efficiency loophole.

For this purpose, we have considered  the possibility of generating a polarisation entangled state via the superposition of the spontaneous fluorescence emitted by two non-linear crystals (rotated for having orthogonal polarisation) driven by the same pumping laser . The crystals are put in cascade along the propagation direction of the pumping laser and the superposition is obtained by using an appropriate optics. If the path between the two crystal is smaller than the coherence length of the laser, the two photons path are indistinguishable and a polarisation entangled state is created. The possibility of easily obtaining a non maximally entangled state is very important, for it has been recognised that in this case the lower limit on the total detection efficiency for eliminating the detection loophole is as small as 0.67 .

The experimental set-up is composed of two crystals of $`LiIO_3`$ placed along the pump laser propagation, 250 mm apart, a distance smaller than the coherence length of the pumping laser. This guarantees indistinguishability in the creation of a couple of photons in the first or in the second crystal. A couple of planoconvex lenses of 120 mm focal length centred in between focalises the spontaneous emission from the first crystal into the second one. A hole of 4 mm diameter is drilled into the centre of the lenses to allow transmission of the pump. A small quartz plate compensates the displacement of the pump due to birefringence in the first crystal. Finally, a half-wavelength plate immediately after the condenser rotates the polarisation of the Argon beam. We have used as photo-detectors two avalanche photodiodes with active quenching (EG&G SPCM-AQ) coupled to an optical fiber.

A very interesting degree of freedom of this configuration is given by the fact that by tuning the pump intensity between the two crystals, one can easily tune the ratio $`f`$ between the two component of the entangled state ($`f=1`$ for a maximally entangled state). Furthermore, every Bell state can be generated (for example for applications to quantum information processing) changing the phase between the two components (moving the second crystal) and/or changing the polarisation on one of the branches after the second crystal. The main difficulty of this configuration is in the alignment, which is of fundamental importance for having a high visibility. This problem has been solved using  an optical amplifier scheme, where a solid state laser is injected into the crystals together with the pumping laser. We think that the proposed scheme is well suited for leading to a further step toward a conclusive experimental test of non-locality in quantum mechanics. The main advantage of the proposed configuration with respect to most of the previous experimental set-ups is that all the entangled pairs are selected (and not only $`<50\%`$ as with beams splitters), furthermore it does not require delicate compensations for the optical paths ofthe ordinary and extraordinary rays after the crystal.

At the moment, the results which we are going to present are still far from a definitive solution of the detection loophole; nevertheless, being the first test of Bell inequalities using a non-maximally entangled state, they represent an important step in this direction. Furthermore, this configuration permits to use any pair of correlated frequencies and not only the degenerate ones. We have thus realised this test using for a first time two different wave lengths (at $`633`$ and $`789`$ nm).

It must be acknowledged that a set-up for generating polarisation entangled pairs of photons, which presents analogies with our, has been realised recently in Ref.  The main difference between the two experiments is that in  the two crystals are very thin and in contact with orthogonal optical axes: this permits a ”partial”  superposition of the two emissions with opposite polarisation. This overlapping is mainly due to the finite dimension of the pump laser beam, which reflects into a finite width of each wavelength emission. A much better superposition can be obtained with the present scheme. More recently , they have also performed a test about local realism using a particular kind of Hardy equalities . Their result is in agreement with quantum mechanics modulo the detection loophole, no discussion concerning the elimination of loopholes for this equality is presented.

As a first check of our apparatus, we have measured the interference fringes, varying the setting of one of the polarisers, leaving the other fixed. We have found a high visibility, $`V=0.973\pm 0.038`$.

Our results are summarised by the value obtained for the Clauser-Horne sum

$$CH=N(\theta _1,\theta _2)N(\theta _1,\theta _2^{})+N(\theta _1^{},\theta _2)+N(\theta _1^{},\theta _2^{})N(\theta _1^{},\mathrm{})N(\mathrm{},\theta _2)$$
(1)
which is strictly negative for local realistic theory. In (1), $`N(\theta _1,\theta _2)`$ is the number of coincidences between channels 1 and 2 when the two polarisers are rotated to an angle $`\theta _1`$ and $`\theta _2`$ respectively ($`\mathrm{}`$ denotes the absence of selection of polarisation for that channel). On the other hand, quantum mechanics predictions for $`CH`$ can be larger than zero. In our case we have generated a state with $`f0.4`$: in this case the largest violation of the inequality is reached for $`\theta _1=72^o.24`$ , $`\theta _2=45^o`$, $`\theta _1^{}=17^o.76`$ and $`\theta _2^{}=0^o`$. Our experimental result is $`CH=512\pm 135`$ coincidences per second, which is almost four standard deviations different from zero and compatible with the theoretical value predicted by quantum mechanics.


text:
# Power correction analyses in e⁺e⁻ annihilation

## 1 INTRODUCTION

QCD studies always have to consider how the observed final states consisting of hadrons, e.g. from the process $`\mathrm{e}^+\mathrm{e}^{}`$ hadrons, map on QCD predictions in terms of partons, the quarks and gluons. The transition from the partons of the hard process to the observed hadrons, commonly referred to as hadronisation, is a fundamental property of QCD.

However, there is no fundamental QCD based theory of the transition from partons to hadrons. The most successful description of the hadronisation process in $`\mathrm{e}^+\mathrm{e}^{}`$ annhilation has been achieved with Monte Carlo simulations based on a parton shower followed by formation of hadrons. Recently, analytic models of the hadronisation have become available which predict how hadronisation effects of QCD observables, e.g. event shapes, scale as a power of the energy scale of the hard process. These models are often referred to as power corrections.

In this report we will briefly review some of the analytic models and present comparisons with $`\mathrm{e}^+\mathrm{e}^{}`$ annihilation data. We will also show a more general test of the consistency of power corrections with the gauge structure of QCD.

## 2 PREDICTIONS

We summarise three related approaches to predict the power correction of event shape observables. In all cases the form of the dependence on the scale of the hard process is extracted from an analysis of infrared renormalon singularities. The differences between the approaches come from the regularisation of the infrared singularities in order to obtain meaningful predictions.

### 2.1  DMW Model

The model of Dokshitzer, Marchesini and Webber (DMW)  makes the assumption that evolution of $`\alpha _s`$ to energies below the Landau pole is possible but the form of $`\alpha _s(\mu )`$ is a priori unknown. A non-perturbative parameter $`\alpha _0`$ is introduced as the 0th moment over $`\alpha _s(\mu )`$:

$$\alpha _0=\frac{1}{\mu _I}_0^{\mu _I}\alpha _s(k)𝑑k.$$
(1)
The quantity $`\mu _I`$ is the infrared matching scale where the non-perturbative and perturbative evolution of $`\alpha _s`$ are merged, usually $`\mu _I=2`$ GeV.

For some event shape observables the power corrections have been calculated up to two loops . The result is that hadronisation is described by a shift of the perturbative prediction $`dR_{PT}/dy`$ inversely proportional to the hard scale $`\sqrt{s}=Q`$:

$$\frac{dR}{dy}=\frac{dR_{PT}}{dy}(yPD_y).$$
(2)
The factor $`P\mu _I\alpha _0(\mu _I)/Q`$ is universal  and includes two loop corrections in the quantity $`=1+(1.575C_A0.104n_f)/\beta _0`$ . The factor $`D_y`$ contains all dependencies on the specific event shape observable in question. The prediction for the 1st moment from equation (2) is $`y=y_{PT}+PD_y`$. We note that the DMW model predicts that $`\alpha _0`$ does not depend on the observable. This property of universality can be tested experimentally.

### 2.2  Single dressed gluon approximation

A theoretical analysis of power corrections of the 1st and higher moments of the thrust distribution was presented in . A different approach compared to the DMW model is taken by adding a resummation based on the so-called single dressed gluon approximation (SDG) to the standard perturbative QCD predictions in $`𝒪(\alpha _s^2)`$. The SDG regularises the infrared renormalon singularities.

For the first moment of thrust a power correction proportional to $`1/Q`$ is found, in agreement with the DMW model. The additional resummation leads to a substantial reduction in the relative magnitude of the power correction. It should be noted that power corrections calculated in this approach cannot be directly compared to the DMW model .

For the 2nd and 3rd moment of thrust the power corrections are predicted to be of the form $`1/Q^3`$. However, other terms of the form $`1/Q^2`$ and $`\alpha _s/Q`$ are also predicted to exist but cannot be quantified.

### 2.3  Shape functions

This approach is based on factorisation of soft and hard processes in event shape distributions . The Ansatz is a convolution of the resummed perturbative QCD calculation ($`𝒪(\alpha _s^2)`$+NLLA) $`R_{PT}(y)`$ for the hard process with a so-called shape function $`f(ϵ)`$ for the non-pertubative or soft processes. A factorisation scale $`\mu `$ separates contributions from gluons with $`k_t<\mu `$ taken account of by $`f(ϵ)`$ from those with $`k_t>\mu `$ included in $`R_{PT}(y)`$. For a cumulative event shape distribution $`R(y)`$ at cms energy $`Q`$ one has:

$$R(y)=_0^{yQ}f(ϵ,\mu )R_{PT}(y\frac{ϵ}{Q},\mu )𝑑ϵ.$$
(3)
It is shown that the shape function $`f(ϵ)`$ resums all power corrections of the form $`1/(yQ)^n`$ and that it is related to the energy flow in the event.

The description of non-perturbative effects by a shift $`1/Q`$ of the perturbative QCD prediction in the DMW model is reproduced as the special case where $`f(ϵ)`$ is a delta function. For the 2nd moments of thrust, heavy jet mass and C-parameter a prediction for the power correction of the form $`y_{PT}\lambda _1/Q+\lambda _2/Q^2`$ is given, consistent with section 2.2. Comparison of the predictions with distributions, 1st and 2nd moments of thrust, heavy jet mass and C-parameter from 35 GeV to 189 GeV shows reasonable agreement with the data.

## 3  STANDARD ANALYSES

All experimental analyses study power corrections in the DMW model. For moments the perturbative QCD prediction is $`𝒪(\alpha _s^2)`$ while for distributions the power correction terms are combined with $`𝒪(\alpha _s^2)`$+NLLA calculations.

### 3.1  First and second moments

Studies of power corrections using the 1st moments of the event shape observables thrust, heavy jet mass, total and wide jet broadening and C-parameter have been performed by several collaborations. The LEP experiments ALEPH , DELPHI  and L3  used data recorded at cms energies up to 205 GeV while the reanalysis of data from the PETRA experiment JADE covers in addition the range 35 GeV to 44 GeV .

Figure 1 presents the DELPHI high energy data together with data from other experiments at $`\sqrt{s}<M_{\mathrm{Z}^0}`$. The fit of the $`𝒪(\alpha _s^2)`$ perturbative QCD prediction combined with a power correction describes the data well.

Figure 2 shows the results by L3 for 2nd moments from analyses of the LEP 2 high energy data as well as from low energy data obtained by selecting events at $`\sqrt{s}=M_{\mathrm{Z}^0}`$ with hard initial or final state photon radiation. The fits of the energy evolution of the 2nd moments are based on a simple extension of the DMW model with $`y^2=y^2_{PT}+2y_{PT}PD_y+A_2/Q^2`$. The quantity $`A_2`$ is fitted with $`\alpha _s`$ and $`\alpha _0`$ fixed to the results from fits to 1st moments and it is found to be consistent with zero only for heavy jet mass and the wide jet broadening.

### 3.2 Differential distributions

In the DMW model power corrections to differential event shape distributions are described by a shift $`1/Q`$ of the perturbative prediction. This prediction has been tested with reanalysed JADE data together with data from SLD and the LEP experiments  and also with LEP data by ALEPH . Figure 3 presents distributions of the C-parameter from $`\sqrt{s}=35`$ GeV to 183 GeV with a fit of the DMW prediction . The fitted prediction within the fit ranges indicated by solid lines describes the data well within the errors.

## 4 STUDY OF QCD COLOUR FACTORS

In this analysis fits of the DMW model to event shape distributions as shown in section 3.2 are generalised to vary one of the QCD colour factors $`n_f`$, $`C_A`$ or $`C_F`$ along with $`\alpha _s`$ and $`\alpha _0`$ . Sensitivity to the colour factors comes mainly from the running of $`\alpha _s`$ in the perturbative prediction. Potential biases from assuming standard QCD in the hadronisation corrections should not be present. It is observed that fits to thrust and C-parameter are stable while fits to total and wide jet broadening don’t converge well in all cases. Results of the fits to individual colour factors, $`\alpha _s`$ and $`\alpha _0`$ are shown in figure 4 and are consistent with standard QCD based on the SU(3) symmetry group. The total uncertainties on the measurements of the colour factors are competitive with traditional analyses using 4-jet final states at $`\sqrt{s}=M_{\mathrm{Z}^0}`$.

## 5 SUMMARY

A summary of all results for $`\alpha _s(M_{\mathrm{Z}^0})`$ and $`\alpha _0`$ of analyses based on the DMW model is shown in figure 5. The label PPSTL refers to the analyses of data from PEP, PETRA, SLC, TRISTAN and LEP performed together with reanalysed JADE data . The results for $`\alpha _s(M_{\mathrm{Z}^0})`$ are consistent with the current world average $`\alpha _s(M_{\mathrm{Z}^0})=0.118\pm 0.003`$  while the results for $`\alpha _0`$ mostly cluster within $`\alpha _00.5\pm 25\%`$. The results for $`\alpha _0`$ support the prediction of universality within 25%.

The other power corrections models are promising new approaches and need experimental testing. From an experimental point of view it would also be desirable to quantify the relations between the different models.

The study of QCD colour factors within the DMW model yields results consistent with standard SU(3) QCD with thrust and C-parameter. The total uncertainties on the colour factors are competitive with other measurements of the colour factors.

The field of power corrections in $`\mathrm{e}^+\mathrm{e}^{}`$ annihilation is rapidly developing and it gives one the opportunity to learn more about soft QCD and the physics of hadronisation.

## QUESTIONS

Is it normal that your analysis \[of the QCD colour factors\] works for both 1-T and C because they are almost 100% correlated? Do you know any other event shape variable for which it works?

It is consistent that the analysis works for $`1T`$ and $`C`$, because they are correlated. However, the observables $`B_T`$ and $`B_W`$ are also highly correlated with $`1T`$ and $`C`$ but the fits are not always stable. This may indicate that the DMW model is not as good an approximation for $`B_T`$ and $`B_W`$ as for $`1T`$ and $`C`$. Other event shapes for which the analysis would be possible are heavy jet mass and energy-energy-correlation.

What are the indications of an insufficient perturbative description of the wide jet broadening \[in the colour factor analysis\]?

The analysis of the colour factors used the latest corrected calculations for the jet broadening observables. It is not known if the perturbative or the power correction part of the prediction for $`B_W`$ is insufficient.


text:
# CERN-TH 2000-232 IFIC/00-49 hep-ph/0009053 Monte Carlo simulation for jet fragmentation in SUSY-QCD

## 1 Introduction

QCD, being an essential part of the Standard Model, successfully describes accelerator data for production of hadrons in $`e^+e^{}`$ annihilation and in deep-inelastic scattering. There are two distinctive parts in these calculations: the perturbative QCD computation of the parton cascade in jets, and the parton hadronization, in which low-virtuality partons are converted non-perturbatively into hadrons.

The QCD parton cascade is usually studied in Modified Leading Logarithmic Approximation (MLLA), where large logarithms, $`\mathrm{ln}(Q^2)`$ and $`\mathrm{ln}(1/x)`$, play a crucial role (here $`Q^2`$ is the maximum of the perpendicular momentum $`k_{}`$, and $`x=k_{}/k_{}^{\mathrm{max}}`$).

This approximation is characterized by remarkable features.

In the MLLA the QCD cascade has a probabilistic interpretation, provided by the absence of interference terms in the tree diagrams. The colour coherence effect is taken into account in the MLLA. It suppresses the emission of soft gluons and results in the Gaussian peak of the parton distribution in terms of $`\xi =\mathrm{ln}(1/x)`$ (hump-backed plateau).

The evolution of parton cascades in the MLLA (as well as in the LLA) is adequately described by the Dokshitzer–Gribov–Lipatov–Altarelli–Parisi (DGLAP) equations .

The parton spectra can be obtained analytically and by Monte Carlo (MC) simulations.

Examples for analytical solutions are the limiting spectrum  and the Gaussian spectrum , in which we include the distorted Gaussian spectrum . The limiting spectrum is the most accurate of them. Since we have a special interest in it, we shall shortly review below the basic assumptions under which this solution is obtained.

The limiting spectrum gives the energy spectrum of partons $`D_{\mathrm{lim}}(\xi ,Y)`$ for a given center-of-mass energy $`\sqrt{s}`$ of an $`e^+e^{}`$-pair. Here $`D`$ is the number of cascade partons, $`Y=\mathrm{ln}(\sqrt{s}/2\mathrm{\Lambda })`$ and $`\mathrm{\Lambda }`$ is the dimensional QCD scale.

The analytical expression for $`D_{\mathrm{lim}}(\xi ,Y)`$ is given in Refs. . There are two fundamental parameters involved in the limiting spectrum solution: the scale $`\mathrm{\Lambda }`$ and the minimal virtuality $`Q_0`$ of partons, down to which the cascade develops perturbatively; $`Q_0`$ can be viewed as the effective mass of the partons. Two assumptions are necessary for the validity of the limiting spectrum. The QCD coupling constant $`\alpha _s(k_{}^2)`$ evolves with $`k_{}^2`$ effectively as in the one-loop approximation with three flavors $`n_f=3`$ for all $`k_{}^2`$,

$$\alpha _s(k_{}^2)=\frac{12\pi }{(332n_f)\mathrm{ln}(k_{}^2/\mathrm{\Lambda }^2)}.$$
(1)
As a matter of fact, $`\mathrm{\Lambda }`$ in Eq. (1) is treated in the limiting spectrum solution as a free parameter to fit the $`e^+e^{}`$-data. The best fit corresponds to $`\mathrm{\Lambda }=`$ 250 – 270 MeV. For this range of $`\mathrm{\Lambda }`$ values, $`\alpha _s(M_Z)`$, given by Eq. (1) with $`n_f=3`$ is in the interval 0.118–0.120, to be compared with the average experimental value $`\alpha _s(M_Z)=0.1184\pm 0.0031`$ . Therefore, the phenomenological parameter $`\mathrm{\Lambda }`$ coincides well with $`\mathrm{\Lambda }_{\mathrm{QCD}}`$, which fits the experimental value of $`\alpha _s(M_Z)`$, in the one-loop approximation with $`n_f=3`$.

The second assumption, necessary for the derivation of the limiting spectrum, is $`Q_0=\mathrm{\Lambda }`$. It gives a reasonable value of $`Q_0`$, but the exact equality of these values has no theoretical justification.

The limiting spectrum solution is valid only for small $`x1`$. In this region, which includes the maximal values of multiplicity (in the Gaussian peak), it describes very accurately experimental data at all available energies $`\sqrt{s}`$ (see, e.g., ). The large $`x`$ up to $`x=1`$ give the dominant contribution to the total momentum of cascade partons. Therefore, the limiting spectrum solution does not guarantee that $`xD_{\mathrm{lim}}(x,s)dx`$ precisely equals 2, although it can be valid up to large $`x`$ .

Monte Carlo (MC) simulations of the QCD cascade give a more precise description of the cascade evolution. They are valid for all $`x`$ including $`x1`$. In contrast to the limiting spectrum, in MC simulations one can use $`\alpha _s(k_{}^2)`$ with the measured value of $`\mathrm{\Lambda }_{\mathrm{QCD}}`$, varying the number of flavours and two-loop corrections beeing taken into account. The assumption $`\mathrm{\Lambda }=Q_0`$, specific to the limiting spectrum, is not needed. MC simulations are based on a probabilistic interpretation of the jet cascade. Parton branching is described by the Altarelli–Parisi functions, and the probability of parton evolution between two values of virtualities without branching is given by the Sudakov form factor. Finally, the coherent effect in the soft gluon emission (destructive interference) is conveniently taken into account by angular ordering $`\theta _1>\theta _2>\theta _3\mathrm{}`$ , where the indices number the generations (non-ordered processes are suppressed ). The first MC simulation with angular ordering was performed in Ref. . At present there are several detailed MC simulations, e.g. , which differ mainly in their description of hadronization. We shall now briefly discuss the problem of hadronization.

The description of parton hadronization is based on the assumption of Local Parton–Hadron Duality (LPHD) . This hypothesis implies that when $`Q_0`$ is small enough (of the order of $`\mathrm{\Lambda }`$) there is a proportionality between the spectra of partons and hadrons, with relations between their momenta, which are local in the phase space. Such an interpretation can be justified by the idea of preconfinement .

As far as spectra are concerned, LPHD implies a proportionality between the hadron and parton spectra. In Refs. , it is emphasized that, most reasonably, this proportionality holds not on a “one parton – one hadron” basis, but for the number of particles averaged over a finite interval $`\mathrm{\Delta }\xi 1`$.

The LPHD hypothesis for limiting spectrum straightforwardly results  in

$$D_{\mathrm{had}}(x,\sqrt{s})=K_h(Q_0)D_{\mathrm{part}}(x,\sqrt{s},Q_0),$$
(2)
where the constant $`K_h`$ is universal, in the sense that it does not depend on $`\sqrt{s}`$.

Equation (2) completes our description of the limiting spectrum, expressing the hadron spectrum through the spectrum of partons. The constant $`K_h`$, which connects the two spectra, is found from a comparison with experimental data as $`K_h1.3`$ for $`\mathrm{\Lambda }=Q_0270`$ MeV , and it does not change with energy unless some new physics (e.g. supersymmetry) appears.

In MC simulations, the parameter $`Q_0`$ is in principle a free parameter found by fitting experimental data. For HERWIG  and PYTHIA , for example, $`Q_01`$ GeV. Several detailed hadronization models are used in simulations, e.g. the independent fragmentation model, the Lund string model , and the cluster fragmentation model . Usually, these models use many free parameters and require to keep track of the four-momentum evolution of all partons.

The calculations described above are valid up to $`\sqrt{s}1`$–10 TeV. At higher energies the production of supersymmetric particles is expected to change the results. One might be interested in much higher energies, being inspired by the production of superheavy particles up to the GUT scale in the Universe. Such particles can be produced by topological defects and by many processes at the post-inflationary stage of the Universe. Recently superheavy particles with masses $`M_X10^{12}`$$`10^{14}`$ GeV attracted much attention as a source of the observed Ultra High Energy Cosmic Rays (UHECR) with energies $`10^{19}`$$`10^{20}`$ eV (for recent reviews, see ).

The limiting spectrum for SUSY QCD was calculated in Ref.  for very high energies $`\sqrt{s}`$, corresponding to masses of superheavy particles $`M_X10^{12}`$$`10^{14}`$ GeV. The supersymmetric partons (squarks and gluinos or jointly spartons) participate in the cascade until the virtualities $`t`$ of the particles drop below the mass scale of SUSY particles, $`tM_{\mathrm{SUSY}}^2`$. Then a SUSY particle decays, producing in the end the Lightest Supersymmetric Particle (LSP), for which the lightest neutralino is usually considered. The role of supersymmetric partners is two-fold: they double the number of parton types in the cascade, and they change the evolution of $`\alpha _s(k_{}^2)`$. Even at small $`tM_{\mathrm{SUSY}}^2`$, the cascade remembers the number of flavours at large $`t`$ because, for example, each squark leaves after its decay a quark, which continues QCD cascading. At large $`t`$ and small $`x1`$ gluons and gluinos dominate and their “children” constitute the dominant part of the cascade at small $`t`$. Therefore, the dominant contribution to the limiting spectrum is given by gluons and gluinos.

The SUSY QCD limiting spectrum solution has two drawbacks with respect to ordinary QCD. First, the number of flavours that determine the evolution of the coupling constant according to Eq. (1) has to be fixed to one value of $`n_f`$ for the whole range of $`k_{}^2`$. Second, the limiting spectrum for ordinary QCD is normalized by experimental data, which are absent in the case of SUSY QCD. Normalization due to the conservation of momentum $`xD_{\mathrm{lim}}(x,s)dx=2`$ is unreliable, since the limiting spectrum is not valid for large $`x`$, which give the main contribution to the integral (see discussion in ).

During the last few years, the production and decays of supersymmetric particles have been included in most MC simulations focusing on LHC studies. Although the LHC will operate above the expected threshold of SUSY particle production, its energy is not large enough for these particles to participate in the QCD cascade. Therefore, all currently available MC simulations consider only on-shell decays of spartons and neglect possible branchings of gluinos and squarks111The future C++ version of HERWIG will include branchings of spartons .. Another obstacle against the use of standard MC simulations at extremely large energies around $`\sqrt{s}10^{12}10^{14}`$ GeV is that the necessary numerical precision and the required amount of memory space and computing time become a challenge for present–day computers.

We have therefore developed a new MC simulation, which includes as cascading particles not only gluons and quarks but also gluinos and squarks. We consider cascades that are initiated by the decay to two jets of superheavy particles with mass $`M_X10^{12}`$$`10^{14}`$ GeV, or by $`e^+e^{}`$ annihilation at $`s=M_X^2`$. SUSY partons, squarks and gluinos, are produced in the fragmentation of ordinary partons and vice versa. All squarks and gluinos are assumed to have equal masses, for which we use $`M_{\mathrm{SUSY}}=200`$ GeV and $`M_{\mathrm{SUSY}}=1000`$ GeV. When the virtuality of the cascading particles drops below $`M_{\mathrm{SUSY}}^2`$, spartons decay to LSPs (neutralinos), which freely escape. The perturbative development of the cascade continues with ordinary partons until their virtualities reach $`Q_0^2`$, for which we use $`Q_0^2=0.625`$ GeV2 to fit the data at small energies $`\sqrt{s}`$. We use a new hadronization procedure. It is based on a model-independent, phenomenological approach, in which hadronization functions for large $`s`$ (or $`M_X`$) are calculated from hadron spectra observed at small $`s`$ ($`M_X`$). This method can be used for any type of hadrons as well as for photons and neutrinos, if their spectra are known with good enough accuracy at small energy.

Following Ref. , we shall use the following notation:
$`\mathrm{\Lambda }`$ is the dimensional QCD scale,
$`Y=\mathrm{ln}(\sqrt{s}/2\mathrm{\Lambda })`$,
$`t=p_\mu ^2`$ is the virtuality of cascade partons,
$`Q^2=t_{\mathrm{max}}`$ is the virtuality of the primary parton,
$`Q_0^2`$ is the minimum virtuality of the perturbative evolution of the QCD cascade,
$`z=E^{}/E`$, where $`E`$ and $`E^{}`$ are the energies of ingoing and outgoing partons at fragmentation,
$`\zeta =1\mathrm{cos}\theta `$, where $`\theta `$ is the angle between two outgoing partons,
$`\stackrel{~}{t}=\zeta E^2`$,
$`k_{},k_{}`$ are the transverse and parallel momenta transferred, respectively,
$`x=k_{}/k_{}^{\mathrm{max}}`$,
$`\xi =\mathrm{ln}(1/x)`$.

## 2 MC Simulation of the Perturbative Phase of SUSY QCD Cascades

The perturbative part of our simulation is very similar to those of MCs for ordinary QCD cascades, except for including spartons and the condition for their exit from the cascade. We consider a superheavy $`X`$ particle with mass $`M_X`$ which decays into two jets with energy $`E_{\mathrm{jet}}=M_X/2`$. We assume that the primary partons produced in the $`X`$ particle decay have the maximum virtuality $`Q^2=m_X^2/4`$ and that the $`X`$ particle has equal branching ratios to all partons. As to the first assumption, in reality, there is a distribution of partons with different $`t`$, but the Sudakov form factors suppress small $`t`$ values. The second assumption is made because of the unspecified interactions of the $`X`$ particles.

Our simulation closely follows the angular ordered parton cascade algorithm developed in Refs. . It is convenient to use in this algorithm the variable $`\stackrel{~}{t}=\zeta E^2`$, where $`E`$ is the energy of the incoming parton, $`\zeta 1\mathrm{cos}\theta `$, and $`\theta `$ is the angle between the two emitted partons. A primary parton with energy $`E_{\mathrm{jet}}`$ ($`=m_X/2`$) and angular variable $`\xi _01`$ initiates a cascade, which proceeds until the ordinary partons reach the minimal virtuality $`\stackrel{~}{t}=4Q_0^2`$. Here the perturbative evolution of the cascade terminates.

In each branching of an incoming parton $`i`$ with $`\stackrel{~}{t}^{}`$, we generate with the veto algorithm  a new $`\stackrel{~}{t}`$ and $`z`$ according to the probability distribution

$$\mathrm{d}𝒫_i(\stackrel{~}{t},z)=\underset{jk}{}\frac{\mathrm{d}\stackrel{~}{t}}{\stackrel{~}{t}}\frac{\mathrm{d}z}{2\pi }\alpha _s\left[z^2(1z)^2\stackrel{~}{t}\right]P_{ijk}(z)\frac{\mathrm{\Delta }_i(\stackrel{~}{t}^{})}{\mathrm{\Delta }_i(\stackrel{~}{t})}.$$
(3)
Here, the sum includes all possible branching channels $`jk`$, $`z^2(1z)^2\stackrel{~}{t}`$ is the parton transverse momentum, and $`\mathrm{\Delta }_i`$ is the product of the individual Sudakov-like form factors $`\mathrm{\Delta }_{ijk}`$  ,

$$\mathrm{\Delta }_{ijk}(\stackrel{~}{t})=\mathrm{exp}\left[_{4\stackrel{~}{t}_{\mathrm{min}}}^{\stackrel{~}{t}}\frac{\mathrm{d}t^{}}{t^{}}f_{ijk}(t^{})\right]$$
(4)
with

$$f_{ijk}(\stackrel{~}{t})=_{z_{\mathrm{min}}}^{z_{\mathrm{max}}}\frac{\mathrm{d}z}{2\pi }\alpha _s\left[z^2(1z)^2\stackrel{~}{t}\right]P_{ijk}(z).$$
(5)
The unregularized Altarelli–Parisi splitting functions $`P_{ijk}(z)`$ of SUSY QCD  are given in Table 1.

The angular ordering $`\zeta _j,\zeta _k<\zeta _i`$ for the branching $`ijk`$, which takes into account colour coherence, is equivalent to $`\stackrel{~}{t}_j<z^2\stackrel{~}{t}_i`$ and $`\stackrel{~}{t}_k<(1z)^2\stackrel{~}{t}_i`$. These conditions result in

$$z_{\mathrm{min}}=\sqrt{\stackrel{~}{t}_{\mathrm{min}}/\stackrel{~}{t}},z_{\mathrm{max}}=1\sqrt{\stackrel{~}{t}_{\mathrm{min}}/\stackrel{~}{t}}.$$
(6)
For the evolution of the running coupling $`\alpha _s`$ as a function of gluon virtuality $`t`$ at small momentum transfer $`t<t_{\mathrm{SUSY}}`$, we use the standard two-loop dependence with variable $`n_f`$ and thresholds, and normalize $`\alpha _s`$ as $`\alpha _s(M_Z)=0.119`$, which corresponds to $`\mathrm{\Lambda }_{\overline{MS}}^{(5)}=222`$ MeV. At large momentum transfer $`t>t_{\mathrm{SUSY}}`$ we use minimal SUSY-SU(5) coupling constant evolution , normalizing the coupling constant at $`\sqrt{t}=M_{\mathrm{GUT}}=110^{16}`$ GeV, as $`\alpha _s(M_{\mathrm{GUT}}^2)1/25.8`$. Explicitly we use

$$\alpha _s(t)=\frac{\alpha (M_{\mathrm{GUT}}^2)}{1+b_s/(4\pi )\mathrm{ln}(t/M_{\mathrm{GUT}}^2)\alpha (M_{\mathrm{GUT}}^2)},$$
(7)
where $`b_s=9n_f`$ is a constant, that governs the evolution of the coupling constant with $`t`$. At $`t>t_{\mathrm{SUSY}}`$, $`n_f=6`$ and $`b_s=3`$. The above assumption means that we introduce, instead of many thresholds corresponding to SUSY particles with different masses, a single threshold at $`t=t_{\mathrm{SUSY}}`$. This is a reasonable thing to do in view of the large uncertainties in our knowledge of mass spectrum of SUSY particles. Equation (7) approximates accurately enough the evolution of $`\alpha _s(t)`$ as calculated in Ref. , when $`t_{\mathrm{SUSY}}210^5`$GeV2. Starting from this value, $`\alpha _s(t)`$ evolves in the regime of Eq. (7). Note that $`t_{\mathrm{SUSY}}`$ does not necessarily coincide with the scale $`M_{\mathrm{SUSY}}`$, the universal mass of squarks and gluino, for which we use as two representative values $`M_{\mathrm{SUSY}}=200`$ GeV and $`M_{\mathrm{SUSY}}=1`$ TeV. In particular, the low value of $`t_{\mathrm{SUSY}}`$ used here is compatible with much larger $`M_{\mathrm{SUSY}}`$, as emphasized in Ref. .

Finally, we have to specify the value of the cut-off $`\stackrel{~}{t}_{\mathrm{min}}`$ for the cascade evolution. We do not distinguish between different quark flavours and use $`\stackrel{~}{t}_{\mathrm{min}}=0.625`$ GeV2 for all branchings in which only ordinary particles are produced and $`\stackrel{~}{t}_{\mathrm{min}}=M_{\mathrm{SUSY}}^2`$, where $`M_{\mathrm{SUSY}}`$ is the typical mass scale of the spartons, for branchings in which SUSY particles are produced, respectively.

Let us now describe a step $`ijk`$ in our simulation. For an incoming parton $`i`$ with $`\stackrel{~}{t}^{}`$ we generate first a new cascade variable $`\stackrel{~}{t}`$, according to the probability distribution given by the ratio $`\mathrm{\Delta }_i(\stackrel{~}{t}^{})/\mathrm{\Delta }_i(\stackrel{~}{t})`$. Then we select the branching channel $`jk`$ using as weight $`f_{ijk}(\stackrel{~}{t})`$ and generate $`z`$ according to the probability distribution $`\alpha _s\left[z^2(1z)^2\stackrel{~}{t}\right]P_{ijk}(z)`$.

The last ingredient in the perturbative part of our simulation is the exit of supersymmetric particles from the cascade. We assume that the neutralino $`\stackrel{~}{\chi }`$ is the LSP and that R-parity is conserved. Reaching $`\stackrel{~}{t}_{\mathrm{min}}=M_{\mathrm{SUSY}}^2`$, squarks and gluinos decay as $`\stackrel{~}{q}q+\stackrel{~}{\chi }`$ and $`\stackrel{~}{g}q+\overline{q}+\stackrel{~}{\chi }`$, thus producing UHE LSPs.

In fact, we are running in this work two Monte Carlo similations: with the ordinary QCD and with SUSY QCD.

In the former case supersymmetric partons are not included, and for perturbative calculations we assume the SM particle content with $`\alpha _s(t)`$ evolution in two-loop approximation with proper thresholds. We fix $`Q_0^2=0.625`$ GeV2. We need these calculations mostly to test our method.

The assumptions of SUSY QCD Monte Carlo are described above. At $`\stackrel{~}{t}_{min}<M_{SUSY}^2`$ cascade develops according to ordinary QCD scheme.

## 3 Hadronization

The Monte Carlo simulation described in the last section is completely determined by perturbative physics. How the spectrum of coloured quarks and gluons $`D_i(x,\sqrt{s})`$ is transformed into the spectrum of hadrons $`D_{\mathrm{had}}(x,\sqrt{s})`$ is still an open problem. Monte Carlo simulations have to use some hadronization model (see Introduction), which describes the non-perturbative evolution of the cascade for $`\stackrel{~}{t}<4\stackrel{~}{t}_{\mathrm{min}}`$. Two hadronization models, the cluster fragmentation model  used in HERWIG  and the Lund string model  used in PYTHIA , require the knowledge of the four-momenta of all the partons. Thus these models need detailed time and memory-consuming computations.

We suggest here a phenomenological, model-independent hadronization scheme based on the knowledge of the hadron spectra at energies $`\sqrt{s}`$ smaller than the energy of interest. This method is valid for any hadron type and can be applied to the secondary particles, such as photons and neutrinos, as well. The application of this method is somewhat restricted (e.g. it cannot give the angular distribution of particles in a jet or correlations), but its use is very efficient for the decay of superheavy particles, where multiplicity, and hence the number of partons to follow in a simulation is very large.

Our hadronization scheme depends on only one theoretical assumption, which is reliable and testable. Namely, we assume that the unknown non-perturbative physics can be factorized into hadronisation functions $`f_i(z)`$ that do not depend on $`\sqrt{s}`$,

$$D_h(x,\sqrt{s})=\underset{i=q,g}{}_x^1\frac{\mathrm{d}z}{z}D_i(x/z,\sqrt{s})f_i^h(z),$$
(8)
where the index $`h`$ runs through different types of hadrons, e.g. $`\pi ^0,\pi ^\pm ,N`$ etc.

The functions $`f_i^h(z)`$ give the probability that a parton $`i`$ with energy $`E`$ is converted into a hadron $`h`$ with energy $`zE`$. It is implicitly assumed in Eq. (8) that the perturbative cut-off $`Q_0`$ is fixed, and $`f_i^h`$ is determined for this value of $`Q_0`$, although in principle for every $`Q_0`$ and $`D_i(x,\sqrt{s},Q_0)`$ one can find $`f_i^h(z,\sqrt{s},Q_0)`$ to fit the observed hadron spectra.

Equation (8) with energy independent hadronization functions follows from basic principles and is confirmed (see below) at energies of $`e^+e^{}`$ colliders. It has the form of a Volterra integral equation of the first kind though in contrast to the standard case, the RHS contains not one but two unknown functions $`f_i`$ for every $`h`$. In principle, the two functions $`f_g(x)`$ and $`f_q(x)`$ can be uniquely determined if $`D_h`$ is known as an analytic function without errors for two different values of $`\sqrt{s}`$. In practice, $`D_h(x)`$ is known only as a discrete set of experimental data and Eq. (8) represents an ill-posed inversion problem222Volterra integral equations of the first kind can be solved normally by linearization, even if the LHS are data. However, the lower integration limit in (8) does not represent a sharp cut-off because the kernels $`D_i(x)`$ vanish for $`x1`$. Therefore, Eq. (8) behaves effectively like a Fredholm equation, and these are known to be extremely ill-conditioned. . Instead of solving Eq. (8) by an inversion method, we prefer to find physically motivated trial functions for $`f_i`$ to fit the experimental data at $`\sqrt{s}=91.2`$ GeV.

In terms of the more convenient variable $`\xi =\mathrm{ln}(1/x)`$, Eq. (8) has the form

$$D_h(\xi ,\sqrt{s})=\underset{i=q,g}{}_0^ld\xi ^{}D_i(\xi \xi ^{},\sqrt{s})f_i(\xi ^{}),$$
(9)
where the index $`h`$ in the hadronization functions is suppressed.

In the limiting spectrum, when $`Q_0=\mathrm{\Lambda }`$, the hadronization functions $`f_i`$ are proportional to delta functions. Inspired by this analytical solution we choose for $`f_i`$ Gaussian functions,

$$f_i(\xi )=a_i\mathrm{exp}\left(\frac{(\xi \xi _{\mathrm{max},i})^2}{\sigma _i^2}\right).$$
(10)
With this hadronization function the approximate proportionality holds between spectra of partons and hadrons as LPHD demands. The position of the peak in the hadronization function determines the shift between the maxima of parton and hadron spectra. While for gluons the hadronization function $`f_g(\xi )`$ should vanish for $`\xi 0`$, because gluons have to split their energy to a $`q\overline{q}`$ pair, for quarks $`f_q(\xi )`$ can be finite at $`\xi =0`$.

The hadronization functions we obtained for $`Q_0^2=0.625`$ GeV2 from a fit to LEP data at $`\sqrt{s}=91.2`$ GeV are shown in Fig. 1.

Our hadronization scheme has been tested by two methods: for relatively small energies, $`\sqrt{s}=58`$ GeV and 133 GeV, we confronted our calculations with LEP data, and for very large $`\sqrt{s}`$ (or $`M_X`$) we compared the calculated spectrum with the limiting spectrum, using a special case when it is correct (see below). In both cases ordinary QCD Monte Carlo was used.

Figures 2-4 display a comparison between the charged hadron spectrum from our MC simulation for ordinary QCD and experimental data  at $`\sqrt{s}=58`$, 91.2 and 133 GeV, respectively.

Let us now discuss whether the hadronization functions $`f_i(\xi )`$ found from the fit to data at $`\sqrt{s}=91.2`$ GeV can really be used at $`M_X=10^{12}`$$`10^{16}`$ GeV.

First of all we note, that a test can be given by LPHD, which demands approximate proportionality between parton and hadron spectra. It implies that the $`\xi ^{}`$ values, that give the dominant contribution to the integral in Eq. (9), are about the same at $`\sqrt{s}=91.2`$ GeV and at large $`M_X`$. Numerical tests show that this is indeed the case for both the quark and the gluon contribution.

As a critical test of our hadronization scheme, we compared the limiting spectrum with the results of our simulation for a special, well-controlled case of ordinary QCD with the number of quark flavours $`n_f=3`$, and with $`\alpha _s(k_{}^2)`$ given by Eq.(1). For the limiting spectrum in this case we can use the normalization constant $`K_h1.3`$ obtained by fitting experimental data . Since we do not introduce new high-energy physics, the limiting spectrum is valid for any initial energy $`M_X`$ (see Introduction). In Fig. 5, we show the ratio of these two (charged) hadron spectra. The agreement between the two spectra is excellent, except for the small $`\xi <6`$ region where it is known that the limiting spectrum is not valid. The disagreement reaches $`50\%`$ at $`\xi 2.1`$ ($`x0.12`$).

In conclusion, we think that our hadronization recipe is a valid alternative to the extrapolation of the Lund string or the cluster fragmentation model to extremely large $`M_X`$.

## 4 Results: Spectra of Hadrons and Secondary Particles

Using the algorithm for the perturbative evolution of the SUSY QCD cascade as described in Section 2 and our hadronization scheme from Section 3, we can now compute the fragmentation spectra of hadrons. As numerical values for $`M_X`$, we choose in the graphs given as examples three values interesting for UHECR physics, $`M_X=10^{12},10^{13},10^{14}`$ GeV, as well as $`M_X=10^6`$ and $`10^{16}`$ GeV as lowest and highest scale of interest. Similarly, we use $`M_{\mathrm{SUSY}}=200`$ GeV and $`M_{\mathrm{SUSY}}=1000`$ GeV as two representative values for the SUSY mass scale.

In Figs. 8 and 9, the hadron spectra $`\mathrm{d}N_{\mathrm{had}}/\mathrm{d}\xi `$ from SUSY QCD MC simulations are displayed as function of $`\xi `$ for $`M_{\mathrm{SUSY}}=200`$ GeV and $`M_{\mathrm{SUSY}}=1000`$ GeV, respectively; in both figures the spectra were calculated for $`M_X=10^{12},10^{13},10^{14}`$ GeV. For the GUT scale $`M_X=10^{16}`$ GeV, the hadron spectra $`\mathrm{d}N_{\mathrm{had}}/\mathrm{d}\xi `$ are shown in Fig. 10 and for the low scale $`M_X=10^6`$ GeV in Fig. 11. The hadron spectra depend only weakly on $`M_{\mathrm{SUSY}}`$, with increasing differences for larger values of $`M_X`$. Both effects are easy to understand: when spartons disappear from the cascade at $`\stackrel{~}{t}M_{\mathrm{SUSY}}^2`$ due to on-shell decays, each of them leaves there an ordinary parton with similar virtuality. Therefore, the cascade proceeds as if nothing had happened, except that some energy is lost through the emission of neutralinos and leptons, which is not large ($`10\%`$). Second, the importance of spartons for the cascade decreases for smaller values of $`M_X`$, thereby reducing also the dependence of the hadron spectra on $`M_{\mathrm{SUSY}}`$ for smaller $`M_X`$.

The weak dependence of the spectra on $`M_{SUSY}`$ justifies our choice of an universal value for the masses of supersymmetric particles. Indeed, if we assume now that supersymmetric particles have different masses in the range 200-1000 GeV, the resulting hadron spectra will differ less than in Figs. 10 and 11.

The signature of supersymmetry in decays of superheavy $`X`$ particles is the production of LSPs, which we assume as stable neutralinos. They are generated in the cascade mostly when the virtuality of the spartons approaches $`M_{\mathrm{SUSY}}^2`$. The calculated neutralino spectra are shown in Figs. 1215 for the same parameters as the hadron spectra in Figs. 811. Like the hadron spectra, they have the characteristic Gaussian form, however with a shifted position of their maxima due to their larger cut-off $`M_{\mathrm{SUSY}}`$ in the shower development. The energy fraction taken away by the neutralinos is typically $`10\%`$ for values of $`M_X`$ interesting for UHECR physics, with a minimum of $`5\%`$ for $`M_X=10^6`$ GeV and $`M_{\mathrm{SUSY}}=1`$ TeV and a maximum of $`12\%`$ for $`M_X=10^{16}`$ GeV and $`M_{\mathrm{SUSY}}=200`$ GeV.

We have only derived a common hadronization function for all hadrons and, consequently, we cannot calculate directly, e.g. pion or nucleon spectra through Eq. (8). Since the fraction of energy $`ϵ_i`$ going into different meson and baryon species is determined by the non-perturbative process of hadronization, these fractions as the hadronization functions themselves do not depend on $`s`$. Thus, we can use the value from $`Z`$ decay, $`ϵ_N0.05`$ and $`ϵ_\pi 0.95`$. Then

$$\frac{\mathrm{d}N_{\mathrm{nucl}}}{\mathrm{d}x}=ϵ_N\frac{\mathrm{d}N_{\mathrm{had}}}{\mathrm{d}x},\frac{\mathrm{d}N_\pi }{\mathrm{d}x}=ϵ_\pi \frac{\mathrm{d}N_{\mathrm{had}}}{\mathrm{d}x}.$$
(11)
Using the hadron spectra obtained in the last Section, it is simple to calculate analytically the spectra of secondary particles, photons and neutrinos. The normalized photon spectrum from a decay of one $`X`$ particle at rest is given by

$$\frac{\mathrm{d}N_\gamma }{\mathrm{d}x}=\frac{2}{3}ϵ_\pi _x^1\frac{\mathrm{d}y}{y}\frac{\mathrm{d}N_{\mathrm{had}}}{\mathrm{d}y}.$$
(12)
The total neutrino spectrum, given by the sum from decays of pions and muons, can be presented in the following form,

$$\frac{\mathrm{d}N_\nu }{\mathrm{d}x}=\frac{2}{3}ϵ_\pi \left(\frac{\mathrm{d}N_{\nu _\mu }}{\mathrm{d}x}(\pi \mu \nu _\mu )+\frac{\mathrm{d}N_{\nu _\mu }}{\mathrm{d}x}(\mu \nu _\mu \nu _ee)+\frac{\mathrm{d}N_{\nu _e}}{\mathrm{d}x}(\mu \nu _\mu \nu _ee)\right),$$
(13)
where for pion decay

$$\frac{dN_{\nu _\mu }}{\mathrm{d}x}(\pi \mu \nu _\mu )=R_{Rx}^1\frac{\mathrm{d}y}{y}\frac{\mathrm{d}N_{\mathrm{had}}}{\mathrm{d}y}$$
(14)
and for muon decay

$$\frac{\mathrm{d}N_{\nu _i}}{\mathrm{d}x}(\mu \nu _\mu \nu _ee)=R_x^1\frac{\mathrm{d}y}{y}_y^{y/r}\frac{\mathrm{d}y^{}}{y^{}}\frac{\mathrm{d}N_{\nu _i}}{\mathrm{d}y}\frac{\mathrm{d}N_{\mathrm{had}}}{\mathrm{d}y^{}},$$
(15)
with

$$\frac{\mathrm{d}N_{\nu _e}}{\mathrm{d}y}=26y^2+4y^3,\frac{\mathrm{d}N_{\nu _\mu }}{\mathrm{d}y}=\frac{5}{3}3y^2+\frac{4}{3}y^3,$$
(16)
and $`r=(m_\mu /m_\pi )^2`$, $`R=1/(1r)`$.

The resulting nucleon, photon and neutrino spectra $`x^3\mathrm{d}N_i/\mathrm{d}x`$ are shown as functions of $`x`$ together with the spectra of neutralinos in Fig. 17 for $`M_{\mathrm{SUSY}}=200`$ GeV and $`M_X=10^{12}`$ GeV and $`M_X=10^{14}`$ GeV, respectively. We have multiplied the spectra by $`x^3`$ in order to facilitate the comparison of our spectra with the energy spectra of observed UHECR. At $`x>0.7`$, the spectra have some uncertainties because of the unknown branching ratios of the $`X`$ particle into (s)partons and fluctuations due to the small number of produced particles. The excess of nucleons over secondary particles, neutrinos and gamma, at largest $`x`$ is a result of kinematical effect and the steep spectra of pions and nucleons in the end of the spectrum.

## 5 Discussion

In this Section we compare for large $`M_X`$ the results of our MC for the two cases, SUSY QCD and ordinary-QCD, with other computations, and most notably with limiting spectrum calculations. The latter case of ordinary QCD is formally a special case of our MC simulation for SUSY QCD in the limit $`M_{\mathrm{SUSY}},t_{\mathrm{SUSY}}\mathrm{}`$, i.e. $`\alpha _s`$ is given by two-loop approximation with variable $`n_f`$, and the probability to produce a sparton is zero.

The validity of our method has been proved by the tests described in Section 3. If no new physics beyond the three light quark flavours is introduced, the $`k_{}`$-dependence of $`\alpha _s`$ is given by Eq. (1) with $`n_f=3`$ and the limiting spectrum with $`K_h=1.3`$ is valid for arbitrary high energies. We can calculate the hadron spectrum in our ordinary QCD MC (hadronization procedure included), introducing there the same assumptions about $`n_f`$ and $`\alpha _s`$. The excellent agreement is illustrated by Fig. 5. The disagreement seen at large $`x`$ is natural, because the limiting spectrum is not valid there.

It is instructive to compare our MC for ordinary QCD with variable $`n_f`$ and the exact behaviour of $`\alpha _s(k_{})`$ with the limiting spectrum with fixed number of flavours $`n_f=3`$ and $`n_f=6`$. It is clear that, in neither case, $`\alpha _s(k_{})`$ from Eq. (1) describes correctly $`\alpha _s`$ in the whole interval of $`k_{}`$, and the MC spectrum should be between these two solutions. Figures. 6 and 7 show that this is indeed the case. The accuracy of each limiting spectrum compared with the MC spectra is better than 30–50%.

In Ref. , HERWIG was used to obtain fragmentation spectra in case of ordinary QCD. The maximal mass $`M_X`$ possible to simulate was $`M_X=10^{11}`$ GeV and even for this not very large value of $`M_X`$ the computations required several months. The spectra were displayed only for large $`x>0.01`$, beyond the Gaussian peak. One of the conclusions of this work was that at $`x0.2`$ the proton yield is higher than the photon and neutrino yield. However, it was later realized that this result is caused by the tendency of HERWIG to overproduce protons at large $`x`$ (Ref. , see also ).

Let us come over to our SUSY QCD MC and compare the simulated spectra with the SUSY limiting spectrum . The spectra disagree both in the position of the Gaussian peak and in its height. To clarify which assumptions of the SUSY QCD limiting spectrum are responsible for this disagreement, we re-run the SUSY QCD MC simulation with a set of assumptions as similar as possible to those used in the derivation of the SUSY QCD limiting spectrum. We found that the main reason for the disagreement is the universal dependence of $`\alpha _s(t)`$, taken as $`\alpha _s^1(t)=(b_s/4\pi )\mathrm{ln}(t/\mathrm{\Lambda }^2)`$, with $`b_s=3`$ for SUSY, together with $`\mathrm{\Lambda }=Q_0=250`$ MeV. It differs from $`\alpha _s`$ with a variable number of flavours, which is used in SUSY QCD MC, by a factor 1.4–3 in the whole $`k_{}^2`$ range, with largest disagreement at small $`k_{}`$. Changing the evolution of $`\alpha _s(k_{})`$, an agreement can be reached between the MC and the SUSY QCD limiting spectrum: we run the SUSY QCD MC including only gluons and (massless) gluinos with fixed $`b_s=3`$ and with frozen $`\alpha _s(\stackrel{~}{t})`$ for $`\stackrel{~}{t}<0.9`$ GeV2, which is a reasonable physical assumption. The comparison with the SUSY QCD limiting spectrum for partons is shown in Fig. 16. The two spectra agree indeed quite well.

An interesting alternative approach to computing the fragmentation spectra produced by decays of superheavy particles was suggested in a recent work . In this method, the event generator SPYTHIA  was used to simulate fragmentation spectra of partons and spartons into protons, photons, and neutrinos at the scale $`M_X=10^4`$ GeV. Then the DGLAP equations were used to evolve the fragmentation functions up to the scale $`10^{12}`$$`10^{13}`$ GeV.

It is premature to compare our results, since in  preliminary results are presented, but the spectra, as displayed in  and , do not agree well with ours. In particular, the Gaussian peak is broader than in our calculations, Fig. 17. Comparing these spectra one should be aware of the differences in methods and assumptions. For example, we treat spartons as cascading particles, while in SPYTHIA spartons are taken as on-shell particles, which decay but do not cascade. On the other hand, the SPYTHIA spectrum is used only as the input, and the evolution to higher energies includes cascading. This difference will be eliminated with the C++ version of HERWIG , which will be available soon.

## 6 Summary

We have developed a new MC simulation for jet fragmentation in ordinary QCD and SUSY QCD, which is valid for initial energies up to the GUT scale. The simulation includes a perturbative part, operating at virtualities higher than the infrared cut-off $`Q_0^2=0.625`$ GeV, and a hadronization part.

The perturbative part for SUSY QCD includes squarks and gluinos as cascade particles with a universal mass $`M_{\mathrm{SUSY}}`$. The evolution of $`\alpha _s(t)`$ takes into account the correct number of active flavours and spartons at a given $`t`$. The influence of the scale $`M_{\mathrm{SUSY}}`$ on the hadron spectrum is rather weak for the studied range $`300M_{\mathrm{SUSY}}1000`$ GeV. It implies that if supersymmetric particles have different masses in the range 200 -1000 GeV, the resulting difference in the hadron spectra remains small (see Section IV and Figs. 10 and 11).

The hadronization scheme is model-independent and based on the well justified and tested assumption that the hadronization function $`f_i^h(z)`$, see Eq. (8), does not depend on $`\sqrt{s}`$. Thus the hadronization function could be calculated from LEP data. Our scheme was tested at $`\sqrt{s}=58`$ and 133 GeV against experimental data and for very large values of $`M_X`$ by a comparison with the limiting spectrum solution for ordinary QCD. For the aim of this comparison, we calculated hadron spectra using the MC for ordinary QCD in the case when the limiting spectrum is known to be correct: $`n_f=3`$ and $`\alpha _s(t)`$ given by Eq. (1). The excellent agreement between both spectra is illustrated by Fig. 5.

The spectra of nucleons and secondary particles, photons, neutrinos, as well as neutralinos, have been calculated and presented in Fig. 17. These spectra can be used for calculations of fluxes of ultra high energy cosmic rays, produced by superheavy dark matter and by topological defects.

## Acknowledgments

We acknowledge with gratitude the participation of R. Sang in the initial stage of this work. We are grateful to Yu.L. Dokshitzer, V.A. Khoze, S.S. Ostapchenko and B.R. Webber for many useful explanations and comments. MK would like to thank the Alexander von Humboldt-Stiftung for a Feodor-Lynen grant and the EC for a Marie-Curie grant. The work of MK was partially also supported by DGICYT grant PB95-1077 and by the EC under the TMR contract ERBFMRX-CT96-0090. The work by VB was performed within INTAS Project No 99-01065.