format: text

text:
# The jet-disk symbiosis model for Gamma Ray Bursts: cosmic ray and neutrino background contribution

## 1 Introduction

More than 30 years after their discovery, thanks to the Burst and Transient Source Experiment (BATSE) and the Italian-Dutch satellite BeppoSAX, the scientific community now knows that Gamma Ray Bursts (GRBs) are isotropically distributed in the sky (Fishman $`\&`$ Meegan 1995) and that at least some of them are at cosmological distances (GRB970228: Djorgovski et al. 1999b, GRB970508: Metzger et al. 1997, GRB971214: Kulkarni et al. 1998, GRB980613: Djorgovski et al. 1999a, GRB980703: Djorgovski et al. 1998, GRB990123: Hjorth et al. 1999, GRB990510: Vreeswijk et al. 1999, GRB990712: Galama et al. 1999). But the present data available for redshift position and host galaxy localization are still too few to give us good statistics to study the evolution of GRBs and their redshift distribution. Before the discovery of GRB afterglows by BeppoSAX, the only way to study their distributions was to compare some GRB properties (like for example the intensity), with some parametric models (Fenimore $`\&`$ Bloom 1995, Cohen $`\&`$ Piran 1995, Kommers et al. 1999). Because of this lack of information, it is still necessary to assume that GRBs follow the statistical distribution of some other better known objects to obtain the GRBs fluence or flux distribution itself.

The origin of GRBs is still controversial. According to different models, their progenitor can be identified with the merging of two neutron stars, or with the collapse of a massive star. In the model presented by Pugliese et al. (1999), GRBs are created inside a pre-existing jet in a binary system formed by a neutron star and an O/B/WR companion, where the input energy comes from the collapse of the neutron star into a black hole and the emission is due to synchrotron radiation from the ultrarelativistic shock waves that propagate along the jet with a low-energy cut-off in the electron distribution. Following this scenario, the birth of GRBs cannot happen too far from the region where the progenitor formed, and this implies that their rate should be connected with the Star Formation Rate (SFR). Already other authors studied the connection between the SFR and GRBs flux distribution. For example, Wijers et al. (1998) showed that the assumption that the GRB rate is proportional to the SFR in the universe is consistent with the GRB flux distribution.

In Sect. 2 we calculate the cumulative distribution of GRB fluences using two SFR distributions as a function of redshift, the one by Miyaji (Miyaji et al. 1998), and the other by Madau (Madau et al. 1996). We compared it with the data from the BATSE catalogue. In Sect. 3 we calculate the maximum energy available in our model to obtain high energy cosmic rays. In Sect. 4 we present our results for the contribution of GRBs to the cosmic rays distribution, both Galactic and extragalactic and in Sect. 5 the eventual contribution from GRBs to the neutrino flux.

## 2 The rate of GRBs

The aim of this section is the calculation of the number of GRBs per year and per 100 $`\mathrm{Mpc}^3`$ as a function of redshift, assuming beamed emission in a jet, (e.g. Pugliese et al., 1999). The following calculations are quite general and we emphasize that their validity does not depend on the particular model used.

To obtain a result as close as possible to the data, we chose the fluence as the quantity that can well represent the characteristics of GRBs. We follow Petrosian $`\&`$ Lloyd (1997), who showed that the fluence is the most appropriate parameter to study the cosmological evolution of GRBs. The relation between the redshift $`z`$ and the fluence $`f`$ is given by (Petrosian $`\&`$ Lloyd 1997) as a function of the specific luminosity L expressed in $`[\mathrm{erg}\mathrm{s}^1\mathrm{Hz}^1]`$:

$$f=\frac{L\mathrm{\Delta }t\mathrm{\Delta }\nu }{4\pi d_\mathrm{L}^2(1+z)^{\alpha 3}}\mathrm{erg}/\mathrm{cm}^2.$$
(1)
Here $`\alpha `$ is the photon flux spectral index, $`d_\mathrm{L}=(2c/h)(1+z\sqrt{1+z})`$ is the cosmological luminosity distance (Weinberg, 1972), and $`H_0=h(100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1)`$ is the Hubble constant. In the transformation from the emitter frame to the observer frame the two corresponding redshift contributions from the frequency and the time dependence are cancelled.

In our model (Pugliese et al. 1999), we used $`\alpha =2`$, and the corresponding value for the fluence is:

$$f=\frac{L\mathrm{\Delta }t\mathrm{\Delta }\nu }{4\pi d_\mathrm{L}^2}(1+z)\mathrm{erg}/\mathrm{cm}^2.$$
(2)
In this way the redshift as a function of the fluence is:

$$1+z(f,f_{})=\left(\frac{1}{2}\sqrt{\frac{f_{}}{f}}+1\right)^2,$$
(3)
where $`f_{}=\frac{L(\nu _1)\mathrm{\Delta }t\mathrm{\Delta }\nu }{4\pi ^2}\frac{h^2}{c^2}`$ $`\mathrm{erg}/\mathrm{cm}^2`$ is the reference fluence and $`L(\nu _1)`$ is expressed in $`[\mathrm{erg}\mathrm{s}^1\mathrm{Hz}^1]`$.

At this point of our calculations it is important to define the role of the parameter $`f_{}`$. In fact a relevant question is whether GRBs are standard candles (i.e. $`f_{}`$ is constant) or whether they are distributed with a Luminosity Function (LF) with, e.g., a power law in $`f_{}`$. Developing our calculations, we arrived at the same results found by Kommers et al. (1999). Here the authors showed that the best model is the one in which the star formation rate is combined with the luminosity function distribution. In agreement with them, also in our model GRBs cannot be considered as standard candles. We use a distribution for the luminosity function with a power law index $`\beta `$, and the corresponding law for the parameter $`f_{}`$ is given by:

$$\varphi (f_{})df_{}=\varphi _0f_{}^\beta df_{},$$
(4)
where $`\varphi _0`$ is the normalization parameter equal to $`1/(f_\mathrm{b}^1f_\mathrm{a}^1)`$. We adopt $`f_\mathrm{a}=10^8\mathrm{erg}/\mathrm{cm}^2`$, and $`f_\mathrm{b}=2.2\times 10^4\mathrm{erg}/\mathrm{cm}^2`$.

Here we do not yet answer the question what the physics of this luminosity function may be. It is plausible in the context of our model, that it is directly connected to the mass flow in the pre-existing jet prior to the GRB explosion. If this were the correct interpretation, the mass accretion rate, and correspondingly, the mass flow rate in the jet may follow a power law distribution, a point which we will pursue elsewhere.

We also calculate the GRB rate using two different star formation rates and compare the corresponding results with the data.

### 2.1 The number count

The number count of GRBs sources is given by the expression $`dN(z)=F(z)(dt/dz)dVdz`$, where $`F(z)`$ is equal to the product of the SFR $`\psi (z)`$ and the luminosity function $`\varphi (z)`$. Obviously the luminosity function may change with redshift $`z`$, but for simplicity and as a first step, we use this ansatz. In term of the fluence $`f`$ of GRBs in the gamma ray band, and of the reference fluence $`f_{}`$, the number count is:

$`dN(f,f_{})`$ $`=`$ $`[\psi (f,f_{})]\left[{\displaystyle \frac{dt}{dz}}(f,f_{})\right]\left[4\pi d_\mathrm{L}^2(f,f_{}){\displaystyle \frac{dd_\mathrm{L}}{dz}}(f,f_{})\right]`$ (5)
$`\times [dz(f,f_{})],`$
where:

The first term $`\psi (z(f,f_{}))`$ represents the star formation rate as a function of the redshift.

The second term represents the temporal interval in which the rate is calculated. It is given by:

$$\frac{dt}{dz}(f,f_{})=\frac{1}{h}\frac{1}{(\frac{1}{2}\sqrt{f_{}/f}+1)^5}.$$
(6)
The third term is the interval of volume as a function of the fluence. It corresponds to the value:

