Flavor Dependent $U(1)$ Symmetric Zee Model with a Vector-like Lepton

We extend the Zee model by introducing a vector-like lepton doublet and a flavor dependent global $U(1)$ symmetry. Flavor changing neutral currents in the quark sector can be naturally forbidden at tree level due to the $U(1)$ symmetry, while sufficient amount of lepton flavor violation is provided to explain current neutrino oscillation data. In our model, additional sources of CP-violation appear in the lepton sector, but their contribution to electric dipole moments is much smaller than current experimental bounds due to the Yukawa structure constrained by the $U(1)$ symmetry. We find that there is a parameter region where the strongly first order electroweak phase transition can be realized, which is necessary for the successful scenario of the electroweak baryogenesis in addition to new CP-violating phases. In the benchmark points satisfying neutrino data, lepton flavor violation data and the strongly first order phase transition, we show that an additional CP-even Higgs boson $H$ mainly decays into a lighter CP-odd Higgs boson $A$, i.e., $H \to AZ$ or $H \to AA$ with a characteristic pattern of lepton flavor violating decays of $A$.


I. INTRODUCTION
Neutrino oscillations are the phenomena which indicate clear evidence for the necessity of physics beyond the standard model (SM). From various measurements, it has been known that masses of neutrinos have to be much smaller than those of charged fermions, e.g., O(10 12 ) times smaller than the top quark mass. This strongly suggests that the neutrino mass is generated by a different mechanism from that for charged fermions, i.e., Dirac masses via Yukawa interactions with a Higgs doublet field. In such context, Majorana masses for neutrinos can be a good candidate, which are effectively described by the dimension five Weinberg operator [1]. The question is then how the Weinberg operator can be written in terms of the ultra-violet physics. The simplest example has been known as the type-I seesaw mechanism [2][3][4], where only right-handed neutrinos are sufficient to be added to the SM. This, however, requires huge Majorana masses of the right-handed neutrinos, typically order O(10 14 ) GeV assuming order one Dirac Yukawa couplings, so that direct detections for such heavy particles are quite challenging.
As an alternative direction, we can consider a scenario where neutrino masses are generated via quantum effects by which new particles are not needed to be super heavy. This idea has originally been realized in the model proposed by A. Zee [5] (Zee model), where neutrino masses are generated at one-loop level 1 . In the Zee model 2 , right-handed neutrinos are not introduced, while the Higgs sector is extended by adding an additional Higgs doublet and charged singlet fields, by which the lepton number is explicitly broken via scalar interactions. The Zee model predicts a characteristic structure of the mass matrix for Majorana neutrinos, e.g., all the diagonal elements to be zero. Such a strong prediction ironically turns out to kill the model itself, because of the contradiction with the observed neutrino oscillations [8][9][10]. Therefore, modifications or extensions of the Zee model are inevitable.
The simplest modification would be just not imposing a Z 2 symmetry which is originally introduced to avoid tree level flavor changing neutral currents (FCNCs) mediated by neutral Higgs bosons. In this case, both Higgs doublet fields can couple to each type of fermions, i.e., up-type quarks, down-type quarks and charged leptons, so that we obtain sufficient sources of lepton flavor violations in order to explain current neutrino data [10][11][12][13]. How-ever, this requires unnatural fine-tunings in the quark sector to avoid various constraints from flavor changing processes such as the B-B mixing. Recently, in Ref. [14] another possibility has been proposed, where a global U (1) symmetry is introduced instead of the Z 2 symmetry 3 . Taking flavor dependent charge assignments for lepton fields, we can obtain additional sources of lepton flavor violations, and at the same time matrices for quark Yukawa interactions are diagonal.
In this paper, we clarify that the minimal Zee model with the U (1) symmetry cannot explain current neutrino oscillation data 4 . We thus add a vector-like lepton doublet, as one of the simplest extensions, to the model, in which we assume a weak mixing between the new vector-like field and the SM leptons to avoid large contributions to charged lepton flavor violating (CLFV) decays. In this extension, an anti-symmetric Yukawa matrix among the lepton doublets and the charged singlet scalar becomes a 4 × 4 form including six independent complex parameters which can be analytically solved in terms of the neutrino parameters driven by experiments. This extended model also provides new sources of CP-violation in the lepton sector, which is well motivated for the explanation of baryon asymmetry of Universe [18]. Interestingly, it is clarified that effects of the CP-violating phases on the electron electric dipole moment (EDM) are negligibly small because of the structure of the Yukawa matrices constrained by the U (1) symmetry. We then study the electroweak phase transition, and find a region of the parameter space where the strongly first order phase transition (FOPT) is realized which is needed for the successful scenario of the electroweak baryogenesis [19,20]. In addition, we discuss collider signatures of our model, particularly focusing on Higgs boson decays into a flavor violating lepton pair in the benchmark parameter points which satisfy neutrino data, CLFV data and the strongly FOPT.
This paper is organized as follows. In Sec. II, we define our model, and give the Yukawa interaction terms and the Higgs potential. Constraints from perturbative unitarity and vacuum stability are also discussed. In Sec. III, we discuss neutrino masses which are generated at one-loop level, and numerically show the necessity of the extension of the minimal Zee model in order to reproduce the current neutrino oscillation data. Sec. IV is devoted for the discussion of various constraints from flavor experiments such as the 3 The Zee model has also been extended by introducing a supersymmetry [15], an A 4 symmetry [16] and an SU (3) c × SU (3) L × U (1) X symmetry [17]. 4 We find an error in the structure of the neutrino mass matrix in Ref. [14]. electron EDM and CLFV decays. Collider phenomenologies are then discussed in Sec. V, particularly focusing on production and decay of additional neutral Higgs bosons at the LHC.
In Sec. VI, we show the electroweak phase transition as a cosmological consequences in our model. Conclusions are given in Sec. VII. In Appendix, we give the approximate formulae for each element of the anti-symmetric matrix F (Appendix A), explicit expressions for the amplitude of the CLFV decays (Appendix B) and those for the effective potential at finite temperature (Appendix C).

