Explaining the ANITA Anomaly with Inelastic Boosted Dark Matter

We propose a new physics scenario in which the decay of a very heavy dark-matter candidate which does not interact with the neutrino sector could explain the two anomalous events recently reported by ANITA. The model is composed of two components of dark matter, an unstable dark-sector state, and a massive dark gauge boson. We assume that the heavier dark-matter particle of EeV-range mass is distributed over the galactic halo and disintegrates into a pair of lighter -- highly boosted -- dark-matter states in the present universe which reach and penetrate the Earth. The latter scatters {\it in}elastically off a nucleon and produces a heavier dark-sector unstable state which subsequently decays back to the lighter dark matter along with hadrons, which induce Extensive Air Showers, via on-/off-shell dark gauge boson. Depending on the mass hierarchy within the dark sector, either the dark gauge boson or the unstable dark-sector particle can be long-lived, hence transmitted significantly through the Earth. We study the angular distribution of the signal and show that our model favors emergence angles in the range $\sim 25^\circ -35^\circ$ if the associated parameter choices bear the situation where the mean free path of the boosted incident particle is much larger than the Earth diameter while its long-lived decay product has a decay length of dimensions comparable to the Earth radius. Our model, in particular, avoids any constraints from complementary neutrino searches such as IceCube or the Auger observatory.


I. INTRODUCTION
The ANtarctic Impulsive Transient Antenna (ANITA) Collaboration recently reported two anomalous upwardmoving cosmic-ray-like events whose inferred energy lies in ∼ 0.5 − 1 EeV [1,2]. Within the framework of the Standard Model (SM), one could envision, for example, an EeV-scale tau-neutrino traversing the Earth and converting to a tau lepton which is responsible for producing an Extensive Air Showers (EAS) in the atmosphere after escaping the surface of the Earth. However, given the typical amount of energy that such a neutrino would have to carry, and the interaction strength of the latter with nuclei, its propagation through the Earth on large distances is extremely unlikely to let a subsequent tau lepton escape the Earth at a sizable emergence angle (see also FIG. 8 for notational conventions). Indeed, a boosted tau neutrino having an energy of O(1) EeV has a transmission probability lower than 10 −6 to induce an EeV-scale tau lepton even if one includes multiple τ − ν τ regenerations [2], and a dedicated study in Ref. [3] has concluded that SM explanations are excluded by (at least) 5σ confidence. Furthermore, Ref. [4] explored the hypothesis of a tau-neutrino origin together with a dedicated acceptance estimate, and reached the conclusion that the ANITA events disfavor a diffuse tau-neutrino interpretation in comparison with the Auger [5] and IceCube [6] upper limits. This suggests that new physics should come into play in order to accommodate the anomaly in a self- * heurtier@email.arizona.edu † doojinkim@email.arizona.edu ‡ jcpark@cnu.ac.kr § seodongshin@yonsei.ac.kr consistent manner, stimulating the theory community to propose various beyond-the-SM (BSM) scenarios.
For completeness, we should mention that the ANITA anomalous events might be explained in pure SM interpretations which were recently proposed. A downwardmoving high-energy cosmic rays might indeed be able to produce the ANITA's inverted polarity signals through coherent transition radiation from the geomagneticinduced current in the associated air showers [7] or could be reflected by sub-layer of the ice-sheet of the Antartic [8]. Such proposals remain for now speculative and still require experimental confirmations, in particular from the ANITA Collaboration itself, and therefore we shall focus in this work on the possibility that such events are a sign of new physics.
Among other attempts, several trials assumed the existence of non-SM neutrino species in order to address the signal. For example, the authors of Ref. [9,10] adopt sterile neutrinos induced by Ultra High Energy (UHE) cosmic rays as "agents" propagating through the Earth in order to explain the ANITA events, while Ref. [11] interprets the signal with lepto-quark-mediated sterile neutrinos. By contrast, Ref. [12] hypothesizes super-heavy right-handed neutrinos captured in the Earth which subsequently decay to active neutrinos in the vicinity of the Earth surface. Some of the latest attempts imagine new physics scenarios within the supersymmetric framework. Examples include interpretations utilizing nextto-lightest supersymmetric tau lepton (i.e, "stau") [3], long-lived bino in R-parity violating supersymmetry [13], and neutrino-initiated supersymmetric sphaleron transitions [14]. In contrast to the approaches taken in preceeding works, the authors of Refs. [15] and [16] explored the possibility that the flux of boosted particles necessary to explain the signal is due to the decay of a super-arXiv:1905.13223v1 [hep-ph] 30 May 2019 heavy dark-matter particle in the Milky Way. In [15] the ANITA anoumalous events are explained with a darkmatter species decaying into a pair of right-handed neutrinos, which later on convert into taus inside the Earth and produce high-energy EAS in the atmosphere. In a complementary manner, also using the decay of a superheavy dark-matter particle, the authors of Ref. [16] proposed an interpretation in terms of the Askaryan emission created by elastic scattering of feebly interacting particles within the ice sheet. The authors of Ref. [17] finally reviewed different topologies of processes which could explain the signal with EAS and proposed a model in which a heavy dark-matter decays into a boosted light component which converts into taus after scattering within the Earth. As another recent BSM interpretation, the authors of Ref. [18] argued that the reflection of radio pulses could reproduce the ANITA signals via axionphoton conversion.
As a generic feature of most of the aforementioned interpretations, we point out that previous works rely on the interaction which might exist between new-physics states and ordinary SM neutrinos in order to induce the relevant experimental signatures. The new states thus effectively "transport" the neutrinos through the Earth to enhance the chance that they convert into leptons in the vicinity of the surface.
While producing leptons inside the Earth from the conversion of neutrinos is probably the most natural SMfriendly interpretation at low energy, the situation is radically different at very high energy since the neutrino is not a good candidate to help propagating through the Earth on large distances. Moreover, there are no compelling reasons that neutrinos are necessary ingredients to explain the ANITA anomaly. Indeed, the only requirement for ANITA to detect a signal is that a significant electric field is locally produced around the detector. Therefore, as far as EAS are concerned, the production of any boosted hadronic final states in the atmosphere, whatever its parent particle is, would suffice to propose an interpretation for the anomalous events reported by the Collaboration.
Furthermore, the connection between the neutrino sector and new-physics states has been thoroughly probed by different experiments in the past and eventually leads to tensions between existing experimental constraints and the parameter choices required in order to fit the data. In particular, the production of highly boosted -and therefore very long-lived -τ leptons close to the Earth surface, which can pass through terrestrial detectors such as IceCube before decaying in the atmosphere, may leave sizable signatures which should already have been observed.
In light of these considerations, we propose a darkmatter scenario not accompanying the production of SM leptons inside the Earth (in contrast to that in [15] and [16]) in order to explain the two ANITA anomalous events. The model that we consider is built upon a non-minimal dark sector scenario, containing • φ, a super-heavy scalar dark matter, • χ 1 , a light fermionic dark-matter component, • χ 2 , an unstable fermionic dark-sector state, and • X, a dark gauge boson interacting with the aforementioned fermionic states and mixing with the SM photon through kinetic mixing.
We assume that the whole signal process sketched in FIG. 1 begins with the decay of φ into a pair of χ 1 's over the galactic halo in the present universe. 1 Each χ 1 then acquires a large Lorentz boost factor due to the mass gap between the two dark matter species. 2 Such boosted χ 1 scatters off material in the Earth to the unstable dark-sector state χ 2 , via a t-channel exchange of the dark gauge boson X. By construction, χ 2 is heavier than χ 1 , and hence decays back to χ 1 along with a quark pair whose hadronization eventually creates an EAS for the ANITA detector as long as the production happens in the atmosphere. In this scenario, the particle χ 1 behaves like a neutrino, as it is not only invisible but also ultra-relativistic, while the model parameters associated with χ 1 are relatively less constrained by existing experiments. Depending on the mass spectrum among χ 1 , χ 2 and X, either X or χ 2 can be naturally long-lived so that the quark pair qq can emerge in the atmosphere with a sizable probability. We shall discuss both of the possibilities in our analysis and demonstrate that they successfully accommodate the ANITA events, yielding the highest event probability at an emergence angle θ em of 25 • − 35 • around which the two reported events are associated.
We additionally provide generic insights in order to understand the angular distribution of the events predicted for ANITA in the general case where a long-lived particle from the dark sector would be produced through scattering of a boosted dark-sector particle in the Earth and decay into hadronic final states in the atmosphere. In particular, we show that the situation where the mean free path of the boosted incoming particle is much larger than the Earth diameter but where its decay product has a mean decay length of dimensions comparable to the Earth radius can lead to an interesting "translucent" case in which the angular distribution can culminate at moderately large emergence angles.
The paper is structured as follows. In Sec. II we start by developing general arguments and demonstrate how the Earth transparency can help predicting events for ANITA with relatively large emergence angles. In Sec. III, we define our benchmark model which may give rise to the signal of interest. Our analysis procedure including relevant Monte Carlo simulation is detailed in Sec. IV. We then report our main findings and results in Sec. V, followed by some discussions and interpretations. Section VI is reserved for conclusions and future prospects.

