Connecting ANITA Anomalous Events to a Non-thermal Dark Matter Scenario

The ANtarctic Impulsive Transient Antenna (ANITA) collaboration has observed two EeV-energy, upward going events originating from below the horizon. As no standard model (SM) particles, propagating through the earth at such energy and exit angles, can give rise to the required survival probability for the observed events, several beyond standard model (BSM) proposals have come up. We propose a scenario where a $Z_2$ odd sector is responsible for such anomalous events. The next to lightest $Z_2$ odd particle or the next to lightest stable particle (NLSP), created from ultra high energy neutrino interactions with nucleons, can pass through the earth and then decay into the lightest $Z_2$ odd particle or the lightest stable particle (LSP) and a tau lepton. The long-lived nature of the NLSP requires its coupling with the LSP to be very small, keeping the LSP out of thermal equilibrium in the early universe. The LSP can then be a non-thermal dark matter while the tau leptons produced from NLSP decay after passing through the earth can explain the ANITA events. We first show that a purely non-thermal dark matter scenario can not give rise to the required event rates while a hybrid scenario where dark matter can have both thermal as a well non-thermal contribution to its relic abundance, serves the purpose.

Here we propose a specific type of dark matter scenario where DM relic abundance is generated partially or fully from a non-thermal mechanism. DM contributes around 27% to the present universe's energy density according to the latest cosmology data provided by the Planck satellite [26]. In terms of density parameter Ω DM and h " Hubble Parameter{p100 km s´1Mpc´1q, the present DM abundance is conventionally reported as [26]: Ω DM h 2 " 0.120˘0.001 at 68% CL. Since none of the SM particles could provide a suitable DM candidate, several BSM proposals have come up out of which the weakly interacting massive particle (WIMP) paradigm has been the most popular one. In this framework, a dark matter candidate typically with electroweak scale mass and interaction rate similar to electroweak interactions can give rise to the correct dark matter relic abundance, a remarkable coincidence often referred to as the WIMP Miracle. However, none of the experiments looking for direct detection of DM has found any signal yet, motivating the community to seek a paradigm shift. One such scenario that has drawn attention in the last few years is non-thermal dark matter [27]. In such a framework, DM particles have so feeble interactions with the remaining thermal bath that it never attains thermal equilibrium at any epoch in the early universe. However, it can be produced from decay of some heavy particles or scattering processes, popularly known as the freeze-in mechanism [27][28][29][30], leading to a new paradigm called freeze-in (or feebly interacting) massive particle (FIMP). For a recent review of this DM paradigm, please see [31]. Typical FIMP or non-thermal DM models involve tiny couplings of DM with SM or other particles in the thermal bath and often lead to long lived particles. Such long lived charged particles have the potential to pass through the earth before decaying into DM and a tau lepton which could explain the anomalous events observed by ANITA [15,16].
Motivated by these, we propose a model which is an extension of the scotogenic model [32] that provides a common origin of DM and light neutrino masses. The original scotogenic model is an extension of the SM by three right handed neutrinos and one additional scalar doublet all of which are odd under an in-built and unbroken Z 2 symmetry. It is worth mentioning at this point that the observation of non-zero neutrino masses and large leptonic mixing [33] has also been another motivation for BSM physics for last few decades. While the addition of singlet right handed neutrinos to the SM content can give rise to the usual seesaw mechanism [34][35][36] for neutrino mass at tree level, the scotogenic framework can explain the origin of neutrino mass and dark matter in a unified manner. The Z 2 odd particles, the lightest of which is the DM candidate, give rise to light neutrino masses at one loop level in this model. We extend this model suitably to explain the AAEs. We first consider an extension of this model to incorporate a purely non-thermal or FIMP type DM and name it as Model-I. Here, the DM is a gauge singlet right handed neutrino whose relic is generated purely from the non-thermal contribution, by virtue of it's small couplings to SM leptons. Although correct DM properties as well as light neutrino masses can be realized in such a scenario, it fails to explain the AAEs. We then consider a hybrid setup where DM relic receives both thermal as well as non-thermal contributions. In this scenario, referred to as Model-II hereafter, DM is the lightest neutral component of the Z 2 odd scalar doublet. Due to electroweak gauge interactions, the DM in this model can not be of purely non-thermal origin. However, it can receive non-thermal contribution from the decay of heavier particles after it freezes out, similar to the recent works [37,38]. In such scenarios, the thermally under-abundant DM can satisfy the correct relic abundance criteria due to the late non-thermal contributions. We show how the correct DM and neutrino phenomenology can be obtained in this hybrid model along with the explanation for the AAEs. This paper is organised as follows. In section 2.1 and 2.2 we discuss the two models along with the corresponding DM phenomenology. In section 3 we discuss the origin of ANITA anomalous events in both the models. We briefly comment upon different ways to probe our models in section 4 and finally conclude in section 5.

