Radiative Model of Neutrino Mass with Neutrino Interacting MeV Dark Matter

We consider the radiative generation of neutrino mass through the interactions of neutrinos with MeV dark matter. We construct a realistic renormalizable model with one scalar doublet and one complex singlet together with three light singlet Majorana fermions, all transforming under a dark U(1) symmetry which breaks softly to Z_2. We study in detail the scalar sector which supports this specific scenario and its rich phenomenology.

The nature of dark matter (DM) is one of the most disputed topics in cosmology. Until one (or two) decade(s) ago, only a few candidates prevailed in the literature, among which were neutralinos (a thermal, cold, DM species) and axions (also cold DM but non-thermal).
Astrophysical and cosmological anomalies since in the last 10-15 years however led many authors to study more exotics scenarios, such as light DM, leptophilic DM, sterile neutrinos [1,2,3,4,5]. So far most DM studies have focused on either almost massless particles (axions), keV particles (sterile neutrinos) or GeV to TeV DM candidates (as provided by supersymmetry and Kaluza-Klein theories 1 ) but the range between a few keV and GeV has been somewhat disregarded.
In cosmology, both keV and GeV-TeV DM candidates are generally assumed to be collisionless, even though their annihilations or decay are invoked to explain the observed DM abundance. About a decade ago, it was pointed out that -even weak-strength -DM interactions could erase the DM primordial interactions and should not be disregarded when the DM is relatively light (a few MeV) and coupled to neutrinos or photons [7,8]. Indeed the damping of the DM primordial fluctuations has two origins, as shown in these references: one is the collisional damping, which suppresses the matter fluctuations until the DM is kinetically decoupled from any other species; this is analogous to the Silk damping. The other source is the DM free-streaming which erases fluctuations that have not been erased yet by the DM collisions.
The resulting linear matter power spectrum associated with light DM candidates coupled to radiation features damped oscillations in addition to an exponential cut-off [9,10,11]. This makes these scenarios interesting alternatives to vanilla CDM and Warm DM candidates. 1 For a review see Ref. [6].
There has been much interest in the DM-neutrino coupling since these first studies but with the twist of DM self-interactions [12,13,14,15]. However, as shown in Refs. [16,17,18], a sole DM-neutrino coupling can solve the missing satellite (which is a deficit of dwarf galaxy haloes in Milky Way-like DM haloes) and the too big to fail problems when the DM elastic scattering off neutrinos is of the order of which is the required value to explain the observed DM abundance in thermal DM scenarios.
The correspondence between Eq. (1) and Eq. (2) thus suggests that current cosmological problems could be related to the current DM abundance.
Even more puzzling is the possibility to explain in some specific models [19,20,21] the existence of small neutrino masses in the presence of such a DM-neutrino coupling.
It is therefore tempting to assume that there exists a framework in which DM-neutrino interactions can explain simultaneously the missing satellite and too big to fail problems, the existence of small neutrinos masses and the observed DM abundance.
In this paper we construct such a framework. We envision a fundamental Yukawa coupling of the formN R ν L ζ 2 where the dark matter candidate, here referred to as N , is a Majorana fermion and both the fermion N and the scalar ζ 2 are light, with masses of order a few MeV. In Section 2, we review the elastic scattering cross section among the neutrino and DM and the related process of DM annihilation into neutrino pair based on this Yukawa coupling. To support this specific scenario, we study an extension of the standard model with one additional scalar doublet and one additional complex singlet, both of which transform nontrivially under a dark global U (1) D that is softly broken into a discrete Z 2 (Section 3). We show how realistic neutrino masses may be obtained with a scalar mass spectrum including the light ζ 2 without conflicting with present data at the Large Hadron Collider (LHC) (Sections 3 and 4). We examine also in detail the scalar sector and obtain theoretical and phenomenological constraints on its parameter space (Section 4). Numerical results are presented in Section 5 and conclusion in Section 6. Some useful formulas are collected in the Appendix.