II. UNDERSTANDING THE ANGULAR DISTRIBUTION OF EAS
The ANITA Collaboration has thus far collected 85.5day data for three different flights, and recorded more than 30 interesting events. Of them, 2 events are identified upward-moving with an energy of ∼ 0.6 EeV, as they do not show a phase reversal characteristic of EAS. Hence, high-energetic upward-propagating cosmic rays may induce these events. Event specifications in terms of reconstructed energy E rec and emergence angle θ em are summarized below, while we refer to [1,2] for more details. EeV and θ em = 35.5 ± 1 • One of the most puzzling features of these anomalous events, besides the propagation of very high-energy particles on large distances, is the angular distribution of such events. While most of the existing interpretations simply rely on the total effective area of the detector in order to compute a total number of events for a certain astrophysical flux, Refs. [9,10,[15][16][17] have attempted to describe the angular repartition of the events predicted by an incoming flux of sterile/right-handed neutrinos. However, the distributions in these studies generically favor relatively small emergence angles ( 10 • ), even though they provide an explanation for the number of events observed by the ANITA Collaboration at angles larger than SM neutrinos would do.
In light of this phenomenological challenge, we discuss what class of scenarios would preferentially give rise to events with a relatively large emergence angle. As detailed in Sec. III, all the ingredients in our benchmark model, the incident particle arriving on the Earth, the particles propagating through the Earth, and the particle inducing an EAS in the low atmosphere, belong to a dark sector. In this section we shall develop some fairly model-independent insight about how the scenario of this sort can give rise to an angular distribution of EAS which is in favor of smaller or larger emergence angles recorded by ANITA. To this end, we simplify the situation at hand, assuming that a particle A scatters off a nucleon and produces a particle B which propagates through the Earth and decays into pairs of quarks: where σ AN and Γ B parameterize the scattering process of A and the decay of B, respectively. As described in more detail in Sec. IV, the effective area of the ANITA detector is a convolution of following different effects: 1. The bigger σ AN is, the more B particles are produced inside the Earth for a given incoming flux of A. Thereafter, the longer-lived the B particle is, the more probable it decays after escaping the surface rather than inside the Earth.
2. Once escaping the Earth, the shorter-lived the B particle is, the more it decays before passing the ANITA detector.
3. Depending on the location of the decay and the opening angle between the ANITA detector and the shower axis, the resulting electric field may initiate detection of an event.
The effective area of the detector per unit emergence angle can then be sketched to be of the form where r ch denotes the spatial direction of the chord of interest and where dΩ E corresponds to the Earth elementary area per unit emergence angle (integrated over the azimuthal angle around the rotation axis of the Earth). The quantity P ex (θ em , E ex ) denotes the probability that a particle B escapes the Earth with an energy of E ex for a given incoming particle A (step 1.). Finally, P dec (θ em , E ex ) and P det ( r ch ) respectively stand for the probability that the particle B decays in the atmosphere (step 2.) and the probability that the shower produced by such a decay is detected by ANITA (step 3.). For simplicity, we assume for now that the detected showers are so boosted that we can consider only the chords crossing the ANITA location and simply look at the variation of these quantities as a function of the emergence angle only. In our scenario, the particles A and B are so feebly interacting that we can neglect any loss of energy which could arise while propagating. Therefore, the exit energy E ex can be uniquely related to the energy of the incoming A particle.
The elementary solid angle contributing to the effective area, integrated over the azimuthal angle, for a given emergence angle (see FIG. 2) is where R E stands for the radius of the Earth. Therefore, neglecting the altitude of the ANITA radio-antenna, we see that a particle -which would ideally propagate freely through the Earth, decay in the low atmosphere with probability being one, and produce a sufficiently energetic EAS -would be detected by ANITA with a total effective area equal roughly to the Earth area, which is a strength of the ANITA detector. Obviously, Eq. (2.3) implies that the angular distribution of the events detected by ANITA would be necessarily suppressed around θ em = 0 • and θ em = 90 • while showing a peak at θ em = 54.7 • . However, in the presence of interactions in the Earth, the interaction probability of the particle A along a chord of length l ch varies as where l A denotes the mean free path of A. This distribution sharply peaks around l l ch in the case of l ch l A , but becomes relatively flat in the case of l ch l A . Therefore, the shorter l A is, the more the scattering processes happen immediately after A enters the Earth surface.
In the former case (i.e., l ch l A ), the efficiency of the scattering process would then be almost independent of the emergence angle. Thereafter, if the particle B produced by such a scattering were long-lived enough, it would escape the Earth with a probability being higher at low emergence angles. This is because the probability that it decays inside the Earth is smaller for small propagation lengths as compared to particles produced deep under the surface. Such a scenario would therefore be disfavored by the fact that ANITA has seen anomalous events at relatively large emergence angles.
By contrast, in the latter case (i.e., l ch l A ), the Earth is almost transparent to the incoming particles so that the probability in Eq. (2.4) can be approximated to   . lA is fixed to 10 5 km for illustration. Shown are the cases with three representative values of the decay length lB, one smaller than, one comparable to, and one larger than the Earth radius by purple, teal, and brown curves, respectively.
This implies that the interaction probability along a chord of l ch is simply l ch /l A . If particle B is long-lived enough, then the exit probability becomes (roughly) proportional to P scat. , i.e., and the relation is valid as far as the decay length of B is larger than the chord length. This shows that the resulting probability distribution is a linearly increasing function in l ch or equivalently in θ em . 3 On the other hand, if particle B is relatively short-lived (i.e., its decay length is shorter than the chord length), the event where A scatters at a depth larger than the B decay length would not allow B to escape the Earth surface since it is likely to decay earlier. The volume of the Earth contributing efficiently to the production of B which may escape is therefore delimited by the surface l = min(l ch , l B ) as depicted in FIG. 3. The exit probability is then scaling like In other words, the proportionality of P ex simply follows l ch for l ch < l B and is frozen to l ch = l B (corresponding to θ em = θ sat ) for l ch > l B . This implies that if the decay product is short-lived as compared to the radius of 3 Note that the effect is accentuated as the mean density along the chord is larger when the chord length increases.

FIG. 5.
Notational conventions for the geometry of the shower opening.
the Earth, the exit probability has to saturate to a constant value, independent of the emergence angle θ em as is schematically shown in FIG. 4 with illustrative values of the decay length l B .
As far as the decay probability is concerned, for a particle B escaping the Earth with an emergence angle θ em , the probability that B decays before reaching an altitude h atm ∼ 15 km is given by Once the shower created in the lower atmosphere, with a maximum angular opening θ sh , we assume that in our analysis ANITA can observe the signal as long as the detector is located inside the cone originating at the point of the decay, as depicted in FIG. 5. Therefore, the capability of detecting a shower, encoded in P det for a given emergence angle, need to be weighted by the surface area of the shower A sh at the detector location which is, in turn, quadratically proportional to the distance between the ANITA detector and the decay point d A : where h A (= 35 km) is the altitude of the ANITA detector above the Earth. However, the reader should note that if showers are emitted too far from the detector, the resultant electric field amplitude may be diminished and as a consequence the detector may not be triggered. A proper treatment for such an effect would require a careful simulation study on the shower electric field, which is beyond the scope of this work. Nevertheless, we expect that it would tend to suppress even more the effective area towards lower angles in favor of large emergence angles, enhancing the effect that we are describing here.  • Opaque case: In the case, similarly to the SM neutrino one, where the Earth is opaque enough to absorb most of the B particles before they can escape the Earth, the exit probability is suppressed at large emergence angles (see also discussions associated with purple curves). Therefore, given the similar suppression of the decaying probability of B exhibited in FIG. 6, the effective area of ANITA will peak at a very small emergence angle.
• Transparent case: In the case where the Earth is transparent, both to the A particle and to the decaying particle B, the exit probability will basically be given by the evolution of the chord length (see also discussions associated with brown curves). As a consequence, it will be maximized at a large angle. The suppression of the decay probability will compete only at large emergence angles, leading to an overall shape of the effective area with a maximum relatively close to the vertical (i.e., θ em = 90 • ).
• emergence angle around 25 • − 35 • , which could potentially explain the value of the emergence angle of the anomalous events observed by the ANITA Collaboration (see also discussions associated with teal curves).
We also visualize this qualitative summary of the three regimes in FIG. 7, not only incorporating the different contributions to the effective area but taking into account the fact that different emergence angles are, in reality, affected by different mean density along the propagation chord ρ ch (we adopted the Preliminary Earth Reference Model of Ref. [21] in our study for the density profile). The latter factor essentially affects the exit probability P ex . We simply re-weight it by ρ ch /ρ max with ρ max being the maximum density inside the Earth as we are here interested in spectral behaviors of A eff . From now on we explore a model of inelastic boosted dark matter, and show that viable regions of parameter space can describe the above-discussed translucent case.

III. BENCHMARK MODEL
We are now in the position to define our benchmark dark-matter model to explain the ANITA anomaly. Once the model definition is established, we discuss production of boosted dark matter followed by its expected flux near the Earth in the first subsection. The next subsection focuses on the scattering of the boosted dark matter in the Earth to create a long-lived particle. Finally, a detailed discussion on the decay of the long-lived particle appears in the last subsection.

A. Benchmark model and production of boosted dark matter
As briefly mentioned in the Introduction, the minimal particle contents include super-heavy dark matter φ, light dark matter χ 1 , dark-sector unstable state χ 2 , and dark gauge boson X. For simplicity, we assume that φ is scalar while χ 1,2 are Dirac-fermionic. The interaction Lagrangian L int contains the following operators: where X µν (F µν ) is the field strength tensor for the dark gauge boson X (the ordinary SM photon) and is the usual kinetic mixing parameter. The term in (3.1) governs the decay of φ to a χ 1 pair and a Yukawa coupling y φ parameterizes the decay strength. On the other hand, the second term in (3.2) takes care of the upscattering process of χ 1 to χ 2 and the decay of χ 2 to χ 1 and X, with g 12 encoding their strength. We emphasize that in the above model the heavier dark matter species φ does not have direct couplings to SM particles, whereas the lighter one χ 1 has direct interactions through a dark gauge boson portal. Therefore, any loop-induced decay modes of φ to SM particles are suppressed enough for φ to predominantly decays to a pair of χ 1 's.
Speaking of dark matter relics, we assume that φ dominates over χ 1 , so the dark-matter halo in our galaxy essentially consists of φ. Be aware that the φ dominance is, in principle, not a necessary condition to explain the ANITA events. The reduction of the φ fraction would decrease the flux of χ 1 linearly, but it can be compensated by lessening the lifetime of φ (as shown in Eq. (3.4)) as far as such a trade-off is consistent with various astrophysical and cosmological observations. We further imagine the situations where the mass of χ 1 , m 1 is of MeV to GeV-range, so it can be thermally produced in the early universe. The present-day relic density of such χ 1 is determined by interrelationships among , m 1 , and the mass of dark gauge boson m X . The relic of χ 1 can be subdominant in a wide range of parameter space as long as m 1 > m X (and in the resonance regime 2m 1 m X ), and therefore our aforementioned assumption is readily satisfied. On the other hand, EeV-scale super-heavy φ under consideration is produced non-thermally [23,24], and it is possible for such φ to take care of almost entire relic abundance of cold dark matter.
Since relic dark matter is mostly given by φ, it is reasonable to assume that the spatial distribution of φ obeys standard dark matter halo density distributions. Our choice in this study is the Navarro-Frenk-White profile [25,26], so the φ energy density ρ φ has the form of where ρ 0 0.3 GeV/cm 3 is the local dark-matter density at r 8.33 kpc and where r s = 24 kpc is the scale radius. The flux of boosted dark matter χ 1 , F 1 , for a given direction is then expressed as follows: where τ φ (or equivalently the inverse of decay rate Γ φ ) is the mean lifetime of φ and where the prefactor 2 takes into account the fact that φ decays to a pair of χ 1 's. Here the incident direction is parameterized by the declination angle θ d and the right-ascension angle φ ra in the International Celestial Reference System (ICRS). We remark that the radial coordinate r in ρ φ is a function of these angles and the line-of-sight s, and so is the dark-matter energy density ρ φ .

B. Scattering-off of boosted dark matter
Once such a boosted χ 1 reaches the Earth, it may scatter off a nucleon via a t-channel exchange of dark gauge boson X and a dark-sector unstable state χ 2 comes out. The nucleon usually breaks up due to a large energy transfer, so the process is essentially initiated with deep inelastic scattering (DIS) between χ 1 and nucleon N . The DIS of this sort is extensively explored in Ref. [27], and we simply quote their final result for σ DIS χ1N : 5) where m N denotes the nucleon mass of interest, α is the fine structure constant, and f (x) describes the associated parton distribution function (PDF). The relations among a dimensionful quantity Q 2 and two dimensionless quan-tities x and y are given by where p 1 and p 2 denote four momenta of the incoming χ 1 and the outgoing χ 2 , respectively. x implies the momentum fraction carried by a parton while y is fractional energy loss. For a given E 1 one can calculate numerically the χ 1nucleon DIS cross sections, convolving relevant PDFs. In our analysis we model the total DIS cross section for a convenient parameter survey, after calculating σ DIS χ1N of several representative parameter points with MSTW2008NNLO [28]. We find that the following model describes functional behaviors of σ DIS χ1N in 0.1 EeV E 1 10 EeV and m X 10 GeV fairly well: where K(x) has the form of (3.10) Note that the above empirical model does not depend on m 1 and m 2 . This is because our numerical study suggests that the variation in σ DIS χ1N be of order (at most) 2 − 3% with m 1,2 5 GeV. As we shall show in Sec. V, some parameter choices belonging to these ranges look highly plausible in generating ANITA-like events. So, we slenderize our parameter survey, relying on the abovegiven model and focusing on the associated parameter regions.
Given the cross section, one can estimate the mean free path of χ 1 . Kinetic theory suggests that the average travel distance of χ 1 (henceforth denoted by L 1 ) should be where n is the mean nucleon number density along the χ 1 propagation chord. Since it depends on Earth layers that χ 1 traverses, we explicitly express its dependence on the emergence angle θ em . Certainly, the functional behavior of n in θ em is highly nontrivial, as the associated density profile inside the Earth shows a drastic change especially at the boundary between adjacent interior layers. Just to build up our intuition, we find that  for θ em = 30 • at which the chord length is the same as the Earth radius.