$$4\pi d_\mathrm{L}^2\frac{d}{dz}d_\mathrm{L}=4\pi \frac{c^3}{h^3}\frac{f_{}}{f}\left(\frac{1}{2}\sqrt{\frac{f_{}}{f}}+1\right)\left(\sqrt{\frac{f_{}}{f}}+1\right).$$
(7)
The fourth term represents the redshift interval and it is given by:

$$dz(f,f_{})=\frac{1}{2\pi }\frac{\sqrt{f}_{}}{(\sqrt{f})^3}\left(\frac{1}{2}\sqrt{\frac{f_{}}{f}}+1\right)df.$$
(8)
The integration over the parameter $`f_{}`$ will be shown in the Eq. 11.

### 2.2 The GRB rate distribution using the Miyaji SFR

Here we use the star formation rate as a function of the redshift presented by Miyaji et al. (1998), and compare the corresponding GRB rate with the data corrected for selection and calibration effects by Petrosian (priv. comm.) for the 4B BATSE catalogue.

In their paper Miyaji et al. (1998), calculated the distribution of Seyfert galaxies as a function of redshift. We considered this same distribution for the SFR, and approximated it with the following function:

$$\psi (z)=A\mathrm{exp}(az),$$
(9)
where $`z`$ is the redshift, $`A10^6[h^3\mathrm{Mpc}^3]`$, $`a=2.5`$, and we used $`\mathrm{\Omega }=1`$, $`\mathrm{\Lambda }=0`$ as a simple reference. This law is valid up to $`z=2`$, and then continues as a constant to high redshifts.

This SFR as a function of the fluence is given by:

$$\psi (f,f_{})=A\mathrm{exp}\left[a\left[\left(\frac{1}{2}\sqrt{\frac{f_{}}{f}}+1\right)^21\right]\right].$$
(10)
Using the LF given in Eq. 4, the rate of GRB fluence is obtained from the product of the terms of the Eqs. 68, and Eq. 10, where we now integrate over $`f_{}`$, that is:

$`N(f)df`$ $`=`$ $`2{\displaystyle \frac{c^3}{h^4}}A\varphi _0(f^\beta df){\displaystyle _{x_{\mathrm{min}}}^{x_{\mathrm{max}}}}x^{1/2}\times `$ (11)
$`\mathrm{exp}\left[\mathrm{a}[({\displaystyle \frac{1}{2}}\sqrt{\mathrm{x}}+1)^21]\right]\times `$
$`\left({\displaystyle \frac{1}{2}}\sqrt{x}+1\right)\left(\sqrt{x}+1\right)^3dx,`$
where $`x=f_{}/f`$, and $`x_{\mathrm{max}}`$ and $`x_{\mathrm{min}}`$ are the limits of integration defined as the maximum and the minimum of the intersection between the interval $`[f_\mathrm{a},f_\mathrm{b}]`$ relative to the LF distribution and the interval in which the ratio $`f_{}/f`$ is defined.

In Fig. 1 we plotted the cumulative curves corresponding to different luminosity function power law indexes ($`\beta =1.5`$ dotted line, $`\beta =2.0`$ dashed line, and $`\beta =3.0`$ solid line) in Eq. 11 and compared them with the corrected data of the 4B catalogue from BATSE. It is evident that even if we change the luminosity function distribution index, we cannot obtain any better fit or calculate any GRB rate using this SFR distribution.

### 2.3 The GRB rate distribution using the Madau SFR

In this paragraph we calculate the GRB rate based on the same procedure as in the last section, using the SFR as a function of redshift presented by Madau et al. (1996) and the luminosity function distribution of Eq. 4.

We approximated the Madau SFR by the combination of two exponential functions. The first part of the curve is:

$$k_1(z)\mathrm{exp}(a_1z),$$
(12)
where $`z`$ is the redshift, and $`a_1=2.9`$ is a constant given by the IR source counts. This law is valid up to $`z^{}=1.7`$, after which it is substituted by the function:

$$k_2(z)\mathrm{exp}(a_2z),$$
(13)
and $`a_2`$ is a constant.

Because of the uncertainties of the dust extinction at high redshift and also the difficulties in the redshift determination for the SCUBA sources (Sanders 1999), the star formation history beyond the redshift $`z=2`$ is still unclear, therefore we considered different slopes for this part of the Madau SFR curve and plotted them in Fig. 2 together with the experimental data given by Steidel (1999).

The cumulative GRB fluence distribution is obtained from Eq. 11, where the SFR function is substituted with Eq. 12 and Eq. 13. The parameters that we can change to fit the data are the redshift $`z^{}`$, the power law index $`a_2`$ in the Eq. 13 and the power law index $`\beta `$ and the upper limit $`f_\mathrm{b}`$ of the interval in which the luminosity function distribution is defined. Following one of the possible curves that fit the data shown by Steidel (1999), we chose $`z^{}=1.7`$, $`\beta `$ determines the slope of the cumulative distribution curve, $`a_2`$ defines curves with different slope, and variations in $`f_\mathrm{b}`$ correspond to little changes in the part of the curve relative to the strongest GRBs, we chose $`f_\mathrm{b}=2.2\times 10^4\mathrm{erg}/\mathrm{cm}^2`$. To reproduce the corrected data by Petrosian of the 4B catalogue, we probed all the different cases changing the slopes of the SFR and the power law index of the fluence distribution. The only values for the power law index of the luminosity distribution that successfully fits the data was $`\beta =1.55`$.

In Fig. 3 we show the cumulative fluence distribution corresponding to the different slopes of Fig. 2 and the power law index $`\beta =1.55`$. We compared these curves with the data corrected for selection and calibration effects by Petrosian for the 4B BATSE catalogue: we can fit these data only if $`a_2`$ is in the range $`0.8÷1.3`$. In Fig. 4 the data are fitted with $`a_2=1.0`$.

Considering the total number of GRBs in the BATSE catalogue, an observing time of 8 years, a volume scale of $`h^310^{10.8}\mathrm{Mpc}^3`$, a beaming factor $`\frac{4\pi }{2\pi \theta ^2}=200\theta _{1\mathrm{j}}^2`$, where $`\theta _{1\mathrm{j}}=\theta _\mathrm{j}/(10^1\mathrm{rad})`$ is the opening angle of the jet, and the factor 22 coming from the integral of the distribution in Eq. 10 calculated with these new SFR and LF, we obtain the following rate of GRBs:

$`10^{5.4}(h^3\theta _{1\mathrm{j}}^2)\mathrm{GRBs}(\mathrm{yr})^1(100\mathrm{Mpc}^3).`$ (14)
This value obtained considering beaming effects is close to the number given by Piran (1999).

The result depends strongly on the power index, while it is not influenced by the interval $`[f_\mathrm{a},f_\mathrm{b}]`$ of the LF distribution we used. We defined a lower limit for the fluence equal to $`10^8\mathrm{erg}/\mathrm{cm}^2`$, and a ratio $`f_\mathrm{b}/f_\mathrm{a}`$ equal to $`10^4`$. A change in the upper limit of this interval corresponds to a small change in the tail of the fluence distribution curve, for $`x>10^4`$. In our model, this interval in the fluence corresponds to an interval in the initial energy deposited in the jet in the range $`[10^{48},10^{52}]`$ ergs.

## 3 Maximum energy available

Another important step is to check if it is possible to produce neutrinos and high energy cosmic rays with our model (Pugliese et al. 1999). Therefore we calculated the maximum energy available in the emission region. Irrespective of any acceleration mechanism, the maximum energy of charged particles, here protons, is that given by the spatial limit in the comoving frame, where the gyromotion just fits the space available:

$`E_{\mathrm{max}}^{(\mathrm{ob})}={\displaystyle \frac{1}{2}}\gamma _{\mathrm{sh}}e(B_\mathrm{j}\mathrm{\Delta }z)^{(\mathrm{sh})}.`$ (15)
Here $`e`$ is the charge of an electron, $`B`$ is the magnetic field in the shock frame given by $`B_\mathrm{j}^{(\mathrm{sf})}(t)=10.24(E_{51}^{1/4}\dot{M}_{5\mathrm{j}}^{3/4}v_{0.3}^{3/4})(ϵ^{1/2}\theta _{1\mathrm{j}}^1\gamma _{\mathrm{m},2}^{1/2})t_5^{3/4}`$, and $`\mathrm{\Delta }z`$ is the thickness of the emission region. For the calculation of the thickness we have two possibilities, and each of them corresponds to one of the following cases:

