Cosmic ray event generator Sibyll 2.1
Eun-Joo Ahn, Ralph Engel, Thomas K. Gaisser, Paolo Lipari, Todor Stanev
aa r X i v : . [ h e p - ph ] N ov FERMILAB-PUB-09-304-A
Cosmic ray interaction event generator Sibyll 2.1
Eun-Joo Ahn , , Ralph Engel , Thomas K. Gaisser , Paolo Lipari , and Todor Stanev Center for Particle Astrophysics, Fermi NationalAccelerator Laboratory, Batavia, IL 60510-0500, USA Bartol Research Institute, Department for Physics and Astronomy,University of Delaware, Newark, DE 19716, USA Forschungszentrum Karlsruhe, Institut f¨ur Kernphysik,Postfach 3640, 76021 Karlsruhe, Germany and INFN sezione Roma “La Sapienza”,Dipartimento di Fisica Universit´a di Roma I,Piazzale Aldo Moro 2, I-00185 Roma, Italy (Dated: November 1, 2018)
Abstract
The cosmic ray interaction event generator sibyll is widely used in extensive air shower sim-ulations. We describe in detail the properties of sibyll sibyll . INTRODUCTION Cosmic ray interactions in the atmosphere can be regarded as high energy fixed targetcollisions involving heavy particles. Because of their low intensity, cosmic rays with energiesabove 10 eV can only be studied indirectly through the extensive air showers (EAS) theyinitiate in the atmosphere. The analysis of EAS relies on air shower Monte Carlo simulationswhich uses hadronic interaction models. At higher energies, where the cosmic ray energyis beyond the reach of man-made accelerators, hadronic interaction properties have to beextrapolated. The difficulties in the extrapolation are augmented by the fact that, whilethe forward region contains most of the energetics and is important for shower development,most of the accelerator measurements are made in the central region.The event generator sibyll [1] is intended for air shower cascade simulations. It is arelatively simple model that is able to reproduce many features of hadronic interactions infixed target and collider experiments. sibyll is based on the dual parton model (DPM) [2–4],the Lund Monte Carlo algorithms [5, 6], and the minijet model [7–10]. The hard interactioncross section is calculated according to the minijet model. For hadron-nucleus interactions,the interaction probability for each nucleon inside the nucleus is calculated based on theimpact parameter distribution. The total interaction cross section is calculated using theGlauber scattering theory [11]. For a nucleus-nucleus interaction the semisuperpositionmodel [12] is used to determine the point of first interaction for the nucleons of the projectilenucleus. The fragmentation region is emphasized as appropriate for air shower simulations.Versions 1.6 and 1.7 of sibyll have been released and used since the early 1990s. The onlydifference between the two is that version 1.7 can have neutral pion interactions, which isimportant only for air showers above 10 eV because at lower energy all neutral pions decaybefore they interact.Several shortcomings of sibyll sibyll b dependence ( b is the impactparameter) as used for hard interactions. While in version 1.7 the cross section for diffractiondissociation is parametrized independently of the eikonal model, a two-channel eikonal modelbased on the Good-Walker model [14, 15] is used in sibyll sibyll to make a referenceof the implemented physics models and ideas available. We will outline the overall structureand improvements made, within details of the soft interactions and diffraction dissociation.We compare sibyll with fixed target and collider data, and we show how it performs in airshower simulations. Finally, we list some remaining shortcomings of sibyll II. HADRON-HADRON INTERACTIONA. Basic DPM picture sibyll q ,color triplet) and diquark ( qq , color antitriplet). Soft gluons are exchanged in an interactionand the color field gets reorganized. The projectile quark (diquark) combines with the targetdiquark (quark) to form two strings. Each string fragments separately following the Lundstring fragmentation model [6].The fractional energy x of the quark f q ( x ) is chosen from a distribution of f q ( x ) = (1 − x ) α ( x + µ /s ) / , (1)where α = 3 . µ = 0 .
35 GeV is the effective quark mass. The diquark energy fraction3s then f qq ( x ) = 1 − f q ( x ). If particles 1, 2 collide to form strings a and b , the energy andmomentum of the strings are as follows E a = √ s x ,q + x ,qq ) , E b = √ s x ,qq + x ,q ) (2) p a = √ s x ,q − x ,qq ) , p b = √ s x ,qq − x ,q ) (3)To fragment the string, a q -¯ q pair or qq - qq pair is generated at one of the randomly chosenends of the string. The new flavor combines with the existing one to form a hadron, andthe remaining (anti)flavor becomes the new end. A primordial p T of equal magnitude andopposite signs is assigned to the pairs, with a Gaussian distribution where the mean is energydependent h p T i = (cid:20) p + 0 .
08 log (cid:18) √ s
30 GeV (cid:19)(cid:21)
GeV / c , (4)with p = 0 . u, d ), 0.45 ( s ), 0.6 ( qq ) and √ s being the c.m. energy of the hadron-hadroninteraction. The energy fraction of each new particle follows the Lund fragmentation function f ( z ) = (1 − z ) a z exp (cid:20) − b m T z (cid:21) , (5)where a = 0 . b = 0 . m T = p m + p T is the transverse mass, and z is the fractionof the new particle energy with respect to the parent quark or diquark. The fragmentationprocess continues until the remaining string mass is less than a “threshold mass.” Thethreshold mass here is defined as the quark masses of the string ends plus the quark/diquarkpair mass plus (1 . ± .
2) GeV. The string finishes the fragmentation by forming two finalhadrons.
B. Hard interactions and minijets
Already in the range of collider energies, at √ s ∼
100 GeV, the original DPM picture ofjust two strings cannot explain what is observed, namely:1. high multiplicity;2. increase of mean p T ;3. high p T jets;4. rise of central rapidity density. 4hese new features can be interpreted as the emergence of hard interactions which becomeprominent as energy increases, in the form of minijets. Minijets are described with pertur-bative QCD, where partons from the colliding hadrons experience hard scattering. Minijetshave a transverse momentum larger than some momentum transfer scale p min T ≫ Λ QCD ,where perturbative calculation holds, but smaller than a typically reconstructed collider jet.The minijet formalism described below is based on Refs. [1, 16] with modifications made inthe new version.Minijets are perceived as part of the hard interaction described by perturbative QCD.The cross section is calculated within the QCD-improved parton model in leading order is σ QCD ( s, p min T ) = K Z dx Z dx Z dp T × X i,j,k,l
11 + δ k,l f a,i ( x , Q ) f b,j ( x , Q ) dσ i,j → k,l QCD (ˆ s, ˆ t ) dp T Θ( p T − p min T ) , (6)where f a,i ( x , Q ) and f b,j ( x , Q ) are the parton distribution functions of parton i ( j ) inparticle a ( b ). The transverse momentum of the scattered partons is denoted by p T . Thecalculation is done for four light flavors. Higher order corrections are accounted for by settingthe factor K =2 and the factorization scale Q = p T . sibyll g ( x ) ∼ /x at small x . Datafrom HERA [18, 19] suggest a steeper increase at low x . sibyll /x with ∆ = 0 . − .
4. As inthe previous version, Eq. (6) has been calculated separately and is included in the code intabular form.The change in the low- x region affects the minijet cross section substantially at high en-ergies. The cross section cannot rise without limit at high energies [22, 23]. If the number ofgluons times the transverse resolution scale of hard interaction ( ∼ /p T ) becomes compara-ble to the proton size, nonlinear effects, possibly saturation, cannot be neglected. Anotherfactor to take into account is the use of collinear factorization approximation in calculatingminijet cross sections, where the transverse momenta of the incoming partons ( i, j ) shouldalways be smaller than the transverse momenta of the scattered partons ( k, l ). This approx-imation is used to sum the parton densities and only the leading term ln( p T ) is considered.The collinear factorization approximation breaks down for ln(1 /x ) ≫ ln( p T ). The ln(1 /x )term becomes important at high energies and needs to be taken into account [24]. In order5o restrict the calculation of the minijet cross section to the region of phase space where theQCD-improved parton model is expected to be reliable, the following transverse momentumcut is applied p min T ( s ) = p T + Λ exp (cid:20) c q ln( s/ GeV ) (cid:21) , (7)where p T = 1 GeV, Λ = 0 .
065 GeV, c = 0 .
9. This parametrization follows from the geometricsaturation condition [23] α s ( p T ) p T · xg ( x, p T ) ≤ πR p , (8)where α s is the strong coupling constant, g ( x, p T ) the gluon density, and R p is the ef-fective radius of a proton in transverse space. The scale Q = p T is assumed. In thelimit ln(1 /x ), ln( Q / Λ ) → ∞ (double leading-logarithmic approximation) the steeply ris-ing gluon density g ( x, Q ) can be written xg ( x, Q ) ∼ exp " − n f ln ln Q Λ ln Q Λ ln 1 x ∼ x . , (9)with Λ being the QCD renormalization scale and n f is the number of quark flavors. Thefunctional form of Eq. (7) follows from inserting Eq. (9) in Eq. (8), however, the parametersin Eq. (7) cannot be derived directly from first principles. sibyll p min T = √ p min T values for both versions are shown in Fig. 1.The minijet cross section quickly rises to exceed the total cross section. This is interpretedas the collision forming more than one minijet. The average number of hard interactions n hard occurring at energy s and at impact parameter b is [9] n hard ( b, s ) = A ( b ) σ QCD ( s ) , (10)where A ( b ) is the profile function for the hadron-hadron collision. The baryon and mesonprofile functions follow those of Refs. [10, 25] and are given in Appendix B. Following theconvention given in Ref. [26], where the basic equations are in Appendix A, the inelasticcross section is σ inel = Z d b (cid:2) − e − χ ( b,s ) (cid:3) , (11)6 p T m i n ( s ) ( G e V ) E CM (GeV) FIG. 1: Minimum transverse momentum ( p min T ) required for the collision to qualify as a hardscattering. sibyll √ √ s in version 2.1. where the eikonal is χ ( b, s ) = χ hard ( b, s ) + χ soft ( b, s ) = 12 n hard ( b, s ) + 12 n soft ( b, s ) . (12)The number of soft interactions is defined analogously to the hard one n soft ( b, s ) = A soft ( b ) σ soft ( s ). The soft part of the eikonal is discussed in the following subsection.The hard part of the eikonal is interpreted as having a probability of exp[ − n hard ( b, s )] forno minijet production at energy s and impact parameter b . Equation (11) can be reorganizedas σ inel = Z d b (cid:8) − e − n hard ( b,s ) + e − n hard ( b,s ) − e − n hard ( b,s ) e − n soft ( b,s ) (cid:9) = ∞ X N =1 σ N + Z d b e − n hard ( b,s ) (cid:2) − e − n soft ( b,s ) (cid:3) , (13)where σ N = Z d b n hard ( b, s ) N N ! e − n hard ( b,s ) (14)is the cross section for production of N pairs of minijet. This interpretation follows fromthe Abramovski-Gribov-Kancheli [27] cutting rules, where σ N is the term with exactly N cut parton ladders, summed over all uncut ladders [28]. The probability distribution for7 -4 -3 -2 -1 E m i n ij e t / E C M E CM (GeV) -2 -1 < N m i n ij e t > E CM (GeV) FIG. 2: Minijet production, with the energy fraction carried by the minijets (left panel) and theaverage number of minijets produced (right panel) over a range of √ s . obtaining N minijet pairs is P N = σ N σ inel , (15)and the mean number of minijet pairs produced per interaction is h N i = ∞ X N =0 N P N = σ QCD σ inel . (16)The contribution of minijets to the overall particle production for a p - p collision is shownin Fig. 2, where the energy fraction carried by the minijets and mean number of minijetsproduced are shown as a function of c.m. energy. Minijets start becoming important at √ s ≈ q -¯ q pair is generated at each end, and a leading particle at each end iscreated. Then the string fragments in the standard way. The fraction of energy going intothe minijet from each hadron 1, 2 ( x an x ) is obtained by selecting x from the effectiveparton density function [29] f ( x ) = g ( x ) + 49 [ q ( x ) + ¯ q ( x )] . (17)The code uses the approximation that for small transfer momentum the cross sections for g - g , q - g and q - q scattering are proportional to t − and are in ratio 1 : 4 / / . The8ransverse momentum follows dσd ˆ t ∝ t , (18)where ˆ t is the four-momentum transferred squared Mandelstam variable, and ˆ t > (cid:0) p min T (cid:1) .We emphasize that the full parton structure functions of the u, d, s, c quarks and gluon areused for the calculation of the hard cross section. The above approximation is made onlywhen sampling partonic final states and the parton density is parametrized in this simpleway by adding quarks and gluons with the approximate weights. C. Soft interactions sibyll χ soft = CA ( b ), having an impact parameter profile function identical to that of the hardcounterpart, and C = 123 GeV − is chosen to reproduce the low energy inelastic cross sectionof 32 mb. Only one soft interaction is permitted, and the hard-soft interaction division wasenergy-independent at p T = √ p min T ( s ) allows a larger range of phase space for softinteractions. The eikonal form of χ soft = A soft ( b ) σ soft ( s ) is kept. We adopt some aspects ofRegge theory in order to accommodate multiple soft interactions. Inspired by Regge theory[30] the energy dependence of σ soft is taken as sum of two power laws, one term for Pomeronexchange and another term for Reggeon exchange [31] σ soft ( s ) = X (cid:18) ss (cid:19) ∆ eff + Y (cid:18) ss (cid:19) − ǫ , (19)The index ǫ for Reggeon exchange at low energy is expected to be very similar to the onefound in fits by Donnachie and Landshoff [31]. The parameter ∆ eff , in contrast, depends onthe subdivision of the Pomeron term into soft and hard contributions and is hence a functionof the transverse momentum cutoff (7). Here we implicitly assume that minijets form thehard part of the Pomeron [32].The parameters X , Y and ǫ and ∆ eff are determined by fitting the measured total, elasticand inelastic cross sections for p - p and p - ¯ p interactions. Based on the GRV parton densities[20, 21] ǫ ≈ . eff ≈ . b b b b FIG. 3: Geometric interpretation of a hard interaction (left) and soft interaction (right).
In the following we describe soft interaction by extrapolating the picture of hard interac-tions into the domain of low momentum transfer, in which this picture cannot be justifiedwithin perturbative QCD. Only comparison of the corresponding model predictions withdata can later prove this description as useful.While hard interactions are approximately pointlike in character, soft interactions involvea larger transverse interaction area. The low p T of the partons in a soft interaction spreadsout the interaction area (from the uncertainty principle ∆ p ∆ b ∼ p T event which localizes the collision to a small region. The geometry of hard and soft inter-actions in impact parameter space is schematically shown in Fig. 3. The energy-dependentincrease of the interaction region of soft interactions has to be taken into account in the softprofile function of the hadron-hadron collision A soft yz ( s, b ) = Z d b d b d b A y ( b ) A z ( b ) A soft ( s, b ) δ (2) ( b − b + b − b ) , (20)where the profile function A y/z is analogous to the hard interaction case. The fuzzy area ofthe right hand diagram in Fig. 3 is represented by A soft ( s, b ), which is parametrized as aGaussian with a energy-dependent width A soft ( s, b ) = 14 π B s ( s ) exp (cid:20) − | b | B s ( s ) (cid:21) , (21)with B s ( s ) = B + α ′ (0) ln (cid:18) ss (cid:19) . (22)In the limit of B s ( s ) → A soft ( s, b ) becomes a delta function and Eq. (20) becomes equalto the hard profile function. 10n order to calculate Eq. (20), an exponential form factor which corresponds to a Gaussianshape in transverse space is used for the proton or meson profile function A y ( b ) as a firstestimate, using data to fit the parameters. For a proton-proton collision, this yields A soft pp ( s, b ) = 14 π (2 B p + B s ( s )) exp (cid:20) − b B p + B s ( s )) (cid:21) , (23)where B p characterizes the transverse size of a proton and is fitted to data.For generating the string configurations in inelastic events, the number of soft ( N s ) andhard ( N h ) interactions is sampled from (see also [33]) σ N s ,N h = Z d b [ n soft ( b, s )] N s N s ! [ n hard ( b, s )] N h N h ! e − n hard ( b,s ) − n soft ( b,s ) , (24)with the inelastic cross section given by σ inel = X N s + N h ≥ σ N s ,N h . (25)The probability distribution Eq. (24) is tabulated during initialization of sibyll and laterused to draw event configurations.In analogy to hard interactions, soft interactions are simulated by a pair of gluons whichare fragmented the same way as a minijet pair. Thus, there is one valence quark string pairand n s − /x distribution. This distribution corresponds to the one expected for a scenario ofsaturated gluon density in central collisions. A minimum mass of m soft = 1 GeV is requiredfor strings between sea quarks to regularize the singular part of the distribution and toensure applicability of string fragmentation. One of the multiple soft interactions alwaysinvolves the valence quarks and the momentum fractions are then sampled from Eq. (1).Implementing multiple soft interactions affects the model predictions at intermediateenergy, which can be seen, for example, in the inelastic and total cross sections between √ s = 50 −
900 GeV in Fig. 4.
D. Diffraction dissociation
Diffraction is a collision where there are no quantum numbers exchanged between thecolliding particles. A characteristic feature is a large rapidity gap in the final state. Unfortu-nately diffraction physics is not satisfactorily understood even on the level of phenomenology.11 comprehensive description of diffraction can be found in Ref. [34]. We give only a briefdescription of the diffraction model used in the new version and put the relevant equationsin Appendix C.In version 1.7, diffraction was considered part of the inelastic, no-minijet event but wasnot otherwise included within the physics framework. The cross sections σ diff were simply 9%each for forward and backward diffraction and 4% for double diffraction of the σ inel (includingminijet production) at 30 GeV, and was assumed to increase with energy as σ diff ∝ ln( s ),with the diffractive event probability P diff = σ diff /σ inel . With this treatment, however, σ diff becomes larger than the cross section for inelastic no-minijet events at high energy. Also,the minijet cross section is a fit to the total inelastic cross section, which includes diffractiveevents as well. This resulted in an underestimation of minijet production. sibyll ⋆ ) which is a superposition of different low-mass diffractive final states. The cross sectionsand relevant equations are given in Appendix C. The high-mass diffraction dissociation iscalculated by extracting the corresponding eikonal function from the soft and hard eikonalparts. Including the excited state of projectile and target in the eikonal formalism gaveconsiderable improvement to the multiplicity distribution (see. Sec. III A).Diffraction dissociation is treated with a strict kinematic cutoff M x /s < .
1, which followsfrom considerations on coherence and diffractive particle production [36]. The net effect isthat the quasielastically scattered protons do not lose more than ∼
20 % at maximumin diffraction dissociation at this energy. The diffractively dissociating particle undergoesa phase-space decay if the mass of the excited system is very low. For higher masses,diffracted particles are divided into two valence components of quark-diquark or quark-antiquark that are connected by a color string which subsequently fragments. The stringcarries the diffracted particle’s momentum and quantum numbers and does not create extra p T . The one-string decay threshold is set to ∆ M = 0 . M is the massdifference of the incoming particle and the excited state of it.Diffractively excited states of a mass of more than 10 GeV are considered as being pro-duced by a Pomeron-hadron interaction. The decay of these states is described with multiple12oft and hard interactions by generating a π - p interaction at √ s = ∆ M , as motivated bydata from the UA4 Collaboration [37]. E. Nucleus interactions
The physics framework for the hadron-nucleus and nucleus-nucleus interactions remainsthe same in the new version. The difference in the cross sections between the versions comesfrom the improvements made to the hadron-hadron interactions. Detailed descriptions canbe found in Ref. [1] for the hadron-nucleus interaction and Ref. [12] for the nucleus-nucleusinteraction. Here, we summarize the basic concepts for the sake of completeness.The interaction length is calculated from the production cross section. Elastic andquasielastic interactions, where no new particles are produced, do not contribute to theair shower development and are not considered. The production cross section, i.e. contribu-tion to particle creation, for a hadron-nucleus ( hA ) interaction where A is the mass numberof the nucleus is σ hA prod = σ hA tot − σ hA el − σ hA qe , (26)with σ hA tot , σ hA el , and σ hA qe being the total, elastic, and quasielastic cross sections, respectively.They are calculated within the Glauber model [11] from the p - p , π - p , and K - p cross sections.The inelastic and total cross sections of p - p , p -air and π -air, π - p collision are shown in Fig. 4.The minijet cross sections are also shown.In a hadron-nucleus interaction, the number of target nucleons directly participating inthe interaction, also known as wounded nucleons, is determined from the production crosssection. The mean number h N w i of wounded nucleon per interaction is given by standardGlauber theory. In analogy with Eqs. (15) and (16), one finds h N w i = 1 σ hA prod X N w N w σ N w = A σ hp inel σ hA prod , (27)where σ N w is the cross section for interaction with N w nucleons and σ hp is the hadron-nucleoncross section.The string model is applied to the fragmentation of the partonic system as well. Thetarget nucleus is seen as N w pairs of valence q - qq , and the projectile hadron is viewed as onevalence q -¯ q or q - qq pair and N w − q -¯ q pairs. The N w color-connected partons undergostring fragmentation. Most of the energy is carried by the valence pair string, and the13 σ ( m b ) E CM (GeV) σ tot,p-air σ inel,p-air σ tot,p-p σ inel,p-p σ ( m b ) E CM (GeV) σ tot, π -air σ inel, π -air σ tot, π -p σ inel, π -p FIG. 4: The total and inelastic cross section of p -air, p - p (left panel) and π -air, π - p (right panel)collision for version 2.1 (red solid lines) and 1.7 (blue dotted lines). The minijet cross section isalso shown in dashed lines, where the higher one is version 2.1. sea strings contribute to giving the proper multiplicity in the target region. This happensbecause there are more valence quarks on the nucleus side for the sea quarks to coupleto. Intranuclear interactions and Fermi momentum inside the target nucleus are ignoredin sibyll . Intranuclear interactions are greatly suppressed at high energy due to the timeneeded for a hadron to form as an independent object (formation time). At high energies,minijets are added to the collision, by generating them for each of the N w wounded nucleonsin the same way as in a hadron-hadron interaction.The nucleus-nucleus interaction is treated with the semisuperposition model [12], whichis between the simple superposition model and full Glauber theory [11]. The superpositionmodel treats each nucleon of the projectile independently and as a consequence the interac-tion lengths of the nucleons have an exponential distribution based on the hadron-nucleuscross section. In reality, the nucleus interaction length is very small and a nucleus will inter-act quickly in the atmosphere. In the semisuperposition model, the number of interactingnucleons in the projectile for each nucleus interaction is determined from Glauber theory,where the remaining spectator nucleons fragment into lighter nuclei. Though the interac-tion and fragmentation is treated as a nucleon-nucleus interaction, the distribution of thesenucleon-subshowers reflects correctly the nucleus-air cross section.14 II. COMPARISON WITH EXPERIMENTAL DATA
Both fixed target and collider experiments give valuable guidance in modeling hadronicinteractions. Fixed target experiments provide data for the forward region which are mostrelevant to cosmic ray interactions, but the energies are relatively low, E lab ∼ several hundredGeV. Collider experiments can probe higher energies ( E lab ∼ GeV) but most of theinformation is collected in the central region. Some collider experiments such as H1 andZEUS are able to detect forward events [38]. The anticipated experiments LHCf [39] andTOTEM [40] in the LHC at CERN are expected to collect information in the forward regionat an energy equivalent to cosmic rays of E lab ∼ GeV.
A. Charged particle multiplicities
Pions are the most numerous particles, followed by kaons and baryons. There is overallgood agreement with experimental data in the forward region at low energies. The differenceof charged particle production between the two versions is due to the improved treatmentof multiple soft interactions, usage of GRV parton densities, and a consistent inclusion ofdiffraction dissociation, which also leads to more minijet production. These improvementsgive a better agreement with data for version 2.1, especially in the central region.The NA49 experiment measured the rapidity y and Feynman x F distribution of chargedparticles for p - p [41] and p - C [42] collisions at E lab = 158 GeV. Figure 5 shows the sibyll results compared to the data for π + and π − , which became available only after the eventgenerator had been released. Good agreement between model predictions and data is found.The excess of π + over π − is due to the flavor content of the proton ( uud ). In version 2.1, thisdiscrepancy is stronger and more particles are produced in the central region which reflectsthe changes made to the soft interaction. However, the difference between the two versionsis small in this respect.Fixed target experiments at FNAL used π + , K + , p as projectiles and p, C for targets.The inclusive cross section Ed σ/dp for each charged particle species has been measured at E lab = 100 GeV at a given p T [43]. The results for π + π − production at p T = 0 . π + s are overproduced while15 -2 -1 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 dn / d x F x Fpp -> π + pC -> π + (x 10) 1.72.1 10 -2 -1 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 dn / d x F x Fpp -> π - pC -> π - (x 10) 1.72.100.20.40.60.811.2 -4 -2 0 2 4 dn / d y y π + 1.72.1 00.20.40.60.811.2 -4 -2 0 2 4 dn / d y y π - 1.72.1 FIG. 5: The Feynman x ( x F ) and rapidity ( y ) distribution of pions plotted against NA49 result of p - p [41] and p - C [42] collision at E lab = 158 GeV. Version 2.1 (1.7) results are shown in red solid(blue dotted) lines. The left (right) panels show the production of π + ( π − ). The upper panelsshow the x F distribution, where the p - C collision results are multiplied by factor 10 in order toshow both interactions on the same plot. The lower panels show the y distribution: the upper(lower) set of lines and data points are from the p - C ( p - p ) collision. π − s are underproduced in the forward region, indicating the role of the valence componentsin determining the leading particle. The kaon projectile shows an overall good agreementfor version 2.1. For the proton projectile, there is an overproduction of π − s compared todata, especially in the forward region.Pseudorapidity ( η ) distributions of charged particles from collider experiments are com-16
0 10 20 30 40 50 60 70 80 90 100 E d σ / dp ( m b / G e V ) p z (GeV) π + p -> π + π + C -> π + -1
0 10 20 30 40 50 60 70 80 90 100 E d σ / dp ( m b / G e V ) p z (GeV) K + p -> π + K + C -> π + -2 -1
0 10 20 30 40 50 60 70 80 90 100 E d σ / dp ( m b / G e V ) p z (GeV) pp -> π + pC -> π + -1
0 10 20 30 40 50 60 70 80 90 100 E d σ / dp ( m b / G e V ) p z (GeV) π + p -> π - π + C -> π - -1
0 10 20 30 40 50 60 70 80 90 100 E d σ / dp ( m b / G e V ) p z (GeV) K + p -> π - K + C -> π - -1
0 10 20 30 40 50 60 70 80 90 100 E d σ / dp ( m b / G e V ) p z (GeV) pp -> π - pC -> π - FIG. 6: Inclusive cross sections of π + - p, C (left panels), K + - p, C (center panels), p - p, C (rightpanels) collisions for π + (upper panels) and π − (lower panels) production, where the events have p T = 0 . E lab = 100 GeV, and are compared with results from a FNALfixed target experiment [43]. Version 2.1 (1.7) results are shown in red solid (blue dotted) lines. pared with sibyll in Fig. 7. It shows the η distribution of charged particles from p - ¯ p collisions at E c . m . = 1800 GeV (CDF [44]), 630 GeV (P238 [45]), 200 GeV (UA5 [46]) and53 GeV (UA5 [47]). The improvements made to version 2.1 most prominently show in thecentral region. The role of the minijets and soft interactions is visible in the central region,where version 1.7 lacks secondary particles especially as the energy increases, while hav-ing more particles in the peripheral region. This trait can be seen at low energies in the pp → π + , π − figures in Fig. 6. Version 2.1 gives an excellent description of P238 data andtends to slightly overestimate the particles at low energies. It should be noted that the η range and trigger condition for 53 GeV is different than for higher energies at UA5. The twoversions are similar for events with large | η | beyond the scope of current collider detectormeasurements.The distributions of charged particle multiplicity at UA5 [48] also give information athigher energies. Figure 8 shows the distribution of charged particle multiplicity for p - ¯ p dn c h / d η η dn c h / d η η FIG. 7: Pseudorapidity ( η ) distribution of p -¯ p collision, for various c.m. energies, for version 1.7(left panel) and 2.1 (right panel) Data points from top are E c . m . = 1800 GeV [44], 630 GeV [45],200 GeV [46], 53 GeV [47]. Note the deficit of charged particles in version 1.7 at high energies inthe central region. -5 -4 -3 -2 -1
0 20 40 60 80 100 P n c h n ch | η | < 0.5 (x0.01)| η | < 1.5| η | < 3.0 (x100) FIG. 8: Probability of charged particle multiplicity ( n ch ) distribution at UA5 detector with E c . m . =900 GeV for three different η ranges; | η | < . × . | η | < . | η | < × K s . E c . m . = 900 GeV, at three different η ranges | η | < . , . , .
5, where theresults for | η | < . | η | < . K s has a very short lifetime and its decay produces charged particlesthat have a non-negligible effect, especially at high multiplities. The treatment of K s isconsidered as one of the uncertainties in the interpretation of the experimental data. Bothstable and unstable cases are plotted, which can be considered as an error band. Theimprovements made in soft interaction and diffraction in version 2.1 are evident in the widerdistribution of n ch as well as in the increase in multiplicity. An underestimation of the crosssection for double diffraction dissociation in sibyll is probably the reason for the lack oflow-multiplicity events satisfying the UA5 trigger condition. B. Leading particle
The leading particle from the fragmentation carries a significant fraction of the totalenergy and becomes the primary particle in the next interaction of the air shower. Theelasticity K = E lead /E proj , the fraction of the leading particle with respect to the collisionenergy, of a collision affects the multiplicity as well as the speed of shower development inthe atmosphere. Thus it is important to get a correct description of the behavior of theleading particles.The NAL bubble chamber experiment has data for p - p interactions at E lab = 102, 205,303, 405 GeV and measured the x F of the leading proton [49]. Figure 9 shows the sibyll results plotted against the NAL data. The sharp dip at x F ≈ . x F = 0 . sibyll cannot have a photon or positron projectile,we simulated a p - p collision at a slightly lower energy of E c . m . = 210 GeV and used eventsfrom one hemisphere, i.e. events with p c . m .z >
0. Figure 10 shows the leading protons of the sibyll results plotted against the ZEUS data. They are plotted as a function of x lab = E/E p ,19 d σ / d x F ( m b ) x FE lab = 405 GeV (x1000)E lab = 303 GeV (x100)E lab = 205 GeV (x10)E lab = 102 GeV 1.72.1 FIG. 9: Distribution of leading protons as a function of x F compared to the p - p collision data fromthe NAL bubble chamber [49]. Version 2.1 (1.7) results are shown in red solid (blue dotted) lines.The results and datapoints have been multiplied by factor of tens (405 GeV by 1000, 303 GeV by100, 205 GeV by 10, 102 GeV by 1). the energy of the proton or neutron divided by the colliding proton energy in the lab frame,which is essentially the elasticity. The leading proton displays similar behavior to that of theNAL bubble chamber. Again, the better diffraction treatment is evident around x lab = 0 . dn / d x l a b x lab p FIG. 10: Leading proton distribution of sibyll compared to the low- Q [beam pipe calorimeter(BPC)] and higher- Q [deep inelastic scattering (DIS)] ZEUS data [50, 51]. Version 2.1 (1.7)outcome is shown in red solid (blue dotted) lines. The transverse momentum of the protons arelimited to p T < . . The c.m. collision energy of sibyll is 210 GeV. C. Strange particle production
A FNAL fixed target experiment measured the production of very forward strange par-ticles produced in p - Be collisions at E lab = 300 GeV [52]. The inclusive cross section of Λ and K s have been plotted for angles in the range θ = 0 . − . p z are likely to be from par-ticle production associated with the projectile, and particles with small p z are from centralproduction. Both sibyll versions give agreements in the forward direction, with a tendencyto slightly overproduce high- p z particles and underproduce low- p z particles. Strange particleproduction directly affects production of high energy muons and neutrinos. IV. AIR SHOWER PERFORMANCE
The development of an air shower depends on a number of factors, some of which arethe production cross section, inelasticity, and multiplicity. For the description of the air21 -3 -2 -1
50 100 150 200 250 300 E d σ / dp ( m b / G e V ) p z (GeV) Λ -3 -2 -1
50 100 150 200 250 300 E d σ / dp ( m b / G e V ) p z (GeV)K s0 1.72.1 FIG. 11: Inclusive cross sections of Λ (left panel) and K s (right panel) production from p - Be collision at E lab = 300 GeV for various angles [52]. Version 2.1 (1.7) results are shown in red solid(blue dotted) lines. From the top, the sets of lines and data are for angles 0.25, 1.5, 2.9, 5.0, 6.9,8.8 mrad. The results and data have been multiplied by factors of tens to show them in one figure(0.25 mrad by 10 , 1.5 mrad by 10 , 2.9 mrad by 10 , 5.0 mrad by 100, 6.9 mrad by 10, 8.8 mradby 1). shower development the hadronic interaction model has to describe correctly the particleinteractions in a wide range of energies. Observables such as depth of shower maximum X max ,electron number N e and muon number N µ at ground will depend on the characteristics ofhadronic interactions . We briefly summarize how air showers are affected by those threeparameters.Increasing the cross section will cause the shower to start earlier in the atmosphere, re-sulting in a smaller X max as well as a smaller fluctuation. The number of electrons measuredis highly dependent on the position of the shower maximum; the closer to X max the larger Extensive studies have been carried out in Ref. [53], where cross sections, elasticity, and multiplicity arevaried to see the effect on observables such as X max , N e and N µ .
400 450 500 550 600 650 700 750 800 850 900 15 16 17 18 19 20 X m ax ( g / c m ) log (E/[eV]) Sibyll 2.1Sibyll 1.6QGSJET II-3EPOS 1.61
FIG. 12: The mean X max of protons and iron nuclei as a function of primary energy at zenithangle θ = 45 ◦ . In the case of sibyll (red solid lines), qgsjet II-3 (orange dot-dashed lines), and epos − − . For sibyll − − . the N e . As most muons are produced from decay of pions and kaons, N µ is expected toremain stable.An increase in the mean multiplicity lessens the energy per particle of the secondaries,which results in a quicker development of the shower with smaller fluctuations. The in-creased multiplicity is expected to increase N µ . The number of electrons are most numerousat X max and decreases away from it. Hence, a quicker development results in a larger dis-tance between X max and ground which will contribute to decreasing the number of electronsdetected on ground.A large elasticity gives a larger fraction of the total energy to the leading particle, andslows the shower development as well as giving a smaller multiplicity. The X max will becloser to the ground and N µ will be smaller. The closer proximity of X max to the groundcauses a larger N e despite the smaller multiplicity [53].23he new version has the shower developing more quickly than the old version. The twoversions are compared with qgsjet II-3 [58] and epos sibyll X max overall but maintains the same shape. qgsjet II-3 has a verylarge multiplicity at high energies, resulting in smaller shower maximum.The muon number is an important indicator for cosmic ray composition studies, as show-ers of heavier nuclei contain more muons than that of protons. Version 2.1 produces moremuons than version 1.7, which cause the new version of sibyll to extract a lighter cosmicray composition from experimental data than the old one. However, both versions producefewer low energy muons than other models such as qgsjet
II-3 or epos N µ decreases as increasing zenith angle.At sufficiently high energies above ∼
10 GeV, most muons do not decay and N µ increases.Figure 13 shows the energy dependence of the average number of muons normalized to theprimary energy for the two versions of sibyll . The average number of muons at sea level( h N µ i ) is plotted with energies above E thr µ =0.3, 1, 3, 10 and 30 GeV in proton-initiatedshowers at zenith angle θ = 0 ◦ . V. SHORTCOMINGS AND FUTURE IMPROVEMENTS
Though many features have been improved, sibyll • The current nucleus-nucleus collision uses the semisuperposition model. Implementingthe full Glauber model will give a more accurate description. • Antibaryon production is not satisfactorily described. There is not enough producedand the distribution of antiprotons in the central region is incorrect compared withdata [61]. The overall normalization can be improved by increasing the diquark produc-tion fraction. The current method of fragmentation suppresses antibaryon formationtogether with other particles in the non-end of strings.24 E [eV]0.0010.0040.0070.010 N µ / E [ G e V ] − E µ thr =0.3 GeV E µ thr =3 GeV E µ thr =30 GeV(c)(b)(a) FIG. 13: Average number of muons at sea level h N µ i , obtained in proton showers with zenithangle θ = 0 ◦ . Each energy represents 5000 showers simulated with the hybrid method. The solid(dotted) line represents the values obtained with sibyll sibyll • The currently used energy-dependent transverse momentum cutoff is independent ofthe relevant gluon density of the interaction with the target nucleus, which varies withthe impact parameter of the collision. A new energy- and impact parameter-dependentcutoff to p T to prevent parton density saturation would improve the modeling. Refer-ence [62] is an example of an attempt to constrain minijet formation. In addition theprofile functions could also be more refined. • A consistent treatment of coherent and incoherent diffraction dissociation in hadron-nucleus and nucleus-nucleus interactions is required. This can be achieved by usinga two-channel model in the Glauber calculation similar to the one presented here for25 - p interactions. • Include charm quark. The current model has u, d s quarks and gluons. This im-provement will be relevant more to muon and neutrino detectors than cosmic rayobservations. A version containing charm will be released soon.
VI. CONCLUSION
In this paper, we described the overall model of sibyll together with the changes madein the version 2.1, and listed the shortcomings and possible ways of improving the model.The 2.1 version still keeps the DPM-minijet structure but with modifications and additions.Results from HERA suggest a steeper parton density for gluons at low x , and the GRVparton densities replaced the EHLQ parton densities which resulted in an increased QCDcross section. Concepts from Regge theory are used to allow multiple soft interactions. Theenergy-dependent transverse momentum cutoff for ensuring perturbative QCD to be validis a better discriminator between soft and hard interactions than the previously appliedenergy-independent cutoff. The two-channel eikonal model is used to describe diffraction,with better if not complete success. The physics framework for the hadron-nucleus andnucleus-nucleus interactions remains unchanged.These improvements produce more particles with a wider distribution in momentumspace, as are evident when comparing with experimental data. Both versions give a goodfit to the rapidity and Feynman x distribution for fixed target experiments such as NA49.However the changes made to the new version are very evident in the central region, in thepseudorapidity and the overall charged particle multiplicity distribution. As a consequence,air showers described with sibyll sibyll sibyll sibyll in comparison with other models would be valuable in analysis of cosmic ray airshower data. Acknowledgments
Work on this project at the University of Delaware is supported in part by a grant fromthe Office of Science of the U.S. Department of Energy, DE-FG02-91ER40626. Fermilab isoperated by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 withthe United States Department of Energy. [1] R. S. Fletcher, T. K. Gaisser, P. Lipari, and T. Stanev, Phys. Rev.
D50 , 5710 (1994).[2] A. Capella and A. Krzywicki, Phys. Rev.
D18 , 3357 (1978).[3] A. Capella and J. Tran Thanh Van, Z. Phys.
C10 , 249 (1981).[4] A. Capella, U. Sukhatme, C.-I. Tan, and J. Tran Thanh Van, Phys. Rept. , 225 (1994).[5] H.-U. Bengtsson and T. Sjostrand, Comput. Phys. Commun. , 43 (1987).[6] T. Sjostrand, Int. J. Mod. Phys. A3 , 751 (1988).[7] T. K. Gaisser and F. Halzen, Phys. Rev. Lett. , 1754 (1985).[8] G. Pancheri and Y. Srivastava, Phys. Lett. B159 , 69 (1985).[9] L. Durand and P. Hong, Phys. Rev. Lett. , 303 (1987).[10] L. Durand and H. Pi, Phys. Rev. D38 , 78 (1988).[11] R. J. Glauber and G. Matthiae, Nucl. Phys.
B21 , 135 (1970).[12] J. Engel, T. K. Gaisser, T. Stanev, and P. Lipari, Phys. Rev.
D46 , 5013 (1992).[13] R. Engel, T. K. Gaisser, T. Stanev, and P. Lipari,
Proceedings of the 26th InternationalCosmic Ray Conference, Salt Lake City, Utah, 1999 (AIP Melville, NY, 2000) Vol. 1, p. 415.[14] M. L. Good and W. D. Walker, Phys. Rev. , 1857 (1960).[15] R. S. Fletcher, Phys. Lett.
B320 , 373 (1994).[16] T. K. Gaisser and T. Stanev, Phys. Lett.
B219 , 375 (1989).[17] E. Eichten, I. Hinchliffe, K. D. Lane, and C. Quigg, Rev. Mod. Phys. , 579 (1984).[18] H1, C. Adloff et al. , Nucl. Phys. B497 , 3 (1997), hep-ex/9703012.
19] ZEUS, J. Breitweg et al. , Phys. Lett.
B407 , 432 (1997), hep-ex/9707025.[20] M. Gluck, E. Reya, and A. Vogt, Z. Phys.
C67 , 433 (1995).[21] M. Gluck, E. Reya, and A. Vogt, Eur. Phys. J. C5 , 461 (1998), hep-ph/9806404.[22] L. V. Gribov, E. M. Levin, and M. G. Ryskin, Phys. Rept. , 1 (1983).[23] E. M. Levin and M. G. Ryskin, Phys. Rept. , 267 (1990).[24] L. N. Lipatov, Phys. Rept. , 131 (1997), hep-ph/9610276.[25] L. Durand and H. Pi, Phys. Rev. D43 , 2125 (1991).[26] M. M. Block and R. N. Cahn, Rev. Mod. Phys. , 563 (1985).[27] V. A. Abramovsky, V. N. Gribov, and O. V. Kancheli, Yad. Fiz. , 595 (1973), [Sov. J.Nucl. Phys. , 308 (1974)].[28] K. A. Ter-Martirosyan, Phys. Lett. B44 , 377 (1973).[29] B. L. Combridge and C. J. Maxwell, Nucl. Phys.
B239 , 429 (1984).[30] P. D. B. Collins,
An Introduction to Regge Theorie & High Energy Physics, (CambridgeUniversity Press,, Cambridge, 1977).[31] A. Donnachie and P. V. Landshoff, Phys. Lett.
B296 , 227 (1992), hep-ph/9209205.[32] A. Capella, J. Tran Thanh Van, and J. Kwiecinski, Phys. Rev. Lett. , 2015 (1987).[33] P. Aurenche et al. , Phys. Rev. D45 , 92 (1992).[34] V. Barrone and E. Predazzi,
High-Energy Particle Diffraction (Springer, Berlin, Germany,2002).[35] A. B. Kaidalov, Phys. Rept. , 157 (1979).[36] K. Goulianos, Phys. Rept. , 169 (1983).[37] UA4, D. Bernard et al. , Phys. Lett. B166 , 459 (1986).[38] H. Abramowicz and A. Caldwell, Rev. Mod. Phys. , 1275 (1999), hep-ex/9903037.[39] LHCf, O. Adriani et al. , CERN-LHCC-2006-004.[40] TOTEM, V. Berardi et al. , CERN-LHCC-2004-002.[41] NA49, C. Alt et al. , Eur. Phys. J. C45 , 343 (2006), hep-ex/0510009.[42] NA49, C. Alt et al. , Eur. Phys. J.
C49 , 897 (2007), hep-ex/0606028.[43] D. S. Barton et al. , Phys. Rev.
D27 , 2580 (1983).[44] CDF, F. Abe et al. , Phys. Rev.
D41 , 2330 (1990).[45] R. Harr et al. , Phys. Lett.
B401 , 176 (1997), hep-ex/9703002.[46] UA5, G. J. Alner et al. , Z. Phys.
C33 , 1 (1986).
47] UA5, K. Alpgard et al. , Phys. Lett.
B112 , 183 (1982).[48] UA5, R. E. Ansorge et al. , Z. Phys.
C43 , 357 (1989).[49] J. Whitmore, Phys. Rept. , 273 (1974).[50] ZEUS, S. Chekanov et al. , Nucl. Phys. B658 , 3 (2003), hep-ex/0210029.[51] ZEUS, S. Chekanov et al. , Nucl. Phys.
B637 , 3 (2002), hep-ex/0205076.[52] P. Skubic et al. , Phys. Rev.
D18 , 3115 (1978).[53] R. Ulrich et al. , (2009), 0906.0418.[54] T. Pierog et al. , Nucl. Phys. Proc. Suppl. , 159 (2006), astro-ph/0411260.[55] T. Bergmann et al. , Astropart. Phys. , 420 (2007), astro-ph/0606564.[56] D. Heck, G. Schatz, T. Thouw, J. Knapp, and J. N. Capdevielle, FZKA-6019.[57] C. L. Pryke, Astropart. Phys. , 319 (2001), astro-ph/0003442.[58] S. Ostapchenko, Phys. Rev. D74 , 014026 (2006), hep-ph/0505259.[59] T. Pierog and K. Werner, Phys. Rev. Lett. , 171101 (2008), astro-ph/0611311.[60] J. Alvarez-Muniz, R. Engel, T. K. Gaisser, J. A. Ortiz, and T. Stanev, Phys. Rev.
D66 ,033011 (2002), astro-ph/0205302.[61] E735, T. Alexopoulos et al. , Phys. Rev.
D48 , 984 (1993).[62] T. C. Rogers, A. M. Stasto, and M. I. Strikman, Phys. Rev.
D77 , 114009 (2008), 0801.0303.
APPENDIX A: AMPLITUDE CONVENTIONS
The conventions and parameters are adopted from Ref. [26]. The Mandelstam variables s , t , u are used, and k is the c.m. momentum. The c.m. scattering angle θ is related to t by t = − k sin ( θ/ , (A1)and a new parameter is defined − q ≡ t . The scattering amplitude in the c.m. frame[ f c . m . ( s, t )] and the Lorentz invariant scattering amplitude [ M ( s, t )] are related by M ( s, t ) = − π √ s f c . m . ( s, t ) . (A2)The scattering amplitude f c . m . ( s, t ) and impact parameter function a ( s, b ) are Fourier29ransforms of each other via f c . m . ( s, t ) = kπ Z d b e i q · b a ( s, b ) , (A3) a ( s, b ) = 14 πk Z d q e − i q · b f c . m . ( s, t ) . (A4)With this scattering amplitude, the differential elastic cross section can be expressed as dσ el d Ω c . m . = | f c . m . | , (A5a) dσ el dt = πk | f c . m . | , (A5b) σ tot = 4 πk Im f c . m . ( θ = 0) , (A5c)where the last relation used the optical theorem. The elastic slope parameter B el is definedfrom an approximation of the elastic scattering cross section in small t region as dσ el dt = (cid:20) dσ el dt (cid:21) t =0 e Bt . (A6)Using Eqs. (A5a) and (A5b), (cid:20) dσ el dt (cid:21) t =0 = πk (cid:20) dσ el d Ω c . m . (cid:21) θ =0 = π (cid:12)(cid:12)(cid:12)(cid:12) ( ρ + i ) Im f c . m . (0) k (cid:12)(cid:12)(cid:12)(cid:12) = π (cid:12)(cid:12)(cid:12)(cid:12) ( ρ + i ) σ tot π (cid:12)(cid:12)(cid:12)(cid:12) , (A7)where ρ is the real to imaginary ratio of f c . m . (0).When using eikonals, the cross sections follow the usual convention: σ tot = 2 π Z db (1 − e − χ ( b,s ) ) (A8) σ el = π Z db (1 − e − χ ( b,s ) ) (A9) σ inel = π Z db (1 − e − χ ( b,s ) ) . (A10)We neglect the real part of the elastic scattering amplitude for calculating the eikonal func-tions. APPENDIX B: HARD INTERACTION MINIJET CROSS SECTION AND PRO-FILE FUNCTIONS.
Each proton or meson is characterized by a (transverse) density profile function A z ( b ),where z can be p, π, K . The probability for the two partons in particles y and z to collide30s found by integrating over all possible impact parameters b and b for a given impactparameter b of the collision A hard yz ( ν y , ν z , b ) = Z d b d b A y ( ν y , b ) A z ( ν z , b ) δ (2) ( b − b − b ) . (B1)The profile function of a proton is given by A p ( ν p , b ) ≈ π ) Z d k T (cid:18) k T ν p (cid:19) − e i k T · b = ν p | ν p b | K ( | ν p b | ) (B2)with ν p ≈ .
71 (GeV/c) . The profile function of a meson is A m ( ν m , b ) = 1(2 π ) Z d k T (cid:18) k T ν m + η (cid:18) k T ν m (cid:19) + (cid:19) − e i k T · b , (B3)where ν m and η are adjustable parameters. For pions, ν π ≈ .
54 (GeV/c) and η ≈
0. Hencethe profile function for a pp interaction is A hard pp ( ν p , b ) = Z d b ′ A p ( | b − b ′ | ) A p ( | b ′ | )= ν p π
18 ( ν p b ) K ( ν p b ) , (B4)and for a pπ interaction A hard pπ ( ν p , ν π , b ) = ν p π − ζ ) (cid:20) ν p b K ( ν p b ) + ζ − ζ [ K ( ν π b ) − K ( ν p b )] (cid:21) , (B5)where ζ = ( ν p /ν π ) . The profile functions are normalized to Z d b A hard yz ( ν p , ν π , b ) = Z d b A y ( ν m , b ) = 1 (B6) APPENDIX C: DIFFRACTION DISSOCIATION
In low-mass diffraction, the two-channel eikonal model is used [15, 35]. Only two statesare considered here; a nondiffractive state and a generic diffractive state denoted ⋆ . The31iffractive scattering of particles Y and Z can be expressed in the following matrix elements h Y Z | M int | Y Z i = M Born (C1a) h Y Z | M int | Y ⋆ Z i = β Y M Born (C1b) h Y Z | M int | Y Z ⋆ i = β Z M Born (C1c) h Y Z | M int | Y ⋆ Z ⋆ i = β Y β Z M Born (C1d) h Y ⋆ Z | M int | Y ⋆ Z i = (1 − α Y ) M Born (C1e) h Y Z ⋆ | M int | Y Z ⋆ i = (1 − α Z ) M Born (C1f) h Y ⋆ Z ⋆ | M int | Y ⋆ Z ⋆ i = (1 − α Y ) (1 − α Z ) M Born . (C1g)The coefficients α and β may depend on energy. A matrix ˆ χ ( s, b ) for the eikonal χ isintroduced to calculate M . The eikonal matrix is diagonalized, and the cross sectionscalculated. The hadronic states Y and Z are defined as | Y, Z i ∼ , | Y ⋆ , Z i ∼ , | Y, Z ⋆ i ∼ , | Y ⋆ , Z ⋆ i ∼ . (C2)The eikonal matrix readsˆ χ ( s, b ) = β Y β Z β Y β Z β Y − α Y β Y β Z β Z (1 − α Y ) β Z β Y β Z − α Z β Y (1 − α Z ) β Y β Z β Z (1 − α Y ) β Y (1 − α Z ) (1 − α Y ) (1 − α Z ) χ ( s, b ) . (C3)After diagonalizing ˆ χ ( s, b ), the cross sections can be calculated. The total cross section isgiven by σ tot Y Z = 2 Z d b h Y Z | (cid:0) − e − ˆ χ ( s, b ) (cid:1) | Y Z i = 2 Z d b ∞ X n =1 f el ,nY f el ,nZ ( − n − ( χ ( s, b )) n n ! , (C4)where f el ,nj = (1 − α j γ j )(1 − α j − γ j ) n + (1 + α j γ j )(1 − α j + γ j ) n , (C5)32nd γ j = q α j + β j j = Y, Z . (C6)Consequently the elastic cross section reads σ el Y Z = Z d b (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∞ X n =1 f el ,nY f el ,nZ ( − n − ( χ ( b , s )) n n ! (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (C7)The cross section for single diffraction dissociation of particle Y follows from σ SD ,YY Z = Z d b (cid:12)(cid:12) h Y ⋆ Z | (cid:0) − e − ˆ χ ( b ,s ) (cid:1) | Y Z i (cid:12)(cid:12) = Z d b (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∞ X n =1 f diff ,nY f el ,nZ ( − n − ( χ ( b , s )) n n ! (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (C8)using f diff ,nj = q γ j − α j γ j [(1 − α j + γ j ) n − (1 − α j − γ j ) n ] . (C9)Finally the expression for double diffraction dissociation is given by σ DD Y Z = Z d b (cid:12)(cid:12) h Y ⋆ Z ⋆ | (cid:0) − e − ˆ χ ( b ,s ) (cid:1) | Y Z i (cid:12)(cid:12) = Z d b (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∞ X n =1 f diff ,nY f diff ,nZ ( − n − ( χ ( b , s )) n n ! (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (C10)Note that after carrying out the sum over n in Eqs. (C4), (C7), (C8), (C10) the cross sectionscan be written as impact parameter integrals over a sum of exponentials.The parameter range for α j and β j is limited by the unitarity constraint that all eikonalfunctions have to be non-negative1 − α j − γ j ≥ α j < / β j > . (C11)A good description of the data is found for α = 0 . β = 0 . N s soft and N h hard interactions follow from σ N s ,N h = Z d b X k =1 Λ k [2 λ k χ soft ( s, b )] N s N s ! [2 λ k χ hard ( s, b )] N h N h ! exp {− λ k ( χ soft ( s, b ) + χ hard ( s, b )) } , (C12)33ith Λ = (cid:18) − α Y γ Y (cid:19) (cid:18) − α Z γ Z (cid:19) Λ = (cid:18) − α Y γ Y (cid:19) (cid:18) α Z γ Z (cid:19) Λ = (cid:18) α Y γ Y (cid:19) (cid:18) − α Z γ Z (cid:19) Λ = (cid:18) α Y γ Y (cid:19) (cid:18) α Z γ Z (cid:19) (C13)and λ = (1 − α Y − γ Y )(1 − α Z − γ Z ) λ = (1 − α Y − γ Y )(1 − α Z + γ Z ) λ = (1 − α Y + γ Y )(1 − α Z − γ Z ) λ = (1 − α Y + γ Y )(1 − α Z + γ Z ) . (C14)For high-mass diffraction, it is assumed that a constant fraction of each cut (soft or hardinteraction) corresponds to an rapidity-gap final state. The corresponding cross section iswritten, see Eq. (C12) σ SDhm = δ ( σ , + σ , ) , (C15a) σ DDhm = δ ( σ , + σ , + σ , ) + β Y σ SD ,Z lm + β Z σ SD ,Y lm . (C15b)The factor δ is estimated by comparing with HERA data: 10% of all deep inelastic scatteringevents at low x correspond to diffraction ( δ ≈ ..