C. Decay of a long-lived particle
Once χ 2 is created, it decays back to χ 1 together with hadrons (through a quark pair) which invoke EAS. As mentioned in the Introduction, the χ 2 decay takes place via either on-shell or off-shell dark gauge boson, depending on the mass relation among χ 1 , χ 2 and X. Labeling them by "on-shell" scenario and "off-shell" scenario respectively, we discuss their characteristic features separately.
On-shell scenario: In this scenario m 2 is greater than m 1 + m X so that an on-shell X comes out as a decay product of χ 2 . Unless the coupling g 12 and/or a mass gap m 2 −m 1 −m X is unusually small, the χ 2 decay is prompt (i.e., it occurs at the scattering point), and X may be long-lived and propagating. The total decay width of X, Γ X is given by where C i denotes the color factor for particle species i which runs over all charged SM fermions lighter than half the mass of X. 4 As a rough estimate we see that the laboratory-frame mean decay length of X, X,lab is X,lab ∼ 400 km (3.14) with E X being the amount of energy that X carries.
One possible issue with this scenario is disappearance of X via X + N → anything, before it decays. To develop some intuition on such a possibility, its cross section should be first estimated. If X is energetic enough to be treated nearly massless, we can get some hint from the photon-nucleon DIS cross section. From the experimental measurements of σ DIS γp at HERA [29,30]  , the center-of-mass energy of the X-p system, is given by 2E X m p . 5 One may estimate the typical distance for a single scattering process using the relation in (3.11), and see that the distance is much larger than the mean decay length above for a given set of E X and . In other words, dark gauge boson is strongly inclined to decay before disappearing.
Off-shell scenario: Contrary to the previous scenario, the mass of dark gauge boson X is larger than the mass difference between χ 2 and χ 1 . Therefore, if the mass gap δm (≡ m 2 −m 1 ) (kinematically) allows for opening decay modes to SM fermion pairs, the associated decay process occurs via an off-shell X. χ 2 is now long-lived and the associated decay width Γ 2 is formulated [31] by the form of under the assumption of m i δm m 2 m X (see Appendix B in Ref. [31] for the exact formula). Our rough estimate shows that the laboratory-frame mean decay length of χ 2 , 2,lab is 2,lab ∼ 12, 500 km where E 2 is the energy of the scattered χ 2 .
We close this section, rendering a few comments. It is summarized that the event topology under consideration begins with an upscattering of an invisible particle and accompanies a long-lived intermediary state decaying to hadrons. Indeed, the structure of upscatterdecay has been adopted in a diverse range of astrophysical phenomena [32][33][34][35][36], search strategies in particle accelerator experiments [37][38][39][40][41][42][43][44], and inelastic boosted dark matter [27,31,45,46]. In particular, Refs. [47][48][49] have considered the upscattering of dark matter to an excited state somewhere outside a detector, followed by its deexcitation in the detector along with visible particles. In the context of the ANITA anomaly, Refs. [13,17] have taken related strategies.

