Prompt atmospheric neutrinos in the quark-gluon string model

We calculate the atmospheric flux of prompt neutrinos, produced in decays of the charmed particles at energies beyond 1 TeV. Cross sections of the $D$-mesons and ${\Lambda}^{+}_{c}$ baryons production in pA and $\pi$A collisions are calculated in the phenomenological quark-gluon string model (QGSM) which is updated using of the recent measurements of cross sections of the charmed meson production in the LHC experiments. A new estimate of the prompt atmospheric neutrino flux is obtained and compared with the limit of the IceCube experiment as well as with predictions of other charm production models.


Introduction
At present time the operating neutrino telescopes focus on the detection of astrophysical high-energy neutrino fluxes: IceCube, a cubic kilometer detector at the South Pole [1,2,3], ANTARES [4,5] located in the Mediterranean Sea, and underwater Baikal Gigaton Volume Detector (Baikal-GVD), a cubic kilometer-scale array, which is currently under construction in lake Baikal [6,7].
The Baikal-GVD has a module structure and consists of functionally independent sub-arrays (clusters) of optical modules (OMs) and is designed to detect astrophysical neutrino fluxes at energies from a few TeV up to 100 PeV. Five clusters have been already installed (every includes 288 OMs), current instrumented volume of the Baikal-GVD is the largest in the Northern Hemisphere (∼ 0.25 km 3 ), and the first highenergy neutrino induced events are reconstructed. The first phase (GVD-1) to be completed by 2021 and will comprise a e-mail: sinegovsky@jinr.ru b e-mail: sorokovikov@jinr.ru 9 clusters (2592 OMs), the full-scale GVD with an instrumented volume about of 2 km 3 will consist of 10 4 light sensors.
The diffuse flux of high-energy astrophysical neutrinos was revealed in 2013 at IceCube detector [8,9], and for 6 years about 100 neutrino events were detected in IceCube experiment [2]. Another important discovery was made recently: on 22 September 2017 IceCube have detected the high-energy neutrino event coincident both in the direction and the time with the gamma-ray flare from the blazar TXS 0506+056 [10]. The event was later confirmed [11] by the archival IceCube data which display an excess of high-energy neutrino events (against the atmospheric neutrino flux) between Sept. 2014 and Mar. 2015, that give 3.5σ -evidence for neutrino flux from the direction of TXS 0506+056. This supports the hypothesis that the blazar TXS 0506+056 is the individual source of high-energy neutrinos and, presumably, the source of high-energy cosmic rays.
Essential progress has been achieved in experimental studies of astrophysical and atmospheric neutrino fluxes, however the prompt atmospheric neutrinos have not yet been detected. High-energy neutrinos arising from decays of mesons and baryons, produced in hadronic collisions of cosmic rays with Earth's atmosphere, compose the background against the neutrinos from distant astrophysical sources. The atmospheric neutrinos comprise two components, which are distinguished by zenith-angle distributions and the energy spectra. The anisotropic component, arising from decays of pions and kaons, has the softer spectrum ("conventional" or π, K-neutrinos). The second component, quasi-isotropic flux produced at higher energies, mainly in decays of short-lived heavy charmed mesons and baryons D, Λ + c , is characterized by harder spectrum. This component ("prompt" neutrinos) is most uncertain because of scarce measurements and wide spread in model predictions of the charm production cross sections at very high-energies.
The high energy interactions of cosmic rays with the Earth's atmosphere are dominated by the soft processes with small momentum transfer, which are beyond the scope of perturbative technique of the quantum chromodynamics (QCD). Perturbative QCD models of the charmed particle production encounter difficulties related to the nonperturbative dynamics contribution. Thus, the elaboration of phenomenological models beyond the pQCD is required for comprehensive research of the charm production in hadronic interactions at high energies.
The quark-gluon string model (QGSM) was developed [12,13,14] to describe the soft and semihard hadronic processes at high energies: it has been applied for successful explanation of characteristics of mesons and baryons production in hadron-nucleon collisions. QGSM was one of the first models to estimate the prompt atmospheric neutrino flux [15]. The recent data on the cross sections of charmed particle production, obtained in experiments at the LHC [16,17,18,19], allow the QGSM free parameters updating. The updated version of QGSM is applied to calculate the prompt neutrino flux in the neutrino energy range 1 TeV -100 PeV. The calculation is based on the hadronic cascade model [15,20,21] and cross sections of D meson and Λ + c baryon production in pA-and πA-collisions which are computed with the updated QGSM. We compare our result with the constraint obtained in the IceCube experiment [1] as well as with predictions of the color dipole model (ERS) [22], SIBYLL 2.3c [23], the NLO pQCD models, BEJKRSS [24] and GRRST [25].