For the first case we used the thickness of the emission region in the shock frame (sh), given by $`z/(4\gamma _{\mathrm{sh}})`$:

$`E_{\mathrm{max}}^{\mathrm{I}(\mathrm{ob})}`$ $`=`$ $`\gamma _{\mathrm{sh}}\left(eB^{(\mathrm{sh})}{\displaystyle \frac{z_\mathrm{j}}{4\gamma _{\mathrm{sh}}}}\right)1.91\times 10^{21}(E_{51}^{1/4}\dot{M}_{5\mathrm{j}}^{1/4}\times `$ (16)
$`v_{0.3}^{1/4})(ϵ_0^{1/2}\theta _{1\mathrm{j}}^1\gamma _{\mathrm{m},2}^{1/2})t^{1/4}\mathrm{eV}.`$
In the second case we considered the width of the shell:

$`E_{\mathrm{max}}^{\mathrm{II}(\mathrm{ob})}`$ $`=`$ $`\gamma _{\mathrm{sh}}(eB^{(\mathrm{sh})}\theta _\mathrm{j}z_\mathrm{j})0.93\times 10^{23}E_{51}^{1/2}\times `$ (17)
$`(ϵ_0^{1/2}\gamma _{m,2}^{1/2})t^{1/2}\mathrm{eV}.`$
At this point it is necessary to check when the thickness of the emission region is lower than the width, which corresponds to $`E_{\mathrm{max}}^{\mathrm{I}(\mathrm{ob})}<E_{\mathrm{max}}^{\mathrm{II}(\mathrm{ob})}`$. To do this we calculated the time at which these two quantities are equal:

$`t_{(E_{\mathrm{max}}^{\mathrm{I}(\mathrm{ob})}=E_{\mathrm{max}}^{\mathrm{II}(\mathrm{ob})})}5.62\times 10^6(E_{51}^{1/4}\dot{M}_{5\mathrm{j}}^{1/4}v_{0.3}^{1/4})\theta _{1\mathrm{j}}\mathrm{s}.`$ (18)
It means that for about two months $`E_{\mathrm{max}}^{\mathrm{I}(\mathrm{ob})}<E_{\mathrm{max}}^{\mathrm{II}(\mathrm{ob})}`$. Obviously, only $`E_{\mathrm{max}}^{\mathrm{I}(\mathrm{ob})}`$ is valid, so after 10 seconds the upper limit of the maximum energy available in our model is about $`1.07\times 10^{21}\mathrm{eV}`$, while after two months the energy to consider is $`E_{\mathrm{max}}^{\mathrm{II}(\mathrm{ob})}`$. Adiabatic losses will diminish these energies for charged particles like protons.

## 4 Cosmic ray contribution

Cosmic rays are ionized nuclei, mainly protons, that extend from low energies (few hundred MeV) up to very high energies (about $`3\times 10^{20}`$ eV). Their spectrum is described by a power law $`(dN/dE)=E^{(\kappa +1)}`$, and shows two breaks in the slope. From low energies up to about $`5\times 10^{15}`$ eV, known as the knee, the spectrum follows a pure power law with $`\kappa 1.7`$. The detailed shape of this break and the precise position are still unknown. Beyond the knee up to about a second break point at $`3\times 10^{18}`$ eV, known as the ankle, the pure power law has an index $`\kappa 2`$.

It has been proposed (Biermann 1993, Stanev et al. 1993), that three components contribute to the cosmic ray spectrum: a) explosion of supernovae into a homogeneous interstellar medium (ISM), accelerating particles up to energies of about $`10^5`$ GeV. The spectrum for these particles is a power law with an index of -2.75, after considering the leakage from our Galaxy. b) Explosion of stars into their former stellar wind (like Wolf Rayet stars), producing particles with energies up to about $`3\times 10^9`$ GeV. The corresponding spectrum switches at the knee from -2.67 to -3.07 and this difference in the spectral index derives from a diminution of the particle curvature drift energy gain. In Biermann’s cosmic ray model for the Galactic component (Biermann 1997), the energetic protons are produced in the shocks of supernova explosions in the interstellar medium, while all the heavier elements are produced in the shock waves propagating in the stellar wind of the progenitor star. c) Production of particles with energies up to $`10^{12}`$ GeV from the hot spots of Fanaroff Riley class II radio galaxies. Their spectrum has an index -2 at the source, and one needs to take the interaction with the cosmological microwave background into account.

We expect the spectrum of cosmic rays above about $`5\times 10^{19}`$ GeV to be strongly attenuated because of the interaction of nuclei and protons with the 2.7 K cosmic microwave background, giving rise to the so-called Greisen-Zatsepin-Kuzmin (GZK) cut-off. The extragalactic sources cannot produce all the total cosmic ray energy density observed at Earth, but they could give a contribution to the ultra-high energy (UHE) part of the spectrum. The origin of the cosmic rays above $`3\times 10^{18}`$ eV is not yet clear, but the common idea is that they are extragalactic and probably connected with the most powerful radio galaxies (Biermann $`\&`$ Strittmatter 1987, Berezinsky $`\&`$ Grigor’eva 1988, Rachen $`\&`$ Biermann 1993, Rachen et al. 1993). At the moment, our knowledge of the sources of the highest energy cosmic rays is limited by the small number of events detected by the present experiments. Considering that at $`10^{20}`$ eV the rate of cosmic rays is about 1 event per $`\mathrm{km}^2`$ per century, it is clear that to detect them it is necessary to have both large aperture detectors and a long exposure time.

In this general context, it is interesting to check the eventual Galactic and extragalactic energetic contribution given by GRBs to the cosmic ray spectrum in our model. In fact GRBs seem to be very powerful explosions and they inject a large amount of energy and elementary particles into the interstellar medium.

We used the GRB rate obtained with the SFR from Madau and two different energetic approaches to calculate the contribution from GRBs to the cosmic rays and the neutrino spectra. First we considered that each GRB gives the same contribution equal to $`10\%`$ of a fixed initial energy of $`10^{51}\mathrm{ergs}`$. Secondly we assume that each GRB contributes proportional to its own fluence, which we assume is distributed with a power law, corresponding to a range of initial energy of $`[10^{48},10^{52}]`$ ergs.

### 4.1 Extragalactic contribution

For GRBs it is important to identify which particles contribute to the cosmic ray flux. As Rachen $`\&`$ Mészáros (1998) showed, during the main burst protons lose most of their energy because of adiabatic expansion, while neutrons can be better candidates to obtain ultra high energy cosmic rays (UHECR) and neutrinos. In fact neutrons carry about $`80\%`$ of the proton energy, and because they are not coupled to the magnetic field, they can escape the fireball and through the $`\beta `$-decay give a cosmic ray proton spectrum. We followed this same logic in our calculations below.

We assumed that the total energy discharged by each GRB in the ISM for hadrons as well as neutrinos is equal to $`10\%`$ of the initial energy $`E_{51}=E/(10^{51}\mathrm{erg})`$ deposited in the jet, that is $`\eta _{10}=(10\%E)/(10^{50}\mathrm{erg})`$. It means that in the Hubble time, the energetic contribution per unit volume inside the whole universe given by all the GRBs is equal to $`10^{21.0}(h^3\theta _{1\mathrm{j}}^2\eta _{10})\mathrm{erg}/\mathrm{cm}^3`$.

To derive the spectrum of GRBs and to compare it with the one of cosmic rays out of our Galaxy, we need to calculate the normalization factor $`N`$ in $`N\left(\frac{E}{E_0}\right)^2dE`$, where $`N`$ is expressed in $`[\mathrm{GeV}^1\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1]`$.

$`N`$ is obtained directly from the integration of this power law, remembering that the result of this integral is equal to the energy per volume produced by all the GRBs:

