text:
# Statistical characterisation of full-disk EUV/XUV solar irradiance and correlation with solar activity

## 1 Introduction

The statistical characteristics of solar irradiance have direct implications for coronal heating and terrestrial climate change; these statistics can be addressed both through flare events and, as in the present paper, through global measures. Many authors have identified transient flaring events and quantified their frequency distributions. The earliest observations of soft X-ray bursts by Drake (drake (1974)), and of hard X-ray bursts by Datlowe et al. (datlowe (1974)), Lin et al. (lin (1984)), and Dennis (dennis (1985)), found that the number of bursts with a given peak photon flux follows an inverse power law with logarithmic slope $`a_E1.8`$. These findings gave rise to the model of the active corona as the superposition of a large number of impulsive phenomena, from the largest *flares* ($`10^{33}`$ ergs) down to *nanoflares* ($`10^{24}`$ ergs) and possibly beyond. The contribution of nanoflares to coronal heating is reviewed by Berghmans berghmans (2002). Kucera et al. (kucera (1997)) present evidence for changes in the upper limit of flare energies with the sizes of active regions, while secular variations are generally not seen. However, Bai (bai (1993)) and Bromund et al. (bromund (1995)) report a 154-day periodicity in hard X-ray burst distributions, and Wheatland (wheatland (2000)) finds a variation in the flaring rate.

In this paper we investigate a different statistic: the distribution of fluctuations in the irradiance from the full solar disk, using EUV/XUV observations by the SOHO CELIAS/SEM instrument (described in Sect. 2). This global data encompasses both flaring and other activity such as waves, slow flows, and CMEs. Importantly, this also includes flare activity that could not be identified by event selection from images. By characterising the distributions with single power-law tails, we find that the exponent thereof differs by a factor of two between periods of high and low solar activity.

## 2 SOHO/SEM Data

The Charge, Element, and Isotope Analysis System/Solar Extreme Ultraviolet Monitor (CELIAS/SEM) transmission grating spectrometer, on board the Solar and Heliospheric Observatory (SOHO) spacecraft, has measured with 15 s resolution the full disk absolute photon flux at 1 AU since the launch of SOHO in 1995 Dec. Full details of the instrumentation and calibration are given by Hove et al. (hove (1995)) and Judge et al. (judge1 (1998)), and the data are available via the internet111http://www.usc.edu/dept/space\_science/sem\_first.htm. Previous EUV/XUV measurements suffered from instrument degradation, but the SEM data is highly stable and reliable with an estimated absolute uncertainty of $`\pm 10\%`$ ($`1\sigma `$) (Judge et al. judge2 (2002)). Owing to the statistical correlation of electron temperature and thermal energy in flaring processes (Aschwanden asch4 (1999)), narrowband detectors provide biased estimates of flare energy distributions (Aschwanden & Charbonneau asch1 (2002)); to overcome this problem, we choose the broadband (0.1–50 nm) channel 2 data. In Fig. 1 we plot the full-disk solar flux averaged over 15 s intervals from 1996 Jan–June (a period of low activity) and 2000 Jan–June (a period of high activity).

The six or so quasiperiodic variations with the largest amplitude are caused by the movement of active regions across the disk with the sun’s 27-day rotation period, plus contributions from the evolution of active regions over days to months (Fröhlich & Pap frohlich1 (1989)) and from the 11-year solar cycle (Willson & Hudson willson (1988)).

Since we wish to study only the short-time fluctuations, it is necessary to remove the long-term variations from the time-series. To achieve this we use a technique known as detrending, and the results are shown in Fig. 2.

Detrending is effected by subtracting from each data point the mean of the measurements over a certain time $`\tau `$ centred on that point. $`\tau `$ is chosen to be much less than the timescales to be removed and much greater than the timescales of interest. Here we set $`\tau =40`$ hours (10000 15 s sampling intervals) so as to remove the dominant 27-day cycle while retaining variations on timescales of up to a few hours.

## 3 Full-disk irradiance distributions

In Fig. 3 we plot the number distributions $`N(x)`$ of the detrended intensities $`x`$ during periods of low and high activity (shown as time-series in Fig. 2), using log-linear axes to magnify the small probabilities in the tails of the distributions.

It proved impossible to fit single distribution functions to the entire range of detrended data in this figure. Gaussian distributions, which are parabolic on these axes, clearly do not apply; log-normal distributions are inappropriate since many detrended values are negative. Fréchet distributions (Chapman et al. chapman (2002); Greenhough et al. greenh (2002)) may be fitted but they fit only $`30\%`$ of the bins within the errors shown in Fig. 3. Figure 4, which plots on log-log axes the ranges to the right of the dashed lines in Fig. 3, shows that the tails of the distributions approximate well to power laws. Single power laws of the form $`N(x)x^{a_T}`$, where $`a_T`$ is the tail exponent, can be fitted to give (i) $`a_T=3.0\pm 0.2`$ for 1996 Jan–June, and (ii) $`a_T=1.5\pm 0.1`$ for 2000 Jan–June, both having an $`R^2`$ goodness-of-fit of 0.99.

To obtain the correct value for $`a_T`$ from the slope $`m`$ of $`\mathrm{log}N(x)`$ against $`\mathrm{log}x`$ (since $`a_T=m`$), the data are placed into bins whose widths increase exponentially with $`x`$. One may also fit composite functions such as broken power laws. If the fitted slopes are allowed to change at, for example, $`x=10^9`$ (upper plot) and $`x=5\times 10^9`$ (lower plot), the best fits are respectively $`a_T=2.8\pm 0.23.3\pm 1.2`$ (overall $`R^2=0.95`$) and $`a_T=1.3\pm 0.21.5\pm 0.2`$ (overall $`R^2=0.97`$). Importantly, regardless of the details of the chosen fits, we find a significant difference in the tails of the detrended full-disk solar irradiance observations between the selected six-month periods of low and high activity.

In Fig. 5 we use single power laws as the simplest description of the distribution tails for intermediate times. We plot the slope $`a_T`$ (with error bars showing $`95\%`$ least-squares confidence bands) for each available month of data from 1996 Jan to 2000 Dec, against the sunspot number for that month (a proxy for solar activity)222Available at http://sidc.oma.be/html/sunspot.html..

Each value is calculated using a least-squares fit from the bin proceeding the mode up to the last non-empty bin in the log-log distribution of detrended data for that month. In three cases – 1998 Oct, 1999 Feb, and 2000 Nov – the error bars extend beyond the range of all the calculated slopes, and these points are not included in the figure. The plot is suggestive of a trend, albeit with a large variance, toward flatter power-law slopes as the sunspot number increases, which is consistent with the six-month averages presented in Fig. 4. To our knowledge, these results represent the first clear activity-related distinction in the statistics of full-disk solar irradiation.

## 4 Conclusions and discussion

By removing the low-frequency (greater than a few hours) variations in full-disk EUV/XUV solar irradiance observations, and calculating the distributions of the remaining detrended fluctuations, we obtain distributions whose tails can be characterised by single power laws. The power-law exponent shows significant variation with activity: for the low-activity period 1996 Jan–June, we find $`a_T=3.0\pm 0.2`$, while $`a_T=1.5\pm 0.1`$ for the high-activity period 2000 Jan–June. Comparing $`a_T`$ and the sunspot number for each month of data suggests that there is in general a flattening of the slope with increasing activity. These results provide evidence for activity-related changes in irradiance fluctuations of the full solar disk. This is particularly intriguing given the diverse range of phenomena that contribute to this integrated emission measure.

Changes in the value of $`a_T`$ may reflect important variations in the proportion or character of emissions from non-flaring regions. However, uncertainty is introduced by the dependence of the absolute photon flux on the relative spectral distribution, which may change during the intense flaring events (Ogawa et al. ogawa (1998)) that are more frequent during high activity. One would have to remove the XUV component from the EUV/XUV total – a process involving much uncertainty (Judge et al. judge2 (2002)) – in order to test for this effect . Furthermore, since the relationship between total quantities and their underlying components is an unsolved problem in the field of non-equilibrium statistical mechanics (discussed by Chapman et al. (chapman (2002)), and by Vekstein & Jain (jain (2003)) in the context of nanoflares), we do not attempt to make a quantitative link between $`a_T`$ and the flare statistic $`a_E`$ (see Sect. 1) in this paper.

###### Acknowledgements.

We are grateful to Bogdan Hnat and Nick Watkins for helpful suggestions. J. G. acknowledges a CASE Research Studentship from the UK Particle Physics and Astronomy Research Council in association with UKAEA. G. R. acknowledges a Leverhulme Emeritus Fellowship. This work was also supported in part by the UK DTI. Data provided by the CELIAS/SEM experiment on the SOHO spacecraft, and the Solar Influences Data analysis Center (SIDC, Royal Observatory of Belgium).


format: text

text:
# The Issue of Time in Quantum Geometrodynamics

## References


text:
# 1 Abstract

## 1 Abstract

The recent success of a proposal from some time ago to explain the spectrum of cosmic rays allows some strong conclusions to be made on the physics of supernovae: In the context of this specific proposal to explain the origin of cosmic rays, the mechanism for exploding supernovae of high mass has to be the one proposed by Bisnovatyi-Kogan more than 30 years ago, which was then based on a broader suggestion by Kardashev: A combination of the effects of rotation and magnetic fields explodes the star. Interestingly, this step then leads inevitably to some further suggestions, useful perhaps for the study of gamma ray bursts and the search of a bright standard candle in cosmology.

## 2 Introduction

Cosmic Ray physics has inspired us for many years, and has given rise to many interesting books and reviews, some of which are .

Several burning questions exist in high energy astrophysics:

* What is the physics of the supernova explosion of massive stars? The most common idea to be explored over the last few decades is that the burst of neutrinos, certainly sufficient in energy, as shown by the supernova 1987A, is doing it. However, the detailed mechanism has yet to be worked out successfully.
* What is the origin of Galactic cosmic rays? What makes it possible for cosmic rays from Galactic sources to reach energies such as $`\mathrm{3\hspace{0.17em}10}^{18}`$ eV?
* What is the mechanism of Gamma Ray Bursts, and what is it’s relation to supernovae?
* Is there a brighter standard candle possibly available for cosmology than the now so famous Supernovae type Ia?

Here, in this review, we wish to combine several recent advances made to give a tentative answer to all these questions, starting with the work on cosmic rays:

## 3 The cosmic ray spectrum

We wish to explain the entire spectrum of cosmic rays, and focus here on Galactic cosmic rays, so the energies up to about $`\mathrm{3\hspace{0.17em}10}^9`$ GeV. Some time ago we made a proposal to explain Galactic and extragalactic cosmic rays in six consecutive papers, . The basic idea was to distinguish the different sites of supernovae into the lower mass stars, that explode into the interstellar medium, and the higher mass stars that explode into their own stellar wind. The basic data to explain are as follows, :

* The overall spectrum is about approx. $`E^{2.7}`$ until the “knee”, which is a bend downwards at around $`\mathrm{3\hspace{0.17em}10}^{15}`$ eV
* approx. $`E^{3.1}`$ beyond the knee
* a slight downward dip from $`\mathrm{3\hspace{0.17em}10}^{17}`$ eV, sometimes referred to as the “second knee”
* a transition near $`\mathrm{3\hspace{0.17em}10}^{18}`$ eV, with then approx. $`E^{2.7}`$ again
* uncertainty beyond $`\mathrm{5\hspace{0.17em}10}^{19}`$ eV, either a mild cutoff (HIRES) or a continuation (AGASA)

