The Zee-Babu Model revisited in the light of new data

We update previous analyses of the Zee-Babu model in the light of new data, e.g., the mixing angle $\theta_{13}$, the rare decay $\mu\to e \gamma$ and the LHC results. We also analyse the possibility of accommodating the deviations in $\Gamma(H\to \gamma\gamma)$ hinted by the LHC experiments, and the stability of the scalar potential. We find that neutrino oscillation data and low energy constraints are still compatible with masses of the extra charged scalars accessible to LHC. Moreover, if any of them is discovered, the model can be falsified by combining the information on the singly and doubly charged scalar decay modes with neutrino data. Conversely, if the neutrino spectrum is found to be inverted and the CP phase $\delta \neq \pi$, the masses of the charged scalars will be well outside the LHC reach.


I. INTRODUCTION
The observed pattern of neutrino masses and mixing remains one of the major puzzles in particle physics. Moreover, massive neutrinos provide irrefutable evidence for physics beyond the Standard Model (SM) and many theoretical possibilities have been proposed to account for the lightness of neutrinos. With the running of the LHC, it is timely to explore neutrino mass models in which the scale of new physics is close to the TeV. In particular, radiative mechanisms are especially appealing, since small neutrino masses are generated naturally due to loop factors. On the other hand, new physics effects can be sizable also in low energy experiments, for instance lepton flavour violating rare decays of charged leptons, α → β γ, providing complementary probes for such models.
In this paper we consider the Zee-Babu model (ZB) of neutrino masses 1 , which just adds two (singly and doubly) charged scalar singlets to the SM. Neutrino masses are generated at two loops and are proportional to the Yukawa couplings of the new scalars and inversely proportional to the square of their masses. This is phenomenologically quite interesting because the new scalars cannot be very heavy or have very small Yukawa couplings, otherwise neutrino masses would be too small. As a consequence, such scalars may be accessible at the LHC, and in principle they could explain the slight excess over the SM prediction found by ATLAS in the diphoton Higgs decay channel H → γγ (currently CMS does not see any excess, see section III for the latest data). They also mediate a variety of lepton flavour violating (LFV) processes, leading to rates measurable in current experiments.
The phenomenology of the ZB model has been widely analyzed: neutrino oscillation data was used to constrain the parameter space of the model, LFV charged lepton decay rates calculated and collider signals discussed [6][7][8]. Non-standard neutrino interactions in the ZB model have also been thoroughly studied, in correlation with possible LHC signals and LFV processes [9]. In [8], some of us performed an exhaustive numerical study of the full parameter space of the model using Monte Carlo Markov Chain (MCMC) techniques, which allow to efficiently explore high-dimensional spaces. The purpose of this work is twofold: (1) to update such analysis, in the light of the recent measurement of the neutrino mixing angle θ 13 , the new MEG limits on µ → eγ, the lower bounds on doubly-charged scalars coming from LHC data, and, of course, the discovery of a 125 GeV Higgs boson by ATLAS and CMS [10,11] and (2) also study the possibility of accommodating deviations from the SM prediction for the Higgs diphoton decay channel, and the effects of the new couplings of the model in the stability of the scalar potential. A possible enhancement of the Higgs diphoton decay rate in the ZB model together with the vacuum stability of the scalar potential has been studied in [12], however a consistent updated analysis including all constraints is lacking.
The outline of the paper is the following. In section II we briefly review the main features of the ZB model, discussing perturbativity and naturality estimates for the allowed ranges of the free parameters of the model. We summarize present constraints from recent neutrino oscillation data, low energy lepton-flavour violating processes, universality and stability of the scalar potential. We also review the collider phenomenology of the ZB model, discussing current limits from LHC, and briefly comment on the prospects for non-standard neutrino interactions. In section III we analyze in detail the contributions of the ZB charged scalars to both, Γ(H → γγ) and Γ(H → Zγ). After some analytic estimates in section IV, we present the results of our MCMC numerical analysis in section V and we conclude in section VI.
Renormalization group equations for the ZB model and relevant loop functions are collected in the appendices.

