format: text

text:
# 1 Introduction

## 1 Introduction

Very high energy ‘cosmic’ neutrinos with energies in excess of about 10 TeV offer a unique source of valuable information about energetic events far away in the Universe; see, for example, the reviews in Refs. . This has led to the development of neutrino telescopes, which use photo-multiplier tubes to detect the Cerenkov radiation emitted from the charged leptons produced in charged-current neutrino interactions in a large volume of water or ice, deep underground; see, for example, the reviews in Refs. . If we consider muon neutrinos, then up to about 100 TeV the spectrum is dominated by atmospheric neutrinos from the decays of pions or kaons produced by cosmic ray interactions in the Earth’s atmosphere. At higher energies the increased lifetime of these mesons means that they interact before they have the opportunity to decay. Above 500 TeV the decays of the much shorter-lived charmed particles become the dominant source of atmospheric muon neutrinos. These are known as ‘prompt’ neutrinos111For tau neutrinos, we will see that prompt production dominates at about 10 TeV and above.. Their energy dependence follows the original cosmic ray spectrum, while the spectrum of ‘conventional’ atmospheric neutrinos falls off by an extra power of the energy because of the competition between the decay and interaction of the parent meson. It is clearly essential to quantify the flux of ‘prompt’ neutrinos as accurately as possible, since they provide the background to the sought-after ‘cosmic’ neutrinos.

There exist many models of the ‘prompt’ neutrino flux in the relevant $`10^4`$$`10^9`$ GeV energy range, which yield predictions which differ by more than two orders of magnitude. Some of the models are purely phenomenological and have arbitrary continuation to high energy from the domain constrained by accelerator data. However, it has been noted that, with the inclusion of the next-to-leading order (NLO) contributions , perturbative QCD gives a satisfactory description of the observed features of the available accelerator data on charm production, see for example . Moreover, the simplified LO calculation reproduces the same behaviour when multiplied by an overall $`K`$ factor, $`K2.3`$, to account for the NLO contribution. Several authors have therefore used perturbative QCD to predict the prompt neutrino flux .

The perturbative approach, however, faces the same problem of extrapolation to high energies. The LO diagram for forward high energy charm production is shown in Fig. 1. The cross section may be written in terms of the Feynman variable $`x_F=p_L/p_L^{\mathrm{max}}`$, where $`p_L`$ is the longitudinal momentum of the charm quark; at high energies $`x_FE_c/E`$, where $`E`$ is the incident proton energy and $`E_c`$ is the energy of the charm quark. Using the notation of Fig. 1, we have

$$\frac{d\sigma }{dx_F}(ppc+\mathrm{})=𝑑x_1𝑑x_2𝑑zg(x_1,\mu _F^2)\frac{d\overline{\sigma }_{ggc\overline{c}}}{dz}g(x_2,\mu _F^2)\delta (zx_1x_F),$$
(1)
where $`d\overline{\sigma }/dz=\overline{s}d\overline{\sigma }/d\overline{t}`$ with $`z=(m_c^2\overline{t})/\overline{s}`$, and where $`g(x)`$ is the gluon density of the proton. The Mandelstam variables $`\overline{s}`$ and $`\overline{t}`$ refer to the $`ggc\overline{c}`$ subprocess.222Throughout we take the mass of the charm quark to be $`m_c=1.25`$ GeV, following Ref. . We know that this value, taken together with the NLO contribution ($`K2.3`$), gives a good description of all the available fixed-target data for $`c\overline{c}`$, or rather D meson, production (which are in the region $`E250`$ GeV) . The problem is that in the high energy domain of interest we sample the gluon density at very small $`x_2`$; to be specific $`x_2M_{c\overline{c}}^2/2x_Fs10^9`$$`10^4`$, where $`\sqrt{s}`$ is the total $`pp`$ c.m. energy. There are no data which determine the gluon for $`x\stackrel{<}{}10^4`$, and, as a rule, parton distributions for $`x<10^5`$ are not available. For example, Ref.  shows a range of predictions for the prompt flux neutrino obtained by continuing the gluon distribution below $`x=10^5`$ using the power law dependence $`xgx^\lambda `$ with $`\lambda `$ in the range 0–0.5. Of course, the prediction depends crucially on the value taken for $`\lambda `$, and at the highest neutrino energy shown, $`10^9`$ GeV, the rates span about two orders of magnitude.

The goal of the present paper is to diminish the uncertainty in the predictions of the prompt neutrino flux. The major problem is to obtain the most reliable method of extrapolation based on the present understanding of the small $`x`$ regime. In order to do this we begin, in Section 2, by comparing different physically-motivated extrapolations of the gluon to very small $`x`$:

* DGLAP gluon with a double leading log (DLL) extrapolation,
* unified DGLAP/BFKL gluon  with $`x^\lambda `$ extrapolation,
* extrapolation with saturation effects included.

In Section 3 we compare the predictions for the $`x_F`$ distribution of charm quarks produced in high energy $`pp`$ collisions ($`E10^5`$ to $`10^9`$ GeV) using the three models for extrapolating the gluon into the small $`x`$ regime, typically $`x10^4`$ to $`10^8`$. We argue that the extrapolation based on model (iii) is the most reliable, and so for the remainder of the paper we show results for this gluon. To determine the prompt neutrino flux we need to extend the calculation to high energy proton–air collisions. This is done in Section 4, where we also consider the uncertainties associated with the extrapolation based on model (iii).

In Section 5 we describe the formalism that we shall use to calculate the prompt lepton spectra. Starting from the production of $`c\overline{c}`$ pairs from the incoming cosmic ray flux, we allow for their fragmentation into charmed ($`D^\pm `$, $`D^0`$, $`D_s`$ and $`\mathrm{\Lambda }_c`$) hadrons, and for their subsequent semileptonic decays. We include the effects of the lifetime of the charmed hadrons and, also, for the lifetime and decay modes of the $`\tau `$ lepton in the $`D_s\tau \nu _\tau `$ decays. The results of the calculation of the prompt $`\nu _\mu `$ and $`\nu _\tau `$ fluxes are presented in Section 6, and compared with the ‘conventional’ atmospheric fluxes. These latter fluxes arise from $`\pi `$, $`K\mathrm{}`$ decays and $`\nu _\mu \nu _\tau `$ oscillations respectively. We find that the prompt $`\nu _\tau `$ spectrum for $`E>10^4`$ GeV lies much above its conventional atmospheric background, whereas for the prompt $`\nu _\mu `$ spectrum this is only achieved for $`E>10^6`$ GeV. The origins of the prompt $`\nu _\tau `$ flux are the $`D_s\tau \nu _\tau `$ decays which occur with a branching fraction of $`6.4\pm 1.5\%`$ . It is relevant to note that high energy $`\nu _\tau `$’s, unlike $`\nu _\mu `$’s, are not depleted in number by absorption in the Earth. They will always penetrate the Earth due to the $`\nu _\tau \tau \nu _\tau \mathrm{}`$ regeneration sequence . This is clearly crucial for upgoing high energy neutrinos through the Earth and can be important for horizontal neutrinos (particularly of high energy) entering a detector deep underground.

At first sight, it appears that the prompt neutrino flux from $`b\overline{b}`$ production will, relative to the flux of $`c\overline{c}`$ origin, be suppressed, first, by a factor of order $`m_c^2/m_b^2`$ and, second, by the fact that the gluon is sampled at larger $`x`$. However, for the prompt $`\nu _\tau `$ flux of $`b\overline{b}`$ origin, the suppression is partly compensated by the existence of significant $`\tau \nu _\tau `$ semileptonic decays of all the beauty hadrons ($`B^\pm `$, $`B^0`$, $`B_s`$ and $`\mathrm{\Lambda }_b`$), in contrast to just the $`D_s\tau \nu _\tau `$ decays for charm. We calculate the prompt flux arising from $`b\overline{b}`$ production, fragmentation and decay in Section 7. Finally, in the concluding section, we discuss the implications of our results for neutrino astronomy and cosmic ray physics. Also, there, we summarize the uncertainties in the predictions of the prompt neutrino fluxes.

