Constraining strongly-coupled new physics from cosmic rays with machine learning techniques
Cosmic rays interacting with the atmosphere allow for the probing of fundamental interactions at ultra-high energies. We thus obtain limits on strongly-coupled new physics models via their imprints on cosmic ray air showers. Using the Monte Carlo event generators Herwig and HERBVI, and the air shower simulator CORSIKA, to simulate such processes, we apply machine learning algorithms to the simulated observables to discriminate the events arising via new physics from the QCD background, before using the signal and background discrimination performance to set potential limits on the cross sections of the new physics models.
The recent discovery of the Higgs boson [1, 2] was the last missing piece to establish the Standard Model of particle physics as an effective theory describing interactions at TeV, thereby confirming the paradigm that nature can be described to a high precision with perturbative quantum field theory in such an energy range. However, many UV completions of the Standard Model predict fundamental modifications to that paradigm, in particular that the theory transitions from a weakly-coupled into a strongly-coupled regime not too far beyond the electroweak scale, e.g. in the range TeV. Examples of such theories***See also Ref.  for selected resonance cross sections and simplified models with mediators to strongly coupled sectors [4, 5] at 100 TeV proton-proton collisions. are composite Higgs models [6, 7, 8, 9], little string theories , Higgsplosion [11, 12, 13] and classicalization [14, 15].
While the former results in the production of strongly-coupled resonances, e.g. or heavy scalar particles, which are usually short-lived and decay into a small number of Standard Model particles, the latter two examples result in the production of a multi-particle final state where the energy of the phenomenon is subsequently distributed over a plethora of particles, not unlike the -violating sphaleron process of the Standard Model. If such processes can be realised with appreciable probabilities, separating signals with a small number of final state objects from large QCD-induced Standard Model backgrounds is a significantly bigger task in a collider environment than for final states with particles.
To access energies of TeV in fundamental interactions protons have to be collided at TeV center-of-mass energies to account for the fact that the individual quarks and gluons in the proton only carry a fraction of the proton’s energy. In the absence of a proton-proton collider that can access such energies, we focus on ultra-high-energy cosmic rays to study whether strongly-coupled new physics can be probed in their interactions with the atmosphere. When a highly energetic proton hits the atmosphere large momentum transfers occur which eventually give rise to an extended air shower of photons, hadrons and leptons. As a whole this air shower is a highly complex object which can arguably obfuscate the hard process that initiates the shower.
In recent years, however, for high-energy events at the LHC, novel analysis techniques have been devised to study jets, complex collimated sprays of hadrons, and their substructure . The remarkable success of these techniques, e.g. in discriminating electroweak scale resonances from QCD-induced backgrounds, makes it plausible that one can apply similar techniques to the study of cosmic ray air showers in separating Standard Model processes from decays of heavy resonances or multi-particle phenomena .
Thus, we use machine learning techniques to analyse the structure of air showers to discriminate the kinematic distributions that heavy resonances would leave compared to Standard Model induced processes. First, we describe the simulation setup, where we use Herwig and HERBVI to generate the hard processes and simulate the air shower using CORSIKA. Then, we show the effects of the new physics models on two air shower observables compared to the background QCD process. Finally, we train machine learning algorithms to classify the events and use this to derive simple estimates of the limits on the cross sections of these processes.
In this section we describe all the steps in our simulation of cosmic ray air showers from models of new physics.
New Physics Processes
To represent possible processes that can arise in non-perturbative solutions to and UV completions of the Standard Model, we consider a -violating sphaleron process, a heavy gauge boson decaying to two Standard Model photons, and a heavy scalar boson decaying to two Standard Model leptons. The masses of the and resonances are 10 TeV, with widths of 100 GeV.
The sphaleron process we study includes a change in baryon and lepton number and is of the form , where and are the numbers of electroweak gauge and Higgs bosons respectively. As it was suggested in [18, 19, 20, 21] that the production cross section for sphalerons is enhanced if produced in association with many gauge bosons, in our simulation we select and . Such sphalerons could also be searched for at IceCube  or at high-energy proton-proton colliders [23, 24], and if observable, they could improve our understanding of the underlying mechanism of electroweak symmetry breaking . At the level of observability of a high-energy collision on the surface of our atmosphere, such a multi-particle production process mimics the kinematic features induced by processes from Higgsplosion or classicalization. Thus, we will take the sphaleron as representative of models with enhanced production mechanisms for elementary scatterings, where .
Hard Interaction Simulation
To simulate the hard interaction for the background QCD and heavy / processes, we use the Herwig 7  Monte Carlo event generator. Herwig collides the two protons, computes the partonic interaction, and simulates the parton shower as well as the hadronic phase transition. To generate the sphaleron processes, we use the HERBVI [27, 28] tool which is implemented in Herwig. The final-state particles after hadronisation are then passed to the air shower simulation.
Air Shower Simulation
A cosmic ray air shower is the phenomenon of observable secondary particles produced by a high-energy cosmic ray colliding with the upper atmosphere. In the following we briefly describe the different stages of such a shower.
The process starts with a cosmic ray heading towards the Earth, which we call the primary particle. In principle any particle could be the primary particle in the collision. However, in this work we focus on nuclear matter, and as representatives of the table of elements we choose a proton, carbon and iron.
Usually ordinary high-energy QCD describes the hard interaction when a primary particle hits an air nucleus in the upper atmosphere. However, the probability for the particular process is determined by its cross section, and in this study we also consider the other processes described above for the hard interaction. Regardless of the physics guiding the hard interaction, there will be a QCD parton shower as well as a hadronic phase transition.
As the interaction located in the upper atmosphere is directed downwards, a cascade of secondary interactions will follow. This is the air shower. The secondary particles will produce bremsstrahlung of any form. Furthermore, they will collide with other air molecules feeding the cascade until the total energy is diluted and the shower dies away.
In experiments like the Pierre Auger Observatory several detectors are used to capture a signal from the air shower. Firstly, there are muon chambers distributed on the ground. These chambers count high-energy muons and establish an estimate of the distribution of the muon density. Furthermore, there are fluorescence detectors , which count the number of charged particles as a function of atmospheric depth by measuring their induced Cherenkov radiation.
To analyse new physics in cosmic ray air showers we need to simulate the whole interaction chain described above. To do so, we process the particles generated from Herwig and HERBVI with the CORSIKA  air shower simulator. We use the GHEISHA  interaction model to treat the low-energy hadronic interactions, and the QGSJET  interaction model to treat high-energy hadronic interactions. A thinning procedure is applied to the shower simulation, which restricts the number of particles in each shower stage as a computational requirement.
The incoming primaries that we simulate have zero inclination and interact at a height of 18 km, with energies ranging from GeV to GeV. The corresponding centre-of-mass (CM) collision energies for the hard interaction, which consists of a proton in the cosmic ray nucleus interacting with a proton in the air nucleus, are given by . Here, is the atomic weight of the primary nucleus: for a proton, for carbon and for iron. For the carbon and iron nuclei, the energy is assumed to be evenly distributed amongst its nucleons. Table 1 shows the values of the collision energies corresponding to the primary particles that we consider.
From the simulation results, we extract the number of muons observed at ground level, having survived through the thinning procedure. We do not apply a dethinning procedure to this observable . In addition, from the distribution of charged particles as a function of the shower depth , we can deduce the shower maximum by performing a -fit of a Gaisser-Hillas function  to the data. This function is given by,
where , , and are to be determined from the fit. In principle there is no reason why one should not include more observables usually studied in air shower experiments, such as the risetime. However, for the purposes of this study we limit it to just these two observables.
As a test of the reliability of using Herwig, with its capability for generating new physics processes, to generate the hard interaction and then processing the events with CORSIKA, we can also generate the full primary-to-air-shower chain for the QCD events with CORSIKA alone by using its own hard process simulation. We find that there is good agreement between them, and in Fig. 1 we show a comparison of the and distributions for the two approaches for a primary proton at both GeV and GeV, which spans the energy range we consider. The differences are small, although the distributions are not identical, but for the purposes of this study we will ignore any small systematic uncertainties that may arise due to the use of Herwig as the hard process generator.
Since we are not only interested in ordinary proton-proton interactions, but actually study nucleus-air collisions we need to model the additional nucleonic complexity. As the air is at rest and its binding energy is low compared to the energies we are interested in we regard it as a stationary proton. However, we cannot use such a simple ansatz for the high-energy primary particle. In principle, we might view the interaction of a nucleus with a proton in the air as a proton-proton interaction. However, we have to take the nucleonic remainder of the now-destroyed primary into account. There are two extremes we can study. We could assume that the impact was so fast that the nucleus stays untouched but with one fewer proton. On the other hand, we could assume that the nucleus is destroyed and completely fragments into its proton and neutron components. A comparison of both approaches is shown in Fig. 2 for a cosmic ray iron at GeV and a cosmic ray carbon at GeV. We find the differences between the two extremes are small, and so for the rest of this study we consider a completely fragmented remainder nucleus.
Results and Limits
In this section we show the effect of the new physics models on the two air shower observables, and train machine learning algorithms to classify the events into signal and background classes. From this, we derive possible limits on the cross sections of the new physics processes.
Classification of new physics events
The and distributions for each new physics model in each energy and primary bin are presented in Fig. 3, along with the background QCD distributions. The distributions shown have each been calculated from 1000 simulated points using a Gaussian kernel density estimate, with the cross showing the maximum of the distribution and the two contours enclosing 68% and 95% of the data. It is clear from the plots for carbon and iron in Fig. 3 that the new physics effects are washed out by the interactions of the remainder nucleus, and thus the parameter distributions are almost identical. Therefore, we only consider the four proton bins in the energy range GeV where the processes are kinematically possible, with the assumption that the energy and primary compositions can be determined independently from these parameters†††We note that there is a relationship between and the composition of the primary. Indeed, one could be tempted to interpret variations of the primary composition as being potential signs of new physics. However, for the sake of this analysis we ignore these effects and their systematics, and assume that the primary compositions and energies are well-determined., so that these parameters can be used for the new physics classification.
In each of these energy and primary bins, we train a machine learning algorithm to independently classify the three new physics models vs the QCD background in the two-dimensional parameter space of and . The machine learning algorithms that we use are a linear discriminant analysis (LDA), a quadratic discriminant analysis (QDA), a support vector machine (SVM) and a multilayer perceptron (MLP), and we use Scikit-learn  for their implementation.
For each new physics model and bin combination, the 2000 data points (1000 signal and 1000 background) are split into training, validation and test sets. We perform hyperparameter scans over the important hyperparameters of each algorithm, and the algorithm which has the highest accuracy on the validation set is used. We then calculate the ROC curve for each new physics model on the test set, which allows one to easily obtain the background efficiency for any chosen signal efficiency . Fig. 4 shows the ROC curves for the vs QCD background classification for a primary proton at GeV for the four machine learning algorithms considered, along with the area-under-curve (AUC) scores for each algorithm.
Following the analysis in Ref. , we can use a simple counting procedure in each proton bin to set a limit for the cross section of each new physics process in terms of the proton-air cross section. The probability for a new physics process to occur in the collision of a proton with the air can be expressed as,
where is the energy-dependent proton-air cross section, and is the average atomic mass of air. For a measured number of events, with a signal efficiency of and a background efficiency of , we can set a 95% confidence limit by requiring that , where and . Assuming that the number of background events is far greater than the number of signal events, this gives the limit,
The efficiencies and can be read off from the ROC curves. Choosing a signal efficiency of , the corresponding background efficiencies are shown in Table 2, with the associated limit factor for a representative number of events in each bin, which reflects the suppression of the cosmic ray flux as a function of energy [36, 37]. In cases where very strong separation is possible, the background efficiency is set to a minimum value of to ensure that the limits are conservative estimates.
Ultra-high-energy cosmic rays interacting with the atoms in the atmosphere are natural high-energy hadron colliders. In comparison with the LHC the event rate recorded through the fly eye at Auger is much smaller. However, the collision energies recorded reach beyond TeV. Thus, Auger might become more sensitive than the LHC in new physics scenarios that are realised at energies outside the kinematic reach of the LHC, and for cross sections that are comparable with QCD interactions. Examples of such scenarios would be potentially unsuppressed sphaleron production or a strongly-coupled dark sector. We find that it is possible to set a model-independent limit on the cross sections of such new physics processes by considering their effects on cosmic ray air showers via the observables and . Using multi-variate data analysis techniques, a strong separation between signal and QCD background interactions can be achieved. However, based on our classification approach, this is only possible for proton primary particles as the effect is washed out for heavier primaries.
We would like to thank Mikael Chala and Christoph Englert for interesting discussions, and Stephen Webster for help with Herwig.
-  Aad G. et al., Phys. Lett.B, 716 (2012) 1.
-  Chatrchyan S. et al., Phys. Lett.B, 716 (2012) 30.
-  Arkani-Hamed N., Han T., Mangano M. Wang L.-T., Phys. Rept., 652 (2016) 1.
-  Englert C., NordstrÃ¶m K. Spannowsky M., Phys. Rev.D, 94 (2016) 055028.
-  Becciolini D., Gillioz M., Nardecchia M., Sannino F. Spannowsky M., Phys. Rev.D, 91 (2015) 015010 [Addendum: Phys. Rev.D92,no.7,079905(2015)].
-  Kaplan D. B., Georgi H. Dimopoulos S., Phys. Lett., 136B (1984) 187.
-  Caracciolo F., Parolini A. Serone M., JHEP, 02 (2013) 066.
-  Barnard J., Gherghetta T. Ray T. S., JHEP, 02 (2014) 002.
-  Ferretti G., JHEP, 06 (2014) 142.
-  Antoniadis I., Arvanitaki A., Dimopoulos S. Giveon A., Phys. Rev. Lett., 108 (2012) 081602.
-  Khoze V. V. Spannowsky M., Nucl. Phys.B, 926 (2018) 95.
-  Khoze V. V. Spannowsky M., Phys. Rev.D, 96 (2017) 075042.
-  Khoze V. V., Reiness J., Spannowsky M. Waite P., Precision measurements for the Higgsploding Standard Model (2017).
-  Dvali G., Giudice G. F., Gomez C. Kehagias A., JHEP, 08 (2011) 108.
-  Dvali G. Gomez C., JCAP, 1207 (2012) 015.
-  Marzani S., Soyez G. Spannowsky M., Looking inside jets: an introduction to jet substructure and boosted-object phenomenology (2019).
-  Brooijmans G., Schichtel P. Spannowsky M., Phys. Lett.B, 761 (2016) 213.
-  Ringwald A., Nucl. Phys.B, 330 (1990) 1.
-  Khoze V. V. Ringwald A., Nucl. Phys.B, 355 (1991) 351.
-  Khoze V. V. Ringwald A., Phys. Lett.B, 259 (1991) 106.
-  Tye S. H. H. Wong S. S. C., Phys. Rev.D, 92 (2015) 045005.
-  Ellis J., Sakurai K. Spannowsky M., JHEP, 05 (2016) 085.
-  Ellis J. Sakurai K., JHEP, 04 (2016) 086.
-  Ringwald A., Sakurai K. Webber B. R., JHEP, 11 (2018) 105.
-  Spannowsky M. Tamarit C., Phys. Rev.D, 95 (2017) 015006.
-  Bellm J. et al., Eur. Phys. J.C, 76 (2016) 196.
-  Gibbs M. J., Ringwald A., Webber B. R. Zadrozny J. T., Z. Phys.C, 66 (1995) 285.
-  Gibbs M. J. Webber B. R., Comput. Phys. Commun., 90 (1995) 369.
-  Abraham J. et al., Nucl. Instrum. Meth.A, 620 (2010) 227.
-  Heck D., Knapp J., Capdevielle J. N., Schatz G. Thouw T., CORSIKA: A Monte Carlo code to simulate extensive air showers (1998).
-  Fesefeldt H., The Simulation of Hadronic Showers: Physics and Applications (1985).
-  Kalmykov N. N., Ostapchenko S. S. Pavlov A. I., Nucl. Phys. Proc. Suppl., 52 (1997) 17.
-  Stokes B. T., Cady R., Ivanov D., Matthews J. N. Thomson G. B., Astroparticle Physics, 35 (2012) 759.
-  Gaisser T. K. Hillas A. M., International Cosmic Ray Conference, 8 (1977) 353.
-  Pedregosa F. et al., J. Machine Learning Res., 12 (2011) 2825.
-  Fenu F., The cosmic ray energy spectrum measured using the Pierre Auger Observatory in The Pierre Auger Observatory: Contributions to the 35th International Cosmic Ray Conference (ICRC 2017) 2017 pp. 9–16 [PoSICRC2017,486(2018)].
-  Gora D., Universe, 4 (2018) 128.
-  Abreu P. et al., Phys. Rev. Lett., 109 (2012) 062002.