IV. ANALYSIS PROCEDURE
In this section, we detail the procedure that we used in order to compute the differential number of events per emergence angle at which the ANITA detector would observe events according to our model prediction, redefining notation and conventions more carefully if needed. Generically, the expected number of signal events N sig is given by the exposure time t exp times the inner product of flux and effective area vectors. From Eq. (3.4) we see that our signal flux varies with the angles θ d and φ ra (see also FIG. 8 for the notations). 6 Labeling the "effective" area projected onto the flux vector by A eff , we remark that A eff potentially depends on θ d and φ ra . Therefore, 6 Note, however, that the flux producing upward-going showers at the south pole corresponds to the decay of DM particle decaying in the outskirt of the galaxy. Therefore the angular distribution of the incoming flux is essentially isotropic in our case. The impact of the choice of DM profile is therefore negligible in our work.
N sig is obtained after integrating over the incoming flux direction θ d and φ ra : where dΩ ICRS = d(cos θ d )dφ ra . Note that in our model the energy of the incoming boosted dark-matter state, E 1 , is uniquely determined by the mass of the heavy darkmatter state φ. Therefore, the flux of χ 1 particles, F 1 , is implicitly depending on E 1 . Likewise, (as we will see later) A eff varies with E 1 , so N sig is essentially a function of the energy E 1 . For brevity, we will omit the E 1 dependence in what follows, unless it is necessary for a better understanding. We further note that for a given energy E 1 the effective area differs according to the emergence angle θ em with which the decaying particle (either X or χ 2 ) escapes the surface of the Earth and according to its energy E ex when it exits the Earth crust. Therefore, for a given direction of the flux arrival (θ d , φ ra ), the effective area A eff (θ d , φ ra ) can be formulated as Since one of the main goals in our work is to obtain an appropriate angular distribution of the signal as seen by ANITA, we seek to obtain the differential number of signal events in θ em , i.e., (4.4) The key quantity in order to compute this angular distribution is the elementary effective area A eff in Eq. (4.3) per emergence angle θ em and exit energy E ex for a given direction d Ω ICRS and energy E 1 of the incoming flux of χ 1 particles. This area essentially encodes all the dynamical details of the detection, from the χ 1 entrance in the Earth, its propagation and conversion into a long-lived state, to the possible decay of the latter in the atmosphere and the EAS detection by the radio-antennas of ANITA. Inspired by the analysis scheme in Ref. [4], we modulize the whole calculational procedure in the following fashion: for a given incoming energy and direction of the flux, we (i) integrate over the entire Earth surface, scanning over all the surface points of impact where a flux of χ 1 particles arriving along parallel chords strikes; (ii) compute the probability that the χ 1 propagates through the Earth, scatters and creates a long-lived particle (either X or χ 2 ) which escapes the Earth surface with a given energy E ex ; (iii) compute the probability that the long-lived particle decays in the atmosphere at a relevant height, creating an EAS; (iv) calculate the probability that such shower is detected by ANITA.
Such decomposition can be formulated as where n E is the Earth-surface unit normal vector while n 1 is a unit vector lying in the flux direction of χ 1 (see also FIG. 8). The solid angle element with respect to the Earth center is denoted by dΩ E = d(cos θ E )dφ E . In the third row we explicitly express the dependence of the differential exit probability dP ex /dE ex on E 1 as well. In the fourth row the decay probability of the long-lived particle is integrated along its propagation chord. Finally, in the last row P det describes the capability for a given EAS to yield a sufficient electric field along a direction which would reach the payload hence result in detection.
We mark its dependence on the shower opening angle θ sh and the energy of EAS E EAS . The ensuing subsections are devoted to discussing in more details the calculation of the steps (ii) through (iv).