II. MODEL
We briefly review the Zee model with a global U (1) symmetry which has been proposed in Ref. [14]. The content of the scalar sector is the same as the original Zee model, which is composed of two isospin doublet fields Φ 1 and Φ 2 and a pair of singly-charged scalar singlets S ± . This model, however, cannot explain current neutrino oscillation data, as it will be shown in the next section. One of the simplest extensions is the introduction of a vector-like lepton doublet L 4 ≡ (ν T , T − ) T to the model. In this section, we first discuss Yukawa interactions with L 4 , and then consider the Higgs potential.

A. Yukawa Interactions
The most general form of the Lagrangian for the lepton sector is given by where L L ( R ) are the left-handed (right-handed) lepton fields. The indices A and B (= 1, . . . , 4) represent the flavor with L 1-3 L and 1-3 R to be identified with the SM lepton fields, and 4 R ≡ T R which is the charged component of L 4 R . The supscript c denotes the charge conjugation. The structure of the Yukawa matricesỸ 1 andỸ 2 is constrained by the U (1) symmetry depending on its charge assignment. Throughout the paper, we take the Class-I assignment defined in Ref. [14], where Φ 1 and the right-handed tau lepton τ R are charged with q( = 0) and −q, respectively, while all the other fields are neutral. 5 This assignment 5 The same structure of the Yukawa interaction can be realized by imposing a Z 2 symmetry instead of the U (1) symmetry, where Φ 1 and τ R are assigned to be odd. In this case, an additional term (Φ † 1 Φ 2 ) 2 appears in the Higgs potential.
provides the largest number of non-zero elements of Yukawa interaction matrices given as where × denotes a non-zero complex value. The fourth column has to be zero due to the gauge invariance. The matrixF is the anti-symmetric 4 × 4 form, so that it is described by six independent parameters. We note that Yukawa interaction terms for quarks are the same form as those in the SM, where only Φ 2 couples to quarks due to the U (1) symmetry.
Thus, FCNCs do not appear in the quark sector at tree level.
In order to separately write fermion mass terms and interaction terms, we introduce the Higgs basis [21,22] defined as where we introduced shorthand notation for the trigonometric functions as s X = sin X and c X = cos X. In addition, we defined tan β = Φ 0 2 / Φ 0 1 and In Eq. (4), G ± and G 0 are the Nambu-Goldstone (NG) bosons which are absorbed into the longitudinal components of W ± and Z bosons, respectively, while H ± , h 1,2 and A are physical charged, CP-even and CP-odd Higgs bosons, respectively. The VEV v is related to . Their mass eigenstates are defined as where the mixing angles α and χ are expressed in terms of the parameters in the Higgs potential, see Eq. (32). We identify h with the discovered Higgs boson with a mass of about 125 GeV.
In the Higgs basis, Eq. (1) is rewritten as whereỸ The mass matrix for the charged leptons is expressed bỹ This matrix can be diagonalized by the unitary rotations The interaction terms are then extracted in the mass eigenstates of the Higgs bosons and the charged fermions as where P L (P R ) is the projection operator for the left (right) handed fermions, and These Yukawa matrices can be rewritten by using Eq. (9) as, where P 123 ≡ diag(1, 1, 1, 0) and P β ≡ diag(cot β, cot β, − tan β, 0). We note that V L does not appear in the above expression due to its unitarity.
In general, the matrices Y Φ and Y Φ contain non-zero off-diagonal elements which can introduce large effects on CLFV decays mediated by scalar bosons. In order to avoid such CLFV decays, we assume that the fourth lepton doublet L 4 is weakly mixed with three generations of the SM leptons. This can naturally be realized by taking the mass parameter M to be much larger than charged lepton masses, e.g., the tau lepton mass.
Let us clarify how the interaction matrices defined above can be expressed in the weak mixing scenario. First, the unitary matrices V R can be expressed as (13) where V ij are the 4 × 4 unitary matrices with i-i and j-j elements to be cos α ij and i-j , while all the other diagonal (off-diagonal) elements to be unity (zero). Second, the weak mixing scenario indicates that the mixing angles α a4 (a = 1, 2, 3) are small, so that we reparametrize them as α a4 → α a4 with 1. In this case, the unitary matrices can be expressed as, where V R,0 is given by with the upper-left (lower-right) block being the 3 × 3 (1 × 1) form. Third, using this expansion, the Yukawa matrices Y Φ and Y Φ are expressed as where with X 3×3 (X 3×1 ) being the 3 × 3 (3 × 1) part of a 4 × 4 matrix X. The expressions for Y Φ,0 and δY Φ are given by replacing P 3×3 Interestingly, the mixing effect on masses of active neutrinos is not suppressed by 2 , but as it will be explained in the next section. It is also important to mention here that F , δY Φ and δY Φ can be complex, so that they can provide new sources of CP-violation, and their effects on EDMs will be discussed in Sec. IV.
It is clear that in the → 0 limit, Y Φ becomes diagonal, while Y Φ takes the block diagonal form as shown in Eq. (17). Furthermore, in the scenario with a softly-broken Z 2 symmetry, the matrix Y Φ also becomes proportional to the diagonalized mass matrix as in the two Higgs doublet models (THDMs): where ζ = cot β (− tan β) for the Type-I and Type-Y (Type-II and Type-X) THDM [23].