Elastic scattering and annihilation cross sections
In (thermal) scenarios where DM can scatter off neutrinos, the collisional damping scale is determined by the integral where a is the scale factor, ρ ν the neutrino energy density, Γ ν the neutrino interaction rate, v the neutrino velocity, (ρ + p) tot is the sum over the energy density and pressure of all the species coupled to the DM while DM still interacts with neutrinos (which includes the DM itself). This length is directly proportional to the neutrino density and velocity (which is equal to c if one assumes that the DM kinetic decoupling from neutrinos happens well before neutrinos become non relativistic) and the neutrino kinetic decoupling time [7,8]. Its magnitude also depends on the period over which the DM is coupled to neutrinos; hence the integral over time, with t dec(DM−ν) (the DM decoupling time from neutrino) as upper bound.
The CMB and linear matter power spectra in the presence of such a DM-neutrino coupling can be easily predicted using the Boltzmann formalism [11,22]. Both agree with the damping length estimate obtained analytically using the above formula (in the absence of mixed damping). But the matter power spectrum ultimately sets the stronger constraint, namely if the cross section is independent of the temperature or if the cross section depends on the neutrino energy [22]. This confirms that a weak strength cross section can erase DM fluctuations at cosmologically relevant scales, if the DM is relatively light. The simplest elastic scattering process N ν → N ν that gives rise to such an effect relies on the exchange of a fermion (scalar) if the DM is a scalar (fermion). The cross section for a Majorana candidate coupled to neutrinos with a coupling g is given by the u and s-channels diagrams, leading to: in the absence of a close degeneracy between the mediator and DM masses. Here we also implicitly assume MeV DM, i.e. that DM is non relativistic at the DM-neutrino decoupling, which occurs slightly below a keV. The annihilation diagrams (t and u-channels) lead to the dominant contribution σv g 4 4 π m 2 where we again assume that there is not a strict degeneracy between the DM and mediator masses and neglect the neutrino mass.
Eqs. (6) and (7) cannot be satisfied simultaneously with the same values of the mass and couplings, unless the DM mass is slightly smaller than a few keV. Yet thermal keV annihilating DM particles into neutrinos are already ruled out [23,24,25,26], as they would change the number of relativistic degrees of freedom at nucleosynthesis and CMB time, by too large an amount. The only possibility for thermal DM candidates coupled to neutrinos is to have a mass above a few MeV [25].
In order to explain both the DM abundance and solve cosmological problems, one thus needs thus to get rid off the temperature dependence of the elastic scattering cross section.
This occurs if the mass splitting between N and ζ 2 are of the order of a few keV or below.
Indeed in this case the elastic scattering cross section reads while the annihilation cross section is given by Therefore, a scenario where the DM is of a few MeVs but the mediator is only slightly heavier than the DM by a few keVs can solve the missing satellite and too big to fail problems (in the absence of baryonic physics) and also explain the DM observed abundance.
Note that the presence of baryonic interactions could alter these values. Depending on the magnitude of the effect, one might either lose the above correspondence or be able to make a temperature dependent elastic scattering and temperature independent annihilation cross section compatible. Given that such studies do not exist yet, we will take the above numbers (see Eqs. (8) and (9)) at face value.
We now investigate whether such a scenario can also give rise to neutrino masses. Note that such a scenario predicts a slightly larger value of N eff than 3.046 and H 0 71 km/s/Mpc [22].