Production of charmed particles in QGSM
The nonperturbative quark-gluon string model (QGSM) gives unified descriptions of the soft hadronic processes. The model is based on the reggeon calculus, the topological 1/N c -expansion of the amplitudes and the color string dynamics (see for more details [12,13,26,27,28] and references therein). The QGSM, having a small number of parameters, succeeded in describing of the multiparticle production in hadron-nucleus collisions at high energies [12,13,14,26,27,28,29,30,31].
The planar and cylinder-type diagrams involve quarkquark, quark-gluon and gluon-gluon scattering. The summation of qq, qg and gg-diagrams leads to the Regge-behavior of scattering amplitude. The s-channel cuttings of cylindric diagrams describe the multiparticle production, and in the t-channel, these diagrams correspond to gluon exchanges. It is argued [31] that cylindric diagrams lead to the pomeron pole as the color singlet made up of sea quarks and soft gluons. The pomeron can be related to a sum of the ladder diagrams with exchange of reggeised gluons. Sum of gluon ladders with possible quark loop insertions may produce the pomeron trajectory. The simplest two-gluons exchange leads to pomeron with the intercept α P (0) = 2S g − 1 = 1. Account of mixing between gg and qq Regge trajectories (glueballs and qq resonances) at the small t, and effects of the small distance perturbative dynamics lead to the supercritical pomeron with ∆ = α P (0) − 1 ∼ 0.15 − 0.25, which ensures the Froissart behaviour of the total cross section.
To calculate inclusive cross sections of charmed hadron production, one needs know the distribution functions of the quarks of the colliding particles and the fragmentation functions of the quarks and diquarks. The inclusive cross sections of charmed hadron production are defined as the convolution of distribution functions of the valence (and the sea) quarks and diquarks of the colliding particles with the functions of quarks (diquarks) fragmentation into a charmed hadron. These functions are expressed in terms of intercept α R (0) of the Regge trajectory (α R (t) ≃ α R (0) + α ′ R t in linear approximation), including the α ψ (t) trajectory of the cc bound states. The complete set of distributions and fragmentation functions can be found in Refs. [13,26,27,28,29].
For a nucleon target, the inclusive cross section of production of a hadron h (h = D + , D − , D 0 ,D 0 , Λ + c ) can be written as a sum over n-pomeron cylinder diagrams: where σ n (s) is the cross section of the 2n-strings (chains) production, corresponding to the s-channel discontinuity of the multipomeron diagrams (n cut pomerons and arbitrary number of external pomerons taking part in the elastic rescattering); ϕ h n (s, x) is the x-distribution of the hadron h produced in the fission of 2n quark-gluon strings: ϕ h 0 (s, x) accounts of the contribution of the diffraction dissociation of colliding hadrons, n = 1 corresponds to the strings formed by valence quarks and diquarks, terms with n > 1 are related to sea quarks and antiquarks; x = 2p / √ s is the Feynman variable, p is the longitudinal momentum of the produced hadron, √ s is the total energy of the two colliding hadrons in the c.m.f., is the mean square of the transverse momentum, and m h is the mass of the hadron h.
The cross sections σ n (s) were calculated [32] in the quasieikonal approximation which accounts of the low-mass diffractive excitations of the colliding particles and corresponds to maximum inelastic diffraction consistent with the unitarity condition. Only nonenhanced graphs were consid-ered with neglect of interactions between pomerons: where Here σ P (s) is the pomeron contribution to the total cross section, z(s) is the function representing a relative contribution of the successive rescatterings, γ P and R 2 are characteristics of the pomeron residue; s 0 = 1 GeV 2 . The term with n = 0 in (1) corresponds to the elastic scattering and the diffrac- The parameter C = 1 + σ DD /σ el > 1 takes into account the the low-mass diffractive dissociation. The values of above parameters are found from experimental data on the total and differential cross sections of elastic pp and pp scattering at high energies [14,27,31,33]: In the case of D meson production in pp interaction, the functions ϕ h n (s, x) can be written [13] as follows: where x ± (s) = 1 2 x 2 + 4m 2 ⊥ /s ± x . For π − p interactions [13,27]: The functions F qq (x ± ), and F D(n) q sea (x ± ) defined as convolution of the quark distributions with the fragmentation functions, take into account contributions of the valence quarks, diquarks, and sea quarks. For example, in pp collisions [26,27,28]: In case of π − p interactions [27,28]: where f j p (x), f j π (x) are the distribution functions of quarks, antiquarks and diquarks in colliding hadrons, j = q,q, qq; G D j (x/x 1 ) are the fragmentation functions. At limits x → 0 and x → 1 these functions are defined by Regge asymptotics, and for the intermediate values of x the interpolation is used [13,26,27]. In particular, More details on the functions ϕ h n (s, x), f j p (x) and G D j (x/x 1 ) can be found in [13,26,27,28,29].
The distribution function ϕ h n (s, x) for the case of Λ c production in pp collisions can be written [26,35] as following: Here F 1qq denotes the distribution at the leading diquark fragmentation with weight a Λ c 1 and F 0qq is the distribution for the nonleading fragmentation of diquarks written with the central density parameter aΛ c 0 : where  [16,17,18,19,37,38,39,40,41,42,43,44,45].
Charmed baryons Λ c have harder spectrum in comparison with charmed mesons in the region x > 0.1. Fragmentation process of the charmed baryon differs from fragmentation of D mesons since Λ c baryon consist of three quarks. The diquark fragmentation functions are divided into two parts which describe different kinematical regions: F n 0qq (x) (central region) and F n 1qq (x) (fragmentation region). The first term in (13) corresponds to direct Λ + c baryon production and the second one is connected to a pair pro- In pp interactions Λ + c baryon can be produced by leading ud diquark, that leads to enhancement of Λ + c spectra in comparison with Λ − c . Accounting the sea diquark contribution is important only for antibaryon production in the forward region, while for baryon spectra sea diquarks contribute small practically at all x [36]. At x → 1 dominate direct Λ + c production in comparison with the pair Λ + c /Λ − c production because of diquark function fragmentations in case of pair production is suppressed in forward region by additional term The distribution functions of charmed particles in (5) contain free parameters that cannot be calculated within the framework of the quark-gluon string model, and their values should be found from a comparison with experiments. The intercept α ψ (0) of the poorly studied cc -trajectory noticeably affects the cross sections of D meson production. Two values of α ψ (0) have been used by QGSM authors [13,31]: α ψ (0) = −2.2, obtained from the mass spectrum on the assumption of linear Regge trajectory, and α ψ (0) = 0 (nonlinear trajectory), derived from the perturbative calculations. Basing on experimental data on Λ c production, they consider [31] nonperturbative value α ψ (0) = −2.2 as preferable one. If the Regge trajectory α ψ (t) is linear (similar to light hadrons), then the x-distributions of charmed particles become softer in comparison with the case of α ψ (0) = 0.
The coefficient a 1 provides an unified description of the kinematic regions x → 0 and x → 1 in the case of leading fragmentation (when the valence quarks take part in the fragmentation). Now there are no clear arguments for choice of its value, and different authors apply various values among which two extreme values may be chosen: a 1 = 2 [30] and a 1 = 30 [13]. New measurements of the total cross sections of charmed meson production at high energies in the experiments ALICE [16,17,18] and ATLAS [19] allow a check of the QGSM predictions for extreme values of the parameter a 1 .
The parameter a h in (5)   collisions at E lab = 500 GeV. Experimental data are from [46]. Same notation for lines as in Fig. 3. The results of calculation of the cross sections of D meson production in pp collisions in comparison with experimental data are shown in Figs. 1-3. The total cross section of D/D mesons production in pp collisions as a function of center-of-mass energy is calculated in the QGSM for four sets of free parameters (Fig. 1). Here the experimental results in a wide energy range [37,38,39,40,41,42,43,44,45] including LHC measurements [16,17,18,19] are presented. Calculation with α ψ (0) = 0 and a 1 = 30 does not agree with experimental data at √ s < 1 TeV, while the calculations with α ψ (0) = −2.2 are in close agreement with the measurements in a wide energy range. At low energies the cross sections calculated with α ψ (0) = −2.2 for extreme values a 1 = 2 and a 1 = 30 differ by a factor 2 − 5, but the influence of the parameter a 1 tends to diminish with energy and becomes negligible at high energies ( √ s > 1 TeV). Figure 2 represents the calculated differential cross sections of D mesons production at the laboratory energies 400 GeV and 800 GeV in comparison with the measurements of experiments NA27 [40] and E743 [41].
As can be seen from Figs. 2, 3, QGSM with the intercept α ψ (0) = −2.2 better describes the experimental data for differential cross sections. Figure 3 presents the comparison of the experimental data (pp collisions, 400 GeV) [40] with calculations of differential cross sections for each sort of D mesons production (D + , D − , D 0 ,D 0 ). The cross sections of D + and D 0 weakly depend on the parameter a 1 , while cross sections of D − andD 0 mesons calculated for a 1 = 2 are smaller (e.g. for α ψ (0) = −2.2 by a factor 2 − 10) as compared to the case of a 1 = 30. The intercept α ψ (0) = −2.2 noticeably better describes the experimental data for D − and D 0 , whereas values α ψ (0) = 0, a 1 = 2 lead to better description of D 0 mesons. In spite of data spread for D + , the intercept α ψ (0) = −2.2 seems preferable.
The factor a 1 amplifies the contribution of the leading fragmentation (D − /D 0 ) with participation of the valence quarks. Production of D − andD 0 in pp interactions has a higher probability because these mesons contain the valence quarks of colliding protons. The contribution of the leading fragmentation functions dominates, and x-distribution of D − and D 0 is harder in comparison with D + and D 0 . The influence of the parameter a 1 on the cross section of all D mesons production is also noticeable (Fig. 2).
The calculations of cross sections of charmed meson production in π − p collisions are compared to experimental data in Figs. 4-7. Figure 4 shows the differential cross sections of neutral D 0 /D 0 mesons production in π − p collisions at energy 500 GeV in comparison to the data of the experiment E791 [46]. The calculation for values α ψ (0) = −2.2 and a 1 = 2 agrees with the measurement data in the small x range, and prediction obtained for α ψ (0) = −2.2 and a 1 = 30 better describes the experimental data in the fragmentation region. The differential cross sections of D/D mesons computed at energy 360 GeV in comparison to the data of experiments WA92 (350 GeV) [47] and NA27 (360 GeV) [48] are presented in Fig. 5. In spite of small differences of the beam energies, data of the two experiments differ at small-x and at x > 0.5. The calculation for α ψ (0) = −2.2 and a 1 = 2 describes well most of experimental points.
The differential cross sections for each type of D mesons in comparison to measurements (WA92) are shown in Fig. 6.   [52]. The same notation for lines as in Fig. 9. The calculations of the leading and nonleading differential cross sections of D mesons production in π − p collisions at energy 360 GeV (LF) are shown in Fig. 7 along with experimental data of NA27 [48].
Comparison of the cross section of Λ + c baryon production in pp collisions with experimental data is shown in Fig. 8. The differential cross section was calculated at √ s = 62 GeV (left panel) the experimental data were obtained for energies √ s = 62 GeV [49] and 63 GeV [50]. There is appreciable difference of the cross section measurements in these two experiments. The calculation with the parameter α ψ (0) = −2.2 agrees with the later experiment.
The right panel of Fig. 8 shows the total cross section of Λ + c production as a function of center-of-mass energy. The experimental points are taken from Ref. [49], the calculation was made for the same parameter sets. The large spread of the total cross section data prevents from making definite choice of the intercept α ψ (0).
In Figs. 9-11, we show the differential cross sections of D mesons production in pp collisions as a function of the rapidity (y) in comparison to LHCb measurements. Experimental data were obtained at energies √ s = 5 TeV [51], 7 TeV [52] and 13 TeV [53] for rapidity range 2 y 4.5 that corresponds to x 10 −3 − 10 −2 . The points plotted in Figs. 9-11 were obtained from original experimental data by summing them over transverse momentum bins (for each bin in y).
There is expected the problem in describing the experimental data on D mesons in the small x range. In the QGSM version under consideration, the inclusive spectra of charmed particles production are averaged over transverse momentum, while the LHCb data were obtained for the transverse momentum interval p ⊥ 15 GeV/c at the rapidity values 2.0 − 4.5. However the calculation with α ψ (0) = −2.2 describes satisfactorily the experimental measurements on the production D ± mesons (unlike D 0 /D 0 ) at energies √ s = 5 TeV and 7 TeV (Figs. 9, 10). The experimental data on neutral D 0 /D 0 and charged D + /D − mesons differ by factor 2 − 3 (at fixed energy), while QGSM predicts close values of dσ /dy. That is the model describes the cross sections of charged mesons much better than D 0 /D 0 . One possible explanation of this discrepancy is that D 0 /D 0 events contain a mixture of vector mesons D * 0 /D * 0 , decays of which might contribute to pseudoscalar D 0 mesons.
It is possible also that in case of the small-x events and large transverse momentum at high energy we encounter with the problem of "enhanced" diagrams, relating to interactions between pomerons, which are neglected in the quasieikonal approximation. Account of enhanced diagrams leads to x-distributions, rising as 1/x at small x [31]. Perhaps it will take a significant revision of the QGSM parameters and calculation technique for this region of kinematics.
Note however that small-x region (10 −4 − 10 −3 ) gives minor contribution to the atmospheric neutrino flux because of the dominating peripheric processes in the cosmic-ray induced hadronic cascade: the small-x values are suppressed under the integral by a factor x γ , where γ is the spectral index of cosmic ray protons (γ ≈ 1.7 − 2.0).

