Neutrino masses and magnetic moments of electron and muon in the Zee Model

We explore parameter space in the Zee Model to resolve the long-standing tension of the electron and muon anomalous magnetic moment (AMM). The model comprises a second Higgs doublet and a charged singlet at electroweak scale and generates Majorana neutrino masses at one-loop level; the neutral partner of the $SU(2)_L$ doublet contributes to the AMM of electron and muon via one loop and two-loop corrections. We propose two minimal flavor structures that can explain these anomalies while fitting the neutrino oscillation data. We find that the neutral Higgs resides in the mass range of roughly 10-300 GeV or 1-30 GeV, depending on the flavor structures. The model is consistent with constraints from colliders, electroweak precision data, and lepton flavor violation. To be comprehensive, we examine the constraints from the electric dipole moment (EDM) and find a region of parameter space that gives a sizable contribution to muon EDM while simultaneously giving corrections to muon AMM. In addition to the light scalar, the two charged scalars with masses as low as 100 GeV can induce nonstandard interactions $\varepsilon_{ee}$ as large as $8\%$, potentially hinting at new physics. We also investigate the projected capability of future lepton colliders to probe the currently allowed parameter space consistent with both electron and muon AMMs via direct searches in the $\ell^{+}\ell^{-}\to \ell^{+}\ell^{-}(H \to \ell^{+}\ell^{-})$ channel.

The flavor-changing nature of the second Higgs doublet that gives rise to large contributions to a also plays a crucial role in neutrino masses and oscillations, which reveals insights into the flavor structure within the model. The opposite signs of the anomalies can be explained by either choosing Yukawa couplings with opposite signs or adjusting the mass splitting between the scalars so that ∆a µ is positive from the dominant one-loop diagram. In contrast, the two-loop Barr-Zee diagram provides negative correction to a e . We propose two minimal Yukawa textures that achieve these while also providing excellent fits to the neutrino oscillation parameters. We find that the simultaneous explanation of the observed disparities in the lepton AMMs requires the scalar mass to be in the range of 10 − 300 GeV for TX-I and 1 − 30 GeV for TX-II while satisfying various LFV measurements such as 1 → 2 γ, including other relevant constraints from low energy physics, collider searches, and fit to the neutrino oscillation data. It is also worth noting that the Yukawa couplings, if complex, can not only generate a but could also have a sizable effect on the electric dipole moment (EDM) d . We find that an extremely small complex phase of O(10 −5 ) is required to satisfy the current limits of electron EDM [143] while simultaneously satisfying a e . However, this is not the case for muon EDM [144], as there exist regions of parameter space that can give a sizable contribution to muon EDM, which can potentially be measured in future experiments [145][146][147]. The charged scalars in the model can induce nonstandard neutrino interactions (NSI) ε ij [148][149][150] (for a recent review on NSI in the context of the Zee Model, see Ref. [17]), which, if observed, could be direct indicators for new physics. We find that diagonal NSI can be as large as 8% for ε ee .
We have also evaluated constraints on the Yukawa couplings for both the textures from direct searches in the e + e − → + − (H → + − ) ( = e, µ) channel at LEP. Moreover, we also explore the future sensitivity of the channel, as mentioned earlier, to the Yukawa couplings at the ILC operating at √ s = 1 TeV with an integrated luminosity L = 500 fb −1 [151][152][153][154][155]. We also study the projected sensitivity of the µ + µ − → µ + µ − (H 2 → µ + µ − ) channel at a muon collider (MuC) with configuration: √ s = 3 TeV, L = 1 ab −1 [156,157]. These searches could preclude simultaneous explanation of both AMMs in the mass range of roughly 10 − 300 GeV depending on the Yukawa coupings for a given texture within the model.
The rest of the paper is organized as follows. In Sec. 2 we present the basic description of the Zee Model, including the Yukawa lagrangian and radiative neutrino mass generation. In Sec. 3 we discuss the anomalous magnetic moments and the possible ways of resolving the two anomalies simultaneously in the model. This section also points out the model predictions for muon EDM (cf. Sec. 3.1). Sec. 4 briefly summarizes the NSI from charged scalars in the model, followed by a detailed study of the various textures of Yukawa coupling matrices which could incorporate both anomalies in Sec. 6. The low energy constraints such as i → j γ, trilepton decays, T-parameter constraints, muonium oscillations, and direct experimental constraints are discussed in Sec. 5. In Sec. 7 we perform a detailed cut-based collider analysis at the detector level to estimate the projected capability of future lepton colliders to probe the Yukawa structure of the additional scalar Higgs boson considered in this work. The results of our analysis on the anomalous magnetic moment of lepton and neutrino oscillation fit are given in Sec. 8, followed by the conclusion in Sec. 9.

