Non-Standard Interactions in Radiative Neutrino Mass Models

Models of radiative Majorana neutrino masses require new scalars and/or fermions to induce lepton number violating interactions. We show that these new particles also generate observable neutrino nonstandard interactions (NSI) with matter. We classify radiative models as type-I or II, with type-I models containing at least one Standard Model (SM) particle inside the loop diagram generating neutrino mass, and type-II models having no SM particle inside the loop. While type-II radiative models do not generate tree-level NSI, popular models which fall under the type-I category are shown, somewhat surprisingly, to generate observable NSI at tree-level, while being consistent with direct and indirect constraints from colliders, electroweak precision data and charged-lepton flavor violation (cLFV). We survey such models where neutrino masses arise at one, two and three loops. In the prototypical Zee model which generates neutrino masses via one-loop diagrams involving charged scalars, we find that diagonal NSI can be as large as ($8\%, 3.8\%, 9.3\%$) for ($\varepsilon_{ee},\varepsilon_{\mu \mu}, \varepsilon_{\tau \tau}$), while off-diagonal NSI can be at most ($10^{-3}\%, 0.56\%, 0.34\%$) for ($\varepsilon_{e\mu},\varepsilon_{e \tau}, \varepsilon_{\mu \tau}$). In one-loop neutrino mass models using leptoquarks (LQs), $(\varepsilon_{\mu\mu},\, \varepsilon_{\tau\tau})$ can be as large as $(21.6\%,\,51.7\%)$, while $\varepsilon_{ee}$ and $(\varepsilon_{e\mu},\, \varepsilon_{e\tau},\,\varepsilon_{\mu\tau})$ can at most be 0.6\%. The most stringent constraints on the diagonal NSI are found to come from neutrino oscillation and scattering experiments, while the off-diagonal NSI are mostly constrained by low-energy processes, such as atomic parity violation and cLFV. While our analysis is focused on radiative neutrino mass models, it essentially covers all NSI possibilities with heavy mediators.


Introduction
The origin of tiny neutrino masses needed to explain the observed neutrino oscillation data is of fundamental importance in particle physics. Most attempts that explain the smallness of these masses assume the neutrinos to be Majorana particles, in which case their masses could arise from effective higher dimensional operators, suppressed by a high energy scale that characterizes lepton number violation. This is the case with the seesaw mechanism, where the dimension-five operator [1] suppressed by an inverse mass scale Λ is induced by integrating out Standard Model (SM) singlet fermions [2][3][4][5][6], SU (2) L triplet scalars [7][8][9][10], or SU (2) L triplet fermions [11] with mass of order Λ. 1 In Eq. (1.1), L stands for the lepton doublet, and H for the Higgs doublet, with i, j, k, l denoting SU (2) L indices, and ik is the SU (2) L antisymmetric tensor. Once the vacuum expectation value (VEV) of the Higgs field, H 0 246 GeV is inserted in Eq. (1.1), Majorana masses for the neutrinos given by m ν = v 2 /Λ will be induced. For light neutrino masses in the observed range, m ν ∼ (10 −3 − 10 −1 ) eV, the scale Λ should be around 10 14 GeV. The mass of the new particle that is integrated out need not be Λ, since it is parametrically different, involving a combination of Yukawa couplings and Λ. For example, in the type-I seesaw model the heavy right-handed neutrino mass goes as M R ∼ y 2 D Λ, which can be near the TeV scale, if the Dirac Yukawa coupling y D ∼ 10 −6 . However, it is also possible that y D ∼ O(1), in which case the new physics involved in neutrino mass generation could not be probed directly in experiments. 2 An alternative explanation for small neutrino masses is that they arise only as quantum corrections [14][15][16] (for a review, see Ref. [17]). In these radiative neutrino mass models, the tree-level Lagrangian does not generate O 1 of Eq. (1.1), owing to the particle content or symmetries present in the model. If such a model has lepton number violation, then small Majorana masses for neutrinos will be induced at the loop level. The leading diagram may arise at one, two, or three loop level, depending on the model details, which will have an appropriate loop suppression factor, and typically a chiral suppression factor involving a light fermion mass as well. 3 For example, in the two-loop neutrino mass model of Refs. [15,16], small and calculable m ν arises from the diagram shown in Fig. 43, which is estimated to be of order , (1.2) assuming normal ordering of neutrino masses and requiring large µ − τ mixing. Here f, h are Yukawa couplings involving new charged scalars with mass of order M . Even with f ∼ h ∼ 1, to obtain m ν ∼ 0.1 eV, one would require the scalar mass M ∼ TeV. This type of new physics can be directly probed at colliders, enabling direct tests of the origin of neutrino mass. When the mediators of neutrino mass generation have masses around or below the TeV scale, they can also induce other non-standard processes. The focus of this paper is neutrino non-standard interactions (NSI) [18] induced by these mediators. These NSI are of great phenomenological interest, as their presence would modify the standard threeneutrino oscillation picture. The NSI will modify scattering experiments, as the production and detection vertices are corrected; they would also modify neutrino oscillations, primarily through new contributions to matter effects. There have been a variety of phenomenological studies of NSI in the context of oscillations, but relatively lesser effort has gone into the ultraviolet (UV) completion of models that yield such NSI (for a recent update, see Ref. [19]). A major challenge in generating observable NSI in any UV-complete model is that there are severe constraints arising from charged lepton flavor violation (cLFV) [20]. One possible way to avoid such constraints is to have light mediators for NSI [21][22][23]. In contrast to these attempts, in this paper we focus on heavy mediators, and study the range of NSI allowed in a class of radiative neutrino mass models. 4 Apart from being consistent with cLFV constraints, these models should also be consistent with direct collider searches for new particles and precision electroweak constraints. We find, somewhat surprisingly, that the strengths of the diagonal NSI can be (20-50)% of the weak interaction strength for the flavor diagonal components in a class of popular models that we term as type-I radiative neutrino mass models, while they are absent at tree-level in another class, termed type-II radiative models.

Type-I and type-II radiative neutrino mass models
We propose a nomenclature that greatly helps the classification of various radiative models of neutrino mass generation. One class of models can be described by lepton number violating effective higher dimensional operators, similar to Eq. (1.1). A prototypical example is the Zee model [14] which introduces a second Higgs doublet and a charged SU (2) L -singlet scalar to the SM. Interactions of these fields violate lepton number, and would lead to the effective lepton number violating (∆L = 2) dimension 7 operator with indices i, j, .. referring to SU (2) L , and e c standing for the SU (2) L singlet let-handed positron state. Neutrino masses arise via the one-loop diagram shown in Fig. 4. The induced neutrino mass has an explicit chiral suppression factor, proportional to the charged lepton mass inside the loop. Operator O 2 can be obtained by cutting the diagram of Fig. 4. We call radiative neutrino mass models of this type, having a loop suppression and a chirality suppression proportional to a light charged fermion mass, and expressible in terms of an effective higher dimensional operator as in Eq. (1.3) as type-I radiative models. A classification of low dimensional operators that violate lepton number by two units has been worked out in Ref. [26]. Each of these operators can generate a finite set of type-I radiative neutrino mass models in a well-defined manner. Lepton number violating phenomenology of these operators has been studied in Ref. [27].
Another well known example in this category is the two-loop neutrino mass model of Refs. [15,16], which induces an effective d = 9 operator Neutrino masses arise in this model via the two-loop diagrams shown in Fig. 43, which has a chiral suppression factor proportional to m 2 , with standing for the charged leptons of the SM. This category of type-I radiative neutrino mass models is populated by one-loop, twoloop, and three-loop models. Popular one-loop type-I models include the Zee model [14] (cf. Sec. 4), and its variant with LQs replacing the charged scalars (cf. Sec. 5). This variant is realized in supersymmetric models with R-parity violation [28]. Other one-loop models include SU (2) L -triplet LQ models (cf. Sec. 7.1.6) wherein the neutrino mass is proportional to the up-type quark masses [29,30]. Ref. [31] has classified simple realizations of all models leading to d = 7 lepton number violating operators, which is summarized in Sec. 2. Popular type-I two-loop models include the Zee-Babu model [15,16] (cf. Sec. 7.2.1), a variant of it using LQs and a diquark [32] (cf. Sec. 7.2.2), a pure LQ extension [33] (cf. Sec. 7.2.3), a model with LQs and vector-like fermions [34] (cf. Sec. 7.2.4), and the Angelic model [35] (cf. Sec. 7.2.5). We also present here a new two-loop model (cf. Sec. 7.2.9) with LQs wherein the neutrino masses arise proportional to the up-type quark masses. Type-I three-loop models include the KNT model [36] (cf. Sec. 7.3.1), an LQ variant of the KNT model [37] (cf. Sec. 7.3.4), the AKS model [38] (cf. Sec. 7.3.2), and the cocktail model [39] (cf. Sec. 7.3.3). For a review of this class of models, see Ref. [17].
A systematic approach to identify type-I radiative models is to start from a given ∆L = 2 effective operators of the type O 2 of Eq. (1.3), open the operator in all possible ways, and identify the mediators that would be needed to generate the operator. Such a study was initiated in Ref. [26], and further developed in Refs. [31,40]. We shall rely on these techniques. In particular, the many models suggested in Ref. [31] have been elaborated on in Sec. 7, and their implications for NSI have been identified. This method has been applied to uncover new models in Ref. [41].
In all these models there are new scalar bosons, which are almost always necessary for neutrino mass generation in type-I radiative models using effective higher dimensional operators. For future reference, we list in Table I all possible new scalar mediators in type-I radiative models that can couple to neutrinos, along with their SU (3) c × SU (2) L × U (1) Y quantum numbers, field components and electric charges (in superscript), and corresponding
Lagragian terms responsible for NSI. We will discuss them in detail in 4, 5 and 7. The models discussed in Sec. 7 contain other particles as well, which are however not relevant for the NSI discussion, so are not shown in Table I. Note that the scalar triplet ∆(1, 3, 1) could induce neutrino mass at tree-level via type-II seesaw mechanism [7][8][9][10], which makes radiative models involving ∆ field somewhat unattractive, and therefore, is not included in our subsequent discussion. There is one exception to the need for having new scalars for type-I radiative models (see Sec. 7.1.1). The Higgs boson and the W, Z bosons of the SM can be the mediators for radiative neutrino mass generation, with the new particles being fermions. In this case, however, there would be tree-level neutrino mass á la type-I seesaw mechanism [2][3][4][5][6], which should be suppressed by some mechanism or symmetry. Such a model has been analyzed in Refs. [42,43], which leads to interesting phenomenology, see Sec. 7.1.
From the perspective of neutrino NSI, these type-I radiative models are the most interesting, as the neutrino couples to a SM fermion and a new scalar directly, with the scalar mass near the TeV scale. We have analyzed the ranges of NSI possible in all these type-I radiative models. Our results are summarized in Fig. 58 and Table XX. A second class of radiative neutrino mass models has entirely new (i.e., non-SM) particles inside the loop diagrams generating the mass. These models cannot be derived from effective ∆L = 2 higher-dimensional operators, as there is no way to cut the loop diagram and generate such operators. We term this class of models type-II radiative neutrino mass models (cf. Sec. 8). The induced neutrino mass may have a chiral suppression, but this is not proportional to any light fermion mass. Effectively, these models generate operator O 1 of Eq. (1.1), but with some loop suppression. From a purely neutrino mass perspective, the scale of new physics could be of order 10 10 GeV in these models. However, there are often other considerations which make the scale near a TeV, a prime example being the identification of a WIMP dark matter with a particle that circulates in the loop diagram generating neutrino mass.
A well-known example of the type-II radiative neutrino mass model is the scotogenic model [44] which assumes a second Higgs doublet and right-handed neutrinos N beyond the SM. A discrete Z 2 symmetry is assumed under which N and the second Higgs doublet are odd. If this Z 2 remains unbroken, the lightest of the Z 2 -odd particles can serve as a dark matter candidate. Neutrino mass arises through the diagram of Fig. 56. Note that this diagram cannot be cut in any way to generate an effective higher dimensional operator of the SM. While the neutrino mass is chirally suppressed by M N , this need not be small, except for the desire for it (or the neutral component of the scalar) to be TeV-scale dark matter. There are a variety of other models that fall into the type-II category [45][46][47][48][49][50].
The type-II radiative neutrino mass models will have negligible neutrino NSI, as the neutrino always couples to non-SM fermions and scalars. Any NSI would be induced at the loop level, which would be too small to be observable in experiments. As a result, in our comprehensive analysis of radiative neutrino mass models for NSI, we can safely ignore type-II models.
One remark is warranted here. Consider an effective operator of the type Such an operator would lead to neutrino masses at the two-loop level, as can be seen in an explicit model shown in Fig. 57. Although this model can be described as arising from an effective ∆L = 2 operator, the neutrino mass has no chiral suppression here. The mass scale of the new scalars could be as large as 10 10 GeV. Such models do belong to type-I radiative models; however, they are more like type-II models due to the lack of a chiral suppression. In any case, the NSI induced by the LQs that go inside the loop diagram for neutrino masses is already covered in other type-I radiative models that we have analyzed. Another example of this type of operator is L i L j H k H l ik jl (H † H), which is realized for instance in the model of Ref. [43] (see Sec. 7.1.1). Such effective operators, which appear as products of lower operators, were treated as trivial in the classification of Ref. [26].

Summary of results
We have mapped out in this paper the allowed ranges for the neutrino NSI parameters ε αβ in radiative neutrino mass models. We present a detailed analysis of the Zee model [14] with light charged scalar bosons. To map out the allowed values of ε αβ , we have analyzed constraints arising from the following experimental and theoretical considerations: i) Contact interaction limits from LEP (cf. Sec. 4.6); ii) Monophoton constraints from LEP (cf. Sec. 4.11); iii) Direct searches for charged scalar pair and single production at LEP (cf. Sec. 4.7.1); iv) Pair production of charged scalars at LHC (cf. Sec. 4.7.2); v) Higgs physics constraints from LHC (cf. Sec. 4.10); vi) Lepton universality in W ± decays (cf. Sec Imposing these constraints, we find that light charged scalars, arising either from the SU (2) L -singlet or doublet field or an admixture, can have a mass near 100 GeV. Neutrino NSI obtained from the pure SU (2) L -singlet component turns out to be unobservably small. However, the SU (2) L -doublet component in the light scalar can have significant Yukawa couplings to the electron and the neutrinos, thus inducing potentially large NSI. The maximum allowed NSI in this model is summarized below: These values are significantly larger than the ones obtained in Ref. [51], where the contributions from the doublet Yukawa couplings of the light charged Higgs were ignored.
We have also analyzed in detail leptoquark (LQ) models of radiative neutrino mass generation. As the base model we analyze the LQ version of the Zee model, the results of which can also be applied to other LQ models with minimal modification. This analysis took into account the following experimental constraints: i) Direct searches for LQ pair and single production at LHC (cf. Sec. 5.3); ii) APV (cf. Sec. 5.1.1); iii) Charged lepton flavor violation (cf. Secs. 5.1.4 and 5.1.5); and iv) rare meson decays (cf. Sec. 5.1.6). Including all these constraints we found the maximum possible NSI induced by the singlet and doublet LQ components, as given below: Our results yield somewhat larger NSI compared to the results of Ref. [52] which analyzed, in part, effective interactions obtained by integrating out the LQ fields. We also analyzed a variant of the LQ model with SU (2) L -triplet LQs, which have couplings to both up and down quarks simultaneously. The maximum NSI in this case are found to be as follows: For completeness, we also list here the maximum possible NSI in the two-loop Zee-Babu model: The NSI predictions in all other models analyzed here will fall into one of the above categories.
The rest of the paper is structured as follows. In Sec. 2, we discuss the classification of low-dimensional lepton-number violating operators and their UV completions. In Sec. 3, we briefly review neutrino NSI and establish our notation. Sec. 4 discusses the Zee model of neutrino masses and derives the various experimental and theoretical constraints on the model. Applying these constraints, we derive the allowed range for the NSI parameters. Here we also show how neutrino oscillation data may be consistently explained with large NSI. In Sec. 5 we turn to the one-loop radiative model for neutrino mass with LQs. Here we delineate the collider and low energy constraints on the model and derive the ranges for neutrino NSI. In Sec. 6, we discuss a variant of the one-loop LQ model with triplet LQ. In Sec. 7 we discuss other type-I models of radiative neutrino mass and obtain the allowed values of ε αβ . We briefly discuss NSI in type-II models in Sec. 8. In Sec. 9 we conclude. Our results are tabulated in Table XX and summarized in Fig. 58.

