ISVHECRI 2016 Multiple production up to √ s = 13 TeV with the generator GHOST adapted to cosmic ray simulation

The model HDPM of CORSIKA has been updated and developed on the base of the recent measurements by LHCb, CMS, TOTEM... The new model GHOST involving a 4-source production reproduces correctly the pseudorapidity distributions of charged secondaries and has been testified with the data in the mid and forward rapidity region, especially in the complex case of TOTEM and also with the recent measurements of CMS up to √ s = 13 TeV (9.2 · 1016 eV). Special calculations have been devoted to the semi-inclusive data playing an important role in the cosmic ray simulation (fluctuations in the earliest collisions, individual cascades and remarkable events). Taking into account the violation of KNO scaling, the negative binomial distribution has been used pointing out an asymptotic behaviour of total charged multiplicities at primary energies exceeding 40 TeV (7 · 1017 eV). The interpretation of the most recent measurements suggests for EAS generated by primary protons a larger production of muons and a cascade maximum at higher altitude than previously assumed.


Introduction
In the sections that follow, we shall discuss the new improvements in cosmic ray simulation suggested by the most recent results of the LHC, especially the data of ATLAS, CMS, TOTEM, LHCb... up to √ s = 13 TeV used as the largest energy now available to elaborate and testify the new Monte Carlo generator suitable for cosmic ray simulations.
Section 2 is devoted to the consequences of the violation of KNO scaling on the semi-inclusive data after considering the limits of the earliest scaling behaviour only valid now in small energy bands.
Section 3 is focused on the basic features of very high energy hadron interactions and guidelines useful for extrapolations of EAS simulations up to primay energies of 10 20 eV.The 4-component model GHOST (Generator Hadronic Oversteering Secondary Treatment) reproduces general characteristics of p-p collisions up to √ s = 13 TeV.
Section 4 evocates some measurements not explained in cosmic ray data above √ s = 1 TeV, especially in high energy spectra of secondaries measured in EAS, spikes and coplanar emission in Gamma-ray families.
Observing that Extensive Air Showers (EAS) are random individual events, we emphasize here the interest of the semi inclusive data delivered by the LHC at ultra high energy.

Consequences of the violation of KNO scaling and semi inclusive data
Each new step in energy reached by the colliders is an opportunity to evaluate the predictions of the interaction models.The simulation of the hadronic cascade in cosmic ray physics needs at least the parameters describing the different cross sections, the total multiplicity of secondaries, their energy and transverse momenta distributions, the energy dependence of those parameters and the statistic laws governing their fluctuations.The violation of Koba, Nielsen, Olesen (KNO) scaling observed in UA5 [1], for instance, introduced the negative binomial distribution well adapted to the prediction of multiplicities of secondaries in cosmic ray cascades.

KNO scaling in central region
Several attempts have been performed to introduce other scaling based on the Feynmann's scaling [2] describing the rapidity distribution as a central plateau with gaussian wings; the link with the pseudorapidity distribution was treating the central dip with dynamical considerations, whereas the extension from ISR energies up to SPS energies introduced both a rising plateau and a width increase derived from UA5 data on pseudorapidity distributions.The KNO scaling (concerning all the distribution) was rejected in 1983 following the contradiction between ISR and UA5 data; taking into account the wider distributions of multiplicity than KNO predicted also by the Dual Parton Model (DPM) [3] and the "fireball model" [1].It was decided to explore the multiplicity distributions in limited pseudorapidity intervals.A tolerable scaling agreement was obtained in central region for η ≤ 1.3 or η ≤ 1.5 in the case of ISR, UA5 and UA1 for instance.However, several uncertainties remain.We have simulated with GHOST such situation for comparison between LHC data ( √ s = 7 TeV) and UA5 ( √ s = 540 GeV).Important discrepancies (Figure 1) of the distributions are distinguished between ATLAS and UA5 for η ≤ 2.5 and η ≤ 3. as well as in the case of ALICE and UA5 respectively for η ≤ 1. and η ≤ 1.5 [4] [5].Since LHC measurements of the average multiplicities n ch are not yet very accurate, a strict comparison versus the KNO variable z = n ch / n ch cannot be shown, but there remains poor chance for this scaling in central region.Our histograms calculated with GHOST confirm the large discrepancies with UA5 data when the energy is rising up to √ s = 7 TeV.Similarly, the attempts made in the fragmentation region of the inclusive pseudorapidity distributions plotted in the beam frame were not conclusive; the measurements were obtained for ISR up to η − y beam = −0.75 for 53 GeV in ISR, but only up to η − y beam = −2.4 for 900 GeV in UA5 [6].The genuine fragmention region cannot be reached because the range of forward directions at |x f | ≥ 0.05 is experimentally unavialable.

