Charm production in SIBYLL

. SIBYLL 2.1 is an event generator for hadron interactions at the highest energies. It is commonly used to analyze and interpret extensive air shower measurements. In light of the ﬁrst detection of PeV neutrinos by the IceCube collaboration the inclusive ﬂuxes of muons and neutrinos in the atmosphere have become very important. Predicting these ﬂuxes requires understanding of the hadronic production of charmed particles since these contribute signiﬁcantly to the ﬂuxes at high energy through their prompt decay. We will present an updated version of SIBYLL that has been tuned to describe LHC data and extended to include the production of charmed hadrons.


Introduction
SIBYLL [1] is a hadronic interaction model that is widely used in air shower simulations.It is available as one of the standard hadronic interaction models for high energy in the simulation packages AIRES, CORSIKA, CONEX and SENECA.SIBYLL is also used for calculating atmospheric lepton fluxes [2].
The current version of the model is SIBYLL 2.1 [3].It is designed to allow simulation of hadronic interactions in the energy range from √ s ≈ 10 GeV up to 400 TeV.At the time of tuning the parameters of this model, TeVatron data ( √ s ∼ 2 TeV) were the highest energy measurements available.In this work we present a new version of SIBYLL tuned to accelerator data including those from LHC.In addition, this version has been extended to include a phenomenological model of the production of charmed hardons.
In the first section we describe the updates of the model motivated by LHC data.This includes the refit of the cross section parameters, the extension of the fragmentation model to increase baryon pair production, and the update of the parton distribution functions.In the second section we describe the model of charm production, how the parameters are adjusted to describe data and what conclusions can be drawn from applying the model in calculations of the inclusive flux of atmospheric leptons.

