Constraining strongly coupled new physics from cosmic rays with machine learning techniques

TL;DR

This paper uses machine learning on cosmic-ray air shower simulations to set limits on strongly coupled new physics models at ultra-high energies.

Contribution

It introduces a novel approach combining Monte Carlo simulations and machine learning to constrain new physics from cosmic-ray data.

Findings

Limits on new physics cross sections established

Machine learning improves event discrimination

Method applicable to future cosmic-ray observations

Abstract

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. We then use the signal and background discrimination performance to set potential limits on the cross sections of the new physics models.

Tables

Table 1. Table 1 : Centre-of-mass collision energies corresponding to the primary particle energies considered.

$E_{\text{lab}}^{\text{primary}}$ [GeV]	$E_{\text{CM}}^{\text{p}}$	$E_{\text{CM}}^{\text{C}}$	$E_{\text{CM}}^{\text{Fe}}$	[TeV]
${10}^7$	4.3	1.3	0.6
${10}^8$	13.7	4.0	1.8
${10}^9$	43.3	12.5	5.8
${10}^{10}$	137.0	39.5	18.3
${10}^{11}$	433.1	125.0	57.8

Equations

$$N(X)=N_{\max }{\left(\frac{X-X_0}{X_{\max }-X_0}\right)}^{\frac{X_{\max }-X_0}{\lambda }}{e}^{\frac{X_{\max }-X}{\lambda }},$$

$${\mathcal{P}}_{\mathrm{new}}=A\frac{{\sigma }_{\mathrm{new}}}{{\sigma }_T(E_{\mathrm{lab}})},$$

$${\sigma }_{\mathrm{new}}\lesssim \sqrt{\frac{4{ϵ}_{\mathrm{B}}}{{ϵ}_{\mathrm{S}}^{2}N{A}^2}}{\sigma }_T\equiv f{\sigma }_T.$$

Full text

1 Introduction

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 $\mathcal{O}(1)$ 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, they predict 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 $10-100$ TeV. Examples of such theories***See also Ref. [3] 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 [10], Higgsplosion [11, 12, 13] and classicalization [14, 15].