B. Higgs Potential
The most general Higgs potential is given by In the above expression, m 2 3 and µ terms softly break the U (1) symmetry, and their complex phases are removed by using the phase redefinition of the scalar fields without the loss of generality. Thus, there is no CP-violating phase in the Higgs potential. We note that the (Φ † 1 Φ 2 ) 2 term is forbidden because of the U (1) symmetry. Thus, a pseudo-NG boson, corresponding to A, appears due to the spontaneous breaking of the U (1) symmetry. The µ term plays a crucial role for the neutrino mass generation, as this term breaks two units of the lepton number when we assign the lepton number of (Φ 1 , Φ 2 , S ± ) to be (0, 0, ±2) 6 .
After imposing the tadpole conditions for two CP-even Higgs bosons, all the masses of 6 By this assignment, the lepton number is conserved in the Yukawa interaction terms. the physical Higgs bosons are expressed as where (M 2 even ) ij and (M 2 ± ) ij (i, j = 1, 2) are the elements of the squared mass matrices for the CP-even and singly-charged Higgs bosons in the basis of (h 1 , h 2 ) and (H ± , S ± ), respectively.
Each element is given as The mixing angles are expressed in terms of these matrix elements: From the above discussion, we can choose the following twelve parameters as inputs: where tan σ ≡ σ 2 /σ 1 . Among the above parameters, v and m h are fixed to be about 246 GeV and 125 GeV by experiments, respectively, and σ 3 is not relevant to the following discussion.
The parameters of the Higgs potential are constrained by considering bounds from perturbative unitarity and vacuum stability. In Refs. [24,25], all the independent eigenvalues of the s-wave amplitude matrix for 2 body to 2 body elastic scatterings (a i ) have been given in the high energy limit. Requiring |a i | ≤ 1/2, the unitarity bound is expressed in our notation as 1 2 where x 1,2,3 are the eigenvalues for the following 3 × 3 matrix The vacuum stability of the potential has also been discussed in Refs. [24,26]. Requiring the potential bounded from below in any direction with large scalar field values, the parameters in the Higgs potential are constrained to be within the following domain [24] where Ω 1 = λ 1 , λ 2 , σ 3 > 0; λ 1 σ 3 + σ 1 > 0; λ 2 σ 3 + σ 2 > 0; with D = min(0, λ 4 ).