LHC updates
Before discussing the update of the model it is worthwhile to mention that SIBYLL 2.1 already describes the general characteristics of hadronic interactions at 7 TeV remarkably well (see dashed blue histogram in Fig. 4 or the review by d'Enterria et al. [4]).
a e-mail: felix.riehn@kit.edulaboratory energy E Lab (eV) Figure 1.Inelastic p-p cross section in SIBYLL.The updated cross section is shown in blue, the old version is in black.The red squares are the measurements by TOTEM [5].The black diamond at the highest energy is the estimate from the Auger Observatory [6,7].The second energy axis shows the equivalent laboratory energy for p-p interactions as applicable to air shower detectors (one proton at rest).The measurement ranges of the IceTop air shower array [8] and the Pierre Auger Observatory [7] are indicated by black lines.

Proton proton cross section
The hadron-proton cross section in SIBYLL follows from unitarizing hard and soft cross section contributions, separated by an energy dependent cutoff in p ⊥ , and terms due to diffraction dissociation.More details on the structure of the model can be found in the publication for version 2.1 [3].Measurements at LHC suggest (see Fig. 1) that SIBYLL 2.1 overestimates the cross section at high energies.The inelastic cross section measured in the TOTEM experiment, which has the highest precision for a measurement of the total cross section at LHC, is 73.5 +1.9 −1.4 mb [5] whereas SIBYLL predicts 80 mb.The rise of the pp cross section beyond 1 TeV is mainly driven by hard parton scattering (hard minijets).
In SIBYLL an eikonal approximation is used to combine the parametrization of soft scatterings with the perturbative calculation of the minijets into an unitary amplitude, which then defines the total and elastic cross sections.The size of the soft and hard contributions in this formalism depends on the size of the particular cross section and the profile function.
In order to make the inelastic cross section compatible with the TOTEM result without changing the hard cross section (calculated within QCD), the profile function of the distribution of the hard partons in transverse (impact parameter) space has been made more narrow so that peripheral collisions are less likely to produce minijets.
The downside of this approach is that central collisions now exhibit very high densities of interacting partons (profile functions are normalized), which means that some events will have a large number of minijets and consequently a large number of final state particles produced.This effect will produce a tail in the multiplicity distribution that is not observed in data.However central collisions are rare so the average multiplicity and most other observables are still compatible with the measurements [4].
Since our goal is a model capable of describing interactions a decade and more higher in center-of-mass energy, the effects of high parton densities have to be considered, even if the mean multiplicity still agrees with current experiments.A microscopic model of parton density saturation could limit the number of scatterings in central collisions and thereby repair the multiplicity.In the current model, saturation is implemented only in an impact parameter independent way as an energy-dependent lower p ⊥ -cutoff for the minijets.
In addition to changing the hard profile function we adjust the parameters of the soft cross section parametrization to fit the p-p and p-p cross section data.
Since the proton profile also enters the meson-nucleon cross sections, we refit the parametrization of the soft contribution there as well.
The resulting cross section is shown in Fig. 1 as a blue line.The old cross section for comparison is shown as a black solid line.The data point of the highest energy is the estimation of the p-p cross section from air shower measurements at the Auger Observatory at energies of about √ s = 57 TeV [6,7].This value has not been used to fit the cross section in SIBYLL and therefore can be seen as an indication that the extrapolation by the model is reasonable.center-of-mass energy s (GeV) laboratory energy E Lab (eV) Figure 2. Average antiproton multiplicity as a function of centerof-mass energy.The low energy data are a compilation of fixed target and ISR experiments that cover the full phase space or were extrapolated to full phase space [9].The CMS data [10] are taken in a phase space region with |y| < 1.0.PHENIX [11] data are taken in the range |η| < 0.35.The prediction by the models are shown for the full and CMS phase spaces only.SIBYLL 2.1 is shown as dashed line, the updated version as solid line.

Baryon production
Particle production in interaction models primarily depends on the implementation of the fragmentation process.
Fragmentation is a non-perturbative process so the rates of particle production cannot be calculated from first principles, which means the parameters in the model have to be set by comparison with experiment.
In SIBYLL string fragmentation [12] is used as the fragmentation model.The string model simplifies hadronization by assuming a uniform energy density in the color field stretched between two partons which eventually is split in two by quark-antiquark pair production.The splitting is continued until the remaining energy is just enough to form two hadrons.Baryons are produced by introducing diquark-antidiquark pairs instead of q q pairs.The probability of producing a diquark pair rather than a quark pair (P q/qq ) in a string breakup is the parameter that controls baryon pair production.In version 2.1 it is set to 0.04 .
For simplicity, only two string classes are distinguished in SIBYLL: the 2 string configuration for the 2 → 2 sea parton scattering and two single strings connecting valence quarks/diquarks.The essential difference between the two is that the latter configuration has valence flavor attached to the string ends, where as the former is in total flavor neutral.This distinction is necessary to describe the differences between the forward/backward regions and the central region of phase space.
The result of this treatment of baryon production in SIBYLL 2.1 for the antiproton multiplicity is shown in Fig. 2 as dashed black lines together with a compilation of data.The multiplicity for full phase space, typically measured in fixed target experiments at low energies, is shown in the upper set of lines whereas the multiplicity in the central region (|η| < 2 ), the region typically accessible in collider experiments, e.g.CMS [10], is shown in the lower set.The current model describes the threshold at low energies well but is not capable of describing the central, high energy data at the same time.
In order to allow for a meaningful extrapolation to high energies, instead of introducing an arbitrary energy dependent parametrization for P q/qq , one can relate the baryon production frequency to soft and semihard interactions, whose energy dependence is different, to increase the baryon production mainly at high energy.With the minijet cross section being derived from perturbative QCD the extrapolation to higher energy is then given by the model and, at low energy, threshold effects due to the large mass of the baryon pairs are important.
Furthermore minijets mostly produce particles in the central region which is exactly where the high energy data by CMS reveal a deficit for SIBYLL 2.1.This assumption is supported by the observation of the ratio of antiprotons to charged pions compared to the central charged multiplicity (see e.g.Fig. 15 in Ref. [10]).
The simplest possible coupling of the diquark production parameter to minijets is to choose a different but fixed value of P q/diq in the fragmentation of minijets in comparison to all other fragmentation processes.The resulting model describes the data much better (solid blue line in Fig. 2), especially in the central region.
Measurements of baryon production at LHC energies that cover the forward phase space could test the assumptions made in this model.

Transverse momentum of minijets
In SIBYLL 2.1 the momentum fractions that determine the kinematics of the minijets are taken from an effective parton density function [13] where g(x) and q(x) are parametrized according to Eichten et al. (EHLQ) [14] and the quark distribution function includes contributions from three light flavors (u, d, and s) and the valence quarks.
In the updated version the same effective parton distribution function (PDF) is used but the quark and gluon contributions are sampled from the same PDF parametrizations (GRV [17,18]) that are used in the calculation of the hard minijet cross section.
The main difference between these parametrizations is the behavior at low x which, in the case of the GRV parametrization, is much steeper.
In combination with the correction of a mistake in the definition of the p min ⊥ the steeper PDFs give a better description of the spectra in the range of intermediate transverse momenta (2 − 5 GeV/c) than in SIBYLL 2.1, see Fig. 3. inv.yield Sibyll 2.3rc1 Sibyll 2.1

Other updates
Other general and more technical aspects of the model that have been updated but are not discussed here are: the transverse momentum acquired in the soft scattering of valence quarks as well as in the string fragmentation is now sampled from an exponential transverse mass distribution rather than a Gaussian as in the previous version.
Another aspect of direct importance to air shower predictions is the enhanced forward production of vector mesons with respect to pseudoscalar mesons (pions) in meson nucleon interactions [19].Since this mechanism has a large influence on muon production in air showers it has been implemented in the new version of SIBYLL.
Furthermore the implemented Glauber model for the calculation of the different cross sections (total, elastic, diffractive, and quasi-elastic) in hadron-nucleus interactions has been extended to include a consistent treatment intermediate low-mass excitations, leading to enhanced screening effects [20].

Comparison to data
To show the compatibility of the updated model with experimental data we look at the charged particle pseudorapidity distribution.The advantage of this observable is that it is very sensitive to the details of the parton level interaction structure and kinematics (n jets , x i ) as well as to the subsequent fragmentation process (dN ch string /dη).The changes introduced in the cross section are expected to increase the central multiplicity at energies beyond 1 TeV whereas the increased baryon production, due to the higher mass of baryonic particles, can be expected to lead to an overall decrease in the multiplicity.Fig. 4 shows that both effects approximately cancel one another at LHC energy and that also the updated model (solid blue histogram) describes the CMS data well.CMS@7TeV CDF@1.8TeVUA5@900GeV UA5@200GeV NA22@20GeV Sibyll 2.3rc1 Sibyll 2.1 Figure 4. Pseudorapidity distribution of charged particles.The data are from NA22 [21],UA5 [22],CDF [23] and CMS [16].The prediction by SIBYLL 2.1 is shown by the dashed line, the one for the updated model by the solid line.

Charm quark extension
SIBYLL 2.1 is limited to the production of particles containing u, d, and s quarks.In version 2.2c of SIBYLL it was shown that a simple phenomenological extension of the fragmentation model, based on the family connection between strange and charmed hadrons, can account for the production of charmed particles at low energy [24].In this approach the normalization is set by the rate at which charm quarks appear relative to strange quarks P c = 0.004 .

New charm model
Due to the high mass of the charm quark the production of charmed hadrons in the fragmentation process is suppressed by a large factor.Instead the dominant channel is the direct production of charm quarks in parton-parton scattering.In the context of QCD the leading contribution gg → cc is often referred to as QCD gluon-gluon fusion [25].The momentum transfer of the reaction due to the charm quark mass Q 2 > Q min ∼ 2m c means that the process can be expected to be calculable within perturbation theory.
The SIBYLL event generator includes only the dominant terms of hard parton-parton scattering at high energy and does not distinguish between the hadronization of the different parton configurations.All parton-parton scattering processes are fragmented into hadrons through an unflavored two string configuration, similar to 2 scattered gluons (usually referred to as hard minijets).
To account for the dominating hard scattering contribution the charm quark fraction is increased in the fragmentation of the hard minijets.To keep the threshold behavior at low energy the charm quark fraction is suppressed exponentially in low mass strings.Specifically where ŝ is the invariant mass of the scattering partons and m eff = 20 GeV 2 is the effective mass scale.To account for string configurations of higher order charm production is not limited to the end of the strings but extends over the whole string.This part of the phenomenological model for charm production is referred to as perturbative component.
Next to the dominant contribution from hard scattering, experiments have shown that there is an asymmetry in charm production in the fragmentation region (i.e. at large x F ) [26,27], which suggests a contribution from charm production in soft interactions.Two models, which can be used to explain this forward production of heavy flavor, are the intrinsic charm model [28] and the flavor excitation model [29].
In SIBYLL we chose a model which could represent either mechanism by adding charm quarks to any string attached to soft scattered partons as well (non-perturbative component).These will include valence strings which, due to the large momentum fraction and the attached flavor of the valence quarks, are able to produce the observed asymmetry at large x F .

Tuning the charm parameters
The values of the parameters in Eq. 2 are adjusted separately for the perturbative and non-perturbative contribution.The perturbative part is tuned to describe the p ⊥spectra of D mesons measured by the ALICE [30] and LHCb [31] experiments in central phase space, since this is where its contribution is expected to be dominant (Fig. 6).The parameters for the soft contribution are set to account for the missing production at low energies (Fig. 7).
In Fig. 5 the cross section for inclusive charm production is shown as a function of the center-of-mass energy.The ALICE data include an extrapolation from central to total phase space.The cross section for D meson production that is measured directly by ALICE is shown by the lower blue points and lines.The dotted line represents the inclusive D meson production cross section without subtracting the decays of resonances of higher mass, e.g.D * .It is shown here because the low energy measurements are not corrected for this either.
The resulting model correctly describes the rise of the inclusive charm cross section with energy and reproduces the spectra at different energies.

Discussion
Charmed particles were introduced to the model because they are expected to contribute significantly to the inclusive flux of atmospheric muons and neutrinos at high energy [39].The inclusive fluxes can be obtained by solving the corresponding cascade equation [40].The results depend on the spectrum weighted moments where x L is the energy fraction of the considered final state particle in the laboratory frame, dn dx L is the hadronic laboratory energy E Lab (eV)  [26,[32][33][34].The measurements at the highest energies are cc from ALICE [30,35].Here data are shown extrapolated to full phase space (red circles) and visible only (blue empty squares).At intermediate energies the data taken at RHIC by the STAR [36] and PHENIX [37] experiments are shown (also extrapolated).The model prediction for the inclusive cc cross section is shown by the solid line, the prediction for the production of D-mesons is shown by the dotted line.production spectrum, and γ is the power law index of the integral all-nucleon spectrum of cosmic rays.With 1.7 < γ < 2.3 it is evident from Eq. 3 that the contribution to the lepton flux is largest for charm production at large x L .  .Feynman-x spectra of charged D mesons in p-p fixed target interactions with P Lab = 400 GeV [32] and 800 GeV [33].
Unfortunately, particle production at large x L is difficult to study at high energy.So far there are no experiments that cover this part of phase space and are capable of particle identification (PID) at the same time.The most forward detector at the LHC with PID capabilities is LHCb (2.5 < y < 4.5).For SIBYLL the contribution from this phase space to Z pD is only about 10 % at √ s = 7 TeV (in-  tegrating the green line in Fig. 8).The prediction of the contribution of charmed particles to the inclusive neutrino and lepton fluxes therefore are not well constrained by the measurements at the LHC.In addition, large-x L production is dominated by soft, non-perturbative processes, so the prediction can not be well constrained by theoretical arguments either.In Fig. 9 a comparison of the weighted energy spectrum, i.e. the integrand in Eq. 3, for D mesons in p-p interactions between the model by Martin, Ryskin and Stasto (MRS) [38] and SIBYLL is shown.The MRS model is a perturbative calculation of the charm production that is extended to low momentum fractions using additional assumptions and accounting for saturation.The energy of the comparison is 7 TeV and the index, with which the spectrum is weighted, is γ = 2.To compare the shapes of the distributions the models are scaled such that MRS and the perturbative component in SIBYLL are equal at x F = 0.19 .One can see that the MRS model and the perturbative component in SIBYLL show a similar behavior.The main difference between the models is that SIBYLL predicts additional charm production from the non-perturbative component that is dominating in the forward direction.
A detailed calculation of the atmospheric lepton fluxes using the model discussed here can be found in our second contribution [41] to this conference.In that paper, the role of the all-nucleon spectrum and the atmosphere are discussed as well.
It should be mentioned that the entire discussion here was focused on proton-proton interactions.What matters for the atmospheric fluxes is the charm production in nucleon-nucleus interactions.In principle, the model is implemented such that central (perturbative) charm production should scale approximately with the number of binary interactions, while forward charm production scales with the number of projectile participants.In practice these scaling expectations are not really satisfied because of additional energy-momentum constraints.Given the strong dependence of the atmospheric fluxes on the forward production nuclear screening effects in this region could have a large effect.For central production, the measurements confirm the binary scaling [42].

Summary and Outlook
An improved version (2.3rc1) of the hadronic interaction model SIBYLL has been presented.The current status of the update of the p-p cross section, the extension of the fragmentation model to describe increased baryon production, and the new charm production model have been described in more detail.The perturbative component of the charm model was found to be compatible with the analytic MRS calculation.It was also shown that the experimental data currently available on charm production do not directly restrict the predictions for the inclusive muon and neutrino fluxes.Only indirectly, by comparing model predictions with charm measurements in phase space regions covered also at colliders, constraints on atmospheric lepton fluxes can be derived.
In the future we plan to estimate how large the uncertainty in the atmospheric fluxes due to the limited phase space coverage of the measurements is.This can be achieved by looking for a set of parameters in the charm model that either minimizes or maximizes the forward charm yield while still being compatible with experimental data.

Figure 3 .
Figure 3. Inclusive cross section for charged particles as function of the transverse momentum.The results obtained with the old and new versions of SIBYLL are compared with CMS data at different c.m. energies [15, 16].

Figure 5 .
Figure 5. Inclusive charm and D-meson cross sections as a function of c.m. energy.The data at low energy are D-meson cross sections in fixed target experiments[26,[32][33][34].The measurements at the highest energies are cc from ALICE[30,35].Here data are shown extrapolated to full phase space (red circles) and visible only (blue empty squares).At intermediate energies the data taken at RHIC by the STAR[36] and PHENIX[37] experiments are shown (also extrapolated).The model prediction for the inclusive cc cross section is shown by the solid line, the prediction for the production of D-mesons is shown by the dotted line.

Figure 6 .
Figure 6.Transverse momentum spectrum of different types of D-mesons in the rapidity interval 4.0 < y < 4.5.Data were taken at √ s = 7 TeV with the LHCb detector [31].

Figure 7
Figure 7. Feynman-x spectra of charged D mesons in p-p fixed target interactions with P Lab = 400 GeV[32] and 800 GeV[33].

Figure 9 .
Figure 9.Comparison of the weighted spectrum of D-mesons in SIBYLL and the MRS model [38].The MRS spectra are shown for different fragmentation assumptions.
Weighted spectrum for D-mesons in SIBYLL at √ s = 7 TeV.The contributions from the perturbative and nonperturbative model components are shown by the blue and red lines, respectively.Note the negligible contribution to the energy spectrum from the phase space covered by the LHCb experiment (2.5 < y < 4.5, green line).