$$_{E_1}^{E_2}\left[\frac{4\pi }{c}N\left(\frac{E}{E_0}\right)^2EdE\right]=10^{21.0}(h^3\theta _{1\mathrm{j}}^2\eta _{10}).$$
(19)
The corresponding value is

$$N10^{8.8}(h^3\theta _{1\mathrm{j}}^2\eta _{10})\mathrm{GeV}^1\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1.$$
(20)
We assumed that the total energy discharged by each GRB in the ISM for hadrons as well as neutrinos is proportional to its own fluence in the $`\gamma `$ band. To obtain the energetic contribution per unit volume inside the whole universe given by all the GRBs in the Hubble time, we integrated the Luminosity Function distribution and had $`10^{23.6}(h^3\theta _{1\mathrm{j}}^2\eta _{10})\mathrm{erg}/\mathrm{cm}^3`$.

The corresponding normalization factor $`N`$ in the curve to plot $`N\left(\frac{E}{E_0}\right)^2dE`$ is

$$N10^{11.4}(h^3\theta _{1\mathrm{j}}^2\eta _{10})\mathrm{GeV}^1\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1.$$
(21)
In Fig. 5 we plot the all particle energy spectrum as measured by different ground-based experiments (Biermann $`\&`$ Wiebel-Sooth 1999), the spectrum in the case in which each GRB gives the same contribution to CR’s (solid line) and the spectrum corresponding to a contribution from each GRB proportional to its own fluence (dashed line). The dotted lines represent unlikely contributions because beyond $`10^{18}`$ eV the interactions with the cosmological microwave background are relevant and they make the curves much flatter, and therefore much lower. From this graphic it is clear that in our model GRBs do not give any energetic contribution to the extragalactic cosmic ray spectrum.

### 4.2 Galactic contribution

To calculate the contribution from GRBs to the Galactic cosmic ray spectrum it is necessary to know the GRBs production rate in our Galaxy. We considered that the ratio between the in/out GRB rates is equal to the ratio of the infrared (IR) luminosity inside and outside our Galaxy. In our Galaxy this luminosity is $`L_{\mathrm{IR}}10^{10}L_{}10^{43.6}\mathrm{erg}/\mathrm{s}`$ at $`60\mu \mathrm{m}`$. Using the equation (1) of Malkan $`\&`$ Stecker (1998), we have an extragalactic IR luminosity $`L_{\mathrm{IR}}10^{44.6}h^3\mathrm{erg}(100\mathrm{Mpc}^3)^1`$. This means that the Galactic GRB rate is equal to $`10^{6.4}(h^3\theta _{1\mathrm{j}}^2\eta _{10})\mathrm{GRBs}\mathrm{per}\mathrm{year}`$.

Following the same procedure used for the extragalactic case, we calculated first the case in which each GRB gives the same contribution to the cosmic ray spectrum. We obtained a diffuse density energy of GRBs in our Galaxy equal to $`10^{16.0}(h^3\theta _{1\mathrm{j}}^2\eta _{10})\mathrm{erg}/\mathrm{cm}^3`$.

In our Galaxy, cosmic rays have a time scale to escape with an energy dependence that goes as $`E^{1/3}`$ (Biermann 1995, Biermann et al. 1995). This term has to be taken in account to calculate the total spectrum of the primary cosmic rays inside our Galaxy. It means that it is necessary to multiply the injection spectrum $`E^2`$ times the leakage term $`E^{1/3}`$ to obtain the final plot of the contribution from GRBs to the Galactic cosmic ray spectrum: $`N\left(\frac{E}{E_0}\right)^{7/3}dE`$. The normalization factor is obtained using the same procedure of the last section:

$$_{E_1}^{E_2}\left[\frac{4\pi }{c}N\left(\frac{E}{E_0}\right)^2EdE\right]=10^{16.0}(h^3\theta _{1\mathrm{j}}^2\eta _{10}).$$
(22)
The corresponding value is

$$N10^{3.9}(h^3\theta _{1\mathrm{j}}^2\eta _{10})\mathrm{GeV}^1\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1.$$
(23)
For the second case we assume that each GRB contribution is proportional to its own fluence in the $`\gamma `$-band, and integrating the luminosity function distribution we had an energy density in our Galaxy equal to $`10^{18.7}(h^3\theta _{1\mathrm{j}}^2\eta _{10})\mathrm{erg}/\mathrm{cm}^3`$.

The corresponding normalization factor $`N`$ in the spectrum $`N\left(\frac{E}{E_0}\right)^{7/3}dE`$ is

$$N10^{6.5}(h^3\theta _{1\mathrm{j}}^2\eta _{10})\mathrm{GeV}^1\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1.$$
(24)
In Fig. 6 we compared the all particle energy spectrum as measured by different ground-based experiment (Biermann $`\&`$ Wiebel-Sooth 1999), with the spectrum from GRBs in the case that each of them gives the same contribution (solid line) and with the one in which the contribution is proportional to the fluence (dashed line). The dotted lines beyond $`10^{18}\mathrm{eV}`$ show where interactions with the microwave background may become relevant. Since massive star formation is highest in the Galactic central region, any contribution is limited to energies for which the Larmor radius becomes as large as the Galactic disk, and the AGASA data suggest only a small anisotropy at the Galactic center. But even if from an energetic point of view, GRBs could give a contribution at these high energies, this is ruled out considering that the temporal interval between two GRBs in our Galaxy is larger than the leakage time of cosmic rays at these energies. Therefore, in our model also in our Galaxy GRBs cannot give any contribution to the cosmic ray spectrum.

## 5 Neutrino production

We also would like to probe the energetic contribution of GRBs to the neutrino flux. Both the very high energy (VHE) neutrinos, with energies in the range $`10^{10}÷10^{17}`$ eV, and the ultra high energy (UHE) neutrinos, with energies $`E10^{17}`$ eV originate in the interactions of protons with photons, through the reactions $`p\gamma n\pi ^+`$, $`\pi ^+\mu ^+\overline{\nu _\mu }`$, $`\mu ^+e^+\nu _\mu \overline{\nu _\mathrm{e}}`$. The efficiency for the neutrino production depends on the fraction of proton energy converted into charged pions and on how much energy pions and muons keep before decaying.

Some authors (Waxman $`\&`$ Bahcall 1997), proposed that GRBs can be associated with neutrinos produced in the photohadronic reaction. Pions come from the interactions between accelerated protons and gamma rays in the fireball, and their decay produces neutrinos and anti-neutrinos together with other elementary particles. Waxman $`\&`$ Bahcall (1997) gave an upper limit for the corresponding neutrino flux.

Other authors (Rachen $`\&`$ Mészáros 1998), argued that the upper limit obtained in this way is optimistic. In fact, protons emitted in the earliest burst do not have enough energy to leave the expansion region because of adiabatic losses. Instead neutrons can easily escape and contribute to the neutrino flux.

To check what is the neutrino flux in our model, we calculated the initial photon density number in the $`p\gamma `$ interaction from the synchrotron photons in our model. We obtained a number of $`\mathrm{photons}/\mathrm{cm}^3`$ equal to:

$`L^{(\mathrm{sh})}/\left(4\pi {\displaystyle \frac{z}{4\gamma _{\mathrm{sh}}}}{\displaystyle \frac{100\mathrm{KeV}}{\gamma _{\mathrm{sh}}}}c\right)`$ $`=`$ $`6.94\times 10^{16}(E_{51}^{1/4}\dot{M}_{5\mathrm{j}}^{5/4}\times `$ (25)
$`v_{0.3}^{5/4})\theta _{1\mathrm{j}}^2t^{7/4}.`$
To obtain the number of hits in the $`p\gamma `$ collisions, we used a cross section $`\sigma _{\mathrm{p}\gamma }=2.70\times 10^{28}\mathrm{cm}^2`$, corresponding to the maximum of the curve describing this reaction. This implies a number of hits per proton equal to $`0.99(E_{51}^{1/4}\dot{M}_{5\mathrm{j}}^{5/4}v_{0.3}^{5/4})\theta _{1\mathrm{j}}^2t^{3/4}`$.

