Can the ANITA anomalous events be due to new physics?

The ANITA collaboration has observed two ultra-high-energy upgoing air shower events that cannot originate from Standard Model neutrinos that have traversed the Earth. Several beyond-the-standard-model physics scenarios have been proposed as explanations for these events. In this paper we present some general arguments making it challenging for new physics to explain the events. One exceptional class of models that could work is pointed out, in which metastable dark matter decays to a highly boosted lighter dark matter particle, that can interact in the Earth to produce the observed events.


I. INTRODUCTION
High energy neutrino astrophysics could provide a window to new physics at energies beyond those accessible on Earth, thanks to the long propagation distances possible for ultra-high-energy (UHE) neutrinos, compared to charged particles. The Auger and IceCube experiments are sensitive to UHE neutrinos, with several events at PeV energies now observed by IceCube. Upper limits on the flux at higher energies are established, up to energies of 10 11 GeV.
The cross section for neutrinos to interact with nucleons through deep inelastic scattering (DIS) [1] grows with energy, and at ultra high energies above ∼ 10 8 GeV the Earth becomes opaque to neutrinos. A promising strategy then is to search for Earth-skimming neutrinos that produce charged leptons with high probability, that do not lose too much energy before they exit [2]. The vast majority of such events would be initiated by ν τ converting to τ leptons via DIS since electrons or muons would be absorbed by the ice over much shorter distances and have a negligible probability to exit compared to τ 's.
ANITA is a balloon-borne instrument that detects polarized radio emission from the electromagnetic component of cosmic ray showers in the atmosphere, reflected from ice at the South Pole, or arriving without reflection from events near the horizon. Because of phase inversion of the signal upon reflection, the difference between such events is distinguishable, and the zenith angle of the τ whose decay produced the shower can be determined to within 0.3 • .
ANITA has made four flights since 2007, and observed two anomalously steep events, henceforth called ANITA anomalous events (AAEs), apparently emerging from the Earth, of energy E ν ∼ 0.6 EeV on the first [3] and third [4] flights, corresponding to particles that traversed chord lengths of 5700 and 7200 km, respectively. 1 The mean free path of neutrinos at this energy is a few hundred km, making it exceedingly unlikely that standard model processes could explain these events [6].
Several models beyond the Standard Model (SM) have been suggested as possible explanations. These fall into three broad categories: (1) SM neutrinos of astrophysical origin could convert to beyond-the-standard-model (BSM) particles through interactions in the Earth, followed by propagation of the BSM particle until it reconverts to SM particles that initiate the observed hadronic air shower [6][7][8][9]; (2) dark matter (DM) that has accumulated within the Earth and decays to BSM particles that reconvert to SM particles shortly below the Antarctic surface could induce the air showers; or (3) an exotic flux of BSM particles, such as sterile neutrinos [10,11] incident upon the Earth interact to produce the observed particles near the Earth's surface. A possible source of such flux is the decay of long-lived DM particles [12][13][14].
A generic challenge for models in the first category is that, despite the increased probability for a BSM particle to traverse the Earth compared to EeV-energy neutrinos, there is still a large reduction in efficiency because of the need for this particle to reconvert within a thin layer near the surface of the Earth before it exits into the atmosphere. A major goal of the present work is to carefully quantify this efficiency and to show that it puts such explanations of the ANITA anomalous events at odds with null searches by IceCube and other experiments such as Auger. Models in the second category face the problem that the amount of DM particles that may have accumulated within the Earth cannot be large enough to explain the two AAEs. On the other hand, third category models have the advantage that the small efficiency factor can be overcome by assuming a larger flux of sterile BSM particles, that is relatively less constrained than the neutrino flux. Nevertheless such a large flux of sterile BSM particles can be problematic if at the same time active neutrinos or other nonsterile particles are produced. We will argue that previously proposed models of this kind are unlikely to be viable, but demonstrate an example of a model that does work. The essential ingredient is an ultraheavy DM particle that decays exclusively to another, much lighter, DM particle, whose flux is unconstrained. This boosted decay product then converts to τ 's within the Earth which in turn induce the hadronic air showers.
The paper is structured as follows: In Section II, we perform a model-independent analysis of scenarios where a neutrino scatters into weakly-interacting states that can traverse the Earth. These states will convert to τ or directly decay to hadrons seen by ANITA. Single conversion and cascade decay processes are considered here, but neither diffuse nor pointlike ν τ fluxes can explain the AAEs. Constraints on models that attempt to explain the AAEs from dark matter decays either within the Earth or in the galactic halo are discussed in Section III. In either case, we find that it is not possible to achieve a large enough flux to explain the ANITA events, while remaining consistent with limits from IceCube and other experiments. In Section IV we give an example of a class of models that can however be consistent, involving decays of a heavy subdominant component of metastable DM into the lighter dominant DM particle, that can interact with matter in the Earth to produce the anomalous events. We conclude in Section V.

