Simulation study on the effects of diffractive collisions on the prediction of the observables in ultra-high-energy cosmic ray experiments

The mass composition of ultra-high-energy cosmic rays is important for understanding their origin. Owing to our limited knowledge of the hadronic interaction, the interpretations of the mass composition from observations have several open problems, such as the inconsistent interpretations of $\langle X_{\mathrm{max}}\rangle $ and $\langle X_{\mathrm{max}}^{\mu}\rangle $ and the large difference between the predictions by the hadronic interaction models. Diffractive collision is one of the proposed sources of the uncertainty. In this paper, we discuss the effect of the detailed characteristics of diffractive collisions on the observables of ultra-high-energy cosmic-ray experiments by focusing on three detailed characteristics. These are the cross-sectional fractions of different collision types, diffractive-mass spectrum, and diffractive-mass-dependent particle productions. We demonstrated that the current level of the uncertainty in the cross-sectional fraction can affect 8.9$\mathrm{g/cm^2}$ of $\langle X_{\mathrm{max}}\rangle $ and 9.4$\mathrm{g/cm^2}$ of $\langle X_{\mathrm{max}}^{\mu}\rangle $, whereas the other details of the diffractive collisions exhibit relatively minor effects.


Introduction
The origin of ultra-high-energy cosmic rays (UHECRs) above 10 18 eV is one of the most important open questions in astrophysics. The mass composition of these cosmic rays is the key information to understand their origin. UHECRs are observed by measuring extended air showers and mass-sensitive observables, such as depth of the maximum shower development X max , depth of the maximum muon productions X µ max , and number of muons at the ground N µ extracted from air shower data. The mass composition is estimated by comparing these observables and their predictions based on Monte Carlo (MC) simulations using experiments, such as the Pierre Auger Observatory [1][2][3][4][5] and the Telescope Array Collaboration [6]. However, owing to our limited knowledge about the hadronic interactions at such a high energy, interpretations of the mass composition have several open problems. The interpretation from X µ max yields a heavier composition than that from X max [4], and c The Author(s) 2012. Published by Oxford University Press on behalf of the Physical Society of Japan. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by-nc/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.  productions, is not well understood yet. Experimentally, these collisions are characterized by a rapidity gap of the produced particles from dissociation of the incident particles [12]. Because the energy of the produced particles in the diffractive collisions is higher than in the other types of collisions, the diffractive collision is expected to be more relevant to the air shower observables [8]. There are three types of diffractive collisions: single diffractive (SD) collisions, double diffractive (DD) collisions, and central diffractive (CD) collisions. One of the initial particles dissociates, whereas the other is intact in the SD collisions. Both the initial particles dissociate in the DD collisions. The CD collisions are diffractive collisions with particle productions without dissociation of the initial particles. Figure 1 displays the Feynman diagrams of these collisions. In this study, the inelastic collisions excluding the diffractive collisions are called as non-diffractive (ND) collisions. The SD collisions with a projectile cosmic-ray dissociation and those with a target air nucleus dissociation induce different effects on an air shower, as displayed in Fig. 2. In this paper, the former collision type is called as the projectile single diffractive (pSD) collision, whereas the latter is called as the target single diffractive (tSD) collision. 3/26 The particle productions in the diffractive collisions are characterized by diffractive mass M X , which is an invariant mass of the particles in a dissociation system defined as, where p i and n are the four-momentum of the i-th particle and the number of particles in the dissociation system, respectively. The modeling of the diffractive collisions can be divided into two parts: the differential cross-section of the diffractive mass (diffractive-mass spectrum), dσ dMX , and the particle productions of the dissociation system as a function of M X . Therefore, there are three characteristics in the diffractive collisions: • cross-section of each collision type, • diffractive-mass spectrum, and • particle productions as a function of the diffractive mass.
Because the diffractive collisions are characterized by a low multiplicity and a high elasticity, they deepen the air shower developments [8]. The collision types in the diffractive collisions also affect the UHECR observables as follows. If the first interaction of an air shower is a pSD or a DD collision, the projectile particle dissociates and produces several particles. However, if the first interaction is a tSD or a CD collision, the projectile particle is intact; therefore, X max becomes one interaction length deeper. The modeling of the diffractive collisions is also important. The diffractive collisions with smaller M X are characterized using a smaller number and higher energy of produced particles, resulting larger X max . There are differences in the modeling of the particle productions as a function of M X between the models, and these differences affect the uncertainty of the energy and the number of produced particles in the diffractive collisions.
Recently, the cross-sections of the diffractive collisions were measured by TOTEM [13,14], ATLAS [15][16][17], CMS [18], and ALICE [19] at the Large Hadron Collider (LHC) in CERN. Since most of the particles are produced only in the very forward regions in the low diffractive-mass cases, measurements of the low diffractive-mass events with M X < 3.4 GeV are limited and there are large uncertainties in the cross-section. Since the total and elastic cross-sections are precisely measured using Roman Pots by TOTEM [20,21] and ATLAS [22], the uncertainty in the inelastic cross-section is very small. However, the relative cross-sections between the collision types are equally important as the absolute cross-sections. Moreover, the particle productions in the diffractive collisions are not well constrained by the collider experiments. Thus, these three quantities have large uncertainties.
There are large differences in the predictions of the diffractive collisions in the latest hadronic interaction models. Figure 3 displays the cross-sectional fraction in the proton-air collisions for 10 19 eV and 10 17 eV projectile protons. The difference in the cross-sectional fractions of the ND collisions in the hadronic interaction models is approximately 10%, and SIBYLL 2.3c [23,24] and EPOS-LHC [25] present the largest and smallest values, respectively. The difference in the cross-sectional fractions is relatively larger for the DD and CD collisions than that for the pSD and tSD collisions. Figure 4 displays the diffractivemass spectra of the proton-air collisions for 10 19 eV and 10 17 eV projectile protons. Here, log 10 (ξ) = log 10 (M 2 X /s), where √ s is the center-of-mass energy of the proton-air collision. Large differences can be noted between SIBYLL 2.3c and the other models, particularly 4/26  Approximately 40,000 events are simulated for each case using CONEX v6.40 [11].
in the lowest diffractive-mass regions. QGSJET II-04 [26] displays a strong peak in the low diffractive-mass region, whereas SIBYLL 2.3c does not presented diffractive-mass dependencies and EPOS-LHC has a bimodal spectrum. Even though the latest hadronic interaction models are updated using the experimental results from the LHC, they do not reproduce the results of the measurements of the diffractive collisions from the LHC [27].
In the following sections, to understand the effect of these quantities in predicting the observables, the effects of the cross-sectional fraction, diffractive-mass spectrum, and particle productions on the UHECR observables are discussed. The focus is on the large differences in the cross-sectional fractions of the models and in the diffractive-mass spectra of SIBYLL 2.3c and the other models.