Model Description
The Zee model [13,158], built on the Two Higgs Doublet Model [159,160], is perhaps the simplest extension of the SM that can generate non-zero radiative neutrino mass at a one-loop level. The model is based on SM gauge symmetry SU (3) C × SU (2) L × U (1) Y and consists of a SU (2) L doublet scalar H 2 in addition to an SM Higgs doublet H 1 and a charged singlet scalar η ± . In Higgs basis [161], where only the neutral component of H 1 takes a vacuum expectation value (VEV), H 0 1 = v 246 GeV, the doublets can be represented as where (G + , G 0 ) are the Goldstone modes, (H 0 1 , H 0 2 ) and A are the neutral CP-even and CP-odd scalars, and H + 2 is a charged scalar field. In the CP conserving limit, where the quartic couplings are real, the field A decouples from {H 0 1 , H 0 2 }. Then one can rotate the CP-even states into a physical basis {h, H}, where the mixing angle parametrized as is given by Here λ 6 is the quartic coupling of the term We base our analysis under the alignment/decoupling limit [162][163][164][165], when α → 0, agreeing with the LHC Higgs data [166,167]  with the mixing angle ϕ defined as where, µ is the coefficient of the cubic coupling H α 1 H β 2 αβ η − in the scalar potential, with αβ being the SU (2) L antisymmetric tensor. This cubic term, along with Eq. (2.6), would break the lepton number by two units.
The leptonic Yukawa interaction in the Higgs basis can be expressed as where, {i, j} are flavor indices, and H a ≡ iτ 2 H a , τ 2 being the second Pauli matrix. l c and L represent the left-handed antileptons and lepton doublets. Note, f ij = −f ji is an antisymmetric matrix and can be made real by a phase redefinitionP fP , whereP is a diagonal phase matrix, whereas {Y, Y } are general complex asymmetric matrices. We assume that H 2 is leptophilic to avoid dangerous flavor violating processes in the quark sector, such as π + → e + ν which would otherwise occur at unacceptably large decay rates for Y ie ∼ O(1). Here, after electroweak symmetry breaking, M = Y v √ 2 is the charged lepton mass matrix, and chosen to be diagonal without loss of generality.
Neutrino masses in this model are zero at the tree level. However, due to explicit lepton number violation, a non-zero Majorana neutrino mass M ν is induced as quantum corrections at the one-loop level, as shown in Fig. 1. The dot in the internal fermion line represents the mass insertion due to the SM Higgs VEV. The Yukawa couplings of Eq. (2.6), together with the H 1 H 2 η trilinear term in the scalar potential, guarantees lepton number violation. These interactions result in an effective (∆L [17,[168][169][170][171]. Note that a companion diagram can be obtained by reversing the arrows of the internal particles. Thus, the neutrino mass matrix is given by where κ is the one-loop factor κ = 1 16π 2 sin 2ϕ log With the other possibility, namely, Y 1, the stringent charge LFV (cLFV) constraints on f Yukawa coupling restrict the maximum NSI to ≤ 10 −8 [172], well below any experimental sensitivity in the foreseeable future. Furthermore, the presence of a singly charged singlet leads to lepton universality violation which, for instance, would alter the decay rate of the muon. The Fermi constant extracted from the modified muon decay would be different from the SM. This new Fermi constant has constraints from the CKM unitarity measurements giving strong limits on the Yukawa couplings, f [16]. The charged current interactions leading to leptonic decays will also be modified so that such interactions are no longer universal, leading to further constraints [173,174]. Moreover, choosing Y 1, one cannot get the required correction to g − 2 of electron and muon, and η + always contributes the wrong sign to the g − 2 of muon (cf. Sec. 3 for details). The Yukawa couplings f and Y on the mass basis of charged leptons can be written as where Y is multiplied by (ν e ,ν µ ,ν τ ) (or (e L , µ L , τ L )) from the left and (e R , µ R , τ R ) T from the right in the charged scalar, H + 2 (or neutral scalar, 1 √ 2 (H 0 2 + iA)) interaction. It is worth mentioning that the Yukawa coupling Y cannot be taken diagonal; otherwise, all the diagonal entries of the neutrino mass matrix would vanish, yielding neutrino mixing angles that are not compatible with the neutrino oscillation data [158,[175][176][177].

Anomalous Magnetic Moments and Related Process
Virtual corrections due to the new scalar fields present in the model can modify the electromagnetic interactions of the charged leptons. The contribution from Y of Eq. (2.6) is part of the SM contribution, a SM , in the decoupling limitα = 0. The contribution from f Yukawa coupling to AMMs is negligible due to the strong limit from cLFV and constrains from tiny neutrino mass, as aforementioned. Thus, the Yukawa Y interactions of leptons with the physical scalars in the alignment limit is given by where, φ = (H + iA)/ √ 2. The charged scalar h + contribution is ignored by choosing a negligible mixing between charged scalars ϕ ∼ 0. Neutral scalar contributions to anomalous magnetic moments at one-loop [178] as shown in Fig. 2 (a) is where, In the above expression, + and − correspond to H and A, respectively, with the second part representing the chiral enhancement. The contribution from charged the Higgs H + from Fig. 2  There are also two-loop Barr-Zee [179,180] diagrams arising from the neutral scalars and charged lepton loop [181,182], as shown in Fig. 2 (c), contributing to the AMM corrections. The two-loop correction is where, with z = m 2 i /m 2 φ , and the coefficients may be obtained from Eq. (3.1) as The analytical expressions for the integrals of Eq. (3.6) are given in Appendix A.2.
One-loop: First, we investigate the one-loop contribution to g − 2 of muon and electron from Eq. (3.2). As seen from Eq. (3.2), for degenerate neutral scalars, i.e., m H = m A , the chiral enhancement vanishes, making the effective contribution to ∆a positive. When the couplings are diagonal, a positive contribution to ∆a can also be achieved by considering CP-odd scalar heavier than CP-even scalar, an ideal scenario to explain ∆a µ . However, such an assumption would contradict with ∆a e , for which an overall negative correction would require a lighter pseudoscalar. It is worth mentioning that the negative contribution from the charged scalar H + can partially cancel out the non-chiral part of the neutral scalar corrections. However, the mass of the charged scalar in our choice of the parameter space is always much larger than the mass of the neutral scalar making this effect negligible. Hence, the real Yukawa couplings Y ii (i = e, µ) cannot explain both the AMMs simultaneously from one-loop corrections. Note that if the lepton mediator is different from the external line Y ij (i = j), Yukawa couplings with opposite signs can explain the both the AMMs, subject to constraints from LFV (cf. Sec. 6 and Sec. 5). However, concurrent solutions do exist for Y µµ ∈ R ∩ Y ee ∈ I for m A > m H and vice versa for m A < m H . Note that any nonzero phase in the Yukawa couplings would be constrained by EDMs (cf. Sec. 3.1). The corrections to a can appear from one or more of the following Yukawa couplings as shown: Two-loop: There are several contributions to corrections for AMMs arising at two-loop, which are studied in detail in [181][182][183][184]. A typical but most relevant two-loop Barr-Zee diagram is shown in Fig. 2 (c). The contribution from the Z boson line in the figure is typically suppressed by a factor of O(10 −2 ) compared to the photon line. This is partly due to the massive Z boson propagator and the smallness of its leptonic vector coupling, g cos θ W (− 1 4 + sin 2 θ W ) ∼ −0.015 [185]. Similar diagrams with charged scalars replacing the lepton loop also exist. We suppress, for simplicity, the contribution from this type by considering a small value of quartic coupling λ 7 (H † 2 H 2 H † 1 H 2 ) in the Higgs potential. Note that this quartic coupling does not contribute to the scalar boson masses. Other diagrams that involve W + W − H 1 and H 1 H 2 H 2 couplings also vanish in the alignment limit, α = 0 (see Fig. 3, 5, 6 in Ref. [181]). Moreover, there is no W + W − H 2 and H 2 H 2 H 2 vertex in the Higgs basis where H 2 does not get VEV. The only other non-vanishing diagram, with charged scalars and W ± respectively replacing the neutral scalar and photon lines, involves neutrinos. In this case, the loop factor, which is a function of masses of the internal particles, becomes negligible. Interestingly, this diagram becomes considerably dominant if vector-like fermions are involved [182] instead of neutrinos.
Hence, Fig. 2 (c) is the only two-loop diagram where diagonal Yukawa couplings Y ii are relevant to AMM corrections, as seen from Eq. (3.5) and Eq. (3.7). Moreover, since there is a chiral enhancement in the lepton loop, Y τ τ contribution becomes relevant. It should be noted that taking Y τ τ ∼ 0, one-loop contribution always dominates two-loop, failing to give concurrent solutions for both AMMs. Here we take m H < m A and choose the same sign for Yukawa couplings such that the scalar contribution dominates the pseudoscalar allowing for an overall negative correction to AMMs, thereby providing the right sign to explain ∆a e . ∆a µ > 0 and ∆a e < 0: The two-loop corrections, even though suppressed by α em /π in comparison to one-loop, have a factor of m i /m l from chiral enhancement, as can be seen from Eq. (3.5). This plays an important role in providing the right signs to the AMMs. In the case of ∆a e , the factor m i /m e enhances the two-loop contribution over the corresponding one-loop correction. Therefore, with a choice of heavier pseudoscalar, the overall correction to ∆a e can be made negative. In doing so, one also gets a negative contribution for ∆a µ enhanced by m i /m µ from two-loop, which can become comparable to a one-loop contribution. However, by choosing Y τ τ small enough, one can suppress the two-loop contribution and get the correct sign for ∆a µ from one-loop (see Sec. 8.1 for more details).
As aforementioned, one-loop by itself can explain both AMMs without the necessity to go to the two-loop level by taking the diagonal couplings Y ii ∼ 0. In this scenario, one of the Yukawa couplings can be chosen to be negative so that, for a heavier pseudoscalar, the non-chiral part would explain ∆a µ while the chiral part, enhanced by log (m i /m φ ) (c.f. (A.1)), explains ∆a e .

Electric Dipole Moments
The electric dipole moment of leptons places stringent constraints on the imaginary part of the Yukawa couplings of the scalar field φ. We study these constraints by turning on the relevant couplings such that the two-loop contribution becomes subdominant. These constraints are only significant when there is a chirality flip in the fermion line inside the loop, depicted in Fig. 2 (a). In such a scenario lepton EDM is given by [186] with +(−) corresponding to A (H), and A tiny complex phase of O(10 −5 ) is required to satisfy the current limits on eEDM |d e | ≤ 1.1 × 10 −29 e-cm from ACME [143] while also satisfying ∆a e . We can simply avoid the electron EDM limit by taking all the relevant couplings that give rise to ∆a e real. However, for the case of muon EDM, there exist regions of parameter space where the phase of the Yukawa couplings can be significant and provide enough correction to satisfy ∆a µ while also remaining compatible with the current upper limits from eEDM measurements |d µ | ≤ 1.9 × 10 −19 e-cm [144] and can potentially be measurable in future experiments [145][146][147]. Fig. 3 shows the µEDM ranges that can be probed at the near future experiments and satisfy ∆a µ within 1σ for different values of Y µµ = |Y µµ |e iθ . Here we set all the other Yukawa couplings small for simplicity and to satisfy the flavor constraints. Different color bands (red, green, yellow, blue) represent different choices of Yukawa couplings |Y µµ | (1, 0.1, 0.01, 0.001) by allowing the phase to take arbitrary values. The pink and purple shaded regions are excluded from e + e − → µ + µ − H searches obtained from BABAR [187]  The collider bounds on Y µµ from e + e − → µ + µ − H searches obtained from BABAR [187](pink) and LHC [188](purple) experiments are projected to the maximum EDM that can be obtained for the coupling. The gray band shows the current experimental bound on µEDM, while the dashed and dotted-dashed gray lines provide the sensitivity reach from future experiments [145][146][147].
and LHC [188], respectively. The bounds were obtained in the Y µµ − m H plane, which was then projected onto the maximum values of µEDM, obtained at θ = π/4. The gray region is excluded from current experiments, and the gray dashed [146] and dotted-dashed lines [145,147] are the projected sensitivities from various proposed experiments.
It is worth noting that ∆a µ cannot be satisfied for θ ∈ [ π 4 , 3π 4 ] since the dominant chirally enhanced term in Eq. (3.2) is ∝ cos 2θ which is ≤ 0 for the said region of complex phases regardless of the value of Y µµ and m H .

Non-standard Neutrino Interaction
In the Zee model, the charged scalars η + and H + 2 can induce charged-current NSI at tree level [17]. Since the model can have leptophilic Yukawa couplings Y ie (i = e, µ, τ ) of order unity, significant NSI can be generated. As mentioned before, f 1 due to strong constraints from LFV as well as to correctly reproduce the neutrino oscillation parameters (see Sec. 8). Thus the contributions from f Yukawa couplings are heavily suppressed. Using the dimension-6 operators for NSI [148], the effective NSI parameters in the model can be expressed as where h + and H + are the physical masses of the charged scalars, ϕ is the mixing angle between the scalars and is define in Eq. (2.4) and Eq. (2.5). Since we wish to make the  neutral scalar light to explain AMM of muon and electron, we take a limit when doublet charged scalar is lighter than singlet charged scalar field. In this limit, the contribution from H + dominates, as can be seen, from Eq. (4.1). Note that there is a strong constraint from the T-parameter on the choice of neutral and charged scalar masses and mixing among scalars (see Sec. 5.3 for details). Due to the strong constraints from LFV (cf. Sec. 5.1 and Sec. 5.2), one cannot get sizable off-diagonal NSI ε ij (i = j) in the model. However, sizable diagonal NSI ε ii can potentially be generated by the matrix elements (Y ee , Y µe , Y τ e ), contingent on satisfying ∆a and reproducing neutrino oscillation parameters, which is discussed in the following section.

Low-energy Constraints
In this section, we summarize various relevant low-energy constraints on the parameter space that can potentially explain the observables being explored here. We can safely ignore charged lepton flavor violation (cLFV) involving the f ij couplings as they are tiny to satisfy the neutrino mass constraint. On the other hand, we require Y ∼ O(1) to explain both AMMs and to induce maximum NSI, which leads to various flavor violating processes that are severely constrained by experimental data.
The l 1 → l 2 γ decays are induced radiatively at one-loop diagrams as shown in Fig. 4, and the general expression for such decays involving the neutral scalar field φ reads as [189]  . For each process, there are diagrams with and without chiral enhancement from different charged lepton mediators, f . The first row in the "Constraints" column shows the bound on the diagram with no chiral enhancement, whereas the second row gives the bound on the chirally enhanced diagram contributed by the charged lepton mediator. Table II. Constraints on Yukawa couplings from radiative decay of charged leptons mediated by charged scalar, H − . The constraints here are similar to that of the neutral scalars with no chiral enhancement, given in Table I. where, for The second term in Eq. (5.1) appears from the chirally enhanced radiative diagrams, whereas the first term has no chirality flip in the fermion line inside the loop. The bounds on the Yukawa couplings as a function of the mediator masses are shown in Table I. Table I has two rows showing the constraints arising from the chiral enhancement of charged lepton mediator in addition to the diagram without the enhancement. Similarly, the decay rate for 1 → 2 γ from charged scalar H + in the Zee model can be expressed as Note that there are no chirally enhanced contributions in these decays. The bounds on the Y Yukawa couplings as a function of the mediator masses are shown in Table II.

Process
Exp. Bound [1,192] Constraints Table III. Constraints on Yukawa couplings as a function of neutral scalar mass from trilepton decays l i → l klj l l of charged leptons. Slightly weaker constraints on tau lepton decays from BABAR, ATLAS and LHC can be found in refs. [195][196][197], respectively.

Trilepton Decays
The flavor-changing nature of the new scalar bosons allows for the processes of form l i → l klj l l to realize at the tree level, and it imparts one of the most stringent constraints on the model that precludes simultaneous explanation of AMMs, reducing the flavor structure down to just two textures of Eq. (6.1). The partial rates for such trilepton decays are obtained in the limit when the masses of the decay products are neglected. The decay rate can be read as [16]: δ lk is the symmetry factor that takes care of identical particles in the final state. This expression is relevant for both CP-even and CP-odd neutral scalar mediated decays. Using the total muon and tau decay widths, Γ tot µ = 3.00 × 10 −19 GeV and Γ tot τ = 2.27 × 10 −12 GeV, we calculate the branching ratios for various processes and summarize the constraints on the Yukawa couplings as a function of the neutral scalar masses in Table III.

T-parameter Constraints
The oblique parameters S, T, and U quantify the deviation of a new physics model from the SM through radiative corrections arising from shifts in gauge boson self energies [198][199][200]. Out of these observables, the T-parameter imposes the most stringent constraint. In the decoupling limitα = 0, T-parameter in the Zee model can be written as [201]  where, The allowed regions of scalar masses m A and m H for the choice of mixing angle between charged scalars sin ϕ = 0.1 under the alignment limit are shown in Fig. 5, where the green and blue regions are excluded at 1σ and 2σ from T= 0.03 ± 0.12 [1]. Here we fix the charged scalar mass heavier than neutral scalar mass by choosing m h + = m H + 200 GeV and m H + = m H + 200 GeV such that collider constraints on light charged scalars are easily satisfied [17]. Notice that the choice of ϕ is arbitrary and can be made small. Such a choice only makes the overall scale of neutrino mass smaller and does not alter the phenomenology for ∆a and NSI, as both can be incorporated with just the second Higgs doublet.

Muonium-antimuonium oscillations
The nontrivial mixing between bound states of muonium (M : e − µ + ) and antimuonium M : e + µ − implies a non-vanishing LFV amplitude for e − µ + → e + µ − [202][203][204][205][206][207][208][209]. These oscillation probabilities were measured by the PSI Collaboration, with P (M ↔M ) < 8.3 × 10 −11 at 95% C.L. [203], while MACE Collaboration at CSNS attempts to improve the sensitivity at the level of O(10 −13 ) [210]. These oscillations place a stringent constraint on the product of Yukawa couplings Y eµ and Y µe , thereby excluding a significant portion of parameter space as shown in Fig. 13 of TX-II. The muonium-antimuonium oscillation probability is given by [211,212] P (M →M ) = 64α 6 m 6 red τ 2 where, m red = m e m µ /(m e + m µ ) is the reduced mass between a muon and an electron, α is the QED fine structure constant, and τ µ is muon lifetime. G MM is the Wilson coefficient [213,214] associated with the dimension-six four fermion operator in the effective Hamiltonian density

Direct experimental constraints
In this section, we analyze various direct experimental constraints on the neutral scalar with mass in the range of 100 MeV to 300 GeV range that explains ∆a . There are various experimental constraints one needs to consider, such as dark photon searches, rare Z−decay constraints, and LEP and LHC constraints. The lack of observation of dark photon A d in searches through the e + e − → γA d , with A d → e + e − channel at KLOE [215,216] and BaBar [217,218] sets strong constraints on the couplings. By recasting the results from BaBar and KLOE, one can put a bound on the mass of the light scalars and the corresponding Yukawa couplings, as depicted by the brown and blue shaded region in Fig. 12. Similarly, for the scalar mass m H > 200 MeV, the dark-boson searches at the BaBar [187] can be recast to limit on mass and the Yukawa couplings via the process e + e − → µ + µ − H [219,220], shown as the pink shaded region in Fig. 12. We observe that the searches at BaBar exclude Y µµ 6 × 10 −3 (1.5 × 10 −2 ) at m H = 0.212 GeV (10 GeV). In the relatively heavier mass regime, m H O(10) GeV, Y µµ is constrained by LHC search limits. For example, direct Z searches at the LHC in the pp → µ + µ − (Z → µ + µ − ) channel leads to upper limits on the production cross-section times branching ratio [188]. These upper limits could, in turn, be recast to upper bounds on the Yukawa couplings of our interest. We illustrate its implication on our parameter space as a purple shaded region in the {Y µµ , m H } plane in Fig. 12 (left panel).
The Yukawa couplings of the leptophilic Higgs boson H are also susceptible to constraints from searches at the LEP. Direct searches in contact interaction processes e + e − → + − ( = e, µ) and Higgs production in association with leptons e + e − → + − (H → + − ) can potentially constrain the Yukawa couplings in both TX-I and TX-II.
In the heavy mass limit, LEP data exclude an effective cutoff scale Λ m H /Y ij [221]. The LEP contact interaction constraints on Λ for a light neutral scalar H are no longer applicable. However, due to the t−channel contribution of H/A that interferes with the SM process, the cross-section of e + e − → + α − α can still be modified. By implementing the model file in the FeynRules package [222], we computed the cross-sections using MadGraph5_aMC@NLO [223], which are then compared with the measured cross-sections [221,224] to get a limit on Yukawa coupling as a function of scalar mass; Y ee < 0.8 and Y eµ < 0.74 for the benchmark value of m H = 130 GeV [17]. This limit is slightly weaker than the bounds from searches in the e + e − → + − (H → + − ) channel, discussed later in great detail.
The most constraining Higgs production processes for Y ee and Y eµ are e + e − → e + e − (H → e + e − ) and e + µ − (H → e + µ − ), respectively. Likewise, Y µµ would be most sensitive to the e + e − → µ + µ − (H → µ + µ − ) process. However, the production cross-section of the latter at LEP is roughly an order of magnitude smaller than its former two counterparts. Therefore, we ignore the constraints on Y µµ from direct searches in the e + e − → µ + µ − (H → µ + µ − ) channel for the present analysis. In order to correctly estimate the LEP bounds on Y ee (Y eµ ), we generate signal events in the e + e − → e + e − (H → e + e − ) (e + e − → e + µ − (H → e + µ − )) channel assuming √ s = 207 GeV for several values of m H at the leading order (LO) using the MadGraph5_aMC@NLO [223] package. We also simulate the respective dominant background processes: 4e (2e2µ) in the same framework. We reconstruct the Higgs boson by identifying the e + e − (e + µ − ) pair with the smallest ∆R = ∆η 2 + ∆φ 2 separation, where ∆φ is the azimuthal angle difference between the opposite sign lepton pairs. Consequently, we compute the signal and background efficiencies and translate them into upper bounds on Y ee (Y eµ ) as a function of m H . We illustrate the current upper limits from LEP on Y ee and Y eµ as grey shaded regions in the right panel of Fig. 12 and Fig. 13, respectively. We would like to note that we do not include any LEP detector effects in these studies, and therefore our "LEP estimations" must be treated cautiously. As such, these upper limits are rather conservative estimations.
In principle, the Higgs boson could have been reconstructed by identifying the opposite sign lepton pairs with invariant mass closest to m H . However, such a reconstruction strategy leads to an almost 100% signal efficiency for all choices of m H since we are restricted to a simplistic truth level analysis without accounting for detector effects. This could potentially entail considerably stronger upper limits on the Yukawa couplings than a realistic scenario leading to misleading implications. We note that we have adopted this reconstruction strategy in Sec. 7, where we perform a detailed cut-based collider analysis at the detector level to study the future sensitivity on Y ee , Y eµ and Y µµ from future lepton colliders.
For completeness, we would like to note that the charged scalar in the model can be pair produced directly at colliders via s-channel off-shell photon or Z exchange with H ± further decaying into ν. The corresponding LEP bound is m H ± 80 GeV [225] for τ ν final state. Moreover, these leptonic final states ν mimic slepton searches in supersymmetric models that can be recast for our scenario and provide a lower bound on its mass m H ± 100 GeV [17]. Similarly at the LHC, the s-channel Drell-Yan process pp → γ * /Z * → H + H − provides a somewhat stronger limit of m H ± 200 GeV [17] for BR(H ± → eν) = 1. Note that the LHC limits can be evaded by lowering the branching ratio, which can be always be done with more than one Yukawa coupling.

Flavor Structures
This section explores the texture of the Yukawa coupling matrix Y for a unified explanation of electron and muon g − 2 while correctly reproducing neutrino oscillation data. The Yukawa coupling f appearing in the neutrino mass formula given in Eq. (2.7) is taken to be arbitrary and small such that it automatically satisfies all the constraints and generates the correct order of the neutrino masses. However, the Yukawa matrix Y has non-trivial structures and its various components explain ∆a e/µ as previously discussed in Sec. 3, for instance, Eq. (3.8) for one-loop texture. It turns out that the Yukawa couplings Y ie (i = e, µ, τ ), which is a part of the texture to explain ∆a , also induce NSI at some level. A cursory glance at the analysis of LFV processes in Sec. 5 reveals two types of Yukawa textures to accommodate both AMMs in conjunction with neutrino observables, which are: These textures are studied in detail in Sec. 8. Here the green, blue, and red color-coded entries respectively represent ∆a µ , ∆a e , and NSI. The same couplings that can explain more than one observables are encircled with the corresponding colors. The ×'s denote the Yukawa couplings required to satisfy the five neutrino oscillation observables (∆m 2 21 , ∆m 2 31 , sin 2 θ 13 , sin 2 θ 23 , sin 2 θ 12 ) while satisfying the flavor constraints such as i → j γ and trilepton decay (c.f. Sec. 5). For more details, see Sec. 8. We also note that the zeros in the matrices of Eq. (6.1) need not be exactly zero, but they need to be sufficiently small so that the flavor-changing processes remain under control (cf. Sec. 5). Note that the flavor structures are proposed under the assumption that m H is smaller than all other new scalar bosons.

TX-I:
Here, we examine the texture where both ∆a e and ∆a µ could be explained from the diagonal Yukawa couplings Y ii . For the choice of real Yukawa couplings, Y µµ can explain ∆a µ at one-loop order parametrized by the chirally enhanced term of Eq. (3.2). The one-loop contribution from Y ee alone always gives the wrong sign to ∆a e . However, the inclusion of a third Yukawa coupling Y τ τ , appearing from the two-loop Barr-Zee diagram, can generate the correct sign for ∆a e [141]. Furthermore, the same Yukawa coupling Y ee induces NSI, as previously discussed in Sec. 4. On the other hand, for complex Yukawa couplings, concurrent solutions do exist, for instance, Y µµ ∈ R ∩ Y ee ∈ I for m A > m H and vice versa for m A < m H . Note that any nonzero phase in the Yukawa couplings would be constrained by EDMs (cf. Sec.3.1).

TX-II:
The only other flavor structure to incorporate ∆a e and ∆a µ while satisfying neutrino oscillation data is given in Eq. (6.1) as TX-II with Yukawa couplings Y eµ and Y µe . With non-zero Y eµ (e ↔ µ), the diagonal couplings Y ii (i = e, µ) are highly constrained from µ → eγ (cf. Sec. 5.1). Thus, the two-loop contributions are highly suppressed and can be safely ignored. To explain both the anomalies, Y eµ and Y µe can take opposite signs to get ∆a e negative. Moreover, considering a hierarchy among the two Yukawa couplings, one can explain ∆a µ from the non-chiral part of Eq. (3.2). Details on the choice of the Yukawa couplings as a function of the mass of scalar field is given in Fig. 13 in Sec. 8. Note that Y µe ∼ O(1) can induce NSI. However, such a choice necessarily requires |Y µe | > |Y eµ | to get the correct order for ∆a , which turns out to be not in favor with the normal hierarchy (NH) solution to the neutrino oscillation data (cf. Sec. 8). For the sake of completeness, we also point out that Y τ e can induce large NSI, τ τ ∼ 9.3% [17], and provide correction to a e from tau mass chiral enhancement involving the coupling Y eτ . As was the case before, the ∆a e emerges from one-loop diagrams alone. However, there exists no choice of Yukawa couplings that can incorporate ∆a µ in this scenario; LFV processes highly suppress all the couplings that provide corrections a µ .

Collider Analysis
The Yukawa structure of the new Higgs boson H required to simultaneously explain ∆a , NSI as well as neutrino oscillation data entails concurrent non-zero values of all three diagonal elements {Y ee , Y µµ , Y τ τ } or the first and second generation off-diagonal entries {Y eµ , Y µe }, as discussed in Sec. 6. At lepton colliders, these couplings could be directly accessed via searches in the e + e − → + − (H → + − ) channel, where = e, µ. In this section, we explore the sensitivity of the aforesaid channels to probe Y ee and Y eµ at the projected ILC configuration: √ s = 1 TeV with an integrated luminosity L = 500 fb −1 [151][152][153][154][155]. We also study the projected sensitivity on Y µµ from direct searches in the µ + µ − → µ + µ − (H → µ + µ − ) channel at the future muon collider (MuC) with configuration: √ s = 3 TeV, L = 1 ab −1 [156,157].
We generate signal and background events at leading order (LO) using the MadGraph5_aMC@NLO [223] framework. Showering and hadronization is performed with Pythia-8 [226,227] while fast detector response is simulated with Delphes-3.5.0 [228]. We utilize the default ILCgen [229] and MuonCollider [230] detector cards to simulate the response of the ILC and the muon collider.

Y ee at ILC
We study the projected sensitivity on Y ee at the ILC using the channel e + e − → e + e − H → e + e − H → e + e − . (7.1) A few typical leading order (LO) Feynman diagrams of this process are illustrated in Fig. 6.
The dominant SM background is e + e − → e + e − e + e − . We select events containing exactly two isolated electrons and two isolated positrons in the final state with p T > 2 GeV and |η| < 3.0. We further impose p T,e 1,2 > 10 GeV and p T,e 3,4 > 5 GeV, where e 1 and e 4 are the highest and lowest p T leptons, respectively. In Fig. 7, we illustrate the p T distribution of the highest and second-highest p T leptons, e 1 and e 2 , respectively, for three signal benchmarks, M H = 10, 130 and 230 GeV, and the 4e background process. In the signal process, the Higgs boson can recoil against an electron or a positron as seen in the first two diagrams in Fig. 6. This leads to an overall improvement in the p T of the recoiling electron with increasing M H . At relatively large m H , this recoiling electron becomes the dominant constituent in e 1 . Correspondingly, we expect to see an upward shift in the peak of p T,e 1 distribution with increasing m H . We observe this behaviour in Fig. 7 where p T,e 1 peaks at ∼ 60 GeV in the M H = 10 GeV scenario and the peak position shifts to p T,e 1 ∼ 90 and ∼ 130 GeV in the m H = 130 and 230 GeV scenarios, respectively. Furthermore, we observe that the overall distribution gets flatter with increasing m H . The background p T distribution of e 1 and e 2 peaks at s ∼ 90 and ∼ 60 GeV, respectively, which roughly coincides with the peak in the M H = 130 GeV scenario. The 4e background process also includes diagrams where an on-shell Z boson recoils against an electron or positron leading to the aforesaid similarity in peak positions. Next, let us focus on the pseudorapidity distributions of the final state leptons. We present the η e 1 and η e 2 distributions in Fig. 7. We observe that e 1 and e 2 in the background are mostly produced in the forwards regions of the detector due to back-to-back production of electron-positron pairs. On the contrary, the leading and sub-leading p T leptons in signal benchmarks with large m H (∼ 230 GeV) are mostly produced in the central regions |η| 1.0 by virtue of their larger transverse momenta. At relatively smaller Higgs masses M H ∼ 10 and 130 GeV, the peaks in pseudorapidity distribution roughly coincides with that of background.   Taking cognizance of the m Hrec distributions in Fig. 8, we require signal and background events to satisfy m H ± 10 GeV. In addition, we also optimize the selection cuts on p T,e 1 ,η e 1 and η e 2 for several signal benchmarks in order to maximize the signal significance σ S = S/ √ S + B, where S and B are the signal and background yields. In Table IV, we present the optimized selection cuts, signal and background yields, along with the signal significance, for several signal benchmarks within the range 10 GeV ≤ M H ≤ 350 GeV. We utilize the optimized signal and background efficiencies obtained in the previous step to derive the projected sensitivity on Y ee as a function of M H at the ILC 1 TeV machine. We present our results as upper limit projections at 2σ and 5σ in the {Y ee , M H } plane in Fig. 8.
We observe that the projected upper limits are weaker at small m H ∼ 10 GeV, and improves by a factor of ∼ 3 until m H ∼ 50 GeV. Afterwards, it gradually improves with increasing m H . At a particular {Y ee , m H }, the projected sensitivity is determined by three factors: the signal production cross-section ∝ Y ee and m H , signal efficiency ∝ m H and background efficiency. At the ILC 1 TeV run, the LO production cross-section of the signal process in Eq. 7.1, σ Yee , peaks at roughly ∼ 100 fb with m H ∼ 220 − 240 GeV for Y ee = 1.0, and falls down to ∼ 7 fb (∼ 70 fb) at m H ∼ 10 GeV (350 GeV), respectively. On the other hand, we observe that the signal efficiency in Table IV improves  GeV. At small m H ( 70 GeV), the small signal production cross-section coupled with a relatively small signal efficiency leads to weaker sensitivity on Y ee . In the 70 GeV m H 230 GeV region, both signal production cross-section and efficiency improves while the background efficiency deteriorates. All these factors contribute towards improving the projected sensitivity. As we move to signal benchmarks with heavier m H , the signal production cross-section registers a decrement. However, the signal efficiency continues to increase until m H = 310 GeV and the background efficiency continues to plummet. The later two counters the decrements in signal cross-section, and the projected upper limits continue improving. At m H = 350 GeV, the signal efficiency is O(8%) smaller than its m H = 310 GeV counterpart. However, this decrement in signal efficiency is not reflected in the projected upper limits in Fig. 8 since the background efficiency falls down by a relatively larger rate O(40%). Overall we observe that the 1 TeV ILC machine would be able to probe Y ee up to Y ee 0.085 (at 2σ) at m H = 100 GeV via direct searches in the e + e − → e + e − (H → e + e − ) channel.

Y eµ at ILC
To explore the projected sensitivity to Y eµ , we consider the channel We consider events containing an OS electron pair and an OS muon pair. All four leptons ( = e, µ) are required to have p T > 2 GeV and |η| < 3.0. The leptons are also required to satisfy p T 1,2 > 10 GeV and p T 3,4 > 5 GeV, where 1 and 4 are the highest and lowest p T leptons, respectively. The kinematic behavior of the final state leptons are similar to that exhibited in the e + e − → e + e − (H → e + e − ) channel. One major difference however is the absence of t-channel electron exchange diagrams, shown in the right panel of Fig. 6, since the H couples to e ± µ ∓ pair instead of e + e − pair. In Fig. 9, we present the p T and η distributions of l 1 and l 2 for three signal benchmarks m H = 10, 130 and 240 GeV, and the 2e2µ background. We observe that the overall features are roughly similar to that in Sec. 7.1. Here as well, the peak of the p T distributions shifts to larger values with increasing m H . Similarly, the relative concentration of 1 and 2 in the central regions of the detector improves with m H .  Table V. Optimized selection cuts on p T,e1 , η e1 and η e2 , signal efficiency and background efficiency, from cut-based analysis in the e + e − → µ + e − (H → µ + e − ) channel at √ s = 1 TeV ILC for several signal benchmark points.
In order to reconstruct the Higgs boson, we perform mass minimization similar to that in Sec. 7.1. We associate an OS electron-muon pair with H that minimizes (m 2 e ± µ ∓ − m 2 H ). In Fig. 10, we present the invariant mass distribution of the reconstructed Higgs boson for several signal benchmarks and the 2e2µ background. We impose |m Hreco ± 10 GeV| to filter the signal from the 2e2µ background continuum. Additionally, we also optimize selection cuts on p T, 1 , η 1 and η 2 . We present the optimized signal and background efficiency for several choices of m H along with the respective set of optimized cuts in Table V. Both signal and background efficiencies follow a similar behaviour to that in Sec. 7.1. We observe that the signal efficiency improves by a factor of ∼ 40 from falling down to ∼ 0.627 at m H = 350 GeV. We translate these efficiencies into upper limit projections in {Y eµ , m H } plane, shown in Fig. 10. We observe that the 1 TeV ILC machine would be able to probe Y eµ up to ∼ 0.1 and ∼ 0.09 at m H = 100 and 200 GeV, respectively, at 2σ. We note that the projected sensitivity on Y eµ is relatively weaker compared to Y ee (see Sec. 7.1) due to relatively smaller signal production cross-section in the present channel.

Y µµ at future muon collider
The e + e − → µ + µ − (H → µ + µ − ) channel offers a direct probe to the Y µµ coupling at electron-positron colliders. However, the production cross-section of this process is at least an order of magnitude smaller relative to e + e − → e + e − (H → e + e − ) process due to the absence of t-channel Z/γ * exchange diagrams, leading to weaker sensitivity. The issue of smaller production cross-section could be, however, circumvented at a muon collider (MuC). Muons produce less synchrotron radiation relative to the electrons, and therefore, a MuC machine could be operated at a much higher center of mass energies compared to an electron-positron collider. Consequently, the most stringent sensitivity on muon Yukawa coupling via direct searches in the + − H channel could perhaps be achieved at a future MuC. Taking this into consideration, in the present section, we explore the potential capability of a muon collider with configuration We require events to contain exactly four isolated muons in the final state. We assume detector coverage up to |η| < 2.5, and also impose p T,µ > 10 GeV. The final state muons exhibit kinematic features that are similar to that of the final state electrons in Sec. 7.1. Consequently, we adopt the mass minimization strategy from Sec. 7.1 to reconstruct the Higgs boson. We do not present the m Hreco distribution in the present scenario due to its close similarity with Fig. 8. The dominant background source is the µ + µ − → 4µ process. We impose additional selection cuts on the p T of the final state muons in order to improve signal-background discrimination: p T,µ 1 > 90 GeV p T,µ 2 > 20 GeV, p T,µ 3 > 15 GeV, (7.3) where, µ 1 represents the highest p T muon. We tabulate the signal and background efficiencies for various signal benchmarks corresponding to different values of m H in Table. VI. In the Higgs mass range of our interest 10 GeV < m H < 350 GeV, we observe that the signal efficiency improves with increasing m H . The background efficiency, on the other hand, falls down with m H except for m H close to m Z where the leptons from Z resonance fill into the m Hreco distribution. We utilize these efficiencies to derive projected upper limits on Y µµ as a function of m H . We illustrate the 2σ and 5σ projections in Fig. 11. The signal efficiency improves by more than one order of magnitude from ∼ 0.006 at m H = 10 GeV to ∼ 0.151 at m H = 40 GeV, while the signal production cross-section improves by a factor of ∼ 6. The background efficiency, on the other hand, remains almost unchanged. This leads to an order of magnitude improvement in the projected upper limits on Y µµ . At m H = 100 GeV, the background efficiency is roughly 4 times higher than at its neighbour signal benchmark points in Table. VI due to its closeness to the Z resonance. This translates into a weakening in the projected upper limits in the vicinity of m Z as seen in Fig. 11. Above m Z , we observe that the signal efficiency continues to grow while the background efficiency keeps falling gently with increasing m H . The signal production cross-section also continues to rise. All these factors lead to a gradual strengthening in the projected upper limits on Y µµ with increasing m H . At m H = 300 GeV, we observe that MuC would be able to probe Y µµ up to ∼ 0.06 at 2σ.

Results and Discussions
In this section, we present numerical analysis for the model parameter space and reconcile electron and muon g − 2 within their 1σ measured values while being consistent with all the low-energy, LHC, and LEP constraints discussed in the previous sections. After exhausting all the possibilities, we find two minimum textures discussed in Eq. (6.1) of Sec. 6 to incorporate both of these anomalies and have a consistent neutrino oscillation fit (∆m 2 21 , ∆m 2 31 , sin 2 θ 13 , sin 2 θ 23 , sin 2 θ 12 ). The neutrino mass matrix given in Eq. (2.7) is diagonalized by a unitary transformation where M ν is the diagonal mass matrix, and U is the 3 × 3 PMNS lepton mixing matrix. We diagonalize the mass matrix numerically by scanning over the input parameters while being consistent with ∆a and LFV constraints. For the ease of satisfying the flavor constraints on the Yukawa coupling f , we factor out f µτ into the overall factor and define a 0 = κf µτ , where κ is the one-loop factor given in Eq. (2.8). Moreover, we perform a constrained minimization where the observables are confined to 3σ of their experimental measured values. The fits to the two textures discussed in Eq. (6.1) are shown in the subsequent sections. It is beyond the scope of this work to explore the entire parameter space of the theory; instead, we find benchmark points to show that the model is consistent with neutrino oscillation data while explaining the anomalies for both textures. As was eluded to earlier, we choose the masses of the scalar bosons such that the major contributions to the AMMs appear from the

Fit to TX-I
The allowed parameter space of the flavor structure of TX-I in Eq. (6.1) is explored here, as shown in Fig. 12. The green and yellow bands in both plots correspond to 1σ and 2σ regions that can simultaneously explain both the AMMs. Here we fix Y τ τ = 0.05 such that it provides the required sign for ∆a e from the Barr-Zee diagram. Note, a smaller value of the Y τ τ requires a larger Y ii (i = τ ) to satisfy AMMs, excluding more parameter ranges. On the other hand, making Y τ τ larger does allow wider parameter space as the plot in Fig. 12 would shift downwards. However, it conflicts with the fit to neutrino oscillation data as there would be a large hierarchy in the elements of the neutrino mass matrix. The various shaded regions in Fig. 12 are excluded from various experimental constraints. The pink and purple shaded regions in Fig. 12 are excluded from e + e − → µ + µ − H searches at BaBar [187] and LHC [188]. Here we considered BR(H → µµ)=1 (dotted line) and 0.5 (purple shaded region) to show that more parameter space is allowed for BR < 1. Blue and brown regions in the right plot are exclusion regions from the dark photon searches through e + e − → γA d channel at BABAR [187] and KLOE [215]. The combination of these constraints on Yukawa couplings Y ee and Y µµ would exclude light scalar mass below 10 GeV in the parameter space of our interest. The black region is obtained from LEP [221] constraints through e + e − → e + e − (H → e + e − ) searches, and the black dash-dotted line pink and purple regions are obtained from e + e − → µ + µ − H searches at BABAR [187] and LHC [188], respectively; blue and brown regions from dark photon searches through e + e − → γH at BABAR [187] and KLOE [215]. Black shaded region is excluded from e + e − → e + e − H searches at LEP [221,224]. Blue and pink dash-dotted are the projected sensitivities from dark-photon searches at Belle-II [231,232] through e + e − → γH and e + e − → µ + µ − H searches, respectively. The black dash-dotted line on the left plot shows the projected sensitivity reach on Y µµ from direct searches in the µ + µ − → µ + µ − (H → µ + µ − ) channel at the 3 TeV Muon Collider assuming L = 1 ab −1 (cf. Fig 11), and the one on the right plot shows the projected upper limits on Y ee from direct searches in the e + e − → e + e − (H → e + e − ) channel at the 1 TeV ILC assuming L = 500 fb −1 (cf. For the Yukawa texture above, the corresponding fit for the neutrino oscillation parameters and ∆a are shown in Table VII  . For comparison, the 3σ allowed range for the oscillation parameters and 1σ for ∆a are also given. Figure 13. The parameter space of the Yukawa couplings Y eµ Y µe (left) and Y eµ (right) as a function of scalar mass (m H ) satisfying both AMMs. The green and yellow band in both the figures respectively represents 1σ and 2σ allowed region. The region shaded in gray is excluded by the muonium oscillation PSI experiment [203], whereas the dash-dotted line represents the future sensitivity bound from the MACE experiment [210]. Black shaded region is excluded from e + e − → e ± µ ∓ H searches at LEP [221,224], whereas blue dash-dotted line is the projected sensitivity from ILC with √ s = 1 TeV with integrated luminosity L = 500 fb − 1.
τ → 3µ. Thus, for a larger value of Y τ τ , say 0.05, Y µµ 0.8 is required to get NH solution.
For such a large Yukawa coupling, AMM can be satisfied by increasing the scalar mass to about 300 GeV. As noted in Sec. 4, order one coupling Y ee can also induce NSI; for the benchmark point Y ee = 0.31, ϕ = 0.1, and m H + = 285 GeV (∆m = 200 GeV), we find ε ee = 4.2%. This can be improved up to 8% [17] with a proper choice of parameters and is limited by direct experimental searches TEXONO [233].

Fit to TX-II
The parameter space explored in TX-II of Eq. (6.1) is provided in this section. As aforementioned, off-diagonal couplings Y eµ and Y µe need to take opposite signs to get a negative ∆a e from the chirally-enhanced part of Eq. (3.2). Furthermore, one of the two couplings needs to be large to provide the required correction to ∆a µ via the non-chiral part of Eq. (3.2). The allowed parameter space as a function of Yukawa couplings and scalar mass is shown in Fig. 13, where the green and yellow bands correspond to 1σ and 2σ regions of ∆a µ/e . The same couplings that give rise to ∆a e/µ also induce muonium oscillations, and the probability of these oscillations at the PSI experiment [203] excludes a considerable region of the allowed parameter space, with a lower bound on the scalar mass of 8 (1) GeV at 1σ (2σ). These bounds are expected to improve at the MACE experiment [210] that can probe/exclude all the allowed parameter space. These couplings are directly accessible at lepton colliders via searches in the e + e − → e ± µ ∓ (H → e ± µ ∓ ) channel and can be used to obtain the bound from LEP [221] as shown by the black shaded region in Fig. 13. The LEP bound imposes an upper limit of about 30 GeV on the scalar mass with the simultaneous explanation of ∆a e/µ . Furthermore, the projected sensitivity of Yukawa coupling Y eµ from ILC is shown in the dash-dotted line, which can probe Y eµ up to ∼ 0.1 at m H = 100 GeV and explore the parameter space for AMMs as low as 20 GeV scalar mass. Our results for the fit to the TX-II of Eq. (6.1) is shown below: The corresponding fit for the neutrino oscillation and ∆a are shown in Table VII as Model Fit II. The fits are in excellent agreement with the observed experimental values. Here the non-zero Yukawa couplings are introduced to get the neutrino oscillation fit, whereas 0's in the texture are to satisfy LFV such as µ → eγ and τ → eµµ. It is worth pointing out that although it is ideal to choose Y µe ∼ O(1) instead of Y eµ such that it can induce diagonal NSI, ε µµ , there is no solution to neutrino oscillation data while simultaneously explaining ∆a and satisfying flavor constraints.

Conclusion
In conclusion, we have proposed a novel scenario in a radiative neutrino mass model, namely the Zee Model, that solves the intriguing deviations observed in muon and electron AMMs simultaneously within 1σ uncertainty, while being consistent with neutrino oscillation data, as well as all the relevant flavor and collider constraints. The neutral scalar in the second Higgs doublet is responsible for explaining both of these anomalies via one-loop and two-loop contributions. The two-loop Barr-Zee diagram is essential for getting the negative correction to a e . Furthermore, by exhausting all the possibilities, we find two minimum textures in the model that can incorporate both of these anomalies and have a consistent neutrino oscillation fit. We observe that the currently allowed parameter space accommodates the scalar mass in the range of roughly 10-300 GeV in TX-I and 1-30 GeV in TX-II. In TX-I, for m H > 200 GeV (m H > 40 GeV), direct searches in the e + e − → e + e − (H → e + e − ) channel at the 1 TeV ILC would be able to probe all allowed parameter space points with Y ee 0.07 (Y ee > 0.1) at 2σ uncertainty. Similarly, through searches in the µ + µ − → µ + µ − (H → µ + µ − ) chanel, the 3 TeV MuC collider would be able to probe the entire region with m H 40 GeV irrespective of Y µµ . In TX-II, we observe that all allowed parameter space points with m H > 20 GeV fall within the projected exclusion reach of 1 TeV ILC via searches in the e + e − → µ + e − (H → µ + e − ) channel. The model also features charged scalars required for generating neutrino masses which induce diagonal NSI ε ee as large as 8%. Furthermore, the model can also give a sizable EDM of muon that can potentially be measured in future experiments. This task of connecting both the AMMs and finding textures to get neutrino oscillation fit is a highly non-trivial. As we have shown in detail, the allowed parameter space in the model can be probed in future experiments, which will either discover NP or rule out the possibility of explaining both AMMs within the context of the Zee Model.

A.1 One-loop
In the limit m i m φ the expression for one-loop contribution to AMM [235] from neutral scalars simplifies to ∆a (1) In the expression above, − and + correspond to H and A, respectively. Similarly, the charged scalar contribution is simply