6
The simplest finite one-loop radiative model of neutrino mass through dark matter is the scotogenic model (from the Greek "scotos" meaning darkness) proposed in 2006 [27]. It assumes an exactly conserved Z 2 symmetry [28] and extends the standard model (SM) of particle interactions with the addition of one scalar doublet (η + , η 0 ) and three singlet Majorana fermions N 1,2,3 which are odd under Z 2 . Many aspects of its phenomenology have been studied in detail [29]. Whereas the masses of η and N are usually considered to be heavy, this mechanism also allows N to be light [30]. With the discovery [31,32]  Higgs interaction, it does not work with just the inert Higgs doublet η, nor with the addition of a real scalar singlet. However, as shown in this paper, it will work with the η doublet plus a complex singlet χ, which is naturally maintained with a dark U (1) D symmetry, softly broken to Z 2 so that N may have a Majorana mass and ν gets a one-loop radiative mass as shown in Fig. 1. If only Z 2 is used, then many more allowed terms appear in the Lagrangian, such as (Φ † Φ)χ 2 + H.c. on top of (Φ † Φ)χ * χ, which are unnecessary complications (and must indeed be small) in the subsequent discussion of the scalar sector.
The scalar potential is given by where the m 2 4 term breaks U (1) D softly to Z 2 . Note that the quartic term (Φ † η) 2 present in the original Z 2 model [27] is now forbidden. The one-loop mechanism for neutrino mass is depicted in Fig. 1, from which we can easily see that the Majorana mass term for N i breaks U (1) D to Z 2 . This diagram is similar to that required in a supersymmetric extension [33]. We note that the m 2 4 and Majorana mass terms are the only two terms in the model that softly break U (1) D into Z 2 . This is a well-known method in symmetry breaking, because soft breaking generates only finite corrections in the renormalizable terms of the Lagrangian and a well-known rationale for setting m 4 small against other mass parameters (because its absence enlarges the symmetry of the theory).
The model however remains renormalizable. These soft terms are analogous to the sfermion mass and the gaugino mass in MSSM which break supersymmetry softly. The origin of these softly breaking terms may be revealed only at a higher theory. In the case of MSSM, these soft terms could arise from supergravity. In our case, we will treat them as phenomenological terms, just like the superpartner mass terms in MSSM, without worrying about the higher theory. For an early application of this idea in neutrino physics, see for example [34].
The mass of the SM-like Higgs h (≡ φ R ) and charged Higgs η ± are given by: The neutral components of η and χ will mix through µ term of the potential. The mass-squared matrices spanning (η R,I , χ R,I ) are given by Let ζ 1R , ζ 2R , ζ 1I , ζ 2I be the mass eigenstates with masses m 1R , m 2R , m 1I , m 2I : ζ 1I = cos θ I η I − sin θ I χ I , ζ 2I = sin θ I η I + cos θ I χ I , then the neutrino mass matrix is given by where M k are the masses of N k . Note that in the limit of m 2 and θ R = θ I . Hence the neutrino mass would be zero.
where s = sin θ 0 and c = cos θ 0 which diagonalize the (η 0 , χ) mass-squared matrix in the absence of m 2 4 with eigenvalues m 2 10 and m 2 20 . The one-loop neutrino mass matrix is then of the form For m 10  In this scenario, N is dark matter with a mass of 3 MeV, ζ 2 has a mass of 4 MeV and interacts withν L N R with strength 0.1. This is thus a possible scenario for neutrino interacting MeV dark matter which obtains the correct relic abundance, as discussed in Section 2. Note that the ζ 2 mass splitting is small, i.e. 3 keV, and both ζ 2R and ζ 2I decay to νN .
If the cosmological missing-satellite problem and the too-big-to-fail problem are solved using elastic N ν scattering, then ζ 2R and ζ 2I should be both only a few keV above M , in which case Eq. (21) is not valid. Let m 2 2R = M 2 (1 + δ R ) and m 2 2I = M 2 (1 + δ I ), with δ R,I of order 10 −3 , then the radiative neutrino mass becomes (s 2 h 2 /32π 2 )M (δ I − δ R ), which is of order 0.1 eV for hs = 0.1 as desired.
However, for the small couplings implied by Eqs. (8) and (9), the induced neutrino mass is negligible. On the other hand, only N 1 needs to be light, whereas N 2,3 can be heavy and the usual acceptable neutrino masses are obtained. The important point of this study (to be justified in the subsequent sections) is that the mass of one scalar, i.e. ζ 2 (ζ 2R or ζ 2I ), can be of order MeV. The other two neutral scalars ζ 1R and ζ 1I can be heavy. We will assume ζ 1R and ζ 1I are heavier than m h /2 so that they do not contribute to the invisible Higgs width.

Phenomenology of the Scalar Sector: Theoretical and Experimental Constraints
The scalar potential Eq. (10) has altogether 12 parameters and 1 vacuum expectation value (vev) v. Two of them (m 2 1 and v) can be eliminated by the minimization condition and W gauge boson mass. At the LHC, both ATLAS and CMS experiments had performed several measurements of the newly discovered scalar particle in different channels. The combined measured mass performed by ATLAS and CMS collaborations based on the data GeV [35]. This measurement if interpreted as the SM Higgs boson allow us to fix λ 1 . We are then left with 10 independent parameters: In our numerical analysis presented in the next Section, these parameters are scanned in the confined domain that fulfill various theoretical and experimental constraints which are discussed below.

Unitarity Constraints
Our scalar potential is similar to the one in the 2 Higgs doublet model except augmented by a complex singlet field χ. We can carefully use the full set of unitarity constraints derived for the 2 Higgs doublet model in [36]. In Appendix A.1, we list the set of unitarity constraints that we will use. Some of the 2 → 2 scattering amplitudes have been modified to take into account the presence of χ. In summary, the requirement that the largest eigenvalues of all the partial wave matrices a 0 s for different channels to respect the unitarity constraints where the definitions of the eigenvalues (a ± , b ± , and so on in the above equation) in terms of the quartic couplings in the scalar potential can be found in Appendix A.1.