5/26
3. Air shower simulations 3.1. Simulation set-up Air shower events were simulated using CONEX v6.40 [11] with three hadronic interaction models, EPOS-LHC, QGSJET II-04, and SIBYLL 2.3c. UrQMD [28,29] was used for the collisions with incident particle energy less than 80 GeV. Primary protons of 10 19 and 10 17 eV with incident zenith angle 60 • were assumed. The following three types of modifications of CONEX were applied: (1) the collision type at the first interaction was added to the output of each shower event. For the definition of the collision type, we used the collision information obtained from each hadronic interaction model, similar to in a previous study [9]. (2) To simulate the events used in Sec. 4.2.2 and 5.2, the fractions of various types of diffractive collisions were modified. (3) When the diffractive-mass dependence, as described in Sec. 4.3 and 5.3, was studied, the collision type at the first interaction was fixed to the pSD collision and the diffractive mass at the first interaction was calculated for each air shower from the momenta of the particles in the dissociation system. The particles from the dissociation system and an atmospheric nucleus were separated using a rapidity threshold, which here is set as 1.5 in the laboratory frame. The number of produced events was 40,000 for each case. Although the calculation is effective for most of the events, less than 0.5 % of the events have masses smaller than that of a proton, which is non-physical, owing to the miss-separation of the dissociation system in the calculations. These non-physical events were not used in the analysis.
CONEX is sufficiently fast for simulating a sufficient number of events because of employing the one-dimensional cascade equation. To confirm the predictions by CONEX, X max calculated by CONEX was compared to the result by a three-dimensional air shower simulation package, COSMOS 8.035 [30], for a 10 17 eV proton primary case. To reduce the time-consuming calculation, a one-dimensional (1D) analytical calculation is conducted for the calculations of electromagnetic cascade showers in COSMOS. In this comparison, QGSJET II-04 was used as a high-energy interaction model, and the zenith angle was 0 • . UrQMD was used for the low-energy interactions below 80 GeV in CONEX, whereas DPM-JET III [31] and PHITS [32] were used in COSMOS for the interactions in 2-80 GeV and less than 2 GeV, respectively. In COSMOS, the number of electrons was calculated every 25 g/cm 2 , and X max was calculated by fitting them to the Gaisser-Hillas function. The number of simulated air showers was 3,000 for each package. The result of X max was 683.4 ± 1.3 g/cm 2 for CONEX and 680.9 ± 1.3 g/cm 2 for COSMOS. The one-dimensional cascade equation in CONEX presents a reasonable agreement with the three-dimensional (3D) cascade simulation by COSMOS.

Longitudinal shower profiles and collision type
Before discussing about X max , X µ max , and N µ , the basic relationship between the diffractive collision type and the longitudinal air shower profile is discussed. To understand the relationship, the events are categorized by the collision type at the first interaction of the air shower, and the mean longitudinal shower profile is calculated for each category. The profiles of the electrons plus positrons (electrons hereafter unless mentioned) in 10 19 eV showers with a zenith angle of 60 • using EPOS-LHC are displayed in Fig. 5 (a). The peak position of the 6/26 diffractive collision categories is 40-50 g/cm 2 deeper than that of the ND collision category, whereas the differences between the diffractive collision categories are small.
The shower profiles of the muons and muon productions using EPOS-LHC are presented in Figures 5 (b) and (c), respectively. Similar to the electron profiles, the peak positions is 40-50 g/cm 2 deeper than that by the ND. Additionally, the pSD and DD collision categories exhibit a smaller peak height than the other categories. For depth X > 1100 g/cm 2 , the number of muons of the tSD and CD collision categories is the largest among the five categories, whereas those of the pSD and DD collision categories are the smallest. The difference in the number of muons of the categories depends on the depth, and it decreases in a deeper atmosphere. The difference between the average and the ND collision categories is small.
These properties are directly related to the observable variables in the air shower experiments. The peak depth of the electron profile defines X max . It is found that the diffractive collisions deepened X max , whereas the collision types in the diffractive collisions are insensitive to X max . The peak depth of the muon production profile defines X µ max , the diffractive collisions deepened X µ max , whereas the collision types in the diffractive collisions are insensitive to X µ max as in case of X max . The number of muons at the ground defines another observable, N µ . The differences in N µ in the diffractive collision categories depends on the depth, whereas the relationship between the categories are the same for depth X > 1100 g/cm 2 , N µ of the tSD and CD collisions categories is large, whereas that of the pSD and DD collisions categories is smaller. These features depend on the choice of the hadronic interaction model. The quantitative discussions of the model dependencies and the effects on the observables are discussed in Sec. 4 and Sec. 5.

