TeV Scale Resonant Leptogenesis with $L_{\mu}-L_{\tau}$ Gauge Symmetry
in the Light of Muon $(g-2)$

Here we consider the popular and minimal model based on the gauged $L_{\mu}-L_{\tau}$ symmetry which is anomaly free He et al. (1991a, b). Apart from the SM fermion content, the minimal version of this model has three heavy right handed neutrinos (RHN) leading to to type I seesaw origin of light neutrino masses Minkowski (1977); Mohapatra and Senjanovic (1980); Yanagida (1979); Gell-Mann et al. (1979); Glashow (1980); Schechter and Valle (1980). The same RHNs can lead to leptogenesis Fukugita and Yanagida (1986); Davidson et al. (2008) via out-of-equilibrium decay into SM leptons. However, with hierarchical RHN spectrum, there exists a lower bound on the scale of leptogenesis, known as the Davidson-Ibarra bound $M_{1}>10^{9}$ GeV Davidson and Ibarra (2002). Several earlier works Adhikary (2006); Chun and Turzynski (2007); Ota and Rodejohann (2006); Asai et al. (2017, 2020) considered different scenarios like supersymmetry, high scale leptogenesis within type I seesaw framework of this model. However, it is also possible to have TeV scale leptogenesis via resonant enhancement due to tiny mass splitting between RHNs, known as the resonant leptogenesis Pilaftsis (1999); Pilaftsis and Underwood (2004); Moffat et al. (2018); Dev et al. (2018a). It should be noted that earlier works on leptogenesis in gauged $L_{\mu}-L_{\tau}$ model considered high scale breaking of such gauge symmetry and compatibility with muon $(g-2)$ explanation from a low scale vector boson was missing. We intend to perform a general analysis as well as to bridge this gap showing the possibility of TeV scale leptogenesis along with muon $(g-2)$ from a light $L_{\mu}-L_{\tau}$ vector boson.

Motivated by the recent measurements of the muon $(g-2)$, in this work we consider the possibility of low scale leptogenesis and constrain the $L_{\mu}-L_{\tau}$ gauge sector from the requirements of successful leptogenesis and $(g-2)_{\mu}$. Since the decaying RHNs can have this new gauge interaction which can keep them in equilibrium for a longer epochs and can also initiate some washout processes, the requirement of successful leptogenesis for a fixed scale of leptogenesis around a TeV can lead to constraints on the gauge sector parameter space. Since leptogenesis is a high scale phenomena and the requirement of $(g-2)_{\mu}$ needs a low scale breaking $L_{\mu}-L_{\tau}$ gauge symmetry, one can not realise both in the minimal version of the model. We first discuss the minimal model from the requirement of satisfying neutrino data and baryon asymmetry from leptogenesis and then consider an extended model which can accommodate muon $(g-2)$ as well. While we do not pursue the study of $R_{K}$ anomalies in this model, one may refer to Biswas and Shaw (2019) for common origin of muon $(g-2)$ and $R_{K}$ anomalies along with dark matter in extensions of minimal $L_{\mu}-L_{\tau}$ model. Explanation of similar flavour anomalies along with muon $(g-2)$ in this model have also been studied Crivellin et al. (2015); Altmannshofer et al. (2016). Recently, a dark matter extension of the $L_{\mu}-L_{\tau}$ model was also found to provide a common origin of muon $(g-2)$ and electron recoil excess reported by the XENON1T collaboration Borah et al. (2020, 2021).

The new field content apart from the SM ones are shown in table 1. Only the second and third generations of leptons are charged under the $L_{\mu}-L_{\tau}$ gauge symmetry. The relevant Lagrangian can be written as

where $H$ is the SM Higgs doublet. The covariant derivatives for RHNs are defined by