To calculate the energy rate of neutrinos in our model we assumed that about $`20\%`$ of the energy of protons $`E_\mathrm{p}=10^{50}\mathrm{erg}`$ goes into neutrinos and that in a first approximation there are no multiple hits for the same proton in the $`p\gamma `$ interaction. Following the same procedure used to calculate the extragalactic cosmic ray contribution from GRBs, we obtained that the normalization factor corresponding to the spectrum $`(E/E_0)^2`$ in the case in which each GRB gives the same contribution to the neutrino spectrum is $`N1.55\times 10^{10}(h^3\theta _{1\mathrm{j}}^2\eta _{10})\mathrm{GeV}^1\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1`$. This spectrum is plotted in Fig. 7 and it shows that in our jet model GRBs give a contribution to the cosmological neutrino flux at a low level.

## 6 Discussion and conclusions

In the first part of this article we calculated the GRB rate and compared the corresponding cumulative distribution in fluence with the observational data. There were two main points to decide on: a) which luminosity function distribution and b) which star formation rate were the best to reproduce the data.

We checked if in our jet model GRBs were standard candles. But we did not obtain any good fit, therefore we tried a power law for the luminosity distribution, $`\varphi (f)f^\beta `$. Together with this function we used the SFR from Miyaji et al. (1998), in which the SFR grows linearly from $`z=0`$ up to about $`z=2`$ and after this redshift is flat up to $`z=6`$. There was not a good agreement between the distribution in fluence we obtained and the corrected data by Petrosian for the 4B BATSE catalogue. Therefore we used the SFR model from Madau et al. (1996), in which the rate follows the same behavior of Miyaji’s up to a redshift $`z^{}`$, and after this redshift the SFR begins to decrease. There are still some uncertainties about the shape of this second part of the distribution, so we considered different curves with different slopes (depending on a constant $`a_2`$) for this decreasing part. We have only two free parameters that we can change to fit the data, the power law index $`\beta `$, and the exponential index $`a_2`$. The redshift $`z^{}`$, and the upper limit $`f_\mathrm{b}`$ of the interval in which the LF is defined are not really free parameter because their range is limited by the observations. We can reproduce the 4B BATSE corrected data by Petrosian using the following values: redshift $`z^{}=1.7`$, $`\beta =1.55`$, $`f_\mathrm{b}=2.2\times 10^4`$ and $`a_2`$ in the range \[0.8,1.3\]. These parameters are remarkably constrained. Thus, given a final model for GRBs and a cosmological model, we may be able to derive strong limits on cosmological parameters.

The key point of the second part of the article is the calculation of the GRB rate inside and outside our Galaxy. The results we obtained depends mainly on three other parameters, the value of the Hubble constant $`H_\mathrm{o}`$, the opening angle of the jet $`\theta _\mathrm{j}`$, and the power law index $`\alpha `$ we assumed for the electron distribution in our model (see Eq. 1).

The choice of $`\alpha =2`$ has been done in the first version of our model, where even if we simplified in many places the physics used, we obtained a good agreement with the data. Any possible small changes in this interval will influence the results obtained in this work in a marginal way, because the dependence on $`\alpha `$ in the expression of $`N`$ is not strong. However, putting $`\alpha =3`$, which is not suggested by the data, would strongly influence our results.

On the other hand, changes in $`H_\mathrm{o}`$ and $`\theta _\mathrm{j}`$, because of the strong dependence in $`N`$ will influence the results. It is interesting to note that both, a lower value of the Hubble constant and a smaller opening angle of the jet, go into the direction of decreasing $`N`$. But these changes cannot influence our results because the interactions of cosmic rays with the microwave background and the large difference between the Galactic GRBs and the leakage time of CR’s rule out that in our jet model GRBs can give any contribution beyond $`10^{18}`$ eV.

In the calculation of the contribution from GRBs to the neutrino flux, it is important to define the parameters that characterize the neutrino production. The number of hits in the $`p\gamma `$ collisions has been obtained considering only the major photohadronic interaction channel, the one that gives a single pion. At higher energies it is possible to have the channels for the production of multi-pions, $`p\gamma n2\pi ^+\pi ^{}`$ and $`p\gamma n3\pi ^+2\pi ^{}`$, but as the energy increases the corresponding cross sections decrease. We used only the first channel with the highest cross section, because the energies involved are not high enough to require secondary channels.

In the context of our jet-disk symbiosis model for GRBs, there is only one way to make the extragalactic contribution to CRs significant at the highest energies, and that is to drastically increase the CR energies per GRB deposited. We consider this implausible in the context of our model because of energy conservation requirements. We use $`10\%`$ of the entire energy available, and so the CR contribution may be increased over our simple calculation by a factor of a few, but not more. This implies that the contributions given by each GRB to CR and neutrino flux cannot be much bigger than the energy emitted in the $`\gamma `$-ray band.

Using a relatively small set of parameters, the jet-disk symbiosis model applied to GRBs, a tested star formation rate and the fundamental physics of the photohadronic interactions we arrive at the conclusion that GRBs are not able to give any significant contribution to the high energy cosmic ray spectrum both inside and outside our Galaxy and predict only a low flux of neutrinos.

A main conclusion of this work is that fitting the corrected fluence distribution of GRBs with the jet-model is well possible. The fit allows strong constraints of the star formation rate as a function of redshift.

###### Acknowledgements.

GP thanks S. S. Larsen and E. Ros for useful discussions. We are grateful to V. Petrosian for the 4B corrected data he gave us, and for discussions. We want also to thank our referee R.A.M.J. Wijers for the helpful advice and comments that he gave us to improve our article. PLB wishes to thank T. Piran for extensive discussions of GRBs. GP is supported by a DESY grant 05 3BN62A 8. HF is supported by a DFG grant 358/1-1&2.


text:
# Auroral field-aligned currents by incoherent scatter plasma line observations in the E region

## 1 Introduction

The incoherent scatter spectrum consists mainly of two lines, the widely used and relatively strong ion line and the very weak easily forgotten broadband electron line. There is also another line present, the plasma line, due to scattering from high frequency electron waves, namely Langmuir waves. From the downgoing and upgoing Langmuir waves, two plasma lines can be detected by the radar. The frequency shift from the transmitted signal is the frequency of the scattered Langmuir wave plus the Doppler shift caused by electron drift. Plasma lines can be used to measure the electron drift and hence the line-of-sight electric current. The problem in ion line analysis with the uncertainty of the radar constant can be solved by the plasma line frequency determination and when that is done the speed of measurement can be significantly increased by including the plasma line in the ion line analysis. However, since the frequency of the plasma lines is not known beforehand, and the frequency is varying with height, it is difficult to measure them with enough resolution.

There have been a number of reports on plasma line measurements and their interpretation. Most of them have discussed the frequency shift from the transmitted pulse and the scattering has mainly been from the F-region peak, e.g. Showen (1979), Kofman et al. (1993) and Nilsson et al. (1996a). The latter two showed also that the simple formula for the Langmuir wave frequency,

$$f^2=f_p^2(1+3k^2\lambda _D^2)+f_c^2\mathrm{sin}^2\alpha $$
(1)
where $`f_p`$ is the plasma frequency, $`k`$ the wave number $`f_c`$ the electron gyro frequency, $`\lambda _D`$ the electron Debye length and $`\alpha `$ the angle between the scattering wave and the magnetic field, is valid to within a few kHz and thus enough to set the radar system constant. To be able to deduce any electron drifts, or current, out of the positions of the lines these authors also show that Eq. (1) is not sufficient, and it is necessary to carry out more accurate calculations. Hagfors and Lehtinen (1981) had also to go to further expansions in deriving the ambient electron temperature from the plasma lines. The fact that so many reports deal with the F region peak is due to the altitude profile shape of the plasma line frequency, which according to Eq. (1) will also show a peak around that height. The measurements can thus be made relatively easily using rather coarse height resolution but good frequency resolution and only detect the peak frequency. Measurements using the same strategy, but at other heights, have been made with a chirped radar by matching the plasma line frequency height gradient and the transmitter frequency gradient (Birkmayer and Hagfors, 1986; Isham and Hagfors, 1993). This technique allows determination of the Langmuir frequency with very high frequency resolution, but do not use the radar optimally, since the chirped pulse cannot be used for anything else than plasma line measurements.