III. NEUTRINO MASSES
Majorana masses for the active neutrinos are generated at one-loop level. In the ν L (= V L ν L ) basis, we obtain the mass matrix as where C diag is the diagonal matrix given by with m A being the charged lepton mass. For m A We note that the mass of the sterile neutrino ν T is simply given by the vector-like mass M . 7 In the weak mixing scenario, the neutrino mass matrix can be expressed as where We note that the contribution from the order 2 term m 2 ν is much smaller than that from m 1 ν due to not only the suppression of the factor but also no enhancement by the large lepton mass m T . Therefore, the neutrino masses can be well approximated by considering the term up to m 1 ν . The neutrino mass matrix becomes the one given in the original Zee model by taking → 0 and P 3×3 β → cot βI 3×3 in Y Φ ,0 , which is consistent with the expression given in Ref. [10]; As aforementioned, this structure cannot accommodate current neutrino oscillation data.
The neutrino mass matrix can be diagonalized by introducing the PMNS matrix U PMNS as follows: where m i (i = 1, 2, 3) are the mass eigenvalues. For the normal ordering (NO) and the inverted ordering (IO) cases, these are m 1 m 2 m 3 and m 3 m 1 m 2 , respectively.
The matrix U PMNS can be parameterized as, where δ CPV is the Dirac CP-phase, and ϕ 1,2 are the Majorana phases. From Eqs. (42) and (47), we obtain By solving this equation, the elements of F are expressed in terms of the following model parameters: with {α AB } = {α 12 , α 13 , α 23 , α 14 , α 24 , α 34 } (similarly to ϕ AB ), and the neutrino parameters appearing in the right-hand side of Eq. (49): with m 0 being the smallest eigenvalue of the neutrino masses. We define the two squared mass differences as ∆m 2 sol ≡ m 2 2 − m 2 1 and ∆m 2 atm ≡ |m 2 3 − m 2 1 | (|m 2 3 − m 2 2 |) for the NO (IO) case. Now, we can reproduce the neutrino parameters shown in Eq. (51) for each fixed value of the model parameters given in Eq. (50). Typical order of the elements of F can be estimated for tan β 1 as follows: where m ν is the typical value of the elements of the neutrino mass matrix, i.e., m ν = O(10 −3 -10 −2 ) eV. In Appendix A, we present better approximated formulae for these elements.

A. Electric Dipole Moments
As we have discussed in the previous section, non-zero CP-violating phases appear in the lepton sector. In this subsection, we study effects of these phases on EDMs, particularly the electron EDM 9 . At one-loop level, contributions from the neutral Higgs bosons (h, H and A) are expressed as where and the loop function A is given by From Eq. (12), we can easily show that the product (Y Φ Y Φ ) 11 becomes a real value: 8 As mentioned in Introduction, if we do not impose any new symmetries such as Z 2 and U (1), we can find a solution to satisfy the current neutrino data. This, however, requires a fine-tuning in the quark Yukawa sector in order to suppress FCNCs mediated by scalar bosons. 9 The CP-violating source in the quark Yukawa sector is the same as in the SM, i.e, only arising from the Kobayashi-Maskawa phase. Thus, its effect on EDMs is negligibly small [29].
A similar argument holds for (Y Φ Y Φ ) 11 just by replacing P 123 → P β in the above expression.
On the other hand, a non-zero imaginary part appears from the product (Y Φ Y Φ ) 11 as follows This combination appears from the h or H loop contribution, which is proportional to Thus, the magnitude of the electron EDM can be estimated for m T m φ as Therefore, the typical value of the electron EDM is much smaller than the current upper limit given from the ACME collaboration, i.e., |d e /e| < 1.1×10 −29 cm at the 90% confidence level [30]. We note that diagrams with charged Higgs boson loops are negligibly small, because these contributions are proportional to the neutrino mass 10 . Thus, the one-loop contributions to the electron EDM can safely be ignored.
One might think that two-loop Barr-Zee (BZ) type diagrams [31] give an important contribution to the electron EDM, rather than the one-loop contribution as in CP-violating THDMs [32][33][34]. However, their contributions also do not give a significant contribution as explained below. For BZ diagrams with a neutral gauge boson (γ or Z) exchange, they are proportional to Im(Y Φ ) 11 or Im(Y Φ ) 11 which is zero in our model. For those with a W boson exchange, its contribution can further be separated into that with H ± and S ± exchanges, where the former is again proportional to Im(Y Φ ) 11 and the latter is negligibly small as it picks us a tiny neutrino mass from an internal neutrino line. 11 Therefore, our model is safe from the constraints of the EDMs even if we consider O(1) CP-phases in the Yukawa couplings.

