FIMP and Muon ($g-2$) in a U$(1)_{L_{\mu}-L_{\tau}}$ Model

The tightening of the constraints on the standard thermal WIMP scenario has forced physicists to propose alternative dark matter (DM) models. One of the most popular alternate explanations of the origin of DM is the non-thermal production of DM via freeze-in. In this scenario the DM never attains thermal equilibrium with the thermal soup because of its feeble coupling strength ($\sim 10^{-12}$) with the other particles in the thermal bath and is generally called the Feebly Interacting Massive Particle (FIMP). In this work, we present a gauged U(1)$_{L_{\mu}-L_{\tau}}$ extension of the Standard Model (SM) which has a scalar FIMP DM candidate and can consistently explain the DM relic density bound. In addition, the spontaneous breaking of the U(1)$_{L_{\mu}-L_{\tau}}$ gauge symmetry gives an extra massive neutral gauge boson $Z_{\mu\tau}$ which can explain the muon ($g-2$) data through its additional one-loop contribution to the process. Lastly, presence of three right-handed neutrinos enable the model to successfully explain the small neutrino masses via the Type-I seesaw mechanism. The presence of the spontaneously broken U(1)$_{L_{\mu}-L_{\tau}}$ gives a particular structure to the light neutrino mass matrix which can explain the peculiar mixing pattern of the light neutrinos.