4.
Effect of diffractive collisions on X max and X µ , respectively. Even in the same category, there are differences in the predictions of the models. SIBYLL 2.3c displays that the largest difference is for X max and X µ max of the pSD and tSD collision categories, which is 20-30 g/cm 2 . Moreover, the size of the difference is as large as the difference between the ND and pSD collisions categories. For EPOS-LHC and QGSJET II-04, the differences between the pSD, tSD, and DD collision categories are smaller than those between the pSD and ND collision categories. These features suggest that the cross-sectional fraction of the diffractive collisions to the ND collisions is more important than that of the pSD and tSD collision for the EPOS-LHC and QGSJET II-04 cases. In comparison, the cross-sectional fractions of the pSD and tSD collision categories may affect the SIBYLL 2.3c case. Details of the effects of the cross-sectional fractions are discussed in the following section.
Focusing on the values in each collision category, the difference between QGSJET II-04 and EPOS-LHC is found to be mostly independent of the collision category, and it is about 10 g/cm and 30 g/cm for X max and X µ max , respectively. The difference in SIBYLL 2.3c and the other models, however, vary with the collision categories. For example, the differences in X µ max of SIBYLL 2.3c and EPOS-LHC for the pSD and the DD collision categories are larger than these for the ND and the tSD collision categories. The projectile protons 7/26   (cosmic rays) dissociate only for the pSD and the DD collisions, therefore, the diffractive mass spectrum affect the shower development only in the pSD and the DD case. Thus, this may be related with a large difference of the diffractive mass spectrum between SIBYLL 2.3c and the other models. It is important to examine the effects of the modeling of the diffractive collisions, including the effect of the diffractive mass spectrum on X max and X µ max .

Effect of cross-sectional fractions of collision types
The effects of the cross-sectional fractions is two-fold. One is the effect of the first interaction in the air showers focusing on all the cross-sectional fractions. The other is the effect including the secondary interactions over the whole air shower, focusing on one cross-sectional fraction that displays the largest effect at the first interaction.

Effect of first interaction.
The effects of the cross-sectional fractions at the first interaction were studied using X max and X µ max , as introduced in Fig. 6. The expected value of X max with an arbitrary mixture of the collision types, X modified max , is defined as the weighted mean of X max of the i-th collision type X i max , where f i is the cross-sectional fraction of the i-th collision type. X µ,modified max is defined similarly as above. Using this reweighting method, the effects of the cross-sectional fractions can be studied by changing f i without repeating the time-consuming MC simulations. Instead of directly changing f i 's, the four ratios of the cross sections, R 1 , R 2 , R 3 , and R 4 , is introduced. As illustrated in Fig.7, the ratios are defined between 0 and 1. R 1 is the ratio of all the diffractive collisions to the inelastic collisions. R 2 is the ratio of the SD (pSD + tSD) Fig. 7 Schematic view of four ratios R 1 , R 2 , R 3 , and R 4 collisions to the sum of the SD and DD collisions. R 3 is the ratio of the tSD collisions to all the SD collisions. R 4 is the ratio of the CD collisions to all the diffractive collisions. The cross-sectional fractions, f i 's, can be calculated from these four ratios as, X modified max and X µ,modified max are calculated using the cross-sectional fractions and Eq. 2 by changing the ratios.
The results of changing the cross-sectional fractions using ratios R 1 , R 2 and R 3 are displayed in Figures 8 and 9. Upper panels display X modified max and X µ,modified max (lines) and the original predictions X original max and X µ,original max (black circles) by three interaction models for 10 19 eV proton primaries. The middle plots are the same as the upper plots except for the 10 17 eV proton primaries. Magenta (solid and dash-dotted), blue (dashed and dash-twosotted), and green (dotted and dash-three-dotted) lines display the results of EPOS-LHC, QGSJET II-04, and SIBYLL 2.3c, respectively. The error bars of these plots display the statistical errors of the MC simulations. The bottom plots present the differences in X modified max and X µ,modif ied max from the original predictions, where ∆ X max is defined as In the bottom plots, the error bars only for 10 19 eV are drawn for visibility. The calculation for R 4 is not performed, because SIBYLL 2.3c does not include CD collisions. Based on the upper plot of Fig. 8 (a), the original R 1 ranges from 0.08 (SIBYLL 2.3c) to 0.18 (EPOS-LHC), which is regarded as the model uncertainty of R 1 . When R 1 is changed in this range, X max changes by 3.7 g/cm 2 for SIBYLL 2.3c and 5.0 g/cm 2 for QGSJET II-04, and X µ max changes by 3.4 g/cm 2 for SIBYLL 2.3c and 4.8 g/cm 2 for QGSJET II-04 for 10 19 eV. However, the shifts when R 2 and R 3 are changed within the model uncertainties are smaller than 0.4 g/cm 2 , which is comparable to the statistical errors. These results signify that the cross-section ratio of the diffractive collisions to the inelastic collisions R 1 has a dominant effects on X max and X µ max . In comparison, the cross-section ratios of the diffractive collision types i. e. R 2 and R 3 , are not sensitive to X max and X µ max .