Classification of ∆L = operators and their UV completions
It is instructive to write down low-dimensional effective operators that carry lepton number of two units (∆L = 2), since all type-I radiative models can be constructed systematically from these operators. Here we present a summary of such operators through d = 7 [26]. We use two component Weyl notation for SM fermions and denote them as The Higgs field of the SM is denoted as H 1, 2, 1 2 . The ∆L = 2 operators in the SM are all odd-dimensional. The full list of operators through d = 7 is given by [26]: Not listed here are products of lower-dimensional operators, such as O 1 × HH, with the SU (2) L contraction of HH being a singlet. Here O 1 is the Weinberg operator [1], while the remaining operators are all d = 7. 5 In this paper, we shall analyze all models of neutrino mass arising from these d = 7 operators for their NSI, as well as the two-loop Zee-Babu model arising from O 9 of Eq. (1.4). A few other models that have been proposed in the literature with higher dimensional operators will also be studied. The full list of d = 9 models is expected to contain a large number, which has not been done to date. Each of these d = 7 operators can lead to finite number of UV complete neutrino mass models. The generic diagrams that induce all of the d = 7 operators are shown in Fig. 1. Take for example the operator O 2 in Eq. (2.2b). There are two classes of models that can generate this operator with the respective mediators obtained from the following contractions (see Table II): In the naming convention of Ref. [26], operators were organized based on how many fermion fields are in them. Operators O5 − O7, which are d = 9 operators, appeared ahead of the d = 7 operator O8.   [31]. Here φ and η generically denote scalars and ψ is a generic vectorlike fermion. The SM quantum numbers of these new fields are as indicated.
Here the pairing of fields suggests the mediator necessary. The (LL) contraction would require a scalar that can be either an SU (2) L singlet, or a triplet. The (e c H) contraction would require a new fermion, which is typically a vectorlike fermion. Thus, O 1 2 has two UV completions, with the addition of a vectorlike lepton ψ 1, 2, − 3 2 to the SM, along with a scalar which is either a singlet η + (1, 1, 1), or a triplet ∆(1, 3, 1). The choice of ∆(1, 3, 1) can lead to the generation of the lower d = 5 operator at tree level via type-II seesaw, and therefore, is usually not employed in radiative models. The model realizing O 1 2 with ψ 1, 2, − 3 2 vectorlike lepton and η + (1, 1, 1) scalar is discussed in Sec. 7.1.2. Similarly operator O 2 2 has a unique UV completion, with two scalars added to the SM -one η + (1, 1, 1) and one Φ 1, 2, 1 2 . This is the Zee model of neutrino mass, discussed at length in Sec. 4. Operators O 3a and O 3b in Eq. (2.2c) can be realized by the UV complete models given in Table. III [31]. Here all possible contraction among the fields are shown, along with the required mediators to achieve these contractions. Fields denoted as φ and η are scalars, while ψ is a vectorlike fermion. The SM quantum numbers for each field are also indicated in the Table. We shall analyze neutrino NSI arising from each of these models in Sec. 7.
The UV completions of operators O 4 and O 8 are shown in Tables IV and V respectively [31]. These models will all be analyzed in Sec. 7 for neutrino NSI. Note that in both O 4 and O 8 , pairing of un-barred and barred fermion fields is not included, as the mediators for such an UV completion will have to be vector bosons which would make such models difficult to realize. As a result, only O 4b can be realized with scalar and fermionic exchange.  [31]. Here the models in the top segment require a new scalar φ and a vectorlike fermion ψ, while those in the lower segment require two scalar fields φ and η.

Classification based on topology of diagrams
Rather than classifying radiative neutrino mass models in terms of effective ∆L = 2 operators, one could also organize them in terms of the topology of the loop diagrams [12,53,54].  Possible one-loop topologies are shown in Fig. 2 [12,53], and the two-loop topologies are shown in Fig. 3 [54]. Note that in the two-loop diagrams, two Higgs particles that are connected to two internal lines in possible ways are not shown. For the purpose of NSI, we find the classification based on type-I and type-II suggested here most convenient. The classification based on the diagram topology does not specify whether the internal particles are SM fermions or not, and the NSI effects arise only when neutrino couples to the SM fermions. Let us also note that the first diagram of Fig. 2 and the first two diagrams of Fig. 3 are the ones that appear most frequently in the explicit type-I radiative models that we discuss in subsequent sections.

Neutrino non-standard interactions
Neutrino NSI can be of two types: Neutral Current (NC) and Charged Current (CC). The CC NSI of neutrinos with the matter fields in general affects the production and detection of neutrinos, while the NC NSI affects the neutrino propagation in matter. In the lowenergy regime, neutrino NSI with matter fields can be formulated in terms of an effective four-fermion Lagrangian as follows [18]: where G F is the Fermi coupling constant, and P X (with X = L, R) denotes the chirality projection operators P L,R = (1 ∓ γ 5 )/2. These projection operators can also be reparameterized into vector and axial components of the interaction. The dimensionless coefficients ε αβ are the NSI parameters that quantify the strength of the NSI between neutrinos of flavors α and β and the matter field f ∈ {e, u, d} (for NC) and f = f ∈ {u, d} (for CC). If ε αβ = 0 for α = β, the NSI violates lepton flavor, while for ε αα = ε ββ , it violates lepton flavor universality. The vector component of NSIs, ε f V αβ = ε f L αβ +ε f R αβ , affects neutrino oscillations by providing a new flavor-dependent matter effect. 6 The effective Hamiltonian for the matter effect is given by where U PMNS is the standard 3 × 3 lepton mixing matrix, E is the neutrino energy, N e (x) is the electron number density as a function of the distance x traveled by the neutrino in matter, and the 1 in the 1 + ε ee term is due to the standard CC matter potential. The Hamiltonian level NSI in Eq. (3.3) is related to the Lagrangian level NSI in Eq. (3.1) as follows: 4) where N f (x) is the number density of fermion f at position x, and N p(n) /N e is the average ratio of the density of protons (neutrons) to the density of electrons along the neutrino propagation path. Note that the coherent forward scattering of neutrinos with nucleons can be thought of as the incoherent sum of the neutrino scattering with the constituent quarks, because the nucleon form factors are equal to one in the limit of zero momentum transfer. Assuming electric charge neutrality of the medium, we can set N p (x)/N e (x) = 1 and define the ratio Y n (x) ≡ N n (x)/N e (x) to rewrite Eq. (3.4) as In the Earth, the ratio Y n which characterizes the matter chemical composition can be taken to be constant to very good approximation. According to the Preliminary Reference Earth Model (PREM) [55], Y n = 1.012 in the mantle and 1.137 in the core, with an average value Y n = 1.051 all over the Earth. On the other hand, for solar neutrinos, Y n (x) depends on the distance to the center of the Sun and drops from about 1/2 in the center to about 1/6 at the border of the solar core [56,57]. In the following sections, we will derive the predictions for the NSI parameters ε αβ in various radiative neutrino mass models, which should then be compared with the experimental and/or global fit constraints [58][59][60][61] on ε αβ using Eq. (3.5). We would like to emphasize two points in this connection: (i) Depending on the model, we might have NSI induced only in the neutrino-electron or neutrino-nucleon interactions, or involving only left-or right-chirality of the matter fields. In such cases, only the relevant terms in Eq. (3.5) should be considered, while comparing with the experimental or global fit constraints.
(ii) Most of the experimental constraints [59] are derived assuming only one NSI parameter at a time, whereas within the framework of a given model, there might exist some non-trivial correlation between NSI involving different neutrino flavors, as we will see below. On the other hand, the global fits [60,61] usually perform a scan over all NSI parameters switched on at the same time in their analyses, whereas for a given model, the cLFV constraints usually force the NSI involving some flavor combinations to be small, in order to allow for those involving some other flavor combination to be sizable. To make a conservative comparison with our model predictions, we will quote the most stringent values from the set of experimental and global fit constraints, as well as the future DUNE sensitivities [62][63][64][65]; see Tables IX and XVII. 4 Observable NSI in the Zee model One of the simplest extensions of the SM that can generate neutrino mass radiatively is the Zee Model [14], wherein small Majorana masses arise through one-loop diagrams. This is a type-I radiative model, as it can be realized by opening up the ∆L = 2 effective d = 7 operator O 2 = L i L j L k e c H l ij kl , and since the induced neutrino mass has a chiral suppression factor proportional to the charged lepton mass. Due to the loop and the chiral suppression factors, the new physics scale responsible for neutrino mass can be at the TeV scale. The model belongs to the classification O 2 2 of Table II. The model assumes the SM gauge symmetry SU (3) c × SU (2) L × U (1) Y , with an extended scalar sector. Two Higgs doublets Φ 1,2 (1, 2, 1/2), and a charged scalar singlet η + (1, 1, 1) are introduced to facilitate lepton number violating interactions and thus nonzero neutrino mass. The leptonic Yukawa Lagrangian of the model is given by: where {α, β} are generation indices, {i, j} are SU (2) L indices, Φ a ≡ iτ 2 Φ a (a = 1, 2) and c denotes the left-handed antilepton fields.
Here and in what follows, a transposition and charge conjugation between two fermion fields is to be understood. Note that due to Fermi statistics, f αβ = −f βα . Expanding the first term of the Lagrangian Eq. (4.1) leads to the following couplings of η + : The presence of two Higgs doublets Φ 1,2 allows for a cubic coupling in the Higgs potential,

Scalar sector
We can start with a general basis, where both Φ 1 and Φ 2 acquire vacuum expectation values (VEVs): (4.4) However, without loss of generality, we can choose to work in the Higgs basis [68] where only one of the doublet fields gets a VEV v given by In the CP-conserving limit where Im(λ 5,6 ) = 0, the CP-odd state will decouple from the CP-even states. One can then rotate the CP-even states into a physical basis {h, H} which would have masses given by [69]: whereas the CP-odd scalar mass is given by The mixing angle between the CP-even eigenstates {H 0 1 , H 0 2 }, defined as is given by We will identify the lightest CP-even eigenstate h as the observed 125 GeV SM-like Higgs and use the LHC Higgs data to obtain constraints on the heavy Higgs sector (see Sec. 4.10). We will work in the alignment/decoupling limit, where β − α → 0 [70][71][72][73], as suggested by the LHC Higgs data [74,75].

Neutrino mass
In the Higgs basis where only the neutral component of H 1 gets a VEV, the Yukawa interaction terms in Eq. (4.1) of fermions with the scalar doublets H 1 and H 2 become where Y and Y are the redefined couplings in terms of the original Yukawa couplings y 1 and y 2 given in Eq. (4.1) and where H a = iτ 2 H a (a = 1, 2) with τ 2 being the second Pauli matrix. After electroweak symmetry breaking, the charged lepton mass matrix reads as (4.20) Without loss of generality, one can work in a basis where M is diagonal, i.e., M = diag (m e , m µ , m τ ). The Yukawa coupling matrix f involving the η + field in Eq. (4.1) is taken to be defined in this basis. The Yukawa couplings in Eq. (4.1), together with the trilinear term in the scalar potential Eq. (4.3), generate neutrino mass at the one-loop level, as shown in Fig. 4. Here the dot (•) on the SM fermion line indicates mass insertion due to the SM Higgs VEV. There is a second diagram obtained by reversing the arrows on the internal particles. Thus, we have a symmetric neutrino mass matrix given by where κ is the one-loop factor given by with ϕ given in Eq. (4.13). The matrix f that couples the left-handed lepton doublets to the charged scalar η + can be made real by a phase redefinitionP fP , whereP is a diagonal phase matrix, while the Yukawa coupling Y in Eq. (4.19) is in general a complex asymmetric matrix: Here the matrix Y is multiplied by (ν e ,ν µ ,ν τ ) from the left and (e R , µ R , τ R ) T from the right, in the interaction with the charged scalar H + . Thus the neutrino NSI will be governed by the matrix elements (Y ee , Y µe , Y τ e ), which parametrize the couplings of ν α with electrons in matter. Note that in the limit Y ∝ M l , as was suggested by Wolfenstein [76] by imposing a discrete Z 2 symmetry to forbid the tree-level flavor changing neutral currents (FCNC) mediated by the neutral Higgs bosons, the diagonal elements of M ν would vanish, yielding neutrino mixing angles that are not compatible with observations [77,78]. For a variant of the Zee-Wolfenstein model with a family-dependent Z 4 symmetry which is consistent with neutrino oscillation data, see Ref. [79].
From Eq. (4.21) it is clear that only the product of the Yukawa couplings f and Y is constrained by the neutrino oscillation data. Therefore, by taking Y ∼ O(1) and f 1 in the neutrino mass matrix Eq. (4.21), we can correctly reproduce the neutrino oscillation parameters (see Sec. 4.13). This choice maximizes the neutrino NSI in the model. We shall adopt this choice.
Since the model has two Higgs doublets, in general both doublets will couple to up and down quarks. If some of the leptonic Yukawa couplings Y αe of Eq. (4.23) are of order unity, so that significant neutrino NSI can be generated, then the quark Yukawa couplings of the second Higgs doublet H 2 will have to be small. Otherwise chirality enhanced meson decays, such as π + → e + ν will occur with unacceptably large rates. Therefore, we assume that the second Higgs doublet H 2 is leptophilic in our analysis.

Charge breaking minima
To have sizable NSI, we need a large mixing ϕ between the singlet and doublet charged scalar fields η + and H + 2 . From Eq. (4.13), this means that we need a large trilinear µ-term. But µ cannot be arbitrarily large, as it leads to charge breaking minima (CBM) of the potential [66,67]. We numerically analyze the scalar potential given by Eq. (4.7) to ensure that it does not develop any CBM deeper than the charge-conserving minimum (CCM).
We take µ 2 2 , µ 2 η > 0. The field H 1 is identified approximately as the SM Higgs doublet, and therefore, the value of λ 1 is fixed by the Higgs mass (cf. Eq. (4.8)), and the corresponding mass-squared term is chosen to be negative to facilitate electroweak symmetry breaking (µ 2 1 > 0 in Eq. (4.7)). Note that the cubic scalar coupling µ can be made real as any phase in it can be absorbed in η − by a field redefinition.
In order to calculate the most general minima of the potential, we assign the following VEVs to the scalar fields: where v η and v 1 can be made real and positive by SU (2) L × U (1) Y rotations. A nonvanishing VEV v η would break electric charge conservation, as does a nonzero value of sin γ. Thus, we must ensure that the CBM of the potential lie above the CCM. The Higgs potential, after inserting Eq. (4.24) in Eq. (4.7), reads as 2 3 cos δ + λ 6 v 2 1 cos (θ 2 + δ ) + λ 7 v 3 2 cos (θ 3 + δ ) + λ 10 v 2 η cos (θ 4 + δ )] +v 1 v 2 cos γ 2 [λ 4 + λ 5 cos (θ 1 + 2δ )] − 2µv η cos δ sin γ}. (4.25) Here θ 1 , θ 2 , θ 3 , and θ 4 are respectively the phases of the quartic couplings λ 5 , λ 6 , λ 7 , and λ 10 . For simplicity, we choose these quartic couplings, as well as λ 9 to be small. This choice does not lead to any run-away behavior of the potential. We keep all diagonal quartic couplings to be nonzero, so that the potential remains bounded. (All boundedness conditions are satisfied if we choose, as we do for the most part, all the quartic couplings to be positive.) We also keep the off-diagonal couplings λ 3 and λ 8 nonzero, as these couplings help in satisfying constraints from the SM Higgs boson properties from the LHC. Eq. (4.25) yields five minimization conditions from which {v 1 , v 2 , v η , δ, γ} can be solved numerically for any given set of masses and quartic couplings. The mass parameters are derived from the physical masses of h + , H + and h in the CCM. We vary m h + from 50 to 500 GeV and choose three benchmark points for m H + : {0.7, 1.6, 2.0} TeV. To get an upper limit on the mixing angle ϕ (cf. Eq. (4.13)] for our subsequent analysis, we keep λ 3 = λ 8 fixed at two benchmark values (3.0 and 2.0) and vary the remaining nonzero quartic couplings λ 2 and λ η in the range [0.0, 3.0]. Our results on the maximum sin ϕ are shown in Fig. 5. We do not consider values of the quartic couplings exceeding 3.0 to be consistent with perturbativity considerations [80]. Each choice of mixing angle ϕ, and the parameters λ 2 , λ η , m h + , and m H + yields different minimization conditions deploying different solutions to the VEVs. We compare the values of the potential for all cases of CBM and CCM. If any one of the CBM is deeper than CCM, we reject the solution and rerun the algorithm with different initial conditions until we meet the requirement of electroweak minimum being deeper than all CBM.
For values of the mixing angle sin ϕ above the curves shown in Fig. 5 for a given m H+ , the potential develops CBM that are deeper than the electroweak minimum, which is unacceptable. This is mainly due to the fact that for these values of ϕ, the trilinear coupling µ becomes too large, which drives the potential to a deeper CBM [66], even for positive µ 2 η . From Fig. 5 it is found that sin ϕ < 0.23 for m H + = 2 TeV, while sin ϕ = 0.707 is allowed when m H + = 0.7 TeV. In all cases the maximum value of |µ| is found to be about 4.1 times the heavier mass m H + . Note that we have taken the maximum value of the mixing ϕ max = π/4 here, because for ϕ > π/4, the roles of h + and H + will be simply reversed, i.e., H + (h + ) will become the lighter (heavier) charged Higgs field. The CBM limits from Fig. 5 will be applied when computing neutrino NSI in the model.

Electroweak precision constraints
The oblique parameters S, T and U can describe a variety of new physics in the electroweak sector parametrized arising through shifts in the gauge boson self-energies [81,82] and impose important constraints from precision data. These parameters have been calculated in the context of the Zee model in Ref. [83]. We find that the T parameter imposes the most stringent constraint, compared to the other oblique parameters. The T parameter in the Zee model can be expressed as [83]: where the symmetric function F is given by In order to generate large NSI effects in the Zee model, the mixing between the singlet and the doublet charged scalar, parametrized by the angle ϕ, should be significant. This mixing contributes to the gauge boson self-energies and will therefore be bounded from the experimental value of the T parameter: T = 0.01 ± 0.12 [84]. For simplicity, we assume no mixing between the neutral CP-even scalars h and H. Furthermore, we take the heavy neutral CP-even (H) and odd (A) scalars to be degenerate in mass. In Fig. 6, we have shown our results from the T parameter constraint, allowing for two standard deviation error bar, in the heavy neutral and charged Higgs mass plane. Here we have fixed the light charged scalar mass m h + = 100 GeV. As shown in the figure, when the masses m H and m H ± are nearly equal (along the diagonal), the T parameter constraint is easily satisfied.
From Fig. 6, we also find that for specific values of m H and m H ± , there is an upper limit on the mixing sin ϕ. This is further illustrated in Fig. 7. Here, the colored regions (both green and red) depict the allowed parameter space in m + H − sin ϕ plane resulting from the T parameter constraint. For example, if we set m H = 0.7 TeV, the maximum mixing that is allowed by T parameter is (sin ϕ) max = 0.63. The mass splitting between We choose λ 5 = −λ 4 , which would maximize the mass splitting, as long as the quartic couplings remain perturbative. The red region in Fig. 7 depicts the scenario where |λ 4 |, |λ 5 | > 3.0, which we discard from perturbativity requirements in a conservative approach. Satisfying this additional requirement that these couplings be less than 3.0, we get an upper limit on sin ϕ < 0.59. For the degenerate case m H ± = m H with λ 4 = λ 5 , the upper limit is stronger: sin ϕ < 0.49.