I. INTRODUCTION
The presence of Dark Matter (DM) is a well established truth and demands the extension of the Standard Model (SM). Starting from the time of Zwicky (who examined the coma-cluster in 1933 and gave the prediction of excess matter) till now, many collaborations have reported the existence of DM. The most prominent examples are, flatness of galaxy rotation curve [1], gravitational lensing [2] and the Bullet cluster observed by NASA's Chandra satellite [3]. In the last one decade, satellite borne experiments like WMAP (Wilkinson Massive Anisotropy Probe) [4] and Planck [5] not only confirmed the existence of DM but also measured the amount of DM present in the Universe with an unprecedented accuracy. However, the particle content of dark matter as well as its interaction strength with the visible world are two questions which remain unresolved. Depending on their production mechanism in the early Universe, the dark matter particles can be classified into two groups. The first group is the thermal dark matter. These class of particles were produced thermally and maintained their thermal as well as chemical equilibrium with the thermal bath. They eventually decoupled from the plasma when their interaction rates became less than the expansion rate of the Universe. Being relic particles, their comoving number density froze to a particular value. This is known as the usual freezeout mechanism [6,7]. The dark matter candidates which were produced through the freeze-out scenario are usually called the Weakly Interacting Massive Particle or WIMP [6,7]. These types of particles have weak scale interaction cross section and their signatures are expected to be seen at various ongoing direct detection experiments like LUX [8,9] and XENON-1T [10]. However, the direct detection experiments are yet to find any real signal of dark matter particles. This forces us to look for other types of dark matter particles which can be an alternative to the WIMP scenario. One such class of dark matter candidates is the Feebly Interacting Massive Particle or FIMP [11][12][13][14][15][16][17][18]. The interaction rates of these particles are so feeble that they never attained thermal equilibrium with the thermal bath. As a result, their initial number densities were negligible compared to others which were in thermal equilibrium. As the Universe cooled down, they were produced predominantly from the decays of heavy particles. If the decaying mother particles are in thermal equilibrium at the early stage of the Universe, the production of dark matter particles from a particular decay mode is expected to be maximum around a temperature equal to the mass of the corresponding decaying particle. In principle, FIMP can also be produced from the annihilation of SM particles and other heavy particles (beyond SM) as well. The production mechanism of FIMP is known as freeze-in [11] which is opposite to the freeze-out scenario. Unlike the freeze-out mechanism, in case of freeze-in the comoving number density and hence the relic density of a FIMP is directly proportional to its couplings with other particles. Since the interaction strength of a FIMP is extremely weak, thus one can evade all the constraints arising from the direct detection experiments.
Besides the dark sector, the other puzzles which require the extension of the SM are the observation of neutrino mass [19], discrepancy in the anomalous magnetic moment of muon (g−2) [20] from its SM prediction and also the baryon asymmetry of the Universe [21]. The existence of extremely small neutrino masses were first confirmed by the Super-Kamiokande collaboration [22] from the observation of neutrino flavor oscillations in atmospheric neutrino data. Thereafter, many outstanding neutrino experiments like SNO [19], KamLand [23], Daya Bay [24], RENO [25], Double Chooz [26], T2K [27,28] and NOνA [29,30] have precisely measured the two mass squared differences and three mixing angles between different neutrino flavours.
In the present work, we try to explain the existence of dark matter, tiny nature of neutrino masses, their intergenerational mixing angles and muon (g − 2) anomaly simultaneously within the framework of our proposed model. Since we already know that the SM is unable to address these issues, we have to think of a beyond Standard Model (BSM) scenario. There are many well motivated extensions of the SM like SUSY, two Higgs doublet model, extension of the SM gauge group by extra U(1) and many more. In this work we have extended the SM gauge group by a local U(1) Lµ−Lτ symmetry [31][32][33]. Therefore, the complete gauge group in this model is SU(2) L × U(1) Y × U(1) Lµ−Lτ . Since we are increasing the SM gauge group by a local U(1) Lµ−Lτ symmetry, we must check the cancellation of axial vector anomaly [34,35] and mixed gravitational-gauge anomaly [36,37]. The advantage of the U(1) Lµ−Lτ extension is that these anomalies cancel automatically between the second and third generations [38][39][40] of SM fermions. We also extend the particle content of the SM and include three right-handed neutrinos and two SM gauge singlet scalars. The new particles are given appropriate U(1) Lµ−Lτ charge. One of the scalars picks up a Vacuum Expectation Value (VEV) breaking U(1) Lµ−Lτ symmetry spontaneously, while the other does not take any VEV. The U(1) Lµ−Lτ charge n µτ of the other scalar φ DM is chosen in such a way that it remains stable even after U(1) Lµ−Lτ breaking and becomes a dark matter candidate. The new gauge boson Z µτ after spontaneous breaking of U(1) Lµ−Lτ becomes massive and gives additional contribution to the muon (g − 2), helping reconcile the observed data with the theoretical prediction in this model. The three righthanded neutrinos carry Majorana masses and give rise to light neutrino mass matrix via the Type-I seesaw mechanism. In our earlier work in Ref. [32], we showed that our model can explain the nonzero neutrino masses and their intergenerational mixing angles. Also in that work, we considered the scalar field φ DM as a WIMP type dark matter candidate and checked its viability in various direct detection experiments. We found that all the existing constraints are satisfied around the resonance regions, where the mediator mass is nearly equal to twice of dark matter mass.
In this work, we show that φ DM could serve as a FIMP type dark matter candidate in this U(1) Lµ−Lτ extension of the SM. All particles of our model including the additional scalar h 2 , Z µτ and the three right-handed neutrinos are in thermal equilibrium in the early Universe except φ DM . In order to make this possible, φ DM must be feebly interacting with other particles. We choose its couplings with the visible sector to be extremely small (see Section II and IV for detailed discussions) such that it remains out of equilibrium throughout its evolution in the early Universe 1 . We compute the relic density of φ DM by solving the Boltzmann equation where we consider all the possible production modes of φ DM from the decays as well as the annihilations of SM and BSM particles. We find that in our model, φ DM is produced not only from the decays of h 1 and h 2 , but also from the pair annihilation of N i (i = 2, 3) mediated by the extra neutral gauge boson Z µτ . Therefore, the dark matter phenomenology is intricately intertwined with the phenomenology of neutrino masses and muon (g − 2).
Rest of the paper is organised in the following way. In Section II we describe the model in detail. In Section III we discuss muon (g − 2) and neutrino masses and mixings. The Boltzmann equation for FIMP considering all possible production channels is given in Section IV. Our main results are presented in Section V. Finally, we summarise the present work in Section VI. All the relevant decay widths and annihilation cross sections are given in Appendices A and B.