II. THE ZEE-BABU MODEL
We follow the notation of [8]. As mentioned above, the Zee-Babu model only contains, in addition to the SM, two charged singlet scalar fields with weak hypercharges ±1 and ±2 respectively (we use the convention Q = T 3 + Y ).
Notice that we can assign lepton number −2 to both scalars, h + and k ++ , in such a way that total lepton number L (or B−L) is conserved in the complete Lagrangian, except for the trilinear coupling µ of the scalar potential; thus, lepton number is explicitly broken by the µcoupling. It is important to remark that lepton number violation requires the simultaneous presence of the four couplings Y , f , g and µ, because if any of them vanishes one can always assign quantum numbers in such a way that there is a global U (1) symmetry. This means that neutrino masses will require the simultaneous presence of the four couplings.
Regarding the physical free parameters in the ZB model, our convention is the following: without loss of generality, we choose the 3 × 3 charged lepton Yukawa matrix Y to be diagonal with real and positive elements. We also use fermion field rephasings to remove three phases from the elements of the matrix g and charged scalar rephasings to set µ real and positive, and to remove one phase from f . In summary we have 12 moduli (3 from Y , 3 from f and 6 from g ab ), 5 phases (3 from g and 2 from f ) and the real and positive parameter µ, plus the rest of real parameters in the scalar potential. As discussed in [8], this choice is compatible with the standard parametrization of neutrino masses and mixings.
After electroweak symmetry breaking, the masses of charged leptons are m a = Y aa v, with v ≡ H 0 = 174 GeV, the VEV of the standard Higgs doublet, while the physical charged scalar masses are given by In principle, the scale of the new mass parameters of the ZB model (m h , m k and µ) is arbitrary. However from the experimental point of view it is interesting to consider new scalars light enough to be produced in the second run of the LHC. Also theoretical arguments suggest that the scalar masses should be relatively light (few TeV), to avoid unnaturally large one-loop corrections to the Higgs mass which would introduce a hierarchy problem. Therefore, in this paper we will focus on the masses of the new scalars, m h , m k , below 2 TeV.
The Yukawa couplings of the new scalars of the model enter in the neutrino mass formula and in several LFV processes, and are strongly bounded for the scalar masses we are considering except in a few corners of the parameter space where we require that the theory remains perturbative. Since one-loop corrections to Yukawa couplings are order one expects from perturbativity f, g 4π, although, as we will see, for the scalar masses considered here, phenomenological constraints are always stronger.
The couplings of the charged scalars in the scalar potential, apart from the stability constraints described in section II E, are essentially free. However, for the theory to make sense as a perturbative theory we also impose the limit λ h,k,kH,hH,hk < 4π.
The trilinear coupling among charged scalars µ, on the other hand, is different, for it has dimensions of mass and it is insensitive to high energy perturbative unitarity constraints.
However, it induces radiative corrections to the masses of the charged scalars of order Requiring that the corrections in absolute value are much smaller than the masses we can derive a naive upper bound for this parameter, µ 4π min(m h , m k ), but it is difficult to fix an exact value of µ for which the contributions to the scalar masses are unacceptably large, leading to a highly fine-tuned scenario.
A large value of µ, as compared with scalar masses, is also disfavoured because it could lead to a deeper minimum of the scalar potential for non-vanishing values of the charged fields, therefore breaking charge conservation. This phenomenon has also been studied in the context of supersymmetric theories (see for instance [13][14][15]). As an example, by looking at the particular direction |H| = |h| = |k| = r, and requiring that the charge breaking minimum is not a global minimum, V (r = 0) > 0, one obtains Assuming no cancellations between the λ's or mass terms, neglecting λ H and m 2 H , and using the perturbative limit for the rest of the couplings λ i < ∼ 4π one finds a very conservative bound on µ Tighter limits can be obtained by looking at all directions in the potential and/or allowing for cancellations.
Given that the neutrino masses depend linearly on the parameter µ, as we will see in the next section, the ability of the model to accommodate all present data is quite sensitive to the upper limit allowed for µ. Thus we choose to implement such limit in terms of a parameter κ, and discuss our results for different values of κ = 1, 5, 4π. Notice that we are using the naturality upper bound (expressed in terms of min(m h , m k )), which in general is much more restrictive than the upper bound obtained by requiring that the minimum of the potential does not break charge conservation (expressed in terms of max(m h , m k )).
The lowest order contribution to neutrino masses involving the four relevant couplings appears at two loops [1,2] and its Feynman diagram is depicted in fig. 1.
The calculation of this diagram gives the following mass matrix for the neutrinos (defined as an effective term in the Lagrangian where I ab is the two-loop integral, which can be calculated analytically [16]. However, since m c , m d are the masses of the charged leptons, necessarily much lighter than the charged scalars, we can neglect them and obtain a much simpler form whereĨ(r) is a function of the ratio of the masses of the scalars r ≡ m 2 k /m 2 h , which is close to one for a wide range of scalar masses. Within this approximation the neutrino mass matrix can be directly written in terms of the Yukawa coupling matrices, f , A very important point is that since f is a 3 × 3 antisymmetric matrix, det f = 0 (for 3 generations), and therefore det M ν = 0. Thus, at least one of the neutrinos is exactly massless at this order.
The neutrino Majorana mass matrix M ν can be written as where D ν is a diagonal matrix with real positive eigenvalues, and U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) leptonic mixing matrix. We are left with only two possibilities for the neutrino masses, m i : • Normal hierarchy (NH): the solar squared mass difference is ∆ S = m 2 2 , the atmospheric mass splitting ∆ A = m 2 3 and m 1 = 0, with m 3 m 2 .
• Inverted hierarchy (IH): The standard parametrization for the PMNS matrix is where c ij ≡ cos θ ij , s ij ≡ sin θ ij and since one of the neutrinos is massless, there is only one physical Majorana phase, φ, in addition to the Dirac phase δ.
Process Experiment (90% CL) Bound (90% CL) In order to provide neutrino masses compatible with experiment, the Yukawa couplings of the charged scalars cannot be too small and their masses cannot be too large. This immediately gives rise to a series of flavour lepton number violating processes, as for instance µ − → e − γ or µ − → e + e − e − , with rates which can be, in some cases, at the verge of the present experimental limits. Therefore, we can use these processes to obtain information about the parameters of the model and hopefully to confirm or to exclude the model in a near future by exploiting the synergies with direct searches for the new scalars at LHC.
In this section we follow the notation of [8], where all the relevant formulae can be found, and update the new bounds. We collect the relevant tree-level lepton flavour violating  [18] for the different couplings. There seems to be a 2σ discrepancy in G exp τ /G exp e , which we interpret as a bound. If confirmed and interpreted within the ZB model, one obtains that |f µτ | 2 − |f eµ | 2 = 0.05 (m h /TeV) 2 . As we will see in section IV, for NH spectrum f eµ ∼ f µτ /2, therefore one needs m h ∼ 4 f µτ TeV, which is easily achieved. For IH spectrum, however, f µτ ∼ 0.2f eµ (f µτ ∼ (0.15 − 0.3) f eµ if we vary the angles in their 3σ range), and therefore, if this measurement is confirmed, the IH scheme in the ZB model  Finally, one-loop level lepton flavour violating constraints coming from − a → − b γ decays 2 and anomalous magnetic moments of electron and muon are collected in table III, including the recent limit on BR(µ → eγ) from the MEG Collaboration [23].
Given that lepton number is not conserved, another interesting low energy process that could arise in the ZB model is neutrinoless double beta decay (0ν2β). However, since the singly and doubly charged scalars do not couple to hadrons and are singlet under the weak SU(2) (therefore, do not couple to W gauge bosons), the 0ν2β rate is dominated by the Majorana neutrino exchange [24] and it is proportional to the |(M ν ) ee | 2 matrix element. In the NH case, Using neutrino oscillation data, one obtains 0.001 < ∼ eV|(M N H ν ) ee | < ∼ 0.004 eV and there-fore it is outside the reach of present and near future 0ν2β decay experiments.
C. Non-standard interactions.
The heavy scalars of the ZB model induce non-standard lepton interactions at tree level, which have been thoroughly analyzed in [9]. In particular, by integrating out the singly charged scalar h + , the following dimension-6 operators are generated: where refer to the charged leptons and the standard NSI parameters ρσ αβ are given by Regarding neutrino propagation in matter, the relevant NSI parameters are m αβ = ee αβ . Since the couplings f σβ are antisymmetric, in the ZB model only m µτ , m µµ and m τ τ are non zero. NSI can also affect the neutrino production in a neutrino factory, via the processes µ → eν β ν α . Source effects in the ν µ → ν τ and ν e → ν τ channels are produced by the NSI respectively. Notice that m µτ = − s * µτ , since both NSI parameters are related to the couplings f eµ and f eτ .
As we discuss in section V, the ratios of Yukawa couplings f eµ /f µτ and f eτ /f µτ are entirely determined by the neutrino mixing angles and Dirac phase of the PMNS matrix U -see eqs. (37) and (38) -, so the impact of the improved bounds on BR(µ → eγ) can be easily estimated: given that the limit is now ∼ 0.05 times smaller than in the study of ref. [9], and the contribution of the singly charged scalar h + to BR(µ → eγ) depends on |f * eτ f µτ | 2 , the current constraints on |f αβ | are roughly a factor 2 tighter than before. Therefore, since the strength of the NSI depends on ρσ αβ ∝ f σβ f * ρα , generically we expect that the allowed size of the NSI is reduced by a factor ∼ 1/4. According to [9] 3 , this implies that in the most favorable case of IH neutrino mass spectrum, s eτ and s µτ are in the range 3 × (10 −5 − 10 −4 ), which is in a range difficult to probe, but it might be in a future neutrino factory with a ν τ near detector [25].