II. ANITA EVENTS FROM NEUTRINO FLUX
In this section, we analyze the possibility that the ANITA events are consistent with BSM explanations, assuming that the BSM particles come from ultra-highenergy (UHE) neutrinos interacting with matter in the Earth. Other sources of BSM particles will be considered in following sections.
We assume that the ANITA anomalous events are air showers initiated by the hadronic decay of energetic τ leptons, or alternatively by the decay of BSM particles directly into hadrons. 2 The predicted number of such events is the convolution of the differential exposure of ANITA to up-going air showers, and the differential diffuse neutrino flux φ(E ν ) dE ν : Here • dA = R 2 sin θ E dθ E dφ E is the differential area of the Earth observed by ANITA, with θ E being the 2 The energetic τ 's which induce air showers above the Antarctic surface can also pass through the IceCube detector, where they could be misidentified as highly energetic muon tracks. At most three such events have been seen [6], while a rough estimate of the number of events is approximately one order of magnitude larger than the number of AAEs [6,11]. It is however possible that a detailed analysis, including in particular instrumental effects, could lead to a lower number of expected events at IceCube, thus resolving this tension and making τ -induced air showers viable (cf. ref. [6]). Hence we discuss both the case of τ -induced air showers and air showers directly induced by a BSM particle in this paper.
angle between the South Pole and a point on the surface as seen from the center of the Earth, φ E being the longitude, and R = 6371 km the Earth's radius; • θ em is the emergence angle of the τ (respectively the BSM particle that causes the hadronic air shower), measured with respect to the horizon-see fig. 1; • P obs is the probability that an UHE ν that enters the Earth with energy E ν leads to a τ (or possibly some BSM particle) that exits the Earth near the south pole at an emergence angle θ em with energȳ E τ and initiates a hadronic air shower within D = 10 km; 3 • dΩ d is the detection solid angle of ANITA. It depends onĒ τ and the τ 's decay position, which requires detailed simulations of Extensive Air Showers (EAS) [15]. Here we neglect this dependence, making the approximation that dΩ d can be integrated independently.
The live time of the two experiments was 17.25 days for ANITA I [16] and 7 days for ANITA III [4]. ANITA II was not sensitive to up-going air showers, and data from ANITA IV have not been released so far.