II. MODEL
We have considered the U(1) Lµ−Lτ extension of the SM, where L µ and L τ are the muon and tau lepton numbers, respectively. Hence, the complete gauge group of the minimally extended SM is SU(2) L × U(1) Y × U(1) Lµ−Lτ . The particle content of our model is given in Tables I and  II.
where Γ is the dark matter production rate and H is the Hubble parameter [41]  The Lagrangian for this model is as follows, where L SM represents the Lagrangian of SM and L N is the Lagrangian for the RH neutrinos which contains its kinetic energy terms, mass terms and Yukawa terms with the SM leptons doublets. The expression of L N is, whereφ h = i σ 2 φ * h and the quantities M µτ , M ee have dimension of mass while the Yukawa couplings h eµ , h eτ and y α s are dimensionless constants. In Eq. (1), L DM part is the dark matter Lagrangian which contains the kinetic term of the dark matter candidate φ DM and the interaction terms of φ DM with the SM-like Higgs (φ h ) and the extra Higgs (φ H ). The expression of the L DM is as follows, The fourth term in Eq. (1) represents the kinetic term of the extra Higgs singlet φ H . The potential V (φ h , φ H ) contains the quadratic and quartic interactions of φ H and the interaction term between φ h and φ H . The potential V (φ h , φ H ) has the following form, The last term in the Eq. (1) is the kinetic energy term of the extra gauge boson Z µτ in terms of the field strength tensor F αβ µτ which has the following form, F αβ µτ = ∂ α Z β µτ − ∂ β Z α µτ . We can write down the generic form of the covariant derivative which are appearing in the Eqs. (1)-(3) is, where X is any SM singlet field which has U(1) Lµ−Lτ charge Q µτ (X) (see Table II) and g µτ is the gauge coupling of the U(1) Lµ−Lτ gauge group. The U(1) Lµ−Lτ symmetry breaks spontaneously once φ H picks up a VEV and consequently the extra gauge field becomes massive with mass term M Zµτ = g µτ v µτ . In the unitary gauge, φ h and φ H take the following form after spontaneous breaking of the SU(2) L × U(1) Y × U(1) Lµ−Lτ gauge symmetry, where v and v µτ are the VEVs of φ h and φ H respectively. The presence of the mutual interaction between φ h and φ H in Eq. (4) induces the mixing between H and H µτ . Hence, the scalar mixing matrix has the following form, We see from the expression of M 2 scalar that if we take λ hH = 0, then H and H µτ can represent the physical states (zero mixing between H and H µτ ). But in our case, λ hH = 0 and therefore, we need to introduce two physical states in the following way, where α is the mixing angle. The above mixing matrix becomes diagonal with respect to the physical states h 1 , h 2 and the eigenvalues of M 2 scalar represent the mass terms M h 1 and M h 2 of the scalar fields h 1 and h 2 respectively. In this work, we choose h 1 as the SM-like Higgs boson. The quartic couplings λ h , λ H and λ hH can be expressed in terms of M h 1 , M h 2 , VEVs (v, v µτ ) and mixing angle (α) in the following way, After symmetry breaking, mass of the dark matter is, In the present scenario, two of the three scalars take VEVs while the third one does not, i.e., In addition, the stability of the vacuum requires [43] the following constraints on the quartic couplings of the scalar fields, Besides the lower limit, there is an upper limit on the quartic couplings and the Yukawa couplings due to the perturbative regime, constraining them to be less than 4π and √ 4π, respectively.