B. Constraints from Charged Lepton Flavor Violation
In this subsection, we discuss CLFV processes. The Yukawa interactions in Eq. (10) Branching ratios 10 Contributions from ν T and the charged Higgs boson loop vanish, because these are proportional to 11 . 11 This argument is not qualitatively changed in the discussion with the mass eigenstates of the charged Higgs bosons H ± 1 and H ± 2 .
(BRs) of these processes are given by where α em is the fine structure constant, and C ij are numerical constants associated with the BR of the charged lepton, i.e., C 21 = BR(µ → eν µνe ) = 1, C 31 = BR(τ → eν τνe ) = 0.1784 and C 32 = BR(τ → µν µ ν e ) = 0.1736. In Eq. (59), a φ L,R denote an amplitude which is estimated from a one-loop diagram with a scalar boson φ = {h, H, A, H ± 1 , H ± 2 } running in the loop; explicit forms of these amplitudes are summarized in Appendix B.
In addition we consider spin-independent µ → e conversion via H exchange; CP-odd scalar exchange induces spin-dependent µ → e conversion process which is less constrained.
We can write the BR for the process such that [35][36][37] where Γ cap is the rate for the muon to transform to a neutrino by capture on the nucleus, S (p,n) is the integral over the nucleus for lepton wavefunctions with corresponding nucleon density, and f N ∼ 0.3 is effective coupling between Higgs and nucleon N defined  [36,39]. The current upper limits on the above BRs with 90% confidence level are found in Refs. [40][41][42][43][44][45] as We note that the Yukawa interactions also induce three body decays such as µ → eeē via the diagram with off-shell H and A at tree level. BRs of such decay modes are found to be much smaller than the current upper bound, e.g. BR(µ → eeē) < 1.0 × 10 −12 [44], because BR(µ → eeē) is proportional to the small lepton masses, i.e., (m e m µ /v 2 ) 2 = O(10 −17 ).
We also note that the Z boson couplings with right-handed charged leptons can contain flavor violation, because 1-3 R and 4 R (= T 4 R ) belong to the different representation under the SU (2) L × U (1) Y group. In the mass eigenbasis, the coefficients of the Z µ¯ a R γ µ b R vertices (a, b = 1, 2, 3 and a = b) are proportional to 2 , so that the off-shell Z boson contribution to the decay rate of µ → eeē is suppressed by 4 . Thus, we can safely avoid the bound by taking < O(10 −3 ). 12 We thus concentrate on the constraints from the i → j γ processes in the following discussion.
In the following, we perform numerical evaluations of the neutrino mass matrix and the BRs of the CLFV decays. First, we scan the input parameters as We fix the mass parameters to be (m T , m H , m H ± 1 , m H ± 2 , m A ) = (1000, 350, 350, 300, 70) GeV, where the scalar boson masses are chosen such that the FOPT can be realized as we see in the next section. In addition, the other parameters are fixed to be s χ = 0.01, = 10 −4 and ϕ 1,2 = ϕ AB = 0. We then solve Eq. (49) to obtain F AB for each set of the input parameters, where we use the best fit values of the neutrino parameters θ ij , |∆m 2 atom | and ∆m 2 sol from the global fit results [27,28] for both the NO and IO cases. Finally, we apply the alignment limit s β−α = 1, which is supported by the current LHC data [49,50], and is favored by the constraint from the electroweak ρ parameter in the case of m H = m H ± 1 13 .
In Fig. 1, we show the correlations between the BRs of i → j γ processes, where the left and right plots correspond to the NO and IO cases. The upper bounds of the BRs are indicated by the dashed lines. We find a strong prediction of BR(τ → µγ) BR(τ → eγ) BR(µ → eγ) in both the NO and IO cases. We see that larger BRs tend to be predicted in the IO case as compared with the NO case, and only little parameter sets are allowed in the IO case. We note that if we adopt smaller s χ such as 10 −3 , most of the parameter sets are excluded since the magnitude of F AB increases by one order resulting in larger BRs.
In Fig. 2, we also show correlation between BR(µ → eγ) and BR(µ → e) Al where the latter is µ-e conversion rate for 27 13 Al. We find that BR(µ → e) Al > BR(µ → eγ) in some 12 Milder bounds can be obtained from the decay of the Z boson into a pair of different lepton flavor. The current upper limits on BR(Z → ) are given to be of order 10 −7 , 10 −6 and 10 −5 for the eµ, eτ and µτ final state [46], respectively. We also check that bounds from the µ-e conversion provide a slightly stronger limit on as compared with those from µ → 3e [44,47,48]. 13 In our numerical calculation, we take the small mixing angle for the charged Higgs bosons, i.e., s χ = 10 −2 , so that the H ± 1 states almost coincide with the charged Higgs bosons from the doublet H ± . Thus, the one-loop corrections to the ρ parameter from the scalar boson loops almost completely vanish due to the approximate custodial symmetry [51,52]. region. These region correspond to the case where (Y Φ ) 11 (Y Φ ) 12 (21) so that µ → eγ is more suppressed than µ-e conversion.
Before closing this section, let us briefly comment on the other observables in the lepton sector under the constraints from the CLFV data. In particular, we consider the effective Majorana mass m ee for the neutrinoless double beta decay 0νββ, which can be non-zero as follows: m ee = |m 1 cos 2 θ 12 cos 2 θ 13 + m 2 sin 2 θ 12 cos 2 θ 13 e iϕ 1 + m 3 sin 2 θ 13 e i(ϕ 2 −2δ CPV ) |.
In our model, however, there is no particular prediction of the neutrino mixing angles, the squared mass differences and δ CPV , so that m ee can be obtained just by inputting the the value of m ν larger than cosmological limit by PLANCK m i < 0.12 eV [54]. In future searches for the neutrinoless double beta decay, m ee would be probed below about 0.05 eV [53], and our model can be indirectly tested.