and in more detail, considering electrons and nuclei:

* electrons $`E^{2.7}`$ as well (from radio, low $`E`$), up to about 20 GeV
* electrons $`E^{3.3}`$ observed to 3 TeV:
loss dominated, so injection $`E^{2.3}`$
* positron fraction a few percent
* abundances enriched
* anti-proton fraction about $`10^4`$
* heavy elements tend to have a slightly flatter spectrum

## 4 The arguments

Already in 1934 the original suggestion by W. Baade & F. Zwicky was that the most energetic particles probably come from the most energetic phenomonea known to us at the time, supernovae, . E. Fermi then provided an argument on how reflection at magnetic irregularities can enhance the energy of a charged particle, 1949 and 1954, . The detailed physics of the Fermi-acceleration was worked out in a simple approximation in a series of papers, starting in 1977 by I. Axford et al., G. Krymsky, , R. Bell in two papers, , and then again by R. Blandford & J. Ostriker,  with a nice extensive review by L. Drury, .

### 4.1 Supernova explosions into the ISM

Then P. Lagage & C. Cesarsky, , worked out the details for the maximum energy of a particle subject to Fermi acceleration in the shock of a supernova, exploding into the interstellar medium. Presciently, they already used the self-similar expansion into a very tenuous medium, now known to exist from X-ray observations. For the details of the acceleration of energetic particles at the shock they use the concept that most of the acceleration happens in a regime, where the magnetic field is parallel to the shock normal - as well its disturbances \- and provides most of the scattering. Here the expansion solution is a similarity or Sedov-solution, a self-similarity, where the outer shock radius $`r`$ behaves as $`rt^{2/5}`$, and the shock speed is $`\dot{r}r^{3/2}`$, and so slows down considerably over time. This slow late evolution makes these supernova remnants dominant in observations, they just live a rather long time.

* $`E_{max}\mathrm{\hspace{0.17em}10}^5ZB_6GeV`$, using already the hot ISM of density $`\mathrm{3\hspace{0.17em}10}^3\mathrm{cm}^3`$, confirmed only in 1997, . The expansion is self-similar, usually referred to as a Sedov expansion, and yields then a maximum particle energy of about $`Z\mathrm{\hspace{0.17em}10}^{14}`$ eV, an order of magnitude below the knee.
* Maybe secondary acceleration is then required in the ISM - maybe by shockwaves inside superbubbles, the motley mixture of powerful stellar winds, and young supernova explosions. This picture has mostly been explored by I. Axford. A quantitative prediction is still outstanding for this scenario.

### 4.2 The knee

The key problem, recognized immediately, is that this maximum energy does not even come to the energy at the knee, and so that feature cannot be explained in such an approach. Therefore many attempts have been made:

#### 4.2.1 Transport?

Cosmic rays get injected with a spectrum, which is subsequently modified by losses from the disk, and so obviously, to get a kink in the spectrum may mean a change in behaviour of these diffusive losses; such a kink could be visible as the “knee”. This is a question of the transport of cosmic rays through a turbulent magnetic ionized gas; general discussions of the magnetic field are in, e.g., . Such a change in transport implies a special scale in the interstellar medium, and would also lead to an increased anisotropy at the kink. However, there is no such scale as yet recognized.

#### 4.2.2 Turbulence?

The evidence from the interstellar medium and from plasma physics suggests - if described by an isotropic approximation with an eye towards transport of cosmic rays - that the spectrum is a Kolmogorov spectrum. There are three arguments in favour of using a Kolmogorov spectrum for the context of cosmic rays: First, the difference in behaviour between nuclei, and electrons, allows the exponent of the diffusive law to be limited by about 0.4, so the energy dependence of the diffusion coefficient $`\kappa `$ is weaker than $`\kappa E^{0.4}`$, . Second, any steeper dependence would produce anisotropies at the higher energies, which are not seen, . And, finally, any steeper dependence would seriously dilute the number of stars contributing at the higher energies, making any reasonable powerlaw fit to the data meaningless, despite the fact, that such fits work well to describe the data, .

#### 4.2.3 Spallation?

Cosmic rays interact, and in fact most isotopes of Lithium, Beryllium, and Bor, as well as the odd-Z elements, and the sub-iron elements arise from cosmic ray interaction, the destruction of nuclei in a spallation process, see M. Garcia-Munoz et al., . Therefore, a change in this interaction may lead to a loss of nuclei at higher energy; the interstellar medium is discussed, with some emphasis on the Galactic Center region, e.g., in . However, when working through this idea, one finds that the column for interaction must be much higher than what is observed; in fact, from the abundances of all the spallation products, as well as positrons and anti-protons we have a fairly good idea what the interaction actually is, and it is not enough to provide such a drastic change. Furthermore, that data suggest today, that the “knee” is a feature in rigidity $`E/Z`$, and so plasma physical processes are implied.

#### 4.2.4 New source?

Obviously, in some proposals we require a new source at $`\mathrm{3\hspace{0.17em}10}^{18}`$ eV, or maybe even at $`10^{14}`$ eV, and so one might realistically ask, whether a new source might not be necessary already at lower energies, if that source is supposed to work so well at higher energies. A new source is easiest to accomodate, if the new source, dominant at higher energy, has a flatter spectrum than the source dominant at lower energy; in such a case the transition is trivial, since the higher energy source just overtakes the lower energy source naturally in flux. However, if the transition is supposed to occur at or near the knee, then the higher energy source must have a steeper spectrum rather than a flatter spectrum, and so we require a multitude of strange coincidences: The higher energy source must have a peaked spectrum, peaking just when the lower energy source peters out, and so have the same flux at that energy, while not contributing at lower energy in any significant way. That is really hard to believe. So, it is then actually simpler, and has been suggested, that we use a new source from the start, and, e.g., adopt the point of view, that Gamma Ray Bursts produce all the cosmic rays, in work by C. Dermer (discussed in a talk at the Aspen meeting, January 2002). However, then we are back where we started, we still need an argument why we have a change in slope at a specific rigidity, the knee. Also, the expected flux of cosmic rays from Gamma Ray Bursts is nowhere near to what is necessary, . However, Gamma Ray Bursts should produce cosmic rays, and the question is whether there is any way to ascertain their contribution, perhaps through high energy neutrinos. Another possible source could be the activity in the Galactic Center, . However, there is no observational evidence at this time, that such a localized activity could be dominant at the location of the Solar system.

#### 4.2.5 Change in acceleration?

Therefore, what remains, is more mundane perhaps, a small modification in acceleration of the existing cosmic rays. The proposal by I. Axford et al.  is that normal supernovae provide cosmic rays to about the knee, and that then in the supernova wind bubbles around the really massive stars, the cosmic rays get further acceleration, and then have a slightly diminished efficiency in acceleration; after all, the detailed physical properties are somewhat different in the environment of many supernovae and their predecessor winds, as compared to the interstellar medium. The only step missing in this very nice argument is a quantitative prediction as to what expect in this picture; one may surmise that the predictions may not turn out to be quite so different from the proposal described below.

## 5 Proposal for the Origin of Galactic Cosmic Rays

The riddle, which provoked a new approach now some time ago, was the key observation, that the radio shells of supernova remnants show a radial direction for the magnetic field, parallel to the shock-normal. This is in exact contradiction to the simple plasma physics expectation that the enhancement of the magnetic field component parallel to the shock surface should make that component dominant behind the shock, see .

Additional support for a new approach came from the observation that radio features could be seen in motion, with an amplitude of velocity quite close to the overall expansion speed of the supernova shock front, and yet measured with respect to the expanding frame of reference.

Third, the thickness of the observed shells both in the radio and also to some degrees in X-rays, is much too large to allow a simple explanation. There must be some reason for these shells to be almost 1/3 of the outer radius oberved.

Finally, there appeared a theoretical argument, that cosmic ray dominated shocks are unstable, in work by R. Ratkiewicz, I. Axford, G. Webb, e.g., .

Therefore, we introduced the key premise that the particle transport in the shock region is dominated by large scale turbulence; in such a case there is no time for any turbulent cascade, and so no time for a Kolmogorov cascade either. Therefore the transport is just given by a specific scale, and a specific velocity, naturally the thickness of the shocked region, and the velocity difference across the shock. This is in the spirit of the early work on turbulence such as by L. Prandtl and Th.v. Karman. This premise allows then to explain all these observed features.

This means that the shock region is strongly turbulent, and the front of the shock, as seen along the line of sight at the edge of some observed supernova remnannt, is actually a superposition of many smaller shocks, advancing and retreating, like the waves on a beach. The geometric scale of these shocks is again given by the radial overall shock thickness, just given by momentum conservation. Therefore, curvature drifts do play a role, and introduce a scale in energy per charge $`E/Z`$.

It can be surmised, that the limit cycle arguments by M. Malkov & L. O’C. Drury will eventually lead to a similar conclusion, . Such a convergence of theoretical concepts will be an important test of the picture.

### 5.1 Supernovae into a stellar wind

The concept then leads to the consideration of supernovae exploding into their own stellar wind, . Stellar winds may be considered in a first approximation as a Parker wind, with $`B_\varphi \mathrm{sin}\theta /r`$ in spherical coordinates. The dominant magnetic field character is an Archimedian spiral, with the field mostly very nearly tangential. The key aspects are i) that the wind is magnetic, with the magnetic field given by the star, and not the interstellar medium, and ii) that the wind has a density gradient $`r^2`$, which means that the shock stays fast, and does not slow down as in an explosion into a medium of constant density. Also, magnetic fields in massive stars exist, they have been detected through nonthermal radio emission and maser emission, .

As the magnetic fields get very weak and also radial towards the pole of a Parker-type wind, one needs to consider shock acceleration in two limiting regimes:

At the pole and near to it we need to use the acceleration dominated by the magnetic field parallel to the shock normal, as done by . However, at the equator, where the shock direction is perpendicular to the dominant magnetic field, we need to use another limit, as shown by R. Jokipii, . Therefore we also need a matching condition, and this leads to a proposal for identifying the knee as arising with this matching energy, the maximum energy near the pole and close to it:

$$E_{max}=ZerB(3/4v_{sh}/c)^2$$
(1)
The maximum energy in the equatorial region, and in fact, most of $`4\pi `$ is

$$E_{max}=ZerB$$
(2)
Since $`B`$ is inversely proportional to radius $`r`$ in the wind, these energies are actually constant throughout the wind, see e.g. . We still need to show, that the particles actually have enough time to reach these maximum energies, but this appears plausible. We also need to show, why this leads to a different spectrum beyond the knee: This is due to the effect that curvature drifts, which provide a small proportion of the particle’s energy gain, when subject to many shock transitions, get slightly weaker when the Larmor radius of the particle moves past the scale corresponding to the maximum energy (really $`E/Z`$) near the pole, see  for a very detailed discussion.

Those stars which have powerful winds, are usually enriched in chemical elements such as Helium, Carbon and Oxygen, and so this line of reasoning actually picks up a theme, that Wolf Rayet stars are a major contributor to cosmic rays, . Therefore, supernova explosions into stellar winds are proposed to provide most of the Galactic cosmic rays, with the notable exception of the element Hydrogen, for which they provide only some fraction.