D. Bounds on the masses of the charged scalars.
Regarding limits on singly-charged bosons decaying to leptons, the best limit still comes from LEP II, m h > 100 GeV.
ATLAS and CMS have placed limits on doubly-charged boson masses from searches of di-lepton final states, using data samples corresponding to √ s = 7 TeV with an integrated luminosity of 4.7 fb −1 and 4.9 fb −1 , respectively [26,27]. The authors of [28] show that, with current data at 8 TeV and 20 fb −1 , all the bounds are expected to become about ∼ 100 GeV more stringent if no significant signal is seen. The main production mechanisms of doubly-charged bosons at hadron colliders are pair production via an s-channel exchange of a photon or a Z-boson, and associated production with a charged boson via the exchange of a W-boson (see [29,30] for a general analysis of the production and detection at LHC of doubly charged scalars belonging to different electroweak representations). In the Zee-Babu model, the associated production is absent, because the new scalars are SU (2) L singlets.
The ATLAS analysis [26] focuses on the ee, µµ, eµ channels and assumes that the rest of the channels can make up to 90% of the total decays. Then, the limits for the Zee-Babu 100% of the times, i.e., BR(k ++ → ) = 1, and also considering several benchmark points with different branching ratios.
Note that whenever the branching ratio to τ τ is less than 30% (see table I and VI of [27]), the bounds are ∼ 280 GeV, provided that there is a significant fraction of decays into light leptons (ee, µµ, eµ).
In the Zee-Babu model the decay width of k ±± into same sign leptons is given by Since the g ab couplings are free parameters, the BRs of the different decay modes are a priori unknown, so we can not apply directly these bounds. As we will see in the numerical analysis, section V, once neutrino oscillation data and low energy constraints are taken into account, the branching ratio to τ τ is very small in the Zee-Babu model, less than about 1%.
Then, a conservative limit is m k > 220 GeV.
Moreover, in the ZB model for m k > 2 m h > 200 GeV, it can happen that the doubly charged scalar decays predominantly into hh, which can easily escape detection. This way the constraints from dilepton searches could be evaded. The relevant decay width is given by Then, even for g ab ∼ 1, for m h = 100 GeV and m k = 200 GeV, we have that Γ(k→hh) Γ(k→ ) ≥ 1 for µ ≥ m k , which is still natural as long it is not very large. Thus, we take m k ≥ 200 GeV in the numerical analysis.
E. Stability of the potential.
In this section we consider further constraints on the ZB model parameter space coming from vacuum stability conditions. The Hamiltonian in quantum mechanics has to be bounded from below, this requires that the quartic part of the scalar potential in eq. (2) should be positive for all values of the fields and for all scales. Then, if two of the fields H, k or h vanish one immediately finds 4 : Moreover the positivity of the potential whenever one of the scalar fields H, h, k is zero where we have defined Eq. (25) constrains only negative mixed couplings, λ xH , λ hk (x = h, k), since for positive ones the potential is definite positive and only the perturbativity limit, λ xH , λ hk < ∼ 4π applies.
Finally, if at least two of the mixed couplings are negative, there is an extra constraint, which can be written as: We have checked that the above conditions, eqs. (24,25,27), are equivalent to the ones derived in [32] for the Zee model, but they differ from the ones used in [12] for the ZB model, which seem not to be symmetric under the exchange of α, β, γ, as they should. Our constraints also agree with the results obtained by using copositive criteria (see for instance ref. [33]). 4 We do not consider the possibility of zero couplings, which can only appear at very specific scales.

