Full parameter scan of the Zee model: exploring Higgs lepton flavor violation

TL;DR

This paper performs a comprehensive numerical analysis of the Zee model, exploring its parameter space to understand neutrino masses, lepton flavor violation, and collider prospects, revealing testable predictions for upcoming experiments.

Contribution

It provides the first full parameter scan of the general Zee model, linking neutrino data, lepton flavor violation signals, and collider accessibility within a unified framework.

Findings

Normal neutrino mass ordering is favored over inverted.

Higgs lepton flavor violating decay BR can reach 10^{-2}.

Future experiments can test large parts of the model's parameter space.

Abstract

We study the general Zee model, which includes an extra Higgs scalar doublet and a new singly-charged scalar singlet. Neutrino masses are generated at one-loop level, and in order to describe leptonic mixing, both the Standard Model and the extra Higgs scalar doublets need to couple to leptons (in a type-III two-Higgs doublet model), which necessarily generates large lepton flavor violating signals, also in Higgs decays. Imposing all relevant phenomenological constraints and performing a full numerical scan of the parameter space, we find that both normal and inverted neutrino mass orderings can be fitted, although the latter is disfavored with respect to the former. In fact, inverted ordering can only be accommodated if $\theta_{23}$ turns out to be in the first octant. A branching ratio for $h \to \tau \mu$ of up to $10^{-2}$ is allowed, but it could be as low as $10^{-6}$. In addition, if future expected sensitivities of $\tau\to \mu\gamma$ are achieved, normal ordering can be almost completely tested. Also, $\mu e$ conversion is expected to probe large parts of the parameter space, excluding completely inverted ordering if no signal is observed. Furthermore, non-standard neutrino interactions are found to be smaller than $10^{-6}$, which is well below future experimental sensitivity. Finally, the results of our scan indicate that the masses of the additional scalars have to be below $2.5$ TeV, and typically they are lower than that and therefore within the reach of the LHC and future colliders.

Figures

Tables

Table 1. Table 1: Experimental 95 % percent 95 95~{}\% C.L. upper bounds on HLFV decays from ATLAS and CMS in the tau sector using the 13 TeV data sets. In our numerical scan, we will use the strongest upper bounds from CMS shown in this table.

HLFV observable	ATLAS	CMS
$\mathrm{Br}\left(h\to \tau \mu \right)$	$1.43\%$  Aad:2016blu	$1.20\%$ CMS:2016qvi
$\mathrm{Br}\left(h\to \tau e\right)$	$1.04\%$ Aad:2016blu	$0.69\%$ Khachatryan:2016rke

Table 2. Table 2: Priors on the 19 free real parameters used in the scan. For Y 2 τ ​ τ superscript subscript 𝑌 2 𝜏 𝜏 Y_{2}^{\tau\tau} , Y 2 τ ​ μ superscript subscript 𝑌 2 𝜏 𝜇 Y_{2}^{\tau\mu} , Y 2 τ ​ e superscript subscript 𝑌 2 𝜏 𝑒 Y_{2}^{\tau e} , and Y 2 μ ​ τ superscript subscript 𝑌 2 𝜇 𝜏 Y_{2}^{\mu\tau} , we scan real and imaginary parts independently, while f μ ​ τ superscript 𝑓 𝜇 𝜏 f^{\mu\tau} , f e ​ τ superscript 𝑓 𝑒 𝜏 f^{e\tau} , and Y 2 e ​ τ superscript subscript 𝑌 2 𝑒 𝜏 Y_{2}^{e\tau} can be assumed to be real without loss of generality. We use logarithmic priors for all parameters except for tan ⁡ β 𝛽 \tan\beta , where uniform priors are used.

Parameter	Range
Complex: $Y_2^{\tau \tau }$, $Y_2^{\tau \mu }$, $Y_2^{\tau e}$, $Y_2^{\mu \tau }$	$[{10}^{-12},{10}^{-1}]$
Real: $f^{\mu \tau }$, $f^{e\tau }$, $Y_2^{e\tau }$	$[{10}^{-12},{10}^{-1}]$
$\tan \beta$	$[0.3,50]$
${\lambda }_1$, ${\lambda }_2$, $\left|{\lambda }_3\right|$, $\left|{\lambda }_5\right|$	$[{10}^{-5},\sqrt{4\pi }]$
${\mu }_h$, ${\mu }_2$ [GeV]	$[1,{10}^7]$
$\mu$ [GeV]	$[1,{10}^7]$

Table 3. Table 3: The current experimental values and 1 ​ σ 1 𝜎 1\sigma errors for electroweak precision tests (see appendix  A ) and universality (see section  3.4.1 ). The limits on the parameters S 𝑆 S , T 𝑇 T , and U 𝑈 U are derived from a fit to different electroweak precision data, see ref.  Baak:2014ora for more details. The correlation coefficients among the obervables S 𝑆 S , T 𝑇 T , and U 𝑈 U are ρ ST = 0.90 subscript 𝜌 ST 0.90 \rho_{\rm ST}=0.90 , ρ SU = − 0.59 subscript 𝜌 SU 0.59 \rho_{\rm SU}=-0.59 , and ρ TU = − 0.83 subscript 𝜌 TU 0.83 \rho_{\rm TU}=-0.83   Baak:2014ora .

Observable	Central value $\pm$ $1\sigma$ error
$S$	$0.05\pm 0.11$ Baak:2014ora
$T$	$0.09\pm 0.13$ Baak:2014ora
$U$	$0.01\pm 0.11$ Baak:2014ora
$\left|g_{\tau }^{\exp }/g_{\mu }^{\exp }\right|$	$1.0011\pm 0.0015$ Pich:2013lsa
$\left|g_{\tau }^{\exp }/g_e^{\exp }\right|$	$1.0030\pm 0.0015$ Pich:2013lsa
$\left|g_{\mu }^{\exp }/g_e^{\exp }\right|$	$1.0018\pm 0.0014$ Pich:2013lsa

Table 4. Table 4: The best-fit values and 1 ​ σ 1 𝜎 1\sigma errors for the leptonic mixing parameters and mass-squared differences from the NuFIT group (NuFIT 3.0, November 2016)  Gonzalez-Garcia:2014bfa ; Esteban:2016qun ; NuFIT3.0 . See also the discussion in section  2.3 . Note that we use symmetric lower and upper errors (choosing the largest of the two when different asymmetric errors are present). Here Δ ​ m 3 ​ ℓ 2 = Δ ​ m 31 2 > 0 Δ subscript superscript 𝑚 2 3 ℓ Δ subscript superscript 𝑚 2 31 0 \Delta m^{2}_{3\ell}=\Delta m^{2}_{31}>0 for NO and Δ ​ m 3 ​ ℓ 2 = Δ ​ m 32 2 < 0 Δ subscript superscript 𝑚 2 3 ℓ Δ subscript superscript 𝑚 2 32 0 \Delta m^{2}_{3\ell}=\Delta m^{2}_{32}<0 for IO. Note that for sin 2 ⁡ θ 23 superscript 2 subscript 𝜃 23 \sin^{2}\theta_{23} there are two minima in their distribution  Esteban:2016qun ; NuFIT3.0 for both orderings, corresponding to the first and second octants. For δ 𝛿 \delta the distribution is also not χ 2 superscript 𝜒 2 \chi^{2} distributed. Therefore, for sin 2 ⁡ θ 23 superscript 2 subscript 𝜃 23 \sin^{2}\theta_{23} and δ 𝛿 \delta , the two-dimensional complete distribution is used.

Observable	NO	IO
${\sin }^2{\theta }_{12}$	$0.306\pm 0.012$	$0.306\pm 0.012$
${\sin }^2{\theta }_{13}$	$0.02166\pm 0.00075$	$0.02179\pm 0.00076$
${\sin }^2{\theta }_{23}$	$0.441\pm 0.027$	$0.587\pm 0.024$
$\Delta m_{21}^2$ $\left[{\mathrm{eV}}^2\right]$	$\left(7.50\pm 0.19\right)\cdot {10}^{-5}$	$\left(7.50\pm 0.19\right)\cdot {10}^{-5}$
$\Delta m_{3\ell }^2$ $\left[{\mathrm{eV}}^2\right]$	$\left(2.524\pm 0.040\right)\cdot {10}^{-3}$	$-\left(2.514\pm 0.041\right)\cdot {10}^{-3}$
$\delta$ $[{}_{}{}^{\circ }]$	$261\pm 59$	$277\pm 46$

Table 5. Table 5: The current experimental 90 % C.L. upper bounds for relevant charged lepton flavor violating processes and electric and magnetic dipole moments.

Observable	Upper bound
$\mathrm{Br}\left({\tau }^{-}\to {\mu }^{-}\gamma \right)$	$4.4\cdot {10}^{-8}$ Patrignani:2016xmw
$\mathrm{Br}\left({\tau }^{-}\to e^{-}\gamma \right)$	$3.3\cdot {10}^{-8}$ Patrignani:2016xmw
$\mathrm{Br}\left({\mu }^{+}\to e^{+}\gamma \right)$	$4.2\cdot {10}^{-13}$ TheMEG:2016wtm
$\mathrm{Br}\left({\tau }^{-}\to {\mu }^{-}{\mu }^{+}{\mu }^{-}\right)$	$2.1\cdot {10}^{-8}$ Patrignani:2016xmw
$\mathrm{Br}\left({\tau }^{-}\to {\mu }^{-}{\mu }^{+}{e}^{-}\right)$	$2.7\cdot {10}^{-8}$ Patrignani:2016xmw
$\mathrm{Cr}{\left(\mu \to e\right)}_{\mathrm{Au}}$	$7\cdot {10}^{-13}$ Bertl:2006up
$\left|\delta a_e\right|$	$2\cdot {10}^{-12}$ Patrignani:2016xmw
$\left|d_{\mu }/e\right|$ [cm]	$1\cdot {10}^{-19}$ Patrignani:2016xmw
$\left|d_e/e\right|$ [cm]	$8.7\cdot {10}^{-29}$ Baron:2013eja

Table 6. Table 6: The current experimental 95 % C.L. upper bounds on parameters related to the neutrino masses, see also section  2.3 . The two upper bounds on m e ​ e subscript 𝑚 𝑒 𝑒 m_{ee} are due to the sensitivity on the nuclear matrix elements, whereas for m ν e subscript 𝑚 subscript 𝜈 𝑒 m_{\nu_{e}} , we show results coming from two different experiments. We use the most stringent values of these two upper bounds in the scan.

Observable	Upper bound
$m_{ee}$ $\left[\mathrm{meV}\right]$	$[190,450]$ Albert:2014awa
$m_{{\nu }_e}$ $\left[\mathrm{eV}\right]$	$[2.05,2.3]$ Lobashev:2003kt ; Kraus:2004zw ; Aseev:2011dq
$\sum m_i$ $\left[\mathrm{eV}\right]$	$0.23$ Ade:2015xua

Table 7. Table 7: Some results of our numerical scan. We show the ranges of the 95 % C.L. regions of different observables and parameters for NO and IO for two values κ = 1 , 10 𝜅 1 10 \kappa=1,10 for the naturality upper bound on the trilinear coupling μ 𝜇 \mu , see eq. ( 46 ).

	NO		IO
Quantity	$\kappa =1$	$\kappa =10$	$\kappa =1$	$\kappa =10$
${\chi }_{\min }^2$	10.7	11.0	21.7	21.5
$\mathrm{Br}\left(h\to \tau \mu \right)$	$[1\cdot {10}^{-6},1\cdot {10}^{-2}]$	$[1\cdot {10}^{-6},1\cdot {10}^{-2}]$	$[2\cdot {10}^{-7},4\cdot {10}^{-3}]$	$[1\cdot {10}^{-7},5\cdot {10}^{-3}]$
$\mathrm{Br}\left(h\to \tau e\right)$	$[1\cdot {10}^{-10},2\cdot {10}^{-4}]$	$[1\cdot {10}^{-10},2\cdot {10}^{-4}]$	$[6\cdot {10}^{-9},3\cdot {10}^{-4}]$	$[3\cdot {10}^{-9},3\cdot {10}^{-4}]$
$\mathrm{Br}\left(\tau \to \mu \gamma \right)$	$[8\cdot {10}^{-10},3\cdot {10}^{-8}]$	$[1\cdot {10}^{-10},3\cdot {10}^{-8}]$	$[3\cdot {10}^{-11},3\cdot {10}^{-8}]$	$[3\cdot {10}^{-11},3\cdot {10}^{-8}]$
$\mathrm{Br}\left(\mu \to e\gamma \right)$	$[1\cdot {10}^{-21},6\cdot {10}^{-13}]$	$[3\cdot {10}^{-22},6\cdot {10}^{-13}]$	$[1\cdot {10}^{-31},1\cdot {10}^{-12}]$	$[1\cdot {10}^{-34},1\cdot {10}^{-12}]$
$\mathrm{Cr}\left(\mu \to e\right)$	$[1\cdot {10}^{-21},4\cdot {10}^{-13}]$	$[1\cdot {10}^{-21},4\cdot {10}^{-13}]$	$[3\cdot {10}^{-17},3\cdot {10}^{-13}]$	$[3\cdot {10}^{-17},3\cdot {10}^{-13}]$
$m_A,m_H$ $\left[\mathrm{TeV}\right]$
$m_{h_1^{+}}$ $\left[\mathrm{TeV}\right]$
$\sin \left(\beta -\alpha \right)$	$[0.98,1.0]$	$[0.98,1.0]$	$[0.97,1.0]$	$[0.97,1.0]$