A. Exit probability
Throughout this work, we perform a Monte Carlo simulation in order to estimate accurately the exit probability P ex appearing in the third row of Eq. (4.5), as a function of the emergence angle θ em . To this end, we adapted the public code provided with Ref. [50] for our BSM scenario. Assuming the Earth surface to be covered by a 4 km-thick ice layer, the simulation accounts for the varying density of the Earth crust along the propagation chord.
The quantity dP ex (E 1 ; θ d , φ ra ; θ em , θ E , φ E )/dE ex is defined as the probability that an incoming boosted darkmatter particle χ 1 -which arrives with an incidence angle θ i (corresponding to an emergence angle θ em = 90 • − θ i ) and an energy E 1 -exits the Earth with an energy E ex . For every emergence angle θ em , each event of the Monte Carlo simulation therefore tracks along the propagation chord the different possible interactions, according to the two aforementioned benchmark scenarios.
On-shell scenario: 1. An incoming χ 1 particle scatters off a nucleon with the cross-section appearing in Eq. (3.9) and converts into a particle χ 2 which promptly decays into a dark photon X and a boosted particle χ 1 , 2. the dark photon X propagates until it reaches the Earth surface, as long as it does not decay inside the Earth, and 3. the (secondary) boosted dark-matter particle χ 1 produced by the χ 2 decay propagates until it rescatters inside the Earth or escapes the surface.
Off-shell scenario: 1. An incoming χ 1 particle scatters off a nucleon with the cross-section shown in Eq. (3.9) and converts into particle χ 2 , 2. the dark-sector particle χ 2 propagates until it decays inside the Earth or escapes the surface, and 3. in the case that χ 2 decays inside the Earth, its decay product (secondary) χ 1 propagates until it rescatters inside the Earth or escapes the surface.
Pursuing the propagation of secondary χ 1 's created in the Earth by the χ 2 decay allows us to keep track of possible regeneration which could potentially increase the exit probability at large emergence angles and low exit energies. However, our simulation study suggests that the regeneration play an insignificant role in the scenarios that we investigate in this work, as we have elaborated in Sec. II. We therefore focus on cases where the Earth is relatively transparent to the flux of incoming particles.