10/26
These differences between the ratios originate from the differences in X max and X µ max of the categories seen in Fig. 6. The impact of the ratio of the diffractive and non-diffractive collisions is large, whereas the impact of the fractions of the different types of the diffractive collisions are small.
The ratio dependences of 10 17 eV for R 1 are similar to those of 10 19 eV, whereas the model uncertainties of the ratios of 10 17 eV (horizontal spread of the black circles) are smaller than those of 10 19 eV. Therefore, the effect of the model uncertainties of R 1 is smaller than that of 10 19 eV. When R 1 prediction is changed from SIBYLL 2.3c to EPOS-LHC, X max becomes 2.8 g/cm 2 larger for SIBYLL 2.3c and 3.8 g/cm 2 larger for QGSJET II-04, and X µ max becomes 3.0 g/cm 2 larger for SIBYLL 2.3c and 3.8 g/cm 2 larger for QGSJET II-04. When R 2 and R 3 are changed within the model uncertainties of the ratios, the shifts are smaller than 1 g/cm 2 .

4.2.2.
Effect of interactions over a whole air shower. The effects of the interactions over a whole air shower on X max and X µ max are discussed in this section, focusing on the crosssectional ratio of the diffractive collisions to inelastic collisions R 1 . It exhibits the largest effect at the first interaction, as presented in Sec.4.2.1. In this study, the ratio for collisions with projectile energy larger than 10 15 eV was modified. This threshold energy is defined by a previous study using CONEX [8]. The procedure of this modification is described below; (1) Reference R 1 ratios as a function of energy are estimated from the simulated ratios of EPOS-LHC and SIBYLL 2.3c, as displayed in Fig. 10 by fitting with a function defined as a log 10 E + b, where a and b are free parameters. The magenta (green) solid and the dashed lines are the fitting results of the ratio for proton-Nitrogen and π + -Nitrogen collisions simulated by EPOS-LHC (SIBYLL 2.3c), respectively. The fitting results of the proton-Nitrogen (π + -Nitrogen) collisions are used for the energy-dependent R 1 ratio for baryon-air nucleus (meson-air nucleus) collisions at a given energy in an air shower. Hereby, the modification with reference R 1 ratio taken from EPOS-LHC and SIBYLL 2.3c results is called as the EPOS-based R 1 case and the SIBYLL-based R 1 case, respectively. (2) For every collision with a projectile energy larger than 10 15 eV, a collision type is selected by following the cross-sectional fraction, f i , after the R 1 modification. (3) During the air shower simulation, a collision is generated according to each model, and the collision type of the produced event is compared to the one selected in (2). If the two types are not the same, another collision is generated repeatedly until those two types become the same.
Results of X max with this modification are summarized in Tab. 1. The difference between the value of EPOS-LHC and the one obtained with the modification using EPOS-based reference R 1 accounts for the systematic uncertainty of the method, and the size of the difference is 2.0 g/cm 2 . The difference for the SIBYLL 2.3c case is 0.9 g/cm 2 . This systematic uncertainty is due to the difference between the original prediction of R 1 and reference R 1 , which is defined by fitting. We discuss the impact of the change in R 1 from the original prediction to the reference one; however, there are some systematic uncertainties. The effect of the cross-sectional ratio is 8.9±0.4 g/cm 2 at the maximum for the EPOS-LHC-based case and 4.2±0.4 g/cm 2 at the maximum for the SIBYLL 2.3c-based case, and they are 1-2 11 Fig. 9 Ratio dependence of X µ max . Same as Fig. 8 times as large as the effects of the cross-section ratio at the first interaction, 3.7-5.0 g/cm 2 . This suggets that the effects of the ratios at the first and secondary interactions are equally important for the X max predictions. There is another important point. The differences in the X max predictions of the models with this modification are 31.9±0.5 g/cm 2 for the EPOS-LHC-based case and 32.0±0.4 g/cm 2 for the SIBYLL 2.3c-based case, as displayed in the bottom row of Tab. 1.   Fig. 10 Energy dependencies of the cross-section ratio, R 1 . The ratios of the proton-Nitrogen collisions (filled markers) and π + -Nitrogen collisions (open markers) for three hadronic interaction models using CRMC v1.7 [33] are displayed. Straight lines are the fitted results of the cross-section ratios with a log 10 E + b. The magenta solid (dashed) line is the fitted result for the proton (π + )-Nitrogen collisions by EPOS-LHC, and the green solid (dashed) line is the fitted result for proton (π + ) Nitrogen collisions by SIBYLL 2.3c. Moreover, they are larger than the difference without any modification: 27.4±0.4 g/cm 2 . This suggests that even if the cross-section ratio, R 1 , is fixed to the same value in three models, the difference in the X max predictions of the models become larger. In the original models, other sources of differences cancel out each other and reduce the difference in the X max predictions by chance. The cross-section ratios affect the mean value of X max predictions; however, they cannot solve the differences in the X max predictions of the models.
X µ max can be similarly discussed based on Tab. 2. The effects of the cross-section ratio, R 1 , are 9.4±0.4 g/cm 2 for the EPOS-LHC-based case and 4.4±0.4 g/cm 2 for the SIBYLL 2.3cbased case at maximum. The differences in the X µ max predictions between the models, with this modification, are smaller by 5.2±0.6 g/cm 2 and 0.4±0.6 g/cm 2 than those of the original predictions for the EPOS-LHC-based case and the SIBYLL 2.3c-based case, respectively. Moreover, they are less than 15 % of those of the original predictions. The cross-section ratios affect the mean value of the X µ max predictions; however, they cannot solve the differences in the X µ max predictions of the models.  Table 2 Results of X µ max with cross-sectional fraction modifications over the whole air shower. Reference R 1 (E) of the modification is set by fitting the ratio of the predictions by EPOS-LHC or SIBYLL 2.3c with energy above 10 15 eV.