While the neutral component of the Higgs doublet $H$ breaks the electroweak gauge symmetry, the singlets break $L_{\mu}-L_{\tau}$ gauge symmetry after acquiring non-zero vacuum expectation values (VEV). Denoting the VEVs of singlets $\Phi_{1,2}$ as $v_{1,2}$, the new gauge boson mass can be found to be $M_{Z_{\mu\tau}}=g_{\mu\tau}\sqrt{(v^{2}_{1}+4v^{2}_{2})}$ with $g_{\mu\tau}$ being the $L_{\mu}-L_{\tau}$ gauge coupling. Clearly the model predicts diagonal charged lepton mass matrix $M_{\ell}$ and diagonal Dirac Yukawa of light neutrinos. Thus, the non-trivial neutrino mixing will arise from the structure of right handed neutrino mass matrix $M_{R}$ only which is generated by the chosen scalar singlet fields. The right handed neutrino mass matrix, Dirac neutrino mass matrix and charged lepton mass matrix are given by

Here $v$ is the VEV of neutral component of SM Higgs doublet $H$. We can find the light neutrino mass matrix can be found by the type I seesaw formula as

The mass matrices for RHN as well as light neutrinos in this model do not possess any specific structure and hence it is straightforward to fit the neutrino oscillation data. In the presence of only one scalar singlet $\Phi_{1}$, the RHN mass matrix has zero entries at $(\mu\mu)$ and $(\tau\tau)$ entries leading to a two-zero minor structure Asai et al. (2017) in $M^{-1}_{\nu}$. The light neutrino mass matrix, although does not contain any zeros, leads to two constraints among its elements. In view of tight constraints on neutrino parameters from global fit data Esteban et al. (2019); de Salas et al. (2021); Zyla et al. (2020) as well as cosmology bounds on sum of light neutrino masses from Planck 2018 data Aghanim et al. (2018), it is difficult to satisfy the data using the constrained structure of mass matrices. Similar conclusion was also arrived at in earlier works Asai et al. (2017, 2020). This becomes more restrictive when we constrain two of the RHNs to be in TeV regime with tiny mass splittings for resonant leptogenesis. Therefore, we have introduced another singlet scalar $\Phi_{2}$ which gives rise to general structures of $M_{R}$ and $M_{\nu}$.

It should be noted that, a kinetic mixing term between $U(1)_{Y}$ of SM and $U(1)_{L_{\mu}-L_{\tau}}$ of the form $\frac{\epsilon}{2}B^{\alpha\beta}Y_{\alpha\beta}$ can exist in the Lagrangian where $B^{\alpha\beta}=\partial^{\alpha}X^{\beta}-\partial^{\beta}X^{\alpha},Y_{\alpha\beta}$ are the field strength tensors of $U(1)_{L_{\mu}-L_{\tau}},U(1)_{Y}$ respectively and $\epsilon$ is the mixing parameter. Even if this mixing is considered to be absent in the Lagrangian, it can arise at one loop level with particles charged under both the gauge sectors in the loop. We consider this mixing to be $\epsilon=g_{\mu\tau}/70$. While the phenomenology of muon $(g-2)$, leptogenesis in our model is not much dependent on this mixing, the experimental constraints on the model parameters can crucially depend upon this mixing. We therefore choose it to be small, around the same order as its one loop value.

where $V_{ij}$ are the elements of the matrix $V$ which diagonalises $M_{R}$.

The CP asymmetry parameter corresponding to the CP violating decay of RHN $N_{i}$ (summing over all lepton flavours) is given by Pilaftsis and Underwood (2004)

Since we are mainly focusing on the parameter space where $M_{1}\backsimeq M_{2}<M_{3}$ such that the leptogenesis is mainly governed by the resonant enhancement between $N_{1}$ and $N_{2}$ and the CP asymmetry coming from the decay of $N_{3}$ is negligible. The resonance condition is satisfied when $M_{2}-M_{1}\simeq\Gamma_{1}/2$ where $\Gamma_{1}$ is the decay width of the lightest RHN $N_{1}$. The relevant asymmetry parameters namely, $\epsilon_{1}$ and $\epsilon_{2}$ crucially depend upon $Im[(h^{\dagger}h)^{2}_{12}]$ and $Im[(h^{\dagger}h)^{2}_{21}]$ respectively. From equation (13) one can write

