Low Scale $U(1)_X$ Gauge Symmetry as an Origin of Dark Matter, Neutrino Mass and Flavour Anomalies

We study a generic leptophilic $U(1)_X$ extension of the standard model with a light gauge boson. The $U(1)_X$ charge assignments for the leptons are guided by lepton universality violating (LUV) observables in semileptonic $b \to s\ell\ell$ decays, muon anomalous magnetic moment and the origin of leptonic masses and mixing. Anomaly cancellation conditions require the addition of new chiral fermions in the model, one of which acts as a dark matter (DM) candidate when it is stabilised by an additional $\mathcal{Z}_2$ symmetry. From our analysis, we show two different possible models with similar particle content that lead to quite contrasting neutrino mass origin and other phenomenology. The proposed models also have the potential to address the anomalous results in $b\to c\ell\nu_{\ell}$ decays like $R(D), R(D^*)$, electron anomalous magnetic moment and the very recent KOTO anomaly in the kaon sector. We also discuss different possible collider signatures of our models which can be tested in future.


I. INTRODUCTION
The large hadron collider (LHC) at CERN is operational for more than last ten years and so far apart from the discovery of Higgs boson, no new particles or interactions have been found. No evidence for the theoretically wellmotivated models like supersymmetry, extra dimension etc have been found. Yet there is a list of unsolved puzzles in particle physics. In the standard model (SM) of particle physics, we do not have explanations for neutrino masses, the existence of dark matter (DM) and the domination of matter over antimatter in the Universe [1]. Nature may still be supersymmetric, or there may be an extra dimension; however, these extensions of SM have failed to show up at the LHC. Even if they are absent, there are a lot of things to learn, at a fundamental level. In principle, one could write down models which are consistent with the present observations at the collider and will show distinct features only at the high luminosity, as for example see [2]. There are indications that we have not yet fully understood the working rules of our Universe at the fundamental level, and that is motivating enough for the particle physics community to keep looking for it. In addition to the ongoing LHC experiment, we already have a few experimental facilities which are operational or will start functioning very soon, and very quickly data will be collected in unprecedented amounts. We can hope that the upcoming data will guide us to establish the more fundamental theory of elementary particles and their interactions.
Apart from the direct searches at the collider, the low energy observables play an essential role for indirect detection of a new particle(s) or interaction(s). In this regard, B-factories have played a significant role in the last couple of decades [3] and will remain productive in near future [4]. In the last couple of years, LHCb has also produced significant results, for a brief review see [5,6]. In the low energy data, new physics (NP) contributions in an observable can be pinpointed through the deviation of its measured value from the respective SM prediction. At the moment there are a few measurements in b → c and b → s decays which show some degree of discrepancies with their respective SM predictions, for very recent updates see [7,8]. Apart from these long standing anomalies, more recently an excess of events have been observed in the rare K L → π 0 νν decay (d → s FCNC process) at the KOTO experiment at J-PARC [9].
The measurements of various angular observables in B → K * µ + µ − [10,11] and B s → φµ + µ − [12] decays are available, and in a few of them there are discrepancies between the theory and experiment. Very recently, with the data collected by the LHCb experiment during the years 2011, 2012 and 2016, a complete set of CP-averaged angular observables has been measured in B → K ( * ) µ + µ − decay [13]. To date, this is the most precise measurement, and the data still shows discrepancies between the theoretical predictions and the measured value in a couple of those angular observables. Note that these angular observables are not free from hadronic uncertainties. However, there are theoretically clean observables like R(K ( * ) ) = B(B→K ( * ) µ + µ − ) B(B→K ( * ) e + e − ) the measured values of which [14,15] are not in good agreement with the corresponding SM expectations. There are new physics explanations of these observations, for a recent update on the model-independent new physics explanation of these data see [16][17][18] and the references therein.
Similar to the observables R(K ( * ) ), we define R(D ( * ) ) = B(B→D ( * ) τ ντ ) B(B→D ( * ) ν ) (with = µ, e) which is associated with the b → c decays. The measured values of these observables [7] have also shown some degree of discrepancies with the respective SM predictions, for details see [19][20][21]. The most recent predictions (in SM) differ from the one obtained using the old Belle data [22,23]. The bounds on the model-independent new physics Wilson-coefficients (WC) can be seen from [21]. It is found that the data still allows sizeable NP contributions in these decays.
Apart from the above mentioned results, the muon anomalous magnetic moment (g − 2) µ is another longstanding puzzle. It has been measured very precisely while it has also been predicted in the SM to a great accuracy. At present the difference between the predicted and the measured value is given by ∆a µ = a exp µ − a SM µ = 26.1(7.9) × 10 −10 , which shows there is still room for NP beyond the SM (for details see [1]). In a recent article, the status of the SM calculation of muon magnetic moment has been updated [24]. According to this study, the difference is given by ∆a µ = 27.9(7.6) × 10 −10 , which is a 3.7σ discrepancy. Analogous to muon magnetic moment, measurements are also available for electron magnetic moment (g − 2) e . The most recent result obtained from measurement of the fine structure constant of QED [25], shows a deviation from the SM. The excess is given by ∆a e = −8.7(3.6) × 10 −13 .
In this study, we look for a NP model which is capable of addressing all the above-mentioned results. At first, we consider a simple model by extending the SM gauge group with an additional U (1) X gauge symmetry 1 . The resulting complete gauge group of the model will be SU (3) c × SU (2) L × U (1) Y × U (1) X which is an extension of SM by an abelian factor. The advantage of such an extension is that it introduces a minimal set of free parameters. The other most important feature of the new gauge symmetry we adopt here is that it is leptophilic in nature i.e. only the leptons will be charged under U (1) X , not the quarks. For an explanation of the above mentioned anomalous results, the lepton generations must have different charges under U (1) X . The degree of fermion non-universality should explain the observed discrepancies in R(K ( * ) ) and muon anomalous magnetic moment. In this minimal model with GeV scale mass of U (1) X gauge boson, we can not explain R(D * ) and the data on electron anomalous magnetic moment.
However, charging the fermions under this new gauge group in the absence of additional chiral fermions generally leads to triangle anomalies which must be cancelled in order to validate the gauge theory at the quantum level. Hence, in order to cancel the gauge anomalies, we need to introduce additional degrees of freedom into our model, in terms of chiral fermions. Here, following the constraints from gauge anomaly cancellation, we discuss only two different possible scenarios in which we can explain the existing data on DM and neutrino oscillation. In extended version of such minimal model with more particle and interactions, there will be additional Feynman diagrams which will contribute to R(D * ) and ∆a e that help us to explain the observed data.
In a similar direction, studies are available in the literature with a heavy U (1) X gauge boson [27][28][29][30][31][32][33]. While such models with heavy U (1) X gauge boson have been extensively studied, there have been very few studies on low mass regions [34][35][36]. However, our working model is very much different compared to the one discussed in the references mentioned above and we also correlate the flavour anomalies with origin of neutrino mass and dark matter. Both the scenarios we discuss here consider the viability of a leptophilic U (1) X gauge symmetry in a way that it is anomaly free, predicts lepton flavour non-universality and the origin of light neutrino masses while the stability of DM candidate is ensured by an additional Z 2 symmetry which also plays a non-trivial role in neutrino mass generation for one of the models.
This paper is organised as follows. In section II we briefly discuss our overall framework followed by the corresponding analysis of flavour anomalies in section III by considering only the SM particle spectrum along with a massive leptophilic and family non-universal U (1) X gauge boson. We then move onto the discussions of the complete models in sections IV, V covering the details of flavour anomalies, dark matter and neutrino mass. In section VI, we discuss the possibility of explaining KOTO anomaly within our toy models. In section VII we discuss about different Higgs invisible and charged lepton flavour violating decays and also comment on other possible ways to probe our model at the LHC and finally summarise our findings in section VIII.

II. OUR FRAMEWORK
As mentioned before, our goal is to extend the SM by an Abelian U (1) X symmetry with a corresponding massive gauge boson X. We restrict our study to only low mass regime (GeV scale) of this additional gauge boson and allow only the leptons to couple to it. The charge assignments of the different SM particles under the different gauge groups are listed in Table. I and the NP interaction Lagrangian is given by where g X is the gauge coupling of the U (1) X group, i represents the lepton generation and n i are the charges of the lepton families under U (1) X which we want to constrain from anomaly cancellation requirements as well as flavour phenomenology. Here, in the above Lagrangian, L i is the left-handed lepton doublet while e R i is the right-handed singlet with same gauge charge n i . Here, while writing the above Lagrangian, we have assumed that the U (1) X charges for the right and left-handed leptons are same, leading to a vector type interaction. In eq. (3), B µν and X µν are the standard U (1) Y and U (1) X field stress tensors, respectively, and the factor represents the kinetic mixing between them. We assume that the leptophilic X mixes kinetically with the SM Z boson with a strength . This mixing will be helpful to get contributions in various low energy observables like R(K), R(K * ) through penguin diagrams with the lepton vertex dominated by the above interaction and the one-loop quark vertex modified by the mixing parameter . In muon or electron anomalous magnetic moments or in other lepton flavour violating (LFV) decays, at leading order, this mixing parameter does not have any specific role.
(1, 2, − 1 2 ) n3 eR (1, 1, −1) n1 µR (1, 1, −1) n2 τR (1, 1, −1) n3 As mentioned before, assigning charges to the SM fermions under a generic U (1) X symmetry leads to non-zero contributions to the one-loop triangle diagrams and makes the model anomalous. Therefore in order to realise a anomaly-free renormalisable model, one needs to put additional chiral fermions into the model which may also provide a natural candidate for DM. At the same time the additional chiral fermions required for anomaly cancelation could be made useful for neutrino mass generation as well. For similar construction of Abelian gauge extended models in the context of DM and neutrino mass generation, see [37][38][39][40][41][42][43][44][45] and references therein.
The equations that govern the anomaly cancellation requirements in our setup are given by : (E) [U(1) X ] : From the above set of conditions (A-E) one can infer that : • n 1 + n 2 + n 3 = 0 ensures anomaly cancellation of all the anomalies except eq. (7).
• In order to ensure eq. (7) is also zero, we can add N extra fermions with U (1) X charges (m 1 , m 2 , ... etc.) such The one way of cancelling the anomaly without adding more fermions is to consider equal and opposite charges for any two generations of leptons and let the charge of the third generation be zero. These are the symmetries like U (1) Le−Lµ , U (1) Lµ−Lτ which has been discussed earlier in the references [29,31,33]. However, if we want to consider non-zero charges for all the three lepton generations, then we need to have additional chiral fermions in our model for anomaly cancellation. So without choosing random charges and adding fermions in an ad-hoc manner, we can try to constrain the possible values of n 1 , n 2 and n 3 from the available low-energy data. Note that n 1 and n 2 will be sensitive to the observables like R(K ( * ) ) as well as electron and muon anomalous magnetic moments. There will not be any contributions to the lepton flavour violating decays and the rare decays like B s → µµ or B s → ee. Also, depending on the lepton in the final state, the b → c ν decays (with = e, µ, τ ) will be sensitive to the charges as mentioned above. However, due to the low mass of X, the new contributions in B → D ( * ) ν decays are much smaller as compared to the corresponding SM counterpart. Therefore, effectively we can get constraints on n 1 and n 2 using the available data on b → s decays (for = µ, e) and anomalous magnetic moments; however, due to unavailability of sufficient data, n 3 can not be constrained. We then look for possible solutions for the charges (n 1 , n 2 , n 3 ) such that n 1 + n 2 + n 3 = 0. Such a prescription also allows us to constrain the mass of X and the kinetic mixing parameter effectively. The detailed analysis is described in the next section.