Charged lepton flavor violation constraints
Charged lepton flavor violation (cLFV) is an integral feature of the Lagrangian Eq. (4.1) of the model. We can safely ignore cLFV processes involving the f αβ couplings which are assumed to be of the order of 10 −8 or so to satisfy the neutrino mass constraint, with Y αβ couplings being order one. Thus, we focus on cLFV proportional to Y αβ . Furthermore, as noted before, NSI arise proportional to (Y ee , Y µe , Y τ e ), where the first index refers to the neutrino flavor and the second to the charged lepton flavor in the coupling of charged scalars h + and H + . After briefly discussing the cLFV constraints arising from other Y αβ , we shall focus on the set (Y ee , Y µe , Y τ e ) relevant for NSI. The neutral scalar bosons H and A will mediate cLFV of the type µ → 3e and τ → µee at tree-level, while these neutral scalars and the charged scalars (h + , H + ) mediate processes of the type µ → eγ via one-loop diagrams. Both of these processes will be analyzed below. We derive limits on the couplings Y αβ as functions of the scalar masses. These limits need to be satisfied in the neutrino oscillation fit, see Sec. 4.13 for details. The constraints derived here will also be used to set upper  Tables VI and VII. We now turn to the derivation of these bounds.

4.5.1
α → β + γ decays The decay α → β +γ arises from one-loop diagrams shown in Fig. 8. The general expression for this decay rate can be found in Ref. [85]. Let us focus on the special case where the FCNC coupling matrix Y of Eq. (4.23) has nonzero entries either in a single row, or in a single column only. In this case, the chirality flip necessary for the radiative decay will occur on the external fermion leg. Suppose that only the right-handed component of fermion f α has nonzero Yukawa couplings with a scalar boson B and fermion F , parametrized as (4.29) The electric charges of fermions F and f are Q F and Q f respectively, while that of the boson B is Q B , which obey the relation Q f = Q F − Q B . The decay rate for f α → f β + γ is then given by Here α = e 2 /4π is the fine-structure constant, t = m 2 F /m 2 B , and the function f F (t) and f B (t) are given by (4.31) These expressions are obtained in the approximation m β m α . Let us apply these results to α → β + γ mediated by the charged scalars (h + , H + ) in the Zee model where the couplings have the form Y αβνα P R β h + sin ϕ, etc. Here Q F = 0, while Q B = +1. Eq. (4.30) then reduces to (with t 1)  Table VI: Constraints on Yukawa couplings as a function of heavy neutral scalar mass from α → β + γ processes.
However, nonzero values of (Y ee , Y µe , Y τ e ), needed for NSI, would lead to α → β + γ mediated by the heavy neutral scalars. Taking H and A to be degenerate, the Yukawa couplings are of the form¯ α P R β H. Thus, Q F = −1 and Q B = 0 in this case, leading to the decay width We show the constraints on these product of Yukawa couplings for a fixed mass of the neutral Higgs m H in Table VI. The severe constraint coming from µ → eγ process prevents the off-diagonal NSI parameter ε eµ from being in the observable range. However, ε eτ and ε µτ can be in the observable range, consistent with these constraints.

Electron anomalous magnetic moment
Another potential constraint comes from anomalous magnetic moment of leptons (g − 2) α , which could get contributions from both charged and neutral scalars in the Zee model. The heavy neutral scalar contribution can be ignored here. For the Yukawa couplings relevant for NSI, the charged scalar contribution to muon g−2 is also absent. The only non-negligible contribution is to the electron g − 2, which can be written at one-loop level as [89] Comparing this with ∆a e ≡ a exp e − a SM e = (−87 ± 36) × 10 −14 (where a e ≡ (g − 2) e /2), based on the difference between the experimental measurements [90] and SM calculations [91] with the updated value of the fine-structure constant [92], we find that the charged scalar contribution (4.34) goes in the right direction. However, for the allowed parameter space in m h + − Y ee sin ϕ plane (see Fig. 18), it turns out to be too small to explain the 2.4σ discrepancy in ∆a e . For example, with |Y τ e | sin ϕ = 0.75 and m h + = 150 GeV, which is a consistent choice (cf. Fig. 18), we would get ∆a e = −2.2 × 10 −14 , an order of magnitude too small to be relevant for experiments.

4.5.3
α →¯ β γ δ decays The Yukawa coupling matrix Y of the second Higgs doublet (cf. Eq. (4.23)) would lead to trilepton decay of charged leptons mediated by the neutral scalars of the theory. The treelevel Feynman diagrams for such decays are shown in Fig. 9. Partial rates for the trilepton decays are obtained in the limit when the masses of the decay products are neglected. The partial decay width for µ →ēee is given as follows: (4.35) The partial decay width for τ →¯ α β γ is given by (4.36) Here S = 1 (2) for β = γ (β = γ) is a symmetry factor. Using the total muon and tau decay widths, Γ tot µ = 3.00 × 10 −19 GeV and Γ tot τ = 2.27 × 10 −12 GeV respectively, we calculate the cLFV branching ratios for the processes µ − → e + e − e − , τ − → e + e − e − and τ − → e + e − µ − using Eqs. (4.35) and (4.36). We summarize in Table VII the current experimental bounds on these branching ratios and the constraints on the Yukawa couplings Y αβ as a function of mass of neutral Higgs boson m H = m A . It is clear from Table VII that these trilepton decays put more stringent bounds on product of Yukawa couplings compared to the bounds arising from loop-level α → β γ decays. This also implies that off-diagonal NSI are severely constrained.
As already noted, the light charged Higgs h + would mediate α → β + γ decay if more than one entry in a given row of Y is large. The heavy neutral Higgs bosons mediate trilepton decays of the leptons when there are more than one nonzero entry in the same column (or same row) of Y . This last statement is however not valid for the third column of Y . For example, nonzero Y τ τ and Y µτ will not lead to tree-level trilepton decay of τ . Table VII: Constraints on Yukawa couplings as a function of heavy neutral scalar mass from α →¯ β γ δ decay (with at least two of the final state leptons of electron flavor to be relevant for NSI).

Exp. bound Constraint
Apart from the first column of Y , we shall allow nonzero entries in the third column as well. In particular, for diagonal NSI ε αα , we need one Y αe entry for some α to be nonzero, and to avoid the trilepton constraints, the only other entry that can be allowed to be large is Y βτ with β = α. On the other hand, for off-diagonal NSI ε αβ (with α = β), we must allow for both Y αe and Y βe to be non-zero. In this case, however, the trilepton decay β → α ee is unavoidable and severely restricts the NSI as we will see in Sec. 4.12. Also, the other entry that can be populated is Y γτ with γ = α, β. This will lead to τ → + γ decays, which, however, do not set stringent limits on the couplings (cf. Table VI). Some benchmark Yukawa textures satisfying all cLFV constraints are considered in Sec. 4.13 to show consistency with neutrino oscillation data.

LEP contact interaction
Electron-positron collisions at center-of-mass energies above the Z-boson mass performed at LEP impose stringent constraints on contact interactions involving e + e − and a pair of fermions [95]. Integrating out new particles in a theory one can express their effect via higher-dimensional (generally dimension-6) operators. An effective Lagrangian, L eff , can parametrize the contact interaction for the process e + e − → ff with the form [96] where δ ef is the Kronecker delta function, f refers to the final sate fermions, g is the coupling strength, Λ is the new physics scale and η f ij = ±1 or 0, depending on the chirality structure. LEP has put 95% confidence level (CL) lower limits on the scale of the contact interaction Λ assuming the coupling g = √ 4π [95]. In the Zee model, the exchange of new neutral scalars (H and A) emerging from the second Higgs doublet will affect the process e + e − → + α − β (with α,β = e, µ, τ ), and therefore, the LEP constraints on Λ can be interpreted as a lower limit on the mass of the heavy neutral scalar, for a given set of Yukawa couplings. Here we assume that H and A are degenerate, and derive limits obtained by integrating out both fields.
In general, for + α − β → + γ − δ via heavy neutral scalar exchange, the effective Lagrangian in the Zee model can be written as By Fierz transformation, we can rewrite it in a form similar to Eq. (4.37): Thus, the only relevant chirality structures in Eq. (4.37) are LR and RL, and the relevant process for deriving the LEP constraints is e + e − → + α − α : Comparing this with Eq. (4.40), we obtain where Λ − corresponds to Λ with η LR = η RL = −1. The LEP constraints on Λ were derived in Ref. [95] for g = √ 4π, which can be translated into a lower limit on m H /|Y ee | using Eq. (4.42), as shown in Table VIII. Similarly, for e + e − → µ + µ − , Eq. (4.39) is Since for NSI, only Y µe (neutrino interaction with electron) is relevant, we can set Y eµ → 0, and compare Eq. (4.43) with Eq. (4.40) to get a constraint on m H /|Y µe |, as shown in Table VIII. Similarly, for e + e − → τ + τ − , we can set Y eτ → 0 and translate the LEP limit on Λ − into a bound on m H /|Y τ e |, as shown in Table VIII.
The LEP constraints from the processes involving qq final states, such as e + e − → cc and e + e − → bb, are not relevant in our case, since the neutral scalars are leptophilic. We will use the limits quoted in Table VIII while deriving the maximum NSI predictions in the Zee model.

LEP constraints on light neutral scalar
The LEP contact interaction constraints discussed in Sec. 4.6 are not applicable if the neutral scalars H and A are light. In this case, however, the cross section of e + e − → + TeV   Table VIII: Constraints on the ratio of heavy neutral scalar mass and the Yukawa couplings from LEP contact interaction bounds.
can still be modified, due to the t-channel contribution of H/A, which interferes with the SM processes. We implement our model file in FeynRules package [97] and compute the e + e − → + α − α cross-sections in the Zee model at the parton-level using MadGraph5 event generator [98]. These numbers are then compared with the measured cross sections [95,99] to derive limits on m H/A as a function of the Yukawa couplings Y αe (for α = e, µ, τ ). For a benchmark value of m H = m A = 130 GeV, we find the following constraints on the Yukawa couplings Y αe relevant for NSI: This implies that the second charged scalar H + can also be light, as long as it is allowed by other constraints (see Fig. 11). We will use this finding to maximize the NSI prediction for the Zee model (see Sec. 4.12.2).

LHC constraints
Most of the LHC searches for heavy neutral scalars are done in the context of either MSSM or 2HDM, which are not directly applicable in our case because H and A do not couple to quarks, and therefore, cannot be produced via gluon fusion. The dominant channel to produce the neutral scalars in our case at the LHC is via an off-shell Z boson: pp → Z → HA → + − + − . 7 Most of the LHC multilepton searches assume a heavy ZZ ( ) resonance [101,102], which is not applicable in this case. The cross section limits from inclusive multilepton searches, mostly performed in the SUSY context with large missing transverse energy [103,104], turn out to be weaker than the LEP constraints derived above.

Collider constraints on light charged scalar
In this section, we discuss the collider constraints on the light charged scalar h ± in the Zee model from various LEP and LHC searches.