## 2 The gluon at small $`x`$ and high energy $`c\overline{c}`$ production

As mentioned above, one possible method of extrapolating the gluon into the $`x<10^5`$ regime is to resum the leading $`\alpha _S\mathrm{ln}Q^2\mathrm{ln}1/x`$ terms within the DGLAP framework, which leads to a small $`x`$ behaviour333In practice, we use the DLL continuation to $`x<10^5`$ in the form
$$g(x,Q^2)=g(x=10^5,Q^2)\mathrm{exp}\left(\sqrt{\frac{16N_C}{b}\mathrm{ln}\frac{\alpha _S(Q)}{\alpha _S(Q_0)}\mathrm{ln}\frac{x}{x_0}}\sqrt{\frac{16N_C}{b}\mathrm{ln}\frac{\alpha _S(Q)}{\alpha _S(Q_0)}\mathrm{ln}\frac{10^5}{x_0}}\right),$$
where the LO coupling $`\alpha _S(Q)=4\pi /(b\mathrm{log}(Q^2/\mathrm{\Lambda }_{\mathrm{QCD}}^2))`$ with $`n_f=4`$ and $`b=25/3`$. We take $`Q_0^2=1`$ GeV2 and $`x_0=0.25`$. We use MRST2001 partons  with $`\mathrm{\Lambda }_{\mathrm{QCD}}=220`$ MeV. The reliability of this form of DLL extrapolation has been checked using GRV partons  which are tabulated down to $`x=10^9`$.

$$xg(x,Q^2)x_0g(x_0,Q_0^2)\mathrm{exp}\left(\sqrt{\frac{16N_C}{b}\mathrm{ln}\frac{\alpha _S(Q)}{\alpha _S(Q_0)}\mathrm{ln}\frac{x}{x_0}}\right).$$
(2)
We denote this extrapolation by MRST on the Figures below.

As far as we fix the scale $`Q^2`$ and extrapolate to much smaller $`x`$, the leading contribution comes from $`\alpha _S\mathrm{ln}1/x`$ terms, which, at leading order, are resummed by the BFKL equation . Therefore a more reliable extrapolation is obtained by solving a unified DGLAP/BFKL equation  for the gluon. This equation is written in terms of the gluon, unintegrated over its transverse momentum, which should be used with the off-mass shell444That is we use the replacement in (1)
$$xg(x,\mu _F^2)\sigma _{ggc\overline{c}}xg(x,Q_0^2)\sigma _{\mathrm{on}}+_{Q_0^2=1\mathrm{GeV}}^{\mathrm{}}f_g(x,Q_t^2)\sigma _{\mathrm{off}}𝑑Q_t^2$$
with $`x=x_2`$; where $`f_g`$ is the unintegrated gluon density as defined in , and where $`\sigma _{\mathrm{on},\mathrm{off}}`$ are the $`ggc\overline{c}`$ cross sections with the small $`x`$ gluon on-, off-mass-shell. The large $`x_1`$ gluon, which is clearly in the DGLAP regime, is always taken to be on-mass-shell. matrix element for the hard $`ggc\overline{c}`$ subprocess amplitude . In this way the result embodies the main part of the NLO DGLAP correction. Besides this, the unified equation includes a kinematic constraint (or consistency condition)  which accounts for the major part of the higher-order corrections to the BFKL evolution.555Explicit expressions for the next-to-leading log BFKL terms can be found in . Indeed, the power behaviour generated corresponds to $`\lambda 0.3`$ which is much less than the LO BFKL behaviour $`x^{\omega _0}`$ with $`\omega _0=12\alpha _S\mathrm{ln}2/\pi `$. Moreover, the charm production cross section $`d\sigma /dx_F`$ calculated in terms of the KMS unintegrated gluons  is found to coincide, within 10% accuracy, with the prediction obtained with conventional DGLAP gluons  for $`10^4`$$`10^5`$ GeV laboratory energies, corresponding to $`x=10^3`$$`10^4`$ for which deep inelastic accelerator data exist. Note that, in this comparison, the prediction based on conventional partons was calculated at LO and multiplied by $`K=2.3`$, which, according to Ref. , accounts well for the NLO corrections. On the other hand the unified DGLAP/BFKL prescription already incorporates the main NLO effect at small $`x`$, and so no $`K`$ factor is present in this $`x`$ domain. We denote the results obtained using this gluon by KMS on the Figures below.

At $`10^5`$ GeV these predictions should be reliable. However, as we proceed to higher energies, we sample gluons in smaller and smaller $`x`$ regimes with increasing gluon density, and so we must account for the effect of saturation. To study this effect we start with the Golec-Biernat, Wüsthoff (GBW) model of deep inelastic scattering (DIS)  (and diffractive DIS ). Let us outline the basis of the model, as applied to $`q\overline{q}`$ production in DIS. The production of $`q\overline{q}`$ pairs is described by the probability of the formation of the pair by the initial photon multiplied by the cross section for the $`q\overline{q}`$–proton interaction, $`\widehat{\sigma }`$. The first stage is given by the effective photon wave functions $`\mathrm{\Psi }_{T,L}`$ (for transverse, longitudinal polarisations), which depend on the momentum fraction $`z`$ carried by the quark and the transverse separation $`\stackrel{}{r}`$ between the $`q`$ and $`\overline{q}`$. The deep inelastic cross sections have the form

$$\sigma _{T,L}(x,Q^2)=d^2r_0^1𝑑z\underset{q}{}\left|\mathrm{\Psi }_{T,L}^q(\stackrel{}{r},z,Q^2)\right|^2\widehat{\sigma }(x,\stackrel{}{r}).$$
(3)
For small $`r`$, the dipole cross section $`\widehat{\sigma }`$ is proportional to $`r^2`$. To allow for multiple interactions, Golec-Biernat and Wüsthoff  parameterize $`\widehat{\sigma }`$ by the form

$$\widehat{\sigma }(x,r)=\sigma _0\left(1\mathrm{exp}\left(\frac{r^2}{4R_0^2(x)}\right)\right),$$
(4)
with an $`x`$-dependent saturation radius

$$R_0(x)=\frac{1}{Q_0}\left(\frac{x}{x_0}\right)^{\lambda /2}.$$
(5)
The parameterization is a simplified version of the well-known Glauber expression for, say, describing the multiple interactions of a pion passing through a nucleus

$$\sigma _{\pi A}=d^2b_t\left[1\mathrm{exp}\left(\sigma _{\pi N}T(b_t)\right)\right].$$
(6)
The integral $`d^2b_t`$ gives the nuclear area $`\pi R_A^2`$, which is replaced by $`\sigma _0`$ in (4), and the mean nucleon density $`T`$ is parameterized by $`1/R_0(x)^2`$, modulo normalisation. $`\sigma _{\pi N}`$ is the $`\pi `$–nucleon cross section, which is equivalent to the $`q\overline{q}`$-dipole cross section in the GBW model. That is, the exponent $`\sigma _{\pi N}T`$ in (6) is equivalent to $`(\pi ^2r^2\alpha _S/3)xg/\sigma _0`$, where the gluon density $`xg/\sigma _0`$ plays the role of the mean nucleon density $`T`$, and where the factor in brackets plays the role of $`\sigma _{\pi N}`$. It is because the gluon density grows as $`x`$ decreases that we have an $`x`$ dependence in the argument of the exponential in (4).

The GBW model has recently been realised  in terms of the gluon density, including the DGLAP $`\mathrm{ln}Q^2`$ evolution of $`g(x,Q^2)`$. Actually in this improved form it should be considered, not as a model, but as a complete perturbative calculation, which in addition to the conventional LO collinear approach also accounts for the rescattering of the incoming $`q\overline{q}`$ pair.