B. Decay probability in the atmosphere
Once a long-lived particle species i (= X or χ 2 ) comes out of the Earth surface, its decay simply obeys the exponential decay law where i,lab is the laboratory-frame mean decay length of the long-lived particle i. In the two cases where i = X and i = χ 2 , the results appearing in Eqs. (3.14) and (3.17) can be substituted to i,lab with E X and E 2 replaced by E ex . 7 As shown in Ref. [15], in the case where an astrophysical flux would produce tau leptons inside the Earth, the exit probability of the tau leptons would spread over a wide range of exit energies. Indeed, taus quickly lose energy while propagating in the Earth. Moreover, during regeneration processes, the particle produced through a decay, which may then re-scatter off a nucleon, will lead to an event with even lower exit energy. In our case, this sort of energy loss of the particles propagating through the Earth can be negligible since their interactions with the visible sector particles are very feeble. We therefore take into consideration only the energy losses (with respect to incoming energy E 1 ), which are caused by the fractional energy release in the scattering and decay processes.
The value of the exit energy has a crucial impact on the decay probability since the more boosted particles are, the less they are expected to decay before reaching the ANITA detector. Our simulation carefully takes care of possible energy diminishment down to E ex = 100 TeV, although most of the events involve a much larger exit energy in our proposed model due to the reasons discussed above.

