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

TL;DR

This paper proposes a non-thermal dark matter scenario involving a Z2 odd sector to explain ANITA's anomalous upward-going events, suggesting a hybrid thermal and non-thermal dark matter relic abundance.

Contribution

It introduces a novel Z2 odd sector model where NLSP decay produces tau leptons, linking dark matter properties to ANITA's observed anomalies.

Findings

Hybrid dark matter scenario explains ANITA events.

Purely non-thermal dark matter cannot account for event rates.

Long-lived NLSP decays produce tau leptons passing through Earth.

Abstract

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.

Tables

Table 1. Table 1: Field content and transformation properties under the symmetry of the model.

Fields	$Q$	$u_R$	$d_R$	$L$	$e_R$	${\Phi }_1$	${\Phi }_2$	$E_{L,R}$	$N$	${\psi }_1$	${\psi }_2$
$SU\left(3\right)$	3	3	3	1	1	1	1	1	1	3	3
$SU\left(2\right)$	2	1	1	2	1	2	2	1	1	1	1
$U{\left(1\right)}_Y$	$\frac{1}{6}$	$\frac{2}{3}$	$-\frac{1}{3}$	-$\frac{1}{2}$	-1	$\frac{1}{2}$	$\frac{1}{2}$	$\pm 1$	0	$\frac{1}{3}$	$\frac{4}{3}$
$Z_2$	1	1	1	1	1	1	-1	-1	-1	-1	-1

Table 2. Table 2: Numerical values of different relevant parameters used for Model-I.

Parameters	Benchmark
${\mu }_2$	1 TeV
${\mu }_3$	6 TeV
${\mu }_4$	6 TeV
$Y_3$	$2.92\times {10}^{-10}$
$Y_2,Y_4,Y_5$	3.5
$M_{{\Phi }_2}^{\pm }$	1 TeV
$M_{{\Phi }_2}^I$	1.045 TeV
$M_{{\Phi }_2}^R$	1.073 TeV
${\lambda }_1$	0.255
${\lambda }_2$	0.1
${\lambda }_3$	$0$
${\lambda }_4$	$4.0$
${\lambda }_5$	$1.0$
$M_{N_1}$	30 MeV
$M_E$	7 TeV

Table 3. Table 3: Field content and transformation properties under the symmetry of the model.

Fields	$Q$	$u_R$	$d_R$	$L$	$e_R$	${\Phi }_1$	${\Phi }_2$	${\Phi }_3$	$E_{L,R}$	$N$	$U_{L,R}$	$D_{L,R}$
$SU\left(3\right)$	3	3	3	1	1	1	1	1	1	1	3	3
$SU\left(2\right)$	2	1	1	2	1	2	2	2	1	1	1	1
$U{\left(1\right)}_Y$	$\frac{1}{6}$	$\frac{2}{3}$	$-\frac{1}{3}$	-$\frac{1}{2}$	-1	$\frac{1}{2}$	$\frac{1}{2}$	$\frac{1}{2}$	$\pm 1$	0	$\pm \frac{2}{3}$	$\mp \frac{1}{3}$
$Z_2$	1	1	1	1	1	1	-1	-1	-1	-1	-1	-1

Table 4. Table 4: Numerical values of different relevant parameters used for Model-II.

Parameters	First benchmark	Second benchmark
${\mu }_2$	400 GeV	300 GeV
$Y_2^{\prime },Y_2^{{}^{\prime \prime }}$	$1.1\times {10}^{-9}$	$4.9\times {10}^{-10}$
$Y_2,Y_4,Y_5$	0.96	1.27
${\lambda }_1$	0.255	0.255
${\lambda }_2$	0.1	0.1
${\lambda }_3$	0.81	0.81
${\lambda }_4$	-0.3	-0.3
${\lambda }_5$	$1\times {10}^{-8}$	$1\times {10}^{-8}$
$M_{{\Phi }_2^{\pm }}$	1.1 TeV	1.1 TeV
$M_{{\Phi }_3^0}$	416.2 GeV	321.4 GeV
$M_{{\Phi }_3^{\pm }}$	427.6 GeV	335.9 GeV
$M_U$	5 TeV	5 TeV
$M_D$	5 TeV	5 TeV
$M_E$	1 TeV	1 TeV
$M_N$	5 TeV	5 TeV

Equations

$${\mathcal{L}}_Y$$