The power of $`x`$ in (5) reflects the power growth of the gluon density in the small $`x`$ region. The parameters $`\sigma _0,x_0`$ and $`\lambda `$ were fitted to describe the small $`x`$ DIS data . It was shown that, up to $`Q^220`$ GeV2, a good description can be achieved, even without accounting for DGLAP evolution. Interestingly, the value of the power, $`\lambda =0.28`$, turns out to be close to the value found in Ref. .

So far we have considered absorption for DIS. Here we are interested in $`ggc\overline{c}`$, and not $`\gamma gc\overline{c}`$. It is therefore necessary to multiply the $`\gamma gc\overline{c}`$ cross section by

$$\frac{1}{6}\left[1\frac{9}{4}z(1z)\right]\frac{\alpha _S}{\alpha e_c^2}$$
(7)
where the first factor is due to colour, and the second term in square brackets accounts for $`gggc\overline{c}`$ production (see, for example, ).

Note that the approach of Golec-Biernat and Wüsthoff includes only the rescattering of the $`c\overline{c}`$ pair and neglects the enhanced Reggeon diagrams which account for the rescattering of the more complicated Fock components of the photon (gluon) like $`c\overline{c}g`$, $`c\overline{c}gg`$, etc. These extra components have a larger absorptive cross section. In other words, when the gluon density becomes sufficiently large, we must allow for $`gg`$ recombination, which diminishes the rate of $`c\overline{c}`$ production. From this $`t`$ channel viewpoint, the absorption is driven by the triple-Pomeron interaction. With the help of the Balitskii–Kovchegov equation , we may sum up the resulting fan diagrams (formed from different networks of Pomeron–Pomeron recombinations into single Pomeron exchanges). This effect has been studied numerically in recent papers . However, the approach is based on LO BFKL, and does not account for the NLO corrections, which are known to be large. In this case we cannot simply rescale the LO prediction by taking a lower Pomeron intercept, $`\omega _0`$. The problem is that the triple-Pomeron vertex is not known at NLO.

There are reasons, both from phenomenology  and from perturbative evolution , to believe that the triple-Pomeron coupling is small. Nevertheless, at very high energy, we expect the resulting absorption to be stronger. From this point of view we may regard the prediction based on the GBW model as the upper limit for $`c\overline{c}`$ production. Later, for a more realistic estimate of the cross section for $`c\overline{c}`$ production, we take account of the triple-Pomeron vertex by replacing (5) by666Such a form is motivated by the results of Ref. . Recall that the combinatorial factor which corresponds to the fan diagrams is $`N!/N!1`$. Thus, contrary to the eikonal form, where we have $`_N(g)^N/N!=\mathrm{exp}(g)`$, here we deal with a geometrical series type of expression $`_N(g)^N=1/(1+g)`$. Therefore we choose form (8) with the constant $`c`$ as evaluated in .

$$R_0^2=\frac{1}{Q_0^2}\left(c+(x/x_0)^\lambda \right)$$
(8)
for $`x<10^310^4`$ with $`c0.05`$$`0.2`$. This is to protect $`R_0`$ becoming too small for small $`x`$.

The only reason why the above upper limit may be exceeded arises because the GBW saturation model  was formulated at fixed impact parameter, and so does not allow for the growth of the proton radius $`R_p`$ with increasing energy. The radius $`R_p`$ can be determined from the slope $`B`$ of the elastic $`pp`$ cross section,

$$B=B_0+2\alpha ^{}\mathrm{ln}E,$$
(9)
where $`\alpha ^{}`$ is the slope of the Pomeron trajectory, and $`E`$ is the proton energy in the laboratory frame. Indeed, for a large-size dipole the GBW model saturates at $`\sigma =\sigma _0=29`$ mb, whereas the normal soft hadronic cross sections, which should be equivalent to large-size dipoles, continue to grow logarithmically with energy. From a physical point of view, the normalisation $`\sigma _0`$ in (4) is related to the proton area $`\pi R_p^2`$. Of course, we only have the inequality $`\sigma _0<\pi R_p^2B`$, since charm production originates mainly from the centre of the proton. However, since $`\pi R_p^2`$ grows with energy, a conservative upper limit is obtained by multiplying the prediction obtained from the GBW model by the factor $`B(E)/B(E_0)`$, with $`E_010^4`$ GeV, typical of the HERA domain where the parameters of the model were tuned.

A lower limit to $`c\overline{c}`$ production may be obtained if we assume a scaling behaviour for $`dn_c/dx_F=(d\sigma (ppc+\mathrm{})/dx)/\sigma _{\mathrm{inel}}`$, where $`\sigma _{\mathrm{inel}}`$ is the total inelastic $`pp`$ cross section – that is if we assume that $`dn_c/dx_F`$ is independent of energy. Hence the lower limit is

$$\frac{d\sigma (E)}{dx_F}=\frac{d\sigma (E_0)}{dx_F}\frac{\sigma _{\mathrm{inel}}(E)}{\sigma _{\mathrm{inel}}(E_0)},$$
(10)
normalised in the region $`E10^5`$ GeV ($`x10^4`$) where the parton distributions are known. To be more precise we should replace $`\sigma _{\mathrm{inel}}`$ in (10) by the cross section corresponding to the Fock component of the proton wave function which contains charm. However, the cross section for this component will grow with energy faster than $`\sigma _{\mathrm{inel}}`$, and thus (10) may be regarded as an extreme lower limit for the charm yield. We consider the Fock charm component to be generated perturbatively. In principle, it would be possible to have also a non-perturbative ‘intrinsic’ charm component , although there is no firm experimental evidence for its existence. Such an intrinsic charm cross section would originate from the non-perturbative large-size domain, controlled by $`\sigma _{\mathrm{inel}}`$, and hence its contribution would become less important, with increasing energy, than the perturbative cross section.

## 3 Predictions for high energy $`c\overline{c}`$ production

In Fig. 2 we compare the predictions for the $`x_F`$ distribution of charm quarks produced in $`pp`$ collisions, as given by the three models for extrapolating the gluon777The gluons in (1) are evaluated at a scale $`\mu _F`$ equal to the transverse mass of the charm quark for the MRST and KMS models; that is $`\mu _F^2=m_c^2+p_{cT}^2`$. For the GBW extrapolation we take $`\mu _F=1/r`$, where $`r`$ is the separation between the $`c`$ and $`\overline{c}`$ quarks. to small $`x`$. For laboratory energies $`E10^3`$$`10^5`$ GeV we sample the gluons in the $`x`$ region $`10^2`$$`10^4`$ where the parton distributions are known from global analyses. Hence, since each model reproduces the same data, they give essentially the same predictions for $`c\overline{c}`$ production. Recall that the LO DGLAP result, based on MRST partons, was multiplied by a $`K`$ factor of 2.3. It was shown  that such a constant $`K`$ factor reproduces well the NLO perturbative QCD prediction and gives a good description of the available fixed-target data for $`c\overline{c}`$, or rather $`D`$ meson, production for $`E250`$ GeV. Recall that, following , we take the mass of the charm quark to be $`m_c=1.25\mathrm{GeV}`$. Although we use a constant $`K`$ factor, $`K=2.3`$, we have confirmed that the use of the parameterization of the $`K`$ factor, $`K(E_c,x_c)`$, given in eq. (3.4) of , does not appreciably alter any predictions.

Up to $`E10^7`$ GeV, the GBW saturation model practically coincides with the DGLAP (MRST) prediction. For higher energies the GBW cross section is lower due to absorptive effects. On the contrary, the prediction based on KMS partons becomes higher, as well as lower at the lower energies. The ‘unified’ KMS evolution includes the BFKL $`\mathrm{ln}1/x`$ resummation and generates a power growth, $`x^\lambda `$, of the gluon density as we extrapolate to small $`x`$. This evolution embodies a kinematic constraint (or consistency condition) which accounts for a major part of the NLO and higher-order BFKL effects. However, the power, $`\lambda 0.3`$, is appreciable, and the growth exceeds the double logarithmic DGLAP growth of the MRST result. Another feature to note is that the shape of the $`x_F`$ distribution becomes a little steeper with increasing energy, as seen in Fig. 2(d), which shows the predictions of the GBW extrapolation for three different energies.