Equations

$$\left(\begin{array}[]{c}H_{1}\\ H_{2}\end{array}\right)=\left(\begin{array}[]{cc}c_{\beta}&s_{\beta}\\ -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^{+}\\ \dfrac{1}{\sqrt{2}}\left(v+\varphi_{1}^{0}+iG^{0}\right)\end{array}\right)\,,\qquad H_{2}=\left(\begin{array}[]{c}H^{+}\\ \dfrac{1}{\sqrt{2}}\left(\varphi_{2}^{0}+iA\right)\end{array}\right)\,,$$

$$V$$

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

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

$$+{\mu }_h^2{\left|h^{+}\right|}^2+{\lambda }_h{\left|h^{+}\right|}^4+{\lambda }_8{\left|h^{+}\right|}^2{H}_1^{†}{H}_1+{\lambda }_9{\left|h^{+}\right|}^2{H}_2^{†}{H}_2$$

$$+{\lambda }_{10}{\left|h^{+}\right|}^2\left(H_1^{†}{H}_2+H.c.\right)+\left(\mu {ϵ}_{\alpha \beta }{H}_1^{\alpha }{H}_2^{\beta }{h}^{-}+H.c.\right),$$

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

$$\mathcal{M}_{c}^{2}=\left(\begin{array}[]{cc}M_{H^{+}}^{2}&-\mu v/\sqrt{2}\\ -\mu v/\sqrt{2}&M_{33}^{2}\end{array}\right),$$

$$M_{H^{+}}^2={\mu }_2^2+\frac{1}{2}{v}^2{\lambda }_3,M_{33}^2={\mu }_h^2+v^2{\lambda }_8.$$

$$\left(\begin{array}{c}{h}_1^{+}\\ h_2^{+}\end{array}\right)=\left(\begin{array}{cc}{s}_{\phi }& c_{\phi }\\ c_{\phi }& -s_{\phi }\end{array}\right)\left(\begin{array}{c}{h}^{+}\\ H^{+}\end{array}\right),$$

$$s_{2\phi }=\frac{\sqrt{2}v\mu }{m_{h_2^{+}}^2-m_{h_1^{+}}^2}$$

$$m_{h_1^{+},h_2^{+}}^2\equiv \frac{1}{2}\left[M_{H^{+}}^2+M_{33}^2\mp \sqrt{{\left(M_{H^{+}}^2-M_{33}^2\right)}^2+2v^2{\mu }^2}\right].$$

$$\mathcal{M}_{h}^{2}=\left(\begin{array}[]{cc}\lambda_{1}v^{2}&\lambda_{6}v^{2}\\ \lambda_{6}v^{2}&\,m_{A}^{2}+\lambda_{5}v^{2}\end{array}\right)\,,$$

$$m_A^2=M_{H^{+}}^2-\frac{1}{2}{v}^2\left({\lambda }_5-{\lambda }_4\right).$$

$$\left(\begin{array}[]{c}h\\ H\end{array}\right)=\left(\begin{array}[]{cc}s_{\beta-\alpha}&c_{\beta-\alpha}\\ c_{\beta-\alpha}&-s_{\beta-\alpha}\\ \end{array}\right)\left(\begin{array}[]{c}\varphi_{1}^{0}\\ \varphi_{2}^{0}\end{array}\right)$$

$$m_{H,h}^2\equiv \frac{1}{2}\left\{m_A^2+v^2\left({\lambda }_1+{\lambda }_5\right)\pm \sqrt{{\left[m_A^2+v^2\left({\lambda }_5-{\lambda }_1\right)\right]}^2+4v^4{\lambda }_6^2}\right\},$$

$$s_{2(\beta -\alpha )}=-\frac{2{\lambda }_6{v}^2}{m_H^2-m_h^2}$$

$$-{\mathcal{L}}_L=\overline{L}(Y_1^{†}{\Phi }_1+Y_2^{†}{\Phi }_2)e_{\mathrm{R}}+\overline{\tilde{L}}fL{h}^{+}+H.c.,$$

$$m_E=\frac{v}{\sqrt{2}}\left(c_{\beta }{Y}_1^{†}+s_{\beta }{Y}_2^{†}\right).$$

$$-{\mathcal{L}}_L=\overline{L}\left[\frac{\sqrt{2}{m}_E}{v}{H}_1+\left(\frac{Y_2^{†}}{c_{\beta }}-\frac{\sqrt{2}{m}_E{t}_{\beta }}{v}\right)H_2\right]e_{\mathrm{R}}+\overline{\tilde{L}}fL{h}^{+}+H.c.$$

$$-{\mathcal{L}}_L$$

$$+\overline{e_{\mathrm{L}}}\left(\frac{-m_E{s}_{\alpha }}{v{c}_{\beta }}+c_{\beta -\alpha }\frac{Y_2^{†}}{\sqrt{2}{c}_{\beta }}\right)e_{\mathrm{R}}h+\overline{e_{\mathrm{L}}}\left(\frac{m_E{c}_{\alpha }}{v{c}_{\beta }}-s_{\beta -\alpha }\frac{Y_2^{†}}{\sqrt{2}{c}_{\beta }}\right)e_{\mathrm{R}}H$$

$$+i\overline{e_{\mathrm{L}}}\left(-\frac{m_E{t}_{\beta }}{v}+\frac{Y_2^{†}}{\sqrt{2}{c}_{\beta }}\right)e_{\mathrm{R}}A+H.c.$$

$$g_{h^0}^1$$

$$g_{h^0}^2$$

$$g_{h_c}^1$$

$$g_{h_c}^2$$

$$g_{h_c}^3$$

$${\mathcal{M}}_{\nu }=U{D}_{\nu }{U}^T,$$

$${\nu }_{\alpha }=\sum _{i=1}^3{U}_{\alpha i}{\nu }_i.$$

$$\displaystyle U=\left(\begin{array}[]{ccc}c_{13}c_{12}&c_{13}s_{12}&s_{13}e^{-i\delta}\\ -c_{23}s_{12}-s_{23}s_{13}c_{12}e^{i\delta}&c_{23}c_{12}-s_{23}s_{13}s_{12}e^{i\delta}&s_{23}c_{13}\\ s_{23}s_{12}-c_{23}s_{13}c_{12}e^{i\delta}&-s_{23}c_{12}-c_{23}s_{13}s_{12}e^{i\delta}&c_{23}c_{13}\end{array}\right)\left(\begin{array}[]{ccc}1&0&0\\ 0&e^{i\phi_{1}/2}&0\\ 0&0&e^{i\phi_{2}/2}\end{array}\right)\,,$$

$$m_{ee}$$

$$m_{{\nu }_e}$$

$$\sum m_i$$

$${\cal M}_{\nu}=A\,\Big{[}f\,m_{E}^{2}+m_{E}^{2}f^{T}-\frac{v}{\sqrt{2}\,s_{\beta}}(f\,m_{E}\,Y_{2}+Y_{2}^{T}\,m_{E}\,f^{T})\Big{]}\,,$$

Full text

††preprint: ADP-17-1/T1007

1 Introduction

In the Standard Model (SM), neutrinos are massless and lepton flavors are exactly conserved to all orders. However, from neutrino oscillation experiments, we know that neutrinos are not massless and that lepton flavor is not conserved in the neutrino sector.

Whether lepton number is a good symmetry of Nature or not remains an open question. If the SM is considered an effective field theory (EFT), the only dimension-5 operator is the Weinberg operator Weinberg:1979sa , where lepton number is broken by two units, and which gives rise to Majorana masses for the neutrinos after electroweak symmetry breaking. There are different realizations of this operator, both at tree level and one-loop level. High-scale realizations of the Weinberg operator, for instance the type-I seesaw mechanism Minkowski:1977sc ; Yanagida:1980 ; Gell-Mann:1980vs ; Glashow:1979vf ; Mohapatra:1979ia , are well motivated by grand unified theories (GUTs), such as $\rm SO(10)$. However, these models are difficult to test, and therefore, other avenues should be explored, in particular in light of new data from the LHC and low-energy experiments.

Indeed, there is hope to test scenarios in which lepton number violation (LNV) occurs close to the electroweak scale. One such scenario is radiative neutrino mass models, where neutrinos are massless at tree level, but acquire mass at one or more loops. Thus, the new degrees of freedom involved in the generation of neutrino mass cannot be too heavy, and therefore, they can be searched for at the LHC. These new particles typically give rise to enhanced lepton flavor violating (LFV) signals in processes like $\mu\to e\gamma$, $\tau\to 3\mu$, or $\mu e$ conversion, which we denote as charged lepton flavor violation (CLFV). Furthermore, with the discovery of the Higgs boson Aad:2012tfa ; Chatrchyan:2012xdj , the ways to search for LFV have increased and one can look for Higgs lepton flavor violating (HLFV) decays, especially in the $\tau$-$\mu$ and $\tau$-$e$ sectors, which are subject to weaker constraints from low-energy probes than the $\mu$-$e$ sector.

The ATLAS and CMS experiments at the LHC have active programs to search for HLFV decays. Considering their 13 TeV data sets, no signal has been observed.111The CMS 8 TeV data showed a $2.4\sigma$ excess in the channel $h\to\tau\mu$ Khachatryan:2015kon , which is translated into a branching ratio ${\rm Br}(h\to\tau\mu)=(0.84^{+0.39}_{-0.37})~{}\%$. This corresponds to an upper limit $\mathrm{Br}(h\to\tau\mu)<1.51~{}\%$ at $95~{}\%$ C.L. Khachatryan:2015kon . Unfortunately, this small excess has disappeared with the CMS 13 TeV data. Of course, an excess at a lower level of $\mathcal{O}(10^{-3})$ could be observed at the LHC in the forthcoming years. In table 1, we summarize the latest 13 TeV upper limits on HLFV decays from ATLAS and CMS. The LHC is sensitive to ${\rm Br}(h\to\tau\mu,\,\tau e)\gtrsim 0.001$, and therefore, these decays can be used to test the models of new physics with HLFV signals in such a range.

In recent years, several studies have analyzed if a $\mathrm{Br}(h\to\tau\mu)\sim 1~{}\%$ is compatible with low-energy constraints, either using an EFT approach Goudelis:2011un ; DiazCruz:1999xe ; Blankenburg:2012ex ; Harnik:2012pb ; Celis:2013xja ; Banerjee:2016foh ; Belusca-Maito:2016axk ; Herrero-Garcia:2016uab or in a type-III two-Higgs doublet model (2HDM) Davidson:2010xv ; Sierra:2014nqa ; Dorsner:2015mja ; Omura:2015nja ; Botella:2015hoa ; Guo:2016ixx (see also ref. Arganda:2015naa for a study of a supersymmetric inverse seesaw scenario). Reference Herrero-Garcia:2016uab is particularly relevant to us, where it was shown that the only tree-level scenarios, which can accommodate the excess, are models with extra scalars. Furthermore, the connection between HLFV decays and neutrino masses was extensively discussed222See also refs. Arganda:2004bz ; Arana-Catania:2013xma ; Arganda:2014dta ; Arganda:2015naa ; Arganda:2015uca ; Arganda:2016qvs for HLFV studies in supersymmetric and seesaw scenarios. and it was found that the most general version of the Zee model Zee:1980ai ; Cheng:1980qt ; Wolfenstein:1980sy (see also refs. Petcov:1982en ; Zee:1985id ; Bertolini:1987kz ; Bertolini:1990vz ; Smirnov:1994wj ; Smirnov:1996hc ; Frampton:1999yn ; Jarlskog:1998uf ; Ghosal:2001ep ; Kanemura:2000bq ; Balaji:2001ex ; Koide:2002uk ; Brahmachari:2001rn ; Frampton:2001eu ; Assamagan:2002kf ; He:2003ih ; Kanemura:2005hr ; AristizabalSierra:2006ri ; He:2011hs ; Babu:2013pma for different works and variations on the Zee model) was the most promising one. In this simple model, an extra Higgs scalar doublet and a new singly-charged scalar singlet are added to the SM and Majorana neutrino masses are generated at one-loop level. In order to describe leptonic mixing correctly, it is necessary that both scalar doublets couple to the charged leptons in a type-III 2HDM, see e.g. refs. Frampton:2001eu ; Koide:2002uk ; He:2003ih . This is precisely the same requirement needed to have large HLFV Herrero-Garcia:2016uab , and therefore, a complete analysis, taking into account all phenomenological constraints and performing a full numerical scan of the parameter space, is of great interest. This is the aim of this work, including a reduction of the allowed parameter space of the model by taking into account recent data, like the discovery of the Higgs boson Aad:2012tfa ; Chatrchyan:2012xdj , the determination of the leptonic mixing angle $\theta_{13}$ Abe:2011fz ; An:2012eh ; Ahn:2012nd , the results from the latest global fits of neutrino parameters including the hint of leptonic CP violation and the uncertainty on the octant of the leptonic mixing angle $\theta_{23}$ (i.e. if $\theta_{23}$ is smaller or larger than $\pi/4$) Gonzalez-Garcia:2014bfa ; Esteban:2016qun ; NuFIT3.0 , the new limits on CLFV processes like the ones on $\mu\to e\gamma$ TheMEG:2016wtm , and HLFV processes like those shown in table 1. Finally, the impact of future expected limits will also be studied, in particular those coming from $\tau\to\mu\gamma$, where Belle II is expected to reach a sensitivity of $10^{-9}$ HerediadelaCruz:2016yqi , and specially $\mu e$ conversion, which is expected to improve by several orders of magnitude in the near future, see e.g. refs. Carey:2008zz ; Kutschke:2011ux ; Donghia:2016lzt ; ProjectX ; Barlow:2011zza ; Witte:2012zza ; Cui:2009zz ; Kuno:2013mha .