The enhancement of the plasma lines, which occurs in the presence of suprathermal electron fluxes (Perkins and Salpeter, 1965), either photoelectrons or secondaries from auroral electrons, has been investigated by Nilsson et al. (1996b), where they also calculate predictions of plasma line strength for different incoherent scatter radars and altitudes. They also show that the power of the plasma line is rather structured with respect to ambient electron density, depending on fine structure in the suprathermal distributions due to excitations of different atmospheric constituents.

Incoherent scatter plasma lines in aurora are more difficult to measure since the variations in the plasma parameters are strong with large time and spatial gradients. Reports of auroral plasma lines in the aurora are also more rare, and most of them are based on too coarse time resolution (Wickwar, 1978; Kofman and Wickwar, 1980; Oran et al., 1981; Valladares et al., 1988), with resolutions ranging from 30 seconds up to 20 minutes. The enhancements over the thermal level were high, but consistent with what could be expected of model calculations of suprathermal electron flux. They also tried to calculate currents and electron temperature from the frequency shifts of the plasma lines but with very large error bars. Kirkwood et al. (1995) used the EISCAT radar and the filter bank technique and recorded much higher intensities of the plasma lines, since they got down to resolutions of 10 seconds, and showed also that the plasma-turbulence model proposed by Mishin and Schlegel (1994) was not consistent with the data, but could be explained by reasonable fluxes of suprathermal electrons.

In this paper we present data obtained with the high resolution alternating code technique (Lehtinen and Häggström, 1987), as was also done for F-region plasma lines by Guio et al. (1996), with even higher intensities due to the time resolutions of 5 seconds. An interesting, but at the time of the experiment not realisable at EISCAT, technique would have been the type of coded long pulses used by Sulzer and Fejer (1994) for HF-induced plasma lines. From relative strengths between up- and downshifted lines we detect a general trend of upgoing field-aligned currents in the diffuse aurora carried by the suprathermal electrons. We propose a generalisation of the theoretical incoherent scatter spectrum, to include multiple shifted electron distributions, and in one example we do a full 7-parameter fit of the incoherent scatter spectrum, including the enhanced plasma lines assuming Maxwellian secondary electrons, resulting in the first radar measurement of its flux and a current carried by the thermal electrons.

## 2 Experiment

The measurements we present were collected by the 930 MHz EISCAT UHF incoherent scatter radar, with its transmitter located at Ramfjordmoen in Norway (69.6 °N, 19.2 °E, L=6.2). The signals scattered from the ionosphere were received at stations in Kiruna, Sweden and Sodankylä, Finland as well as at the transmitting site. General descriptions of the radar facility are given by Folkestad et al. (1983) and Baron (1984). Local magnetic midnight at Ramfjordmoen is at about 2130 UT. The aim, in the Swedish-Japanese EISCAT campaign in February 1999, was to measure the ionospheric parameters inside and outside the auroral arcs.

For this a 3 channel ion line alternating code (Lehtinen and Häggström, 1987) experiment, optimised to probe the E-region and lower F-region with as high a speed as possible, was developed. The 16 bit strong condition alternating code with bitlengths of 22 $`\mu `$s was used, giving 3 km range resolution, and with a sample rate of 11 $`\mu `$s the range separations in consecutive spectra were 1.65 km. Fig. 1 shows the transmission/reception scheme of the first 20 ms of the radar cycle. The whole alternating code sequence takes about 0.3 s to complete. During this period the incoherent scatter autocorrelation functions (ACF) at the probed heights should not change significantly for the alternating codes to work. In order to keep this as short as possible, the short pulses, normally used for zerolag estimation, were dropped and instead a pseudo zero lag, obtained from decoding the power profiles of the different codes in the alternating code sequence, was used. Fig. 2 shows the range-lag ambiguity function for this lag centred around 0 $`\mu `$s in the lag direction, but the main contribution to the signal comes from around 7 $`\mu `$s. The more normal lag centred at 22 $`\mu `$s is also shown for comparison. The range extents are rather similar but the power is, of course, considerably lower for the pseudo zero lag. Nevertheless, this is taken care of in the analysis and this lag is rather important in events with high temperatures giving broad ion line spectra or narrow ACFs. The transmitting frequencies were chosen to give maximum radiated power for a given high voltage setting.

The monostatic plasma line part of the experiment used two channels covering the same ranges as the ion line but with 3.3 km range separations between the spectra. The frequency setup for the experiment, illustrated in Fig. 3, gives the possibility for 3 upshifted and 3 downshifted bands. For 3 MHz and 6.5 MHz frequency shifts, both the up- and downshifted plasma lines were measured. In addition there was a 4 MHz upshifted and a 5.5 MHz downshifted band. The width of these bands should have been set to match the bit length of the codes used, 50 kHz, but unfortunately this was not the case and 25 kHz wide filters were erroneously used. This gave naturally less signal throughput, and in Fig. 2 the corresponding range/lag ambiguity functions for the plasma line channels are included for the first two lags. The decoding still works, giving just slightly increased unwanted ambiguities, but above all the pseudo zero lag is moved out to a larger lag value. This fact, with one exception, almost ruled out the possibilities to measure plasma line spectra because, as will be shown here later, they are likely to be rather wide, due to large time and height gradients of the plasma parameters. This lead the analysis to use mainly the undecoded zerolag, which after integration over the different codes, almost resembles the shape of a 352 $`\mu `$s (16$`\times 22\mu `$s) long pulse.

The experiment contained a large number of antenna pointings in order to follow the auroral arcs, but as this day was cloudy over northern Scandinavia the transmitting antenna was kept fixed along the local geomagnetic field line. The remote sites, receiving only ion lines, were monitoring the same pulses as the transmitting site, and were used to measure the drifts in the F-region, to derive the electric field. Thus, these antennas intersected the transmitted beam at the F-region altitude giving the best signal, for this day mostly at 170 km.

## 3 Measurements

### 3.1 Ion line

The experiment started at 1900 UT on 14 February 1999 and continued until 2300 UT. Fig. 4 shows an overview of the parameters deduced from the ion line measurements, which were analysed using the on-line integration time, 5 seconds, in order to be compared to the plasma line data. The analysis was done using the GUISDAP package (Lehtinen and Huuskonen, 1996), but a correction of 45% of the radar system constant used in the package had to be invoked to fit the plasma line measurements according to Eq. (1). This short integration time was possible due to the highly optimised mode used, with all the transmitter power concentrated to the E-region. In range, some integration was done, so that at lower heights 2 range gates were added together and with increasing height the number of gates added together increases to 15 in the F-region.

There is a rather strong E-region from the start, but no real arcs, and we interpret this as diffuse aurora. The peak electron density shows some variation, but as time goes the E-region ionisation decreases until 2050 UT, where it is almost gone. The density peak during this time was at around 120 km altitude and the lower edge of the E-region at 110 km, but at times the ionisation reaches down to 100 km. At 2050 UT and onwards until 2240 UT the ionosphere above Tromsø became more active and several auroral arcs passed the beam. Around some of the arcs there are short-lived enhancements of electric field, seen as F-region ion and E-region electron temperature increases. In the last 10 minutes of the experiment the arc activity disappeared and again there was diffuse aurora. From the field aligned ion drifts it is evident that there is a rather strong wave activity in the diffuse aurora until 2050 UT, while it is not so clear in the continuation of the experiment. The last panel with the inferred electron density from the pseudo zero lag shows the same features as the fitted density panel but with highest possible resolution since no height integration is made.

### 3.2 Plasma line

Since the analysis of the plasma lines was forced to handle the undecoded zerolags of the alternating codes, it was necessary to analyse their profile shape. Fig. 5 shows how this analysis was performed. From the fitted parameters of the ion line, electron density and temperature, a profile of the approximate Langmuir frequency can be calculated using Eq. (1). When probing at fixed frequency, there will be scattered signal only from heights where the probing and Langmuir frequencies match each other. The effect of the undecoded zerolag is similar to the one where a normal long pulse is used, but with a lower signal strength. So, the profile shape will be a square pulse centred on the corresponding altitude, since no gating is performed. Because the signal strengths are rather weak compared to the system temperature, there is a great deal of noise in the profile shape. To be able to extract the altitude and signal power, a fitting procedure need to be performed. For this a continous function and a good first guess is needed and the measured plasma line profile was convolved with the pulse shape to have a triangular shape and also showing a peak close to the matching height.