C. Detection probability
The hadronic decay products can create an EAS in the atmosphere which can eventually be detected by the radio-antennas of ANITA at an altitude of h A ∼ 35 km. Experimentally, it is not sufficient that the shower reaches the payload but the local electric field measured by the ANITA detector should be large enough to overcome the associated threshold. Here the probability P det 7 Be aware again that Eq. (3.17) is an approximation valid under restricted mass spectra. Our simulation study was based on the exact formula in Ref. [31].
parameterizes the likelihood of recording an event with such effects considered.
The authors of Ref. [4] have performed a dedicated analysis regarding this probability, drawing a few important observations: (a) it gets gradually challenging for a shower created at altitudes above O(10 km) to induce a fully-developed EAS, (b) the electric field peak, depending on the altitude where the shower is produced, lies at an opening angle θ sh in-between ∼ 1 • and ∼ 2 • , (c) the lower E ex is, the less the ANITA detector is triggered, and (d) the further from the ANITA detector the EAS is produced, the less likely it is to be detected.
Unfortunately, although Ref. [4] provides an approximate formula in order to evaluate the electric field produced by an EAS at the ANITA location, with an emergence angle of 30 • and arbitrary distance from the point of decay and shower angle, the exact functional dependence of the parameters involved in this expression with the emergence angle has not been exhibited. This renders the evaluation of the precise probability P det challenging from a theoretical point of view. For simplicity, we rather use the main conclusions of Ref. [4], following the assumptions made in Ref. [15], that is, we assume that P det = 1 as long as the shower is produced at an altitude lower than 15 km, and whenever the detector is contained in a cone of angular opening θ sh < 1.5 • with respect to the shower axis, and P det = 0 otherwise.
Finally, the E EAS dependence of P det must be taken with care, as the energy of the EAS could potentially be related to E ex in a non-trivial manner. In the "onshell" scenario, the whole energy of the escaping X is transferred to decay products as it decays fully visibly. Therefore, we are allowed to set E EAS = E ex . On the contrary, in the "off-shell" scenario, a certain fraction of E ex is taken away by the outgoing χ 1 . We find that if χ 1,2 are sufficiently heavier than hadronic decay products, χ 1 is produced nearly at rest in the χ 2 rest frame. Therefore, in the laboratory frame, χ 1 takes away ∼ (m 1 /m 2 )E ex , and our simulation study with PYTHIA 8 [51] confirms that a large majority of events conforms to this estimate. In our analysis with off-shell scenario, we survey parameter space satisfying such mass relations, and thus adopt the approximation that E EAS ≈ (δm/m 2 )E ex .

