Planck Scale Origin of Non-zero $\theta_{13}$ and Super-WIMP Dark Matter

We study a discrete flavour symmetric scenario for neutrino mass and dark matter under the circumstances where such global discrete symmetries can be explicitly broken at Planck scale, possibly by gravitational effects. Such explicit breaking of discrete symmetries mimic as Planck suppressed operators in the model which can have non-trivial consequences for neutrino and dark matter sectors. In particular, we study a flavour symmetric model which, at renormalisable level, gives rise to tri-bimaximal type neutrino mixing with vanishing reactor mixing angle $\theta_{13}=0$, a stable inert scalar doublet behaving like a weakly interacting massive particle (WIMP) and a stable singlet inert fermion which does not interact with any other particles. The introduction of Planck suppressed operators which explicitly break the discrete symmetries, can give rise to the generation of non-zero $\theta_{13}$ in agreement with neutrino data and also open up decay channels of inert scalar doublet into singlet neutral inert fermion leading to the realisation of the super-WIMP dark matter scenario. We show that the correct neutrino phenomenology can be obtained in this model while discussing three distinct realisation of the super-WIMP dark matter scenario.


I. INTRODUCTION
Non-zero but tiny neutrino mass and large leptonic mixing have become a well established fact by now, thanks to several experimental efforts in the last two decades [1][2][3][4][5][6][7][8][9]. For a review of neutrino mass and mixing, please see [10]. Among these experiments, the relatively new ones namely, the T2K [5], the Double Chooz [6], the Daya Bay [7], the RENO [8] and the MINOS [9] have not only reconfirmed the fact that neutrinos oscillate by making the measurements of mass squared differences and mixing angles more precise, but also discovered a non-zero value of of reactor mixing angle θ 13 which was thought to be almost vanishing earlier. For a recent global fit of neutrino oscillation data, we refer to [11,12].
Apart from neutrino oscillation experiments, cosmology experiment like Planck also constrains the neutrino sector by putting an upper bound on the sum of absolute neutrino masses i |m i | < 0.11 eV [13]. As indicated by these global fits, in spite of these experimental developments, there still remains several unknowns in the neutrino sector like the nature of neutrinos: Dirac or Majorana, mass hierarchy: normal or inverted, CP violating phases etc. Apart from these, the origin of light neutrino mass also remains unknown in the standard model (SM) as the Higgs field do not have any renormalisable coupling to the neutrinos due to the absence of right handed neutrinos. If we go beyond renormalisable level, then it is possible to generate light neutrino mass of Majorana type via a dimension five Weinberg operator [14] of type (LLHH)/Λ where Λ is an unknown cutoff scale somewhere above the electroweak scale. Usual seesaw models [15][16][17][18] for light neutrino mass attempt to provide a dynamical origin of the Weinberg operator by introducing new fields like heavy right handed neutrinos. Apart from the sub-eV mass, well below the electroweak scale, another puzzling feature related to neutrinos is their large mixing, in sharp contrast with small mixing angles in the quark sector. This has also motivated the study of different flavour symmetry models that can predict such large mixing patterns. One of the very popular flavour symmetric scenarios is the one that predicts a µ − τ symmetric light neutrino mass matrix. However, such a scheme which predicts θ 13 = 0, θ 23 = π 4 and different values of θ 12 depending upon the particular realisation of this symmetry [19], has already been ruled out by recent neutrino experiments. Among such realisations, the tri-bimaximal (TBM) mixing [20][21][22][23][24] which predicts θ 12 = 35.3 o has been the most popular one. This mixing, which was consistent with neutrino data before the discovery of non-zero θ 13 can be realised naturally within several flavour symmetric models based on non-Abelian discrete groups [25][26][27][28]. Among such flavour symmetric models, the discrete group A 4 , group of even permutations of four objects, can reproduce the TBM mixing in the most economical way [29][30][31][32][33][34][35][36][37][38][39]. However, in order to be consistent with the present neutrino data, such TBM or µ − τ symmetric scenarios have to be corrected to generate non-zero value of θ 13 . There have been several attempts in that direction, some of which can be found in  and references therein.
Apart from non-zero neutrino mass and large leptonic mixing, another observed phenomena that has propelled serious hunt for beyond standard model (BSM) physics is the presence of non-baryonic form of matter, or the so called dark matter (DM) in large amount in the present universe. Apart from the longstanding astrophysical evidences [63][64][65], the recent cosmology experiment Planck suggests that almost 26% of the present Universe's energy density is in the form of DM while only around 5% is the usual baryonic matter leading the rest of the energy budget to mysterious dark energy [13]. In terms of density parameter Ω DM and h = Hubble Parameter/(100 km s −1 Mpc −1 ), the present DM abundance is conventionally reported as [13]: Ω DM h 2 = 0.120 ± 0.001 at 68% CL. Since none of the SM particles can satisfy the requirements for being a DM candidate, several BSM proposals have been put forward out of which the weakly interacting massive particle (WIMP) paradigm is perhaps the most popular or the most widely studied one. The coincidence that a stable or cosmologically long lived particle having electroweak scale mass and interactions fulfilling the criteria for observed DM abundance is often referred to as the WIMP Miracle. However, if such particles having electroweak scale mass and interactions really exist in the universe with such a large density, they are expected to pass through the detectors of several DM direct detection experiments giving rise to nuclear recoils. However, there have been no detection of particle DM at any experiments. The direct detection experiments like LUX [66], PandaX-II [67,68] and Xenon1T [69,70] have continued to produce null results so far.
Similar null results follow from collider searches at the large hadron collider (LHC) [71] as well as the indirect detection frontiers [72]. The typical indirect detection experiments excess of antimatter, gamma rays or neutrinos, originating perhaps from dark matter annihilations (for stable DM) or decay (for long lived DM) and no convincing signal has been observed at any experiment operating in this frontier. Such null results from WIMP searches have led to the proposals of several alternative frameworks of DM, specially to scenarios where the interaction between DM and visible sector could be much more weaker than what it is in the WIMP paradigm. In fact, if the coupling of DM to visible sector is sufficiently weak, then DM can never be produced thermally in the universe requiring a non-thermal origin for its relic [73]. In such scenarios, the DM has negligible initial number density in the early universe and is produced from out of equilibrium decay or scattering of visible sector particles. The production mechanism for non-thermal DM is known as freeze-in and the candidates of non-thermal DM produced via freeze-in are often classified into a group called Freeze-in (Feebly interacting) massive particle (FIMP). For a recent review of this DM paradigm, please see [74]. The tiny couplings between DM and visible sector can be naturally realised either by higher dimensional operators [73,75,76] or through some UV complete renormalisable theories [77]. While typical FIMP does not have much direct detection prospects and WIMP direct detection has so far failed, there exists a scenario which is a combination of both and have interesting detection prospects in terms of secondary particles. This is known as super-WIMP scenario [78] where a metastable WIMP decays into a super-weakly interacting dark matter at late epochs. Although DM still has feeble interactions with the visible sector, the metastable WIMP has sizeable interactions and can be detected as, for example, long lived BSM particle at collider experiments [79]. This scenario was also adopted in the context of neutrino mass models in several works including [80][81][82]. In most of such models there exists some exact or approximate discrete symmetries to stabilise the DM or to give rise to a long lived metastable WIMP.
In this work, we try to find a common thread linking the discrete symmetries in neutrino sector and DM sector. We consider the neutrino sector to have an exact µ − τ symmetry while the dark sector has an exact Z 2 symmetry. We find a common source of breaking of these discrete symmetries through Planck scale suppressed operators by using the well known argument that any generic theories of quantum gravity should not respect global symmetries: both discrete and continuous [83][84][85]. Recently, this argument was used to realise light neutrino mass from Planck scale lepton number breaking [86]. Effects of such Planck scale breaking of discrete symmetries on light neutrino parameters was studied within the context of left right symmetric model a few years back in [87]. Here we study the consequences of such Planck scale breaking of discrete symmetries on neutrino sector which generates non-zero θ 13 , as required by present neutrino data. At the same time, such breaking also makes it possible to realise the super-WIMP scenario by assisting the decay of a metastable WIMP to a super-weakly interacting dark matter sector. We consider the presence of dis- We however, restrict ourselves to discussing the origin of non-zero θ 13 in this work, without assuming a global lepton number symmetry at renormalisable level of the theory. This paper is organised as follows. In section II we discuss our flavour symmetric model and possible Planck suppressed operators. In section III we discuss neutrino mass and mixing followed by super-WIMP dark matter phenomenology in section IV. We finally conclude in section V.