The paper is structured as follows. In section 2, we describe the Zee model and its relevant parameters. In section 3, we discuss the phenomenological constraints of the model. Then, in section 4, we perform a numerical scan and present our results for three different scenarios: (i) without neutrino masses (just a type-III 2HDM with an extra charged singlet) and with neutrino masses for both (ii) normal and (iii) inverted neutrino mass orderings. Finally, in section 5, we summarize our results and give our conclusions. In addition, in appendices A and B, we present the contributions of the model to the electroweak precision test parameters $S$, $T$, and $U$ and derive explicit analytical expressions for various loop functions that these parameters are constructed from and which can be used for any model.

2 The general Zee model

In addition to the SM content with a Higgs scalar doublet $\Phi_{1}$, the Zee model Zee:1980ai ; Cheng:1980qt ; Wolfenstein:1980sy contains an extra Higgs scalar doublet $\Phi_{2}$ and a singly-charged scalar singlet $h^{+}$. We start by discussing the most general scalar potential.

2.1 The scalar potential

The following analysis is similar to the ones performed for 2HDMs, see e.g. refs. Gunion:2002zf ; Davidson:2005cw ; Branco:2011iw ; Djouadi:2005gj . One can start in a generic basis, where both Higgs scalar doublets $\Phi_{1}$ and $\Phi_{2}$ take VEVs denoted by $v_{1}$ and $v_{2}$, respectively. Then, one can rotate to the Higgs basis, where only $H_{1}$ takes a VEV denoted as usual as $v=\sqrt{v^{2}_{1}+v^{2}_{2}}\simeq 246$ GeV. The rotation is given by the following transformation Davidson:2005cw :

[TABLE]

where $\tan\beta\equiv v_{2}/v_{1}$ and the short-hand notations $s_{x}\equiv\sin x$ and $c_{x}\equiv\cos x$.333In type-III 2HDM, $\tan\beta$ is an unphysical parameter Davidson:2005cw . For the lepton sector, it can be defined as the ratio of the tau Yukawa coupling (times the vev) and its mass. In general, the definition of $\tan\beta$ will be different for up and down quarks. In our numerical scan (see section 4.1), we will treat it as an arbitrary free parameter. We will also use $t_{x}\equiv\tan x$. In the Higgs basis, the doublets take the form:

[TABLE]

where $\varphi_{1}^{0}$ and $\varphi_{2}^{0}$ are CP-even neutral Higgs fields, $A$ is a CP-odd neutral Higgs field, $H^{+}$ is a charged Higgs field, and $G^{+}$ and $G^{0}$ are the would-be Goldstone bosons, which are eaten by the $W^{+}$ and the $Z$. The most general potential for the Zee model (see e.g. ref. Davidson:2005cw ) is given in the Higgs basis by

[TABLE]

where $\lambda_{i}$ ($i=1,2,\ldots,10,h$) are the quartic couplings, $\mu^{2}_{i}$ ($i=1,2,3,h$) are bare mass-squared parameters, and $\mu$ is a trilinear coupling. In addition, $\epsilon_{\alpha\beta}$ is the rank two antisymmetric Levi-Civita tensor. In general, $\lambda_{5}$, $\lambda_{6}$, $\lambda_{7}$, $\lambda_{10}$, $\mu_{3}$, and $\mu$ can be complex. Note that one can choose $\lambda_{5}$ to be real by redefining $H_{1}$ and $H_{2}$ Davidson:2005cw . Furthermore, without loss of generality, $\mu$ can be chosen to be real and positive by redefining the singlet $h^{-}$. In addition, we choose $\lambda_{6}$ to be real for simplicity. In the numerical scan (see section 4.1), we will treat the three quantities $\mu$, $\mu_{2}$, and $\mu_{h}$ as free real parameters, except for the case when we will set $\mu=0$. In section 3.1, we will comment on the usage of the quartic couplings $\lambda_{i}$ in the numerical scan.

Since only $H_{1}$ takes a VEV, differentiating eq. (2.1) with respect to $H_{1}$ and $H_{2}$, gives the following minimization conditions

[TABLE]

which can be used to eliminate $\mu_{1}^{2}$ and $\mu_{3}^{2}$ as independent variables. Equation (4) applies to both the real and imaginary parts. Inserting $\left\langle H_{1}\right\rangle=(0,v/\sqrt{2})^{T}$ into eq. (2.1), we obtain the squared mass matrices of the charged and neutral CP-even Higgs states. For the charged ones, in the Higgs basis $\left(H^{+},h^{+}\right)$, we have

[TABLE]

where

[TABLE]

The mass eigenstates $h_{1}^{+}$ and $h_{2}^{+}$ are a mixing of $h^{+}$ and $H^{+}$ given by

[TABLE]

where

[TABLE]

and the masses are defined as

[TABLE]

Similarly, the CP-even mass matrix is

[TABLE]

where the mass of the CP-odd Higgs state enters as

[TABLE]

Thus, in the Higgs basis, the mass eigenstates $h$ and $H$ are a mixture of the CP-even states $\varphi_{1}$ and $\varphi_{2}$

[TABLE]

with the masses defined as

[TABLE]

where the CP-even mixing is given by

[TABLE]

that needs to be sufficiently close to zero (i.e. the alignment limit) to give rise to a SM-like Higgs boson Gunion:2002zf .

2.2 The lepton sector

As we will see, in order to describe leptonic mixing, both Higgs scalar doublets must couple to the charged leptons, and thus, we are considering a type-III 2HDM, see e.g. ref. He:2003ih . The most general Yukawa Lagrangian in the generic basis, where both Higgs fields take VEVs, reads

[TABLE]

where $L=(\nu_{\rm L},\,e_{\rm L})^{T}$ and $e_{\rm R}$ are the SU(2) lepton doublets and singlets, respectively, and $\tilde{L}\equiv i\sigma_{2}L^{c}=i\sigma_{2}C\overline{L}^{T}$ with $\sigma_{2}$ being the second Pauli matrix. Due to Fermi statistics, $f$ is an antisymmetric Yukawa matrix in flavor space (i.e. $f^{\alpha\beta}=-f^{\beta\alpha}$), while $Y_{1}$ and $Y_{2}$ are completely general complex Yukawa matrices. Furthermore, the charged-lepton masses are given by

[TABLE]

Note that we will work in the basis where $m_{E}$ is diagonal with real and positive elements $m_{e}$, $m_{\mu}$, and $m_{\tau}$. Moreover, $Y_{2}$ will be a general complex matrix and $Y_{1}$ can be expressed completely in terms of $m_{E}$ and $Y_{2}$ using eq. (16).

In the Higgs basis, we can rewrite eq. (15) using eq. (1) as

[TABLE]

Without loss of generality, rotating the lepton doublets $L^{\alpha}$ and the lepton singlets $e^{\alpha}_{\rm R}$ by the same phase (so that $m_{E}$ remains diagonal and positive), three phases from $f$ can be removed. However, note that the phases from $Y_{2}$ cannot be removed by lepton field redefinitions.

In the mass basis (also for massive neutrinos), using eq. (17), the most general leptonic Lagrangian reads

[TABLE]

We define the following effective couplings for the neutral Higgs fields $h^{0}=(h,\;H,\;A)$, which will turn out to be useful for CLFV processes:

[TABLE]

and for the charged Higgs fields $h_{c}=(h_{1}^{+},\,h_{2}^{+})$:

[TABLE]

One can observe that $g_{h^{0}}^{1}$ is flavor conserving and proportional to $m_{E}/v$.

2.3 Neutrino parameters

A general $3\times 3$ Majorana neutrino mass matrix $\mathcal{M}_{\nu}$, which is defined as an effective mass term in the Lagrangian $\mathcal{L}_{\nu}\equiv-1/2\,\overline{\nu_{\rm L}^{c}}\mathcal{M}_{\nu}\nu_{\rm L}+\mathrm{H.c.}$, can be written as

[TABLE]

where $\nu_{\rm L}$ is the left-handed neutrino flavor eigenfield with three lepton flavors, $D_{\nu}$ is a $3\times 3$ diagonal matrix with positive real eigenvalues and $U$ is the $3\times 3$ unitary leptonic mixing matrix, which relates the neutrino mass eigenfields $\nu_{i}$ ($i=1,2,3$) with definite masses $m_{i}$ and the neutrino flavor eigenfields $\nu_{\alpha}$ ($\alpha=e,\mu,\tau$):

[TABLE]

The standard parametrization for $U$ is Patrignani:2016xmw

[TABLE]

where $c_{ij}\equiv\cos\theta_{ij}$ and $s_{ij}\equiv\sin\theta_{ij}$ ($\theta_{12}$, $\theta_{13}$, and $\theta_{23}$ being the three leptonic mixing angles and $\theta_{12},\theta_{13},\theta_{23}\in[0,\pi/2)$). Furthermore, in eq. (32), $\delta$ is the Dirac CP-violating phase ($\delta\in[0,2\pi)$) and $\phi_{1}$ and $\phi_{2}$ are two Majorana CP-violating phases ($\phi_{1},\phi_{2}\in[0,4\pi)$). Using the three definite neutrino masses $m_{1}$, $m_{2}$, and $m_{3}$, we also define the two linearly-independent neutrino mass squared-differences $\Delta m_{21}^{2}\equiv m_{2}^{2}-m_{1}^{2}$ and $\Delta m_{31}^{2}\equiv m_{3}^{2}-m_{1}^{2}$, known as the small and large mass squared-differences, respectively, where the sign of $\Delta m_{31}^{2}$ is still unknown. The case $\Delta m_{31}^{2}>0$ is generally referred to as ‘normal neutrino mass ordering’ (NO), whereas the case $\Delta m_{31}^{2}<0$ is known as the ‘inverted neutrino mass ordering’ (IO). Note that using neutrino oscillation experiments, it is not possible to determine $\phi_{1}$ and $\phi_{2}$ nor the absolute neutrino mass scale. The most up-to-date best-fit values from global analyses of the ordinary neutrino oscillation parameters (i.e. the leptonic mixing parameters and the neutrino mass-squared differences) are $\theta_{12}\simeq 34^{\circ}$, $\theta_{13}\simeq 8.5^{\circ}$, $\theta_{23}\simeq 42^{\circ}$ for NO and $\theta_{23}\simeq 50^{\circ}$ for IO, $\delta\simeq 1.5\pi$, $\Delta m_{21}^{2}\simeq 7.5\cdot 10^{-5}\,{\rm eV}^{2}$, and $\Delta m_{31}^{2}\simeq 2.5\cdot 10^{-3}\,{\rm eV}^{2}$ for NO and $\Delta m_{31}^{2}\simeq-2.5\cdot 10^{-3}\,{\rm eV}^{2}$ for IO Gonzalez-Garcia:2014bfa ; Esteban:2016qun ; NuFIT3.0 . Similar values can also be found in refs. Capozzi:2016rtj ; Forero:2014bxa . Note that both the first and second octants of $\theta_{23}$ are allowed from the global analyses Esteban:2016qun ; NuFIT3.0 for both orderings with a mild preference of the first (second) octant for NO (IO).