The plasma line part of the experiment is overviewed in Fig. 6. The signal strengths shown should be compared with the UHF system temperature of about 90 K. At first glance there is almost nothing in the upshifted part, but a more careful look shows weak signals between 1940 and 2030 UT and after 2250 UT, corresponding to the occurrence of diffuse aurora. Similar echoes can also be seen in the downshifted channel, and are due to plasma lines at 3 MHz offset from the transmitted frequencies, according to the ion line measurements. In the downshifted part after 2100 UT, frequent events of rather strong signals in phase with auroral arcs pass the beam. Most of these are from the 5.5 MHz shift, but some of them also are due to plasma lines at 6.5 MHz. Such events are less frequent and for most of them there are also signals in the upshifted part.

Due to the fact that the same channel is used for several frequencies, there can be several altitudes that fulfill the matching condition between Langmuir and probing frequency. This complicates the analysis somewhat, and there have to be several fits with different numbers of triangles superposed on each other. Some examples of this analysis are shown in Fig. 7, from only one plasma line to several both up- and downshifted lines. Using a lower limit of 2 K signal power, the total number of 5 second integration events with enhanced plasma lines for this evening was 256, and a total of 468 plasma line echoes were detected, divided into 220 on 3 MHz, 19 on 4 MHz, 157 on 5.5 MHz and 72 on 6.5 MHz.

## 4 Theory

In order to relate the plasma line measurements to physical quantities it is necessary to investigate the spectrum of the incoherent scatter process. The Nyquist theorem approach, derived in a long series of papers by Dougherty and Farley (1960, 1963), Farley et al. (1961), Farley (1966) and finally Swartz and Farley (1979), arrives at

$$\sigma (\omega )=\frac{N_er_e^2\mathrm{sin}^2\delta }{\pi }\frac{\left|y_e\right|^2_i\frac{\eta _i\mathrm{}(y_i)}{\omega 𝒌𝒗_i}+\left|jk^2\lambda _D^2+_i\mu _iy_i\right|^2\frac{\mathrm{}(y_e)}{\omega 𝒌𝒗_e}}{\left|y_e+jk^2\lambda _D^2+_i\mu _iy_i\right|^2},$$
(2)
where

$$\eta _i=\frac{n_iq_i^2}{N_ee^2},$$
(3)
$$\mu _i=\frac{\eta _iT_e}{T_i},$$
(4)
and the index $`e`$ stands for electrons and $`i`$ for the different ion species. $`N`$ and $`n`$ are the densities, $`r_e`$ the classical electron radius, $`v`$ the bulk velocity, $`q`$ the charge, $`e`$ the electron charge and $`T`$ the temperature. The complex normalised admittance function, $`y`$, contains most of the physics with the plasma dispersion function and have as main arguments the collision frequency and magnetic field. The same result was also reached by Rosenbluth and Rostoker (1962), using the dressed particle approach. In Fig. 8 there is an example of the spectrum, showing clearly the strong ion line around zero offset frequency and the rather weak plasma lines at rather large offsets. The plasma lines become enhanced by a photo electron or auroral electron produced suprathermal electron distribution, and in order to simulate what this extra distribution does to the spectrum a modification to the formula has to be made. First, a rewriting of Eq. (2) following Swartz (1978) has to be performed in order to separate the electron and ion contributions to the spectrum:

$$S(f)=\frac{1}{\pi }\frac{\left|\frac{N_ey_e}{T_e}\right|^2_i\frac{n_iZ_i^2\mathrm{}(y_i)}{f+kv_i/2\pi }+\left|jC_D+_i\frac{n_iZ_i^2y_i}{T_i}\right|^2\frac{N_e\mathrm{}(y_e)}{f+kv_e/2\pi }}{\left|\frac{N_ey_e}{T_e}+jC_D+_i\frac{n_iZ_i^2y_i}{T_i}\right|^2},$$
(5)
where

$$C_D=\frac{k^2ϵK}{e^2},$$
(6)
$`f`$ is the frequency shift, $`Z=q/e`$, $`ϵ`$ is the dielectricity and $`K`$ is the Boltzmann constant. Here, a normalisation of the spectrum has also been performed, so that the zero lag of the corresponding ACF is the raw electron density and the vector velocity is replaced by the line-of-sight velocity. Using a treatment in analogy to the ion contribution, it is now possible to rewrite the spectrum to support a number of Maxwellian electron distributions as

$$S(f)=\frac{1}{\pi }\frac{\left|_e\frac{N_ey_e}{T_e}\right|^2_i\frac{n_iZ_i^2\mathrm{}(y_i)}{f+kv_i/2\pi }+\left|jC_D+_i\frac{n_iZ_i^2y_i}{T_i}\right|^2_e\frac{N_e\mathrm{}(y_e)}{f+kv_e/2\pi }}{\left|_e\frac{N_ey_e}{T_e}+jC_D+_i\frac{n_iZ_i^2y_i}{T_i}\right|^2}.$$
(7)
In Fig. 9 the effect on the plasma lines of a suprathermal distribution with a reasonable density of $`10^7`$ m<sup>-3</sup> and width of 10 eV is shown, being lower than the ionisation energy for most ions. The plasma lines grow considerable and the integrated power over the bandwidth used in the experiment becomes comparable to the power of the ion line. It may also be noted that there is no effect at all seen in the ion line. Perkins and Salpeter (1965) have shown similar calculations, but when their method was based on large expansions to allow non-Maxwellian distributions one can here more directly superpose a few Maxwellian distributions to explain the measurements and even make fits of the spectra taken to get estimates of the suprathermal distributions, which is of great importance in auroral measurements. A consequence of Eq. (7) is that it is also possible to derive the spectrum assuming currents carried by the suprathermal distribution, since it allows different drift velocities on the various distributions. Indeed, Fig. 10, shows differential strengths on the two plasma lines, with upgoing electrons enhancing mainly the downshifted line and downgoing ones the upshifted line.

## 5 Discussion

Plasma line measurements in the active auroral ionosphere are not an easy task, due to the large variations in the ionospheric parameters. The Langmuir frequencies are largely dependent on the ambient electron density, making the line move considerably as the density changes, which it does on time scales of seconds. Moreover, the density height gradient makes the lines very broad when measured over a specific height interval, and at times even broader than the receiver band. A chirped radar would solve only a part of the problem at the cost of transmitter power. These complications make it hard to draw any conclusions on the power in the lines, as one does not know for sure the scattering volume or the time duration of the scattering. Bearing this in mind and to at least minimise these effects, one can nevertheless look at the different distributions in height and power for the different lines to get an idea of their nature, using the on-line integration time of 5 seconds.

The altitude distributions for the different frequencies are shown in Fig. 11. One must note that the signal levels shown are not corrected for range, as the scattering volumes are not known, so signals from a higher altitude are in fact stronger than the corresponding signal from a lower height. It is clearly seen that the strongest signal is the 5.5 MHz line, followed by 6.5, 4 and 3 MHz. The altitude distribution shows more or less the expected dependence on range, but there are some exceptions: In one point at 188 km in the 3 MHz band and for the 5.5 MHz band the strong values between 130 and 140 km seem to be stronger than the others, even taking into account the range effect. However, the number of points are too few to be used as evidence on altitude effects. Most of the echoes are coming from around 120-150 km altitude and Fig. 12 shows a simulation of the expected strength of the plasma line for given background electron density. It shows a peak at around 5.5 MHz and this is also what the experiment shows. A more realistic suprathermal distribution will decrease the returned power for a number of frequencies and one should see the figure as an upper-limit estimate. Indeed, Nilsson et al. (1996b), have made predictions of the expected strength of Langmuir waves for different heights and carefully derived distributions. These predictions are in rather good agreement with the present measurements showing a strong peak between 5 and 6.5 MHz. Although there is some uncertainty on the scattering volume, the observed strengths of the plasma lines can be used to get some estimates of the distribution of the suprathermal population. Assuming a width of 6 eV, one need, to get to the measured strengths for the 3 MHz case, a suprathermal density of about $`510^8`$ m<sup>-3</sup> which is 0.5 % of the thermal population. For the 5.5 and 6.5 MHz cases, it is enough with only $`110^8`$ m<sup>-3</sup>, but the uncertainty of volume is even larger here due to the active environment and the values should maybe be the same as for the lower frequency offsets.