III. MUON (g − 2) AND NEUTRINO MASS
In our previous work [32] we studied in detail the muon (g−2) and neutrino mass phenomenology in the framework of present model. The difference here with the previous work is that here we wish to produced the dark matter via the freeze-in mechanism instead of freeze-out, as in [32]. Our aim is to simultaneously explain the dark matter relic abundance as well as the muon (g − 2). The region of the g µτ -M Zµτ parameter space that can explain the muon (g − 2) has been shown in [44]. This allowed region is seen to be very small and mostly constrained by the neutrino trident process experiments such as CHARM-II, CCFR [45,46] and 4-leptons decay [44]. In determining the allowed parameter space for the present model, we have considered the relic density constraint [5] and also the bound on invisible decay width of SM-like Higgs [47]. The one loop contribution [48,49] to muon (g − 2) for the U(1) Lµ−Lτ gauge boson Z µτ mediated diagram is given by, where, r = (M Zµτ /m µ ) 2 , m µ being the muon mass. The discrepancy between the observed and SM predicted value [42] of muon (g − 2) is, In Fig. 1 we show the allowed and disfavoured regions in the g µτ -M Zµτ plane. The region above the green dashed line is ruled out by the neutrino trident experiment CCFR [44]. The grey region inside the blue dashed-dotted line and black dashed line can explain the (g − 2) anomaly in ±2σ range [44]. As will be discussed in much detail later, the red dots in this figure span the parameter region which can satisfy the dark matter relic abundance (cf. Eq. (24)). We see that for g µτ ≥ 3 × 10 −3 no red points exist because the contribution from the Z µτ mediated diagram to the relic abundance becomes too large. See Appendix A, for the expression of cross section of The region of the parameter space compatible with both the dark matter relic abundance and muon (g − 2) lies in the narrow overlapping zone. The benchmark point (values of g µτ and M Zµτ ) used in all further results shown in this work is marked by the star in Fig. 1 and corresponds to M Zµτ = 100 MeV and g µτ = 9 × 10 −4 . Such low mass Z µτ gauge boson can be searched by looking 2µ + / E T final states in LHC or future collider experiments [50,51]. For these values of the parameters, the contribution to muon (g − 2) from Eq. (13) comes out to be which lies within the ±2σ range of the observed value [44]. The neutrino mass generation via the Type-I seesaw has been discussed in detail in [32]. From the neutrino part of the Lagrangian given in the Eq. (2), we can write down the Majorana mass matrix for the RH neutrinos after spontaneous symmetry breaking as, In general M R 2 can be complex but by proper phase rotation all the components can be made real except one and we choose to take M µτ as the complex component [52]. The light neutrino mass matrix is given as, where M D is the Dirac mass term, and here it comes as diagonal. The expression for the neutrino mass matrix is, where p = h eµ h eτ v 2 µτ − M ee M µτ e iξ and f e , f µ and f τ are the diagonal components of M D . In [32], we have shown that for the chosen (g µτ , M Zµτ ) value, we can satisfy all the constraints of the mixing angles, mass square differences [53] and cosmological bound on the sum of the neutrino masses [5]. We also showed the region in the model parameter space which satisfy all the neutrino parameter constraints.

IV. FIMP DARK MATTER AND BOLTZMANN EQUATION
As can be seen from the Eqs. (1-4), the dark matter particle can interact with the thermal bath containing both SM as well as BSM particles only via h 1 , h 2 and Z µτ . The Feynman diagrams relevant for the production of dark matter are shown in Fig. 2. The corresponding couplings are listed in Table III. From Table III, one can see that the couplings of φ DM with scalar bosons h 1 and h 2 depend on the parameters λ Dh and λ DH , while that with extra neutral

