Non-Standard Interactions in Radiative Neutrino Mass Models

TL;DR

This paper investigates how radiative Majorana neutrino mass models inherently produce observable nonstandard neutrino interactions with matter, analyzing their magnitude, constraints, and implications across various loop-induced models.

Contribution

It classifies radiative neutrino mass models into two types and demonstrates that type-I models can generate observable tree-level NSI consistent with experimental constraints.

Findings

Type-I models can produce observable NSI at tree-level.

Diagonal NSI can reach up to 9.3% in certain models.

Off-diagonal NSI are generally below 1%, constrained by low-energy experiments.

Abstract

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.

Tables

Table 1. Table I: Summary of new particles, their S ​ U ​ ( 3 ) c × S ​ U ​ ( 2 ) L × U ​ ( 1 ) Y 𝑆 𝑈 subscript 3 𝑐 𝑆 𝑈 subscript 2 𝐿 𝑈 subscript 1 𝑌 SU(3)_{c}\times SU(2)_{L}\times U(1)_{Y} quantum numbers (with the non-Abelian charges in boldface), field components and electric charges (in superscript), and corresponding Lagrangian terms responsible for NSI in various type-I radiative neutrino mass models discussed in Secs.  4 , 5 and 7 . Here Φ ~ = i ​ τ 2 ​ Φ ⋆ ~ Φ 𝑖 subscript 𝜏 2 superscript Φ ⋆ \widetilde{\Phi}=i\tau_{2}\Phi^{\star} , with τ 2 subscript 𝜏 2 \tau_{2} being the second Pauli matrix. For a singly charged scalar, η + superscript 𝜂 \eta^{+} and h + superscript ℎ h^{+} are used interchangeably, to be consistent with literature.

Particle Content	Lagrangian term
${\eta }^{+}(\mathrm{𝟏},\mathrm{𝟏},1)$ or $h^{+}(\mathrm{𝟏},\mathrm{𝟏},1)$	$f_{\alpha \beta }{L}_{\alpha }{L}_{\beta }{\eta }^{+}$ or $f_{\alpha \beta }{L}_{\alpha }{L}_{\beta }{h}^{+}$
$\Phi (\mathrm{𝟏},\mathrm{𝟐},\frac{1}{2})=({\varphi }^{+},{\varphi }^0)$	$Y_{\alpha \beta }{L}_{\alpha }{\ell }_{\beta }^c\tilde{\Phi }$
$\Omega (\mathrm{𝟑},\mathrm{𝟐},\frac{1}{6})=\left(\begin{array}{c}{\omega }^{2/3},{\omega }^{-1/3}\end{array}\right)$	${\lambda }_{\alpha \beta }{L}_{\alpha }{d}_{\beta }^c\Omega$
$\chi (\mathrm{𝟑},\mathrm{𝟏},-\frac{1}{3})$	${\lambda }_{\alpha \beta }^{\prime }{L}_{\alpha }{Q}_{\beta }{\chi }^{\star }$
$\overline{\rho }(\overline{\mathrm{𝟑}},\mathrm{𝟑},\frac{1}{3})=({\overline{\rho }}^{4/3},{\overline{\rho }}^{1/3},{\overline{\rho }}^{-2/3})$	${\lambda }_{\alpha \beta }^{\prime \prime }{L}_{\alpha }{Q}_{\beta }\overline{\rho }$
$\delta (\mathrm{𝟑},\mathrm{𝟐},\frac{7}{6})=({\delta }^{5/3},{\delta }^{2/3})$	${\lambda }_{\alpha \beta }^{\prime \prime \prime }{L}_{\alpha }{u}_{\beta }^c\delta$
$\Delta (\mathrm{𝟏},\mathrm{𝟑},1)=({\Delta }^{++},{\Delta }^{+},{\Delta }^0)$	$f_{\alpha \beta }^{\prime }{L}_{\alpha }{L}_{\beta }\Delta$

Table 2. Table II: Minimal UV completions of operator 𝒪 2 subscript 𝒪 2 \mathcal{O}_{2}   Cai:2014kra . Here ϕ italic-ϕ \phi and η 𝜂 \eta generically denote scalars and ψ 𝜓 \psi is a generic vectorlike fermion. The SM quantum numbers of these new fields are as indicated.

	${𝒪}_2^1$
	$L\left(LL\right)\left(e^{c}H\right)$
$\varphi$	$(\mathrm{𝟏},\mathrm{𝟏},1)$
$\psi$	$(\mathrm{𝟏},\mathrm{𝟐},-\frac{3}{2})$

Table 3. Table III: Minimal UV completions of operators 𝒪 3 ​ a subscript 𝒪 3 𝑎 \mathcal{O}_{3a} and 𝒪 3 ​ b subscript 𝒪 3 𝑏 \mathcal{O}_{3b} Cai:2014kra . Here the models in the top segment require a new scalar ϕ italic-ϕ \phi and a vectorlike fermion ψ 𝜓 \psi , while those in the lower segment require two scalar fields ϕ italic-ϕ \phi and η 𝜂 \eta .

	${𝒪}_3^1$	${𝒪}_3^2$	${𝒪}_3^3$	${𝒪}_3^4$	${𝒪}_3^5$	${𝒪}_3^6$
	$Q\left(LL\right)\left(d^{c}H\right)$	$d^c\left(LL\right)\left(QH\right)$	$L\left(L{d}^c\right)\left(QH\right)$	$L\left(LQ\right)\left(d^{c}H\right)$	$L\left(LQ\right)\left(d^{c}H\right)$	$L\left(L{d}^c\right)\left(QH\right)$
$\varphi$	$(\mathrm{𝟏},\mathrm{𝟏},1)$	$(\mathrm{𝟏},\mathrm{𝟏},1)$	$(\mathrm{𝟑},\mathrm{𝟐},\frac{1}{6})$	$(\mathrm{𝟑},\mathrm{𝟏},-\frac{1}{3})$	$(\mathrm{𝟑},\mathrm{𝟑},-\frac{1}{3})$	$(\mathrm{𝟑},\mathrm{𝟐},\frac{1}{6})$
$\psi$	$(\mathrm{𝟑},\mathrm{𝟐},-\frac{5}{6})$	$(\mathrm{𝟑},\mathrm{𝟏},\frac{2}{3})$	$(\mathrm{𝟑},\mathrm{𝟏},\frac{2}{3})$	$(\mathrm{𝟑},\mathrm{𝟐},-\frac{5}{6})$	$(\mathrm{𝟑},\mathrm{𝟐},-\frac{5}{6})$	$(\mathrm{𝟑},\mathrm{𝟑},\frac{2}{3})$
	${𝒪}_{3a}$	${𝒪}_{3a}$	${𝒪}_{3a}$	${𝒪}_{3b}$	${𝒪}_{3a},{𝒪}_{3b}$	${𝒪}_{3a},{𝒪}_{3b}$

Table 4. Table IV: Minimal UV completions of the operators 𝒪 4 ​ a subscript 𝒪 4 𝑎 \mathcal{O}_{4a} and 𝒪 4 ​ b subscript 𝒪 4 𝑏 \mathcal{O}_{4b} . Note that only the operator 𝒪 4 ​ b subscript 𝒪 4 𝑏 \mathcal{O}_{4b} is generated. Fields ϕ italic-ϕ \phi and η 𝜂 \eta are scalars, while the ψ 𝜓 \psi fields are vectorlike fermions.

	${𝒪}_4^1$	${𝒪}_4^2$
	$\overline{Q}\left(LL\right)\left(\overline{u^c}H\right)$	$\overline{u^c}\left(LL\right)\left(\overline{Q}H\right)$
$\varphi$	$(\mathrm{𝟏},\mathrm{𝟏},1)$	$(\mathrm{𝟏},\mathrm{𝟏},1)$
$\psi$	$(\mathrm{𝟑},\mathrm{𝟐},\frac{7}{6})$	$(\mathrm{𝟑},\mathrm{𝟏},-\frac{1}{3})$
	${𝒪}_{4b}$	${𝒪}_{4b}$

Table 5. Table V: Minimal UV completions of operator 𝒪 8 subscript 𝒪 8 \mathcal{O}_{8} . Fields ϕ italic-ϕ \phi and η 𝜂 \eta are scalars, while the ψ 𝜓 \psi fields are vectorlike fermions.

	${𝒪}_8^1$	${𝒪}_8^2$	${𝒪}_8^3$
	$L\left(\overline{e^c}\overline{u^c}\right)\left(d^{c}H\right)$	$\overline{u^c}\left(L{d}^c\right)\left(\overline{e^c}H\right)$	$\overline{e^c}\left(L{d}^c\right)\left(\overline{u^c}H\right)$
$\varphi$	$(\mathrm{𝟑},\mathrm{𝟏},-\frac{1}{3})$	$(\mathrm{𝟑},\mathrm{𝟐},\frac{1}{6})$	$(\mathrm{𝟑},\mathrm{𝟐},\frac{1}{6})$
$\psi$	$(\mathrm{𝟑},\mathrm{𝟐},-\frac{5}{6})$	$(\mathrm{𝟏},\mathrm{𝟐},-\frac{1}{2})$	$(\mathrm{𝟑},\mathrm{𝟐},\frac{7}{6})$

Table 6. Table VI: Constraints on Yukawa couplings as a function of heavy neutral scalar mass from ℓ α → ℓ β + γ → subscript ℓ 𝛼 subscript ℓ 𝛽 𝛾 \ell_{\alpha}\to\ell_{\beta}+\gamma processes.

Process	Exp. bound	Constraint
$\mu \to e\gamma$	BR <  4.2 $\times {10}^{-13}$ TheMEG:2016wtm	$\left|Y_{\mu e}^{\star }{Y}_{ee}\right|$ <  $1.05\times {10}^{-3}{\left(\frac{m_H}{700\text{GeV}}\right)}^2$
$\tau \to e\gamma$	BR <  3.3 $\times {10}^{-8}$ Aubert:2009ag
$\tau \to \mu \gamma$	BR <  4.4 $\times {10}^{-8}$ Aubert:2009ag

Table 7. Table VII: Constraints on Yukawa couplings as a function of heavy neutral scalar mass from ℓ α → ℓ ¯ β ​ ℓ γ ​ ℓ δ → subscript ℓ 𝛼 subscript ¯ ℓ 𝛽 subscript ℓ 𝛾 subscript ℓ 𝛿 \ell_{\alpha}\to\bar{\ell}_{\beta}\ell_{\gamma}\ell_{\delta} decay (with at least two of the final state leptons of electron flavor to be relevant for NSI).

Process	Exp. bound	Constraint
${\mu }^{-}\to e^{+}{e}^{-}{e}^{-}$	BR <  1.0 $\times {10}^{-12}$  Bertl:1985mw	$\left|Y_{\mu e}^{\star }{Y}_{ee}\right|$ <  $3.28\times {10}^{-5}{\left(\frac{m_H}{700\text{GeV}}\right)}^2$
${\tau }^{-}\to e^{+}{e}^{-}{e}^{-}$	BR <  1.4 $\times {10}^{-8}$  Amhis:2016xyh	$\left|Y_{\tau e}^{\star }{Y}_{ee}\right|$ <  $9.05\times {10}^{-3}{\left(\frac{m_H}{700\text{GeV}}\right)}^2$
${\tau }^{-}\to e^{+}{e}^{-}{\mu }^{-}$	BR <  1.1 $\times {10}^{-8}$   Amhis:2016xyh	$\left|Y_{\tau e}^{\star }{Y}_{\mu e}\right|$ <  $5.68\times {10}^{-3}{\left(\frac{m_H}{700\text{GeV}}\right)}^2$

Table 8. Table VIII: Constraints on the ratio of heavy neutral scalar mass and the Yukawa couplings from LEP contact interaction bounds.

Process	LEP bound LEP:2003aa	Constraint
$e^{+}{e}^{-}\to e^{+}{e}^{-}$	${\Lambda }_{LR/RL}^{-}>10$ TeV	$\frac{m_H}{\left|Y_{ee}\right|}>1.99$ TeV
$e^{+}{e}^{-}\to {\mu }^{+}{\mu }^{-}$	${\Lambda }_{LR/RL}^{-}>7.9$ TeV	$\frac{m_H}{\left|Y_{\mu e}\right|}>1.58$ TeV
$e^{+}{e}^{-}\to {\tau }^{+}{\tau }^{-}$	${\Lambda }_{LR/RL}^{-}>2.2$ TeV	$\frac{m_H}{\left|Y_{\tau e}\right|}>0.44$ TeV

Table 9. Table IX: Maximum allowed NSI (with electrons) in the Zee model, after imposing constraints from CBM (Sec.  4.3 ), T 𝑇 T -parameter (Sec.  4.4 ), cLFV searches (Sec.  4.5 ), LEP contact interaction (Sec.  4.6 ), direct collider searches (Sec.  4.7 ), lepton universality (LU) in W 𝑊 W decays (Sec.  4.8 ), LHC Higgs data (Sec.  4.10 ), and LEP monophoton searches (Sec.  4.11 ). We also impose the individual constraints, taking one NSI parameter at a time, from either neutrino-electron scattering or neutrino oscillation experiments (as shown in the third column), like CHARM-II  Barranco:2007ej , TEXONO  Deniz:2010mp and BOREXINO  Agarwalla:2019smc (only ϵ α ​ β e ​ R superscript subscript italic-ϵ 𝛼 𝛽 𝑒 𝑅 \epsilon_{\alpha\beta}^{eR} are considered, cf. Eq. ( 4.88 )) or IceCube  Esmaili:2013fva as well as the global-fit constraints (as shown in the fourth column), taking all NSI parameters simultaneously, from neutrino oscillation+COHERENT data  Esteban:2018ppq (only ε α ​ β p superscript subscript 𝜀 𝛼 𝛽 𝑝 \varepsilon_{\alpha\beta}^{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)  dev_pondd .

NSI	Zee Model	Individual	Global-fit	DUNE
	Prediction (Max.)	constraints	constraints Esteban:2018ppq	sensitivity dev_pondd
${\epsilon }_{ee}$	0.08	$[-0.07,0.08]$ Deniz:2010mp	$[-0.010,2.039]$	$[-0.185,0.380]$
	(TEXONO)			($[-0.130,0.185]$)
${\epsilon }_{\mu \mu }$	0.038	$[-0.03,0.03]$ Barranco:2007ej	$[-0.364,1.387]$	$[-0.290,0.390]$
	(CHARM)	$[-0.017,0.038]$ (ours)		($[-0.192,0.240]$)
${\epsilon }_{\tau \tau }$	0.093	$[-0.093,0.093]$ Esmaili:2013fva	$[-0.350,1.400]$	$[-0.360,0.145]$
	(IceCube)			($[-0.120,0.095]$)
${\epsilon }_{e\mu }$	$1.5\times {10}^{-5}$	$[-0.13,0.13]$ Barranco:2007ej	$[-0.179,0.146]$	$[-0.025,0.052]$
	(LEP + LU + cLFV + $T$-param.)			( $[-0.017,0.040]$)
${\epsilon }_{e\tau }$	0.0056	$[-0.19,0.19]$ Deniz:2010mp	$[-0.860,0.350]$	$[-0.055,0.023]$
	(LEP + LU + cLFV + $T$-param.)			($[-0.042,0.012]$)
${\epsilon }_{\mu \tau }$	0.0034	$[-0.10,0.10]$ Barranco:2007ej	$[-0.035,0.028]$	$[-\mathrm{0.0.015},0.013]$
	(LEP + LU + cLFV + $T$-param)			($[-0.010,0.010]$)

Table 10. Table X: Values of parameters chosen for different sets of Yukawa structure given in Eqs. ( 4.105 )-( 4.113 ) to fit the neutrino oscillation data.

BP	$x_1$	$x_2$	$y_{11}$	$y_{12}$	$y_{13}$	$y_{21}$	$y_{22}$	$y_{23}$	$y_{31}$	$y_{32}$	$y_{33}$	$a_0\left({10}^{-9}\right)$
BP I (IH)	$-7950$	34	$-1.0$	0	$-0.01$	0	0.001	0.08	0	0.05	0.70	0.017
BP II (IH)	14	4.7	0	0.05	0.01	1.0	0	0.02	0	0.06	0.03	0.19
BP III (NH)	$-9.9$	0.27	0.01	0	0.07	0	0.13	$-0.007$	$-1.0$	0	$-0.036$	0.6

Table 11. Table XI: 3 ​ σ 3 𝜎 3\sigma allowed ranges of the neutrino oscillation parameters from a recent global-fit  Esteban:2018azc (without NSI), along with the model predictions for each BP.

Oscillation	$3\sigma$ allowed range	Model prediction
parameters	from NuFit4 Esteban:2018azc	BP I (IH)	BP II (IH)	BP III (NH)
$\Delta m_{21}^2({10}^{-5}{\mathrm{eV}}^2$)	6.79 - 8.01	7.388	7.392	7.390
$\Delta m_{23}^2\left({10}^{-3}{\mathrm{eV}}^2\right)$(IH)	2.412 - 2.611	2.541	2.488	-
$\Delta m_{31}^2\left({10}^{-3}{\mathrm{eV}}^2\right)$(NH)	2.427 - 2.625	-	-	2.505
${\sin }^2{\theta }_{12}$	0.275 - 0.350	0.295	0.334	0.316
${\sin }^2{\theta }_{23}$ (IH)	0.423 - 0.629	0.614	0.467	-
${\sin }^2{\theta }_{23}$ (NH)	0.418 - 0.627	-	-	0.577
${\sin }^2{\theta }_{13}$ (IH)	0.02068 - 0.02463	0.0219	0.0232	-
${\sin }^2{\theta }_{13}$(NH)	0.02045 - 0.02439	-	-	0.0229

Table 12. Table XII: Constraints on Yukawa couplings and LQ masses from μ − e 𝜇 𝑒 \mu-e conversion in different nuclei. For | λ e ​ d ′ ⁣ ⋆ ​ λ μ ​ d ′ | superscript subscript 𝜆 𝑒 𝑑 ′ ⋆ subscript superscript 𝜆 ′ 𝜇 𝑑 |\lambda_{ed}^{\prime\star}\lambda^{\prime}_{\mu d}| , the same constraints apply, with m ω subscript 𝑚 𝜔 m_{\omega} replaced by m χ subscript 𝑚 𝜒 m_{\chi} .

Nucleus	Experimental	$Z_{\mathrm{eff}}$	$F_p$	${\Gamma }_N$ Suzuki:1987jf	Constraint
	Limit	Kitano:2002mt	Kitano:2002mt	$\left({10}^6{s}^{-1}\right)$	on $\left|{\lambda }_{ed}^{\star }{\lambda }_{\mu d}\right|$
${}_{}{}^{48}{}_{22}{}^{}$Ti	BR  Kaulard:1998rb	17.6	0.54	2.59
${}_{}{}^{197}{}_{79}{}^{}$Au	BR  Bertl:2006up	33.5	0.16	13.07
${}_{}{}^{208}{}_{82}{}^{}$Pb	BR  Honecker:1996zf	34.0	0.15	13.45

Table 13. Table XIII: Constraints on the Yukawa couplings λ ′ superscript 𝜆 ′ \lambda^{\prime} as a function of the singlet LQ mass from ℓ α → ℓ β ​ γ → subscript ℓ 𝛼 subscript ℓ 𝛽 𝛾 \ell_{\alpha}\to\ell_{\beta}\gamma processes.

Process	Exp. limit	Constraint
$\mu \to e\gamma$	BR < 4.2 $\times {10}^{-13}$ TheMEG:2016wtm	$\left|{\lambda }_{ed}^{\prime }{\lambda }_{\mu d}^{\prime \star }\right|$ <  $2.4\times {10}^{-3}{\left(\frac{m_{\chi }}{\text{TeV}}\right)}^2$
$\tau \to e\gamma$	BR < 3.3 $\times {10}^{-8}$ Aubert:2009ag	$\left|{\lambda }_{ed}^{\prime }{\lambda }_{\tau d}^{\prime \star }\right|$ <  $1.6{\left(\frac{m_{\chi }}{\text{TeV}}\right)}^2$
$\tau \to \mu \gamma$	BR < 4.4 $\times {10}^{-8}$ Aubert:2009ag	$\left|{\lambda }_{\mu d}^{\prime \star }{\lambda }_{\tau d}^{\prime }\right|$ <  $1.9{\left(\frac{m_{\chi }}{\text{TeV}}\right)}^2$

Table 14. Table XIV: Constraints on couplings and the LQ mass from semileptonic tau decays. Exactly the same constraints apply to λ ′ superscript 𝜆 ′ \lambda^{\prime} couplings, with m ω subscript 𝑚 𝜔 m_{\omega} replaced by m χ subscript 𝑚 𝜒 m_{\chi} .

Process	Exp. limit Tanabashi:2018oca	Constraint
$\tau \to \mu {\pi }^0$	BR < $1.1\times {10}^{-7}$
$\tau \to e{\pi }^0$	BR < $8\times {10}^{-8}$
$\tau \to \mu \eta$	BR < $6.5\times {10}^{-8}$
$\tau \to e\eta$	BR < $9.2\times {10}^{-8}$
$\tau \to \mu {\eta }^{\prime }$	BR < $1.3\times {10}^{-7}$
$\tau \to e{\eta }^{\prime }$	BR < $1.6\times {10}^{-7}$

Table 15. Table XV: Constraints on the χ 𝜒 \chi LQ Yukawa couplings from D 0 → ℓ + ​ ℓ − → superscript 𝐷 0 superscript ℓ superscript ℓ D^{0}\to\ell^{+}\ell^{-} and D + → π + ​ ℓ + ​ ℓ − → superscript 𝐷 superscript 𝜋 superscript ℓ superscript ℓ D^{+}\to\pi^{+}\ell^{+}\ell^{-} decays.

Process	Exp. limit Tanabashi:2018oca	Constraint
$D^0\to e^{+}{e}^{-}$	BR < $7.9\times {10}^{-8}$
$D^0\to {\mu }^{+}{\mu }^{-}$	BR < $6.2\times {10}^{-9}$
$D^{+}\to {\pi }^{+}{e}^{+}{e}^{-}$	BR < $1.1\times {10}^{-6}$
$D^{+}\to {\pi }^{+}{\mu }^{+}{\mu }^{-}$	BR < $7.3\times {10}^{-8}$
$D^{+}\to {\pi }^{+}{e}^{+}{\mu }^{-}$	BR < $3.6\times {10}^{-6}$

Table 16. Table XVI: Constraints on the ratio of LQ mass and the Yukawa coupling from LEP  LEP:2003aa and HERA  Abramowicz:2019uti contact interaction bounds.

LQ	LEP		HERA
type	Exp. bound LEP:2003aa	Constraint	Exp. bound Abramowicz:2019uti	Constraint
${\omega }^{2/3}$	${\Lambda }_{LR}^{-}>5.1$ TeV	$\frac{m_{\omega }}{\left|{\lambda }_{ed}\right|}>1.017$ TeV	${\Lambda }_{LR}^{-}>4.7$ TeV	$\frac{m_{\omega }}{\left|{\lambda }_{ed}\right|}>0.937$ TeV
${\chi }^{-1/3}$	${\Lambda }_{LL}^{-}>3.7$ TeV	$\frac{m_{\chi }}{\left|{\lambda }_{ed}\right|}>0.738$ TeV	${\Lambda }_{LL}^{-}>12.8$ TeV	$\frac{m_{\chi }}{\left|{\lambda }_{ed}\right|}>2.553$ TeV

Table 17. Table XVII: Maximum allowed NSI (with d 𝑑 d -quarks) in the one-loop LQ model, after imposing the constraints from APV (Sec.  5.1.1 ), cLFV (Secs.  5.1.2 , 5.1.5 , 5.1.6 ), LEP and HERA contact interaction (Sec.  5.2 ), perturbative unitarity and collider (Sec.  5.3 ) constraints. We also impose the constraints from neutrino-nucleon scattering experiments, like CHARM II  Davidson:2003ha , NuTeV  Davidson:2003ha , COHERENT  Coloma:2017ncl and IceCube  Salvado:2016uqu , as well as the global-fit constraints from neutrino oscillation+COHERENT data  Esteban:2018ppq , whichever is stronger. The scattering and global-fit constraints are on ε α ​ β d superscript subscript 𝜀 𝛼 𝛽 𝑑 \varepsilon_{\alpha\beta}^{d} , so it has been scaled by a factor of 3 for the constraint on ε α ​ β subscript 𝜀 𝛼 𝛽 \varepsilon_{\alpha\beta} 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 ε e ​ e subscript 𝜀 𝑒 𝑒 \varepsilon_{ee} where the APV constraint is weaker than HERA (Fig.  33 (a))) and for ε μ ​ μ subscript 𝜀 𝜇 𝜇 \varepsilon_{\mu\mu} which has an additional constraint from D + → π + ​ μ + ​ μ − → superscript 𝐷 superscript 𝜋 superscript 𝜇 superscript 𝜇 D^{+}\to\pi^{+}\mu^{+}\mu^{-} decay (see Fig.  33 (b)). In the last column, we also show the future DUNE sensitivity  dev_pondd for 300 kt.MW.yr exposure (and 850 kt.MW.yr in parentheses).

	LQ model prediction (Max.)		Individual	Global-fit	DUNE
NSI	Doublet	Singlet	constraints	constraints Esteban:2018ppq	sensitivity dev_pondd
${\epsilon }_{ee}$	0.004	0.0069	$[-1.8,1.5]$ Davidson:2003ha	$[-0.036,1.695]$	$[-0.185,0.380]$
	(LHC + APV)	(LHC+HERA)			($[-0.130,0.185]$)
${\epsilon }_{\mu \mu }$	0.216	0.0086	$[-0.024,0.045]$ Davidson:2003ha	$[-0.309,1.083]$	$[-0.290,0.390]$
	(LHC+PU)	($D\to \pi \mu \mu$)	$[0.0277,0.0857]$ (ours)		($[-0.192,0.240]$)
${\epsilon }_{\tau \tau }$	0.343		$[-0.225,0.99]$ Coloma:2017ncl	$[-0.306,1.083]$	$[-0.360,0.145]$
	(LHC + Unitarity)				($[-0.120,0.095]$)
${\epsilon }_{e\mu }$	$1.5\times {10}^{-7}$		$[-0.21,0.12]$ Coloma:2017ncl	$[-0.174,0.147]$	$[-0.025,0.052]$
	(LHC + $\mu -e$ conv.)				($[-0.017,0.040]$)
${\epsilon }_{e\tau }$	0.0036		$[-0.39,0.36]$ Coloma:2017ncl	$[-0.618,0.330]$	$[-0.055,0.023]$
	(LHC + $\tau \to e{\pi }^0)$				($[-0.042,0.012]$)
${\epsilon }_{\mu \tau }$	0.0043		$[-0.018,0.0162]$ Salvado:2016uqu	$[-0.033,0.027]$	$[-0.015,0.013]$
	(LHC + $\tau \to \mu {\pi }^0$)				($[-0.010,0.010]$)

Table 18. Table XVIII: Constraints on couplings and the LQ mass from semileptonic tau decays in the triplet LQ case. Here we have assumed all the triplet fields ( ρ ¯ 4 / 3 superscript ¯ 𝜌 4 3 \bar{\rho}^{4/3} , r ​ h ​ o ¯ 1 / 3 superscript ¯ 𝑟 ℎ 𝑜 1 3 \bar{rho}^{1/3} , ρ ¯ − 2 / 3 superscript ¯ 𝜌 2 3 \bar{\rho}^{-2/3} ) to have the same mass m ρ subscript 𝑚 𝜌 m_{\rho} .