Empirical scaling and semi-inclusive data
Fortunately an "empirical scaling" was pointed out in the central region by UA5 as a relation between the normalized central density and the normalized multiplicity [6].The ratio ρ(0)/ ρ(0) = f (z) (instead of f (z, s) is here independent of √ s), ρ(0) being the central density for a given charged multiplicity n, ρ(0) being the average central density (for a more simple legibility n, n ch , N and N ch are used interchangeable).The KNO variable z = n/ n is the normalized multiplicity.The general dependence on z is shown in the Fig. 2. Following our adjustment, we propose for f (z) the numerical representation: This kind of scaling is not yet confirmed in LHC data as the semi-inclusive measurements are performed only in small windows of pseudorapidity, but it remains present in models extrapolations at larger energies.Following UA5 and our considerations in CORSIKA user's guide [7], the fluctuations of < N ch > for a fixed primary energy can be reproduced at very high energy by a more simple relation corresponding to the negative binomial distribution (NBD) at very high energy; this simple function [8] represents in terms of KNO reduced variablesz = n n the fluctuations as: The values of the parameters correspond here to n = < N ch > for the charged NSD component, P n being the probability of an incoming multiplicity N ch .The evolution of those fluctuations of N ch is plotted in the Fig. 3 at ISR and UA5 energies, completed by the curve and the histogram illustrating the calculation with GHOST at √ s = 8 TeV.The asymptotic form at very high energy (40 TeV) in the unfulfilled project of the Superconducting Super Collider (SSC 1986) is the most enlarged variation with the simple expression (z) = 4z e −2z .The probability in p-p collision to receive randomly very small N ch or very large N ch increases at very high energy as shown in the Fig. 3.
Such contrast changes the statistical treatment of the earliest interactions in hadronic cascades at UHE and we underline the neighbourhood of the histogram at √ s = 8 TeV with the asymptotic expression of the SSC.Thanks to the function ζ the semi-inclusive data can be governed by the integro-differential system: ) m r is the ratio of central mean rapidity density and mean central pseudorapidity density derived from the "dip" existing in the centre of the pseudorapidity distribution, resulting from the mass m and the transverse mass m T of the secondaries as m r = 1 − m 2 /m 2 T .In the case of the 4 gaussian generation (one pair of functions in each hemisphere), symmetrics around the center of mass, a i (e −0.5u i + e −0.5v i ) ( 7) it is possible to use the opportunity of the scaling of function ζ in the relation between the center y i and the width σ i of each gaussian function as After introducing one proportion χ of the multiplicity distributed to the pair of gaussian centered in central region and in mid-rapidity region, it is possible to obtain with a minimal Monte Carlo generation of random deviates the total original rapidity distribution source of the distribution of pseudorapidity shown in the Fig. 4 for various range of multiplicity (respectively χ = 20% and 80% or for the total NSD multiplicity in the Fig. 7).The couples rapiditytransverse momentum (y, p t ) are generated using the same p t distribution as in CORSIKA [7].The total average inelasticity of the NSD collisions calculated with GHOST is 0.71 indicating a larger proportion of the primary energy distributed to the secondaries.The generator GHOST has been used here to appreciate the elements provided by semi-inclusive measurements which will be distributed soon by the LHC.The NSD behaviour is shown here in the Fig. 4, corresponding to the histogram of the Fig. 3 for √ s = 8 TeV.

Theoretical models and phenomenological collision generators
High energy cosmic ray interactions are mainly understood with the "soft physics" data provided by the colliders.Therefore, the approach of microscopic models starts from the parton distribution functions, distinguished for valence quarks and diquarks, sea quarks and gluons.In the case of p-p collisions, there are common features between models used in accelerators or in cosmic rays simulations, respectively models such as Dual Parton Model, Phojet, Pythia on one side, models like QGSJet, Sybill on the other side [9].The geometrical and statistical models such as the cluster model GENCL of UA5 [10] are more connected with the phenomenological features of the experimental data and built on analogy with quantum optics.Similarly, the Hybrid Dual Parton Model (HDPM) was inspired by the most simple double pomeron structure of DPM.It represents the pair of quark-diquark chains stretched between the valence quark and diquark of the projectile and the target by one pair of gaussian functions with a correct overlap to reproduce the expected rapidity distribution with the profile of one plateau with gaussianshaped wings [11][12][13].
Unfortunately, in the interpretation of cosmic ray data, the approach has a limit due to the difficulties of measurement in the very forward region of rapidity and the determination of the properties of the leading cluster.In other words, the role of the pilot particle fundamental in hadronic cosmic ray cascades remains dependent on adequate assumptions.
All the models have been extended to p-A and A-A collisions in order to describe the collision of cosmic ray hadrons with air constituents and how the comparison with the LHC data becomes also possible.