On the other hand, at low energies, we see from Fig. 3 that the KMS prediction falls away. Indeed, it is about a factor of two below the GBW/MRST predictions for $`E_c=10^3`$ GeV. Here, we are sampling the gluon at $`x`$ values above 0.01, where the $`\mathrm{ln}1/x`$ resummation is not effective and where LO DGLAP evolution takes over. To put it another way, it was observed in the KMS unified BFKL/DGLAP approach , that the cross section for an off-shell gluon is enhanced by $`\mathrm{ln}1/x`$ effects for $`x\stackrel{<}{}\mathrm{\hspace{0.33em}0.01}`$, whereas it rapidly tends to the LO on-shell DGLAP formula as $`x`$ increases above this value. Hence in the low energy regime a factor $`K2`$ should be included in the KMS prediction, with the factor dying away with increasing energy, as we enter the BFKL regime.

A convenient way to summarize the relevant energy behaviour of the $`d\sigma ^c/dx_F`$ cross section is to plot the ‘$`Z`$ moment’  of the $`x_F`$ distribution, see for example . For high energies ($`E>10^6`$ GeV) the incoming primary cosmic ray flux falls down as $`E^{(\gamma +1)}`$ with $`\gamma =2.02`$. Therefore the charm flux is proportional to the moment

$$\sigma Z_c\frac{d\sigma ^c}{dx}x^{2.02}𝑑x.$$
(11)
This moment is shown in Fig. 3,

where the difference between the saturation model and the other two models becomes apparent for $`E_c>10^6`$ GeV. Note that here we fix the energy $`E_c`$ of the outgoing $`c`$ quark, rather than that of the incident proton which was used in Fig. 2.

Although the GBW model predicts the smallest cross section of the three models, it should be considered as the upper limit888Modulo a possible growth of the proton radius, see (9). for $`c\overline{c}`$ production as it only accounts for part of the absorptive effects. For the reasons mentioned earlier, the GBW model is more than a model – rather it should be regarded as a full leading-order QCD prediction with $`c\overline{c}`$ rescattering effects included. Of course there is, in addition, absorption of the gluons in the evolution process. It appears likely that the consequent reduction of the $`c\overline{c}`$ cross section due to this additional absorption may be partially compensated by the growth of the proton radius with increasing energy. We investigate this, and other effects, in the next Section; see Fig. 4 later. Therefore, from now on, we will base our study on the GBW approach and its variations.

## 4 $`c\overline{c}`$ production in proton–air collisions

For a precise study we need the charm yield in $`p`$-nuclear (air) collisions. An advantage of the GBW saturation model is that it may be straightforwardly extended from $`pp`$ to $`p`$-nuclear collisions. Within the GBW framework we may account for the eikonal rescattering of the $`c\overline{c}`$ pair within the nucleus by replacing $`\widehat{\sigma }`$ in (4) by

$$A\widehat{\sigma }\underset{n=1}{\overset{\mathrm{}}{}}\frac{((A1)\widehat{\sigma }/8\pi B_A)^{n1}}{nn!}$$
(12)
where, for air, the mean atomic number $`A=14.5`$ and the slope $`B_A=r^2/6=29\mathrm{GeV}^2`$. Note that in the numerator we have $`(A1)`$ and not $`A`$. In this way we exclude rescattering on the nucleon where the $`c\overline{c}`$ pair is created. This rescattering is controlled by a different slope ($`B_A`$), and is already accounted for in the cross section given in (4). We have taken the root-mean-square nuclear radius999For example, Ref.  gives the r.m.s. radius of oxygen as 2.7 fm. To obtain the result for air we take the usual $`rA^{\frac{1}{3}}`$ dependence. $`\sqrt{r^2}=2.6`$ fm, and assumed a gaussian distribution of the nucleon density for light nuclei. Note that the replacement of $`\widehat{\sigma }`$ occurs before the integration over the $`c`$$`\overline{c}`$ separation in (3). In summary, the inclusive $`c\overline{c}`$ cross section for proton–air interactions is given by the sum of proton–nucleon cross sections, with the only nuclear effect being the enhanced absorption of the produced $`c\overline{c}`$ pair. In the Appendix we give a simple parameterization which reproduces the proton–air $`c\overline{c}`$ production cross section in the energy interval $`10^4<E<10^{11}\mathrm{GeV}`$ to within $`\pm 5\%`$.

The $`\gamma =2.02`$ moment of charm production in high energy proton–air collisions is shown in Fig. 4.

The dashed curve is the prediction of the extrapolation based on the original GBW saturation model, while the upper dotted curve (marked by $`B`$) includes the possible effect of the growth of the proton radius with energy, as discussed at the end of Section 2, see in particular (9). Instead of using the Pomeron slope, $`\alpha ^{}=0.25\mathrm{GeV}^2`$, measured in elastic $`pp`$ scattering, in (9) we use the value $`0.11\mathrm{GeV}^2`$ together with slope $`B=4.4\mathrm{GeV}^2`$ at $`W=90\mathrm{GeV}`$, which were deduced from elastic $`J/\psi `$ production data at HERA . These values are appropriate for charm production. The lower dotted line (marked $`g_{3P}`$) reflects the inclusion of an additional absorptive effect—the absorption of gluons. We take $`c=0.2`$ in (8), which corresponds to the largest estimate101010Triple-Pomeron phenomenology gives $`g_{3P}0.5\mathrm{GeV}^1`$. However, the bare triple-Pomeron vertex may be a few times larger, since the phenomenological value already accounts for some screening effects. It was argued in Ref.  that in perturbative QCD we expect $`g_{3P}=(0.5`$$`2)\mathrm{GeV}^1`$., $`g_{3P}=2\mathrm{GeV}^1`$, of the triple-Pomeron vertex. This choice is made to show the full extent of the uncertainty in $`c\overline{c}`$ production. If both the above effects are included (the radius growth and $`g_{3P}`$), then the solid curve (marked $`B+g_{3P}`$) is obtained. Conservatively, we predict that $`Z_c`$ will lie within the shaded region in Fig. 4; the most likely behaviour is that it will follow the cross-hatched region. Even the conservative prediction has much less uncertainty than previous estimates. However, for completeness, the dot-dashed curve shows the extreme lower limit, given by scaling formula (10), but where now $`\sigma _{\mathrm{inel}}`$ is the proton–air cross section. All the variations of the original GBW model were normalized in the region $`E10^5\mathrm{GeV}`$, where the partons sampled in the hard subprocess are known, and the model tuned to the data.

## 5 Prompt neutrinos: development of the air shower

Our aim is to predict the spectra of prompt $`\nu _\mu `$ and $`\nu _\tau `$ neutrinos produced in the atmosphere by cosmic rays. Prompt leptons originate from the following sequence: the production of $`c\overline{c}`$ pairs, their fragmentation into charm hadrons which then undergo semi-leptonic decay. In the lower energy range ($`E<10^7\mathrm{GeV}`$) it is possible to estimate the leptonic spectra by simply taking a product of moments of the various distributions, see, for example, Ref. . However, for $`E\stackrel{>}{}\mathrm{\hspace{0.33em}10}^7\mathrm{GeV}`$ the decay length of $`D`$ mesons becomes comparable with the depth of the atmosphere, and so it is necessary to follow the development of the air shower in more detail. It is described by a set of equations in terms of the ‘depth’ $`X`$ of the atmosphere transversed by a particle in the vertical direction. $`X`$ is related to the height $`h`$ by