II.1. Conversion probability
In the following we will compute upper limits on P obs for several "process topologies," independent of a specific BSM realization. In the SM the process topology is which is understood to include τ regeneration and energy loss. The BSM process will occur in parallel with any new physics processes, and so its effects should be included alongside the latter for quantitative results, even if the SM contribution is quite small. Assuming the presence of BSM physics, several topologies are conceivable. In the first, the UHE ν converts to a BSM particle X, which subsequently converts or decays to a τ : A second possibility is to have cascade decays of several BSM particles, This includes 1 as the special case n = 1. As we will show below, for n > 1 one can increase P obs somewhat. It is also possible that X decays directly into hadrons that initiate the observed air showers. We denote this by The initial scattering process for BSM physics that converts ν to X is parametrized by the interaction length where σ BSM is the BSM cross section for scattering of ν on nucleons (or possibly electrons) and n ⊕ is the number density of target particles in the Earth. In case of more than one BSM conversion, such as cascade decays, the ith interaction length is denoted as i . Similarly, the mean free path for the standard model DIS interaction is denoted by These interaction lengths generally depend upon E ν .
0 The SM process ν → τ → hadrons. The total conversion probability P SM can be expressed as where B h 0.65 is the branching ratio for a τ to decay into hadrons (we count only the hadronic decays for initiating the kind of air showers observed by ANITA) and the chord length between the entry point of the UHE ν and the exit point of the τ is c = 2R ⊕ sin θ em . For illustration see the line labeled 0 in fig. 1. The probability distribution for the incoming neutrino to scatter into a τ at position z 1 is and the probability for τ to exit with energyĒ τ , starting from initial energy E ν , is given by [2] dP exit in terms of the τ decay length τ (E τ ) = cτ τ E τ /m τ . (Here we neglected the fact that the τ on average receives only 80% of E ν in the interaction [17].) The Earth density is taken as a constant for simplicity, ρ ⊕ = 2.7 g/cm 3 , and β τ ∼ (0.4 − 0.8) × 10 −6 cm 2 /g in the energy range of interest. We take β τ = 0.6 × 10 −6 cm 2 /g for our estimates. The probability for the τ to decay in the atmosphere at position z 3 is dP τ The probability distribution (5) is integrated overĒ τ between the initial value E ν and a minimum value E min ∼ 0.1 EeV that could still be associated with the observed AAEs.
With these analytic formulae, we can reproduce the conversion probability given by the Monte Carlo simulations in ref. [18] with vanishing ice thickness for θ em 3 • . Unlike the simulations, we neglect τ regeneration [19], which is important for larger θ em and lowĒ τ . After exiting the Earth, the τ does not suffer further energy loss, and it will decay in the atmosphere. The resulting probabilities are shown as a function of θ em in fig. 2, in which P SM drops steeply at large θ em . Even though τ regeneration increases the probabilities for large θ em , this falls far short of being able to explain the AAEs within the SM [6].
1 The BSM process ν → X → τ → hadrons. To obtain a higher conversion probability at large θ em , we consider the possibility that the UHE ν can convert to a weakly interacting particle X with a larger mean free path than neutrinos, followed by X decaying or rescattering into a τ . The conversion probability can be written as with superscript A denoting that we neglect the SM contribution to the total conversion probability here (see eq. (12) below). The probability distributions are given by dP X→τ dP ν→X is the differential probability that ν converts to X, determined by the mean free path for BSM interactions BSM (in addition to the damping by the SM contribution to the ν conversions). dP X→τ quantifies the decay  or rescattering of X (produced at z 1 ) into τ at z 2 , with mean free path X . In the rescattering scenario, there are two possible subcases: the interaction may simply be the inverse of the one that produced X (if the original interaction involved the full third generation SU(2) L doublet (ν τ , τ ) and not just the right-handed component of τ ), or it may be a different interaction. In the first case, one has X ∼ = BSM , so this can be considered as a special case of the more general situation. The general situation is like the second case, where X and BSM are independent parameters. dP exit /dE τ is to account the energy loss of τ as in eq. (7) with z 1 → z 2 .
We must also include the SM component, similar to 0 , which is however modified by the new physics because the mean free path for ν → τ is decreased by the BSM contribution, We omit the number of the process in the superscript of P (B) SM , since this same modification will apply to all three BSM scenarios. The total probability P (1A) + P

(B)
SM can be expressed in terms of dimensionless variables x and y, For an incident 3 EeV neutrino to generate a τ decaying in the lower atmosphere, at a fixed value of θ em = 27 • representative of an AAE, we find the largest possible probability 3.2 × 10 −4 , around x ∼ 0, y ∼ 0.906. The value of y is close to the chord length c (27 • ) 0.906R , which makes X-decays close to the surface of the Earth likely. At 27 • , the BSM probability P (1A) (which is an increasing function of 1/x) dominates over P (B) SM by far. Fig. 2 illustrates the dependence of P obs on the emergence angle θ em for the two choices (x, y) = (0, 0.9) and (2, 0.9), respectively, and E ν = 3 EeV. The green (solid) line is the total probability, while the other curves show the individual contributions P (1,A) and P (B) SM , as well as the SM contribution (red dashed curve). (Since the latter neglects neutrino regeneration effects, the true SM probability is higher for θ em 3 • .) 2 The BSM process ν → X 1 → · · · → X n → τ → hadrons. Rather than a single BSM particle X in the conversion, it is conceivable to have a cascade of decays or several consecutive conversions.
One expects that increasing the number of steps increases the conversion probability, since the final step can then have a higher likelihood of producing a τ near the surface of the Earth. Going to a very large number of steps is theoretically unlikely, and part of the initial E ν is lost in each decay. Therefore we give the numerical results up to five steps, fixing θ em = 27 • and c = 5785 km, by neglecting the energy loss in every conversion, except for the energy loss of τ , taking the initial ν energy to be 3 EeV. In table I, second column, we show the maximum probabilities as a function of the number of steps in the cascade. They are achieved when the conversion rate of ν → X 1 is high (small x) and all the decay lengths of the X i are roughly i ∼ c /n. The probabilities include the hadronic branching ratio B h .
3 The BSM process ν → X 1 → · · · → X n → hadrons. The AAEs do not necessarily originate from τ decays. A BSM particle decaying in the atmosphere to hadrons could also give rise to the observed events, for example a dark photon decay A → h. As in case 2, we take θ em = 27 • and E ν = 3 EeV, and neglect the energy loss of BSM particles while they traverse the Earth. The maximum probabilities are shown in column 3 of table I, where we take the branching ratio X n → hadrons as 1. These can be approximated by an analytic formula, P n n/2π D/( c sin θ em ), which works well if c /n D/ sin θ em .
In the following, we will refer to the modeldependent probability defined in this section generically as P obs (θ em ), considering it as a function of the emergence angle.

