text:
# On the effects of the stellar magnetic field on the structure of T Tauri accretion discs

## 1 Introduction

It is now quite widely accepted that in the formation process of low-mass ($`M2M_{}`$) stars, a disc accretion phase is typically present (see Hartmann 1998 for an extensive review and references). In particular, the classical T Tauri stars (CTTS) are thought to be still surrounded by discs and the emission properties of some typical objects suggest that these discs are not merely “passive” dust bodies, irradiated by the central star, but are rather bona-fide energy producing, viscous accretion discs, similar to the ones present in mass transferring close binary systems, and whose modeling dates back to the classical works of Shakura & Sunyaev (1973) and Lynden-Bell & Pringle (1974). Because of their ubiquity and of some interesting theoretical challenges they pose, accretion discs have remained at the centre stage of astrophysical research (see e.g. Frank, King & Raine 1992; Papaloizou & Lin 1995; Lin & Papaloizou 1996 for reviews).

One particular issue in this context is the fact that in many CTTS the measured spectrum presents an unusually high emission in the infra red (called IR excess), far more than expected from disc reprocessing of stellar light or viscous heating of the disc (see Bertout 1989, Beckwith et al. 1990, Hartmann et al. 1994). Our wish to investigate this problem, as well as some other facets of CTTS discs, has been the main motivation for the work reported on in this paper. It consists of setting up a numerical calculation of the properties of an accretion disc around magnetic CTTS, using a scheme of magnetic interaction, originally proposed in the context of diamagnetic accretion by King \[King 1993\].

This description of the magnetic interaction assumes that as material moves through the magnetosphere it interacts with the local magnetic field via a velocity dependent acceleration (force per unit mass, $`f_{\mathrm{mag}}`$) of the form:

$$f_{\mathrm{mag}}=K[𝐮𝐮_f]_{},$$
(1)
where $`𝐮`$ and $`𝐮_f`$ are the velocities of the material and magnetic field lines respectively, $`K`$ is a suitable ”magnetic drag” coefficient (see below) and the suffix $``$ refers to the vector component perpendicular to the field line. The magnetic acceleration, as expressed in equation 1, is intended to represent the dominant term of the magnetic interaction, with $`K`$ containing the relevant parameters determining the effective magnetic time-scale.

Within the above description (1) there are still a number of different possibilities to model the inner disc - magnetosphere interaction. In this paper we shall use the diamagnetic blob accretion (DBA) model, although it is possible to formulate the diametrically opposite case (complete magnetic penetration of the disc) in the general form (1) as well. Wynn, Leach & King \[Wynn, Leach & King 2001\] show how to calculate the appropriate coefficient $`K`$ for the latter case and we plan to repeat our calculations in the future, using this prescription.

In the DBA approach one assumes that the fluid constituting the accretion flow, i.e. the disc and its surroundings, is composed of blobs immersed in a dilute interblob plasma (see Aly & Kuijpers 1990). The blobs behave diamagnetically in the presence of the stellar magnetic field and thus suffer a surface drag force acting on them. The model has since its introduction been applied to various systems, notably to the intermediate polar class of cataclysmic variables (CV). Wynn & King \[Wynn & King 1995\] and Wynn, King & Horne \[Wynn, King & Horne 1997\] incorporated the diamagnetic drag force into a particle hydrodynamical numerical code (HYDISC), originally developed by Whitehurst \[Whitehurst 1988\] to simulate the accretion disc in non-magnetic close binary systems. They found important properties of the white dwarf’s spin evolution and its effect on the disc and explicitly applied their findings to the moderately magnetic CV AE Aqr. A more recent study using this approach was its application to the system EX Hya \[King & Wynn 1999\].

The DBA approach was extended to the study of accreting T Tauri stars for the first time by King & Regev \[King & Regev 1994\], hereafter KR. They performed calculations of individual blob (that is, ballistic) orbits, including the interaction with magnetic loops of the central star and have shown that this mechanism can eject blobs from the system, in directions pointing away from the disc plane. By estimating the angular momentum thus expelled and including it in the overall angular momentum budget KR found that a stellar spin equilibrium value is compatible with the observations (i.e. $``$ an order of magnitude less than the breakup value) if the magnetic loops extend to a few stellar radii. Pearson & King \[Pearson & King 1995\], hereafter PK, generalised the above work into a $`N`$-particle simulation, in a similar manner to the one mentioned above for CVs. This work also focused on the issue of the slow observed T Tauri spin rates and confirmed the idea that such equilibrium spin rates can be achieved by the expulsion of matter from the disc until the corotation radius (see below) coincides with the edge of the magnetic loop. The additional new finding, not anticipated by KR, of this work was that the ejection of material comes to a halt as the star approaches its spin equilibrium value.

The problem of the slow spin rates and rotational evolution of T Tauri stars as a result of the star’s magnetic interaction with its accretion disc has been also approached, in a number of works, by analytical and semi-analytical methods, that is, stopping short of multi-dimensional numerical simulations. Königl (1991), Cameron & Campbell (1993) and Armitage & Clarke (1996) used quite different approaches and all found that a stellar dipole magnetic field of strength $`1\mathrm{kG}`$ is able to regulate the stellar spin to a quasi-static value, in the observed range. In the process of the magnetic interaction the inner accretion disc gets disrupted and the inner (abrupt) termination radius of the disc was estimated in the above papers.

KR and PK used the DBA approach focusing on the ejection of mass from the disc, in an effort to link the low spin rates of T Tauri stars with outflows from young stellar objects (YSO) within a unified scenario. In the present work we utilise the DBA model in the calculation of the properties of the steady (or quasi-steady) accretion disc itself, that is, the material that remains close to the equatorial plane while spiraling in due to viscous torques, after the spin period has already stabilised. Some pecularities of T Tauri light-curves (i.e. the temporal luminosity variations) have already been treated using a model based on single blob ballistic orbits, which leave the the disc plane but ultimately return and impinge on the disc surface \[Ultchin et al. 1997\]. In order to investigate whether the DBA model is consistent with some of the basic properties of T Tauri spectra (see Bertout 1989 and references therein) and in particular the above mentioned IR excess, we have performed numerical simulations using a code based on HYDISC, suitably modified so as to adapt it for the problem at hand.

This paper is organised as follows. In the next section we discuss the basics of blob dynamics and estimate the relevant time-scales. §3 describes the numerical code used in the accretion flow simulation and the procedure for finding the spectrum emitted by such flows. In §4 the results of our simulations for a score of parameter values are described in some detail and finally, §5 summarises this work in comparison to the results of other approaches.

## 2 Diamagnetic accretion flow dynamics

Following KR and PK we model the material of a T Tauri disc as being composed of gas blobs immersed in a tenuous interblob plasma. These blobs are assumed to behave diamagnetically in the presence of the magnetic field of the central star, which threads the plasma. Drell, Folley & Ruderman (1965) showed that the typical time-scale on which a diamagnetic object of mass $`m`$ (i.e. a blob) loses energy, when it moves in a magnetic field $`B`$, is

$$t_\mathrm{d}=\frac{c_Am}{B^2l_b^2}$$
(2)
where $`c_A=\sqrt{B^2/(4\pi \rho _p)}`$ is the Alfvén speed of the interblob plasma, whose density is $`\rho _p`$, and $`l_b`$ is the blob size. There is a minimum condition on the conductivity for this equation to hold and, as it was shown by KR, this condition is satisfied for all reasonable parameters in the conditions appropriate for T Tauri accretion.

This allows to cast the diamagnetic force per unit mass, acting on a blob, into one simple formula, whose form is in accordance with equation 1:

$$𝐟_d=K_d[𝐮_b𝐮_f]_{}\widehat{𝐧},$$
(3)
where $`K_d=1/t_\mathrm{d}`$ (see below) is the diamagnetic drag coefficient. The magnetic drag force obviously vanishes if the blob velocity $`𝐮_b`$ relative to the velocity of the local field line $`𝐮_f`$ is parallel to the field lines, it is proportional to the perpendicular (to the local field) component of this relative velocity and is directed in the normal direction $`\widehat{𝐧}`$ to the field. This is the meaning of the shorthand vector notation in equation 3. Its explicit component form will be given below.

Assuming that the magnetic field rotates with the star and does not change on the time-scale of interest, the above equation is best explicitly written in the frame of reference rotating with the central star and had the form

$$𝐟_d=K_d[𝐯(𝐯𝐛)𝐛],$$
(4)
where $`𝐯`$ is the blob velocity in the rotating frame, that is, relatively to the field lines, and $`𝐛`$ is a unit vector in the local field direction (see KR).

Substituting explicitly the above form of $`c_A`$ and $`m=l_b^3\rho _b`$, where $`\rho _b`$ is the blob density, in the drag coefficient one gets

$$K_d(𝐱)=2\sqrt{\pi }\rho _p^{1/2}(\rho _bl_b)^1BKf(𝐱),$$
(5)
where the second equality parametrizes $`K_d`$ in terms of a constant coefficient $`K`$ and a spatially dependent ($`𝐱`$ denotes the position vector) function, which is of order $`1`$. The blob parameters and their space dependence, as well as the plasma density, are rather uncertain. The only variable which can be explicitly written is the magnetic field (in the approximation that it is fixed in the rotating frame, that is, anchored to the star), because it can be chosen a priori.

In this work we shall mainly use a dipolar magnetic field configuration and group, for the sake of convenience, all the other variables into a single fixed parameter, incorporating in it most of our ignorance of the complex plasma instabilities, which are supposed to create and shape the blobs. Within this framework we have

$$K_d=K(R/R_{})^3[1+3(z/R)]^{1/2},$$
(6)
in cylindrical coordinates $`(R,z)`$ whose $`z`$ axis coincides with the dipole axis. $`R_{}`$ is the stellar radius and