BR(H→XX)
where θ W is the Weinberg angle and λ HAA = (m 2 H − m 2 A ) cot 2β/v. Differently from the leptonic decays of the Higgs bosons, decay rates into a quark pair are the same as those given in the Type-I THDM. Thus, among the hadronic decays, the decay rates of H/A → bb We here comment on the decay of the SM-like Higgs boson h. In the THDMs, the couplings of h coincide with those of the SM Higgs boson in the alignment limit at tree level. In our model, however, the Yukawa matrix Y Φ is not exactly the same as the SM one even in the alignment limit, i.e., non-zero off-diagonal elements appear in the fermion mass eigenbasis. Such off-diagonal component is highly suppressed by the factor of 2 , so that the size of BRs of LFV decays of h can be estimated by O( 4 ). In fact, we numerically checked  [55]. The shaded region is excluded by the constraint from B d → µµ [56,57]. Let us move on to the discussion of production processes of the additional Higgs bosons and constraints from current LHC data. At the LHC, the additional Higgs bosons can mainly be produced via the gluon fusion process gg → H/A. 14 The pair productions pp → HA/H ± H/H ± A can also be used for smaller mass cases. So far, no discovery of significant signatures has been reported, and it has taken lower limits on their masses and upper limits on relevant coupling constants. See, e.g., Ref. [58] for the recent analysis of the constraints from direct searches at the LHC in the THDM. We see in [58] that in the alignment limit s β−α = 1 with degenerate masses of the additional Higgs bosons, most of the parameter region has not been excluded yet, e.g., tan β 2 is allowed in the Type-I THDM. It is important to mention here that these searches typically have sensitivities to relatively larger mass regions, e.g., above 100 GeV, so that they cannot be simply applied to cases for lighter additional Higgs bosons and/or the case with mass differences.
From the above analysis, we see that H can mainly decay into ZA and AA, and A further decays into a lepton pair. Thus, clear signals of H can be obtained by the decay chain of H → ZA/AA → 4 , and its cross section can be estimated by where σ(gg → H) is the gluon fusion cross section for H, which can be estimated to be 4 pb × cot 2 β for m H = 350 GeV [59]. A similar process with the same final states has been searched at the LHC, i.e., pp → Z → φφ(φ → 2 ) process ( = e, µ) with Z and φ being an extra neutral gauge boson and a neutral scalar boson, respectively. For m Z = 350 GeV which can be replaced by m H in our case, the upper limit has been set to be σ 4 ∼ 2 fb [55]. In the left panel of Fig. 5, we show the sum of the BRs of A → ( = e or µ) and that including at least one τ . We see that the former quantity rapidly increases as 10 −3 , 10 −2 and 10 −1 at tan β 2, 4 and 6, respectively, so that at high tan β some parameter points can be excluded by the current LHC data. In fact, as seen in the right panel of where the value of cross section times BR has been constrained to be smaller than about 0.1 pb with the 95% confidence level. In our model, the gluon fusion production of A with the mass of 70 GeV is given to be about 40 pb× cot 2 β [59], while the BR of A → γγ is given of order 10 −4 . Therefore, the cross section times BR is much below the current upper limit. At future collider experiments such as the High-Luminosity LHC [62,63] and lepton colliders 15 , e.g., the International Linear Collider (ILC) [64][65][66], the Circular Electron-Positron Collider (CEPC) [67] and the Future Circular Collider (FCC-ee) [68], our scenario could be tested via the LFV decays of the Higgs bosons with the characteristic decay pattern. It goes without saying that dedicated studies are required to clarify the feasibility of such signatures, and such analyses are beyond the scope of the present paper.

