A model for lepton flavor violating non-standard neutrino interactions

We present a model for Lepton Flavor Violating (LFV) neutral current non-standard interactions of neutrinos with matter fields parameterized by $\epsilon_{\alpha \beta}^f$ with $\alpha \ne \beta$. Here, unlike the previous models, the ratios of the off-diagonal LFV elements of the effective NSI coupling to the diagonal lepton flavor conserving ones ({\it i.e.,} $(\epsilon_{\alpha \beta}^f)^2/(\epsilon_{\alpha \alpha}^f \epsilon_{\beta \beta}^f)$ ) are arbitrary. The model enjoys rich phenomenology, predicting invisible Higgs decay and new meson decay modes observable in upcoming experiments. The model for $\epsilon_{\mu e}^f$ also predicts a $\mu^-$ to $e^-$ conversion rate on nuclei accessible in the planned experiments.


I. INTRODUCTION
Since the establishment of lepton flavor violation in solar and atmospheric neutrino data, a wide program for the measurement of the parameters of neutrino mass matrix has been started and is vehemently going on. We are now entering neutrino precision era with upcoming experiments being sensitive to the small subdominant effects in the neutrino oscillation. These experiments aim to extract the yet unknown neutrino oscillation parameters, especially, the Dirac CP-violating phase, δ CP . This has also instilled wide interest in the neutral current Non-Standard neutrino Interaction (NSI) with matter fields, f ∈ {e, quarks} parameterized as the following effective potential where ǫ f αβ are dimensionless parameters. In the limit ǫ f αβ → 0, the standard model is recovered. In the presence of NSI, the propagation of neutrinos in a medium will be affected. There is rich literature studying the effects of NSI on different neutrino oscillation experiments [1]. It has been shown that if ǫ f αβ are relatively large, NSI will induce degeneracies in the parameter space which should be taken into account when the values of neutrino mixing parameters are extracted from observation. In particular, it has been shown that NSI can even mimic the effects of δ CP [2]. Moreover, it has been shown that neglecting NSI may lead to a wrong determination of the θ 23 octant [3].
Of course, if |ǫ f αβ | ≪ 1, its effects will be negligible. The degeneracies that we mentioned above appear only for large values of ǫ f αβ . The natural question is that whether we can make a viable SU (2) × U (1) invariant model which gives rise to NSI with such large ǫ f αβ without violating the myriad of bounds that already exist. As shown in [4][5][6][7][8], invoking a U ′ (1) gauge symmetry with a light gauge boson, Z ′ with a mass below ∼ 100 MeV coupled to both quarks and neutrinos, we can build such models. Building models for lepton flavor conserving NSI (i.e., ǫ f αα ) is relatively easy as we can just gauge a linear combination of lepton flavors and the Baryon number. However, obtaining sizable lepton flavor violating (LFV) ǫ f αβ with α = β is more challenging [9]. If the NSI does not break SU (2) × U (1), the corresponding charged leptons also receive LFV couplings to Z ′ , leading to fast l − α → l − β Z ′ at a tree level with a rate enhanced by (m l − α /m Z ′ ) 2 due to the longitudinal component of Z ′ [5]. This problem is overcome in [6] by introducing a new fermion singlet under SM but charged under the new U ′ (1) which mixes with active neutrinos ν α and ν β through a new Higgs doublet whose Vacuum Expectation Value (VEV) breaks both SU (2) × U (1) and U ′ (1). Within this model, LFV NSI can be achieved; however, a relation between the diagonal and off-diagonal elements holds As discussed in [6,7], if we generalize the model to include more than one new fermion mixed with ν α and ν β , the Schwartz inequality still implies From the phenomenological point of view, the off-diagonal ǫ f αβ elements are distinguishable from diagonal elements. Ref [10] shows that non-zero ǫ eµ and ǫ eτ provide a better fit to solar neutrino data than the standard ν oscillation scheme. Obviously, this solution does not respect the relation in Eq. (3).
The aim of the present paper is to build a model in which this inequality is violated: In sec II, we present the underlying model for the LFV NSI and discuss the bounds that already exist on the parameters of the model. In sec. III, we summarize our results and suggest strategies to test the model.