From equation (LABEL:eq:complex1) the importance of the phases appearing in the elements of the matrix $V$ can be seen. These phases are appearing in the diagonalising matrix $V$ from the complex parameters in $M_{R}$.

The relevant Boltzmann equations for our setup can be written as

where $n_{N_{i}}$ and $n_{B-L}$ denote the comoving number densities of $N_{i}$ and $B-L$ respectievly. The equilibrium no densities of $N_{i}$’s are defined by $n_{N_{i}}^{\rm eq}=\frac{z^{2}}{2}\kappa_{2}(z)$, (with $\kappa_{i}(z)$ being the modified Bessel function of $i$-th kind) and $z=M_{1}/T=M_{2}/T$ (for resonant leptogenesis). In $\langle\sigma v\rangle_{N_{i}N_{i}\rightarrow XX}$, $X$ denotes any final state particle to which $N_{i}$’s can annihilate into. $D_{1,2}$ are the decay terms and $WID_{1,2}$ are the inverse decay terms for $N_{1}$ and $N_{2}$ decays respectively, which are measures of the rate of decay and inverse decay with respect to the background expansion of the universe parametrised by Hubble expansion rate $\bf H$. They are defined as

with $K_{i}=\Gamma_{i}/{\bf H}(z=1)$ being the decay parameter. The term $\Delta W$ on the right hand side of the equation (20) takes account of all the scattering processes that can act as possible washouts for the generated $B-L$ asymmetry. We identify the following scattering washouts in our model, $lW^{\pm}(Z)\longrightarrow N_{1,2}H$, $lZ_{\mu\tau}\longrightarrow N_{1}H$, $ql\longrightarrow qN_{1,2}$, $lN_{1,2}\longrightarrow qq^{c}$, $lH\longrightarrow l^{c}H^{*}$, $lH\longrightarrow N_{1,2}W^{\pm}(Z)$ and $lN_{1,2}\longrightarrow Z_{\mu\tau},H$. We take all of them into account in our numerical analysis and defined $\Delta W$ as,

Thus, apart from the usual Yukawa or SM gauge coupling related processes, we also have washout processes involving $Z_{\mu\tau}$ gauge boson. Since interactions involving $Z_{\mu\tau}$ can cause dilution of $N_{1}$ abundance as well as wash out the generated lepton asymmetry, one can tightly constrain the $L_{\mu}-L_{\tau}$ gauge sector couplings from the requirement of successful leptogenesis at low scale. Similar discussions on impact of such Abelian gauge sector on leptogenesis can be found in Iso et al. (2011); Okada et al. (2012); Heeck and Teresi (2016); Dev et al. (2018b); Mahanta and Borah (2021); Fileviez Pérez et al. (2021).

After solving the above Boltzmann equations, we convert the final $B-L$ asymmetry $n_{B-L}^{f}$ just before electroweak sphaleron freeze-out into the observed baryon to photon ratio by the standard formula

where $a_{\rm sph}=\frac{8}{23}$ is the sphaleron conversion factor (taking into account two Higgs doublets). Here $g_{*},g_{*}^{0}$ denote the relativistic degrees of freedom at scale of leptogenesis $T=M_{1}$ and scale of recombination respectively. The final baryon to photon ratio is then compared with Planck 2018 data Zyla et al. (2020); Aghanim et al. (2018)

and constraints on the model parameters are obtained. In figure 1, we show the evolution of lepton asymmetries for different benchmark values of gauge coupling $g_{\mu\tau}$ (left panel) and gauge boson mass $M_{Z_{\mu\tau}}$ (right panel) while keeping other parameters fixed. We can see that the asymmetry decreases with the increase in $g_{\mu\tau}$. It is because with the increase in $g_{\mu\tau}$ the annihilations of $N_{1,2}$ through $Z_{\mu\tau}$ and to a pair of $Z_{\mu\tau}$ increases, which tries to bring the $N_{1,2}$ number densities close to their equilibrium densities. Thus, depletion in number densities of $N_{1,2}$ leads to a decrease in the final asymmetry.