II.2. Limits on BSM assuming a diffuse ν-flux
The predicted number of ANITA events is in general given by eq. (2), where conversion probabilities for different models are given in table I. We use the exposure estimate obtained from the reflective events observed by the ANITA-I and III flights [6], For the incoming neutrino flux, we take the limits from IceCube [20] and Auger [21]: the energy times diffuse isotropic flux is 6 km −2 sr −1 yr −1 at E ν = 3 EeV. Assuming that ν τ constitutes 1/3 of the total diffuse flux, due to oscillations, this implies E ντ φ(E ντ ) 2 km −2 sr −1 yr −1 . Our conclusions below do not change significantly on varying the energy in the range of [0.1 EeV, 10 EeV].
H ref AN does not include the conversion probability P obs . The latter is a function of θ em , so we average over it in order to obtain the estimated total exposure. As mentioned above we assume that dΩ d does not depend on θ em . Then the predicted number of events from a diffuse flux in the region visible to ANITA can be estimated as To check whether the AAEs can be compatible with the BSM scenarios discussed above we perform an extended likelihood analysis. The extended likelihood takes into account the total number of observed AAEs n as well as the measured emergence angles θ i em of these events, and is defined as where f (θ em ) = P obs (θ em ) sin θ em sin θ E dθ E /dθ em dθ E sin θ E P obs (θ em ) sin θ em , is the probability distribution function to observe the emergence angle θ em , and L(0, θ em ) ≡ 1. From L one can compute the p-value for ANITA to observe at least two events with the respective emergence angles 27 • , 35 • . The angular integration region is defined to be such that in terms of the Poisson distribution P (k, µ) = µ k e −µ /k!.
In practice, the important contributions to (19) come from n ≥ 2, and for n ≥ 3 one can neglect the angular dependence and take L(n) = P (n; µ). For BSM model 1 , we find that the parameters x and y that maximize the extended likelihood are x ∼ 0, y ∼ 1.13. As expected y is slightly different from the one given in the previous section since we are considering two events with emergence angles 27 • and 35 • , respectively. The corresponding p-value is 3.7 × 10 −8 , which implies that the BSM topology 1 is inconsistent with the data at the level of 5.5 σ. As shown in table I, for the model 2 , the discrepancies are roughly the same, ∼ 5.5 σ; for the model 3 , the discrepancy can be reduced to the level of 4.6 σ. It turns out that increasing the number of steps does not correspond to an decrease of the p-value, which depends upon the total acceptance and not just on P obs .
We conclude that if a diffuse neutrino flux is assumed to be the source, BSM physics is insufficient to explain the ANITA anomalous events. The minimum incompatibility, tension at the level of 4.6 σ, is obtained for a cascade of several BSM particles X i in which the last one directly decays to hadrons with a branching fraction close to unity.