Constraints from LEP searches
At LEP, h ± can be pair-produced through the s-channel Drell-Yan process mediated by either γ or Z boson (see Fig. 10 (a)). It can also be pair-produced through the t-channel to parity [100]. Figure 10: Feynman diagrams for pair-and single-production of singly-charged scalars h ± at e + e − collider.
processes mediated by a light neutrino (see Fig. 10 (b)). In addition, it can be singlyproduced either in association with a W boson (see Fig. 10 (c)) or via the Drell-Yan channel in association with leptons (see Fig. 10 (d)).
It is instructive to write down the explicit formula for the pair-production (Figs. 10 (a) and 10 (b) cross section: where β = 1 − 4m 2 h + /s, s is the squared center-of-mass energy, e and g are the electromagnetic and SU (2) L coupling strengths, respectively, and c w ≡ cos θ w , s w ≡ sin θ w (θ w being the weak mixing angle). Note that the t-channel cross section depends on the Yukawa coupling Y αe , and it turns out there is a destructive interference between the s and t-channel processes. Similarly, the differential cross section for the production of h ± W ∓ ( Fig. 10 (c)) is given by where θ is the angle made by the outgoing h ± with respect to the initial e − -beam direction, (4.49) The analytic cross section formula for the single-production of charged Higgs via Drell-Yan process ( Fig. 10 (d)) is more involved due to the three-body phase space and is not given here. We implement our model file in FeynRules package [97] and compute all the cross-sections at the parton-level using MadGraph5 event generator [98].
Once produced on-shell, the charged scalar will decay into the leptonic final states ν α β through the Yukawa coupling Y αβ . Since we are interested in potentially large NSI effects, the charged scalar must couple to the electron. Due to stringent constraints from cLFV processes, especially the trilepton cLFV decays (see Table VII), which is equally applicable for the product of two Yukawa entries either along a row or column, both Y αe and Y αµ (or Y αe and Y βe ) cannot be large simultaneously. So we consider the case where BR eν + BR τ ν = 100% and BR µν is negligible, in order to avoid more stringent limits from muon decay. 8 Electron channel: For a given charged scalar decay branching ratio to electrons, BR eν , we can reinterpret the LEP selectron searches [105] to put a constraint on the charged scalar mass as a function of BR eν . In particular, the right-handed selectron pair-production e + e − → e R e R , followed by the decay of each selectron to electron and neutralino, e R → e R + χ 0 , will mimic the e + e − νν final state of our case in the massless neutralino limit. So we use the 95% CL observed upper limits on the e R e R production cross section [105] for m χ = 0 as an experimental upper limit on the quantity and derive the LEP exclusion region in the plane of charged scalar mass and BR eν , as shown in Fig. 11 (a) by the orange-shaded region. Here we have chosen Y ee sin ϕ = 0.1 and varied Y τ α (with α = µ or τ ) to get the desired branching ratios. We find that for BR eν = 1, charged scalar masses less than 100 GeV are excluded. For BR eν < 1, these limits are weaker, as expected, and the charged scalar could be as light as 97 GeV (for BR eν = 0.33), if we just consider the LEP selectron (as well as stau, see below) searches. Fig. 11 (b) shows the same constraints as in Fig. 11 (a), but for the case of Y ee sin ϕ = 0.2. The LEP selectron constraints become stronger as we increase Y ee and extend to smaller BR eν . However, the mass limit of 100 GeV for BR eν = 1 from Fig. 11 (a) still holds here. This is because the charged scalar pair-production cross section drops rapidly for m h + > 100 GeV due to the kinematic threshold of LEP II with √ s = 209 GeV and is already below the experimental cross section limit even for Y ee sin ϕ = 0.2. In this regime, the single-production channel in Fig. 10 (d) starts becoming important, despite having a three-body phase space suppression.
Figs. 11 (c) and 11 (d) show the same constraints as in Fig. 11 (a) and 11 (b) respectively, but for the Y ee = 0 case. Here we have fixed Y τ e sin ϕ and varied Y τ α (with α = e or µ) to get the desired branching ratios. In this case, the single-production channel in association with the W boson (cf. Fig. 10 (c)) goes away, and therefore, the limits from selectron and stau searches become slightly weaker. Note that for the NSI purpose, we must have a non-zero Y αe (for α = e, µ or τ ). Therefore, the t-channel contribution to the pair-production (cf. Fig. 10 (b)), as well as the Drell-Yan single-production channel are always present. 9 Tau channel: In the same way, we can also use the LEP stau searches [105] to derive an upper limit on and the corresponding LEP exclusion region in the plane of charged scalar mass and BR τ ν , as shown in Fig. 11 by the blue shaded region. We find that for BR τ ν = 1, charged scalar masses less than 104 (105) GeV are excluded for Y ee sin φ = 0.1 (0.2).
For BR τ ν = 0, a slightly stronger limit can be obtained from the LEP searches for the charged Higgs boson pairs in the 2HDM [107]. Their analysis focused on three kinds of final states, namely, τ ντ ν, csτ ν (orcsτ ν) and cscs, under the assumption that BR τ ν +BR cs = 1, which is valid in the 2HDM as the couplings of the charged Higgs boson to the SM fermions are proportional to the fermion masses. In our case, the observed LEP upper limit on σ(e + e − → h + h − )BR 2 τ ν for BR τ ν = 1 can be recast into an upper limit on and the corresponding exclusion region is shown in Fig. 11 by the green shaded region. We can also use the LEP cross section limit on csτ ν for BR τ ν = 1 as an upper limit on σ(e + e − → h ± W ∓ )BR τ ν BR W →cs and the corresponding exclusion region is shown in Fig. 11 by the cyan shaded region, which is found to be weaker than the τ ντ ν mode. Figure 11: Collider constraints on light charged scalar h ± in the Zee model for (a) We plot the h ± branching ratios to τ ν and eν (with the sum being equal to one) as a function of its mass. All shaded regions are excluded: Blue and orange shaded regions from stau and selectron searches at LEP (see Sec. 4.7.1); purple region from selectron searches at LHC (see Sec. 4.7.2); yellow, brown, and pink regions from W universality tests in LEP data for µ/e, τ /e, and τ /µ sectors respectively (see Sec. 4.8); light green and gray regions from tau decay universality and lifetime constraints respectively (see Sec. 4.9). The W universality constraints do not apply in panels (b) and (c), because the h ± W ∓ production channel in Fig. 10 (c) vanishes in the Y ee → 0 limit. Figure 12: Feynman diagrams for pair production and single production of singly-charged scalars h ± at LHC.

Constraints from LHC searches
As for the LHC constraints, there is no t-channel contribution to the singlet charged scalar production. The only possible channel for pair-production is the s-channel Drell-Yan process pp → γ /Z → h + h − (see Fig. 12 (a)), followed by the leptonic decay of h ± → ν. There are also single-production processes as shown in Fig. 12 (b)-(d), which turn out to be less important. The relevant LHC searches are those for right-handed selectrons/staus: , which will mimic the + ν − ν final states from h + h − decay in the massless neutralino limit. The √ s = 13 TeV LHC stau searches focus on the stau mass range above 100 GeV and it turns out that the current limits [108] on the stau pairproduction cross section are still a factor of five larger than the h + h − pair-production cross section in our case; therefore, there are no LHC limits from the tau sector. A √ s = 8 TeV ATLAS analysis considered the mass range down to 80 GeV [109]; however, the observed cross section is still found to be larger than the theoretical prediction in our case even for BR τ ν = 1.
As for the selectron case, we take the √ s = 13 TeV CMS search [110], which focuses on the selectron masses above 120 GeV, and use the observed cross section limit on σ(pp → e + R χ 0 e − R χ 0 ) to derive an upper limit on σ(pp → h + h − )BR 2 eν , which can be translated into a bound on the charged scalar mass, as shown in Fig. 11 by the purple shaded regions. It is evident that the LHC limits can be evaded by going to larger BR τ ν 0.4, which can always be done for any given Yukawa coupling Y αe by choosing an appropriate Y βτ . This however may not be the optimal choice for NSI, especially for Y ee = 0, where the lepton universality constraints restrict us from having a larger BR τ ν . Thus, the LHC constraints will be most relevant for ε ee , as we will see in Fig. 18 (a).

Constraints from lepton universality in W ± decays
The presence of a light charged Higgs can also be constrained from precision measurements of W boson decay rates. The topology of the charged Higgs pair production h + h − (Fig. 10 (a) and 10 (b)) and the associated production h ± W ∓ (Fig. 10 (c)) is very similar to the W + W − pair production at colliders, if the charged Higgs mass is within about 20 GeV of the W boson mass. Thus, the leptonic decays of the charged Higgs which are not necessarily flavor-universal can be significantly constrained from the measurements of lepton universality in W decays. From the combined LEP results [111], the constraints on the ratio of W branching ratios to leptons of different flavors are as follows: Note that while the measured value of R µ/e agrees with the lepton universality prediction of the SM, R SM µ/e = 1, within 1.1σ CL, the W branching ratio to tau with respect to electron is about 1.8σ and to muon is about 2.7σ away from the SM prediction: R SM τ / = 0.9993 (with = e, µ), using the one-loop calculation of Ref. [112].
The best LEP limits on lepton universality in W decays come from the W + W − pair-production channel, where one W decays leptonically, and the other W hadronically, i.e., e + e − → W + W − → νqq [111]. However, due to the leptophilic nature of the charged Higgs h ± in our model, neither the e + e − → h + h − channel (Figs. 10 (a) and 10 (b)) nor the Drell-Yan single-production channel ( Fig. 10 (d)) will lead to νqq final state. So the only relevant contribution to the W universality violation could come from the h ± W ∓ production channel (Fig. 10 (c)), with the W decaying hadronically and h ± decaying leptonically. The pure leptonic channels (eνeν and µνµν) have ∼ 40% uncertainties in the measurement and are therefore not considered here.
Including the h ± W ∓ contribution, the modified ratios R / can be calculated as follows: where σ(W + W − ) and σ(h ± W ∓ ) are the production cross sections for e + e − → W + W − and e + e − → h ± W ∓ respectively, BR W ν denotes the branching ratio of W → ν (with = e, µ, τ ), whereas BR ν denotes the branching ratio of h ± → ν as before (with = e, τ ). At LEP experiment, the W + W − pair production cross section σ W + W − is computed to be 17.17 pb at √ s = 209 GeV [111]. Within the SM, W ± decays equally to each generation of leptons with branching ratio of 10.83% and decays hadronically with branching ratio of 67.41% [84]. We numerically compute using MadGraph5 [98] the h ± W ∓ cross section at √ s = 209 GeV as a function of m h ± and BR ν , and compare Eq.  This is shown in Figs. 11 (a) and 11 (b) by yellow, brown, and pink shaded regions for µ/e, τ /e, and τ /µ universality tests, respectively. Note that these constraints are absent in Figs. 11 (c) and 11 (d), because when Y ee = 0, there is no W ± h ∓ production at LEP (cf. Fig. 10 (c) in the Zee model. But when Y ee is relatively large, these constraints turn out to be some of the most stringent ones in the m h + -BR ν plane shown in Figs. 11 (a) and 11 (b), and rule out charged scalars below 110 GeV (129 GeV) for Y ee sin ϕ = 0.1 (0.2). These constraints are not applicable for m h ± > 129 GeV, because h ± W ∓ can no longer be produced on-shell at LEP II with maximum √ s = 209 GeV. As mentioned before, the measured W branching ratio to tau with respect to muon is 2.7σ above the SM prediction. Since in our case, h ± decays to either eν or τ ν, but not µν, this contributes to R τ µ only in the numerator, but not in the denominator. Therefore, the 2.7σ discrepancy can be explained in this model, as shown by the allowed region between the upper and lower pink-dashed curves in Fig. 11 (a) with Y ee sin ϕ = 0.1. 10 The upper pink-shaded region with larger BR τ ν gives R τ µ > 1.122, which is above the allowed 2σ range given in Eq. (4.55). On the other hand, the lower pink-shaded region with smaller BR τ ν gives R τ µ < 1.018, which is below the allowed 2σ range given in Eq. (4.55). For larger Yukawa coupling Y ee , as illustrated in Fig. 11 (b) with Y ee sin ϕ = 0.2, the whole allowed range of parameter space from R τ /µ shifts to lower values of BR τ ν . This is because the h ± W ∓ production cross section σ(h ± W ∓ ) in Eq. (4.56) is directly proportional to |Y ee | 2 , and therefore, for a large Y ee , a smaller BR τ ν would still be compatible with the R τ /µ -preferred range.

Constraints from tau decay lifetime and universality
In order to realize a light charged scalar h − consistent with LEP searches, we have assumed that the decay h − → τν β proceeds with a significant branching ratio. h − also has coupling with eν α , so that non-negligible NSI is generated. When these two channels are combined, we would get new decay modes for the τ lepton, as shown in Fig. 13. This will lead to deviation in τ -lifetime compared to the SM expectation. The new decay modes will also lead to universality violation in τ decays, as the new modes preferentially lead to electron final states. Here we analyze these constraints and evaluate the limitations these pose for NSI.
The effective four-fermion Lagrangian relevant for the new τ decay mode is given by This can be recast, after a Fierz transformation, as This can be directly compared with the SM τ decay Lagrangian, given by It is clear from here that the new decay mode will not interfere with the SM model (in the limit of ignoring the lepton mass), since the final state leptons have opposite helicity in the two decay channels. The width of the τ lepton is now increased from its SM value by a factor 1 + ∆, with ∆ given by [114] The global fit result on τ lifetime is τ τ = (290.75 ± 0.36) × 10 −15 s, while the SM prediction is τ SM τ = (290.39 ± 2.17) × 10 −15 s [84]. Allowing for 2σ error, we find ∆ ≤ 1.5%. If the only decay modes of h − are h − →ν α e − and h − →ν β τ − , then we can express |Y βτ | 2 in terms of |Y αe | 2 as Using this relation, we obtain where ε αα is the diagonal NSI parameter for which the expression is derived later in Eq. (4.79). Therefore, a constraint on ∆ from the tau lifetime can be directly translated into a constraint on ε αα : (4.64) An even stronger limit is obtained from e−µ universality in τ decays. The experimental central value prefers a slightly larger width for τ → µνν compared to τ → eνν. In our scenario, h − mediation enhances τ → eνν relative to τ → µνν. We have in this scenario Figure 14: (a) New contribution to h → γγ decay mediated by charged scalar loop.
which constrains ∆ ≤ 0.002, obtained by using the measured ratio Γ(τ →µνν) Γ(τ →eνν) = 0.9762 ± 0.0028 [84], and allowing 2σ error. This leads to a limit In deriving the limits on a light charged Higgs mass from LHC constraints, we have imposed the τ decay constraint as well as the universality constraint on ∆, see Fig. 11. Avoiding the universality constraint by opening up the τ → µνν channel will not work, since that will be in conflict with µ → eνν constraints, which are more stringent. The Michel parameters in τ decay will now be modified [115]. While the ρ and δ parameters are unchanged compared to their SM value of 3/4, ξ is modified from its SM value of 1 to However, the experimental value is ξ = 0.985 ± 0.030 [84], which allows for significant room for the new decay. Again, our choice of Yukawa couplings does not modify the µ → eνν decay, and is therefore, safe from the Michel parameter constraints in the muon sector, which are much more stringent.

Constraints from Higgs precision data
In this subsection, we analyze the constraints on light charged scalar from LHC Higgs precision data. Both ATLAS and CMS collaborations have performed several measurements of the 125 GeV Higgs boson production cross sections and branching fractions at the LHC, both in Run I [116] and Run II [117,118]. Since all the measurements are in good agreement with the SM expectations, any exotic contributions to either production or decay of the SM-like Higgs boson will be strongly constrained. In the Zee model, since the light charged scalar is leptophilic, it will not affect the production rate of the SM-like Higgs h (which is dominated by gluon fusion via top-quark loop). However, it gives new contributions to the loop-induced h → γγ decay (see Fig. 14 (a)) and mimics the tree-level h → W W → 2 2ν channel via the exotic decay mode h → h ± h ∓ → h ± ν → 2 2ν (see Fig. 14 (b)). Both these contributions are governed by the effective hh + h − coupling given by Therefore, the Higgs precision data from the LHC can be used to set independent constraints on these Higgs potential parameters, as we show below. The Higgs boson yield at the LHC is characterized by the signal strength, defined as the ratio of the measured Higgs boson rate to its SM prediction. For a specific production channel i and decay into specific final states f , the signal strength of the Higgs boson h can be expressed as where µ i (with i = ggF, VBF, V h, and tth) and µ f (with f = ZZ , W W , γγ, τ + τ − , bb) are the production and branching rates relative to the SM predictions in the relevant channels. As mentioned above, the production rate does not get modified in our case, so we will set µ i = 1 in the following. As for the decay rates, the addition of the two new channels shown in Fig. 14 will increase the total Higgs decay width, and therefore, modify the partial widths in all the channels.
To derive the Higgs signal strength constraints on the model parameter space, we have followed the procedure outlined in Ref. [69,119], using the updated constraints on signal strengths reported by ATLAS and CMS collaboration for all individual production and decay modes at 95% CL, based on the √ s = 13 TeV LHC data. The individual analysis by each experiment examines a specific Higgs boson decay mode corresponding to various production processes. We use the measured signal strengths in the following dominant decay modes for our numerical analysis: 130] and h → bb [131][132][133].
We formulate the modified h → γγ decay rate as where the scaling factor κ γ is given by , (4.71) where N f c = 3 (1) is the color factor for quark (lepton), f is the sum over the SM fermions f with charge Q f , and the loop functions are given by [134] A 0 (τ ) = −τ + τ 2 f (τ ), (4.72) The parameters τ i = 4m 2 i /m 2 h are defined by the corresponding masses of the heavy particles in the loop. For the fermion loop, only the top quark contribution is significant, with the next leading contribution coming from the bottom quark which is an 8% effect. Note that the new contribution in Eq. (4.71) due to the charged scalar can interfere with the SM part either constructively or destructively, depending on the sign of the effective coupling λ hh + h − in Eq. (4.68).
As for the new three-body decay mode h → h ± h ∓ → h ± ν, the partial decay rate is given by where Y is the Yukawa coupling defined in Eq. (4.19), The partial decay widths of h in other channels will be the same as in the SM, but their partial widths will now be smaller, due to the enhancement of the total decay width. A comparison with the measured signal strengths therefore imposes an upper bound on the effective coupling λ hh ± h ∓ which is a function of the cubic coupling µ, quartic couplings λ 3 and λ 8 , and the mixing angle sin ϕ (cf. Eq. (4.69)). For suppressed effective coupling λ hh ± h ∓ to be consistent with the Higgs observables, we need some cancellation between the cubic and quartic terms. In order to have large NSI effect, we need sufficiently large mixing sin ϕ, which implies large value of µ (cf. Eq. (4.13)). In order to find the maximum allowed value of sin ϕ, we take λ 3 = λ 8 in Eq. (4.69) and show in Fig. 15 the Higgs signal strength constraints in the λ 8 −sin ϕ plane. The red, blue, yellow, cyan, and green shaded regions are excluded by the signal strength limits γγ, W W , ZZ , τ τ , and bb decay modes, respectively. We have fixed the light charged Higgs mass at 100 GeV, and the different panels are for different benchmark values of the heavy charged Higgs mass: m H + = 700 GeV (upper left), 2 TeV (upper right), 1.6 TeV (lower left) and 450 GeV (lower right). The first choice is the benchmark value we will later use for NSI studies, while the other three values correspond to the minimum allowed values for the heavy neutral Higgs mass (assuming it to be degenerate with the heavy charged Higgs to easily satisfy the T -parameter constraint (cf. Sec. 4.4)) consistent with the LEP contact interaction bounds for O(1) Yukawa couplings (cf. Sec. 4.6). From Fig. 15, we see that the h → γγ signal strength gives the most stringent constraint. If we allow λ 8 to be as large as 3, then we can get maximum value of sin ϕ up to 0.67 (0.2) for m H + = 0.7 (2) TeV.
In addition to the modified signal strengths, the total Higgs width is enhanced due to the new decay modes. Both ATLAS [102] and CMS [135] collaborations have put 95% CL upper limits on the Higgs boson total width Γ h from measurement of off-shell production in the ZZ → 4 channel. Given the SM expectation Γ SM h ∼ 4.1 MeV, we use the CMS upper limit on Γ h < 9.16 MeV [135] to demand that the new contribution (mostly from h → h ± h ∓ , because the h → γγ branching fraction is much smaller) must be less than

5.1
MeV. This is shown in Fig. 15 by the grey shaded region (only visible in the upper right panel), which turns out to be much weaker than the signal strength constraints in the individual channels. Figure 16: Feynman diagrams for charged scalar contributions to monophoton signal at LEP.

Monophoton constraint from LEP
Large neutrino NSI with electrons inevitably leads to a new contribution to the monophoton process e + e − → ννγ that can be constrained using LEP data [136]. In the SM, this process occurs via s-channel Z-boson exchange and t-channel W -boson exchange, with the photon being emitted from either the initial state electron or positron or the intermediate state W boson. In the Zee model, we get additional contributions from t-channel charged scalar exchange (see Fig. 16). Both light and heavy charged scalars will contribute, but given the mass bound on the heavy states from LEP contact interaction, the dominant contribution will come from the light charged scalar. The total cross section for the process e + e − → ν ανβ γ can be expressed as σ = σ SM + σ NS , where σ SM is the SM cross section (for α = β) and σ NS represents the sum of the pure non-standard contribution due to the charged scalar and its interference with the SM contribution. Note that since the charged scalar only couples to right-handed fermions, there is no interference with the W -mediated process (for α = β = e). Moreover, for either α or β not equal to e, the W contribution is absent. For α = β, the Z contribution is also absent.
The monophoton process has been investigated carefully by all four LEP experiments [84], but the most stringent limits on the cross section come from the L3 experiment, both on [137] and off [138] Z-pole. We use these results to derive constraints on the charged scalar mass and Yukawa coupling. The constraint |σ − σ exp | ≤ δσ exp , where σ exp ± δσ exp is the experimental result, can be expressed in the following form: We evaluate the ratio σ exp /σ SM by combining the L3 results [137,138] with an accurate computation of the SM cross section, both at Z-pole and off Z-pole. Similarly, we compute the ratio σ NS /σ SM numerically as a function of the charged scalar mass m h + and the Yukawa coupling Y αβ sin ϕ. For comparison of cross sections at Z-pole, we adopt the same event acceptance criteria as in Ref. [137], i.e., we allow photon energy within the range 1 GeV < E γ < 10 GeV and the angular acceptance 45°< θ γ < 135°. Similarly, for the off Z-pole  analysis, we adopt the same event topology as described in Ref. [138]: i.e., 14°< θ γ < 166°, 1 GeV < E γ , and p γ T > 0.02 √ s. We find that the off Z-pole measurement imposes more stringent bound than the Z-pole measurement bound. As we will see in the next section (see Fig. 18), the monophoton constraints are important especially for the NSI involving tau-neutrinos. We also note that our monophoton constraints are somewhat weaker than those derived in Ref. [139] using an effective four-fermion approximation.

NSI predictions
The new singly-charged scalars η + and H + 2 in the Zee Model induce NSI at tree level as shown in Fig. 17. Diagrams (a) and (d) are induced by the pure singlet and doublet components of the charged scalar fields and depend on the Yukawa couplings f and Y respectively (cf. Eqs. (4.1) and (4.19)). On the other hand, diagrams (b) and (c) are induced by the mixing between the singlet and doublet fields, and depend on the combination of Yukawa couplings and the mixing angle ϕ (cf. Eq. (4.13)). As mentioned in Sec. 4.2, satisfying the neutrino mass requires the product f · Y to be small. For Y ∼ O(1), we must have f ∼ 10 −8 to get m ν ∼ 0.1 eV (cf. Eq. (4.21)). In this case, the NSI from Fig. 17 (a) and (c) are heavily suppressed. So we will only consider diagrams (b) and (d) for the following discussion and work in the mass basis for the charged scalars, where η + and H + 2 are replaced by h + and H + (cf. Eq. (4.12)).
The effective NSI Lagrangian for the contribution from Fig. 17 (b) is given by Thus, the diagonal NSI parameters ε αα depend on the Yukawa couplings |Y αe | 2 , and are always positive in this model, whereas the off-diagonal ones ε αβ (with α = β) involve the product Y αe Y βe and can be of either sign, or even complex. Also, we have a correlation between the diagonal and off-diagonal NSI: which is a distinguishing feature of the model. Fig. 17 (d) gives a sub-dominant NSI contribution as follows: Hence, the total matter NSI induced by the charged scalars in the Zee model can be expressed as To get an idea of the size of NSI induced by Eq. (4.82), let us take the diagonal NSI parameters from the light charged scalar contribution in Eq. (4.79): Thus, for a given value of m h + , the NSI are maximized for maximum allowed values of |Y αe | and sin ϕ. Following Eq. (4.68), we set the trilinear coupling λ hh + h − → 0, thus minimizing the constraints from Higgs signal strength. We also assume λ 3 = λ 8 to get Now substituting this into Eq. (4.13), we obtain .

(4.85)
Furthermore, assuming the heavy charged and neutral scalars to be mass-degenerate, from LEP contact interaction constraints (cf. Sec. 4.6), we have where Λ α = 10 TeV, 7.9 TeV and 2.2 TeV for α = e, µ, τ , respectively [95]. Combining Eqs. (4.83), (4.85) and (4.86), we obtain Using benchmark values of m h + = 100 GeV and λ 8 = 3, we obtain:  Fig. 18 by the light purple shaded region. The model predictions for NSI are then compared with the current experimental constraints from neutrino-electron scattering experiments (red shaded), as well as the global fit results from neutrino oscillation plus COHERENT data (orange shaded); see Table IX for more details. 11 For neutrino-electron scattering constraints, we only considered the constraints on ε eR αβ [140][141][142][143], since the dominant NSI in the Zee model always involves right-handed electrons (cf. Eq. (4.78)). These scattering experiments impose the strongest limits for ε µµ and ε τ τ , restricting them to be less than 3.8% and 43%, respectively, although the model allows for much larger NSI (cf. Fig. 18). For ε µµ , we have rederived the CHARM II limit following Ref. [140], but using the latest PDG value for s 2 w = 0.22343 (on-shell) [84]. Specifically, we used the CHARM II measurement of the Z-coupling to right-handed electrons g e R = 0.234 ± 0.017 obtained from their ν µ e → νe data [144] and compared with the SM value of (g e R ) SM = s 2 w to obtain a 90% CL limit on ε µµ < 0.038, which is slightly weaker than the limit of 0.03 quoted in Ref. [141].
As for the global-fit constraints, we use the constraints on ε p αβ from Ref. [60], assuming that these will be similar for ε e αβ due to charge-neutrality in matter. Also shown (blue solid lines) are the future sensitivity at long baseline experiments, such as DUNE with 300 kt.MW.yr and 850 kt.MW.yr of exposure, derived at 90% CL using GloBES3.0 [145] with the DUNE CDR simulation configurations [146]. Here we have used δ (true) = −π/2 for the true value of the Dirac CP phase and marginalized over all other oscillation parameters [65].
Taking into account all existing constraints and this possibility of light h + and H + , the maximum possible allowed values of the NSI parameters in the Zee model are shown in the second column of Table IX, along with the combination of the relevant constraints limiting each NSI parameter (shown in parentheses). Thus, we find that for the diagonal NSI, one cannot get significantly large ε ee and ε µµ , but ε τ τ as large as 43% can be allowed in this model and there is a good portion of the allowed region for ε τ τ within reach of DUNE sensitivity. As for the off-diagonal NSI, they require the presence of at least two non-zero Yukawa couplings Y αe , and their products are all heavily constrained from cLFV; therefore, one cannot get sizable off-diagonal NSI in the Zee model that can be probed at DUNE or any other neutrino experiment.  [141], TEX-ONO [142] and BOREXINO [143]. Orange shaded region in (c) is excluded by global fit constraints from neutrino oscillation+COHERENT data [60]. We also show the future DUNE sensitivity in blue solid lines, for both 300 kt.MW.yr and 850 kt.MW.yr exposure [65].    . We also impose the constraints from neutrino-electron scattering experiments (as shown in the third column), like CHARM-II [141], TEXONO [142] and BOREXINO [143] (only eR αβ are considered, cf. Eq. (4.78)) as well as the global fit constraints (as shown in the fourth column) from neutrino oscillation+COHERENT data [60] (only ε p αβ are considered), whichever is stronger. The maximum allowed value for each NSI parameter is obtained after scanning over the light charged Higgs mass (see Figs. 18 and 19) and the combination of all relevant constraints limiting the NSI are shown in parentheses in the second column. In the last column, we also show the future DUNE sensitivity for 300 kt.MW.yr exposure (and 850 kt.MW.yr in parentheses) [65].

Light neutral scalar case
Now we consider the case where the neutral scalars H and A are light, so that the LEP contact interaction constraints (cf. 4.6) are not applicable. In this case, both h + and H + contributions to the NSI in Eq. (4.82) should be kept. For concreteness, we fix m H + = 130 GeV to allow for the maximum H + contribution to NSI while avoiding the lepton universality constraints on H + (cf. Sec. 4.8). We also choose the neutral scalars H and A to be nearly mass-degenerate with the charged scalar H + , so that the T -parameter and CBM constraints are easily satisfied. The Higgs decay constraints can also be significantly relaxed in this case by making λ hh + h − → 0 in Eq. (4.68). The NSI predictions for this special choice of parameters are shown in Fig. 20. Note that for higher m h + , the NSI numbers are almost constant, because of the m H + contribution which starts dominating. We do not show the off-diagonal NSI plots for this scenario, because the cLFV constraints still cannot be overcome (cf. Fig. 19).

Consistency with neutrino oscillation data
In this section, we show that the choice of the Yukawa coupling matrix used to maximize our NSI parameter values is consistent with the neutrino oscillation data. The neutrino where c ij ≡ cos θ ij , s ij ≡ sin θ ij , θ ij being the mixing angle between different flavor eigenstates i and j, and δ is the Dirac CP phase. We diagonalize the neutrino mass matrix (4.21) numerically, assuming certain forms of the Yukawa coupling matrices given below. The unitary matrix thus obtained is converted to the mixing angles θ ij using the following relations from Eq. (4.90): (4.91) Since the NSI expressions in Eq. (4.82) depend on Y αe (the first column of the Yukawa matrix), we choose the following three sets of benchmark points (BPs) for Yukawa textures to satisfy all the cLFV constraints, see Tables VI and VII. For simplicity, we also take all the elements of Yukawa matrix to be real.     Yee . Similarly, for BPs II and III, one can absorb Y µµ and Y τ τ respectively in the overall factor a 0 to get the mass matrix parameters in terms of the ratios x i and y ij .
For each set of Yukawa structure, we show in Table X the best-fit values of the parameters x i , y ij and a 0 . For BP I and II, we obtain inverted hierarchy (IH) and for BP III, we get normal hierarchy (NH) of neutrino masses. The model predictions for the neutrino oscillation parameters in each case are shown in Table XI, along with the 3σ allowed range from a recent NuFit4 global analysis [147]. It is clear that the fits for all the three sets are in very good agreement with the observed experimental values. We note here that the NuFit4 analysis does not include any NSI effects, which might affect the fit results; however, it is sufficient for the consistency check of our benchmark points. A full global analysis of the oscillation data in presence of NSI to compare with our benchmark points is beyond the scope of this work.
In addition to the best fit results in the tabulated format, we also display them in Fig. 21 in the two-dimensional projections of 1σ, 2σ and 3σ confidence regions of the global fit results [147] (without inclusion of the Super-K atmospheric ∆χ 2 -data). Colored regions (grey, magenta, cyan) are for normal hierarchy, whereas regions enclosed by solid, dashed, dotted lines are for inverted hierarchy. The global-fit best-fit points, along with the Figure 21: Global oscillation analysis obtained from NuFit4 [147] for both Normal hierarchy (NH) and Inverted hierarchy (IH) compared with our model benchmark points (BP1, BP2, BP3). Gray, Magenta, and Cyan colored contours represent 1σ, 2σ, and 3σ CL contours for NH, whereas solid, dashed, and dotted lines respectively correspond to 1σ, 2σ, and 3σ CL contours for IH. Red, purple, and (blue, black, brown) markers are respectively best-fit from NuFit for IH and NH, and benchmark points I, II and III for Yukawa structures given in Eqs.

NSI in one-loop leptoquark model
There are only four kinds of scalar leptoquarks that can interact with the neutrinos at the renormalizable level in the SM (see Table I): Ld c Ω, LQχ , LQρ and Lu c δ. 12 In this section and next, we discuss neutrino mass models with various combinations of these LQs. Our focus is again the range of neutrino NSI that is possible in these models. We note in passing that all these scalar LQ scenarios have gained recent interest in the context of semileptonic  [149]. We start with a LQ variant of the Zee model that generates small neutrino masses at one-loop level, via the operator is O 3b (cf. Eq. (2.2c)). It turns out that O 3b will induce neutrino masses at one-loop, while O 3a , owing to the SU (2) L index structure, will induce m ν at the two-loop level. A UV complete model of O 3a will be presented in Sec. 7.2.3. More precisely, the model of this section corresponds to O 8 3 of Table III, which involves two LQ fields and no new fermions. All other realizations of O 3 will be analyzed in subsequent sections.
The phenomenology of the basic LQ model generating O 8 3 will be analyzed in detail in this section, and the resulting maximum neutrino NSI will be obtained. The constraints that we derive here on the model parameters can also be applied, with some modifications, to the other O 3 models, as well as other one-loop, two-loop and three-loop LQ models discussed in subsequent sections.
To realize operator O 3b the SU (2) L doublet and singlet scalars of the Zee model [14] are replaced by SU (2) L doublet and singlet LQ fields. This model has been widely studied in the context of R-parity breaking supersymmetry, where the LQ fields are identified as the Q and d c fields of the MSSM [28,150,151]. For a non-supersymmetric description and analysis of the model, see Ref. [30].
Once the neutral component of the SM Higgs doublet acquires a VEV, the cubic term in the scalar potential (5.2) will generate mixing between the ω −1/3 and χ −1/3 fields, with the mass matrix given by: where m 2 ω and m 2 χ include the bare mass terms plus a piece of the type λv 2 arising from the SM Higgs VEV. The physical states are denoted as {X  Table III. In SUSY models with R-parity violation, ω −1/3 is identified asd and χ 1/3 asd c .
with the mixing angle given by The squared mass eigenvalues of these states are: Neutrino masses are induced via the one-loop diagram shown in Fig. 22. The mass matrix is given by: Here M d is the diagonal down-type quark mass matrix. Acceptable neutrino masses and mixings can arise in the model for a variety of parameters. Note that the induced M ν is proportional to the down-quark masses, the largest being m b . In the spirit of maximizing neutrino NSI, which are induced by either the ω −1/3 or the χ −1/3 field, without relying on their mixing, we shall adopt a scenario where the couplings λ αβ are of order one, while λ αβ 1. Such a choice would realize small neutrino masses. One could also consider λ ∼ O(1), with λ 1 as well. However, in the former case, there is a GIM-like suppression in the decay rate for α → β + γ [33], which makes the model with λ ∼ O(1), λ 1 somewhat less constrained from cLFV, and therefore we focus on this scenario. The reason for this suppression will be elaborated in Sec. 5.1.4.

Low-energy constraints
One interesting feature of the LQ model presented in this section is that the radiative decay α → β + γ is suppressed in the model due to a GIM-like cancellation. On the other hand, µ−e conversion in nuclei gives a stringent constraint on the Yukawa couplings of the model, as do the trilepton decays of the lepton to some extent. Since the product |λλ | 1 in order to generate the correct magnitude of the neutrino masses (cf. Eq. (5.7)), we shall primarily consider the case where |λ | 1 with |λ| being of order one. This is the case where the constraints from radiative decays are nonexistent. If on the other hand, |λ| 1 and |λ | is of order unity, then these radiative decays do provide significant constraints. This situation will be realized in other LQ models as well; so we present constraints on the model of this section in this limit as well. The processes that are considered are: α → β + γ, µ − e conversion in nuclei, α →¯ β γ δ (with at least two of the final state leptons being of same flavor), τ → π, τ → η, τ → η (where = e or µ), and APV.

Atomic parity violation
The strongest constraints on the λ ed and λ ed couplings come from atomic parity violation (APV) [152], analogous to the R-parity violating supersymmetric case [153]. The diagrams shown in Fig. 23 lead to the following effective couplings between up/down quarks and electrons: where we have used the Fierz transformation in the second step. The parity-violating parts of these interactions are given by On the other hand, the parity-violating SM interactions at tree-level are given by with Correspondingly, the weak charge of an atomic nucleus with Z protons and N neutrons is given by where (2Z + N ) and (Z + 2N ) are respectively the number of up and down quarks in the nucleus. The presence of the new PV couplings in Eq. (5.9) will shift the weak charge to There are precise experiments measuring APV in cesium, thallium, lead and bismuth [154]. The most precise measurement comes from cesium (at the 0.4% level [155]), so we will use this to derive constraints on LQ. For 133 55 Cs, Eq. (5.13) becomes (5.14) Taking into account the recent atomic structure calculation [152], the experimental value of the weak charge of 133 55 Cs is given by [84] Q exp Comparing this with Eq. (5.14), we obtain the corresponding 2σ bounds on λ ed and λ ed as a function of the LQ mass as follows: The APV constraint on down-quark coupling of the LQ is stronger than the up-quark coupling constraint due to the fact that the experimental value of Q w (cf. Eq. (5.15)) is 1.5σ larger than the SM prediction (cf. Eq. (5.16)), while the doublet LQ contribution to Q w goes in the opposite direction (cf. Eq. (5.14)).

µ − e conversion
Another constraint on the LQ model being discussed comes from the cLFV process of coherent µ − e conversion in nuclei (µN → eN ). We will only consider the tree-level contribution as shown in Fig. 24, since the loop-level contributions are sub-dominant. Following the general procedure described in Ref. [114], we can write down the branching ratio for this process as [33] BR where p e and E e are the momentum and energy of the outgoing electron respectively, Z and A are the atomic number and mass number of the nucleus respectively, Z eff is the effective atomic number, F p is the nuclear matrix element, and Γ N is the muon capture rate of the nucleus. Here we take | p e | E e m µ and use the values of Z eff and F p from Ref. [156], and the value of Γ N from Ref. [157]. Comparing the model predictions from Eq. (5.19) with the experimental limits for different nuclei [158][159][160], we obtain the constraints on the Yukawa couplings (either λ or λ ) and LQ mass as shown in Table XII.

5.1.3
α →¯ β γ δ decay Leptoquarks do not induce trilepton decays of the type µ → 3e at the tree-level. However, they do induce such processes at the loop level. There are LQ mediated Z and photon penguin diagrams, as well as box diagrams. These contributions have been evaluated for the LQ model of this section in Ref. [33]. With the Yukawa couplings λ being of order one, but with |λ | 1, the branching ratio for µ − → e + e − e − decay is given by [33]  The constraint on |λ ed λ µd | from the trilepton decay (cf. Eq. (5.22)) turns out to be weaker than those from µ − e conversion (cf . Table XII). Similarly, the constraints on |λ ed λ τ d | and |λ µd λ τ d | from the trilepton decay (cf. Eqs. (5.23) and (5.24)) turn out to be weaker than those from semileptonic tau decays (cf . Table XIV).

5.1.4
α → β γ constraint The lepton flavor violating radiative decay α → β +γ arises via one-loop diagrams with the exchange of LQ fields (see Fig. 25). These diagrams are analogous to Fig. 8, but with the charged and neutral scalars replaced by LQ scalars. Note that the photon can be emitted from either the LQ line, or the internal fermion line. It turns out that the LQ Yukawa coupling matrix λ leads to suppressed decay rates for α → β + γ, owing to a GIM-like cancellation. The coupling of the ω 2/3 LQ has the form αL d c βR ω 2/3 , which implies that Q B = 2/3 and Q F = −1/3 in Eq. (4.30). Consequently, the rate becomes proportional to a factor which is at most of order (m 2 b /m 2 ω ) 2 . Thus, the off-diagonal couplings of λ are unconstrained by these decays.  On the other hand, the χ −1/3 LQ field does mediate α → β + γ decays, proportional to the Yukawa coupling matrix λ . The relevant couplings have the formū L L χ , which implies that Q F = −2/3 and Q B = 1/3 in Eq. (4.30). We find the decay rate to be where 9 = 3 2 is a color factor. Here we have assumed t = m 2 F /m 2 B → 0, since the LQ is expected to be much heavier than the SM charged leptons to satisfy the experimental constraints. The limits on the products of Yukawa couplings from these decays are listed in Table XIII.

Rare D-meson decays
The coupling matrix λ of Eq. (5.1) contains, even with only diagonal entries, flavor violating couplings in the quark sector. To see this, we write the interaction terms in a basis where the down quark mass matrix is diagonal. Such a choice of basis is always available and conveniently takes care of the stringent constraints in the down-quark sector, such as from rare kaon decays. The χ leptoquark interactions with the physical quarks, in this basis, read as Here V is the CKM mixing matrix. In particular, the Lagrangian contains the following terms: − L Y ⊃ −λ αd (V ud α uχ + V cd α cχ ) + H.c. The presence of these terms will result in the rare decays D 0 → + − as well as D → π + − where = e, µ. The partial width for the decay D 0 → + − is given by Here we have used the effective Lagrangian arising from integrating out the χ field to be and the hadronic matrix element Using f D = 200 MeV, we list the constraint arising from this decay in Table XV. It will turn out that the NSI parameter ε µµ will be most constrained by the limit D 0 → µ + µ − , in cases where χ leptoquark is the mediator. Note that this limit only applies to SU (2) L singlet and triplet LQ fields, and not to the doublet LQ field Ω. The doublet LQ field always has couplings to a SU (2) L singlet quark field, which does not involve the CKM matrix, and thus has not quark flavor violation arising from V . The semileptonic decay D + → π + + − is mediated by the same effective Lagrangian as in Eq. (5.39). The hadronic matrix element is now given by term is proportional to the final state lepton mass, it can be ignored. For the form factor F + (q 2 ) we use (5.42) For the D → Dπ decay constant we use g D Dπ = 0.59 [166]. Vector meson dominance hypothesis gives very similar results [167]. With these matrix elements, the decay rate is given by    [95] and HERA [168] contact interaction bounds.
The function F is defined as Note that in the limit of infinite D mass, this function F reduces to m 6 D /24. The numerical value of the function is F 2.98 GeV 6 . Using f D = 200 MeV, f π = 130 MeV, g D Dπ = 0.59 and the experimental upper limits on the corresponding branching ratios [84], we obtain bounds on the λ couplings as shown in Table XV. These semileptonic D decays have a mild effect on the maximal allowed NSI. Note that the experimental limits on D 0 → π 0 + − are somewhat weaker than the D + decay limits and are automatically satisfied when the D + semileptonic rates are satisfied.

Contact interaction constraints
High-precision measurements of inclusive e ± p → e ± p scattering cross sections at HERA with maximum √ s = 320 GeV [168] and e + e − → qq scattering cross sections at LEP II with maximum √ s = 209 GeV [95] can be used in an effective four-fermion interaction theory to set limits on the new physics scale Λ > √ s that can be translated into a bound in the LQ mass-coupling plane. This is analogous to the LEP contact interaction bounds derived in the Zee model 4.6. Comparing the effective LQ Lagrangian (5.8) with Eq. (4.37) (for f = u, d), we see that for the doublet LQ, the only relevant chirality structure is LR, whereas for the singlet LQ, it is LL, with η d LR = η u LL = −1. The corresponding experimental bounds on Λ − and the resulting constraints on LQ mass and Yukawa coupling are given in Table XVI. In principle, one could also derive an indirect bound on LQs from the inclusive dilepton measurements at the LHC, because the LQ will give an additional t-channel contribution to the process pp → + − . However, for a TeV-scale LQ as in our case, the LHC contact interaction bounds [169,170] with √ s = 13 TeV are not applicable. Recasting the LHC dilepton searches in the fully inclusive category following Ref. [171] yields constraints weaker than those coming from direct LQ searches shown in Fig. 29.

LHC constraints
In this section, we derive the LHC constraints on the LQ mass and Yukawa couplings which will be used in the next section for NSI studies.

Pair production
At hadron colliders, LQs can be pair-produced through either gg or qq fusion, as shown in Fig. 28 (a), (b) and (c). Since LQs are charged under SU (3) c , LQ pair production at LHC is a QCD driven process, solely determined by the LQ mass and strong coupling constant, irrespective of their Yukawa couplings. Although there is a t-channel diagram [cf. Fig. (28) (c)] via charged lepton exchange through which LQ can be pair-produced via quark fusion process, this cross-section is highly suppressed compared to the s-channel pair production cross-section.
There are dedicated searches for pair production of first [172,173], second [173][174][175] and third generation [175][176][177] LQs at the LHC. Given the model Lagrangian 5.1, we are interested in the final states containing either two charged leptons and two jets ( jj), or two neutrinos and two jets (ννjj). Note that for the doublet LQ Ω = (ω 2/3 , ω −1/3 ), the jets Figure 29: LHC constraints on scalar LQ in the LQ mass and branching ratio plane. For a given channel, the branching ratio is varied from 0 to 1, without specifying the other decay modes which compensate for the missing branching ratios to add up to one. Black, red, green, blue, brown and purple solid lines represent present bounds from the pair production process at the LHC, i.e., looking for e + e − jj, µ + µ − jj, τ + τ − bb, τ + τ − tt, τ + τ − jj and ννjj signatures respectively. These limits are independent of the LQ Yukawa coupling. On the other hand, black (red) dashed, dotted and dot-dashed lines indicate the bounds on LQ mass from the single production in association with one charged lepton for LQ couplings λ ed (µd) = 2, 1.5 and 1 respectively for first (second) generation LQ.
will consist of down-type quarks, while for the singlet LQ χ −1/3 , the jets will be of up-type quarks. For the light quarks u, d, c, s, there is no distinction made in the LHC LQ searches; therefore, the same limits on the corresponding LQ masses will apply to both doublet and singlet LQs. The only difference is for the third-generation LQs, where the limit from τ + τ − bb final state is somewhat stronger than that from τ + τ − tt final state [175,177].
In Fig. 29, we have shown the LHC limits on LQ mass as a function of the corresponding branching ratios for each channel. For a given channel, the branching ratio is varied from 0 to 1, without specifying the other decay modes which compensate for the missing branching ratios to add up to one. For matter NSI, the relevant LQ couplings must involve either up or down quark. Thus, for first and second generation LQs giving rise to NSI, we can use e + e − jj and µ + µ − jj final states from LQ pair-production at LHC to impose stringent bounds on the λ αd and λ αd couplings (with α = e, µ) which are relevant for NSI involving electron and muon flavors. There is no dedicated search for LQs in the τ + τ − jj channel to impose similar constraints on λ τ d and λ τ d relevant for tau-flavor NSI. There are searches for third generation LQ [176,177] looking at τ + τ − bb and τ + τ − tt signatures which are not relevant for NSI, since we do not require λ τ t (for χ −1/3 ) or λ τ b (for ω 2/3 ) couplings. For constraints on λ τ d , we recast the τ + τ − bb search limits [175][176][177] taking into account the b-jet misidentification as light jets, with an average rate of 1.5% (for a b-tagging efficiency of 70%) [178]. As expected, this bound is much weaker, as shown in Fig. 29.
However, a stronger bound on NSI involving the tau-sector comes from ννjj final state. From the Lagrangian (5.1), we see that the same λ τ d coupling that leads to τ + τ − dd final state from the pair-production of ω 2/3 also leads to ν τντ dd final state from the pairproduction of the SU (2) L partner LQ ω −1/3 , whose mass cannot be very different from that of ω 2/3 due to electroweak precision data constraints (similar to the Zee model case, cf. Sec. 4.4). Since the final state neutrino flavors are indistinguishable at the LHC, the ννjj constraint will equally apply to all λ αd (with α = e, µ, τ ) couplings which ultimately restrict the strength of tau-sector NSI, as we will see in the next subsection. The same applies to the λ τ d couplings of the singlet LQ χ −1/3 , which are also restricted by the ννjj constraint.

Single production
LQs can also be singly produced at the collider in association with charged leptons via sand t-channel quark-gluon fusion processes, as shown in Fig. 28 (d) and (e). The single production limits, like the indirect low-energy constraints, are necessarily in the masscoupling plane. This signature is applicable to LQs of all generations. In Fig. 29, we have shown the collider constraints in the single-production channel for some benchmark values of the first and second generation LQ couplings λ ed and λ µd (since d jets cannot be distinguished from s jets) equal to 1, 1.5 and 2 by dot-dashed, dotted and dashed curves respectively. The single-production limits are more stringent than the pair-production limits only for large λ ed , but not for λ µd . There is no constraint in the τ j channel, and the derived constraint from τ b channel is too weak to appear in this plot.

How light can the leptoquark be?
There is a way to relax the ννjj constraint and allow for smaller LQ masses for the doublet components. This is due to a new decay channel ω −1/3 → ω 2/3 +W − which, if kinematically allowed, can be used to suppress the branching ratio of ω −1/3 → νd decay for relatively smaller values of λ αd couplings, thereby reducing the impact of the ννjj constraint. The partial decay widths for ω −1/3 → ω 2/3 + W − and ω −1/3 → ν α d β are respectively given by In deriving Eq. (5.44), we have used the Goldstone boson equivalence theorem, and in Eq. (5.45), the factor in the denominator is not 8π (unlike the SM h → bb case, for instance), because only one helicity state contributes. The lighter LQ ω 2/3 in this case can only decay to α d β with 100% branching ratio. Using the fact that constraints from τ + τ − jj channel are weaker, one can allow for ω 2/3 as low as 522 GeV, as shown in Fig. 29 by the solid brown curve, when considering the λ τ d coupling alone. This is, however, not applicable to the scenario when either λ ed or λ µd coupling is present, because of the severe constraints from e + e − jj and µ + µ − jj final states.

NSI prediction
The LQs ω −1/3 and χ −1/3 in the model have couplings with neutrinos and down-quark (cf. Eq. (5.1)), and therefore, induce NSI at tree level as shown in Fig. 30 via either λ or λ couplings. From Fig. 30, we can write down the effective four-fermion Lagrangian as where we have used Fierz transformation in the second step. Comparing Eq. (5.46) with Eq. (3.1), we obtain the NSI parameters Np(x) = 1, one can obtain the effective NSI parameters from Eq. (3.5) as To satisfy the neutrino mass constraint [cf. Eq. (5.7)], we can have either λ or λ couplings of O(1), but not both simultaneously. As mentioned in Sec. 5.1, the choice λ 1 and λ ∼ O(1) is less constrained from cLFV.

Doublet leptoquark
First, taking the λ-couplings only and ignoring the λ contributions, we show in Figs. 31 and 32 the predictions for diagonal (ε ee , ε µµ , ε τ τ ) and off-diagonal (ε eµ , ε µτ , ε eτ ) NSI parameters respectively from Eq. (5.48) by black dotted contours. Colored shaded regions in each plot are excluded by various theoretical and experimental constraints. In Figs. 31 (b) and (c), the yellow colored regions are excluded by perturbativity constraint, which requires the LQ coupling λ αd < 4π √ 3 [179]. Red shaded region in Fig. 31 (a) is excluded by the APV bound (cf. Sec. 5.1.1), while the brown and cyan regions are excluded by HERA and LEP contact interaction bounds, respectively (cf. Table XVI). Red shaded region in Fig. 31 (c) is excluded by the global fit constraint from neutrino oscillation+COHERENT data [60]. Blue shaded regions in Figs. 31 (a) and (b) are excluded by LHC LQ searches (cf. Fig. 29) in the pair-production mode for small λ αd (which is independent of λ αd ) and single-production mode for large λ αd ) with α = e, µ. Here we have assumed 50% branching ratio to ej or µj, and the other 50% to τ d in order to relax the LHC constraints and allow for larger NSI. Blue shaded region in Fig. 31 (c) is excluded by the LHC constraint from the ννjj channel, where the vertical dashed line indicates the limit assuming BR(ω −1/3 → νd) = 100%, and the unshaded region to the left of this line for small λ τ d is allowed by opening up the ω −1/3 → ω 2/3 W − channel (cf. Sec. 5.3.3). Note that we cannot completely switch off the ω −1/3 → νd channel, because that would require λ τ d → 0 and in this limit, the NSI will also vanish.
The red line in Fig. 31 (b) is the suggestive limit on ε dR αβ from NuTeV data [140] (cf. Table XVII). This is not shaded because there is a 2.7σ discrepancy of their s 2 w measurement with the PDG average [84] and a possible resolution of this might affect the NSI constraint obtained from the same data. Here we have rederived the NuTeV limit following Ref. [140], but using the latest value of s 2 w (on-shell) [84] (without including NuTeV). Specifically, we have used the NuTeV measurement of the effective coupling g µ R 2 = 0.0310 ± 0.0011 from ν µ q → νq scatterings [180] which is consistent with the SM prediction of g µ R 2 SM = 0.0297. Here g µ R 2 is defined as where g u R = − 2 3 s 2 w and g d R = 1 3 s 2 w are the Z couplings to right-handed up and down quarks respectively. Only the right-handed couplings are relevant here, since the effective NSI Lagrangian (5.46) involves right-handed down-quarks for the doublet LQ component ω 2/3 . In Eq. (5.49), setting ε uR µµ = 0 for this LQ model and comparing g µ R 2 with the measured value, we obtain a 90% CL on ε dR µµ < 0.029, which should be multiplied by 3 (since ε αβ ≡ 3ε dR αβ ) to get the desired constraint on ε αβ shown in Fig. 31 (b). Figure 31: Predictions for diagonal NSI (ε ee , ε µµ , ε τ τ ) induced by doublet LQ in the oneloop LQ model are shown by black dotted contours. Colored shaded regions are excluded by various theoretical and experimental constraints. Yellow colored region is excluded by perturbativity constraint on LQ coupling λ αd [179]. Blue shaded region is excluded by LHC LQ searches (Fig. 29) in subfigure (a) by e+jets channel (pair production for small λ ed and single-production for large λ ed ), in subfigure (b) by µ+jets channel, and in subfigure (c) by ν+jet channel. In (a), the red, brown and cyan shaded regions are excluded by the APV bound (cf. Eq. 5.18), HERA and LEP contact interaction bounds (cf. Table XVI) respectively. In (b), the red line is the suggestive limit from NuTeV [140]. In (c), the red shaded region is excluded by the global fit constraint from neutrino oscilla-tion+COHERENT data [60]. We also show the future DUNE sensitivity in blue solid lines for both 300 kt.MW.yr and 850 kt.MW.yr [65].  Fig. 29). In (a) and (b), the brown and green shaded regions are excluded by τ → π 0 and τ → η (with = e, µ) constraints (cf . Table XIV). In (a), the red shaded region is excluded by the global fit constraint on NSI from neutrino oscilla-tion+COHERENT data [60]. In (b), the yellow shaded region is excluded by perturbativity constraint on LQ coupling λ αd [179] combined with APV constraint (cf. Eq. (5.18)). In (c), the red shaded region is excluded by µ → e conversion constraint. Also shown in (b) are the future DUNE sensitivity in blue solid lines for both 300 kt.MW.yr and 850 kt.MW.yr [65]. (c) (d) Figure 33: Additional low-energy constraints on NSI induced by singlet LQ. Subfigure (a) has the same APV and LHC constraints as in Fig. 18 (a), the modified HERA and LEP contact interaction bounds (cf . Table XVI), plus the D + → π + e + e − constraint, shown by green shaded region (cf. Sec. 5.1.6). Subfigure (b) has the same constraints as in Fig. 18 (b), plus the D + → π + µ + µ − constraint, shown by light-green shaded region, and D 0 → µ + µ − constraint shown by brown shaded region (cf. Sec. 5.1.6). Subfigure (c) has the same constraints as in Fig. 19 (a), plus the τ → µγ constraint, shown by purple shaded region. Subfigure (d) has the same constraints as in Fig. 19 (b), plus the τ → eγ constraint, shown by purple shaded region.  2), perturbative unitarity and collider (Secs. 5.3) constraints. We also impose the constraints from neutrino-nucleon scattering experiments, like CHARM II [140], NuTeV [140], COHERENT [181] and IceCube [182], as well as the global fit constraints from neutrino oscillation+COHERENT data [60], whichever is stronger. The scattering and global fit constraints are on ε d αβ , so it has been scaled by a factor of 3 for the constraint on ε αβ in the Table. The maximum allowed value for each NSI parameter is obtained after scanning over the LQ mass (see Figs. 31 and 32) and the combination of the relevant constraints limiting the NSI are shown in parentheses in the second column. The same numbers are applicable for the doublet and singlet LQ exchange, except for ε ee where the APV constraint is weaker than HERA (Fig. 33 (a))) and for ε µµ which has an additional constraint from D + → π + µ + µ − decay (see Fig. 33 (b)). In the last column, we also show the future DUNE sensitivity [65] for 300 kt.MW.yr exposure (and 850 kt.MW.yr in parentheses).
For ε ee , the most stringent constraint comes from APV (Sec. 5.1.1), as shown by the red shaded region in Fig. 31 (a) which, when combined with the LHC constraints on the mass of LQ, rules out the possibility of any observable NSI in this sector. Similarly, for ε µµ , the most stringent limit comes from NuTeV. However, if this constraint is not considered, ε µµ can be as large as 21.6%. On the other hand, ε τ τ can be as large as 34.3%, constrained only by the LHC constraint on the LQ mass and perturbative unitarity constraint on the Yukawa coupling (cf. Fig. 31 (c)). This is within the future DUNE sensitivity reach, at least for the 850 kt.MW.yr (if not 300 kt.MW.yr) exposure [65], as shown in Fig. 31 (c). Note that from oscillation data alone, ε τ τ − ε µµ is constrained to be less than 9.5% [60].
As for the off-diagonal NSI in Fig. 19, the LHC constraints (cf. Sec. 5.3) are again shown by blue shaded regions. The yellow shaded region in Fig. 19 (b) is from the combination of APV and perturbative unitarity constraints. However, the most stringent limits for all the off-diagonal NSI come from cLFV processes. In particular, τ → π 0 and τ → η (with = e, µ) impose strong constraints (cf. Sec. 5.1.5) on ε µτ and ε eτ , as shown in Figs. 32 (a) and (b). For ε eµ , the most stringent limit comes from µ−e conversion (cf. Sec. 5.1.2), as shown in Fig. 32 (c). The maximum allowed NSI in each case is tabulated in Table XVII, along with the current constraints from neutrino-nucleon scattering experiments, like CHARM [140], COHERENT [181] and IceCube [182], as well as the global fit constraints from neutrino oscillation+COHERENT data [60] and future DUNE sensitivity [65]. It turns out that the cLFV constraints have essentially ruled out the prospects of observing any off-diagonal NSI in this LQ model in future neutrino experiments. This is consistent with general arguments based on SU (2) L gauge-invariance [20].