## 6 Stellar evolution

For the physics of cosmic rays and supernovae we need to consider stellar evolution in some detail. What happens to stars of different zero-age main sequence masses $`M`$:

* Single stars with $`M<\mathrm{\hspace{0.33em}8}M_{}`$ give no supernovae.
* Stars in the mass range $`8<M<\mathrm{\hspace{0.33em}15}M_{}`$ explode as supernovae, but the explosion goes into the interstellar medium. Such explosions lead a Sedov type expansion, which becomes slow and has a long-lasting remnant.
* Stars with $`15<M<\mathrm{\hspace{0.33em}30}M_{}`$ explode as supernovae with a substantial wind, but enriched only in Helium. The mass in the wind-shell is moderate. We refer to these stars as Red Super Giant stars (RSG).
* Stars $`>\mathrm{\hspace{0.33em}30}M_{}`$ explode as supernovae into a strong wind; the chemical abundances are strongly enriched, since the mass loss eats deeply down back into the star. The mass in the wind-shell is large. We refer to these stars as Wolf-Rayet stars (WR).

Obviously, these numbers for the masses are approximations.

## 7 Predictions

The predictions, all made some time ago, here focussing on wind-SNe, were described in .

They gave the following quantitative spectra:

* $`E^{2.74\pm 0.04}`$ for the interstellar medium supernovae, with the abundances that correspond to the medium through which the supernova shock races. This spectrum is predicted to have a cutoff in the 100 TeV range.
* $`E^{2.670.02\pm 0.02}`$ below knee for the wind supernovae, for which the shock races through the predecessor wind, and correspondingly has enriched abundances.
* $`E^{3.070.07\pm 0.07}`$ above knee for those same stars winds.
* A bending of the wind spectra at the “knee” near $`Z\mathrm{\hspace{0.17em}10}^{15}`$ eV, and a final cutoff at the “ankle” near $`Z\mathrm{\hspace{0.17em}3\hspace{0.17em}10}^{17}`$ eV. These specific rigidities have an uncertainty of about a factor of 2.
* When the supernova explodes, a powerful shock wave races through the wind, and then smashes into the shell; this wind-shell is made up of both wind-material and interstellar medium material, from the environment of the predecessor star. Such a shell has two shocks at its boundaries, one slowing down the wind, and one speeding up the interstellar material. If there is no convective instability, then a contact discontinuity separates wind material from environmental interstellar medium material.
* The shock is loaded with energetic particles from cosmic ray acceleration, and then these cosmic rays suffer from spallation in the shell, . Both cosmic ray particles as well as wind material particles get broken up. This, in the case of WR shells, gives the Be, B, Li nuclei in cosmic rays, which serve then as tracers of cosmic ray lifetimes, and interaction. It also gives the odd-$`Z`$ elements, as well as the sub-Fe elements in cosmic rays.
* As a new prediction following from the work in  we find a small spectral flattening is expected due to differential spallation for Fe-like nuclei. The fraction of nuclei removed by spallation from the set of Fe-like nuclei is so large, that across the range of energy per particle measured, this fraction is larger for lower energy, and so the spectrum flattened. This proposal is quantitatively consistent with the data.

At that stage we did not differentiate the cosmic ray acceleration scenario for red supergiant winds (RSG), and Wolf Rayet (WR) winds. We now do differentiate those two kinds of stars with winds in their interaction in the wind shell, . The RSG stars have a weaker wind, therefore a lower mass shell as compared to the WR stars, and so we argued that the cosmic ray interaction, when the supernova shock finally hits the shell, is convective for RSG stars, and diffusive for WR stars. In this picture the RSG stars produce the gamma ray emission from the Galactic plane, and the WR stars produce the spallation products such as Beryllium, and give only some small contribution to the Galactic gamma ray emission at moderate photon energies.

### 7.1 Implications

Such a model implies, as had been predicted already by B. Peters in the fifties, , that the turnover at the knee, and at the final cutoff is gradual, as the elements roll off in sequence of $`Z`$. Furthermore, since the knee and final cutoff energy are given by stellar parameters, the magnetic field, implicitly the rotation (the wind must be asymptotic already at the surface, which implies slow rotation at the surface) and the explosion energy, these parameters should be very nearly the same for all stars.

This corresponds exactly to what has been argued in a very different context already by G. Bisnovatyi-Kogan in (1970), who based his argument in turn on a broader suggestion by N. Kardashev (1964), . In this picture, the predecessor star is rotating and has a magnetic field. When the core of the star collapses, the core is spun up to a fast rotating disk, and the collapse along the plane of symmetry stops. The potential energy at that stage then is transmitted to the envelope by the torque of the magnetic fields (i.e. their angular momentum transport), and that potential energy constitutes the explosion energy.

Therefore, if the data are really well fit by the model, the conclusion is strongly suggested that not only the concept of magneto-rotational explosion works, but that in fact it works with very similar numbers in all stars.

## 8 New data

Since the original detailed quantitative predictions were made, new data have appeared, which provide very serious tests, e.g. :

* At present the data set from the KASCADE experiment is the best across the knee region, and provides a first challenging test. As we show below they are very well fit by the data. Other data also show consistency, but with large error bars, and with unknown systematics still, so further conclusions can not be drawn from them.
* In paper  we predicted the energy dependence of the B/C ratio, from interaction in the thick WR-star shells, as $`E^{5/9}`$. A refined treatment was given in . As V. Ptuskin then showed some time later, at the ICRC in Salt Lake City, , the best fit to the data gives an energy dependence of $`E^{0.54}`$, which is quite consistent.
* In 1997 St. Hunter et al., , and M. Mori, , showed that the gamma ray emission from the Galactic plane could only be fitted by a cosmic ray spectrum, which was substantially flatter than the one observed. In  we suggested that the spectrum could be explained as interaction in the source region of cosmic rays, specifically the RSG star shells. The interaction happens then with a spectrum of $`E^{7/3}`$. We will discuss the details of this fit elsewhere.

### 8.1 Specific tests

The tasks ahead of us are the following, to be quite specific:

* We normalize all abundances at TeV, using the data collection of .
* We need to remember, that Hydrogen has two components, only one of which may have a knee (wind-SNe, RSG), the other one has a cutoff below the knee (ISM-SNe); on the other hand Hydrogen is extremely low in WR winds. We need to note here, that the data do not force us at this time to assume that the knee of the RSG stars and the WR stars is the same; in fact, it could be, that both knee and cutoff for RSG stars are quite different, much lower in rigidity than for WR stars. In our first data fits, described below, we do not yet distinguish RSG and WR stars as regards the knee and final cutoff. The KASCADE data are forceful only for the abundances from the WR stars.
* Also Helium may have two different observable knees or cutoffs, corresponding again to RSG supergiants, and WR supergiants, but they also might be same; we have to keep this uncertainty in mind.
* We then extrapolate with the predicted spectrum to the knee in rigidity, bend at the knee in rigidity, and extrapolate to the cutoff. As the cutoff is also in rigidity, there is a rolloff of the overall spectrum, starting from the cutoff in Helium, which we can actually identify with the “second knee”.
* We then compare to KASCADE, Akeno, HIRES, … One important point is to remember, that any fit to the data across the knee has to “reach” the data at higher energy, such as from Akeno, etc. This restricts any possible fits quite severely, see .

### 8.2 Work to be done

The work is spread among many of our partners, as follows:

* The detailed abundances are being considered with A. Popescu (Bukarest), and N. Langer (Utrecht).
* The energy in explosions is considered with G. Pavalas (Bukarest).
* The positron spectra are done with W. Rhode et al.(Wuppertal).
* The gamma ray spectrum of the Galaxy is done with S. Casanova (Bonn), R. Engel (Karlsruhe), W. Rhode (Wuppertal), and many others.
* The anti-proton spectra are done with E.-S. Seo, R. Sina (Univ. Maryland), and R. Engel (Karlsruhe).
* The detailed fit to the KASCADE data is first described in this paper, in the following main section.
* Later we will also consider a more detailed fit to the Akeno, AGASA and HIRES near the ankle with S. Ter-Antonyan (coauthor here), A. Vasile (coauthor here), Zh. Cao (Utah) and St. Westerhoff (New York).

## 9 Cosmic ray properties from air shower data

The following section was written by S. Ter-Antonyan.

In general, the relation between energy spectra ($`\mathrm{}_A/E_0`$) of primary nuclei ($`A`$) and measured EAS size spectra at observation level ($`\mathrm{\Delta }I/\mathrm{\Delta }N_{e,\mu }^{}`$) is determined by an integral equation

$$\frac{\mathrm{\Delta }I(\theta )}{\mathrm{\Delta }N_{e,\mu }^{}}=\underset{A}{}_{E_{min}}^{\mathrm{}}\frac{\mathrm{}_A}{E_0}\frac{W(E_0,A,\theta )}{N_{e,\mu }^{}}𝑑E_0$$
(3)
where the kernel function ($`W/N_{e,\mu }^{}`$) of the equation depends on primary energy, type of primary nucleus ($`A=1,4,\mathrm{},59`$), zenith angle ($`\theta `$), $`AA_{Air}`$ interaction model and response functions of measurements.
Eq. (3) is a typical ill-posed problem and has an infinite set of solutions for unknown primary energy spectra. However, the integral Eq. (3) turns to a Fredholm equation, if the type of a primary nucleus is defined directly in the experiment (as it is in the balloon and satellite measurements  where energy spectra for different nuclei are obtained up to $`10^{15}`$ eV).
Similarly to  we reconstruct primary energy spectra based on a transformation of the integral Eq. (3) to a parametric equation with unknown spectral parameters. For that, instead of unknown primary energy spectra of the integral Eq. (3) we consider energy spectra according to the predictions of a multi-component model . In this case, unknowns are such model parameters as cut-off and knee energies, spectral indices, scale factors and fractions of different components. EAS muon and electron size spectra on the left side of the Eq. (3) are taken from recent KASCADE publications .
The kernel function of Eq. (3) $`W/N_{e,\mu }`$ is preliminarily calculated by a Monte-Carlo method provided a given ($`AA_{Air}`$) interaction model of primary nuclei and atmosphere.
Solutions of Eq. (3) for unknown spectral parameters at *a priori* known primary spectra are easy to obtain by means of $`\chi ^2`$-minimization method . It is clear that solutions found that way are partial solutions of the integral Eq. (3). However, if the values of the spectral parameters are consistent with the predictions of a primary spectrum model calculations , it will confirm the use of this model in the investigation energy region.
We have already done such investigations  by testing a multi-component model of primary cosmic ray origin with observed EAS size spectra from KASCADE (1020 g/cm<sup>2</sup>) and ANI (700 g/cm<sup>2</sup>) using experiments in 5 different zenith angular intervals.
Here we consider the inverse problem based on KASCADE EAS electron and ”truncated” muon size spectra  in 3 zenith angular intervals. The combined analysis of electron and muon EAS size spectra allows us to determine the influence of a nuclear interaction model on the quality of the inverse problem solution more accurately.

### 9.1 Parametrization of the primary energy spectra

The energy spectra of primary nuclei according to the multi-component model of primary cosmic ray origin  at energies $`10^{12}10^{18}`$ eV are presented here in a 3-component form:

$$\frac{\mathrm{}}{E_A}=\beta \mathrm{\Phi }_A\left(\delta _{A,1}\frac{d\mathrm{}_1}{dE_A}+\delta _{A,2}\frac{d\mathrm{}_2}{dE_A}\right)+\mathrm{\Phi }_A^{EG}\frac{d\mathrm{}_3}{dE_A}$$
(4)
where the $`\beta `$ is a dimensionless normalization parameter, $`\mathrm{\Phi }_A`$ are scale spectral factors and model parameters $`\delta _{A,i=1,2}`$ are the fractions of each component ($`\delta _{A,1}+\delta _{A,2}=1`$) in a primary flux of nuclei (A).
The first component (ISM) is derived from explosions of a normal supernova into an interstellar medium with expected rigidity-dependent power law spectra

$$\frac{d\mathrm{}_1}{dE_A}=\{\begin{array}{cc}E_A^{\gamma _1}:\hfill & E_A<E_{cut}Z\hfill \\ 0:\hfill & E_A>E_{cut}Z\hfill \end{array}$$
(5)
where the model parameter $`E_{cut}Z`$ is a cut-off energy of ISM component at $`Z`$ nuclear charge.
The second component (SW) is a result of explosions of stars into their former stellar winds with expected rigidity-dependent power law spectra

$$\frac{d\mathrm{}_2}{dE_A}=\{\begin{array}{cc}E_A^{\gamma _2}:\hfill & E_A<E_kZ\hfill \\ E_k^{\gamma _2}(E_A/E_k)^{\gamma _3}:\hfill & E_A>E_kZ\hfill \\ 0:\hfill & E_A>E_{cut}^{SW}Z\hfill \end{array}$$
(6)
where the model parameter $`E_kZ`$ is a knee energy of SW component and $`E_{cut}^{SW}2.210^5`$ TeV is a corresponding cut-off energy .
The third, extragalactic (EG) component is approximated also by rigidity-dependent power law spectra

$$\frac{d\mathrm{}_3}{dE_A}=\{\begin{array}{cc}E_{ank}^{2.75}(E_A/E_{ank})^2:\hfill & E_A<E_{ank}Z\hfill \\ E_A^{2.75}:\hfill & E_A>E_{ank}Z\hfill \end{array}$$
(7)
The values of model predictions  for spectral parameters are:
$`\gamma _1=2.75\pm 0.04`$, $`\gamma _2=2.67\pm 0.03`$, $`\gamma _3=3.07\pm 0.1`$, $`E_{cut}120`$ TV, $`E_k700`$ TV, $`\mathrm{\Phi }_{A=1}^{EG}0.032`$ $`(m^2secsterTeV)^1`$, $`E_{ank}6.510^5`$ TV at factors of uncertainty $`2`$.
In  we have already obtained the evaluations of spectral indices of ISM and SW components: $`\gamma _1=2.78\pm 0.03`$, $`\gamma _2=2.65\pm 0.03`$ as solutions of parametric Eq. (3) using KASCADE  and  EAS size spectra at 5 zenith angular intervals. This is in quite good agreement with a similar early test . Since these values agreed with the model predictions here we set them fixed.

### 9.2 EAS size spectra

In general, the kernel function of Eq. (3) in EAS inverse problems is determined by

$$\frac{W}{N_{e,\mu }^{}}_{\theta _1}^{\theta _2}_0^{\mathrm{}}\frac{G(E_0,A,\theta )}{N_{e,\mu }}\frac{P(N_{e,\mu })}{N_{e,\mu }^{}}\frac{\mathrm{sin}\theta }{\mathrm{\Delta }_\theta }𝑑\theta 𝑑N_{e,\mu }$$
(8)
where $`G/N_{e,\mu }`$ is an expected EAS size electron (truncated muon) spectrum at a given observation level for $`E_0,A,\theta `$ parameters of a primary nucleus and depends on $`AA_{Air}`$ interaction model, $`\mathrm{\Delta }_\theta =\mathrm{cos}\theta _1\mathrm{cos}\theta _2`$, and $`P/N_{e,\mu }^{}`$ is the error function of measurements .
We calculated the shower spectra $`G/N_{e,\mu }`$ at KASCADE observation level (1020 g/cm<sup>2</sup>) using the CORSIKA6016(NKG) EAS simulation code  with the QGSJET01  and SIBYLL2.1  interaction models. Input parameters in simulations were: primary energies $`E_A32,100,\mathrm{}3.210^5`$ TeV, 4 groups of primary nuclei $`A1,4,16,56`$ and 3 zenith angular intervals similarly to KASCADE experimental data . Intermediate values are calculated using log-linear interpolations. Relative statistical errors of calculated shower spectra were less than $`3\%`$.
The expected detectable EAS size spectra $`W/N_{e,\mu }^{}`$ were calculated using log-Gaussian approximations for $`G/N_{e,\mu }`$ and error functions $`P/N_{e,\mu }^{}`$ which were obtained with high accuracy in the investigated energy range.
We also investigated the error function of measurements using CORSIKA code and for reconstructions of truncated muon sizes we obtained accuracy approximation: $`\mathrm{\Delta }N_\mu /N_\mu 15/\sqrt{N}_\mu +0.05`$ that exceeds the corresponding values of KASCADE approximation . In our calculations we used a more precise standard option of CORSIKA6016.
For the left side of Eq. (3) we used KASCADE EAS electron and truncated muon size spectra from  (Fig. 2, symbols). In order to make a comparison and subsequent normalization of our data with direct balloon and satellite measurements at about $`10^{14}`$ energy range we also considered the data from early publications of KASCADE  with corresponding zenith angular corrections (the first 7 symbols for each angular interval in Fig. 2a).

### 9.3 Results

Minimization of $`\chi ^2(𝐈,𝐏)`$-functional with a measurement vector

$$𝐈\{\mathrm{\Delta }I/\mathrm{\Delta }N_{e,i,k}^{},\mathrm{\Delta }I/\mathrm{\Delta }N_{\mu ,j,k}^{}\}$$
(9)
and a corresponding prediction vector $`𝐏`$ from the right-hand part of Eq. (3) was carried out at $`i=1,\mathrm{},42`$ energy intervals of EAS electron size spectra, $`j=1,\mathrm{},26`$ energy intervals of EAS truncated muon size spectra and $`k=1,2,3`$ zenith angular bins from KASCADE data . However, the combined analysis of electron and muon size spectra at $`\chi ^2`$-minimization requires to include in the expected shower spectra 2 additional unknown dimensionless parameters $`\eta _e`$ and $`\eta _\mu `$ which define a constant bias of each spectrum due to peculiarities of interaction model and systematic measurement errors .
In the first instance, inner spectral parameters ($`\gamma _{1,2,3}`$, $`E_{cut}`$, $`E_k`$, $`\delta _A`$, $`\mathrm{\Phi }_A`$) of primary spectra (4-6) have been evaluated using normalized dimensionless EAS size spectra in the measurement vector

$$𝐈_𝐧\{\frac{1}{U_e}\frac{\mathrm{\Delta }I}{\mathrm{\Delta }N_{e,i,k}^{}},\frac{1}{U_\mu }\frac{\mathrm{\Delta }I}{\mathrm{\Delta }N_{m,j,k}^{}}\}$$
(10)
where the scalar parameters $`U_e=_{i,k}(\mathrm{\Delta }I/\mathrm{\Delta }N_{e,i,k})`$ and $`U_\mu =_{j,k}(\mathrm{\Delta }I/\mathrm{\Delta }N_{\mu ,j,k})`$. Corresponding normalizations of spectra were performed for all components of the prediction vector $`𝐏_𝐧`$ at the minimization of the $`\chi ^2(𝐈_𝐧,𝐏_𝐧)`$-functional. After evaluation of the inner spectral parameters the values of $`\eta _e`$ and $`\eta _\mu `$ were determined using the $`\chi ^2(𝐈,𝐏)`$ functional by definition (8) at fixed inner spectral parameters.
Finally, the dimensionless parameter $`\beta `$ in expression (4) was determined by a normalization of the all-particle spectrum obtained $`(\mathrm{}/E_A)`$ with JACEE data at 100 TeV energy.

Table 1 contains the values of all spectral parameters which were obtained by the method above at QGSJET01  and SIBYLL2.1  interaction models and KASCADE EAS data . The normalization factors in Table 1 are: $`\beta =1.13\pm 0.05`$ at QGSJET model and $`\beta =1.0\pm 0.05`$ SIBYLL model.
The energy spectra obtained of different nuclei and the corresponding all-particle spectrum in comparison with JACEE  and RUNJOB  direct measurements and EAS data  up to energy $`310^{18}`$ eV are presented in Fig. 1.

The thick line corresponds to the expected energy spectra of the extragalactic component . This flux is also included in the all-particle spectrum in Fig. 1
The expected EAS electron and truncated muon size spectra corresponding to our primary spectra obtained with the QGSJET interaction model in comparison with the KASCADE data  are shown in Fig. 2a,b.

The value of $`\chi ^21`$ with 2$`\%`$ uncertainty (the width of a shaded area) of expected data (and 4$`\%`$ uncertainty at SIBYLL model). For biases of shower spectra (Fig. 2a,b) we obtained the values: $`\eta _e=1.23\pm 0.04`$, $`\eta _\mu =1.12`$ at QGSJET model and $`\eta _e=1.0\pm 0.01`$, $`\eta _\mu =1.3\pm 0.02`$ at SIBYLL model. These values point out the possible existence of systematic biases both in the EAS measurements and in the interaction models.
In Fig. 3 the expected muon lateral distribution functions (symbols) and corresponding KASCADE data shaded area) from  at different intervals of detected truncated muon size are shown.

In order to test the behavior of primary energy spectra obtained in a larger energy range we also calculated the expected hadron energy spectra at mountain level for comparison with data of the PION experiment (3200m a.s.l.)  at 1-7 TeV hadron energies at observation level (700 g/cm<sup>2</sup>). The calculation was carried out using the QGSJET and SIBYLL interaction models and using the method introduced above without normalizations ($`\beta =1`$). The best agreement of the expected hadron flux at mountain level in the framework of the 3-component primary model predictions and PION data was achieved at $`\delta _{A=1,2}=0.35\pm 0.05`$ and $`\delta _{A>1,2}=0.9\pm 0.1`$ fractions of component. These results are shown in Fig. 4. The effective primary energy range corresponding to these data is $`E_010100`$ TeV/nucleon.

Thus, predictions of the multi-component model of the cosmic ray origin  explain the measured KASCADE EAS data in the knee region ($`E_010^{14}10^{17}`$ eV) and hadron spectra at mountain level ($`E_010^{13}10^{14}`$ eV) with an accuracy of 10-15$`\%`$ in the framework of the QGSJET and SIBYLL interaction models respectively.
The rigidity-dependent behavior of spectra for different primary nuclei is the same for the two interaction models.
The agreement of the expected all particle spectrum and world data in $`E_010^{17}10^{18}`$ eV primary energy range displays the presence of the extragalactic component of primary cosmic rays in accordance with the 3-component model predictions.

## 10 The fit and consequences

Here we have described the actual fit to the data; we did this in two steps, first just a direct fit to the KASCADE elemental data as published, as done by A. Vasile, and then by an extensive Monte Carlo run to consider all published detailed distributions from KASCADE directly, described in the section above written by S. Ter-Antonyan. We refer here to the ICRC publications, such as . An independent, but similar approach is described by G. Schatz, ; recent KASCADE papers are . It is gratifying, that a first direct fit, as done by A. Vasile, and the extensive fitting done by S. Ter-Antonyan, not only agree with each other, but actually agree with the prediction, and with the early tests, .

We summarize, that the KASCADE data as well as a first approximation to the HiRes, Akeno and Yakutsk data can be fully fitted by the model with parameters very close to what had been predicted. We need to emphasize here, that we checked the consistency of the predicted model with the data; we did not check the uniqueness of the solution found. We do note that no other model has been either inferred nor proposed, which fits all the same data, but it is quite conceivable that there is such a model, yet to be worked out or published. The approaches by M. Malkov, J.-P. Meyer and L. O’C. Drury, E. Berezho et al., and I. Axford et al. all provide a different approach to the quest for the origin of cosmic rays. The various approaches may yet converge to a common theory.

Focussing now on the concept originally started some time ago, there are a number of important conseqences from the success of the model proposed:

### 10.1 A common final state for massive stars

The are several very interesting consequences from the success of this approach. Some of these consequences follow by necessity, which entails if these consequences can be falsified, then the original model fails. We will identify these critical steps below.

* The abundances in stellar winds near to and beyond 20 solar masses, RSG and WR supergiants dominate the cosmic ray chemical composition – this is in contrast to the superbubble concept; however, since the superbubbles in fact also have a mixture of wind abundances, these two approaches may be consistent in this point with each other. The argument, that the abundances of cosmic rays derive from accelerated dust particles is starkly different, .
* The energy in the cosmic ray population of Helium, Carbon, and Oxygen is so large as to require a much larger energy content in cosmic rays per supernova; note that the RSG and WR stars are the dominant contributors, stars which are quite rare. The energy requirement in cosmic rays per supernova of this kind is then about $`10^{51}`$ ergs.
* If we allow for an overall efficiency of 10 % of putting explosion energy into cosmic rays, the explosion energy is then required to be $`10^{52}`$ ergs for such supernovae.
* The most important conclusion from the approach outlined here is that stars near to or above 25 solar masses explode due to the magneto-rotational mechanism. We cannot say for sure at this stage at which exact mass we have a transition in supernova mechanisms to the magneto-rotational mechanism, since the cosmic ray argument is really clear only for the WR stars.
* Knee energy and cutoff energy need to be the same for all supernovae contributing. This is easiest if those stars form a “common final state”, as regards rotation, magnetic fields, mass, explosive energy.
* There is first evidence for an explosion energy of order $`10^{52}`$ ergs from a very massive star from SN1998bw, . That specific supernova was also highly asymmetric, with a rotational symmetry quite possible, consistent with the mechanism proposed by Bisnovatyi-Kogan.

### 10.2 The magneto-rotational mechanism

The magnetorotational supernova mechanism was suggested by Bisnovatyi-Kogan, . The idea was to get the energy for the explosion from the rotational (gravitational) energy using the magnetic field. Most stars are rotating differentially. The magnetic field is “frozen” into the ionized gas of the star. Then differential rotation leads to an amplification of the toroidal component of the magnetic field. When magnetic pressure becomes comparable with the gas pressure a compression wave forms and moves out to the envelope of the star. When moving along a steeply decreasing density profile this wave transforms to a strong fast MHD shock which pushes part of the envelope of the star to the outside, and so forms an explosion. Simulations of the magnetorotational mechanism for a magnetized cloud  show that this mechanism is rather effective. The simulations of the collapse of a rotating pre-supernova star show that after the short collapse stage the star consists of a rapidly rotating core and a slowly rotating envelope. At the transition region between the core and the envelope the rotation is very differential and in such a situation the inclusion of even a weak initial magnetic field could lead to amplification of the toroidal component and so to a magnetorotational explosion. Further work is in .

### 10.3 Implications for Gamma Ray Bursts?

There is an interesting resemblance of this explosion mechanism to the concept of a hypernova developed in the context of Gamma Ray Bursts, . Using the concept, , does indeed suggest an energy scale of also $`10^{52}`$ erg, derived in that specific model as an upper energy limit implied by a fit to the fluence distribution, allowing for observational selection effects. The recent observations by Schaefer et al., , support the concept that at least some GRBs blow up into a stellar wind.

This leads then one more time to the question of what the relation is between Gamma Ray Burst explosions and Supernova explosions, a question explored by many.

There is much recent work also by C. Wheeler, et al., .

### 10.4 A bright standard candle for cosmology?

There is a corollary also for cosmology: If the explosion energy is the same for all supernovae above about a zero age main sequence mass of near 25 solar masses, as suggested by the arguments above, one may speculate even further: Could it be that also the maximum luminosity, or the integral over the luminosity time curve is the same for all such supernovae? There is the obvious problem, that these explosions are rotationally symmetric only, and so present vastly different views to observers at different angles to the symmetry axis. This may be resolved by infrared observations, and also by polarization observations, just as for active galactic nuclei, for which infrared observations give the best approximation for an isotropic emission, and polarization observation do give clues for the angle between symmetry axis and line of sight.

If we could thus derive from observations the maximum luminosity or its time integral, we might have a much brighter standard candle for use in cosmology than the supernovae of type Ia, which derive from the collapse of a white dwarf in a stellar binary system, and have a fairly low explosion energy.

### 10.5 Transition: Beyond GZK cutoff

Let us finish with a brief note on the events beyond the Galaxy, the events with an energy beyond $`\mathrm{3\hspace{0.17em}10}^{18}`$ eV. Recent reviews are in, e.g., . The latest data and their discussions are in  for HiRes, and in  for AGASA.

The main dispute at the present time is the apparent discrepancy between the AGASA data and the HiRes data, with the HiRes data suggesting a downturn of the overall spectrum in qualitative agreement with the GZK-cutoff concept, although at a somewhat higher energy, as would be expected from a highly inhomogeneous source distribution, such as any subpopulation of galaxies. On the other hand, the AGASA data suggest a continuation of the spectrum to higher energies, which would require a new component to appear from under the lower flux contributors, in perfect agreement with a scenario involving the decay of topological defects. The real discrepancy between the two data sets is only about 2 sigma, as discussed at the meeting in Paris June 2002. It is widely expected that AUGER data will resolve this issue, and then maybe the further experiments will be required, such as EUSO and OWL, the European and American space missions to detect airshowers at extremely high energies.

## 11 Conclusions

The successful fit of the new KASCADE data, under the constraint to also fit the higher energy data from HiRes, AGASA and Yakutsk, with the cosmic ray source model originally proposed some time ago, and developed further since then, leads to a number of important conclusions. The most important points are:

* The origin of Galactic cosmic rays may be much closer to a full understanding of their origin. Sofar all quantitative tests, which have been made, show consistency with the model proposed earlier. But much work remains to be done.
* All stars above a zero age main sequence mass of about 25 solar masses converge to a “common final state”.
* Those stars explode with about $`10^{52}`$ ergs in a highly anisotropic explosion, following the mechanism originally proposed by G. Bisnovatyi-Kogan more than 30 years ago, which was then based on a broader suggestion by Kardashev, involving rotation, magnetic fields and potential energy.
* These supernovae may constitute a new very bright standard candle, useful in cosmology, provided we could determine their luminosity integrated over $`4\pi `$ from infrared, polarization or other observations.

## 12 Acknowledgements

P.L. Biermann would like to acknowledge the hospitality, first at the University of Maryland, offered by his colleague Eun-Suk Seo in the early spring of 2002, and then during his prolonged stay at the University of Paris VII, in the late spring 2002, offered by his colleagues Norma Sanchez and Hector de Vega, as well as their hospitality at many meetings in Erice, Paris and Palermo. PLB would especially like to acknowledge a day-long debate with I. Axford at the Taiwan meeting in April 2002, arguing about the agreements and also differences between the superbubble concept and the stellar wind concept. PLB would also like to express his appreciation to N. Langer for a similarly long discussion at Utrecht in May 2002, arguing about supernova physics. And, PLB would also like to acknowledge many discussions with G. Schatz on the intricacies of trying to elucidate cosmic ray spectra from shower data. The material in this paper was the topic of a long helpful debate with T. Stanev in Paris in the spring of 2002, for which PLB is also very grateful. Help with some of the data was generously supplied by Zh. Cao from Utah and D. Bergman from Rutgers. The authors would like to acknowledge with gratitude discussions with many friends and colleagues, especially E.J. Ahn, I. Axford, Zh. Cao, S. Casanova, M. Chirvasa, A. Donea, R. Engel, T. Enßlin, H. Falcke, C. Galea, R. Engel, A. Haungs, H. Kang, M. Kaufman, P.P. Kronberg, N. Langer, A. Lazarian, H. Lee, K. Mannheim, S. Markoff, G. Medina-Tanco, A. Meli, F. Munyaneza, A. Olinto, B. Nath, G. Pavalas, A. Popescu, R. Protheroe, G. Pugliese, J. Rachen, W. Rhode, G. Romero, E. Roulet, D. Ryu, N. Sanchez, K. Sato, G. Schäfer, G. Schatz, H. Seemann, E.-S. Seo, R. Sina, T. Stanev, V. Tudose, H. de Vega, Y.P. Wang, A. Watson, T. Weiler, St. Westerhoff, and Ch. Zier. Work with PLB is being supported through the AUGER theory and membership grant 05 CU1ERA/3 through DESY/BMBF (Germany); further support for the work with PLB comes from the DFG, DAAD, Humboldt Foundation (all Germany), grant 2000/06695-0 from FAPESP (Brasil) through G. Medina-Tanco, a grant from KOSEF (Korea) through H. Kang and D. Ryu, a grant from ARC (Australia) through R.J. Protheroe, and European INTAS/ Erasmus/ Sokrates/ Phare grants with partners V. Berezinsky, M.V. Rusu, and V. Ureche. All these sources of support are gratefully acknowledged. The work reported here was also specifically supported by the DAAD through a fellowship to S. Ter-Antonyan.


text:
# Does the Red Queen reign in the kingdom of digital organisms?

## 1 Introduction

Two major principles of evolutionary biology have been observed in competition experiments between variants of RNA viruses with identical or almost identical fitness: competitive exclusion and Red Queen dynamic . The competitive exclusion principle refers to the ultimate outcome of these competition experiments, and states that when two or more species live on the same resource, eventually all but one will die out . The Red Queen dynamic refers to the initial phase of these competition experiments, where the two competing virus variants both increase in fitness while they remain in roughly equal concentrations. Van Valen  had originally proposed the Red Queen dynamic as a metaphor for the struggle for existence of species in complex ecosystems. These species, just like the Red Queen in Lewis Carroll’s *Through the Looking Glass*, would constantly have to run (that is, to adapt to changing conditions) just to remain where they were (see also Ref. ).

Solé et al.  studied the competition of neutral virus variants in a simple bitstring model with additive fitness, related to Kauffman’s $`NK`$ model , and could confirm both competitive exclusion and Red Queen dynamic in their model. Moreover, Solé et al. showed that the competitive exclusion principle follows immediately from the quasispecies equations  that describe virus evolution.

The Red Queen effect, on the other hand, is not a necessary conclusion of the quasispecies dynamic, as we can see from a simple thought experiment: Assume that we allow two viruses to compete that are both situated on top of the highest peak in the fitness landscape. None of the two viruses can accumulate any further beneficial mutations, and the outcome of a competition between them will be determined by genetic drift. Clearly, the Red Queen dynamic can occur only if beneficial mutations are sufficiently abundant, so that adaptation happens on a faster time scale than genetic drift. The additive model of Solé et al. has a fairly high rate of positive mutations for all but the very best sequences in the fitness landscape, which explains why they observed the Red Queen dynamic. Such simple additive or multiplicative fitness landscapes lead to a smooth increase in average fitness over time , and such increase has been reported for RNA viruses . However, in many cases beneficial mutations are rare, which leads to the frequently-oberseved stepwise increase in average fitness .

Here, I present a study of the influence of the immediate neighborhood in the fitness space on the competition dynamics of identical organisms. I carried out this study with digital organisms (self-replicating computer programs) , using the Avida platform. Avida is similar to Tom Ray’s Tierra , and has been described extensively in the literature . In Avida, the self-replicating computer programs are rewarded with increased CPU speed when they perform certain logical computations. By changing the bonus structure, that is, by changing which computations are rewarded with what bonuses, the researcher can shape the fitness landscape in which the digital organisms evolve. I studied three different strains of digital organisms that were located on the top (A), at the base (B), and some distance away from the nearest fitness peak (C) (Fig. 1).

## 2 Materials and methods

### 2.1 Computer experiments

I used Avida version 1.99, which is available from http://sourceforge.net/
projects/avida. I used the default setup, apart from the following modifications: I switched off length changes, and set insertion and deletion mutations to zero. As the original ancestor of all organisms in this study, I used the handwritten organism organism.heads.100, which comes with the Avida distribution. I used two different environment files (the environment determines the logical operations for which organisms are rewarded, and thus defines the fitness landscape). The first one, environment 1, did not reward any logical computations, so that fitness was a direct function of an organism’s replication efficiency. The second one, environment 2, rewarded all possible one-, two-, and three-input logical operations. In this environment, replication efficiency was only a minor component of fitness, whereas the successful completion of as many logical computations as possible resulted in high fitness values. The mutation rate in all experiments was 0.75 miscopied instructions per offspring organism, unless explicitly stated otherwise.

I created the strains A, B, and C (as illustrated in Fig. 1) as follows. For strain A, I inocculated a population of $`N=5000`$ organisms with the handwritten ancestor, and propagated this population for 14,000 generations in environment 1. Then, I extracted the most abundant organism, which replicated roughly two times faster than the original ancestor. By generation 14,000, the population had not experienced any significant fitness improvements for over 10,000 generations, which showed that it had indeed reached the top of a fitness peak. For strain B, I first evolved from the handwritten ancestor an organism that could perform all possible one- and two-input logical computations. Then, I inocculated a population of size $`N=200`$ with this evolved organism, in order to let the organism accumulate mutations. After 50 generations, I extracted from the population a variant with drastically reduced (by a factor of fifteen) fitness. This variant differed in six instructions (out of 100) from its evolved ancestor, and had lost five logical operations. It was therefore at the base of a fitness peak in environment 2, where many additional logical computations could be learned. For strain C, I simply used the handwritten ancestor (also in environment 2). The handwritten ancestor did not perform any logical computations. In Avida, the evolution of the first logical computation is harder than the evolution of additional ones, which guaranteed that C was further away from the fitness peak than B.

I carried out all competition experiments with the following protocol. For each strain, I used three different population sizes, $`N=100`$, $`N=1000`$, and $`N=10,000`$. I did 100 replicates per population size and competing strain. In each competition experiment, I inocculated the population with $`N`$ identical copies of the strain under study, and marked half of them with an inheritable neutral marker. Then, I propagated the population in the respective environment (environment 1 for strain A, environment 2 for strains B and C) at fixed size $`N`$ until the population consisted of either all marked or all unmarked organisms. Each generation, I recorded the average and maximum fitness of the marked and unmarked organisms, as well as their relative frequencies.

### 2.2 Coalescent theory

For the competition of two strains in a completely flat landscape, we can calculate the distribution of extinction times from coalescent theory. Let $`p_j(t)`$ be the probability that after $`t`$ generations, each of the $`N`$ organisms in the population is a direct descendant of one of a group of only $`j`$ organisms that were present at $`t=0`$. Further, let $`m_j=[1(1/2)^{j1}]`$ be the probability that this group of $`j`$ contained both marked and unmarked organisms. Then, the probability that neither the marked nor the unmarked strain is extinct after $`t`$ generations is

$$P(t_{\mathrm{extinct}}>t)=\underset{j=2}{\overset{N}{}}m_jp_j(t).$$
(1)
The quantity $`p_j(t)`$ has been given in various forms in the literature . The following version is based on the result by Tavaré :

$$p_j(t)=\underset{k=j}{\overset{N}{}}\rho _k(t)(1)^{kj}(2k1)C_{jk},$$
(2)
where $`\rho _k(t)=\mathrm{exp}[k(k1)t\sigma ^2/(2N)]`$, and $`C_{jk}`$ satisfies the recursion (for $`kj`$):

$`C_{11}`$ $`=1,`$ (3a)
$`C_{jk+1}`$ $`={\displaystyle \frac{(Nk)(k+j1)}{(N+k)(kj+1)}}C_{jk},`$ (3b)
$`C_{j+1k}`$ $`={\displaystyle \frac{(k+j1)(kj)}{j(j+1)}}C_{jk}.`$ (3c)
The quantity $`\sigma ^2`$ is the variance of the offspring distribution. In Avida, it has the value $`\sigma ^2=2/3`$, because organisms in Avida can have either zero, one, or two offspring organisms per generation, with equal probabilities (if all organisms have equal fitness).

## 3 Results and Discussion

I measured the time to extinction of the losing clone in all competition experiments and derived the distributions of extinction times. The nine distributions (three strains at three different population sizes) are shown in Fig. 2. The extinction times are the longest for competitions of strain A, intermediate for competitions of strain C, and shortest for competitions of strain B. For strains A and B, the extinction times grow systematically with increasing population size, whereas for strain C, the extinction times grow from $`N=100`$ to $`N=1000`$, but then decrease again for $`N=10,000`$. (Notice that the line with solid diamonds lies on the left of the line with solid squares in Fig. 2b.)

Typical competition dynamics for the three strains are displayed in Figs. 3-5. In all three figures, I show the average and maximum fitness of the winning and the losing clone as a function of time, as well as the relative concentration of the winning clone as a function of time.

For strain A, average and maximum fitness of the two competing clones remained typically almost unchanged until the losing clone disappeared (Fig. 3). Changes in relative clone size were mostly due to genetic drift. However, the theoretical prediction Eq. (1) works only moderately well for the smallest population size, and overestimates the extinction times for the larger population sizes substantially (Fig. 2). I verified that this effect is not a shortcomming of the theoretical prediction by carrying out additional competition experiments with all mutations switched off (data not shown). For these experiments, the extinction data were in very good agreement with the theoretical prediction, which implies that the reduced extinction times in the competitions with mutations must be due to the accumulation of deleterious mutations. Indeed, in Fig. 3, we see that the maximum fitness of the losing clone is—after approximately 500 generations—consistently lower than the maximum fitness of the winning clone. Even though the difference in maximum fitness between the two clones is very small, it is sufficient to accelerate the extinction process. And the smaller the losing clone becomes, the more likely it is to experience even further reductions in its maximum fitness. The final stage of the competition is mutational meltdown : Decreasing clone size accelerates loss of the highest fitness mutants, which in turn results in even further reduction of clone size. The clone decreases in size and loses fitness at an ever accelerating pace, until it has disappeared.

Fig. 3 clearly shows that mutational meltdown takes place towards the end of the competition and leads to a reduced extinction time. At first glance it is surprising that mutational meltdown should be responsible for the increasing deviations between theory and measured extinction times as the population size increases. After all, mutational meltdown is commonly associated with small population sizes. The reason why this effect here becomes more pronounced at larger population sizes is the following: When the relative difference in fitness between two mutants is smaller than the inverse of the population size, then these two mutants are effectively neutral, in the sense that they are equally affected by genetic drift. Therefore, larger population sizes can resolve finer fitness differences between mutants. In the case of strain A, the fitness difference between the winning and the losing clone is miniscule, so that at small population sizes, drift is the dominant factor. Once the population size is sufficiently large, however, this small fitness difference turns the population dynamic deterministic, and the clone that loses the fittest mutant first will slowly but surely disappear.

A competition between two clones of strain C typically started out just like one between two clones of strain A. However, often one of the two clones managed eventually to acquire a beneficial mutation with substantial selective advantage, and would then quickly exclude the other clone. Since beneficial mutations were fairly rare for strain A, the second clone did not have the chance to pick up an even better mutation in the short time that remained before it died out. Clearly, the time to the first beneficial mutation determined therefore the distribution of extinction times, unless it was much larger than the typical time to extinction by genetic drift. The time to the first beneficial mutation grows with decreasing population size, while the time to extinction by drift grows with increasing population size. The distributions of extinction times for strain C are determined by these two constraints: For $`N=100`$, beneficial mutations are very rare, and the extinction times are dominated by the effects of drift. For $`N=1000`$, beneficial mutations are still rare, but nevertheless sufficiently abundant, so that the extinction times are clearly shorter than the ones for drift alone. Finally, for $`N=10,000`$, beneficial mutations are so frequent that the time to extinction is on average even shorter than for $`N=1000`$.

Strain B showed a competition dynamic very similar to the one described by Solé et al. . Both the marked and the unmarked clone gained substantially in fitness during the competition, and both clones would alternatingly take the lead in fitness gains (Fig. 5). However, this Red-Queen dynamic came at a price: The time to extinction of either clone was consistently shorter than for strains A or C. Apparently, the constantly changing relative growth rates of the two competing clones introduced increased fluctuations in the clone sizes, so that one of the clones was quickly reduced to a size at which it became prone to mutational meltdown, or was at least substantially impaired in its ability to acquire further beneficial mutations.

## 4 Conclusions

The location in fitness space from where a strain is taken has a strong influence on its competition dynamic. A clear arms race of mutants with ever increasing fitness can only be observed when beneficial mutations are abundant. When beneficial mutations are rare or completely absent, then either the clone that finds a beneficial mutation first wins or the clone that loses the highest-fitness mutant first loses. In general, it seems that a positive mutation rate will always reduce the competition time, so that the loser dies out earlier than it would in the absence of mutations. The Red-Queen dynamic, where both clones acquire mutants of ever increasing fitness, is particularly unstable. In this case, competitions last the shortest.

The results that I have presented here were obtained in computer experiments with digital organisms. Therefore, it is not a priori clear that my conclusions apply directly to RNA viruses. Nevertheless, I believe it is very likely that they do. In particular, the way in which RNA virus strains are typically prepared for competition experiments (frequent transfers at small population sizes, which leads to accumulation of deleterious mutations, and transfers to new environments, where many advantageous mutations can be acquired) are similar to my preparation of strain B. The fact that they show competition dynamics very similar to that of strain B is therefore reassuring. To test my other predictions in a virus system, one would need strains similar to A or C. In principle, it should be possible to prepare a virus strain which is located at or near the top of a fitness peak, by propagating the virus for a long time in a constant environment at a large effective population size. With such a strain, the competition dynamic should be more similar to my strain A, or maybe C, than B. If this turns out to be true, then competition experiments at various population sizes can be used as a reliable tool to map out the neighborhood in fitness space of a particular virus strain.

## Acknowledgments

This work was supported by the NSF under contract No. DEB-9981397. I would like to thank D. A. Drummond for interesting discussions on coalescent theory, and C. Adami for helpful comments on the manuscript.


text:
# Anomalous magnetotransport in (Y1-xGdx)Co2 alloys: interplay of disorder and itinerant metamagnetism.

## I Introduction

Interplay of structural disorder and magnetic interactions opens a rich field of new physical phenomena. Among them are the actively discussed possibility of disorder-induced Non-Fermi Liquid (NFL) behavior near a magnetic Quantum Critical Point (QCP) as well as a broader scope of effects of disorder on magnetotransport. Hertz76 ; Millis93 ; Varma2001  Structurally disordered alloys Y1-xGdxCo2 are quasi-binary solid solutions of Laves phase compounds YCo2 and GdCo2. The compounds belong to a large family of isostructural composites RCo2. YCo2 is an enhanced Pauli paramagnet whose itinerant Co-3d electron system is close to magnetic instability. In external magnetic field of about 70 T this system undergoes a metamagnetic transition into ferromagnetic (FM) ground state.Goto89  GdCo2 is, on the other hand, a ferrimagnet with a Curie temperature of 400 K in which the spontaneous magnetization of 4f moments is anti-parallel to the induced magnetization of the Co-3d band. Compounds of this family and their alloys provide a convenient ground for experimental studies of magnetotransport phenomena. The electronic structure in the important for the transport vicinity of the Fermi energy is composed mainly of Co-3d states and is, to the first approximation, the same for all compounds of RCo2 family. It has been found that the main contribution to the resistivity of RCo2 comes from the scattering of conduction electrons on magnetic fluctuations due to strong s–d exchange coupling, Gratz95  therefore the transport properties are expected to be especially sensitive to the magnetic state of the sample. GdCo2 occupies a special place in RCo2 family since Gd 4f magnetic moment has no orbital contribution and, therefore crystal-field effects are not important for this compound.

The experimental results on the transport properties of Y1-xGdxCo2 alloys has been published partly in our previous article. Nakama2001  Here we analyze these and new experimental results in order to reveal the physical mechanism of anomalous megnetotransport properties observed in the alloys. In this paper we will discuss the magnetotransport properties of the FM alloys. The properties of paramagnetic alloys will be published elsewhere.

## II Experimental

Samples of Y1-xGdxCo2 were prepared from pure components by melting in an arc furnace under a protective Ar atmosphere and were subsequently annealed in vacuum at 1100 K for about one week. An X-ray analysis showed no traces of impurity phases. A four–probe dc method was used for electrical resistivity measurements. Magnetoresistivity (MR) was measured with longitudinal orientation of electrical current with respect to the magnetic field. The size of the samples was typically about 1$`\times `$ 1$`\times `$10 mm3. Magnetization was measured by a SQUID magnetometer for samples from the same ingot as that used for the resistivity and AC susceptibility measurements.

## III Experimental results

The magnetic phase diagram of the Y1-xGdxCo2 system inferred from the transport and magnetic measurements Nakama2001  is shown in Fig. 1.

Curie temperature $`T_\mathrm{c}`$ decreases with increasing content of Y and eventually drops to zero. A precise determination of the critical concentration $`x_c`$ which separates the magnetically ordered ground state and the paramagnetic region is difficult, since on the onset of the long range order its signatures in the magnetic and transport properties are very weak. The first firm evidence of the long range order are found for alloy with x=0.14 in ac susceptibility at T=27 K, Fig. 2.

Quantum critical scaling theory predicts that when Curie temperature of a FM system continuously depends on an external parameter $`x`$, this dependence is expressed as Millis93 :

$$T_\mathrm{c}\left|xx_\mathrm{c}\right|^{\frac{z}{d+z2}}$$

with critical index $`z=3`$ for a FM system of spatial dimension $`d=3.`$ The experimental Tc vs. $`x`$ dependency does follow this relation, but with additional kink at $`x=x_\mathrm{t}`$. A possible origin of this kink will be discussed later. Linear extrapolation of the phase separation line on the phase diagram Fig. 1 to $`T_\mathrm{c}=0`$ gives as the critical concentration $`x_c=0.12`$. We do not claim however that QCP exists in this alloy system. Direct experimental verification that Tc $``$ 0 as $`x`$ approaches $`x_\mathrm{c}`$ from magnetically ordered state is difficult for a disordered alloy system.

The very surprising result is the positive MR in the FM phase at low temperatures, see Figs. 1 and 3. The well known theoretical result for MR of a localized moment ferromagnet was derived long ago by Kasuya and De Gennes. Kasuya56  As it follows from their theory, MR of a metallic ferromagnet should be negative, having a maximum absolute value at Curie temperature, and approaching zero as T $``$ 0, and in the limit of high temperatures. Qualitatively this behavior has been supported by experiment, as well as by later more detailed theoretical calculations. The present experimental results are in a qualitative agreement with this theoretical behavior only for alloys with large Gd content ($`x0.4`$) (Fig. 3). MR of the FM alloys with composition $`0.3>x>0.14`$ fundamentally differs from this theoretical behavior. Let us note that this composition range falls into the region of the phase diagram between the paramagnetic phase and the additional phase boundary indicated by the kink in the Tc vs $`x`$ dependency, see Fig. 1. MR of these alloys is positive below Curie temperature and is very large.

The known mechanisms of a positive MR can not explain the experimental data. A rough estimate of Lorenz force-driven MR one can get from a comparison with the MR of pure YCo2. Burkov98  In the most pure samples of YCo2 (with residual resistivity of about 2 $`\mu \mathrm{\Omega }`$cm) the Lorenz force-driven positive contribution to the total MR does not exceed 5%. On the other hand the resistivity of the FM alloys at low temperatures falls into the region from 30 to 100 $`\mu \mathrm{\Omega }`$ cm, i.e. at least one order of magnitude larger than the resistivity of pure YCo2. Therefore, according to Koehler’s rule, this mechanism can give MR of only about 0.5%, this has to be compared with the experimental MR of almost 40%.

Weak localization effect is known to give positive MR. However our estimates show that this mechanism can give a contribution which is at least two order of magnitude smaller than the observed MR.

## IV Discussion

The key for understanding the mechanism of the positive MR is a combination of strong dependence of the magnetic susceptibility $`\chi `$ of metamagnetic Co-3d subsystem on the effective magnetic field and of the structural disorder in the R - sublattice of the alloys. In case of GdCo2 as well as in the case of other magnetic RCo2 compounds with heavy R-elements, the 4f–3d exchange interaction is described by introducing an effective field which acts on 3d electrons as:

$$B_{\mathrm{eff}}=n_{\mathrm{fd}}M_\mathrm{f}B$$

where B is external field, Mf is the uniform magnetization of R - sublattice, and $`n_{\mathrm{fd}}`$ \- is the f–d coupling constant (in case of GdCo2 $`n_{\mathrm{fd}}50`$ T/f.u.$`\mu _\mathrm{B}`$Goto2001 ). In Y1-xGdxCo2 alloys the Gd moments are randomly distributed over the R-sites of the crystal lattice. Therefore, the effective field acting on 3d electrons depends on the local distribution of Gd moments and is therefore a random function of coordinate. This random field can be characterized by a distribution function $`P\{B_{\mathrm{eff}}(r)\}`$. The spatially fluctuating effective field induces an inhomogeneous magnetization of Co-3d electron system:

$$m(r)=\chi \left(B_{\mathrm{eff}}\right)B_{\mathrm{eff}}\left(r\right).$$

Therefore even at zero temperature in the ferromagnetic ground state, there are two kind of static magnetic fluctuations in the system: i. $`M_\mathrm{f}(r)`$; ii. $`m(r)`$. These fluctuations give an additional contribution to the resistivity. Since at T=0 K in ferromagnetic phase the 4f magnetic moments are saturated the corresponding contribution to the resistivity does not depend on external magnetic field. On the other hand, the 3d magnetic moment, as we will see later, is not saturated even in the ferromagnetic ground state. Therefore the 3d magnetization does depend on the external field. However, as long as 3d susceptibility is field independent and uniform, the external magnetic field will change only the mean (non-fluctuating part) value of the magnetization, whereas the magnitude of the fluctuations of $`m`$ and corresponding contribution to the resistivity remain unchanged. However actually, the 3d system is close to the metamagnetic instability. Therefore the 3d susceptibility $`\chi `$ is field dependent and gives rise to static magnetic fluctuations which are dependent on magnetic field resulting in non-zero magnetoresistivity.

For a qualitative analysis of this new mechanism of MR we make the following assumptions.
1.We consider the system only in its ground state.
2.Correlations between potential and spin-dependent scattering are neglected.
3.Only the contributions which may be strongly dependent on external magnetic field are retained, e.g. we do not include the effects related to a change of d-density of states in magnetic field, and contributions due to potential scattering and scattering on 4f moments. In fact only the contribution related to the scattering on fluctuations of Co-3d magnetization will be considered.

Figure 4 shows schematically the dependency of the Co-3d magnetization of RCo2 compounds on effective magnetic field Goto2001  and the distribution function $`P\{B_{\mathrm{eff}}(r)\}`$. The metamagnetic transition is indicated by the rapid increase of the magnetization around 70 T.

For the further discussion it is important to have an estimate of the magnitude of the fluctuations of the effective field, i.e. the width of the distribution function $`P\{B_{\mathrm{eff}}(r)\}`$. We can get a hint considering experimental data on field (0 to 7 T) dependency of the magnetization, Fig.5.

Two important points follow from these data: i. there is a non-saturating para-process in these field dependencies above about 2 T; ii. the estimated susceptibility of this para-process ($``$ 0.016 $`\mu _B`$/T) is larger than the susceptibility of 3d system below and above metamagnetic transition ($``$ 0.002 $`\mu _B`$/T), however it is smaller than the susceptibility in the transition region ($``$ 0.04 $`\mu _B`$/T).Goto2001  This indicates, that the scale of the fluctuations is larger than the width of the metamagnetic transition. Therefore we can not treat the fluctuating part of the effective field as a perturbation and have to resort to a phenomenological analysis.

The distribution function $`P(B_{\mathrm{𝐞𝐟𝐟}})`$ describing the fluctuating effective field depends on the alloy composition. For diluted alloys ($`x0`$) the most probable value $`B_{\mathrm{av}}`$ of $`B_{\mathrm{eff}}`$ is close to zero. As $`x`$ increases $`B_{\mathrm{av}}`$ shifts to higher values and in a certain range of the concentrations, the function will have essentially non-zero weight for both $`B_{\mathrm{𝐞𝐟𝐟}}<B_0`$ and $`B_{\mathrm{𝐞𝐟𝐟}}>B_0`$, see Fig. 4. In this case there shall be regions with low and with high 3d magnetization in the sample. The resistivity resulting from this static disorder in 3d magnetization can be expressed as:

$$\rho _m=\rho _{sd}y\left(1y\right).$$

Where $`y=_{B_0}^{\mathrm{}}P\left(B_{\mathrm{eff}}\right)𝑑B`$ is the volume fraction of the high magnetization component. Parameter $`y`$ depends on the alloy composition $`x`$, and on external magnetic field $`B`$. In zero field, point $`y=0`$ corresponds to $`x=0.`$ As the content of Gd increases $`B_{\mathrm{av}}`$ shifts to larger effective fields, and finally, at some alloy composition $`x_t`$, it becomes larger than $`B_0`$, i.e. at this composition almost whole volume is occupied by the high magnetization component, i.e. $`y1`$. With a further increase of Gd content the mean magnetization of Co should increase with a smaller rate, determined by the slope of the $`m(B)`$ dependency above $`B_0`$, however $`y=1`$ in this region. According to this scenario $`\rho _m`$ will increase with $`x`$ at first, reach a maximum value at $`x`$ which corresponds to $`y=0.5`$ and will decrease with further increase of $`x`$ approaching to zero at $`xx_t`$. We believe that $`x_t`$ corresponds to the kink on the phase diagram, Fig. 1. The expected variation of $`\rho _m`$ with the alloy composition is schematically shown in Fig. 6. The total experimental resistivity includes additionally contributions coming from potential scattering and from scattering on 4f moments ($`\rho _0`$), which both are proportional to $`x(1x)`$ and also are depicted in Fig. 6. For comparison the experimental resistivity is shown in this picture too. We present here the resistivity measured at T=2 K in field of 15 T to exclude the contributions related to spin-flip scattering.

The experimental low-temperature resistivity shows the expected variation with $`x`$ (one needs to keep in mind that the relation between $`x`$ and $`y`$ is in general non-linear, especially around magnetic phase boundary $`x_c`$). The resistivity attains the maximum value of about 100 $`\mu \mathrm{\Omega }cm`$ in the region, which corresponds to maximum static magnetic disorder at $`y0.5`$ (room temperature resistivity of the alloys weakly depends on $`x`$ and is about 150 $`\mu \mathrm{\Omega }cm`$). About the same value was obtained for the high temperature limit (maximum magnetic disorder) of magnetic part of the resistivity, arising from scattering on 3d temperature-induced magnetic fluctuations in YCo2.Gratz95  This suggests that the main part of the experimental low temperature resistivity at $`x<x_t`$ originates from the scattering on the magnetic fluctuations of 3dmagnetization, i.e. is identical to $`\rho _\mathrm{m}`$.

Let us show that $`\rho _m`$ is of order of the spin-fluctuation contribution to resistivity at large temperatures, i.e. is of order of 100 $`\mu \mathrm{\Omega }`$cm. In case of YCo2 the spin-fluctuation contribution is the most important and is nearly independent of temperature above about 200 K. Hamiltonian of s-d exchange interaction is given by

$$H_{sd}=G𝑑\mathrm{𝐫𝐬}(𝐫)𝐒_d(𝐫)$$

here $`𝐬(𝐫)`$ and $`𝐒_d(𝐫)`$ are spin density of s- and d- electrons, correspondingly.111Hereafter we take $`\mu _B1`$.

Spin fluctuation contribution to resistivity has form Ueda75 :

$$\rho =\frac{3m}{4ne^2}G^2N_s\frac{1}{T}\underset{0}{\overset{1}{}}𝑑qq^3\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{d\omega \omega }{\mathrm{sinh}^2\left(\omega /2T\right)}Im\chi (q\frac{2k_F}{k_F^{}},\omega )$$
(1)
where $`N_s`$ is density of states of s-electrons, $`\frac{k_F}{k_F^{}}`$ is ratio of Fermi momentum of s- and d- electrons. Dynamic susceptibility $`\chi (q,\omega )`$ is given by the equation:

$$\chi ^1(q,\omega )=\chi ^1\left(q\right)\left(1i\omega /\mathrm{\Gamma }_q\right)$$

here $`\mathrm{\Gamma }_q`$ is damping of the spin-fluctuations, whereas static nonlocal susceptibility $`\chi ^1\left(q\right)`$ is given by

$$\chi ^1\left(q\right)=\chi ^1+Aq^2/N_d$$

with $`\chi =N_d/\zeta (T)`$ , where $`\zeta (T)`$ is inverse Stoner factor, renormalized by spin fluctuations, $`N_d`$ is density of states of d-electrons, and $`A<1`$ is a dimensionless constant.

At large temperatures $`T>\mathrm{\Gamma }_q`$ the expression (1) reduces to:

$$\rho =\frac{3\pi m}{ne^2}G^2N_sT\underset{0}{\overset{1}{}}𝑑qq^3\chi \left(q\frac{2k_F}{k_F^{}}\right)\frac{3\pi m}{4ne^2}G^2N_sT\chi $$
(2)
The last equality is valid when $`2Ak_F/k_F^{}<\zeta (T)<1`$ which must be the case for YCo2. Note that neglecting momentum dependence of susceptibility means that s-electrons see the d-spin fluctuations as point scatterers.

Scattering of conducting s – electrons by the static random distribution of spin density of d-electrons we can estimate in the following way.

Due to random distribution of magnetic moments of Gd, s-electrons experience scattering by

$$\frac{G}{2}S_d\left(𝐫\right)=\frac{G}{2}\chi _m(𝐫𝐫^{})\delta B_{eff}(𝐫^{})$$

random potential. Here $`\chi _m(𝐫𝐫^{})`$ is nonlocal susceptibility near metamagnetic transition. We estimate correlation function of fluctuating effective field as

$$\delta B_{eff}(𝐫)\delta B_{eff}(𝐫^{})=\left(2S_{\mathrm{Gd}}n_{\mathrm{fd}}\right)^2\delta \left(𝐫𝐫^{}\right)x\left(1x\right)a^3.$$

Here $`a^3`$ is volume of the formula unit. In Born approximation the corresponding contribution to the resistivity is given by the expression:

$$\rho _m=\frac{m}{ne^2}G^2N_s\left(2S_{\mathrm{Gd}}n_{\mathrm{fd}}\right)^2x\left(1x\right)a^3\underset{0}{\overset{1}{}}𝑑qq^3\chi _m^2\left(q\frac{2k_F}{k_F^{}}\right)$$
(3)
Both expressions (1) and (3) give the contributions to the resistivity due to scattering on spin fluctuations, however in the first case they are of thermal origin, whereas in the second case the fluctuations are due to randomness of the effective field.

Assuming that nonlocality of the susceptibility is not important we obtain from (2) and (3):

$$\frac{\rho _m}{\rho }=\frac{\left(2S_{\mathrm{Gd}}n_{\mathrm{fd}}\right)^2x\left(1x\right)a^3\chi _m^2}{3\pi T\chi }$$

Using experimental results for $`\chi _\mathrm{m}`$ and $`\chi `$ Goto89 ; Goto2001 ; Burzo72  we find $`\frac{\rho _m}{\rho }`$ in the range from 0.5 to 3, i.e. the resistivity caused by the static magnetic fluctuations is of the same order as the temperature–induced spin fluctuation resistivity. The uncertainty is mainly due to determination of $`\chi _\mathrm{m}`$ near the metamagnetic transition. Taking $`\chi _\mathrm{m}`$ as the susceptibility of YCo2 at B=70 T (at the field of metamagnetic transition)Goto89  gives the upper bound for $`\frac{\rho _m}{\rho }`$. Whereas $`\chi _\mathrm{m}`$ for the disordered alloy of $`x`$=0.18, estimated from our results on M(B), Fig. 5 gives the lower bound.

In external magnetic field, the effective field $`B_{\mathrm{𝐞𝐟𝐟}}`$ decreases 222The effective field decreases for antiferromagnetic 4f-3d exchange, for ferromagnetic exchange (like that in RCo2 compounds with light R-elements) the effective field increases with the external field., therefore $`y`$ also decreases. Depending on the value of $`y_0`$ – the volume fraction in zero field, $`\rho _\mathrm{m}`$ will either increase or decrease, resulting in positive or negative magnetoresistivity: for $`0.5<y_0<1`$ it will be positive, whereas for $`0<y_0<0.5`$ we will have a negative MR. In agreement with the model, the experimental MR is positive at $`0.15<x<x_t`$ and quickly decreases at $`x>x_t=0.3`$ where $`y1`$. Nearly linear field dependencies of MR, observed for $`x<0.3`$, see Fig. 3, implies that the width of $`P(B_{\mathrm{eff}})`$ in this composition range is larger than our experimental field limit of 15 T. The region $`y<0.5`$ almost coincides with the paramagnetic region of the phase diagram. The model predicts negative MR for this region, and this prediction agrees with the experimental result. The model also gives a satisfactory description of the residual resistivity behaviour in this region, see Fig. 6. This agreement suggests that at $`x<x_c`$ the system is actually in spin-glass state. In this region additional essential contributions to MR are present. First, in the paramagnetic region there is a negative MR due to suppression of magnetic disorder in 4f magnetic moment system by external magnetic field. This negative contribution may be the reason why the cross-over point from the positive to negative magnetoresistivity does not coincide with the maximum of resistivity, see Fig. 1. Secondly, there can be additional contributions, both positive and negative, near to zero-temperature magnetic phase boundary due to closeness to QCP.

An independent test of the model is based on the observation that the critical magnetic field $`B_0`$ of the metamagnetic transition in 3d subsystem increases under external hydrostatic pressure (P). Saito99 ; Yamada95  Therefore, basing on the model, we expect that the resistivity of the alloys with the composition left of the resistivity maximum (see Fig. 1) will decrease with pressure, whereas for the alloys right of the maximum it will increase with the increasing pressure. The experimental results for three alloy compositions are shown in Fig. 7.

The sign of the pressure effect is in agreement with the model prediction: the resistivity decreases with pressure for $`x`$=0.1, whereas it increases with P for $`x`$=0.18 and $`x`$=0.3. Moreover, there is a good scaling of pressure and magnetic field dependencies of the resistivity, Figs. 3 and 7 (P$`\alpha `$B) with scaling parameter $`\alpha `$ which is close to literature data on the pressure dependence of B0: $`\frac{dB_0}{dP}1.5`$ T/kBar.Saito99 ; Yamada95

## V Conclusion

It has been found that at low temperatures there is a large contribution to the resistivity related to scattering on magnetic fluctuations in metamagnetic itinerant 3d system, induced by fluctuating effective field of 4f moments. Large positive MR, found in the FM Y1-xGdxCo2 alloys and strong pressure dependence of the resistivity are explained as arising from a combination of static magnetic disorder and strong magnetic field dependence of magnetic susceptibility.

We want to emphasize that this mechanism of resistivity (and of MR) is not material specific, rather it should be common for a broad class of disordered itinerant metamagnets with strong coupling of conduction electrons to the magnetic fluctuations. In Y1-xGdxCo2 the relevant disorder originates from random distribution of d–f exchange fields, however a similar effect should arise when there is a random distribution of local susceptibilities (a corresponding treatment for Kondo systems was recently developed in Wilhelm2001 ). Positive MR observed in FM alloys Y(Co1-xAlx)2 Nakama2000  and in ErxY1-xCo2 Hauser2001  may be explained by this mechanism.

## ACKNOWLEDGMENT

This work is supported by grants 02-02-17671 and 01-03-17794 of Russian Science Foundation, Russia. We thank Dr. P. Konstantinov for stimulating discussions.