Process	Exp. limit Tanabashi:2018oca	Constraint
$\tau \to \mu {\pi }^0$	BR < $1.1\times {10}^{-7}$
$\tau \to e{\pi }^0$	BR < $8\times {10}^{-8}$
$\tau \to \mu \eta$	BR < $6.5\times {10}^{-8}$
$\tau \to e\eta$	BR < $9.2\times {10}^{-8}$
$\tau \to \mu {\eta }^{\prime }$	BR < $1.3\times {10}^{-7}$
$\tau \to e{\eta }^{\prime }$	BR < $1.6\times {10}^{-7}$

Table 19. Table XIX: Constraints on the singly-charged scalar Yukawa couplings  Nebot:2007bc . Here g α exp subscript superscript 𝑔 exp 𝛼 g^{\rm exp}_{\alpha} stands for the effective gauge coupling extracted from muon and tau decays in the different leptonic channels.

Observable	Exp. limit	Constraint
$\mu \to e\gamma$	BR < $4.2\times {10}^{-13}$ TheMEG:2016wtm	$\left|f_{e\tau }^{\star }{f}_{\mu \tau }\right|$ < $1.09\times {10}^{-3}{\left(\frac{m_{h^{+}}}{\mathrm{TeV}}\right)}^2$
$\tau \to e\gamma$	BR < $3.3\times {10}^{-8}$ Aubert:2009ag	$\left|f_{e\mu }^{\star }{f}_{\mu \tau }\right|$ < $0.71{\left(\frac{m_{h^{+}}}{\mathrm{TeV}}\right)}^2$
$\tau \to \mu \gamma$	BR < $4.4\times {10}^{-8}$ Aubert:2009ag	$\left|f_{e\mu }^{\star }{f}_{e\tau }\right|$ < $0.82{\left(\frac{m_{h^{+}}}{\mathrm{TeV}}\right)}^2$
lep./had. univ.	$\sum _{q=d,s,b}{\left|V_{uq}^{\exp }\right|}^2=0.9992\pm 0.0011$ Tanabashi:2018oca
$\mu /e$ univ.	$g_{\mu }^{\exp }/g_e^{\exp }=1.0001\pm 0.0020$ Tanabashi:2018oca
$\tau /\mu$ univ.	$g_{\tau }^{\exp }/g_{\mu }^{\exp }=1.0004\pm 0.0022$ Tanabashi:2018oca
$\tau /e$ univ.	$g_{\tau }^{\exp }/g_e^{\exp }=1.0004\pm 0.0023$  Tanabashi:2018oca

Table 20. Table XX: A comprehensive summary of type-I radiative neutrino mass models, with the new particle content and their ( S ​ U ​ ( 3 ) c , S ​ U ​ ( 2 ) L , U ​ ( 1 ) Y ) 𝑆 𝑈 subscript 3 𝑐 𝑆 𝑈 subscript 2 𝐿 𝑈 subscript 1 𝑌 (SU(3)_{c},~{}SU(2)_{L},~{}U(1)_{Y}) charges, and the maximum tree-level NSI allowed in each model. Red-colored exotic particles are odd under a Z 2 subscript 𝑍 2 Z_{2} symmetry. 𝒮 𝒮 {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\mathcal{S}} and ℱ ℱ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\mathcal{F}} represent scalar and fermion fields respectively.

Term	$𝒪$	Model	Loop	$𝒮/$	New particles	Max NSI @ tree-level
			level	$ℱ$		$\left|{\epsilon }_{ee}\right|$	$\left|{\epsilon }_{\mu \mu }\right|$	$\left|{\epsilon }_{\tau \tau }\right|$	$\left|{\epsilon }_{e\mu }\right|$	$\left|{\epsilon }_{e\tau }\right|$	$\left|{\epsilon }_{\mu \tau }\right|$
$L{\ell }^c{\Phi }^{\star }$	${𝒪}_2^2$	Zee Zee:1980ai	1	$𝒮$	${\eta }^{+}(\mathrm{𝟏},\mathrm{𝟏},1)$, ${\Phi }_2(\mathrm{𝟏},\mathrm{𝟐},1/2)$	0.08	0.038	0.093	$𝒪\left({10}^{-5}\right)$	0.0056	0.0034
	${𝒪}_9$	Zee-Babu Zee:1985id ; Babu:1988ki	2	$𝒮$	$h^{+}(\mathrm{𝟏},\mathrm{𝟏},1)$, $k^{++}(\mathrm{𝟏},\mathrm{𝟏},2)$
	${𝒪}_9$	KNT Krauss:2002px	3	$𝒮$	${\eta }_1^{+}(\mathrm{𝟏},\mathrm{𝟏},1)$, ${\eta }_2^{+}(\mathrm{𝟏},\mathrm{𝟏},1)$
				$ℱ$	$N(\mathrm{𝟏},\mathrm{𝟏},0)$	0	$0.0009$	$0.003$	0	0	0.003
$LL\eta$	${𝒪}_9$	1S-1S-1F Cepedello:2018rfh	3	$𝒮$	${\eta }_1(\mathrm{𝟏},\mathrm{𝟏},1)$, ${\eta }_2(\mathrm{𝟏},\mathrm{𝟏},3)$
				$ℱ$	$F(\mathrm{𝟏},\mathrm{𝟏},2)$
	${𝒪}_2^1$	1S-2VLL Cai:2014kra	1	$𝒮$	$\eta (\mathrm{𝟏},\mathrm{𝟏},1)$
				$ℱ$	$\Psi (\mathrm{𝟏},\mathrm{𝟐},-3/2)$
	${𝒪}_3^{\prime }$	AKS Aoki:2008av	3	$𝒮$	${\Phi }_2(\mathrm{𝟏},\mathrm{𝟐},1/2)$,  ${\eta }^{+}(\mathrm{𝟏},\mathrm{𝟏},1)$,  ${\eta }^0(\mathrm{𝟏},\mathrm{𝟏},0)$	$𝒪\left({10}^{-10}\right)$	$𝒪\left({10}^{-10}\right)$	$𝒪\left({10}^{-10}\right)$	$𝒪\left({10}^{-10}\right)$	$𝒪\left({10}^{-10}\right)$	$𝒪\left({10}^{-10}\right)$
$L{\ell }^c{\varphi }^{\star }$				$ℱ$	$N(\mathrm{𝟏},\mathrm{𝟏},0)$
—	${𝒪}_{d=15}$	Cocktail Gustafsson:2012vj	3	$𝒮$	${\eta }^{+}(\mathrm{𝟏},\mathrm{𝟏},1)$, $k^{++}(\mathrm{𝟏},\mathrm{𝟏},2),{\Phi }_2(\mathrm{𝟏},\mathrm{𝟐},1/2)$	0	0	0	0	0	0
$W/Z$	${𝒪}_2^{\prime }$	MRIS Dev:2012sg	1	$ℱ$	$N(\mathrm{𝟏},\mathrm{𝟏},0)$, $S(\mathrm{𝟏},\mathrm{𝟏},0)$	0.0013	$𝒪\left({10}^{-4}\right)$	0.0028	$𝒪\left({10}^{-5}\right)$	$𝒪\left({10}^{-4}\right)$	0.0012
$L\Omega d^c$	${𝒪}_3^8$	LQ variant of Zee AristizabalSierra:2007nf	1	$𝒮$	$\Omega (\mathrm{𝟑},\mathrm{𝟐},1/6)$, $\chi (\mathrm{𝟑},\mathrm{𝟏},-1/3)$	0.004	0.216	0.343	$𝒪\left({10}^{-7}\right)$	0.0036	0.0043
$\left(LQ{\chi }^{\star }\right)$	${𝒪}_8^4$	2LQ-1LQ Babu:2010vp	2	$𝒮$	$\Omega (\mathrm{𝟑},\mathrm{𝟐},1/6)$, $\chi (\mathrm{𝟑},\mathrm{𝟏},-1/3)$	(0.0069)	(0.0086)
	${𝒪}_3^3$	2LQ-1VLQ Babu:2011vb	2	$𝒮$	$\Omega (\mathrm{𝟑},\mathrm{𝟐},1/6)$
				$ℱ$	$U(\mathrm{𝟑},\mathrm{𝟏},2/3)$
$L\Omega d^c$	${𝒪}_3^6$	2LQ-3VLQ Cai:2014kra	1	$𝒮$	$\Omega (\mathrm{𝟑},\mathrm{𝟐},1/6)$
				$ℱ$	$\Sigma (\mathrm{𝟑},\mathrm{𝟑},2/3)$	0.004	0.093	0.093	$𝒪\left({10}^{-7}\right)$	0.0036	0.0043
	${𝒪}_8^2$	2LQ-2VLL Cai:2014kra	2	$𝒮$	$\Omega (\mathrm{𝟑},\mathrm{𝟐},1/6)$
				$ℱ$	$\psi (\mathrm{𝟏},\mathrm{𝟐},-1/2)$
	${𝒪}_8^3$	2LQ-2VLQ Cai:2014kra	2	$𝒮$	$\Omega (\mathrm{𝟑},\mathrm{𝟐},1/6)$
				$ℱ$	$\xi (\mathrm{𝟑},\mathrm{𝟐},7/6)$
$L\Omega d^c$	${𝒪}_3^9$	Triplet-Doublet LQ Cai:2014kra	1	$𝒮$	$\rho (\mathrm{𝟑},\mathrm{𝟑},-1/3)$, $\Omega (\mathrm{𝟑},\mathrm{𝟐},1/6)$	0.0059	0.0249	0.517	$𝒪\left({10}^{-8}\right)$	0.0050	0.0038
$\left(LQ\overline{\rho }\right)$
	${𝒪}_{11}$	LQ/DQ variant Zee-Babu Kohda:2012sr	2	$𝒮$	$\chi (\mathrm{𝟑},\mathrm{𝟏},-1/3)$ , $\Delta (\mathrm{𝟔},\mathrm{𝟏},-2/3)$
	${𝒪}_{11}$	Angelic Angel:2013hla	2	$𝒮$	$\chi (\mathrm{𝟑},\mathrm{𝟏},1/3)$
				$ℱ$	$F(\mathrm{𝟖},\mathrm{𝟏},0)$
$LQ{\chi }^{\star }$	${𝒪}_{11}$	LQ variant of KNT Nomura:2016ezz	3	$𝒮$	$\chi (\mathrm{𝟑},\mathrm{𝟏},-1/3)$,  ${\chi }_2(\mathrm{𝟑},\mathrm{𝟏},-1/3)$	0.0069	0.0086	0.093	$𝒪\left({10}^{-7}\right)$	0.0036	0.0043
				$ℱ$	$N(\mathrm{𝟏},\mathrm{𝟏},0)$
	${𝒪}_3^4$	1LQ-2VLQ Cai:2014kra	1	$𝒮$	$\chi (\mathrm{𝟑},\mathrm{𝟏},-1/3)$
				$ℱ$	$𝒬(\mathrm{𝟑},\mathrm{𝟐},-5/6)$
$L{u}^c\delta$	${\widetilde{𝒪}}_1$	3LQ-2LQ-1LQ (New)	1	$𝒮$	$\overline{\rho }(\overline{\mathrm{𝟑}},\mathrm{𝟑},1/3)$, $\delta (\mathrm{𝟑},\mathrm{𝟐},7/6)$, $\xi (\mathrm{𝟑},\mathrm{𝟏},2/3)$	0.004	0.216	0.343	$𝒪\left({10}^{-7}\right)$	0.0036	0.0043
$\left(LQ\overline{\rho }\right)$						(0.0059)	(0.007)	(0.517)		(0.005)	(0.0038)
$L{u}^c\delta$	${𝒪}_{d=13}$	3LQ-2LQ-2LQ(New)	2	$𝒮$	$\delta (\mathrm{𝟑},\mathrm{𝟐},7/6)$, $\Omega (\mathrm{𝟑},\mathrm{𝟐},1/6)$,$\hat{\Delta }({\mathrm{𝟔}}^{\star },\mathrm{𝟑},-1/3)$	0.004	0.216	0.343	$𝒪\left({10}^{-7}\right)$	0.0036	0.0043
$LQ\overline{\rho }$	${𝒪}_3^5$	3LQ-2VLQ Cai:2014kra	1	$𝒮$	$\overline{\rho }(\overline{\mathrm{𝟑}},\mathrm{𝟑},-1/3)$
				$ℱ$	$𝒬(\mathrm{𝟑},\mathrm{𝟐},-5/6)$	0.0059	0.0007	0.517	$𝒪\left({10}^{-7}\right)$	0.005	0.0038
	All Type-II Radiative models					0	0	0	0	0	0

Equations

$${\mathcal{O}}_1=L^i{L}^j{H}^k{H}^l{ϵ}_{ik}{ϵ}_{jl}$$

$$m_{\nu }\approx \frac{f^{2}h}{(16{\pi }^2{)}^2}\frac{m_{\mu }^2}{M},$$

$${\mathcal{O}}_2=L^i{L}^j{L}^k{e}^c{H}^l{ϵ}_{ij}{ϵ}_{kl}$$

$${\mathcal{O}}_9=L^i{L}^j{L}^k{e}^c{L}^l{e}^c{ϵ}_{ij}{ϵ}_{kl}.$$

$${\mathcal{O}}_1^{\prime }=L^i{L}^j{H}^k{H}^l{ϵ}_{ik}{ϵ}_{jl}(u^c{u}^c)(u^c{u}^c{)}^{\star }.$$

$$L\left(1,2,-\frac{1}{2}\right),e^c(1,1,1),Q\left(3,2,\frac{1}{6}\right),d^c\left(\overline{3},1,\frac{1}{3}\right),u^c\left(\overline{3},1,-\frac{2}{3}\right).$$

$${\mathcal{O}}_1$$

$${\mathcal{O}}_2$$

$${\mathcal{O}}_3$$

$${\mathcal{O}}_4$$

$${\mathcal{O}}_8$$

$${\mathcal{O}}_2^1=L(LL)(e^{c}H),{\mathcal{O}}_2^2=H(LL)(L{e}^c).$$

$$\displaystyle\mathcal{L}_{\rm NSI}^{\rm NC}\$$

$$\displaystyle\mathcal{L}_{\rm NSI}^{\rm CC}\$$

$$H=\frac{1}{2E}{U}_{\mathrm{PMNS}}\left(\begin{array}{ccc}0& 0& 0\\ 0& \Delta m_{21}^2& 0\\ 0& 0& \Delta m_{31}^2\end{array}\right)U_{\mathrm{PMNS}}^{†}+\sqrt{2}{G}_F{N}_e(x)\left(\begin{array}{ccc}1+{\epsilon }_{ee}& {\epsilon }_{e\mu }& {\epsilon }_{e\tau }\\ {\epsilon }_{e\mu }^{\star }& {\epsilon }_{\mu \mu }& {\epsilon }_{\mu \tau }\\ {\epsilon }_{e\tau }^{\star }& {\epsilon }_{\mu \tau }^{\star }& {\epsilon }_{\tau \tau }\end{array}\right),$$

$$\displaystyle\varepsilon_{\alpha\beta}\$$

$$={\epsilon }_{\alpha \beta }^{eV}+⟨\frac{N_p(x)}{N_e(x)}⟩(2{\epsilon }_{\alpha \beta }^{uV}+{\epsilon }_{\alpha \beta }^{dV})+⟨\frac{N_n(x)}{N_e(x)}⟩({\epsilon }_{\alpha \beta }^{uV}+2{\epsilon }_{\alpha \beta }^{dV}),$$

$${\epsilon }_{\alpha \beta }={\epsilon }_{\alpha \beta }^{eV}+\left[2+Y_n(x)\right]{\epsilon }_{\alpha \beta }^{uV}+\left[1+2Y_n(x)\right]{\epsilon }_{\alpha \beta }^{dV}.$$

$$-{\mathcal{L}}_Y\supset f_{\alpha \beta }{L}_{\alpha }^i{L}_{\beta }^j{ϵ}_{ij}{\eta }^{+}+(y_1{)}_{\alpha \beta }{\tilde{\Phi }}_1^i{L}_{\alpha }^j{\ell }_{\beta }^c{ϵ}_{ij}+(y_2{)}_{\alpha \beta }{\tilde{\Phi }}_2^i{L}_{\alpha }^j{\ell }_{\beta }^c{ϵ}_{ij}+H.c.,$$

$$-{\mathcal{L}}_Y\supset 2{\eta }^{+}\left[f_{e\mu }({\nu }_e\mu -{\nu }_{\mu }e)+f_{e\tau }({\nu }_e\tau -{\nu }_{\tau }e)+f_{\mu \tau }({\nu }_{\mu }\tau -{\nu }_{\tau }\mu )\right]+H.c.$$

$$V\supset \mu {\Phi }_1^i{\Phi }_2^j{ϵ}_{ij}{\eta }^{-}+H.c.,$$

$$\langle\Phi_{1}\rangle\ =\ \frac{1}{\sqrt{2}}\left(\begin{array}[]{c}0\\ v_{1}\\ \end{array}\right),\hskip 14.22636pt\langle\Phi_{2}\rangle\ =\ \frac{1}{\sqrt{2}}\left(\begin{array}[]{c}0\\ v_{2}e^{i\xi}\\ \end{array}\right).$$

$$\left(\begin{array}[]{c}H_{1}\\ H_{2}\\ \end{array}\right)\ =\ \left(\begin{array}[]{cc}c_{\beta}&e^{-i\xi}s_{\beta}\\ -e^{i\xi}s_{\beta}&c_{\beta}\\ \end{array}\right)\left(\begin{array}[]{c}\Phi_{1}\\ \Phi_{2}\\ \end{array}\right)\,,$$

$$H_{1}=\left(\begin{array}[]{c}G^{+}\\ \frac{1}{\sqrt{2}}(v+H_{1}^{0}+iG^{0})\\ \end{array}\right),\hskip 28.45274ptH_{2}\ =\ \left(\begin{array}[]{c}H_{2}^{+}\\ \frac{1}{\sqrt{2}}(H_{2}^{0}+iA)\\ \end{array}\right)\,,$$

$$\displaystyle V(H_{1},H_{2},\eta)\ =\$$

$$+\frac{1}{2}{\lambda }_1(H_1^{†}{H}_1{)}^2+\frac{1}{2}{\lambda }_2(H_2^{†}{H}_2{)}^2+{\lambda }_3(H_1^{†}{H}_1)(H_2^{†}{H}_2)+{\lambda }_4(H_1^{†}{H}_2)(H_2^{†}{H}_1)$$

$$+\left[\frac{1}{2}{\lambda }_5(H_1^{†}{H}_2{)}^2+\left\{{\lambda }_6(H_1^{†}{H}_1)+{\lambda }_7(H_2^{†}{H}_2)\right\}{H}_1^{†}{H}_2+H.c.\right]$$

$$+{\mu }_{\eta }^2|\eta {|}^2+{\lambda }_{\eta }|\eta {|}^4+{\lambda }_8|\eta {|}^2{H}_1^{†}{H}_1+{\lambda }_9|\eta {|}^2{H}_2^{†}{H}_2$$

$$+({\lambda }_{10}|\eta {|}^2{H}_1^{†}{H}_2+H.c.)+(\mu {ϵ}_{ij}{H}_1^i{H}_2^j{\eta }^{-}+H.c.)$$

$${\mu }_1^2=\frac{1}{2}{\lambda }_1{v}^2,{\mu }_3^2=\frac{1}{2}{\lambda }_6{v}^2,$$

$$\displaystyle M_{\rm charged}^{2}\$$

$$\displaystyle\quad M_{2}^{2}\$$

$$m_{h^{+},H^{+}}^2=\frac{1}{2}\left\{M_2^2+M_3^2\mp \sqrt{(M_2^2-M_3^2{)}^2+2v^2{\mu }^2}\right\},$$

$$h^{+}$$

$$H^{+}$$

$$\sin 2\phi =\frac{-\sqrt{2}v\mu }{m_{H^{+}}^2-m_{h^{+}}^2}.$$

Full text

††preprint: FERMILAB-PUB-19-304-T, OSU-HEP-19-04††arxiv: 1907.09498

1 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 Weinberg:1979sa

[TABLE]

suppressed by an inverse mass scale $\Lambda$ is induced by integrating out Standard Model (SM) singlet fermions Minkowski:1977sc ; Mohapatra:1979ia ; Yanagida:1979as ; GellMann:1980vs ; Glashow:1979nm , $SU(2)_{L}$ triplet scalars Schechter:1980gr ; Cheng:1980qt ; Mohapatra:1980yp ; Lazarides:1980nt , or $SU(2)_{L}$ triplet fermions Foot:1988aq with mass of order $\Lambda$.111For a clear discussion of the classification of seesaw types see Ref. Ma:1998dn . 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 $\epsilon_{ik}$ is the $SU(2)_{L}$ antisymmetric tensor. Once the vacuum expectation value (VEV) of the Higgs field, $\langle H^{0}\rangle\simeq 246$ GeV is inserted in Eq. (1.1), Majorana masses for the neutrinos given by $m_{\nu}=v^{2}/\Lambda$ will be induced. For light neutrino masses in the observed range, $m_{\nu}\sim(10^{-3}-10^{-1})$ eV, the scale $\Lambda$ should be around $10^{14}$ GeV. The mass of the new particle that is integrated out need not be $\Lambda$, since it is parametrically different, involving a combination of Yukawa couplings and $\Lambda$. For example, in the type-I seesaw model the heavy right-handed neutrino mass goes as $M_{R}\sim y_{D}^{2}\Lambda$, which can be near the TeV scale, if the Dirac Yukawa coupling $y_{D}\sim 10^{-6}$. However, it is also possible that $y_{D}\sim{\cal O}(1)$, in which case the new physics involved in neutrino mass generation could not be probed directly in experiments.222This is strictly true for one generation case. For more than one generation, the scale could be lower Kersten:2007vk .

An alternative explanation for small neutrino masses is that they arise only as quantum corrections Zee:1980ai ; Zee:1985id ; Babu:1988ki (for a review, see Ref. Cai:2017jrq ). In these radiative neutrino mass models, the tree-level Lagrangian does not generate ${\cal 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.333The magnitude of $m_{\nu}$ would be too small if it is induced at four or higher loops, assuming that the diagrams have chiral suppression factors proportional to the SM charged fermion masses; see Sec. 7.4. For example, in the two-loop neutrino mass model of Refs. Zee:1985id ; Babu:1988ki , small and calculable $m_{\nu}$ arises from the diagram shown in Fig. 43, which is estimated to be of order

[TABLE]

assuming normal ordering of neutrino masses and requiring large $\mu-\tau$ mixing. Here $f,h$ are Yukawa couplings involving new charged scalars with mass of order $M$. Even with $f\sim h\sim 1$, to obtain $m_{\nu}\sim 0.1$ eV, one would require the scalar mass $M\sim$ 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) Wolfenstein:1977ue induced by these mediators. These NSI are of great phenomenological interest, as their presence would modify the standard three-neutrino 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. Dev:2019anc ). 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) Gavela:2008ra . One possible way to avoid such constraints is to have light mediators for NSI Farzan:2015hkd ; Babu:2017olk ; Denton:2018xmq . 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.444Analysis of Ref. Forero:2016ghr ; Dey:2018yht of neutrino NSI in a model with charged singlet and/or doublet scalars, although not in the context of a neutrino mass model, is analogous to one model we analyze. 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.

1.1 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 Zee:1980ai 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 $(\Delta L=2)$ dimension 7 operator

[TABLE]

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 ${\cal 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. Babu:2001ex . 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. deGouvea:2007qla .

Another well known example in this category is the two-loop neutrino mass model of Refs. Zee:1985id ; Babu:1988ki , which induces an effective $d=9$ operator

[TABLE]

Neutrino masses arise in this model via the two-loop diagrams shown in Fig. 43, which has a chiral suppression factor proportional to $m_{\ell}^{2}$, with $\ell$ standing for the charged leptons of the SM.

This category of type-I radiative neutrino mass models is populated by one-loop, two-loop, and three-loop models. Popular one-loop type-I models include the Zee model Zee:1980ai (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 Hall:1983id . 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 Dorsner:2017wwn ; AristizabalSierra:2007nf . Ref. Cai:2014kra 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 Zee:1985id ; Babu:1988ki (cf. Sec. 7.2.1), a variant of it using LQs and a diquark (DQ) Kohda:2012sr (cf. Sec. 7.2.2), a pure LQ extension Babu:2010vp (cf. Sec. 7.2.3), a model with LQs and vector-like fermions Babu:2011vb (cf. Sec. 7.2.4), and the Angelic model Angel:2013hla (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 are proportional to the up-type quark masses. Type-I three-loop models include the KNT model Krauss:2002px (cf. Sec. 7.3.1), an LQ variant of the KNT model Nomura:2016ezz (cf. Sec. 7.3.4), the AKS model Aoki:2008av (cf. Sec. 7.3.2), and the cocktail model Gustafsson:2012vj (cf. Sec. 7.3.3). For a review of this class of models, see Ref. Cai:2017jrq .

A systematic approach to identify type-I radiative models is to start from a given $\Delta L=2$ effective operators of the type ${\cal 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. Babu:2001ex , and further developed in Refs. Angel:2012ug ; Cai:2014kra . We shall rely on these techniques. In particular, the many models suggested in Ref. Cai:2014kra 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. Klein:2019jgb .

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}\times SU(2)_{L}\times U(1)_{Y}$ quantum numbers, field components and electric charges (in superscript), and corresponding Lagrangian 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 $\Delta({\bf 1},{\bf 3},1)$ could induce neutrino mass at tree-level via type-II seesaw mechanism Schechter:1980gr ; Cheng:1980qt ; Mohapatra:1980yp ; Lazarides:1980nt , which makes radiative models involving $\Delta$ 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 Minkowski:1977sc ; Mohapatra:1979ia ; Yanagida:1979as ; GellMann:1980vs ; Glashow:1979nm , which should be suppressed by some mechanism or symmetry. Such a model has been analyzed in Refs. Pilaftsis:1991ug ; Dev:2012sg , 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. 59 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 $\Delta 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 ${\cal 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 Ma:2006km 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. 57. 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 Kubo:2006yx ; FileviezPerez:2009ud ; Law:2013saa ; Restrepo:2013aga ; Baek:2015mna ; Dutta:2018qei .

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

[TABLE]

Such an operator would lead to neutrino masses at the two-loop level, as can be seen in an explicit model shown in Fig. 58. Although this model can be described as arising from an effective $\Delta 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}\epsilon_{ik}\epsilon_{jl}(H^{\dagger}H)$, which is realized for instance in the minimal radiative inverse seesaw model (MRISM) of Ref. Dev:2012sg (see Sec. 7.1.1). Such effective operators, which appear as products of lower operators, were treated as trivial in the classification of Ref. Babu:2001ex .

1.2 Summary of results

We have mapped out in this paper the allowed ranges for the neutrino NSI parameters $\varepsilon_{\alpha\beta}$ (cf. Sec. 3) in radiative neutrino mass models. We present a detailed analysis of the Zee model Zee:1980ai with light charged scalar bosons (cf. Sec. 4). To map out the allowed values of $\varepsilon_{\alpha\beta}$, 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^{\pm}$ decays (cf. Sec. 4.8); vii) Lepton universality in $\tau$ decays (cf. Sec. 4.9); viii) Electroweak precision data (cf. Sec. 4.4); ix) charged-lepton flavor violation (cf. Sec. 4.5); x) Perturbative unitarity of Yukawa and quartic couplings; and xi) charge-breaking minima of the Higgs potential (cf. Sec. 4.3).

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 (cf. Table IX):

[TABLE]

