Reconciling dark matter, $R_{K^{(*)}}$ anomalies and $(g-2)_{\mu}$ in an ${L_{\mu}-L_{\tau}}$ scenario

We propose an anomaly free unified scenario by invocation of an extra local ${\rm U(1)}_{L_{\mu}-L_{\tau}}$ gauge symmetry. This scenario simultaneously resolves the $R_{K^{(*)}}$ anomalies, the dark matter puzzle and the long-standing discrepancy in muon's anomalous magnetic moment. A complex scalar ($\eta$) having nonzero ${L_{\mu}-L_{\tau}}$ charge has been introduced to break this new U(1) symmetry spontaneously. Moreover, for the purpose of studying dark matter phenomenology and $R_{K^{(*)}}$ anomalies in a correlated manner, we introduce an inert ${\rm SU(2)}_L$ scalar doublet ($\Phi$), a $\mathbb{Z}_2$-odd real singlet scalar ($S$) and a $\mathbb{Z}_2$-odd coloured fermion ($\chi$) which transforms vectorially under the ${\rm U(1)}_{L_{\mu}-L_{\tau}}$ symmetry. This extra gauge symmetry provides a new gauge boson $Z_{\mu\tau}$ which not only gives additional contribution to both $b\to s \ell\ell$ transition and $(g-2)_{\mu}$ but also provides a crucial annihilation channel for dark matter candidate $\rho_1$ of the present scenario. This $\rho_1$ is an admixture of CP-even neutral component of $\Phi$ and $S$. Our analysis shows that the low mass dark matter regime ($M_{\rho_1}\lesssim 60$ GeV) is still allowed by the experiments like XENON1T, LHC (via Higgs invisible branching) and Fermi-LAT, making the dark matter phenomenology drastically different from the standard Inert Doublet and the Scalar Singlet models. Furthermore, the present model is also fairly consistent with the observed branching ratio of $B\to X_s\gamma$ in $3\sigma$ range and is quite capable of explaining neutrino masses and mixings via Type-I seesaw mechanism if we add three right handed neutrinos in the particle spectrum. Finally, we use the latest ATLAS data of non-observation of a resonant $\ell^+\ell^-$ signal at the 13 TeV LHC to constrain the mass-coupling plane of $Z_{\mu\tau}$.