Again, from figure 1, we can see that the asymmetry decreases with increase in $M_{Z_{\mu\tau}}$. This is because we have two types of annihilation of $N_{1,2}$ involving $M_{Z_{\mu\tau}}$, the first one being the annihilation into a pair of leptons and the second one is the annihilation into a pair of $M_{Z_{\mu\tau}}$. Out of these two processes, the dominant one is the annihilation into a pair of leptons as can be seen in figure 2 where the two annihilation cross-sections are shown as a function of center of mass energy squared $s$ for comparison. The annihilation of $N_{1}$ into a pair of lepton is possible through a s-channel diagram mediated by $Z_{\mu\tau}$ and through a t-channel diagram mediated by scalar doublet. The s-channel diagram mediated by singlet scalar due to mixing with the SM like Higgs is ignored as such mixing is generated only after electroweak symmetry breaking. As we are working in the region of parameter space where $M_{Z_{\mu\tau}}<2M_{1}$, therefore with increase in $M_{Z_{\mu\tau}}$ the cross sections for the $Z_{\mu\tau}$ mediated processes increase as seen from the right panel plot of figure 2. This tends to bring the $N_{1,2}$ number density closer to its equilibrium density and therefore leading to a decrease in $B-L$ asymmetry. We show this effect in figure 3. Clearly, increase in gauge coupling takes the number density of $N_{1}$ closer to its equilibrium number density, as seen from the left panel plot of figure 3 . Increase in $M_{Z_{\mu\tau}}$ also has a similar but milder effect, as seen from the right panel plot of figure 3.

In figure 4, we summarise the results in $g_{\mu\tau}-M_{Z_{\mu\tau}}$ parameter space considering $Z_{\mu\tau}$ mass larger than a GeV. The pink coloured band is the favoured region in $g_{\mu\tau}-M_{Z_{\mu\tau}}$ plane from the Fermilab’s result on muon (g-2). The yellow coloured region is excluded by the upper bound on neutrino trident process measured by the CCFR collaboration Altmannshofer et al. (2014). Therefore one can clearly see that the muon $(g-2)$ favoured region in the high mass regime of $Z_{\mu\tau}$ is completely excluded from the bound on neutrino trident process. Additionally, the LHC bounds searches for multi-lepton final states signatures also rule out some part of the parameter space . We show the exclusion region from LHC measurements of $Z\rightarrow 4\mu$ Sirunyan et al. (2019) in terms of the dotted line. The region below the blue dashed line is allowed by HFAG lepton universality test at $2\sigma$ level Amhis et al. (2021). The leptogenesis favoured parameter space for two different scales of leptogenesis $M_{1}=1.38$ TeV, $M_{1}=1.95$ TeV are shown in blue and green coloured bands respectively. As $M_{1}$ increases, the annihilation rates of $N_{1}$ decreases, allowing slightly larger values of $g_{\mu\tau}$ for fixed $Z_{\mu\tau}$ mass. On the other hand, as $Z_{\mu\tau}$ mass is increased towards $2M_{1}$, the annihilation rate of $N_{1}$ mediated by $Z_{\mu\tau}$ increases, requiring smaller values of $g_{\mu\tau}$ in order to take $N_{1}$ out of equilibrium resulting in generation of lepton asymmetry.

If sub-GeV regime of $Z_{\mu\tau}$ is explored, there exists an allowed region in $g_{\mu\tau}-M_{Z_{\mu\tau}}$ parameter space consistent with muon $(g-2)$ and other experimental bounds as discussed recently in Borah et al. (2020); Zu et al. (2021); Amaral et al. (2021); Zhou (2021); Borah et al. (2021). As we will see below, this allowed region correspond to $M_{Z_{\mu\tau}}\sim 10-100$ MeV, $g_{\mu\tau}\sim 10^{-4}-10^{-3}$. However, such sub-GeV $M_{Z_{\mu\tau}}$ will correspond to low scale breaking of $L_{\mu}-L_{\tau}$ gauge symmetry around 100 GeV, even below the sphaleron transition temperature. Prior to the symmetry breaking scale, the RHN mass matrix has very restrictive structure as can be realised by considering vanishing singlet scalar VEVs in $M_{R}$ given in (7). The fact that the model is inconsistent with successful leptogenesis in the $\mu-\tau$ symmetric limit was also noted in earlier works Adhikary (2006). Thus, one needs to go beyond the minimal setup in order to find a consistent picture accommodating both leptogenesis and muon $(g-2)$ while agreeing with experimental data including light neutrino mass and mixing.