These values are significantly larger than the ones obtained in Ref. Herrero-Garcia:2017xdu , where the contributions from the doublet Yukawa couplings of the light charged Higgs were ignored.

We have also analyzed in detail LQ models of radiative neutrino mass generation. As the base model we analyze the LQ version of the Zee model (cf. Sec. 5), the results of which can also be applied to other LQ models with minimal modifications. 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 (cf. Table XVII):

[TABLE]

[TABLE]

Our results yield somewhat larger NSI compared to the results of Ref. Wise:2014oea 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 (cf. Eq. (6.186)):

[TABLE]

For completeness, we also list here the maximum possible tree-level NSI in the two-loop Zee-Babu model (cf. Eq. (7.196)):

[TABLE]

The NSI predictions in all other models analyzed here will fall into one of the above categories (except for the MRISM discussed in Sec. 7.1.1). Our results for the base models mentioned above are summarized in Fig. 59, and the results for all the models analyzed in this paper are tabulated in Table XX. We emphasize that while our analysis is focused on radiative neutrino mass models, it essentially covers all NSI possibilities with heavy mediators, and thus is more general.

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 $\varepsilon_{\alpha\beta}$. 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. 59. In Appendix A, we present the analytic expressions for the charged-scalar production cross sections in electron-positron collisions.

2 Classification of $\Delta L=2$ operators and their UV completions