In addition to the ordinary neutrino oscillation parameters, the following effective neutrino parameters appear naturally in different contexts Patrignani:2016xmw

[TABLE]

where $m_{ee}$ (or $m_{0\nu 2\beta}$) is the effective electron neutrino mass parameter that could be measured in neutrinoless double beta decay ($0\nu 2\beta$) experiments Wolfenstein:1981rk ; Bilenky:1987ty (see also ref. DellOro:2016tmg for a recent review), $m_{\nu_{e}}$ (or $m_{\beta}$) is the effective neutrino mass parameter measured in (single) beta decay experiments McKellar:1980cn , and finally, $\sum m_{i}$ is the sum of the three neutrino masses, which, in the future, could be determined by cosmology, but at present it is only restricted by an upper bound, see e.g. refs. Hannestad:2010yi ; GonzalezGarcia:2010un .

2.4 Neutrino masses in the Zee model

As can be seen from the potential and the Yukawa Lagrangian of the Zee model, eqs. (2.1) and (2.2), respectively, in order to have LNV and therefore neutrino masses, we need the simultaneous presence of $Y_{1}$, $Y_{2}$, $f$, and $\mu$. In the Zee model, the one-loop diagram shown in figure 1, where the charged scalars $h^{+}_{1}$ and $h^{+}_{2}$ run in the loop, generates neutrino masses. The complete neutrino mass matrix is then given by (see e.g. ref. He:2011hs )

[TABLE]

where we have defined

[TABLE]

with $\varphi$ being the mixing angle for the charged scalars given in eq. (8). Therefore, in the Zee model, due to loop and chiral suppression, the new physics scale can be light.

Assuming $f^{e\mu}=0$, neglecting $m_{e}\ll m_{\mu},\,m_{\tau}$, and keeping only the term proportional to $m_{\mu}$ in the 3-3 element,444Keeping the 3-3 element to order $m_{\mu}/m_{\tau}$ is phenomenologically relevant for the following two reasons. First, to have that all neutrinos are massive. Second, to obtain a constrain on $Y_{2}^{\mu\tau}$, which enters in $h\to\tau\mu$. In our numerical analysis, we keep all terms proportional to $m_{\mu}$. we obtain the following (symmetric) Majorana mass matrix

[TABLE]

Note that in a simpler scenario, where one neglects terms proportional to $m_{\mu}$, one neutrino will be massless. On the other hand, taking terms proportional to $m_{\mu}$ into account, all neutrinos will obtain masses.

In our analysis, we assume zero Yukawa couplings in the $e$-$\mu$ sector, i.e. we assume $f^{e\mu}=0$ and $Y_{2}^{\mu\mu}=Y_{2}^{\mu e}=Y_{2}^{e\mu}=Y_{2}^{ee}=0$, which means that the non-zero Yukawa couplings are $f^{e\tau}$, $f^{\mu\tau}$, $Y_{2}^{\tau\tau}$, $Y_{2}^{\tau\mu}$, $Y_{2}^{\tau e}$, $Y_{2}^{\mu\tau}$, and $Y_{2}^{e\tau}$. We assume all these Yukawa couplings to be complex except for $f^{e\tau},\,f^{\mu\tau}$, and $Y_{2}^{e\tau}$ (which does not enter in neutrino masses), see the discussion in section 2.2. Thus, the Yukawa couplings will constitute eleven free real parameters in the numerical scan that will be described in section 4.1.

In order to obtain correct mixing angles, we need both $Y_{2}^{\tau\mu}$ and $Y_{2}^{\tau e}$ different from zero, as they enter in the 1-2 submatrix of eq. (38). Therefore, it is clear that reproducing the leptonic mixing angles correctly will imply restrictions on $\mathrm{Br}(h\to\tau\mu)$, $\mathrm{Br}(h\to\tau e)$, and other LFV processes. In fact, from this argumentation, a lower bound on the product $\mathrm{Br}(h\to\tau\mu)\cdot\mathrm{Br}(h\to\tau e)$ (in addition to an upper bound from other CLFV processes) is expected.

2.5 The (minimal) quark sector

Although the Zee model only deals with the lepton sector, the SM Higgs scalar doublet needs to couple to the SM quarks, like tops and bottoms, in order to be observed via its production and decay modes at the LHC Khachatryan:2016vau . In the generic basis, the most general Lagrangian in the quark sector is given by

[TABLE]

where $Q=(u_{\rm L},\,d_{\rm L})^{T}$ are the SU(2) quark doublets, $u_{\rm R}$ and $d_{\rm R}$ are the SU(2) quark singlets, and $\tilde{\Phi}_{i}\equiv i\sigma_{2}\Phi_{i}^{*}$ ($i=1,2$). However, flavor violation in the quark sector is severely constrained (see e.g. ref. Harnik:2012pb ), so we will assume the simplest scenario in which $Y_{d2}=Y_{u2}=0$. Then, we can use the basis, where the up-type quark mass matrix is diagonal. Furthermore, we assume the Yukawa couplings $Y_{d1}$ and $Y_{u1}$ to be Hermitian. The masses for the third generation quarks are given by $m_{b}=(Y^{33}_{d1})^{*}c_{\beta}v/\sqrt{2}$ and $m_{t}=(Y^{33}_{u1})^{*}c_{\beta}v/\sqrt{2}$. Therefore, the interactions of the physical neutral Higgs fields with quarks are given by

[TABLE]

where the corresponding Feynman rule for the CP-odd scalar $A$ includes a $\gamma_{5}$.

Note that if we had taken the couplings to quarks as general as possible, including the first generation, there would have been other phenomenological implications. In particular, related to neutrino masses, there would have been be new contributions to neutrinoless double beta decay and new universality and non-standard neutrino interactions with matter, stemming from interactions of the charged scalars $h^{+}_{1}$ and $h^{+}_{2}$. However, when naturality constraints are imposed on the Yukawa couplings to the leptons and the up and down quarks, see eq. (44), these contributions are subdominant. We will therefore only discuss the universality constraints and the non-standard neutrino interactions generated through leptonic interactions, see section 3.4, and we will only consider the contributions to neutrinoless double beta decay mediated by $W$ bosons, i.e. the contributions from the light neutrinos.

3 Phenomenology

3.1 Stability of the potential

A Hamiltonian in quantum mechanics has to be bounded from below, which requires the quartic part of the scalar potential in eq. (2.1) to be positive for all values of the fields and for all scales. Then, if two of the three fields $H_{1},\,H_{2}$, and $h$ vanish, one immediately finds

[TABLE]

For a general 2HDM potential with $\lambda_{6}=\lambda_{7}=0$, it has been shown in ref. Ivanov:2008er that the additional necessary conditions are

[TABLE]

However, when $\lambda_{6},\lambda_{7}\neq 0$, it has been shown that in addition to the previous conditions, the following condition Ferreira:2009jb

[TABLE]

is both necessary and sufficient to ensure stability of the potential. Other stability conditions for similar potentials are discussed in refs. Kanemura:2000bq ; Kannike:2012pe ; Herrero-Garcia:2014hfa .

In this work, due to the large number of parameters, we will set $\lambda_{4}=\lambda_{7}=\lambda_{8}=\lambda_{9}=\lambda_{10}=\lambda_{h}=0$, since they do not significantly impact phenomenology, even though their presence is expected to somewhat open the allowed parameter space. Thus, the four free Higgs couplings are $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$, and $\lambda_{5}$, which we will treat as free real parameters, while $\lambda_{6}$ is a derived parameter that can be computed from eq. (13). In the numerical scan (see section 4.1), we impose the conditions from eqs. (41)–(43).

3.2 Naturality and perturbativity

There are naturality and perturbativity constraints on the Yukawa couplings and on the quartic and trilinear couplings of the potential. In order not to have large fine-tuned cancellations between the different Yukawa couplings (see e.g. refs. Dorsner:2015mja ; Herrero-Garcia:2016uab ), we demand that

[TABLE]

One can also obtain an upper bound on $\mu$, which contributes to the scalar masses. In fact, using eq. (9), it is clear that naturality demands that $\mu\lesssim\sqrt{2}m^{2}_{h^{+}_{2}}/v$. Interestingly, we can also derive a natural upper bound using the $125$ GeV Higgs boson Aad:2012tfa ; Chatrchyan:2012xdj , due to the fact that $\mu$ contributes at one-loop level to its mass. The relevant coupling of the light Higgs boson to the charged scalars induced by $\mu$ is

[TABLE]

We demand that the one-loop contribution to the Higgs mass fulfills $\delta m_{h}/m_{h}\lesssim\kappa$, where we choose $\kappa=1\,(10)$, which corresponds to no (10 %) fine-tuning.555The fine-tunings in the Higgs mass squared, which is the relevant parameter in the Lagrangian, would be 1 % (100 %) for $\kappa=1\,(10)$. Neglecting logarithms and factors of two in the Higgs self-energies, we obtain

[TABLE]

Taking $s_{\beta-\alpha}\sim 1$, we find an upper bound of $1.5\,(15)$ TeV for $\kappa=1\,(10)$. In addition, we impose that all the quartic couplings are perturbative:

[TABLE]

3.3 Charged lepton flavor violation and electric and magnetic moments

3.3.1 Trilepton decays

The presence of the second Higgs doublet gives rise to tree-level trilepton decays $\ell_{i}\to\ell_{j}\overline{\ell_{k}}\ell_{l}$.666At one loop and two loops, there are dipole contributions which dominate the rate. However, these are strongly bounded by $\tau\to\mu\gamma$ and $\mu\to e\gamma$, see section 3.3.2. Also, box diagrams are very suppressed, see ref. Mitsuda:2001vh . The ratio of branching ratio reads

[TABLE]

where $G_{F}=g^{2}/(\sqrt{2}v^{2})\simeq 1.166\cdot 10^{-5}\,{\rm GeV}^{-2}$ is the Fermi coupling constant and the Wilson coefficients $D_{\rm PP^{\prime}}$ (${\rm P,P^{\prime}=L,R}$) are given by the coherent sum of the contributions from the neutral Higgs fields. In addition, $\xi=1/2\,(1)$ when there are two (no) indistinguishable particles in the final state. We are interested in tau decays.777Other processes, like $\mu\to e\overline{e}e$, are absent at tree level, since we assume $Y_{2}^{\mu e}=Y_{2}^{e\mu}=0$. Also, tree level contributions to $\tau\to\mu e\bar{e}$ and $\tau\to ee\bar{e}$ are suppressed by $m_{e}$. In this case, for $\tau\to\mu\mu\bar{\mu}$, the Wilson coefficient $D_{\rm LL}$ is given by

[TABLE]

where $h^{0}=(h,\,H,\,A)$. Similarly, for $\tau\to e\mu\bar{\mu}$, one can simply substitute $Y_{2}^{\mu\tau}\to Y_{2}^{e\tau}$. Furthermore, $D_{\rm RL}$ is obtained from $D_{\rm LL}$ by changing $(g^{1}_{h^{0}})^{*}\to g^{1}_{h^{0}}$. Finally, $D_{\rm RR}$ ($D_{\rm LR}$) is obtained from $D_{\rm LL}$ ($D_{\rm RL}$) by making the replacement $Y_{2}^{\alpha\beta}\to(Y_{2}^{\beta\alpha})^{*}$ and conjugating the vertex factors.

As expected, these tree-level processes do not restrict the parameter space as much as $\tau\to\mu\gamma$, $\tau\to e\gamma$, or $\mu\to e\gamma$ does, since they always involve a muon mass suppression (squared) and are therefore irrelevant. In table 5, the upper bounds used in the numerical scan for the various observables are presented.

3.3.2 $\ell_{i}\to\ell_{j}\gamma$ decays

One of the most constrained CLFV process is the radiative process $\ell_{i}\to\ell_{j}\gamma$ with $\ell_{i}$ being the physical charged leptons $e$, $\mu$, and $\tau$. This process always arises at loop level and it can be viewed as stemming from an effective operator of the form (Davidson:2010xv, )

[TABLE]

where $\Lambda$ is the scale of new physics, $\sigma^{\mu\nu}=i[\gamma^{\mu},\gamma^{\nu}]/2$, and $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is the electromagnetic field strength tensor with $A_{\mu}$ being the photon field. It is useful to define

[TABLE]

and similarly, $C_{\rm L}=C^{\dagger}_{\rm R}$. Then, it follows that

[TABLE]

We use the expressions for the Wilson coefficients given in refs. Goudelis:2011un ; Blankenburg:2012ex ; Harnik:2012pb , adapted to the Zee model. We also include the two-loop Barr–Zee contributions as given in ref. Chang:1993kw . At one-loop level, the dominant contribution for $\tau\to\mu\gamma$ reads

[TABLE]