II.3. Limits on BSM from point sources
The IceCube diffuse flux upper limit for ν µ (that we presume to also hold for ν τ ), E ν φ(E ν ) 2 km −2 sr −1 yr −1 , rules out both SM and BSM explanations of the AAEs. A loophole could in principle be bright neutrino point sources. If there exists a large number of point sources with an isotropic distribution, their total flux cannot exceed the diffuse flux limit. Thus for an isotropic distribution of many point sources, one reaches the same conclusion as for the diffuse flux. On the other hand, when the number of point sources N ps is low, the single source flux limit is bounded by total diffuse flux divided by N ps ; see fig. 8 of [22]. The optimal situation for explaining the AAEs is then to have only two points sources, both in the Northern hemisphere, that give rise to the two AAEs. For a single point source of muon neutrinos, assuming a spectrum of the form φ ps ∼ E −2 , the sensitivity of IceCube is φ ps E ν = 0.13 km −2 yr −1 at E ν = 1 EeV [22]. The IceCube sensitivity to ν τ is lower because of poorer angular resolution, but we will assume that the two fluxes are equalized by oscillations. Since the angle is fixed for point sources, the unit of flux does not include sr −1 , and this changes the estimated number of observed events from eq. (16) to the form where dA ps = (h/ sin θ em ) 2 dΩ d is the differential area over which ANITA observes an air shower originating from a point source at an angle θ em and h ∼ = 35 km is its altitude. Taking dA ps sin θ em 9 × 10 4 km 2 , we can estimate in the BSM model 1 , the maximum µ ps = 1×10 −5 for two point sources in the northern hemisphere, excluding the model at the 6.5 σ confidence level. The reason for this stronger exclusion is the assumption of the spectrum φ(E ν ) ∼ E −2 ν [22]. It leads to fewer events at E ν ∼ 1 EeV from anisotropic point sources than does the diffuse flux.
This suggests another possible loophole (apparently taken in ref. [7]), namely to have a non-monotonic spectrum with a peak at energies ∼ 1 EeV. Ref. [22] computes the sensitivity and discovery potential of IceCube to such sources, as a function of E ν and declination δ (see fig. 4 in that reference). For δ = 30 • , similar to the angles of the AAEs, the differential flux limits are φ ps E ν ∼ 6 × 10 3 km −2 yr −1 at E ν = 0.5 EeV and φ ps E ν ∼ 1.6 × 10 5 km −2 yr −1 at E ν = 1 EeV. Taking the former flux we find that BSM model 1 is in tension with the ANITA data at only 1.6 σ, while with the latter, all the models including the SM are consistent with the data.
However at declination δ = 30 • and E ν ∼ 1 EeV, the Auger experiment is more sensitive than IceCube, since events from this direction are down-going or Earthskimming [23]. A limit of E 2 ν φ ps < 140 EeV km −2 yr −1 is reported for an E −2 ν flux from a point source. For a flux concentrated at E ν ∼ 1 EeV, the limit is presumably weaker, but probably not greatly so since Auger's sensitivity to differential spectra is peaked around 1 EeV. It is reasonable to assume the differential limit at ∼ 1 EeV is a factor of 5 weaker than the E 2 ν flux. Taking the differential limit of E ν φ ν < 700 km −2 s −1 at E ν ∼ 1 EeV, this would be in 3.1 σ tension for being able to explain the AAEs. We conclude that nonpower-law flux point sources at high declination, while perhaps not ruled out, can only be marginally consistent with the AAEs.
Finally, intense transient sources lasting 24 hours, coming from the northern hemisphere and peaked in the EeV range, might have been missed by both Auger and IceCube, so that we are unable exclude this as a possible explanation, either within the SM or coming from new physics. A study of this possibility is beyond the scope of the present work.

III. ANITA EVENTS FROM DECAYING DARK MATTER
Since very large fluxes are required to overcome the small transmission probability through the Earth, one is motivated to look for an exotic source of BSM particles that could produce AAEs, whose flux could be much larger than the experimental limits on the neutrino diffuse or pointlike fluxes. Decaying dark matter has been suggested. In one scenario, the signal is dominated by decays of DM that has been trapped in the Earth [13], while in another it is DM decaying in the galactic halo that dominates [14]. We will argue that both of these are in varying degrees of tension with complementary constraints. However it is possible that decaying DM gives rise only to highly boosted BSM particles, that interact in the Earth to produce the AAEs. We point out in the following Section that this third scenario, heretofore unconsidered, can be viable.