III. ANALYSIS
In the following subsection, we will discuss different observables which will be useful to constrain various model parameters like U (1) X charges of leptons n i , new gauge coupling g X , new gauge boson mass M X , and the kinetic mixing parameter .
A. Exclusive b → s (with = e, µ) decays As mentioned earlier, the measured values of R(K ( * ) ) in the semi-leptonic B-meson decay reported by the experimental collaborations provide an indication of lepton flavour universality violation (LFUV). The measured value of R(K) by LHCb is given by [14] R(K) = 0.846 +0.060 +0.016 −0.054 −0.014 , for q 2 ∈ [1.1, 6] GeV 2 , where q 2 is the squared momentum of the leptons in the final state. This result has a deviation from the SM prediction by 2.5σ. Similar measurements are available for R(K * ) by LHCb and Belle collaborations. While the LHCb Collaboration has reported [46] R( Belle presented their first measurement [15] of R(K * ) in B 0 and B + decays in April 2019 which reports Although the measurements from Belle are compatible with the SM expectations, they have comparatively large uncertainties. Thus, considering the more precise results from LHCb, the anomaly in R(K * ) stands at ∼ 2.4σ. The effective Hamiltonian for the b → s transitions is given by [47]: where O i and O i 's are the dimension six effective operators and C i 's are the corresponding Wilson coefficients (WC). Although the semi-leptonic operators O ( ) 9 ∝ (sγ µ P L(R) b)(μγ µ µ) and O ( ) 10 ∝ (sγ µ P L(R) b)(μγ µ γ 5 µ) are relevant for the decay b → s + − . The analysis with the very recent data suggests that O 9 is the only one operator scenario that can simultaneously explain all the data in b → s decays [16,17]. However, there are a few two or three operator scenarios which can best explain the data at the moment, and that includes the combination O 9 and O 10 [16,17]. In our model, the leading contributions to the Wilson coefficients will come from the diagrams shown in Fig. 1, and we will have contributions only in C 9 due to the vectorial coupling of the X to the leptons. There will be contributions in both b → sµ + µ − and b → se + e − decays. As can be seen from eq. (3), due to the absence of axial-vector coupling of X to the leptons, we do not have contributions to C 10 . Therefore, at the leading order, the new Wilson coefficient (WC) is given by with Here, m t and M W are the top quark and W -boson masses, respectively. The sine of the Weinberg angle is defined as S W = sin θ W and C W = 1 − S 2 W . Also, n ≡ (n 1 , n 2 ) depending upon the lepton flavour it contributes to. Note that we are working in a model with the mass of X in the GeV or sub-GeV range, in particular, we are focusing in the region M X > 2m µ . On the other hand for B → K ( * ) decays, the allowed values of q 2 lie in the range 4m 2 < q 2 < (M B − M K ( * ) ) 2 . In such a situation, one cannot Taylor expand the X propagator in powers of q 2 /M 2 X . Therefore, the new WC, as shown in eq. (13) will have explicit q 2 dependence and in general, could be complex. Note that for the X-boson, we have introduced the Breit Wigner (BW) propagator. In this form of the propagator, we will get a finite analytic expression for the amplitude at the resonance region. This is because, around the mass of X, the zeroth-order propagator vanishes and the higher-order effects are leading, which is given by the imaginary part proportional to the X decay width. The imaginary part will receive contributions from every particle into which X can decay. In general, without a priori knowledge of all the decay channels of X, it is hard to predict its total decay  width. However, we have considered a leptophilic X, and its primary decay channels are the dilepton final states, like + − and νν with = e, µ. Hence, one needs to estimate the decay width Γ X ≈ Γ(X → ) + Γ(X → νν).
In this model, there are free parameters which need to be constrained using the existing data. In particular, the constraints from low energy experiments, like neutrino trident production (NTP) bound, rare kaon decay K + → ν µ µ + X(→ νν), BaBar 4µ channel search 2 etc. along with cosmological observations of Big Bang nucleosynthesis (BBN) are important. As can be seen from [49], the current data allow a gauge coupling g X ∼ 0.001 for M X ∼ 0.5 GeV and it could be > ∼ 0.002 for M X > ∼ 1.0 GeV. On the other hand, the kinetic mixing parameter is constrained from neutrino-electron scattering experiments like CHARM-II, GEMMA and TEXONO, for details see [50]. Mixing strength > ∼ 10 −3 is ruled out for gauge bosons of mass around the electroweak (EW) scale. For keV scale bosons, the bound is even tighter O(10 −6 ). LEP II has put a lower bound on the ratio of new gauge boson mass to the new gauge coupling to be M X /g X ≥ 7 TeV [51]. However, since we are interested in the low mass of the gauge boson, bounds from hadron colliders like ATLAS and CMS will not be very relevant. Similarly, LEP bound is also not applicable in such low mass regime 3 . With all these inputs, the X decay width as mentioned above, will be of order O(10 −9 -10 −7 ) GeV for g X ∈ (10 −4 , 10 −3 ), which is much smaller than M X . In the limit Γ X M X → 0 (narrow-width approximation (NWA)), the BW becomes a delta distribution: δ(q 2 − M 2 X ). LHCb has done a dedicated search for light hidden-sector bosons by measuring the branching fraction B(B 0 → K * 0 χ(µ + µ − )) = B(B 0 → K * 0 χ) × B(χ → µ + µ − ). Here, χ is the light boson in the hidden sector similar to X in our case. Depending on the lifetime τ (χ), LHCb has put bounds on the above mentioned branching fraction for a given mass range of χ [53]. One can refer to Fig.7 of the supplemental material of reference [53] in which the ratio B(B 0 → K * 0 χ(µ + µ − ))/B(B 0 → K * 0 µ + µ − ) has been plotted as a function of m(χ) (with 214 ≤ m(χ) ≤ 4350 MeV) for different values of τ (χ) including τ (χ) = 0. Here, the branching fraction B(B 0 → K * 0 µ + µ − ) is defined for 1.1 < m 2 (µ + µ − ) < 6.0 GeV 2 . Note that if we choose the lifetime τ (χ) = 1000 ps, which corresponds to a very small decay width of X, the ratio as mentioned above could be of order one. However, the bounds on the same ratio will be ≤ O(10 −2 ) (at 95% Confidence Level (CL)) for τ (χ) = 10 ps, which is even the case in the limit τ (χ) → 0. The width Γ X ≈ 10 −9 GeV corresponds to a lifetime ≈ 10 −4 ps which is close to zero.
In our model, we have estimated B(B 0 → K * 0 X(µ + µ − ))/B(B 0 → K * 0 µ + µ − ) within the accessible ranges of M X . The normalisation B(B 0 → K * 0 µ + µ − ) has been measured by LHCb for 1.1 < q 2 < 6.0 GeV 2 [54], which is given by (1.6 ± 0.3) × 10 −7 . The dependences of this ratio on different model parameters like , the charge n 2 and the coupling g X are shown in Fig. 2. A close inspection of Figs. 2b and 2a suggests that if we choose < ∼ 10 −4 , for values of n 2 as large as 5, the constraints from LHCb will be satisfied within the accessible ranges of M X . In the rest of our analysis, we will consider ≈ 10 −4 . Note that even for a relatively large gauge coupling (∼ 10 −3 ) we can still be able to satisfy the upper bound provided by LHCb. However, in such cases, the allowed values of n 2 will be very much restricted (as shown in Fig. 2c). (1) and (2), respectively.

B. Anomalous Magnetic Moments
Another important observable which could be useful to put tight constraints on the model parameters is the anomalous magnetic moments of muon or electron. As one can see from eq. (3), since we do not have lepton-flavour violating couplings of X, the gauge boson mediated diagram will not contribute to decays like τ → µγ, µ → eγ etc.
The effective vertex of photon with any charged particle is given by: The factor g µ ≡ 2(F 1 (0) + F 2 (0)), and the anomalous magnetic moment is given as a µ ≡ F 2 (0) = 0 (since F 1 (0) = 1 at all order). In our model, the diagram that will contribute to muon and electron anomalous magnetic moments is given in Fig. 3. In our model, the contribution to ∆a is given by where , ≡ e, µ and n (= n i ) denotes the U (1) X charge of the lepton. Our analytical expression can be compared with the one obtained in [55]. Note that the contributions in ∆a for both = µ and e are positive; however, in the case of electron magnetic moment, the expectation is negative. Also, as compared to the requirement, the contribution in electron anomalous magnetic moment is negligibly small. The dependences of ∆a µ on various model parameters are shown in Fig. 4. We can easily explain the excess in ∆a µ for values of g X of order O(10 −3 ), and the data prefers a value of M X < ∼ 1 GeV. In such situation, the value of n 2 need not be >> 1. However, if we choose g X ≈ 10 −4 then in order to explain the excess in muon (g − 2), we need relatively larger values of n 2 (>> 1).
We have already pointed out in the introduction that the current measurement of the fine structure constant poses a negative ∼ 2.4σ deviation in anomalous magnetic moment of electron from its theory prediction [25]: In electron anomalous magnetic moment, the contribution from the diagram in Fig. 3 will be positive and is given by The allowed parameter spaces for n 1 , n 2 , g X and M X which passes the constraints from R(K ( * ) ), the anomalous magnetic moment of the muon and B(B → K ( * ) µµ). Note that the observables in b → s decays are defined in high ( 1.1 < q 2 < 6.0 (GeV 2 )) as well as in the low-q 2 (0.045 < q 2 < 1.1 (GeV 2 )) regions.
for M X = 1 GeV, g X = 0.001 and n 1 = 1. Hence, we can not explain the current trend of data in ∆a e with only an additional U (1) X gauge boson.

C. Combined parameter spaces
In this subsection, we discuss the constraints obtained on the model parameters from a simultaneous analysis of the observables in b → s decays and ∆a µ . As we can see from eqs. (10) and (11), data are available in two different q 2 regions, one for q 2 ∈ [0.045, 1.1] GeV 2 (low-q 2 ) and the other for q 2 ∈ [1.1, 6] GeV 2 (high-q 2 ). In our analysis, we have considered the inputs from R(K ( * ) ), B(B → K ( * ) µ + µ − ) (in both the q 2 regions) and ∆a µ . For R(K ( * ) ), we have not considered the Belle data since it has significant errors, and we have considered the LHCb data on it at their 2σ CL. Similarly, ∆a µ has been considered in its 3σ CL interval.
As mentioned earlier, the mixing parameter plays a crucial role in constraining the other relevant new parameters. We have noted that when ≈ 10 −4 the allowed regions of the other parameters are more relaxed than the one obtained for ≈ 5 × 10 −4 . We scan the parameters over the following intervals: 0.5 ≤ M X ≤ 1.5 (in GeV), −5 ≤ n 1 ≤ 5, −5 ≤ n 2 ≤ 5, 0.1 ≤ g X (×10 3 ) ≤ 3. Here, we would like to mention that the low mass regions 0.22 < M X < 0.5 (GeV) are also allowed by the data as discussed above, in the next section, we will show it in a specific scenario. The allowed parameter spaces for = 1 × 10 −4 are shown in Fig. 5. Note that n 2 and n 1 have a nice correlation, higher positive values of n 2 prefers higher negative values of n 1 . Also, within our chosen parameter values, only negative values of n 1 are allowed. For a fixed value of n 2 , a wide range of values of n 1 is allowed. However, as expected, the scenario n 1 = n 2 is excluded. Here, we have shown only the positive values of n 2 , which are allowed by the data. The allowed values of n 2 are symmetrically distributed about the origin along the n 2 -axis. In addition, we see that for ≈ 1 × 10 −4 , within the given range of M X , the allowed values of the coupling g X lies in between 0.5 × 10 −3 and 2 × 10 −3 . To be conservative, we have not considered values of g X larger than 2 × 10 −3 since other low energy observables constrain higher values, for details see [49].

IV. THE EXTENSION OF U (1)X WITH ADDITIONAL DEGREES OF FREEDOM
A certain combination of n 1 , n 2 would lead to a particular extension of the SM with chiral fermions. While such extension is not unique, we stick to minimal possible extensions in order to address the problems discussed earlier. Therefore we can now proceed towards making a specific choice for these parameters in order to complete our model in a way that the extension is minimal. We have already been able to constrain n 1 and n 2 from low energy data while n 3 remains unconstrained. One can easily see that the minimal way to cancel the anomalies would be to add three chiral singlet fermions with U (1) X charges n 1 = −n 1 , n 2 = −n 2 and n 3 = −n 3 . This will make the sum of the charges as well as sum of the cubes of the charges equal to zero. The three fermions can be considered to be 3 right-handed neutrinos (RHNs). As a benchmark scenario we choose (n 1 , n 2 , n 3 ) = (−1, 2, −1) which is in good agreement with the flavour data as shown in Fig. 5. For n 2 = 2,n 1 = −1, the correlation between M X and g X is shown in Fig. 6, here we have shown the region 0.25 < ∼ M X < ∼ 1.0 (GeV). For these values of [n 2 , n 1 ] and for g X = 0.001, within the allowed ranges of q 2 and M X the numerical values of the WCs ∆C µ 9 and ∆C e 9 will lie in between [−0.827, −1.83] and [0.413, 0.91], respectively. These values of the WCs are consistent with the result (within 2σ CI) obtained from a global fit to all the available data in b → s decays considering NP effects in both the muon and electron final states [18]. One can make the fermion content richer by adding more chiral fermions with appropriate charges that satisfy the anomaly cancellation requirements. However, we would like to have a plausible explanation for the neutrino masses, and at the same time, we want to keep our model minimal. Therefore, we extend our model with only three RHNs. All these fermions couple directly to SM leptons via SM Higgs (due to equal and opposite U (1) X charges of right and left handed leptons while SM Higgs remains chargeless under it), and therefore we cannot consider one of them to be our DM candidate. One can, of course, add a Dirac fermion on top of this which will not contribute to any anomaly and assign this to be the DM. But such a scenario will be ad-hoc and less motivating since the DM does not arise naturally from the anomaly cancellation requirements. Also, its mass remains a free parameter without being connected to the scale of U (1) X symmetry breaking. Thus we need to look beyond this minimal solution by extending the particle content further 4 . On the other hand, imposing a discrete Z 2 symmetry on the new chiral fermions (at least in one of them) will help us in forbidding their direct coupling to SM fermions and SM Higgs. In non-minimal or UV complete version of such minimal scenarios, it is possible to realise such Z 2 symmetry as a remnant after spontaneous symmetry breaking of U (1) X [39][40][41][43][44][45]. Our minimal setup here will enable us to have a DM candidate without adding new fermions apart from the RHNs. Under such a scenario, there are two possibilities with the different origin of light neutrino masses but with almost the same DM phenomenology, which we discuss in the following section.

V. TOY MODELS
In the following subsections, we discuss the toy models which have been built considering the U (1) X charge assignments of the SM leptons and new chiral fermions as described in the previous section i.e. (n 1 , n 2 , n 3 ) = (−1, 2, −1) and (n 1 , n 2 , n 3 ) = (1, −2, 1). We consider the additional fermions (namely, N 1 , N 2 and N 3 ) to be right-handed, hence, their U (1) X charges will be the sign-flipped version of (n 1 , n 2 , n 3 ) i.e. (−1, 2, −1). Based on how we are imposing the Z 2 symmetry on the new chiral fermions, one can come up with different models, and here we will discuss two such toy models.
A. Toy Model I

Particle Content
In this scenario, we consider that all generations of the RHNs, N i (i = 1, 2, 3), to be odd under a discrete Z 2 symmetry while all the SM particles are even. Thus to write a Yukawa term for the RHNs with the SM leptons, we would require an additional Higgs doublet (H 2 ) which is also odd under this discrete symmetry. The unbroken Z 2 (1, 1, 0) 4 + symmetry prevents H 2 from acquiring a non-zero vacuum expectation value (vev), and it remains inert. However, it plays a crucial role in neutrino mass generation by the radiative seesaw mechanism [57], which has been described later. To give mass to the chiral fermions, we require at least two singlet neutral scalars with non-zero U (1) X charges which we choose to be 2 and 4, respectively. The lightest singlet neutral fermion can be a suitable DM candidate since the Z 2 symmetry protects its decay into other lighter particles. However, since N 2 is unique from the other singlet fermions in terms of its U (1) X charge, so we consider this to be our DM candidate and ensure that it is the lightest among all the Z 2 odd fermions. In Table II, we have shown the entire particle content alongside with their respective charges with respect to different symmetries of the model.

Lagrangian and Scalar Mass Spectrum
In the set up given above, the total Lagrangian can be written as where n i , n i are the U (1) X charges of the SM lepton generations and RHN generations, respectively. The relevant Yukawa interactions of the RHNs is given by : The scalar Lagrangian L S can be written as: where the covariant derivative is given by Here, G X = (g X X + g Y ) and (Y, X) are the hypercharges related to U (1) Y and U (1) X gauge groups respectively. In the above mentioned scenario, X = 4 and 2 for ϕ 2 and ϕ 1 , respectively, while X = 0 for H 1 and H 2 . As defined earlier, is the kinetic mixing parameter. The scalar potential V (H 1 , H 2 , ϕ 1 , ϕ 2 ) is defined as with the doublet and singlet scalars after the electroweak symmetry breaking (EWSB) are defined as After spontaneous symmetry breaking, all the scalars apart from H 2 acquires a vev and is responsible for giving mass to other particles. In order to spontaneously break the electroweak symmetry as well as U (1) X , we must have µ 2 1 < 0, µ 2 3 < 0 and µ 2 4 < 0. Also since the inert doublet does not acquire a vev, µ 2 2 > 0. Here, the term proportional to δ in the scalar potential (23) will play an important role in determining the mass of pseudo-scalars like A 2 . The potential minimization conditions are given by The gauge boson mass term can be obtained from the kinetic terms in eq. (21) which is given by with Note that H 2 does not acquire a vev; hence, it does not play any role in the mass generation of the gauge bosons or fermions. We obtain the masses of W -boson, Z-boson and photon as in case of SM. The X boson mass has been obtained as a combination of the vevs of the singlet scalars and the vev of H 1 ; the contributions from H 1 is suppressed by the factor 2 .
In eq. (26), to obtain the masses of the neutral gauge bosons, we need to carry out the standard electroweak rotation as given below where A µ is the photon field. Note that after the symmetry breaking there will be a remaining mixing between Z 0 µ and X µ , which can be written as: Here, we have neglected the mixing between the photon and the new gauge boson. The masses of the physical heavy gauge bosons (Z, Z ) can be obtained after diagonalising the above matrix by a rotation, and the masses are given by In the limit that M X << M Z 0 and mixing parameter << 1, we obtain the masses as and the mixing angle is given by : On the other hand, the mixing mass matrix for the CP even and Z 2 even neutral scalars (h , s 1 , s 2 ) is given by : The physical scalars (h, s 1 , s 2 ) are obtained after diagonalising the above mass mixing matrix and they are related to the unphysical ones by an orthogonal transformation. We consider a general real orthogonal 3 × 3 rotation matrix O with three mixing angles (no phase) for diagonalising the above mentioned mass mixing matrix as where c αij ≡ cos α ij and s αij ≡ sin α ij . In order to make the notation simpler, we redefine the angles as α 12 ≡ α 1 , α 13 ≡ α 2 and α 23 ≡ α 3 . Another important variable is the ratio between the vevs v and v 1 which we have defined as tan β = v1 v . In general, to keep the analysis simple, we can assume that the mixing of s 2 with s 1 and h are negligibly small. In such situation, we need to focus only on the mixing between s 1 and h, i.e s α1 or c α1 . There are studies on the singlet scalar extension of the SM, and bounds are available on the respective model parameters like tan β and s α1 ; for example, see [58][59][60]. These studies took into account the bounds from various experimental measurements like the precision observables S, T and U parameters, W -mass, LEP and LHC bounds. Alongside, they have considered various theory inputs, like perturbative unitary constraints on scalar self-interactions, vacuum stability, etc. All these studies suggest that for 300 ≤ M s1 ≤ 800 (in GeV), one can safely assume | sin α 1 | ≤ 0.3 and tan β > 1. Note that in our model, we have two singlet scalars, and as discussed above we have more free parameters. In general, we can expect that the bounds as mentioned above will be little more relaxed in case of our model parameters. However, to be on the safe side, we have used these bounds in our analysis. This will help us to constrain a few of the other model parameters. In our analysis, we have considered |sin α 1 | < ∼ 0.2 and tan β = 2.0 and M s1 = 500 GeV [59,60] which is even more conservative. The corresponding values of M s2 can be obtained from eq. (27) after the evaluation of v 2 .
The mass mixing matrix for the CP odd and Z 2 even neutral scalars (A 1 , A 2 ) is given by After diagonalizing this matrix with an orthogonal transformation we will obtain one massless goldstone (A 1 ) corresponding to the gauge boson of U (1) X and another massive physical CP odd scalar (A 2 ) of mass − 2 expression that the dimensionful coupling δ has to be negative. Also here, for simplicity, we can limit our discussion to the value s γ << 1. The masses of the neutral and charged inert scalars are given by : where λ L = (λ 3 +λ 4 +λ 5 ) and λ A = (λ 3 +λ 4 −λ 5 ). To summarise, in Appendix B, we have presented various couplings in terms of the relevant physical masses, vevs and the mixing angles. These are the most general relations from which one can obtain the approximate relations for small mixing angle. The coupling strength of the interaction between H 1 and H 2 is defined by λ L = λ 3 + λ 4 + λ 5 . In an inert two Higgs doublet model (2HDM) where H 0 is considered as a suitable DM candidate, the bound on this type of coupling is given by λ L < ∼ 6 × 10 −3 [61]. We have not explored this possibility. In our study, the doublet H 2 is relevant for the neutrino mass generation and the required coupling is λ 5 which we have treated as free parameter.

DM Phenomenology
We adopt the thermal DM paradigm where DM gets produced in the early Universe thermally followed by its freeze-out from the thermal bath which decides its present day abundance. The relic abundance of DM can be computed by solving the appropriate Boltzmann equation and the model parameters can be constrained by comparing the calculated relic with observed abundance which, in terms of density parameter Ω and h = Hubble Parameter/(100 km s −1 Mpc −1 ), is conventionally reported as [62]: Ωh 2 = 0.120 ± 0.001 at 68% C.L. We solve the Boltzmann equation numerically using micrOMEGAs [63] where the model information has been supplied to micrOMEGAs using FeynRules [64].
In our model, we choose N 2 as the DM candidate which is supposed to be the lightest RHN. Note that its U (1) X charge is different from the other two RHNs. The dominant contributions to the relic abundance of DM will come from the annihilation diagrams shown in Fig. 7. There are a few other diagrams which are shown in Fig. 33 in Appendix A whose contributions in the DM relic abundance will be sub-leading 5  In this model, apart from M X (≡ M Z ) and g X (≡ g Z ) the other parameters that are relevant for DM phenomenology , respectively. The co-annihilation diagrams in Fig. 34 are sensitive to Y ij (with i, j = 1 or 3). Therefore, the relic density is almost insensitive to these parameters since the contributions from these diagrams are suppressed. Considering the bounds from the low energy data in the rest of our analysis, we have fixed the mass M Z at 1 GeV; also, we have set g Z ≈ 10 −3 . As mentioned earlier, with a particular choice of tan β one can fix the value of M s1 since M h is known. Once this is done the allowed values M s2 can be fixed from eq. (27). In this regard, the perturbativity of the scalar couplings will also play an important role. Since we have chosen tan β ≈ 2 and M s1 ≈ 500 GeV, the corresponding values of M s2 and M A2 will be limited to < ∼ 200 GeV. Accordingly, the mass of DM will be restricted because s-channel annihilations are the dominant annihilation process for the DM.
In Note that for M s2 < M h the bound on relic density is satisfied for DM masses close to, but less than M s2 /2. In such scenario, the relic density is under-abundant at or near M DM ∼ M h /2. On the contrary, when M s2 > M h the relic density is satisfied at a DM mass close to M h /2, provided we assume that there is a mixing between h and s 2 , i.e s α2 = 0. As discussed earlier, we have restricted our analysis to the small values of s α2 (≈ 0.01). Later we will see that such a restriction will be useful to evade stringent bounds on the direct detection of Higgs portal DM from XENON-1T experiment [66]. However, we have shown the dependences of the relic density on s α2 in Fig. 8c. As expected, in the no-mixing scenario, the resonance peak due to the Higgs mass vanishes. In such situation, when M s2 > M h , the relic density will be satisfied for values of DM mass close to M s2 /2 instead at M h /2, which is the case for non-zero s α2 . Note that there will be another resonance peak at M DM ∼ M A2 /2, at or around which the relic density will be much lower than the existing bound. One can also see from the Figs. 8a and 8d that the co-annihilations with inert scalars do not play much role since the relic abundance do not change with change in the value of coupling Y 22 and the masses M H0 or M H ± . This is precisely due to the fact that the scalar masses are much larger than DM mass making the coannihilation processes inefficient.
In Fig. 9a, we have shown the allowed ranges of the DM mass obtained from a scan with the constraints from relic density projected against the upper limit on direct detection cross-section σ SI (spin-independent) of DM from XENON 1T experiment [66]. To generate this plot we consider tan β = 2, and the values of the other relevant parameters are the following: s α1 ∼ 0.01, 0 < s α3 < 0.01, 0 < s α2 < 0.01, 20 ≤ M s2 ≤ 200 GeV and M s1 ∼ 500 GeV. In this model, the spin-dependent direct detection cross-section is highly suppressed; hence we have not considered it for a numerical study. Note that the allowed values of the DM mass lies in between 5 and 90 GeV. It is evident that there will be an allowed region near M Z /2 which is not shown in this plot. The dependencies of σ SI on s α2 and s α3 are shown in Fig. 9b. As expected, the large values of s α3 allows the relatively larger values of s α2 , however, we will stick to the low values of both like s α2 ≈ s α3 = 0.01.
FIG. 10: Diagrams contributing to flavour changing charge current process b → c ν . predictions and the relevant measurements the reader can look at [7,[19][20][21]. A certain degree of discrepancy has been found between the predictions and their respective measurements. Also, in both R(D) and R(D * ), the measured values are higher than the respective predictions. In a model-independent effective theory approach, the Hamiltonian describing the b → c ν transitions with all possible four-fermion operators in the lowest dimension is given by where the operator bases are defined as and the corresponding Wilson coefficients (WC) are given by C W ( W = V 1 , V 2 , S 1 , S 2 , T ). In the above-mentioned basis neutrinos are assumed to be left handed. The other theory details and the results of the model-independent new physics analysis on these modes can be seen in [21,[67][68][69] and the references therein.
We have noticed that due to the coupling of the Z to the lepton families, there would be a vertex correction diagram contribution to the channel b → cτν τ as shown in Fig. 10a. However, that contribution is not sufficient in order to explain the anomalies in R(D ( * ) ). With the addition of the inert scalar doublet and RHNs, we will have additional diagrams as shown in Fig. 10b. The lepton vertex is modified due to the loop corrections coming from the scalars H ± , H 0 and N i , (i = 1, 2, 3). Hence one can obtain a bound on the Yukawa couplings Y ij of eq. (19) and masses of H ± , H 0 and RHN from the semi-leptonic b → c decays.
The diagram given in Fig. 10b will contribute to C V1 of eq. (40) which is the WC of the four-fermion operator O V1 in eq. (41). The following is the corresponding mathematical expression : and Here (i, j) denote the generation index of the lepton and RHN respectively and M H denotes the mass of H 0 or A 0 running in the loop. Hence depending on the generation of RHN running in the loop, we would have different contributions to C V1 corresponding to each lepton flavour. As shown earlier, N 2 is our DM candidate, hence, the corresponding diagram will contribute only toB → D ( * ) µν µ decays and will be proportional to Y 2 22 . There will be other diagrams with N 1 and N 3 which can contribute simultaneously to bothB → D ( * ) τν τ andB → D ( * ) eν e decays. Since the loop factor as mentioned above is sensitive to M Ni , if we consider the masses of N 1 and N 3 equal, for simplicity, then the total NP WC for the tau mode would be proportional to R 2 3 ≡ (Y 2 31 + Y 2 33 ) while that for the electron mode would be proportional to R 2 1 ≡ (Y 2 11 + Y 2 13 ). If we assume that Y 13 << Y 11 or Y 31 << Y 33 then we can write down the following approximate relations: R 2 3 ≈ Y 2 33 and R 2 1 ≈ Y 2 11 . Since the semi-leptonic branching fractions of the B meson into the light lepton channels are precisely measured [70,71] and SM consistent, global fits of the NP Wilson coefficients to R(D), R(D * ) data [21] is done by considering NP only in the τ decay modes. However, since we have contributions to all semi-leptonic decay channels, we consider NP in both the numerator and denominator of R(D), R(D * ) while also ensuring that the contribution to the light lepton modes does not overshoot the experimental limits on their branching fractions.
We perform a parameter space scan of R 1 , R 3 and M N (≡ M N1 = M N3 ) by fixing Y 22 = 0.1 and the masses of the inert scalars as shown in Fig. 11. The blue and the orange points are the allowed regions for M H ± = 500 GeV and 1000 GeV respectively, when the RHN masses are varied between 100 − 500 GeV. All these allowed points satisfy the experimental constraints on R(D), R(D * ) and the branching fraction ofB → D ( * ) ν at their respective 2σ confidence interval (CI). These parameter spaces also satisfy the bound B(B c → τ ν) < 30% and the corresponding expression in terms of C V1 can be seen in [72]. From Fig. 11a it is evident that the data prefers R 3 > R 1 . Also, as expected, athough it is allowed, we don't necessarily need large values for R 1 and the data allows a solution like R 1 ≈ 0 while R 3 > ∼ 1. In our model, C V1 is positive; hence, the new contribution will interfere constructively with the SM and increase the relevant branching fractions from their SM predictions. Notice that there are minimal dependencies of R D ( * ) on M H ± and M H 0 . However, as can be seen from Fig. 11b, it is almost independent of M N .
The similar type of diagrams as given in Figure 23 with the replacement c → u will contribute to b → uτν τ processes like B ± → τ ± ν τ ,B 0 → π + τ −ν τ decays. We have checked that the required values of the parameters for an explanation of the data in b → cτν τ decays can accommodate the current observation B(B ± → τ ± ν τ ) = (1.09 ± 0.24) × 10 −4 [1]. Similarly, our model will contribute to s → uτν τ , c → sτν τ decays which will lead to semileptonic and purely leptonic decays of K and D/D s -mesons, respectively.

Anomalous Magnetic Moment and LFV
a. Magnetic moments :-In toy model-I, the additionl contributions to anomalous magnetic moments will come from the type of diagram given in Fig. 12 (for i = j) with charged Higgs and RHNs in the loop. The contribution is given by where (i, j) denotes the generation of the lepton and RHN respectively, λ j = which are suppressed compared to the gauge boson mediated diagram for it, which is shown in Fig. 3 (with X in the loop), by two or three orders in magnitude. On the contrary, as shown in eq. (18), the contribution to electron anomalous moment from the same diagram is negligibly small. Therefore, by extending the symmetry of the SM by an abelian U (1) X gauge group without additional degrees of freedom, we cannot explain the observed discrepancy in the electron magnetic moment. Note that ∆a e (eq. (17)) has a significant error, and at 3-σ CI it is consistent with zero. Therefore, it would be too early to prejudge the potential impact of new physics on this observable. Here, we will show that our model has the potential to predict negative values ∆a e although it is difficult to explain the data in one or two-σ CI. From eq. (44) it is clear that ∆a e is sensitive to our predefined variable R 1 ≡ Y 2 11 + Y 2 13 . In Fig. 13, we show the variation of ∆a e with M N for two different values of R 1 and three values of M H ± . The chosen values of R 1 , M N and M H ± explain the data on R(D), R(D * ) as previously shown in Fig. 11. From the plot we see that the contribution is negative and small and we can not reach the present experimental limit within its 1 or 2-σ CI. To get a large negative contribution we need small values of M H ± and M N , and a large value of R 1 (R 1 > 1). However, R 1 >> 1 is not allowed by B(B → D ( * ) ν) data. Therefore, within the allowed parameter space of R 1 and M N , we have a contribution to the electron anomalous magnetic moment which is negative and of order O(10 −14 ) and within the 3σ range of the experimental data.
b. Lepton flavour violation:-The same one loop diagram given in Fig. 12 will also contribute to the LFV process τ → eγ. Therefore, one must ensure that the contribution is within the current experimental limit B(τ → eγ) < 3.3 × 10 −8 [1]. However, in our model there will not be any contribution to τ → µγ or µ → eγ. The expression for the partial decay width Γ( i → j γ) for the diagram in Fig. 12 is given by [73] : where, and r ≡ In the above expression for the decay width we have a combination |Y 11 Y 31 | or |Y 33 Y 13 | depending on whether N 1 or N 3 runs in the loop. So we can constrain the allowed values of these product couplings from the experimental upper limit on the branching fraction of τ → eγ. As we have seen earlier, if we assume that the off-diagonal Yukawas are much smaller in value than the diagonal ones, the data on B → D ( * ) ν allow R 3 ≡ Y 33 ∼ 1 and R 1 ≡ Y 11 ∼ 1 for M N in the range (100 − 500) GeV or more. Therefore, in general, the magnitude of the product couplings as mentioned above could be small even if we assume Y 33 ∼ Y 11 ≈ 1.
In Figs

Neutrino Mass Generation
The neutrino mass will be generated by radiative scotogenic mechanism in a way similar to the original proposal of [57] as depicted in Fig. 15 and will be mainly moderated by the mass splitting between H 0 and A 0 . The one-loop contribution is given by : where, M k is the mass of the RHN N k running in the loop. The Majorana mass matrix, however, has the following texture :  [62], as well as oscillation data on the neutrino mass squared differences and mixing angles [74,75]. Hence it is convenient to rewrite the Yukawa couplings in terms of the light neutrino parameters in order to automatically incorporate the above constraints on the couplings. One useful way of achieving this is through the Casas-Ibarra (CI) parametrisation [76] extended to the radiative seesaw model [77] which enables us to express the Yukawa coupling matrix as where, U is the usual Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix, D ν is the diagonal light neutrino mass matrix, R is an arbitrary complex orthogonal matrix satisfying RR T = 1 and Λ is a diagonal matrix with elements and where L i (m 2 ) is the mass function defined as Note that we are working in a basis where the charged lepton mass matrix is not diagonal. The PMNS mixing matrix can be parametrised as where c ij = cos θ ij , s ij = sin θ ij and δ is the leptonic Dirac CP phase. The diagonal matrix P = diag(1, e iα , e i(β+δ) ) contains the Majorana CP phases α, β that appears when ν is Majorana and are not constrained by neutrino oscillation data but has to be probed by alternative experiments. This leptonic mixing matrix is related to the diagonalising matrices of charged lepton and neutrino mass matrices as U = V † L U ν and as mentioned above, V L is not a unit matrix in our model. It consists of a rotation in (1 − 3) plane which can be parametrised as where c l 13 = cos θ l 13 , s l 13 = sin θ l 13 and δ l is an arbitrary phase which we assume to be zero for simplicity. Using this and the above parametric form of PMNS mixing matrix U , one can parametrise U ν which can then be used to parametrise the light neutrino mass matrix as In the above expression for M ν , the diagonal light neutrino mass matrix is denoted by M ν (diag) = diag(m 1 , m 2 , m 3 ) where the light neutrino masses can follow either normal ordering (NO) or inverted ordering (IO). For NO, the three neutrino mass eigenvalues can be written as M diag ν = diag(m 1 , m 2 1 + ∆m 2 21 , m 2 1 + ∆m 2 31 ) while for IO, they can be written as 2 23 , m 3 ) Structure of this parametric form of light neutrino mass matrix can now be compared with the structure of light neutrino mass matrix predicted by the model. Note that the model not only predicts a specific structure of right handed neutrino mass matrix given by eq. (49), but also predicts the Dirac Yukawa coupling matrix to have a similar structure Using the formula for light neutrino masses given in eq. (48), it can be shown that the above mentioned textures of Dirac Yukawa coupling matrix Y and right handed neutrino mass matrix M R lead to a very specific structure of light neutrino mass matrix with two independent zeros namely, (M ν ) eµ = (M ν ) µe = 0, (M ν ) µτ = (M ν ) τ µ = 0 where the equality (M ν ) αβ = (M ν ) βα results due to Majorana nature of light neutrinos giving rise to a complex symmetric structure of mass matrix. We numerically solve these two texture zero complex equations in order to evaluate the unknowns namely, the lightest neutrino mass m 1 (NO), m 3 (IO), leptonic Dirac CP phase δ as well as two Majorana CP phases α, β. The additional rotation angle in charged lepton sector θ l 13 is considered as a free parameter which can lie anywhere in (0, π/2). The other known parameters namely, three mixing angles, two mass squared differences are varied in 3σ range [75]. We find that these textures in light neutrino mass matrix predict a large value of the lightest neutrino mass, which is in tension with Planck 2018 bound on sum of absolute neutrino masses i m i ≤ 0.12 eV [62] as well as bounds on absolute neutrino mass scale from laboratory based experiments like KATRIN [78]. Even if we consider a non-zero CP phase in charged lepton correction matrix V L , this conclusion does not change. This is not surprising, given the fact that almost all possible two-zero textures in diagonal charged lepton basis are ruled out by latest experimental data [79].
One possible way to make it consistent with neutrino data without changing the model significantly is to change the U (1) X charge of the singlet scalar φ 2 from 4 to 1. This results in a right handed neutrino mass matrix having only one zero at (22) entry. While the lightest eigenstate of singlet fermion mass matrix can still be a DM candidate, no zeros appear in the light neutrino mass matrix even with the same Dirac Yukawa (57). Such a general structure of light neutrino mass matrix can be fitted with light neutrino data as there are sufficient free parameters, unlike in the previous case with two texture zeros. It is very unlikely that such a setup will change our DM and flavour physics results significantly. In the following subsection, we have added a discussion on this modified scenario.

Modified Setup for Toy Model I
As mentioned in the previous section, the light neutrino mass matrix that we obtain in this scenario violates the Planck 2018 bound on the sum of absolute neutrino masses. We also identified that a possible way out of this issue is by choosing the U (1) X charge of the singlet scalar φ 2 to be 1 instead of 4. In this subsection, we will briefly point out the changes that will occur in our theoretical setup and how it might affect the other observables. First of all, the Yukawa interactions given in eq. (20) will be modified as given below in eq. (58). Note that the first three terms of the Yukawa lagrangian remain unchanged, however, the interaction term involving N 2 and φ 2 has changed. Also, there will be a little change in the scalar potential, the trilinear term in eq. (23) now becomes δ ϕ 2 ϕ 2 ϕ † 1 + h.c. . Hence, the pseudoscalar mass, which primarily depended on this trilinear term, modifies to M A2 = − v 2 2 δ √ 2v1s 2 γ 1/2 . Recall that the gauge boson mass M X mass and gauge coupling g X are related to the singlet vevs (eq. (27)). With the change in the U (1) X charge of ϕ 2 , the aboove relation changes to and for << 1, we obtain M Z M X = g X (4v 2 1 + v 2 2 ). For simplicity, if we consider v 1 = v 2 , then from eq. (59), v 1 ≈ 450 GeV for M X = 1 GeV and g X 0.001. Therefore the masses of s 1 and s 2 will be restricted to be < ∼ 450 GeV for the Yukawa and quartic couplings to remain perturbative.
The analysis of relic abundance and the direct detection cross section will be in a similar line as discussed in subsection V A 3. The annihilation via Z remains the same. However, we can not consider a pure N 2 state as our DM candidaye, since the the Yukawa lagrangian does not have a majorana mass term for N 2 . Therefore, in principle, the lightest particle of N 1 and N 3 can be our DM candiate, and the dominating contributions will come from the annihilation diagrams shown in Fig. 16. In such situation, as before, depending on the mass of s 1 , the relic abundance will once again be satisfied near the resonances ie.e near M DM ∼ M s1 /2. In the presence of s-channel annihilation, the role of co-annihilations are expected to be sub-dominant as in the previous setup.
The possibility of mixing of the pure states N 1 , N 2 and N 3 can be considered by rotating the interaction basis N i to a new basis N i by using a general unitary transformation as which will result in a mass matrix of the form In the rotated basis, the lowest mass eigenstate can be considered as the DM candidate which will contribute via the annihilation diagram as given in Fig. 16. The Yukawa lagrangian responsible for the RHN masses also gets modified such that the Majorana mass mixing matrix now becomes This is in contrast to the mass matrix we obtained before in eq. (49). Since the mixing angles (s ν ij ) of O ν R are completely arbitrary, we have full freedom of choosing them in a way such that M 2 < M 1 , M 3 and the Yukawa couplings are also perturbative.
It is important to note that the contributions to the other observables like anomalous magnetic moments, LFV decays and R(D ( * ) ) remain unaltered. We have already seen that a charged Higgs and RHN mediated diagram contributes to the magnetic moments of the leptons (cf. Fig. 12). In the modified set-up, the changes occur in the Majorana-Yukawa interactions, which involves the coupling Y φ ij , and they do not contribute to ∆a µ,e , LFV decays or R(D ( * ) ).

Particle Content
In this toy model, we have the same particle content as in the previous case, with the only difference being that all particles except N 2 are even under the discrete Z 2 symmetry. This will once again prevent it from interacting directly with SM leptons. However, in this scenario, the neutrino mass generation mechanism will be different from the previous one. The particle content, along with their respective gauge quantum number and charges, has been described in Table III.

Lagrangian and Scalar Mass Spectrum
In this scenario, the successful generation of charged lepton and light neutrino masses require H 2 to be charged under U (1) X . The relevant Yukawa interactions are given by: where both i and j can take values (1, 3) and k = 2. Thus only the second generation of lepton doublet couples to N 1,3 via the second Higgs doublet H 2 . The scalar lagrangian will be similar to the one defined in eq. (21) with the scalar potential as given below: In this case, all the scalars acquire a vev and are given by : Under such a scenario, electroweak symmetry breaking of the scalars require µ 2 i < 0 (i = 1, 2, 3, 4) and the minimization conditions are given by: where, λ = c Λ 2 , the usefulness of this term will be discussed later in this subsection. The covariant derivative can be defined in the same way as in the previous case eq. (22). From the kinetic part of the scalar Lagrangian, we obtain the mass of the W-boson as : One can rewrite the mass of W as M 2 We also express the ratio of the two vevs as v u = tan β. The neutral gauge bosons (W 3 µ , B µ , X µ ) on the other hand mix and the mixing matrix is given by : After the usual Weinberg rotation as given in eq. (29), we obtain the masses of the physical neutral gauge bosons as : One can immediately see that in the limit << 1, M 2 Z = M 2 X g 2 X (9u 2 + 4(v 2 1 + 4v 2 2 ) + 9u 2 ). In this model, none of the scalars are Z 2 odd, therefore, in principle, both the CP even and CP odd neutral components mix to give two (4 × 4) mixing mass matrices; one for (h , s 1 , s 2 , H 0 ) and the other for (z , A 1 , A 2 , A 0 ) as given below in eqs. (72) and (73), respectively. We also have a (2 × 2) mixing matrix for the charged scalars (w ± , h ± ) as given in eq. (74).
We therefore require two (4 × 4) rotation matrices (cf. Appendix C) to diagonalize the CP even and CP odd Higgs which we denote by R α and R θ respectively (as shown in eq. (75)) and an orthogonal rotation by angle γ for the charged scalars.
where G z , G z and G ± are the massless Goldstones corresponding to the physical vector bosons Z, Z and W respectively. Notice that in the scalar potential in eq. (64), we have added a higher dimensional symmetry breaking term proportional to λ(= c Λ 2 ) apart from the trilinear term. The relevance of this term can be easily understood from the scalar mass matrix given in eq. (73). In this matrix, the elements of the first and fourth rows and columns are proportional to λ; hence, if we set λ = 0, the resulting mass matrix will be a 2 × 2 matrix with determinant zero, which results in zero-mass pseudo-scalar fields (not allowed). Also, in absence of this term, the U (1) X symmetry can be broken by the vev of H 2 alone, and we don't need the additional singlet scalars.
We denote the angles in R α(θ) by α ij (θ ij ). Thus, we have many unconstrained terms in the rotation matrices R α,θ with at least 6 mixing angles in each, we make the following assumptions to simplify the analysis: i. The mixing angles of h with the singlet scalars are α 12 ≡ α 2 and α 13 ≡ α 3 , respectively. Also, we have not considered very large mixing scenarios.
ii. For simplicity, the mixing angles of H 0 with the singlet scalars are set to zero, i.e α 24 = α 34 ≡ 0. Also, the possibility of mixing between the two singlet scalars has been neglected, i.e α 23 0.
iii. We denote the mixing of the h and H 0 by α 14 ≡ α.
A similar approximation is also considered for the rotation matrix R θ . This helps us to eliminate some of the mixing angles for each of the matrices. Therefore, we are left with the following free parameters: Mixing angles (α, α 2 , α 3 , θ, θ 2 , θ 3 , γ).
The couplings expressed in terms of masses and mixing angles can be found in Appendix D. These model parameters are constrained from both theoretical requirements of unitarity, vacuum stability, perturbativity etc. and experimental data on electroweak observables, Higgs decays and so on. We have assumed small values of α 2 and α 3 so that we can utilize the existing bound on the parameters like α and β of a two Higgs doublet model (2HDM) scenario with and without an additional singlet. For recent analyses of extended 2HDM see [80][81][82][83][84]. It has been shown that large singlet doublet admixture is allowed by the LEP and LHC data [82]. However, a large admixture does not allow a large value for | cos(β − α)| [80]. In our analysis, we have considered the scenarios with tan β ≤ 5 and |cos(β − α)| < 0.1, also, we have assumed sin α 2 ≤ 0.1 and sin α 3 ≤ 0.1. We identify h to be the 125 GeV Higgs boson discovered at the LHC and restrict the parameters in the following range : For the above range of masses, the scale Λ ∼ (300 − 600) GeV for c = −1. Note that tan β > 1 allows only the scenario M h > M H 0 otherwise λ H2 will pick up a very large value. One can have the scenario M h < M H 0 when tan β < 1, however, these choices will lead to the large values of λ H1 , and at the same time top-Yukawa y t >> 1.

DM Phenomenology
In this case, at the leading order, the contributions to the relic abundance and the direct detection cross-section will come from a similar type of annihilation diagrams, as shown in Fig. 7. Hence the true relic abundance is expected to be satisfied only around the resonances of the different scalar and vector mediators. There will be no coannihilations in this case. Apart from M Z and g Z , the other model parameters which will have a dominant role in DM searches . The other parameters which will have a subdominant role are given by M s1 , M A2 , tan β and s α2 . Therefore, we have fixed their values at M s1 = 400 GeV, M A2 = 200 GeV and s α2 = 0.1, respectively. In Fig. 17a, we have shown the variation of the dark matter relic abundance with DM mass for two different values of tan β. The nature of the curve is similar to the one observed in our toy model 1 (see Fig. 8a). When the DM mass is in the sub-GeV range, the Z mediated annihilation will be dominant similar to the previous case. As expected, the current bound on relic density will be satisfied at the DM masses close to the value M Z /2, and at a value M DM < M s2 /2. There are different peaks for M DM > M s2 /2 which correspond to the different resonance annihilation of the DM through the Higgs portal. In all the resonances for M DM > M s2 /2, the relic is much below the present observed abundance. The allowed values of DM mass are mostly limited in the sub-GeV to less than 50 GeV mass. Note that the relic is almost insensitive to the value of tan β. Also, as shown in Fig. 17b   In Fig. 18a, we have shown the regions of M DM allowed by relic density bound and the current experimental limit on the DM direct detection cross section σ SI from XENON 1T. To generate this plot we consider tan β = 2, and the values of the other relevant parameters are the following: 0 < s α3 < 0.01, and 10 ≤ M s2 ≤ 100 GeV. All the other relatively less relevant parameters are fixed at the values as mentioned above. The maximum value of M DM allowed by the data on the relic and σ SI is ∼ 40 GeV. Note from Fig. 18b that the current limit on σ SI put stringent bound FIG. 19: The diagrams which will contribute to muon and electron magnetic moments. on s α3 . For example, for M DM ≈ 30 GeV the allowed value of s α3 can not be larger than 0.01. Here, we have shown the plot for tan β = 2; however, as shown above, the results will be similar for other allowed values of tan β. Like the toy model-I, in this model, the contributions to spin-dependent direct detection cross sections are negligibly small.

Electron Anomalous Magnetic Moment and LFV
a. Magnetic Moments :-In this model, apart from the contribution from a U (1) X gauge particle as has been discussed in sub-section III B, the contributions to the muon and electron magnetic moment will come from the respective diagrams shown in Fig. 19. The contributions from these diagrams from left to the right, respectively, are summarised in the following equations: Here, we have defined |R 2 | 2 ≡ (Y ν 21 ) 2 + (Y ν 23 ) 2 in the same way as we defined R 1 , R 3 in the previous toy model. To do so, we have assumed the same masses for N 1 and N 3 . Note that the contributions in ∆a (H) e is sensitive to the Yukawa coupling |Y 12 |, and the contributions in ∆a µ are coming from the diagrams with ν τ and N 1 /N 3 in the loop, respectively.
The variations of ∆a e with |Y 12 | for different values of M H 0 and M A 0 are shown in Fig. 20a. Note that the contribution to ∆a e is highly suppressed and the values like Y 12 > ∼ 0.01 are allowed. As can be seen from Figs. 20b, the contribution to ∆a µ from the diagram with ν τ in the loop is highly suppressed. In Fig. 21, we have shown the variation of ∆a µ with |R 2 | for different values of M H 0 and M N . Note that the contribution in ∆a µ from the diagram with right-handed neutrinos are significant and have negative values. We have already shown earlier that the contribution from the diagram with X with a mass M X ≈ 0.5 GeV can accommodate the current discrepancy in ∆a µ . The large negative contribution from diagrams with N 1 or N 3 (in the loop) will reduce the value of ∆a µ obtained from a diagram with X. However, note that, for lower values of R 2 i.e. R 2 < ∼ 0.5, the effects are not that significant. Therefore, to explain ∆a µ we can restrict R 2 to a low value. In Fig. 21b, we have shown the variation of the total contribution to ∆a µ with R 2 for the ultimate choices of the other relevant parameters.
b. Lepton Flavour Violation :-In our second model, there won't be any contribution to the processes like τ → µγ, µ → eγ or τ → eγ. However, from eq. (63), one can see that the charge lepton mass matrix is not diagonal and is given by Due to the presence of off-diagonal terms in the charged lepton mass matrix, we have a contribution to the lepton flavour violating decay τ → 3µ as shown in Fig. 22a. Although the contribution is mixing suppressed, the stringent (b) The total contributions in ∆a µ from the diagrams in Fig. 19 and 3. The grey horizontal line is the allowed 3-σ lower limit [1], while the dashed grey line represents the 3-σ lower limit of a very recent estimate [24].
limit on the branching fraction will put a direct constraint on the Yukawa coupling Y 32 since the process occurs at tree level. The upper bound on the branching fraction from the Belle Collaboration [85] is at 90% CL. The amplitude for the process can be written in the form where and the branching fraction is given by [86]  where T τ is the lifetime of the τ lepton, x = 2Eμ/m τ is the reduced energy of the antimuon, and θ is angle between the polarization of the τ and the momentum of the antimuon. In Fig. 22b, we show the variation of the branching fraction of τ → 3µ with Y 32 for three different values of M H 0 , and in each of these cases, we have chosen two different values of tan β. It is evident from the plot that experimental upper limit on B(τ → 3µ) restricts the allowed regions of Y 32 and M H 0 , and the preferable choice is Y 32 The decay width is also sensitive to the value of tan β. Therefore, for all practical purposes it is convenient to set Y 32 to a very small value, say 0.001, in order to evade this strong bound.

Additional contribution to b → c ν
In this case, the diagrams that will contribute to b → c ν (with = e, µ and τ ) decays are given in Figs. 23. The diagrams in Fig. 23a will contribute to b → cµν µ decays, whereas those in Figs. 23b and 23c will contribute to b → ceν e and b → cτν τ decays, respectively. The resulting Wilson coefficient contributing to b → cµν µ can be written as: region is the experimental 1σ allowed region while the dotted lines denote the same in the 2σ CI. The shaded and unshaded colored bands (red,blue,green) signify the theoretical uncertainty in the 1σ and 2σ CI respectively.
where the coupling g HH − W + between the charged Higgs, W-boson and the different neutral scalars are given by : where, R αij is the (ij) element of the rotation matrix given that diagonalises the CP even mass matrix given in Appendix C, and similarly for R θij . The factor ∆ j W HH is given in eq. (43). However, it is evident from the Lagrangian that the dominant contributions would come from the diagrams containing H 0 and A 0 in the loop since their coupling with the neutrinos are not mixing suppressed unlike the other Higgses.
The contributions to b → ceν e and b → cτν τ from the diagrams as mentioned above are negligibly small. The contributions come from the second term of the Yukawa Lagrangian in eq. (63). The Wilson coefficients in these cases are sensitive to the off-diagonal Yukawas of the charged lepton mass mixing matrix. Also, we have already obtained a direct bound on the coupling Y 32 from the branching fraction of τ → 3µ in the previous section and fixed it to 0.001. For such a small Yukawa, the resulting WC will have a value of order O(10 −9 ) which is negligibly small and hence it can safely be neglected. Note that even if we choose Y 32 ≈ 0.01, the contributions to the WC will be of order O(10 −7 ). Similar conclusion holds for the electron final state. Therefore, we will only focus on the contribution to the muon channel.
If we assume, for simplicity, the masses of N 1 and N 3 are equal, then the Wilson coefficient will be proportional to the parameter R 2 which we have already defined in the previous sub-section. We can then constrain the parameter space of R 2 and M N ≡ M N1 = M N3 from the data on R(D), R(D * ) andB → D ( * ) µν µ . In Fig. 24, we show the variation of the branching fractions ofB → D ( * ) µν µ with R 2 for two different values of M N and M H ± each, keeping the masses of the neutral scalars H 0 , A 0 fixed. Also we keep tan β fixed at 2 and choose α such that cos(β − α) = 0.1. We see that even for a large (∼ O(1)) R 2 , the branching fraction ofB → D ( * ) µν µ remains within the 2σ experimental range. Also, we have noted that the change in the branching fraction with the mass of the charged scalar or the heavy neutrino is insignificant (the legends overlap in the figure). We have also studied the impact on R(D * ) and R(D), and the results are shown in Figs 25a and 25b, respectively. Note that it is hard to explain R(D * ) even if we take both the theory and the measured errors within their 2σ confidence interval (CI) 6 . However, we can conveniently explain the observation in R(D) even if we consider the data and theory at their 1σ CI. Also, even though large values of R 2 is allowed by the branching fractions to the muon mode, the data on R(D ( * ) ) restricts the value of R 2 to < ∼ 0.5 for the entire range of RHN or charged scalar mass.

Neutrino Mass Generation
In this case the minimal seesaw mechanism will help us give rise to the neutrino mass. From the Yukawa interactions eq. (63), one can obtain the Dirac neutrino mass matrix M D and the Majorana mass matrix M R to be : Thus the light active neutrino masses can be obtained from the seesaw formula given by: By using the structure of M D , M R in the type I seesaw [87][88][89][90] formula for light neutrino masses mentioned above, we find a general structure of light neutrino mass matrix without any textures unlike that found in toy model I. However, the light neutrino mass matrix has rank 2 predicting the lightest neutrino mass to be vanishing. While neutrino oscillation experiments can not constrain such a scenario, other experiments like neutrinoless double beta decay which is sensitive to absolute neutrino mass scale can shed more light on such scenario in future.

VI. KOTO ANOMALY
Very recently, the KOTO experiment at J-PARC reported an excess of events over the SM expectation for the rare decay process K L → π 0 νν [9]. They have reported four candidate events in the relevant signal region, whereas the SM expectation is only 0.1 ± 0.02 events. The corresponding measured value is given by where the quoted errors are given at 68% and 95% (within the parenthesis) CI, respectively. The SM prediction, on the other hand, is [91,92]  which is about two orders of magnitude smaller than the one measured at KOTO. An upper bound on the branching fraction of K + → π + νν decay reported by the NA62 Collaboration [93,94]: Here again, we have quoted the 90%(95%) confidence level (CL) limit and the measured value is consistent with the respective SM prediction of (9.11 ± 0.72) × 10 −11 . The explanation of these excess events in K L → π 0 νν require new contribution in d → s + invisible decays beyond the SM. However, the same NP will contribute to K + → π + νν decay as well. Also, the above mentioned branching fractions follow a model-independent bound [95], which is totally based on isospin symmetry. From the measured values, it looks difficult to explain both the branching fractions simultaneously. However, it is important to note that the experimental measurement at NA62 [93] excludes the following kinematic regions of the missing mass : 0.01 < M 2 miss < 0.026 GeV 2 /c 2 and M 2 miss > 0.068 GeV 2 /c 2 , respectively. Hence, for the missing masses within these kinematic regions, one can avoid the bound given in eq. (93). Following this, there are different NP explanations available in the literature, for example, see [96][97][98][99][100][101].
The presence of a low scale U (1) X gauge boson in our model allows us to study the prospect of explaining the KOTO excess via K L → π 0 Z (Z → νν) decays which is possible in both of our toy models. The relevant Feynman diagram is given in Fig. 26a. The corresponding branching fraction can be expressed as i.e via the resonance production of Z and then by a subsequent Z → νν decay. The branching fraction of the rare K L → π 0 Z decay is given by [99] B The effective vertex for the interaction [d L γ µ s L ]Z is given by which is obtained after calculating the loop diagram given in Fig. 26a. Here, f K L π 0 + is the meson decay form factor, λ(x, y, z) = x 2 + y 2 + z 2 − 2(xy + yz + zx), and C(x t ) can be found in Eq. (14). Note that in our model Z plays the role of a mediator of DM (N 2 ) interactions, and the bound on the relic abundance is satisfied at the DM mass M N2 ≈ M Z /2. Since in our model Z primarily decays to µ + µ − , e + e − and νν, therefore the resonance production of Z and its subsequent on-shell decay to N 2 N 2 will be kinematically forbidden. As we have seen earlier the other low energy data like ∆a µ , data in b → sµµ, b → see etc. preferred a mass region 0.22 < ∼ M Z < ∼ 1 (in GeV) for a given coupling 0.0005 < ∼ g Z < ∼ 0.001, for example see Fig. 6. In Fig. 26b, we have shown the variation of B(K L → π 0 νν) with M Z for different values of the mixing parameter . Note that the coupling g Z cancels in the ratio Γ(Z →νν) Γ Z , hence the branching fraction is insensitive to the variation of g Z . As can be seen from the plot, for mass range 0.30 < M Z < 0.35 GeV, we are able to explain the observed branching excess in B(K L → π 0 νν). Note that in this region we can avoid the model-independent bound given in eq. (93). Though we have discussed the phenomenology of our toy models with M Z = 1 GeV, the observations made are equally valid for the mass window 0.3 < ∼ M Z (GeV) < 1.

A. Higgs Invisible decays
In any NP model it is quite an exciting prospect to look into the non-standard or undetected decays of the SM Higgs as a complementary search for BSM particles. The toy models that we discussed above constitute of a dark matter particle that couple to the SM Higgs boson through its mixing with singlet scalars. Also, there is a viable parameter space in both the models where the dark matter mass is lighter than the Higgs. Under such a scenario, it would be a useful exercise to study the model contribution to Higgs invisible decay and use the available data on it to further constrain the model parameters. All the relevant diagrams contributing to the Higgs invisible decay are shown in Fig. 27. From the figure it is almost clear that the dominating invisible decay of Higgs would be to the DM N 2 for M N2 < M h /2 and the decay to 4ν via gauge bosons would in general be suppressed compared to this tree level decay. However, since the new gauge boson in our model, Z , is light (sub-GeV), the contribution mediated by Z could also be quite significant. The third decay, involving one heavy boson Z and a light Z , is in general suppressed by the very small value of the gauge mixing parameter (∼ 10 −4 ) that we considered. Hence, one can safely neglect the contribution from this decay mode. Therefore we will only consider the first two decay modes in our calculation. Both ATLAS and CMS have looked into such invisibly decaying Higgs mainly through its inclusive production in the vector boson fusion mode, as well as in the associated production of a Higgs with a Z boson. The constraint on the Higgs invisible decay branching fraction from the ATLAS experiment at LHC is [102] while the recent ATLAST announcement [103] puts a more stringent constraint at 13%. The Higgs decay to SM particles is known to be around 4 MeV. In the following two subsections, we will discuss the impact of this upper limit on the model parameters, in particular the mixing angles.

Toy Model 1
The invisible decay width of Higgs to dark matter is given by : On the other hand, the decay of Higgs to SM neutrinos via XX, or equivalently via Z Z is given by : where g hZ Z = 8g 2 X (c α2 s α1 v 1 + 4s α2 v 2 ) is the effective coupling of SM Higgs with Z via mixing with the singlet scalars. For our model, [B(Z → 2ν )] 2 ≈ 0.14 for M Z = 1 GeV and g Z = 10 −3 . Therefore, the total invisible decay width of Higgs is given by the sum of the decay widths as mentioned above. In Fig. 28, we have shown the variation of the Higgs invisible decay with the mass of dark matter for different values of the mixing angles. Note that only small mixing like s α2 = 0.01, s α1 = 0.01 are allowed by the current limit for the entire mass range of N 2 . However, s α3 could be as large as 0.1.

Toy Model 2
In our second model, the contribution to the Higgs invisible decay width is given by Also, the decay of Higgs to SM neutrinos via Z Z is given by : where g hZ Z = 8g 2 X (R 12 v 1 + 4R 13 v 2 ) is the effective coupling of SM Higgs with Z via mixing with the singlet scalars. Again, the sum of squares of branching fraction, [B(Z → 2ν )] 2 ≈ 0.14 for M Z = 1 GeV as mentioned before. In Fig. 29, we have shown the dependencies of the B(h → invisible) with the DM mass and other relevant parameters in toy model 2, like sine of the mixing angles and tan β. Note that our chosen benchmark values like s α2 = 0.01, s α3 = 0.01 and tan β = 2 or 4 are allowed by the current bound on B(h → invisible).

B. LFV decays of Higgs
It is evident from the Yukawa interactions in eqs. (20) and (63) that there exists lepton flavour violating decays of the Higgs (h) for both the Toy models. However, there are notable differences between the allowed LFV channels in the two models. The U (1) X charge assignments of the charged leptons are such that in both the models we will get the h → τ e decay, for example, see   So far no excess have been observed in these channels at the LHC searches and the most recent upper limits on the lepton flavour violating branching fractions of the Higgs boson by CMS [104] reads These limits will be helpful to constrain the lepton flavour violating Higgs couplings Y ij where (i, j = 1, 3) & i = j. The general expression for the branching fraction for the LFV Higgs decay is given by Here, Γ tot h = Γ(h → SM ) + Γ(h → Invisible) and y denotes the effective LFV coupling. In most of the allowed parameter spaces, we can expect Γ(h → Invisible) << Γ(h → SM ); however, there are regions where it might be relevant to consider. In both the models, for h → τ e decays the effective LFV coupling is given by y ≡ Y 13 2 + Y 31 2 .
As mentioned earlier, there won't be any contribution to h → µτ or h → µe decays in Toy Model I. In Toy Model II, the expression for the braching fraction for h → µτ decay is given by which in the limit α → 0, gives us the corresponding expression for h → eτ decay. We will obtain the expression for B(h → µe) after replacing Y 32 by Y 12 in eq. (104). In Fig. 31a we have shown the variation of B(h → eτ ) with the effective coupling y. Since the contribution will be similar for both the models, we have not shown it separately for the two. It can be clearly understood from the plot that the coupling cannot be larger than ∼ 0.005 irrespective of the value of tan β or other angles. There is very little dependence on s α2 or s α3 which is coming from the contributions in Γ(h → Invisible) (see eq. (98)) in the denominator. For illustrative purpose, we have shown the variation for M N2 = 20 GeV; however, we have checked that the variation does not change significantly on changing the DM mass.
In Fig. 31b we show a similar variation of the branching ratio to the µτ mode with Y 32 for two different values of tan β and other mixing angles. As mentioned earlier, this decay mode is specific for Toy Model II only. Since this process is mixing induced, both tan β and Y 32 are tightly constrained from the data. As expected, the branching fraction is sensitive to both the mixing parameters β and α. Note that for tan β = 2, the allowed values of Y 32 is Y 32 < 0.01. However, for tan β > 2, more higher values of the Yukawa coupling are allowed. In general, higher values of the tan β prefers higher values of the Yukawa coupling. This is expected since the constraint |cos(β − α)| = 0.1 implies smaller sin α for large tan β. Once again the conclusions are not affected significantly by the DM mass. Also, we have noted that for values like tan β = 2 and s α2 = s α3 = 0.01, the branching fraction B(h → µe) ≤ 0.4% for Y 12 ≤ 0.01.  In this subsection, we would like to briefly discuss some other possible collider signatures of our Toy models at the LHC, in addition to the specific ones mentioned above. Here, we will only mention a few exciting channels for the search at the LHC; a detailed analysis is beyond the scope of this paper. It is quite evident that the most notable collider signals would be the multi-lepton final states with or without an associated missing energy ( / E T ). More importantly, the dilepton (2 ) + / E T channel which can be probed with excellent precision in the high luminosity     Table. IV. Please note that the cross-section for the dilepton + / E T channel for Toy model I quoted here excludes the contribution from the charged Higgs mediated diagram (which we have separately shown in Table. V).
colliders in the near future, has the potential to discriminate the two Toy models. A few of the dominating production channels are mentioned in Table IV. Note that the search for 4 final states will not be a unique test of our toy models. This is because, in all the cases, the primary decay channels are via the production of Z which is very common in NP models with an additional U (1) X gauge bosons. The details of the collider searches for such an extension with a sub-GeV M Z can be seen from [52]. However, the search of (2 + / E T ) could be helpful to probe Toy model I since the process is also mediated by the production and decay of the new scalars. In Table IV, we have mentioned only the dominating production channels. We have noted that for the allowed values of the model parameters, as discussed earlier, the production cross-section σ 2 + / E T in Toy Model-I is much larger than that in Toy Model-II. This is due to the presence of an additional channel pp → H + H − → µ + µ − / E T in Toy Model-I. The corresponding production crosssection of this specific channel for some benchmark scenarios are given in Table V. Here, the events are generated in MADGRAPH [105] at √ s = 14 TeV. Depending on the charged Higgs mass and associated coupling Y 22 , the cross-section can be quite large. For the rest of three channels as mentioned in Table IV, the estimated production cross-sections for a few benchmark scenarios are given in Table VI. As one can see, the production cross-sections are very small, which is expected since the diagrams mentioned here are mostly mixing induced which we have assumed to be small. Note that the cross-sections are insensitive to the mass of Z . We have checked that in both the toy models, among these three channels the dominating contribution will come from pp → Z → Z s 2 → + − / E T . Therefore, this channel will not be helpful to discriminate the signatures of Toy Model-I from that of Toy Model-II. However, as one can see that σ(pp → H + H − → µ + µ − / E T ) is much larger than the production cross-sections for the rest of three channels, therefore, at the colliders a dedicated search for µ + µ − + / E T signature could be helpful to probe Toy Model-I. Please note than in order to obtain long-lived charged scalars whose decay length (cτ ) > ∼ 0.1 mm, we need Yukawa coupling of the order of ∼ (10 −6 − 10 −5 ). This will give very low production cross-section of the final state and therefore techniques like displaced muon and kink vertex will not be applicable. So we do not discuss it further.
Another interesting collider signature could be the production of (µ ± + / E T ) which is an exclusive feature of Toy Model I and can be a smoking gun signal. It can be observed at a pp collider like the LHC where the intermediate particles leading to such a final state are the inert charged Higgs (H ± ) and inert neutral scalars (H 0 , A 0 ) as shown in Fig. 32. The readers may recall that the inert Higgs couples to muon along with the dark matter. Due to electroweak interaction, it is possible to have sizeable production of H + -H 0 which then decay to give a mono-muon plus missing energy final state. There are no other contributing diagrams to this muon specific signal. The other mono-lepton channels (say the mono-electon for example) will be kinematically suppressed due to the associated heavy neutrinos in the final state. The mono-muon signal is cleaner than the mono-jet searches and therefore it is possible to tag the muon. The major background is the W ( ν) process but one can expect a clean signal away from the W -boson mass window. The other minor backgrounds include t, tt, Z/γ * ( ), γ + jets and V V (where V stands for the SM vector bosons W, Z). As can be seen from Table. VII, for a few suitable benchmark values of the scalar and DM masses and coupling Y 22 , it is possible to obtain a few femtobarns of production cross-section. In Toy Model-II, the U (1) X charge of H 2 forbids its coupling with muon and N 2 simultaneously, instead it couples with the other RHNs (N 1 , N 3 ). Also the scalars are in general lighter than N 1 , N 3 . Therefore, once again, the µ ± + / E T production at LHC will be    kinematically suppressed. The probable collider signatures of Toy Model-II will be the productions of τ + µ − µ − e + , τ − µ + µ + e − and τ τ µe events at the LHC via the production and decay of H 0 H 0 . This is possible only in Toy Model-II since in this model, H 0 takes part in LFV interactions, which is not allowed in Toy Model-I. In a few benchmark scenarios, the corresponding production cross-sections are given in Table VIII. As expected, the cross-sections are highly sensitive to the mass of H 0 . Note that for the above mentioned four-lepton final states the SM background will be highly suppressed. Therefore, a dedicated search of these four lepton states with specific flavour and charge could be useful to test our Toy Model-II.

VIII. SUMMARY
We have extended the SM by an Abelian U (1) X gauge group which results in a massive gauge boson (X) that couples only to leptons and has a small kinetic mixing with the SM Z boson. We have considered only the low masses of X (M X < ∼ 1 GeV). In this kind of extension, we will get new contributions to flavour changing processes like b → s + − decays, and the new contribution will be in ∆C 9 which is the WC of the operator O 9 . Here, O 9 is a left-handed quark current operator with vector muon/electron coupling. Also, in this model, the contributions such flavour changing processes will be in both the electron and muon final states. At the same time, we will get new contributions in anomalous magnetic moments of the muon. We use the present data on R(K), R(K * ), the ratio of branching fraction B(B 0 → K * 0 χ(µ + µ − ))/B(B 0 → K * 0 µ + µ − ), B(B → K ( * ) e + e − ), and muon anomalous magnetic moment to constrain U (1) X charges of the SM leptons. Also, the values ∆C µ 9 and ∆C e 9 which obtained from the analysis are consistent with the global fit results of the data in b → s + − decays including various angular observables. Additionally, we consider all upper bounds from different experimental data on such light Abelian gauge boson mass and its couplings. Now charging the SM fermion under a generic U (1) X symmetry makes the theory anomalous. To get an anomaly-free renormalisable model, we have incorporated additional chiral fermions into the model. In order to fit our requirements with a minimal particle content, we have considered a scenario where the three generation of leptons having vector type interactions U (1) X have corresponding charges (n 1 , n 2 , n 3 ) = (−1, 2, −1) respectively. Such a choice is consistent with the data and also ensures anomaly cancellations after adding one right-handed neutrino per fermion generation having equal, and opposite U (1) X charges as that of SM lepton in that generation. Also, we have added additional singlet and doublet Higgs fields to get the desired mass spectrum. To prevent a direct coupling of the RHNs with the lepton doublets via SM Higgs, we impose a discrete Z 2 symmetry on the particles in two different ways which lead to two distinct models and phenomenology. This kind of symmetry restrictions will provide a natural candidate for DM in our extended models. At the same time, the chosen particle content of the models can also generate light neutrino masses, in agreement with neutrino oscillation data.
We are able to successfully study the DM phenomenology which is almost similar for the two Toy models but they have different neutrino mass generation mechanisms. The scalar content is very rich with an additional scalar doublet and two scalar singlets apart from the usual SM Higgs doublet. However, the second scalar doublet has very distinct features and plays different roles in each of the two Toy models. The low gauge bisin X mass allows us to evade stringent constraints from LHC, while facing tight constraints from other low energy experiments. We ensure that our analysis is consistent with the LEP II bounds on U (1) X gauge boson mass and coupling and the bounds from other light boson search experiments. Few preliminary results have also been shown and discussed.
The two toy models lead to phenomenological implications that can be tested at the collider experiments. The promising channels are the 4-lepton final states, dilepton (2 ) + / E T etc. Also, in the low energy experiments the potential signatures may come from the FCNC processes, like b → s(d)+invisible, s → d+invisible, c → u+invisible which will lead to rare decays of B q /K/D mesons to a relatively lighter meson final state with invisible particles. For example, the recently observed excess of events in the K L → π 0 + invisible decay at the KOTO experiment can be explained for a mass window 0.30 < M Z (M X ) < 0.35 GeV. For these values of M Z , all the other observables discussed in this article are consistent with the corresponding data. Also, both the models will contribute to LFV h → τ e decays.
We have discussed some distinct features of both the models which could be helpful to discriminate the signatures of the two models at different experiments. At the LHC, the production of µ ± + / E T and lepton-specific µ + µ − + / E T events could be the possible signatures of Toy Model-I which is not possible to get in Toy Model-II. On the other hand, the search for the specific multi-leptonic states like τ + µ − µ − e + , τ − µ + µ + e − and τ τ µe could be useful to identify the potential signatures of Toy Model-II. In the context of Higgs LFV, Toy Model-II contributes to h → τ µ and h → µe decays while Toy Model-I does not. Similarly, there are a few examples of the potential observables in the low energy sector: Toy Model-I contributes significantly to the semileptonic or purely leptonic decays B q /K/D mesons via the following quark level transitions: b → c(u)τν τ , s → uτν τ , c → sτν τ . However, Toy Model-II does not have significant contributions to these decays with τ in the final states. In Toy Model-II, we do not have any contributions in τ → eγ, τ → µγ and µ → eγ LFV decays, however, this model contributes to τ → 3µ decay. On the other hand, Toy Model-I contributes only in τ → eγ, not in all the other LFV decays as considered above. More precise data from future experiments will be able to discriminate between such toy models while confirming or ruling out some part of the available parameter space. where s ij ≡ sin α ij and c ij ≡ cos α ij .