Vertex
Vertex Factor gauge boson Z µτ involves g µτ and n µτ . For the dark matter to be a suitable FIMP candidate, the cross section of the diagrams listed in Fig. 2 should be very small. The complete expressions for the contribution from each of the diagrams is given in Appendix A and B. The processes involving h 1 and h 2 can be easily made feeble enough by taking λ Dh and λ DH ∼ 10 −12 . As we will see later, the other important production mechanism of φ DM is shown by the Feynman diagram where N 2 and N 3 annihilate to φ DM via the new gauge boson Z µτ . The expression for the cross section of this process is given in Eq. (A10) in Appendix A. We see that the cross section for this process is proportional to ∼ g 4 µτ n 2 µτ /10 2 . Since we fix g µτ = 9 × 10 −4 to explain the anomalous muon (g − 2), we take n µτ ∼ 10 −5 to keep σ N j N j →φ † DM φ DM small enough so that φ DM stays out of chemical equilibrium. This choice for n µτ also ensures that there is a remnant Z 2 symmetry even when U(1) Lµ−Lτ symmetry is broken spontaneously, which enables φ DM to be stable. Thus, φ DM behaves as a FIMP dark matter, it stays out of thermal equilibrium at all times, and is produced by the freeze-in mechanism.
The evolution of comoving number density of FIMP produced from the decays as well as annihilations of the SM and BSM particles is governed by the Boltzmann equation. The Boltzmann equation in terms of the comoving number density of φ DM is given below. This equation contains both decay as well as annihilation terms. While deriving the Boltzmann equation for the FIMP φ DM , we have taken all the particles except φ DM in thermal equilibrium and hence their number densities follow the Maxwell-Boltzmann distribution function.
In the above equation Y φ DM = n φ DM s is the comoving number density, n φ DM represents the actual number density of the dark matter candidate φ DM while s is the entropy of the Universe. The GeV is the Planck mass. The function g (z) is related to the degrees of freedom and has the following expression, In the above equation g s (z) and g ρ (z) are the effective degrees of freedom corresponding to the entropy and energy densities of the Universe. If the decaying particles (h 1 , h 2 ) are in thermal equilibrium 3 then the thermal average of the decay width can be expressed in a simpler form as follows, In the above expression K 1 (z) and K 2 (z) are the modified Bessel's functions of order one and two, respectively. Next we explain the Boltzmann equation given in Eq. (19) term by term. The first term on the R.H.S represents the contribution to Y φ DM arising from the decays of h 1 and h 2 and it is proportional to the equilibrium comoving number density of the decaying particle. On the other hand, the inverse decay term proportional to the comoving number density of φ DM occurs with a negative sign as it washes out the φ DM number density. However, as we have mentioned before that the initial number density of φ DM is extremely small due to its non-thermal origin and hence the negative feedback coming from the inverse processes can be safely neglected.
Similarly, the second, third and fourth terms in the R.H.S of Eq. (19) indicate the net contribution to Y φ DM coming from the annihilations of SM and BSM particles at the early stage of the Universe. Unlike the first term (decay term) of Eq. (19) these terms are proportional to the second power of the comoving number densities of relevant particles. The term proportional to Y 2 p (p = W, Z, h 1 , h 2 , f ) (in the second term of the Boltzmann equation) represents the increment of Y φ DM from the pair annihilation of p and its antiparticle while the feedback arising from the inverse process φ DM φ † DM → pp is proportional to Y 2 φ DM and comes with a negative sign. The rest of the annihilation terms represent the production processes of φ DM from the annihilations of two different particles such as (N i , N j ), (h 1 , h 2 ) and consequently these terms are proportional to the comoving number densities of two different initial state particles. The terms indicating the inverse processes In the above Eq. (22), s is the Mandelstam variable, M A and M B are the masses of the initial state particles, σ A B→ φ † DM φ DM is the annihilation cross section, while K i (x) is the modified Bessel function of order i. In order to get the comoving number density (Y φ DM ) of the DM particle φ DM , we have to solve the Boltzmann equation numerically. After determining the comoving number density Y φ DM of dark matter particle φ DM at the present epoch, we can determine the relic density [54,55] from the following relation, where M DM is in GeV. The observed value of the dark matter relic density given by the Planck collaboration [5], is In this work we have used the above relic density bound.