Singlet leptoquark
Now if we take the λ couplings instead of λ in Eq. (5.48), the NSI predictions, as well as the constraints, can be analyzed in a similar way as in Figs. 31 and 32. Here the APV (cf. Eq. (5.18)), as well as the LEP and HERA contact interaction constraints on ε ee (cf. Table XVI) are somewhat modified. In addition, there are new constraints from D + → π + + − and D 0 → + − (cf. Sec. 5.1.6) for ε ee and ε µµ , as shown in Fig. 33 (a) and (b). For ε ee , the D + → π + e + e − constraint turns out to be much weaker than the APV constraint. The D 0 → e + e − constraint is even weaker and does not appear in Fig. 33 (a). However, for ε µµ , the D + → π + µ + µ − constraint turns out to be the strongest, limiting the maximum allowed value of ε µµ to a mere 0.8%, as shown in Fig. 33 (b) and in Table XVII.
The NuTeV constraint also becomes more stringent here due to the fact that the singlet LQ χ couples to left-handed quarks (cf. Eq. (5.46)). So it will affect the effective coupling g L . For ε µµ , we use the NuTeV measurement of g µ L 2 = 0.3005 ± 0.0014 from ν µ q → νq scatterings [180] which is 2.7σ smaller than the SM prediction of g µ L 2 SM = 0.3043. Here g µ L 2 is defined as where g u L = 1 2 − 2 3 s 2 w and g d L = − 1 2 + 1 3 s 2 w . For the SM prediction, we have used the latest PDG value for on-shell s 2 w = 0.22343 from a global fit to electroweak data (without NuTeV) [84] and comparing g µ L 2 with the measured value, derive a 90% CL constraint on 0.0018 < ε µµ < 0.8493. Note that this prefers a non-zero ε µµ at 90% CL (1.64σ) because the SM with ε µµ = 0 is 2.7σ away and also because there is a cancellation between g d L (which is negative) and ε µµ (which is positive) in Eq. (5.50) to lower the value of g µ L 2 to within 1.64σ of the measured value.
For the off-diagonal sector, there are new constraints from τ → γ relevant for ε µτ and ε eτ , as shown in Figs. 33 (c) and (d). However, these are less stringent than the τ → π 0 and τ → η constraints discussed before. There are no new constraints for ε τ τ and ε eµ that are stronger than those shown in Figs. 31 (c) and 32 (c) respectively, so we do not repeat these plots again in Fig. 33.  Table III [31].