$$K=\sqrt{\pi }\rho _p^{1/2}(\rho _bl_b)^1B_0,$$
(7)
where $`B_0`$ is the polar field strength. This $`K`$ is assumed to be constant and will serve as a parameter of our problem. For a discussion of the different options to chose $`K_d`$ see e.g. Wynn, King & Horne (1997).

We focus in this work on the case of an aligned dipole, that is, when the dipole axis coincides with the disc axis, but will briefly discuss (in §4.3) the inclined dipole case as well.

The magnetic drag term has to be added to the blob equations of motion, which include also gravitational, centrifugal and coriolis terms (in the rotating frame). Scaling the independent variables by their natural values, that is, lengths by the corotation radius, $`R_{\mathrm{co}}(GM_{}/\mathrm{\Omega }^2)^{1/3}`$ (with $`M_{}`$ the central star’s mass and $`\mathrm{\Omega }`$ its rotational angular velocity) and time by $`1/\mathrm{\Omega }`$, the following blob equations of motion are found \[Ultchin et al. 1997\]:

$$\ddot{x}=2y+x\frac{x}{R^3}k_d(𝐱)[\dot{x}(\dot{𝐱}𝐛)b_x],$$
(8)
$$\ddot{y}=2x+y\frac{y}{R^3}k_d(𝐱)[\dot{y}(\dot{𝐱}𝐛)b_y],$$
(9)
$$\ddot{z}=\frac{z}{R^3}k_d(𝐱)[\dot{z}(\dot{𝐱}𝐛)b_z],$$
(10)
These constitute the nondimensional equations of motion for the cartesian components, in the rotating frame, of the blob trajectory $`𝐱(t)`$ and are thus ready to be incorporated into HYDISC (see in the next section). $`R`$, the cylindrical radius, is obviously $`R(x^2+y^2)^{1/2}`$ and the calculation of $`𝐛`$ is quite simple for an aligned dipole. For a tilted dipole, a somewhat more involved (but straight forward) calculation is needed.

The drag term has been written as $`k_d`$, the nondimensional version of $`K_d`$ (of equation 6). The explicit form of this function, following from equations 6 and 7, brings out the physical meaning of the constant coefficient $`k`$ appearing in front of the nondimensional spatially dependent part of the magnetic field

$$k_d(𝐱)=k(R/R_{})^3[1+3(z/R)]^{1/2},$$
(11)
which is

$$k=2.82\left(\frac{B_0}{100\mathrm{G}}\right)\left(\frac{\rho _p}{10^{10}\mathrm{g}\mathrm{cm}^3}\right)^{1/2}\left(\frac{l_b\rho _b}{100\mathrm{g}\mathrm{cm}^2}\right)^1\left(\frac{P_{}}{10^6\mathrm{s}}\right)\frac{t_{\mathrm{rot}}}{t_{\mathrm{mag}}(R_{},100\mathrm{G})}.$$
(12)
Here reasonable (order of magnitude) estimates of the parameters for accreting T Tauri systems have been used to scale the relevant variables. $`P_{}t_{\mathrm{rot}}`$ is the stellar spin period and $`B_0`$ the polar magnetic field strength. The disc density scaling value is estimated from models of standard accretion discs around T Tauri stars, using a representative accretion rates of $`\dot{M}510^8M_{}\mathrm{yr}^1`$ \[Basri & Bertout 1989\]. The blob parameters are rather uncertain, but it is reasonable that the blob size is probably of the order of the vertical scale height in the disc while its density should be several orders of magnitude larger than that of the disc. This gives the value with which we scale $`l_b\rho _b`$. In any case, our results will be seen to depend only on the parameter $`k`$ (and several other quantities related to the numerical model) and hence on the particular combination given in equation 12.

Time-scale estimates can now be easily seen to demonstrate that for $`k1`$ the magnetic drag time-scale near the star (where $`RR_{}`$ and thus the contribution of the space dependent part, arising from eq. 11, is very close to $`1`$) for a polar field of $`100\mathrm{G}`$, denoted here by $`t_{\mathrm{mag}}(R_{},100\mathrm{G})`$, is much larger than the dynamical time-scale there, since the latter is much shorter than $`t_{\mathrm{rot}}`$ (T Tauri stars are slow rotators). Thus we do not expect a significant effect near the star due to the magnetic drag in this case. Since the magnetic field decays (and thus the magnetic drag time scale grows) faster with $`R`$ than the dynamical time-scale, the effect of the magnetic field on the motion will be small everywhere for $`k1`$

Let us examine now what happens at the corotation radius. From the functional dependence of the dipole field (eq. 11) we have the term $`(R_{\mathrm{co}}/R_{})^3=(\mathrm{\Omega }/\mathrm{\Omega }_K)^2`$, with the ratio of the $`\mathrm{\Omega }`$-s evaluated at the stellar surface. This term is of the order of $`10^2`$ in T Tauri stars, but $`t_{\mathrm{rot}}`$, by definition, is equal to the dynamical time scale at the corotation radius. If one chooses $`k100`$, a significant effect of the magnetic force at corotation is expected and obviously more so inside corotation. Thus we expect a significant modification of the inner accretion disc in this case.

These time-scale estimates provide the rationale for our choice of the $`k`$ parameter range (or $`B_0`$ with all the other parameters constant) in the simulations (described in §4). Note that the viscous time-scale did not enter our considerations here. In accretion disc theory it is always assumed that this time-scale is much larger than the dynamical time-scale and consequently if the the magnetic time-scale approaches the dynamical one they are both much shorter than the viscous time scale. However, in regions where the magnetic time-scale is extremely long (far out in the disc, or even quite close in, if the field is weak enough) it might be even longer than the viscous time-scale and thus the magnetic drag would not interfere with the viscous evolution of the disc.

## 3 The numerical method

### 3.1 General considerations

To facilitate an accretion flow simulation in our system we adapt HYDISC, a numerical particle code originally written by Whitehurst \[Whitehurst 1988\] to study accretion in CV, to the problem of a single central magnetic accretor (a T Tauri star). The introduction of the magnetic interaction into the particles’ equations of motion, using the DBA approach, has already been explained in the previous section. Our blobs are being considered as the particles of the particle scheme, utilized in HYDISC to simulate the hydrodynamical aspects of the flow, like pressure and viscosity. As it was already mentioned in the Introduction, similar schemes based on HYDISC have been designed for simulating accretion flows in magnetic CV \[Wynn & King 1995, Wynn, King & Horne 1997\] and also for investigating the outflows and spin evolution of an accreting T Tauri star \[Pearson & King 1995\].

Our calculation is, in principle, similar to the latter work because the gravitational effects of the companion star are removed from the code, however the purpose of this work and therefore the simulation details are quite different. We are interested in obtaining a steady-state structure of the disc itself and the spectral shape of the emitted radiation after the spin period of the central star has settled onto its equilibrium value. This calls for the determination of an appropriate spatial computational domain, a suitable particle injection scheme (initializing the calculation) and the need for a method to assess the temperature (and from it, the integrated spectrum) of the accreting matter.

The inner and outer boundaries of our computational region are the stellar radius, $`R_{}`$, and the ”escape” radius, $`R_{\mathrm{out}}`$ (which is a parameter typically set at $`10R_{\mathrm{co}}`$), respectively. Particles are injected into the computational domain (in a manner described below) and if during the course of the evolution they are found to leave this domain, they are removed from the calculation altogether. Once in the computational zone, each particle is subject to the field forces formulated explicitly in the rotating frame by the equations of motion 8 through 10. In addition, short-range forces, reflecting the interaction between adjacent particles and thus simulating the hydrodynamical aspect of the dynamics, are employed using the HYDISC version of the particle hydrodynamics approach \[Whitehurst 1988\]. These forces are controlled by three parameters: $`C`$, determining the pressure force; $`Q`$, a constant of order unity and $`r_𝐯`$, the ”chaining” mesh size. Physically $`C`$ is the sound speed, $`Q`$ is in fact the coefficient of restitution of the particle-particle collisions and is less than $`1`$ for inelastic collisions, which must be assumed to simulate viscous flows, and $`r_𝐯`$ is the viscous scale-length (for details, see Whitehurst 1988).

The fact that HYDISC separates the treatment of the long-range (field) forces from the short-range (hydrodynamical) ones makes this code convenient for application to our problem. Most changes are introduced in the field part (as explained in the previous section on the equations of motion) and in the initialization schemes, that is, particle injection (see below), while the treatment of particle interactions remains essentially the same as in HYDISC, with an appropriate choice of the parameters.

The structure of HYDISC, as explained above, allows also to estimate, relatively easily, the energy dissipation (see below).

### 3.2 Particle injection and disc initialization

Our problem involves an accretion disc around a single star and thus the injection scheme has to be chosen accordingly. The disc around a T Tauri star is not being continuously replenished by an external stream, as is the case in discs fed by the secondary star in a close binary system, like a CV. This said, we still may chose between two distinct possibilities: a ”burst” injection, in which all particles are injected simultaneously (or almost so) and a continuous injection. In both cases we have also to select the spatial distribution of the injected particles. Clearly, in order to obtain meaningful results the ultimate configuration should not depend on the initial conditions. Computation time limitations dictate that the computational zone must be limited (as well as the total number of particles). Thus, the seemingly best option, i.e. to continuously inject particles at a radius much larger than $`R_{\mathrm{co}}`$ and wait for the disc to be built up by viscosity, is prohibitive. We have adopted instead what we call a ”modified burst” injection scheme.