Vacuum Stability
A necessary condition for the stability of the vacuum comes from requiring the potential given in Eq. (10) to be bounded from below when the scalar fields become large in any direction of the field space. At large field values, the scalar potential is dominated by quartic couplings, the bounded from below constraints will depend only on the quartic couplings.
The constraints ensuring tree level vacuum stability are: • If λ 6 > 0 and λ 7 > 0, • If λ 6 < 0 or λ 7 < 0, in additional to the above constraints, we also have Additional constraints also related to the stability issue come from the requirement of the absence of tachyonic masses. They are Details of derivation of these constraints can be found in Appendix A.2.

Invisible Decay of the Higgs
Our neutrino model requires MeV warm dark matter particle which can be identified as the lightest Majorana neutrino state of N 1,2,3 . The SM Higgs h → N i N i can occur only through one loop. Hence its branching ratio is small and we will ignore this invisible mode in our analysis. On the other hand, due to the mixing of complex field χ with the inert doublet η, we have 2 light dark Higgses ζ 2R and ζ 2I , one is CP even and one is CP odd. These states are not stable since they can decay via χ D → N ν where χ D is the lighter state of ζ 2R and ζ 2I and N is the DM, the lightest of N 1,2,3 . Thus the tree level decay h → χ D χ D → N N νν will be invisible. The SM Higgs couplings to these dark Higgses ζ 2R and ζ 2I are given in Table 1 in Appendix A.3.
The openings of one of the non-standard decays of the Higgs boson such as h → ζ i ζ j can modify the total width of the Higgs boson and can have significant impact on LHC results. The couplings of the SM Higgs to the neutral dark Higgses η R,I , χ R,I is given by As one can see from Eq. (13) If these conditions are fulfilled, we have M 2 R,I = vF. Once M 2 R,I are diagonalized by some orthogonal matrices, the coupling matrix F will be also diagonal in the mass eigenstate basis. Therefore the couplings hζ 2R ζ 2R and hζ 2I ζ 2I will be proportional to the mass of the dark Higgses and are therefore negligible for MeV dark Higgses.
The decay rate for h → ζ a ζ b can be found in Appendix A.3. In our case only h → ζ iR ζ jR and h → ζ iI ζ jI exists. Furthermore, provided that the alignment condition of M 2 R,I = vF can be satisfied, only h → ζ 2R ζ 2R and h → ζ 2I ζ 2I will contribute to the SM Higgs invisible width. The other diagonal decays h → ζ 1R ζ 1R and h → ζ 1I ζ 1I will be kinematically not accessible if we assume m ζ 1R and m ζ 1I are larger than m h /2.

Z Decay Width
The measurement of Z-boson total width Γ Z at LEP set stringent bounds on any extra contribution ∆Γ Z from new decay channels. In our case, Z can decay to ζ 2R ζ 2I through the mixing of the neutral component of inert doublet with the complex singlet.
The Zζ a ζ b couplings are listed in Table 1 sin θ R ≈ sin θ I ≤ 0.47 .

S and T Parameters
If the scale of new physics is much larger than the electroweak scale, virtual effect of the new particles in the loops are expected to contribute through vacuum polarization corrections to the electroweak precision observables. These corrections are known as oblique corrections and are parameterized by S, T and U parameters [41]. In our case, the inert Higgs doublet couples to the W and Z gauge bosons via the covariant derivative. Due to mixing effects, the complex singlet χ will couple to the weak gauge bosons as well. Therefore, both η and χ will contribute to S and T parameters which are very well constrained by electroweak precision data. Analytic formulas for ∆S and ∆T modified by the mixing angles as compared with the IHDM formulas are collected in Appendix A.5 for convenience. Thus our model will remain viable as long as ∆S and ∆T are compatible with the fitted values [42] which are given by: ∆S = 0.06 ± 0.09 and ∆T = 0.10 ± 0.07 .