Effects of diffractive collision modeling
As we discussed in Sec. 4.1, there are differences between the model predictions even for the same collision type at the first interaction, and these differences vary with the collision type. This suggests that the modeling of the diffractive collisions affect X max and X µ max . The modeling of the diffractive collisions can be divided into two parts; the diffractive-mass spectrum and the diffractive-mass dependent particle productions. Firstly, to understand the overview of the effect of the modeling, the diffractive-mass dependencies of X max and X µ max is discussed. Figures 11 (a) and 12 (a) display the diffractive-mass dependent X max and X µ max , respectively, for events whose collision type at the first interaction is fixed to the pSD collisions and for 10 19 eV incident protons. X max and X µ max become larger as the diffractive mass becomes smaller, and the diffractive-mass dependencies are small for log 10 (ξ) > −6. The dependencies of X max and X µ max are similar. The difference between the maximum and the minimum X max ( X µ max ) in the diffractive mass range is less than 30 g/cm 2 in each model. This suggests that even if we assume an extreme diffractive-mass spectrum, X max and X µ max originated by the pSD collisions are only affected by 30 g/cm 2 . Since projectile protons dissociate only for the pSD and the DD collisions, the diffractive-mass spectrum only affects the pSD and the DD collisions. Their cross-section are 5-12% in inelastic collisions, thus, the effect of the diffractive-mass spectrum at the first interaction is expected to be less than 3 g/cm 2 .
It is difficult to separate the effect of the particle production mechanisms, which is another important part of the diffractive collision modeling, and also the effect of the secondary interactions because the difference between models for each diffractive-mass bin displayed in Figures 11 (a) and 12 (a) reflects both the effects. However, the differences in the shape may indicate the differences in the particle production mechanisms, because it is hard to explain the change of the shape of the diffractive-mass dependencies only using the effect of the secondary interactions; if we assume the same particle production mechanism between SIBYLL 2.3c and EPOS-LHC, the secondary interactions should make a few g/cm 2 differences at log 10 (ξ) = −8, whereas more than 15 g/cm 2 differences at log 10 (ξ) = −2. The shape of diffractive-mass dependencies of EPOS-LHC and QGSJET II-04 are similar, whereas a difference can be seen between SIBYLL 2.3c and the other models; SIBYLL 2.3c only displays a small dip at log 10 (ξ) = −8. These features suggest that the modeling of particle production mechanisms in SIBYLL 2.3c is different between the low and the middle diffractive-mass regions, and that may affect X max and X µ max . These diffractive-mass dependencies are 15 Fig. 12 Diffractive-mass dependencies of X µ max for events whose collision type at the first interaction is the pSD collisions for (a) 10 19 eV and (b) 10 17 eV. same for 10 17 eV proton primaries as displayed in Figures 11 (b) and 12 (b). It is interesting that the models exhibit a good agreement in X max predictions only around log 10 (ξ) = −6 as found in Fig. 11 (b). The reason of such differences in diffractive-mass dependencies and energy dependence of the differences between models is worth studying in future.
To estimate the effect of the diffractive-mass spectrum quantitatively, the diffractive-mass spectrum at the first interaction was artificially changed. X max ( X µ max ) after modification of the diffrative mass spectrum is defined as, where P i MX is a probability of the i-th diffractive-mass bin of an arbitrary model in Fig. 4, and X (µ)i max is X (µ) max taken from the i-th diffractive-mass bin of Fig. 11 (Fig. 12). The effects of replacing P i MX between the three models on X max and X µ max , E Xmax , are defined as follows, where X The results of X max are displayed in Tab. 3 for 10 19 eV and in Tab. 4 for 10 17 eV. When the diffractive-mass spectrum of QGSJET II-04 is replaced by a spectrum of another model, X max becomes smaller. The size of the effects E Xmax is (4.1±0.7) to (5.1±0.7) g/cm 2 for 10 19 eV and (4.3±0.7) to (7.1±1.2) g/cm 2 for 10 17 eV. Differences between the model predictions when the same diffractive-mass spectrum was assumed at the first interaction are summarized in the bottom row of Tab. 3 and Tab. 4. The size of the differences is (18.6±0.7) to (21.0±1.2) g/cm 2 for 10 19 eV, which is larger by 4.5-6.9 g/cm 2 than the original difference. Therefore, the differences between the model predictions become larger even when the same diffractive-mass spectrum is used in the first pSD collisions. It is noted that these differences include the effect of the particle productions and the secondary interactions, and these two effects cannot be separated. For 10 17 eV, the size of the differences in predictions with the same diffractive-mass spectrum is (9.0±0.7) to (15.3±1.2) g/cm 2 , which is larger by 3.8-10.1 g/cm 2 than the original difference.
The results of X µ max are displayed in Tab. 5 for 10 19 eV and in Tab. 6 for 10 17 eV. The tendencies of the effect of the diffractive-mass spectrum are similar as the X max case. When the diffractive-mass spectrum of QGSJET II-04 is replaced by that in another model, X µ max becomes smaller. E Xmax is (2.8±0.7) to (4.8±0.7) g/cm 2 for 10 19 eV and (4.4±0.8) to (7.1±0.8) g/cm 2 for 10 17 eV. The sizes of differences between the model predictions by assuming the same diffractive-mass spectrum is larger by 3.5-4.9 g/cm 2 and by 4.8-7.2 g/cm 2 than the original differences for 10 19 eV and 10 17 eV, respectively. The differences between the model predictions become larger even when the same diffractive-mass spectrum is used in the first pSD collisions.
It must be noted that these estimations are performed by fixing the first interactions. The diffractive-mass spectrum only affects the pSD and the DD collisions. Because the crosssection of the pSD and the DD collisions are 5 to 12 % of the inelastic collisions, the effect of the diffractive collision modeling at the first interactions is 10 times smaller than the estimations above. Therefore, the total effects of the diffractive-mass spectrum on X max and X µ max are expected to be 0.5 g/cm 2 at maximum. These effects of diffractive-mass spectrum are much smaller than the effect of the cross-sectional fractions. The effects of the diffractive-mass spectrum are similar in X max and X µ max . From the results in Sec. 4.2 and 4.3, the effects of the diffractive collisions on X max are similar to that on X µ max . Therefore, the effects of the cross-sectional fraction and the diffractive mass spectrum do not solve the problem of inconsistent interpretations of the cosmic-ray mass composition from X max and X µ max . It must be noted that the collisions of low energy projectiles below 10 15 eV also affect X max and X µ max . According to Ref. [34], collisions with lower energy projectiles affect more on X µ max than X max .