The particles were injected at an equatorial ring around a radius $`R_{\mathrm{inj}}`$ with a Gaussian spread in the individual particles injection radii and a small Gaussian displacement out of the disc plane as well. Each particle was initially assigned an azimuthal Keplerian velocity plus a random contribution, appropriate for the local disc temperature. The temperature of the disc, used to calculate the above mentioned random contributions to the injected particles’ velocities and the sound speed (and thus the constant $`C`$), was estimated using the analytical standard disc model. The actual formula which we have used is identical to formula 6 of PK:

$$T(r)=T(3)\left[\frac{27}{13^{1/2}}r^3(1\sqrt{r})\right]^{1/4},$$
(13)
where $`r`$ is the radius in units of $`R_{}`$ and T(3), is the temperature at 3 stellar radii, where true estimates are available \[Basri & Bertout 1989\]. Although the temperature in our disc is somewhat different, this estimate actually influences only the pressure in the disc (in addition to the initial condition) and therefore the error it introduces should be negligible (see PK).

In order not to lose too quickly significant numbers of particles from the simulations, the ”modified burst” scheme consisted of extending the burst duration to a finite (typically a few stellar rotation periods) time $`\tau _\mathrm{b}`$, with the injection rate decaying with time in a Gaussian manner. In addition, particles were added at a constant (but rather slow as compared to the burst) rate $`\dot{N}_{\mathrm{cont}}`$ so as to compensate for particles lost from the computational zone (accreted or expelled). The spatial distribution of the added particles was kept fixed in time. Thus the modified burst scheme actually gave rise to the following number of injected particles as a function of time (expressed here in units of the rotation period)

$$N_{\mathrm{inj}}(t)=N_\mathrm{b}(t)+\dot{N}_{\mathrm{cont}}t,$$
(14)
where $`N_\mathrm{b}(t)=_0^t\dot{N}_0\mathrm{exp}(t_1^2/2\tau _\mathrm{b}^2)𝑑t_1`$ . Typically, in our simulations the total number of particles in the burst (i.e. $`N_\mathrm{b}(t=3\tau _b)`$ say) was between 5000 and 10000 and $`\dot{N}_{\mathrm{cont}}=150`$ particles per rotational period (see below).

### 3.3 Calculation of the emitted spectrum

With the present DBA scheme it is impossible to obtain a detailed thermal structure of the accretion flow and therefore we aim at approximating just the effect of the magnetic field on the emitted spectrum, relatively to the non-magnetic ($`k=0`$) case. We assume that the disc is optically thick in the region of interest and that the energy dissipated in it by the viscous interaction is immediately locally radiated out in the vertical direction. These assumptions are essentially the same as in the standard Shakura & Sunyaev treatment of geometrically thin and optically thick accretion discs. What is different here is the fact that the viscous dissipation rate must be computed from the details of the numerical particle scheme and can not be found directly from the disc’s viscous flow. As explained above, a part of the kinetic energy of the particles (blobs) is lost in collisions as these are inelastic ($`Q<1`$). After dividing the computational zone into cells, by a two-dimensional grid in the disc plane, one can calculate the energy released in each cell $`\mathrm{\Delta }E_{ij}`$ during a given time $`\mathrm{\Delta }t`$, provided the masses of the particles are known. The composite spectrum can then be found by adding up the contributions of all the cells in the computational domain, with each cell assumed to contribute a black-body emission spectrum, corresponding to its effective temperature

$$T_{ij}^{\mathrm{eff}}=\left(\frac{\mathrm{\Delta }E_{ij}}{\mathrm{\Delta }tA_{ij}\sigma }\right)^{1/4},$$
(15)
where $`A_{ij}`$ is the cell surface area and $`\sigma `$ is the Stefan radiation constant.

As explained in §2 (see equations 11 and 12) we have lumped the unknown values of the blob properties into a single parameter and chosen a specific spatial dependence for it. Thus for a given value of the coefficient $`k`$, different combinations of the constituent parameters (and among them the blob’s mass) is allowed. Assuming some fixed average blob mass can thus give a spectral shape, which should be considered as a qualitative estimate and have a relative meaning, that is, the spectral shapes for different values of $`k`$ can be perceived as an assessment of the effect of the magnetic field strength on the shape of the spectrum. In addition, since our computational zone does not include (due to computing time limitations) the extended outer portion of the disc, the results must be interpreted remembering this fact. In particular, since the outer cool parts of accretion discs contribute mainly to the long wavelength part of the emitted spectrum, our computed spectra clearly underestimate the spectrum strength for large $`\lambda `$ (see below in Fig. 3).

## 4 The simulations

### 4.1 Simulation parameters

In choosing the parameters of the simulations we have tried to optimize the sometimes conflicting demands to obtain, on the one hand, meaningful results and on the other hand to limit, as much as possible, the CPU demands of the numerical calculations.

We have set the stellar parameters in all the simulations to be $`M_{}=M_{\mathrm{}}`$, $`R_{}=1.5R_{\mathrm{}}`$, $`P_{}=10^6\mathrm{s}11\mathrm{days}`$, giving for the corotation radius $`R_{\mathrm{co}}\mathrm{1.5\hspace{0.33em}10}^{12}\mathrm{cm}14R_{}`$. As explained before, the corotation radius serves as our length unit.

Experimenting with different injection radii (the location of the centre of the injection ring), we have found that with the injection scheme (14) an injection radius of $`R_{\mathrm{inj}}=1.4R_{\mathrm{co}}`$ is effective in building up rather quickly an accretion disc, which does not differ essentially from the one resulting (after a longer calculation) from a larger injection radius. In these tests we have shut off the magnetic interaction, but remembered from previous works that the injection radius should not be too close to corotation. Consequently in all the simulations we have used this value of the injection radius, with a Gaussian spread of up to $`\pm 0.1R_{\mathrm{co}}`$.

As discussed in the previous section, the parameters related to the particle scheme treatment of the pressure and viscosity include $`C`$, $`Q`$ and $`r_\mathrm{v}`$. We chose $`C`$ to be equal to the isothermal sound speed, derived from the standard Shakura & Sunyaev disc temperature distribution. $`Q`$ and $`r_\mathrm{v}`$ determine the viscous effects and have to be chosen with care. $`Q`$, the restitution coefficient has a global effect on the simulation and since its reasonable values should lie between $`0.5`$ and $`1`$ ($`Q<0.5`$ allows particle interpenetration and $`Q=1`$ corresponds to elastic collisions). We have found that changing $`Q`$ within the above mentioned range has only a small effect on the final quasi steady state, but it influences the length of the time until such a state is reached. Whitehurst (1988) reached a similar conclusion by experimenting with $`Q=0.6`$ and $`0.7`$ values in the original HYDISC code. Thus, wishing to shorten the simulation time, without, however, making the viscosity unphysically large, we have opted for $`Q=0.7`$ and used this value in all our simulations, save for a few test cases.

The parameter $`r_\mathrm{v}`$, setting the viscous interaction length scale, has a local effect on the simulation. When deciding on the range of variation for this parameter we have taken into account the fact that too small a value of $`r_\mathrm{v}`$ causes the viscous evolution time to be unacceptably long. In addition, because the number of cells in the computational range becomes large, one has to increase the number of particles in order to keep the particle treatment meaningful. Both effects increase the CPU time and we have found that for $`r_\mathrm{v}0.02`$ we were unable to efficiently compute the disc evolution. On the other hand for $`r_\mathrm{v}`$ approaching the value of $`0.1`$ the increase in viscosity causes the particles to clump together and the evolution becomes unrealistically viscosity dominated. Performing test calculations with no magnetic interaction ($`k=0`$) we found that an initial ring injected around $`1.4R_{\mathrm{co}}`$ viscously spreads to a structure very similar to a standard thin disc in a time of $`100300`$ rotational periods when $`r_\mathrm{v}`$ is kept within the range $`0.08r_\mathrm{v}0.03`$. Using such a calculation it is also possible to find what is the effective $`\alpha `$ (the viscosity coefficient in the Shakura & Sunyaev prescription) with a given choice of $`Q`$ and $`r_\mathrm{v}`$. A test calculation with $`r_\mathrm{v}=0.05`$ (in units of $`R_{\mathrm{co}}`$, see below) and $`Q=0.7`$ was very similar to the evolution of a ring of matter to as it is found from a standard exercise with an $`\alpha `$ viscosity (see e.g. Pringle 1981) value of $`0.1`$.

The parameters $`k`$ and $`r_\mathrm{v}`$ were varied in the different simulations. The range of the variation of $`k`$ was chosen so as to correspond to polar stellar magnetic value of up to several kiloGauss (when the other parameters in $`k`$ are set to their representative values, see equation 7) and $`r_\mathrm{v}`$ was chosen within the range mentioned above.

### 4.2 Results of representative cases

We have performed 12 full simulations of an accretion disc in an aligned dipole magnetic field for different choices of the pair of parameters $`k`$ and $`r_\mathrm{v}`$. The simulations were initialized using the injection scheme described in equation 14 with parameters chosen so as to insure that the number of particles does not drop below $`10^4`$ during the simulation. The simulations were continued until a quasi-steady configuration was achieved, that is, no significant changes of the structure were observed over a significant time (see below). We have chosen to describe here in some detail the simulations with $`r_\mathrm{v}=0.05`$ (simulation group 3 in our list), that is, with the viscosity strength chosen in the middle of its feasible range. As discussed in the previous section this choice corresponds to an effective $`\alpha 0.1`$. Three runs were performed for different values of the magnetic interaction strength ($`k=0,10,100`$), where the non-magnetic case ($`k=0`$) was run to serve as reference, in particular for the spectra, which have quantitative meaning only in comparison to each other (as explained before).