III. H → γγ AND H → Zγ
It remains an open question whether the 125 GeV Higgs boson discovered by ATLAS [10] and CMS [11] where we have defined δR(m x , λ xH ) for the scalar x with mass m x and coupling to the Higgs λ xH as: with τ i ≡ A 0 (τ h,k ) > 0, therefore in order to obtain a constructive interference we need to consider negative couplings λ hH , λ kH .
As discussed in sec. II E, stability of the potential imposes that 2 √ λ H λ x + λ xH > 0, for Since M H ∼ 125 GeV fixes the value of the Higgs self-coupling to λ H ∼ 0.13, it is immediately apparent that large and negative λ xH couplings are going to be in conflict with stability of the potential, unless we push λ x close to the naive perturbative limit (λ x < 4π), for which −3 < ∼ λ hH , λ kH . Notice that this fact is not a special feature of the ZB model, but a generic problem of any scenario in which the enhancement of the Higgs diphoton decay rate is due to a virtual charged scalar.
We can consider three different cases: • If m k m h , For the same masses and couplings of both singlets, the doubly charged produces a larger enhancement/suppression than the singly-charged, due to its greater charge.
The largest enhancement can happen when both charged scalars are about the same mass and these masses are low enough. We show in fig. 2 the prediction of the ratio R γγ when the doubly charged scalar k dominates, for different values of the coupling with the Higgs, λ kH . Both an enhancement (as seen by ATLAS [35]) or a suppression (as seen by CMS [36]), can be accommodated. In fact, deviations from the SM value are expected, i.e., R γγ = 1, in particular for below the TeV scale singlets and sizeable scalar couplings. Of course, even for light singlets it is possible that R γγ ≈ 1, either because the relevant scalar couplings are tiny or due to a cancellation between the contributions of the singly charged and the doubly charged scalars.
In principle, the enhancement R γγ induced by a singly charged scalar h of similar mass and coupling to the Higgs λ hH ∼ λ kH is smaller; however since the lower limit on m h from LEP II direct searches is weaker m h > 100 GeV, as discussed in the previous section, and the largest contribution occurs for lower masses, the resulting values of R γγ for the allowed range of m h are comparable to the doubly charged case.
In summary, to obtain R γγ ∼ 1.5 we need m h < 350 GeV and/or m k < 500 GeV. As it will be shown in the numerical analysis section, these scalar masses are in tension with describing neutrino oscillation data and being compatible with current low-energy bounds in the ZB model if naturality is required at the level of κ = 1, especially for the NH Figure 2: R γγ in the presence of a doubly charged particle. Both an enhancement (as seen by ATLAS [35]) or a suppression (as seen by CMS [36]), can be accommodated. For the same masses and couplings, the singly-charged produces a smaller enhancement/suppression than the doublycharged, due to its smaller charge. spectrum. Moreover, the large negative values of the couplings λ xH ∼ −2 required to obtain such enhancement also induce vacuum instability of the ZB scalar potential, unless the corresponding coupling λ x is close to the perturbative limit, λ x ∼ 8.
decay rate in the ZB model with respect to the SM one is: where A Zγ SM is the SM H → Zγ decay amplitude, with λ i ≡ In fact, to have an enhancement in the H → γγ channel, we need negative couplings of the singlets with the Higgs, which in turn implies that the H → Zγ channel is reduced with respect to SM prediction, as can be seen in fig. 3.
We show in fig. 4 the contours of R γγ = 1.55 (0.78) motivated by the experimental results of ATLAS and CMS [35,36] in terms of singly and doubly charged masses, for various negative (positive) couplings.