and similarly, we obtain $C_{\rm R}^{0}$ with the replacement $Y_{2}^{\alpha\beta}\to(Y_{2}^{\beta\alpha})^{*}$. For the charged scalars, we find (using $\sum_{i}|U_{\tau i}|^{2}\approx 1$)

[TABLE]

where the contribution to $C_{\rm R}^{+}$ is zero, since we assume $f^{e\mu}=0$.

The total contribution to the Wilson coefficient is $C_{\rm L}=C_{\rm L}^{0}+C_{\rm L}^{+}+C_{\rm L}^{\mbox{\scriptsize 2-loops}}$, and similarly for $C_{\rm R}$. For $\mu\to e\gamma$, the dominant contributions, proportional to $m_{\tau}$, are given by the neutral Higgs fields and $\tau$ running in the loop

[TABLE]

Note that we also add the contributions from $f$ to $\mu\to e\gamma$, which are proportional to $\,m_{\mu}$, with $\nu_{\tau}$ running in the loop:

[TABLE]

The contribution in eq. (59) strongly constraints the antisymmetric Yukawa coupling $f$ of the singly-charged scalar singlet.

3.3.3 Electron and muon electric dipole and anomalous magnetic moments

When flavor is conserved, anomalous magnetic moments (AMMs) are generated. The dominant contributions are given by loops with neutral Higgs fields and tau leptons, similar to those in eqs. (56) and (57). It can be defined as Blankenburg:2012ex

[TABLE]

Then, using eq. (51), the muon AMM is given by

[TABLE]

Similarly, the electron AMM is obtained by replacing the indices $\mu\to e$ everywhere.

For the electron AMM, there is no disagreement between theory and experiment, and in fact, it represents one of the most precisely measured quantities in all of physics with an experimental 90 % C.L. upper bound on possible new physics contributions of $2\cdot 10^{-12}$ Patrignani:2016xmw . On the other hand, for the muon AMM, there is an experimental deviation from the theoretical prediction of the SM, i.e. $\left|\delta a_{\mu}\right|=(2.88\pm 0.80)\cdot 10^{-9}$ Patrignani:2016xmw . Unfortunately, neither a 2HDM nor the Zee model can accommodate this discrepancy Dorsner:2015mja .

If there is CP violation, electric dipole moments (EDMs) are generated, which can be defined as

[TABLE]

which gives

[TABLE]

The experimental 90 % C.L. upper bounds for the muon and electron EDMs are $1\cdot 10^{-19}\,e$ cm Patrignani:2016xmw and $8.7\cdot 10^{-29}\,e$ cm Baron:2013eja , respectively (see also table 5). Therefore, for the muon EDM, this bound is almost non-existing, whereas the electron one strongly constrains the imaginary parts of $Y_{2}^{\tau e}$ and $Y_{2}^{e\tau}$ (in our scenario, $Y_{2}^{e\tau}$ is assumed to be real for simplicity).

3.3.4 $\mu e$ conversion in nuclei

Finally, $\mu e$ conversion is a very interesting process because the sensitivity in the next generation experiments is expected to increase by around four orders of magnitude. As the dipole contribution is already heavily constrained by $\mu\to e\gamma$, we can use the current limits to restrict the monopole photonic contribution. In the limit, where the transferred momentum is zero, the conversion rate, relative to the muon capture rate, can be expressed in gold as Kitano:2002mt

[TABLE]

where $V^{p}_{\rm Au}=0.0974\,m_{\mu}^{5/2}$, $\Gamma_{\rm Au}=8.6\cdot 10^{-18}\,{\rm GeV}$ Kitano:2002mt , and the vector current coefficient is given by Dorsner:2015mja

[TABLE]

Similarly, $g_{\rm RV}$ is obtained by replacing $Y_{2}^{\alpha\beta}\to(Y_{2}^{\beta\alpha})^{*}$ in eq. (65). The currently best experimental limit is $\mathrm{Cr}(\mu\to e)_{\rm Au}<7\cdot 10^{-13}$ at 90 % C.L. by the SINDRUM II experiment at PSI Bertl:2006up .

For titanium, the best experimental limit is $\mathrm{Cr}(\mu\to e)_{\rm Ti}<4.3\cdot 10^{-12}$ at 90 % C.L. Bertl:2006up . In the future, PRISM/PRIME Barlow:2011zza ; Witte:2012zza expects to achieve a sensitivity of $\mathcal{O}(10^{-18})$. For aluminium, the future experimental sensitivity is $\mathrm{Cr}(\mu\to e)_{\rm Al}<6\cdot 10^{-17}$ at 90 % C.L. in the Mu2e experiment at Fermilab Carey:2008zz ; Kutschke:2011ux ; Donghia:2016lzt and of $\mathcal{O}(10^{-17})$ in the COMET experiment at J-PARC Cui:2009zz ; Kuno:2013mha . Mu2e may also run at Project X using either Al or Ti ProjectX with an expected sensitivity of $\mathcal{O}(10^{-19})$. In the numerical analysis, we will discuss how the Zee model will be constrained if the expected sensitivity of this process is achieved in the future.

3.4 Leptonic interactions of the charged scalars

Now, we study four-lepton interactions of the singly-charged scalar mass eigenstates $h_{1}^{+}$ and $h_{2}^{+}$ that involve two charged leptons and two neutrinos.888Since the new Higgs scalar doublet also couples to quarks (to the third generation in our scenario), there are four-fermion interactions between quarks and leptons, which we do not analyze any further, see the discussion at the end of section 2.5. These give rise to muon and tau decays into lighter charged leptons and neutrinos as well as to non-standard neutrino interactions (NSIs), see e.g. refs. Bilenky:1993bt ; Davidson:2003ha .

At tree level, we can integrate out $h_{1}^{+}$ and $h_{2}^{+}$, see e.g. ref. Bilenky:1993bt . Therefore, using eq. (2.2) and expanding to first order in $p^{2}/m^{2}_{h_{1,2}^{+}}$, we obtain

[TABLE]

where we have defined the effective Yukawa couplings

[TABLE]

and the effective masses

[TABLE]

Note that for the definition of $\tilde{M}_{12}^{2}$ we have used the charged-scalars mixing $s_{2\varphi}$ as defined in eq. (8).

3.4.1 Universality

The second operator in eq. (3.4), which is second order in $f$, i.e. $|\mathcal{Y}_{2}^{\rm eff}|^{2}\propto f^{\dagger}f$, couples to the left-handed leptons, like charged currents in the SM. This implies that it interferes constructively with the $W$ boson. In the SM, the Fermi constant extracted from muon decay $G_{\mu}^{\rm SM}$ and the one extracted from hadronic decays $G_{\beta}^{\rm SM}$ are tested to be equal with great precision. The presence of the charged scalars modifies the muon decay rate Santamaria1987 ; Nebot:2007bc , which implies that $G_{\beta}^{\rm SM}=G_{\mu}^{\rm SM}\neq G_{\mu}^{\rm Zee}$, where $G_{\mu}^{\rm Zee}$ is the Fermi constant from muon decay in the Zee model. Therefore, we have

[TABLE]

where $\tilde{M}_{2}$ is defined in eq. (68).

In the SM, unitarity of the quark mixing matrix $V$ holds to great precision. In our scenario, as we assume $f^{e\mu}=0$, we also have that $V$ is unitary up to order $\mathcal{O}(1/m^{4}_{h_{1,2}^{+}})$:

[TABLE]

On the hand, other leptonic decays may not be universal (in the SM, they are mediated by gauge interactions and are therefore universal). The ratio of flavor violating decays can be tested among the different generations via the effective couplings given by

[TABLE]

In our scenario the expressions (71)–(73) will generally, but not necessarily, be different from one, i.e. they deviate from the SM prediction.

3.4.2 Non-standard neutrino interactions

Apart from standard neutrino interactions (including neutrino oscillations), the new singly-charged scalar fields $h_{1}^{+}$ and $h_{2}^{+}$ introduced in the Zee model will induce NSIs at tree level. These NSIs are new LFV processes that are not allowed in the SM, but could be probed in future neutrino oscillation experiments, and are usually treated using an effective four-fermion operator.

Interestingly, the operators in the second line of eq. (3.4) violate lepton number.999Their gauge invariant EFT operators are, of course, dimension 7 deGouvea:2014lva , $LLL\overline{e_{\rm R}}\Phi_{1,2}$. Indeed, they involve the same combination of four leptons that appears inside the neutrino mass diagram, see figure 1, and thus, their coefficients are proportional to the same lepton-number combination appearing in the neutrino mass formula, see eq. (38). Hence, they are subject to constraints from neutrino masses and therefore suppressed.

The operators in the first line of eq. (3.4) do not violate lepton number and, in principle, they give rise to NSIs that are not suppressed by neutrino masses. Applying Fierz identities one can express them in various ways. Using ref. Nieves:2003in , they can be written in the flavor basis with the usual NSI language as Antusch:2008tz ; Ohlsson:2009vk

[TABLE]

where $\chi^{\rho\sigma}_{\alpha\beta}$ and $\varepsilon^{\rho\sigma}_{\alpha\beta}$ are the canonical NSI parameters given by

[TABLE]

where $\tilde{M}_{1}$ and $\tilde{M}_{2}$ are defined in eq. (68).

For neutrinos propagating in ordinary matter, these operators induce the following matter NSI parameters

[TABLE]

In our scenario, the only relevant matter NSI parameters are $\chi^{\rm m}_{\tau\tau}$ and $\varepsilon^{\rm m}_{\tau\tau}$ (cf. ref. Ohlsson:2009vk ), since we assume $Y_{2}^{\mu e}=Y_{2}^{e\mu}=0$ and $f^{e\mu}=0$. Now, we can derive an upper bound on $\chi^{\rm m}_{\tau\tau}$ applying CLFV and HLFV limits. For illustration, let us impose the limits from ${\rm Br}(h\rightarrow\tau e)$ from table 1. Using a similar equation to eq. (88), but for the $\tau e$ channel, which depends on $Y_{2}^{e\tau}$, we obtain

[TABLE]

which is below present and most probably future experimental sensitivity. Reproducing small neutrino masses and fulfilling other stronger CLFV constraints, e.g. constraints on $\mu\to e\gamma$ and $\tau\to\mu\gamma$, imply that $\chi^{\rm m}_{\tau\tau}$ and $\varepsilon^{\rm m}_{\tau\tau}$ are much smaller than $10^{-8}$, which is beyond any future experimental sensitivity. This can be seen in our numerical scan. Other limits and future prospects on NSIs can be found in refs. Biggio:2009nt ; Blennow:2016etl .

Finally, at a future neutrino factory, there could be source NSIs in the process $\mu\to e\overline{\nu_{\beta}}\nu_{\alpha}$. In our scenario, the only relevant source NSI parameters $\chi^{\mu e}_{\alpha\beta}$ and $\varepsilon^{\mu e}_{\alpha\beta}$ are those with tau neutrinos, i.e. $\chi^{\mu e}_{\tau\tau}$ and $\varepsilon^{\mu e}_{\tau\tau}$. However, these are also very small, at least below $10^{-6}$.

3.5 Higgs signals

In the Zee model, the couplings to SM particles are modified with respect to their SM values. For instance, the couplings to gauge bosons are

[TABLE]

and similarly, for $g_{hZZ}$, $g_{HZZ}$, and $g_{AZZ}$, changing $m^{2}_{W}\to m^{2}_{Z}$. Clearly, close to the decoupling limit, $s_{\beta-\alpha}\to 1$, the light Higgs interactions are sufficiently SM-like Gunion:2002zf . As we will see in section 3.5.2, in order to have HLFV, we cannot be exactly at the decoupling limit, but it needs to be close enough to fulfill the bounds on the Higgs decays measured at the LHC Khachatryan:2016vau . The other Higgs couplings in the Zee model are modified as in eq. (2.2) (see also eq. (88)) for leptons, eq. (40) for quarks, and eq. (82) for photons.

The Higgs results at the LHC are usually given in terms of the global signal strength defined as

[TABLE]

where $\sigma_{X}(h)$ is the cross section for the production mode $X$ and ${\rm Br}(h\to Y)$ is the Higgs branching ratio for the decay mode $Y$. By definition, in the SM, $\mu_{if}^{\rm SM}=1$ for all production modes $i$ and decay channels $f$. At the LHC, there are four production modes available for the Higgs boson, where the dominant one is gluon-gluon fusion (ggF), mainly through a top loop. The subdominant ones are vector boson fusion (VBF), associated production with a vector boson $Vh$ (where $V=W,Z$), and the associated production with a top-quark pair $t\bar{t}h$. The production modes are usually grouped into two effective modes according to ${\rm ggF}+t\bar{t}h$ and ${\rm VBF}+Vh$. We consider the five decay channels, where a signal has been detected, namely $\gamma\gamma$, $WW^{*}$, $ZZ^{*}$, $b\bar{b}$, and $\tau\bar{\tau}$. For instance, for $gg\to h\to b\overline{b}$, the signal strength is given by