In order to achieve successful leptogenesis with sub-GeV $Z_{\mu\tau}$, we extend the minimal model with two additional Higgs doublets as shown in table 2. The presence of these additional Higgs doublets $H_{2,3}$ are required to allow corresponding Dirac Yukawa couplings with $N_{\mu,\tau}$ whose $L_{\mu}-L_{\tau}$ charges are chosen as $\pm x$ with $x<1$. While the Higgs doublet $H_{1}$ is neutral under this additional gauge symmetry, $H_{2,3}$ are charged. Singlet scalar charges are chosen to be $x,2x$ respectively. Since $L_{\mu}-L_{\tau}$ gauge boson mass will be proportional to $x$ as well, in addition to $g_{\mu\tau}$ and singlet VEVs, one can have symmetry breaking scale above the scale of leptogenesis while having $g_{\mu\tau}-M_{Z_{\mu\tau}}$ in the desired range for explaining muon $(g-2)$ by choosing small $x$.

With this new particle content in table 2, the singlet fermion Lagrangian remain same as before (except the new gauge charges of $N_{\mu,\tau}$) while the neutrino Dirac Yukawa terms change to the following

While the right handed neutrino mass matrix $M_{R}$ has the same structure as before, the Dirac neutrino mass matrix is changed to

Here we define $\tan\beta=v_{2}/v_{1}=v_{3}/v_{1}$, where $v_{1}$, $v_{2}$ and $v_{3}$ are the VEV of neutral components of Higgs doublets $H_{1}$, $H_{2}$ and $H_{3}$ respectively. We are assuming $v_{2}=v_{3}$ for simplicity so that $v_{1}\sqrt{1+2\tan^{2}\beta}=246$ GeV. While deviation from this $v_{2}=v_{3}$ assumption can lead to different phenomenology of such multi-Higgs doublet models, we do not expect leptogenesis results to change significantly and hence stick to this simple limit. Similarly, assuming the singlet VEVs to be identical to $u$, the $Z_{\mu\tau}$ mass can be derived as

The model we adopt here has two SM-singlet scalars and three Higgs doublets. Since leptogenesis occurs below the scale of $U(1)_{L_{\mu}-L_{\tau}}$ symmetry breaking, the singlet scalars are massive at the scale of leptogenesis. While they do not play any role in leptogenesis, their masses can still be constrained from experimental data. Bounds on such singlet scalars arise primarily due to its mixing with the SM like Higgs boson Robens and Stefaniak (2015); Chalons et al. (2016). The strongest bound on singlet scalar-SM like Higgs mixing angle ($\theta_{m}$) comes form $W$ boson mass correction López-Val and Robens (2014) at NLO for $250{\rm~{}GeV}\lesssim m_{s_{i}}\lesssim 850$ GeV as ($0.2\lesssim\sin\theta_{m}\lesssim 0.3$) where $m_{s_{i}}$ is the mass of singlet scalar $s_{i}$. On the other hand, for $m_{s_{i}}>850$ GeV, the bounds from the requirement of perturbativity and unitarity of the theory turn dominant which gives $\sin\theta_{m}\lesssim 0.2$. For lower values singlet masses $m_{s_{i}}<250$ GeV, the LHC and LEP direct search Khachatryan et al. (2015); Strassler and Zurek (2008) and measured Higgs signal strength Strassler and Zurek (2008) restrict the mixing angle $\sin\theta_{m}$ dominantly ($\lesssim 0.25$). The bounds from the measured value of EW precision parameter are mild for $m_{s_{i}}<1$ TeV. Since singlet scalars do not play any role in leptogenesis, we can tune their couplings with the SM like Higgs so that the mixing angle remains as small as required. On the other hand, the physical scalars arising from Higgs doublets are also constrained from experiments. While they remain massless at the scale of leptogenesis due to unbroken electroweak symmetry, several studies like Mühlleitner et al. (2017); Haller et al. (2018); Misiak and Steinhauser (2017); Arhrib et al. (2018) and references therein, have studied the phenomenology of such multi-Higgs doublet models. While the bounds from existing experiments depend upon the type of such models, the lower bounds on additional scalar masses, arising out of such additional Higgs doublets can range from around 100 GeV to a few hundred GeVs. Since we can satisfy these bounds by appropriate tuning of scalar potential parameters and without affecting the results related to leptogenesis, we do not elaborate them further in this work.