While the former results in the production of strongly coupled resonances (such as $Z^{\prime}$ 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 $(B+L)$-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 $\mathcal{O}(100)$ particles.

To access energies of $\mathcal{O}(10)$ TeV in fundamental interactions, protons have to be collided at $\mathcal{O}(100)$ 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 instead 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 [16]. 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 [17, 18]. Previous work aimed at setting limits on new physics using cosmic-ray interactions has either predominantly focused on exploiting primary and secondary neutrinos [19, 20, 21, 22], hadronic shower particles [23], or very light resonances [24]. Here, instead, we study whether the detailed interactions of the hard process involving very heavy particles could leave an imprint strong enough to discriminate new physics from Standard Model QCD-induced backgrounds as measured at the Pierre Auger Observatory.

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 QCD-induced processes. First, we describe the simulation setup, where we use Herwig and HERBVI to generate the hard processes, followed by the simulation of 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.

2 Simulation Setup

In this section we describe all the steps in our simulation of cosmic-ray air showers from models of new physics.

2.1 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 $(B+L)$-violating sphaleron process, a heavy gauge boson $Z^{\prime}$ decaying to two Standard Model photons, and a heavy scalar boson $h^{\prime}$ decaying to two Standard Model leptons. The masses of the $Z^{\prime}$ and $h^{\prime}$ resonances are 10 TeV, with widths of 100 GeV.

The sphaleron process we study includes a change in baryon and lepton numbers of $\Delta B=\Delta L=-3$ and is of the form $qq\to 7\bar{q}+3\bar{l}+n_{V}W/Z+n_{h}h$, where $n_{V}$ and $n_{h}$ are the numbers of electroweak gauge bosons and Higgs bosons, respectively. Since it was suggested in Refs. [25, 26, 27, 28] that the production cross section for sphalerons is enhanced if produced in association with many gauge bosons, in our simulation we select $n_{V}=24$ and $n_{h}=0$. Such sphalerons could also be searched for at IceCube [29] or at high-energy proton-proton colliders [30, 31], and if observable, they could improve our understanding of the underlying mechanism of electroweak symmetry breaking [32]. 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 $2\to n$ scatterings, where $n\gg 1$.

2.2 Hard Interaction Simulation

To simulate the hard interaction for the background QCD and heavy $Z^{\prime}$/$h^{\prime}$ processes, we use the Herwig 7 [33] 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 [34, 35] tool which is implemented in Herwig. The final-state particles after hadronisation are then passed to the air shower simulation.

2.3 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 1660 water Cherenkov counters on the ground, which measure both muons, and the electromagnetic shower component. These detectors count high-energy muons and establish an estimate of the distribution of the muon density at ground level. Furthermore, there are fluorescence detectors [36], which measure the fluorescence emission of air molecules in the ultraviolet range as a function of atmospheric depth.

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 [37] air shower simulator. We use the GHEISHA [38] interaction model to treat the low-energy hadronic interactions, and the QGSJET [39] 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 $E_{\mathrm{lab}}=10^{7}$ GeV to $E_{\mathrm{lab}}=10^{11}$ 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 $\sqrt{s}\simeq\sqrt{2m_{\mathrm{p}}E_{\mathrm{lab}}/A_{\mathrm{N}}}$. Here, $A_{\mathrm{N}}$ is the atomic weight of the primary nucleus: $A_{\mathrm{N}}=1$ for a proton, $A_{\mathrm{N}}=12$ for carbon and $A_{\mathrm{N}}=56$ 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 $\rho_{\mu}$ observed at ground level, having survived through the thinning procedure. We do not apply a dethinning procedure to this observable [40]. In addition, from the distribution $N(X)$ of charged particles as a function of the shower depth $X$, we can deduce the shower maximum $X_{\text{max}}$ by performing a $\chi^{2}$-fit of a Gaisser-Hillas function [41] to the data. This function is given by,

[TABLE]

where $N_{\mathrm{max}}$, $X_{\mathrm{max}}$, $X_{0}$ and $\lambda$ 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 to determine whether these are sufficient for a meaningful discrimination between signal and background, and leave a more complicated analysis with more observables to future studies.

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 $\rho_{\mu}$ and $X_{\text{max}}$ distributions for the two approaches for a primary proton at both $10^{7}$ GeV and $10^{11}$ 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 as well, 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 except 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 $10^{7}$ GeV and a cosmic-ray carbon at $10^{11}$ 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.

3 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.

3.1 Classification of new physics events

The $X_{\text{max}}$ and $\rho_{\mu}$ 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. We also show the effect on the average $X_{\text{max}}$ and $\rho_{\mu}$ values as a function of the mass of the $Z^{\prime}$ in Fig. 4 for a primary proton at $10^{9}$ GeV, and we expect this behaviour to be representative of variations in the mass for other bins and new physics models.

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 $10^{8}-10^{11}$ GeV where the processes are kinematically possible, with the assumption that the energy and primary compositions can be determined independently of these parameters†††We note that there is a relationship between $X_{\text{max}}$ 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 $X_{\text{max}}$ and $\rho_{\mu}$. 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 [42] 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 in order to prevent overfitting the model to the training set. We then calculate the ROC curve for each new physics model on the test set, which allows one to easily obtain the background efficiency $\epsilon_{\mathrm{B}}$ for any chosen signal efficiency $\epsilon_{\mathrm{S}}$. Fig. 5 shows the ROC curves for the $Z^{\prime}$ vs QCD background classification for a primary proton at $10^{9}$ GeV for the four machine learning algorithms considered, along with the area-under-curve (AUC) scores for each algorithm. For this particular energy and primary bin, we find that the MLP has the highest AUC score. In Fig. 6, we show the output of the MLP on the signal and background test sets, where a larger output corresponds to a higher confidence from the MLP model that the particular event is a signal event. We see that the MLP is able to discriminate most events correctly.

3.2 Limits

Following the analysis in Ref. [17], 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,

[TABLE]

where $\sigma_{T}(E_{\mathrm{lab}})$ is the energy-dependent proton-air cross section, and $A=14.6$ is the average atomic mass of air. For a measured number of $N$ events, with a signal efficiency of $\epsilon_{\mathrm{S}}$ and a background efficiency of $\epsilon_{\mathrm{B}}$, we can set a 95% confidence limit by requiring that $S/\sqrt{S+B}\gtrsim 2$, where $B=\epsilon_{\mathrm{B}}N$ and $S=\epsilon_{\mathrm{S}}NA\sigma_{\mathrm{new}}/\sigma_{T}$. Assuming that the number of background events is far greater than the number of signal events, this gives the limit,

[TABLE]

The efficiencies $\epsilon_{\mathrm{S}}$ and $\epsilon_{\mathrm{B}}$ can be read off from the ROC curves. Choosing a signal efficiency of $\epsilon_{\mathrm{S}}=0.8$, the corresponding background efficiencies are shown in Table 2, with the associated limit factor $f$ for a representative number of $N$ events in each bin, which reflects the suppression of the cosmic-ray flux as a function of energy [43, 44]. In cases where very strong separation is possible, the background efficiency is set to a minimum value of $\epsilon_{\mathrm{B}}=0.05$ to ensure that the limits are conservative estimates. We also show in Table 2 the derivable limits for the case where a systematic uncertainty of 5% on the background is assumed. In this case, the 95% confidence limit is obtained from requiring $S/\sqrt{S+B+\delta^{2}B^{2}}\gtrsim 2$, where $\delta=0.05$ is the systematic uncertainty.

For proton energies in the range $10^{8}-10^{11}$ GeV, the proton-air cross section ranges from $\sim 450$ mb to $600$ mb [45]. Thus the limits on the new physics processes in Table 2 range from $\sim 80$ $\mu$b to $8$ mb. In Fig. 7 we show the $95\%$ confidence limits on the new physics cross sections as a function of the number of events for a proton at $10^{9}$ GeV where no systematic uncertainty on the background is assumed.

4 Conclusions

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 $\mathcal{O}(100)$ TeV. Thus, Auger might become more sensitive than the LHC for 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 $\rho_{\mu}$ and $X_{\mathrm{max}}$. 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.

Acknowledgements.

We would like to thank Mikael Chala and Christoph Englert for interesting discussions, and Stephen Webster for help with Herwig.