[TABLE]

For the contribution to the $\chi^{2}$ function (to be discussed in section 4) from the Higgs decay channels, we need to take into account correlations between different production modes. Thus, for each of the decay modes $f=\gamma\gamma,WW^{*},ZZ^{*},b\bar{b},\tau\bar{\tau}$, the contribution to the $\chi^{2}$ function is defined as

[TABLE]

where $\mu^{f}_{1,2}$ are the results in the Zee model, $\hat{\mu}^{f}_{1(2)}$ are the measured Higgs signal strengths, $\hat{\sigma}^{f}_{1(2)}$ are the standard deviations, and $\rho$ is the correlation. The index 1 stands for the combination ${\rm ggF}+t\bar{t}h$ and the index 2 for the combination ${\rm VBF}+Vh$. The numerical values are given in refs. Aad:2013wqa ; ATLAS-CONF-2014-009 ; ATLAS-CONF-2013-079 ; Aad:2015vsa .

3.5.1 $h\to\gamma\gamma$ decays

In the Zee model, the decay of the Higgs boson to two photons is modified by two factors. First, the couplings to gauge bosons and top quarks are changed, since we have two Higgs doublets. Second, there are new extra charged scalars couplings to the Higgs boson. In ref. Kanemura:2000bq , a study of $h\to\gamma\gamma$ in the Zee model has been performed. However, $\lambda_{7}$ and $\lambda_{10}$ were set to zero. Thus, in the following, we will analyze this decay in our scenario.

The value of the $h\to\gamma\gamma$ decay width in the Zee model with respect to the SM one is given by Ellis:1975ap ; Shifman:1979eb ; Carena:2012xa

[TABLE]

where $\lambda_{{\cal S}H}$ is the coupling of a charged scalar ${\cal S}$ with mass $m_{\cal S}$ to the Higgs field, which is coming from a term in the potential of the form $(H^{\dagger}H)({\cal S}^{\dagger}{\cal S})$. Note that we have used the modified couplings to tops and $WW$ given in eqs. (40) and (78), respectively. Here, $\tau_{i}\equiv 4m_{i}^{2}/m_{H}^{2}$ and $A_{i}(x)$ ($i=0,1/2,1$) are loop functions:

[TABLE]

We need to compute the couplings to charged scalars $\lambda_{{\cal S}H}$ in the Zee model. Since we are in the Higgs basis, the terms of the potential in eq. (2.1), which are coupling the Higgs boson to the charged scalars, are those involving $\lambda_{3}$, $\lambda_{7}$, $\lambda_{8}$, $\lambda_{10}$, and $\mu$. Using the rotations to the mass basis, i.e. eqs. (7) and (12), the relevant interactions $hh_{1}^{+}h_{1}^{-}$ and $hh_{2}^{+}h_{2}^{-}$ read

[TABLE]

In the case when there is no mixing $\mu\,(\varphi)\to 0$, we obtain

[TABLE]

where the first equation agrees with eq. (F1) in ref. Gunion:2002zf . Note that, in our scenario, only the terms proportional to $\lambda_{3}$ and $\mu$ in eq. (86) will contribute, as we set the other couplings to zero.

3.5.2 Higgs lepton flavor violation

The Zee model predicts HLFV interactions that can be sizable. From the leptonic Lagrangian, i.e. eq. (2.2), the branching ratio of $h\to\tau\mu$ is given by

[TABLE]

and similarly, $\mathrm{Br}(h\to\tau e)\propto(|Y_{2}^{\tau e}|^{2}+|Y_{2}^{e\tau}|^{2})$. We can expand around the decoupling limit, i.e. $\beta-\alpha\approx\pi/2$, by using eq. (14), to obtain Bizot:2015qqo

[TABLE]

Thus, in order to have large HLFV, i.e. $\mathrm{Br}(h\to\tau\mu)\sim 1\,\%$, we need

[TABLE]

In principle, this can be achieved quite easily. For instance, choosing $c_{\beta}\sim 0.5$, $\lambda_{6}\sim 0.05$, and $m_{H}\sim 2v$, we can obtain the desired branching ratio for $\sqrt{|Y_{2}^{\tau\mu}|^{2}+|Y_{2}^{\mu\tau}|^{2}}\sim 0.002$. In order to have a sizable $\mathrm{Br}(h\to\tau\mu)$, the correct neutrino mass scale can be obtained with very small singly-charged Yukawa couplings $f^{e\tau}$ and $f^{\mu\tau}$.

Note that for a type-III 2HDM there is an upper bound on $\mathrm{Br}(h\to\tau\mu)\cdot\mathrm{Br}(h\to\tau e)$ from combining the rates of $\mu\to e\gamma$ and $\mu e$ conversion  Dorsner:2015mja (which currently saturates the bound), as all combinations of couplings relevant to these HLFV processes enter in CLFV with tau leptons running in the loop, see secs. 3.3.2 and 3.3.4:

[TABLE]

In the Zee model, as we will see, reproducing the leptonic mixing angles correctly implies that there are also lower bounds on the HLFV processes.

4 Numerical analysis

4.1 Scan of the parameter space

In order to study the large parameter space of our scenario of the Zee model and to be able to investigate how large CLFV and HLFV processes can be, we perform a full numerical scan using the software MultiNest Feroz:2007kg ; Feroz:2008xx ; Feroz:2013hea . MultiNest is a Bayesian inference tool that uses so-called nested sampling and especially suitable when there are possibly several maxima in the parameter space. It is designed to determine the Bayesian evidence, but as a byproduct, it also yields the posterior distribution that is relevant for a Bayesian analysis. Nevertheless, it also maximizes the likelihood, which is relevant for a frequentist analysis. We are interested in the maximization of the likelihood and we perform a fully frequentist analysis. We scan over all free parameters in the model, in total 19 real parameters, which are given in table 2 together with their chosen allowed parameter ranges. The plots, including best-fit points and $1\sigma$ and $2\sigma$ confidence regions, are produced using the graphical interface Superplot Fowlie:2016hew .

We impose the stability conditions on the scalar couplings given in eqs. (42), (43), and (47). Direct searches on singly-charged scalars from LEP II imply $m_{h^{+}_{1}},m_{h^{+}_{2}}>80$ GeV Abbiendi:2013hk . In the scan, we assume that all Higgs bosons, except the light one (with the mass fixed to $m_{h}=125.5$ GeV), are heavier than $100$ GeV, i.e. not only the charged ones. This means that $m_{A},m_{H},m_{h^{+}_{1}},m_{h^{+}_{2}}>100$ GeV. The scan over the free parameters is performed for three cases: (i) $\mu=0$, i.e. no neutrino masses,101010In this case, we have only 18 free real parameters, since $\mu$ is set to zero. (ii) $\mu\neq 0$ with neutrino masses in NO, and (iii) $\mu\neq 0$ with neutrino masses in IO.

The quantity that is maximized is the likelihood $L$, which is equivalent to minimizing the $\chi^{2}$ function: $\chi^{2}=-2\ln L$. We assume Gaussian likelihoods, and thus, the contribution to the $\chi^{2}$ function from the $N$ signals is given by

[TABLE]

where $\mathcal{O}^{i}_{\rm th}$ ($\mathcal{O}^{i}_{\rm obs}$) is the theoretical prediction (experimental measurement) of the observable $i$ with the respective standard deviation $\sigma_{i}$. For $s_{23}^{2}$ and $\delta$, we use their combined two-dimensional distribution function (therefore including their correlation) given in refs. Esteban:2016qun ; NuFIT3.0 . For the $M$ upper bounds, we use

[TABLE]

where $B^{j}$ is the experimental upper limit at $1\sigma$ significance level of the observable $j$. According to the definitions for the standard Gaussian distributions assumed for the observables, the 90 % C.L. and 95 % C.L. upper limits are normalized (i.e. divided) by the factors 1.645 and 1.949, respectively. The limits for the signals are rescaled in the same way. The choice of eq. (93) for the $\chi^{2}$ function is made so that we do not penalize deviation from zero when this is not supported by data.

We take into account all the relevant bounds and observables described in section 3, except for the muon AMM discrepancy with respect to the SM, which cannot be accommodated in our present scenario of the Zee model. Thus, adding eqs. (92) and (93), the total $\chi^{2}$ function reads

[TABLE]

In tabs. 1 and 3, we present different observables that yield positive signals (see also appendix A). In table 4, we show the values for the neutrino oscillation parameters from a global fit Esteban:2016qun ; NuFIT3.0 (see also section 2.3). Furthermore, in table 5, observables for which there are only upper bounds are presented, and in table 6, we show the upper bounds for different neutrino mass parameters (see also section 2.3). Finally, we also treat the naturality constraints of eqs. (44) and (46) as bounds and include them in $\chi_{\rm bounds}^{2}$.

4.2 Results of scan

In this section, we discuss the main results of our numerical scan, performed for $\kappa=1$ with the naturality upper limit of $\mu\lesssim 1.5$ TeV, see eq. (46), unless otherwise stated. At the end of the section, we discuss how the results would change for $\kappa=10$ ($\mu\lesssim 15$ TeV). The numerical scan is performed for three different cases of our scenario of the Zee model, i.e. for $\mu=0$, NO, and IO. For these cases, the values for the minima of the $\chi^{2}$ function, $\chi_{\rm min}^{2}$, at the respective best-fit points are

1. (i)

$\mu=0$, i.e. massless neutrinos: $\chi_{\rm min}^{2}\simeq 5.1$, 2. (ii)

$\mu\neq 0$, massive neutrinos in NO: $\chi_{\rm min}^{2}\simeq 10.7\,(11.0)$ for $\kappa=1\,(10)$, 3. (iii)

$\mu\neq 0$, massive neutrinos in IO: $\chi_{\rm min}^{2}\simeq 21.7\,(21.5)$ for $\kappa=1\,(10)$.

Thus, in our scenario of the Zee model, neutrino masses and leptonic mixing can be accommodated in both NO and IO. However, IO is disfavored compared to NO. We also note that if $\theta_{23}$ will turn out to be in the second octant, IO cannot be accommodated.

In figure 2, we present the contributions to $\chi_{\rm min}^{2}$ from the observables in both NO and IO. Note that we do not show the contributions from the upper bounds, since these are always satisfied, and thus, the corresponding contribution to $\chi_{\rm min}^{2}$ is exactly zero. Note that we present the combined contribution from $s^{2}_{23}$ and $\delta$ to $\chi^{2}_{\rm min}$. In IO, the largest contributions stem from $s^{2}_{12}$ and $s^{2}_{23}+\delta$, although they are all within $3\sigma$ of their experimental values. In NO, the corresponding contributions are small. These observables account for the fact that the fit in IO is much worse than in NO. All other observables in IO are of the same size as those in NO and within $2\sigma$. In addition, we perform a run assuming $Y_{2}^{\mu\tau}=0$, which renders one neutrino massless, see eq. (38). In this case, we find that IO is significantly better than NO, which is in agreement with the results of ref. He:2011hs , and the value of $\chi_{\rm min}^{2}$ in IO is of the same size as for $Y_{2}^{\mu\tau}\neq 0$ (i.e. quite large), whereas in NO, the difference is more than one order of magnitude with $\chi_{\rm min}^{2}\sim\mathcal{O}(100)$. Thus, in this scenario, both NO and IO are basically excluded. Moreover, there are non-negligible contributions to $\chi_{\rm min}^{2}$ from the Higgs signals, especially $h\to ZZ$. Regarding $h\to\gamma\gamma$, we find that it is around the SM value. Note that in this case we choose some of the couplings of the scalar potential to be zero, see eqs. (86), so its value can be modified by turning them on. The values for the universality parameters $g_{ij}$, see eqs. (71)–(73), are always very close to one, even closer than the experimental values, and compatible with observations at $2\sigma$. However, these also give non-negligible contributions to $\chi_{\rm min}^{2}$.

In the following, we present the allowed regions for the most interesting parameters and observables. In figure 3, we plot the leptonic mixing parameters $s^{2}_{12}$ and $s^{2}_{23}$ for NO (left panel) and IO (right panel), where one can clearly see that the fit is very good for NO, while for IO it crucially depends on the octant of $\theta_{23}$. In fact, the fitted value of $s^{2}_{12}$ would be $5\sigma$ away from its experimental best-fit value if $\theta_{23}$ lies in the first octant. For NO, it is clear that both octants are viable.