LEP Limits
Due to the Z 2 symmetry, all interactions that involve the dark Higgses must contain a pair of them. The precise measurements of W and Z widths at LEP [40] can be used to set a limit on the mass of the inert Higgses. In order not to significantly modify the decay widths of W and Z, we request that the channels W ± → {ζ iR η ± , ζ iI η ± } and/or Z → {ζ iR ζ jI , η + η − } are kinematically closed. This leads to the following constraints: At e + e − colliders, the production mechanism for inert Higgs is while at hadron machines we have Because of Z 2 , the inert Higgs can not decay to SM fermions. Thus the LEPII and Tevatron searches for charged Higgs and neutral Higgs can not be applied to our model. The inert charged Higgs can decay via η ± → W ± ζ 2R , W ± ζ 2I or through cascade decay via η ± → Similarly, the neutral dark Higgses ζ 1R and ζ 1I can decay into Zζ 2I and Zζ 2R respectively, or through cascade decays like ζ 1R → In all cases the final states of the these production mechanisms both at lepton or hadron colliders would be multi-leptons or multi-jets, depending on the decay products of W ± and Z, plus missing energies carried by the dark Higgses.
To certain extent, the signatures for the inert charged or neutral Higgses would be similar to the supersymmetry searches for charginos and neutralinos at the e + e − or hadron colliders.
Detailed phenomenological implications of this model at the LHC are interesting to explore but it is beyond the scope of this present work.