$$X=_h^{\mathrm{}}\rho (h^{})𝑑h^{}$$
(13)
where $`\rho (h)`$ is the density of the atmosphere at vertical height $`h`$. We take the same exponential profile of the atmosphere111111We have checked that the numerical results essentially do not depend on the precise parameterization of the density of the atmosphere. as was used in Ref. . The sequence of equations determine $`\varphi _a(E,X)`$, which are the fluxes of the corresponding particles with energy $`E`$ at depth $`X`$, where $`a=N,c,i,l`$ (that is nucleon, $`c`$ quark, charmed hadron, lepton). The initial flux $`\varphi _N(E,0)`$ is the known primary cosmic ray flux. All other initial fluxes are zero, that is $`\varphi _a(E,0)=0`$ for $`a=c,i,l`$. The set of equations which determine $`\varphi _N\varphi _c\varphi _i\varphi _l`$ are

$$\frac{\varphi _N(E,X)}{X}=\frac{1}{\mathrm{\Lambda }_N(E)}\varphi _N(E,X);$$
(14)
$$\varphi _c(E,X)=_E^{\mathrm{}}𝑑E^{}𝑑x_c\varphi _N(E^{},X)\frac{n_A}{\rho }\frac{d\sigma }{dx_c}^{pc}\delta (Ex_cE^{});$$
(15)
$$\frac{\varphi _i(E,X)}{X}=\frac{1}{\mathrm{\Lambda }_i(E,X)}\varphi _i(E,X)+_E^{\mathrm{}}𝑑E^{}𝑑x\varphi _c(E^{},X)\frac{dn}{dx}^{ci}\delta (ExE^{})$$
(16)
with $`i=D^\pm ,D^0,\overline{D}^0,D_s^\pm ,\mathrm{\Lambda }_c`$; and finally

$$\frac{\varphi _l(E,X)}{X}=\underset{i}{}_E^{\mathrm{}}𝑑E^{}𝑑x\varphi _i(E^{},X)\frac{1}{\lambda _i^{\mathrm{dec}}(E^{},X)}B(il)\frac{dn}{dx}^{il}\delta (ExE^{})$$
(17)
where $`B(il)`$ is the branching fraction of the decay of the charmed hadron $`i`$ to lepton $`l`$. The nucleon attenuation length is

$$\mathrm{\Lambda }_N(E)\frac{\lambda _N(E)}{1Z_{NN}},$$
(18)
where $`Z_{NN}`$ is the spectrum-weighted moment for nucleon regeneration and $`\lambda _N`$ is the interaction thickness

$$\lambda _N=\rho (h)/\underset{A}{}\sigma _{NA}(E)n_A(h).$$
(19)
$`n_A(h)`$ is the number density of air nuclei of atomic number $`A`$ at height $`h`$ and $`\sigma _{NA}`$ is total $`NA`$ inelastic cross section. Instead of the sum over $`A`$ in (19), we take the mean value $`A=14.5`$ for air. Note that the factor $`n_A/\rho `$ in (15) arises from (19). For the nucleon–air cross section, $`\sigma _{N\mathrm{air}}`$, we take the parameterization of Bugaev et al. . For the incoming cosmic ray flux we take the parametrisation given in  denoted as TIG with knee. Also from  we take parametrisation of $`Z_{NN}`$ which depends on energy and takes into account the knee, which is consistent with . This $`Z`$ factor includes the regeneration by the $`p,n,N^{}\mathrm{}`$ particles.

The attenuation length $`\mathrm{\Lambda }_i(E,X)`$ of the charmed hadrons consists of two parts: the decay length $`\lambda _i^{\mathrm{dec}}`$ and attenuation due to their strong interactions with air nuclei. The decay lengths, $`\lambda _i^{\mathrm{dec}}`$, of the various charmed hadrons depend, via their $`X`$ dependence, on the density of the atmosphere,

$$\lambda _i^{\mathrm{dec}}=c\tau _i\frac{E}{m_i}\rho (X),$$
(20)
where $`m_i`$ and $`\tau _i`$, the mass and lifetime of the $`i^{\mathrm{th}}`$ charmed hadron, are taken from . The attenuation due to the strong interactions of the produced charmed hadrons with the air has a form similar to (18), namely $`\lambda _i/(1Z_{c\overline{c}})`$. To calculate $`\lambda _i`$ we assume an absorptive cross section equal to half the absorptive $`p`$–air cross section (based on quark counting), and we take a charm regeneration factor $`Z_{c\overline{c}}=(0.8)^\gamma `$. That is, we estimate that the leading charm quark carries a fraction $`x=m_c/(m_c+m_q)0.8`$ of the incoming energy, where $`m_q0.3\mathrm{GeV}`$ is the mass of a light constituent quark. For simplicity, we take the same $`Z_{c\overline{c}}`$ and $`\lambda _i`$ for all charm hadrons; for $`\mathrm{\Lambda }_c`$ we expect the larger cross section to be approximately compensated by the larger $`Z_{c\overline{c}}`$. Thus, finally, $`\mathrm{\Lambda }_i`$ is given by

$$\frac{1}{\mathrm{\Lambda }_i}=\frac{1}{\lambda _i^{\mathrm{dec}}}+\frac{1Z_{c\overline{c}}}{\lambda _i}.$$
(21)
From (14) the light baryon flux is given by

$$\varphi _N(E,X)=\varphi _N(E,0)\mathrm{exp}(X/\mathrm{\Lambda }_N(E)).$$
(22)
If we insert (15) into (16), then the individual charm hadron fluxes are

$$\varphi _i(E,X)=_0^X𝑑X^{}\mathrm{exp}\left(_X^{}^X\frac{dX^{\prime \prime }}{\mathrm{\Lambda }_i(E,X^{\prime \prime })}\right)S_{Ni}(E,X^{})$$
(23)
where

$$S_{Ni}(E,X)=_E^{\mathrm{}}𝑑E^{}\frac{n_A}{\rho }\frac{1}{E^{}}\frac{d\sigma }{dx}^{Ni}(E,E^{})\varphi _N(E^{},X)$$
(24)
and $`x=E/E^{}`$.

The spectra of charmed hadrons $`d\sigma ^{Ni}/dx`$ are calculated using the three models shown in Fig. 2. The ratios of the different charm yields (after the hadronization of the $`c`$ quark) are given in Ref.  to be

$$\frac{\sigma (D_s^+)}{\sigma (D^+,D^0)}=0.2,\frac{\sigma (\mathrm{\Lambda }_c)}{\sigma (D^+,D^0)}=0.3,\frac{\sigma (D^+)}{\sigma (D^0)}=0.5.$$
(25)
More recent data  favour a smaller value121212For isolated $`c`$-quark fragmentation, PYTHIA  gives a ratio of 0.08. of the first ratio quoted in (25), namely

$$\frac{\sigma (D_s^+)}{\sigma (D^+,D^0)}0.1.$$
(26)
We take this value in our analysis.

Note that the $`\mathrm{\Lambda }_c`$ baryon is produced by the recombination of a $`c`$ quark with a spectator diquark of the incoming nucleon131313The PYTHIA/JETSET Monte Carlo gives only 2.5% $`\mathrm{\Lambda }_c`$ baryons in the fragmentation of an isolated $`c`$ jet, that is in the absence of spectator diquarks.. It is not produced from a $`\overline{c}`$ quark. Therefore the parton momentum fraction $`x_\mathrm{\Lambda }`$ carried by the $`\mathrm{\Lambda }_c`$ is $`x_\mathrm{\Lambda }=x_d+x_c`$. The diquark momentum fraction $`x_d`$ will be less than $`\frac{2}{3}(1x_c)`$, as part of the energy is carried away by the third valence quark (the factor $`\frac{2}{3}`$), by the $`\overline{c}`$ quark and by gluons. To allow for this, we therefore take $`x_d=\frac{1}{2}(1x_c)`$, which leads to