Collision type and N µ
As discussed in Sec. 3.2, the number of muons of the tSD and the CD collisions categories is larger than that of the ND collisions and that of the pSD and the DD collisions is smaller than the ND collisions for the depth X > 1100 g/cm 2 . To understand the overview of the effects of the diffractive collisions, N µ is calculated for the simulated events by categorizing a collision type at the first interaction. Hereby, the zenith angle and the muon counting altitude 17/26 Table 3 The effect of the diffractive-mass spectrum on X max for the 10 19 eV proton primary case. E Xmax is the size of effects estimated by three diffractive-mass spectra as defined in Equation 10.  Table 4 The effect of the diffractive-mass spectrum on X max for the 10 17 eV proton primary case.  Table 5 The effect of the diffractive-mass spectrum on X µ max for the 10 19 eV proton primary case.  Table 6 The effect of diffractive-mass spectrum on X µ max for the 10 17 eV proton primary case. are fixed at 60 • and 0 m above the sea level, respectively, following the previous works that discussed N µ for inclined air showers [5,9]. The results of N µ are displayed in Fig. 13 for three hadronic interaction models. N µ of the tSD and the CD collision categories is approximately 5 % larger than that of the ND collision category, and the sizes of differences are the same level for three interaction models. N µ of the pSD and the DD collision categories is smaller than that of the ND collision category, however, the size of the differences depend on the hadronic interaction model; SIBYLL 2.3c displays 5 % differences between the pSD and the ND collision categories, whereas QGSJET II-04 exhibits a very small difference. Moreover, even in the same category, there are differences in predictions between the models. These features suggest both the cross-sectional fractions and the diffractive collision modeling affect N µ , those are discussed in Sec. 5.2 and 5.3, respectively. The same method used in Sec. 4 is followed, and the details of diffractive collisions affect N µ and the problem of muon excess is focused.