II. A 4 MODEL WITH TBM MIXING
In this section we outline the flavour model based on non-Abelian discrete group A 4 augmented by a Z 4 symmetry for dark sector. A brief summary of A 4 group, its representations and product rules are given in appendix A. The particle content of the model, relevant for discussing the lepton sector and dark matter sector is shown in table I. Apart from the SM lepton doublets and charged lepton singlets, there are three right handed neutrinos required for implementing the type I seesaw mechanism [15][16][17][18]. In addition to the usual SM Higgs doublet, there exist three additional Higgs doublets required to generate Dirac mass term for neutrinos. There also exists two more scalar fields responsible for generating the desired right handed neutrino mass matrix and also to break the A 4 flavour symmetry spontaneously. Two more fields one scalar and one fermion, both charged under the Z 4 symmetry are included in order to achieve the desired DM phenomenology.
Fields L e , L µ , L τ e R , µ R , τ R H H T φ N ξ N R ψ η A 4 1, 1 , 1 1,1 ,1 1 3 3 1 3 1 1 whereH T = τ 2 H * T . Due to the additional Z 4 symmetry, the fields ψ, η do not have any renormalisable coupling with SM leptons or right handed neutrinos. We assume ψ be vector like so that a bare mass term M ψψ ψ is allowed in the Lagrangian. As mentioned earlier, gravity is not supposed to obey such global symmetries and hence we can write down Planck suppressed terms which explicitly break both A 4 and Z 4 . The dimension five terms involving one of the A 4 triplet flavon φ N which explicitly break the discrete symmetries but preserve the gauge symmetries of the standard model can be written as We can also write down the dimension five terms simultaneously involving two of the A 4 triplet flavons H T , φ N which explicitly break the discrete symmetries but preserve the gauge symmetries of the standard model. They can be written as One can prevent the coupling of η, ψ with standard model leptons at dimension five level by considering an additional U (1) X gauge symmetry under which ψ is a Dirac fermion with charge n 1 while η has a charge 2n 1 . Thus, η can behave like a next to lightest particle which decays only to ψ due to the Planck scale effects on Z 2 symmetry breaking.
As discussed below, we choose a hierarchy of flavon vacuum expectation value (VEV)'s in order to achieve the desired phenomenology. While the SU (2) L singlet flavon VEV's can be arbitrary, the same does not apply to that of the triplet flavon which is also a doublet under SU (2) L . In order to be in agreement with precesion electroweak measurements [10], the vev of the neutral components of A 4 singlet Higgs doublet H(v 1 ) and A 4 triplet Higgs GeV. The scalar potential involving H, H T only can be written as which is similar to the scalar potential of a multi Higgs doublet model with additional and detection prospects of the components of H T will be similar to two or multi Higgs doublet models [90,91] as the physical masses of all the components can be kept around the TeV scale.