These three runs were all allowed to continue for a total simulation time of $`250P_{}`$. A quasi-steady state was achieved already after approximately half of this time. The CPU time requirement for each run, with the number of particles kept at $`10^4`$ (the accreted and ejected particles were replenished as explained above), was about $`100`$ hours on a Sun ULTRA 2 workstation, when the program was coded in C, to increase its efficiency. In Fig. 1 we show the surface number density of the particles (a quantity proportional to the disc surface mass density, customarily called $`\mathrm{\Sigma }`$) as a function of radius for the three values of $`k`$. Note that in the innermost and outermost parts of the disc, regions in which the surface density curves are dashed, the number of particles is too small to render the result to be quantitatively reliable. We also plot on this graph the analytical result, arising from a standard thin disc formula with inner cutoff (see below). This is represented by the fully dashed curve, which which follows the $`k=100`$ case for most of the range.

From this figure it is evident that the effect of the magnetic field on the disc structure is important, as it significantly modifies the functional dependence of the surface density on the radius. For $`k=0`$ (no magnetic field) the particles form a rather ”collapsed” with its density growing sharply towards the central object. This is so because our simulations were tailored for the high $`k`$ case and therefore the particle injection rate was too small for the buildup of a steady disc in a small $`k`$ case. A growing magnetic field strength has a dynamical effect on the disc. As could be expected from a previous study, in which individual orbits of particles were calculated \[Ultchin et al. 1997\], and from the DBA simulations of PK, the magnetic interaction drags blobs residing inside corotation towards the star and pushes out blobs whose orbits are outside corotation. These effects are modified here by the fact that the dipole magnetic field decays outwards and by the interblob interactions. The combined effect is a significant depletion of blobs close to the star (a ”hole”) and an increase of density at large radii (a ”flattening” of the functional dependence of the surface density on radius). The size and depth of the hole and the extent of the flattening of the surface density obviously depend on $`k`$, That is, on a combination of the disc density (which depends on the accretion rate) and the strength of the magnetic field.

Analytic estimates of the hole radii at spin equilibrium, as given by Cameron and Campbell (1993), predict that the disc should end abruptly close to the corotation radius (but depending also on the fastness parameter of the star). Our numerical results do not yield an abrupt disc termination, but the value of $`R_{\mathrm{co}}`$ is certainly consistent with these results. In the following comparisons of our results to analytical solutions we have used this value of the inner cutoff radius, that is, the analytical formulae giving rise to the fully dashed curves in Figs. 1 and 2 (see below) result from

$$\mathrm{\Sigma }(R)=\frac{\dot{M}}{3\pi \nu }\left[1\sqrt{R_{\mathrm{co}}/R}\right]$$
(16)
and

$$T^4(R)=\frac{3GM\dot{M}}{8\pi R^3\sigma }\left[1\sqrt{R_{\mathrm{co}}/R}\right],$$
(17)
with the values of $`\nu `$ and $`\dot{M}`$ chosen appropriately, from our calculation, so as to enable the comparison.

If Fig. 2 we show the comparison of the analytical (effective) temperature distribution as given in (17) with the numerical result for the case of $`k=100`$, when the outer parts of the disc resemble most the analytical result (see Fig. 1). The lack of the inner, hottest, part in the numerical results is apparent. In these calculations the method described in §3.3 was used to monitor the energy dissipation during the simulations and calculate the discs temperature structures and their spectra. The spectral flux distributions ($`\lambda F_\lambda `$) in the three cases of $`k`$ depicted in Fig. 1 are displayed in Fig. 3.

It is apparent from the figure that with increasing $`k`$ the spectral distribution of the emitted radiation is shifted towards longer wavelengths. This trend is quite clear since as we have seen in Fig. 1 and 2, the hottest material, which resides close to the star, is significantly depleted by the magnetic interaction. In addition, as a result of the increase of the surface density in the outer cooler parts, the contribution to the long-wavelength part of the spectrum is more prominent. We stress once again that the temperature calculation and hence that of the spectrum suffers from the fact that the mass of the blobs is essentially unknown, so that the results displayed in Fig. 3 are of a qualitative nature and their value is mainly in uncovering the trend of the spectral behaviour. In addition, computer time limitations do not allow us to extend the calculation to the outer parts of the disc (here, in simulation 3 we have $`R_{\mathrm{out}}=10R_{\mathrm{co}}`$). Thus the spectral energy distributions in Fig. 3 fail to account for the cool material lying outside the outer limit of the computational zone and which should have an effect of increasing the strength at long wavelengths. Since the disc is more extended for large $`k`$ this effect should matter the most in the strongly magnetic case. The actual structure of the accretion disc for the three values of $`k`$ can be seen in Fig. 4, in which maps of the the surface density and temperature are shown.

From the examination of Fig. 3 and the above considerations it is thus reasonable to conclude that the magnetic interaction, at least as it is dealt with in this work, tends to flatten the spectrum for long wavelengths, consistently with the observed IR excess mentioned in the Introduction.