Effects of the cross-sectional fractions of the collision types
The effect of the first interaction by changing the cross-sectional fractions is estimated as in Sec. 4.2. The estimation of N µ for each type of the first interaction is displayed in Fig. 13. The fractions after changing one of the ratios R 1 , R 2 , and R 3 are calculated using Eq. 3, and an arbitrary mixture of N µ with modified fractions, N modified µ , is defined as, where N i µ is the average number of muons of the i-th type of the first interaction. The results are displayed in Fig. 14. The upper and middle panels display N modified µ (lines) and the original predictions N original µ (black circles) for 10 19 eV and 10 17 eV proton primaries, respectively. Magenta (solid and dash-dotted), blue (dashed and dash-two-dotted), and green (dotted and dash-three-dotted) lines display the results of EPOS-LHC, QGSJET II-04, and SIBYLL 2.3c, respectively. The error bars of these plots display the statistical errors of the MC simulations. The bottom plots exhibit the percentile differences, ∆N µ , defined as For 10 19 eV proton primaries, when R 1 is changed within the model uncertainty, N µ change 0.09 %, 0.08 % and 0.16 % for EPOS-LHC, QGSJET II-04 and SIBYLL 2.3c, respectively, whereas statistical errors are ±0.1 %. The effect of the model uncertainty in the R 2 case is 0.25 % at maximum, where the larger R 2 makes the N µ larger. The effect of the model uncertainty in the R 3 case is 0.10 % at maximum, where a larger R 3 makes N µ larger. N µ of the tSD and the CD collision categories is approximately 5 % larger than that of the ND collision category as seen in Fig. 13. However, owing to the small cross-sectional fraction of these collisions, the total effects on N µ are small. For 10 17 eV proton primaries, the ratio 19/26 dependence of N µ are similar to the 10 19 eV cases, whereas the size of effects is slightly larger as displayed in the bottom plots in Fig. 14. The largest effects of the R 1 , R 2 , and R 3 cases are 0.11 % in SIBYLL 2.3c, 0.30 % in EPOS-LHC, and 0.18 % in SIBYLL 2.3c, respectively, whereas the statistical errors are ±0.1 %. It can be conclude that the effect of the cross-sectional fractions at the first interaction is small. The modification of the cross-section ratio R 1 for collisions of projectiles above 10 15 eV over the whole air shower is performed as in Sec. 4.2.2, and the results are provided in Tab. 7. The result of EPOS-LHC (SIBYLL 2.3c) with modification to the EPOS (SIBYLL) based reference R 1 (E) exhibits systematic uncertainty of this method, amounting to 0.46±0.11 % (0.09±0.12 %). The size of the effect of this modification is 0.83±0.12 % (from (2.162±0.002)×10 7 to (2.144±0.002)×10 7 ) for SIBYLL 2.3c with modification using the EPOS based R 1 (E) and 0.15±0.11 % (from (1.999±0.002)×10 7 to (2.002±0.002)×10 7 ) for EPOS-LHC with modification using the SIBYLL based R 1 (E). Because the effect of the first interaction is 0.16 % for SIBYLL 2.3c when R 1 is changed to the prediction by EPOS-LHC, and 0.09 % for EPOS-LHC when R 1 is changed to the prediction by SIBYLL 2.3c, the modifications of secondary collisions above 10 15 eV make approximately 1.5-5 times larger effect than that of the first interaction. However, they are still smaller than 1 %. The cross-sectional fractions have a small effect on N µ even after considering the secondary interactions.
Muons are produced by the decay of pions and kaons in hadronic cascade showers and the effect of the diffractive collisions continue over the shower. From the discussion above, the secondary collisions of projectiles above 10 15 eV can enhance the effect 1.5-5 times compared to the first interaction. A previous study of 10 15 eV primary showers also suggested this effect is typically 4-7 [35]. Combining these studies through a 10 19 eV primary shower, enhancement of more than 10 due to the secondary interactions is possible. If 0.25 % effect of the crosssection ratio R 2 at the first interaction is enhanced by 10 times, the total effect on N µ is expected to be 2.5 %. This effect is still a minor source of the problem of muon excess, where 30 % discrepancies on N µ are reported [5,7]. Table 7 The results of N µ with cross-section ratio R 1 modifications for collisions with projectiles above 10 15 eV in a whole air shower. The reference R 1 (E) of modification is set by fitting the ratio of predictions by EPOSLHC or SIBYLL 2.3c with energy above 10 15 eV.   effects of modeling, the diffractive-mass dependencies of N µ are calculated as displayed in Fig. 15 (a) and 15 (b) for 10 19 eV and 10 17 eV, respectively. The simulations are performed by fixing the first interaction to the pSD collisions. A large difference can be seen in the low diffractive-mass regions around log 10 (ξ) = −9 (log 10 (ξ) = −7) for 10 19 eV (10 17 eV). SIBYLL 2.3c displays a dip structure in the low diffractive-mass region, whereas EPOS-LHC and QGSJET II-04 do not exhibit any such structure. The difference between the maximum and the minimum values in the diffractive-mass dependencies is approximately 5 % for SIBYLL 2.3c. Therefore, even if we assume an extreme diffractive-mass spectrum at the first interaction, N µ of the pSD collision category changes by 5 %. Since the projectile protons dissociate only for the pSD and the DD collisions and the diffractive-mass spectrum affects only these collisions, and the cross-sectional fraction of the pSD and the DD collisions in all inelastic collisions is 5-12 %. Thus, the effects at the first interaction is expected to be 0.6 % at maximum.
It is difficult to separate the effects of the particle productions and the secondary interactions in this analysis. However, a large difference in the shape of diffractive-mass dependencies between the models suggest the effects of particle production mechanisms, because it is hard to explain these differences in the shape only with the secondary interactions. Only SIBYLL 2.3c exhibits a dip structure in the low diffractive-mass region and small differences can be seen between EPOS-LHC and QGSJET II-04, and that may affect N µ .
The effect of the diffractive-mass spectrum is estimated by artificially changing the spectrum. N µ after modification of the diffractive-mass spectrum is defined as, where P i MX is a probability of the i-th diffractive-mass bin of an arbitrary model in Fig. 4, and N i µ is N µ taken from the i-th diffractive-mass bin of Fig. 15. The results for 10 19 eV are listed in Tab. 8. The largest effect of changing the diffractive-mass spectrum is 1.7±0.2 % (from (2.078±0.002)×10 7 to (2.042±0.004)×10 7 ) when the diffractive-mass spectrum of SIBYLL 2.3c is replaced with the spectrum of QGSJET II-04. The effect of replacing P i MX between the three models, E Nµ , is defined as follows, where N modified,max. when the diffractive-mass spectrum of QGSJET II-04 is assumed. This is caused by a combination of the diffractive-mass dependence of SIBYLL 2.3c and the diffractive-mass spectrum of QGSJET II-04; owing to the dip structure of SIBYLL 2.3c found in Fig. 15 and the low diffractive-mass peak of QGSJET II-04 found in Fig. 4, N µ is reduced when SIBYLL 2.3c meets the diffractive-mass spectrum of QGSJET II-04, whereas N µ is not changed so much for the other cases.
In this calculation, we focus mainly on the air showers with fixing the first interaction to the pSD collisions. The diffractive-mass spectrum only affects the pSD and the DD collisions. Since the cross-section of the pSD and the DD collisions are only 5-12 % of all inelastic collisions, the effects of the first interaction is 10 times smaller. Thus, the effect of the diffractive-mass spectrum at the first interaction is 0.1 % of N µ at maximum and as large as the effect of cross-sectional fractions. As discussed in the previous section, the effect of diffractive collisions can be enhanced by more than a factor 10 due to the secondary interactions. If 0.1 % effect of the diffractive-mass spectrum at the first interaction is enhanced by 10 times, the total effects on N µ is expected to be less than a few % and the effect is one of the minor sources of the problem of muon excess. The effect of the modeling of the diffractive collisions over the whole air shower is already discussed for SIBYLL 2.3c in ref. [10]. By changing the parameters of the diffractive collisions that affect the diffractive-mass spectrum of SIBYLL 2.3c, they concluded that the number 23/26 Table 9 The effect of diffractive-mass spectrum on the number of muon at the ground for the 10 17 eV proton primary case. of muons becomes 5 % smaller than the original SIBYLL 2.3c predictions. The diffractivemass spectrum used in this previous study is different from the current model predictions and the parameters may affect other parts of the hadronic interaction model, i. e. inelastic cross-section. Therefore, we cannot directly compare this previous work with the calculations in this section. However, the shape of the diffractive-mass dependencies in Fig. 15 suggests that if the diffractive-mass spectrum over the whole air shower is changed, whereas the other parts are simulated with EPOS-LHC and QGSJET II-04, the effect of the diffractive-mass spectrum will be much smaller than that in the previous work, because these two models show very small diffractive-mass dependencies.