NSI in a triplet leptoquark model
Ω 3, 2, 1 6 = ω 2/3 , ω −1/3 . The relevant Lagrangian for the neutrino mass generation can be written as These interactions, along with the potential term whereρ is related to ρ by charge conjugation as ρ 3, 3, − 1 3 = ρ 2/3 , −ρ −1/3 , ρ −4/3 , induce neutrino mass at one-loop level via the O 9 3 operator in the notation of Ref. [31], as shown in Fig. 34. The neutrino mass matrix can be estimated as where M d is the diagonal down-type quark mass matrix and M ≡ max(m ω , m ρ ). The NSI parameters read as Note that both λ and λ cannot be large at the same time due to neutrino mass constraints (cf. Eq. (6.3)). For λ λ , this expression is exactly the same as the doublet LQ contribution derived in Eq. (5.48) and the corresponding maximum NSI can be read off from Table XVII for the doublet component. On the other hand, for λ λ, the third term in Eq. (6.4) is analogous to the downquark induced singlet LQ NSI given in Eq. (5.48) (except for the Clebsch-Gordan factor of (1/ √ 2) 2 ), whereas the second term is a new contribution from the up-quark sector. Note that both terms depend on the same Yukawa coupling λ αu = λ αd in the Lagrangian (6.1). This is unique to the triplet LQ model, where neutrinos can have sizable couplings to both up and down quarks simultaneously, without being in conflict with the neutrino mass constraint. As a result, some of the experimental constraints quoted in Sec. 5 which assumed the presence of only down-quark couplings of LQ will be modified in the triplet case, as discussed below:

Atomic parity violation
The shift in the weak charge given by Eq. (5.13) is modified to Assuming m ρ 1/3 = m ρ 4/3 ≡ m ρ and noting that λ αu = λ αd in Eq. (6.1), we obtain Comparing this with the 2σ allowed range (5.17), we obtain the modified constraint which is weaker (stronger) than that given by Eq. (5.18) for the SU (2) L -doublet (singlet) LQ alone.

µ − e conversion
From Eq. (5.19), we see that for the triplet case, the rate of µ − e conversion will be given by For degenerate ρ-mass and λ d = λ u , we obtain the rate to be (3/2) 2 times larger than that given in Eq. (5.19). Therefore, the constraints on |λ ed λ µd | given in Table XII will be a factor of 3/2 stronger.

D-meson decays
The α u βρ 1/3 and α d βρ 4/3 terms in Eq. (7.22) induce flavor violating quark decays. Following the discussion in Sec. 5.1.6, we work in a basis where the down quark mass matrix is diagonal, so there are no constraints from rare kaon decays. However, the α u βρ 1/3 term in Eq. (7.22) now becomes α V id u iρ 1/3 which induces D 0 → + − and D + → π + + − decays. The analysis will be the same as in Sec. 5.1.6, except that the λ αd couplings will now be replaced by λ αd / √ 2. Correspondingly, the constraints on |λ αd | given in Table XV will be √ 2 times weaker. For instance, (6.13)

Contact interaction constraints
The LEP and HERA contact interaction bounds discussed in Sec. 5.2 will also be modified in the triplet LQ case. Here, the interactions are only of LL type, but the effective Yukawa coupling is 3/2 times that of the singlet case in Table XVI. The modified constraint is given by .904 TeV from LEP 3.127 TeV from HERA . (6.14)

LHC constraints
The LHC constraints on theρ fields will be similar to the discussion in Sec. 5.3. Comparing the Lagrangians (5.1) and (7.22), we see thatρ 1/3 will have the same decay modes to νj and j, and therefore, the same constraints as the singlet χ −1/3 discussed in Sec. 5.4.2. In our analysis, we have assumed degenerate mass spectrum for all the triplet LQ fields. But we note here that theρ −2/3 component can in principle be lighter, since it can only decay to νj for which the constraints are weaker (cf. Fig. 29). However, the mass splitting between ρ −2/3 andρ 1/3 cannot be more than ∼ 100 GeV from T -parameter constraints, analogous to the charged scalar case discussed in Sec. 4.4 (cf. Fig. 7). In that case, the limit on m ρ 1/3 for 50% branching ratio to νj and j channels (since they are governed by the same λ αd coupling), one can allow for m ρ −2/3 as low as 800 GeV or so.

NSI prediction
Taking into account all the constraints listed above, we show in Figs. 35 and 36 the predictions for diagonal (ε ee , ε µµ , ε τ τ ) and off-diagonal (ε eµ , ε µτ , ε eτ ) NSI parameters respectively from Eq. (6.4) by black dotted contours. Colored shaded regions in each plot are excluded by various theoretical and experimental constraints, as in Figs. 31 and 32. The main difference is in the NuTeV constraint shown in Fig. 35 (b), which is more stringent than those shown in Figs. 31 (b) and 33 (b). The reason is that in presence of both ε uL µµ and ε dL µµ as in this LQ model (cf. (6.1)), the total contribution to g µ L 2 in Eq. (5.50) is always positive, and therefore, any nonzero ε µµ will make the discrepancy worse than the SM case of 2.7σ. Therefore, we cannot impose a 90% CL (1.64σ) constraint from NuTeV in this scenario. The line shown in Fig. 35 (b) corresponds to the 3σ constraint on ε µµ < 0.0007, which is subject to the same criticism as the discrepancy with the SM, and therefore, we have not shaded the NuTeV exclusion region and do not consider it while quoting the maximum allowed NSI.
From Figs This is also summarized in Fig. 58 and in Table XX. 7 Other type-I radiative models In this section, we briefly discuss the NSI predictions in other type-I radiative models at one-, two-and three-loops. In each case, we present the new particle content, model Lagrangian, Feynman diagrams for neutrino mass generation and expressions for neutrino mass, followed by the expression for NSI parameters. The maximum NSI allowed in each model is summarized in Table XX. 7.1 One-loop models

Minimal radiative inverse seesaw model
This is an exception to the general class of type-I radiative models, where the new particles running in the loop will always involve a scalar boson. In this model, the SM Higgs and Z bosons are the mediators, with the new particles being SM-singlet fermions. 14 The low-energy effective operator that leads to neutrino mass in this model is the dimension-7 operator However, this mechanism is only relevant when the dimension-5 operator given by Eq. (1.1) that leads to the tree-level neutrino mass through the seesaw mechanism is forbidden due to some symmetry. This happens in the minimal radiative inverse seesaw model [43]. In the usual inverse seesaw model [183], one adds two sets of SM-singlet fermions, N and S, with opposite lepton numbers. The presence of a Majorana mass term for the S-field, i.e., µ S SS leads to a tree-level neutrino mass via the standard inverse seesaw mechanism [183]. However, if one imposes a global U (1) symmetry under which the S-field is charged, then the breaking term that is allowed is the Majorana mass term for the N -field, i.e., µ R N N . It can be shown that this term by itself does not give rise to neutrino mass at tree-level, but a non-zero neutrino mass is inevitably induced at one-loop through the diagram shown in Fig. 37 involving the SM Higgs doublet (which gives rise to two diagrams involving the SM Higgs and Z-boson after electroweak symmetry breaking [43]). One can see that the N σ N in Eq. (7.2), but not the Majorana mass terms Sσ ( ) S.
After electroweak symmetry breaking, evaluating the self-energy diagrams that involve the Z-boson and Higgs boson (cf. Fig. 37), the neutrino mass reads as (in the limit µ R M N ) [43,184]: H /m 2 W and x Z = m 2 Z /m 2 W , and we have assumed M N = m N 1 for simplicity.
The NSI in this model arise due to the fact that the light SU (2) L -doublet neutrinos ν mix with the singlet fermions N and S, due to which the 3 × 3 lepton mixing matrix is no longer unitary. The neutrino-nucleon and neutrino-electron interactions proceed as in the SM via t-channel exchange of W and Z bosons, but now with modified strength because of the non-unitarity effect, that leads to NSI [185]. If only one extra Dirac state mixes with the three light states with mixing parameters U α4 (with α = e, µ, τ ), we can write the NSI parameters as where Y n = N n /N e is the ratio of the average number density of neutrons and electrons in matter. Note that for Y n → 1 which is approximately true for neutrino propagation in earth matter, we get vanishing ε eµ and ε eτ up to second order in U α4 . 16 Taking into account all  Table II [31].
the experimental constraints on U α4 U β4 from neutrino oscillation data in the averaged-out regimes, beta decay, rare meson decay, beam dump experiments, cLFV searches, collider constraints from LEP and LHC, as well as electroweak precision constraints [185][186][187][188][189][190], the maximum NSI parameters allowed in this model are summarized in Table XX. We find that For ε eµ and ε eτ , we have used Y n = 1.051 (for average value all over the earth) in Eq. (7.4), in addition to the cLFV constraints on U e4 U µ4 and U e4 U τ 4 . The maximum NSI values listed above (and also summarized in Table XX) are obtained for a relatively light (∼MeV-scale) sterile neutrino, where the experimental constraints are weaker than at higher masses. The NSI expressions (7.4) also apply to two-loop radiative models with two W -boson exchange [191][192][193]. However, the maximum NSI obtainable in these models will be much smaller than the estimate in Eq. (7.5) because the sterile neutrino in this case is required to be heavier for successful neutrino mass generation at two-loop.

One-loop model with vectorlike leptons
This model [31] utilizes the same d = 7 operator O 2 = L i L j L k e c H l ij kl (cf. Eq. (2.2b)), as in the Zee model to generate a one-loop neutrino mass. The new particles added are a scalar singlet η + (1, 1, 1) and a vectorlike lepton ψ 1, 2, − 3 2 = (E, F −− ), which give rise to the O 1 2 operator L(LL)(e c H) (cf. Table II). Neutrino mass is generated via the one-loop diagram shown in Fig. 38. The relevant Lagrangian for the neutrino mass generation reads: where ψ c = (F ++ , −E c ) and H 1, 2, 1 2 is the SM Higgs doublet. Expanding the first two terms, we get   [194]. Here g exp α stands for the effective gauge coupling extracted from muon and tau decays in the different leptonic channels.
The neutrino mass matrix can be estimated as where M is the diagonal mass matrix for the SM charged leptons, M E is the diagonal mass matrix for the vector-like leptons with eigenvalues m E i , and M ≡ max(m η , m E i ). Note that just one flavor of ψ is not sufficient, because in this case, the neutrino mass matrix (7.8) would have a flavor structure given by (f M − M f ), which has all the diagonal entries zero, similar to the Zee-Wolfenstein model [76]. Such a structure is ruled out by observed neutrino oscillation data. Thus, we require at least two flavors of ψ, in which case the diagonal entries of M ν are nonzero, and the model is consistent with experiments. NSI in this model are induced by the f -type couplings in Eq. (7.7), similar to the f -couplings in the Zee model Lagrangian (4.2). The NSI parameters read as Due to the antisymmetric nature of the f couplings, the only relevent NSI parameters in this case are ε µτ , ε µµ , and ε τ τ . These are severely constrained by cLFV searches and universality of charged currents [194], as shown in Table XIX. This is similar to the case of Zee-Babu model discussed later in Sec. 7.2.1. Since the singly-charged scalar mass has to be above ∼ 100 GeV to satisfy the LEP constraints (cf. Sec. 4.7), we obtain from Eq. (7.9) and Table XIX the following maximum values: This is also summarized in Table XX.   Table III). The new particles introduced are a scalar LQ singlet χ 3, 1, − 1 3 and a vectorlike quark doublet Q 3, 2, − 5 6 = D −1/3 , X −4/3 . Neutrino mass is generated at one-loop level as shown in Fig. 39. The QQχ and d c u c χ interaction terms, allowed by gauge invariance, are forbidden by demanding baryon-number conservation in order to avoid rapid proton decay. The relevant Lagrangian for the neutrino mass generation reads as Expanding the first two terms, we get The neutrino mass matrix can be estimated as where M d is the diagonal down-type quark mass matrix, M D is the mass matrix for the down-type VQ with eigenvalues m D i , and M ≡ max(m χ , m D i ). With a single copy of VQ quarks, the rank of M ν is two, implying that the lightest neutrino has zero mass at the one-loop order. This model can lead to consistent neutrino oscillation phenomenology. NSI in this model are induced by the λ-type interactions in Eq. (7.12): This is exactly same as the singlet LQ contribution in Eq. (5.48) and the corresponding maximum NSI can be read off from This is also summarized in Table XX. This is referred to as O 6 3 in Table III. The model has an SU (2) L -doublet LQ Ω 3, 2, 1 6 = ω 2/3 , ω −1/3 and an SU (2) L -triplet vectorlike quark Σ 3, 3, 2 3 = Y 5/3 , U 2/3 , D −1/3 . Neutrino mass is generated at one-loop level via the Feynman diagram shown in Fig. 40. The relevant Lagrangian for the neutrino mass generation can be written as where Ω = iτ 2 Ω is the isospin conjugate field. Expanding the terms in Eq. (7.16), we obtain The neutrino mass can be estimated as where M d and M D are the diagonal down quark mass matrix and vectorlike quark mass matrix respectively, and M ≡ max(m ω , m D i ), with m D i being the eigenvalues of M D . As in previous models with one copy of vectorlike fermion, the rank of M ν is two in this model, implying that the lightest neutrino is massless at the one-loop level. NSI in this model are induced by the doublet LQ component ω −1/3 . The NSI parameters read as This expression is exactly the same as the doublet LQ contribution in Eq. (5.48) and the corresponding maximum NSI can be read off from This is also summarized in Table XX.

Model with SU (2) L -triplet leptoquark and vectorlike quark
This is based on the operator O 5 3 (see Table III) which is realized by adding an SU (2) Ltripletρ 3 , 3, 1 3 = ρ 4/3 ,ρ 1/3 ,ρ −2/3 and a vectorlike quark doublet Q 3, 2, − 5 6 = D −1/3 , X −4/3 . Neutrino mass is generated at one-loop level, as shown as Fig. 41. There is also a two-loop diagram involving ρ 2/3 , which is not considered here, as that would be sub-dominant to the one-loop diagram. The interaction term QQρ is forbidden by demanding baryon-number conservation to avoid proton decay. The relevant Lagrangian for the neutrino mass generation can be written as whereρ is related to ρ by charge conjugation as ρ 3, 3, − 1 3 = ρ 2/3 , −ρ −1/3 , ρ −4/3 . Expanding the terms in Eq. (7.21), we get The neutrino mass can be estimated as where M d and M D are the diagonal mass matrices for down-type quark and vectorlike quark fields, and M = max(m D i , m ρ ), with m D i being the eigenvalues of M D . With a single copy of the vectorlike quark, the matrices y and λ are 3 × 1 dimensional. Consequently the rank of M ν is two, which would imply that the lightest neutrino mass m 1 = 0 at the one-loop level. Realistic neutrino mixing can however be generated, analogous to the model of Ref. [15,16]. NSI in this model are induced by bothρ −2/3 andρ 1/3 fields, which couple to up and down quarks respectively (cf. Eq. (7.22)). The NSI parameters read as This is same as the triplet contribution in Eq. (6.4) and the maximum allowed values are given in Eq. (6.15).

A new extended one-loop leptoquark model
Here we present a variation of the one-loop LQ model of Sec. 5 wherein the neutrino mass is generated with up-quark chiral suppression (see Fig. 42), rather than down-quark mass suppression (as in Fig. 22). The effective operator of the model is of dimension nine, given by O 1 = (LQ)(Lu c )(HH)H , (7.25) which may appear to be a product of O 1 of Eq. (1.1) and the SM operator (Qu c H); but the SU (2) L contractions mix the two sub-operators. To realize this operator at the one-loop level, three SU (3) c -triplet LQ fields are introduced: δ 3, 2, 7 6 = δ 5/3 , δ 2/3 ,ρ 3 , 3, 1 3 = ρ 4/3 ,ρ 1/3 ,ρ −2/3 , ξ 3, 1, 2 3 . Since three new fields are introduced, this model may be viewed as non-minimal, and does not fit into the classification of The corresponding Lagrangian for the neutrino mass generation reads as Neutrino mass is generated by the diagram shown in Fig. 42 where m 1 and m 2 are the masses of the heaviest two LQs among the δ,ρ and ξ fields, and M u is the diagonal mass matrix in the up-quark sector. To get small neutrino masses, we need the product λλ 1. We may take λ ∼ O(1) and λ λ which is preferable to the other case of λ λ , since the λ couplings are constrained by D-meson decays (see Sec. 5.1.6).
After integrating out the heavy LQ fields, Eq. (7.26) leads to an effective NSI Lagrangian with up-quarks in the neutrino propagation through matter. The NSI parameters read as For λ λ , this expression is exactly the same as the doublet LQ contribution derived in Eq. (5.48) and the corresponding maximum NSI can be read off from There are other variations of one-loop LQ models with more exotic particles [29,30], where the neutrino mass is proportional to up-type quark mass. The NSI predictions in these models are the same as in Eq. (7.29).