The results of all the other simulations (i.e. the ones with a different value of the viscosity scale parameter $`r_\mathrm{v}`$ are similar to the ones described above and thus the conclusions are similar. The actual value of $`r_\mathrm{v}`$ has only a little effect on the final (i.e. quasi-steady) surface density and temperature distribution, as long as this value is within the calculation feasibility range (see above). The difference is in the time that has to pass until the configuration settles down. This is understandable as the evolution time is essentially the viscous time scale and the latter decreases with increasing $`r_\mathrm{v}`$.

### 4.3 Effects of special magnetic field geometry

We have performed, in addition to the simulations described above in which the magnetic field was that of an aligned dipole, two similations with different magnetic configurations. In the first one we have used a magnetic dipole whose axis was tilted to the disc axis and in the second one a ”dipole slice” (see below) was introduced in order to estimate the effects of a localised magnetic field configuration.

#### 4.3.1 Tilted dipole

A single simulation was performed with $`r_\mathrm{v}=0.05;k=100`$, that is, similar to one of the simulations described above, but with the magnetic dipole axis tilted at an angle of $`45^{}`$ to the disc axis. The results show relatively few changes as compared to the aligned dipole case. There was almost no perturbation of the disc matter, relatively to the aligned dipole case, far out of the corotations radius. This is understandable since the blobs orbiting at these positions experience in this case a magnetic field whose direction oscillates with the period of the stellar rotation (which is considerably shorter than the blob’s Keplerian period). Thus, it is expected that the blobs actually react to the average magnetic field, i.e. behave in a similar way as in the aligned dipole case.

This expectation is wrong for blobs close to and inside corotation. Indeed, their motion should be channeled there by the magnetic field. This effect could not be observed in our simulations since the density of the particles in the magnetosphere, above the disc plane, is very small. A very high number of particles or alternatively an extremely strong magnetic field is needed in order to resolve, in our type of calculation, the particle flow along field lines. This remark obviously applies to the aligned dipole case as well.

#### 4.3.2 Dipole slice

To simulate the case in which the actual stellar magnetic field resembles a localized loop, similar to the case studied by PK, we have used a ”slice” of the aligned dipole, with Gaussian fall-off in strength in both the radial coordinate $`R`$ and the azimuthal angle $`\theta `$. Thus we have performed a simulation with $`r_\mathrm{v}=0.05;k=500`$ and with the magnetic field centered on the position $`R_0,\theta _0`$

$$B_{\mathrm{slice}}=B_{\mathrm{dipole}}\mathrm{exp}[(\theta \theta _0)^2/2d_\theta ^2]\mathrm{exp}[(RR_0)^2/2d_r^2],$$
(18)
where $`B_{\mathrm{dipole}}`$ is a dipole field (see the functional dependence in equation 6 and $`d_R`$, $`d_\theta `$ are the radial and angular Gaussian fall-off scales in the radial and azimuthal directions. In the simulation we have chosen $`R_0R_{\mathrm{co}}`$, consistently with the findings of KR that this should be the equilibrium configuration.

The results clearly depend on $`d_R`$ and $`d_\theta `$. When we chose $`d_RR_{\mathrm{co}}`$, the disc beyond $`R_{\mathrm{co}}`$ resembled the non-magnetic case ($`k=0`$), since the field is practically zero there. Closer to corotation and inside it, the density distribution looked like that of a full but weaker dipole, since the dipole slice field effectively acts on the orbiting blobs only along a fraction $`d_\theta /2\pi `$ of the blob’s orbit.

## 5 Summary and conclusions

The most prominent feature of our simulation results is the formation of a low density region (a hole) in the inner part of the disc as a result of the magnetic interaction. In addition, the density distribution becomes more spread out, as the slope of $`\mathrm{log}\mathrm{\Sigma }(R)`$ decreases with increasing $`k`$. This behaviour is typical and essentially qualitatively independent of all the other parameters. Although the ”hole” is the more noticeable feature (see e.g. the density maps in Fig. 4), it is less significant to overall observable features than the global change in the density distribution.

We have also seen, by comparing our calculations with what can be expected from analytical results, that we uncover some features (like the absence of hot material close to the inner hole) which can not be obtained from just classical disc models with a hole inside.

Due to computer power restrictions we were unable to achieve a high enough resolution, so as to see the accretion flow along the field lines and outflows from the system. In addition we were able to simulate only a limited portion of the disc and thus the calculations of the spectrum reveal only the general trend. This trend, resulting from the flattening of the surface density distribution, was always to shift emission power toward longer wavelengths. The flattening of $`\lambda F_\lambda `$ is most apparent near its maximum and it is reasonable that were it not for the close (computational) cutoff in $`R`$, the shape would remain flat to longer wavelengths. We propose therefore that the IR excess in T Tauri stars can be attributed to magnetic interaction, which modifies the functional dependence of the surface density in the surrounding discs.

Existing models attempting to explain the IR excess in T Tauri stars fall into two distinct categories: those invoking geometrical factors and others, proposing energy dissipation mechanisms operating preferably in the outer parts of the disc. Our model suggests a correlation between the spectrum flatness and the magnetic field strength, appears to be quite robust (practically any shape of the magnetic field would do) and does not require assumptions about flared shapes of discs or unusual energy dissipation modes. At this stage, the results of our calculation provide little more than support for a qualitative promising idea. As it was mentioned in the Introduction, it is possible to apply a similar prescription to the case in which the magnetic field penetrates the disc. We can reasonably expect that the results of such a calculation should not be too different, at least qualitatively, from the ones presented in this paper. It appears that all that is required is magnetic field lines imparting a torque on the gas, which changes in sign as we cross the corotation radius.

Significant results for both models of the magnetic interaction can be achieved in high resolution (significantly larger number of particles). However, to extend the idea into a reliable quantitative model, full multidimensional MHD simulations, including radiative transfer, have to be ultimately performed.

## Acknowledgments

We thank Giora Shaviv for helpful remarks and suggestions. Yigal Ultchin acknowledges the hospitality of the Astronomy Group at Leicester University and the financial support extended to him at the Technion. This research was supported by grants from the Israel Science Foundation, the Helen & Robert Asher Research Fund and by the Fund for the Promotion of Research at the Technion. We would also like to thank the referee, whose insightful comments helped considerably in improving this paper.


format: text

text:
# 1 ( a,b, taken from [] ). Mean logarithmic mass ⟨𝑙𝑛𝐴⟩ from the analysis of different sets of observables vs. 𝑙𝑔𝑁_𝜇^{𝑡𝑟} ( QGSJET prediction ). The two dashed lines indicate alternative predictions of the Single Source Model [], with the upper line being our latest and ’best estimate’. The sets displayed in (b) do not include the observable 𝑁_𝑒. It is seen that omitting 𝑁_𝑒 results in a heavier mass composition. (c) The ratio of mean values of ⟨𝑙𝑛𝐴⟩ from sets (b) and (a). We argue, in the text, that the ’true’ value of ⟨𝑙𝑛𝐴⟩ is a little below the crosses in Figure 1b.

PROPERTIES OF COSMIC RAY INTERACTIONS AT PEV ENERGIES

A.D.Erlykin 111Corresponding author: Erlykin A.D., e-mail: A.D.Erlykin@durham.ac.uk,2, A.W.Wolfendale 2

(1) Department of Physics, University of Durham, South Road, Durham DH1 3LE, UK

(2) P. N. Lebedev Physical Institute, Leninsky Prosp., Moscow 117924, Russia

## Abstract

An analysis has been made of the present situation with the high energy hadron-nucleus and nucleus-nucleus interaction models. As is already known there are inconsistencies in the interpretation of experimental data on the primary mass composition, which appear when different EAS components are used for the analyses, even for the same experiment. In the absence of obvious experimental defects, there is a cleari need for an improvement to the exising models; we argue that the most promising way is to introduce two effects which should be present in nucleus-nucleus collisions and have not been allowed for before. These are: a few percent energy transfer into the EAS electromagnetic component due to electron-positron pair production or electromagnetic radiation of quark-gluon plasma and a small slow-down of the cascading process in its initial stages associated with the extended lifetime of excited nuclear fragments. The latter process displaces the shower maximum deeper into the atmosphere.

1. Introduction

Among the most popular interaction models which are used for the analysis of experimental data in the PeV region and beyond are QGSJET, SYBILL, NEXUS, DPMJET, VENUS and HDPM, coupled with CORSIKA and partly with the AIRES simulation codes ( see  and references therein for a description of the models ). Early versions of these models differed significantly in their description of the development of the atmospheric cascades. One of the possible ways to evaluate the quality of a model is to derive the primary mass composition from the analysis of different observables: a good model should obviously give the same composition for any set of observables. We used this approach to interpret the data on the maximum development depth of the electromagnetic cascades $`X_{max}`$ and on the ratio of muon to electron sizes $`N_\mu /N_e`$ in terms of the mean logarithm of the primary mass $`lnA`$; the best consistency was obtained for the QGSJET model . A later analysis of the EAS hadron component confirmed this conclusion .

Since then the models have been improved and the differences between their predictions reduced . The precision of experimental data has been also improved and this offers the hope of further testing the models and of determining the actual value of $`lnA`$. Some of the new measurements in which Cherenkov light in the atmosphere is detected have given a deeper $`X_{max}`$ and hence a lighter mass composition than before . The extensive set of measurements on the ground and relevant simulations made by the KASCADE collaboration indicated that there is a difference in the estimate of $`lnA`$ derived from different sets of observables . At lower energies it was also difficult to get a consistent explanation of the measured muon and hadron trigger rates . All these facts provide evidence that further improvements to the models are urgently needed.

Meanwhile, a few calls for a quite radical change of the interaction model have appeared, inspired mostly by the existence of the well known knee in the primary cosmic ray ( CR ) spectrum at PeV energies. We discuss first the need for such a radical change and follow this with our own suggestion for modifications to the model for a better description of the modern experimental data ( i.e. KASCADE, DICE, BLANCA et al. ).

2. The possibility of a radical change in the interaction model

2.1. The possibilities and general remarks

A number of possibilities arise, and they will be considered in turn. Essentially they are of two types. The first considers that the primary spectrum ( represented in the usual way by a power law ) has an energy-independent exponent valid both below and above the PeV region where there is the knee in the spectrum inferred from the detected particles ’low’ in the atmosphere ( electrons, muons, hadrons ). The change in the interaction mechanism is then sudden ( at $``$ 3 PeV ) for CR-air nucleus interactions. The second is of the same philosophy i.e. that the original production spectrum has a unique exponent, but the interactions are in the interstellar medium ( ISM ) or even in the intergalactic medium ( IGM ) and not with air nuclei so that the spectrum of CR incident on the atmosphere already has the knee in it; there is no further radical change in the interaction mechanism for CR - air nucleus collisions.

The general argument against all these ideas is based on the observed anisotropy of the arrival directions of cosmic rays. The measured amplitude and phase of the anisotropy change in the vicinity of the knee and this gives support for its astrophysical origin , i.e. that the spectral shape arises because several sources, or types of source, contribute.

Of relevance, too, is the cosmological aspect, viz the effect of a proposed radical interaction change in the early Universe when the average energy of all ’particles’ was in the PeV region. As we have pointed out already  there are sound cosmological arguments against a radical change for P-P reactions in the PeV region. Some 10-14 sec after the Big Bang the average temperature of the primordial plasma was about 1016K and the mean energy per ’particle’ was of the order of a TeV in the centre of mass system, which corresponds to PeV energies in the laboratory system. If a radical change of the interaction occured at this energy changes to the subsequent evolution of the Universe would have taken place . Now the contemporary model ( i.e. without a radical change ) fits the observed data on the relative abundance of the light elements and it is difficult to imagine this being the case were there to have been a radical change of the interaction at 10-14 s.

Despite ’the case against’ the possibility of such a radical change it would have such important implications for particle physics that we make a detailed examination of some of the recent radical interaction models put forward to explain the knee.

2.2. A radical change for the CR - air nucleus interaction

Nikolsky has argued extensively that there is such a radical change. The latest idea  is close to that of his paper of 1995 . The author assumes that at a primary energy of about 3-5 PeV the multiplicity of secondary particles increases dramatically. The attenuation of the cascade speeds up and the intensity of showers in the lower half of the atmosphere falls, thereby forming the knee. There is actually no detailed model to be checked, because no numerical estimates have been given of the basic interaction parameters and their energy dependence, therefore it is difficult to assess the idea quantitatively.

One remark should be made, however. There is a difference between the previous and present pictures of the interaction advocated by the author. The 1995 model had the need to withdraw a part of the cascade energy and to put it into some unobservable components ( high energy muons, neutrinos etc. ) in order to form the knee. Now, in 2001, on the other hand, more energy is transferred to low energy particles of the hadronic cascade in the initial stages of its development. In our mind this is a backward step. The suggested increase of the energy transferred into the hadronic cascade in the atmosphere should be seen in the spectra of muons and Cherenkov light emitted by the charged particles in the atmosphere as an ’inverse knee or ankle’, i.e. as a decrease of the spectral slope. Nothing like that is actually observed. EAS muon size spectra at mountains: Tien-Shan , EAS-TOP  and at sea level: Ohya , KASCADE  and MSU  all show the ’ordinary’ knee, i.e. an increase of the spectral slope. To the existing 5 Cherenkov light spectra: TUNKA , HEGRA , DICE , BLANCA  and CACTI  has been added another one which is based on the combined use of Cherenkov light and EAS electron size measurement . All 6 demonstrate clearly the existence of the knee just like that in the particle size spectrum with a sharp increase of slope, which rules out the possibility of radical increasing the energy fraction deposited by the cascade in the atmosphere.

The possibility of having no knee in the primary spectrum and hiding the missing energy in high energy muons and neutrinos ( an apparently attractive idea ) has also been discussed by Petrukhin . The proposed model suggests the production at threshold of one or a few hypothetical heavy particles in each interaction, which then decay into muons and neutrinos. The energy flux contained in his muon spectra is less than the missing energy flux, needed to form the knee, by a factor of about 20 even for the case of the maximum flux of high energy muons. This means that the bulk of the missing energy should be carried away by neutrinos. The author succeeds in making the energy spectrum of muons from the hypothetical decay much harder than that of the background muons from $`\pi `$,$`K`$ decays and prompt muons. In this way it is possible to make the missing energy visible only as a very low flux of extremely high energy muons. However, such a model of the ’new physics’ ignores everything we know about the multiple production at high energies ( small cross-section for the production of heavy particles, conservation of lepton number in weak decays, etc.) and seems hardly credible.

Another model of ’new physics’ discussed by Kazanas and Nicolaidis , in spite of some theoretical arguments, is actually purely phenomenological. In it one opens a new channel for multiple production and puts part of the primary energy into some unobserved particles ( technihadrons, gravitons ), then, if one fits the energy dependence of the cross-section for the new channel and the energy fraction contained in the new particles, one can get a reasonable fit to the observed primary cosmic ray energy spectrum, including the knee. It is difficult to argue against the fit, but there seems to us to be no support for such speculation.

2.3. Radical changes in the interactions in the interstellar medium

It is, in principle, possible that the knee is due to interactions of some form beyond the atmosphere - in the ISM or IGM. If, as might be assumed by analogy with the ’X-particles’ produced by cosmic strings, the electromagnetic component is a prominent end-product of this interaction, then this hypothesis can be ruled out because the ensuing cascade in the Universe ( inverse compton scattering on the cosmic microwave background ) would produce a gamma ray intensity in excess of that measured ( see, e.g.  for the related cosmic string arguments ). However, there are other possibilities for the radical change of the interaction, which can now be considered.

The possibility of finite neutrino mass has given birth to models in which the knee is caused by the interaction of the primary protons with massive relic neutrinos. Assuming that the mass of the electron antineutrino is 0.4 eV/c2, Wigmans  converts some of the primary protons by inverse $`\beta `$-decay reaction into neutrons and positrons, starting from an energy of about 4 PeV. Indeed, protons will disappear from the primary flux, but we ourselves argue that even if the shower from a secondary positron misses the detector, the shower produced by the neutron will be practically indistinguishable from a proton-initiated shower, because neutrons preserve more than 99% of the primary proton energy; thus EAS measurements would not see a knee. Taking into account the problem of the very small cross-section for the neutrino induced reactions, which requires us to assume a very dense halo of neutrinos gravitationally attracted to the cosmic ray source region ( be it the Galaxy or source within the Galaxy ), one should remark that the PeV neutron will decay again into a proton and a positron within tens of parsecs and the model loses its value.

Another model of the intermediate type was proposed by Ehrlich , He used some theoretical ideas that the electron neutrino is a tachyon, which in fact does not contradict the existence of a finite neutrino mass. In this case, starting from the same energy of 4 PeV, the primary proton has the possibility to decay into a neutron, a positron and an electron neutrino. This decay might create a narrow peak of neutrons in the primary cosmic ray flux in the region close to 4.5 PeV, which can perhaps imitate the knee. All the scepticism put forward in the previous paragraph is applicable also in this case, however.

The common deficiency of all the models on the radical change of the interaction is that no attempts have been made to reconcile it with the multitude of other experimental data, different from those used for the motivation of the model. The conclusion which can be drawn from this survey is that at the present time there is neither a reasonable model nor solid argument which justifies a radical change of the interaction. We conclude that such a change as is necessary for the interaction model should not be radical and not responsible for the origin of the knee, but only responsible for restoring consistency for all cosmic ray measurables.

3. What kind of the modification of the interaction model do we need ?

The variety of observed EAS components and precise measurements of their characteristics by the KASCADE experiment allowed this collaboration to derive the primary mass composition using the multivariate analysis of data. Remarkably, as has been pointed out already, it was found that not only do different methods applied to the same set of observables give different results, but that there is a systematic difference between the results obtained using the same method, but different components. In Figure 1 the mean logarithm of primary mass $`lnA`$ derived from different sets of observables is shown vs. the so called truncated muon number $`N_\mu ^{tr}`$ ( the figure is taken from  ). $`N_\mu ^{tr}`$ is the number of low energy muons ( $`>0.25GeV`$ ) contained between 40 and 200 m from the shower axis and is an observable adopted by the KASCADE collaboration as a measure of the primary energy, independent of the primary mass. They also use the number of higher energy muons ( $`>2.4GeV`$ ) $`N_\mu ^{}`$ as well as some observables related to the hadron component for the analysis of the mass composition . The fundamental problem is, of course, that all the $`lnA`$ values should be the same, at the same value of $`N_\mu ^{tr}`$, and they are not.

A number of features are of interest and, presumably, of some importance.

(i) The dispersion of the individual $`lnA`$ values about smooth lines through them indicates that non-systematic errors are $`\pm `$0.05.

(ii) The use of the electron size $`N_e`$ results in an average lighter composition ( Fig.1a ). On the other hand, omitting electrons and using just muons and hadrons results in a heavier composition ( Fig.1b ).

(iii) Considering the mean values from the two sets of data ( filled-in circles ), the ratio of $`lnA`$ for (b) to (a) increases smoothly with $`lgN_\mu ^{tr}`$, i.e. the discrepancy between (b) and (a) rises with increasing primary energy ( Fig.1c ).

(iv) Within each of (a) and (b) the differences between the $`lnA`$ values using different combinations of data ( $`N_e,N_\mu ^{}`$ etc. ) are systematic and real ( see (i) ).

(v) In both (a) and (b), the highest estimated $`lnA`$ values are for the cases where $`N_\mu ^{}`$ is absent - i.e. not used in the analysis. The lowest values of $`lnA`$ are where $`N_\mu ^{}`$ has the greatest weight in the analysis, i.e. the number of parameters ( 2 in 1a and 3 in 1b ) is a minimum. It is an indication that the higher energy muons ( $`>2.4GeV`$ ) in the observed showers are in the greater deficit compared with model predictions .

(vi) With respect to (v), taking only the two sets ( one in 1a and one in 1b ) which did not use $`N_\mu ^{}`$, it is interesting to note that the former points are close to our lower variant ( lower dashed curve ) in 1a and close to our ’best estimate’ ( upper dashed curve ) in 1b. The differences are, therefore, significant.

The difference in the $`lnA`$ values derived using different observables ( Fig.1 ) points to inadequacies in the models used for the analysis of the experimental data and provides an impetus to correct them. However, as remarked already, it is sure that the corrections should not be radical, because the difference between the $`lnA`$ is not large, typically $`\delta lnA0.20.4`$ or, in terms of the ratio of the mean $`lnA`$ values from Figures 1b and 1a, an increase from 1.12 to 1.37. Below, we show that the mentioned systematic difference between the $`lnA`$ values can be an indication that the energy distribution between the different shower components is slightly different from that in the models: specifically, the actual mean number of electrons $`N_e`$ in EAS appears to be higher, that of muons $`N_\mu `$ slightly lower and that of hadrons $`N_h`$ lower still than in the models.

It is well known from the present models that for the same primary energy the number of muons in nuclei-induced showers is higher than in proton-induced ones. On the other hand the number of electrons and hadrons in nuclei-induced showers observed in the lower half of the atmosphere is lower than in proton ones . If one has an opportunity to measure the primary energy, by muons, Cherenkov light or another technique, and finds that the shower has a low $`N_\mu `$ or a high $`N_e`$ ( actually the ratio $`\frac{N_\mu }{N_e}`$ is important ) the conclusion will be that this shower is initiated by a proton or light nucleus. On the contrary, if one finds that the shower has a low $`N_h`$ the conclusion will be the opposite, i.e. that the shower is initiated by a heavy nucleus. This is exactly what is observed in the showers at sea level. In order to get a better consistency between the $`lnA`$ values obtained using $`N_e`$, $`N_\mu ^{tr}`$ and $`N_h`$ ( Fig.1 ) the model should have higher $`N_e/N_\mu ^{tr}`$ and smaller $`N_h/N_\mu ^{tr}`$ ratios. This situation could be achieved if the corrections give a reduction of $`N_h`$, a smaller reduction of $`N_\mu ^{tr}`$ and an increase in $`N_e`$.

The same conclusion on the improvement of the model can be drawn from an analysis of KASCADE event rates . Both the muon and hadron trigger rates, observed by KASCADE, are lower than expected from the model calculations. This discrepancy indicates again that the actual numbers of muons and hadrons in EAS are lower than in the models, although the energy region responsible for the trigger rates is lower than that analysed for the mass composition around the knee. The argument that the lower trigger rates might be due to a lower primary intensity is disproved by the fact that the ratio of hadron to muon trigger rates, which is independent of the absolute intensity, is also lower than in the calculations . To get agreement the predicted muon trigger rate should be reduced by a factor 0.89 $`\pm `$ 0.06, and the hadron rate by 0.54 $`\pm `$ 0.08 compared with the rates given by the QGSJET model. This analysis is of value because it indicates that the needed reduction of the number of muons in the model should be about 6%, for hadrons it is bigger - about 29%. Because muons and hadrons are the products of hadronic cascades, it is evident from the energy balance that to reduce the energy contained in the hadronic cascade one has to increase the energy transferred to the electromagnetic cascade. At this stage it can be remarked that an increase in this energy is predicted by theory ( although it appears not to have been included into the models hitherto ).

4. Numerical estimates

4.1. The balance between the EAS components

To evaluate the effect of the proposed change of the balance between different cascade components we applied the semi-quantitative analytical approach. We assumed that there are only nucleons (N), pions ($`\pi `$), muons+neutrinos ($`\mu \nu `$) and the electromagnetic component ($`e\gamma `$) in the cascade. The total and partial inelasticity for nucleons: $`K_{tot}^N`$ and $`K_\gamma ^N`$ and for pions $`K_{tot}^\pi =1`$ and $`K_\gamma ^\pi `$, interaction mean free paths for nucleons $`\lambda _N`$ and for pions $`\lambda _\pi `$ were taken to be energy independent, i.e. constant. The system of kinetic equations, describing the longitudinal development of the nucleon, pion, muon+neutrino and the electromagnetic energy of the cascade has been taken as
$`\frac{dE_N}{dX}=\frac{E_N}{\mathrm{\Lambda }_N}`$
$`\frac{dE_\pi }{dX}=(\frac{1}{\mathrm{\Lambda }_\pi }+\frac{E_{cr}}{E_\pi X})E_\pi +\frac{K_{tot}^NK_\gamma ^N}{\lambda _N}E_N`$
$`\frac{dE_{\mu \nu }}{dX}=\frac{E_{cr}}{E_\pi X}E_\pi `$
$`\frac{dE_{e\gamma }}{dX}=\frac{E_\pi }{\mathrm{\Lambda }_\pi }+\frac{K_\gamma ^N}{\lambda _N}E_N`$
Here X is the atmospheric depth, $`\mathrm{\Lambda }_N=\frac{\lambda _N}{K_{tot}^N}`$, $`\mathrm{\Lambda }_\pi =\frac{\lambda _\pi }{K_\gamma ^\pi }`$, $`E_{cr}`$ \- the critical energy for pion decay in the air, $`E_\pi `$ \- the mean energy for pions. The boundary condition at $`X=0`$ is $`E_\pi +E_{\mu \nu }+E_{e\gamma }=0`$
$`E_N=E_0`$
where $`E_0`$ is the total energy of the cascade. The solution of the system of equations is
$`E_N=E_0exp(\frac{X}{\mathrm{\Lambda }_N})`$
$`E_\pi =E_0\frac{K_{tot}^NK_\gamma ^N}{\lambda _N}exp(\frac{X}{\mathrm{\Lambda }_\pi })_0^X(\frac{y}{X})^{\frac{E_{cr}}{E_\pi }}exp(\frac{y}{\mathrm{\Lambda }_{\pi N}})𝑑y`$
$`E_{\mu \nu }=_0^X\frac{E_{cr}E_\pi (y)}{E_\pi y}𝑑y`$
$`E_{e\gamma }=_0^X\frac{E_\pi (y)}{\mathrm{\Lambda }_\pi }𝑑y+\frac{E_0K_\gamma ^N}{K_{tot}^N}(1exp(\frac{X}{\mathrm{\Lambda }_N}))`$
Here $`\mathrm{\Lambda }_{\pi N}=\frac{\mathrm{\Lambda }_\pi \mathrm{\Lambda }_N}{\mathrm{\Lambda }_\pi \mathrm{\Lambda }_N}`$.

We calculated also the longitudinal development of the electron size of the shower $`N_e`$ as
$`N_e(X)=_0^X\frac{dE_{e\gamma }}{dy}\stackrel{~}{N}_e(Xy)𝑑y`$
The longitudinal profile of the electromagnetic cascade has been taken as
$`\stackrel{~}{N}_e(t)=\frac{E_{e\gamma }}{1GeV}exp(tt_{max}2tlnS)`$, where the atmospheric depth $`y`$ is measured in units of the radiation length $`t_0`$: $`t=\frac{y}{t_0}`$, $`t_{max}=1.7+0.76ln(\frac{E_{e\gamma }}{\beta })`$, $`S=\frac{2t}{t+t_{max}}`$ . Here $`E_{e\gamma }`$ is the total energy of the electromagnetic cascade, $`t_0`$ = 37.1 gcm-2 and $`\beta `$ is the critical energy for the electromagnetic cascade in air, equal to 0.081 GeV. Although this approximation is recommended for nuclear cascades, its use for the electromagnetic cascades is more appropriate for our energy-balance approach, since it gives the correct estimate of the depth of the shower maximum ( $`X_{max}610gcm^2`$ ) and its elongation rate ( $`ER65gcm^2`$ ).

The integrals for $`E_\pi `$, $`E_{\mu \nu }`$ and $`E_{e\gamma }`$ were calculated numerically. The basic parameters for the calculation were taken as $`K_{tot}^N`$=0.60, $`K_\gamma ^N`$=0.20, $`K_\gamma ^\pi `$=0.33, $`\lambda _N`$=90 gcm-2, $`\lambda _\pi `$=120 gcm-2, $`E_{cr}`$ = 90 GeV, $`E_\pi `$ = 300 GeV.

The results of the calculation are shown in Figure 2 by the full line. According to our suggestion, we increased the energy fraction transferred by nucleons into the electromagnetic component, leaving all the other characteristics untouched. We increased $`K_\gamma ^N`$ from 0.20 to 0.26 ( the reason for this particular value will be clear later ). The result is shown in Figure 2 by the dashed line: the muon energy at sea level ( 1000 gcm-2 ) decreased by $``$6%, the hadron energy decreased by $``$23%, the energy transferred into the electromagnetic component increased by $``$2%.

4.2. The triangle diagrams

The balance of the energy contained in the major components of the shower is convenient to analyse using the so-called “triangle diagrams” . If the height of an equilateral triangle is equal to 1, then for each point inside this triangle the sum of the distances to its sides is equal to 1. If we know the energy fractions carried by the electromagnetic ($`\delta _{e\gamma }`$), muon ($`\delta _\mu `$) and hadron ($`\delta _h`$) components of the shower at the observation level, so that $`\delta _{e\gamma }+\delta _\mu +\delta _h=1`$, then each shower can be presented by a single point inside the triangle. Our basic shower ( Figure 2, full line ) is shown in Figure 3 by a full circle ( $`\delta _{e\gamma }`$=0.394, $`\delta _\mu `$=0.505, $`\delta _h`$=0.101 ). The desired direction for the shift of the energy balance in the modified model is shown by the straight arrow in Figure 3b.

However, despite the increase of the energy transferred into the electromagnetic component, the preserved electromagnetic energy and the electron size of the shower $`N_e`$ at sea level decreased by 14% due to its faster development and then faster attenuation of the cascade. The point in the triangle diagram moved in the different direction ( open circle ). The shift is similar to the case when we increase the total inelasticity $`K_{tot}^N`$ and $`K_\gamma ^N`$ with $`K_\gamma ^N=\frac{1}{3}K_{tot}^N`$.

If we were to assume a larger fluctuation $`\sigma (K_\gamma ^N)`$ compared with those in the model this would result in a wider distributions of $`N_e`$, $`N_\mu `$ and $`N_h`$ and in a bigger value of $`\sigma (lnA)`$. Within our approach we cannot give the numerical estimates, but we expect that the mean energy balance would not move in the desired direction and the mean value $`lnA`$ would not change significantly. Therefore, an increase of the mean value or increased fluctuations of $`K_\gamma ^N`$ cannot give the required result.

4.3. Data from mountain laboratories

We should remark that the conclusion about the decrease of $`N_e`$ relates only to measurements at sea level. At mountain altitudes the increase of $`K_\gamma ^N`$ results in an increase in $`N_e`$ as seen in Figure 2, however it is not enough here, either. The best complex EAS array at mountain altitude was at Tien-Shan and the data there also indicated the need for the model to be improved . The problem was with the so-called ’observed’ mass composition, i.e. the mass composition of primary particles responsible for showers at the observation level selected by one of their parameters. The Tien-Shan array is at a depth of 690 gcm-2, i.e. PeV-showers are definitely observed beyond their maximum. The $`N_e`$-selected showers have to be enriched by showers initiated by protons and light nuclei, while in contrast the $`N_\mu `$-selected showers have to be enriched by heavy nuclei initiated showers. However, in fact, the ’observed’ composition at Tien-Shan was about the same for these two selections . The consistency can be restored by the same way as in KASCADE, i.e. by changing the balance between the muon and electromagnetic components in the models, in favor of the latter.

4.4. Slowing down the cascade

We have concluded that the mere increase of the energy fraction $`K_\gamma ^N`$ transferred into the electromagnetic component does not change the balance properly, although it does help in some respects and, as we have stated, it must be present. We now argue that another effect that should be present will bring about the desired effect: the slowing down of the development of the cascade in its initial stages. The physics behind this modification will be discussed later. Here we just demonstrate the validity of this argument in a semi-qualitative way. For illustration purposes we slow down the development of the hadronic and electromagnetic cascade by increasing the elongation rate $`ER`$ from 65 gcm-2 to 71 gcm-2, preserving $`K_\gamma ^N`$ = 0.26. The result is shown in Figure 2 by dotted lines and in Figure 3 by an open triangle. The direction of the shift is correct and its magnitude can be adjusted to achieve consistency in the mass composition.

Thus we conclude that the increase of the energy transferred into the electromagnetic component ( $`K_\gamma ^N`$ ) combined with the slowing down of the development of cascades in their initial stages is the most realistic way to improve the particle interaction model and to achieve a consistent estimate of the primary mass composition.

5. Theoretical arguments

5.1. General remarks

All the arguments hitherto were purely phenomenological. However, there are theoretical arguments too which lend support to the phenomenological consideration as will be demonstrated. At the moment nearly all the models for hadron-nucleus (hA) interactions implemented in CORSIKA are well developed and theoretically justified. Nucleus-nucleus (AA) interactions are considered within the framework of the Glauber approach, which reduces the $`AA`$-collision to the sum of a few $`hA`$-collisions. No new processes are assumed to occur in $`AA`$-collisions. However, there are processes which must be present but which have not been included hitherto, and we think that both effects: the increase of $`K_\gamma `$ and the slowing down of the cascade development refer to processes, intrinsic for $`AA`$-collisions, which will cause at least near-equality of the $`lnA`$ values in Figures 1a and 1b.

5.2. The increase of $`K_\gamma `$

At PeV energies, when the Lorenz-factor of the projectile nucleus approaches $`10^6`$ and its charge $`Z`$ is $`1030`$, the density of virtual photons in the contracted coulomb field is so high that when the projectile collides with a nitrogen ( $`Z`$ = 7 ) nucleus of the air and the impact parameter is small ( central collision ) the probability of electron-positron pair production including multiple pair production becomes quite appreciable . Using the numerical estimates made by Erlykin  for p-air collisions, we can extrapolate them to Fe-air collisions and expect an increase of $`K_\gamma ^{Feair}`$ for the average collision at PeV energies by 0.05-0.13. For specific central collisions the effect will be even higher. Our 0.06 of §4.1 was thus rather conservative. We also emphasize the fact that ordinary $`K_\gamma ^{AA}`$ value in AA-collisions is not very high. Typically for Fe-air collisions the average number of wounded nucleons $`n_w`$ is about 7-8 and $`K_\gamma ^{Feair}K_\gamma ^{PA}\frac{n_w}{A_{Fe}}0.03`$. Therefore even a few percent increase of $`K_\gamma ^{AA}`$ in the first interaction can be essential for the subsequent development of the atmospheric cascade.

Experimental results confirm the intense production of electron-positron pairs in $`AA`$ collisions. In the CERES/NA45 experiment the number of $`e^+e^{}`$ pairs in the low invariant mass region $`m_{e^+e^{}}<200MeV/c^2`$ exceeds the number in the higher mass region of $`m_{e^+e^{}}>200MeV/c^2`$ by a factor of $`10^4`$. For the CERES/NA45 collaboration it is just a background, because they are interested in dileptons with high invariant mass  as a signature of an elusive quark-gluon plasma. We believe that to check our assumption the study of the dilepton production should be concentrated preferentially in the region of invariant mass below $`10MeV/c^2`$.

An additional energy transfer into an electromagnetic component can arise also from an excess of direct photons, which has been predicted theoretically as a signature of the quark-gluon plasma  and is now observed in AA-collisions both at low and high transverse momenta .

Another theoretical possibility is to increase the transfer of the energy into the electromagnetic component, by way of $`hA`$-interactions; this could be connected with the breaking of the diquark and its subsequent recombination into 3 $`\pi ^0`$-mesons .

5.3. Slowing down of the cascade development

The theoretical motivation of the second assumption is also connected with the properties of AA-interactions, but relates mainly to peripheral collisions. It is known that even at high energies a projectile nucleus does not fully disintegrate into constituent nucleons in the first interaction but fragment into few pieces of different mass . Some of these fragments are excited and, after the de-excitation, if it occurs in space, give rise to MeV gamma-ray lines, observed from ’discrete’ sources and the interstellar medium . The lifetime of the excited fragments varies from a ’nuclear’ time $`10^{23}`$ sec to millions of years. For AA-interaction at PeV energies both the lifetime before de-excitation and the energy of emitted gamma-quanta are extended by the factor of $`10^510^6`$ due to relativistic effects. As a consequence one can expect an additional sub-PeV electromagnetic cascade to be initiated a few hundred meters below the point of the first interaction. This effect will slow down the development of the electromagnetic cascade and shift its maximum.

In general the shift is limited to no more than 1-2 mean free paths for the inelastic interaction of the fragments. However, there might be an additional effect which is, it must be admitted, speculative. It is based on the so-called “sling” effect . and is connected with the deformation of the shape of the nuclear fragment, which rotates with a high spin . As a result of such deformation one can expect an increase in the mean free path for inelastic interaction which enables the excited nuclear fragments to penetrate deeper into the atmosphere, slowing down the cascade development and resulting in a further shift of $`X_{max}`$. There will be an additional increase of $`N_e`$ for nucleus-induced cascades like that shown in Figure 2b by the dotted line.

All the consequences of the ’sling’ effect: polarisation of the secondary nuclear fragments, deformed shape of the nuclei rotating with a high spin and cross-section fluctuations of the excited nuclei are true effects observed at low energies . The problem is whether they still hold at higher energies. However, even if they don’t, the de-exitation of the nuclear fragments itself will give the needed shift of $`X_{max}`$ if their lifetime is not less than $`10^{12}`$ sec.

It is difficult to make reliable numerical estimates of the shift since there is not enough experimental data and no good theory of the effect. We think that to check our assumption the study of the properties of the secondary nuclear fragments should be made in high energy AA-collisions with different nuclear beams.

6. Discussion

6.1. Application of the new model to PeV energies

6.1.1. The mass composition problem

In order to get better consistency of the results on the primary cosmic ray mass composition, the combination of two effects is needed. The increase of only $`K_\gamma ^{AA}`$ is not enough because although it results in a decrease of $`N_\mu `$ and $`N_h`$ it also leads to a fall in $`N_e`$ if the measurements are made at sea level, whereas we need an increase in $`N_e`$. The slowing down of the electromagnetic cascade alone does not change the ratio between muons and hadrons. The slowing down of just the nuclear and electromagnetic cascades ( e.g. by increasing the mean free path for nucleon interaction $`\lambda _N`$ without increasing $`K_\gamma ^N`$ ) results in an increase in both $`N_e`$ and $`N_h`$. Only the combination of both effects: the increase of $`K_\gamma ^N`$ and the slowing down of the development of the electromagnetic cascade changes the energy balance in the needed direction: $`N_\mu `$ and $`N_h`$ decreased, $`N_e`$ and $`X_{max}`$ increased ( Fig.3 ). The required changes of $`K_\gamma ^{AA}`$ and $`X_{max}`$ are small, because the sensitivity of the energy balance to these changes is rather high.

The result should be a coming together of the $`lnA`$ values in Figures 1a and 1b. There will be an increase in the values in Figure 1a and a small reduction in those in Figure 1b. In view of the problems associated with $`N_\mu ^{}`$ ( see §3 ) most of the predictions in Figure 1b should be disregarded; the ’best estimate’ is then just a little below the crosses in Figure 1b - these values are reassuringly close to our latest estimate indicated by the upper dashed curve in Figure 1.

6.1.2. The $`X_{max}`$ controversy

Although the $`X_{max}`$ aspect for EAS is one of some complexity we discuss it briefly. The assumption of a slowing down of the cascade development and the consequent shift of $`X_{max}`$ in $`AA`$-interactions helps us to understand the striking difference between the mass composition derived recently by Swordy and Kieda  and Fowler et al. , using Cherenkov light measurements and those by Roth et al.  and others, without. Those with Cherenkov light give a significantly lighter mass composition and lower intensity of the primary energy spectrum than those obtained without the use of $`X_{max}`$ . If there is indeed a shift of $`X_{max}`$ and a consequent increase of $`N_e`$ the shower looks more like that initiated by proton and its primary energy, based on the on-ground measurements is underestimated. Insofar as the $`N_e`$ values are most in doubt, from the prediction standpoint the real mass composition should be closer to that derived by using direct measurements in the stratosphere or using EAS muons and hadrons on the ground ( Fig.1b ).

6.2. Application of the new model to higher energies

If the increase of $`K_\gamma ^{AA}`$ and the slowing down of the development of the cascade are due to electron-positron pair production, and de-excitation of nuclear fragments, there will be significant further consequencies of the model because of its charge and energy dependent QED effect and the displacement of the depth of shower maximum. Both these effects will grow with increasing collision energy. In this case it is not surprising that the inconsistency in the estimates of the primary mass composition grows with energy, as noticed by Roth et al.  and indicated in Figure 1c.

Of greatest importance is the mass composition at the highest energies, where, due to the low anisotropy and the rapid growth of $`X_{max}`$ with energy, particles are generally considered to be extragalactic (EG). There is controversy, here, between those who subscribe to the common view that the EG particles are protons ( eg  and ourselves ( eg , following much earlier work by Tkaczyk et al. ). A later publication will deal with this aspect in detail but a few brief remarks here are in order.

(i) The change of $`K_\gamma ^{AA}`$ will be even bigger at the highest energies ( say above 1010 GeV ) and the result will be an increase in the mean primary mass: this follows from the analysis of Ave et.al., in which inclined muon showers were detected - the corresponding primary energy will now be higher and the rates more in accord with expectation for a significant fraction of primary iron nuclei, a recent recalibration of Haverah Park energies  leads to a reduction in primary energy and an even further increase in the iron fraction.

(ii) The displacement of the depth of maximum to greater values will, again, lead to an increase in the mean primary mass as indicated by the $`X_{max}`$ values from the experimental data ( see Fig.1b ).

7. Conclusions

Our analysis of the present situation with the high energy interaction models indicates that there is no support for the introduction of a radical change of the models, sometimes proposed to explain the origin of the knee in the primary cosmic ray energy spectrum, However, the inconsistencies in the interpretation of the experimental data on the primary mass composition, obtained when different EAS components are used for the analysis, indicate the need for some improvement to the models which were used hitherto. We propose that the most promising way is to introduce an additional ( a few percent ) energy transfer into the EAS electromagnetic component combined with a slowing down of the cascade development at its initial stages, which is followed by a small ( $`2030gcm^2`$ ) shift of $`X_{max}`$ into the deeper atmosphere and the consequent increase of $`N_e`$. The most likely processes which can be responsible for such changes ( e.g. electron-positron pair production including multiple pairs, direct photons from the hypothetical quark-gluon plasma, excitation of the secondary nuclear fragments ) are those which occur in nucleus-nucleus collisions and they should indeed be present at some level. The importance of these processes is expected to grow with energy and offer the hope of resolving some controversies at very high energies.

Acknowledgements

The authors are greatful to the UK’s Particle Physics and Astronomy Research Council and to The Royal Society for financial support. Also we thank Professors G.Schatz, H.Rebel, K.H.Kampert, J.Kempa and R.Chapman for useful discussions.