$$+Y_4^{\prime }{u}_R{\psi }_1{E}_R+Y_5^{\prime }{d}_R{\psi }_2{E}_R+Y_6{d}_R{\psi }_{1}N.$$

$$V({\Phi }_1,{\Phi }_2,{\psi }_1,{\psi }_2)$$

$$+{\lambda }_4|{\Phi }_1^{†}{\Phi }_2{|}^2+\left\{\frac{{\lambda }_5}{2}({\Phi }_1^{†}{\Phi }_2{)}^2+\text{h.c.}\right\}+{\mu }_3^2|{\psi }_1{|}^2+{\mu }_4^2|{\psi }_2{|}^2$$

$$+|{\Phi }_1{|}^2({\lambda }_6|{\psi }_1{|}^2+{\lambda }_7|{\psi }_2{|}^2)+|{\Phi }_2{|}^2({\lambda }_8|{\psi }_1{|}^2+{\lambda }_9|{\psi }_2{|}^2)$$

$$+\left\{\frac{{\lambda }_{10}}{2}({\Phi }_1^{†}{)}^2({\psi }_1^{†}{\psi }_2)+\text{h.c.}\right\}+\left\{\frac{{\lambda }_{11}}{2}({\Phi }_2^{†}{)}^2({\psi }_1^{†}{\psi }_2)+\text{h.c.}\right\}$$

$$+{\lambda }_{12}|{\psi }_1{|}^4+{\lambda }_{13}|{\psi }_2{|}^4.$$

$${\Phi }_1=\left(\begin{array}{c}0\\ \frac{1}{\sqrt{2}}(v+h)\end{array}\right),{\Phi }_2=\left(\begin{array}{c}{\Phi }_2^{+}\\ \frac{1}{\sqrt{2}}({\Phi }_2^R+i{\Phi }_2^I)\end{array}\right),$$

$$m_h^2$$

$$(M_{{\Phi }_2}^R{)}^2$$

$$(M_{{\Phi }_2}^I{)}^2$$

$$(M_{{\Phi }_2}^{\pm }{)}^2$$

$$\frac{d{Y}_{{\Phi }_2^{\pm }}}{dz}=-\frac{⟨\sigma v⟩s}{Hz}\left(Y_{{\Phi }_2^{\pm }}^2-(Y_{{\Phi }_2^{\pm }}^{\mathrm{eq}}{)}^2\right)-\frac{⟨{\Gamma }_{{\Phi }_2^{\pm }}⟩}{Hz}{Y}_{{\Phi }_2^{\pm }}$$

$$\frac{d{Y}_{N_1}}{dz}=\frac{⟨{\Gamma }_{{\Phi }_2^{\pm }}⟩}{Hz}{Y}_{{\Phi }_2^{\pm }},$$

$$⟨{\Gamma }_{{\Phi }_2^{\pm }}⟩$$

$$={\Gamma }_{{\Phi }_2^{\pm }}\frac{K_1(M_{{\Phi }_2^{\pm }}/T)}{K_2(M_{{\Phi }_2^{\pm }/}T)},$$

$$Y_{N_1}$$

$$\simeq 4.71239\left[\frac{{\Gamma }_{{\Phi }_2^{\pm }}}{16{\pi }^{4}H(M_{{\Phi }_2^{\pm }})}\frac{45}{g_{\ast }}\right]$$

$${\Omega }_{DM}{h}^2$$

$${\mathcal{L}}_Y$$

$$+Y_2^{\prime }\overline{L}{\Phi }_3{E}_R+\widetilde{Y_3^{\prime }}\overline{L}\widetilde{{\Phi }_3}N+Y_4^{\prime }\overline{Q}\widetilde{{\Phi }_3}{U}_R+Y_5^{\prime }\overline{Q}{\Phi }_3{D}_R$$

$$+Y_2^{\prime \prime }{L}^T\widetilde{{\Phi }_2}{E}_L++Y_3^{\prime \prime }{L}^T\widetilde{{\Phi }_3}{E}_L.$$

$$V({\Phi }_1,{\Phi }_2,{\Phi }_3)$$

$$+{\lambda }_4|{\Phi }_1{|}^2|{\Phi }_2{|}^2+{\lambda }_5|{\Phi }_1^{†}{\Phi }_2{|}^2+\left\{\frac{{\lambda }_6}{2}({\Phi }_1^{†}{\Phi }_2{)}^2+\text{h.c.}\right\}$$