$$x_\mathrm{\Lambda }=\frac{1}{2}(1+x_c)$$
(27)
This is found to be in good agreement with the distribution generated by the PYTHIA Monte Carlo , which has a maximum in the region $`x_\mathrm{\Lambda }0.5`$$`0.6`$. Of course a very slow $`c`$ quark is unlikely to combine with a fast light diquark (with $`x_d0.5`$). Therefore we introduce a cut-off, $`x_c>x_c^0`$, in (27). Assuming a mean velocity of the $`c`$ quark to be $`v^20.25`$, we estimate $`x_c^00.1`$. For $`\mathrm{\Lambda }_b`$ production, which we discuss below, the heavier $`b`$ quark will carry a larger fraction of the $`\mathrm{\Lambda }_b`$ momentum. In this case we take the cut-off to be $`x_b^0=0.25`$. As a check, we also compute the prompt flux, arising from $`c\mathrm{\Lambda }_c`$, using an alternative to prescription (27). We assume that for the diquark $`x_d=(m_d/m_c)/x_c`$ with a constituent diquark of mass $`m_d=2m_q=0.7`$ GeV and $`m_c=1.5`$ GeV. With this assumption

$$x_\mathrm{\Lambda }=1.47x_c,$$
(28)
or, in the case of beauty, $`x_\mathrm{\Lambda }=1.16x_b`$, using $`m_b=4.5`$ GeV.

Rather than using a fragmentation function for $`dn^{cD}/dx`$, for $`D`$ mesons we take $`x_D=0.75x_c`$. This is sufficient for our purposes141414That is the distributions $`dn^{ci}/dx`$ were taken to be proportional to $`\delta (x_D0.75x_c)`$ for $`D`$ mesons, whereas for $`\mathrm{\Lambda }_c`$ we assume that they are proportional to $`\delta (x_\mathrm{\Lambda }\frac{1}{2}(x_c+1))`$ for $`x_c>0.1`$, or $`\delta (x_\mathrm{\Lambda }1.47x_c)`$. We note that PYTHIA gives a harder $`x=E_D/E_c`$ distribution than that shown for the Peterson et al. function  in Ref. .. For illustration we compare, in Fig. 5, the prompt $`\nu _\mu +\overline{\nu }_\mu `$ flux at ground level, $`\varphi _{\nu _\mu +\overline{\nu }_\mu }(E)`$, obtained from (17) using different forms of fragmentation. Clearly the upper curve, corresponding to no fragmentation, gives an overestimate of the flux. Moreover, due to the presence of additional light “sea” quarks, we expect a harder distribution for the fragmentation in $`pp`$ production than that obtained in $`e^+e^{}`$ collisions (lower curve). Hence our use of fragmentation corresponding to one of the middle curves. We see that both the $`\mathrm{\Lambda }_c`$ hadronization prescriptions (27) and (28) give very similar fluxes. At the highest energies the fraction of neutrinos coming from $`\mathrm{\Lambda }_c`$ (relative to those from $`D`$) increases due to the short $`\mathrm{\Lambda }_c`$ lifetime. Therefore the choice $`x_\mathrm{\Lambda }=\frac{1}{2}(1+x_c)`$ of (27), which corresponds to a larger $`x_\mathrm{\Lambda }`$ than (28), gives a larger neutrino flux for $`E_\nu \stackrel{>}{}\mathrm{\hspace{0.33em}10}^6`$ GeV. The results below correspond to using prescription (27).

For the leptonic decay of each charm hadron $`i`$, the distribution $`dn^{il}/dx_l`$ was generated by the PYTHIA Monte Carlo . Note that in PYTHIA 6.2 the $`D_s\tau \nu _\tau `$ branching ratio was set to be 1%, whereas the latest value is $`6.4\pm 1.5\%`$ . Since $`D_s\tau \nu _\tau `$ is almost the only source of prompt $`\nu _\tau `$ neutrinos, it is important to renormalise the yield using the updated branching ratio. It is interesting to note that the $`D_s^+`$ decay produces more prompt $`\overline{\nu }_\tau `$ neutrinos than $`\nu _\tau `$, since the $`\overline{\nu }_\tau `$ from $`\tau ^+`$ decay has the large $`x_l`$. Of course, the reverse is true for $`D_s^{}`$ decay.

## 6 Prompt $`\nu _\mu `$ and $`\nu _\tau `$ fluxes

We present the predicted yields of prompt $`\nu _\mu `$ and $`\nu _\tau `$ neutrinos from charmed hadrons produced by cosmic ray collisions in the atmosphere. The flux of $`\nu _e`$ neutrinos is essentially equal to that of $`\nu _\mu `$ neutrinos. High energy prompt electrons are completely degraded in the electromagnetic cascade. For prompt muons the electromagnetic interaction is much weaker; it was demonstrated in Ref.  that the prompt $`\mu `$ flux is only about 10% smaller than the prompt $`\nu _\mu `$ flux at the surface of the earth.

In Fig. 6 we plot the prompt $`\nu _\mu `$ and $`\nu _\tau `$ fluxes predicted by the three models for the extrapolation of the gluon distribution to very small $`x`$. Although in Figs. 3 and 4 we compared different models of extrapolation using a fixed $`\gamma =2.02`$, in the actual calculation of the neutrino fluxes we used the observed primary cosmic ray flux, which corresponds to different values of $`\gamma `$ above and below the ‘knee’. The sharp fall-off for neutrino energies $`E>10^8\mathrm{GeV}`$ is due to the increase in the decay length of the charmed hadrons, arising from Lorentz time dilation. Clearly in the gluon kinematic domain of interest ($`x<10^5`$ and $`Q^2\stackrel{<}{}\mathrm{\hspace{0.33em}10}\mathrm{GeV}^2`$) saturation effects become important. The reason for the behaviour of the KMS prediction—below at low energies, above at high energies—was given in the discussion concerning Fig. 3. We have argued that extrapolations based on the GBW model and its variations, as shown in Fig. 4, give the most reliable predictions for $`E>10^6\mathrm{GeV}`$. Thus in Fig. 7 we show the spread of predictions of the neutrino fluxes based on the shaded domain in Fig. 4. In Fig. 7 we also show the conventional atmospheric flux of $`\nu _\mu `$ (from $`\pi `$ and $`K`$ decays, etc.). Moreover, there is a small probability that atmospheric $`\nu _\mu `$ neutrinos may oscillate into $`\nu _\tau `$ neutrinos and so provide an ‘atmospheric’ $`\nu _\tau `$ flux. We calculate this flux using the $`3\sigma `$ ranges of the $`\mathrm{sin}^2\theta _{\mathrm{ATM}}`$ and $`\mathrm{\Delta }m_{\mathrm{ATM}}^2`$ neutrino mixing parameters found in an analysis  of the Super-Kamiokande  and MACRO  data. The resulting atmospheric $`\nu _\tau `$ flux is shown by a band in the lower plot of Fig. 7.

We discuss the neutrino fluxes of Fig. 7 in the concluding section, after we have included the contribution to the prompt $`\nu _\tau `$ spectrum arising from $`b\overline{b}`$ production, fragmentation and decay. However, first, we show in Fig. 8 the effect of charmed hadron interactions with the atmosphere.

The effect is illustrated by the difference between the dashed curve, for which the interactions are suppressed (that is $`\lambda _i=0`$), and the continuous curve with the interactions present. For this comparison we use the GBW extrapolation, the one given by the solid curve in Fig. 4.

We emphasize that the predictions for the prompt neutrino fluxes depend strongly on the nuclear content of the primary cosmic rays. So far we have assumed that the cosmic rays are composed entirely of protons. Suppose that the protons were to be replaced by nuclei of the same energy $`E`$ and of atomic number $`A`$. Then we have to scale the energy of the primary interacting nucleon to $`E/A`$. Roughly speaking, this reduces the neutrino flux $`\varphi _\nu `$ in the plateau region ($`10^6`$$`10^7\mathrm{GeV}`$) of the $`E^3\varphi _\nu `$ plot by $`A^3`$. On the other hand, the number of incoming nucleons is $`A`$, so as a consequence we expect an $`A^2`$ suppression. A detailed calculation for incoming nuclei with $`A=7`$  gives the dotted curve in Fig. 8. The suppression of the original GBW continuous curve (for $`A=1`$) is apparent.