V. RESULTS
We have implemented our model in LanHEP [56] to generate all the vertex factors which are required to calculate the relevant annihilation cross sections and decay widths. Corresponding Feynman diagrams are shown in Fig. 2. After putting all the cross sections and decay widths (as listed in Appendix A and Appendix B) in Eq. (19), we solve the Boltzmann equation numerically and study the related phenomenology of FIMP dark matter. Throughout this analysis we keep the extra Higgs mass fixed at M h 2 = 500 GeV. In Fig. 3, we show the variation of the dark matter relic density with z for different choices of the initial temperature. From the figure, it is clear that if the initial temperature is greater T in ≥ 1 TeV, then there is no dependence of the final relic density on the value of the initial temperature. This can be explained in the following way. As the heavy Higgs h 2 is in thermal equilibrium with the cosmic soup, hence the In the left panel of Fig. 4, we show the relative contributions of two different types of production processes (decay and annihilation) to Ωh 2 . The red dotted line represents the contribution from the decay of SM-like Higgs boson h 1 and extra Higgs boson h 2 , while the black dashed line corresponds to the contribution from the all possible annihilation channels of SM and BSM particles (see Fig. 2 for the corresponding Feynman diagrams). The total contribution towards the relic density of φ DM coming from the decay as well as annihilation of different particles is represented by blue dashed-dotted line. The horizontal magenta line indicates the observed value of dark matter relic density (Ωh 2 ∼ 0.12 [5]) at the present epoch. From this plot, it can be seen clearly that for our chosen set of model parameters (M Zµτ = 0.1 GeV, g µτ = 9.0 × 10 −4 , α = 0.01, λ Dh = 9.8 × 10 −13 , λ DH = 1.3 × 10 −11 , M DM = 50.0 GeV, n µτ = 5.5 × 10 −5 ), the decay processes contribute ∼ 67% of dark matter production while rest of the dark matter particles are produced from the annihilations of different SM as well as BSM particles. In the right panel of Fig. 4, the variation of relic density with z (i.e. with respect to the inverse of temperature T ) has been shown for four different values of dark matter mass M DM . For M DM = 10 GeV, 30 GeV and 50 GeV, the relic density is seen to rise with M DM . This agrees with the expression for relic density given in Eq. (23). However, if we take a slightly higher value of dark matter mass, M DM = 60 GeV (blue dotted line), the relic density decreases instead of increasing. This is because M DM = 60 GeV is very close to half of the SM like Higgs boson mass (∼ M h 1 /2) and the decay mode h 1 → φ DM φ † DM becomes phase space suppressed. Therefore, it reduces the contribution arising from h 1 decay and hence the final relic density of dark matter.
The contributions to Ωh 2 arising from the decays of h 1 , h 2 and the annihilations of SM as well as BSM particles are shown respectively in left and right panels of Fig. 5. Here we define a quantity Ω Γ Ω ( Ω σv Ω ) which represents the fractional contribution of a particular decay (annihilation) channel to dark matter relic density. In the left-panel of Fig. 5, the contribution from h 2 decay has been shown by the green dashed-dotted line and that from h 1 decay has been shown by the red dashed line, while the total decay contribution to the dark matter relic density is represented by the black solid line. From this plot one can see that, initially for low values of z (z < 10, corresponding to higher temperatures), the extra Higgs contribution to Ωh 2 is more because of its high mass. On the other hand, for higher values of z (z > 10), the SM-like Higgs decay contribution starts dominating. In the right panel of Fig. 5, we show the contribution coming from different annihilation channels. The total contribution from all the annihilation channels is represented by the black solid line while the other lines show the contribution of individual channels. From this plot it is clearly seen that, the two dominating  annihilation channels are N 2 N 2 and N 3 N 3 . Annihilation channels h 1 h 1 and h 2 h 2 also have significant role in the production processes of dark matter, while the effect of other channels are sub dominant. In the left-panel of Fig. 6, variation of relic density with z for three different value of λ Dh have been shown. The red dashed line for λ Dh = 9.8 × 10 −13 gives the correct relic density. In the right-panel of Fig. 6 we show the variation of relic density for different values of the other quartic coupling λ DH . It is clear from Fig. 6 that the relic density increases with both λ Dh and λ DH as the production modes of φ DM are proportional to these quartic couplings. However, the increment of Ωh 2 with respect to increasing λ Dh or λ DH is not uniform. When we decrease λ Dh (λ DH ) from 9.8 × 10 −13 (1.3 × 10 −11 ) by one order of magnitude, the decrease in Ωh 2 is very small because in this regime we have dominant contribution from the Z µτ mediated right-handed neutrino annihilation channel. However, if we increase λ Dh (λ DH ) from 9.8 × 10 −13 (1.3 × 10 −11 ) by one order of magnitude, Ωh 2 increases by more than order of magnitude since in this case the contribution from decay channels become dominant. In the left-panel of Fig. 7, we have shown the allowed regions in the λ Dh -λ DH plane. The red points in the plane satisfy the relic density bound. Both the parameters λ Dh and λ DH have been varied from 10 −14 to 10 −8 . We see from the figure that for λ Dh ≥ 8 × 10 −11 and λ Dh ≥ 3 × 10 −10 no red points exist, and therefore these regions are disallowed by the relic density bound. One can also notice that there is no lower bound on λ Dh and λ DH . This is because for lower values of λ Dh and λ DH , even though the Higgs mediated annihilation and decay contributions become very less, the Z µτ mediated annihilation channels (N i N i → φ † DM φ DM , i = 2, 3 see Appendix A) contribute fully and hence can explain the relic density bound. In the right-panel of Fig. 7, we show the allowed regions in the M DM -g µτ plane. Here we have varied dark matter mass M DM from 1 GeV to 100 GeV and the U(1) Lµ−Lτ gauge coupling g µτ from 10 −6 to 0.1. The figure shows that the whole range of dark matter mass M DM can satisfy the relic density bound. However, the gauge coupling g µτ ≥ 3 × 10 −3 does not satisfy the relic density bound as over production of φ DM occurs through the annihilation channel N i N i → φ † DM φ DM , i = 2, 3. In Fig. 8, we show the allowed parameter space in the M DM − λ Dh and M DM − λ DH planes in the left and right panels respectively. As it was seen earlier, here too the whole range of considered dark matter mass 4 is allowed. However, there is an upper limit on both λ Dh and λ DH . In both the panels, there exist an anti-correlation between the dark matter mass and the quartic coupling λ Dh(H) since the relic density is proportional to the dark matter mass M DM as well as the coupling constant. Hence, if we increase M DM then to satisfy the relic density, λ Dh(H) must decrease. In the left-panel we observe that around M DM ∼ 62 GeV, there is a rise in λ Dh . This happens because this is the resonance region for the SM-like Higgs and as a result there is little contribution from the decay of SM-like Higgs (phase space suppression). Hence, to satisfy the relic density bound, there is a sudden rise in λ Dh to increase the contribution arising from the decay as well as annihilation processes involving h 1 . Beyond M DM ∼ 62.5 GeV, there is no contribution from the SM-like Higgs. Hence, the coupling λ Dh again starts behaving in the normal way (anti-correlation). No such peculiar behaviour is seen for λ DH because this parameter is important for the decay of h 2 and the chosen mass range of the dark matter is not in the resonance region of the extra Higgs h 2 (M h 2 ∼ 500 GeV). Therefore, for the quartic coupling λ DH , the anti-correlation exists for the entire range of dark matter mass M DM (1-100 GeV).