Discussion and Conclusion
In this work, the effects of the diffractive collisions on X max , X µ max and N µ were discussed with focus on three detail characteristics of the diffractive collisions, (1) cross-sectional fractions among collision types, (2) diffractive-mass spectra, and (3) particle productions. The diffractive collisions make X max and X µ max deeper, and the cross-sectional fraction of the diffractive collisions in the inelastic collisions displays the largest effect among the detail characteristics. If we assume the current difference between the model predictions as an uncertainty, the maximum effect over a whole air shower is 8.9±0.4 g/cm 2 and 9.4±0.4 g/cm 2 on X max and X µ max , respectively. If the same cross-sectional fraction of the diffractive collisions is used for collisions of > 10 15 eV projectiles, the differences between the model predictions become larger by approximately 4.5 g/cm 2 for X max and smaller by 0.4-5.2 g/cm 2 for X µ max . The effect of the cross-sectional fraction of the diffractive collisions at the first interaction on X max and X µ max is approximately 5 g/cm 2 for 10 19 eV, whereas the effects of other cross-sectional fractions and the diffractive-mass spectrum at the first interaction are less than 1 g/cm 2 for 10 19 eV. Therefore, the effects of cross-sectional fractions between the pSD, tSD, and DD collisions and diffractive-mass spectrum are negligible for both X max and X µ max . The cross-section ratio of the diffractive collisions exhibit the large effect on the mean value of X max and X µ max predictions. However, the other details of diffractive collisions exhibit negligible effects and any details of diffractive collisions discussed in this paper cannot solve the differences of X max and X µ max predictions between the models.
We found that the sizes of the effect on X max and X µ max by any details of diffractive collisions are similar. This suggests that the discussions in this paper cannot explain the discrepancy between the interpretations of the cosmic-ray mass composition from X max and X µ max . For N µ , N µ originated from the tSD and the CD collisions is larger than that 24/26 from the ND collisions, and that from the pSD and the DD collisions is smaller than that by the ND collisions. N µ originated from the pSD and the DD collisions depend on the modeling of the diffractive collisions. If the current differences between the model predictions are used as the uncertainties, the effects of cross-sectional fractions and the diffractive-mass spectrum are expected to be a few % and less than a few % of N µ , respectively. The diffractive collisions may be a minor source of the muon excess problem. The effect of low energy collisions of 10 15 eV projectiles are not discussed in this paper, however, they will affect more on X µ max than on X max [34]. Moreover, these low energy collisions are also essential for N µ [35]. Thus, it is important to study the effect of low energy collisions carefully in the future. The effects of particle production mechanisms were discussed when a difference was found after assuming the same collision type and the same diffractive-mass spectrum at the first interaction. If the same diffractive-mass spectrum is used for the pSD collisions at the first interaction, the differences between the model predictions for X max , X max µ , and N µ are 21.0±1.2 g/cm 2 , 30.9±1.3 g/cm 2 , and (0.178±0.007)×10 7 at maximum for 10 19 eV, respectively. These differences are caused by two effects, the effects of particle production mechanisms and the secondary interactions, which cannot be separated in this study. There is a large difference in the diffractive-mass dependencies of N µ and this difference can be connected with the particle production mechanisms. Therefore, the effect of the particle production mechanisms is worth studying in the future.
Finally, the importance of improvements of hadronic interaction models and the experiments are discussed. The improvements of the cross-sectional fraction of the diffractive collisions in the inelastic collisions are most important since it exhibits the large effect on the mean value of X max and X µ max . There are large differences in the cross-sectional fractions of the DD and the CD collisions between the models, and they are the major sources of the differences of the cross-sectional fraction between the models. The latest models are updated using the parts of results in the LHC, but not integrating the latest results on the diffractive collisions at the LHC. It is important to update these models using the latest experimental results of the DD collisions. Experimentally, there is a difficulty in the measurements for the low diffractive-mass collisions, because most of the particles are produced in very forward regions around the beam pipes. There is a large peak of the cross-section at very low diffractive-mass regions in predictions by several models. Therefore, the experimental uncertainty of the cross-section of diffractive collisions is large due to that difficulty. Moreover, the cross sections of heavy ions collisions, i.e. lead and xenon collisions, are measured in the LHC, whereas that of the light ion collisions, which emulate the interactions between a cosmic-ray particle and an atmospheric nuclei are not measured so far. Thus, the measurements for low-diffractive mass regions and light ion collisions are required. Recently, the ATLAS and the LHCf collaborations displayed the first results of their joint analysis for very forward photons produced by the diffractive collisions [36], thanks to their complementary pseudorapidity coverage sharing the same interaction point at the LHC. The measurements by the ATLAS and the LHCf collaborations are helpful in improving the treatments of the low diffractive-mass collisions in the hadronic interaction models.