In figure 4, we plot the allowed regions of the branching ratios ${\rm Br}(h\to\tau\mu)$ and ${\rm Br}(\tau\to\mu\gamma)$. For $\mu=0$ (left panel), only upper bounds on the CLFV and HLFV processes exist and these are compatible with the results of refs. Sierra:2014nqa ; Dorsner:2015mja . For $\mu\neq 0$, one can observe that there are lower bounds on the CLFV and HLFV processes in both NO (middle panel) and IO (right panel). That is, reproducing neutrino masses implies that CLFV and HLFV cannot be arbitrarily small. In particular, ${\rm Br}(h\to\tau\mu)\gtrsim 10^{-6}\,(10^{-7})$ in NO (IO). For NO, the upper bound on ${\rm Br}(h\to\tau\mu)$ saturates the experimental one, while for IO, we obtain ${\rm Br}(h\to\tau\mu)\lesssim 5\cdot 10^{-3}$.

In the future, Belle II is expected to reach a sensitivity on ${\rm Br}(\tau\to\mu\gamma)$ of $\mathcal{O}(10^{-9})$ HerediadelaCruz:2016yqi , which would significantly probe a substantial part of the allowed parameter space of IO and almost the complete allowed region of NO. This is one of the most interesting results of our work. Similarly, if sensitivities of ${\rm Br}(h\to\tau\mu)$ at future colliders reach $10^{-4}\,(10^{-5})$, NO (IO) will be tested at 68 % C.L.

In figure 5, we show the allowed regions of the branching ratios ${\rm Br}(h\to\tau e)$ and ${\rm Br}(h\to\tau\mu)$ for NO (left panel) and IO (right panel). No correlation exists for $\mu=0$, while there is a strong correlation for $\mu\neq 0$, stronger for IO than for NO. In general, we find that ${\rm Br}(h\to\tau e)\lesssim 10^{-2}\,{\rm Br}(h\to\tau\mu)$ or even lower for both orderings. Therefore, observations of $h\to\tau e$ will be considerably more challenging than for $h\to\tau\mu$.

In figure 6, we display the allowed regions of the $\mu e$ conversion rate in gold and the branching ratio ${\rm Br}(h\to\tau\mu)$ for NO (left panel) and IO (right panel). As discussed in section 3.3.4, next generation experiments, using aluminium and titanium, are expected to achieve an improved sensitivity of up to about four orders of magnitude, maybe reaching $\mathcal{O}(10^{-19})$ Donghia:2016lzt . The $\mu e$ conversion rates of these materials are of the same order of magnitude as that of gold,111111Quantitatively, the $\mu e$ conversion rates for Al and Ti scale as $\rm{Cr}(\mu\to e)_{\rm Al}\simeq 0.5\cdot\rm{Cr}(\mu\to e)_{\rm Au}$ and $\rm{Cr}(\mu\to e)_{\rm Ti}\simeq 0.8\cdot\rm{Cr}(\mu\to e)_{\rm Au}$, respectively. and therefore, if a negative result is obtained, IO would be excluded, while there would still be a considerable allowed region for NO. Thus, there is a complementarity between ${\rm Br}(\tau\to\mu\gamma)$ and $\mu e$ conversion. On the other hand, the sensitivity of ${\rm Br}(\tau\to e\gamma)$ is expected to reach about $3\cdot 10^{-9}$ HerediadelaCruz:2016yqi and ${\rm Br}(\mu\to e\gamma)$ is expected to be improved by one order of magnitude, but these are not able to test the model as thoroughly as ${\rm Br}(\tau\to\mu\gamma)$ and $\mu e$ conversion.

In figure 7, we present the allowed regions for the absolute values of the Yukawa couplings $Y_{2}^{\tau e}$ and $Y_{2}^{\tau\mu}$ in both NO (left panel) and IO (right panel). We find that the regions are quite well defined, although they are larger in IO than in NO, and the allowed values are larger for $Y_{2}^{\tau\mu}$ than for $Y_{2}^{\tau e}$. Note that the value of $Y_{2}^{\tau\mu}$, which is a very relevant parameter since it controls the decays of the scalars into $\tau\mu$ (together with $Y_{2}^{\mu\tau}$), is always larger than around $10^{-3}\,(2.5\cdot 10^{-3})$ in NO (IO).

The scale of neutrino masses is controlled by the charged scalar mixing angle $\varphi$, which is proportional to the trilinear coupling $\mu$, see eq. (8), and the Yukawa couplings $f$ and $Y_{2}$. Therefore, in figure 8, we plot the allowed regions of $|f^{\mu\tau}|$ and $s_{2\varphi}$ for both NO (left panel) and IO (right panel). The range of $|f^{\mu\tau}|$ is similar in both orderings. One can clearly see that $s_{2\varphi}$ is close to zero for $|f^{\mu\tau}|\gtrsim 10^{-5}$, while it grows very fast for $|f^{\mu\tau}|\lesssim 10^{-5}$, reaching values of $0.7\,(0.4)$ for NO (IO). As expected, the naturality condition of eq. (46), which is added to the $\chi^{2}$ function, restricts $\mu$ to be smaller than about $3\,(30)$ TeV for $\kappa=1\,(10)$ at $2\sigma$.

The Yukawa couplings $f^{\mu\tau}$ and $f^{e\tau}$ are always in the range $[10^{-7},0.1]$ and have similar allowed regions in both NO and IO. In fact, they are highly correlated with $f^{\mu\tau}\simeq f^{e\tau}\,(0.1f^{e\tau})$ in NO (IO). Their allowed $1\sigma$ C.L. regions lie roughly below $10^{-5}$, suppressing all interactions mediated by the antisymmetric Yukawa coupling $f$ of the singly-charged scalar singlet. Therefore, in order to describe neutrino masses and leptonic mixing, $Y_{2}$ should be much larger than $f$.

In figure 9, we show the allowed regions of the Higgs scalar mass differences $m_{H}-m_{A}$ and $m_{H}-m_{h_{1}^{+}}$. These affect the size of the parameter $T$, see appendix A. We find similar results to those in ref. Benbrik:2015evd . At $1\sigma$ C.L. and at low scalar masses, $m_{H}-m_{h_{1}^{+}}$ is roughly equal to $m_{H}-m_{A}$, thus canceling the contributions to $T$. At $2\,\sigma$ C.L., for the two cases $\mu=0$ and IO, $m_{H}$ and $m_{A}$ can be close to each other and still fulfill $m_{H}\geq m_{h_{1}^{+}}$. The allowed ranges of $m_{A}$ and $m_{H}$ are the same, but the two are not completely correlated, especially for IO. In fact, it is possible to have $m_{A}\simeq 100$ GeV, while $100~{}{\rm GeV}\lesssim m_{H}\lesssim 500~{}{\rm GeV}$.

In addition, we find that the mass of the lightest charged scalar is $m_{h^{+}_{1}}\lesssim 0.9\,(1.7)$ TeV at $1\,(2)\,\sigma$ C.L. in NO, while $m_{h^{+}_{1}}\lesssim 0.7\,(1.1)$ TeV at $1\,(2)\,\sigma$ C.L. in IO. The other charged scalar of the model, i.e. $h_{2}^{+}$, is always heavier than $h_{1}^{+}$, see eq. (9), and can reach values of $\mathcal{O}(100)$ TeV in both orderings. The larger the $m_{h_{2}^{+}}$, the smaller the $s_{2\varphi}$, see eq. (8), and therefore, the smaller the neutrino masses.121212The heavy scalars contribute to $m_{h}$ at one-loop level and could, in principle, pose a problem for naturality, especially $h^{+}_{2}$ as it is the heaviest one. However, their contributions to the Higgs boson self-energy are suppressed by both the loop factor and the quartic couplings of the scalar potential and $m_{h}$ is therefore natural for the values used.

From eq. (89), we know that ${\rm Br}(h\to\tau\mu)$ is proportional to $1/m^{4}_{H}$, i.e., it decouples with the CP-even scalar mass. In figure 10, we display the allowed regions of the mass $m_{H}$ and the branching ratio ${\rm Br}(h\to\tau\mu)$. For NO (left panel), $m_{H}\lesssim 0.9\,(1.7)$ TeV at $1\,(2)\,\sigma$ C.L., while for IO (right panel), $m_{H}\lesssim 0.7\,(1.1)$ TeV at $1\,(2)\,\sigma$ C.L. Therefore, if an extra CP-even scalar (and close by CP-odd and charged scalars, see figure 9) is observed below $0.9\,(0.7)$ TeV, then one expects ${\rm Br}(h\to\tau\mu)\gtrsim 10^{-4}\,(10^{-5})$ for NO (IO) at $1\sigma$ C.L. The heavy CP-even scalar $H$ could also have sizable decays into $\tau\mu$, depending on the scalar spectrum.

In figure 11, we show $\tan{\beta}$ as a function of $s_{\alpha}$ for $\mu=0$ (left panel), NO (middle panel), and IO (right panel). Having Higgs boson decays close to the observed ones implies being close to the decoupling limit, i.e. $s_{\beta-\alpha}\to 1$, and therefore, $\tan{\beta}$ and $s_{\alpha}$ are strongly correlated. For massless neutrinos, $\tan{\beta}$ can reach values up to 15 and $s_{\alpha}$ can approach zero. However, if neutrino masses are introduced, the value of $\tan{\beta}$ is severely constrained to smaller values in both orderings and the allowed range of $s_{\alpha}$ is also reduced. The upper bound on $\tan{\beta}$ is clearly more severe for NO, where $\tan{\beta}$ is driven to the smallest possible values (i.e. below 0.5), while for IO, it can reach values up to 1.4. We have also performed a run forcing the value of the unphysical parameter $\tan{\beta}$ to be large, i.e. $40\lesssim\tan{\beta}\lesssim 50$, and we find that the fit becomes significantly worse, leading to a value for the minimum of the $\chi^{2}$ function of $\mathcal{O}(1000)$, and thus excluding this scenario.

In figure 12, we present the allowed regions of the effective mass parameter $m_{ee}$ that appears in neutrinoless double beta decay and the smallest neutrino mass ($m_{1}$ for NO and $m_{3}$ for IO) for NO (left panel) and IO (right panel). We can see that the smallest neutrino mass is several orders of magnitude smaller in IO than in NO. In fact, in IO, it can be massless, whereas this is not the case in NO. This is consistent with the fact that the fit in NO is bad when $Y_{2}^{\mu\tau}=0$, rendering one neutrino massless. In NO, we obtain $|m_{ee}|\simeq(4-5)\,\mathrm{meV}$, while in IO, $|m_{ee}|$ is one order of magnitude larger, i.e. about 50 meV, and thus, it will be possibly probed in planned neutrinoless double beta decay experiments. Furthermore, the mass of the lightest neutrino mass eigenstate is less than $1.4\cdot 10^{-3}$ eV ($1.3\cdot 10^{-4}$ eV) in NO (IO). Note that in NO there is a lower bound on this mass, $6\cdot 10^{-4}$ eV, while IO is compatible with a massless neutrino.

In addition, we mention the phases of the leptonic mixing matrix $U$. We find that the preferred value for the Dirac CP-violating phase $\delta$ is close to $2\pi$ in NO and close to $\pi$ in IO, thus implying no leptonic CP violation in our scenario of the Zee model. The hint of $\delta\sim 3\pi/2$ from global fits to neutrino oscillation data, if confirmed, can therefore not be accommodated in any of the two orderings. We also find that the values of both Majorana CP-violating phases $\phi_{1}$ and $\phi_{2}$ are around $2\pi$ in both orderings.

As expected, the strong limits from other LFV processes imply that the NSI parameters $\chi^{\rm m}_{\tau\tau}$ and $\varepsilon^{\rm m}_{\tau\tau}$ defined in eq. (76) are very suppressed in both orderings, i.e. $\chi^{\rm m}_{\tau\tau},\varepsilon^{\rm m}_{\tau\tau}\lesssim 10^{-8}$, and therefore well below future experimental sensitivity. The NSI parameters $\varepsilon^{\rm\rho\sigma}_{\alpha\beta}$, given in eq. (75), are also generated in the Zee–Babu model, where they reach values of about $10^{-4}$ Ohlsson:2009vk ; Herrero-Garcia:2014hfa . However, in the Zee model, NSIs turn out to be smaller, since neutrino masses are generated at one-loop level, while in the Zee–Babu model, the latter arise at two loops.