VI. CONCLUSION
In this work, we propose a framework for non-thermal production of a FIMP type dark matter in the gauged U(1) Lµ−Lτ extension of the SM. The particle content of our model includes two SM gauge singlet scalars and three right-handed neutrinos. Both the SM singlet scalars carry U(1) Lµ−Lτ charge. One of them picks up a VEV, breaking U(1) Lµ−Lτ spontaneously, thereby giving the new gauge boson Z µτ mass. The new boson Z µτ gives additional contribution to muon (g − 2) which can explain its measured value. The U(1) Lµ−Lτ breaking also leads to additional terms in the light neutrino mass matrix which has its origin via the Type-I seesaw mechanism. This model can therefore explain consistently the neutrino masses and mixing patterns observed in neutrino oscillation experiments. Finally, the other SM singlet scalar φ DM does not acquire any VEV and remains stable due to a remnant Z 2 symmetry even after the breaking of U(1) Lµ−Lτ due to our choice of its L µ − L τ charge (n µτ ). This scalar can therefore serve as a dark matter candidate.
In this model, φ DM can interact with the visible sector (including both SM as well as BSM particles) only through the scalar bosons h 1 , h 2 and L µ − L τ gauge boson Z µτ . Hence, in order to keep φ DM out of thermal equilibrium we have considered the corresponding couplings of φ DM to be extremely feeble. All other particles (both SM and BSM) in the early Universe, however remain in equilibrium with the thermal soup. We have numerically solved the Boltzmann equation containing all possible production modes of φ DM from decays and annihilations of SM and BSM particles. We have found that, among all of its production modes, φ DM is mainly produced from the decays of h 1 , h 2 and the pair annihilations of right handed neutrinos mediated by Z µτ . The latter process beautifully sets the interplay between the dark matter sector, neutrino masses and muon (g − 2). The solution of Boltzmann equation also indicates that, for our chosen set of model parameters, decay processes of h 1 and h 2 contribute ∼ 67% of the total dark matter relic abundance while rest ∼ 33% contribution is coming from the annihilation of particles (mainly right-handed neutrinos) in the thermal bath. We have shown that, the relic abundance of our FIMP type dark matter candidate φ DM lies within the Planck Limit (0.1172 ≤ Ωh 2 ≤ 0.1224) only when the relevant parameters satisfy the following bounds: λ Dh < ∼ 10 −11 , λ DH < ∼ 10 −10 and g µτ ≤ 3 × 10 −3 for the dark matter mass range of 1 GeV to 100 GeV and n µτ = 5.5 × 10 −5 . Consequently, due to the low coupling strength there does not exist any bound on dark matter spin independent scattering cross section from direct detection experiments.
In conclusion, our proposed U(1) Lµ−Lτ extension of the SM can explain the three main puzzles that demand beyond Standard Model physics. First, it can successfully explain the smallness of neutrino masses via the Type-I seesaw mechanism as well as the peculiar mixing pattern of the neutrinos via the U(1) Lµ−Lτ gauge symmetry that also acts on the lepton flavours, thereby giving a pattern to the light neutrino mass matrix. Second, the additional one loop contribution of the extra neutral gauge boson Z µτ can successfully satisfy the muon (g − 2) data. And finally, the model has a SM singlet scalar with non-zero U(1) Lµ−Lτ , that makes it stable and which acts as a non-thermal dark matter candidate, thereby satisfying constraint on the relic abundance and at the same time evading all bounds coming from direct and indirect dark matter detection experiments. The annihilation of the vector bosons W + and W − into the dark matter φ DM and φ † DM are occurred by the two s channel diagrams which are mediated by the SM-like Higgs h 1 and the SM extra Higgs boson h 2 . Expression of this annihilation cross section is, where |M W W | 2 represents the square of the absolute value of M W W . The couplings g h 1 φ † DM φ DM and g h 2 φ † DM φ DM are given in Eq. (A1) while Γ h 1 and Γ h 2 are the total decay widths of the SM-like Higgs and the extra Higgs respectively (see Appendix B for more details). In the expression of M W W the extra factor 1/9 comes due to average of the polarisation vectors of two initial state W bosons.
The self annihilation of Z boson to dark matter particles φ DM and φ † DM has two s channel diagrams, one of them is mediated by the SM-like Higgs h 1 while another one is mediated by the extra Higgs h 2 . Cross section structure for the Z Z annihilation is similar to W + W − annihilation and has the following form, The annihilation of two h 1 to φ † DM φ DM occurs through one four points interaction process as well as two s channel processes 5 mediated by the SM-like Higgs h 1 and extra Higgs h 2 . The expression of the cross section is the following, Similarly, the annihilation of two h 2 to φ † DM φ DM also occurs through one four points interaction process as well as two s channel processes mediated by the SM like Higgs h 1 and extra Higgs h 2 . The expression of annihilation cross section for h 2 h 2 → φ † DM φ DM is given by Like the previous two cases, here we also have three processes (one four points interaction and two s channels) contributing to the annihilation of h 1 h 2 → φ † DM φ DM . The expression of h 1 h 2 → φ † DM φ DM has the following form Here, the annihilating particles are t,t and the final particles are φ DM , φ † DM . This annihilation is occurred only by the two s channel diagrams mediated via h 1 and h 2 respectively. The expression of the cross section for this process is, where M t is the mass of the top quark and n c = 3 is its color charge.
• N 1 N j → φ † DM φ DM (j = 2, 3) : The annihilation of N 1 and N j to φ † DM φ DM has two s channel diagrams mediated by h 1 and h 2 respectively. The corresponding expression of the cross section is, Unlike the previous cases, here the annihilation of N j N j occurs only through a s channel process mediated by the extra neutral gauge boson Z µτ . The expression of the cross section is, .