QGSM in comparison with different charm production models
Before comparing of the prompt atmospheric neutrino fluxes predictions it would be useful to confront cross sections of charmed particles production of different models. The comparison of the differential cross sections of charmed mesons production in pp collisions for proton energies in the laboratory frame (10 3 and 10 6 TeV) is shown in the top panel of Fig. 12: QGSM (solid and dashed lines), SIBYLL 2.3c [23] (dotted line), perturbative QCD model (BEJKRSS) [24] (the line with symbols), and the dipole model (DM) [22] (dashdotted). The pQCD calculation is rather close to the present work results obtained with parameters α ψ (0) = −2.2, a 1 = 2 (10 3 TeV) and a 1 = 30 (10 6 TeV). Our calculation lies below the DM result for most of the x range, that should lead to the lowered prompt neutrino flux as compared with the result of Enberg et al. [22].
The total cross sections of D mesons production in pp collisions as a function of center-of-mass energy are shown in the bottom panel of Fig. 12 for QGSM, SIBYLL 2.3c [23] and pQCD model (BEJKRSS) [24]. Predictions of the QGSM (for a 1 = 2) and pQCD model are almost the same in a wide energy range, with the exception of energies √ s < 100 GeV (calculation for a 1 = 30 gives large cross sections at √ s < 10 3 GeV).
To calculate the prompt neutrino flux one needs know the cross sections of charmed particles production in col-lisions of hadrons with atmospheric nuclei. The cross sections are recalculated for a nuclear target with average mass number A according to the formula dσ pA /dx = A α dσ pp /dx (for the air we take A = 14.5). The index α depends on x: α ≈ 1 at x → 0 and monotonically decreases with rise of x (α ≈ 0.5 at x → 1) [13]. In [15] the prompt neutrino flux has been calculated for constant α ≈ 0.72 (authors argued that uncertainty due this choice does not exceed 15%), and we use the same value.