III.1. DM decays in the Earth
DM particles χ can accumulate inside the Earth if they have a sufficiently large cross section σ χN for scattering on nucleons. The current upper limit from XENON1T, extrapolated to DM mass m χ = 1 EeV, is σ χN < 10 −39 cm 2 [24], giving a mean free path in the Earth of λ > 0.3 × 10 10 km, which is 5 × 10 5 times greater than the radius of the Earth. The collection efficiency is expected to be = πR ⊕ /(2λ) 3 × 10 −6 (accounting for the average distance through the Earth traversed by a DM particle). The total mass of DM which the Earth can collect is where ρ local DM 0.3 GeV/cm 3 is the local DM density, v rel 220 km/s is the DM velocity and t 4.5 Gyr is the age of the Earth.
On the other hand, in order to explain the N A =2 anomalous events within t A = 34 days, in a solid angle Ω A 4 km 2 /(4πR 2 ), and with DM lifetime τ χ , one needs a total DM mass accumulation of where P A is the probability for a DM decay product to propagate through the Earth and produce a τ that initiates a hadronic shower in the lower atmosphere. We take an optimistic value of P A ∼ 10 −3 based on our previous estimates. The lifetime is constrained by Ice-Cube searches for neutrinos from decaying DM [25] to be τ χ 2 × 10 27 s. This is the weakest limit (corresponding to χ → νν) which in general depends upon the final state particles in the decay. The required mass is five orders of magnitude greater than the amount of DM that could collect inside the Earth, so even with an unrealistically large value for the conversion probability P A ∼ 1, the scenario is still ruled out.
One might ask whether this negative conclusion could be evaded by considering strongly interacting dark matter, that is stopped in the Earth before reaching underground detectors. However the larger cross sections needed have been excluded by other experiments and astrophysical constraints, as summarized in ref. [26]. That reference further closes any remaining loophole by showing that too much heat would be deposited in the Earth relative to the measured heat flow, if the direct detection bound is significantly violated; hence our upper bound on is robust.

III.2. DM decays in the galactic halo
Even without collecting in the Earth, decays in the DM halo of our galaxy can provide a flux of sterile BSM particles that is far bigger than the experimental limit on the diffuse neutrino flux. This could potentially overcome the small survival probability for traversing the Earth, if these BSM particles can convert to ν τ , τ or some other particle that produces a hadronic shower in the atmosphere. This is the strategy of ref. [14], in which a bosonic DM particle χ is presumed to decay into sterile right-handed neutrinos ν R , that have a small mixing angle θ ∼ 0.01 with active neutrinos. Through this mixing, ν R can interact with nucleons in the Earth to produce ν τ , leading to the observed AAEs.
However because of the mixing, there are also direct decays χ → ν R ν τ , with a branching ratio of θ 2 , rendering the decays themselves detectable by IceCube, regardless of ν R interactions in the Earth. Ref. [14] explains the observed AAEs with m χ = 20 EeV and τ χ = θ 2 × 10 27 s. The effective lifetime for the subdominant decays χ → ν R ν τ is therefore simply 10 27 s, independent of θ. This is in tension with IceCube limits on decaying DM in the χ → νν channel [25]. Although the latter reference stops slightly short of m χ = 1 EeV, the limits have been extended to higher masses [27], giving τ χ > 2 × 10 28 s, at m χ = 20 EeV. As a consequence, the expected number of AAEs is at least an order of magnitude too low. Ref. [14] also identifies a region of parameter space around m χ = 5 × 10 4 EeV, τ χ = θ 2 × 10 25 s, corresponding to an effective lifetime of 10 25 s for the χ → ν R ν τ channel. This too is disfavored, since the limits on the lifetime from [27] are two orders of magnitude stronger than what one would need to obtain an order one number of AAEs.