Guidelines and extrapolations at ultra high energy
Conversely, there exist also models elaborated directly for collisions of heavy nuclei which can be tuned to describe p-p collisions; this is the case, for instance, of the non-equilibrium statistical relativistic diffusion model (RDM) based on a Fokker-Planck type transport equation [14].Hence, the earliest and most simple data produced by colliders and accelerators is the charged Non Single Diffraction (NSD) distribution.The characteristics of this distribution: average central rapidity density ρ(0), total average multiplicity of charged secondaries n ch (as far as it can be approached by the integration of the densities collected till the smallest angles of emission), the maximum and the full width at half maximum (FWHM), transverse momentum distribution, are the normal tests of the preliminary validity of the models recognized.
A next step of quality is reached with additive criteria adapted to the same pseudorapidity distribution data considered in small bands of multiplicity i.e. in a semiinclusive analysis as it was shown in the Fig. 4.

Energy dependence of total charged multiplicity and central pseudorapidity densities at LHC
Up to now the average central pseudorapidity density has been measured up to √ s = 13 TeV (only for the inelastic component as measured recently by CMS), the other components being limited at √ s = 8 TeV.For instance the widening of the full pseudorapidity distribution at half maximum (FWHM) can be used as an estimator up to √ s = 8 TeV thanks to TOTEM data giving a part of the pseudorapidity densities near the forward region, but this is not yet the case at √ s = 13 TeV as shown in the Fig. 5.We can observe also in the Fig. 5 the different behaviour of the full width at half maximum (FWHM), σ H W , on the one hand, and the rise of the available rapidity (or rapidity Y 0 of the center of mass in the laboratory system) on the other hand.Both dependencies on s reflect a logarithmic increase σ H W = 0.3182 ln(s) + 0.5 ( 9 ) Y 0 = 0.502 ln(s) + 0.3 (10) We begin however to have new elements in the interval √ s = [1.− 13.] TeV useful for a lever arms of two guidelines in cosmic ray simulation, the average multiplicity dependence and average central pseudorapidity dependence to describe in better condition the cosmic ray data up to about 10 17 eV.New extrapolations risked up to 10 20 eV could be evaluated.The energy dependence in the Fig. 6 of the total mean charged multiplicity (NSD component) is shown hereafter adding to the S p pS collider measurements [1] the data of CMS and TOTEM at √ s = 8 TeV.An important enhancement is thus appearing when we take into account the multiplicity estimated using GHOST after integration of the density of both rapidity and pseudorapidity in the best agreement with the pseudorapidity measured at √ s = 8 TeV by CMS and TOTEM [15].We note a large discrepancy of the asymptotic behaviour (solid line) with the previous tendency (dashed line) expected from UA5 data [1].The first estimation of the total average multiplicity < N ch−U A5 > by a power law has to be updated as we suggest above √ s = 1 TeV introducing another power law for < N ch > as follows: There is also an important rise of multiplicity in the central region characterized by the average central pseudorapidity density ρ(0) with a power law dependence on s such as 0.708s 0.118 .However we use hereafter the parabolic logarithmic fit of CMS [16], the different behaviour with the logarithmic dependence ρ1 (0) of UA5 replaced by ρ(0) as follows is shown in the Fig. 6b ρ(0) = 3.17 − 0.372ln(s) + 0.0291ln(s 2 ) ( 13) ρ1 (0) = 0.24 ln(s) + 0.1 (14) The dependence on energy evaluated in LHC up to √ s = 13 TeV exhibits multiplicity enhancement in central rapidity region as well as in mid-rapidity region and a total charged mean multiplicity close to 100 particles can be expected at √ s = 15 TeV.

Average rapidity and pseudorapidity distributions up to √ s = 13 TeV
At present the charged density measurements at √ s = 8 TeV cover central and mid-rapidity region limited to 7 units of pseudorapidity.A comparison of a profile generated with the double gaussian generator GHOST is shown in the Fig. 7 and compared with the NSD measurements of CMS and TOTEM [15].This result takes into account for both rapidity and pseudorapidity the generation of individual collisions containing for instance the semi-inclusive data shown in the Fig. 4. The simple integration of the simulated distribution in agreement with the data on average densities of pseudorapidity < N ch > = TeV.An attempt with GHOST generator (dotted histogram) suggests one possible pseudorapidity distribution among several tendencies; it is only one example taking account the enhancement of total multiplicity suggested by the Fig. 6, a measurement of central NSD pseudorapidity being necessary.