Energy spectra of the prompt atmospheric neutrinos
In the present work, the calculation of prompt neutrino fluxes is performed with the method [15,20] for QGSM for parameterization of cosmic ray spectrum by Nikolsky, Stamenov, Ushev (NSU) [54], the toy model by Thunman, Ingelman, Gondolo (TIG) [55], and the recent model for cosmic ray spectrum by Hillas and Gaisser (H3a) [56]. The NSU spectrum which takes into account an elemental composition of primary cosmic rays was chosen in order to compare new result with the old one [15,21]. We use also the toy spec- trum TIG (else called the broken power law, BPL) only in order to compare results of different calculations, including the dipole model prediction [22] which was used by IceCube as a benchmark model. The three-component model with mixed extragalactic population, H3a was chosen to compare the QGSM calculation with SIBYLL 2.3c [23], the NLO pQCD predictions BEJKRSS [24] and GRRST [25], as well as with the IceCube experiment limitation on the prompt neutrino flux [1].
Nucleon-nuclear interactions give the main contribution to the prompt neutrino flux, while reactions π ± A add less than 5-10% to the prompt neutrino flux in the energy range 10 − 10 5 TeV. The Fig. 13 shows relative contributions of NA and πA interactions to the prompt fluxes calculated for NSU spectrum and two values of the free parameter of the quark-gluon string model: a 1 = 2 (solid line) and a 1 = 30 (dashed one). Contributions of D mesons and Λ + c baryons are presented in Fig. 14. Figure 15 shows the calculation of the prompt atmospheric neutrinos flux (scaled by E 2 ν ): the band represents this work calculation for the NSU spectrum and the QGSM with the intercept α ψ (0) = −2.2. The band shows uncertainty due to change of the parameter a 1 which ensures unified behavior of the leading fragmentation functions at x → 0 and x → 1. Extreme values of a 1 lead to change of the neutrino flux by a factor 1.4: a 1 = 2 corresponds to lower bound, and a 1 = 30 to upper one. However, influence of the intercept of Regge trajectory α ψ (0) appears to be more substantial: the replacement of α ψ (0) = 0 by α ψ (0) = −2.2 reduces the flux by a factor 3 as compared to the result [15] (solid line) obtained for the same scheme (QGSM+NSU) with intercept α ψ (0) = 0. The dipole model result [22] for the TIG cosmic ray spectrum is also shown in Fig. 15 (dashed line).
The QGSM flux [15] performed for α ψ (0) = 0 and NSU spectrum was considered by IceCube collaborators as too optimistic prediction [57,58]. At the energies E ν > 10 6 GeV it exceeds the ERS result [22] by about 30%, however part E E Fig. 15 Vertical flux of prompt atmospheric neutrinos calculated with QGSM for the NSU cosmic ray spectrum. Shaded band: this work with α ψ (0) = −2.2. The band width corresponds to the variation of a 1 from 2 (lower bound) to 30 (upper bound). The solid line presents the result from Refs. [15,21] obtained with α ψ (0) = 0. Dash line plots the dipole model calculation [22] for the TIG spectrum.
of this excess is related to the difference of the cosmic ray spectra used.
An influence of charm production models on neutrino fluxes is seen in Fig. 16. All results are obtained for the same cosmic ray spectrum TIG. In the energy range beyond 1 PeV, where atmospheric neutrinos from charmed particles dominate, QGSM (shaded band) leads to appreciably lower flux as compared to the dipole model result [22] (dashed line). The predictions of the pQCD models, BEJKRSS [24] and GRRST [25] are compatible in the whole energy range. The QGSM flux at energies above 200 TeV is close to upper bound of the BEJKRSS band.
The uncertainty of the QGSM neutrino flux prediction due to variation of the intercept α ψ (0) is shown in Fig. 17. Uncertainty band is computed for α ψ (0) = −2.2 with a D = 7.0 ·10 −4 (lower bound) and α ψ (0) = 0 with a D = 10 −3 (upper bound). The band calculation is performed for H3a cosmic ray spectrum with use the value a 1 = 2. The bulk of prompt neutrino flux uncertainties is induced by intercept of the charm Regge trajectory (as Fig.17 shows). Rest parameters, a D , a 1 slightly disturb bounds of the uncertainty span. In fact, the uncertainty due to the α ψ (0)-variation absorbs that for the a 1 -variation. α ψ (0) = −2.2 seems to be preferable from the viewpoint of the most experimental measurements on charm production. However, α ψ (0) = −2.2 leads to underestimation of the neutral D meson production at high energy. On the other hand α ψ (0) = 0 predicts harder charm spectra at all x. The largest differences in dσ /dx appear at x → 1, where the charm spectra calculated for α ψ (0) = 0 are larger by several orders of magnitude in comparison with case of α ψ (0) = −2.2. The x ∼ 1 is the region of peripheric processes which E E Fig. 16 Comparison of the prompt neutrino flux calculated with use of the models: QGSM with α ψ (0) = −2.2 (shaded band), the color dipole model (ERS) [22] (dash line), NLO pQCD models, BEJKRSS [24] (hatched area) and GRRST [25] (dash-dotted line). All calculations are performed for the TIG cosmic ray spectrum. Solid line denotes the IceCube upper limit and the dotted one its extrapolation [1] (see text for more details). Fig. 17 Uncertainties of the prompt neutrino flux calculations. The blue band corresponding to QGSM calculations shows the effect of varying crucial parameter: the intercept α ψ (0) = 0 (upper bound of the band) and α ψ (0) = −2.2 (lower bound). The hatched area presents the scale of uncertainties of the flux calculation with the NLO pQCD model BEJKRSS [24]. The SIBYLL 2.3c prediction [23] is compatible with QGSM (the line inside band). All calculations are performed for the H3a cosmic ray spectrum. The solid line with arrows represents the IceCube limitation and dotted line labels its extrapolation [1]. dominate in the atmospheric hadronic cascade and atmospheric neutrino fluxes turn out to be sensitive to the intercept value: the replacement of α ψ (0) = −2.2 by α ψ (0) = 0 increases flux by a factor 2 − 3.