IV. A VIABLE MODEL: DM DECAYS TO DM
The previous example suggests that one might achieve a viable explanation if DM decays solely into a BSM particle whose interactions in the Earth can produce an exiting τ , with sufficiently high probability. Here we present a phenomenological model with a heavy metastable DM particle Ψ decaying to a highly boosted, light DM particle χ, that can interact in the Earth to produce τ leptons.
The heavy component Ψ must have a mass m Ψ ∼ 1 EeV to produce the AAEs. This is above the mass limit for thermal production of DM, but there exist several mechanisms for producing super-heavy DM, for example by gravitational particle production [28][29][30] or reheating/preheating [31][32][33][34] at the end of inflation, freeze-in [35], dilution by entropy production [36], or bubble collisions during a phase transition [32]. We suggest a scenario where Ψ may be a subdominant component of the total dark matter, while the lighter χ particle constitutes most of the DM and gets its relic density from thermal freezeout.
There is a Z 4 charge carried by the new particles that guarantees the stability of χ. The scale of the dimension- (If desired, a generalization of the model that brings Λ below the Planck scale is to introduce a Z 2n discrete symmetry under which so that Ψ → (2n − 1)χ through an operator of dimension 3n,ψχ(χ c χ) n−1 , with n ≥ 3.) For superheavy DM with m Ψ ∼ 3 EeV, the decay product χs are highly energetic, and can interact with nucleons in the Earth by the diagram shown in fig. 3. A novel feature is that the internal L τ can go on shell, so the cross section is logarithmically sensitive to the width of ν τ , Γ ν . Since the decay rate of ν τ is negligible, this width is dominated by the interaction rate of ν τ in the Earth, where n N ∼ = 3.3×10 24 /cm 3 is the density of nucleons and the energy dependence of σ νN , valid around E ν ∼EeV, is taken from ref. [18]. At high energies, the scattering is dominated by lowvirtuality W or Z exchange. We find that the chargedcurrent cross section is (see appendix A for details) whereŝ = xs is the parton-level invariant, with momentum fraction x. The function ln( √ŝ /Γ ν ) denotes the leading behavior at large values of the argument, but we use the more exact expression (A7) for the following estimates.
Averaging (28) over the parton distribution functions of the nucleons gives a result that depends upon the inert Higgs doublet mass m φ , because the cross section is dominated by small x, whose minimum value is x min = m 2 φ /s. The enhancement factor E takes the place of ln( √ŝ /Γ ν ), and E varies from 12, 000 for m φ = 100 GeV to 3, 000 for m φ = 1 TeV, as shown in fig. 4. We find that it is accurately described by the analytic fit

IV.2. ANITA anomalous events
Our model can be considered to be of the type 1 in our classification of section II.1, in the limit where x → 0 (see eq. (14)), which indicates that the BSM particle X = Ψ is already present upon entry of the Earth. The scattering process actually produces two sources of τ : first from the primary vertex giving the doublet L τ shown in fig. 3, and second through the fast decays of the produced φ → χL τ . In section II.1 we showed that the probability for these τ s to be observed is maximized when y = x /R ⊕ ∼ = 1. In the remainder we will assume this condition is satisfied, to partly constrain the parameter space of the model.
The cross section (29) leads to a scattering length in the Earth of taking the density of nucleons to be n N = 3.3×10 24 /cm 3 . Imposing y = 1 fixes y τ as a function of m φ via the relation (30). We achieve lengths of order the Earth radius with reasonable values of y τ 1.2, as shown in fig. 5. The procedure of section II.1, with x = 0 and y = 1 in process 1 gives a probability of P ∼ 5 × 10 −4 at large θ em , like for the AAEs. Then using the likelihood method of section II.2, we find that the best-fit flux of χ particles is Φ χ = 2230 km −2 yr −1 sr −1 . For a rough estimate, we can compute the flux coming from DM decays by assuming that there is a constant mean density of χ, n χ = 3n Ψ (τ u /τ Ψ ) from the decays, where  n Ψ = ρ Ψ /m Ψ = f Ψ Ω cdm ρ crit /m Ψ and τ u is the age of the universe. Here f Ψ 1 is the abundance relative to that if Ψ constituted all of the DM.
The isotropic flux is then which determines the Ψ lifetime when equated with the best-fit flux: A more quantitative estimate using the J-factor for decays in the galactic halo gives a similar estimate. One needs a very high scale Λ ∼ f 1/4 Ψ × 10 24 GeV in the dimension-6 operator, that could be lowered by taking a larger value of n in the generalized version of the model. We do not concern ourselves here with trying to build a UV complete model, but instead emphasize that the general framework may be promising for understanding the AAEs via new physics.