V. RESULTS AND INTERPRETATIONS
We now report our simulation results in this section, followed by discussing their phenomenological implications. In order to demonstrate that the inelastic boosted dark matter scenario can accommodate the ANITA anomalous events, we choose the following two benchmark sets of parameters for the on-shell and off-On-shell X Off-shell X Note that in the on-shell case the masses of χ 1 and χ 2 do not play any role either in the simulation of the effective area or in the decay width of the dark gauge boson, as long as they satisfy the mass relation with X, they are light enough to be considered relativistic, and χ 2 decays instantly. By contrast, in the off-shell scenario, all the mass values affect the propagation and the decay, although m 1 does not play a role in the propagation code. The m φ values are selected in such a way that the associated hadronic decay products collectively carry away an energy within 0.5 − 1 EeV. In FIG. 9 we exhibit the differential exit probability in the emergence angle θ em , integrating over all possible exit energy E ex . The left and right panels correspond to the on-shell and off-shell scenarios, respectively, for which we vary the value of the kinetic mixing parameter with the value of the dark-sector coupling g 12 fixed to unity. We see that the red curves in both scenarios show a different behavior from the others. As we elaborated in Sec. II, this is because the mean decay length of X (for the on-shell scenario in the left panel) or χ 2 (for the off-shell scenario in the right panel) resulting from the associated parameter choices is much shorter than most of the chord lengths so that the exit probability quickly saturates in increasing θ em .
We remark that the variation in g 12 would alter the results in FIG. 9 in different manners in the on-shell and off-shell scenarios. In the former case, raising g 12 would result in enhancing the scattering cross section, while the decay length of X remains unaffected. As far as the mean free path of χ 1 keeps much larger than the dimension of the Earth, this would simply lead the effect of rescaling the total exit probability by g 2 12 . On the contrary, in the latter case, both the scattering cross section and the decay width of χ 2 scale as g 2 12 2 . Therefore, increasing or decreasing the dark-sector coupling is equivalent to inversely varying the kinetic mixing parameter by the same magnitude.
Following the procedure described in Sec. IV, we next integrate over the incoming χ 1 directions defined by the (θ d , φ ra ) pair and the Earth surface dΩ E , in order to find the effective area of the detector. In both scenarios we adapt the lifetimes of the decaying particles such that the total expected number of events is of O(1), given the exposure time t exp = 85.5 days. As can be seen in both the on-shell and off-shell cases, varying 20 40 60 80 On-shell scenario panels) and differential number of events (right panels) as a function of the emergence angle θem anticipated for ANITA-85.5 days, in the cases of "on-shell" scenario (upper panels) and "off-shell" scenario (lower panels). Parameter choices appear in the respective legends.
the value of the kinetic mixing parameter clearly exhibits the transition from the opaque case (large ) to the transparent case (small ). Furthermore, we observe that the maximum value of the angular distribution changes according to what was described in Sec. II. As expected previously, the intermediate case (or translucent case) shows a maximum at emergence angles of O(25 • − 35 • ), which is the range in which the two anomalous events observed by ANITA are lying.
Finally, we make a brief comment on the relic abundance of χ 1 . Since our benchmark dark-sector model assumes that χ 1 is a negligible relic component, it is crucial to check whether or not our parameter choices indeed advocate the assumption. We have explicitly calculated the χ 1 relic density for all sets of parameter choices, employing MicrOMEGAs [52]. The numerical outcomes suggest that m 1 > ∼ 0.45 GeV in the on-shell scenario and m 1 > ∼ 1.9 GeV in the off-shell scenario be good enough to have negligible ( < ∼ 0.1%) abundance for our coupling choices, in addition to the resonance regime satisfying m X 2m 1 .

VI. CONCLUSIONS AND OUTLOOK
The upward-moving anomalous events reported recently by the ANITA Collaboration [1,2] have brought increasing attention in the particle physics community, as it seems disfavored to explain them under the SM framework, e.g., tau-neutrino-induced processes. Therefore, the phenomenon has been considered as a hint to new physics beyond the SM, stimulating a growing number of phenomenological studies in the context of new physics scenarios [3,[9][10][11][12][13][14][15][16][17][18]. In this paper, we conducted a dedicated study on interpreting the ANITA anomaly with scenarios of inelastic boosted dark matter [45].
The proposed model in this paper includes an EeVrange super-heavy (scalar) dark matter φ which predominantly decays to a pair of MeV − GeV-range (fermionic) dark matter particles χ 1 in the galactic halo in the present universe. Due to the large mass difference between the two dark-matter states, the lighter dark-matter particles χ 1 produced by the decay of the super-heavy candidate φ are highly boosted and reach the Earth. An incident χ 1 may therefore upscatter to a heavier, unstable dark-sector state χ 2 via an exchange of dark gauge boson X. χ 2 then decays to χ 1 and hadronic final states through on-shell (off-shell) intermediary X, while X (χ 2 ) often becomes considerably long-lived. Such a long-lived particle may escape from the Earth, carrying an EeVscale energy. Once it decays in the atmosphere at a proper altitude, there arises an EAS which can be detected by ANITA.
We then conducted a parameter survey according to our analysis strategy elaborated in Sec. IV. We examined opaque, translucent, and transparent cases (see Sec. II for their definitions) with different parameter choices, and found that the translucent case can accommodate the ANITA events best in the sense that the expected number of signal events is maximized at an emergence angle in-between 25 • and 35 • . The selected parameter values are consistent with existing constraints. More importantly, the dark-matter relic abundance associated with our parameter choices supports the benchmark darksector model assumption that φ and χ 1 are the dominant and negligible relic components, respectively. Therefore, the anomalous events observed by the ANITA Collaboration would be considered as a sign of the inelastic boosted dark matter scenario adopted in this study.
We would like to emphasize that in our scenario, the very feeble interaction of the particles propagating through the Earth with nuclei renders detectors such as the Auger observatory or IceCube far less competitive than ANITA in terms of signal detection, given their small fiducial volume as compared to the volume of the South Pole atmosphere that ANITA utilizes to observe EAS. Indeed, as opposed to scenarios which require to produce boosted tau leptons inside the Earth crust in order to produce EAS in the atmosphere, our scenario simply requires that boosted X or χ 2 particles reach the Earth surface before decaying, hence very unlikely to leave any significant traces in Earth-based detector.
As a final remark, we point out that our analysis strategy is straightforwardly applicable to other dark-sector scenarios. Obviously, a vector mediator (dark gauge boson in this work) can be replaced by a scalar, e.g., dark Higgs portal models where the hadronic interactions are larger. Also, although our approach is predicated upon a non-minimal boosted dark matter model containing a dark-sector unstable state (χ 2 in this study), scenarios similar to our "on-shell" scenario are available in the minimal boosted dark matter model, by allowing for the following term in the Lagrangian: where g 11 is the coupling constant governing a χ 1 -flavorconserving interaction. The existence of this operator, in turn, allows for the process by radiating a dark gauge boson X off initial-/final-state χ 1 (called "dark-strahlung") [53]. This possibility certainly looks plausible, and we leave the exploration in this direction for a future work.