## 7 Prompt $`\nu _\tau `$ flux from $`b\overline{b}`$ production

At first sight it appears that if we also allow for $`b\overline{b}`$ production, with the same cross section, fragmentation and decay then we will approximately double the $`\nu _\tau `$ flux. However note that, in comparison with charm production, the cross section for beauty production is about 30 times smaller due to the factor $`m_c^2/m_b^2`$ and to the larger value of $`x_2`$ of the gluon structure function which is sampled, see eq. (1), see Fig. 9. Nevertheless, the high energy production of $`\nu _\tau `$ neutrinos from beauty decays is not negligible.

For charm, only the decay of the $`D_s`$ may produce $`\nu _\tau `$ neutrinos, whereas now $`B^\pm `$, $`B^0`$, $`B_s`$ and $`\mathrm{\Lambda }_b`$ semileptonic decays also give rise to a significant $`\nu _\tau `$ flux. Indeed including $`B`$ and $`\mathrm{\Lambda }_b`$ decays enlarges the predicted total prompt $`\nu _\tau `$ flux by about 40% for $`E_\nu 10^5\mathrm{GeV}`$, and even more at higher energies, see Fig. 10(b). This is considerably larger than the $`\nu _\tau `$ flux calculated in , where the beauty induced contribution grows from 1% at $`E_\nu =10^2`$ GeV to about 10% at $`E_\nu =10^6`$ GeV. The reasons why we obtain a larger fraction of $`\nu _\tau `$ neutrinos from beauty are, first, that we use (26), rather than (25), for $`cD_s`$ hadronization and, second, our cross section for charm production is more suppressed at high energies by absorptive corrections, than the more compact $`b\overline{b}`$ production process.

Fig. 10(a) shows the breakdown of the total prompt $`\nu _\tau `$ flux. We take the ratios of the $`B^\pm `$, $`B^0`$, $`B_s`$ and $`\mathrm{\Lambda }_b`$ beauty yields (after hadronization of the $`b`$ quark) to be given by exactly analogous relations to (25) and (26) for charm. Recall that $`\nu _\tau `$’s of $`c\overline{c}`$ and $`b\overline{b}`$ origin come, respectively, from $`D_s\tau \nu _\tau `$ and from the semileptonic $`\tau \nu _\tau `$ decay modes of $`B^\pm `$, $`B^0`$, $`B_s`$ and $`\mathrm{\Lambda }_b`$. Thus we have a direct $`\nu _\tau `$ component and an indirect component coming from the $`\tau \nu _\tau `$ decays. For charm the former component is very small since the direct $`\nu _\tau `$ carries away a small energy fraction. Thus the $`\tau \nu _\tau `$ decay gives the dominant component, which is truncated at an energy when the $`\tau `$ has insufficient time to decay. For beauty the direct and indirect components are comparable until the energy regime is entered at which the $`\tau `$ has not enough time to decay.

To calculate the $`b`$ quark cross section we have used the same GBW model as described in Sections 2–4 with a quark mass $`m_b=4.5\mathrm{GeV}`$. We have approximated the $`bB`$ fragmentation function by a delta function at $`x=E_B/E_b=0.9`$. Also we have assumed that the hadronization and the cross-section of the beauty hadron–air interaction do not depend appreciably on the nature of the heavy quark. So we have used the same attenuation lengths $`\mathrm{\Lambda }_i`$ and the same ratios of the different beauty hadron yields (eqs. (25) and (26)), as were used for the charm case.

We neglect the $`b\overline{b}`$ contribution to the $`\nu _\mu `$ (and $`\nu _e`$) fluxes. They never exceed about 2–3% of that arising from $`c\overline{c}`$ production and decay.

## 8 Concluding discussion

### 8.1 Implications for Physics

The $`\nu _\mu `$ flux shown in Fig. 7(a), and the $`\nu _\tau `$ flux shown in Fig. 10(b), have important implications for neutrino astronomy. These neutrino fluxes arise from cosmic ray interactions with the atmosphere. They therefore provide the background to searches for cosmic neutrinos (for which there are many exciting New Physics scenarios). They also provide a possibility to calibrate the detectors of the new neutrino telescopes. A particularly interesting energy domain is where the prompt neutrino flux has emerged from the sharply falling ‘conventional’ atmospheric neutrino flux. Here $`\nu _\tau `$ appears to have an enormous advantage in searching for a signal for New Physics, see Figs. 7 and 10. Due to mixing, the fluxes of $`\nu _\mu `$ and $`\nu _\tau `$ are equal for incoming cosmic neutrinos. The flux of prompt $`\nu _\tau `$ neutrinos (which arises mainly from $`D_s\tau \nu _\tau `$ decays, but with a significant 40% contribution from $`B`$, $`B_s`$, $`\mathrm{\Lambda }_b`$ semileptonic decays) is about ten times less than that for prompt $`\nu _\mu `$, and the number of $`\nu _\tau `$ neutrinos produced from the conventional atmospheric flux (via $`\nu _\mu \nu _\tau `$ oscillations) is negligible for $`E>10^4\mathrm{GeV}`$. There is only a little time for neutrino mixing in the atmosphere. Thus the tau neutrino flux originating from atmospheric neutrinos is greatly suppressed. As a consequence, tau neutrinos offer an ideal means of identifying neutrinos of cosmic origin, and for searching for New Physics.

There is another reason to concentrate on the $`\nu _\tau `$ flux. For values of $`E_\nu `$ above $`10^4`$ GeV the $`\nu _\mu `$ flux is significantly depleted in the passage of the neutrinos through the Earth. The absorption increases rapidly as $`E_\nu `$ is increased . On the other hand, high energy $`\nu _\tau `$ neutrinos have the unique advantage that they are not depleted in number no matter how much of the Earth that they pass through. If they suffer a charged-current interaction, they produce a $`\tau `$ lepton which subsequently decays, regenerating a $`\nu _\tau `$ neutrino with degraded energy . It is, moreover, interesting to note that tau neutrinos of energies about $`10^7`$ GeV may produce spectacular distinctive signatures in the planned 1 km3 high energy neutrino detectors made of strings of photomultiplier tubes distributed throughout a naturally occurring Cerenkov medium, such as water or ice, deep in the ocean or ice cap. A charged-current $`\nu _\tau `$ interaction may produce a contained ‘double bang’ signature  or a ‘lollipop’ signature . In the former, the first bang corresponds to the hadronic shower produced along with the subsequent $`\tau `$ lepton and the second bang is the shower associated with the $`\tau `$ decay. A lollipop event is when only the second shower occurs within the detector and the $`\tau `$ lepton which initiates the shower is identified by the relatively weak ionization that it causes.

Another interesting possibility follows from the observation that electron neutrinos may be distinguished experimentally from muon and tau neutrinos. Electron neutrinos give rise to distinctive showers, which result from their charged-current interactions, which characteristically are contained within the large detector volume, and hence identified. The conventional atmospheric flux of $`\nu _e`$ neutrinos comes mostly from kaon decays and so the relative flux $`\nu _e/\nu _\mu \stackrel{<}{}\mathrm{\hspace{0.33em}0.1}`$ . On the other hand, the prompt $`\nu _e`$ and $`\nu _\mu `$ fluxes are practically equal to each other. Hence the prompt $`\nu _e`$ flux should be more visible over its conventional background, than the prompt $`\nu _\tau `$ flux151515We thank John Beacom for drawing this possibility to our attention.. Of course, the background reduction is much less than for $`\nu _\tau `$ neutrinos, but this will be partially offset since it is likely that $`\nu _e`$ neutrinos will be easier to identify than $`\nu _\tau `$ neutrinos.