E E
Relative to the BEJKRSS result, the QGSM gives larger flux at low and middle energies. However, our calculation comes near to the BEJKRSS one with energy growth, and the QGSM lower bound is compatible with upper bound of the pQCD prediction at E ν > 10 PeV. On the other hand, our Fig. 18 The comparison of the prompt neutrino flux calculated for the H3a cosmic ray spectrum [56] with use of different models: QGSM (α ψ (0) = −2.2), SIBYLL 2.3c [23], the pQCD models BEJKRSS [24] and GRRST [25]. Solid line is the IceCube constraint, dotted line plots the limit extrapolation [1]. uncertainty band completely covers the SIBYLL 2.3c flux in the whole energy range. Figure 18 presents our preferred result: the prompt neutrino spectrum calculated with QGSM for the magnitudes of parameters α ψ (0) = −2.2, a D = 7 · 10 −4 . The brown band marks the QGSM uncertainty related to the parameter a 1 which rises from 2 (lower bound) to 30 (upper bound). Also shown here are the SIBYLL 2.3c [23] calculation and the NLO pQCD ones (BEJKRSS [24] and GRRST [25] models). All these spectra were obtained for H3a cosmic ray spectrum. The QGSM flux for a 1 = 30 and α ψ (0) = −2.2 is close to the SIBYLL 2.3c calculation beyond 100 TeV and greater than the BEJKRSS flux by a factor ∼ 2 − 4 at 10 5 − 10 6 GeV. With the energy rise these differences decrease, and our prediction at E ν 10 PeV is close to the BEJKRSS flux.
In the IceCube analysis [1], upper limit on the prompt atmospheric neutrino flux was obtained using high statistics collected over six years. The prompt neutrino flux was constrained using the color dipole model [22] corrected for the cosmic ray spectrum parametrization H3p [56]. The solid line with arrows in Figs. 16,17,18 shows the constraint specified by IceCube for the sensitive energy range 9 − 69 TeV. The extrapolation beyond 69 TeV was made by IceCube using again the prompt neutrino spectrum [22]: the original result was also brought into accord with the cosmic-ray model H3p [56]. The dotted line in Figs. 17,18 just shows the dipole model behavior beyond of the sensitivity range with the best-fit parameters describing the systematic uncertainties of the IceCube analysis, e.g., the optical efficiency of the telescope, the Antarctic ice properties, the uncertainties of the cosmic ray spectrum.
Our result is evidently compatible with the IceCube limitation, the same is true for the rest models under discussion: SIBYLL 2.3c [23], the NLO pQCD models BEJKRSS [24] and GRRST [25].