It is instructive to write down low-dimensional effective operators that carry lepton number of two units ($\Delta 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$ Babu:2001ex . We use two component Weyl notation for SM fermions and denote them as

[TABLE]

The Higgs field of the SM is denoted as $H\left({\bf 1},{\bf 2},\frac{1}{2}\right)$. The $\Delta L=2$ operators in the SM are all odd-dimensional. The full list of operators through $d=7$ is given by Babu:2001ex :

[TABLE]

Not listed here are products of lower-dimensional operators, such as $\mathcal{O}_{1}\times\overline{H}H$, with the $SU(2)_{L}$ contraction of $\overline{H}H$ being a singlet. Here $\mathcal{O}_{1}$ is the Weinberg operator Weinberg:1979sa , while the remaining operators are all $d=7$.555In the naming convention of Ref. Babu:2001ex , operators were organized based on how many fermion fields are in them. Operators $\mathcal{O}_{5}-\mathcal{O}_{7}$, which are $d=9$ operators, appeared ahead of the $d=7$ operator $\mathcal{O}_{8}$. 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 $\mathcal{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 $\mathcal{O}_{2}$ in Eq. (2.7b). There are two classes of models that can generate this operator with the respective mediators obtained from the following contractions (see Table II):

[TABLE]

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.666There is a third contraction allowed in principle, $e^{c}(LL)(LH)$. However, the mediator needed to realize this would generate $d=5$ operator $LLLH$ either via type-I or type-II seesaw at tree-level, and hence this contraction is not used in radiative neutrino mass models. Thus, $\mathcal{O}_{2}^{1}$ has two UV completions, with the addition of a vectorlike lepton $\psi\left({\bf 1},{\bf 2},-\frac{3}{2}\right)$ to the SM, along with a scalar which is either a singlet $\eta^{+}({\bf 1},{\bf 1},1)$, or a triplet $\Delta({\bf 1},{\bf 3},1)$. The choice of $\Delta({\bf 1},{\bf 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 $\mathcal{O}_{2}^{1}$ with $\psi\left({\bf 1},{\bf 2},-\frac{3}{2}\right)$ vectorlike lepton and $\eta^{+}({\bf 1},{\bf 1},1)$ scalar is discussed in Sec. 7.1.2. Similarly operator $\mathcal{O}_{2}^{2}$ has a unique UV completion, with two scalars added to the SM – one $\eta^{+}({\bf 1},{\bf 1},1)$ and one $\Phi\left({\bf 1},{\bf 2},\frac{1}{2}\right)$. This is the Zee model of neutrino mass, discussed at length in Sec. 4.

Operators $\mathcal{O}_{3a}$ and $\mathcal{O}_{3b}$ in Eq. (2.7c) can be realized by the UV complete models given in Table. III Cai:2014kra . Here all possible contraction among the fields are shown, along with the required mediators to achieve these contractions. Fields denoted as $\phi$ and $\eta$ are scalars, while $\psi$ 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 $\mathcal{O}_{4}$ and $\mathcal{O}_{8}$ are shown in Tables IV and V respectively Cai:2014kra . These models will all be analyzed in Sec. 7 for neutrino NSI. Note that in both $\mathcal{O}_{4}$ and $\mathcal{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 ${\cal O}_{4b}$ can be realized with scalar and fermionic exchange.

Classification based on topology of diagrams

Rather than classifying radiative neutrino mass models in terms of effective $\Delta L=2$ operators, one could also organize them in terms of the topology of the loop diagrams Ma:1998dn ; Bonnet:2012kz ; Sierra:2014rxa . Possible one-loop topologies are shown in Fig. 2 Ma:1998dn ; Bonnet:2012kz , and the two-loop topologies are shown in Fig. 3 Sierra:2014rxa . Note that in the two-loop diagrams, two Higgs particles that are connected to two internal lines in possible ways are not shown. Recently the three-loop topologies that generate operator ${\cal O}_{1}$ has been classified in Ref. Cepedello:2018rfh .

For the purpose of NSI, we find the classification based on type-I and type-II suggested here more 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.

3 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 low-energy regime, neutrino NSI with matter fields can be formulated in terms of an effective four-fermion Lagrangian as follows Wolfenstein:1977ue :

[TABLE]

where $G_{F}$ is the Fermi coupling constant, and $P_{X}$ (with $X=L,R$) denotes the chirality projection operators $P_{L,R}=(1\mp\gamma^{5})/2$. These projection operators can also be reparametrized into vector and axial components of the interaction. The dimensionless coefficients $\varepsilon_{\alpha\beta}$ are the NSI parameters that quantify the strength of the NSI between neutrinos of flavors $\alpha$ and $\beta$ and the matter fields $f,f^{\prime}\in\{e,u,d\}$. If $\varepsilon_{\alpha\beta}\neq 0$ for $\alpha\neq\beta$, the NSI violates lepton flavor, while for $\varepsilon_{\alpha\alpha}\neq\varepsilon_{\beta\beta}$, it violates lepton flavor universality.

The vector component of NSI, $\varepsilon_{\alpha\beta}^{fV}=\varepsilon_{\alpha\beta}^{fL}+\varepsilon_{\alpha\beta}^{fR}$, affects neutrino oscillations by providing a new flavor-dependent matter effect.777The axial-vector part of the weak interaction gives a nuclear spin-dependent contribution that averages to zero in the non-relativistic limit for the nucleus. The effective Hamiltonian for the matter effect is given by

[TABLE]

where $U_{\rm PMNS}$ is the standard $3\times 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+\varepsilon_{ee}$ term is due to the standard CC matter potential. The Hamiltonian level NSI in Eq. (3.11) is related to the Lagrangian level NSI in Eq. (3.9) as follows:

[TABLE]

where $N_{f}(x)$ is the number density of fermion $f$ at position $x$, and $\langle N_{p(n)}/N_{e}\rangle$ 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 $\langle N_{p}(x)/N_{e}(x)\rangle=1$ and define the ratio $Y_{n}(x)\equiv\langle N_{n}(x)/N_{e}(x)\rangle$ to rewrite Eq. (3.12) as

[TABLE]

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) Dziewonski:1981xy , $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 Serenelli:2009yc ; Coloma:2016gei .

In the following sections, we will derive the predictions for the NSI parameters $\varepsilon_{\alpha\beta}$ in various radiative neutrino mass models, which should then be compared with the experimental and/or global-fit constraints Coloma:2017egw ; Farzan:2017xzy ; Esteban:2018ppq ; Esteban:2019lfo on $\varepsilon_{\alpha\beta}$ using Eq. (3.13). We would like to emphasize two points in this connection:

1. (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.13) should be considered, while comparing with the experimental or global-fit constraints. 2. (ii)

Most of the experimental constraints Farzan:2017xzy 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 Esteban:2018ppq ; Esteban:2019lfo 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 both, as well as the future DUNE sensitivities deGouvea:2015ndi ; Coloma:2015kiu ; Blennow:2016etl ; dev_pondd (cf. 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 Zee:1980ai , 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 $\Delta L=2$ effective $d=7$ operator ${\cal O}_{2}=L^{i}L^{j}L^{k}e^{c}H^{l}\epsilon_{ij}\epsilon_{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 $\mathcal{O}_{2}^{2}$ of Table II.

The model assumes the SM gauge symmetry $SU(3)_{c}\times SU(2)_{L}\times U(1)_{Y}$, with an extended scalar sector. Two Higgs doublets $\Phi_{1,2}({\bf 1},{\bf 2},1/2)$, and a charged scalar singlet $\eta^{+}({\bf 1},{\bf 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:

[TABLE]

where $\{\alpha,\beta\}$ are generation indices, $\{i,j\}$ are $SU(2)_{L}$ indices, $\widetilde{\Phi}_{a}\equiv i\tau_{2}\Phi_{a}^{\star}$ ($a=1,2$) and $\ell^{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_{\alpha\beta}=-f_{\beta\alpha}$. Expanding the first term of the Lagrangian Eq. (4.14) leads to the following couplings of $\eta^{+}$:

[TABLE]

The presence of two Higgs doublets $\Phi_{1,2}$ allows for a cubic coupling in the Higgs potential,

[TABLE]

which, along with the Yukawa couplings of Eq. (4.14), would lead to lepton number violation. The magnitude of the parameter $\mu$ in Eq. (4.16) will determine the range of NSI allowed in the model. Interestingly, $\mu$ cannot be arbitrarily large, as it would lead to charge-breaking minima of the Higgs potential which are deeper than the charge conserving minimum  Barroso:2005hc ; Babu:2014kca (see Sec. 4.3).

4.1 Scalar sector

We can start with a general basis, where both $\Phi_{1}$ and $\Phi_{2}$ acquire vacuum expectation values (VEVs):

[TABLE]

However, without loss of generality, we can choose to work in the Higgs basis Davidson:2005cw where only one of the doublet fields gets a VEV $v$ given by $v=\sqrt{v_{1}^{2}+v_{2}^{2}}\simeq 246$ GeV. The transformation to the new basis $\{H_{1},H_{2}\}$ is given by:

[TABLE]

where $s_{\beta}\equiv\sin{\beta}$ and $c_{\beta}\equiv\cos{\beta}$, with $\tan\beta=v_{2}/v_{1}$. In this new basis, we can parametrize the two doublets as

[TABLE]

where $(G^{+},G^{0})$ are the Goldstone bosons, $(H_{1}^{0}$, $H_{2}^{0})$, $A$, and $H_{2}^{+}$ are the neutral $\mathcal{CP}$-even and odd, and charged scalar fields, respectively. We shall work in the $\mathcal{CP}$ conserving limit, since phases such as $\xi$ in Eq. (4.17) will not have a significant impact on NSI phenomenology which is our main focus here.

The most general renormalizable scalar potential involving the doublet fields $H_{1},H_{2}$ and the singlet field $\eta^{+}$ can be written as

[TABLE]

Differentiating $V$ with respect to $H_{1}$ and $H_{2}$, we obtain the following minimization conditions:

[TABLE]

where, for simplicity, we have chosen $\mu_{3}^{2}$ to be real. The mass matrix for the charged scalars in the basis $\{H_{2}^{+},\eta^{+}\}$ becomes

[TABLE]

where

[TABLE]

The physical masses of the charged scalars $\{h^{+},H^{+}\}$ are given by:

[TABLE]

where

[TABLE]

with the mixing angle $\varphi$ given by

[TABLE]

As we shall see later, this mixing parameter $\varphi$, which is proportional to $\mu$ will play a crucial role in the NSI phenomenology of the model.

Similarly, the matrix for the $\mathcal{CP}$-even and odd neutral scalars in the basis $\{H_{1}^{0},H_{2}^{0},A\}$ can written as Babu:2018uik :

[TABLE]

In the $\mathcal{CP}$-conserving limit where Im($\lambda_{5,6})=0$, the $\mathcal{CP}$-odd state will decouple from the $\mathcal{CP}$-even states. One can then rotate the $\mathcal{CP}$-even states into a physical basis $\{h,H\}$ which would have masses given by Babu:2018uik :

[TABLE]

whereas the $\mathcal{CP}$-odd scalar mass is given by

[TABLE]

The mixing angle between the $\mathcal{CP}$-even eigenstates $\{H_{1}^{0},H_{2}^{0}\}$, defined as

[TABLE]

is given by

[TABLE]

We will identify the lightest $\mathcal{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 $\beta-\alpha\rightarrow 0$ Gunion:2002zf ; Carena:2013ooa ; Dev:2014yca ; Das:2015mwa , as suggested by the LHC Higgs data Bernon:2015qea ; Chowdhury:2017aav .

4.2 Neutrino mass

In the Higgs basis where only the neutral component of $H_{1}$ gets a VEV, the Yukawa interaction terms in Eq. (4.14) of fermions with the scalar doublets $H_{1}$ and $H_{2}$ become

[TABLE]

where $Y$ and $\widetilde{Y}$ are the redefined couplings in terms of the original Yukawa couplings $y_{1}$ and $y_{2}$ given in Eq. (4.14) and where $\widetilde{H}_{a}=i\tau_{2}H^{\star}_{a}$ ($a=1,2)$ with $\tau_{2}$ being the second Pauli matrix. After electroweak symmetry breaking, the charged lepton mass matrix reads as

[TABLE]

Without loss of generality, one can work in a basis where $M_{\ell}$ is diagonal, i.e., $M_{\ell}=\text{diag}\,(m_{e},\,m_{\mu},\,m_{\tau})$. The Yukawa coupling matrix $f$ involving the $\eta^{+}$ field in Eq. (4.14) is taken to be defined in this basis.

The Yukawa couplings in Eq. (4.14), together with the trilinear term in the scalar potential Eq. (4.16), generate neutrino mass at the one-loop level, as shown in Fig. 4. Here the dot ($\bullet$) 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

[TABLE]

where $\kappa$ is the one-loop factor given by

[TABLE]

with $\varphi$ given in Eq. (4.28). From Eq. (4.36) it is clear that only the product of the Yukawa couplings $f$ and $Y$ is constrained by the neutrino oscillation data. Therefore, by taking some of the $Y$ couplings to be of $\sim{\cal O}(1)$ and all $f$ couplings very small in the neutrino mass matrix of Eq. (4.36), 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. With the other possibility, namely, $Y\ll 1$, the stringent cLFV constraints on $f$ couplings (cf. Table XIX) restrict the maximum NSI to $\lesssim 10^{-8}$ Herrero-Garcia:2017xdu , well below any foreseeable future experimental sensitivity.

The matrix $f$ that couples the left-handed lepton doublets to the charged scalar $\eta^{+}$ can be made real by a phase redefinition $\hat{P}f\hat{P}$, where $\hat{P}$ is a diagonal phase matrix, while the Yukawa coupling $Y$ in Eq. (4.34) is in general a complex asymmetric matrix:

[TABLE]

Here the matrix $Y$ is multiplied by $(\bar{\nu}_{e},\bar{\nu}_{\mu},\,\bar{\nu}_{\tau})$ from the left and $(e_{R},\,\mu_{R},\,\tau_{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_{\mu e},\,Y_{\tau e})$, which parametrize the couplings of $\nu_{\alpha}$ with electrons in matter.

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_{\alpha e}$ of Eq. (4.38) 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 $\pi^{+}\rightarrow e^{+}\nu$ will occur with unacceptably large rates. Therefore, we assume that the second Higgs doublet $H_{2}$ is leptophilic in our analysis.

Note that in the limit $Y\propto M_{l}$, as was suggested by Wolfenstein Wolfenstein:1980sy 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_{\nu}$ would vanish, yielding neutrino mixing angles that are not compatible with observations Koide:2001xy ; He:2003ih . For a variant of the Zee-Wolfenstein model with a family-dependent $Z_{4}$ symmetry which is consistent with neutrino oscillation data, see Ref. Babu:2013pma .

4.3 Charge-breaking minima

To have sizable NSI, we need a large mixing $\varphi$ between the singlet and doublet charged scalar fields $\eta^{+}$ and $H_{2}^{+}$. From Eq. (4.28), this means that we need a large trilinear $\mu$-term. But $\mu$ cannot be arbitrarily large, as it leads to charge-breaking minima (CBM) of the potential Barroso:2005hc ; Babu:2014kca . We numerically analyze the scalar potential given by Eq. (4.20) to ensure that it does not develop any CBM deeper than the charge-conserving minimum (CCM).

We take $\mu_{2}^{2},\,\mu_{\eta}^{2}>0$. The field $H_{1}$ is identified approximately as the SM Higgs doublet, and therefore, the value of $\lambda_{1}$ is fixed by the Higgs mass (cf. Eq. (4.21)), and the corresponding mass-squared term is chosen to be negative to facilitate electroweak symmetry breaking ($\mu_{1}^{2}>0$ in Eq. (4.20)). Note that the cubic scalar coupling $\mu$ can be made real as any phase in it can be absorbed in $\eta^{-}$ by a field redefinition.

In order to calculate the most general minima of the potential, we assign the following VEVs to the scalar fields:

[TABLE]

where $v_{\eta}$ and $v_{1}$ can be made real and positive by $SU(2)_{L}\times U(1)_{Y}$ rotations. A non-vanishing VEV $v_{\eta}$ would break electric charge conservation, as does a nonzero value of $\sin\gamma$. Thus, we must ensure that the CBM of the potential lie above the CCM. The Higgs potential, after inserting Eq. (4.39) in Eq. (4.20), reads as

[TABLE]

Here $\theta_{1},\theta_{2},\theta_{3},$ and $\theta_{4}$ are respectively the phases of the quartic couplings $\lambda_{5},\lambda_{6},\lambda_{7},$ and $\lambda_{10}$. For simplicity, we choose these quartic couplings, as well as $\lambda_{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 $\lambda_{3}$ and $\lambda_{8}$ nonzero, as these couplings help in satisfying constraints from the SM Higgs boson properties from the LHC.

Eq. (4.40) yields five minimization conditions from which $\{v_{1},v_{2},v_{\eta},\delta,\gamma\}$ 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 $\varphi$ (cf. Eq. (4.28)] for our subsequent analysis, we keep $\lambda_{3}=\lambda_{8}$ fixed at two benchmark values (3.0 and 2.0) and vary the remaining nonzero quartic couplings $\lambda_{2}$ and $\lambda_{\eta}$ in the range $[0.0,3.0]$. Our results on the maximum $\sin\varphi$ are shown in Fig. 5. We do not consider values of the quartic couplings exceeding 3.0 to be consistent with perturbativity considerations Kang:2013zba . Each choice of mixing angle $\varphi$, and the parameters $\lambda_{2}$, $\lambda_{\eta}$, $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\varphi$ 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 $\varphi$, the trilinear coupling $\mu$ becomes too large, which drives the potential to a deeper CBM Barroso:2005hc , even for positive $\mu_{\eta}^{2}$. From Fig. 5 it is found that $\sin\varphi<0.23$ for $m_{H^{+}}=2~{}{\rm TeV}$, while $\sin\varphi=0.707$ is allowed when $m_{H^{+}}=0.7$ TeV. In all cases the maximum value of $|\mu|$ is found to be about 4.1 times the heavier mass $m_{H^{+}}$. Note that we have taken the maximum value of the mixing $\varphi_{\rm max}=\pi/4$ here, because for $\varphi>\pi/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.

4.4 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 Peskin:1990zt ; Peskin:1991sw and impose important constraints from precision data. These parameters have been calculated in the context of the Zee model in Ref. Grimus:2008nb . 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 Grimus:2008nb :

[TABLE]

where the symmetric function $\mathcal{F}$ is given by

[TABLE]

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 $\varphi$, 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\pm 0.12$ Tanabashi:2018oca . For simplicity, we assume no mixing between the neutral $\mathcal{CP}$-even scalars $h$ and $H$. Furthermore, we take the heavy neutral $\mathcal{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^{\pm}}$ 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^{\pm}}$, there is an upper limit on the mixing $\sin{\varphi}$. This is further illustrated in Fig. 7. Here, the colored regions (both green and red) depict the allowed parameter space in $m_{H}^{+}-\sin{\varphi}$ 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{\varphi})_{\rm max}=0.63$. The mass splitting between the heavy neutral and the charged Higgs bosons is governed by the relation (cf. Eqs. (4.26) and (4.30)):

[TABLE]

We choose $\lambda_{5}=-\lambda_{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 $|\lambda_{4}|,|\lambda_{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{\varphi}<0.59$. For the degenerate case $m_{H^{\pm}}=m_{H}$ with $\lambda_{4}=\lambda_{5}$, the upper limit is stronger: $\sin{\varphi}<0.49$.

4.5 Charged-lepton flavor violation constraints

Charged-lepton flavor violation is an integral feature of the Lagrangian Eq. (4.14) of the model. We can safely ignore cLFV processes involving the $f_{\alpha\beta}$ couplings which are assumed to be of the order of $10^{-8}$ or so to satisfy the neutrino mass constraint, with $Y_{\alpha\beta}$ couplings being order one. Thus, we focus on cLFV proportional to $Y_{\alpha\beta}$. Furthermore, as noted before, NSI arise proportional to $(Y_{ee},\,Y_{\mu e},\,Y_{\tau 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_{\alpha\beta}$, we shall focus on the set $(Y_{ee},\,Y_{\mu e},\,Y_{\tau e})$ relevant for NSI. The neutral scalar bosons $H$ and $A$ will mediate cLFV of the type $\mu\rightarrow 3e$ and $\tau\rightarrow\mu ee$ at tree-level, while these neutral scalars and the charged scalars $(h^{+},\,H^{+})$ mediate processes of the type $\mu\rightarrow e\gamma$ via one-loop diagrams. Both of these processes will be analyzed below. We derive limits on the couplings $Y_{\alpha\beta}$ 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 limits of possible off-diagonal NSI. The various processes considered and the limits derived are summarized in Tables VI and VII. We now turn to the derivation of these bounds.

4.5.1 $\ell_{\alpha}\to\ell_{\beta}+\gamma$ decays

The decay $\ell_{\alpha}\to\ell_{\beta}+\gamma$ arises from one-loop diagrams shown in Fig. 8. The general expression for this decay rate can be found in Ref. Lavoura:2003xp . Let us focus on the special case where the FCNC coupling matrix $Y$ of Eq. (4.38) 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_{\alpha}$ has nonzero Yukawa couplings with a scalar boson $B$ and fermion $F$, parametrized as

[TABLE]

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_{\alpha}\rightarrow f_{\beta}+\gamma$ is then given by

[TABLE]

Here $\alpha=e^{2}/4\pi$ is the fine-structure constant, $t=m_{F}^{2}/m_{B}^{2}$, and the function $f_{F}(t)$ and $f_{B}(t)$ are given by

[TABLE]

These expressions are obtained in the approximation $m_{\beta}\ll m_{\alpha}$.

Let us apply these results to $\ell_{\alpha}\rightarrow\ell_{\beta}+\gamma$ mediated by the charged scalars ($h^{+},\,H^{+})$ in the Zee model where the couplings have the form $Y_{\alpha\beta}\bar{\nu}_{\alpha}P_{R}\ell_{\beta}h^{+}\sin\varphi$, etc. Here $Q_{F}=0$, while $Q_{B}=+1$. Eq. (4.45) then reduces to (with $t\ll 1$)

[TABLE]

If we set $m_{h^{+}}=100$ GeV, $m_{H^{+}}=700$ GeV and $\sin\varphi=0.7$, then the experimental limit BR($\mu\rightarrow e\gamma)\leq 4.2\times 10^{-13}$ Cei:2017puj implies $|Y_{\alpha e}Y_{\alpha\mu}^{\star}|\leq 6\times 10^{-5}$. Similarly, the limit BR($\tau\rightarrow e\gamma)\leq 3.3\times 10^{-8}$ Aubert:2009ag implies $|Y_{\alpha\tau}Y_{\alpha e}^{\star}|\leq 4\times 10^{-2}$, and the limit BR($\tau\rightarrow\mu\gamma)\leq 4.4\times 10^{-8}$ Aubert:2009ag implies $|Y_{\alpha\tau}Y_{\alpha\mu}^{\star}|\leq 4.6\times 10^{-2}$. These are rather stringent constraints, which suggest that no more than one entry in a given row of $Y$ can be large. Such a choice does not however affect the maximum NSI, as the elements of $Y$ that generate them are in the first column of $Y$. Keeping only the entries $(Y_{ee},\,Y_{\mu e},\,Y_{\tau e})$ nonzero does not lead to $\ell_{\alpha}\rightarrow\ell_{\beta}+\gamma$ decay mediated by the charged scalars ($h^{+},\,H^{+})$.

However, nonzero values of $(Y_{ee},\,Y_{\mu e},\,Y_{\tau e})$, needed for NSI, would lead to $\ell_{\alpha}\rightarrow\ell_{\beta}+\gamma$ mediated by the heavy neutral scalars. Taking $H$ and $A$ to be degenerate, the Yukawa couplings are of the form $\bar{\ell}_{\alpha}P_{R}\ell_{\beta}H$. Thus, $Q_{F}=-1$ and $Q_{B}=0$ in this case, leading to the decay width

[TABLE]

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 $\mu\to e\gamma$ process prevents the off-diagonal NSI parameter $\varepsilon_{e\mu}$ from being in the observable range. However, $\varepsilon_{e\tau}$ and $\varepsilon_{\mu\tau}$ can be in the observable range, consistent with these constraints.

4.5.2 Electron anomalous magnetic moment

Another potential constraint comes from anomalous magnetic moment of leptons $(g-2)_{\alpha}$, 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 Babu:2002uu

[TABLE]

Comparing this with $\Delta a_{e}\equiv a_{e}^{\rm exp}-a_{e}^{\rm SM}=(-87\pm 36)\times 10^{-14}$ (where $a_{e}\equiv(g-2)_{e}/2$), based on the difference between the experimental measurements Hanneke:2010au and SM calculations Aoyama:2017uqe with the updated value of the fine-structure constant Parker:2018vye , we find that the charged scalar contribution (4.49) goes in the right direction. However, for the allowed parameter space in $m_{h^{+}}-Y_{ee}\sin\varphi$ plane (see Fig. 18), it turns out to be too small to explain the $2.4\sigma$ discrepancy in $\Delta a_{e}$. For example, with $|Y_{\tau e}|\sin\varphi=0.75$ and $m_{h^{+}}=150$ GeV, which is a consistent choice (cf. Fig. 18), we would get $\Delta a_{e}=-2.2\times 10^{-14}$, an order of magnitude too small to be relevant for experiments.

4.5.3 $\ell_{\alpha}\to\bar{\ell}_{\beta}\ell_{\gamma}\ell_{\delta}$ decays

The Yukawa coupling matrix $Y$ of the second Higgs doublet (cf. Eq. (4.38)) would lead to trilepton decay of charged leptons mediated by the neutral scalars of the theory. The tree-level 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 $\mu\rightarrow\bar{e}ee$ is given as follows:

[TABLE]

The partial decay width for $\tau\rightarrow\bar{\ell}_{\alpha}\ell_{\beta}\ell_{\gamma}$ is given by

[TABLE]

Here $S$ = 1 (2) for $\beta=\gamma\,(\beta\neq\gamma)$ is a symmetry factor. Using the total muon and tau decay widths, $\Gamma^{\rm tot}_{\mu}=3.00\times 10^{-19}$ GeV and $\Gamma^{\rm tot}_{\tau}=2.27\times 10^{-12}$ GeV respectively, we calculate the cLFV branching ratios for the processes $\mu^{-}\to e^{+}e^{-}e^{-}$, $\tau^{-}\to e^{+}e^{-}e^{-}$ and $\tau^{-}\to e^{+}e^{-}\mu^{-}$ using Eqs. (4.50) and (4.51). We summarize in Table VII the current experimental bounds on these branching ratios and the constraints on the Yukawa couplings $Y_{\alpha\beta}$ 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 $\ell_{\alpha}\to\ell_{\beta}\gamma$ decays. This also implies that off-diagonal NSI are severely constrained.

As already noted, the light charged Higgs $h^{+}$ would mediate $\ell_{\alpha}\rightarrow\ell_{\beta}+\gamma$ 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_{\tau\tau}$ and $Y_{\mu\tau}$ will not lead to tree-level trilepton decay of $\tau$. Apart from the first column of $Y$, we shall allow nonzero entries in the third column as well. In particular, for diagonal NSI $\varepsilon_{\alpha\alpha}$, we need one $Y_{\alpha e}$ entry for some $\alpha$ to be nonzero, and to avoid the trilepton constraints, the only other entry that can be allowed to be large is $Y_{\beta\tau}$ with $\beta\neq\alpha$. On the other hand, for off-diagonal NSI $\varepsilon_{\alpha\beta}$ (with $\alpha\neq\beta$), we must allow for both $Y_{\alpha e}$ and $Y_{\beta e}$ to be non-zero. In this case, however, the trilepton decay $\ell_{\beta}\to\ell_{\alpha}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_{\gamma\tau}$ with $\gamma\neq\alpha,\beta$. This will lead to $\tau\rightarrow\ell+\gamma$ 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.

4.6 Collider constraints on neutral scalar mass

In this section, we discuss the collider constraints on the neutral scalars $H$ and $A$ in the Zee model from various LEP and LHC searches.

4.6.1 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 LEP:2003aa . Integrating out new particles in a theory one can express their effect via higher-dimensional (generally dimension-6) operators. An effective Lagrangian, $\mathcal{L}_{\rm eff}$, can parametrize the contact interaction for the process $e^{+}e^{-}\to f\bar{f}$ with the form Eichten:1983hw

[TABLE]

where $\delta_{ef}$ is the Kronecker delta function, $f$ refers to the final sate fermions, $g$ is the coupling strength, $\Lambda$ is the new physics scale and $\eta^{f}_{ij}=\pm 1$ or 0, depending on the chirality structure. LEP has put 95% confidence level (CL) lower limits on the scale of the contact interaction $\Lambda$ assuming the coupling $g=\sqrt{4\pi}$ LEP:2003aa . 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^{-}\to\ell_{\alpha}^{+}\ell_{\beta}^{-}$ (with $\ell_{\alpha,\beta}=e,\mu,\tau$), and therefore, the LEP constraints on $\Lambda$ 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 $\ell_{\alpha}^{+}\ell_{\beta}^{-}\to\ell_{\gamma}^{+}\ell_{\delta}^{-}$ via heavy neutral scalar exchange, the effective Lagrangian in the Zee model can be written as

[TABLE]

By Fierz transformation, we can rewrite it in a form similar to Eq. (4.52):

[TABLE]

Thus, the only relevant chirality structures in Eq. (4.52) are $LR$ and $RL$, and the relevant process for deriving the LEP constraints is $e^{+}e^{-}\to\ell_{\alpha}^{+}\ell_{\alpha}^{-}$:

[TABLE]

with $\eta^{\ell}_{LR}=\eta^{\ell}_{RL}=-1$.

Now for $e^{+}e^{-}\to e^{+}e^{-}$, Eq. (4.54) becomes

[TABLE]

Comparing this with Eq. (4.55), we obtain

[TABLE]

where $\Lambda^{-}$ corresponds to $\Lambda$ with $\eta^{\ell}_{LR}=\eta^{\ell}_{RL}=-1$. The LEP constraints on $\Lambda$ were derived in Ref. LEP:2003aa for $g=\sqrt{4\pi}$, which can be translated into a lower limit on $m_{H}/|Y_{ee}|$ using Eq. (4.57), as shown in Table VIII. Similarly, for $e^{+}e^{-}\to\mu^{+}\mu^{-}$, Eq. (4.54) is

[TABLE]

Since for NSI, only $Y_{\mu e}$ (neutrino interaction with electron) is relevant, we can set $Y_{e\mu}\to 0$, and compare Eq. (4.58) with Eq. (4.55) to get a constraint on $m_{H}/|Y_{\mu e}|$, as shown in Table VIII. Similarly, for $e^{+}e^{-}\to\tau^{+}\tau^{-}$, we can set $Y_{e\tau}\to 0$ and translate the LEP limit on $\Lambda^{-}$ into a bound on $m_{H}/|Y_{\tau e}|$, as shown in Table VIII.

The LEP constraints from the processes involving $q\bar{q}$ final states, such as $e^{+}e^{-}\to c\bar{c}$ and $e^{+}e^{-}\to b\bar{b}$, 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.

4.6.2 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^{-}\to\ell_{\alpha}^{+}\ell_{\alpha}^{-}$ 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 Christensen:2008py and compute the $e^{+}e^{-}\to\ell_{\alpha}^{+}\ell_{\alpha}^{-}$ cross-sections in the Zee model at the parton-level using MadGraph5 event generator Alwall:2014hca . These numbers are then compared with the measured cross sections LEP:2003aa ; Abbiendi:2003dh to derive limits on $m_{H/A}$ as a function of the Yukawa couplings $Y_{\alpha e}$ (for $\alpha=e,\mu,\tau$). For a benchmark value of $m_{H}=m_{A}=130$ GeV, we find the following constraints on the Yukawa couplings $Y_{\alpha e}$ relevant for NSI:

[TABLE]

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).

4.6.3 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\to Z^{\star}\to HA\to\ell^{+}\ell^{-}\ell^{+}\ell^{-}$.888Only the $(H\overset{\leftrightarrow}{\partial}_{\mu}A)Z^{\mu}$ coupling is nonzero, while the $(H\overset{\leftrightarrow}{\partial}_{\mu}H)Z^{\mu}$ and $(A\overset{\leftrightarrow}{\partial}_{\mu}A)Z^{\mu}$ couplings vanish due to parity Grzadkowski:2018ohf . Most of the LHC multilepton searches assume a heavy $ZZ^{(\star)}$ resonance Sirunyan:2018qlb ; Aaboud:2018puo , 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 Chatrchyan:2014aea ; Aaboud:2018zeb , turn out to be weaker than the LEP constraints derived above.

4.7 Collider constraints on light charged scalar

In this section, we discuss the collider constraints on the light charged scalar $h^{\pm}$ in the Zee model from various LEP and LHC searches.

4.7.1 Constraints from LEP searches

At LEP, $h^{\pm}$ can be pair-produced through the $s$-channel Drell-Yan process mediated by either $\gamma$ or $Z$ boson (see Fig. 10 (a)). It can also be pair-produced through the $t$-channel processes mediated by a light neutrino (see Fig. 10 (b)). In addition, it can be singly-produced 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)). The analytic expressions for the relevant cross sections can be found in Appendix A. For our numerical study, we implement our model file in FeynRules package Christensen:2008py and compute all the cross-sections at the parton-level using MadGraph5 event generator Alwall:2014hca . We find good agreement between the numerically computed values and the analytic results presented in Appendix A.

Once produced on-shell, the charged scalar will decay into the leptonic final states $\nu_{\alpha}\ell_{\beta}$ through the Yukawa coupling $Y_{\alpha\beta}$. 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_{\alpha e}$ and $Y_{\alpha\mu}$ (or $Y_{\alpha e}$ and $Y_{\beta e}$) cannot be large simultaneously. So we consider the case where ${\rm BR}_{e\nu}+{\rm BR}_{\tau\nu}=100\%$ and ${\rm BR}_{\mu\nu}$ is negligible, in order to avoid more stringent limits from muon decay.999This choice is consistent with the observed neutrino oscillation data (see Sec. 4.13).

Electron channel: For a given charged scalar decay branching ratio to electrons, ${\rm BR}_{e\nu}$, we can reinterpret the LEP selectron searches lepsusy to put a constraint on the charged scalar mass as a function of ${\rm BR}_{e\nu}$. In particular, the right-handed selectron pair-production $e^{+}e^{-}\to\widetilde{e}_{R}\widetilde{e}_{R}$, followed by the decay of each selectron to electron and neutralino, $\widetilde{e}_{R}\to e_{R}+\widetilde{\chi}^{0}$, will mimic the $e^{+}e^{-}\nu\bar{\nu}$ final state of our case in the massless neutralino limit. So we use the 95% CL observed upper limits on the $\widetilde{e}_{R}\widetilde{e}_{R}$ production cross section lepsusy for $m_{\widetilde{\chi}}=0$ as an experimental upper limit on the quantity

[TABLE]

and derive the LEP exclusion region in the plane of charged scalar mass and ${\rm BR}_{e\nu}$, as shown in Fig. 11 (a) by the orange-shaded region. Here we have chosen $Y_{ee}\sin\varphi=0.1$ and varied $Y_{\tau\alpha}$ (with $\alpha=\mu$ or $\tau$) to get the desired branching ratios. We find that for ${\rm BR}_{e\nu}=1$, charged scalar masses less than 100 GeV are excluded. For ${\rm BR}_{e\nu}<1$, these limits are weaker, as expected, and the charged scalar could be as light as 97 GeV (for ${\rm BR}_{e\nu}=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\varphi=0.2$. The LEP selectron constraints become stronger as we increase $Y_{ee}$ and extend to smaller ${\rm BR}_{e\nu}$. However, the mass limit of 100 GeV for ${\rm BR}_{e\nu}=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 $\sqrt{s}=209$ GeV and is already below the experimental cross section limit even for $Y_{ee}\sin\varphi=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_{\tau e}\sin\varphi$ and varied $Y_{\tau\alpha}$ (with $\alpha=e$ or $\mu$) 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_{\alpha e}$ (for $\alpha=e,\mu$ or $\tau$). 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.101010This might be the reason why the LEP limits derived here are somewhat more stringent than those reported in Ref. Cao:2017ffm , which presumably only considered the $s$-channel contribution.

Tau channel: In the same way, we can also use the LEP stau searches lepsusy to derive an upper limit on

[TABLE]

and the corresponding LEP exclusion region in the plane of charged scalar mass and ${\rm BR}_{\tau\nu}$, as shown in Fig. 11 by the blue-shaded region. We find that for ${\rm BR}_{\tau\nu}=1$, charged scalar masses less than 104 (105) GeV are excluded for $Y_{ee}\sin\phi=0.1~{}(0.2)$.

For ${\rm BR}_{\tau\nu}\neq 0$, a slightly stronger limit can be obtained from the LEP searches for the charged Higgs boson pairs in the 2HDM Abbiendi:2013hk . Their analysis focused on three kinds of final states, namely, $\tau\nu\tau\nu$, $c\bar{s}\tau\nu$ (or $\bar{c}s\tau\nu$) and $c\bar{s}\bar{c}s$, under the assumption that ${\rm BR}_{\tau\nu}+{\rm BR}_{c\bar{s}}=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 $\sigma(e^{+}e^{-}\to h^{+}h^{-}){\rm BR}^{2}_{\tau\nu}$ for ${\rm BR}_{\tau\nu}=1$ can be recast into an upper limit on

[TABLE]

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 $c\bar{s}\tau\nu$ for ${\rm BR}_{\tau\nu}\neq 1$ as an upper limit on $\sigma(e^{+}e^{-}\to h^{\pm}W^{\mp}){\rm BR}_{\tau\nu}{\rm BR}_{W\to c\bar{s}}$ and the corresponding exclusion region is shown in Fig. 11 by the cyan-shaded region, which is found to be weaker than the $\tau\nu\tau\nu$ mode.

4.7.2 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\to\gamma^{\star}/Z^{\star}\to h^{+}h^{-}$ (see Fig. 12 (a)), followed by the leptonic decay of $h^{\pm}\to\ell\nu$. There are also single-production processes as shown in Fig. 12 (b)-(d), which are less important. The relevant LHC searches are those for right-handed selectrons/staus: $pp\to\widetilde{\ell}^{+}_{R}\widetilde{\ell}^{-}_{R}\to\ell^{+}_{R}\widetilde{\chi}^{0}\ell^{-}_{R}\widetilde{\chi}^{0}$, which will mimic the $\ell^{+}\nu\ell^{-}\nu$ final states from $h^{+}h^{-}$ decay in the massless neutralino limit. The $\sqrt{s}=13$ TeV LHC stau searches focus on the stau mass range above 100 GeV and it turns out that the current limits Sirunyan:2018vig on the stau pair-production 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 $\sqrt{s}=8$ TeV ATLAS analysis considered the mass range down to 80 GeV Aad:2014yka ; however, the observed cross section is still found to be larger than the theoretical prediction in our case even for ${\rm BR}_{\tau\nu}=1$.

As for the selectron case, we take the $\sqrt{s}=13$ TeV CMS search Sirunyan:2018nwe , which focuses on the selectron masses above 120 GeV, and use the observed cross section limit on $\sigma(pp\to e^{+}_{R}\widetilde{\chi}^{0}e^{-}_{R}\widetilde{\chi}^{0})$ to derive an upper limit on $\sigma(pp\to h^{+}h^{-}){\rm BR}^{2}_{e\nu}$, 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 ${\rm BR}_{\tau\nu}\gtrsim 0.4$, which can always be done for any given Yukawa coupling $Y_{\alpha e}$ by choosing an appropriate $Y_{\beta\tau}$. This however may not be the optimal choice for NSI, especially for $Y_{ee}\neq 0$, where the lepton universality constraints restrict us from having a larger ${\rm BR}_{\tau\nu}$. Thus, the LHC constraints will be most relevant for $\varepsilon_{ee}$, as we will see in Fig. 18 (a).

4.8 Constraints from lepton universality in $W^{\pm}$ 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^{\pm}W^{\mp}$ (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 ALEPH:2004aa , the constraints on the ratio of $W$ branching ratios to leptons of different flavors are as follows:

[TABLE]

Note that while the measured value of $R_{\mu/e}$ agrees with the lepton universality prediction of the SM, $R_{\mu/e}^{\rm SM}=1$, within $1.1\sigma$ CL, the $W$ branching ratio to tau with respect to electron is about $1.8\sigma$ and to muon is about $2.7\sigma$ away from the SM prediction: $R_{\tau/\ell}^{\rm SM}=0.9993$ (with $\ell=e,\mu$), using the one-loop calculation of Ref. Kniehl:2000rb .

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^{-}\to W^{+}W^{-}\to\ell\nu q\bar{q}^{\prime}$ ALEPH:2004aa . However, due to the leptophilic nature of the charged Higgs $h^{\pm}$ in our model, neither the $e^{+}e^{-}\to h^{+}h^{-}$ channel (Figs. 10 (a) and 10 (b)) nor the Drell-Yan single-production channel (Fig. 10 (d)) will lead to $\ell\nu q\bar{q}$ final state. So the only relevant contribution to the $W$ universality violation could come from the $h^{\pm}W^{\mp}$ production channel (Fig. 10 (c)), with the $W$ decaying hadronically and $h^{\pm}$ decaying leptonically. The pure leptonic channels ($e\nu e\nu$ and $\mu\nu\mu\nu$) have $\sim 40\%$ uncertainties in the measurement and are therefore not considered here.

Including the $h^{\pm}W^{\mp}$ contribution, the modified ratios $R_{\ell/\ell^{\prime}}$ can be calculated as follows:

[TABLE]

where $\sigma(W^{+}W^{-})$ and $\sigma(h^{\pm}W^{\mp})$ are the production cross sections for $e^{+}e^{-}\to W^{+}W^{-}$ and $e^{+}e^{-}\to h^{\pm}W^{\mp}$ respectively, ${\rm BR}^{W}_{\ell\nu}$ denotes the branching ratio of $W\to\ell\nu$ (with $\ell=e,\mu,\tau$), whereas ${\rm BR}_{\ell\nu}$ denotes the branching ratio of $h^{\pm}\to\ell\nu$ as before (with $\ell=e,\tau$). At LEP experiment, the $W^{+}W^{-}$ pair production cross section $\sigma_{W^{+}W^{-}}$ is computed to be 17.17 pb at $\sqrt{s}=209$ GeV ALEPH:2004aa . Within the SM, $W^{\pm}$ decays equally to each generation of leptons with branching ratio of $10.83\%$ and decays hadronically with branching ratio of $67.41\%$ Tanabashi:2018oca . We numerically compute using MadGraph5 Alwall:2014hca the $h^{\pm}W^{\mp}$ cross section at $\sqrt{s}=209$ GeV as a function of $m_{h^{\pm}}$ and ${\rm BR}_{\ell\nu}$, and compare Eq. (4.66) with the measured values given in Eqs. (4.63)-(4.65) to derive the $2\sigma$ exclusion limits in the $m_{h^{+}}$-${\rm BR}_{\ell\nu}$ plane. This is shown in Figs. 11 (a) and 11 (b) by yellow, brown, and pink-shaded regions for $\mu/e$, $\tau/e$, and $\tau/\mu$ 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^{\pm}h^{\mp}$ 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^{+}}$-${\rm BR}_{\ell\nu}$ plane shown in Figs. 11 (a) and 11 (b), and rule out charged scalars below 110 GeV (129 GeV) for $Y_{ee}\sin\varphi=0.1$ (0.2). These constraints are not applicable for $m_{h^{\pm}}>129$ GeV, because $h^{\pm}W^{\mp}$ can no longer be produced on-shell at LEP II with maximum $\sqrt{s}=209$ GeV.

As mentioned before, the measured $W$ branching ratio to tau with respect to muon is $2.7\sigma$ above the SM prediction. Since in our case, $h^{\pm}$ decays to either $e\nu$ or $\tau\nu$, but not $\mu\nu$, this contributes to $R_{\tau\mu}$ only in the numerator, but not in the denominator. Therefore, the $2.7\sigma$ 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\varphi=0.1$.111111Light charged scalar has been used to address the lepton universality issue in $W$ decays in Ref. Dermisek:2008dq . The upper pink-shaded region with larger ${\rm BR}_{\tau\nu}$ gives $R_{\tau\mu}>1.122$, which is above the allowed $2\sigma$ range given in Eq. (4.65). On the other hand, the lower pink-shaded region with smaller ${\rm BR}_{\tau\nu}$ gives $R_{\tau\mu}<1.018$, which is below the allowed $2\sigma$ range given in Eq. (4.65). For larger Yukawa coupling $Y_{ee}$, as illustrated in Fig. 11 (b) with $Y_{ee}\sin\varphi=0.2$, the whole allowed range of parameter space from $R_{\tau/\mu}$ shifts to lower values of ${\rm BR}_{\tau\nu}$. This is because the $h^{\pm}W^{\mp}$ production cross section $\sigma(h^{\pm}W^{\mp})$ in Eq. (4.66) is directly proportional to $|Y_{ee}|^{2}$, and therefore, for a large $Y_{ee}$, a smaller ${\rm BR}_{\tau\nu}$ would still be compatible with the $R_{\tau/\mu}$-preferred range.

4.9 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^{-}\rightarrow\tau\bar{\nu}_{\beta}$ proceeds with a significant branching ratio. $h^{-}$ also has coupling with $e\bar{\nu}_{\alpha}$, so that non-negligible NSI is generated. When these two channels are combined, we would get new decay modes for the $\tau$ lepton, as shown in Fig. 13. This will lead to deviation in $\tau$-lifetime compared to the SM expectation. The new decay modes will also lead to universality violation in $\tau$ 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 $\tau$ decay mode is given by

[TABLE]

This can be recast, after a Fierz transformation, as

[TABLE]

This can be directly compared with the SM $\tau$ decay Lagrangian, given by

[TABLE]

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 $\tau$ lepton is now increased from its SM value by a factor $1+\Delta$, with $\Delta$ given by Kuno:1999jp

[TABLE]

where

[TABLE]

The global-fit result on $\tau$ lifetime is $\tau_{\tau}=(290.75\pm 0.36)\times 10^{-15}$ s, while the SM prediction is $\tau_{\tau}^{\rm SM}=(290.39\pm 2.17)\times 10^{-15}$ s Tanabashi:2018oca . Allowing for 2$\sigma$ error, we find $\Delta\leq 1.5\%.$ If the only decay modes of $h^{-}$ are $h^{-}\to\bar{\nu}_{\alpha}e^{-}$ and $h^{-}\to\bar{\nu}_{\beta}\tau^{-}$, then we can express $|Y_{\beta\tau}|^{2}$ in terms of $|Y_{\alpha e}|^{2}$ as

[TABLE]

Using this relation, we obtain

[TABLE]

where $\varepsilon_{\alpha\alpha}$ is the diagonal NSI parameter for which the expression is derived later in Eq. (4.89). Therefore, a constraint on $\Delta$ from the tau lifetime can be directly translated into a constraint on $\varepsilon_{\alpha\alpha}$:

[TABLE]

An even stronger limit is obtained from $e-\mu$ universality in $\tau$ decays. The experimental central value prefers a slightly larger width for $\tau\to\mu\nu\nu$ compared to $\tau\to e\nu\nu$. In our scenario, $h^{-}$ mediation enhances $\tau\to e\nu\nu$ relative to $\tau\to\mu\nu\nu$. We have in this scenario

[TABLE]

which constrains $\Delta\leq 0.002$, obtained by using the measured ratio $\frac{\Gamma(\tau\to\mu\nu\nu)}{\Gamma(\tau\to e\nu\nu)}=0.9762\pm 0.0028$ Tanabashi:2018oca , and allowing 2$\sigma$ error. This leads to a limit

[TABLE]

In deriving the limits on a light charged Higgs mass from LHC constraints, we have imposed the $\tau$ decay constraint as well as the universality constraint on $\Delta$, see Fig. 11. Avoiding the universality constraint by opening up the $\tau\to\mu\nu\nu$ channel will not work, since that will be in conflict with $\mu\to e\nu\nu$ constraints, which are more stringent.

The Michel parameters in $\tau$ decay will now be modified Babu:2016fdt . While the $\rho$ and $\delta$ parameters are unchanged compared to their SM value of 3/4, $\xi$ is modified from its SM value of 1 to

[TABLE]

However, the experimental value is $\xi=0.985\pm 0.030$ Tanabashi:2018oca , which allows for significant room for the new decay. Again, our choice of Yukawa couplings does not modify the $\mu\to e\nu\nu$ decay, and is therefore, safe from the Michel parameter constraints in the muon sector, which are much more stringent.

4.10 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 Khachatryan:2016vau and Run II Sirunyan:2018koj ; ATLAS:2019slw . 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\to\gamma\gamma$ decay (see Fig. 14 (a)) and mimics the tree-level $h\to WW^{\star}\to 2\ell 2\nu$ channel via the exotic decay mode $h\to h^{\pm}h^{\mp\star}\to h^{\pm}\ell\nu\to 2\ell 2\nu$ (see Fig. 14 (b)). Both these contributions are governed by the effective $hh^{+}h^{-}$ coupling given by

[TABLE]

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

[TABLE]

where $\mu^{i}$ (with $i=$ ggF, VBF, $Vh$, and $t\bar{t}h$) and $\mu_{f}$ (with $f=ZZ^{\star},WW^{\star},\gamma\gamma,\tau^{+}\tau^{-},b\bar{b}$) 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 $\mu^{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. Babu:2018uik ; Jana:2017hqg , 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 $\sqrt{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: $h\to\gamma\gamma$ CMS:1900lgv ; CMS:2018rbc ; ATLAS:2018uso ; ATLAS:2019aqa , $h\to ZZ^{\star}$ CMS:2019chr ; ATLAS:2018bsg , $h\to WW^{\star}$ Aaboud:2018jqu ; Aad:2019lpq ; Sirunyan:2018egh , $h\to\tau\tau$ CMS:2019pyn ; Aaboud:2018pen and $h\to b\bar{b}$ Sirunyan:2018kst ; Aaboud:2018zhk ; Aaboud:2018gay .

We formulate the modified $h\to\gamma\gamma$ decay rate as

[TABLE]

where the scaling factor $\kappa_{\gamma}$ is given by

[TABLE]

where $N^{f}_{c}=3~{}(1)$ is the color factor for quark (lepton), $\sum_{f}$ is the sum over the SM fermions $f$ with charge $Q_{f}$, and the loop functions are given by Djouadi:2005gi

[TABLE]

The parameters $\tau_{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.81) due to the charged scalar can interfere with the SM part either constructively or destructively, depending on the sign of the effective coupling $\lambda_{hh^{+}h^{-}}$ in Eq. (4.78).

As for the new three-body decay mode $h\rightarrow h^{\pm}h^{\mp^{\star}}\to h^{\pm}\ell\nu$, the partial decay rate is given by

[TABLE]

where $Y$ is the Yukawa coupling defined in Eq. (4.34), $\Gamma_{h^{+}}={\rm Tr}(Y^{\dagger}Y)m_{h^{+}}/8\pi$ is the total decay width of $h^{+}$, and $r=m_{h^{+}}^{2}/m_{h}^{2}$. With this new decay mode, the signal strength in the $h\to 2\ell 2\nu$ channel will be modified to include $\Gamma(h\to h^{\pm}\ell\nu\to 2\ell 2\nu)$ along with the SM contribution from $\Gamma(h\to WW^{\star}\to 2\ell 2\nu)$, and to some extent, from $\Gamma(h\to ZZ^{\star}\to 2\ell 2\nu)$.

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 $\lambda_{hh^{\pm}h^{\mp}}$ which is a function of the cubic coupling $\mu$, quartic couplings $\lambda_{3}$ and $\lambda_{8}$, and the mixing angle $\sin{\varphi}$ (cf. Eq. (4.79)). For suppressed effective coupling $\lambda_{hh^{\pm}h^{\mp}}$ 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{\varphi}$, which implies large value of $\mu$ (cf. Eq. (4.28)). In order to find the maximum allowed value of $\sin\varphi$, we take $\lambda_{3}=\lambda_{8}$ in Eq. (4.79) and show in Fig. 15 the Higgs signal strength constraints in the $\lambda_{8}-\sin\varphi$ plane. The red, blue, yellow, cyan, and green-shaded regions are excluded by the signal strength limits $\gamma\gamma,WW^{\star},ZZ^{\star},\tau\tau$, and $b\bar{b}$ 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 ${\cal O}(1)$ Yukawa couplings (cf. Sec. 4.6). From Fig. 15, we see that the $h\to\gamma\gamma$ signal strength gives the most stringent constraint. If we allow $\lambda_{8}$ to be as large as 3, then we can get maximum value of $\sin\varphi$ 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 Aaboud:2018puo and CMS Sirunyan:2019twz collaborations have put 95% CL upper limits on the Higgs boson total width $\Gamma_{h}$ from measurement of off-shell production in the $ZZ\to 4\ell$ channel. Given the SM expectation $\Gamma_{h}^{\rm SM}\sim 4.1$ MeV, we use the CMS upper limit on $\Gamma_{h}<9.16$ MeV Sirunyan:2019twz to demand that the new contribution (mostly from $h\to h^{\pm}h^{\mp\star}$, because the $h\to\gamma\gamma$ 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.

4.11 Monophoton constraint from LEP

Large neutrino NSI with electrons inevitably leads to a new contribution to the monophoton process $e^{+}e^{-}\to\nu\bar{\nu}\gamma$ that can be constrained using LEP data Berezhiani:2001rs . 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^{-}\to\nu_{\alpha}\bar{\nu}_{\beta}\gamma$ can be expressed as $\sigma=\sigma^{\rm SM}+\sigma^{\rm NS}$, where $\sigma^{\rm SM}$ is the SM cross section (for $\alpha=\beta$) and $\sigma^{\rm 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 $\alpha=\beta=e$). Moreover, for either $\alpha$ or $\beta$ not equal to $e$, the $W$ contribution is absent. For $\alpha\neq\beta$, the $Z$ contribution is also absent.

The monophoton process has been investigated carefully by all four LEP experiments Tanabashi:2018oca , but the most stringent limits on the cross section come from the L3 experiment, both on Acciarri:1998vf and off Achard:2003tx $Z$-pole. We use these results to derive constraints on the charged scalar mass and Yukawa coupling. The constraint $|\sigma-\sigma^{\rm exp}|\leq\delta\sigma^{\rm exp}$, where $\sigma^{\rm exp}\pm\delta\sigma^{\rm exp}$ is the experimental result, can be expressed in the following form:

[TABLE]

We evaluate the ratio $\sigma^{\rm exp}/\sigma^{\rm SM}$ by combining the L3 results Acciarri:1998vf ; Achard:2003tx with an accurate computation of the SM cross section, both at $Z$-pole and off $Z$-pole. Similarly, we compute the ratio $\sigma^{\rm NS}/\sigma^{\rm SM}$ numerically as a function of the charged scalar mass $m_{h^{+}}$ and the Yukawa coupling $Y_{\alpha\beta}\sin\varphi$. For comparison of cross sections at $Z$-pole, we adopt the same event acceptance criteria as in Ref. Acciarri:1998vf , i.e., we allow photon energy within the range 1 GeV $<E_{\gamma}<$ 10 GeV and the angular acceptance 45$\degree<\theta_{\gamma}<135\degree$. Similarly, for the off $Z$-pole analysis, we adopt the same event topology as described in Ref. Achard:2003tx : i.e., 14$\degree<\theta_{\gamma}<166\degree$, 1 GeV $<E_{\gamma}$, and $p_{T}^{\gamma}>0.02\sqrt{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. Fox:2011fx using an effective four-fermion approximation.

4.12 NSI predictions

The new singly-charged scalars $\eta^{+}$ 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.14) and (4.34)). 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 $\varphi$ (cf. Eq. (4.28)). As mentioned in Sec. 4.2, satisfying the neutrino mass requires the product $f\cdot Y$ to be small. For $Y\sim{\cal O}(1)$, we must have $f\sim 10^{-8}$ to get $m_{\nu}\sim 0.1$ eV (cf. Eq. (4.36)). In this case, the NSI from Figs. 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 $\eta^{+}$ and $H_{2}^{+}$ are replaced by $h^{+}$ and $H^{+}$ respectively (cf. Eq. (4.27)).

The effective NSI Lagrangian for the contribution from Fig. 17 (b) is given by

[TABLE]

where in the second step, we have used the Fierz transformation. Comparing Eq. (4.88) with Eq. (3.9), we obtain the $h^{+}$-induced matter NSI parameters (setting $\rho=\sigma=e$)

[TABLE]

Thus, the diagonal NSI parameters $\varepsilon_{\alpha\alpha}$ depend on the Yukawa couplings $|Y_{\alpha e}|^{2}$, and are always positive in this model, whereas the off-diagonal ones $\varepsilon_{\alpha\beta}$ (with $\alpha\neq\beta$) involve the product $Y_{\alpha e}Y_{\beta e}^{\star}$ and can be of either sign, or even complex. Also, we have a correlation between the diagonal and off-diagonal NSI:

[TABLE]

which is a distinguishing feature of the model.

Similarly, Fig. 17 (d) gives the $H^{+}$-induced matter NSI contribution:

[TABLE]

Hence, the total matter NSI induced by the charged scalars in the Zee model can be expressed as

$\displaystyle\varepsilon_{\alpha\beta}\ \equiv\ \varepsilon_{\alpha\beta}^{(h^{+})}+\varepsilon_{\alpha\beta}^{(H^{+})}\ =\ \frac{1}{4\sqrt{2}G_{F}}Y_{\alpha e}Y_{\beta e}^{\star}\left(\frac{\sin^{2}{\varphi}}{m_{h^{+}}^{2}}+\frac{\cos^{2}{\varphi}}{m_{H^{+}}^{2}}\right)\,.$

(4.92)

To get an idea of the size of NSI induced by Eq. (4.92), let us take the diagonal NSI parameters from the light charged scalar contribution in Eq. (4.89):

[TABLE]

Thus, for a given value of $m_{h^{+}}$, the NSI are maximized for maximum allowed values of $|Y_{\alpha e}|$ and $\sin{\varphi}$. Following Eq. (4.78), we set the trilinear coupling $\lambda_{hh^{+}h^{-}}\to 0$, thus minimizing the constraints from Higgs signal strength. We also assume $\lambda_{3}=\lambda_{8}$ to get

[TABLE]

Now substituting this into Eq. (4.28), we obtain

[TABLE]

Furthermore, assuming the heavy charged and neutral scalars to be mass-degenerate, the LEP contact interaction constraints (cf. Sec. 4.6) require

[TABLE]

where $\Lambda_{\alpha}=10$ TeV, 7.9 TeV and 2.2 TeV for $\alpha=e,\mu,\tau$, respectively LEP:2003aa . Combining Eqs. (4.93), (4.95) and (4.96), we obtain

[TABLE]

Using benchmark values of $m_{h^{+}}=100\,\text{GeV}$ and $\lambda_{8}=3$, we obtain:

[TABLE]

Although a rough estimate, this tells us that observable NSI can be obtained in the Zee model, especially in the $\tau$ sector. To get a more accurate prediction of the NSI in the Zee model and to reconcile large NSI with all relevant theoretical and experimental constraints, we use Eq. (4.92) to numerically calculate the NSI predictions, as discussed below.

4.12.1 Heavy neutral scalar case

First, we consider the case with heavy neutral and charged scalars, so that the LEP contact interaction constraints (cf. Sec. 4.6) are valid. To be concrete, we have fixed the heavy charged scalar mass $m_{H^{+}}=700$ GeV and the quartic couplings $\lambda_{3}=\lambda_{8}=3$. In this case, the heavy charged scalar contribution to NSI in Eq. (4.92) can be ignored. The NSI predictions in the light charged scalar mass versus Yukawa coupling plane are shown by black dotted contours in Fig. 18 for diagonal NSI and Fig. 19 for off-diagonal NSI. The theoretical constraints on $\sin\varphi$ from charge-breaking minima (cf. Sec. 4.3) and $T$-parameter (cf. Sec. 4.4) constraints are shown by the light and dark green-shaded regions, respectively. Similarly, the Higgs precision data constraint (cf. Sec. 4.10) on $\sin\varphi$ is shown by the brown-shaded region. To cast these constraints into limits on $Y_{\alpha e}\sin\varphi$, we have used the LEP contact interaction limits on $Y_{\alpha e}$ (cf. Sec. 4.6) for diagonal NSI, and similarly, the cLFV constraints (cf. Sec. 4.5) for off-diagonal NSI, and combined these with the CBM, $T$-parameter and Higgs constraints, which are all independent of the light charged scalar mass. Also shown in Figs. 18 and 19 are the LEP and/or LHC constraints on light charged scalar (cf. Sec. 4.7) combined with the lepton universality constraints from $W$ and $\tau$ decays (cf. Secs. 4.8 and 4.9), which exclude the blue-shaded region below $m_{h^{+}}\sim 100$ GeV. In addition, the LEP monophoton constraints from off $Z$-pole search (cf. Sec. 4.11) are shown in Fig. 18 by the light purple-shaded region. The corresponding limit from LEP on $Z$-pole search (shown by the purple dashed line in Fig. 18 (c) turns out to be weaker.

The model predictions for NSI are then compared with the current direct experimental constraints from neutrino-electron scattering experiments (red/yellow-shaded), and the global-fit constraints (orange-shaded) Esteban:2018ppq which include the neutrino oscillation data Tanabashi:2018oca , as well as the recent results from COHERENT experiment Akimov:2017ade ;121212For related NSI studies using the COHERENT data, see e.g. Refs. Coloma:2017egw ; Coloma:2017ncl ; Liao:2017uzy ; AristizabalSierra:2018eqm ; Altmannshofer:2018xyo . see Table IX for more details. For neutrino-electron scattering constraints, we only considered the constraints on $\varepsilon_{\alpha\beta}^{eR}$ Davidson:2003ha ; Barranco:2007ej ; Deniz:2010mp ; Agarwalla:2019smc , since the dominant NSI in the Zee model always involves right-handed electrons (cf. Eq. (4.88)). For $\varepsilon_{\mu\mu}$, we have rederived the CHARM II limit following Ref. Davidson:2003ha , but using the latest PDG value for $s_{w}^{2}=0.22343$ (on-shell) Tanabashi:2018oca . Specifically, we used the CHARM II measurement of the $Z$-coupling to right-handed electrons $g_{R}^{e}=0.234\pm 0.017$ obtained from their $\nu_{\mu}e\to\nu e$ data Vilain:1994qy and compared with the SM value of $(g_{R}^{e})_{\rm SM}=s_{w}^{2}$ to obtain a $90\%$ CL limit on $\varepsilon_{\mu\mu}<0.038$, which is slightly weaker than the limit of 0.03 quoted in Ref. Barranco:2007ej . Nevertheless, the CHARM limit turns out to be the strongest in realizing maximum $\varepsilon_{\mu\mu}$ in the Zee model, as shown in Fig. 18 (b).

There is a stronger constraint on $|\varepsilon_{\tau\tau}-\varepsilon_{\mu\mu}|<9.3\%$ from the IceCube atmospheric neutrino oscillation data Fornengo:2001pm ; Esmaili:2013fva ; Day:2016shw . In general, this bound can be evaded even for large NSI, if e.g. both $\varepsilon_{\mu\mu}$ and $\varepsilon_{\tau\tau}$ are large and there is a cancellation between them. However, in the Zee model, such cancellation cannot be realized, because we can only allow for one large diagonal NSI at a time, otherwise there will be stringent constraints from cLFV (cf. Sec. 4.5). For instance, making both $\varepsilon_{\mu\mu}$ and $\varepsilon_{\tau\tau}$ large necessarily implies a large $\varepsilon_{\mu\tau}$ (due to the relation given by Eq. (4.90)), which is severely constrained by $\tau^{-}\to\mu^{-}e^{-}e^{+}$ (cf. Table VII and Fig. 19 (a)) and also by IceCube itself Esmaili:2013fva ; Salvado:2016uqu ; Aartsen:2017xtt . Therefore, the bound on $\varepsilon_{\tau\tau}-\varepsilon_{\mu\mu}$ is equally applicable to both $\varepsilon_{\mu\mu}$ and $\varepsilon_{\tau\tau}$. This is shown by the brown-shaded regions in Fig. 18 (b) and (c), respectively. This turns out to be the most stringent constraint for $\varepsilon_{\tau\tau}$, although the model allows for much larger NSI, as shown by the black dotted contours in Fig. 18 (c).

For completeness, we also include in Fig. 18 global-fit constraints from neutrino oscillation plus scattering experiments Esteban:2018ppq .131313We use the constraints on $\varepsilon_{\alpha\beta}^{p}$ from Ref. Esteban:2018ppq , assuming that these will be similar for $\varepsilon_{\alpha\beta}^{e}$ due to charge-neutrality in matter. The global-fit analysis assumes the simultaneous presence of all $\varepsilon_{\alpha\beta}$’s, and therefore, the corresponding limits on each $\varepsilon_{\alpha\beta}$ are much weaker than the ones derived from oscillation or scattering data alone, due to parameter degeneracies. For instance, the global-fit constraint on $\varepsilon_{\tau\tau}\in[-35\%,140\%]$ (cf. Table IX) is significantly affected by the presence of nonzero $\varepsilon_{ee}$ and $\varepsilon_{e\tau}$ Friedland:2004ah , which were set to zero in the IceCube analysis of Ref. Esmaili:2013fva .

Also shown in Fig. 18 (blue solid lines) are the future sensitivity at long-baseline neutrino oscillation experiments, such as DUNE with 300 kt.MW.yr and 850 kt.MW.yr of exposure, derived at 90% CL using GloBES3.0 Huber:2007ji with the DUNE CDR simulation configurations Alion:2016uaj . Here we have used $\delta~{}({\rm true})=-\pi/2$ for the true value of the Dirac CP phase and marginalized over all other oscillation parameters dev_pondd . We find that even the most futuristic DUNE sensitivity will not be able to surpass the current constraints on the Zee model. On the other hand, the current neutrino scattering experiments like COHERENT and atmospheric neutrino experiments such as IceCube should be able to probe a portion of the allowed parameter space for $\varepsilon_{\mu\mu}$ and $\varepsilon_{\tau\tau}$, respectively.

4.12.2 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.92) 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 $\lambda_{hh^{+}h^{-}}\to 0$ in Eq. (4.78). 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).

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 can get maximum $\varepsilon_{ee}$ of 8%, $\varepsilon_{\mu\mu}$ of 3.8%, and $\varepsilon_{\tau\tau}$ of 9.3%, only limited by direct experimental searches (TEXONO, CHARM and IceCube, respectively). Thus, the future neutrino experiments could probe diagonal NSI in the Zee model. As for the off-diagonal NSI, they require the presence of at least two non-zero Yukawa couplings $Y_{\alpha 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 by any neutrino scattering or oscillation experiment in the foreseeable future.

4.13 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 mass matrix in the Zee model is given by Eq. (4.36) which is diagonalized by the unitary transformation

[TABLE]

where $\widehat{M}_{\nu}={\rm diag}(m_{1},m_{2},m_{3})$ is the diagonal mass matrix with the eigenvalues $m_{1,2,3}$ and $U_{\rm PMNS}$ is the 3 $\times$ 3 lepton mixing matrix. In the standard parametrization Tanabashi:2018oca ,

[TABLE]

where $c_{ij}\equiv\cos{\theta_{ij}}$, $s_{ij}\equiv\sin{\theta_{ij}}$, $\theta_{ij}$ being the mixing angle between different flavor eigenstates $i$ and $j$, and $\delta$ is the Dirac $\mathcal{CP}$ phase. We diagonalize the neutrino mass matrix (4.36) numerically, assuming certain forms of the Yukawa coupling matrices given below. The unitary matrix thus obtained is converted to the mixing angles $\theta_{ij}$ using the following relations from Eq. (4.100):

[TABLE]

Since the NSI expressions in Eq. (4.92) depend on $Y_{\alpha 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.

[TABLE]

For BP I, substituting $Y$ from Eq. (4.105) in Eq. (4.36), we get a symmetric neutrino mass matrix as follows:

[TABLE]

where $a_{0}=\kappa f_{\mu\tau}Y_{ee}$ fixes the overall scale, and the entries in $M_{\nu}$ are given by

[TABLE]

and we have defined the ratios $x_{1}=\frac{f_{e\mu}}{f_{\mu\tau}}$ , $x_{2}=\frac{f_{e\tau}}{f_{\mu\tau}}$, $y_{13}=\frac{Y_{e\tau}}{Y_{ee}}$, $y_{22}=\frac{Y_{\mu\mu}}{Y_{ee}}$, $y_{23}=\frac{Y_{\mu\tau}}{Y_{ee}}$, $y_{32}=\frac{Y_{\tau\mu}}{Y_{ee}}$, and $y_{33}=\frac{Y_{\tau\tau}}{Y_{ee}}$. Similarly, for BPs II and III, one can absorb $Y_{\mu\mu}$ and $Y_{\tau\tau}$ 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$\sigma$ allowed range from a recent NuFit4 global analysis  Esteban:2018azc . 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\sigma$, $2\sigma$ and $3\sigma$ confidence regions of the global-fit results Esteban:2018azc (without inclusion of the Super-K atmospheric $\Delta\chi^{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 model predictions for each benchmark point, are shown for comparison. It is clear that the theoretical predictions are within the observed 3$\sigma$ range in each case.

5 NSI in one-loop leptoquark model

There are only four kinds of scalar LQs that can interact with the neutrinos at the renormalizable level in the SM (see Table I): $Ld^{c}\Omega$, $LQ\chi^{\star}$, $LQ\bar{\rho}$ and $Lu^{c}\delta$.141414The LQ fields $\Omega,\,\chi^{\star},\,\bar{\rho},\,\delta$ are often denoted as $S_{1},\,S_{3},\,R_{2},\,\widetilde{R}_{2}$ respectively Buchmuller:1986zs . 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 $B$-decay anomalies, viz., $R_{D^{(\star)}}$ and $R_{K^{(\star)}}$ (see e.g., Aebischer:2019mlg ). But it turns out that none of these scalar LQ models can simultaneously explain both $R_{D^{(\star)}}$ and $R_{K^{(\star)}}$ Angelescu:2018tyl .

We start with a LQ variant of the Zee model that generates small neutrino masses at one-loop level, via the operator is ${\cal O}_{3b}$ (cf. Eq. (2.7c)). It turns out that ${\cal O}_{3b}$ will induce neutrino masses at one-loop, while $\mathcal{O}_{3a}$, owing to the $SU(2)_{L}$ index structure, will induce $m_{\nu}$ at the two-loop level. A UV complete model of $\mathcal{O}_{3a}$ will be presented in Sec. 7.2.3. More precisely, the model of this section corresponds to $\mathcal{O}_{3}^{8}$ of Table III, which involves two LQ fields and no new fermions. All other realizations of $\mathcal{O}_{3}$ will be analyzed in subsequent sections.

The phenomenology of the basic LQ model generating $\mathcal{O}_{3}^{8}$ 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 $\mathcal{O}_{3}$ models, as well as other one-loop, two-loop and three-loop LQ models discussed in subsequent sections.

To realize operator $\mathcal{O}_{3b}$ the $SU(2)_{L}$ doublet and singlet scalars of the Zee model Zee:1980ai 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 $\widetilde{Q}$ and $\widetilde{d^{c}}$ fields of the MSSM Hall:1983id ; Barbier:2004ez ; Dreiner:1997uz . For a non-supersymmetric description and analysis of the model, see Ref. AristizabalSierra:2007nf .

The gauge symmetry of the model denoted as $\mathcal{O}_{3}^{8}$ is the same as the SM: $SU(3)_{c}\times SU(2)_{L}\times U(1)_{Y}$. In addition to the SM Higgs doublet $H\left({\bf 1},{\bf 2},\frac{1}{2}\right)$, two $SU(3)_{c}$ triplet LQ fields $\Omega\left({\bf 3},{\bf 2},\frac{1}{6}\right)=\left(\omega^{2/3},\omega^{-1/3}\right)$ and $\chi^{-1/3}\left({\bf 3},{\bf 1},-\frac{1}{3}\right)$ are introduced. The Yukawa Lagrangian relevant for neutrino mass generation in the model is given by

[TABLE]

Here $\{\alpha,\beta\}$ are family indices and {$i,j$} are $SU(2)_{L}$ indices as before. As in the Zee model, a cubic scalar coupling is permitted, given by

[TABLE]

which ensures lepton number violation.

Once the neutral component of the SM Higgs doublet acquires a VEV, the cubic term in the scalar potential (5.119) will generate mixing between the $\omega^{-1/3}$ and $\chi^{-1/3}$ fields, with the mass matrix given by:

[TABLE]

where $m_{\omega}^{2}$ and $m_{\chi}^{2}$ include the bare mass terms plus a piece of the type $\lambda v^{2}$ arising from the SM Higgs VEV. The physical states are denoted as $\{X_{1}^{-1/3},X_{2}^{-1/3}\}$, defined as

[TABLE]

with the mixing angle given by

[TABLE]

The squared mass eigenvalues of these states are:

[TABLE]

Neutrino masses are induced via the one-loop diagram shown in Fig. 22. The mass matrix is given by:

[TABLE]

Here $M_{d}$ is the diagonal down-type quark mass matrix. Acceptable neutrino masses and mixing can arise in the model for a variety of parameters. Note that the induced $M_{\nu}$ 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 $\omega^{-1/3}$ or the $\chi^{-1/3}$ field, without relying on their mixing, we shall adopt a scenario where the couplings $\lambda_{\alpha\beta}$ are of order one, while $\lambda^{\prime}_{\alpha\beta}\ll 1$. Such a choice would realize small neutrino masses. One could also consider $\lambda^{\prime}\sim{\cal O}(1)$, with $\lambda\ll 1$ as well. However, in the former case, there is a GIM-like suppression in the decay rate for $\ell_{\alpha}\rightarrow\ell_{\beta}+\gamma$ Babu:2010vp , which makes the model with $\lambda\sim{\cal O}(1),\lambda^{\prime}\ll 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.

5.1 Low-energy constraints

One interesting feature of the LQ model presented in this section is that the radiative decay $\ell_{\alpha}\rightarrow\ell_{\beta}+\gamma$ is suppressed in the model due to a GIM-like cancellation. On the other hand, $\mu-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 $|\lambda\lambda^{\prime}|\ll 1$ in order to generate the correct magnitude of the neutrino masses (cf. Eq. (5.124)), we shall primarily consider the case where $|\lambda^{\prime}|\ll 1$ with $|\lambda|$ being of order one. This is the case where the constraints from radiative decays are nonexistent. If on the other hand, $|\lambda|\ll 1$ and $|\lambda^{\prime}|$ 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: $\ell_{\alpha}\to\ell_{\beta}+\gamma$, $\mu-e$ conversion in nuclei, $\ell_{\alpha}\to\bar{\ell}_{\beta}\ell_{\gamma}\ell_{\delta}$ (with at least two of the final state leptons being of same flavor), $\tau\to\ell\pi$, $\tau\to\ell\eta$, $\tau\to\ell\eta^{\prime}$ (where $\ell=e$ or $\mu)$, and APV.

5.1.1 Atomic parity violation

The strongest constraints on the $\lambda_{ed}$ and $\lambda^{\prime}_{ed}$ couplings come from atomic parity violation (APV) Roberts:2014bka , analogous to the $R$-parity violating supersymmetric case Kao:2009fg . The diagrams shown in Fig. 23 lead to the following effective couplings between up/down quarks and electrons:

[TABLE]

where we have used the Fierz transformation in the second step. The parity-violating parts of these interactions are given by

[TABLE]

On the other hand, the parity-violating SM interactions at tree-level are given by

[TABLE]

with

[TABLE]

Correspondingly, the weak charge of an atomic nucleus with $Z$ protons and $N$ neutrons is given by

[TABLE]

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.126) will shift the weak charge to

[TABLE]

There are precise experiments measuring APV in cesium, thallium, lead and bismuth Safronova:2017xyt . The most precise measurement comes from cesium (at the 0.4% level Porsev:2009pr ), so we will use this to derive constraints on LQ. For ${}_{55}^{133}$Cs, Eq. (5.130) becomes

[TABLE]

Taking into account the recent atomic structure calculation Roberts:2014bka , the experimental value of the weak charge of ${}_{55}^{133}$Cs is given by Tanabashi:2018oca

[TABLE]

whereas the SM prediction is  Roberts:2014bka ; Tanabashi:2018oca

[TABLE]

based on a global-fit to all electroweak observables with radiative corrections. Assuming new radiative corrections from LQ are small and saturating the difference between Eqs. (5.132) and (5.133), we obtain a $2\sigma$ allowed range of $\delta Q_{w}$:

[TABLE]

Comparing this with Eq. (5.131), we obtain the corresponding $2\sigma$ bounds on $\lambda_{ed}$ and $\lambda^{\prime}_{ed}$ as a function of the LQ mass as follows:

[TABLE]

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.132)) is $1.5\sigma$ larger than the SM prediction (cf. Eq. (5.133)), while the doublet LQ contribution to $Q_{w}$ goes in the opposite direction (cf. Eq. (5.131)).

5.1.2 $\mu-e$ conversion

Another constraint on the LQ model being discussed comes from the cLFV process of coherent $\mu-e$ conversion in nuclei ($\mu N\to 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. Kuno:1999jp , we can write down the branching ratio for this process as Babu:2010vp

[TABLE]

where $\vec{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_{\rm eff}$ is the effective atomic number, $F_{p}$ is the nuclear matrix element, and $\Gamma_{N}$ is the muon capture rate of the nucleus. Here we take $|\vec{p}_{e}|\simeq E_{e}\simeq m_{\mu}$ and use the values of $Z_{\rm eff}$ and $F_{p}$ from Ref. Kitano:2002mt , and the value of $\Gamma_{N}$ from Ref. Suzuki:1987jf . Comparing the model predictions from Eq. (5.136) with the experimental limits for different nuclei Kaulard:1998rb ; Bertl:2006up ; Honecker:1996zf , we obtain the constraints on the Yukawa couplings (either $\lambda$ or $\lambda^{\prime}$) and LQ mass as shown in Table XII.

5.1.3 $\ell_{\alpha}\to\bar{\ell}_{\beta}\ell_{\gamma}\ell_{\delta}$ decay

LQs do not induce trilepton decays of the type $\mu\to 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. Babu:2010vp . With the Yukawa couplings $\lambda$ being of order one, but with $|\lambda^{\prime}|\ll 1$, the branching ratio for $\mu^{-}\to e^{+}e^{-}e^{-}$ decay is given by Babu:2010vp

[TABLE]

where

[TABLE]

Here we have kept only those couplings that are relevant for neutrino NSI, and we have assumed that there are no accidental cancellations among various contributions. Using BR($\mu\rightarrow 3e)<1.0\times 10^{-12}$ Bertl:1985mw , we obtain

[TABLE]

Analogous constraints from $\tau\to 3e$ and $\tau\to 3\mu$ are less stringent. For example, from BR($\tau\to 3e)<1.4\times 10^{-8}$ Amhis:2016xyh , and using Eq. (5.137) with a multiplicative factor of BR($\tau\to\bar{\nu}_{\ell}\ell\nu_{\tau})=0.174$, we obtain

[TABLE]

Similarly, from BR($\tau\to 3\mu)<1.2\times 10^{-8}$ Amhis:2016xyh we obtain

[TABLE]

The constraint on $|\lambda_{ed}\lambda_{\mu d}^{\star}|$ from the trilepton decay (cf. Eq. (5.139)) turns out to be weaker than those from $\mu-e$ conversion (cf. Table XII). Similarly, the constraints on $|\lambda_{ed}\lambda_{\tau d}^{\star}|$ and $|\lambda_{\mu d}\lambda_{\tau d}^{\star}|$ from the trilepton decay (cf. Eqs. (5.140) and (5.141)) turn out to be weaker than those from semileptonic tau decays (cf. Table XIV).

5.1.4 $\ell_{\alpha}\to\ell_{\beta}\gamma$ constraint

The lepton flavor violating radiative decay $\ell_{\alpha}\to\ell_{\beta}+\gamma$ 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 $\lambda$ leads to suppressed decay rates for $\ell_{\alpha}\to\ell_{\beta}+\gamma$, owing to a GIM-like cancellation. The coupling of the $\omega^{2/3}$ LQ has the form $\ell_{\alpha L}\overline{d_{\beta R}^{c}}\omega^{2/3}$, which implies that $Q_{B}=2/3$ and $Q_{F}=-1/3$ in Eq. (4.45). Consequently, the rate becomes proportional to a factor which is at most of order $(m_{b}^{2}/m_{\omega}^{2})^{2}$. Thus, the off-diagonal couplings of $\lambda$ are unconstrained by these decays.

On the other hand, the $\chi^{-1/3}$ LQ field does mediate $\ell_{\alpha}\to\ell_{\beta}+\gamma$ decays, proportional to the Yukawa coupling matrix $\lambda^{\prime}$. The relevant couplings have the form $\bar{u}_{L}\ell_{L}\chi^{\star}$, which implies that $Q_{F}=-2/3$ and $Q_{B}=1/3$ in Eq. (4.45). We find the decay rate to be

[TABLE]

where $9=3^{2}$ is a color factor. Here we have assumed $t=m_{F}^{2}/m_{B}^{2}\to 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.

5.1.5 Semileptonic tau decays

The decays $\tau^{-}\to\ell^{-}\pi^{0},\,\ell^{-}\eta,\,\ell^{-}\eta^{\prime}$, with $\ell=e$ or $\mu$ will occur at tree level mediated by the doublet LQ $\omega^{2/3}$ or the singlet LQ $\chi^{-1/3}$. The relevant Feynman diagrams are shown in Fig. 26. The decay rate for $\tau^{-}\to\ell^{-}\pi^{0}$ mediated by $\omega$ LQ is given by

[TABLE]

where

[TABLE]

If this decay is mediated by the $\chi$ LQ, the same relation will hold, up to a factor of $|V_{ud}|^{2}$, with the replacement $(\lambda,\,m_{\omega})\rightarrow(\lambda^{\prime},\,m_{\chi})$. The rates for $\tau^{-}\to\ell^{-}\eta$ and $\tau^{-}\to\ell^{-}\eta^{\prime}$ can be obtained from Eq. (5.143) by the replacement $(f_{\pi},\,m_{\pi})\rightarrow(m_{\eta},\,f_{\eta}^{q})$ and $(m_{\eta^{\prime}},\,f_{\eta^{\prime}}^{q})$ respectively. Here we have defined the matrix elements to be

[TABLE]

The sign difference in Eq. (5.145) is due to the fact that the state $|\pi^{0}\rangle=(u\bar{u}-d\bar{d})/\sqrt{2}$. As for $|\eta\rangle$ and $|\eta^{\prime}\rangle$ states, these are obtained from the mixing of the flavor states $|\eta_{q}\rangle=(\bar{u}u+\bar{d}d)/\sqrt{2}$ and $|\eta_{s}\rangle=\bar{s}s$:

[TABLE]

The matrix elements entering semileptonic $\tau$ decays are then related as

[TABLE]

where $f_{q}$ is defined through

[TABLE]

The mixing angle $\phi$ and the decay parameter $f_{q}$ have been determined to be Beneke:2002jn

[TABLE]

Using these relations, and with $f_{\pi}\simeq 130$ MeV, we have $f_{\eta}^{q}\simeq 108$ MeV and $f_{\eta^{\prime}}^{q}\simeq 89$ MeV Li:2005rr . Using these values and the experimental limits on the semileptonic branching ratios Tanabashi:2018oca , we obtain limits on products of Yukawa couplings as functions of the LQ mass, which are listed in Table XIV. It turns out that these limits are the most constraining for off-diagonal NSI mediated by LQs.

We should mention here that similar diagrams as in Fig. 26 will also induce alternative pion and $\eta$-meson decays: $\pi^{0}\to e^{+}e^{-}$ and $\eta\to\ell^{+}\ell^{-}$ (with $\ell=e$ or $\mu$). In the SM, ${\rm BR}(\pi^{0}\to e^{+}e^{-})=6.46\times 10^{-8}$ Tanabashi:2018oca , compared to ${\rm BR}(\pi^{0}\to\gamma\gamma)\simeq 0.99$. Specifically, the absorptive part of $\pi^{0}\to e^{+}e^{-}$ decay rate151515The dispersive part of $\pi^{0}\to e^{+}e^{-}$ decay rate is found to be 32% smaller than the absorptive part in the vector meson dominance Babu:1982yz . is given by Berman:1960zz ; Babu:1982yz

[TABLE]

where $\beta=\sqrt{1-4m_{e}^{2}/m_{\pi}^{2}}$. For LQ mediation, the suppression factor $(m_{e}/m_{\pi})^{2}\sim 1.4\times 10^{-5}$ is replaced by the factor $(m_{\pi}/m_{\omega})^{4}\sim 3.3\times 10^{-16}$ for a TeV-scale LQ. Similar suppression occurs for the $\eta$ decay processes $\eta\to\ell^{+}\ell^{-}$ (with $\ell=e$ or $\mu$) Geffen:1965zz ; Babu:1982yz . Therefore, both pion and $\eta$ decay constraints turn out to be much weaker than those from $\tau$ decay given in Table XIV.

5.1.6 Rare $D$-meson decays

The coupling matrix $\lambda^{\prime}$ of Eq. (5.118) 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 $\chi$ LQ interactions with the physical quarks, in this basis, read as

[TABLE]

Here $V$ is the CKM mixing matrix. In particular, the Lagrangian contains the following terms:

[TABLE]

The presence of these terms will result in the rare decays $D^{0}\rightarrow\ell^{+}\ell^{-}$ as well as $D\rightarrow\pi\ell^{+}\ell^{-}$ where $\ell=e,\,\mu$. The partial width for the decay $D^{0}\rightarrow\ell^{+}\ell^{-}$ is given by

[TABLE]

Here we have used the effective Lagrangian arising from integrating out the $\chi$ field to be

[TABLE]

and the hadronic matrix element

[TABLE]

Using $f_{D}=200$ MeV, we list the constraint arising from this decay in Table XV. It will turn out that the NSI parameter $\varepsilon_{\mu\mu}$ will be most constrained by the limit $D^{0}\to\mu^{+}\mu^{-}$, in cases where $\chi$ LQ 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 $\Omega$. 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^{+}\rightarrow\pi^{+}\ell^{+}\ell^{-}$ is mediated by the same effective Lagrangian as in Eq. (5.156). The hadronic matrix element is now given by

[TABLE]

with $q^{2}=(p_{1}-p_{2})^{2}$. Since the $F_{-}(q^{2})$ term is proportional to the final state lepton mass, it can be ignored. For the form factor $F_{+}(q^{2})$ we use

[TABLE]

For the $D^{\star}\to D\pi$ decay constant we use $g_{D^{\star}D\pi}=0.59$ Burdman:2003rs . Vector meson dominance hypothesis gives very similar results Babu:1987xe . With these matrix elements, the decay rate is given by

[TABLE]

The function ${\cal F}$ is defined as

[TABLE]

Note that in the limit of infinite $D^{\star}$ mass, this function ${\cal F}$ reduces to $m_{D}^{6}/24$. The numerical value of the function is ${\cal F}\simeq 2.98$ GeV6. Using $f_{D}=200$ MeV, $f_{\pi}=130$ MeV, $g_{D^{\star}D\pi}=0.59$ and the experimental upper limits on the corresponding branching ratios  Tanabashi:2018oca , we obtain bounds on the $\lambda^{\prime}$ 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}\rightarrow\pi^{0}\ell^{+}\ell^{-}$ are somewhat weaker than the $D^{+}$ decay limits and are automatically satisfied when the $D^{+}$ semileptonic rates are satisfied.

5.2 Contact interaction constraints

High-precision measurements of inclusive $e^{\pm}p\to e^{\pm}p$ scattering cross sections at HERA with maximum $\sqrt{s}=320$ GeV Abramowicz:2019uti and $e^{+}e^{-}\to q\bar{q}$ scattering cross sections at LEP II with maximum $\sqrt{s}=209$ GeV LEP:2003aa can be used in an effective four-fermion interaction theory to set limits on the new physics scale $\Lambda>\sqrt{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.125) with Eq. (4.52) (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 $\eta^{d}_{LR}=\eta^{u}_{LL}=-1$. The corresponding experimental bounds on $\Lambda^{-}$ 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\to\ell^{+}\ell^{-}$. However, for a TeV-scale LQ as in our case, the LHC contact interaction bounds Aaboud:2017buh ; Sirunyan:2018ipj with $\sqrt{s}=13$ TeV are not applicable. Recasting the LHC dilepton searches in the fully inclusive category following Ref. Faroughy:2016osc yields constraints weaker than those coming from direct LQ searches shown in Fig. 29.

5.3 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.

5.3.1 Pair production

At hadron colliders, LQs can be pair-produced through either $gg$ or $q\bar{q}$ 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 CMS:2018sxp ; Khachatryan:2015vaa , second Aaboud:2019jcc ; Khachatryan:2015vaa ; Sirunyan:2018kzh and third generation Sirunyan:2018vhk ; Aaboud:2019bye ; Sirunyan:2018kzh LQs at the LHC. Given the model Lagrangian 5.118, we are interested in the final states containing either two charged leptons and two jets ($\ell\ell jj$), or two neutrinos and two jets ($\nu\nu jj$). Note that for the doublet LQ $\Omega=(\omega^{2/3},\omega^{-1/3})$, the jets will consist of down-type quarks, while for the singlet LQ $\chi^{-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 $\tau^{+}\tau^{-}b\bar{b}$ final state is somewhat stronger than that from $\tau^{+}\tau^{-}t\bar{t}$ final state Aaboud:2019bye ; Sirunyan:2018kzh .

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 $\mu^{+}\mu^{-}jj$ final states from LQ pair-production at LHC to impose stringent bounds on the $\lambda_{\alpha d}$ and $\lambda^{\prime}_{\alpha d}$ couplings (with $\alpha=e,\mu$) which are relevant for NSI involving electron and muon flavors. There is no dedicated search for LQs in the $\tau^{+}\tau^{-}jj$ channel to impose similar constraints on $\lambda_{\tau d}$ and $\lambda^{\prime}_{\tau d}$ relevant for tau-flavor NSI. There are searches for third generation LQ Sirunyan:2018vhk ; Aaboud:2019bye looking at $\tau^{+}\tau^{-}b\bar{b}$ and $\tau^{+}\tau^{-}t\bar{t}$ signatures which are not relevant for NSI, since we do not require $\lambda^{\prime}_{\tau t}$ (for $\chi^{-1/3}$) or $\lambda_{\tau b}$ (for $\omega^{2/3}$) couplings. For constraints on $\lambda_{\tau d}$, we recast the $\tau^{+}\tau^{-}b\bar{b}$ search limits Sirunyan:2018vhk ; Aaboud:2019bye ; Sirunyan:2018kzh taking into account the $b$-jet misidentification as light jets, with an average rate of $1.5\%$ (for a $b$-tagging efficiency of 70%) Chatrchyan:2012jua . As expected, this bound is much weaker, as shown in Fig. 29.

However, a stronger bound on NSI involving the tau-sector comes from $\nu\bar{\nu}jj$ final state. From the Lagrangian (5.118), we see that the same $\lambda_{\tau d}$ coupling that leads to $\tau^{+}\tau^{-}dd$ final state from the pair-production of $\omega^{2/3}$ also leads to $\nu_{\tau}\bar{\nu}_{\tau}dd$ final state from the pair-production of the $SU(2)_{L}$ partner LQ $\omega^{-1/3}$, whose mass cannot be very different from that of $\omega^{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 $\nu\bar{\nu}jj$ constraint will equally apply to all $\lambda_{\alpha d}$ (with $\alpha=e,\mu,\tau$) couplings which ultimately restrict the strength of tau-sector NSI, as we will see in the next subsection. The same applies to the $\lambda^{\prime}_{\tau d}$ couplings of the singlet LQ $\chi^{-1/3}$, which are also restricted by the $\nu\bar{\nu}jj$ constraint.

5.3.2 Single production

LQs can also be singly produced at the collider in association with charged leptons via $s$- and $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 mass-coupling 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 $\lambda_{ed}$ and $\lambda_{\mu 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 $\lambda_{ed}$, but not for $\lambda_{\mu d}$. There is no constraint in the $\tau j$ channel, and the derived constraint from $\tau b$ channel is too weak to appear in this plot.

5.3.3 How light can the leptoquark be?

There is a way to relax the $\nu\bar{\nu}jj$ constraint and allow for smaller LQ masses for the doublet components. This is due to a new decay channel $\omega^{-1/3}\to\omega^{2/3}+W^{-}$ which, if kinematically allowed, can be used to suppress the branching ratio of $\omega^{-1/3}\to\nu d$ decay for relatively smaller values of $\lambda_{\alpha d}$ couplings, thereby reducing the impact of the $\nu\bar{\nu}jj$ constraint. The partial decay widths for $\omega^{-1/3}\to\omega^{2/3}+W^{-}$ and $\omega^{-1/3}\to\nu_{\alpha}d_{\beta}$ are respectively given by

[TABLE]

In deriving Eq. (5.161), we have used the Goldstone boson equivalence theorem, and in Eq. (5.162), the factor in the denominator is not $8\pi$ (unlike the SM $h\to b\bar{b}$ case, for instance), because only one helicity state contributes.

The lighter LQ $\omega^{2/3}$ in this case can only decay to $\ell_{\alpha}d_{\beta}$ with 100% branching ratio. Using the fact that constraints from $\tau^{+}\tau^{-}jj$ channel are weaker, one can allow for $\omega^{2/3}$ as low as 522 GeV, as shown in Fig. 29 by the solid brown curve, when considering the $\lambda_{\tau d}$ coupling alone. This is, however, not applicable to the scenario when either $\lambda_{ed}$ or $\lambda_{\mu d}$ coupling is present, because of the severe constraints from $e^{+}e^{-}jj$ and $\mu^{+}\mu^{-}jj$ final states.

5.4 NSI prediction

The LQs $\omega^{-1/3}$ and $\chi^{-1/3}$ in the model have couplings with neutrinos and down-quark (cf. Eq. (5.118)), and therefore, induce NSI at tree level as shown in Fig. 30 via either $\lambda$ or $\lambda^{\prime}$ couplings. From Fig. 30, we can write down the effective four-fermion Lagrangian as

[TABLE]

where we have used Fierz transformation in the second step. Comparing Eq. (5.163) with Eq. (3.9), we obtain the NSI parameters

[TABLE]

For $Y_{n}(x)\equiv\frac{N_{n}(x)}{N_{p}(x)}=1$, one can obtain the effective NSI parameters from Eq. (3.13) as

$\displaystyle\varepsilon_{\alpha\beta}\ \equiv\ 3\varepsilon_{\alpha\beta}^{d}\ =\ \frac{3}{4\sqrt{2}\ G_{F}}\left(\frac{\lambda_{\alpha d}^{\star}\lambda_{\beta d}}{m_{\omega}^{2}}+\frac{\lambda_{\alpha d}^{\prime\star}\lambda^{\prime}_{\beta d}}{m_{\chi}^{2}}\right)\,.$

(5.165)

To satisfy the neutrino mass constraint [cf. Eq. (5.124)], we can have either $\lambda_{\alpha d}^{\star}\lambda_{\beta d}$ or $\lambda_{\alpha d}^{\prime\star}\lambda^{\prime}_{\beta d}$ of ${\cal O}(1)$, but not both simultaneously, for a given flavor combination $(\alpha,\beta)$. But we can allow for $\lambda_{\alpha d}^{\star}\lambda_{\beta d}$ and $\lambda_{\alpha^{\prime}d}^{\prime\star}\lambda^{\prime}_{\beta^{\prime}d}$ simultaneously to be of ${\cal O}(1)$ for either $\alpha\neq\alpha^{\prime}$ or $\beta\neq\beta^{\prime}$, which will be used below to avoid some experimental constraints for the maximum NSI predictions.

5.4.1 Doublet leptoquark

First, let us consider the doublet LQ contribution by focusing on the $\lambda$-couplings only. We show in Figs. 31 and 32 the predictions for diagonal $(\varepsilon_{ee},\,\varepsilon_{\mu\mu},\,\varepsilon_{\tau\tau})$ and off-diagonal $(\varepsilon_{e\mu},\,\varepsilon_{\mu\tau},\,\varepsilon_{e\tau})$ NSI parameters respectively from Eq. (5.165) by black dotted contours. Color-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 $\lambda_{\alpha d}<\sqrt{\frac{4\pi}{\sqrt{3}}}$ DiLuzio:2017chi . 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 Esteban:2018ppq . 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 $\lambda_{\alpha d}$ (which is independent of $\lambda_{\alpha d}$) and single-production mode for large $\lambda_{\alpha d}$) with $\alpha=e,\mu$. Here we have assumed 50% branching ratio to $ej$ or $\mu j$, and the other 50% to $\tau 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 $\nu\bar{\nu}jj$ channel, where the vertical dashed line indicates the limit assuming ${\rm BR}(\omega^{-1/3}\to\nu d)=100\%$, and the unshaded region to the left of this line for small $\lambda_{\tau d}$ is allowed by opening up the $\omega^{-1/3}\to\omega^{2/3}W^{-}$ channel (cf. Sec. 5.3.3). Note that we cannot completely switch off the $\omega^{-1/3}\to\nu d$ channel, because that would require $\lambda_{\tau d}\to 0$ and in this limit, the NSI will also vanish.

The red line in Fig. 31 (b) is the suggestive limit on $\varepsilon_{\alpha\beta}^{dR}$ from NuTeV data Davidson:2003ha (cf. Table XVII). This is not shaded because there is a $2.7\sigma$ discrepancy of their $s_{w}^{2}$ measurement with the PDG average Tanabashi:2018oca 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. Davidson:2003ha , but using the latest value of $s_{w}^{2}$ (on-shell) Tanabashi:2018oca (without including NuTeV). Specifically, we have used the NuTeV measurement of the effective coupling $\left(\tilde{g}^{\mu}_{R}\right)^{2}=0.0310\pm 0.0011$ from $\nu_{\mu}q\to\nu q$ scatterings Zeller:2001hh which is consistent with the SM prediction of $\left(\tilde{g}^{\mu}_{R}\right)_{\rm SM}^{2}=0.0297$. Here $\left(\tilde{g}^{\mu}_{R}\right)^{2}$ is defined as

[TABLE]

where $g_{R}^{u}=-\frac{2}{3}s_{w}^{2}$ and $g_{R}^{d}=\frac{1}{3}s_{w}^{2}$ 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.163) involves right-handed down-quarks for the doublet LQ component $\omega^{2/3}$. In Eq. (5.166), setting $\varepsilon_{\mu\mu}^{uR}=0$ for this LQ model and comparing $\left(\tilde{g}^{\mu}_{R}\right)^{2}$ with the measured value, we obtain a 90% CL on $\varepsilon_{\mu\mu}^{dR}<0.029$, which should be multiplied by 3 (since $\varepsilon_{\alpha\beta}\equiv 3\varepsilon_{\alpha\beta}^{dR}$) to get the desired constraint on $\varepsilon_{\alpha\beta}$ shown in Fig. 31 (b).

Also note that unlike in the Zee model case discussed earlier, the IceCube limit on $|\varepsilon_{\tau\tau}-\varepsilon_{\mu\mu}|$ Esmaili:2013fva is not shown in Figs. 31 (b) and (c). This is because the NSI parameters in the LQ model under consideration receive two contributions as shown in Eq. (5.165). Although we cannot have both $\lambda$ and $\lambda^{\prime}$ contributions large for the same $\varepsilon_{\alpha\beta}$, it is possible to have a large $\lambda$ contribution to $\varepsilon_{\alpha\beta}$ and a large $\lambda^{\prime}$ contribution to $\varepsilon_{\alpha^{\prime}\beta^{\prime}}$ (with either $\alpha\neq\beta$ or $\beta\neq\beta^{\prime}$), thus evading the cLFV constraints (which are only applicable to either $\lambda$ or $\lambda^{\prime}$ sectors), as well as the IceCube constraint on $|\varepsilon_{\tau\tau}-\varepsilon_{\mu\mu}|$, which is strictly applicable only in the limit of all $\varepsilon_{e\alpha}\to 0$. This argument can be applied to all the LQ models discussed in subsequent sections, with a few exceptions, when the NSI arises from only one type of couplings; see e.g. Eq. (7.200) and (7.205)). So we will not consider the IceCube limit on $\varepsilon_{\mu\mu}$ and $\varepsilon_{\tau\tau}|$ for our LQ NSI analysis, unless otherwise specified.

For $\varepsilon_{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 $\varepsilon_{\mu\mu}$, the most stringent limit of 8.6% comes from NuTeV. However, if this constraint is not considered, $\varepsilon_{\mu\mu}$ can be as large as 21.6%. Similarly, $\varepsilon_{\tau\tau}$ 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 dev_pondd , as shown in Fig. 31 (c).

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, $\tau\to\ell\pi^{0}$ and $\tau\to\ell\eta$ (with $\ell=e,\mu$) impose strong constraints (cf. Sec. 5.1.5) on $\varepsilon_{\mu\tau}$ and $\varepsilon_{e\tau}$, as shown in Figs. 32 (a) and (b). For $\varepsilon_{e\mu}$, the most stringent limit comes from $\mu-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 Davidson:2003ha , COHERENT Coloma:2017ncl and IceCube Salvado:2016uqu , as well as the global-fit constraints from neutrino oscillation+COHERENT data Esteban:2018ppq and future DUNE sensitivity dev_pondd . 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 Gavela:2008ra .

5.4.2 Singlet leptoquark

Now if we take the $\lambda^{\prime}$ couplings instead of $\lambda$ in Eq. (5.165), 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.135)), as well as the LEP and HERA contact interaction constraints on $\varepsilon_{ee}$ (cf. Table XVI) are somewhat modified. In addition, there are new constraints from $D^{+}\to\pi^{+}\ell^{+}\ell^{-}$ and $D^{0}\to\ell^{+}\ell^{-}$ (cf. Sec. 5.1.6) for $\varepsilon_{ee}$ and $\varepsilon_{\mu\mu}$, as shown in Fig. 33 (a) and (b). For $\varepsilon_{ee}$, the $D^{+}\to\pi^{+}e^{+}e^{-}$ constraint turns out to be much weaker than the APV constraint. The $D^{0}\to e^{+}e^{-}$ constraint is even weaker and does not appear in Fig. 33 (a). However, for $\varepsilon_{\mu\mu}$, the $D^{+}\to\pi^{+}\mu^{+}\mu^{-}$ constraint turns out to be the strongest, limiting the maximum allowed value of $\varepsilon_{\mu\mu}$ 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 $\chi$ couples to left-handed quarks (cf. Eq. (5.163)). So it will affect the effective coupling $\left(\tilde{g}^{\ell}_{L}\right)$. For $\varepsilon_{\mu\mu}$, we use the NuTeV measurement of $\left(\tilde{g}^{\mu}_{L}\right)^{2}=0.3005\pm 0.0014$ from $\nu_{\mu}q\to\nu q$ scatterings Zeller:2001hh which is $2.7\sigma$ smaller than the SM prediction of $\left(\tilde{g}^{\mu}_{L}\right)_{\rm SM}^{2}=0.3043$. Here $\left(\tilde{g}^{\mu}_{L}\right)^{2}$ is defined as

[TABLE]

where $g_{L}^{u}=\frac{1}{2}-\frac{2}{3}s_{w}^{2}$ and $g_{L}^{d}=-\frac{1}{2}+\frac{1}{3}s_{w}^{2}$. For the SM prediction, we have used the latest PDG value for on-shell $s_{w}^{2}=0.22343$ from a global-fit to electroweak data (without NuTeV) Tanabashi:2018oca . Comparing $\left(\tilde{g}^{\mu}_{L}\right)^{2}$ with the measured value, we derive a 90% CL constraint

on $0.0018<\varepsilon_{\mu\mu}<0.8493$. Note that this prefers a non-zero $\varepsilon_{\mu\mu}$ at 90% CL ($1.64\sigma$) because the SM with $\varepsilon_{\mu\mu}=0$ is $2.7\sigma$ away and also because there is a cancellation between $g_{L}^{d}$ (which is negative) and $\varepsilon_{\mu\mu}$ (which is positive) in Eq. (5.167) to lower the value of $\left(\tilde{g}^{\mu}_{L}\right)^{2}$ to within $1.64\sigma$ of the measured value.

For the off-diagonal sector, there are new constraints from $\tau\to\ell\gamma$ relevant for $\varepsilon_{\mu\tau}$ and $\varepsilon_{e\tau}$, as shown in Figs. 33 (c) and (d). However, these are less stringent than the $\tau\to\ell\pi^{0}$ and $\tau\to\ell\eta$ constraints discussed before. There are no new constraints for $\varepsilon_{\tau\tau}$ and $\varepsilon_{e\mu}$ 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.

6 NSI in a triplet leptoquark model

This is the $\mathcal{O}_{3}^{9}$ model of Table III Cai:2014kra . In this model, two new fields are introduced – an $SU(2)_{L}$-triplet scalar LQ $\bar{\rho}\left({\bf\bar{3}},{\bf 3},\frac{1}{3}\right)=\left(\bar{\rho}^{4/3},\,\bar{\rho}^{1/3},\,\bar{\rho}^{-2/3}\right)$ and an $SU(2)_{L}$-doublet LQ $\Omega\left({\bf 3},{\bf 2},\frac{1}{6}\right)=\left(\omega^{2/3},\,\omega^{-1/3}\right)$. The relevant Lagrangian for the neutrino mass generation can be written as

[TABLE]

These interactions, along with the potential term

[TABLE]

where $\bar{\rho}$ is related to $\rho$ by charge conjugation as $\rho\left({\bf{3}},{\bf 3},-\frac{1}{3}\right)=\left(\rho^{2/3},\,-\rho^{-1/3},\,\rho^{-4/3}\right)$, induce neutrino mass at one-loop level via the ${\cal O}_{3}^{9}$ operator in the notation of Ref. Cai:2014kra , as shown in Fig. 34. The neutrino mass matrix can be estimated as

[TABLE]

where $M_{d}$ is the diagonal down-type quark mass matrix and $M\equiv{\rm max}(m_{\omega},m_{\rho})$. The NSI parameters read as

$\displaystyle\varepsilon_{\alpha\beta}\ =\ \frac{3}{4\sqrt{2}G_{F}}\left(\frac{\lambda_{\alpha d}^{\star}\lambda_{\beta d}}{m_{\omega}^{2}}+\frac{\lambda_{\alpha u}^{\prime\star}\lambda^{\prime}_{\beta u}}{m_{\rho^{-2/3}}^{2}}+\frac{\lambda_{\alpha d}^{\prime\star}\lambda^{\prime}_{\beta d}}{2m_{\rho^{1/3}}^{2}}\right)\,.$

(6.171)

Note that both $\lambda$ and $\lambda^{\prime}$ cannot be large at the same time due to neutrino mass constraints (cf. Eq. (6.170)). For $\lambda\gg\lambda^{\prime}$, this expression is exactly the same as the doublet LQ contribution derived in Eq. (5.165) and the corresponding maximum NSI can be read off from Table XVII for the doublet component.

On the other hand, for $\lambda^{\prime}\gg\lambda$, the third term in Eq. (6.171) is analogous to the down-quark induced singlet LQ NSI given in Eq. (5.165) (except for the Clebsch-Gordan factor of $(1/\sqrt{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 $\lambda^{\prime}_{\alpha u}=\lambda^{\prime}_{\alpha d}$ in the Lagrangian (6.168). 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:

6.1 Atomic parity violation

The shift in the weak charge given by Eq. (5.130) is modified to

[TABLE]

Assuming $m_{\rho^{1/3}}=m_{\rho^{4/3}}\equiv m_{\rho}$ and noting that $\lambda^{\prime}_{\alpha u}=\lambda^{\prime}_{\alpha d}$ in Eq. (6.168), we obtain

[TABLE]

Comparing this with the $2\sigma$ allowed range (5.134), we obtain the modified constraint

[TABLE]

which is weaker (stronger) than that given by Eq. (5.135) for the $SU(2)_{L}$-doublet (singlet) LQ alone.

6.2 $\mu-e$ conversion

From Eq. (5.136), we see that for the triplet case, the rate of $\mu-e$ conversion will be given by

[TABLE]

For degenerate $\rho$-mass and $\lambda^{\prime}_{\ell d}=\lambda^{\prime}_{\ell u}$, we obtain the rate to be $(3/2)^{2}$ times larger than that given in Eq. (5.136). Therefore, the constraints on $|\lambda^{\prime\star}_{ed}\lambda^{\prime}_{\mu d}|$ given in Table XII will be a factor of 3/2 stronger.

6.3 Semileptonic tau decays

The semileptonic tau decays such as $\tau^{-}\to\ell^{-}\pi^{0},~{}\ell^{-}\eta,~{}\ell^{-}\eta^{\prime}$ will have two contributions from $\bar{\rho}^{1/3}$ and $\bar{\rho}^{4/3}$. The relevant terms in the Lagrangian (7.208) are

[TABLE]

where we have assumed a basis with diagonal down-type quark sector. Using the matrix element (5.145), we find the modified decay rate for $\tau^{-}\to\ell^{-}\pi^{0}$ from Eq. (5.143):

[TABLE]

Thus, for $m_{\rho^{-1/3}}=m_{\rho^{4/3}}$, the $\tau^{-}\to\ell^{-}\pi^{0}$ decay rate is suppressed by a factor of 1/4, compared to the doublet or singlet LQ case (cf. Eq. (5.143)). So the constraints on $\lambda^{\prime}_{\ell d}\lambda_{\tau d}^{\star}$ from $\tau\to\ell\pi^{0}$ shown in Table XIV will be a factor of 2 weaker in the triplet LQ case.

On the other hand, using the matrix element (5.146), we find that the modified decay rate for $\tau^{-}\to\ell^{-}\eta$ becomes

[TABLE]

which is enhanced by a factor of 9/4 for $m_{\rho^{-1/3}}=m_{\rho^{4/3}}$, compared to the doublet or singlet LQ case. So the constraints on $\lambda_{\ell d}\lambda_{\tau d}^{\star}$ from $\tau\to\ell\eta$ shown in Table XIV will be a factor of 3/2 stronger in the triplet LQ case. The same scaling behavior applies to $\tau\to\ell\eta^{\prime}$ constraints. These modified constraints are summarized in Table XVIII.

6.4 $\ell_{\alpha}\to\ell_{\beta}+\gamma$

The cLFV decay $\ell_{\alpha}\to\ell_{\beta}+\gamma$ arises via one-loop diagrams with the exchange of $\bar{\rho}$ LQ fields, analogous to Fig. 25. The relevant couplings in Eq. (6.168) have the form $\ell u\bar{\rho}^{1/3}=\overline{u^{c}}P_{L}\ell\bar{\rho}^{1/3}$ for which $Q_{F}=-2/3$ and $Q_{B}=1/3$ in the general formula (4.45), whereas for the couplings $\ell d\bar{\rho}^{4/3}=\overline{d^{c}}P_{L}\ell\bar{\rho}^{4/3}$, we have $Q_{F}=1/3$ and $Q_{B}=4/3$. Substituting these charges in Eq. (4.45) and taking the limit $t=m_{F}^{2}/m_{B}^{2}\to 0$ (since the LQs are expected to be much heavier than the SM charged leptons), we obtain

[TABLE]

where $9=3^{2}$ is a color factor and we have assumed $m_{\rho^{-1/3}}=m_{\rho^{4/3}}$. The rate in Eq. (6.179) is 9/4 times larger than that given in Eq. (5.142) for the singlet LQ case. Therefore, the constraints on $|\lambda^{\prime}_{\beta d}\lambda_{\alpha d}^{\prime\star}|$ derived in Table XIII will be weakened by a factor of 3/2.

6.5 $D$-meson decays

The $\ell_{\alpha}u_{\beta}\bar{\rho}^{1/3}$ and $\ell_{\alpha}d_{\beta}\bar{\rho}^{4/3}$ terms in Eq. (7.208) 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 $\ell_{\alpha}u_{\beta}\bar{\rho}^{1/3}$ term in Eq. (7.208) now becomes $\ell_{\alpha}V^{\star}_{id}u_{i}\bar{\rho}^{1/3}$ which induces $D^{0}\to\ell^{+}\ell^{-}$ and $D^{+}\to\pi^{+}\ell^{+}\ell^{-}$ decays. The analysis will be the same as in Sec. 5.1.6, except that the $\lambda^{\prime}_{\alpha d}$ couplings will now be replaced by $\lambda^{\prime}_{\alpha d}/\sqrt{2}$. Correspondingly, the constraints on $|\lambda^{\prime}_{\alpha d}|$ given in Table XV will be $\sqrt{2}$ times weaker. For instance,

[TABLE]

6.6 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 $\sqrt{3/2}$ times that of the singlet case in Table XVI. The modified constraint is given by

[TABLE]

6.7 LHC constraints

The LHC constraints on the $\bar{\rho}$ fields will be similar to the discussion in Sec. 5.3. Comparing the Lagrangians (5.118) and (7.208), we see that $\bar{\rho}^{1/3}$ will have the same decay modes to $\nu j$ and $\ell j$, and therefore, the same constraints as the singlet $\chi^{-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 $\bar{\rho}^{-2/3}$ component can in principle be lighter, since it can only decay to $\nu j$ for which the constraints are weaker (cf. Fig. 29). However, the mass splitting between $\bar{\rho}^{-2/3}$ and $\bar{\rho}^{1/3}$ cannot be more than $\sim 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_{\rho^{1/3}}$ for 50% branching ratio to $\nu j$ and $\ell j$ channels (since they are governed by the same $\lambda^{\prime}_{\alpha d}$ coupling), one can allow for $m_{\rho^{-2/3}}$ as low as 800 GeV or so.

6.8 NSI prediction

Taking into account all the constraints listed above, we show in Figs. 35 and 36 the predictions for diagonal $(\varepsilon_{ee},\,\varepsilon_{\mu\mu},\,\varepsilon_{\tau\tau})$ and off-diagonal $(\varepsilon_{e\mu},\,\varepsilon_{\mu\tau},\,\varepsilon_{e\tau})$ NSI parameters respectively from Eq. (6.171) by black dotted contours. Color-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 $\varepsilon^{uL}_{\mu\mu}$ and $\varepsilon^{dL}_{\mu\mu}$ as in this LQ model (cf. (6.168)), the total contribution to $\left(\tilde{g}^{\mu}_{L}\right)^{2}$ in Eq. (5.167) is always positive, and therefore, any nonzero $\varepsilon_{\mu\mu}$ will make the discrepancy worse than the SM case of $2.7\sigma$. Therefore, we cannot impose a 90% CL ($1.64\sigma$) constraint from NuTeV in this scenario. The line shown in Fig. 35 (b) corresponds to the $3\sigma$ constraint on $\varepsilon_{\mu\mu}<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. 35 and 36, we find the maximum allowed values of the NSI parameters in the triplet LQ model to be

[TABLE]

This is also summarized in Fig. 59 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

7.1.1 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.161616There is yet another possibility where the mediators could be new vector bosons; however, this necessarily requires some new gauge symmetry and other associated Goldstone bosons to cancel the UV divergences. The low-energy effective operator that leads to neutrino mass in this model is the dimension-7 operator

[TABLE]

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 Dev:2012sg . In the usual inverse seesaw model Mohapatra:1986bd , 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., $\mu_{S}SS$ leads to a tree-level neutrino mass via the standard inverse seesaw mechanism Mohapatra:1986bd . However, if one imposes a global $U(1)$ symmetry under which the $S$-field is charged, then the $\mu_{S}SS$ term can be explicitly forbidden at tree-level.171717This can be done, for instance, by adding a singlet scalar field $\sigma$ with a global $U(1)$ charge of +2, and by making $N$ and $S$ oppositely charged under this $U(1)$, viz., $N(-1)$ and $S(+1)$, so that the $S\sigma S$ term is forbidden, but $N\sigma N$ and $\overline{S}\sigma N$ are allowed. Furthermore, this global $U(1)$ symmetry can be gauged, e.g., in an $E_{6}$ GUT embedding, where the fundamental representation 27 breaks into ${\bf 16}_{1}+{\bf 10}_{-2}+{\bf 1}_{4}$ under $SO(10)\times U(1)$. The $\nu$ and $N$ belong to the ${\bf 16}_{1}$ subgroup, while the $S$ belongs to ${\bf 1}_{4}$. Adding two scalars $\sigma,\sigma^{\prime}$ with $U(1)$ charges $-2$ and $-5$ respectively allows the Dirac mass term $\overline{N}\sigma S$ and Majorana mass term $N\sigma^{\prime}N$ in Eq. (7.188), but not the Majorana mass terms $S\sigma^{(^{\prime})}S$. In this case, the only lepton number breaking term that is allowed is the Majorana mass term for the $N$-field, i.e., $\mu_{R}NN$. 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 Dev:2012sg ). One can see that the low-energy effective operator that leads to neutrino mass in this model is the $d=7$ operator ${\cal O}^{\prime}_{1}$ of Eq. (1.5) by cutting Fig. 37 at one of the $H$-legs in the loop.

The relevant part of the Yukawa Lagrangian of this model is given by

[TABLE]

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 $\mu_{R}\ll M_{N}$) Dev:2012sg ; Dev:2012bd :

[TABLE]

where $M_{D}\equiv Yv/\sqrt{2}$, $\alpha_{w}\equiv g^{2}/4\pi$, $x_{N}=m_{N}^{2}/m_{W}^{2}$, $x_{H}=m_{H}^{2}/m_{W}^{2}$ and $x_{Z}=m_{Z}^{2}/m_{W}^{2}$, and we have assumed $M_{N}=m_{N}{\bf 1}$ for simplicity.

The NSI in this model arise due to the fact that the light $SU(2)_{L}$-doublet neutrinos $\nu$ mix with the singlet fermions $N$ and $S$, due to which the $3\times 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 Blennow:2016jkn . If only one extra Dirac state with mass larger than $\sim$ GeV (such that it cannot be produced in accelerator neutrino oscillation experiments, such as DUNE) mixes with the three light states with mixing parameters $U_{\alpha 4}$ (with $\alpha=e,\mu,\tau$), we can write the NSI parameters as

$\displaystyle\varepsilon_{ee}\ =$ $\displaystyle\ \left(\frac{Y_{n}}{2}-1\right)|U_{e4}|^{2},\qquad\varepsilon_{\mu\mu}\ =\ \frac{Y_{n}}{2}|U_{\mu 4}|^{2},\qquad\varepsilon_{\tau\tau}\ =\ \frac{Y_{n}}{2}|U_{\tau 4}|^{2},$

$\displaystyle\varepsilon_{e\mu}\ =$ $\displaystyle\ \frac{1}{2}\left(Y_{n}-1\right)U_{e4}U^{\star}_{\mu 4},\qquad\varepsilon_{e\tau}\ =\ \frac{1}{2}\left(Y_{n}-1\right)U_{e4}U^{\star}_{\tau 4},\qquad\varepsilon_{\mu\tau}\ =\ \frac{Y_{n}}{2}U_{\mu 4}U^{\star}_{\tau 4}\,.$

(7.190)

Here $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}\to 1$ which is approximately true for neutrino propagation in earth matter, we get vanishing $\varepsilon_{e\mu}$ and $\varepsilon_{e\tau}$ up to second order in $U_{\alpha 4}$.181818This result is in disagreement with Ref. Blennow:2016jkn , where they have $\varepsilon_{\alpha\beta}=\frac{1}{2}U_{\alpha 4}U^{\star}_{\beta 4}$ for all the off-diagonal NSI parameters, which cannot be the case, because for $\alpha=e$, both CC and NC contributions are present, whereas for $\alpha\neq e$, only the NC contribution matters. Taking into account all the experimental constraints on $U_{\alpha 4}U_{\beta 4}^{\star}$ 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 Atre:2009rg ; Deppisch:2015qwa ; deGouvea:2015euy ; Alonso:2012ji ; Blennow:2016jkn ; Antusch:2014woa , the maximum NSI parameters allowed in this model are summarized in Table XX. We find that Blennow:2016etl

[TABLE]

For $\varepsilon_{e\mu}$ and $\varepsilon_{e\tau}$, we have used $Y_{n}=1.051$ (for average value all over the earth) in Eq. (7.190), in addition to the cLFV constraints on $U_{e4}U_{\mu 4}^{\star}$ and $U_{e4}U_{\tau 4}^{\star}$. The maximum NSI values listed above (and also summarized in Table XX) are obtained for a relatively heavy sterile neutrino (with mass larger than the electroweak-scale), so that the stringent low-energy constraints from beam dump and meson decays can be avoided, and the only relevant constraint comes from the electroweak precision data Antusch:2014woa .

The NSI expressions (7.190) also apply to two-loop radiative models with two $W$-boson exchange Petcov:1984nz ; Babu:1988ig ; Davidson:2006tg . However, the maximum NSI obtainable in these models will be much smaller than the estimate in Eq. (7.191) because the sterile neutrino in this case is required to be heavier for successful neutrino mass generation at two-loop.

7.1.2 One-loop model with vectorlike leptons

This model Cai:2014kra utilizes the same $d=7$ operator ${\cal O}_{2}=L^{i}L^{j}L^{k}e^{c}H^{l}\epsilon_{ij}\epsilon_{kl}$ (cf. Eq. (2.7b)), as in the Zee model to generate a one-loop neutrino mass. The new particles added are a scalar singlet $\eta^{+}({\bf 1},{\bf 1},1)$ and a vectorlike lepton $\psi\left({\bf 1},{\bf 2},-\frac{3}{2}\right)=(E,\,F^{--})$, which give rise to the ${\cal O}_{2}^{1}$ 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:

[TABLE]

where $\psi^{c}=(F^{++},\,-E^{c})$ and $H\left({\bf 1},{\bf 2},\frac{1}{2}\right)$ is the SM Higgs doublet. Expanding the first two terms, we get

[TABLE]

The neutrino mass matrix can be estimated as

[TABLE]

where $M_{\ell}$ 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\equiv{\rm max}(m_{\eta},m_{E_{i}})$. Note that just one flavor of $\psi$ is not sufficient, because in this case, the neutrino mass matrix (7.194) would have a flavor structure given by $(fM_{\ell}-M_{\ell}f)$, which has all the diagonal entries zero, similar to the Zee-Wolfenstein model Wolfenstein:1980sy . Such a structure is ruled out by observed neutrino oscillation data. Thus, we require at least two flavors of $\psi$, in which case the diagonal entries of $M_{\nu}$ are nonzero, and the model is consistent with experiments.

NSI in this model are induced by the $f$-type couplings in Eq. (7.193), similar to the $f$-couplings in the Zee model Lagrangian (4.15). The NSI parameters read as

$\displaystyle\varepsilon_{\alpha\beta}\ \equiv\ \varepsilon_{\alpha\beta}^{ee}\ =\ \frac{1}{\sqrt{2}G_{F}}\frac{f_{e\alpha}^{\star}f_{e\beta}}{m_{\eta^{+}}^{2}}\,.$

(7.195)

Due to the antisymmetric nature of the $f$ couplings, the only relevant NSI parameters in this case are $\varepsilon_{\mu\tau}$, $\varepsilon_{\mu\mu}$, and $\varepsilon_{\tau\tau}$. These are severely constrained by cLFV searches and universality of charged currents Nebot:2007bc , 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 $\sim 100$ GeV to satisfy the LEP constraints (cf. Sec. 4.7), we obtain from Eq. (7.195) and Table XIX the following maximum values:

[TABLE]

This is also summarized in Table XX.

7.1.3 $SU(2)_{L}$-singlet leptoquark model with vectorlike quark

This model Cai:2014kra is the $\mathcal{O}_{3}^{4}$ realization of the dimension-7 operator ${\cal O}_{3}$ (cf. Table III). The new particles introduced are a scalar LQ singlet $\chi\left({\bf 3},{\bf 1},-\frac{1}{3}\right)$ and a vectorlike quark doublet ${\cal Q}\left({\bf 3},{\bf 2},-\frac{5}{6}\right)=\left(D^{-1/3},\,X^{-4/3}\right)$. Neutrino mass is generated at one-loop level as shown in Fig. 39. The $QQ\chi^{\star}$ and $d^{c}u^{c}\chi$ 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

[TABLE]

Expanding the first two terms, we get

[TABLE]

The neutrino mass matrix can be estimated as

[TABLE]

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\equiv{\rm max}(m_{\chi},m_{D_{i}})$. With a single copy of VQ quarks, the rank of $M_{\nu}$ 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 $\lambda$-type interactions in Eq. (7.198):

$\displaystyle\varepsilon_{\alpha\beta}\ =\ \frac{3}{4\sqrt{2}G_{F}}\frac{\lambda_{\alpha d}^{\star}\lambda_{\beta d}}{m_{\chi}^{2}}\,.$

(7.200)

This is similar to the singlet LQ contribution in Eq. (5.165), with the important exception that the NSI get contribution only from the $\lambda$-couplings, and therefore, the IceCube limits on $|\varepsilon_{\mu\mu}-\varepsilon_{\tau\tau}|<9.3\%$ cannot be avoided, just like in the Zee model case. The corresponding maximum NSI can be read off from Table XVII, except for $\varepsilon_{\tau\tau}$:

[TABLE]

This is also summarized in Table XX.

7.1.4 $SU(2)_{L}$ -doublet leptoquark model with vectorlike quark

This is referred to as $\mathcal{O}_{3}^{6}$ in Table III. The model has an $SU(2)_{L}$-doublet LQ $\Omega\left({\bf 3},{\bf 2},\frac{1}{6}\right)=\left(\omega^{2/3},\,\omega^{-1/3}\right)$ and an $SU(2)_{L}$-triplet vectorlike quark $\Sigma\left({\bf 3},{\bf 3},\frac{2}{3}\right)=\left(Y^{5/3},\,U^{2/3},\,D^{-1/3}\right)$. 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

[TABLE]

where $\widetilde{\Omega}=i\tau_{2}\Omega^{\star}$ is the isospin conjugate field. Expanding the terms in Eq. (7.202), we obtain

[TABLE]

The neutrino mass can be estimated as

[TABLE]

where $M_{d}$ and $M_{D}$ are the diagonal down quark mass matrix and vectorlike quark mass matrix respectively, and $M\equiv{\rm max}(m_{\omega},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_{\nu}$ 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 $\omega^{-1/3}$. The NSI parameters read as

$\displaystyle\varepsilon_{\alpha\beta}\ =\ \frac{3}{4\sqrt{2}G_{F}}\frac{\lambda_{\alpha d}^{\star}\lambda_{\beta d}}{m_{\omega}^{2}}\,.$

(7.205)

This expression is similar to the doublet LQ contribution in Eq. (5.165), with the exception that the IceCube limits on $|\varepsilon_{\mu\mu}-\varepsilon_{\tau\tau}|<9.3\%$ cannot be avoided. The corresponding maximum NSI can be read off from Table XVII, except for $\varepsilon_{\mu\mu}$ and $\varepsilon_{\tau\tau}$:

[TABLE]

This is also summarized in Table XX.

7.1.5 Model with $SU(2)_{L}$-triplet leptoquark and vectorlike quark

This is based on the operator ${\cal O}_{3}^{5}$ (see Table III) which is realized by adding an $SU(2)_{L}$-triplet $\bar{\rho}\left({\bf\bar{3}},{\bf 3},\frac{1}{3}\right)=\left(\bar{\rho}^{4/3},\,\bar{\rho}^{1/3},\,\bar{\rho}^{-2/3}\right)$ and a vectorlike quark doublet $\mathcal{Q}\left({\bf 3},{\bf 2},-\frac{5}{6}\right)=\left(D^{-1/3},\,X^{-4/3}\right)$. Neutrino mass is generated at one-loop level, as shown as Fig. 41. There is also a two-loop diagram involving $\rho^{2/3}$, which is not considered here, as that would be sub-dominant to the one-loop diagram. The interaction term $QQ\rho$ is forbidden by demanding baryon-number conservation to avoid proton decay. The relevant Lagrangian for the neutrino mass generation can be written as

[TABLE]

where $\bar{\rho}$ is related to $\rho$ by charge conjugation as $\rho\left({\bf{3}},{\bf 3},-\frac{1}{3}\right)=\left(\rho^{2/3},\,-\rho^{-1/3},\,\rho^{-4/3}\right)$. Expanding the terms in Eq. (7.207), we get

[TABLE]

The neutrino mass can be estimated as

[TABLE]

where $M_{d}$ and $M_{D}$ are the diagonal mass matrices for down-type quark and vectorlike quark fields, and $M={\rm max}(m_{D_{i}},m_{\rho})$, with $m_{D_{i}}$ being the eigenvalues of $M_{D}$. With a single copy of the vectorlike quark, the matrices $y$ and $\lambda^{\prime}$ are $3\times 1$ dimensional. Consequently the rank of $M_{\nu}$ 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. Zee:1985id ; Babu:1988ki .

NSI in this model are induced by both $\bar{\rho}^{-2/3}$ and $\bar{\rho}^{1/3}$ fields, which couple to up and down quarks respectively (cf. Eq. (7.208)). The NSI parameters read as

$\displaystyle\varepsilon_{\alpha\beta}\ =\ \frac{3}{4\sqrt{2}G_{F}}\left(\frac{\lambda_{\alpha u}^{\star}\lambda_{\beta u}}{m_{\rho^{-2/3}}^{2}}+\frac{\lambda_{\alpha d}^{\star}\lambda_{\beta d}}{2m_{\rho^{1/3}}^{2}}\right)\,.$

(7.210)

This is same as the triplet contribution in Eq. (6.171) and the maximum allowed values are given in Eq. (6.186).

7.1.6 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

[TABLE]

which may appear to be a product of $\mathcal{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: $\delta\left({\bf 3},{\bf 2},\frac{7}{6}\right)=\left(\delta^{5/3},\,\delta^{2/3}\right),\,\bar{\rho}\left(\bar{\bf 3},{\bf 3},\frac{1}{3}\right)=\left(\bar{\rho}^{4/3},\bar{\rho}^{1/3},\bar{\rho}^{-2/3}\right)\,,\xi\left({\bf 3},{\bf 1},\frac{2}{3}\right)$. 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

[TABLE]

Neutrino mass is generated by the diagram shown in Fig. 42 using the Lagrangian (7.212), together with the potential terms

[TABLE]

where $\widetilde{H}=(H^{0},\,-H^{-})$ represents the SM Higgs doublet. The neutrino mass matrix can be estimated as

[TABLE]

where $m_{1}$ and $m_{2}$ are the masses of the heaviest two LQs among the $\delta$, $\bar{\rho}$ and $\xi$ fields, and $M_{u}$ is the diagonal mass matrix in the up-quark sector. To get small neutrino masses, we need the product $\lambda\lambda^{\prime}\ll 1$. We may take $\lambda\sim{\cal O}(1)$ and $\lambda^{\prime}\ll\lambda$ which is preferable to the other case of $\lambda\ll\lambda^{\prime}$, since the $\lambda^{\prime}$ couplings are constrained by $D$-meson decays (see Sec. 5.1.6).

After integrating out the heavy LQ fields, Eq. (7.212) leads to an effective NSI Lagrangian with up-quarks in the neutrino propagation through matter. The NSI parameters read as

$\displaystyle\varepsilon_{\alpha\beta}\ =\ \frac{3}{4\sqrt{2}G_{F}}\left(\frac{\lambda_{\alpha u}^{\star}\lambda_{\beta u}}{m_{\delta}^{2}}+\frac{\lambda_{\alpha u}^{\prime\star}\lambda^{\prime}_{\beta u}}{m_{\rho^{-2/3}}^{2}}+\frac{\lambda_{\alpha d}^{\prime\star}\lambda^{\prime}_{\beta d}}{2m_{\rho^{1/3}}^{2}}\right)\,.$

(7.215)

For $\lambda\gg\lambda^{\prime}$, this expression is exactly the same as the doublet LQ contribution derived in Eq. (5.165) and the corresponding maximum NSI can be read off from Table XVII for the doublet component. For $\lambda^{\prime}\gg\lambda$, Eq. (7.215) is the same as Eq. (7.210). This latter choice maximizes NSI in this model and is summarized in Table XX.

There are other variations of one-loop LQ models with more exotic particles Dorsner:2017wwn ; AristizabalSierra:2007nf , where the neutrino mass is proportional to up-type quark mass. The NSI predictions in these models are the same as in Eq. (7.215).

7.2 Two-loop models

7.2.1 Zee-Babu model

This model realizes the operator $\mathcal{O}_{9}$ of Eq. (1.4). In this model Zee:1985id ; Babu:1988ki , two $SU(2)_{L}$-singlet Higgs fields, $h^{+}({\bf 1},{\bf 1},1)$ and $k^{++}({\bf 1},{\bf 1},2)$, are introduced. The corresponding Lagrangian for the generation of neutrino mass reads:

[TABLE]

Majorana neutrino masses are induced at two-loop as shown in Fig. 43 by the Lagrangian (7.216), together with the potential term

[TABLE]

The neutrino mass matrix reads:

[TABLE]

where $M={\rm max}(m_{k^{++}},m_{h^{+}})$ and $\cal{I}$ is a dimensionless function that depends on the ratio of the masses of the two new scalars Babu:2002uu ; Nebot:2007bc ; Schmidt:2014zoa ; Herrero-Garcia:2014hfa . The singly charged scalar $h^{+}$ induces NSI at tree-level through the $f$-type Yukawa coupling in Eq. (7.216). After integrating out the heavy scalars, NSI induced in neutrino propagation through normal matter can be written as in Eq. (7.195), with the replacement $m_{\eta^{+}}\to m_{h^{+}}$, and the maximum NSI are given by Eq. (7.196). These are severely constrained by cLFV searches and universality of charged currents Nebot:2007bc (cf. Table XIX), restricting the maximum NSI to ${\cal O}(10^{-3})$ level Ohlsson:2009vk . These numbers are summarized in Table XX.

7.2.2 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. Here the effective operator is of dimension nine, given by

[TABLE]

In addition to the SM fields, this model Kohda:2012sr employs a scalar LQ $\chi\left({\bf 3},{\bf 1},-\frac{1}{3}\right)$ and a scalar DQ $\Delta\left({\bf 6},{\bf 1},-\frac{2}{3}\right)$. The $\chi~{}(\Delta)$ 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

[TABLE]

Neutrino mass is generated at two-loop via the Lagrangian (7.220) in combination with the potential term

[TABLE]

The neutrino mass matrix can be calculated as

[TABLE]

where $M\equiv$ max($m_{\chi},m_{\Delta}$), $M_{d}$ is the diagonal down-type quark mass matrix, and ${\cal I}$ is a dimensionless two-loop integral defined in terms of the ratio of $m_{\Delta}^{2}$ and $m_{\chi}^{2}$ Babu:2002uu . After integrating out the heavy scalars, the NSI parameters in this model are given by Eq. (7.200) and the corresponding maximum NSI are given in Eq. (7.201). This is also summarized in Table XX.

There are a few variants of this LQ/DQ version of the Zee-Babu model. First, one could replace the color sextet field $\Delta\left({\bf 6},{\bf 1},-\frac{2}{3}\right)$ by a color triplet scalar $\Delta\left({\bf 3},{\bf 1},-\frac{2}{3}\right)$ in Fig. 44. The cubic term $\chi^{\star}\chi^{\star}\Delta$ will not be allowed by Bose symmetry in this case. By assuming two copies of the $\chi$ field, namely, $\chi_{1}$ and $\chi_{2}$, one could restore this coupling from $\chi_{1}^{\star}\chi_{2}^{\star}\Delta$, in which case the diagram of Fig. 44 can be connected Klein:2019jgb . 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 $\chi\left({\bf 3},{\bf 1},-\frac{1}{3}\right)$ by $\rho\left({\bf 3},{\bf 3},-\frac{1}{3}\right)$ and $\Delta\left({\bf 6},{\bf 1},-\frac{2}{3}\right)$ by $\Delta\left({\bf 6},{\bf 1},\frac{4}{3}\right)$. Neutrino NSI will then follow the $\rho$ NSI predictions as in Sec. 7.1.5. In this up-quark variant, one could replace the DQ $\Delta\left({\bf 6},{\bf 1},\frac{4}{3}\right)$ by a color triplet field $\Delta\left({\bf 3},{\bf 1},\frac{4}{3}\right)$ as well Klein:2019jgb .

7.2.3 Model with $SU(2)_{L}$-doublet and singlet leptoquarks

Operator $\mathcal{O}_{8}^{4}$ of Table. (V) does not induce neutrino mass via one-loop diagrams owing to the $SU(2)_{L}$ index structure. This operator will, however, lead to generation of neutrino masses at the two-loop level. A simple realization of $\mathcal{O}_{8}^{4}$ is given in Ref. Babu:2010vp . This model uses the same gauge symmetry and particle content as in the LQ variant of the Zee model (cf. Sec. 5), i.e., $\Omega\left({\bf 3},{\bf 2},\frac{1}{6}\right)=\left(\omega^{2/3},\omega^{-1/3}\right)$ and $\chi\left({\bf 3},{\bf 1},-\frac{1}{3}\right)$, with $\chi$ coupling modified as follows:

[TABLE]

Note that these Yukawa couplings conserve both baryon and lepton number as can be seen by assigning $(B,L)$ charges of $\left(\frac{1}{3},-1\right)$ to $\Omega$ and $\left(\frac{1}{3},1\right)$ to $\chi$. The couplings $\tilde{\lambda}_{\alpha\beta}u^{c}_{\alpha}d^{c}_{\beta}\chi^{\star}$, allowed by the gauge symmetry are forbidden by $B$, and the couplings $\lambda^{\prime}_{\alpha\beta}L_{\alpha}Q_{\beta}\chi^{\star}$ (as in Eq. (5.118)), allowed by gauge symmetry as well as $B$ are forbidden by $L$.191919The simultaneous presence of the $f$ and $\lambda^{\prime}$ couplings will drastically alter the successful $V-A$ structure of the SM Buchmuller:1986iq , and therefore, the $\lambda^{\prime}$ terms must be forbidden in this model by $L$ symmetry. The $L$ symmetry is softly broken by the cubic term in the scalar potential (5.119).

The simultaneous presence of Eqs. (7.223) and (5.119) would lead to neutrino mass generation at two-loop level as shown in Fig. 45. Here $X_{a}$ (with $a=1,2$) are the mass eigenstates obtained from the mixture of the $\omega^{-1/3}$ and $\chi^{-1/3}$ states (cf. Eq. (5.121)). Evaluation of the LQ-$W$ exchange diagrams in Fig. 45 (a) give the neutrino mass matrix as

[TABLE]

where 3 is a color factor, $\alpha$ is the $\omega-\chi$ mixing angle (cf. Eq. (5.122)), $M_{u,d,\ell}$ are diagonal mass matrices for the up- and down-type quarks, and charged leptons, respectively, $V$ is the CKM mixing matrix, $M\equiv{\rm min}(m_{1},m_{2})$ (with $m_{1,2}$ given by Eq. (5.123)), and $\mathcal{I}$ is a dimensionless two-loop integral that depends on $m_{1,2}$, $m_{W}$ and $M_{u,d,\ell}$ Babu:2010vp .

NSI induced in this LQ model has the same features as the LQ variant of the Zee model discussed in Sec. 5.4. Note that the $f_{\alpha\beta}$-couplings in Eq. (7.223) do not lead to neutrino NSI. The expression for the NSI parameters is given by Eq. (7.205), with the maximum allowed values given in Eq. (7.206) and also summarized in Table XX.

7.2.4 Leptoquark model with $SU(2)_{L}$-singlet vectorlike quark

This model utilizes the dimension-7 operator $L^{i}L^{j}\epsilon_{ij}Q^{k}H^{l}\epsilon_{kl}d^{c}$ to generate two-loop neutrino mass Babu:2011vb . This specific realization corresponds to the model $\mathcal{O}_{3}^{3}$ of Table III Cai:2014kra . In addition to the SM fields, an $SU(2)_{L}$-singlet vector-like quarks $U\left({\bf 3},{\bf 1},\frac{2}{3}\right)$ and $U^{c}\left({\bf 3}^{\star},{\bf 1},-\frac{2}{3}\right)$, and a scalar doublet LQ $\Omega\left({\bf 3},{\bf 2},\frac{1}{6}\right)=\left(\omega^{2/3},\omega^{-1/3}\right)$ are added to the SM spectrum. Addition of these fields leads to the following new Yukawa Lagrangian:

[TABLE]

where $\widetilde{\Omega}\equiv i\tau_{2}\Omega^{\star}$. The presence of all three Yukawa terms implies that lepton number is not conserved. Together with the quartic coupling term in the potential

[TABLE]

the Lagrangian (7.2.4) leads to neutrino mass generation at two-loop as shown in Fig. 46. This can be estimated as

[TABLE]

where $M_{d}$ and $M_{U}$ are the diagonal down quark and vectorlike quark mass matrices respectively, and $M=$max($m_{\omega},m_{U_{i}}$), with $m_{U_{i}}$ being the eigenvalues of $M_{U}$.

NSI in this model are induced by the $\omega^{-1/3}$ LQ and are given by Eq. (7.205). The maximum NSI that can be obtained in this model are given in Eq. (7.206) and are also summarized in Table XX.

7.2.5 Angelic model

This model induces operator $\mathcal{O}_{11}$ of Ref. Babu:2001ex :

[TABLE]

In this model Angel:2013hla , one adds two scalar LQs $\chi_{a}\left({\bf 3},{\bf 1},-\frac{1}{3}\right)$ (with $a=1,2$) and a color-octet Majorana fermion $F({\bf 8},{\bf 1},0)$. The relevant Yukawa Lagrangian is written as

[TABLE]

Expanding the first term, we get

[TABLE]

Within this framework, neutrino mass is induced at two-loop level as shown in Fig. 47 which can be estimated as

[TABLE]

where $V$ is the CKM-matrix, $M_{d}$ is the diagonal down-quark mass matrix, $M\equiv{\rm max}(m_{F},m_{\chi_{a}})$, and $\mathcal{I}$ is a loop function containing $m_{\chi_{a}},m_{F}$ and $M_{d}$ Angel:2013hla .

NSI in this model are induced by the singlet LQ $\chi$ and are given by

$\displaystyle\varepsilon_{\alpha\beta}\ =\ \frac{3}{4\sqrt{2}G_{F}}\frac{\lambda_{\alpha da}^{\star}\lambda_{\beta da}}{m_{\chi_{a}}^{2}}\,.$

(7.232)

This is similar to the singlet LQ contribution in Eq. (7.200). The maximum NSI in this model are the same as in Eq. (7.201). This is tabulated in Table XX.

7.2.6 Model with singlet scalar and vectorlike quark

This model realizes the ${\cal O}_{3}^{1}$ operator (cf. Table III) by adding a singlet scalar $\eta^{+}({\bf 1},{\bf 1},1)$ and vectorlike quark $\mathcal{Q}\left({\bf 3},{\bf 2},-\frac{5}{6}\right)=\left(D^{-1/3},X^{-4/3}\right)$. Neutrino mass is generated at two-loop level as shown in the Fig. 48. The relevant Lagrangian for the neutrino mass generation can be read as:

[TABLE]

The neutrino mass can be estimated as

[TABLE]

where $\sin\varphi$ represents the mixing between $W^{+}$ and $\eta^{+}$. The role of the vectorlike quarks in this model is to achieve such a mixing, which requires lepton number violation. Note that only the longitudinal component of $W$ mixes with $\eta^{+}$, which brings in two powers of lepton mass suppression in the neutrino mass estimate – one from the Yukawa coupling of the longitudinal $W$ and the other from a required chirality-flip inside the loop. It is to be noted that Eq. (7.234) does not fit the neutrino oscillation data as it has all diagonal entries zero, owing to the anti-symmetric nature of the $f$-couplings.

Other operators which lead to similar inconsistency with the neutrino oscillation data are $\mathcal{O}_{3}^{2},\>\mathcal{O}_{4}^{1}$ and $\mathcal{O}_{4}^{2}$ (cf. Tables III and IV). Therefore, we do not discuss the NSI prospects in these models.

7.2.7 Leptoquark model with vectorlike lepton

This model is a realization of ${\cal O}_{8}^{2}$ in Table V. This is achieved by adding an $SU(2)_{L}$-doublet LQ $\Omega\left({\bf 3},{\bf 2},\frac{1}{6}\right)$ and a vectorlike lepton $\psi\left({\bf 1},{\bf 2},-\frac{1}{2}\right)=(N,\,E)$. The Lagrangian responsible for neutrino mass generation can be written as

[TABLE]

Neutrino masses are generated at two-loop level via diagrams shown in Fig. 49 and can be estimated as:

[TABLE]

where $M_{d},\>M_{u},\>M_{\ell}$ and $M_{E}$ are the diagonal mass matrices for down quark, up quark, charged leptons and vectorlike leptons, respectively, and $m_{E}$ is the largest eigenvalue of $M_{E}$. The NSI parameters can be written as in Eq. (7.205), with the maximum values given in Eq. (7.206) and also summarized in Table XX.

7.2.8 Leptoquark model with $SU(2)_{L}$–doublet vectorlike quark

This model realizes the ${\cal O}_{8}^{3}$ operator (cf. Table V) by adding an $SU(2)_{L}$-doublet LQ $\Omega\left({\bf 3},{\bf 2},\frac{1}{6}\right)$ and an $SU(2)_{L}$-doublet vectorlike quark $\xi\left({\bf 3},{\bf 2},\frac{7}{6}\right)=\left(V^{5/3},U^{2/3}\right)$. The corresponding Lagrangian for the neutrino mass generation is given by

[TABLE]

Neutrino mass is generated at two-loop level as shown in Fig. 50 and can be estimated as

[TABLE]

where $M_{d},\>M_{u},\>M_{\ell}$ 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 in Eq. (7.205), with the maximum values given in Eq. (7.206).

7.2.9 A new two-loop leptoquark model

Here we propose a new two-loop LQ model for neutrino mass, where one can get NSI with up-quark. The effective $\Delta L=2$ operator is $d=13$, and is given by

[TABLE]

This model utilizes two scalar LQs – $\delta\left({\bf 3},{\bf 2},\frac{7}{6}\right)=\left(\delta^{5/3},\delta^{2/3}\right)$ and $\Omega\left({\bf 3},{\bf 2},\frac{\bf 1}{6}\right)=\left(\omega^{2/3},\omega^{-1/3}\right)$, and a scalar DQ $\hat{\Delta}\left({\bf 6}^{\star},{\bf 3},-\frac{1}{3}\right)=\left(\hat{\Delta}^{-4/3},\hat{\Delta}^{-1/3},\hat{\Delta}^{2/3}\right)$. The relevant Yukawa Lagrangian for the neutrino mass generation reads as

[TABLE]

The relevant terms in the potential that leads to neutrino mass generation read as

[TABLE]

The neutrino mass is induced at two-loop level as shown in Fig. 51 and can be estimated as

[TABLE]

where $M_{u}$ is the diagonal up-type quark mass matrix. Note that $M_{\nu}$ is a symmetric matrix, as it should be, since $h=h^{T}$.

After integrating out the heavy scalars, NSI induced in this model can be written as

$\displaystyle\varepsilon_{\alpha\beta}\ =\ \frac{3}{4\sqrt{2}G_{F}}\frac{f_{\alpha u}^{\star}f_{\beta u}}{m_{\delta}^{2}}\,.$

(7.243)

This is same as the extended one-loop LQ model prediction in Eq. (7.215) for $\lambda\gg\lambda^{\prime}$ with the exception that $\varepsilon_{\mu\mu}$ and $\varepsilon_{\tau\tau}$ are now constrained by IceCube. The maximum allowed values are given in Eq. (7.206). This is also summarized in Table XX.

7.3 Three-loop models

7.3.1 KNT Model

The Krauss-Nasri-Trodden (KNT) model Krauss:2002px generates the $d=9$ operator $\mathcal{O}_{9}$ of Eq. (1.4). SM-singlet fermions $N_{\alpha}({\bf 1},{\bf 1},0)$ and two SM-singlet scalars $\eta_{1}^{+}$ and $\eta_{2}^{+}$ with SM charges $({\bf 1},{\bf 1},1)$ are introduced. The relevant Yukawa Lagrangian is written as

[TABLE]

Tree level mass is prevented by imposing a $Z_{2}$ symmetry under which the fields $\eta_{2}^{+}$ and $N$ are odd, while the other fields are even. The Majorana mass term for $N$ as shown in Eq. (7.244) explicitly breaks lepton number. Neutrino masses are generated at three-loop as shown in Fig. 52 by the Lagrangian (7.244), together with the quartic term in the potential

[TABLE]

The estimated neutrino mass matrix reads as

[TABLE]

where $M_{\ell}$ is the diagonal charged lepton mass matrix, $M_{N}={\rm diag}(m_{N_{\alpha}})$ is the diagonal Majorana mass matrix for $N_{\alpha}$ fermions, $M\equiv{\rm max}(m_{N_{\alpha}},m_{\eta_{1}},m_{\eta_{2}})$, and ${\cal I}$ is a three-loop function obtained in general by numerical integration Cheung:2016ypw .

NSI in the KNT model arise from singly-charged scalar $\eta_{1}^{+}$ that has the same structure as in the Zee-Babu model (cf. Sec. 7.2.1) and are given by Eq. (7.195). The maximum NSI one can get in this model are same as in Eq. (7.196) and also summarized in Table XX.

7.3.2 AKS model

In the Aoki-Kanemura-Seto (AKS) model Aoki:2008av an effective $\Delta L=2$ operator of dimension 11 is induced:

[TABLE]

Note that there is a chiral suppression in this model unlike generic operators of type ${\cal O}_{1}^{\prime}$ given in Eq. (1.5). In addition to the SM fields, the following particles are added: an isospin doublet scalar $\Phi_{2}\left({\bf 1},{\bf 2},\frac{1}{2}\right)$, a singly-charged scalar singlet $\eta^{+}({\bf 1},{\bf 1},1)$, a real scalar singlet $\eta^{0}({\bf 1},{\bf 1},0)$, and two isospin-singlet right-handed neutrinos $N_{\alpha}({\bf 1},{\bf 1},0)$ (with $\alpha=1,2$). The relevant Yukawa Lagrangian for the neutrino mass generation reads as

[TABLE]

where $\Phi_{1}\left({\bf 1},{\bf 2},\frac{1}{2}\right)$ is the SM Higgs doublet. Tree-level neutrino mass is forbidden by imposing a $Z_{2}$ symmetry under which $\eta^{\pm}$, $\eta^{0}$ and $N_{\alpha 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.248) with the quartic term in the potential

[TABLE]

In Fig. 53 $H^{\pm}$ are the physical charged scalars from a linear combination of $\Phi_{1}$ and $\Phi_{2}$. The neutrino mass matrix reads as follows:

[TABLE]

where $\tan\beta\equiv\langle\Phi_{2}^{0}\rangle/\langle\Phi_{1}^{0}\rangle$ and $\mathcal{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

$\displaystyle\varepsilon_{\alpha\beta}=\frac{1}{4\sqrt{2}G_{F}}\frac{y_{e\alpha a}^{\star}y_{e\beta a}}{m_{H^{-}}^{2}}\,.$

(7.251)

This is similar to the heavy charged scalar contribution in Eq. (4.91). However, since the same Yukawa couplings $y_{e\alpha a}$ contribute to the electron mass in Eq. (7.248), we expect

[TABLE]

where $y_{e}$ is the electron Yukawa coupling in the SM. Thus, the maximum NSI in this model are of order of ${\cal O}\left(10^{-10}\right)$, as summarized in Table XX.

7.3.3 Cocktail Model

This model Gustafsson:2012vj induces operator $\mathcal{O}_{d=15}$ at the three-loop level:

[TABLE]

with $\Psi=L$ or $e^{c}$. The model includes two $SU(2)_{L}$-singlet scalars $\eta^{+}({\bf 1},{\bf 1},1)$ and $k^{++}({\bf 1},{\bf 1},2)$, and a second scalar doublet $\Phi_{2}\left({\bf 1},{\bf 2},\frac{1}{2}\right)$, in addition to the SM Higgs doublet $\Phi_{1}\left({\bf 1},{\bf 2},\frac{1}{2}\right)$. The fields $\eta^{+}$ and $\Phi_{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

[TABLE]

which breaks lepton number when combined with the following cubic and quartic terms in the potential:

[TABLE]

The $\Phi_{2}$ field is inert and does not get a VEV. After electroweak symmetry breaking, it can be written as

[TABLE]

For $\kappa_{1}\neq 0$, the singly-charged state $\phi_{2}^{+}$ mixes with $\eta^{+}$ (with mixing angle $\beta$), giving rise to two singly-charged scalar mass eigenstates:

[TABLE]

where $s_{\beta}\equiv\sin\beta$ and $c_{\beta}\equiv\cos\beta$.

The neutrino mass matrix is obtained from the three-loop diagram as shown in Fig. 54 and reads as Gustafsson:2012vj

[TABLE]

where $M_{\ell}$ stands for the diagonal charged lepton mass matrix.

As for the NSI, since both $\Phi_{2}$ and $\eta^{+}$ 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.

7.3.4 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 LQs. The effective operator induced in this model remains as $\mathcal{O}_{11}$ or Eq. (7.219). To achieve this, two isospin-singlet scalar LQs $\chi_{a}^{-1/3}\left({\bf 3},{\bf 1},-\frac{1}{3}\right)$ (with $a=1,2$) and at least two SM-singlet right-handed neutrinos $N_{\alpha}({\bf 1},{\bf 1},0)$ (with $\alpha=1,2$) are supplemented to the SM fields. A $Z_{2}$ symmetry is invoked under which $\chi_{2}^{-1/3}$ and $N$ are odd, while the rest of the fields are even. The relevant Yukawa Lagrangian is as follows:

[TABLE]

Here the first term expands to give $\lambda_{\alpha\beta}\left(\nu_{\alpha}d_{\beta}-\ell_{\alpha}u_{\beta}\right)\chi_{1}^{\star 1/3}$. These interactions, along with the quartic term in the potential

[TABLE]

generate neutrino masses at three-loop level, as shown in Fig. 55. The neutrino mass matrix reads as

[TABLE]

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 $\cal{I}$ is a dimensionless three-loop integral that depends on the ratio of the masses of particles inside the loop Nomura:2016ezz .

NSI in this model arise from the $\chi_{1}^{-1/3}$ interactions with neutrinos and down-quarks. The expression for NSI parameters is given as in Eq. (7.200), with the replacement $m_{\chi}\to m_{\chi_{1}}$. The maximum NSI for this model are the same as those given in Eq. (7.201) and are summarized in Table XX.

7.3.5 $SU(2)_{L}$–singlet three-loop model

This model Cepedello:2018rfh introduces two $SU(2)_{L}$-singlet scalars $\eta_{1}({\bf 1},{\bf 1},1)$ and $\eta_{2}({\bf 1},{\bf 1},3)$, and a singlet fermion $F({\bf 1},{\bf 1},2)$, in addition to the SM fields. The effective operator induced in this model is $\mathcal{O}_{9}$ in Eq. (1.4). The relevant Lagrangian term for the neutrino mass generation can be read as:

[TABLE]

With the potential term

[TABLE]

the Lagrangian (7.262) generates the neutrino mass at three-loop level, as shown in Fig. 56. The neutrino mass matrix can be written as

[TABLE]

where $M_{\ell}$ is the diagonal charged lepton mass matrix and $M\equiv\text{max}(m_{F},m_{\eta_{1}},m_{\eta_{2}})$. NSI in this model arise from singly-charged $\eta_{1}$ that has the same structure as in the Zee-Babu (cf. Sec. 7.2.1) and KNT (cf. Sec. 7.3.1) models and and are given by Eq. (7.195). The maximum NSI one can get in this model are same as in Eq. (7.196) and also summarized in Table XX. Other three-loop models of this type discussed in Ref. Cepedello:2018rfh will have similar NSI predictions.

7.4 Four- and higher-loop models

As noted in the introduction, it is very unlikely that neutrino masses and mixing of the right order can be induced in type-I radiative models at four or higher loops. The magnitude of $m_{\nu}$ in such models would be much smaller than needed to explain neutrino oscillation data, provided that the loop diagrams have chiral suppression proportional to a SM fermion mass. We illustrate below the difficulties with higher loop models with a four loop model presented in Ref. deGouvea:2019xzm .

In Ref. deGouvea:2019xzm an effective $d=9$ operator involving only $SU(2)_{L}$-singlet fermions of the SM was studied. The operator has the form

[TABLE]

Various UV completions are possible to induce this operator, with differing fermion contractions. All these models will induce light neutrino mass only at the four-loop level, since each fermion in ${\cal O}_{\rm s}$ has to be annihilated. A rough (and optimistic) estimate of the four-loop induced neutrino mass is deGouvea:2019xzm

[TABLE]

where $\Lambda$ is the UV cut-off scale. If the other Yukawa couplings involved are all of order one, $\Lambda=(100~{}{\rm MeV}-1~{}{\rm GeV})$ is needed to generate $m_{\nu}\sim 0.05$ eV. However, such a low value of $\Lambda$ will be inconsistent with experimental data on search for new particles, since the mediators needed to induce ${\cal O}_{\rm s}$ are either colored or electrically charged, with lower limits of order TeV on their masses from collider searches.

Models with such higher dimensional operators are nevertheless very interesting, as they can lead to lepton flavor and lepton number violating processes, without being constrained by neutrino masses, as emphasized in Ref. deGouvea:2019xzm . For example, neutrinoless double beta decay may occur at an observable level purely from ${\cal O}_{\rm s}$, which would be unrelated to the neutrino mass.

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_{\nu}$, and therefore, do not generally contribute to tree-level NSI, although small loop-level NSI effects are possible Bischer:2018zbd . To illustrate this point, let us take the scotogenic model Ma:2006km as a prototypical example. The new particles introduced in this model are SM-singlet fermions $N_{\alpha}({\bf 1},{\bf 1},0)$ (with $\alpha=1,2,3$) and an $SU(2)_{L}$ doublet scalar $\eta\left({\bf 1},{\bf 2},\frac{1}{2}\right):(\eta^{+},\eta^{0})$. A $Z_{2}$ symmetry is imposed under which the new fields $N_{\alpha}$ and $\eta$ are odd, while all the SM fields are even. The new Yukawa interactions in this model are given by

[TABLE]

Together with the scalar quartic term

[TABLE]

where $\Phi$ is the SM Higgs doublet, the Lagrangian (8.267) gives rise to neutrino mass at one-loop, as shown in Fig. 57. 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

[TABLE]

where we have assumed $M_{N}$ to be diagonal, and $m_{0}^{2}$ is the average squared mass of the real and imaginary parts of $\eta^{0}$. It is clear from Eq. (8.269) 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. 58, where the new particles added are as follows: one color-sextet DQ $\Delta\left({\bf 6},{\bf 1},\frac{4}{3}\right)$, one $SU(2)_{L}$ doublet scalar LQ $\delta\left({\bf 3},{\bf 2},\frac{7}{6}\right)=(\delta^{5/3},\delta^{2/3})$, and an $SU(2)_{L}$ singlet scalar LQ $\xi\left({\bf 3},{\bf 1},\frac{2}{3}\right)$. The relevant Yukawa Lagrangian is given by

[TABLE]

Together with the scalar potential terms

[TABLE]

where $\Phi$ is the SM Higgs doublet, the Lagrangian (8.270) gives rise to neutrino mass at two-loop level, as shown in Fig. 58. The neutrino mass can be approximated as follows:

[TABLE]

where $m_{1}$ and $m_{2}$ are the masses of the heaviest two LQs among the $\delta$, $\xi$ and $\Delta$ fields that run in the loop. Thus, although this model can be described as arising from an effective $\Delta L=2$ operator $\mathcal{O}^{\prime}_{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 $\delta$ 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. Liu:2016mpf .

9 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. 59 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 $\sim$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^{\pm}$ (cf. Sec. 4.8) and $\tau$ (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 ($\varepsilon_{ee},\,\varepsilon_{\mu\mu},\,\varepsilon_{\tau\tau})$ can be as large as $(8\%,\,3.8\%,\,9.3\%)$, while the off-diagonal NSI parameters $(\varepsilon_{e\mu},\,\varepsilon_{e\tau},\,\varepsilon_{\mu\tau})$ can be at most $(10^{-3}\%,\,0.56\%,\,0.34\%)$, as summarized in Fig. 59 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 ($\varepsilon_{ee},\,\varepsilon_{\mu\mu},\,\varepsilon_{\tau\tau})$ can be as large as $(0.4\%,\,21.6\%,\,34.3\%)$, while off-diagonal NSI $(\varepsilon_{e\mu},\,\varepsilon_{e\tau}\,\varepsilon_{\mu\tau})$ can only be as large as $(10^{-5}\%,\,0.36\%,\,0.43\%)$, as summarized in Fig. 59 and Table XVII. A variant of the LQ model with triplet LQs (cf. Sec. 6) allows for larger $\varepsilon_{\tau\tau})$ which can be as large as $51.7\%$. Neutrino scattering experiments are found to be the most constraining for the diagonal NSI parameters $\varepsilon_{ee}$ and $\varepsilon_{\mu\mu}$, while the cLFV searches are the most constraining for the off-diagonal NSI. $\varepsilon_{\tau\tau}$ 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.

Acknowledgments

We thank Sanjib Agarwalla, Sabya Chatterjee, Peter Denton, Radovan Dermišek, Arman Esmaili, Tao Han, Chris Kolda, Pedro Machado, Michele Maltoni, Ivan Martinez-Soler, Jordi Salvado, Yongchao Zhang and Yue Zhang for useful discussions. The work of KB, SJ, and AT was supported in part by the US Department of Energy Grant Number DE-SC 0016013. The work of BD was supported in part by the US Department of Energy under Grant No. DE-SC0017987 and by the MCSS. This work was also supported by the US Neutrino Theory Network Program under Grant No. DE-AC02-07CH11359. KB is supported in part by a Fermilab Distinguished Scholar program. We thank the Fermilab Theory Group for warm hospitality during the completion of this work. In addition, BD thanks the Department of Physics at Oklahoma State University for warm hospitality during the completion of this work. SJ and AT thank the Department of Physics at Washington University in St. Louis for warm hospitality, where part of this work was done.

[FIGURE:]

Appendix

Appendix A Analytic Expressions for Charged Higgs Cross Sections

It is instructive to write down the explicit formula for the charged-Higgs pair-production (Figs. 10 (a) and 10 (b) cross section:

[TABLE]

where $\beta=\sqrt{1-4m_{h^{+}}^{2}/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}\equiv\cos\theta_{w},~{}s_{w}\equiv\sin\theta_{w}$ ($\theta_{w}$ being the weak mixing angle). Note that the $t$-channel cross section depends on the Yukawa coupling $Y_{\alpha 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^{\pm}W^{\mp}$ (Fig. 10 (c)) is given by

[TABLE]

where $\theta$ is the angle made by the outgoing $h^{\pm}$ with respect to the initial $e^{-}$-beam direction, $\lambda(x,y,z)=x^{2}+y^{2}+z^{2}-2xy-2xz-2yz$, and

[TABLE]

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.