I. INTRODUCTION
With the discovery of the missing piece, the Higgs boson at the Large Hadron Collider (LHC) [1,2] at CERN the Standard Model (SM) of particle physics has been turned into a complete theory. From the last several decades it has been a well known fact that most of the theoretical predictions of this theory are in good agreement with various experimental results. However, at the same time, different experimental results in various directions compelling us to formulate physics beyond the SM (BSM). For example, dark matter relic density has been measured with a great precision from the temperature and polarization anisotropies of the cosmic microwave background (CMB) radiation by experiments like WMAP [3] and Planck [4]. On top of that various indirect evidence such as rotation curve [5], gravitational lensing of distant objects [6], collision between galaxy clusters (such as Bullet cluster [7] etc.) etc. have strongly support for the existence of dark matter. However, in the SM there is no such candidate for dark matter. On the other hand neutrino oscillation experiments [8][9][10][11] firmly established massive nature of at least two neutrinos and have accurately measured three intergenerational mixing angles, both of which are missing in the SM due to non-existence of the right handed counterparts of left handed neutrinos. Besides, the CP-violation in the quark sector is not at all sufficient to explain the observed baryon asymmetry of the Universe [12]. Furthermore, there is an enduring ∼ 3.5σ discrepancy [12] between experimentally measured value of the anomalous magnetic moment of muon [(g − 2) µ ] and its SM predictions, which strongly indicates the presence of a new physics (NP) beyond the SM.
Apart from the above mentioned facts, over the last few years different flavour physics experiments like LHCb, Belle and Babar have been consistently shown that experimental data for different observables are in significant disagreement with respect to the corresponding SM predictions. Indeed this situation demands the invocation of NP effects. Recently the LHCb collaboration has reported additional hints for violation of Lepton Flavour Universality (LFU) between b → sµ + µ − and b → se + e − processes. The LFU violation 1 (LFUV) can be measured with the help of following observables R K and R K * Summary of the corresponding experimental results with their SM predictions for different di-lepton invariant mass squared (q 2 ) ranges are given in Table I.
For the purpose of explaining b → sµ + µ − anomaly, this type of U(1) Lµ−Lτ model has also been modified from its minimal version, albeit in a different approach [31,43,50,51,[96][97][98][99][100][101]. In the present article, we introduce a Z 2 -odd bottom quark like non-standard fermion field χ which is vectorial in nature under the U(1) Lµ−Lτ symmetry. It couples to all generations of the down-type SM quarks via Yukawa like interaction involving a Z 2 -odd scalar doublet Φ. Moreover, we introduce a Z 2 -odd singlet scalar S which helps us to explain the flavour anomaly, dark matter and (g − 2) µ anomaly simultaneously. A Z 2 -even complex scalar singlet field η with a nonzero L µ −L τ charge has been introduced for the purpose of breaking of U(1) Lµ−Lτ symmetry spontaneously. Apart from these fields we have the usual Higgs doublet field H which breaks the SU(2) L ×U(1) Y symmetry. Therefore, in the broken phase of both electroweak (SU(2) L ×U(1) Y ) and U(1) Lµ−Lτ symmetries, we have three physical Z 2 -odd neutral scalars emerge from the mixing between Φ and S. The lightest field among the three physical Z 2 -odd neutral scalars can be considered as a potentially viable dark matter candidate. This is an admixture of both doublet and singlet scalar representations and have distinct phenomenology compared to the standard Inert Doublet [102][103][104] and the Scalar Singlet models [105][106][107][108], where the low mass dark matter regime is almost ruled by the latest bound on spin independent scattering cross section from XENON1T [109] as well as by the upper limit on Higgs invisible branching fraction from LHC [110]. This is mainly due to the fact that in these models in the low mass regime (M DM ≤ 62.5 GeV), dark matter candidate predominantly annihilates into bb final state.
On the contrary, in the present scenario, the dark matter candidate in the low mass regime can annihilate into a pair of L µ − L τ gauge boson Z µτ and the branching fraction of this annihilation channel is controlled by dark sector mixing angle θ D . This actually makes the dark matter freeze-out process extremely correlated with the flavour physics anomalies and (g − 2) µ anomaly, where an O(MeV) light Z µτ plays a pivotal role. Since, Z µτ does not have direct couplings to the first generation leptons and quarks, constraints from the LEP and more recently from the LHC on the g Zµτ − M Zµτ plane are relatively relaxed. Particularly, light gauge boson with M Zµτ < ∼ 100 MeV and also with moderate gauge coupling g Zµτ < ∼ 10 −3 is still allowed from the experiments measuring neutrino trident processes namely CCFR [111] CHARM-II [112]. Moreover, apart from (g − 2) µ anomaly and flavour physics related issues, such a light gauge boson has excellent cosmological implication. The reason is that it can relax the ∼ 3σ tension between the measurements of Hubble constant (H 0 ) from two different epochs 3 by providing extra contribution to the radiation energy density (∆N ef f ∼ 0.2−0.5) through the alteration of neutrino decoupling temperature [114]. In the present scenario the non-standard neutral gauge boson Z µτ emerge from all three neutral gauge bosons associated with SU(2) L , U(1) Y and U(1) Lµ−Lτ gauge groups by diagonalising a 3 × 3 mixing matrix. The additional contribution to the anomalous magnetic moment of muon comes from an effective µ + µ − γ vertex which has been generated from one loop penguin diagram involving Z µτ . Moreover, we also have one loop contribution from a diagram involving other BSM scalar (an orthogonal state of the SM-like Higgs boson arises from the mixing between H and η in the broken phase of the theory). However, its effect on (g − 2) µ is negligibly small.
To this end, we would like to mention another novel signature of the present scenario. The 3 At two different redshifts (z), one is from the CMB experiment Planck [4] at high z while another one is from the local measurement using Hubble Space Telescope [113] at low z. correlation between dark sector and flavour physics sector is not only due to L µ − L τ gauge boson but also due to all the Z 2 -odd neutral particles (including dark matter candidate of the present scenario) along with the coloured Z 2 -odd fermion χ generate non-standard one loop contributions to produce b → sµ + µ − transition. In the present scenario, one can produce nonstandard contributions to both the WCs C µ 9 and C µ 10 respectively, however, the contribution of the latter is insignificant and hence our analysis will be furnished with C NP,µ 9 only. The NP contribution to C µ 9 is obtained from non-standard penguin and self-energy diagrams and there is no further NP contribution from box-diagram at one loop level. Moreover, we consider the constraint from the branching ratio of another flavour changing neutral current (FCNC) process B → X s γ. Hence, we have computed the branching ratio of this decay in the present scenario. Further, neutrino masses and mixings can easily be addressed in these class of L µ − L τ models via Type-I seesaw mechanism by adding three right handed neutrinos, which are singlet under the SM gauge groups and two of them have equal and opposite L µ − L τ charges for anomaly cancellation. Since a detailed analysis on neutrino masses and mixings in the present scenario is beyond the scope of this article, hence for the sake of completeness, we just have added three right handed neutrinos in the Lagrangian and find the Majorana mass matrix for the light neutrinos. A more comprehensive analysis on diagonalisation of the light neutrino mass matrix and thereby finding the mass eigenvalues and mixing angles in the L µ − L τ scenario has already been done in [84].
Finally, in order to impose the constraints on the parameter space of the present scenario from the LHC experiment, we use the latest ATLAS data [115] of non-observation of a resonant + − signal at the LHC running at 13 TeV for the high mass range of Z µτ . Hence, we will estimate the cross section for the process pp → Z µτ → + − at the 13 TeV LHC in the present scenario. Consequently, it will be an interesting part of this exercise that, how the LHC data can constrain the values of non-standard gauge coupling constant as well as the Z − Z µτ mixing angle.
The article is organised as follows. In Sec. II we introduce the model with possible field content and interactions as well as we set our notations. Then in Sec. III, we show the calculational details of flavour physics observables and after that we will discuss (g − 2) µ anomaly in Sec. IV. In Sec. V, we show the viability of our dark matter candidate of the present scenario considering all possible bounds from ongoing experiments and explain how can we correlate the dark matter with the flavour physics anomalies. We briefly discuss neutrino mass generation via Type-I seesaw mechanism in Sec. VI. Sec. VII deals with constraint that are obtained from non-observation of a resonant + − signal at the LHC running at 13 TeV. Finally, we summarize our results and conclude in Sec. VIII.

II. THE U(1) Lµ−Lτ MODEL
In order to facilitate our motivations (discussed in Section I), we propose an anomaly free U(1) Lµ−Lτ gauge extension of the SM. This scenario is free from mixed gauge-gravitational and axial vector gauge anomalies because these anomalies cancel between second and third generations of charged leptons and their corresponding neutrinos due to their equal and opposite L µ − L τ charges. The Lagrangian which remains invariant under the SU(3) C ×SU(2) L ×U(1) Y × U(1) Lµ−Lτ × Z 2 symmetry is given by, are field strength tensors for the two U(1) gauge fields 4B α andX α respectively while the Lorentz indices α, β ≡ 0, 1 . . . 3. The term contains both field strength tensors is the kinetic mixing term betweenB α andX α , which is not forbidden by any of the symmetries of the present model. Full list of particle contents and their quantum numbers under various symmetry groups are given in Table II.

Gauge groups
Fermion fields Scalar fields Quark fields Lepton fields Z 2 symmetry + + + + + + + + + + + + -+ + -- As has been discussed in earlier that the L µ − L τ extension of the SM is anomaly free, however, for the purpose of neutrino mass generation via Type-I seesaw mechanism we invoke three SM gauge singlet right handed neutrinos (N Ri ) having nonzero L µ − L τ charge in such a manner so that their inclusion does not introduce any further anomaly. The Lagrangian of right handed neutrinos is denoted by L N which contains kinetic energy terms, mass terms and Yukawa terms associated with the SM lepton doublets (L Li ) allowed by the symmetries of the present model.
whereH = i σ 2 H * . M ee , M µτ are the bare mass parameters while y eµ , y eτ and y i are the dimensionless Yukawa couplings. In order to generate b → s transition at one loop level involving Z 2 -odd scalars a non-standard SU(2) L singlet fermionic field χ with a colour charge has been introduced in this scenario. This fermion is also Z 2 -odd and has an electric charge identical to SM down-type quarks. Furthermore, both left and right chiral parts of χ field have same L µ − L τ charge making it a vector like fermion under U(1) Lµ−Lτ symmetry. The Lagrangian of this field is given by where M χ is the bare mass parameter for the χ field and f j s are couplings of the Yukawa type interactions among the SM quark doublets (Q Lj ), Z 2 -odd scalar doublet Φ and the right chiral part of χ. The above Yukawa type interactions terms involving s and b quarks have significant roles in b → s transition and hence in the explanation of R K ( * ) anomalies. The covariant derivative D α for the field χ is defined as where g 1 , g Zµτ and g 3 are the U(1) Y , U(1) Lµ−Lτ and SU(3) C gauge coupling constants respectively. n χ is the L µ − L τ charge of χ. Further, Λ a s (a = 1, 2 . . . 8) are the eight Gell-Mann matrices representing the generators for SU(3) C while the corresponding gauge fields are denoted by G a α . The 4 th , 5 th and 6 th terms of the Eq. (2) represent the kinetic terms for all the non-standard scalar representations (η, Φ and S) introduced in the present model for specific purposes. Particularly, the complex singlet (under the SM gauge group) scalar η is necessary to break the U(1) Lµ−Lτ symmetry spontaneously as it is the only scalar field which has not only a U(1) Lµ−Lτ charge but also has a nonzero vacuum expectation value (VEV) v 2 . Consequently, after L µ − L τ symmetry breaking one obtains a massive non-standard neutral gauge boson. It has played crucial roles in different aspects: e.g., (g − 2) µ anomaly explanation, amelioration of the anomalies that are related to b → sµµ transition and most importantly it provides new annihilation channels for the dark matter candidate, which alters its dynamics from the standard case. Moreover, a Z 2 -odd SU(2) L scalar doublet Φ having both U(1) Y and U(1) Lµ−Lτ charges which are required to get the NP contribution to b → s transition via the Yukawa like interaction given in Eq. (5). Although, one of the neutral components of Φ (lightest one) is stable, but for the simultaneous explanation of the dark matter enigma, (g − 2) µ anomaly and R K ( * ) anomalies we include another real singlet scalar field S which is also odd under Z 2 symmetry. Covariant derivatives for the scalar fields η and Φ are given as follows where σ a are the three Pauli's spin matrices with a runs from 1 to 3. n X denotes the L µ − L τ charge of the corresponding scalar fields X = Φ, η. Further, g 2 is the SU(2) L gauge coupling constant and W a α s are the corresponding gauge bosons. Finally, the scalar potential V (H, η, Φ, S) in Eq. (2) contains those interactions terms among the scalar fields which remain invariant under all the symmetries of the present model, has the following form, where m H , m η , m Φ and m S are real parameters having dimension of mass and λ i s (i = H, η, S, 1, 2 . . . 7) are dimension less, real quartic coupling constants because the corresponding operators are self-conjugate in nature. However, the quartic coupling λ 8 can in general be a complex parameter and thus can act as an extra source of CP-violation. Since in this work we are not studying any CP-violating effects, we have taken λ 8 as a real parameter and this assumption will not alter our conclusions. Although, the term proportional to λ 8 has important significance in this model as it generates mixing between Φ and S. Later we will discuss more elaborately on this issue. The component wise structure of the scalar fields are given in the following where v 1 and v 2 are the VEVs of the scalar fields 5 H and η respectively. After breaking of both electroweak and L µ − L τ symmetries by the respective VEVs v 1 and v 2 , one can have mixing between the real components h 1 and h 2 due to the presence of an interaction term proportional to λ 1 in V (H, η, Φ, S). The mixing matrix in the basis 1 scalar ρ 3 exactly coincides with a 0 . In matrix notation, the basis transformation can be shown as where the mixing angle θ D can be expressed in terms of parameters of the Lagrangian as, Among the three states (ρ 1 , ρ 1 and ρ 3 ), we choose ρ 1 as the lightest odd particle (LOP) which is regarded as the stable dark matter candidate in this scenario. Thus, the dark matter candidate in this scenario is an admixture of singlet and doublet states. The expressions for the masses of these Z 2 -odd scalar fields are given below where Further using Eqs. (19)(20)(21), one can establish a relation between M ρ 1 , M ρ 2 , M ρ 3 and θ D , which has the following form Therefore, the mass of the CP-odd scalar ρ 3 is not an independent quantity in the present scenario and it becomes fixed though the above relation once we know other parameters like M ρ 1 , M ρ 2 and θ D . This is a consequence of that, the 2 × 2 and 3 × 3 elements of the dark sector mixing matrix M 2 DM are identical. From the symmetry argument this can be understood as follows. The splitting between the coefficients of φ 0 2 (∝ 2 × 2 element of M 2 DM ) and a 0 2 (∝ 3 × 3 element of M 2 DM ) of a Z 2 -odd doublet Φ is obtained from a term like (H † Φ) 2 (usual λ 5 term in the Inert Doublet Model [102]), which is forbidden here by the U(1) Lµ−Lτ symmetry invariance. Additionally, in the dark sector we also have a charged scalar φ ± and its mass term is given by Let us now find out the effects of the extra U(1) Lµ−Lτ local gauge symmetry on the gauge sector and generate the physical states of the gauge bosons with their proper mass terms. In the Eq. (3),B α andX α are denoted as gauge fields corresponding to gauge groups U(1) Y and U(1) Lµ−Lτ respectively. As mentioned earlier, the kinetic terms for the two U(1) gauge fields with hat notation are not diagonal and it is clearly evident from the presence of a mixing term between two U(1) gauge fields proportional . The kinetic mixing parameters is severely constrained from the electroweak precision data (sensitive mainly in the low mass regime of the extra gauge boson) [116,117] and also from di-lepton searches at the LHC (for relatively high mass regime i.e., few hundred GeV to few TeV range). Now, one can perform a basis transformation from "hat" states to "un-hat" states, due to which the off-diagonal kinetic term vanishes. This can be achieved by applying a following transformation 6 , The above matrix has a special symmetry. If we rotate W α 3 and B α by the Weinberg angle tan θ W = g 1 g 2 , the matrix M 2 gauge reduces to a 2 × 2 block diagonal structure with respect to an intermediate state Z α ≡ cos θ W W α 3 − sin θ W B α and X α while the other orthogonal state i.e. A α = sin θ W W α 3 +cos θ W B α having zero mass eigenvalue becomes completely decoupled. This is possible due to the special choice of the transformation matrix we have considered in Eq. (24). Now, once we reduce a 3 × 3 matrix to a 2 × 2 block diagonal form, we already have made our life very simple and next task is to perform another orthogonal transformation between the states Z α and X α to finally get the physical Z and Z µτ bosons. This is mathematically demonstrated below for both mass matrix as well as eigenstates, 6 This transformation matrix is not a unique one. For a general 2 × 2 real matrix, we have four independent elements. However, using c 1 = c 2 = 1 and c 3 = 0, we have only three independent equations to solve for four variables. Here, c 1 , c 2 and c 3 are coefficients of 1 4 B µν B µν , 1 4 X µν X µν and 2 B µν X µν respectively. Thus, one can express three elements in terms of the fourth one and for each real value of that element, we will have a different transformation matrix which eventually cancels the kinetic mixing term. For the particular matrix that we have used here is obtained by setting 2 × 1 element of the transformation matrix equal to zero. Such a special choice easily reproduces all the phenomena of electromagnetism. 7 In the present scenario, after EWSB one can readily determine the mass of the W ± gauge boson which is exactly equal to that of the SM, i.e., M W = 1 2 g 2 v 1 .
where, the masses of two massive neutral gauge bosons (Z and Z µτ ) are respectively given as and the two orthogonal transformation matrices are given by, Finally, the gauge basis and the mass basis of the neutral gauge bosons are related the following orthogonal transformation with where θ W , as mentioned above, is the familiar Weinberg angle and θ µτ is the mixing angle between two neutral gauge bosons Z and Z µτ . These mixing angles can be expressed in terms gauge coupling constants, VEVs and the kinetic mixing parameters as follows, Before we proceed any further, it is worthwhile to mention about the independent parameters. In this model, in addition to the SM parameters, we have fourteen new parameters in the scalar sector (excluding SM-Like Higgs boson mass and VEV v 1 ), three additional Yukawa like coupling constants and one mass term in the extended quark sector 8 and two more couplings in the gauge sector in the form of new gauge coupling g Zµτ and kinetic mixing parameter . These twenty independent parameters are: In terms of these independent parameters the other parameters appearing in the Lagrangian (Eq. (2)) can be obtained using Eqs. (13)(14)(15), Eqs. (18,19), Eqs. (21,23) and Eqs. (30, 34) 9 .
In the present scenario the NP part of the effective Hamiltonian H eff (≡ H SM eff + H NP eff ) that describes the b → s transitions is given by where G F is the Fermi constant, V ij are the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements. Here we neglect other dimension-six operators for example, C 7 can not give significant contributions to the processes, because it corresponds to the dipole operator that is strictly constrained by branching ratio of B → X s γ [118]. Also four-quark operators [119] cannot play any significant role for the violation of LFU, hence they are irrelevant in this work. Moreover, four-fermion contact interactions with scalar currents could be a natural source of LFU violation, although they are highly constrained by existing measurements of the B s → µ + µ − and B s → e + e − branching ratios [120,121]. The NP contribution to the WC C NP, Here, we are not considering Yukawa like coupling constants and bare mass terms in the extended neutrino sector. 9 Additionally, one needs to use minimization conditions of the scalar potential V (H, η, Φ, S).
can be obtained from while the NP contribution to the WC C NP, although we have found that the contribution of C NP, 10 ( ≡ µ) is insignificant 10 and we will focus only on C 9Z(Zµτ ) ( ≡ µ) in rest of the analysis 11 . α em is the fine structure constant. Here are given in the Appendix A. In Fig. 1 we have shown relevant Feynman diagrams responsible for the additional contribution to the b → sµµ transition. It is clearly evident from these Feynman diagrams that the NP contribution to the WC C NP, 9 is provided by the non-standard bottom like fermion field χ and the dark matter candidate ρ 1 with its partners ρ 2 and ρ 3 . Later we provide the dark matter phenomenology of a weakly interacting massive particle (WIMP) type dark matter candidate ρ 1 and related issues by considering the constraints of flavour physics observables that we 10 Due to this reason there is no significant NP contribution to the decay B s → µ + µ − . Therefore, there is no stringent constraint from the branching ratio of this process to our analysis. 11 Therefore, the present scenario can be considered as a typical scenario which can provide the NP contribution to C 9 ( ≡ µ) only. Although, there is a NP contribution to C e 9 but practically it has no significance due to very small mixing between Z and Z µτ . Hence, the coupling between Z µτ and e + e − pair is effectively vanishing in nature. have considered in this article. To ameliorate the tension between the SM prediction and experimental data for For the purpose of notational simplicity, from now and onwards, we use ∆C 9 for the total NP contributions to the WC C 9 for = µ, i.e., C NP,µ 9 = C µ 9Z + C µ 9Zµτ = ∆C 9 . In order to understand the dependence of ∆C 9 on the model parameters we have shown the variation of ∆C 9 in Fig. 2. In this figure there are four panels which represent the variation of ∆C 9 with respect to four important parameters namely M ρ 1 , M χ , M Zµτ and θ D . In Fig. 2a, we have shown the variation of ∆C 9 with mass of ρ 1 for three different values of the product of Yukawa couplings f 2 and f 3 . Here, one can see that the magnitude of ∆C 9 increases with decreasing values of mass of ρ 1 which enters into loop diagrams (see Feynman diagrams shown in Fig. 1). Consequently, the loop functions are enhanced which in turn increase the magnitude of ∆C 9 . Moreover, as the NP contributions to the WC C 9 (Eq. (36)) is proportional to Yukawa  It is clearly seen from Eq. (36), where ∆C 9 is inversely proportional to M 2 Zµτ . On the other hand, in this plot ∆C 9 increases significantly with the gauge coupling g Zµτ for the considered mass range of M Zµτ (0.01 ≤ M Zµτ (GeV) ≤ 0.1). Finally, in Fig. 2d we have demonstrated the variation of ∆C 9 with respect to the dark sector mixing angle θ D for three different choices of M ρ 1 . In this plot, we have varied θ D in range 0 to π/2. The oscillatory behaviour of ∆C 9 with respect to θ D is due the combined effects of two factors. One is the direct involvement of sine and cosine functions within the expressions of ∆C 9 . Another one is the indirect effect due to the change of M ρ 3 with θ D , where the former undergoes a full oscillation between M ρ 2 to M ρ 1 via Eq. (22) when θ D changes from 0 to π. The morphology of ∆C 9 with respect to θ D fits pretty well with a function like −A sin 2 (2θ D ), where the exact value of the normalisation constant A depends on the values of other parameters namely, M ρ 1 , M ρ 2 , g Zµτ , M Zµτ and M χ . Moreover, the oscillatory behaviour of ∆C 9 vanishes if we set M ρ 1 = M ρ 2 . Under this condition, the dependence of θ D disappears from the expression of M ρ 3 and consequently ∆C 9 becomes independent of θ D . Furthermore, in all the four plots of Fig. 2, the grey coloured band represents 2σ range allowed range of fit value of ∆C 9 for explaining R K ( * ) anomalies [27].
The measurement of inclusive radiative B decay process like B → X s γ has also been shown deviation from the corresponding SM prediction. The world average experimental value of the branching ratio of this process is [122] Br for photon energy E γ > 1.6 GeV in the B-meson rest frame. Under the same conditions the corresponding SM prediction with higher order corrections is [123] Br It is quite evident that the theoretical prediction is in good agreement with the experimental value. Hence this small difference can tightly constrain any NP which contributes to this process. Keeping this in mind we have evaluated the NP contributions to this decay process in the present scenario. Consequently, we use the branching ratio of this process as one of the constraints in our analysis.
At quark level B → X s γ decay is indicated by b → sγ transition. The effective Hamiltonian for this transition at the bottom quark mass (µ b = m b ) scale is given by (see ref. [124,125]) At first the WCs (C i ) have been calculated at electroweak scale (µ W ) and using renormalisation group (RG) equations [124][125][126] they are evolved down to µ b = m b scale. The local operators O 1 ....O 6 represent four quark interactions and the explicit form of these operators can be found in [127]. The remaining operators O 7γ (electromagnetic dipole) and O 8G (chromomagnetic dipole) which are the most important for this decay and the expressions for these operators at the leading order are given by The expressions of the WCs at µ b scale is given by and Apart from these other WCs vanish at the electroweak scale µ W . The superscript "0" indicates the leading logarithmic (LO) approximation. The values of a i , h i andh i can be obtained from [126]. The total (SM+NP) contribution at the LO is represented by the functions The functions corresponding to electromagnetic and chromomagnetic dipole operators due to the NP particles (generated form Fig. III B) are given in the following respectively In SM the branching ratio of B → X s γ has been estimated at a very high level of accuracy including higher order QED and QCD corrections. For example in refs. [129,130] one can find the full next-to leading order (NLO) QCD and QED corrections for this process in two different ways. The present precision level of experimental data requires that one should also include next-to-next-to leading order (NNLO) QCD corrections in this analysis. In this regard the first effort to measure NNLO QCD corrections for this process in SM was described in ref. [131]. Finally in a recent article [123] one can find an updated and more complete NNLO QCD corrections to this process. Using the [123] one can calculate the branching ratio of B → X s γ incorporating NNLO QCD corrections in NP scenario. Therefore, in the current article we also follow the same approach 12 (as given in [123]) to measure NP contribution for this process with NNLO QCD corrections Here ∆C 7 and ∆C 8 represent for the NP contributions to WCs for electromagnetic and chromomagnetic dipole operators. In our convention, Using Dirac equation one can define the magnetic moment M of muon in terms of its spin S and gyromagnetic ratio (g µ ) in the following way which is one of the most accurately measured physical quantities. Ideally the value of g µ is equal to "2". In SM one can easily calculate the one loop correction to this quantity and that gives marginal shift from "2". Hence, to measure the deviation of g µ from its tree level value one can define a quantity namely This quantity has been precisely measured by the CERN experiments and later on by the E821 experiment. The current average experimental value is [12] a exp µ = 116592091.0 ± 54 ± 33 × 10 −11 .
On the other hand total theoretical prediction of this quantity considering all kinds of source of contributions in SM is [12] a th µ = 116591823.1 ± 34 ± 26 × 10 −11 .
It is quite evident from the above Eqs. 56 and 57 that both the experimentally measured and theoretically predicted values of a µ are close to each other, however there still exists some disagreement between these two quantities at the 3.5σ significance which is [12], Therefore, this anomaly with respect to the SM expectation requires the interference of BSM theories where one obtains extra contributions from some NP particles. In the present model 13 , apart from the SM contribution, we have two additional one loop diagrams (see Fig. 4) in which the extra neutral gauge boson Z µτ and extra CP-even scalar H 2 are involved. The additional contribution from Fig. 4a is given by [80,135], with R Zµτ ≡ M 2 Zµτ /m 2 µ and Furthermore, the contribution from the extra CP-even scalar H 2 is given by [136,137] δa with However, we have checked that the contribution of CP-even scalar H 2 is insignificant with respect to Z µτ in the allowed parameter space.
In Fig. 5, we have shown the allowed region of M Zµτ and g Zµτ in g Zµτ − M Zµτ plane by red coloured points, which can explain the discrepancy between theoretical prediction (SM) and experimentally measurable value of the anomalous magnetic moment of muon in 2σ range. The corresponding 1σ allowed region is also indicated by green coloured points. We will come back to this parameter space (g Zµτ − M Zµτ plane) with a detailed analysis, which includes constraints like dark matter relic density, direct detection, observables related to rare B-meson decays (R K ( * ) , Br(B → X s γ)) and also bounds from ongoing and future experiments like CCFR, LHC, DUNE, Borexino etc. in the next section (see Fig. 12 and related discussions).

V. DARK MATTER
We are in a stage, where we can discuss dark matter phenomenology. The scalar sector of the present scenario contains two Z 2 -odd scalar representations, one of them is an SU(2) L doublet Φ having a nonzero L µ − L τ charge while the rest is a gauge singlet scalar S. As we have seen earlier in the Section II, the term proportional to λ 8 in the scalar potential (Eq. 10) enforces a mixing between the CP-even component φ 0 of the doublet Φ and the singlet S. Therefore, in the odd sector we have three physical neutral scalars namely, ρ 1 , ρ 2 and ρ 3 , out of which ρ 1 and ρ 2 are two mutually orthogonal linear combinations of S and φ 0 while ρ 3 coincides with the CP-odd component a 0 as the latter does not have any mixing with others. Being Z 2 -odd, the lightest one among the neutral scalars ρ 1 , ρ 2 and ρ 3 is automatically stable and can be an excellent dark matter candidate of the Universe. In this work, we consider ρ 1 as the potential dark matter candidate and depending upon the dark sector mixing angle θ D , ρ 1 will be either "singlet-like" or "doublet-like" or a mixed state. Later in this Section, we will show that although the combined effects of both dark matter relic density bound and flavour physics anomalies (including (g − 2) µ ) considering in this work dictates that the dark matter candidate ρ 1 to be mostly a "single-like" state, its freeze-out process involves extra annihilation channels involving L µ − L τ gauge boson Z µτ , making this scenario significantly different from the case of standard Scalar Singlet dark matter [105][106][107][108].
The viability of the proposed dark matter candidate ρ 1 has been investigated first by computing its relic density 14 Ω DM h 2 . This requires comoving number density Y at the present epoch (T = T 0 , T 0 is the present temperature of the Universe), which is a solution of the Boltzmann equation involving all relevant annihilation and co-annihilation processes in the collision term. The Boltzmann equation in terms of Y is given by [138][139][140], where Y = i Y i with Y i = n i s being the comoving number density of Z 2 -odd particle i having number density n i and s stands for the entropy density of the Universe. Moreover, x = M ρ 1 T is a dimensionless variable and G is the Newton's gravitational constant. The function g [138] depends on degrees of freedom for entropy and energy densities of the Universe. The quantity σv eff has been defined as [139] σv eff = i j where, σ i j v i j is the thermal averaged annihilation cross section between particle i and j having relative velocity v i j . σ i j v i j has the following expression in terms of cross section σ i j , with where, K i is the i th order Modified Bessel function of second kind and s is the Mandelstam variable. Further, Y eq i and n eq i are the equilibrium values of Y i and n i respectively while n = 14 Here, DM represents the short form of dark matter.
i n i is the total number density of all the odd sector particles. This is the most relevant quantity instead of individual n i s, since all heavier particles, which survive annihilation, will eventually decay into the LOP (ρ 1 ). This is the actual reason of expressing the Boltzmann equation in terms of total comoving number density Y instead of individual Y i s. In the above, , represents the mass splitting between LOP and other heavier Z 2 -odd particles.
After implementation of the present model in FeynRules [141] we have solved Boltzmann equation at T = T 0 using micrOMEGAs [142]. Finally, we have obtained Y (T 0 ) which is related to the relic density of LOP through the following relation [140] Relic density Ω DM h 2 of dark matter has been measured precisely by satellite borne experiments like Planck and WMAP and its present acceptable range is 0.1172 ≤ Ω DM h 2 ≤ 0.1226 at 67% confidence level (C.L.) [4]. Apart from this, one has to take into account the latest bound on dark matter nucleon scattering cross section from the "ton-scale" direct detection experiment namely XENON1T [109], which till now provides the most stringent upper bound on dark matter nucleon spin independent scattering cross section (σ SI ) for dark matter mass ranging from 6 GeV to 1 TeV. Since the dark matter candidate of the present scenario is a scalar, it has only spin independent scattering with nucleon and such scattering is possible only though scalar bosons H 1 and H 2 . Feynman diagrams of such elastic scattering ρ 1 + N → ρ 1 + N are shown in Fig. 6. The corresponding expression of σ SI is given by where g H 1 (H 2 )ρ 1 ρ 1 is the coupling between H 1 (H 2 ) and a pair of ρ 1 . Expressions of these couplings are listed in Appendix B. Moreover, f N and M N are nuclear form factor and nucleon mass respectively. For dark matter scattering mediated by scalars f N ∼ 0.3 [108]. We already know that non-observation of any dark matter signal at direct detection experiments impose severe upper bound on σ SI with respect to dark matter mass. From, the above expression of σ SI , it can be seen clearly that such exclusion limit on σ SI in turn puts an upper bound on the involved couplings like g H 1 ρ 1 ρ 1 and g H 2 ρ 1 ρ 1 .
Moreover, SM Higgs to ρ 1 ρ 1 coupling for M H 1 > 2 M ρ 1 case is also constrained from the maximum allowed limit of Higgs invisible decay width . At present, the upper limit on invisible branching fraction of the SM Higgs boson is 0.24 at 95% C.L. [110]. In the present model, the SM like Higgs boson in addition to its "standard decay modes", can also decay into Z µτ Z µτ , ZZ µτ , ρ 1 ρ 1 and ρ 1 ρ 2 final states 15 . Decay widths of such processes are given below, and Expressions of all the coupling involved in the above decay widths are given in Appendix B. According to the latest results from LHC, Γ Inv H 1 ≤ 0.24 Γ SM Higgs , where Γ SM Higgs = 4.13 MeV, is total decay width of the SM Higgs boson [143].
The dark matter candidate (ρ 1 ) of the present scenario is a thermal WIMP, which remains in equilibrium with the thermal bath until its freeze-out though annihilations and co-annihilations into various final states allowed by the symmetries of the model. In this work, we have considered M ρ 1 between 10 GeV to 1 TeV. For low dark matter masses (i.e. M ρ 1 < 100 GeV), ρ 1 predominantly annihilates into a pair of Z µτ . In some cases, depending upon the relevant couplings, bb, cc and ττ final states are also possible. Moreover, co-annihilations among the Z 2 -odd particles in the low mass regime are insignificant as we have considered all heavier Z 2 -odd particles masses larger than 100 GeV throughout this analysis to evade experimental bounds [104]. Alternatively, for the heavier mass range of ρ 1 , there are various possibilities. First of all depending upon the mass splitting between ρ 1 and other Z 2 -odd particles (parametrised by a quantity ∆ i , defined earlier) there can either be annihilation or   ∆ i ≤ 0.2 [139] for any Z 2 -odd particle i (i = ρ 2 , ρ 3 , φ ± ). In this circumstances, various coannihilations among these dark sector particles like In Fig.9, we have plotted spin independent scattering cross section σ SI of ρ 1 with its mass M ρ 1 , varying between 10 GeV to 1 TeV. In this plot all red coloured points satisfy relic density constraint i.e., 0.1172 ≤ Ω DM h 2 ≤ 0.1226 and bound related to Higgs invisible decay modes as well. The blue dashed-dot line represents the latest bound on σ SI from XENON1T experiment. All the parameter space below the blue dashed-dot line are still allowed and can be probed in near future by "multi-ton-scale" direct detection experiments like XENONnT. Therefore, if we consider direct constrains like relic density, direct detection and Higgs invisible decay only, there are still enough parameter space left (although few portion mostly in the low mass dark matter regime have already been ruled-out) for the entire considered mass range of ρ 1 . However, the situation does not remain same when one tries to explain flavour physics anomalies and (g − 2) µ anomaly within this framework. The allowed parameter space in σ SI − M ρ 1 plane gets severely restricted when we impose bound on the NP contribution to the WC C 9 (i.e. −1.26 ≤ ∆C 9 ≤ −0.63 in 2σ range [27]) to explain R K ( * ) anomalies. This has been indicated by green coloured points in the above plot where one can notice that the low dark matter mass regime (i.e. M ρ 1 < ∼ 100 GeV) is the most favourable to address R K ( * ) anomalies. This can be understood from the behaviour of ∆C 9 (Eq. (36)) with respect to the mass of ρ 1 as illustrated in Fig. 2a, where the magnitude of ∆C 9 sharply decreases with the increase of M ρ 1 . Furthermore, in this framework, we have also tried to explain both Br(B → X s γ) and (g − 2) µ anomaly, the two long-standing anomalies of the SM from their experimental counterparts. These are indicated by cyan and yellow coloured points respectively in σ SI − M ρ 1 plane. For the branching ratio of B → X s γ, we have used 3σ (2.84 ≤ Br(B → X s γ) × 10 4 ≤ 3.80 [122]) while the 2σ band i.e. 115.44 ≤ ∆a µ × 10 9 ≤ 420.56 [12] for (g − 2) µ has been taken into account 16 . We 16 Here we would like to mention that, another constraint e.g.,  have checked that in the low dark matter mass region (M ρ 1 ≤ 100 GeV), ρ 1 predominantly annihilates into the Z µτ pair. This actually makes dark matter physics strongly correlated with the physics of rare B-decays and anomalous magnetic moment of µ, where the role of new gauge boson Z µτ is extremely crucial. Moreover, it also helps us to evade the strong bound coming from the experiments of direct detection [109], indirect detection [145] and also from the collider on Higgs invisible branching [110] for the low mass scalar dark matter [146][147][148], where bb final state is the principal annihilation channel. Therefore, in spite of being a gauge singlet Z 2 -odd scalar field, the mixing with another Z 2 -odd field (part of an SU(2) L doublet) having nonzero L µ − L τ charge, makes the entire dynamics of our dark matter candidate ρ 1 strikingly different from the standard Scalar Singlet dark matter scenario [105][106][107][108]. Finally, for the completeness we would like to mention here that the yellow coloured points in σ SI − M ρ 1 plane are those which satisfy all the experimental results we have considered in this work.
In the left panel of Fig. 10, we have shown ranges of M ρ 1 and M ρ 2 allowed by various experimental results. The allowed region in M ρ 2 − M ρ 1 plane from both relic density as well as direct detection bounds are indicated by the green coloured points. Similar to the previous plot in diagrams and these are negligibly small. The reason is that, apart from the dark matter particle, all nonstandard particles which generate box diagrams are sufficiently massive (especially the non-standard fermion χ, whose mass that we have taken ≥ 1 TeV throughout the analysis). At this point it is relevant to mention that, from the recent 13 TeV LHC data [144], a down-type quark (B) with charge (-1/3) is excluded for masses below 1.22 TeV for the decay channels B → Zb/W t/SM Higgs b. However, this bound is not applicable in our case, as in our model the field χ is odd under Z 2 symmetry, therefore such decays are restricted by the Z 2 symmetry. Although, for the sake of conservative approach we use M χ ≥ 1 TeV in our analysis. Hence, the loop functions are substantially suppressed. Consequently, the NP contribution to B 0 s −B 0 s mixing would not put any stringent constraint in our scenario. Fig. 9, here also when we have imposed various flavour physics constraints, the allowed parameter space shrinks to a smaller region concentrated mainly in the low mass regime of ρ 1 . The parameter space which reproduces ∆C 9 in 2σ range for explaining R K ( * ) anomalies has been shown by the blue coloured points. On the other hand, the red coloured points are indicating those values of M ρ 1 and M ρ 2 which in addition to above mentioned experimental results also satisfy Br(B → X s γ) in 3σ range. Moreover, as we have already known that the dark matter candidate ρ 1 is an admixture of a real scalar singlet S and a CP-even neutral component (φ 0 ) of a doublet Φ. While both S and Φ are Z 2 -odd but only Φ has nonzero L µ −L τ charge. Therefore, the interaction of ρ 1 with L µ − L τ gauge boson Z µτ (e.g. annihilation of ρ 1 into a pair of Z µτ ) is governed by the mixing angle θ D . Larger the mixing angle, larger is the annihilation rate into Z µτ Z µτ final state, making ρ 1 less abundant at the present epoch. Therefore, the relic density bound puts an upper limit on the maximum allowed value of θ D , which is more stringent in the low dark matter mass region where Z µτ Z µτ is the principal annihilation mode. This feature is clearly visible in the right panel of Fig. 10, where we have shown the allowed range of θ D with respect to M ρ 1 . However, in the high mass regime (M ρ 1 ≥ 500 GeV), large values of θ D > ∼ 0.3 rad are still allowed because for such large θ D , ρ 1 is mostly an SU(2) L doublet like state (similar to the Inert Doublet dark matter in high mass range [148][149][150]) which attains the present abundance of dark matter through co-annihilations with other Z 2 -odd fields into various bosonic final states (both vector and scalar). Moreover, we have also seen from the Fig. 9 that the magnitude of ∆C 9 (Eq. (36)) decreases with the increasing values of masses of the particles ρ 1 , ρ 2 and ρ 3 involving within b → s transition loops. Now, although the masses of ρ 1 and ρ 2 are indeed free parameters of the present model, the mass of the remaining scalar ρ 3 becomes fixed for a particular choice of M ρ 1 M ρ 2 and θ D via Eq. (22). Here, M ρ 3 actually oscillates between M ρ 2 and M ρ 1 as we vary θ D from 0 to π/2. As we are working in the limit M ρ 1 < M ρ 2 (since ρ 1 is our dark matter candidate), large θ D ensures low mass for ρ 3 (using Eq. (22)) and hence enlarge loop contribution to ∆C 9 . Thus, R K ( * ) anomalies prefer larger values of θ D , which is a contrasting situation compared to the low mass regime of ρ 1 , where relic density bound favours relatively smaller values of mixing angle to suppress large annihilation into Z µτ Z µτ . As a result, both dark matter relic density bound and R K ( * ) anomalies are simultaneously addressable for 0.01 < θ D (rad) < 0.3, when M ρ 1 is mostly concentrated below 100 GeV range. This has been demonstrated by the blue coloured points in θ D − M ρ 1 plane. Similar to the left panel, here also red colour points represent the portion in the parameter space which has been satisfied by the constraint of Br(B → X s γ) as well.
Since, the allowed values of θ D which satisfy all the experimental results considered in this work fall in the range 0.01 < θ D (rad) < 0.3 (red coloured points in the right panel of Fig. 10), this makes M ρ 2 and M ρ 3 almost degenerate and this has been demonstrated in Fig. 11, where the colour bar is indicating corresponding values of mass of the dark matter candidate ρ 1 .
Finally, in Fig. 12 we have shown our results in g µτ −M Zµτ plane, which is at the present moment extremely constrained by various experimental results. In this figure (Fig. 12), the red coloured points represent those values of g µτ and M Zµτ which explain (g − 2) µ in 2σ range. Here, in the g µτ − M µτ plane, most strongest constraint till now comes from neutrino trident production. Neutrino trident production is a process of producing µ + µ − pair via neutrino scattering in the Coulomb field of a target nucleus (N ), i.e. ν µ (ν µ ) + N → ν µ (ν µ ) + µ + µ − + N . In the SM, this process is possible via W ± and Z bosons only. Moreover, if there exists any new neutral gauge boson (similar to Z µτ in the present work) which couples to both muons and muon-neutrinos then that gauge boson can also contribute significantly to the trident production cross section. However, all the experimental collaborations namely, CCFR [111], CHARM-II [112] and NuTeV [151] have measured neutrino trident events and their measured cross sections are in good agreement with that of the SM prediction i.e. σ CCFR σ SM = 0.82 ± 0.28, σ CHARM−II σ SM = 1.58 ± 0.57 and σ NuTeV σ SM = 0.72 +1.73 −0.72 respectively. These results therefore put strong constraint in the masscoupling plane of the new gauge boson. In Fig. 12, the crossed region above the black dashed line represents 95% C.L. upper bound [152] on g µτ as a function of M Zµτ using neutrino trident cross section measured by the CCFR collaboration 17 . Consequently, all the crossed regions above black dashed line are excluded by neutrino trident production. Besides, there is a further constraint from the measurement of the SM Z boson decay to 4µ final state at the LHC. This 17 Furthermore, it is clearly evident from the Fig. 12, that, due to the consideration of CCFR experimental data we naturally incorporate the constraint of the branching ratio of τ → µν τνµ . The reason is that the parameter space (in g µτ − M Zµτ plane) which describes all the concerned observables simultaneously does not overlap with the portion that has already been ruled out from the branching ratio of τ → µν τνµ [96]. Moreover, we have explicitly checked that the NP contribution for the decay τ → µν τνµ due to Z µτ is practically vanishing in nature in the allowed parameter space of the present scenario.

D U N E
Ω DM h 2 + XENON1T + ΔC 9 + B→ X s γ Ω DM h 2 + XENON1T + ΔC 9 (g-2) μ in 2σ Ω DM h 2 + XENON1T has also been indicated by the grey region at the topmost corner of right side of this plot. Cyan coloured points represent those values of g Zµτ and M Zµτ which satisfy bounds related to dark matter physics namely, relic density, direct detection and Higgs invisible branching ratio. On top of the existing dark matter constraints, the effects of flavour physics observables like R K ( * ) anomalies (2σ bound on ∆C 9 ) and R K ( * ) + Br(B → X s γ) on the mass as well as the coupling of Z µτ have been shown by green and yellow coloured points respectively. Therefore, from this plot it can be easily seen that although maximum portions of g µτ − M Zµτ plane have already been excluded by the results of CCFR collaboration, there is still a small but interesting region left in this parameter space which is 0.01 ≤ M Zµτ (GeV) ≤ 0.1 and 3×10 −4 ≤ g µτ ≤ 10 −3 . This region of the parameter space of the present model can address dark matter, (g − 2) µ anomaly, R K ( * ) anomalies and Br(B → X s γ) simultaneously and more exciting thing is that this parameter space can be probed within a few years by the DUNE experiment [153] measuring neutrino trident events (shown by black dashed line) [154]. This will surely be the test of our model, at least the benchmark points (if not the full model) in the low mass dark matter region which are compatible to both dark matter and flavour physics issues. For completeness in Table III

VI. NEUTRINO MASSES AND MIXINGS
In this section, we will discuss briefly about neutrino masses and mixings. It has now been firmly established from the phenomena of neutrino oscillations that there exist two tiny mass square differences between three neutrino mass eigenstates i.e. ∆m 2 21 = 7.39 +0.21 −0.20 × 10 −5 eV 2 18 and ∆m 2 31 = 2.525 +0.033 −0.032 (−2.512 +0.034 −0.032 ) × 10 −3 eV 2 for the normal(inverted) hierarchy [155] in 3σ range. This also indicates that to explain solar, atmospheric and rector neutrino anomalies though three flavour neutrino oscillation we need at least two neutrino mass eigenstates having nonzero masses corresponding to mass squared differences as mentioned above. Moreover, there are also precise measurements of three intergenerational mixing angles namely the atmospheric mixing angle (40.3 19 , the solar mixing angle (31.61

36.27
• ) and the reactor mixing angle (8.22 [155]. The latter one is the most recent entry in that list. In the present model, although we do not need any extra fermionic degrees of freedom to cancel L µ − L τ anomaly which actually cancels between µ and τ generations of charged leptons and corresponding neutrinos, one can still introduce three right handed neutrinos N Ri (i = e , µ, , τ ) in an anomaly free manner, in the model, to address neutrino masses and mixings via Type I seesaw mechanism. The Lagrangian for right handed neutrinos are given in Eq. (4). The light neutrino mass matrix m ν after spontaneous breaking of both SU(2) L × U(1) Y and U(1) Lµ−Lτ symmetries has the following structure while the mass matrix for the heavy neutrinos coincides with M R . In the above, p = y eµ y eτ v 2 2 − M ee M µτ e iξ . Majorana mass matrix M R and Dirac mass matrix M D are given by, In the present case, due to L µ − L τ flavour symmetry, the Dirac mass matrix is exactly diagonal while before U(1) Lµ−Lτ symmetry breaking only three elements (only two are independent) are 19 Where numbers without(within) brackets are for the normal(inverted) hierarchical scenario.
there in the Majorana mass matrix M R . Only after symmetry breaking we get additional elements proportional to v 2 . Therefore, L µ − L τ symmetry breaking plays a crucial role here to get desire structure of m ν matrix. Also, looking at both M D and M R matrices, one can easily notice that there can only be one complex element. Phases of other elements can be absorbed by redefining both SM leptons and right handed neutrinos. Now, one can calculate mass eigenvalues and mixing angles of light neutrinos by diagonalising this m ν matrix, which is a complex symmetric matrix, indicating the Majorana nature of the light neutrinos. If we consider, v 2 ∼ 10 2 GeV (in the right ballpark to produce desired contribution to (g − 2) µ ), 0.1 < ∼ y eµ , y eτ < ∼ 1.0 and 10 GeV < ∼ M ee , M µτ < ∼ 1 TeV (100 GeV to TeV scale right handed neutrinos) then we need Dirac couplings 10 −7 < ∼ y e , y µ , y τ < ∼ 10 −5 to reproduce neutrino oscillation parameters. Detail analysis of mass matrix diagonalisation and comparison with latest 3σ range of oscillation parameters have been done in Ref. [84]. Moreover, we would like to mention here that although only two right handed neutrinos (N µ R and N τ R ) are sufficient to make the present model anomaly free, such scenario is unable to reproduce all neutrino oscillation parameters due to special flavour structure of the Dirac mass matrix.

VII. CONSTRAINT FROM DI-LEPTON RESONANCE SEARCH AT 13 TEV LHC
Depending on the mass ranges, the non-standard Z boson (which we designate as Z µτ in this article) confronts constraints from collider searches. For example, if the mass of Z µτ is less then SM Z boson then there exists some viable parameter region for the favorable kind among the various NP models that exist in the literature. Furthermore, as the Z µτ has no direct coupling with electron 20 , hence LEP searches cannot provide direct constraint on the light Z µτ . On the other hand, the Tevatron [156,157] and LHC [115,158,159] searches for Z µτ to di-lepton final state only apply, however in this case M Zµτ > 100 GeV. Moreover, only relevant limit to the light Z µτ case obtained from the LHC searches at pp → Z → 4µ [152]. At this point we remark in passing that, in our present article even though in the low mass limit of Z µτ we have obtained certain region of parameter space (depicted in Fig. 12) which has been satisfied by some flavour physics data, dark matter constraints and (g − 2) µ anomaly, however, cross section for a process like pp → Z µτ → + − in that region of parameter space is extremely tiny at the 13 TeV LHC.
On the other hand in the high mass region of Z µτ , the LHC searches put the tightest bound on its mass (3 ∼ 5 TeV [115,158,159]) in the di-muon final states. Thus, in the present article we use the exclusion data obtained by ATLAS collaboration [115] for a di-lepton resonance search at the LHC experiment to constraint parameter space of the present scenario. In order to embed this limit in the present scenario, we first implement the model using FeynRules [141]. Then we generate the cross section for the process pp → Z µτ → + − using Madgraph5 [160] with the default parton distribution functions NNPDF3.0 [161] at 13 TeV LHC 21 . Here (≡ e, µ), 20 Only possible via Z-Z µτ mixing. Therefore, the interaction strength is insignificant. 21 Production of Z µτ at the LHC in the present model is possible due to the couplings of q iqi Z µτ which are however, significant contribution has been generated from µ + µ − final state. Finally, for a specific combination of coupling g Zµτ and Z-Z µτ mixing angle θ µτ we compare the theoretical prediction of cross section for any particular value of mass (confined within the range [0. 5,5] TeV) of Z µτ with the corresponding experimental data given by ATLAS collaboration [115].
In Fig. 13 we show the exclusion curves at 95% C.L. in the M Zµτ − g Zµτ plane for four different values of Z-Z µτ mixing angle θ µτ using the ATLAS data [115] for non-observation of a resonant + − signal at the LHC running at 13 TeV with integrated luminosity 139 fb −1 . In this case the region above a particular curve has been ruled out at 95% C.L. from the non-observation of a resonant + − signal in the 13 TeV run of LHC by ATLAS data [115]. If we focus on a particular curve fixed by a particular value of mixing angle θ µτ then we observe that for the lower values of mass the coupling g Zµτ rapidly falls with the increasing values of mass M Zµτ . This phenomena can be explained in the following way. In the lower mass range if we vary the mass then the cross section does not fall rapidly as desired by the ATLAS data. Hence, to acquire the proper cross section for a particular mass one should decrease the value of the coupling g Zµτ . Once the lower mass range is over then with the increasing values of mass the generated via Z-Z µτ mixing.
curve exactly replicates the exclusion plot as given in [115]. At this point, we would like to mention another notable feature of the exclusion curves (which is true for all over the mass range) that for a fixed value of mass if the mixing angle increases then to satisfy ATLAS data [115] one requires decreasing values of coupling constant g Zµτ . Furthermore, it is clearly evident form the Fig. 13 that as the mixing angle θ µτ increases large amount of area in the M Zµτ − g Zµτ plane has been ruled out by the ATLAS data. Both of the features can be explained, if we analyse the structure of the coupling 22 between Z µτ and + − . If we decompose the coupling then we can find that there is one vectorial part and other is axial vectorial in nature. The latter one has no significant role in the concerned process but totally controlled by the vectorial part. We have also checked that, one can control the coupling (which in turn the vectorial part of the coupling) that satisfy the exclusion data with lower values of mixing angle θ µτ . However, as the mixing increases then one looses the control over the coupling, i.e., there is no variation of coupling with larger mixing angle. Therefore, with larger mixing angle one can not vary the cross section properly, hence one can not have the required cross section for a particular mass. For example, if the mixing is set at 4.5 × 10 −4 rad, then one can not go beyond 1500 GeV mass of M Zµτ . Since in this situation after 1500 GeV mass we can not have the desired cross section by changing the value of g Zµτ . Therefore, in order to translate the exclusion limit obtained by ATLAS data [115] for non-observation of di-lepton resonance search at the LHC experiment in our model we have restricted ourselves within the relatively smaller values of mixings angle θ µτ .

VIII. CONCLUSIONS
In order to simultaneously resolving R K ( * ) anomalies and dark matter enigma, we have proposed a unified scenario by introducing an extra local U(1) Lµ−Lτ gauge symmetry to the Standard Model. This U(1) Lµ−Lτ gauge symmetry provides a neutral non-standard gauge boson Z µτ which has versatile effects on different phenomenological aspects that have been considered in this article. For the purpose of breaking of the U(1) Lµ−Lτ symmetry spontaneously a complex scalar field η has been invoked to the scalar sector in addition to the usual Standard Model Higgs doublet H. Three singlet right handed neutrinos have also been introduced in order to explain the observed oscillation data by incorporating neutrino masses and mixings via Type-I seesaw mechanism. Furthermore, for the proper establishment of correlation between R K ( * ) anomalies and dark matter puzzle, a bottom quark like coloured fermion field χ has been included in this scenario. This non-standard fermion field χ is transformed vectorially under the U(1) Lµ−Lτ symmetry and further it is odd under the Z 2 parity. Apart from these, an SU(2) L scalar doublet Φ with nonzero U(1) Lµ−Lτ charge and a real scalar singlet S have also been incorporated in the present scenario. Both of these non-standard scalar fields are odd under Z 2 symmetry. The mixing (which is parametrised by a mixing angle θ D ) between these two Z 2 -odd scalar fields gives a potential dark matter candidate ρ 1 and also two heavier Z 2 -odd physical particles ρ 2 and ρ 3 . All of these three scalar fields provide significant contributions not only in dark matter phenomenology but also in rare B-meson decay processes.
Existence of lepton flavour universality violation in neutral current sector has been measured by R K ( * ) in which b → s + − ( ≡ e, µ) transition is involved. This type of flavour changing neutral current is highly suppressed in the Standard Model and therefore, even for a small deviation between the experimental data and the Standard Model could play significant role for finding of new physics effects. In this work, the introduced new physics particles have played crucial role in the concerned b → s transition processes which are in general loop induced 23 . Particularly, the dark matter particle ρ 1 with two heavier Z 2 -odd neutral scalar fields ρ 2 , ρ 3 and the nonstandard fermion χ generate extra loop contributions. Furthermore, the extra non-standard gauge boson Z µτ behaves as a propagator (in addition to the SM Z boson) for the process b → s + − . Now, due to the very basic structure of our model, the process b → sµ + µ − is more favourable with respect to b → se + e − , consequently one obtains the significant non-standard contribution to the Wilson coefficients C NP 9 for "µ" but not for "e". Therefore, in our work, we have easily satisfied the current fit result for C NP,µ 9 ∈ [−1.26, −0.63] in 2σ interval to explain the R K ( * ) anomalies and thereby we have constrained the parameter space of the proposed scenario. On top of that, we have also calculated another rare decay process B → X s γ which has also been a class of processes that characterised by b → s transition. We have estimated the branching ratio for B → X s γ process, and have used the corresponding experimental data within 3σ interval as one of the constraints in our analysis. Moreover, we have calculated the contribution of non-standard gauge boson Z µτ to the (g − 2) µ and considering the recent experimental data with some error bars (1σ and 2σ) we have further constrained the parameter space allowed by dark matter and flavour physics observables.
In the present scenario, we have extensively studied the dark matter phenomenology by choosing ρ 1 as a WIMP type dark matter candidate. This ρ 1 is an admixture of a real scalar singlet S and the CP-even neutral component (φ 0 ) of the doublet Φ. In our work, first we have calculated dark matter relic abundance by considering all possible annihilation and co-annihilation channels for a wide range (10 GeV ≤ 1 TeV) of the mass of ρ 1 . Thereafter, we have imposed necessary constraints like Planck limit on relic density (0.1172 ≤ Ω DM h 2 ≤ 0.1226), latest direct detection bounds on σ SI from XENON1T and also the bound on Higgs invisible branching ratio from LHC to find the allowed parameter space. We have found that in the case of low mass region (M ρ 1 < 100 GeV), our dark matter candidate ρ 1 predominantly annihilates into Z µτ pair while co-annihilations among other Z 2 -odd particles are insignificant as we have considered all heavier Z 2 -odd particles masses larger than 100 GeV throughout this analysis to respect the experimental bounds form LEP collider. Due to this primary annihilation channel (ρ 1 ρ 1 → Z µτ Z µτ ), in spite of being a gauge singlet Z 2 -odd scalar field, the mixing with another Z 2 -odd field (part of an SU(2) L doublet) having nonzero L µ − L τ charge, makes the entire dynamics of our dark matter candidate ρ 1 remarkably different from the standard Scalar Singlet dark 23 Apart from leptoquark scenarios where b → s transition is possible at tree level. matter scenario where bb final state is in general the principal annihilation channel and low mass region has already been ruled out by direct detection, indirect detection and also by the upper limit on Higgs invisible decay branching ratio. On the other hand for the higher values of M ρ 1 , depending upon the mass splitting between ρ 1 and other Z 2 -odd particles several annihilation or co-annihilation channels may appear and have contributed significantly to the relic density. Since, one of our prime motivations of this article is to correlate dark matter puzzle with some specific flavour physics anomalies associated with FCNC processes, therefore, we have used experimental data of some flavour physics observables (e.g., R K ( * ) anomalies and Br(B → X s γ)) as further constraints on the parameter space which is already allowed by experiments related to dark matter physics. As a consequence, both the effects of R K ( * ) anomalies and dark matter phenomenology allow only a very restrictive values of dark sector mixing angle θ D which remains confined within a certain range (0.01< θ D (rad) < 0.3) when M ρ 1 ≤ 100 GeV. This is a unique feature of our proposed model. Additionally, we have used some other constraints which have been relevant to our present scenario. For example, we have imposed constraint from neutrino trident production and for that purpose we have used the CCFR experimental data which is currently the most stringent one for the neutrino trident production process. Furthermore, we have imposed constraint from the measurement of the Standard Model Z boson decay to 4µ final state at the LHC. As a consequence there is a substantial amount of reduction in the parameters space due to the inclusion of such constraints. However, there still exists a few portion of the parameter space of the present model which can address dark matter, R K ( * ) anomalies, (g − 2) µ and Br(B → X s γ) simultaneously. Most importantly our predicted parameter space and hence our model can be tested within a few years by neutrino trident processes at DUNE. Therefore, in view of the above discussion we can readily conclude that our proposed scenario can reasonably connect the dark matter puzzle with some of the flavour physics anomalies. Besides, within the scope of our proposed model, we have also briefly discussed the origin of neutrino masses and mixing angles via Type-I seesaw mechanism, which is a common feature of most of the L µ − L τ models.
Finally, for the purpose of constraining the parameter space of the present scenario from the LHC, we have used the latest ATLAS data of non-observation of a resonant + − signal at the LHC running at 13 TeV with an integrated luminosity 139 fb −1 . For this purpose we have estimated the cross section for the process pp → Z µτ → + − at the 13 TeV LHC for the mass range M Zµτ ∈ [0. 5,5] TeV in the present scenario. By comparing the theoretical predictions of the cross section with corresponding ATLAS data of cross section for non-observation of a resonant + − signal at the 13 TeV LHC one yields some specific combination of coupling g Zµτ and Z-Z µτ mixing angle θ µτ . Consequently, with those combinations we have excluded some portion of the parameter space of the present scenario at 95% C.L. From our analysis it has been observed that, for a larger values of mixing angle one can exclude larger region of parameter space in the M Zµτ − g Zµτ plane. For example if the mixing angle is 4.5 × 10 −5 rad then one can maximally exclude the region of parameter space in the M Zµτ − g Zµτ plane.