Furthermore, we check how our results change when imposing the fine-tuning parameter $\kappa$ to be 10 instead of 1, see eq. (46). For both values of $\kappa$, the $1\sigma$ C.L. region is close to the upper limit on $\mu$, even though the entire range down to $\mu=1$ GeV is allowed at $2\sigma$ C.L. The value of $\mu$ is not significantly correlated to the value of $\tan{\beta}$. We find that the allowed ranges for the scalar masses (where the upper bounds determine the neutrino masses) depend critically on $\kappa$, having larger allowed mass ranges the larger the value of $\kappa$. Quantitatively, for $\kappa=1$ at $1\,(2)\,\sigma$ C.L., the upper bounds on the masses are $m_{A},m_{H},m_{h_{1}^{+}}=0.9\,(1.7)$ TeV for NO and $m_{A},m_{H},m_{h_{1}^{+}}=0.7\,(1.1)$ TeV for IO, whereas for $\kappa=10$ at $1\,(2)\,\sigma$ C.L., the upper bounds are $m_{A},m_{H},m_{h_{1}^{+}}=1.6\,(2.5)$ TeV for NO and $m_{A},m_{H},m_{h_{1}^{+}}=0.9\,(1.4)$ TeV for IO.

Finally, we summarize some of our main results in table 7, where we display the $2\sigma$ C.L. regions for some of the most interesting observables and parameters for NO (for $\kappa=1,10$), and IO (for $\kappa=1,\,10$). We emphasize once more that for $\mu=0$ only upper bounds on the CLFV processes (and lower bounds on the scalar masses) exist, whereas for $\mu\neq 0$, there are lower bounds for the CLFV processes (and upper bounds on the scalar masses). This means that the scalar sector cannot be arbitrarily heavy if neutrino masses are to be reproduced. The precise upper bound depends crucially on $\mu$, see eq. (46). For $\mu=0$, the masses could be arbitrarily large, unlike the case of having $\mu\neq 0$ and reproducing neutrino masses, which imposes that they are below about 2 TeV. We do not display the ranges for all observables, such as neutrino masses and leptonic mixing parameters, as their contributions to $\chi_{\rm min}^{2}$ are shown in figure 2. In the Zee model (including $\mu=0$), the value of the muon AMM, which has not been included in the fit, is several orders of magnitude smaller than the experimental one. This implies a deviation of about $3.5\sigma$. The allowed ranges for the scalar masses $m_{A}$ and $m_{H}$ are the same for both values of $\kappa$, but the two are not completely correlated (especially in IO). It is possible to simultaneously have $m_{A}\simeq 100$ GeV and $m_{H}$ varying in the range $100~{}{\rm GeV}\lesssim m_{H}\lesssim 500~{}{\rm GeV}$. We conclude by stating that NO will be tested in the next generation CLFV searches, specially with $\tau\to\mu\gamma$ and $\mu e$ conversion, as well as searches for $h\to\tau\mu$ and the other new scalars at colliders.

5 Summary and conclusions

It is well known that there is LFV in the neutrino sector and this is also expected in the charged-lepton sector. In this work, we have studied the Zee model in detail, which is a simple extension of the SM that can accommodate neutrino masses and leptonic mixing if at the same time sizable signals in LFV processes are generated. We have performed a full numerical scan of the parameter space for three different cases (i) $\mu=0$, which implies massless neutrinos, (ii) NO, and (iii) IO. We have found that neutrino masses and leptonic mixing can be easily accommodated in NO, whereas IO is disfavored in comparison to NO due to the difficulty to fit the leptonic mixing angles $\theta_{12}$ and $\theta_{23}$ as well as the Dirac CP-violating phase $\delta$. In fact, if $\theta_{23}$ turns out to be in the first octant, only NO would be allowed in the Zee model. Note also that none of the orderings can reproduce the hint of $\delta\simeq 3\pi/2$ from global fits to neutrino oscillation data Esteban:2016qun ; NuFIT3.0 .

If expected sensitivities are achieved in $\tau\to\mu\gamma$ and no signal is observed, a significant portion of the allowed parameter space for NO would be ruled out. This would put the Zee model under severe pressure, requiring an extension, e.g. involving the Yukawa couplings that give rise to terms proportional to $m_{e}$ and possibly large hierarchies among them, and relaxing the naturality demands on the trilinear coupling $\mu$ considerably, as ${\rm Br}(\tau\to\mu\gamma)\sim 10^{-9}$ would still rule out NO for $\kappa\sim 10$. If no signals are observed in future $\mu e$ conversion experiments, which are expected to increase their sensitivities by several orders of magnitude, the allowed regions in the parameter space will be strongly reduced, and IO will be basically excluded.

We have analyzed if the predicted rates of the Zee model for HLFV decays are observable at the LHC and future colliders. We have found that $\mathrm{Br}(h\to\tau\mu)$ can be at the percent level, whereas $\mathrm{Br}(h\to\tau e)$ is at least two orders of magnitude smaller. If no signals are observed for ${\rm Br}(h\to\tau\mu)$ at future colliders at the $10^{-5}$ level, both orderings will be excluded at $1\sigma$ C.L. In general, we find that the expected sensitivities for CLFV will be more constraining in the near future. However, both CLFV and HLFV will have significant impact on the allowed parameter space of the Zee model.

In the model, neutrinoless double beta is due to only light neutrinos. Therefore, as usual, current experiments will be only sensitive to IO, while NO will only be tested if a further-order-of-magnitude improvement is achieved. In the Zee model, NSIs are always very suppressed, the strong limits from other CLFV processes and the fact that neutrino masses need to be generated at one-loop level.

In general, we have found that the masses of the new scalars should be at most a few TeV, which implies that they, especially the charged scalars that are pair-produced via Drell–Yan processes, can be searched for at the LHC. In particular, the masses of the neutral scalars and the charged scalar $h_{1}^{+}$ are below $2.5$ TeV, and typically they are lower than that, for both NO and IO. The phenomenology of the scalar sector is very model-dependent, like in general for 2HDMs, although it is possible to have sizable decays of the heavy neutral scalars into $\tau\mu$, correlated with the light Higgs ones.

We conclude by emphasizing that the general Zee model studied in this work is fully testable in the near future by combining different LFV processes. In particular, both orderings should be completely tested by CLFV and HLFV processes in the forthcoming years, as well as collider searches of the new scalars. Furthermore, if a signal in $h\to\tau\mu$ is observed at the LHC or in a future collider, the Zee model will be one of the best-motivated scenarios to accommodate it and at the same time describe neutrino masses and leptonic mixing.

Acknowledgements

We thank Nuria Rius, Arcadi Santamaria, and Thomas Schwetz for useful discussions and comments on the manuscript, Mikael Twengström for technical advice on numerical computations, and Andrew Fowlie for providing help with Superplot. We also acknowledge Thomas Schwetz for suggesting to show the individual contributions to the $\chi^{2}$ function and for discussions regarding global fits of neutrino oscillation data. Numerical computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC Center for High Performance Computing (PDC-HPC) at KTH Royal Institute of Technology in Stockholm, Sweden under project numbers PDC-2016-8, PDC-2016-27, PDC-2016-36, PDC-2016-60, and PDC-2016-82. This work was supported by the University of Adelaide and the Australian Research Council through the ARC Centre of Excellence for Particle Physics at the Terascale (CoEPP) (CE110001104).

Appendix A Electroweak precision tests

The Peskin–Takeuchi parameters $S$, $T$, and $U$ give a parametrization of the new physics contributions to electroweak radiative quantities, in particular to gauge boson self-energies Peskin:1990zt ; Peskin:1991sw . We follow closely ref. Haber:2010bw , extending their results for the 2HDM to the Zee model by adding the extra contributions stemming from the singly-charged scalar singlet $h^{+}$. We refer the reader to ref. Haber:2010bw for additional details on the evaluations of the one-loop self-energies.

In the Zee model, the parameter $T$ is given by131313Using ref. Grimus:2007if , one would obtain a different factor in front of the first term on the third line: $-2s^{2}_{\varphi}c^{2}_{\varphi}\mathcal{F}(m_{h_{1}^{+}}^{2},m_{h_{2}^{+}}^{2})$.

[TABLE]

where $\alpha_{\rm em}\equiv e^{2}/(4\pi)$ is Sommerfeld’s fine-structure constant141414Note that $\pi\alpha_{\rm em}v^{2}=m_{W}^{2}s_{W}^{2}$, where $s_{W}=\sin\theta_{W}$ and $\theta_{W}$ being the Weinberg angle. and the symmetric auxiliary function $\mathcal{F}$ is defined as

[TABLE]

Similarly, for the parameter $S$, following ref. Haber:2010bw and adding the singly-charged contributions to the different gauge boson self-energies, we obtain

[TABLE]

and for the combination $S+U$, we find151515We believe there are typos in the last two terms of ref. Haber:2010bw , which should have opposite signs.

[TABLE]

where the renormalized auxiliary functions $\mathcal{B}_{22}$ and $\mathcal{B}_{0}$ are defined as

[TABLE]

with the Passarino–Veltman functions $B_{22}$ and $B_{0}$ Passarino:1978jh , arising from two-point self-energies. Using eqs. (97) and (98), one can readily obtain an expression for the parameter $U$. In the limit $s_{\varphi}=0$, it can be easily checked that one recovers the 2HDM results of the electroweak precision tests of ref. Haber:2010bw . Finally, in dimensional regularization, the two functions $B_{22}$ and $B_{0}$ read Haber:2010bw

[TABLE]

with

[TABLE]

in $d$ space-time dimensions, where $\gamma\simeq 0.577$ is the Euler–Mascheroni constant. Note that $B_{22}$ and $B_{0}$ are symmetric in their last two arguments. In appendix B, we derive explicit analytical expressions for $B_{22}$ and $B_{0}$ as well as $\mathcal{B}_{22}$ and $\mathcal{B}_{0}$.

Appendix B Explicit analytical expressions for the Passarino–Veltman functions $B_{0}$ and $B_{22}$ and the renomalized auxiliary functions $\mathcal{B}_{0}$ and $\mathcal{B}_{22}$

Following the seminal work by Passarino and Veltman Passarino:1978jh closely, it holds that161616Note that we use a different sign convention for $B_{22}$ than ref. Passarino:1978jh .

[TABLE]

where the additional Passarino–Veltman function $A_{0}$ is given by

[TABLE]

Furthermore, we obtain

[TABLE]

where $x_{1}$ and $x_{2}$ are the roots of the equation $q^{2}x^{2}+(m_{1}^{2}-m_{2}^{2}-q^{2})x+m_{2}^{2}=0$. Inserting eqs. (105), (106), and (107) into eq. (104), it follows after some tedious calculations that the functions $B_{22}$ and $B_{0}$ can be written as

[TABLE]

where

[TABLE]

In fact, it is possible to find closed-form expressions even for the auxiliary functions $\mathcal{B}_{22}$ and $\mathcal{B}_{0}$. Inserting eqs. (108) and (109) into eqs. (99) and (100), respectively, and using the functions $B_{22}$ and $B_{0}$ evaluated at $q^{2}=0$, i.e.

[TABLE]

we obtain for $m_{1},m_{2}>q>0$ and $(m_{1}-m_{2})^{2}>q^{2}$

[TABLE]

whereas for $m_{1},m_{2}>q>0$ and $(m_{1}-m_{2})^{2}<q^{2}$, we have to make the replacements $Y\to iY^{\prime}$ and

[TABLE]

in eqs. (113) and (114), respectively, where $Y^{\prime}\equiv\sqrt{[(m_{1}+m_{2})^{2}-q^{2}][q^{2}-(m_{1}-m_{2})^{2}]}$. In the case when $m_{1}=q$ and $m_{2}=m$ (which is useful for computing eqs. (97) and (98)), we have for $0<m/2<q<m$

[TABLE]

where

[TABLE]

whereas for $0<q<m/2$, we have to make the replacements $Z\to iZ^{\prime}$ and

[TABLE]

in eqs. (115) and (116), respectively, where $Z^{\prime}\equiv\sqrt{m^{2}-4q^{2}}$. Similarly, in the case when $m_{1}=m_{2}=m>q>0$, we find

[TABLE]

where

[TABLE]

B.1 Comments on cancellation of dimensionful logarithms and divergent terms in $S$ and $S+U$

Note that the terms proportional to the dimensionful logarithm $\ln(m_{1}m_{2})$ in eq. (113) as well as in eqs. (115) and (118) sum up to contributions for both $S$ in eq. (97) (contribution of seven $\mathcal{B}_{22}$ functions) and $S+U$ in eq. (98) (contribution of ten $\mathcal{B}_{22}$ functions) that are proportional to $c_{\varphi}^{2}[\ln(m_{A}/m_{h_{1}^{+}})+\ln(m_{H}/m_{h_{1}^{+}})]+s_{\varphi}^{2}[\ln(m_{A}/m_{h_{2}^{+}})+\ln(m_{H}/m_{h_{2}^{+}})]$, which is indeed a result of dimensionless logarithms. Furthermore, it should be noted that all terms on the form $-(\Delta+1)q^{2}/12$ (which contains a divergency) in eq. (99) with eq. (101) cancel exactly for both $S$ and $S+U$ in eqs. (97) and (98), respectively, as they must. This also holds true for the divergent terms in $T$ in eq. (95).