New Constraints on Millicharged Particles from Cosmic-ray Production

We study the production of exotic millicharged particles (MCPs) from cosmic ray-atmosphere collisions which constitutes a permanent MCP production source for all terrestrial experiments Our calculation of the MCP flux can be used to reinterpret existing limits from experiments such as MACRO and Majorana on an ambient flux of ionizing particles. Large-scale underground neutrino detectors are particularly favorable targets for the resulting MCPs. Using available data from the Super-K experiment, we set new limits on MCPs, which are the best in sensitivity reach for the mass range $0.1 \lesssim m_{\chi} \lesssim 0.5$ GeV, and which are competitive with accelerator-based searches for masses up to 1.5 GeV. Applying these constraints to models where a sub-dominant component of dark matter (DM) is fractionally charged allows us to probe parts of the parameter space that are challenging for conventional direct-detection DM experiments, independently of any assumptions about the DM abundance. These results can be further improved with the next generation of large-scale neutrino detectors.


I. INTRODUCTION
The remarkable success of the Standard Model (SM), along with null results for new physics at the LHC, strongly suggests that if new physics exists below the TeV scale it can only be weakly coupled to SM degrees of freedom. While nearly decoupled from the SM, such a dark sector would likely leave its strongest imprint on SM degrees of freedom commensurate with its own dynamical energy scales [1]. It is interesting to note that the MeV -GeV regime both contains many SM particles (e.g. muons, mesons, and nucleons) and hosts a number of persistent anomalies, including the anomalous magnetic moment of the muon [2][3][4]. Furthermore, this energy range is interesting from a phenomenological point of view as it allows for many novel and complementary search strategies that can be used to probe the dark sector. For instance, new physics can often be efficiently probed by fixed target experiments [5][6][7][8][9] with high intensity electron [10,11] and proton beams [12][13][14][15][16][17][18][19][20] where dark sector decays in the upper atmosphere for m χ in the few MeV to few GeV regime. For this, we adopt a minimal MCP model that is based on only two assumptions: 1. The new particle χ couples to the SM photon with a strength Q χ = × e; we remain agnostic as to the origin of this charge.
2. The new particle is stable; this is a natural consequence if Q χ is the smallest (non-zero) charge in the dark sector.
As these features are relatively generic, MCPs can be thought of as a useful representative example of a stable dark sector particle with which to benchmark the impact of neutrino telescopes. In particular, since we consider only primary production in what follows, our constraints apply (possibly conservatively) to any model that satisfies the above two assumptions.
In addition to being a useful benchmark model, MCPs are of interest because of their potential impact on 21cm cosmology (potentially explaining the EDGES anomaly [45][46][47][48]), and their natural appearance in models of light dark matter (DM) interacting with the SM via a massless dark photon [49][50][51][52]. Boosted millicharged DM can also potentially explain a reported excess in direct detection experiments [53,54]. Running parallel to these more cosmological motivations, the lack of constraints in the few MeV -few GeV regime has also motivated the proposal of dedicated detectors such as MilliQan [22,23] and FerMINI [18].
One immediate consequence of our calculation of the cosmic ray-induced MCP flux is that existing bounds on a naturally occurring flux of MCPs can be converted into constraints on Lagrangian parameters and m χ , where the MCP charge is Q χ = × e. In fact, multiple such bounds already exist in the literature, but have never been translated into the − m χ plane because the relation between , m χ and the flux from cosmic rays had never been made explicit. Examples include constraints from MACRO [56,57], Kamiokande-II [64], LSD [65], CDMS [66] and Majorana [58]. Interestingly we find the resulting constraints to be roughly competitive with those from existing collider experiments, but sub-dominant to reported bounds from neutrino experiments [16,59,60].
Here we point out that neutrino telescopes can set new leading bounds on MCP couplings in the 100 MeV -500 MeV regime based on existing data, surpassing the reach of fixed target experiments with neutrino detectors. We demonstrate this point explicitly by providing novel constraints based on published analyses by the Super-Kamiokande (Super-K, SK) collaboration searching for the diffuse supernova neutrino background (DSNB) [55]. Our results, summarized in Fig. 1, suggest that future neutrino telescopes could be able to act as the leading probe of MCPs in this mass regime. Furthermore, our results can be recast as a study of millicharged strongly interacting dark matter (SIDM) [67,68], allowing us to Exclusion limits for MCPs from cosmic-ray interactions (SK, red solid), obtained using analysis results of the diffuse supernova neutrino background search in Super-K [55], as well as sensitivity projections for an improved SK analysis (SK+, red dotted) and near-future Hyper-K (HK, red dashed). We also display new limits (blue) from recasting data of MACRO [56,57] and Majorana [58]. Previous limits from fixed target (SLAC MilliQ [10,11], MiniBooNE [16,59], ArgoNeuT [60]) and collider experiments [6,[61][62][63] (as compiled in Ref. [60]) are shown for comparison. explore a region of interesting parameter space that cannot be easily studied by conventional underground directdetection experiments (see Section V for a more detailed discussion).
Our study establishes neutrino telescopes as an important probe of the same MCP parameter space that motivated the proposal of MilliQan [22,23] (and the similarly designed FerMINI [18]), studies on MCP bounds from neutrino experiments [16,60,69], and the proposed MCP DM explanation of the EDGES 21 cm anomaly [46,47,70]. In the context of MCP DM we emphasize that our constraints are independent of the fractional composition of composition of DM [71][72][73], for other searches see e.g. Ref. [70,[74][75][76]. Finally, the explicit calculation of the MCP flux from cosmic rays presented here will enable the use of neutrino telescopes as a robust platform for studying MCPs, free from cosmological assumptions. This has connections to charge quantization, which is itself connected to, but does not necessarily preclude [49], the existence of magnetic monopoles [77], Grand Unification [78][79][80][81], and quantum gravity [82].
This paper is organized as follows. In Section II, we discuss the production of mesons from cosmic-ray collisions in the upper atmosphere. In Section III, the MCP flux from meson decays is calculated. We then discuss the detection of MCPs in neutrino telescopes in Section IV, and outline the kinematics of detecting MCPs. In Section V, we discuss the millicharged SIDM and the constraints and projections that we can place based on our analysis.

II. COSMIC-RAY MESON PRODUCTION
Cosmic rays produce a sizeable number of mesons from interactions in the upper atmosphere, whose subsequent decay produces a continuous flux of MCPs. While the problem can be studied numerically with Monte Carlo simulations, we present here a semi-analytic treatment, allowing us to transparently illustrate the role of key ingredients.
Incoming cosmic rays are isotropically distributed on the sky, with the associated flux typically quoted in terms of intensity, [I CR ] = GeV −1 cm −2 s −1 str −1 [83]. In our analysis we take this quantity as implemented in DarkSUSY [84], based on Ref. [85], and focus on the dominant component of cosmic rays, free protons. For convenience, we will instead express the intensity in terms of the center-ofmass boost γ cm for CR protons impinging on atmospheric protons at rest and thus introduce I CR (γ cm ) = s/m p , s is the Mandelstam variable for the pp collision, and m p is the proton mass.
Taking into account that all incoming cosmic rays are eventually absorbed by the atmosphere, the amount of primary mesons m produced in these collisions is approximately determined by the ratio of the inclusive cross section σ m for pp → mX with other particles X to the total inelastic cross section for protons passing through atmospheric matter. We note that this is a rather conservative estimate for the total production of mesons given that all final states in these primary interactions tend to trigger further cascades when interacting with the atmosphere, resulting, among others, in a large multiplicity of (lowerenergy) meson states. Here we neglect these contributions, which could be studied with a dedicated Monte Carlo simulations of air showers. We model all interactions in the upper atmosphere as pp collisions and therefore take the elastic cross section to be σ in (pp), whose dependence on γ cm is given in Ref. [83]. The resulting meson flux from cosmic-ray collisions in the upper atmosphere is then given by where Ω eff ≈ 2π is the effective solid angle from which MCPs can arrive at the detector 2 , and P (γ m |γ cm ) represents the probability to get a meson with boost γ m in 2 By rescaling the muon's stopping power [26,83], we estimate that the energy loss of MCPs in the Earth's crust (standard rock) is the lab frame. The latter can be conveniently estimated (see Appendix A 1) from the differential production cross section with respect to x F ≡ p L /p max , where p L is the longitudinal momentum and p max is the maximum possible momentum: Here α = ± denotes the two different possible contributions, see Eq. (A3), and dσ m /dx F is a function of γ cm and x F (γ m ).
The meson-production energy spectrum thus depends on both the total meson cross section, σ m (γ cm ), and the differential cross section with respect to x F , or equivalently on P (γ m |γ cm ). These quantities must be specified across a significant range of γ cm to reflect the large range of cosmic-ray energies, and we do so by interpolating between existing data for selected values of fixed γ cm (see Appendix B). Although both σ m and P (γ m |γ cm ) influence the final resulting MCP flux, we find that the production cross section (which is also better measured) has a much stronger effect than the differential distribution.
In principle, all possible mesons originating from pp interactions and leading to MCPs (i.e. those with substantial electromagnetic decay modes: π, η, ω, ρ, J/ψ, Υ etc.) as well as direct production via Drell-Yan should be considered. For light MCPs produced via π 0 → γχχ a combination of SLAC's milliQ experiment [10] and LSND's search for electron-like scattering events [86] already strongly restricts the MCP parameter space [16]. We therefore restrict our discussion to the case of heavier roughly 50 MeV/km for ∼ 10 −2 . While for the range of and energies that we are interested in here MCPs interact too strongly to penetrate the entire Earth, they are not significantly impeded to reach the detector when originating from the upper hemisphere.
To keep our discussion of meson production tractable, we focus on the dominant η, light vector, and J/ψ mesons. While we have also quantitatively considered Υ meson as well as direct Drell-Yan production, we found these contributions to be negligibly small (six orders of magnitude smaller than J/ψ) since in addition to smaller cross sections these processes require more energetic cosmic rays (with correspondingly much smallerfluxes).
In Appendix B we analyze and fit the available experimental data for η, ρ, ω, φ and J/ψ, finding the total production cross section σ pp→ηX (γ cm ) as well as the standard spectrum "shape parameterization" dσ m /dx F . In Fig. 2 we display the resulting differential cosmic ray intensity I CR (γ cm ) multiplied by σ m (γ cm ). The shape of these curves is determined by the competition between a rising inclusive cross section and a sharply falling cosmicray flux, and illustrates which parts of the cosmic ray spectrum predominantly contributes to a given meson species.

III. MCP FLUX FROM MESON DECAYS
Upon constructing Φ m as outlined above, we can find the associated flux of MCPs from meson decays by folding the meson flux with the unit-normalized spectrum of MCPs in the lab frame, P (γ χ |γ m ), and weighting by the decay branching ratio where the factor of 2 accounts for the contribution from bothχ and χ. The quantity P (γ χ |γ m ) can be calculated from first principles, at leading order in , as where Γ is the decay rate for m → χχ and dΓ/dγ χ is the differential rate with respect to the MCP boost, both evaluated in the lab frame (see Appendix A 2). Anticipating MCP detection, we define the integrated "fast-flux" of MCPs satisfying γ χ ≥ γ cut as where γ cut is set by the relevant experimental threshold. In Fig. 3 we display the mass-dependence of this quantity for several choices of γ cut . The choice γ cut = 1 corresponds to the full integrated MCP flux, as relevant for low-threshold ionization experiments, while γ cut = 6 is adequate for experiments with an electron recoil threshold of T min = 16 MeV (as relevant for the physics analysis of Super-K discussed below).
Fast-flux of MCPs Φcut due to meson decays as a function of MCP mass, mχ, for three different choices of γcut.
The spectrum for γcut = 1 is the full integrated MCP flux. The meson mass thresholds are clearly visible, stemming from η, ω/ρ, φ, and finally J/ψ (sequentially from left to right).

IV. DETECTING MCPS IN LABORATORIES
As alluded to in the introduction, MCPs can deposit ionization energy directly within detectors, which can be used as a probe of MCP couplings [56][57][58][64][65][66]. Lacking an explicit calculation of cosmic rays as MCP source, previous searches have avoided discussing the mass of the incident MCPs, and instead presented constraints on an ambient MCP flux as a function of the fractional charge . Our study allows us to directly translate these results (and future searches) into limits on as a function of m χ , thus making direct contact with Lagrangian parameters. We discuss the details of this translation in Appendix C 2, and show our results in Fig. 1. We find that ionization experiments are competitive with constraints from colliders around the 100 MeV regime, but quickly become subdominant as m χ is increased. Significant improvement in detector exposure for ionization searches is expected in future experiments, and our results establish a quantitative baseline that can be used to estimate the potential future impact of upcoming projects such as LEGEND [87]. We note that for MCPs with large charges of 10 −1 , as relevant for ionization searches, effects of attenuation when passing through Earth to reach typical detector depths of ∼ 1 km of standard rock (i.e. few km water-equivalent) become significant (see e.g. Ref. [26]). Since we do not attempt a detailed translation of ionization bounds in this work, and this region is already well constrained by collider searches (which are not sensitive to attenuation), we do not consider these effects here.
Electron scattering inside Cherenkov detectors, with recoils in the 10 MeV range, is a powerful probe of MCPs [15] (see also [26]). Counting electron-like events with recoil energies, T e = 2m e (E e − m e ), between T min and T max naturally introduces a windowed cross-section Eq. (C1) which can be well approximated (see Fig. 9) as Here, α is the fine-structure constant, Θ is the Heaviside step function and γ cut ≈ 0.6 2T min /m e + 0.4 2T max /m e . The total resulting number of χ − e scattering events N eχ for a given experiment is where N e is the number of electrons within the detector's fiducial volume and t is the data collection period.
Using data sample and analysis results from the DSNB search in Super-K [55] we can place stringent new limits on MCPs from cosmic-ray production. This search looked at inverse beta decaysν e p → ne + with a positron recoil energy 16 MeV < T e + < 88 MeV, corresponding to γ cut ≈ 6, effectively reducing the background from cosmic ray muon spallation at lower energies (note that Cherenkov detectors do not directly differentiate between electrons and positrons). As detailed in Appendix C 2, the results of the likelihood analysis performed in Ref. [55] (dashed curve in their Fig. 19) can directly be employed to constrain the recoil electron spectrum from MCPs.
We show our resulting bounds on MCPs in Fig. 1. The limits are competitive with accelerator-driven searches across the 0.1 m χ 1.5 GeV range. In the range 0.1 m χ 0.5 neutrino telescopes exceed the leading constraints from both MiniBooNE and ArgoNeuT, demonstrating the potential of neutrino telescopes as a "downstream" detector. The results are quoted in terms of events per year per 22.5 kt of water, corresponding to the fiducial volume of Super-K as employed in Ref. [55].
Upcoming large neutrino experiments will be able to further improve on these results. Additional background suppression due to improved neutron tagging will be possible in an upcoming Super-K upgrade with gadolinium doping [88], which we denote as SK+ and assume a reach of ∼ 0.6 events/22.5 kt-yr in Fig. 1. With a fiducial volume of 190 kt the near future Hyper-K water Cherenkov experiment [35] can further improve significantly on results of Super-K. In Fig. 1 we indicate this by assuming a year-long exposure and a sensitivity (in terms of SK's fiducial volume) of ∼ 0.1 events/22.5 kt-yr.
Other near future experiments with sizable fiducial volumes, such as DUNE (40 kton, liquid argon) [37] and JUNO (20 kton, liquid scintillator) [36], will complement water-based Cherenkov detectors as probes of the DSNB [89], and hence can also serve as probes of atmospherically-produced MCPs. The solar and spallation backgrounds at low energies are expected to be present in DUNE's DSNB search [90], with an expected resulting energy cut-off of ∼ 20 MeV for a search as in SK. Due to favorable detector configuration and application of pulse-shape discrimination techniques JUNO can perform DSNB search over a wider energy range [91], down to ∼ 10 MeV, albeit with overall statistics still considerably lower than that of Hyper-K.

V. STRONGLY-INTERACTING DM
A major focus of DM studies is the direct detection of DM particles using terrestrial detectors [95], typically placed underground. However, these searches depend on the local flux of the DM particles of interest that could reach the experiment. It has long been noted that when the DM-SM particle (mostly nuclei and electrons) cross section is large enough, this flux would be significantly attenuated [67,96]. The class of DM models featuring such large interactions with ordinary matter is often referred to as strongly interacting DM (SIDM 3 ).
In Ref. [68], DM-SM interactions through a dark photon kinetically mixing with U (1) Y are studied, focusing on the terrestrial effects on direct detection experiments. In this case, DM scattering with electrons becomes more important than scattering with nuclei, so this is what we focus on in the following. Millicharged DM with cross sections larger than a critical value would have its average energy attenuated and be unable to trigger a detectable signature in ground-based direct-detection experiments [68]. Above this critical cross section, there is a window of available parameter space where MCPs could constitute a sub-dominant component of DM ( 0.4% to avoid cosmological constraints [71][72][73]98]), from hereon referred to as the millicharged SIDM window. New balloon and satellite experiments have been recently proposed [68] to further explore this window, which could accommodate interesting DM models that could potentially explain the EDGES anomaly [45][46][47][48]99].
In this section, we recast our bounds and projections on MCPs to explore this SIDM window. The results in Fig. 4 are shown in terms of a "reference cross section" typically employed for direct detection experiments to compare the sensitivity reaches of different experiments. For millicharged DM, this is given by 4 where µ χe is the reduced mass of the electron and χ and FIG. 4. Constraints and sensitivity reaches that can cover the millicharged SIDM window, including our new bounds from recasting Super-K data [55] (red) and the projection for Super-K (red, dashed) and Hyper-K (red, dotted). We also show the new bounds (blue) from recasting data of MACRO [56,57] and Majorana [58]. Existing acceleratorbased constraints [6,10,11,16,[60][61][62][63] and direct-detection limits [68,[92][93][94] are also shown.
for semiconductor or noble-liquid targets from the local DM flux, taken to be αm e [68]. Above the "ground-based (GB) direct detection critical cross section", one can see a regime enclosed by bounds from accelerator-based experiments [6,10,11,16,[60][61][62][63], constraints from the above-atmosphere detector (RRS) [92], a rocket experiment (XQC) [93,94], and underground direct detection experiments [68]. We plot the bound of Super-K and sensitivity reaches for Super-K+, and Hyper-K. We do not consider bounds based on the MCP acceleration from astrophysical sources [26,38,39,100], since they rely on additional assumptions beyond local DM abundance. The constraint on the ultralight dark-photon mediator is also not shown since it is not an essential ingredient for minimal MCPs.
Our results in Fig. 4 establish new constraints on the millicharged SIDM window. It is important to note that our bounds and projections are independent of any assumption about which fraction of the DM is millicharged. Further, for reference cross sections below approximately 10 −17 cm 2 our results are insensitive to attenuation in the Earth, given that cosmic-ray produced MCPs have much higher energy than that of the local DM flux.

VI. SUMMARY
We have considered MCP production from standard cosmic rays interacting with the atmosphere. This closes a gap in the MCP literature and constitutes a permanent MCP production source for all terrestrial experiments. We presented the first translation of longstanding bounds on an ambient MCP flux into bounds on the MCP charge as a function of its mass m χ , and demonstrated that large-scale underground neutrino experiments are particularly well suited for probing previously inaccessible parameter space. Using existing limits from Super-K's DSNB search we have placed new limits on MCPs for 0.1 m χ 1.5 GeV, which for m χ 0.5 GeV exceed the sensitivity of fixed target experiments such as MiniBooNE and ArgoNeuT. These new limits are highly relevant also in scenarios where MCPs constitute an SIDM component because they are i) independent of the DM fraction made of such MCPs and ii) probe a part of the parameter space that cannot be readily tested with conventional direct-detection experiments. The results presented here will be further improved with upcoming large-scale neutrino experiments, and, since we only consider primary meson production, can likely be further strengthened by a more detailed modeling of cosmic-ray showers.

VII. ACKNOWLEDGEMENTS
We would like to thank Dr. Appendix A: Cosmic-ray production kinematics

Boost of produced mesons
The lab-frame energy of a meson produced in a collision, E m , can be written as E m = γ cm E + γ cm β cm P , where curly script variables refer to center of mass frame quantities. We can re-write this expression in terms of x F = P /p max ("Feynman-x"), where p max = 1 2 √ s(1 − m 2 m /s) is the largest possible longitudinal momentum allowed by kinematic constraints; x F therefore varies from −1 (backwards pointing) to +1 (forward pointing). Written in terms of x F our formula is given by (A1) Since p T = P T p max we can neglect it in our analysis. Therefore we can obtain γ m (x F ) from This equation can be inverted to yield two branches x corresponding to the two solutions of the quadratic equation.

Meson decay to millicharged particles
Most of the decay modes we consider involve two-body final states. For example, in the case of the J/ψ the differential decay in the rest frame of the parent meson is mono-energetic dΓ/dE χ ∝ δ(E χ − 1 2 m J/ψ ) (in this subsection, curly letters refer to meson rest-frame quantities). Upon boosting to the lab frame this becomes a box distribution, Box(E χ |γ m ), of width E (+) Equivalently, in terms of the MCP's lab frame boosts, we have whereγ χ andβ χ are the boost and velocity of the MCP in the meson rest frame. In the case of ρ 0 and φ, the dominant decay mode is also a two body final state (e.g. ρ 0 → χχ). For ω the SM branching ratio for ω → π 0 + − is roughly ten times larger than ω → + − [83], but this decay mode is only accessible for m χ ≤ 1 2 (m ω − m π ) ≈ 325 MeV, as opposed to m χ ≤ 1 2 m ω ≈ 390 MeV for the direct two body decay. We therefore neglect this decay mode 5 which will underestimate the MCP flux by a factor of ∼ O(few) in the window 275 MeV m χ 325 MeV. and focus instead on ω → χχ. The branching ratio for 5 Including ω → π 0 χχ would involve a chiral perturbation theory calculation analogous to the one performed for η → γχχ.

Appendix B: Atmospheric meson production rate
Our treatment of meson production in the upper atmosphere is data driven and centers mostly around the ratio of σ(pp → mX)/σ inel (pp) which varies as a function of center of mass energy. Although we have tried to inform our fits using data across a wide range of center of mass energies (or equivalently γ cm ) there is a limited window of "important" center of mass boosts that is determined by the competition between a rising inclusive cross section and a sharply falling cosmic-ray flux as a function of γ cm (the typical I CR ∼ E −2.7 scaling translates to roughly I CR ∼ γ −4.5 cm ). This is illustrated in Fig. 2 where we see that the relevant ranges are γ cm between 1.5-5 for η mesons, between 1.5-10 for ρ (and ω) and φ, and between 3-25 for J/ψ.
The rest of this section is devoted to our parameterization of the available inclusive cross section data, which we separate into a discussion of σ m (γ cm ) and P (γ m |γ cm ). It is important to note that although P (γ m |γ cm ) is poorly constrained by the data we were able to find, its impact on our sensitivity curves is marginal; this is because the total number of MCPs produced is independent of this quantity. In contrast, although it has a relatively comprehensive dataset, the production cross section σ m (γ cm ) in the window of maximal production (as shown in Fig. 2) can have a substantial impact on the MCP signal (bounds on scale as 4 √ signal) because it alters the total number of MCPs produced. We therefore anticipate that the uncertainties in the production cross section are the dominant source of error in our analysis (at the level of ∼ O(few). √ s/mp. The data is taken from Refs. [103][104][105][106][107][108] and fitted using the piece-wise procedure described in the text; the smooth curve is Eq. (B1). .

η mesons
Eta meson production in pp collisions has been most extensively measured in the near-threshold regime for the exclusive process pp → ηpp [103,108]. Near threshold this is the only available channel, such that this cross section can be taken as a reasonable estimate of the total inclusive cross section. Further away from threshold bona fide measurements of the inclusive cross section are scarcer but we have identified four measurements in the literature at √ s =3.17, 27.45, 38.8, and 53 GeV [103][104][105][106][107]. We split the available data into two subsets, near-threshold exclusive production (defined as pp → ppη measurements for √ s ≤ 3 GeV) and far-fromthreshold inclusive data (defined as pp → ηX for √ s > 3 GeV). We fit the near-threshold data for σ η ( √ s) with the function f (x) = a(x − 2.42) b x c where x ≡ s pp /GeV. For the far-from-threshold data we instead use g(x) = a(1 + |b|/(x − 2.42) 2 ) log 2 (x). In both cases a weighted linear regression to the data was performed. Using the best fit values for both fits, and demanding that the function is continuous we find where the numerical factor comes from the relationship √ s = 2m p γ cm = (1.876 GeV) γ cm . The functions f (x) and g(x), with their best fit values, are given by and γ = 1.59 is chosen such that Eq. (B1) is continuous; the fit is shown vs. the data (with error bars when available) in Fig. 5.
For the differential cross section dσ m /dx F , measurements at NA27 [104] strongly suggest an exponential dis-tribution, where c η depends on γ cm . Measurements from NA27 at √ s = 27.5 GeV (corresponding to γ cm = 14.6) fix c η ≈ 9.5 [104]. One generally expects that c η will be a monotonically increasing function of γ cm , and that c η > 0. The simplest functional form that satisfies these expectations, and agrees with the measurement of [104] is c η (γ cm ) = 9.5 + (slope) × (γ cm − 14.6) ; (B5) we take slope≈ 1 2 . We checked that our sensitivity to MCPs from experiments such as SK are relatively insensitive to the value of the slope parameter.

Light vector mesons
The production cross section for the ρ 0 meson is relatively well measured [103,107,108]. Like the η meson we perform a best fit analysis with the function g(x), but without weighted errors. We find the data to be reasonably well described by A comparison between the available data and our smooth fit is shown in Fig. 6.
We found that the data for σ(pp → ρX) had a much better coverage than the corresponding ω production cross section, and where there are measurements of the ω cross section it is nearly identical to the ρ cross section. We therefore estimated the ω cross section σ ω (γ ω ) σ ρ (γ ρ ).
For the φ meson we find that the functional form h(x) = a(1 + |b|/(x − 2.896) 2 ) −1 x c gives a reasonable fit to the data [103,107,108]. After an unweighted regression we find that σ φ is well described, c.f. Fig. 7, by σ φ (γ cm ) = (0.01 mb)(1.876γ cm ) 1.23 Like the η meson, the longitudinal momentum distributions for the vector mesons were more difficult to find in the literature, and we rely on a single measurement at √ s = 27.5 GeV [104] which shows the x F dependence to be described by Eq. (B4) with c V = c ρ = c ω = c φ ≈ 7.7. We expect this value to be smaller at lower center of mass energies and so take c V = 7.7 + 5.7 13 (γ cm − 14.6) , which, just like the c η , should be viewed as a cartoon of the behaviour of dσ/dx F as a function of γ cm rather than a faithful representation.

J/ψ mesons
For the J/ψ mesons we found two convenient summaries of the available data: one from E-739 ( Fig. 7 in Ref. [109]) and one from HERA-B (Fig. 8 of Ref. [110]). The HERA-B compilation includes measurements at significantly higher center of mass energies. For comparison we plot both sets of data in Fig. 8 where we see that the HERA-B compiled data is roughly consistent with that from the E739 paper, but suggests a steeper growth with rising center of mass energy. We use the best fit to the former to calculate σ m (γ cm ).
For the differential distribution we used the standard parameterization of dσ J/ψ /dx F [111] dσ J/ψ Like c η , the fit parameter c J/ψ depends on the center of mass energy, and like c η the precise value of c J/ψ has a relatively mild effect on the fast-flux of MCPs.  [109] and HERA-B's compilation [110]. The solid curve is digitized from the fit presented in Ref. [110] while the dashed curve is taken from Ref. [109]. For our sensitivity analysis we use the solid curve from [110]. [111] whereas experiments at higher energies find larger values such as c J/Ψ ≈ 6 for √ s ≈ 40 GeV [112]; we did not find a robust set of measurements of c J/ψ spanning the entire range of γ cm relevant for cosmic ray pp collisions. For simplicity, and because our final results are relatively insensitive to the details of the x F distribution, we take c J/ψ to vary linearly with γ cm , In this case, there is data at lower center of mass energies that suggests this formula is a reasonable interpolation.
Appendix C: MCP signal in experiments

MCP-electron scattering
The detection of MCPs is dominated by soft scattering from electrons as can be readily understood by considering the differential scattering cross section which, being mediated by photon exchange, scales as dσ/dQ 2 ∼ 1/Q 4 . For elastic scattering from a target of mass M , the momentum transfer is given by Q 2 = 2M (E − M ), where E is the total recoil energy of the target. The cross section is therefore maximized by scattering off the lightest target possible with the lowest possible recoil energy. In practice, experimental considerations such as detection efficiency and background reduction will set a minimum electron recoil energy which will, in turn, dictate the detection cross section for that given experiment. We will therefore consider a windowed cross section for electron recoils with kinetic energy, T e = (E e − m e ) between T min and T max , or equivalently with momentum transfers between Q 2 min and Q 2 max , Since the four-momentum transfer is directly related to the recoil energy in the lab frame, T e = E e − m e , via Q 2 = −2(p e − p e ) 2 = (2m e E e − 2m 2 e ) = 2m e T e , this is equivalent to demanding that Q 2 ≥ 2m e T min . In the center of mass frame the maximal momentum transfer is given when the scattering is back-to-back such that where P e is the electron's momentum in the center of mass frame, In terms of lab frame variables this implies that where the approximation holds provided m χ m e γ χ . The main consequence of Eq. (C3) is that the lower bound of integration in Eq. (C1) is given (at leading order in m e /E χ ) 6 by Q 2 min = max(2m e T min , 4m 2 e (β χ γ χ ) 2 ), and the upper bound is given by Q 2 max = min(2m e T max , 4m 2 e (β χ γ χ ) 2 ). The effect of this approximation onσ eχ is stated as Eq. (6) in the main text and illustrated in Fig. 9. In summary, the primary driver of the windowed cross section is whether or not the incident MCP is sufficiently boosted to kick the electron above the detection threshold. In principle, the thresholds of large neutrino detectors can be rather low, a few MeV in case of SK, and as low as 200 keV for Borexino. We choose, however, a much higher threshold of ∼ 15 − 16 MeV, that removes all the events generated by solar neutrinos, so that background counting rates reduce to O(few) per year.
In the case of eχ scattering the event shape spectrum is determined by the differential cross section with respect to recoil energy, dσ/dT e ∝ 1/T 2 e , and the incident flux of MCPs. We have confirmed that this shape is very similar to the case of a neutrino spectrum described by temperature of T ν 5 MeV [55], allowing us to readily employ those results.
We emphasize that the DSNB limits from Super-K are given in terms of limits on the scattered positron event rate as a function of the effective neutrino temperature T ν from supernova emission. We note, for the reader's convenience, that in [55] there are two bounds quoted: one for an ensemble of supernovae of different temperatures and one for a single supernova temperature. We use the latter, because it more closely mimics our signal as is clearly shown in Fig. 10 (the diffuse ensemble would be relatively flat as a function of energy).

Ionization experiments
Ionization is a very low threshold process and so we use the full flux (integrated over all boosts γ χ ) of MCPs for ionization experiments; this corresponds to γ cut = 1 as shown in Fig. 3; we denote this total flux by Φ(m χ ). To translate existing bounds on an ambient MCP flux in the literature we demand that where −2 Φ ( m χ ) corresponds to the γ cut = 1 curve in Fig. 3 (i.e. the integrated MCP flux generated in the upper atmosphere), and Φ ion ( ) is the joint exclusion curve obtained by combining data from MACRO [56,57] and Majorana [58] as shown in Fig. 7 of [58]. We then solve for for each value of m χ which determines a critical value of, c (m χ ), above which MCPs are excluded. FIG. 10. Comparison of event shapes for MCP elastic scattering off electrons and inverse beta decay from supernova background neutrinos. The MCP signal was obtained by folding the differential scattering cross section dσeχ/dTe against the cosmic-ray induced MCP flux. The supernova background curves correspond to E 2 ν /(e Eν /Tν + 1) where Eν = Te + 1.3 MeV; these correspond to the fixed temperature profiles used in Ref. [55] as can be readily verified by reproducing their Fig.  19.