Constraints from LHC
Both ATLAS and CMS experiments of the LHC run at 7 ⊕ 8 TeV confirmed the discovery of a scalar particle with mass around 125 GeV identified to be the Higgs field h in our model.
Both groups performed several measurements on this scalar particle couplings to the SM particles such as W + W − , ZZ, γγ and τ + τ − with 20 − 30% uncertainties, while for bb it suffers from larger uncertainty of 40 − 50%. Recently ATLAS [43] published an updated analysis of 7 ⊕ 8 TeV data in which the signal strengths 2.7 +4.6 −4.5 for h → γZ and −0.7 +3.7 −3.7 for h → µ + µ − were reported. Basically, all the LHC data collected so far indicates that the 125 GeV boson couplings to the SM particles are very much SM-like. One of the main tasks of the new LHC run at 13 TeV would be to improve all the aforementioned measurements and probe for new ones, such as h → γZ, µ + µ − and perhaps the trilinear self-coupling of the Higgs. It is expected that the new run of LHC will narrow down the uncertainties of hbb and hτ + τ − measurements to 10 − 13% and 6 − 8% respectively. In the future, if the high luminosity option for LHC (HL-LHC) is available, it can do much to ameliorate the uncertainties to 4 − 7% (hbb) and 2 − 5% (hτ + τ − ) [44]; while for the e + e − Linear Collider (LC), these uncertainties can be cut down further to 0.6% (hbb) and 1.3% (hτ + τ − ) [44,45].
While the tree level SM Higgs couplings to fermions and to weak gauge bosons in our model are identical to the SM one, the loop mediated processes such as h → γγ and h → γZ will receive additional contributions from inert charged Higgs loop that can either enhance or suppress their partial widths [46]. On the other hand, the invisible decay of the SM Higgs into dark Higgs pair is very much suppressed in our model. As a consequence the total width of the SM Higgs will be modified slightly through the additional charged Higgs contributions in the h → γγ and h → γZ modes. ATLAS and CMS collaborations usually present their results in terms of the so-called signal strengths. For a given production channel and a given decay mode of the SM Higgs, the signal strength is defined as where the Higgs mass is evaluated to be the same in both numerator and denominator.
In our analysis for the signal strengths, we will use the following ATLAS results [43]: • h → γγ: R γγ = 1.17 ± 0.27 We now present our numerical results with the implementation of all the theoretical and experimental constraints on the parameter space discussed in the previous Section. Let us classify the dimensionless parameters λ i in the scalar potential into two different sets according to the following two types of constraints: • First set of constraints includes the unitarity constraints in Eq. and R τ + τ − listed in the end of Section (4.2.5). We also require the masses of ζ 1R , ζ 1I and η ± to be heavier than 100 GeV. We refer this set of constraints as C 2 .
Since λ 1 is fixed by the SM Higgs mass, we scan over the other λ i ∈ P in the following range For the dimensional mass parameters in the scalar potential, m 2 1 is fixed by the SM Higgs mass, and m 2 2 and m 2 3 are fixed by Eq. (35) in order to suppress the invisible decay of the SM Higgs. For the m 2 4 and µ parameters, they will chosen in such a way to allow for MeV dark Higgses ζ 2R and ζ 2I . Recall that their masses are provided by the smaller eigenvalues of the two mass matrices in Eq. (13), where A, B and C are given by 2 To obtain very light dark Higgses, we fine tune A + C and (A − C) 2 + 4B 2 to be almost the same size. Define Assuming is small and dropping the 2 term, we have For given , A and C, the µ parameter is determined. If we further drop the term in Eq. (53), it would give an additional constraint on the sign of the product AC > 0: where Eq. (35) has been applied in the last equality. Thus λ 4 + λ 5 and λ 6 should have the same sign. On the other hand, neglecting m 2 4 in C, Eqs. (35), (49) and (51) imply Combining Eqs. (54) and (55), one can conclude λ 4 +λ 5 and λ 6 should be positive, at least for small and m 2 4 . Numerically, we set /GeV 2 = 10 −4 and m 2 4 /GeV 2 = 10 −5 in our analysis.
A systematic scan on λ i in the range defined in Eq. (47) indicated that λ 2 and λ 3 are not very much restricted by all the above constraints. In Fig. (2) we illustrate the allowed range while green points pass both C 1 and C 2 .
We have checked that all the red points in the left panel of Fig. (2) fall in the following domain:  From our previous discussion we demonstrated that under our assumptions λ 6 is positive.
It is clear from the plot at the right panel of Fig. (2) that when λ 6 and λ 7 are both positive, the C 1 constraint set does not restrain λ 6 and λ 7 too much. Even when both C 1 and C 2 are imposed, λ 7 is not very much constrained while the range of λ 6 has shrunk significantly.
This is due to the fact that λ 7 does not contribute to the masses of dark Higgses while λ 6 does. Figure 3: Scatter plot for R γγ in (λ 4 , m η ± ) plane (left) and in (λ 4 , λ 5 ) plane (right). All points pass both C 1 and C 2 sets.
In the left and right panels of Fig. (3) we present the scatter plots of the signal strength R γγ , represented by the color palettes located at the right sides of both panels, on the (λ 4 , m η ± ) and (λ 4 , λ 5 ) planes respectively. In these two plots, both C 1 and C 2 constraint sets are imposed. In our model, since the SM Higgs is produced exactly the same way as in the SM, the production cross sections in the numerator and denominator of R γγ cancel, and the signal strength is simply given by the ratio of branching fractions. Thus R γγ is independent of the LHC energy at Run 1 or 2.
As is well known the loop contributions in h → γγ is fully dominated by W ± with some subleading contribution from top quark which interferes destructively with the W ± .
As alluded earlier, h → γγ receives additional contribution from charged Higgs η ± in this model [46]. The coupling of the SM Higgs to the η ± pair is proportional to λ 4 . If λ 4 is negative (positive) then the η ± loop is constructively (destructively) interference with the W ± s, resulting in an enhanced (suppressed) h → γγ rate with respect to SM one. By comparing with the color palettes for R γγ on the right side of both panels of Fig. (3), it is evident that R γγ is enhanced for negative λ 4 but suppressed for positive λ 4 . Note that λ 4 is restricted only to a small range of negative λ 4 [−0.65, 0] which could enhance h → γγ rate with respect to SM. This range of negative λ 4 corresponds to λ 5 in the range [0.5, 3.7]. These two ranges for λ 4 and λ 5 imply that the charged Higgs η ± is in [100, 325] GeV range where R γγ > 1. It is clear from the left panel of Fig. (3) that larger η ± mass (and so as the two other neutral dark Higgses ζ 1R,1I ) say 500 GeV is also possible. But then the signal strength of R γγ would be very close to its SM value. Correlation between R γZ and R γγ with λ 4 shown in palette at the right.
In Fig. (4) we illustrate R γZ and its correlation with R γγ . In the left panel, we show R γZ as a function of λ 4 while scanning all other parameters. As in the R γγ case, R γZ is enhanced for negative λ 4 but suppressed for positive λ 4 with respect to SM. In the right panel we show the correlation between R γγ and R γZ . At the point λ 4 = 0, the charged Higgs contribution vanishes and both R γγ and R γZ reduces to their SM values. It is interesting to note that for R γγ > 1 we have 1 < R γZ < R γγ while for R γγ < 1 we have R γZ > R γγ . These predictions can be tested at LHC Run II.

Conclusion
In this work, we have presented a realistic renormalizable model with one