IV. ANALYTICAL ESTIMATES
In this section we give some order of magnitude estimates of the free parameters in the ZB model, which complement and help to understand our full numerical analysis. In particular, we want to estimate to which extent light charged scalar masses, for instance like those required to fit an enhanced Higgs diphoton decay rate or to have a chance of being discovered at the LHC, are consistent with neutrino oscillation data and low-energy constraints.
As discussed in sec. II, with respect to the SM the ZB model has 17 extra parameters relevant for neutrino masses (9 moduli and 5 phases from the Yukawa couplings f, g, and 3 mass parameters from the charged scalar sector, m h , m k and µ), plus 5 quartic couplings in the scalar potential. However, some of the free parameters can be traded by the measured neutrino masses and mixings, ensuring in this way that the experimental data is reproduced and reducing the number of free variables as follows.
Since det f = 0, there is an eigenvector a = (f µτ , −f eτ , f eµ ) which corresponds to the zero eigenvalue, f a = 0 [6]. Then, by exploiting the fact that a is also an eigenvector of M ν , we which leads to three equations, one of which is trivially satisfied because one element of D ν is zero. The other two equations allow to write the ratios of Yukawa couplings f ij in terms of the neutrino mixing angles and Dirac phase as follows: f eτ f µτ = tan θ 12 cos θ 23 cos θ 13 + tan θ 13 sin θ 23 e −iδ , in the NH case, and for IH spectrum. Therefore, we choose f µτ as a free, real, parameter and obtain (complex) f eµ and f eτ from the above equations. Notice that the measured values, s 2 12 ∼ 0.3, s 2 23 ∼ 0.4 and s 2 13 ∼ 0.02 imply that, for NH, the first term on the right-hand side of eqs. (37) dominates and leads to f eµ ∼ f µτ /2 ∼ f eτ . Conversely, for IH it is clear that f eτ /f eµ = − tan θ 23 ∼ −1 and |f eµ /f µτ | ∼ |f eτ /f µτ | ∼ 4. Of course, to explain such fine-tuned relations of Yukawa couplings a complete theory of flavour would be needed, which is beyond the scope of this work.
Regarding the Yukawa couplings g, we keep g ee , g eµ and g eτ as free complex parameters and fix the remaining ones (g µµ , g µτ , g τ τ ) by imposing the equality of the three elements m 22 , m 23 and m 33 of the neutrino mass matrix M ν , written in terms of the parameters of the ZB model in eq. (13), and in terms of the masses and mixings measured in neutrino oscillation experiments in eq. (14), i.e., where we have defined ω ab ≡ m a g * ab m b , and ζ = µ 48π 2 M 2Ĩ (r), being r the ratio of the scalar masses, r ≡ m 2 k /m 2 h . Because of the hierarchy among the charged lepton masses, m e m µ , m τ , it is natural to assume that ω ee , ω eµ , ω eτ ω µµ , ω µτ , ω τ τ . Within the approximation ω ea = 0, the equation (39) for neutrino masses is simplified, and we can easily estimate the ranges of parameters consistent with neutrino oscillation data. Thus in this section we neglect them, although we keep all ω ab in the full numerical analysis 5 We then have From the large atmospheric angle we expect which leads to a definite hierarchy among the corresponding g ab couplings: g τ τ : g µτ : g µµ ∼ m 2 µ /m 2 τ : m µ /m τ : 1.
It is now convenient to write the mass matrix elements m ij in terms of the neutrino masses and mixings. In the normal hierarchy case this gives and therefore the hierarchy of couplings in eq. (42) is also obtained. However in the IH case there is a strong cancellation for Majorana and Dirac phases close to π, so we can obtain smaller values of ω ab . In particular, for φ = δ = π and the best fit values of the masses and mixing angles we find which allows for a smaller g µµ and as a consequence a lighter m k still consistent with the experimental limits. Therefore, although in the following analytic approximations we assume the hierarchy of couplings in eq. (42), one has to keep in mind that a larger parameter space is expected to be allowed when φ δ π. Indeed we will confirm in the full numerical analysis of section V that this region is specially favoured for light m k . Now we can estimate the lowest scalar masses able to reproduce current neutrino data.
Using the neutrino mass equation we can write 6 The upper bound on τ → 3µ decay implies that |g µµ | < ∼ 0.4(m k /TeV), while the new MEG limits on µ → eγ lead to |f µτ | 2 < ∼ 1.3 · 10 −3 (m h /TeV) 2 , where ≡ |f eτ /f µτ | ∼ 1/2 (4) for NH (IH). Substituting these constraints in eq. (48) we obtain m 33 0.05 eV which can be translated into a lower bound on the scalar masses. Using that m 33 ∼ 0.025 eV from neutrino oscillation data, if m h > m k then µ ≤ κ m k andĨ(r) ∼ 1, so eq. (49) implies that On the contrary, if m h < m k , we find From the above results 7 , we conclude that: 1. It is easier to reconcile an enhanced Higgs diphoton decay rate with neutrino oscillation data if the former is due to the doubly charged scalar loop contribution, since the lower bounds from neutrino masses are similar, while the BR(H → γγ) can be accounted for by a heavier m k . Moreover, if the enhancement is due to a light m h , then m k can not be very heavy, because otherwise neutrino masses are too small.
2. For a NH neutrino mass spectrum, it is possible to fit simultaneously neutrino oscillation data, lepton flavour violation constraints and an enhanced BR(H → γγ) only if the trilinear coupling µ is large, namely κ > ∼ 4(10) for min(m h , m k ) = 500 (300) GeV, respectively.
3. In general, the case of IH neutrino masses is in conflict with an enhanced Higgs diphoton rate unless κ ∼ O (30). However if we take into account the strong cancellations in ω µµ when φ δ π, and allow for a smaller m 33 ∼ 0.003 eV, it is also possible to fit all data with κ ∼ 4. 7 Our limits in the IH case differ from those in [7]. We traced this difference to the fact that in the estimates of [7] the perturbativity bound |g µµ | < 1 is imposed, but for low masses, m k < 2 TeV, such bound is always satisfied, and the relevant bound is |g µµ | < ∼ 0.4(m k /TeV), which depends on m k and changes the scaling with , leading to a weaker lower bound on the charged scalar masses in our case. We thank Martin Hirsch for discussions about this point.