IV.3. Relic density and other constraints
It is interesting that the same interaction that induces the scattering of χ in the Earth can also determine its thermal abundance through the annihilations χχ → L τLτ . The cross section at threshold is As a rough estimate, we equate this to the nominal level (σv) 0 ∼ = 2 × 10 −26 cm 3 /s that leads to the observed relic density. In conjunction with the determined value of y τ for fitting the AAEs, this imposes a relation between m χ and m φ that is plotted as the solid curve in fig. 5. It is approximately fit by m χ ∼ = 0.3 m φ − 22 GeV, consistent with the requirement that m φ > m χ so that φ → χL τ decays occur. The inert doublet φ could be pair-produced though electroweak interactions at the LHC, leading to τ pairs with missing energy. This is the same signature as for stau pairs in the MSSM that decay to τ s and neutralinos, τ →χ 0 τ . ATLAS has searched for this signal in the 8 TeV data [37] and more recently CMS has searched in the 13 TeV run [38]. So far the sensitivity does not reach the expected production cross section, though in a future run at the HL-LHC it is projected that m φ (in the guise ofτ ) will be excluded below 650 GeV [39].

V. CONCLUSIONS
The origin of the two anomalous events observed by the ANITA balloon experiment remains unexplained. In this paper we explored if and how these AAEs might be explained by new physics, within the assumption that BSM states transverse the Earth and converts to τ or hadrons in the atmosphere. Our conclusions are: a) A diffuse flux of ν τ cannot explain the AAEs, and it is disfavored by ∼ 5σ by considering single conversion and cascade decay models. For point sources, the tension is reduced to ∼ 3σ, since the flux limit is weaker. b) DM decaying inside the Earth cannot account for the AAEs because not enough DM could have been accumulated inside the Earth during its history. c) A large flux of "sterile" BSM particles whose interactions in the Earth can produce a τ near the Antarctic surface could in principle explain the AAEs. Such a flux of sterile particles could originate from the decay of DM in the galactic halo. One must ensure however that the DM decay does not also induce a flux of non-sterile particles, such as active neutrinos, in excess of the limits from IceCube. As an illustration of a viable model, we discuss the case of a metastable EeV-scale dark matter particle that decays exclusively to a lighter dark matter particle χ. The latter can in turn scatter within the Earth to produce a τ air shower, consistently with all constraints. The relic density of the stable DM can be explained by the same interaction that induces the air showers, which can also lead to events at the LHC resembling supersymmetricτ production and decay. and WX are supported by the European Research Council grant NEO-NAT. JC is supported by the Natural Sciences and Engineering Research Council of Canada.

Appendix A: DM-nucleon cross section
Here we give some details of the computation of the charged-current scattering process shown in fig. 3. The spin-averaged squared matrix element takes the simple form |M| 2 = y 2 τ g 4 2 (p 1 · p 3 )(p 2 · p 3 )(p 4 · p 5 ) (2p 2 · p 5 + m 2 W ) 2 (p 1 · p 3 + δ 2 ) 2 (A1) in the limit of massless particles, where for the diagram with virtual ν. The parton-level cross section in terms of the three-body phase space can be written as (A3) where we choose coordinates such that (in the center-ofmass frame) p 1 = E 1 (sin α cos φ, sin α sin φ, cos α) p 2 = − p 1 p 3 = E 3 (0, 0, 1) p 4 = E 4 (sin γ, 0, cos γ) p 5 = − p 3 − p 4 (A4) with E 1 = E 2 = √ŝ /2 and cos γ = 1 +ŝ − 2 √ŝ (E 3 + E 4 ) 2E 3 E 4 being fixed by energy-momentum conservation, since Hereŝ is the partonic center-of-mass energy squared, where x is the quark momentum fraction and m N the nucleon mass. Only for x close to x min = m 2 φ /s does it matter that the external particles are massive; otherwise the masses can be neglected since they are 1 EeV, the energy scale of interest.
The momentum squared in the ν τ propagator is 2p 1 ·p 3 = 2E 1 E 3 (1 − cos α), so the integral over cos α is dominated by cos α ∼ = 1 at high energies. One finds that Since the most sensitive cos α dependence is in this factor, we take cos α → 1 in the remaining part of |M| 2 , which depends upon E 4 . The integral over E 4 can be done analytically, and in the limit E 3 m W it gives The remaining integral over E 3 can also be done analytically, where we evaluated them at scale Q = m W in accordance with the low virtuality of the internal W boson indicated by (A6).