Appendix A A.1 Perturbative Unitarity Constraints
To constrain the scalar potential parameters, one can demand that tree-level unitarity is preserved in a variety of 2 → 2 scattering processes.
Since our model is a 2 Higgs doublet extended with a singlet field, we can use the same procedure developed in [36] to derive the unitarity constraints. According to [36], one computes the S matrix in the non-physical fields basis where the computation is much easier.
The crucial point is that the S matrix expressed in terms of the physical fields (i.e. the mass eigenstate fields) can be transformed into an S matrix for the non-physical fields by making a unitary transformation. The eigenvalues for the S matrix should be unchanged under such a unitary transformation.
The first submatrix M 1 , corresponding to scatterings whose initial and final states being one of the following combinations (w + 1 w − 2 ,w + 2 w − 1 , φ R η I , η R φ I , φ I η I , φ R η R ), is given by Its eigenvalues are determined as The second submatrix M 2 corresponds to scattering with initial and final states belonged to one of the following states (w where λ 45 = λ 4 + λ 5 . This matrix has 8 eigenvalues. Five of them are The other 3 eigenvalues b ± and s 2 are solutions of the following polynomial as eigenvalues as defined previously.
With the two singlet components, more states such as (φ R χ R,I , φ I χ R,I ), (η R χ R,I , η I χ R,I ), etc. can be constructed. But their corresponding scattering matrices will be diagonal and lead to either λ 6 or λ 7 as eigenvalues. These scattering states will not lead to any nontrivial constraints among λ i since they are required to be perturbative, namely |λ i | ≤ 4π for all i.
In our analysis we also include the following two body scattering processes among the 8 . This submatrix only lead to one additional constraint which is The others are duplicated with the previous cases.
With the two singlet components, we can also construct charged states like (χ R w + 1 , χ I w + 1 ) and (χ R η + , χ I η + ) which decouple from the previous charged scattering processes. Again, the scattering matrices in these cases are diagonal with eigenvalues λ 6 and λ 7 respectively.
All the eigenvalues shown in this appendix are required to satisfy the perturbative unitarity constraints as given by Eq. (23).

A.2 Vacuum Stability Constraints on Scalar Potential
At large field values the potential Eq. (10) is dominated only by the part containing the terms that are quartic in the fields The study of V quartic will thus be sufficient to obtain the main constraints from vacuum stability considerations.
One can rewrite the quartic terms using the new parameterization as where we have used x = sin φ. In this form, V quartic /r 4 is a second degree polynomial in x ∈ [0, 1]. One can show that V quartic /r 4 is positive if and only if 3 where we used y = cos 2 θ. The first condition (Eq. (78)) is nothing but a scalar potential without the singlet field. This condition will give us the boundedness from below for 2 Higgs doublet model. We note that A is a second degree polynomial in y. It is positive, if and only if The last condition of the above equation gives the following 2 conditions is positive if and only if a > 0, b > 0 and c > −2 √ ab.
In our analysis, we impose all the constraints derived in this appendix for the quartic couplings λ i . Table 1: General coupling coefficients of hζ a ζ b and Zζ a ζ b vertices.

A.4 Zζ a ζ b Couplings
From the covariant derivative (D µ η) † (D µ η), we have the following derivative couplings where the fields η R,I are related to the physical fields ζ 1R , ζ 2R , ζ 1I and ζ 2I as η R = cos θ R ζ 1R + sin θ R ζ 2R and η I = cos θ I ζ 1I + sin θ I ζ 2I . From the last term in Eq. (94), we get the vertex for Z( µ (k)) → ζ a (p)ζ b (p ) as +(g/2c θw )c ab (p − p ) µ with c ab defined in the last column of Table 1. The decay rate for Z → ζ a ζ b is given by A.5 Formulas for the ∆S and ∆T The analytic expressions for ∆S and ∆T can be given in terms of Passarino-Veltman functions which have been calculated using the software packages FormCalc [48] and Loop-Tools [49]. The SM expressions for S and T have been subtracted properly, we give hereafter only the extra contributions ∆S and ∆T . We take as reference point the Higgs mass m h = 125 GeV, m t = 173 GeV and assume ∆U = 0.
We have checked both analytically and numerically that ∆S and ∆T are UV finite and also independent of the renormalisation scale. In terms of the Passarino-Veltman functions A 0 and B 00 , they are given by  − cos 2 θ R B 00 [0, m 2 ζ 1R , m 2 η ± ] + sin 2 θ I sin 2 θ R B 00 [0, m 2 ζ 2R , m 2 ζ 2I ]) .