(A10)
Appendix B: Expressions of decay widths of h 2 , h 1 and Z µτ Total decay width of h 2 : Decay width for a particular process is generally calculated in the rest frame of the corresponding decaying particle. In the present work, the BSM Higgs h 2 can decay to all the SM particles and also to other BSM particles like Z µτ , φ DM and N j . Here, we have given the expressions of all partial decay widths as well as the total decay width of h 2 .
• h 2 → V V (V = W ± , Z): The width of this decay process is given as follows In the above decay expression S V is the statistical factor. It is 1 for W ± boson and 2 for Z boson.
Finally, the total decay width of the BSM Higgs h 2 is thus the sum of all partial decay widths 6 mentioned above , which is Decay width of h 1 : Besides decaying to SM particles, h 1 has extra decay modes to a pair of Z µτ and φ DM respectively. The expressions of these decay channels are, where the expression of g h 1 φ † DM φ DM is given in Eq. (A1). Therefore, total decay width of the SM-like Higgs h 1 can be written as, Decay width of Z µτ : Since in this work we have considered the low mass of Z µτ (M Zµτ ∼ 100 MeV) and also it has no couplings to quarks, hence it can only decay to neutrinos. Therefore expression of total decay width of Z µτ is given by [1] Y. Sofue and V. Rubin, "Rotation curves of spiral galaxies" Ann. Rev. Astron. Astrophys. 39, 137 (2001) [astro-ph/0010594]. 6 Here we have taken only those decay modes of h 2 which occur in tree level.