III. LEPTON MASSES AND MIXINGS
Up to leading order, owing to the A 4 symmetry, one can write the charged lepton mass matrix in the diagonal form as where v 1 is the VEV of the neutral component of the Higgs doublet H. We use the A 4 product rules given in appendix A in order to evaluate the forms of different mass matrices at renormalisable level. Considering the vacuum alignment of A 4 triplet scalar field H T = (v 2 , 0, 0), the Dirac neutrino mass matrix at leading order can be written as, where, without any loss of generality, one can consider Y ν1 = Y ν2 = Y ν3 .
Considering triplet flavon vacuum alignment as φ N = u(1, 1, 1) and ξ = u ξ the right handed neutrino mass matrix at leading order is given by where a = 2Y ξ u ξ and b = 2Y N u. Hence the light neutrino mass matrix, using the type I seesaw formula is This µ − τ symmetric light neutrino mass matrix can be diagonalised by the usual tribimaximal mixing matrix given by [92] predicting θ 13 = 0 and maximal value of θ 23 . Present experimental observation, however, rules out θ 13 = 0 scenario and has a preference towards higher octant for θ 23 (i.e.θ 23 > 45 • ) [12]. We therefore consider the Planck scale suppressed corrections to such scenarios so that the correct light neutrino phenomenology can be obtained along with the implications for dark sector.
The contribution to the light neutrino mass matrix and neutrino mixing from the Planck suppressed terms which break the discrete symmetries explicitly can arise in a variety of ways. For example, • there can be contributions to light neutrino mass matrix of the type • there can be new contributions to charged lepton mass matrix of the type in principle, lead to deviations from diagonal charged lepton mass matrix. Usually, the leptonic mixing matrix is given in terms of the charged lepton diagonalising matrix (U l ) and light neutrino diagonalising matrix U ν as U = U † l U ν . In the simple case where the charged lepton mass matrix is diagonal which is true in our model at tree level, we can have U l = 1. Therefore we can write U = U ν . But U l may become non-trivial after the corrections are added and it will mimic as a correction in the leading order leptonic mixing which is TBM type.
• there can be corrections to the Dirac neutrino mass matrix of the type which can again be suppressed by an additional factor u/M Pl compared to the leading order contributions. Such corrections will propagate to the light neutrino mass matrix through the type I seesaw formula.
• there can be corrections to the right handed neutrino mass matrix of the type (Y ξ ξ 2 + Y N φ 2 N )N R N R /M Pl which can change the structure of M R from the µ − τ symmetric leading order form mentioned before. Once again this correction will propagate to the light neutrino mass matrix via type I seesaw formula.
All these corrections can be simultaneously or separately sufficient to generate correct neutrino mixing. For representative purpose, we will consider corrections to the Dirac neutrino mass matrix only and that too via Higgs doublet H. Addition of other corrections will make the calculations complicated without any new insights. A more general analysis is beyond the scope of this work and left for future studies. Now, following equation (2), the correction to the Dirac neutrino mass matrix now can be written as, In the present construction, following A 4 multiplication rules given in the appendix, we find With such correction the effective Dirac neutrino mass matrix can be written as (assuming where . Now, along with the Planck suppressed operator contribution in Dirac neutrino mass, the light neutrino mass originating from type I seesaw can be written as After a rotation by U TBM , the light neutrino mass matrix can be written as where Clearly an additional rotation in the 13 plane is required to diagonalise the above matrix via the relation with where m 1 , m 2 , m 3 are the real positive neutrino mass eigenvalues and γ 1,2,3 are the respective Majorana phases. Therefore the effective neutrino mixing matrix (with diagonal charged lepton sector), can be written as where U m is the diagonal Majorana phase matrix given by U m = diag(e iγ 1 , e iγ 2 , e iγ 3 ). Comparing this with the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix one can obtain the correlation between the model parameters and neutrino oscillation parameters. In the above PMNS mixing matrix, P (= diag(1, e iα 21 /2 , e iα 31 /2 )) is the Majorana phase matrix where α 21 = (γ 1 − γ 2 ) and α 31 = (γ 1 − γ 3 ) are the physical Majorana phases after rotating away one common phase which is irrelevant. The rotation parameter θ in the rotation matrix (19) can be obtained as with ψ = 0, π respectively yielding maximal value for the Dirac CP phase δ. The construction of the present set-up is such that we only have considered the Planck scale suppressed terms only for the Dirac neutrino mass matrix. Once we allow correction terms for the charged lepton and/or right handed neutrino mass matrix, non-maximal values for the Dirac CP phase δ becomes allowed [54,93,94]. Comparing the neutrino mixing matrix of our scenario given in (20) with the standard parametric form of mixing matrix given in (21), the neutrino mixing angles can be evaluated as [54] sin 2 θ 13 = 2 3 | sin θ| 2 (23) Thus, the neutrino mixing angles are function of the angle θ which depends upon x 1,2 and they further depend upon the coupling constants y 1,2,3 as described before. As mentioned earlier, in the present set-up we have considered the deviation from TBM mixing (which usually predicts sin 2 θ 12 = 1/3, sin 2 θ 23 = 1/2 and θ 13 = 0) via Planck scale suppressed operators appearing solely in the neutrino sector. The present construction (based on A 4 discrete symmetry) is such that observed θ 13 can be successfully generated in presence of such operators and deviations in θ 12 and θ 23 from their respective TBM values also arise at the same time. Such deviation from TBM mixing can be obtained via a unitary rotation (in the 13 plane) matrix, U 13 , parameterised only by two parameters θ and ψ which also satisfy the condition |(U 13 ) 11 | 2 + |(U 13 ) 13 | 2 = 1. The correlations among θ, ψ and the mixing angles are generic features of this class of models [42,53,95] considered here and can be realised very easily using A 4 discrete symmetry [54,93,94,96]. However in contrast to all other models mentioned above, the present scenario predicts the Dirac CP phase δ to be maximal (as ψ acquires values 0 and π). Such predictions naturally make the model testable as well as distinguishable from other scenarios. Now, the neutrino mass eigenvalues can be written as where α = |b|/|a| and φ ba = φ b −φ a is the phase difference between b and a. Using these, one can now define a ratio of solar to atmospheric mass-squared differences (∆m 2 = ∆m 2 21 = m 2 2 − m 2 1 and |∆m 2 Here, depending on the sign of denominator, the parameter space is divided into two parts: '+' in the denominator yields normal hierarchy of neutrino mass whereas for '-' in the denominator yields inverted neutrino mass hierarchy. Therefore, the neutrino mixing angles and the ratio of solar to atmospheric mass-squared differences are altogether functions of five parameters, namely three coupling constants y 1 , y 2 , y 3 and α, φ ba . Using the observed neutrino oscillation data summarised in latest global fit [12], these parameters can be constrained for both normal and inverted neutrino mass hierarchy, as we discuss below. A. Normal Hierarchy Here we first explore neutrino phenomenology for normal neutrino mass hierarchy. In left panel (right panel) of figure 1 we have shown the allowed points for y 1 vs y 2 (y 3 ).
All these points satisfy 3σ allowed range for the three neutrino mixing angles, ratio of solar to atmospheric mass-squared differences [12]. From the equation (26) and equation (27), one can easily obtain the common factor |a| appearing in the absolute neutrino mass using ∆m 2 21 = m 2 2 − m 2 1 = 2.43 × 10 −5 eV 2 . Therefore the parameter space can be further scale should be close to GUT scale. As far as fine-tuning is concerned, we set a tolerance same as electron Yukawa in the standard model. Different A 4 breaking scale will also have different implications for the super-WIMP dark matter sector, requiring some amount of fine tuning in relevant Yukawa couplings related to mother particle and dark matter, as we discuss below.

IV. SUPER-WIMP DARK MATTER
The dark matter in our setup is similar to the super-WIMP scenario [78] where a metastable WIMP decays into super-weakly interacting dark matter at late epochs. As can be seen from the dimension five Lagrangian (2), there are several terms which can give rise to decay of either ψ or η, the particles of the Z 4 sector. At renormalisable level however, both of them are stable due to the unbroken Z 4 symmetry. While η can be produced thermally in the early universe due to its electroweak gauge interactions, the singlet fermion ψ has negligible thermal abundance. In a general super-WIMP formalism where η is the metastable WIMP and ψ is the DM candidate, one can write down the the decay width of η into two dark matter particles (ψ) as where Y ηψ is the effective Yukawa coupling, m η and m ψ are the masses of the mother particle η and ψ respectively.
However, this doublet scalar, apart from having electroweak gauge interactions, can have sizeable quartic interactions with other scalars like the standard model Higgs doublet and hence can be thermally produced in the early universe. Now, considering the mother particle η to be in thermal equilibrium in the early universe which also decays into the dark matter particle ψ, we can write down the relevant Boltzmann equations for comoving number densities of η, ψ as where x = M sc T , is a dimensionless variable while M sc is some arbitrary mass scale which we choose equal to the mass of η and M Pl is the Planck mass. Moreover, g s (x) is the number of effective degrees of freedom associated to the entropy density of the universe and the quantity g (x) is defined as Here, g ρ (x) denotes the effective number of degrees of freedom related to the energy density of the universe at x = M sc T . The first term on the right hand side of the Boltzmann equation (31) corresponds to the self annihilation of η into standard model particles and vice versa which play the role in its freeze-out. The second term on the right hand side of this equation corresponds to the dilution of η due to its decay into dark matter ψ. Let us denote the freezeout temperature of η as T F and its decay temperature as T D . If we assume that the mother particle freezes out first followed by its decay into dark matter particles, we can consider In such a case, we can first solve the Boltzmann equation for η considering only the self-annihilation part to calculate its freeze-out abundance.
Then we solve the following two equations for temperature T < T F We stick to this simplified assumption T F > T D in this work and postpone a more general analysis without any assumption to an upcoming work. The assumption T F > T D allows us to solve the Boltzmann equation (34) Since flavon VEV is much larger than electroweak one, the first two decay modes will be more dominant than the last one. Out of these two, the mode η → lψ will contribute to the abundance of DM ψ. One can also think of some non-minimal 2 Recently another scenario was proposed where the dark matter freezes out first with underproduced freeze-out abundance followed by the decay of a long lived particle into dark matter, filling the deficit [80,97].
scenarios where there exists some local symmetries which forbids η → ll or even η → lψ depending upon how η, ψ are charged under the local symmetry. In the former case, DM will be produced dominantly from η → lψ mode while in the second case the only possible mode that remains is the η → ψψ one. Since these two modes have very different effective couplings, the DM phenomenology can be very different. We now discuss all these three distinguishing cases one by one.
We note that our second and third cases will involve some local symmetries which are not explicitly broken by Planck scale physics. Here we do not outline a complete model for such cases but point out the distinguishing DM phenomenology without affecting the neutrino phenomenology. The additional gauge interactions of ψ, η which can prevent coupling of ψ, η to SM leptons can also produce ψ thermally in the early universe, in principle. However, if such mediator gauge bosons are very heavy, with masses greater than the reheat temperature after inflation, then the thermal production of such DM remains suppressed [99] and non thermal contribution from the metastable WIMP will become more relevant. Another way is to assume very feeble gauge interactions so that they never get thermalised.

A. Case I
In this subsection we discuss the DM phenomenology in the most general way, corresponding to the symmetries and particle content of the A 4 × Z 4 model discussed before. As where, λ 1 and λ 2 can be written as However in this case, DM is also not perfectly stable as the Z 4 symmetry which protects its stability gets explicitly broken by Planck suppressed terms. As can be seen from (2), ψ can decay into the Higgs boson and leptons due to the term Y 1 L αH ψ. To forbid this decay at tree level, we consider the DM mass to be below 1 MeV. This term also can give rise to a mixing of ψ with light neutrinos of the order which can be tuned appropriately in order to satisfy X-ray or gamma ray bounds.
In this case, the metastable WIMP η not only decays to DM but also the charged lepton pairs. Though DM remains out of thermal equilibrium throughout, the charged leptons are part of the thermal bath and hence producing them from η decay can, in principle, release entropy. In that case, solving the coupled Boltzmann equations for η and ψ described earlier is insufficient and one has to consider a third equation for radiation energy density of the universe. Before going to the actual DM calculations, we first check the amount of entropy release by solving the following three coupled Boltzmann equations namely, where Γ η = Γ η→l + l − + Γ η→ψν . Please note that here we are writing the Boltzmann equations in terms of ordinary number density (n), energy density (ρ) and time (t) unlike the earlier equations in terms of comoving quantities. The first equation above shows the time evolution of the number density of η where the usual dilution due to the expansion of the universe is given by the first term on the right hand side (RHS) and dilution due to the decay is given by the second term on the RHS. There should, in principle, be one more term on RHS namely, − σv ηη→pp n 2 η 0 − (n eq η 0 ) 2 with p denoting any SM particles to which η can self annihilate into. However, we have dropped this term assuming that η has already frozen out before it starts decaying which is a reasonable assumption (like super-WIMP framework) for the couplings governing η decay in our model. The second and third equations are for the evolution of radiation energy density and the DM number density respectively. The figure 7 shows that the contribution to the radiation or entropy release due to the decay of η to the lepton pairs is negligible all the way up to the epochs where η decays to reduce its number density while yielding ψ. The radiation energy density remains dominant during this and does not get any significant contribution from η decay. Since the entropy release is negligible, it thereby justifies the use of coupled Boltzmann equations only for metastable WIMP and DM in our analysis. We now write down the coupled Boltzmann equations for η, ψ in terms of their comoving number densities as We solve these two equations (41) and (42) figure 8, we vary the mother particle's mass while keeping DM mass and mother particle's freeze-out abundance fixed. We see that the final abundance of DM does not depend much upon mother particle's mass as long as the freezeout abundance remains fixed. From these plots it is clear that the final abundance of DM Ω DM h 2 = m DM Y DM s 0 strictly depends on DM mass whereas it is almost independent of the mass of the mother particle. This is due to the fact that η has the same branching ratio to two different decay channels (charged lepton pairs and DM plus lepton) therefore half of it's abundance going into that of DM irrespective of its mass. We summarise our results for case I as a scan plot shown in figure 9. The plot shows freeze-out abundance of η and mass of η from two different directions and requirements. The orange-blue coloured band comes from the relic density requirement of DM. Since DM is being produced from η decay after thermal freeze-out of η, the orange-blue coloured points give us the freeze-out abundance of η and its mass so that DM with a particular mass in that coloured band has the correct relic abundance. For this plot, the effective Yukawa coupling between η and DM has been kept fixed at 10 −12 . As can be seen, the dependence on η mass is marginal which was also observed in the benchmark plots shown in figure 8. This is expected as DM abundance should strongly depend upon the freeze-out abundance of η as well as DM mass. We then check whether the desired freeze-out abundance of η can actually be obtained for some choices of η parameters. We can write down the components of η as Here we consider η R as the lightest component of η. Calculating the thermal freeze-out abundance of η R is similar to the DM relic calculation in the inert doublet model discussed extensively in2 the literature [100][101][102]. Typically there exists two distinct mass regions, M η ≤ 80 GeV and M η ≥ 500 GeV, where correct relic abundance criteria can be satisfied. In both regions, depending on the mass differences m η ± −m η R , m η I −m η R , the coannihilations of η R , η ± and η R , η I can also contribute to the DM relic abundance [103,104]. The parameters which crucially affect the thermal freeze-out abundance of η are the mass splitting ∆M η within the components of η and η couplings with the SM Higgs λ L . It should be noted that there are four different Higgs doublets which can mediate η interactions with the SM particles. Here, for simplicity, we consider the scalar interactions to be mediated only by the SM like Higgs, governed by the coupling λ L . We fix the mass splitting as ∆M η = 50 GeV and scan over λ L to find the thermal freeze-out abundance of η R . The corresponding region of parameter space is shown as red-black colour coded region of figure 9. The overlapping region of orange-blue and red-black colour coded regions of figure 9 contain the points which satisfy the correct DM abundance criteria in our model. As can be seen from the colour codes, smaller the DM mass, more overproduced η R has to be, as expected. And as it is well known for freeze-out of scalar doublet, it is the low mass regime of η R where it is more likely to be overproduced due to the absence of sufficient annihilation channels resulting from kinematic suppression as well as departure from s-channel resonance. One can also have such overproduced regime of η R in the high mass limit m η ≥ 1 TeV which we have not shown in this plot.
It should be noted that, there exists strong bounds on the masses of different components of η as well as their couplings with the SM Higgs. For example, precision measurements of W, Z decay widths at LEP experiment lead to the bounds The region defined by the intersection of the following conditions m η R < 80 GeV, m η I < 100 GeV, m η I − m η R > 8 GeV is also excluded from the non-observation of dijet or dilepton signals at LEP [105,106].
The charged scalars are constrained as m η ± > 70 GeV by reinterpreting the LEP bounds on charginos. We apply these bounds here simply by taking a conservative mass splitting figure 9. Another bound comes from the LHC from the measurements of invisible decay width of the SM Higgs as well as SM Higgs decay into diphoton [106]. We check that the range of parameters we scan through satisfy these bounds.

B. Case II
In this case η can decay only through the interaction φ N i M P l ηνψ as the other interaction is forbidden here. As there is no decay terms present for ψ, it is not required to consider sub-MeV mass of ψ for forbid its decay kinematically. We consider DM mass in the GeV regime here. The relevant coupled Boltzmann equations in this case are almost same as the ones discussed in previous subsection namely, (41) and (42) except the fact that Γ η is now replaced by Γ η→ψν as there is no other decay channel present for η. They can be written as In figure 10 we have shown the variation of DM abundance as a function of (x = mη T ) for different benchmark values of model parameters. The left panel of upper row in figure 10 shows the variation caused by different freeze-out abundance of mother particle while the right panel plot shows the variation with DM mass. In the lower panel plot of figure 10 we show the variation due to change in mother particle's mass. As in case I, the final abundance of DM does not depend much upon mother particle's mass as long as its freeze-out abundance is kept fixed. One interesting feature observed in right panel plot in the upper row of this figure (not noticed in case I) is the change in evolution of DM density as DM mass becomes closer to mother particle's mass. Final abundance of DM is always proportional to DM mass but when m DM becomes very close to the m η then η decays slowly and that can be seen from the brown line corresponding to m DM = 90GeV which grows slower with x compared to the line corresponding to lower values of DM mass m DM = 10GeV . Similar to case I, here also we make a parameter scan that can give rise to the correct relic abundance of super-WIMP DM. The plot is shown in figure 11 where the triangles correspond to DM masses in the white-red colour code while dots correspond to ∆M η with blue-red colour code. The triangular points correspond to the parameter space that gives rise to correct DM abundance for a particular set of (M DM , M η , Ω η h 2 ). As before, the mass of η does not play much role on DM abundance. However, ∆M η plays crucial role in generating the required freeze-out abundance of η. As can be seen from this plot, there exists very small overlaps between triangles and dots which correspond to correct DM abundance as well as realistic freeze-out abundance of mother particle for particular choice of ∆M η . Thus, case II gets more constrained compared to case I in terms of final parameter space.

C. Case III
We finally consider the case where the mother particle can decay into DM only through interactions of typeỸ 1 H M Pl η † ψψ. This differs from case II by the facts (i) mother particle decays into two DM particles instead of one, (ii) the effective strength of the vertex η −ψ −ψ is much weaker Y eff ≈Ỹ 1 10 −17 . The relevant Boltzmann equations are where the decay width Γ η→ψψ is given by equation (30). The corresponding results are shown in figure 12 where we show the evolution of DM relic with temperature for different benchmark choices of model parameters. The overall evolution remains similar to that in case II (shown in figure 10) except for the fact that due to smaller effective Yukawa coupling between mother particle and DM, the yield in DM abundance happens at relatively lower temperatures or higher x. The parameter space scan for case III remains very similar to that for case II (shown in figure 11) and we skip it showing again.
It should be noted that in case I and II, we considered the effective Yukawa coupling between η and ψ to be 10 −12 which will require some fine tuning in the Yukawa couplings appearing with the Planck suppressed operators if the scale of A 4 symmetry that is, u lies close to GUT scale. In case III, the effective Yukawa coupling is very small 10 −16 which arises naturally from the ratio of electroweak scale to Planck scale. In case III, the effective Yukawa coupling of super-WIMP is independent of A 4 breaking scale and requires no finetuning. The smallness of effective Yukawa couplings chosen here also justifies the validity of super-WIMP formalism where the major non-thermal contribution to DM relic comes after metastable WIMP freezes out. In general non-thermal or FIMP dark matter scenario, contribution to DM can arise while mother particle is in thermal equilibrium as well as after the mother particle freezes out [73]. Equilibrium contribution dominates for larger Yukawa coupling Y eff ∼ 10 −9 − 10 −8 and for chosen Yukawa couplings here, the equilibrium contribution remains suppressed below 5%. This also has helped us to solve the coupled PLANCK 2018 FIG. 13: Lifetime of mother particle η as function of effective Yukawa coupling Y eff for its two body decay into dark matter.
Boltzmann equations at two stages: first solving only the mother particle's equation for its thermal freeze-out and then solving the coupled Boltzmann equations for mother particle and DM as discussed above.
Finally, we check whether the three scenarios discussed in our work could affect big bang  [110,111]. However a detailed study in this direction is beyond the scope of present study. We also outline the super-WIMP dark matter phenomenology by considering three distinct scenario: (i) η decays to SM particles as well as ψ and ψ is kinematically long lived, (ii) η decays to ψ and a SM neutrino while ψ is perfectly stable, (iii) η decays into a pair of ψ while ψ is perfectly stable. Out of these, the first scenario correspond to the model which we have discussed in our work while the latter two scenario can be realised if the discrete Z 4 symmetry in the dark sector can be uplifted to a gauge symmetry which does not get broken by gravity effects. While we do not discuss such UV complete gauge symmetric realisation of Z 4 symmetry we outline the interesting differences for super-WIMP phenomenology. The analysis for neutrino sector in all three DM scenarios remain same, however.
We show that correct dark matter relic abundance can be obtained in all three distinct scenarios. While direct detection of super-WIMP dark matter itself may not be very optimistic, the mother particle η can be probed at ongoing experiments. Since the lightest component of η decays to DM as well as SM leptons (in case I) with very feeble couplings one can probe them in colliders either as missing energy or displaced vertex, long charged tracks depending upon the lifetime. If charged component of η is the lightest component of η then for typical super-WIMP couplings discussed in this work, it will give rise to a long charged track in colliders as its decay length will be much larger than typical displaced vertex ones searched for. For more discussions on displaced vertex and disappearing charged track signatures of a similar model, please refer to a recent work [112]. Apart from signatures of mother particle, the DM itself can leave some detectable signatures, specially in case I where it is not perfectly stable but long lived. For example, a 7 keV long-lived fermion DM 3 can decay into a photon and light neutrino at radiative level with W boson and charged leptons of the SM in loop. The corresponding decay width is given by [114] Γ(ψ → νγ) ≈ 1.38 × 10 −29 s −1 sin 2 2θ 1 × 10 −7 where θ denotes the mixing between N 1 and ν. From the observation of the 3.55 keV line, which can arise from the decay of a 7.1 keV sterile neutrino DM, the mixing angle which is in agreement with the observed flux is sin 2 2θ ≈ 7 × 10 −11 [115]. Such a mixing angle can be naturally obtained in the model, as seen from the expression for mixing angle given in equation (39). Although the analysis of the preliminary data collected by the Hitomi satellite (before its unfortunate crash) do not confirm such a monochromatic line [116], one still needs to wait for a more sensitive observation with future experiments to have a final word on it. We leave a more detailed study of detection prospects as well as UV complete realisations of case II, III to future works. The non-Abelian discrete group A 4 , also the symmetry group of a tetrahedron, is a group of even permutations of four objects. This group has four irreducible representations out of which three are one-dimensional and one three dimensional, denoted by 1, 1 , 1 and 3 respectively, being consistent with the sum of square of the dimensions i n 2 i = 12. We denote a generic permutation (1, 2, 3, 4) → (n 1 , n 2 , n 3 , n 4 ) simply by (n 1 n 2 n 3 n 4 ). The group A 4 can be generated by two basic permutations S and T given by S = (4321), T = (2314).

This satisfies
which is called a presentation of the group. Their product rules of the irreducible representations are given as where a and s in the subscript corresponds to anti-symmetric and symmetric parts respectively. Denoting two triplets as (a 1 , b 1 , c 1 ) and (a 2 , b 2 , c 2 ) respectively, their direct product can be decomposed into the direct sum mentioned above. In the S diagonal basis, the products are given as 1 a 1 a 2 + b 1 b 2 + c 1 c 2 1 a 1 a 2 + ω 2 b 1 b 2 + ωc 1 c 2 1 a 1 a 2 + ωb 1 b 2 + ω 2 c 1 c 2 3 s (b 1 c 2 + c 1 b 2 , c 1 a 2 + a 1 c 2 , a 1 b 2 + b 1 a 2 ) In the T diagonal basis on the other hand, they can be written as 1 a 1 a 2 + b 1 c 2 + c 1 b 2