VI. ELECTROWEAK PHASE TRANSITION
In this section, we consider the cosmological consequences of our scenario, particularly focusing on the electroweak phase transition. It has been known that in the electroweak baryogenesis scenario [19,20], the strongly FOPT is required to realize sufficient departure from thermal equilibrium in order to maintain non-zero baryon asymmetry of the Universe.
The criteria for the strongly FOPT can be expressed by ϕ c /T c 1 or ϕ n /T n 1, where T c (T n ) is the critical temperature providing degenerate vacua (nucleation temperature of electroweak bubbles from the bounce solution) and ϕ c (ϕ n ) is the order parameter at T c (T n ).
Although we should use ϕ n /T n rather than ϕ c /T c for the estimation of the strength of the FOPT, we show both of them for comparison. Typically, these two valuables are almost the same with each other unless considering a mechanism of the supercooling scenario in which the thermal barrier does not disappear until much lower temperature T n T c [69,70].
It has been known that additional bosons can enhance the FOPT [71], because their loop effects provide a positive contribution to a cubic field term of the effective potential at finite temperature, which makes a potential barrier higher at around the critical temperature.
Indeed, it has been clarified that additional Higgs bosons in the THDM can make the FOPT stronger as a consequence of the non-decoupling effect if their masses mainly come from the VEV [72]. This can also be described in such a way that a new dimensionful parameter 15 At the LEP experiments, the light A can be produced associated with a fermion pair, i.e., e + e − → ff A.
We have checked that its cross section is of order 10 (1)  M which is irrelevant to the VEV and appears in the mass of Higgs bosons is taken to be smaller than the physical mass parameter of the additional Higgs bosons. In our model, such M parameter directly corresponds to the mass of A, see Eq. (21), so that m A has to be smaller than the other masses of the additional Higgs bosons in order to enhance the FOPT. In Appendix C, we present detailed analytic expressions for the effective potential at finite temperature. We note that the loop effect of the charged singlet fields S ± can not be dominant in our model, because of the smaller degrees of freedom (N = 2) as compared with those of the doublets (N = 8) 16 . We numerically find that the FOPT can be maximal when the mass of singlet-like Higgs boson m H ± 2 is taken to be around 300 GeV in our scenario. In the following, we numerically evaluate the strength of the FOPT, i.e., ϕ c /T c or ϕ n /T n by using CosmoTransitions [76]. We have checked that the behavior of the FOPT in the case of the THDM is consistent with the previous work [77]. Our model predictions can then be simply obtained by taking M = m A in the THDM and modifying the charged Higgs sector by that with H ± 1 and H ± 2 . As mentioned above, we take a smaller value of m A but larger than m h /2, e.g., 70 GeV, to avoid the constraint from h → AA at the LHC, see previous section. In addition, we take into account the constraints from perturbative unitarity and vacuum stability discussed in Sec. II and flavor experiments [56,57] such as In Fig. 6 (left), we show ϕ c /T c and ϕ n /T n as a function of tan β for the case with m A = 70 GeV, m H = m H ± 1 = 350 GeV and m H ± 2 = 300 GeV. We here neglect the effect of the small mixing angle between two charged Higgs bosons. For concreteness, we take s χ = 0. In this plot, we do not impose the constraints from the flavor experiments and the unitarity bound.
We find that in some parameter points indicated by the triangle a two-step phase transition happens, where the empty (filled) triangles represent the value of ϕ c /T c or ϕ n /T n at the first (second) transition, while the circle points represent the case with the one-step phase transition. We see that both ϕ c /T c or ϕ n /T n become maximal at around tan β = 1.2, and ϕ c /T c (ϕ n /T n ) is greater than unity at 0.7 tan β 1.9 (0.7 tan β 2.0). In addition, the value of ϕ n /T n is slightly larger than ϕ c /T c , which can be more clearly seen from the right panel of Fig. 6, in which we show the correlation between ϕ c /T c and ϕ n /T n . We see that ϕ n /T n is typically about 10% larger than ϕ c /T c . This can mainly be explained by the 16 In Refs. [73][74][75], it has been shown that a larger number of singlet scalar fields significantly enhance their loop effects on the effective potential, and makes the FOPT stronger. difference between ϕ c and ϕ n (ϕ n > ϕ c ) whose values are quite sensitive to the temperature at around the phase transition. In fact, around the 10% difference between ϕ c and ϕ n arises from the slight difference between the temperatures, i.e., (T c − T n )/T c = O(10 −3 ).
In Fig. 7, we show the contour plots for ϕ/T at T = T c (left) and T = T n (right) in the m H (= m H ± 1 )-tan β plane using the same parameter set taken in the right panel of Fig. 6. Under the constraints from the flavor experiments and the unitarity bound, we find that there is a region of the parameter space which realizes the strongly FOPT (ϕ n /T n 1) shown as the white area in the right panel. It is clear that such region requires the quite limited parameter choice, i.e., 1.9 tan β 2.3 and 250 m H 400 GeV, which provides a good benchmark scenario for the collider phenomenology, as we discussed it in the previous 1. section.
It is worth discussing detectability of the gravitational waves (GWs) originated from the first-order phase transition since the typical frequency of GWs from the first-order phase transition at the electroweak scale is expected to explore by future planned space-based interferometers such as LISA [78], DECIGO [79] and BBO [80] which survey GWs in the millihertz to decihertz range. The GW spectrum can be parameterized by two dimensionless parameters, related to the released energy density and the change rate of the three dimensional bounce action S 3 , defined as where ρ rad is the radiation energy density. We employ the analytic formula provided in Ref. [81] to estimate the spectrum of the GWs, Ω GW h 2 (f ). Parameters have one-to-one correspondence (tan β, m H ) ↔ (α, β) ↔ (f, Ω GW h 2 ) and correlations with the phase transition strength ϕ n /T n ∝ Ω GW h 2 ∝ α ∝ β −1 .
In Fig. 8, the predicted GW signals are plotted taking parameters on Fig. 7. Red colors on each panel are allowed by the constraint from B d → µµ. We find that the GW signals do not reach the future planned sensitivities of observations. Because the exact non-decoupling limit does not allow due to the global U(1) and the constraint from B d → µµ excluds most of the parameter space, the GW signals correlated with the strength of the first-order phase transition cannot be enough strong.
Finally, we give a comment on the possibility for the generation of baryon asymmetry of the Universe. As it is known that the strongly FOPT is one of the necessary conditions to obtain the non-zero baryon number in the electroweak phase transition. In order to estimate the generated baryon number, we also need to calculate the reflection rate of chiral fermions against the bubble, in which non-zero charges, e.g., the hypercharge and isospin, can be accumulated at the vicinity of the bubble if CP-violation happens in the interaction among the bubble and fermions. Such non-zero charges can be converted into the baryon number from the sphaleron process, which is frozen in the broken phase because of the decoupling of the sphaleron process. See e.g., [83] for more detailed discussions. In Ref. [84], the electroweak baryogenesis has been discussed in a two Higgs doublet model by introducing CP-violating phases in lepton Yukawa couplings, particularly τ -τ and τ -µ couplings to additional Higgs bosons. It has been shown that in order to explain the observed baryon asymmetry the magnitude of these Yukawa couplings has to be of order 0.1-1. We thus plot the values of corresponding Yukawa couplings (Y Φ ) 33 and (Y Φ ) 32 in our model in Fig. 9. It is seen that in the region where the FOPT is realized, i.e., tan β 2, the value of (Y Φ ) 33 and (Y Φ ) 32 can be maximally of order 0.01. Therefore, qualitatively our model can generate baryon asymmetry which is one or two orders of magnitude smaller than the required value, so that additional sources of baryon asymmetry would be needed.