The most interesting thing with incoherent plasma lines is, of course, the possibility to derive differential drifts between ions and electrons, and from these deduce ionospheric currents. For this, one needs to measure the up- and downshifted lines simultaneously. Although the time of the measurement for both lines were not exactly the same in this experiment, the time shifts between them are so small (3-6 ms) compared to the total cycle time (300 ms) of the codes, that this effect is of minor importance. In Fig. 13, the strengths and altitudes of the two concurrently recorded up- and downshifted plasma lines are shown. In general, there are stronger up- than downshifted lines for the 3 MHz case, whereas no such trend can be seen in the 6.5 MHz band. However, there are exceptions to these overall trends and on occasions there are large differences in the signal strengths between the lines. Almost all of the 3 MHz plasma lines were recorded in diffuse aurora and this evident difference in signal power needs a closer examination. When there is a drift of the thermal electrons, the plasma lines shift, and when probing at a fixed frequency, the scatter may not come from the same altitude for the up- and downshifted lines respectively. The strength of the plasma line is also rather altitude dependent due to the damping by the collisions of electrons with ions and neutrals. But Fig. 13 shows no general height difference between the up- and downshifted 3 MHz lines, so this difference in strength cannot be explained by thermal electron bulk drifts. To simulate the effect of current carried by suprathermal electrons, Fig. 14 illustrates the strength of the 3 MHz plasma lines for Maxwellian electron beams of different energies. With no net current the lines are of almost the same strength, and the difference is mainly dependent on where in the receiver band the lines are. However the experiment shows stronger upshifted lines, thus it is evident that the diffuse aurora this night contained fluxes of downgoing suprathermal electrons or, in other words, there was an upgoing current carried by suprathermal electrons. The average power ratio between up- and downshifted lines, 1.4, can be used to estimate the amount of current carried by the suprathermals. As before, using a density of $`510^8`$ m<sup>-3</sup> and temperature of 6 eV, and shifting this population to a current density of 15 $`\mu `$Am<sup>-2</sup> brings to the observed ratios. This current seems somewhat high, but not unrealistic.

For the 6.5 MHz bands there may be a slight difference with respect to the altitude, and that is most likely due to thermal currents causing frequency shifts of the plasma lines. This effect is not very clear, but as lower heights have higher density, or Langmuir frequency, and as the downshifted line is at a slightly lower height, this is most probably an effect of a downgoing current carried by thermal electrons. The correction due to heat-flow in the plasma dispersion function, discussed by Kofman et al. (1993) and confirmed later also by Nilsson et al. (1996a) and Guio et al. (1996), but not taken into account here, would also show the same effect in altitude difference between the lines. The band widths used here, 25 kHz, are much wider than the effect of heat-flow, which is less than or around 1 kHz in the F-region and much lower in the E-region, so that cannot explain the 6.5 MHz height differences. A recent paper by Guio et al. (1998) with proper calculation of the dispersion equation investigates the heat-flow and finds that it is not necessary to invoke the effect at all, but as this paper don’t go to the same extreme the heat-flow has to be considered.

Most of the above discussions on currents are only on directions, but to get any quantitative numbers it is necessary to look at the individual spectra. Of the 256 events of enhanced plasma lines, there is only one that is good enough to investigate. All the other are too broad either due to the Langmuir frequency gradient smearing or the fact that the Langmuir frequency at the height in question is varying during the 5 second time slots. A further reason is that the stationarity condition for the alternating code technique is not fulfilled, due to changes of Langmuir frequency, and the spectra change too much within the 300 ms cycle. Anyway, there is one, and the up- and downshifted bands, together with the ion line band, are shown in Fig. 15. The shifts in this case are around 3 MHz as deduced from the ion line analysis. These spectra were then fitted to the theoretical spectrum in Eq. (7) using 7 ionospheric parameters, the thermal parameters for both ions and electrons and the parameters of a non-shifted suprathermal electron distribution. At the first glance the fit on the ion line seems rather poor, but the fit was done as usual in the time domain with proper weighting on the different lags of the ACF, and in the FFT process to produce Fig. 15 these statistical properties are lost. The fit looks actually better in the time domain, but the frequency domain was chosen in the figure to be more informative. The current density $`j=N_ee(v_iv_e)`$, and with the fitted parameters, the field aligned current carried by the thermal electrons, amount to 12 $`\mu `$Am<sup>-2</sup> downward. This is opposite the general current seen in the plasma line strength for the 3 MHz band, but the magnitudes are comparable. For this case no current could be seen on the suprathermals, and it is evident that the currents are rather structured. Of course, there are many more cases of no plasma lines than enhanced plasma lines in the diffuse aurora, and some cases of only one plasma line either up- or downshifted. It is not surprising that this single example of thermal currents do not follow the general trend. Again, heat-flow is not taken into account, but that should increase the current somewhat. The deduced suprathermal distribution is in good agreement with those Kirkwood et al. (1995) derived from the precipitating flux of primaries; of course in the present case the distribution is Maxwellian and without the fine structure due to the atmospheric constituents, but the numbers are comparable.

## 6 Conclusions

We have made measurements of plasma lines in the active auroral E-region. During 256 periods of 5 second integration we found a total of 468 plasma line echoes, divided into 220 on 3 MHz, 19 on 4 MHz, 157 on 5.5 MHz and 72 on 6.5 MHz. It may seem strange to try to measure plasma lines at such a low frequency as 3 MHz giving low signal levels, but in fact most of the echoes and the most interesting results came from this frequency offset.

The strongest echoes were found at the 5.5 MHz line, and somewhat weaker ones at 6.5 MHz, inside auroral arcs. One must, however, note that the strength measured inside the arcs is mostly a low-limit estimate due to the active environment. The integration period used, 5 s, is rather long in auroral arc conditions, and changes typically occur on shorter time scales. Therefore, effects of gradients in the Langmuir frequency profile, and hence scattering volumes, have not been taken into account.

The Holy Grail in incoherent scatter plasma lines is the possibility to measure currents, and in this case, the field-aligned currents in aurora. The simulations carried out here, the extended full incoherent scatter spectrum with multi-Maxwellian distributions of electrons, show that the strength of the lines is determined by the suprathermal part of the electron distribution, and the frequency mainly by the thermal part. For simultaneous up- and downshifted plasma line differences in intensity we can deduce currents carried by suprathermals, and for differences in frequency, currents carried by thermals.

The simultaneous up- and downshifted frequencies of the 3 MHz line in the diffuse aurora show, on average, an upward field-aligned suprathermal current during the two main periods when they were detected, 1940-2030 UT and 2250-2300 UT. In the arcs in general, there is an indication of downward thermal current as seen from the altitudes of the 6.5 MHz echoes. Of course, no rule is without exceptions, and there are cases where one line is much stronger than the other or the other line is not at all enhanced, indicating strong currents.

In the full 7-parameter fit of the incoherent scatter spectrum with the ion line and the both enhanced plasma lines, we obtained a thermal current with a suprathermal distribution of electrons consistent with distributions derived from precipitating fluxes.

###### Acknowledgements.

One of the authors (I.H.) was working under a contract from NIPR and is grateful to the Director-General of NIPR for the support. We are indebted to the Director and staff of EISCAT for operating the facility and supplying the data. EISCAT is an International Association supported by Finland (SA), France (CNRS), the Federal Republic of Germany (MPG), Japan (NIPR), Norway (NFR), Sweden (NFR) and the United Kingdom (PPARC).