Prompt neutrinos from atmospheric c-cbar and b-bbar production and the gluon at very small x

We improve the accuracy of the extrapolation of the gluon distribution of the proton to very small x, and show that the charm production cross section, needed to calculate the cosmic ray-induced `atmospheric' flux of ultrahigh energy prompt muon and tau neutrinos, may be predicted within perturbative QCD to within about a factor of three. We follow the sequence of interactions and decays in order to calculate the neutrino fluxes as a function of energy up to 10^9 GeV. We also compute the prompt neutrino tau flux from b-bbar production, hadronization and decay. New cosmic sources of neutrinos will be indicated if more prompt neutrinos are observed than predicted. If fewer neutrinos are observed than predicted, then constraints will be imposed on the nuclear composition of cosmic rays. The advantages of studying tau neutrinos are emphasized. We provide a simple parameterization of the prediction for the inclusive cross section for c quark production in high energy proton--air collisions.


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. [1,2]. This has led to the development of neutrino telescopes, which use photomultiplier 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. [2,3]. 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' neutrinos 1 . 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 [4], perturbative QCD gives a satisfactory description of the observed features of the available accelerator data on charm production, see for example [5]- [11]. Moreover, the simplified LO calculation reproduces the same behaviour when multiplied by an overall K factor, K ≃ 2.3, to account for the NLO contribution. Several authors have therefore used perturbative QCD to predict the prompt neutrino flux [7]- [11].
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 max L , where p L is the longitudinal momentum of the charm quark; at high energies x F ≃ E 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, where dσ/dz =s dσ/dt with z = (m 2 c −t)/s, and where g(x) is the gluon density of the proton. The Mandelstam variabless andt refer to the gg → cc subprocess. 2 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 2 ≃ M 2 cc /2x F s ∼ 10 −9 -10 −4 , where √ s is the total pp c.m. energy. There are no data which determine the gluon for x < ∼ 10 −4 , and, as a rule, parton distributions for x < 10 −5 are not available. For example, Ref. [11] 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 xg ∼ x −λ with λ in the range 0-0.5. Of course, the prediction depends crucially on the value taken for λ, 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: (i) DGLAP gluon with a double leading log (DLL) extrapolation, (ii) unified DGLAP/BFKL gluon [12] with x −λ extrapolation, (iii) 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 (E ∼ 10 5 to 10 9 GeV) using the three models for extrapolating the gluon into the small x regime, typically x ∼ 10 −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 cc pairs from the incoming cosmic ray flux, we allow for their fragmentation into charmed (D ± , D 0 , D s and Λ 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 τ lepton in the D s → τ ν τ decays. The results of the calculation of the prompt ν µ and ν τ fluxes are presented in Section 6, and compared with the 'conventional' atmospheric fluxes. These latter fluxes arise from π, K . . . decays and ν µ − ν τ oscillations respectively. We find that the prompt ν τ spectrum for E > 10 4 GeV lies much above its conventional atmospheric background, whereas for the prompt ν µ spectrum this is only achieved for E > 10 6 GeV. The origins of the prompt ν τ flux are the D s → τ ν τ decays which occur with a branching fraction of 6.4 ± 1.5% [13]. It is relevant to note that high energy ν τ 's, unlike ν µ 's, are not depleted in number by absorption in the Earth. They will always penetrate the Earth due to the ν τ → τ → ν τ . . . regeneration sequence [14]. 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 bb production will, relative to the flux of cc origin, be suppressed, first, by a factor of order m 2 c /m 2 b and, second, by the fact that the gluon is sampled at larger x. However, for the prompt ν τ flux of bb origin, the suppression is partly compensated by the existence of significant τ ν τ semileptonic decays of all the beauty hadrons (B ± , B 0 , B s and Λ b ), in contrast to just the D s → τ ν τ decays for charm. We calculate the prompt flux arising from bb 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 cc production As mentioned above, one possible method of extrapolating the gluon into the x < 10 −5 regime is to resum the leading α S ln Q 2 ln 1/x terms within the DGLAP framework, which leads to a small x behaviour 3 3 In practice, we use the DLL continuation to x < 10 −5 in the form where the LO coupling α S (Q) = 4π/(b log(Q 2 /Λ 2 QCD )) with n f = 4 and b = 25/3. We take Q 2 0 = 1 GeV 2 and x 0 = 0.25. We use MRST2001 partons [15] with Λ QCD = 220 MeV. The reliability of this form of DLL extrapolation has been checked using GRV partons [16] which are tabulated down to x = 10 −9 .
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 α S ln 1/x terms, which, at leading order, are resummed by the BFKL equation [17]. Therefore a more reliable extrapolation is obtained by solving a unified DGLAP/BFKL equation [12] 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 shell 4 matrix element for the hard gg → cc subprocess amplitude [18]. 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) [19] which accounts for the major part of the higher-order corrections to the BFKL evolution. 5 Indeed, the power behaviour generated corresponds to λ ≃ 0.3 which is much less than the LO BFKL behaviour x −ω 0 with ω 0 = 12α S ln 2/π. Moreover, the charm production cross section dσ/dx F calculated in terms of the KMS unintegrated gluons [12] is found to coincide, within 10% accuracy, with the prediction obtained with conventional DGLAP gluons [15] 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. [9], 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) [21] (and diffractive DIS [22]). Let us outline the basis of the model, as applied to qq production in DIS. The production of qq pairs is described by the probability of the formation of the pair by the initial photon multiplied by the cross section for the qq-proton interaction, σ. The first stage is given by the effective photon wave functions Ψ T,L (for transverse, longitudinal polarisations), which depend on the momentum fraction z carried by the quark and the transverse separation r between the q andq. The deep inelastic cross sections have the form [23] 4 That is we use the replacement in (1) where f g is the unintegrated gluon density as defined in [12], and where σ on, off are the gg → cc 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. 5 Explicit expressions for the next-to-leading log BFKL terms can be found in [20].
For small r, the dipole cross section σ is proportional to r 2 . To allow for multiple interactions, Golec-Biernat and Wüsthoff [21] parameterize σ by the form with an x-dependent saturation radius 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 The integral d 2 b t gives the nuclear area πR 2 A , which is replaced by σ 0 in (4), and the mean nucleon density T is parameterized by 1/R 0 (x) 2 , modulo normalisation. σ πN is the π-nucleon cross section, which is equivalent to the qq-dipole cross section in the GBW model. That is, the exponent σ πN T in (6) is equivalent to (π 2 r 2 α S /3)xg/σ 0 , where the gluon density xg/σ 0 plays the role of the mean nucleon density T , and where the factor in brackets plays the role of σ π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 [24] in terms of the gluon density, including the DGLAP 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 qq pair.
The power of x in (5) reflects the power growth of the gluon density in the small x region. The parameters σ 0 , x 0 and λ were fitted to describe the small x DIS data [21]. It was shown that, up to Q 2 ≃ 20 GeV 2 , a good description can be achieved, even without accounting for DGLAP evolution. Interestingly, the value of the power, λ = 0.28, turns out to be close to the value found in Ref. [12].
So far we have considered absorption for DIS. Here we are interested in gg → cc, and not γg → cc. It is therefore necessary to multiply the γg → cc cross section by where the first factor is due to colour, and the second term in square brackets accounts for gg → g → cc production (see, for example, [25]).
Note that the approach of Golec-Biernat and Wüsthoff includes only the rescattering of the cc pair and neglects the enhanced Reggeon diagrams which account for the rescattering of the more complicated Fock components of the photon (gluon) like ccg, ccgg, 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 cc production. From this t channel viewpoint, the absorption is driven by the triple-Pomeron interaction. With the help of the Balitskii-Kovchegov equation [26,27], 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 [28,29,30]. 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, ω 0 . The problem is that the triple-Pomeron vertex is not known at NLO.
There are reasons, both from phenomenology [31] and from perturbative evolution [32], 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 cc production. Later, for a more realistic estimate of the cross section for cc production, we take account of the triple-Pomeron vertex by replacing (5) by 6 The only reason why the above upper limit may be exceeded arises because the GBW saturation model [21] 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, where α ′ 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 σ = σ 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 σ 0 in (4) is related to the proton area πR 2 p . Of course, we only have the inequality σ 0 < πR 2 p ∝ B, since charm production originates mainly from the centre of the proton. However, since πR 2 p 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 0 ≃ 10 4 GeV, typical of the HERA domain where the parameters of the model were tuned.
A lower limit to cc production may be obtained if we assume a scaling behaviour for dn c /dx F = (dσ(pp → c + . . .)/dx)/σ inel , where σ 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 normalised in the region E ∼ 10 5 GeV (x ∼ 10 −4 ) where the parton distributions are known.
To be more precise we should replace σ 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 σ 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 [34], 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 σ inel , and hence its contribution would become less important, with increasing energy, than the perturbative cross section.

Predictions for high energy cc 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 gluon 7 to small x. For laboratory energies E ∼ 10 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 cc production. Recall that the LO DGLAP result, based on MRST partons, was multiplied by a K factor of 2.3. It was shown [9] 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 cc, or rather D meson, production for E ∼ 250 GeV. Recall that, following [9], we take the mass of the charm quark to be m c = 1.25 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 [8], does not appreciably alter any predictions.
Up to E ∼ 10 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 ln 1/x resummation and generates a power growth, x −λ , 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, λ ∼ 0.3, is appreciable, and the growth exceeds the double logarithmic DGLAP growth of the MRST result. Another feature 7 The gluons in (1) are evaluated at a scale µ F equal to the transverse mass of the charm quark for the MRST and KMS models; that is µ 2 F = m 2 c + p 2 cT . For the GBW extrapolation we take µ F = 1/r , where r is the separation between the c andc quarks.
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 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 [12], that the cross section for an off-shell gluon is enhanced by ln 1/x effects for x < ∼ 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 K ∼ 2 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σ c /dx F cross section is to plot the 'Z moment' [35] of the x F distribution, see for example [36]. For high energies (E > 10 6 GeV) the incoming primary cosmic ray flux falls down as E −(γ+1) with γ = 2.02. Therefore the charm flux is proportional to the moment 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 limit 8 for cc 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 cc 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 cc 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.

cc 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 cc pair within the nucleus by replacing σ in (4) by where, for air, the mean atomic number A = 14.5 and the slope B A = r 2 /6 = 29 GeV −2 .
Note that in the numerator we have (A − 1) and not A. In this way we exclude rescattering on the nucleon where the cc 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 rootmean-square nuclear radius 9 r 2 = 2.6 fm, and assumed a gaussian distribution of the nucleon density for light nuclei. Note that the replacement of σ occurs before the integration over the c-c separation in (3). In summary, the inclusive cc 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 cc pair. In the Appendix we give a simple parameterization which reproduces the proton-air cc production cross section in the energy interval 10 4 < E < 10 11 GeV to within ±5%.
The γ = 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, α ′ = 0.25 GeV −2 , measured in elastic pp scattering, in (9) we use the value 0.11 GeV −2 together with slope B = 4.4 GeV −2 at W = 90 GeV, which were deduced from elastic J/ψ production data at HERA [38]. 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 estimate 10 , g 3P = 2 GeV −1 , of the triple-Pomeron vertex. This choice is made to show the full extent of the uncertainty in cc 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 σ inel is the proton-air cross section. All the variations of the original GBW model were normalized in the region E ∼ 10 5 GeV, where the partons sampled in the hard subprocess are known, and the model tuned to the data.

Prompt neutrinos: development of the air shower
Our aim is to predict the spectra of prompt ν µ and ν τ neutrinos produced in the atmosphere by cosmic rays. Prompt leptons originate from the following sequence: the production of cc pairs, their fragmentation into charm hadrons which then undergo semi-leptonic decay. In the lower energy range (E < 10 7 GeV) it is possible to estimate the leptonic spectra by simply taking a product of moments of the various distributions, see, for example, Ref. [36,39]. However, for E > ∼ 10 7 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 where ρ(h) is the density of the atmosphere at vertical height h. We take the same exponential profile of the atmosphere 11 as was used in Ref. [7]. The sequence of equations determine φ 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 φ N (E, 0) is the known primary cosmic ray flux. All other initial fluxes are zero, that is φ a (E, 0) = 0 for with i = D ± , D 0 ,D 0 , D ± s , Λ c ; and finally where B(i → l) is the branching fraction of the decay of the charmed hadron i to lepton l. The nucleon attenuation length is where Z N N is the spectrum-weighted moment for nucleon regeneration and λ N is the interaction thickness n A (h) is the number density of air nuclei of atomic number A at height h and σ N A 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 /ρ in (15) arises from (19). For the nucleon-air cross section, σ N −air , we take the parameterization of Bugaev et al. [40]. For the incoming cosmic ray flux we take the parametrisation given in [36] denoted as TIG with knee. Also from [36] we take parametrisation of Z N N which depends on energy and takes into account the knee, which is consistent with [7]. This Z factor includes the regeneration by the p, n, N * . . . particles.
The attenuation length Λ i (E, X) of the charmed hadrons consists of two parts: the decay length λ dec i and attenuation due to their strong interactions with air nuclei. The decay lengths, λ dec i , of the various charmed hadrons depend, via their X dependence, on the density of the atmosphere, where m i and τ i , the mass and lifetime of the i th charmed hadron, are taken from [13]. The attenuation due to the strong interactions of the produced charmed hadrons with the air has a form similar to (18), namely λ i /(1 − Z cc ). To calculate λ 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 cc = (0.8) γ . 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 q ≃ 0.3 GeV is the mass of a light constituent quark. For simplicity, we take the same Z cc and λ i for all charm hadrons; for Λ c we expect the larger cross section to be approximately compensated by the larger Z cc . Thus, finally, Λ i is given by From (14) the light baryon flux is given by If we insert (15) into (16), then the individual charm hadron fluxes are where and x = E/E ′ .
The spectra of charmed hadrons dσ N →i /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. [6] to be More recent data [41] favour a smaller value 12 of the first ratio quoted in (25), namely We take this value in our analysis.
Note that the Λ c baryon is produced by the recombination of a c quark with a spectator diquark of the incoming nucleon 13 . It is not produced from ac quark. Therefore the parton momentum fraction x Λ carried by the Λ c is x Λ = x d + x c . The diquark momentum fraction x d will be less than 2 3 (1 − x c ), as part of the energy is carried away by the third valence quark (the factor 2 3 ), by thec quark and by gluons. To allow for this, we therefore take This is found to be in good agreement with the distribution generated by the PYTHIA Monte Carlo [42], which has a maximum in the region x Λ ≃ 0.5-0.6. Of course a very slow c quark is unlikely to combine with a fast light diquark (with x d ∼ 0.5). Therefore we introduce a cut-off, x c > x 0 c , in (27). Assuming a mean velocity of the c quark to be v 2 ∼ 0.25, we estimate x 0 c ∼ 0.1. For Λ b production, which we discuss below, the heavier b quark will carry a larger fraction of the Λ b momentum. In this case we take the cut-off to be x 0 b = 0.25. As a check, we also compute the prompt flux, arising from c → Λ 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 or, in the case of beauty, x Λ = 1.16x b , using m b = 4.5 GeV.
Rather than using a fragmentation function for dn c→D /dx, for D mesons we take x D = 0.75x c . This is sufficient for our purposes 14 . For illustration we compare, in Fig. 5, the prompt ν µ +ν µ flux at ground level, φ νµ+νµ (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 Λ c hadronization prescriptions (27) and (28) give very similar fluxes. At the highest energies the fraction of neutrinos coming from Λ c (relative to those from D) increases due to 12 For isolated c-quark fragmentation, PYTHIA [42] gives a ratio of 0.08. 13 The PYTHIA/JETSET Monte Carlo gives only 2.5% Λ c baryons in the fragmentation of an isolated c jet, that is in the absence of spectator diquarks.
14 That is the distributions dn c→i /dx were taken to be proportional to δ(x D − 0.75x c ) for D mesons, whereas for Λ c we assume that they are proportional to δ(x Λ − 1 2 (x c + 1)) for x c > 0.1, or δ(x Λ − 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 [43] in Ref. [13]. the short Λ c lifetime. Therefore the choice x Λ = 1 2 (1 + x c ) of (27), which corresponds to a larger x Λ than (28), gives a larger neutrino flux for E ν > ∼ 10 6 GeV. The results below correspond to using prescription (27).
For the leptonic decay of each charm hadron i, the distribution dn i→l /dx l was generated by the PYTHIA Monte Carlo [42]. Note that in PYTHIA 6.2 the D s → τ ν τ branching ratio was set to be 1%, whereas the latest value is 6.4 ± 1.5% [13]. Since D s → τ ν τ is almost the only source of prompt ν τ 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ν τ neutrinos than ν τ , since theν τ from τ + decay has the large x l . Of course, the reverse is true for D − s decay.

Prompt ν µ and ν τ fluxes
We present the predicted yields of prompt ν µ and ν τ neutrinos from charmed hadrons produced by cosmic ray collisions in the atmosphere. The flux of ν e neutrinos is essentially equal to that of ν µ 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. [44] that the prompt µ flux is only about 10% smaller than the prompt ν µ flux at the surface of the earth.
In Fig. 6 we plot the prompt ν µ and ν τ 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 γ = 2.02, in the actual calculation of the neutrino fluxes we used the observed primary cosmic ray flux, which corresponds to different values of γ above and below the 'knee'. The sharp fall-off for neutrino energies E > 10 8 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 < ∼ 10 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 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 ν µ (from π and K decays, etc.). Moreover, there is a small probability that atmospheric ν µ neutrinos may oscillate into ν τ neutrinos and so provide an 'atmospheric' ν τ flux. We calculate this flux using the 3σ ranges of the sin 2 θ ATM and ∆m 2 ATM neutrino mixing parameters found in an analysis [45] of the Super-Kamiokande [46] and MACRO [47] data. The resulting atmospheric ν τ 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 ν τ spectrum arising from bb 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 λ 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 φ ν in the plateau region (10 6 -10 7 GeV) of the E 3 φ ν 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 [48] gives the dotted curve in Fig. 8. The suppression of the original GBW continuous curve (for A = 1) is apparent.
7 Prompt ν τ flux from bb production At first sight it appears that if we also allow for bb production, with the same cross section, fragmentation and decay then we will approximately double the ν τ 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 2 c /m 2 b 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 ν τ neutrinos from beauty decays is not negligible. For charm, only the decay of the D s may produce ν τ neutrinos, whereas now B ± , B 0 , B s and Λ b semileptonic decays also give rise to a significant ν τ flux. Indeed including B and Λ b decays enlarges the predicted total prompt ν τ flux by about 40% for E ν ∼ 10 5 GeV, and even more at higher energies, see Fig. 10(b). This is considerably larger than the ν τ flux calculated in [49], where the beauty induced contribution grows from 1% at E ν = 10 2 GeV to about 10% at E ν = 10 6 GeV. The reasons why we obtain a larger fraction of ν τ neutrinos from beauty are, first, that we use (26), rather than (25), for c → D s hadronization and, second, our cross section for charm production is more suppressed at high energies by absorptive corrections, than the more compact bb production process. Fig. 10(a) shows the breakdown of the total prompt ν τ flux. We take the ratios of the B ± , B 0 , B s and Λ b beauty yields (after hadronization of the b quark) to be given by exactly analogous relations to (25) and (26) for charm. Recall that ν τ 's of cc and bb origin come, respectively, from D s → τ ν τ and from the semileptonic τ ν τ decay modes of B ± , B 0 , B s and Λ b . Thus we have a direct ν τ component and an indirect component coming from the τ → ν τ decays. For charm the former component is very small since the direct ν τ carries away a small energy fraction. Thus the τ → ν τ decay gives the dominant component, which is truncated at an energy when the τ has insufficient time to decay. For beauty the direct and indirect components are comparable until the energy regime is entered at which the τ 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 GeV. We have approximated the b → B 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 Λ 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 bb contribution to the ν µ (and ν e ) fluxes. They never exceed about 2-3% of that arising from cc production and decay.

Implications for Physics
The ν µ flux shown in Fig. 7(a), and the ν τ 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 ν τ 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 ν µ and ν τ are equal for incoming cosmic neutrinos. The flux of prompt ν τ neutrinos (which arises mainly from D s → τ ν τ decays, but with a significant 40% contribution from B, B s , Λ b semileptonic decays) is about ten times less than that for prompt ν µ , and the number of ν τ neutrinos produced from the conventional atmospheric flux (via ν µ → ν τ oscillations) is negligible for E > 10 4 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 ν τ flux. For values of E ν above 10 4 GeV the ν µ flux is significantly depleted in the passage of the neutrinos through the Earth. The absorption increases rapidly as E ν is increased [50]. On the other hand, high energy ν τ 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 τ lepton which subsequently decays, regenerating a ν τ neutrino with degraded energy [14]. It is, moreover, interesting to note that tau neutrinos of energies about 10 7 GeV may produce spectacular distinctive signatures in the planned 1 km 3 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 ν τ interaction may produce a contained 'double bang' signature [51] or a 'lollipop' signature [52]. In the former, the first bang corresponds to the hadronic shower produced along with the subsequent τ lepton and the second bang is the shower associated with the τ decay. A lollipop event is when only the second shower occurs within the detector and the τ 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 ν e neutrinos comes mostly from kaon decays and so the relative flux ν e /ν µ < ∼ 0.1 [53]. On the other hand, the prompt ν e and ν µ fluxes are practically equal to each other. Hence the prompt ν e flux should be more visible over its conventional background, than the prompt ν τ flux 15 . Of course, the background reduction is much less than for ν τ neutrinos, but this will be partially offset since it is likely that ν e neutrinos will be easier to identify than ν τ 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 ν > 10 8 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.

Uncertainties in the predicted prompt ν µ and ν τ fluxes
We have emphasized the importance to neutrino astronomy of reliable predictions of the prompt muon and tau neutrino fluxes, which arise from cc and bb 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 bb production 16 , and should also hold for cc production since the factorization scale µ F ∼ m c .
It is important to note that at fixed target energies, E ∼ 250 GeV, the NLO predictions for dσ/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 x ∼ 10 −4 , predictions, which use partonic structure determined from these data and which agree with fixed target charm data, should be reliable up to about E ∼ 10 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 cc 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 Λ c hyperons. The recombination of the c quark with the ud diquark of the incoming proton gives the Λ c a rather large x. We use the simplified formula x Λ = 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 ν flux is proportional to the x γ=2.02 moment, (11), of the x distribution. Thus the larger value of x carried by the Λ c (as compared to the D mesons) compensates for the small probability of Λ c formation, see (25). As a consequence, the contribution of the Λ c to the prompt ν µ flux increases from about 10% to 40% as the neutrino energy increase from 10 4 to 10 9 GeV. Unfortunately the x distribution of Λ 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 Λ = 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 ν τ flux from charm, which relies entirely on the D s → τ ν τ 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 ν τ flux from beauty, since it originates more uniformly from B ± , B 0 , B s and Λ b semileptonic τ ν τ 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 ν µ and ν τ 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][3][4][5]. To ±5% accuracy for 0.05 < x < 0.6 we have where β = 0.05 − 0.016 ln(E/10 TeV) and    (11), of charm production in proton-nucleon collisions, as a function of the energy of the produced charm quark. The models are as in Fig. 2. The reason why the KMS prediction falls below the other predictions at the smaller values of E c is explained in Section 3. The γ = 2.02 moments are shown for illustration; in the calculation of the neutrino fluxes, the differential cross section x F dσ/dx F is convoluted with the observed primary cosmic ray flux. curves respectively include the growth of the proton radius, triple-Pomeron effects and the combination of these two effects. The dot-dashed curve is the scaling prediction of (10), but with σ inel corresponding to p-air collisions. The γ = 2.02 moments are shown for illustration; in the calculation of the neutrino fluxes, the differential cross section x F dσ/dx F is convoluted with the observed primary cosmic ray flux.  Figure 6: The prompt ν µ +ν µ and ν τ +ν τ fluxes, arising from cc production, fragmentation and decay, obtained using the three different extrapolations of the gluon to very small x (described in the caption to Fig. 2). For the MRST and KMS models, the A dependence is taken to be dσ/dx(p + air → c + . . .) = Adσ/dx(pN → c + . . .).  Figure 7: The prompt (a) ν µ +ν µ and (b) ν τ +ν τ fluxes calculated using the charm production cross sections corresponding to the shaded band in Fig. 4. Also shown are the conventional muon and tau neutrino atmospheric fluxes, where the latter originates, via neutrino mixing transitions, from the former. There is also a contribution to the prompt ν τ +ν τ flux from beauty production, which is not included here, but is shown in Fig. 10(b). The prompt ν e +ν e flux is equal to that for ν µ +ν µ , but the atmospheric flux is a factor of 10 or so less, see the discussion in Section 8.1. with charm-air rescatt.  Figure 8: The dashed and dotted curves correspond, respectively, to the prompt ν µ +ν µ flux obtained by switching off the charmed hadron-air interactions and assuming that the incoming cosmic rays have A = 7 rather than A = 1. The default continuous curve (with air interactions and A = 1) is based on the GBW gluon extrapolation. In this work we assume that the c → Λ c hadronization is given by x Λ = 1 2 (1 + x c ) for x c > 0.1.   Figure 10: The prompt ν τ +ν τ fluxes originating from cc and bb production and decay, which respectively arise from the D s → τ ν τ decay, and from the B ± , B 0 , B s and Λ b semileptonic τ ν τ decay modes. The upper plot shows the breakdown into the direct ν τ contribution (continuous curves) and the indirect τ → ν τ contribution (dashed curves). The lower plot shows the total prompt ν τ +ν τ flux, together with its components of cc and bb origin. Also shown is the nonprompt ν τ +ν τ flux arising from ν µ → ν τ oscillations from the conventional atmospheric ν µ flux.  Figure 11: The continuous curve is the prompt ν µ +ν µ flux from horizontally-incident cosmic rays. The dashed curve, which corresponds to vertically-incident cosmic rays, is taken from Fig. 7.