$$+{\lambda }_4^{\prime }|{\Phi }_1{|}^2|{\Phi }_3{|}^2+{\lambda }_5^{\prime }|{\Phi }_1^{†}{\Phi }_3{|}^2+\left\{\frac{{\lambda }_6^{\prime }}{2}({\Phi }_1^{†}{\Phi }_3{)}^2+\text{h.c.}\right\}$$

$$+\{{\mu }_{23}^2{\Phi }_2^{†}{\Phi }_3+\text{h.c}\}+{\lambda }_7|{\Phi }_2{|}^2|{\Phi }_3{|}^2+{\lambda }_8|{\Phi }_2^{†}{\Phi }_3{|}^2$$

$$+\{{\lambda }_7^{\prime }|{\Phi }_1{|}^2{\Phi }_2^{†}{\Phi }_3+{\lambda }_8^{\prime }{\Phi }_2^{†}{\Phi }_3|{\Phi }_3{|}^2+{\lambda }_9^{\prime }{\Phi }_3^{†}{\Phi }_2|{\Phi }_2{|}^2+\text{h.c.}\}$$

$$+\left\{\frac{{\lambda }_9}{2}({\Phi }_2^{†}{\Phi }_3{)}^2+\text{h.c.}\right\}.$$

$$\frac{d{n}_{{\Phi }_3^0}}{dt}+3H{n}_{{\Phi }_3^0}=-⟨\sigma v{⟩}_{{\Phi }_3^0{\Phi }_3^0\to \mathrm{SM}\mathrm{SM}}(n_{{\Phi }_3^0}^2-(n_{{\Phi }_3^0}^{\mathrm{eq}}{)}^2)+N_{{\Phi }_3^0}{\Gamma }_E{n}_E$$

$$\frac{d{n}_E}{dt}+3H{n}_E=-⟨\sigma v{⟩}_{EE\to \mathrm{SM}\mathrm{SM}}(n_E^2-(n_E^{\mathrm{eq}}{)}^2)-{\Gamma }_E{n}_E.$$

$$\frac{d{Y}_{{\Phi }_3^0}}{dz}=-\frac{⟨\sigma v⟩s}{Hz}\left(Y_{{\Phi }_3^0}^2-(Y_{{\Phi }_3^0}^{\mathrm{eq}}{)}^2\right)+\frac{N_{{\Phi }_3^0}{\Gamma }_N}{Hz}{Y}_E$$

$$\frac{d{Y}_E}{dz}=-\frac{{\Gamma }_E}{Hz}{Y}_E.$$

$$Y_E(z)=Y_E(z_F)\text{exp}\left(-\frac{r}{2}\left(z^2-z_F^2\right)\right).$$

$$\frac{d{Y}_{{\Phi }_3^0}}{dz}=-\frac{⟨\sigma v⟩s}{Hz}\left(Y_{{\Phi }_3^0}^2-(Y_{{\Phi }_3^0}^{\mathrm{eq}}{)}^2\right)+N_{{\Phi }_3^0}rz{Y}_E(z_0)\text{exp}\left(-\frac{r}{2}\left(z^2-z_0^2\right)\right).$$

$${\Gamma }_E=\frac{Y_2^{\prime 2}}{16\pi }{M}_E{\left(1-\frac{M_{{\Phi }_3^0}^2}{M_E^2}\right)}^2$$

$${\tau }_E=\frac{1}{{\Gamma }_E}=10\left(\frac{M_E}{500\mathrm{GeV}}\right)ns.$$

$$\mathcal{N}=A_{\mathrm{eff}}T\Delta \Omega {\int }_{E_i}^{E_f}d{E}_{\nu }\Phi (E_{\nu })ϵ(E_{\nu }),$$