Zee-Babu model
This model realizes the operator O 9 of Eq. (1.4). In this model [15,16], two SU (2) L -singlet Higgs fields, h + (1, 1, 1) and k ++ (1, 1, 2), are introduced. The corresponding Lagrangian for the generation of neutrino mass reads: Majorana neutrino masses are induced at two-loop as shown in Fig. 43 by the Lagrangian (7.30), together with the potential term The neutrino mass matrix reads: where M = max(m k ++ , m h + ) and I is a dimensionless function that depends on the ratio of the masses of the two new scalars [89,[194][195][196]. The singly charged scalar h + induces NSI Figure 43: Neutrino mass generation at two-loop in the Zee-Babu model [15,16]. This model generates operator O 9 of Eq. (1.4).
at tree-level through the f -type Yukawa coupling in Eq. (7.30). After integrating out the heavy scalars, NSI induced in neutrino propagation through normal matter can be written as This is exactly the same as Eq. (7.9) for which the maximum NSI are given by Eq. (7.10). These are severely constrained by cLFV searches and universality of charged currents [194] (cf . Table XIX), restricting the maximum NSI to O(10 −3 ) level [197]. These numbers are summarized in Table XX.

Leptoquark/diquark variant of the Zee-Babu model
One can also generate neutrino mass at two-loop by replacing leptons with quarks in the Zee-Babu model as shown in Fig. 44. In addition to the SM fields, this model [32] employs a scalar LQ χ 3, 1, − 1 3 and a scalar diquark ∆ 6, 1, − 2 3 . The χ (∆) field plays the role of singly (doubly)-charged scalar in the Zee-Babu model. The relevant Yukawa Lagrangian for the neutrino mass generation is written as Neutrino mass is generated at two-loop via the Lagrangian (7.34) in combination with the potential term The neutrino mass matrix can be calculated as where M ≡ max(m χ , m ∆ ), M d is the diagonal down-type quark mass matrix, and I is a dimensionless two-loop integral defined in terms of the ratio of m 2 ∆ and m 2 χ [89]. After integrating out the heavy scalars, the NSI parameters in this model are given by in Fig. 44. The cubic term χ χ ∆ will not be allowed by Bose symmetry in this case. By assuming two copies of the χ field, namely, χ 1 and χ 2 , one could restore this coupling from χ 1 χ 2 ∆, in which case the diagram of Fig. 44 can be connected [41]. The NSI in such a model is identical to the model described in this section. Second, one could replace the internal down quarks of Fig. 44 by up-type quarks, with a simultaneous replacement of χ 3, 1, − 1 3 by ρ 3, 3, − 1 3 and ∆ 6, 1, − 2 3 by ∆ 6, 1, 4 3 . Neutrino NSI will then follow the ρ NSI predictions as in Sec. 7.1.5. In this up-quark variant, one could replace the diquark ∆ 6, 1, 4 3 by a color triplet field ∆ 3, 1, 4 3 as well [41]. Figure 47: Two-loop neutrino mass generation in the Angelic model [35]. This model induces operator O 11 of Ref. [26].
In this model [35], one adds two scalar LQs χ a 3, 1, − 1 3 (with a = 1, 2) and a color-octet Majorana fermion F (8, 1, 0). The relevant Yukawa Lagrangian is written as Expanding the first term, we get Within this framework, neutrino mass is induced at two-loop level as shown in Fig. 47 which can be estimated as   Table V [31].
The neutrino mass can be estimated as  Tables III and IV). Therefore, we do not discuss the NSI prospects in these models.
Neutrino mass is generated at two-loop level as shown in Fig. 50 and can be estimated as where M d , M u , M and M U are the diagonal mass matrices for down quark, up quark, charged leptons and vectorlike quarks, respectively, and m U is the largest eigenvalue of M U . The NSI parameters can be written as This is exactly the same expression as the doublet contribution in Eq. (5.46), with the maximum values given in Table XVII.
The neutrino mass is induced at two-loop level as shown in Fig. 51 and can be estimated as where M u is the diagonal up-type quark mass matrix. Note that M ν is a symmetric matirx, as it should be, since h = h T . After integrating out the heavy scalars, NSI induced in this model can be written as This is same as the extended one-loop LQ model prediction in Eq. (7.29) for λ λ . The maximum allowed values can be read off from Table XVII for the doublet component. This is also summarized in Table XX. 7.3 Three-loop models

KNT Model
The Krauss-Nasri-Trodden (KNT) model [36] generates the d = 9 operator O 9 of Eq. (1.4). SM-singlet fermions N α (1, 1, 0) and two SM-singlet scalars η + 1 and η + 2 with SM charges (1, 1, 1) are introduced. The relevant Yukawa Lagrangian is written as Tree level mass is prevented by imposing a Z 2 symmetry under which the fields η + 2 and N are odd, while the other fields are even. The Majorana mass term for N as shown in Eq. (7.64) explicitly breaks lepton number. Neutrino masses are generated at three-loop as shown in Fig. 52 by the Lagrangian (7.64), together with the quartic term in the potential The estimated neutrino mass matrix reads as where M is the diagonal charged lepton mass matrix, M N = diag(m Nα ) is the diagonal Majorana mass matrix for N α fermions, M ≡ max(m Nα , m η 1 , m η 2 ), and I is a three-loop function obtained in general by numerical integration [199]. Figure 53: Three-loop neutrino mass generation in the AKS model [38]. The model induces operator O 3 of Eq. (7.68).
NSI in the KNT model arise from singly-charged scalar η + 1 that has the same structure as in the Zee-Babu model (cf. Sec. 7.2.1) and are given by The maximum NSI one can get in this model are same as in Eq. (7.10) and also summarized in Table XX.

AKS model
In the Aoki-Kanemura-Seto (AKS) model [38] an effective ∆L = 2 operator of dimension 11 is induced: Note that there is a chiral suppression in this model unlike generic operators of type O 1 given in Eq. (1.5). In addition to the SM fields, the following particles are added: an isospin doublet scalar Φ 2 1, 2, 1 2 , a singly-charged scalar singlet η + (1, 1, 1), a real scalar singlet η 0 (1, 1, 0), and two isospin-singlet right-handed neutrinos N α (1, 1, 0) (with α = 1, 2). The relevant Yukawa Lagrangian for the neutrino mass generation reads as where Φ 1 1, 2, 1 2 is the SM Higgs doublet. Tree-level neutrino mass is forbidden by imposing a Z 2 symmetry under which η ± , η 0 and N αR are odd, while the remaining fields are even. Neutrino masses are generated at three-loop, as shown in Fig. 53, by combining Eq. (7.69) with the quartic term in the potential where tan β ≡ Φ 0 2 / Φ 0 1 and I is a dimensionless three-loop integral function that depends on the masses present inside the loop. NSI in this model are induced by the charged scalar H − . After integrating out the heavy scalars, the NSI expression can be written as This is similar to the heavy charged scalar contribution in Eq. (4.81). However, since the same Yukawa couplings y eαa contribute to the electron mass in Eq. (7.69), we expect ε αβ ∝ y 2 e tan 2 β ∼ O 10 −10 , (7.73) where y e is the electron Yukawa coupling in the SM. Thus, the maximum NSI in this model are of order of O 10 −10 , as summarized in Table XX.

Cocktail Model
This model [39] induces operator O 9 of Eq. (1.4) at the three-loop level. The model includes two SU (2) L -singlet scalars η + (1, 1, 1) and k ++ (1, 1, 2), and a second scalar doublet Φ 2 1, 2, 1 2 , in addition to the SM Higgs doublet Φ 1 1, 2, 1 2 . The fields η + and Φ 2 are odd under a Z 2 symmetry, while k ++ and all SM fields are even. With this particle content, the relevant term in the Lagrangian reads as −L Y ⊃ y αβ Φ 1 L α c β + Y αβ c α β k ++ + H.c. , (7.74) which breaks lepton number when combined with the following cubic and quartic terms in the potential: The Φ 2 field is inert and does not get a VEV. After electroweak symmetry breaking, it can be written as For κ 1 = 0, the singly-charged state φ + 2 mixes with η + (with mixing angle β), giving rise to two singly-charged scalar mass eigenstates: where s β ≡ sin β and c β ≡ cos β.
The neutrino mass matrix is obtained from the three-loop diagram as shown in Fig. 54 and reads as [39] M ν ∼ g 2 (16π 2 ) 3 M (Y + Y T )M , (7.78) where M stands for the diagonal charged lepton mass matrix. As for the NSI, since both Φ 2 and η + are odd under Z 2 and the SM fields are even, there is no tree-level NSI in this model. Note that neutrino mass generation utilizes the W boson couplings, thus the neutrino matter effects in this model are the same as in the SM.

Leptoquark variant of the KNT model
One can replace the charged leptons in the KNT model (cf. Sec. 7.3.1) by quarks, and the charged scalars by leptoquarks. The effective operator induced in this model remains as O 9 or Eq. (1.4). To achieve this, two isospin-singlet scalar LQs χ −1/3 a 3, 1, − 1 3 (with a = 1, 2) and at least two SM-singlet right-handed neutrinos N α (1, 1, 0) (with α = 1, 2) are supplemented to the SM fields. A Z 2 symmetry is invoked under which χ −1/3 2 and N are odd, while the rest of the fields are even. The relevant Yukawa Lagrangian is as follows: Here the first term expands to give λ αβ (ν α d β − α u β ) χ generate neutrino masses at three-loop level, as shown in Fig. 55. The neutrino mass matrix reads as where the factor 15 comes from total color-degrees of freedom, M d and M N are the diagonal down-type quark and right-handed neutrino mass matrices, respectively, and I is a dimensionless three-loop integral that depends on the ratio of the masses of particles inside the loop [37]. NSI in this model arise from the χ −1/3 1 interactions with neutrinos and down-quarks. The expression for NSI parameters is given by which is the same as the singlet contribution in Eq. (5.48). The maximum NSI for this model are the same as those given in Eq. (7.38) and are summarized in Table XX. 8 Type II radiative models As discussed in the introduction (cf. Sec. 1.1), type-II radiative neutrino mass models in our nomenclature contain no SM particle inside the loop diagrams generating m ν , and therefore, do not generally contribute to tree-level NSI, although small loop-level NSI effects are possible [200]. To illustrate this point, let us take the scotogenic model [44] as a prototypical example. The new particles introduced in this model are SM-singlet fermions N α (1, 1, 0) (with α = 1, 2, 3) and an SU (2) L doublet scalar η 1, 2, 1 2 : (η + , η 0 ). A Z 2 symmetry is imposed under which the new fields N α and η are odd, while all the SM fields are even. The new Yukawa interactions in this model are given by −L Y ⊃ h αβ (ν α η 0 − α η + )N β + 1 2 (M N ) αβ N α N β + H.c. where Φ is the SM Higgs doublet, the Lagrangian (8.1) gives rise to neutrino mass at oneloop, as shown in Fig. 56. Since this diagram does not contain any SM fields inside the loop, it cannot be cut to generate an effective higher-dimensional operator of the SM. Therefore, we label it as a type-II radiative model. The neutrino mass in this model is given by where we have assumed M N to be diagonal, and m 2 0 is the average squared mass of the real and imaginary parts of η 0 . It is clear from Eq. (8.3) that the neutrino mass is not chirally suppressed by any SM particle mass.
A new example of type-II-like radiative model is shown in Fig. 57, where the new particles added are as follows: one color-sextet diquark ∆ 6, 1, 4 3 , one SU (2) L doublet scalar LQ δ 3, 2, 7 6 = (δ 5/3 , δ 2/3 ), and an SU (2) L singlet scalar LQ ξ 3, 1, 2 3 . The relevant Yukawa Lagrangian is given by −L Y ⊃ f αβ (ν α δ 2/3 − α δ 5/3 )u c β + λ αβ u c α ∆u c β + H.c. (8.4) Together with the scalar potential terms V ⊃ µδ † Φδ + µ δ 2 ∆ + H.c. , (8.5) where Φ is the SM Higgs doublet, the Lagrangian (8.4) gives rise to neutrino mass at two-loop level, as shown in Fig. 57. The neutrino mass can be approximated as follows: where m 1 and m 2 are the masses of the heaviest two LQs among the δ, ξ and ∆ fields that run in the loop. Thus, although this model can be described as arising from an effective ∆L = 2 operator O 1 of Eq. (1.5), the neutrino mass has no chiral suppression here. In this sense, this can be put in the type-II radiative model category, although it leads to tree-level NSI induced by the δ LQs, as in the one-loop type-I model discussed in Sec. 7.1.6. A similar two-loop radiative model without the chiral suppression can be found in Ref. [201].

Conclusion
We have made a comprehensive analysis of neutrino non-standard interactions generated by new scalars in radiative neutrino mass models. For this purpose, we have proposed a new nomenclature to classify radiative neutrino mass models, viz., the class of models with at least one SM particle in the loop are dubbed as type-I radiative models, whereas those models with no SM particles in the loop are called type-II radiative models. From NSI perspective, the type-I radiative models are most interesting, as the neutrino couples to a SM fermion (matter field) and a new scalar directly, thus generating NSI at tree-level, unlike type-II radiative models. After taking into account various theoretical and experimental constraints, we have derived the maximum possible NSI in all the type-I radiative models.
Our results are summarized in Fig. 58 and Table XX.
We have specifically analyzed two popular type-I radiative models, namely, the Zee model and its variant with LQs replacing the charged scalars, in great detail. In the Zee model with SU (2) L singlet and doublet scalar fields, we find that large NSI can be obtained via the exchange of a light charged scalar, arising primarily from the SU (2) L -singlet field but with some admixture of the doublet field. A light charged scalar with mass as low as ∼100 GeV is found to be consistent with various experimental constraints, including charged lepton flavor violation (cf. Sec. 4.5), monophoton constraints from LEP (cf. Sec. 4.11), direct searches for charged scalar pair and single production at LEP (cf. Sec. 4.7.1) and LHC (cf. Sec. 4.7.2), Higgs physics constraints from LHC (cf. Sec. 4.10), and lepton universality in W ± (cf. Sec. 4.8) and τ (cf. Sec. 4.9) decays. In addition, for the Yukawa couplings and the mixing between singlet and doublet scalars, we have considered the contact interaction limits from LEP (cf. Sec. 4.6), electroweak precision constraints from T -parameter (cf. Sec. 4.4), charge breaking minima of the Higgs potential (cf. Sec. 4.3), as well as perturbative unitarity of Yukawa and quartic couplings. After imposing all these constraints, we find diagonal values of the NSI parameters (ε ee , ε µµ , ε τ τ ) can be as large as (8%, 3.8%, 43%), while the off-diagonal NSI parameters (ε eµ , ε eτ , ε µτ ) can be at most (10 −3 %, 0.56%, 0.34%), as summarized in Fig. 58 and Table IX. Most of these NSI values are still allowed by the global fit constraints from neutrino oscillation and scattering experiments, and some of these parameters can be probed at future long-baseline neutrino oscillation experiments, such as DUNE.
We have also analyzed in detail the LQ version of the Zee model, the results of which can be applied to other LQ models with minimal modification. This analysis took into account the experimental constraints from direct searches for LQ pair and single production at LHC (cf. Sec. 5.3), as well as the low-energy constraints from APV (cf. Sec. 5.1.1), charged lepton flavor violation (cf. Secs. 5.1.4 and 5.1.5) and rare meson decays (cf. Sec. 5.1.6), apart from the theoretical constraints from perturbative unitarity of the Yukawa couplings. Including all these constraints we found that diagonal NSI (ε ee , ε µµ , ε τ τ ) can be as large as (0.4%, 21.6%, 34.3%), while off-diagonal NSI (ε eµ , ε eτ ε µτ ) can be as large as (10 −5 %, 0.36%, 0.43%), as summarized in Fig. 58 and Table XVII. A variant of the LQ model with triplet LQs (cf. Sec. 6) allows for larger ε τ τ ) which can be as large as 51.7%. Neutrino scattering experiments are found to be the most constraining for the diagonal Figure 58: Summary of maximum NSI strength |ε αβ | allowed in different classes of radiative neutrino mass models discussed here. Red, yellow, green, cyan, blue and purple bars correspond to the Zee model, minimal radiative inverse seesaw model, leptoquark model with singlet, doublet and triplet leptoquarks, and Zee-Babu model respectively.
NSI parameters ε ee and ε µµ , while the cLFV searches are the most constraining for the off-diagonal NSI. ε τ τ is the least constrained and can be probed at future long-baseline neutrino oscillation experiments, such as DUNE, whereas the other NSI parameters are constrained to be below the DUNE sensitivity reach.      [118] ATLAS Collaboration, "Combined measurements of Higgs boson production and decay using up to 80 fb −1 of proton-proton collision data at √ s = 13 TeV collected with the ATLAS experiment," Tech. Rep. ATLAS-CONF-2019-005, 2019.   [129] CMS Collaboration, "Measurement of Higgs boson production and decay to the τ τ final state," Tech. Rep. CMS-PAS-HIG-18-032, 2019.