So far we have shown the prompt neutrino fluxes arising from vertically-incident cosmic rays. However, the predicted fluxes are sensitive to the zenith angle. In particular, as we depart from the vertical incident direction, more atmosphere is encountered, which gives more time for the charmed hadrons to decay and hence allows more high energy neutrinos to be produced. This is well illustrated by Fig. 11, which shows the prompt neutrino fluxes produced by horizontally-incident cosmic rays.

Indeed, for $`E_\nu >10^8\mathrm{GeV}`$, the flux in the horizontal direction is noticeably larger. The observation of the zenith angle dependence will allow the ‘atmospheric’ background neutrinos to be identified and hence cosmic neutrinos (and therefore New Physics) to be isolated.

### 8.2 Uncertainties in the predicted prompt $`\nu _\mu `$ and $`\nu _\tau `$ fluxes

We have emphasized the importance to neutrino astronomy of reliable predictions of the prompt muon and tau neutrino fluxes, which arise from $`c\overline{c}`$ and $`b\overline{b}`$ production, hadronization and decay in the atmosphere. They provide the background to the search for cosmic neutrinos originating from Active Galactic Nuclei or elsewhere. We have argued that the predictions based on perturbative QCD, given in Figs. 7(a) and 10(b), have much less uncertainty than those that already exist in the literature. How reliable is perturbative QCD and the extrapolation of the gluon into the very small $`x`$ domain? Perturbative QCD is expected to be valid for $`b\overline{b}`$ production161616Unlike $`c\overline{c}`$ production, there are no accelerator $`b\overline{b}`$ data in the fragmentation region which is relevant to this study. However, there exist Tevatron data for central $`b\overline{b}`$ production with, typically, $`p_{bT}>5`$ GeV  which appear to exceed the NLO QCD prediction. Recently it has been shown  that the discrepancy may be reduced to an acceptable level by using a fragmentation function with $`x`$ larger than that of the conventional Peterson et al. parameterization ., and should also hold for $`c\overline{c}`$ production since the factorization scale $`\mu _Fm_c`$. It is important to note that at fixed target energies, $`E250`$ GeV, the NLO predictions for $`d\sigma /dx_F`$ are in agreement with the available data, while for very high energies the mean scale, $`1/r`$, increases due to saturation effects. Since the gluon distribution is determined by HERA deep inelastic data down to $`x10^4`$, predictions, which use partonic structure determined from these data and which agree with fixed target charm data, should be reliable up to about $`E10^5`$ GeV, and perhaps an order of magnitude or so above. Fig. 3 shows examples of three very different models which illustrate that this is indeed the case. (Recall that in Section 3 we explained why the KMS prediction is not applicable at the lower energies shown in Fig. 3.) However, as we proceed to higher energies, and sample smaller and smaller $`x`$ values and hence increasing gluon density, we must include the effects of saturation. Though we have labelled the predictions by GBW, we have not simply used the saturation model of Golec-Biernat and Wüsthoff, but rather we investigated the uncertainties in including saturation effects in some detail, see Fig. 4. We concluded that the most reliable predictions would be obtained using the continuous curve in Fig. 4, with conservative errors shown by the shaded band, which are reflected in the shaded bands on the prompt neutrino fluxes shown in Fig. 7. These bands represent the uncertainty in $`c\overline{c}`$ production.

In addition, there are also uncertainties associated with the fragmentation of the charm quark (see Fig. 5) and with the $`D`$ meson attenuation due to its strong interaction with the atmosphere (see Fig. 8). Another source of uncertainty is the $`x`$ distribution of the $`\mathrm{\Lambda }_c`$ hyperons. The recombination of the $`c`$ quark with the $`ud`$ diquark of the incoming proton gives the $`\mathrm{\Lambda }_c`$ a rather large $`x`$. We use the simplified formula $`x_\mathrm{\Lambda }=\frac{1}{2}(1+x_c)`$, where a fast diquark with $`x_{ud}0.5`$ recombines with any $`c`$ quark with $`x_c>0.1`$. The relatively large uncertainty arises because the $`\nu `$ flux is proportional to the $`x^{\gamma =2.02}`$ moment, (11), of the $`x`$ distribution. Thus the larger value of $`x`$ carried by the $`\mathrm{\Lambda }_c`$ (as compared to the $`D`$ mesons) compensates for the small probability of $`\mathrm{\Lambda }_c`$ formation, see (25). As a consequence, the contribution of the $`\mathrm{\Lambda }_c`$ to the prompt $`\nu _\mu `$ flux increases from about 10% to 40% as the neutrino energy increase from $`10^4`$ to $`10^9`$ GeV. Unfortunately the $`x`$ distribution of $`\mathrm{\Lambda }_c`$ hyperons at high energies ($`E>10`$ TeV) is not measured yet. However, the $`x^{2.02}`$ moment given by PYTHIA is in agreement with the above prescription, (27), that $`x_\mathrm{\Lambda }=\frac{1}{2}(1+x_c)`$ for $`x_c>0.1`$.

Overall, we see that the uncertainty in the prompt neutrino flux predictions is about a factor of 3. For the $`\nu _\tau `$ flux from charm, which relies entirely on the $`D_s\tau \nu _\tau `$ decay, there is also an extra uncertainty associated with the branching ratio and with the production ratio of (26). There is less uncertainty in the $`\nu _\tau `$ flux from beauty, since it originates more uniformly from $`B^\pm `$, $`B^0`$, $`B_s`$ and $`\mathrm{\Lambda }_b`$ semileptonic $`\tau \nu _\tau `$ decays.

Finally we note that the predictions have been made assuming that the incoming cosmic rays are predominantly composed of protons. If, however, it should happen that the observed $`\nu _\mu `$ and $`\nu _\tau `$ prompt fluxes lie below the predictions of Figs. 7(a) and 10(b), then the measurements may give important constraints on the nuclear composition of cosmic rays.

## Acknowledgements

It is a special pleasure to acknowledge that over many years we have benefited from many valuable and illuminating discussions with Jan Kwieciński, on this and many other topics. We thank Francis Halzen for drawing the problem of prompt neutrinos to our attention. A.M.S. thanks Leonard Leśniak for discussions. This work was partially supported by the UK Particle Physics and Astronomy Research Council, the Russian Fund for Fundamental Research (grant 01-02-17095), the British Council–Polish KBN joint research programme, and Polish KBN grants 2P03B 05119, 5P03B 14420.

## Appendix: charm production in proton–air collisions

Here we present a simple parameterization which may be used to reproduce the inclusive cross section of $`c`$-quark production in proton–air collisions, based on the GBW model (which yields the continuous curves in Figs. 2–5). To $`\pm 5\%`$ accuracy for $`0.05<x<0.6`$ we have

$$x\frac{d\sigma }{dx}(p+\mathrm{air}c+\mathrm{})=Ax^\beta (1x^{1.2})^n$$
(29)
where $`\beta =0.050.016\mathrm{ln}(E/10\mathrm{TeV})`$ and

$$\begin{array}{cc}\hfill n=& 7.6+0.025\mathrm{ln}(E/10\mathrm{TeV})\hfill \\ \hfill A=& 140+(11\mathrm{ln}(E/0.1\mathrm{TeV}))^{1.65}\mu \mathrm{b}\hfill \end{array}\}\mathrm{for}10^4<E<10^8\mathrm{GeV},$$

$$\begin{array}{cc}\hfill n=& 7.6+0.012\mathrm{ln}(E/10\mathrm{TeV})\hfill \\ \hfill A=& 4100+245\mathrm{ln}(E/10^5\mathrm{TeV})\mu \mathrm{b}\hfill \end{array}\}\mathrm{for}10^8<E<10^{11}\mathrm{GeV}.$$

For $`E=10^{12}\mathrm{GeV}`$, the parameterization overestimates the cross section by about 10% for $`x<0.2`$.