Conclusion
The new calculation of the atmospheric neutrino flux from decays of the charmed particles is performed with updated version of the quark-gluon string model. The QGSM parameters α ψ (0) and a 1 are examined by comparison of calculated cross sections for the charmed meson production with the data of measurements obtained at LHC and in other experiments. Though the data of the LHCb experiment do not allow unique choice of the varying parameters magnitudes, we consider the intercept α ψ (0) = −2.2 as a preferred value against α ψ (0) = 0. Note, that α ψ (0) = −2.2 appears proper magnitude, because it is in accordance also with observable (and natural) pattern: heavier quarks have lower intercept of the Regge trajectory.
At high energy the differential cross sections of charm production are more sensitive (as compared with the total cross section) to the change of the parameter a 1 , which brings the neutrino flux deviation about (20 − 40)% for extreme values of a 1 . The analysis shows that intercept of Regge trajectory α ψ (0) causes more noticeable effect on the charm production and therefore on the prompt atmospheric neutrino flux. Updated version of QGSM with α ψ (0) = −2.2 leads to decrease of the prompt neutrino flux by a factor ∼ 2 − 3 as compared to the former QGSM prediction [15,21] obtained with α ψ (0) = 0.
In the energy range beyond 1 PeV, where atmospheric neutrinos from the decay of charmed particles dominate, as it is expected, the new QGSM flux is significantly lower in comparison with the color dipole model [22]. The QGSM flux obtained for intercept α ψ (0) = −2.2 and H3a cosmic ray spectrum is compatible with the NLO pQCD predictions at E ν > 10 PeV, and upper bound of our calculation does not differ practically from the SIBYLL 2.3c result. The updated QGSM calculation of the prompt atmospheric neutrino flux is consistent with the IceCube limitation.
The performed calculations confirm viability of QGSM as appropriate phenomenological model of the high-energy hadronic interactions which allows account of effects beyond the pQCD. The updated version of QGSM is the suitable approach to provide reasonable prediction of the atmospheric prompt neutrino flux. Undoubtedly, the current version QGSM must be developed taking into account the large-mass diffraction dissociation and resonances production. Also all parameters of the model need in comprehensive revision and the thorough fitting by use of recent and future results of high-energy hadronic experiments.