II. THE MODEL
In this section, we build an underlying model for effective Lagrangian in Eq. (1). Our model contains a new U (1) gauge boson Z ′ µ which couples both to the matter fields f (f = e, u or d) and to neutrinos as follows and For the low energy-momentum transfer, we can then write In sect. II A, we introduce the underlying scenario that leads to the off-diagonal coupling of (5) in two versions: First in minimal version that violates the lepton number and then in the lepton number conserving case. We then discuss the bounds on the parameters of this sector of the model and briefly comment on the prospects for testing the model in the future. In sect II B, we discuss the model for the interaction of Eq. (4) and the phenomenological consequences of this model.
A. An electroweak invariant model for Z ′ µνα γ µ ν β with α = β and without Z ′ µlα γ µ l β In sec. II A 1, we introduce the minimal version of the model which breaks lepton number and may therefore induce large contribution to neutrinoless double beta decay and/or lepton number violating processes such as B − → π + µ − µ − . In sec. II A 2, we show that with a slight change in the content of the model, a lepton number conserving model can be built without large contribution to neutrinoless double beta decay and other lepton number violating processes. In sec. II A 3, we shall review the experimental bounds on the parameters of model. We will discuss which part of the parameter space gives us large enough g ν .

Lepton number violating version of the model
Let us first introduce two right-handed Weyl fermions which are singlets of the standard model gauge group but under the new gauge U ′ (1) transform as ψ 1 → e iα ψ 1 and ψ 2 → e −iα ψ 2 . (In sec II A 2, we shall discuss another version of the model in which ψ 1 and ψ 2 are promoted to be Dirac fermions.) Having opposite charges, the U ′ (1) gauge anomaly will be canceled. Moreover, we can write mass term for them as where in the right hand side of the equation, we have replaced so N 1 and N 2 have equal masses. The gauge interaction can be then written as Notice that Z ′ converts N 1 and N 2 to each other. If we construct a mechanism that mixes N 1 only with ν α and N 2 only with ν β , we obtain an interaction of type in Eq. (5) without diagonal couplings of formν α γ µ ν α Z ′ µ orν β γ µ ν β Z ′ µ and without corresponding LFV couplings for charged leptons,l β γ µ l α Z ′ µ . In the following, we show that it is possible to build a model for such mixing after spontaneous gauge symmetry breaking. To do so, we shall invoke the famous inverse seesaw mechanism [11] for neutrino mass generation involvingÑ 1L ,Ñ 2L ,Ñ 1R andÑ 2R where (Ñ 1L ,Ñ 1R ) and (Ñ 2L ,Ñ 2R ) form Dirac fermions. In our model, these Weyl fermions are all singlets both under the standard model gauge group and the new U ′ (1) gauge group. In our model, the U ′ (1) is broken by the VEV of a new scalar φ which is a singlet of the SM gauge group and its U ′ (1) charge is equal to that of ψ 1 and opposite to that of ψ 2 ; i.e., under Let us introduce a Z 2 symmetry under which, The parities ofÑ i , L α and L β are shown in table I. The other standard model particles are parity even under this Z 2 . As we mentioned above, our model invokes inverse seesaw mechanism with the following potential which respects Let us definem Assigning lepton number equal to 1 toÑ i , it is easy to show that the potential in Eq. (10) is lepton number conserving and does not induce a mass for light neutrinos. To obtain mass for active neutrinos, the following lepton number violating masses should be added: plus Including these terms, SM neutrinos will obtain a mass proportional to the µ terms. This is the basis of the so-called inverse seesaw mechanism. Notice that µ 21 in Eq (12) breaks Z 2 making it possible to obtain a mixing between the α and β flavors in the neutrino mass matrix. We shall not elaborate further on neutrino masses as the inverse seesaw mechanism is widely studied in the literature.
To obtain a mixing between N i and the neutral leptons, we assign lepton number of −1 to N 1 and N 2 (or equivalently to and introduce the following Z 2 and gauge invariant coupling: The vacuum expectation value of φ, v φ = φ , breaks the U ′ (1) symmetry and therefore induces mass mixing terms as Notice that the M N mass term in Eq. (7) is an explicit source of lepton number violation. Despite this source of lepton number violation, in the limit that the µ terms in Eqs. (11) and (12) vanish, SM neutrinos will remain massless. This can be understood because, in the limit of µ → 0, the symmetric mass matrix for (ν β ,Ñ 2L , cÑ * 2R , cN * 2 ) can be written as whose determinant vanishes so the lightest mass eigenvalue will vanish independent of the values of m 2 , M 2 or M N . The light neutrino will be a linear combination as where β and β ′ are small mixing angles. Similar consideration holds valid for the mass matrix of ν α ,Ñ 1 and N 1 . We can similarly write We then obtain We demand the masses ofÑ i and N i to be larger than ∼ 500 MeV to avoid the bounds from supernova type II cooling and meson decay. This can be achieved if M 2 , M 1 , M N > 500 MeV. The bounds from low energy experiments (such as meson decays) set an upper bound on sin 2 β + sin 2 β ′ and on sin 2 α + sin 2 α ′ . In the limit M 1 ≫ m 1 (M 2 ≫ m 2 ), we find sin α ′ ≫ sin α (sin β ′ ≫ sin β). In order for sin α and sin β to saturate the bounds and therefore to obtain largest possible g ν , we focus on the range that m 2 ∼ M 2 and m 1 ∼ M 1 . Remembering that m 1 and m 2 are given by φ which also contribute to the Z ′ mass, we find The following remarks are in order.
• Notice that the Z 2 symmetry guarantees that N i mixes with only one of ν α or ν β and as a result while we obtain LFV coupling Z ′ µν α γ µ ν β , we do not obtain LFC couplings of Z ′ µν α γ µ ν α and Z ′ µν β γ µ ν β . As a result, the bound |ǫ αβ | 2 ≤ |ǫ αα ǫ ββ | within the model(s) in Ref [6] does not apply here. Without the Z 2 symmetry, N 1 and N 2 could mix simultaneously with ν α and ν β leading to Z ′ µν α γ µ ν α and Z ′ µν β γ µ ν β along with Z ′ µν α γ µ ν β . • We could obtain a mixing between active neutrinos and N i in a more economic version of the model without introducingÑ iL (with a lepton number violating mass of the formÑ T iR cÑ iR ) but in this case the active neutrinos would obtain a mass of order m N sin 2 α ∼ m N sin 2 β which for m N ∼ 500 MeV would imply sin α, sin β < 10 −5 rendering ǫ αβ too small.

Lepton number conserving version of the model
Promoting ψ 1 and ψ 2 to Dirac fermions, we can impose lepton number conservation up to small effects induced by the µ terms in Eqs. (11) and (12). The Z 2 symmetry then implies the mass terms for ψ 1 and ψ 2 to be of form The rest of features of the model will be similar to the model described in sec II A 1.

Bounds on the model parameters
In this section, we discuss the bounds on the mixing of sterile neutrinos with active neutrinos and other parameters of the model. We also comment on the possibility to test the model. For sterile neutrinos lighter than a few MeV, there are strong bounds on the mixing from cosmology [12]. For masses below ∼ 100 MeV, there are strong constraints from the supernova cooling. For masses below that of Kaon (≃ 500 MeV), strong bounds on the mixing come from the Kaon and pion decay. We therefore assume the masses of sterile neutrinos, determined by m N and M i , to be heavier than 500 MeV.
In the literature, strong bounds on the mixing of sterile neutrinos with mass heavier than 500 MeV with active neutrinos have also been reported from NuTeV [13], WA66 [14], CHARM II [15], BELLE [16], Higgs decay [17], NA62 [18], L3 [19], DELPHI [20] and ATLAS+CMS [21]. Ref. [22] gives an updated compilation of the bounds and a forecast for future searches (see also, [23]). A sterile neutrino mixed with ν α can decay into ν α ν γνγ and, if the kinematics allows, into l αlγ ν γ (where γ may or may not correspond to α) via electroweak interaction with a rate suppressed by the square of the mixing. The bounds that we enumerated are all based on searches for the signature of the final charged particles. In our model, there is a possibility of faster two body decay into active neutrinos and Z ′ . If Z ′ decays into a neutrino pair, it will not show up in these experiments so all these bounds can be avoided. Let us formulate the condition for avoiding the bounds. Notice that in our model, in addition to N 1 and N 2 , we havẽ N 1 andÑ 2 which also mix with the active neutrinos. In order to open the decay mode into Z ′ forÑ 1 andÑ 2 (or to be more precise for mass eigenstates composed mainly ofÑ i and N i ), we need a large mixing between them. In sect II A 1, we already showed that m i ∼ M i . If we further impose the condition m i ∼ M i ∼ M N , this condition will be fulfilled and the decay rates of all these sterile neutrinos into Z ′ and active neutrinos will be of the same order. In order for the two body decay into Z ′ ν a to dominate over the electroweak three body decay, we just need g 2 ψ ≫ (g 4 SU(2) /16π 2 )(m N /m W ) 4 which can be readily satisfied. A more challenging requirement is that Z ′ dominantly decays into neutrinos. This requires g 2 ν ≫ g 2 e , for m Z ′ > 2m π g 2 ν ≫ g 2 q , for m Z ′ > 2m µ g 2 ν ≫ g 2 µ and for m Z ′ > 2m τ g 2 ν ≫ g 2 τ . In summary, in order to avoid the enumerated bounds on the mixing of sterile neutrinos with mass above 500 MeV, the decay mode of the sterile neutrinos into Z ′ ν a and then Z ′ → ν aνa must dominate. This in turn requires m N ∼ m i ∼ M i and g 2 ν ≫ g 2 f . In case that the mass of sterile neutrino is of Majorana type, the null results from neutrinoless double beta decay searches set a strong bound on the mixing with ν e , ranging from few×10 −5 to few×10 −2 for 500 MeV < M N < ∼ 10 TeV.
If we are interested in ǫ eτ or ǫ eµ close to the present bounds, we then need to adopt the lepton number conserving version of the model as described in sec. II A 2 to avoid the bounds from 0νββ. Moreover, for the Majorana type sterile neutrinos, there is also a strong bound on the mixing with ν µ of order of 10 −2 from searches for lepton number violating decay mode B − → π + µ − µ − from the LHCb [42].
Since active neutrinos mix with sterile neutrinos, the PMNS 3 × 3 matrix will not be unitary. There are strong bounds on the violation of the unitarity of the PMNS matrix [24]. Some of these bounds do not however apply to our case. Most notably despite the deviation of the PMNS matrix from unitarity, in our model we do not obtain a significant contribution to l − α → l − β γ at one loop level as ν α and ν β mix with different sets of sterile neutrinos. In other words, the heavier mass eigenstate either have a contribution from ν α or from ν β but not from both, making one loop contribution absent. As a result, the bounds on (U † P MN S · U P MN S ) αβ | α =β from l α → l β γ discussed in the literature [24] does not apply here. There will be a two-loop contribution in which both W and Z ′ propagate but the effect will be both GIM and two-loop suppressed and therefore negligible.
As long as the heaviest sterile neutrino mixed with the active neutrinos is much lighter than m Z /2, the bounds on the invisible decay width of Z and those from the leptonic decay modes of W can also be relaxed in our model. In other words, the deviation from standard model prediction for (Z →invisibles) and for (W → l+missing energy) will be suppressed not only by the square of mixing but also by O(m 2 N /m 2 W ). However for m N > m K , K + (π + ) → l + α ν and K + (π + ) → l + β ν will be suppressed by cos 2 α and cos 2 β, respectively. Moreover, the rate of the muon decay which is used to extract G F will be affected. To be on the safe side, we take sin α sin β < 10 −3 to satisfy the bounds from the violation of the unitarity of the PMNS matrix [24]. We then find There are also direct bounds on ( α |(g ν ) eα | 2 ) 1/2 and ( α |(g ν ) µα | 2 ) 1/2 from the K + and π + decays into e + and µ + plus missing energy which are again around 10 −3 [25]. As shown in [26], the DUNE near detector can probe small values of g ν and even determine its flavor structure. As discussed before, the contribution to l α → l β γ in our model is two-loop suppressed but at one loop, we obtain l α → l β Z ′ with a rate estimated as where m 2 ψ /m 2 W comes from the GIM suppression and (m 2 α /m 2 Z ′ ) is the enhancement factor due to the longitudinal component of Z ′ . In PDG [27], there are explicit bounds on such exotic decay modes of τ : Br(τ → eZ ′ ) < 2.7 × 10 −3 and Br(τ → µZ ′ ) < 5 × 10 −3 .
These bounds can be easily satisfied for g ν < 10 −3 and m Z ′ > 10 MeV. The present bound on Br(µ → eZ ′ ) is of order of 10 −5 [28]. Taking (g ν ) µe < 10 −3 , m ψ ∼ GeV and m Z ′ ∼ 10 MeV, we find that Br(µ → eZ ′ ) is of similar order. Future muon decay experiments can easily test this prediction [28]. Similarly, at one loop level and through Z ′ exchange, µ to e conversion can take place with The present bound is set by SINDRUM II collaboration which is 7 × 10 −13 [29] so the present bound can be easily satisfied. The next generation Mu2e and COMET experiments can probe R down to 5 × 10 −17 [30] so they can be sensitive down to (g ν ) eµ g q ∼ 10 −8 . In our model, the Higgs will have new invisible decay modes intoν βÑ2 andν αÑ1 with rates of λ 2 β m H /(4π) and λ 2 α m H /(4π), respectively. Taking λ α,β H =m α ,m β < 2 GeV, the present bound on the branching ratio of the invisibles decay mode, Br(H → invisibles) < 0.2 [27], can be satisfied. Future searches for H → invisibles can test the model. In fact, if Br(H → invisibles) down to O(1%) is measured, the entire parameter space of interest to us can be probed. In principle, we can consider the masses of the sterile neutrinos to be heavier than the Higgs mass but then the mixing parameters sin α and sin β will be suppressed by m 1,2 /m N ∼ (m Z ′ /m N )(Y 1,2 /g ψ ). (Remember that we require relatively light Z ′ to obtain sizable NSI.)

B.
Couplings of matter fields to Z ′ In this section, we discuss two mechanisms for coupling Z ′ to matter fields and then discuss the bounds on the coupling. We also evaluate the maximum values of ǫ f αβ that can be achieved, combining these bounds with the bounds on g ν discussed in sec II A 3. We then propose ideas to test the model.
• Another possibility is a kinetic mixing between Z ′ and the hypercharge gauge where Z ′ µν and B µν are field strength of these gauge bosons. Going to the canonical basis where both gauge boson mass terms and the kinetic terms are diagonal and properly normalized, we find that Z ′ obtains a vector-like coupling to charged fermions which in the leading order in δ is given by As expected, the axial part of the coupling to charged fermions as well as the coupling to neutrinos, which originate from the contribution of Z µ to B µ (rather than that of A µ to B µ ), are further suppressed by a factor of (m 2 Z ′ /m 2 Z )δ sin θ W . Remembering that the Deuteron dissociation process that is invoked by SNO to measure the rate of the neutral current interaction of solar neutrinos is only sensitive to the axial current, no significant bound from this measurement applies to our model. Up to O(δ), N 1 and N 2 will not obtain couplings to the photon so we do not need to worry about the bounds on the millicharged particles. 1 Notice that within this scenario, the effective coupling, ǫ f will be proportional to the electric charge of the matter field, Q f . This means that for neutrinos propagating in neutral media (i.e., all the media such as Earth, Sun, supernova and etc that we know of), there will be no NSI effect on the propagation. However, for the coherent elastic neutrino nucleus scattering experiments such as COHERENT [32] or CONUS [33], the effects of NSI will be present and proportional to the square of the atomic number. The appearance of NSI effects in these experiments and the lack of evidence for them in the neutrino propagation experiments such as DUNE and atmospheric or solar neutrino experiments can be indicative of this particular scenario in which the coupling of the mediator to the matter field originates from the kinetic mixing.
Let us now discuss the experimental bounds. For Z ′ lighter than ∼ 100 MeV, there are strong bounds on g e from beam dump experiments [34] and supernova cooling. As a result, for the case that a e = 0 or in the kinetic mixing scenario, we should take Z ′ to be heavier than ∼ 100 MeV. Still we will have strong bounds from neutrino electron scattering experiments such as BOREXINO, GEMMA and CHARM but these bounds are already encoded in the bounds on NSI parameter, ǫ e αβ . As discussed before, the condition g e ≪ g ν has to be fulfilled to guarantee that the invisible decay mode of Z ′ dominates. Taking (g ν ) αβ = 10 −3 g ψ , m Z ′ = 500 MeVg ψ (see Eq. 16) and g e = 0.1g ν , we find ǫ e αβ = 0.01((g e /g ν )/0.1) which is close to the upper bounds on ǫ e αβ [7]. In the kinetic model, ǫ p αβ = −ǫ e αβ and in the Baryon number gauging model, ǫ q αβ = −ǫ e αβ /(3a e ) which again is close to the present upper bounds [7,[36][37][38][39]. In the scenario where we gauge a linear combination of lepton and Baryon numbers, if we set a e = 0, a wider mass range for Z ′ will open. Z ′ lighter than ∼ 10 MeV is disfavored by cosmology [40]. For 130 MeV > m Z ′ > 10 MeV we can have g q ∼ g τ ∼ g µ ∼ g ν < 10 −3 without violating any bound [41]. Notice that for this mass range even if g q , g µ , g τ > g ν , Z ′ → νν will be still the dominant decay mode. Taking g q = 10 −3 (which saturates the bound from π 0 → γZ ′ [35]), (g ν ) αβ = 10 −3 g ψ and m Z ′ ∼ 500 MeVg ψ (with g ψ > 0.02 to make m Z ′ heavier than 10 MeV and therefore avoid the cosmological bounds), we obtain ǫ q αβ = 0.12 g ψ g q 10 −3 so the values of ǫ q αβ can easily saturate the present bounds on them. In this model, we also obtain diagonal NSI couplings as To have anomaly cancellation we should have a µ + a τ = 3 so a µ and a τ cannot be simultaneously zero unless new chiral doublets are added to the model to cancel the anomalies. We can however have a µ = 0 and a τ = 3 (or a τ = 0 and a µ = 3). In this case, it is obvious that for any α = β, |ǫ q αβ | 2 > |ǫ q αα ǫ q ββ |. More interestingly, since we have a freedom in the choice of a µ /a τ and g q /(g ν ) αβ , we can obtain an arbitrary ratio of ǫ q αβ /ǫ q αα . For example, we can have ǫ q ee = ǫ q µµ = 0, a value of ǫ q eµ saturating the bound and ǫ q τ τ = 0 with an arbitrary value of ǫ q eµ /ǫ q τ τ For heavier Z ′ , visible decay modes Z ′ → µμ and Z ′ → π + π − open so as we discussed in sec II A 3, we should impose g f ≪ g ν . Like the case of nonzero g e , we can still obtain large ǫ q αβ with ǫ q αβ ≫ ǫ q αα , ǫ q ββ . Notice that in our model, all the new particles, including the new scalar φ, have a mass above ∼ 10 MeV and decay away well before the neutrino decoupling in the early universe so the model is not constrained by the bound on the extra relativistic degrees of freedom from cosmological constraints.

III. SUMMARY AND OUTLOOK
We have built a SU (2) × U (1) invariant UV complete model which leads to NSI with LFV couplings ǫ f αβ with α = β and arbitrary ratios of ǫ f αβ /ǫ f αα and ǫ f αβ /ǫ f ββ . The effects of such NSI can show up in various neutrino oscillation experiments and in different coherent elastic neutrino nucleus scattering experiments. The model is based on a pair of fermions which are singlets of the standard model but have opposite charges under a new U ′ (1) gauge symmetry with a relatively light gauge boson of mass 10 MeV < m Z ′ < few 100 MeV. The lower bounds on m Z ′ is imposed by BBN and CMB [40]. After electroweak and U ′ (1) symmetry breaking, the new fermions mix with the ordinary neutrinos of different flavors leading to LFV couplings that we require. The model incorporates the inverse seesaw mechanism which helps to decouple the scales of the masses of the new sterile neutrinos from the masses of active neutrinos. The values of ǫ f αβ are given by the mixing between the sterile neutrinos and the active ones. We have presented the model in two versions: lepton number violating version and lepton number conserving one. We have also proposed two different mechanisms to couple Z ′ to the matter fields based on gauging a linear combination of the Baryon and lepton numbers and based on a kinetic mixing between the field strengths of the new gauge boson and that of the hypercharge. We have discussed how we can distinguish them by combining the results from the neutrino oscillation experiments and from the coherent elastic neutrino nucleus scattering experiments.
We have taken the masses of the sterile neutrinos to be heavier than 500 MeV (i.e., m K + ) to avoid the strong bounds on their mixing from the Kaon and pion decay. In our model, the sterile neutrinos promptly decay into Z ′ and SM neutrinos and appear as missing energy in colliders so the bounds that exist on the mixing of the sterile neutrinos with mass above 0.5 GeV from searching for their charged decay products do not apply here. On the other hand, we discuss that since the active sterile mixing, sin α and sin β are expected to be of order of the ratio of m Z ′ to the masses of sterile neutrinos, large ǫ αβ (being proportional to sin α sin β/m 2 Z ′ ) implies that the sterile neutrino masses are not much heavier than a few GeV. The invisible decay of the Higgs also imposes a similar bound. The masses of sterile neutrinos should be in the range between 500 MeV to a few GeV. The entire range can be probed by improving the bounds on the invisible decay rate of the Higgs. In principle, it can also be probed by studying the neutral current neutrino nucleus scattering experiments in which the neutrino beam is energetic enough to produce the sterile neutrinos.
In the case of ǫ eµ and ǫ eτ , we require a mixing between the sterile neutrinos with ν e . To avoid the strong bounds from neutrinoless double beta decay searches on the mixing of a Majorana sterile neutrino with ν e , we must adopt the lepton number conserving version of the model. In the case of ǫ τ µ , it is also favorable to adopt the lepton number conserving version to avoid the constraint from searches for B + → π + µ − µ − by the LHCb [42].
The couplings of the active neutrinos to Z ′ , g ν , can be as large as 10 −3 constrained by the bounds on K + → µ + Z ′ ν [25]. As shown in [26] smaller values of g ν can be tested by the near detector of DUNE. On the other hand, the couplings to quarks can also be as large as 10 −3 constrained by π 0 → Z ′ γ for m Z ′ < m π 0 [35]. Within this part of the parameter space of our model, we can obtain values of ǫ q αβ saturating the present bounds on it. We have shown that despite LFV, at one loop level, there is no contribution to l − α → l − β γ so no significant bound comes even from the very stringent bound on Br(µ − → e − γ) on our model parameter space. However, at one loop level, there is a contribution to µ − to e − conversion on nuclei through virtual Z ′ exchange which can be probed by COMET and mu2e experiments in future. Moreover, for m Z ′ < m µ at one loop level, we obtain µ − → e − Z ′ whose effect may be probed by studying the energy spectrum of e − emitted in the muon decay. A signal for µ − − e − conversion in the experiments without a signal for µ − → e − γ in the future searches combined with the traces of Z ′ in the muon and meson decays along with observable ǫ f eµ and a deviation of the Higgs invisible decay rate from the SM prediction can be considered smoking gun signatures for the present model.