We then follow the same procedure as above by diagonalising $M_{R}$ and rewrite $M_{D}$ in the diagonal $M_{R}$ basis followed by solving the Boltzmann equations numerically to find the asymmetry. We keep the values of $g_{\mu\tau}$ and $M_{Z_{\mu\tau}}$ in the region favoured by muon $(g-2)$. With such small values of $g_{\mu\tau}$ and $M_{Z_{\mu\tau}}$ while keeping the singlet VEVs above the masses of the RHNs, the required value of parameter $x$ becomes less than unity reducing the interaction strength of the RHNs with $Z_{\mu\tau}$ compared to the earlier model. As a consequence, the annihilation of RHNs into a pair of $Z_{\mu\tau}$ or other annihilation mediated by $Z_{\mu\tau}$ remain sub-dominant and hence do not play any significant role in our analysis. Therefore, $N_{1,2}$ abundances are primarily determined by their decays.

Similarly, in the muon $(g-2)$ favoured regime of $g_{\mu\tau}-M_{Z_{\mu\tau}}$, the washouts involving $Z_{\mu\tau}$ are very feeble and do not produce any significant effect on asymmetry. The dominant washout processes in this regime are the inverse decays and the scattering $ql\longrightarrow qN_{1,2}$. In figure 5 we show the evolution of $N_{1}$ abundance for different benchmark parameters. Since $N_{1}$ abundance is primarily governed by decays in this regime, we do not see much differences due to change in model parameters like mass splitting $\Delta M$ or $\tan\beta$. Similarly, in figure 6, we show the evolution of $n_{B-L}$ for different combination of mass splitting $\Delta M$ (left panel) and $tan\beta$ (right panel). Clearly, the asymmetry is maximum for smallest mass splitting. It is because of the resonance enhancement of the asymmetry. The decay widths $\Gamma_{1,2}$ are of the order of $10^{-10}-10^{-9}$ GeV, within the parameter space we are interested in and as the mass splitting reaches near that value the CP asymmetry increases sharply to $\mathcal{O}(1)$ due to the resonance enhancement condition $\Delta M\sim\Gamma/2$ mentioned earlier. On the right panel of figure 6 it is seen that asymmetry increases for larger values of $v_{1}\tan\beta$. Once we fix the entries of the Dirac mass matrix from neutrino oscillation data, the Yukawa couplings ($Y_{D\mu}$ and $Y_{D\tau}$) relevant for leptogensis are very tightly constrained to $v_{1}tan\beta$. Since $v_{1}\tan\beta=v_{2}=v_{3}$, one requires smaller values of Yukawa couplings ($Y_{D\mu}$ and $Y_{D\tau}$) for larger values of $v_{1}tan\beta$ to satisfy light neutrino data, for a fixed scale of leptogenesis. As a consequence of that the decay widths of $N_{1,2}$ decrease and also the corresponding inverse decay rates. The small difference in the decay rates of $N_{1}$ is marginally visible in the evolution plots of $n_{N_{1,2}}$ of figure 5 (right panel) as well. Since lepton asymmetry is generated primarily due to resonant enhancement, decrease in Yukawa coupling (or increase in $v_{1}tan\beta$) decreases the inverse decay rates which increases the $B-L$ asymmetry as can be seen from the right panel plot of figure 6.