Model-I: Pure Freeze-in
As mentioned earlier, we first discuss a model where DM relic is generated purely from the freeze-in mechanism. The particle content of Model-I are shown in table 1. There are five different types of BSM fields all of which are odd under an in-built and unbroken Z 2 symmetry. We need three copies of singlet right handed neutrino N in order to account for non-thermal DM and light neutrino masses simultaneously. The lightest of the right handed neutrino N 1 is the DM candidate having tiny Yukawa couplings with leptons while the heavier right handed neutrinos generate the light neutrino masses at radiative level. While Φ 2 , N are same as those in the minimal scotogenic model and are sufficient to realize N 1 as freeze-in DM and generate light neutrino masses and mixing, the other three types of particles E L,R , ψ 1 , ψ 2 are required in order to explain the AAEs. It is worth mentioning that due to very small Yukawa couplings associated with one of the right handed neutrinos N 1 , the lightest neutrino mass becomes vanishingly small, predicting a hierarchical light neutrino mass spectrum.
The scalar potential is, As mentioned, the lightest right handed neutrino N 1 is the dark matter candidate in this model and being gauge singlet, it can be prevented from being produced thermally in the early universe, if the corresponding Yukawa couplings are very small. We consider such a possibility where N 1 is produced non-thermally from the decay of Φ 2 . The decay of Φ 2 can occur either in equilibrium or after its freezeout, the latter is similar to the superWIMP formalism [39]. It should be noted that all the components of Φ 2 can contribute to the non-thermal production of N 1 namely, Φ2 Ñ N 1 l˘, Φ 0 2 Ñ N 1 ν, l " e, µ, τ . Since the Z 2 symmetry remains unbroken, the neutral components of Φ 2 do not acquire any non-zero vacuum expectation value (VEV) which can be ensured by choosing µ 2 2 ą 0. We parametrize the two scalar doublets Φ 1,2 as where Φ 1 is identified with the SM Higgs doublet whose neutral component acquires a VEV denoted by v, responsible for electroweak symmetry breaking (EWSB). As follows from the scalar potential, after EWSB, the physical scalars (originating from the two scalar doublets) have the following masses, The relevant parameters as well as physical masses of different fields are listed in table 2. The other parameters which do not find mention in this table can be chosen as per phenomenological requirements and we do not list them in this table as they do not affect our numerical calculations.
In order to compute the abundance of N 1 , we first write down two Boltzmann equations corresponding to the evolution of Φ2 as well as N 1 . They can be written in terms of their comoving number densities pY " n{s, n " number density, s " entropy densityq as, with z " M Φ2 {T and thermally averaged decay width of Φ2 where K i 's are the modified Bessel functions. Now, for the most part of the parameter space of the Yukawa Y 3 the freeze-in process happens while Φ2 is still in equilibrium. So, we can safely integrate the equation (6). So, the solution of the equation (6) is given as Table 2: Numerical values of different relevant parameters used for Model-I.
and the relic is given as For a chosen benchmark point, the temperature evolution of comoving number densities for DM and its mother particle are shown as function of temperature in figure 1. We choose the benchmark point as M N 1 " 30 MeV, M Φ2 " 1 TeV, Y 3 " 2.92ˆ10´1 0 and find that the correct relic density is satisfied. It can be clearly seen from figure 1 that the number density of DM that is N 1 increases as temperature cools down due to decay of different components of Φ 2 . The equilibrium contribution of Φ2 decay dominates over its post freeze-out contribution to the abundance of N 1 for the chosen benchmark. While N 1 receives non-thermal contribution from all the components of Φ 2 , the coannihilations among different components of Φ 2 have been incorporated following Ref. [40]. While the plot shown in figure  1 corresponds to a particular benchmark point of DM mass and Yukawa coupling, we show the allowed parameter space in terms of these two parameters in figure 2 which gives rise to correct DM relic. While Y 3 , M N 1 are being varied to generate this plot, the other parameters are kept fixed at values mentioned in table 2. The observed correlation between Y 3 and M N 1 from this plot agrees with the approximate analytical expression for DM relic density derived above. Since freeze-in abundance increases with increase in DM coupling [27], the corresponding DM mass needs to be smaller in order to be in agreement with observed abundance. It should also be noted that the requirement of such small Yukawa coupling of DM N 1 for its non-thermal nature effectively decouples one of the singlet neutrinos from the light neutrino mass generation mechanism at one loop level. The other two heavy singlet neutrinos N 2,3 can however, have order one couplings with leptons and can generate the required light neutrino mass. Therefore, the model predicts vanishingly small lightest neutrino mass. Although DM relic is satisfied for this model, but it can not explain the AAEs as we discuss in the upcoming section. Before that, we outline Model-II which satisfies all relevant DM constraints apart from explaining the AAEs.

Model-II: Freeze-out + Freeze-in
In this subsection, we discuss our Model-II where DM relic is generated from a hybrid setup consisting of both freeze-out and freeze-in contributions. The particle content of this model is shown in table 3. There exist six different types of fields, apart from the SM fields, all of which are odd under an in-built Z 2 symmetry. The SM fields, on the other hand, are even under Z 2 symmetry. In order to generate correct neutrino mass, we require at least two copies of right handed singlet neutrinos N .
While the minimal scotogenic model has only two different types of Z 2 odd particles, the roles of other four types of Z 2 odd particles in this model will become clear when we discuss the explanation of the AAEs in upcoming section. Table 3: Field content and transformation properties under the symmetry of the model.
The relevant part of the Yukawa Lagrangian is given by, The scalar potential is, We assume the neutral component Φ 0 3 of Φ 3 to be the lightest Z 2 odd particle and hence the DM candidate while E is the NLSP. Due to the long lived singlet fermion E which decays at late times to dark matter and a tau lepton, we can have a non thermal contribution to the freeze-out abundance of DM. The physical masses of the components of scalar doublets Φ 2 , Φ 3 after EWSB can be found in a way similar to Model-I. However, due to the presence of these two Z 2 odd scalar doublets, there can be mixing between them as well. For simplicity, we ignore such mixing, which amounts to tuning the respective bilinear and quartic couplings like µ 23 , λ 7 , λ 8 , λ 9 in the scalar potential (11).
The coupled Boltzmann equations can be written as, Here, N Φ 0 3 is the average number of DM particles produced from a single decay of E which is one in this case, Γ E is the decay width of E, and H is the Hubble parameter. Assuming that the singlet fermion E do not contribute dominantly to the total energy budget, we can take the comoving entropy density (g ‹s ) and the comoving energy density (g ‹ ) to be approximately constant. Further, we assume that almost all of E decays during the radiation dominated epoch, which is satisfied for the benchmark values of Yukawa couplings we use here. Also, if we assume E decays only after its thermal freeze-out, the first term on the right hand side of equation (13) can be ignored 1 and one can analytically solve the Boltzmann equation for n E above. Writing the above equations in terms of with s " 2π 2 45 g˚sT 3 being the entropy density and changing the variable from time t to z " to the epoch after the freeze-out of E from where we start integrating the differential equation for DM) however, depends on the initial abundance of E and can be found by calculating its freeze-out abundance. Using this solution for Y E , the equation for Y Φ 0 3 can be rewritten as, The decay width of fermion E (assuming , ignoring the mass of leptons is given by, where M E is the mass of fermion E and Y 1 2 , its coupling with the tau lepton and DM. The lifetime requirement of E from ANITA anomaly point of view is [15], For M E " 5 TeV, the required lifetime of 100 ns, can be obtained for Y 1 2 « 10´1 0´1 0´9, typical coupling for non-thermal dark matter. It should be noted that E is produced by interactions of high energy neutrinos with nucleons in a process mediated by the charged component of Φ 2 . The production cross section of E in this process can be large by suitable tuning of the corresponding Yukawa couplings Y 2 , Y 4 , Y 5 , Y 2 2 . By keeping Φ 2 heavier than E, in order to forbid the latter's decay into the former, we can have long lived E which decays after traveling through earth into dark matter and tau lepton. We will discuss the details of AAEs in this model in the upcoming section.
Parameters First benchmark Second benchmark µ 2 400 GeV 300 GeV  In order to show the DM relic density results, we choose two benchmark points, the details of which are given in table 4. For the chosen benchmark points, we then solve the coupled Boltzmann equations for Φ 0 3 , E and show their comoving number densities as function of temperature in figure 3. As can be seen from the plots, the comoving number densities of both Φ 0 3 and E decreases as temperature falls below their respective masses, owing to the Boltzmann suppression for non-relativistic species. At some epochs, the individual rate of interactions of E as well as Φ 0 3 with the SM particles freezes out giving rise to the plateau shaped regions. However, since E is unstable and decays into Φ 0 3 , the comoving number density of E drops very fast at some epoch due to its decay with a corresponding rise in the number density of Φ 0 3 . In the absence of E decay, the freeze-out abundance of Φ 0 3 remains under-abundant compared to the Planck bound. The mass of Φ 0 3 in both the benchmark points given in table 4 falls in the mass regime 80 GeV ď M Φ 0 3 ď 500 GeV where scalar doublet DM fails to give rise to correct thermal relic as is well known from earlier studies on Z 2 odd scalar doublet dark matter. The role of non-thermal contribution in filling this deficit was discussed in earlier work [37].

Explanation of ANITA Events
The estimation of number of AAEs can be facilitated from the evaluation of survival probability which basically represents the fraction of incoming ultra-high energy neutrino flux that ultimately gives rise to τ -leptons near the earth surface. The total number of events can be estimated as, where the effective area A eff « 4 km 2 , T is the exposure time which we take « 25 days by combining the exposures of ANITA-I and ANITA-III, ∆Ω is the acceptance angle which corresponding to isotropic and anisotropic neutrino flux source is given as, ∆Ω « # 2π sr for isotropic case, 2πp1´cos δ θ q « 0.0022 sr for anisotropic case, where δ θ « 1.5˝is the angular uncertainty relative to parent neutrino direction. Note that we have not considered the exposure time for ANITA-II as it was not sensitive to such events. The range of integration in equation (20) should be such that the correct range of shower energy can be obtained.
As we consider τ lepton as the origin of observed events, with E τ " E ν {4 and observed shower energy 0.2-1 EeV, we take E i " 0. 8 EeV, E f " 4 EeV. As far as the flux ΦpE ν q is concerned, for the isotropic case if one assumes the source of the EeV neutrinos to be the Greisen-Zatspin-Kuzmin (GZK) mechanism then in that case, following [41], for the concerned shower energy, Φ iso « 10´2 5 (GeV cm 2 s sr)´1. Otherwise, if some localized source is the origin of the EeV neutrinos, then the upper limit on such anisotropic flux would be Φ aniso « 3.2ˆ10´2 0 (GeV cm 2 s sr)´1 [42][43][44][45][46]. Thus by considering the mentioned parameters in the isotropic and anisotropic cases to get two AAEs, one should have, N « # 2.0ˆ10 2 for isotropic case, 2.0ˆ10 4 for anisotropic case.
This implies that for the observed two AAEs in the isotropic case the should be " 10´2, and for the anisotropic case should be " 10´4. By considering SM interactions, authors of [15] predicted the survival probability SM " 10´7, giving N " 2ˆ10´5 for the isotropic flux. Clearly in the SM explanation of these two events is very unlikely. While using the SM interaction and considering the anisotropic flux one gets N " 2ˆ10´3 which again makes it very implausible to get explanation of the observed events in the SM. Thus in any BSM scenario trying to explain the two AAEs, one should strive to increase the pE ν q and we will show that this can be achieved in one of the SM extensions that we are considering.
In the following section we discuss the shower events, generated from τ decay, which is interpreted as the origin of the AAEs. Basically, τ which induces the air shower above Antarctic surface could also pass through the IceCube detector giving similar neutrino events, but no such event is observed by IceCube collaboration. In Ref. [15], the authors have interpreted 3 events of PeV energies in IceCube data [41] as sub-EeV earth emergent cosmic rays similar to ANITA. The exposure of IceCube is estimated at 54.0 km 2 sr yr which is 20 times that of ANITA exposure, 2.7 km 2 sr yr [15]. By using the relative exposures, the number of events N at ANITA can be estimated as N « N IC {20 giving N " 0.15 for N IC " 3 events at IceCube which is roughly an order of magnitude away from the observed events by ANITA. As pointed out in [15], ANITA estimated event rate serves as an upper limit which can be reduced for a better angular coverage and a detailed analysis including particular instrumental effects is needed to resolve it. In the following section, while remaining consistent with IceCube, we will estimate the required survival probability for N " 0.15. In the estimation of survival probability we have closely followed the approach of [19].
As discussed in the considered scenarios, the NLSP decays to LSP and a τ -lepton which is assumed to be the origin of the observed AAEs. While propagating through the earth, the NLSP may loose its energy through the electromagnetic processes. The average energy loss of a particle traveling a distance l can be estimated from the equation,´dE{dl " βpEqEρplq where ρplq is the density of the earth, and β is the radiative energy loss. While using the parametrization for βpEq from [47] we found β " 6.5ˆ10´1 0 cm 2 {gm for NLSP initial energy E " 1 EeV and mass M N LSP " 1 TeV. Then with ρ " 4 gm{cm 3 [48] for l " 7000 km one can find the energy loss due to radiation to be " 17%, implying the NLSP energy E " 0.8 EeV. The corresponding final shower energy E{2 " 0.4 EeV is consistent with the ANITA observations.

Model-I
The basic idea is to consider the ultra-high energy (UHE) neutrino interactions with the matter inside earth. At the elementary level the incoming UHE neutrino with energy E ν interacts with the quarks mediated by the appropriate new physics particle and produce the BSM particle which travels through the earth and decays to τ -lepton which actually gives rise to the AAEs. Given the particle content shown in table 1 and interactions in equation (1) the relevant processes in this case are νd Ñ Φ2 ψ 2 and νu Ñ Φ2 ψ 1 which are mediated by the particle E, where Φ2 is the charged scalar arising from the doublet Φ 2 . In this case the final τ -lepton is obtained from the decay of the Φ2 , which decays to N 1 and τ . In this model, N 1 is identified as DM candidate. In terms of Bjorken scaling variables x " Q 2 {2m N E 1 ν , y " E 1 ν {E ν , where m N " pm p`mn q{2, the average mass of proton and neutron, E 1 ν " E ν´E Φ2 is the energy loss in the lab frame,´Q 2 the invariant momentum transfer between neutrino (antineutrino) and particle Φ2 , the differential cross-section of neutrino (antineutrino) nucleon scattering can be written as, where s " 2m N E ν is the center of mass energy, f q px, Q 2 q and f q px, Q 2 q are the parton distribution functions (PDFs) of quarks and antiquarks with q " u, d. In our analysis, we consider the standard p1 : 1 : 1q flavor ratio of neutrino (antineutrino) on earth and use CTEQ6 leading order PDF sets [49] while calculating the cross-sections. By using the interaction cross-section obtained from equation (23), the interaction length, l Φ 2 , can be written as, where N A " 6.022ˆ10 23 cm´3 is the Avogadro number and ρ « 4 is the earth density in water equivalent unit. A convenient representation of the density profile of the earth is given in [48] which is used to calculate the approximate value of ρ for the relevant chord lengths for the AAEs. The decay length of Φ2 , l D , in the earth rest frame can be written as, where Γ Φ2 " ΓpΦ2 Ñ τ´N 1 q is the rest frame decay width of Φ2 into τ lepton and DM N 1 which is obtained from equation (7). Here we used the approximation E Φ2 " p1´ξqE ν {2, where ξ accounts for the energy loss of Φ2 due to radiation (which is " 17% as estimated before), to get the observed shower energy E τ " E Φ2 {2 " 0.4 EeV by taking into account the incident neutrino energy E ν " 2 EeV. As mentioned before, by considering only the SM interactions it is not possible to get the desired survival probability for any SM particle to pass a chord distance of " 7000 km. In considered BSM scenarios, as a result of additional interactions, the survival probability of the neutrino flux can be estimated as [19], where l C is the chord length, l 0 is the mean interaction length which is 290 and 265 km for ANITA-I and ANITA-III events respectively [15], and h « 10 km is the window of τ production near the surface of earth.
In figure 4, we plotted the survival probability BSM of Φ2 by normalizing it with respect to SM , as a function of M Φ2 , and it's decay coupling Y 3 with DM and tau lepton, for ANITA-I (left) and ANITA-III (right) events. While scanning the parameter space, we considered the benchmark relevant for the correct relic density which are quoted in table 2 and shown as a star in figure 4 for M Φ2 " 1 TeV and Y 3 " 2.92ˆ10´1 0 . The color scheme corresponds to different regions of parameter space where BSM { I,III SM ě 1. In our analysis, we considered I SM " 4.4ˆ10´7 and III SM " 3.2ˆ10´8 for ANITA-I and ANITA-III events respectively [15]. By using the benchmark in this model we get BSM " 2.45ˆ10´7 for ANITA-I event and BSM " 2.17ˆ10´7 for ANITA-III event. When these values are used in equation (22), we found N " 0.00009 for isotropic flux and N " 0.009 for anisotropic flux. Clearly, in this model, it is not possible to get the desired number of events for perturbative Yukawa couplings while being consistent with DM relic density estimation.

Model-II
The basic tenet of explaining the AAEs in this model follows a similar path as the Model-I, mentioned in the previous subsection. In the Model-II, outlined in subsection 2.1, the necessary interactions can be obtained from the Yukawa Lagrangian defined in equation (10). In this case, the typical relevant processes would be νd Ñ EU and νu Ñ ED mediated by the charged scalar arising from the doublet Φ 2 . The produced E subsequently decays to τ and the DM candidate which is the neutral component of Φ 3 . Now for the mentioned processes, the total differential cross-section of neutrino (antineutrino) scattering in terms of Bjorken scaling variables x " Q 2 {2m N E 1 ν , y " E 1 ν {E ν , with E 1 ν " E ν´EE the energy loss in the lab frame,´Q 2 the invariant momentum transfer between neutrino (antineutrino) and particle E, can be written as,  where s " 2m N E ν is the center of mass energy. Similar to the previous case the interaction length, l E , can be given as, where σ is the interaction cross-section obtained from Eq. (27). The decay length of E in the earth rest frame can be written as, where Γ E " ΓpE Ñ τ´Φ 0 3 q is the rest frame decay width of E obtained from equation (18). As before, we used the approximation E E " p1´ξqE ν {2, where ξ accounts particle E energy loss due to radiation (which is " 17% as estimated before), to get the observed shower energy E τ " E E {2 " 0.4 EeV by taking into account the incident neutrino energy E ν " 2 EeV. In this case, the survival probability of the neutrino flux, BSM , can also be written following equation (26) by making the substitution l Φ2 Ñ l E and the appropriate l D defined in equation (29).
In figure 5, we plotted the survival probability BSM of E normalized to SM , as a function of it's mass, M E , and it's decay coupling Y 1 2 to DM and tau lepton for ANITA-I (left) and ANITA-III (right) events. While scanning the parameter space, we considered the parameters relevant for correct relic density which are quoted in table 4 and specifically first benchmark point is shown as a star in figure 5 for M E " 1 TeV and Y 1 2 " 1.1ˆ10´9. The color scheme corresponds to different regions of parameter space where BSM { I,III SM ě 1. While using the first benchmark point, we get BSM " 4.33ˆ10´6 for ANITA-I event and BSM " 3.41ˆ10´6 for ANITA-III event. These values are then used in equation (22) and we found N " 0.0015 for isotropic flux, and N " 0.154 for anisotropic flux. We repeat the similar exercise for second benchmark point by taking into account the relic density estimations and get N " 0.17 for anisotropic flux. Our result for the second benchmark point is shown in figure 6. Clearly for both the benchmarks, it is not possible to get the desired survival probability for isotropic flux and so anisotropic flux is required while being consistent with DM relic density estimation.

Complementary Probes of the Model
Apart from providing an explanation to AAEs along with dark matter, the model we propose here can also offer tantalizing signatures at different experiments operating at energy as well as intensity frontiers. We briefly comment on such possibilities related to model-II which is successful in explaining the observed AAEs discussed before. Clearly, the model is an extension of the widely studied inert Higgs doublet model or scotogenic model of radiative neutrino mass with additional Z 2 odd particles. The scalar sector of model-II contains two Z 2 odd scalar doublets Φ 2,3 . By virtue of their electroweak gauge interactions, the components of these scalar doublets can be produced significantly in protonproton collisions of the Large Hadron Collider (LHC). Depending upon the mass spectrum of its components, the heavier ones can decay into the lighter ones and a gauge boson, which finally decays into a pair of leptons or quarks. Therefore, we can have either pure leptonic final states plus missing transverse energy (MET), hadronic final states plus MET or a mixture of both. The MET corresponds to DM or light neutrinos. In several earlier works [50][51][52], the possibility of opposite sign dileptons plus MET was discussed. In [53], the possibility of dijet plus MET was investigated with the finding that inert scalar masses up to 400 GeV can be probed at high luminosity LHC. In another work [54], trilepton plus MET final states was also discussed whereas mono-jet signatures have been studied by the authors of [55,56]. The enhancement in dilepton plus MET signal in the presence of additional vector like singlet charged leptons (like we have in model-II) was also discussed in [57]. Exotic signatures like displaced vertex and disappearing or long-lived charged track for compressed mass spectrum of inert scalars and singlet fermion DM was studied recently by the authors of [58].
To summarise, the search strategy and the bounds on Z 2 odd scalars depend crucially on the model and spectrum of lighter particles to which they can decay. Due to identical gauge interactions, the bounds are somewhat similar to the ones on sleptons in supersymmetric models. As mentioned in [33], collider bounds on charged sleptons are also dependent upon lightest neutralino mass among other details related to the mass spectrum. Since we are not performing a detailed collider analysis in our work, we consider a conservative lower bound on such Z 2 odd charged scalar masses (ą 100 GeV) in agreement with LEP bounds [33,59]. Also, the presence of Z 2 odd vector like singlet charged leptons can lead to dilepton plus MET signatures at colliders. As noted from the collider studies in [57,60] within a similar model, the chosen values of vector like charged lepton masses are safe from existing collider limits. As can be seen from the above analysis, such bounds are satisfied by the parameter space shown in figure 4, 5, and 6.
Apart from collider signatures of inert scalars in the model, as discussed in several earlier works mentioned above, the model can also have interesting signatures at lepton flavor violating decays like µ Ñ eγ, µ Ñ eee and µ Ñ e conversions [61,62]. With inclusion of vector-like quarks and vector-like charged leptons, like we have in model-II, one can also explain flavor anomalies in b Ñ s decays as well as muon g´2 [60]. Thus, the model we have studied here, which successfully explains the ANITA anomalous events also have rich phenomenology that can be accessed at different experiments operating at cosmic, energy as well as intensity frontiers.

Summary and Conclusions
We have proposed a beyond standard model framework to explain the two anomalous upward going ultra high energy cosmic ray air shower events reported by the ANITA collaboration. The novel feature of the model is the way it relates the origin of AAEs to non-thermal origin of dark matter in the universe as well as the origin of light neutrino masses. Sticking to minimality we extend the standard model by few additional particles, all of which are odd under an in-built and unbroken Z 2 symmetry. While two types of such Z 2 odd fields namely, a scalar doublet and gauge singlet neutral fermions play the role of generating light neutrino masses at one loop level similar to the scotogenic scenarios, the other fields are responsible for generating the AAEs. The non-thermal nature of dark matter is related to its tiny coupling with the next to lightest Z 2 odd particle whose long-livedness makes its passage through earth possible followed by its decay into DM and a tau lepton required to explain the AAEs.
We first consider a dark matter scenario where its relic is generated purely from non-thermal production. In this case, one of the three Z 2 odd gauge singlet neutral fermions is the DM candidate. Being gauge singlet, its non-thermal nature is dictated by its tiny coupling with the SM leptons and Z 2 odd scalar doublet. While such tiny coupling results in vanishingly small lightest neutrino mass [37], the abundance of DM is generated from decay of the Z 2 odd scalar doublet, the NLSP in this case. While the correct DM relic is obtained by suitable choices of DM, NLSP masses and their coupling, it was found that the model can not explain the AAEs. We then consider a hybrid setup where DM relic is obtained from both thermal as well as non-thermal contribution. The neutral component of a second Z 2 odd scalar doublet is considered to be the DM candidate which annihilates into the SM particles at a rate larger than the one required for correct thermal relic abundance. This leads to an under-abundant thermal dark matter and the deficit can be filled in by a non-thermal contribution from a singlet fermion, considered to be a charged vector like singlet fermion E. The light neutrino masses arise as usual from gauge singlet neutral fermions and another Z 2 odd scalar doublet, thereby not getting affected by the smallness of DM coupling with E. The model is found to explain the AAEs successfully.
Although we have confined our discussions here to DM relic and anomalous ANITA events, the Model-II of our study, successful in explaining both of these, can have other consequences which can be tested at experiments operating at different frontiers. For example, the components of Z 2 odd scalar doublet can be produced at a significant rate at collider experiments like the LHC due to its electroweak gauge interactions. Our model also has additional colored vector like Z 2 odd fermions which can also give rise to similar signatures at colliders. While we briefly mention such additional aspects of probing this model at variety of experiments, the detailed phenomenological study of such scenarios is beyond the scope of this present work and left for future works.