$$\Delta \Omega \approx \{\begin{array}{l}2\pi \text{sr for isotropic case,}\\ 2\pi (1-\cos {\delta }_{\theta })\approx 0.0022\text{sr for anisotropic case,}\end{array}$$

$$\mathcal{N}\approx \{\begin{array}{l}2.0\times 10^2ϵ\text{for isotropic case,}\\ 2.0\times 10^4ϵ\text{for anisotropic case.}\end{array}$$

Full text

1 Introduction

Recently, the ANtarctic Impulsive Transient Antenna (ANITA) collaboration has reported two anomalous upward going ultra high energy cosmic ray (UHECR) air shower events corresponding to deposited energies of $0.6\pm 0.4$ EeV and $0.56^{+0.3}_{-0.2}$ EeV respectively [1, 2]. These two events, which originate from well below the horizon, with angles of elevation $-27.4^{\circ}\pm 0.3^{\circ}$ and $-35.0^{\circ}\pm 0.3^{\circ}$ respectively, have $\geq 3\sigma$ evidence for occurring due to direct upward moving earth emergent UHECR like air showers above the Antarctic ice surface [2]. However, no particle in the standard model (SM) can survive such a passage through earth, a distance of several thousand kilometers corresponding to the observed zenith angles for ANITA anomalous events (AAEs). Therefore, several beyond standard model (BSM) proposals have been put forward to explain the observed events. They include sterile neutrino mixing [3, 4], heavy dark matter (DM) [5, 6, 7, 8, 9, 10, 11, 12], inelastic boosted DM [13, 14], long lived charged particles like stau [15, 16], R-parity violating supersymmetry (SUSY) [17], SUSY sphalerons configurations [18], leptoquarks [19], and radio pulses from axion-photon conversion [20]. Additionally, in Ref. [21], $\nu_{\tau}$ neutrino origin of AAEs is studied in compliance with Auger [22] and IceCube [23] upper limits. Note that explanations of AAEs in SM scenario are also proposed [24, 25], but they remain speculative and require experimental confirmations.

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 $\Omega_{\rm DM}$ and $h=\text{Hubble Parameter}/(100\;\text{km}~{}\text{s}^{-1}\text{Mpc}^{-1})$, the present DM abundance is conventionally reported as [26]: $\Omega_{\text{DM}}h^{2}=0.120\pm 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.

2 Model

2.1 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 $\Phi_{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},\psi_{1},\psi_{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 relevant part of the Yukawa Lagrangian is given by,

[TABLE]

The scalar potential is,

[TABLE]

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 $\Phi_{2}$. The decay of $\Phi_{2}$ can occur either in equilibrium or after its freeze-out, the latter is similar to the superWIMP formalism [39]. It should be noted that all the components of $\Phi_{2}$ can contribute to the non-thermal production of $N_{1}$ namely, $\Phi^{\pm}_{2}\rightarrow N_{1}l^{\pm},\Phi^{0}_{2}\rightarrow N_{1}\nu,\;l\equiv e,\mu,\tau$. Since the $Z_{2}$ symmetry remains unbroken, the neutral components of $\Phi_{2}$ do not acquire any non-zero vacuum expectation value (VEV) which can be ensured by choosing $\mu_{2}^{2}>0$. We parametrize the two scalar doublets $\Phi_{1,2}$ as

[TABLE]

where $\Phi_{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,

[TABLE]

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 $\Phi^{\pm}_{2}$ as well as $N_{1}$. They can be written in terms of their comoving number densities $(Y=n/s,n\equiv{\rm number\;density},s\equiv{\rm entropy\;density})$ as,

[TABLE]

[TABLE]

with $z=M_{\Phi^{\pm}_{2}}/T$ and thermally averaged decay width of $\Phi^{\pm}_{2}$

[TABLE]

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 $\Phi^{\pm}_{2}$ is still in equilibrium. So, we can safely integrate the equation (6). So, the solution of the equation (6) is given as

[TABLE]

and the relic is given as

[TABLE]

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\;{\rm MeV},M_{\Phi^{\pm}_{2}}=1\;{\rm TeV},Y_{3}=2.92\times 10^{-10}$ 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 $\Phi_{2}$. The equilibrium contribution of $\Phi^{\pm}_{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 $\Phi_{2}$, the coannihilations among different components of $\Phi_{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.

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

The relevant part of the Yukawa Lagrangian is given by,

[TABLE]

The scalar potential is,

[TABLE]

We assume the neutral component $\Phi^{0}_{3}$ of $\Phi_{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 $\Phi_{2},\Phi_{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 $\mu_{23},\lambda_{7},\lambda_{8},\lambda_{9}$ in the scalar potential (11).

The coupled Boltzmann equations can be written as,

[TABLE]

[TABLE]

Here, $N_{\Phi^{0}_{3}}$ is the average number of DM particles produced from a single decay of $E$ which is one in this case, $\Gamma_{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_{\star s}$) and the comoving energy density ($g_{\star}$) 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 555Note that the first term on the right hand side of equation (13) can also be ignored if $E$ decays while in thermal equilibrium $n_{E}=n^{\rm eq}_{E}$, we consider the general case in our numerical calculations. and one can analytically solve the Boltzmann equation for $n_{E}$ above. Writing the above equations in terms of $Y_{\Phi^{0}_{3}}=\frac{n_{\Phi^{0}_{3}}}{s},Y_{E}=\frac{n_{E}}{s}$ with $s=\frac{2\pi^{2}}{45}g_{*s}T^{3}$ being the entropy density and changing the variable from time $t$ to $z=M_{\Phi^{0}_{3}}/T$, where $M_{\Phi^{0}_{3}}$ is the mass of DM $\Phi^{0}_{3}$, we get,

[TABLE]

[TABLE]

The equation for $Y_{E}$ can be solved analytically giving,

[TABLE]

Here, $z_{F}=M_{\Phi^{0}_{3}}/T_{F}$ is the point of freeze-out and usually takes a value of $\mathcal{O}(20)$. Also, $r=\dfrac{\Gamma_{E}}{H\,z^{2}}=\dfrac{\Gamma_{E}\,M_{\rm Pl}}{\pi M_{\Phi^{0}_{3}}^{2}}\sqrt{90/g_{\star}}$, depends on the decay width of the mother particle $E$. $Y_{E}(z_{0})$ (where $z_{0}$ corresponds 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_{\Phi^{0}_{3}}$ can be rewritten as,

[TABLE]

The decay width of fermion $E$ (assuming $Y^{\prime}_{3}=Y^{{}^{\prime\prime}}_{3}=Y^{\prime}_{2}/2$), ignoring the mass of leptons is given by,

[TABLE]

where $M_{E}$ is the mass of fermion $E$ and $Y^{\prime}_{2}$, its coupling with the tau lepton and DM. The lifetime requirement of $E$ from ANITA anomaly point of view is [15],

[TABLE]

For $M_{E}\sim 5$ TeV, the required lifetime of 100 ns, can be obtained for $Y^{\prime}_{2}\approx 10^{-10}-10^{-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 $\Phi_{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^{{}^{\prime\prime}}_{2}$. By keeping $\Phi_{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.

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 $\Phi^{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 $\Phi^{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 $\Phi^{0}_{3}$ with the SM particles freezes out giving rise to the plateau shaped regions. However, since $E$ is unstable and decays into $\Phi^{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 $\Phi^{0}_{3}$. In the absence of $E$ decay, the freeze-out abundance of $\Phi^{0}_{3}$ remains under-abundant compared to the Planck bound. The mass of $\Phi^{0}_{3}$ in both the benchmark points given in table 4 falls in the mass regime $80\;\text{GeV}\leq M_{\Phi^{0}_{3}}\leq 500\;\text{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].

3 Explanation of ANITA Events

The estimation of number of AAEs can be facilitated from the evaluation of survival probability $\epsilon$ which basically represents the fraction of incoming ultra-high energy neutrino flux that ultimately gives rise to $\tau$-leptons near the earth surface. The total number of events can be estimated as,

[TABLE]

where the effective area $A_{\rm eff}\approx 4$ km2, $T$ is the exposure time which we take $\approx 25$ days by combining the exposures of ANITA-I and ANITA-III, $\Delta\Omega$ is the acceptance angle which corresponding to isotropic and anisotropic neutrino flux source is given as,

[TABLE]

where $\delta_{\theta}\approx 1.5^{\circ}$ 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 $\tau$ lepton as the origin of observed events, with $E_{\tau}=E_{\nu}/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 $\Phi(E_{\nu})$ 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, $\Phi_{\rm iso}\approx 10^{-25}$ (GeV cm2 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 $\Phi_{\rm aniso}\approx 3.2\times 10^{-20}$ (GeV cm2 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,

[TABLE]

This implies that for the observed two AAEs in the isotropic case the $\epsilon$ should be $\sim 10^{-2}$, and for the anisotropic case $\epsilon$ should be $\sim 10^{-4}$. By considering SM interactions, authors of [15] predicted the survival probability $\epsilon_{\rm SM}\sim 10^{-7}$, giving $\mathcal{N}\sim 2\times 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 $\mathcal{N}\sim 2\times 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 $\epsilon(E_{\nu})$ 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 $\tau$ decay, which is interpreted as the origin of the AAEs. Basically, $\tau$ 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~{}\rm km^{2}~{}sr~{}yr$ which is 20 times that of ANITA exposure, $2.7~{}\rm km^{2}~{}sr~{}yr$ [15]. By using the relative exposures, the number of events $\mathcal{N}$ at ANITA can be estimated as $\mathcal{N}\approx\mathcal{N_{IC}}/20$ giving $\mathcal{N}=0.15$ for $\mathcal{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 $\mathcal{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 $\tau$-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=\beta(E)E\rho(l)$ where $\rho(l)$ is the density of the earth, and $\beta$ is the radiative energy loss. While using the parametrization for $\beta(E)$ from [47] we found $\beta=6.5\times 10^{-10}$ cm${}^{2}/$gm for NLSP initial energy $E=1$ EeV and mass $M_{NLSP}=1$ TeV. Then with $\rho=4~{}\rm{gm/cm^{3}}$ [48] for $l\sim 7000$ km one can find the energy loss due to radiation to be $\sim 17\%$, implying the NLSP energy $E\sim 0.8$ EeV. The corresponding final shower energy $E/2\sim 0.4$ EeV is consistent with the ANITA observations.

3.1 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_{\nu}$ interacts with the quarks mediated by the appropriate new physics particle and produce the BSM particle which travels through the earth and decays to $\tau$-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 $\nu\overline{d}\to\Phi_{2}^{-}\psi_{2}$ and $\nu\overline{u}\to\Phi_{2}^{-}\psi_{1}$ which are mediated by the particle $E$, where $\Phi_{2}^{-}$ is the charged scalar arising from the doublet $\Phi_{2}$. In this case the final $\tau$-lepton is obtained from the decay of the $\Phi_{2}^{-}$, which decays to $N_{1}$ and $\tau$. In this model, $N_{1}$ is identified as DM candidate.

In terms of Bjorken scaling variables $x=Q^{2}/2m_{N}E^{\prime}_{\nu}$, $y=E^{\prime}_{\nu}/E_{\nu}$, where $m_{N}=(m_{p}+m_{n})/2$, the average mass of proton and neutron, $E^{\prime}_{\nu}=E_{\nu}-E_{\Phi^{-}_{2}}$ is the energy loss in the lab frame, $-Q^{2}$ the invariant momentum transfer between neutrino (antineutrino) and particle $\Phi^{-}_{2}$, the differential cross-section of neutrino (antineutrino) nucleon scattering can be written as,

[TABLE]

where $s=2m_{N}E_{\nu}$ is the center of mass energy, $f_{q}(x,Q^{2})$ and $f_{\overline{q}}(x,Q^{2})$ are the parton distribution functions (PDFs) of quarks and antiquarks with $q=u,d$. In our analysis, we consider the standard $(1:1:1)$ 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_{\Phi_{2}}$, can be written as,

[TABLE]

where $N_{A}=6.022\times 10^{23}~{}\rm cm^{-3}$ is the Avogadro number and $\rho\approx 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 $\rho$ for the relevant chord lengths for the AAEs. The decay length of $\Phi^{-}_{2}$, $l_{D}$, in the earth rest frame can be written as,

[TABLE]

where $\Gamma_{\Phi^{-}_{2}}\equiv\Gamma({\Phi^{-}_{2}\rightarrow\tau^{-}N_{1}})$ is the rest frame decay width of $\Phi^{-}_{2}$ into $\tau$ lepton and DM $N_{1}$ which is obtained from equation (7). Here we used the approximation $E_{\Phi^{-}_{2}}=(1-\xi)E_{\nu}/2$, where $\xi$ accounts for the energy loss of ${\Phi^{-}_{2}}$ due to radiation (which is $\sim 17\%$ as estimated before), to get the observed shower energy $E_{\tau}=E_{\Phi^{-}_{2}}/2\sim 0.4$ EeV by taking into account the incident neutrino energy $E_{\nu}=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 $\sim 7000$ km. In considered BSM scenarios, as a result of additional interactions, the survival probability of the neutrino flux can be estimated as [19],

[TABLE]

where $l_{\Earth}$ 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\approx 10$ km is the window of $\tau$ production near the surface of earth.

In figure 4, we plotted the survival probability $\epsilon_{BSM}$ of $\Phi^{-}_{2}$ by normalizing it with respect to $\epsilon_{SM}$, as a function of $M_{\Phi^{-}_{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_{\Phi^{-}_{2}}=1$ TeV and $Y_{3}=2.92\times 10^{-10}$. The color scheme corresponds to different regions of parameter space where $\epsilon_{BSM}/\epsilon^{I,III}_{SM}\geq 1$. In our analysis, we considered $\epsilon^{I}_{SM}=4.4\times 10^{-7}$ and $\epsilon^{III}_{SM}=3.2\times 10^{-8}$ for ANITA-I and ANITA-III events respectively [15]. By using the benchmark in this model we get $\epsilon_{BSM}=2.45\times 10^{-7}$ for ANITA-I event and $\epsilon_{BSM}=2.17\times 10^{-7}$ for ANITA-III event. When these values are used in equation (22), we found $\mathcal{N}=0.00009$ for isotropic flux and $\mathcal{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.

3.2 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 (2.2). In this case, the typical relevant processes would be $\nu d\to EU$ and $\nu\overline{u}\to E\overline{D}$ mediated by the charged scalar arising from the doublet $\Phi_{2}$. The produced $E$ subsequently decays to $\tau$ and the DM candidate which is the neutral component of $\Phi_{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^{\prime}_{\nu}$, $y=E^{\prime}_{\nu}/E_{\nu}$, with $E^{\prime}_{\nu}=E_{\nu}-E_{E}$ the energy loss in the lab frame, $-Q^{2}$ the invariant momentum transfer between neutrino (antineutrino) and particle $E$, can be written as,

[TABLE]

where $s=2m_{N}E_{\nu}$ is the center of mass energy. Similar to the previous case the interaction length, $l_{E}$, can be given as,

[TABLE]

where $\sigma$ is the interaction cross-section obtained from Eq. (27). The decay length of $E$ in the earth rest frame can be written as,

[TABLE]

where $\Gamma_{E}\equiv\Gamma({E\rightarrow\tau^{-}\Phi^{0}_{3}})$ is the rest frame decay width of $E$ obtained from equation (18). As before, we used the approximation $E_{E}=(1-\xi)E_{\nu}/2$, where $\xi$ accounts particle $E$ energy loss due to radiation (which is $\sim 17\%$ as estimated before), to get the observed shower energy $E_{\tau}=E_{E}/2\sim 0.4$ EeV by taking into account the incident neutrino energy $E_{\nu}=2$ EeV.

In this case, the survival probability of the neutrino flux, $\epsilon_{BSM}$, can also be written following equation (26) by making the substitution $l_{\Phi_{2}^{-}}\to l_{E}$ and the appropriate $l_{D}$ defined in equation (29).

In figure 5, we plotted the survival probability $\epsilon_{BSM}$ of $E$ normalized to $\epsilon_{SM}$, as a function of it’s mass, $M_{E}$, and it’s decay coupling $Y^{\prime}_{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^{\prime}_{2}=1.1\times 10^{-9}$. The color scheme corresponds to different regions of parameter space where $\epsilon_{BSM}/\epsilon^{I,III}_{SM}\geq 1$. While using the first benchmark point, we get $\epsilon_{BSM}=4.33\times 10^{-6}$ for ANITA-I event and $\epsilon_{BSM}=3.41\times 10^{-6}$ for ANITA-III event. These values are then used in equation (22) and we found $\mathcal{N}=0.0015$ for isotropic flux, and $\mathcal{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 $\mathcal{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.

4 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 $\Phi_{2,3}$. By virtue of their electroweak gauge interactions, the components of these scalar doublets can be produced significantly in proton-proton 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], tri-lepton 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 [59, 33]. 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 $\mu\rightarrow e\gamma,\mu\rightarrow eee$ and $\mu\rightarrow 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\rightarrow 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.

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

Acknowledgements

DB acknowledges the support from IIT Guwahati start-up grant (reference number: xPHYSUGI-ITG01152xxDB001), Early Career Research Award from DST-SERB, Government of India (reference number: ECR/2017/001873) and Associateship Programme of Inter University Centre for Astronomy and Astrophysics (IUCAA), Pune. The research of UKD is supported by the Ministry of Science, ICT & Future Planning of Korea, the Pohang City Government, and the Gyeongsangbuk-do Provincial Government. The research of GT is supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, grant number 2019R1F1A1052231. DB and UKD thanks the organizers of WHEPP XVI where important discussions are carried out. We would like to thank Bhavesh Chauhan for useful discussion.