VII. CONCLUSIONS
We have investigated the Zee  Finally, as a cosmological consequence, we have studied the electroweak phase transition by using the effective potential at finite temperature. The content of the Higgs sector in our model is the same as that of the original Zee model, i.e., two Higgs doublets and a pair of charged singlet scalars. Due to the U (1) symmetry, the coefficient of the (Φ † 1 Φ 2 ) 2 , denoted as λ 5 , is forbidden, and thus the CP-odd Higgs boson A becomes a pseudo-NG boson whose mass arises from the (Φ † 1 Φ 2 ) term. Because of this property, we need a light A in order to make non-decoupling effects of the additional Higgs bosons stronger, which plays an important role to realize the strongly FOPT. We have found that the strongly FOPT can be realized in the case where 1.9 tan β 2.3 and 250 m H 400 GeV for m A = 70 GeV, which is compatible with the constraints from the flavor experiments and LHC data. Therefore, we have shown that our model contains two important ingredients to realize the successful electroweak baryogenesis scenario, i.e., the additional CP-violating phases and the strongly FOPT. According to the qualitative discussion, our model can generate baryon asymmetry which is one or two orders of magnitude smaller than its observed amount, so that additional sources of baryon asymmetry would be required. In Sec. III, we show that the matrix elements F can be expressed in terms of the neutrino mass matrix by solving Eq. (49). Neglecting the electron mass except for terms proportional to 1/m e , we find the following approximate formulae for the elements of the F matrix: where y ij ≡ (Y 0 Φ ) ij , δ ij ≡ (δY Φ ) ij and r ≡ m µ /m τ . The matrix elements m ij dat denote those of the right-hand side of Eq. (47). These expressions agree with the exact ones typically within a 10% error.

Appendix B: Decay Amplitudes for i → j γ Processes
Here, we summarize the analytic formulae for the amplitudes of the i → j γ processes denoted by (a L,R ) ij in Eq. (59). The contributions from diagrams with the charged Higgs bosons inside loop are obtained as The loop functions are defined by In order to discuss the one-loop effective potential at finite temperature, we introduce the classical constant field configurations ϕ = (ϕ 1 , ϕ 2 ) where (ϕ 1 , ϕ 2 )| T =0 ≡ (v 1 , v 2 ) at the VEVs of our present vacuum. We here do not consider the possibility of CP-breaking and/or charge-breaking vacua for simplicity.
The effective potential receives additional contributions from thermal loop diagrams, and is modified to [71] V eff ( ϕ, T ) = V tree ( ϕ) with V tree ( ϕ) ≡ m 2 where n i and M i ( ϕ) denote the degrees of the freedom and the field-dependent masses for particles i = h, G 0 , H, A, G ± , H ± 1 , H ± 2 , γ L,T , Z L,T , W ± L,T , t and b, respectively. Namely, n h,G 0 ,H,A = 1, n G ± ,H ± 1 ,H ± 2 = 2, n γ L ,Z L = 1, n W ± L = 2, n γ T ,Z T = 2, n W ± T = 4, n t,b = −12. (C5) The renormalization scale Q is set at v in our analysis. We take the MS scheme, where the numerical constants c i are determined to be 3/2 (5/6) for scalars and fermions (gauge bosons). The thermal correction is given by I B,F (a 2 i ) = ∞ 0 dx x 2 ln 1 ∓ exp(− x 2 + a 2 i ) for bosons (−) and fermions (+), respectively.
In order to take ring-diagram contributions into account, we have introduced the fielddependent masses depending on the temperature in the effective potential by [85] where Π i (T ) denote the finite temperature contributions to the self energies of the fields i.