Cosmic ray data related with the LHC energy domain
A large part of cosmic ray remarkable events such as coplanar emission [17] or spikes [9] have been observed in different conditions (mountain altitude or stratosphere).They appear mainly between 1 PeV and 100 PeV in laboratory frame being covered now by the LHC data up to √ s = 13 TeV.Very small emulsion X-ray chambers of 40 × 50 cm have collected data in the stratosphere with balloons (i.e.JACEE, siberian flights above 30 km altitude), Concorde supersonic flights at 17 km altitude and even with japanese air liners near 10 km altitude).In parallel Emulsion-X-ray chambers or and X-ray chambers were also installed above 4 km altitude in Pamir, Tibet with very large area in order to collect events at energies exceeding 10 17 eV.
As expected more A-A collisions were collected with balloons, when more p-A collsions were collected during airplane flights with a statistics exhausted above 5 • 10 16 eV interesting with "spikes" in pseudorapidity distribution.This data contains events with local multiplicity fluctuations which were also observed at lower energy in UA5 data suggesting the examination of different hypothesis such as formation of quark-gluon plasma, but the Monte Carlo simulations of the CERN suggest also that the fluctuations occurred with a frequency similar to the data sample.New comparison of those statistics in LHC and in altitude flights appears necessary.
Similarly near and above √ s = 5 TeV, the unexpected alignment of secondaries has been observed in 2 independent stratospheric flights, in parallel with the collection of Pamir registrations.Additional simulations are now necessary to understand the coplanar emission in a new context; for instance diquark breaking of valence diquark could generate very large transverse momenta from the extreme string tension between quarks contained in each diquark.There could be also other investigations in relation with the ridges seen in p-p and p-A collisions at the LHC.

Conclusion
The 4 sources model GHOST is a very fast Monte Carlo collision generator able to reproduce both inelastic and NSD data up to √ s = 13 TeV.This circumstance allows a more refined analysis of cosmic ray data at energies ISVHECRI 2016 above the "knee" at 3 − 5 10 15 eV.The enhancement of multiplicity suggests in EAS simulation a larger abundance of muons and a maximum of development at higher altitude than with previous models.Altogether an abundance of heavy primaries lower than with previous models can be expected above the knee.Furthermore, the energy guide lines derived from the complete data of the LHC shall provide new elememts in the polemics on GZK cut-off.
The data collected at √ s = 7 − 8 TeV is confirmed and will be completed soon for the NSD, inelastic components at √ s = 13 TeV and hopefully for cosmic ray data in the leading cluster with LHCf.This could help the interpretation of cosmic ray anomalies in the "knee" region above 2 TeV.The high energy photons in EAS in Tian Shan and the gamma ray families above 5 TeV (near 10 PeV in laboratory frame) exhibit steeper integral energy spectra than at lower energies.Such signal could be in relation with an important part of energy lost in the leading cluster.

Figure 3 .
Figure 3. Fluctuations of total charged NSD multiplicity.The full line represents the asymptotic behaviour at √ s = 40 TeV, the histogram and the dashed line correspond both to the calculation of GHOST and the function (z, k) at √ s = 8 TeV, whereas the dotted and dashed curves correspond respectively to √ s =900 GeV and √ s = 53 GeV.

Figure 5 .
Figure 5. Energy dependence of the available rapidity Yo and of the FWHM, σ H W , respectively, open squares, solid line and dashed line from GHOST, dark squares from CMS-TOTEM, triangles from QGSJET01 model, horizontal line from the relativistic diffusion model RDM [14].

Figure 6 .
Figure 6. Figure 6a (upper) Total charged multiplicity (NSD) dependence on s for E 0 = 10 13 −10 17 eV following UA5 (solid line) and the asymptotic tendency using GHOST (dashed line).3 lower points are from UA5, the upper point at 13 TeV is estimated with GHOST. Figure 6b (bottom) Central pseudorapidity density dependence on s for E 0 = 10 13 −10 17 eV (solid line for NSD component with one point from GHOST at √ s = 13 TeV, dashed line with diamonds devoted to the inelastic component is aproximately 0.8695• ρ(0) NSD).

Figure 7 .
Figure 7. Average pseudorapidity density simulated with GHOST at √ s = 8 TeV compared with CMS and TOTEM measurements taking into account respective limits of Pt.

8
TeV an average total multiplicity < N ch > = 80 ± 2. The results obtained by GHOST are also compared with the recent data concerning the inelastic pseudorapidity respectively at √ s = 8 TeV and √ s = 13 TeV [16] as shown in the Fig. 8.There are not measurements published up to now concerning the NSD component at √ s = 13

Figure 8 .
Figure 8.Average pseudorapidity density simulated with GHOST compared with CMS measurements for the inelastic data of CMS at √ s = 8 TeV (lower histogram) and at √ s = 13 TeV (middle histogram).The upper histogram (dotted) concerns the NSD data and is just a preliminary attempt.
This work was thanks to the French-Polish Collaboration Agreement IN2P3-COPIN.Part of Z. Plebaniak work was supported by NCN grant UMO-2016/20/T/ST9/00589.