V. NUMERICAL ANALYSIS
In order to explore exhaustively the highly multi-dimensional parameter space of the ZB model, naive grid scans are completely inappropriate, the method of choice is resorting to Monte Carlo driven Markov Chains (MCMC) that incorporate all the current experimental information described in precedence. As parameters we will use {s 2 ij , ∆ A , ∆ S , δ, φ, f µτ , m h , m k , µ, g ee , g eµ , g eτ }, and we allow them to vary within the ranges showed in table IV.
Had we tried to use our MCMC to obtain a posteriori probability distribution functions with a canonical Bayesian meaning, the choice of priors would have had a significant role.
Nevertheless, since our aim is to explore where in parameter space could the ZB model adequately reproduce experimental data without weighting in the available parameter space volume (that is, the "metric" in parameter space given by the priors), we will represent instead profiles of highest likelihood (equivalently profiles of minimal χ 2 ≡ −2 ln L with L the likelihood) which, on the contrary, can be interpreted on a frequentist basis. This is not a choice that we make because of the merits or demerits of either statistical school: our goal remains to understand if and where the ZB "works well", i.e. could fit experimental data. The interpretation of the results/plots will be clear: they show the regions where the model is in agreement with data without regard to their size when the remaining information (parameters and observables) is marginalized over 8 . In this case, exploring the parameter space in a uniform, logarithmic or other manner, in some given parameter will not affect our results (only the computational efficiency required to reach them will be, of course, affected).
For the modelling of experimental data we typically resort to individual Gaussian likelihoods for measured quantities. Bounds are implemented through smooth likelihood functions that include, piecewise, a constant and a Gaussian-like behaviour. For the sake of clarity: if the experimental bound for a given observable O is B O [90%CL] at 90% CL (1.64σ in one dimension), the χ 2 contribution associated to the model prediction O th for this observ-able is  To compare our results with the analysis presented a few years ago by some of us [8] some remarks are in order: first, here we have updated the experimental input on LFV and neutrino oscillation parameters, especially important are the determination of sin θ 13 and the new limits on µ → eγ. This tends to reduce the allowed regions but not dramatically.
Second, although the scanning of parameters is performed like in [8], we have chosen here to present results in terms of profiles of highest likelihood, which are insensitive to the volume of the parameter space and the priors used to scan it. This allows us to explore regions where parameters are fine tuned (after all, Yukawa couplings always require a certain degree of fine tuning). This is important since, as we have seen, the model is highly constrained at present and less conservative assumptions could exclude it before time, at least in the region of low masses. Moreover, we focus only on the region of masses with phenomenological interest (m h,k < 2 TeV) precisely to explore better the region of low masses.
In fig. 5 we depict the points allowed by neutrino oscillation data and all low energy constraints in the plane (m h , m k ) for the two mass orderings (NH and IH) and different values of the fine-tuning parameter in eq. (9) (κ = 1 darker, κ = 5 dark, κ = 4π light). The results of the numerical analysis imply that in general the indirect lower bounds on m h and m k from neutrino oscillation data and low energy constraints are stronger than the current limits from direct searches, except when cancellations occur for δ, φ ∼ π, especially in the IH case, and/or when naturality assumptions on µ are relaxed, allowing for κ = 4π. In table V we summarize the lower bounds on the scalar masses obtained for the three values of the naturality parameter κ, and two illustrative values of the Dirac phase, δ = 0, π.
The correlation between the CP phase δ of the neutrino mixing matrix and the scalar masses is illustrated in fig. 6, where we plot δ versus the doubly charged scalar mass, m k . 9 9 The correlation of δ with m h is entirely analogous.  Such correlation is especially relevant in the IH case, where scalar masses lower than ∼ 1 TeV are only allowed if δ ∼ π. A similar correlation with the phase φ was already found in [8] for IH spectrum, so we do not show it here.
Regarding the singly charged scalar h ± , the width of its decay modes (eν, µν, τ ν) is fixed by the f ia couplings to leptons (see for instance [7,8] for the relevant formulae). Therefore, after the measurement of θ 13 , present neutrino oscillation data determine completely the BRs of h from eqs. (37) and (38), up to a residual dependence on the CP phase δ in the case of NH spectrum. In this case, a very precise measurement of the branchings in the µν or τ ν channels (probably in a next generation collider) will predict the CP phase δ, and viceversa.
We show the ranges attainable by the different BRs in fig. 7, as a function of δ, splitting the two currently allowed octants of θ 23 . The most significant change between octants is the interchange of the µν and τ ν for the IH case. Clearly, the best option to discriminate between hierarchies is the eν channel.
An important point of the ZB model is that the doubly charged scalar can decay to two singly charged scalars, which are difficult to detect at the LHC. However, in fig. 5 we see that for a NH neutrino mass spectrum m h > 200 GeV, and the channel k → hh is closed for m k < 400 GeV. Therefore, present bounds on m k from dilepton searches at LHC discussed in II D apply. For the IH case, the k → hh channel is always open and can be dominant, unless κ = 1, for which we obtain that it is closed in the region m k < 440 GeV. Thus in general current direct bounds from LHC are weaker.
Let us now turn to the g ab couplings. We find always g τ τ g µτ , both for the NH and IH cases, in agreement with the analytic estimates in eq. (42); however the expected ratio g µµ /g µτ ∼ m τ /m µ is only fulfilled for the NH spectrum, since in the IH case large cancellations when the phases of the PMNS matrix U are δ ∼ φ ∼ π lead to smaller g µµ ∼ g µτ . This can be seen in fig. 8, where we show the ratios g τ τ /g µτ and g µµ /g µτ obtained in the numerical simulation as a function of δ, together with the expectation based on the analytic approximations, which is just a constant fixed by the charged lepton masses (red horizontal line) 10 .
To set the absolute scale of the couplings we present in fig. 9 the value of the largest couplings against m k , namely g µµ in the NH case, and g µτ in the IH case. We see that in both cases the couplings are always in the range from 10 −2 to 1 and therefore they tend to dominate the decays of the k ++ . 10 In the NH case there can also be cancellations with the g eτ terms, which have been neglected in eq. (43), that would allow much smaller values of g τ τ and g µτ , but those only occur for κ = 4π and in a tiny region of the parameter space.  Regarding the couplings g ea , which are not determined by the neutrino mass matrix, bounds from LFV charged lepton decays strongly constrain g eτ and g eµ to be less than O(0.01), while g ee can be larger, O(1). The constraint on |g ee g eµ | from µ → 3e implies that |g ee g eµ | < 2.3 × 10 −5 (m k /TeV) 2 and it is illustrated in fig. 10.
Since the widths of the k ±± leptonic decay modes are directly related to these couplings, from the above results we can readily infer the corresponding BRs. We find that the probability of k → eµ, eτ, τ τ is always negligible (even in the IH case, g eµ can be at most 0.1 and only when δ = π). For m k < ∼ 400 GeV, and NH neutrino spectrum, BR(k → ee) + BR(k → µµ) ∼ 1, since the k → hh decay channel is closed; therefore k ±± can not evade current LHC bounds on doubly-charged scalar searches and the limit m k > 310 GeV applies (400 GeV if no signal is found at 8 TeV with 20 fb −1 [28]

VI. CONCLUSIONS
We have analyzed the ZB model in the light of recent data: the measured neutrino mixing angle θ 13 , limits from the rare decay µ → eγ and LHC results. Although the model contains many free parameters, neutrino oscillation data and low energy constraints are powerful enough to rule out sizeable regions of the parameter space. A large source of uncertainty comes from the mass scale of the new physics, which is unknown. Since we are interested on possible signatures at the LHC, we present results for the masses of the extra scalar fields below 2 TeV. Previous analyses [7,8] have shown that larger mass scales are always allowed, given the absence of significant deviations from the SM besides neutrino masses.
Even within this reduced scenario, there is still a free mass parameter, the trilinear coupling between the charged scalars, µ, which remains mainly unconstrained. Naturality arguments together with perturbativity and vacuum stability bounds, indicate that µ can not be much larger than the physical scalar masses, m k , m h , but it is not possible to determine a precise theoretical limit. Because the neutrino masses depend linearly on the parameter µ, the ability of the model to accommodate all present data is quite sensitive to the upper limit allowed for it, so we have considered three limiting values, µ < κ min(m k , m h ), with κ = 1, 5, 4π. Within the above ranges for the mass parameters of the ZB model, we have performed an exhaustive numerical analysis using Monte Carlo Markov Chains (MCMC), incorporating all the current experimental information available, both for NH and IH neutrino masses. The results of the analysis are presented in sec. V and summarized in figs. 5 -10.
We have addressed the possibility that the slight excess in the Higgs diphoton decay observed by the ATLAS collaboration is due to virtual loops of the extra charged scalars of the ZB model, h ± and k ±± . Note that in the Zee-Babu model, as the new particles are singlets, there is a negative correlation between H → γγ and H → γZ. Although a similar study has been performed in [12], it was limited to the scalar sector parameters of the model, and neutrino data, which we find crucial to determine the allowed charged scalar masses, was not included in the analysis. In agreement with [12], we find that in order to accommodate an enhanced H → γγ decay rate, large and negative λ hH , λ kH couplings are needed, together with light scalar masses m h < 350 GeV, m k < 500 GeV. Such couplings are in conflict with the stability of the potential, unless the self-couplings λ h,k are pushed close to the naive perturbative limit, ∼ 4π. As a consequence, even if vacuum stability and perturbativity constraints are satisfied at the electroweak scale, RGE running leads to non-perturbative couplings at scales < ∼ O (100 TeV).
When neutrino data and low energy constraints are taken into account, we still find regions of the parameter space in which such enhancement is compatible with all current experimental data; in particular, it seems easier if the enhancement is due to the doublycharged scalar loop contribution. As can be seen in fig. 5, in the NH case, the trilinear coupling µ should be near its upper limit, while in the IH case lower masses can be achieved in the region δ ∼ φ ∼ π due to cancellations.
Regarding LHC bounds on the doubly-charged scalar mass, they are largely dependent on the BRs of the k ±± decay modes, namely same sign leptons ± a ± b and h ± h ± . The leptonic decay widths are controlled by the g ab couplings to the right-handed leptons, which are in principle unknown. By imposing that the measured neutrino mass matrix is reproduced, within the approximation m e = 0 one obtains analytically that g τ τ : g µτ : g µµ ∼ m 2 µ /m 2 τ : m µ /m τ : 1, while there is no information on the g ea couplings. Our numerical analysis confirms the above ratio of couplings in the case of NH, but for the IH spectrum there can be large cancellations if the PMNS matrix phases δ, φ are close to π, leading to g τ τ g µτ ∼ g µµ .
Moreover, in NH, if m k < 400 GeV for κ = 4π (m k < 600 GeV if κ = 5), m h < m k /2 is ruled out, therefore the decay channel k → hh is kinematically closed and the LHC bounds from doubly-charged scalar searches can not be evaded. In IH, however, for δ ∼ φ ∼ π the k → hh channel is open unless κ = 1, while if δ = π, indirect bounds on m k set a much stronger constraint than direct LHC searches.
As a consequence, if the light neutrino spectrum is NH, k decays mainly to ee, µµ, and the current bound from LHC is m k > 310 GeV, while if the spectrum is IH, k may also decay to µτ and hh, so the present bound is weaker, about 200 GeV. Were a doubly-charged boson discovered at LHC, the measurement of its leptonic BRs could rule out the ZB model, or predict a definite neutrino mass spectrum. Conversely, if a CP phase δ = π is measured in future neutrino oscillation experiments together with an IH spectrum, the mass of the charged scalars of the ZB model will be pushed up well outside the LHC reach.
Note: During the final stages of this work we became aware of [40], where an analysis of the Zee-Babu model was performed. Our bounds on the scalar masses are comparable to theirs taking into account the slightly different procedures, in particular that they fix the neutrino oscillation parameters to their best fit values and we allow them to vary in their two sigma range. While in our work we focus on prospects for the LHC, in [40] the possibility of detecting the doubly charged singlet in a future linear collider is studied.