Due to the feeble gauge portal interactions in the low mass regime of $Z_{\mu\tau}$, we do not find any strong correlation in $g_{\mu\tau}-M_{\mu\tau}$ plane from leptogenesis unlike in the previous model. Nevertheless, it is observed that low scale leptogenesis is possible in the region of $g_{\mu\tau}-M_{Z_{\mu\tau}}$ plane favoured by muon ($g-2$). We show a summary plot in figure 7 showing the allowed parameter space in $g_{\mu\tau}-M_{Z_{\mu\tau}}$ plane and show three possible benchmark points for TeV scale leptogenesis which also satisfy the light neutrino data.

In equation (29) we show the numerical value of the Yukawa matrix $h$ for first benchmark point of figure 7.

Since the viability of resonant leptogenesis crucially depends upon quasi-degenerate right handed neutrinos, we show the structure of the RHN mass matrix too corresponding to the first benchmark point of figure 7 as

in the units of TeV. Accordingly, the bare mass terms $M_{ee},M_{\mu\tau}$ as well as the singlet scalar Yukawa couplings with the RHNs can be tuned. The above mass matrix gives rise to two quasi-degenerate RHNs with masses $M_{1}=M_{2}-\Delta M=0.53$ TeV with mass splitting $\Delta M=1$ keV and $M_{3}=1.6$ TeV.

In summary plot of figure 7, the parameter space corresponding to Fermilab’s muon $(g-2)$ data is shown between the dashed pink coloured lines. Several experimental constraints are also shown in the summary plot of figure 7. The cyan coloured exclusion band labelled as CCFR correspond to the upper bound on neutrino trident process measured by the CCFR collaboration Altmannshofer et al. (2014). In the high mass regime for $Z_{\mu\tau}$ the exclusion region labelled as BABAR corresponds to the limits imposed on four muon final states by the BABAR collaboration Lees et al. (2016). The astrophysical bounds from cooling of white dwarf (WD) Bauer et al. (2020); Kamada et al. (2018) excludes the upper left triangular region. Very light $Z_{\mu\tau}$ is ruled out from cosmological constraints on effective relativistic degrees of freedom Aghanim et al. (2018); Kamada et al. (2018); Ibe et al. (2020); Escudero et al. (2019). This arises due to the late decay of such light gauge bosons into SM leptons, after standard neutrino decoupling temperatures thereby enhancing $N_{\rm eff}$. More stringent constraints apply if DM is in the sub-GeV regime and its thermal relic is dictated by annihilation mediated by $Z_{\mu\tau}$. However, in our case, such additional constraints do not arise as DM is much heavier and uncharged under $L_{\mu}-L_{\tau}$ gauge symmetry. The observation of coherent elastic neutrino-nucleus cross section in liquid argon (LAr) and cesium-iodide (CsI) performed by the COHERENT Collaboration Akimov et al. (2017, 2021) also leads to constraint on $g_{\mu\tau}-M_{Z_{\mu\tau}}$ parameter space. Adopting the analysis of Cadeddu et al. (2021); Banerjee et al. (2021), we impose these bounds, labelled as COHERENT LAr and COHERENT CsI respectively in figure 7. Clearly, even after incorporating all existing experimental bounds, there still exists a small parameter space between a few MeV to around 100 MeV consistent with all bounds and the requirement of explaining muon $(g-2)$. The three points denoted by stars correspond to benchmark values of parameters for which the correct baryon asymmetry can be generated via resonant leptogenesis. Therefore, we can have successful TeV scale resonant leptogenesis while also satisfying the Fermilab’s data on muon (g-2). This allowed region from muon (g-2) and also TeV scale leptogenesis remains within the reach of future experiments like NA62 at CERN Krnjaic et al. (2020) (Orange dashed line in figure 7), the NA64 experiment at CERN (dot-dashed line of green colour in figure 7) Gninenko et al. (2015); Gninenko and Krasnikov (2018). Possible future confirmation of this muon $g-2$ favoured parameter space will also indicate the possible scale of leptogenesis in this model.