An explanation for the muon and electron $g-2$ anomalies and dark matter

We propose simple models with a flavor-dependent global $U(1)_\ell$ and a discrete $\mathbb{Z}_2$ symmetries to explain the anomalies in the measured anomalous magnetic dipole moments of muon and electron, $(g-2)_{\mu,e}$, while simultaneously accommodating a dark matter candidate. These new symmetries are introduced not only to avoid the dangerous lepton flavor-violating decays of charged leptons, but also to ensure the stability of the dark matter. Our models can realize the opposite-sign contributions to the muon and electron $g-2$ via one-loop diagrams involving new vector-like leptons. Under the vacuum stability and perturbative unitarity bounds as well as the constraints from the dark matter direct searches and related LHC data, we find suitable parameter space to simultaneously explain $(g-2)_{\mu,e}$ and the relic density.


I. INTRODUCTION
The Standard Model (SM) for elementary particles has successfully explained a plethora of phenomena in various experiments. Despite its tremendous success, physics beyond the SM (BSM) is strongly called for to explain neutrino oscillations, dark matter (DM) and baryon asymmetry of the Universe that cannot be accommodated within the SM. The question is then how we can experimentally show the existence of such a new physics model. A discovery of new particles, of course, would provide a direct proof. However, no report of such discoveries has been given so far, though there is still a possibility for their detection in future collider experiments, such as the High-Luminosity LHC [1] and the Future Circular Colliders (FCCs) [2]. In addition to the direct searches, precision measurements of certain observables can also offer good opportunities to probe new physics (NP). Deviations in measured values of the observable from their SM predictions can be attributed to the effects of new particles.
Among various observables, the anomalous magnetic dipole moment of the muon, dubbed the muon g − 2, has long been thought to be a harbinger for NP [3] and attracted a lot of attention for almost two decades because of the discrepancy between its experimental value measured at Brookhaven National Laboratory (BNL) [4] and the SM expectation. According to recent studies about the hadronic vacuum polarization contributions [5][6][7][8] to the muon g − 2, the discrepancy is at about 3.3σ level [9], with the experimental value higher than the SM prediction. See also the recent review on the muon g − 2 [10]. On the other hand, the experimental value of the electron g − 2 has been updated in 2018 [11] from a precision determination of the fine-structure constant α em . Interestingly, this measurement also shows a possible disagreement between the data and theory, with the measured value lower than the SM prediction by about 2.4σ [11]. These tantalizing opposite deviations have invited many studies to explore suitable NP models [12][13][14][15][16][17][18][19][20].
In order to accommodate both g − 2 anomalies simultaneously, a characteristic flavordependent structure is called for. In this paper, we propose a new model with a set of new particles whose interactions are constrained by a flavor-dependent global U (1) symmetry and a Z 2 symmetry, and demonstrate its capabilities to simultaneously accommodate both anomalies and, at the same time, offer a DM candidate. These new symmetries do not only play an important role in explaining both anomalies, but also forbid dangerous flavor-violating decays of the charged leptons, such as µ → eγ. Furthermore, they also guarantee the stability of the DM candidate, which is the lightest neutral particle among the new particles. We find regions in the parameter space that can satisfy the relic density and the direct search constraint of the DM while successfully explaining both g − 2 anomalies. This paper is organized as follows. In Sec. II, we define our model and give the Yukawa interactions and the scalar potential that are compliant with the symmetries. In Sec. III, we discuss the new contributions to the muon and electron g − 2, and scan the parameter space for regions that can explain both anomalies. Sec. IV is devoted to the discussion on DM physics and the collider phenomenology. Our conclusion is summarized in Sec. V.

II. MODEL
In addition to the SM gauge symmetry SU (2) L ⊗ U (1) Y , our model has an additional global U (1) and an exact Z 2 symmetries. The particle content in the lepton and scalar sectors is given in Table I 1 . The lepton sector is comprised of new vector-like isospin singlets χ a (a = e, µ) in addition to the SM left-(right-) handed lepton doublets (singlets) L L ( R ) with = e, µ, τ . The scalar sector is also extended from the SM one by introducing additional scalar isospin doublet η D and singlet η S fields, with the SM Higgs doublet field denoted by Φ. All of the new fields (χ a and η D,S ) are assigned to be odd under the Z 2 symmetry. In Table I, the hypercharge Y D is chosen to be either 0 or 1 in order to include at least one neutral particle in the Z 2 -odd sector to be a DM candidate, provided it is the lightest among all the Z 2 -odd particles. For simplicity, we assume η S to be a real field for the scenario with Y D = 1.
The Z 2 -even scalar doublet field is parameterized as usual as The U (1) charges depend on the lepton flavor with q e = q µ . The parameter Y D appearing in the hypercharges for Z 2 -odd particles can be either 0 or 1.
while the Z 2 -odd scalar doublet can be parameterized as where G ± and G 0 are the Nambu-Goldstone bosons that are absorbed as the longitudinal components of the W ± and Z bosons, respectively. The vacuum expectation value (VEV) v is fixed by The lepton Yukawa interactions and the mass term for χ a are given by where ( e R , µ R , τ R ) = (e R , µ R , τ R ). Because of the U (1) symmetry, we can naturally realize the flavor-diagonal couplings f L and f R , so that contributions from the new particles to lepton flavor-violating processes such as µ → eγ can be avoided at all orders. It should be emphasized here that analogous to the GIM mechanism, this structure cannot be achieved in a model with only one vector-like lepton, where it is impossible to accommodate both muon and electron g − 2 while suppressing the µ → eγ decay to the level consistent with the current experimental bound. In general, the new Yukawa couplings f a L,R can be complex, but we assume them to be real for simplicity in the following discussions. The Lagrangian for the quark and gauge sectors are the same as in the SM.
The most general form of the scalar potential consistent with all the symmetries is given by where with τ 2 being the second Pauli matrix. The phases of λ 5 and κ parameters can be removed by a redefinition of the scalar fields without loss of generality. Therefore, CP symmetry is preserved in the scalar potential. We require µ 2 Φ , µ 2 D , µ 2 S > 0 in order to preserve the stability of the SM vacuum.
The squared mass of the Higgs boson h is given by m 2 h = v 2 λ 1 in both scenarios of Y D = 1 and Y D = 0. On the other hand, the mass formulas for the Z 2 -odd scalar bosons are different in the two scenarios. For the scenario with Y D = 1, the singlet field η S is neutral (η 0 S ≡ η S ), so that the η 0 H and η 0 S fields can mix with each other. By introducing a mixing angle θ, the mass eigenstates of these neutral scalar fields can be defined through  where s θ ≡ sin θ and c θ ≡ cos θ. The mixing angle can be expressed as where M 2 H is the mass matrix in the basis of (η 0 H , η 0 S ): The squared masses of the scalar bosons are then given by From the above expressions, we can write the parameters in the scalar potential in terms of the physical parameters as follows: , After fixing m h and v to their experimental values, the remaining ten independent parameters in the scalar potential are then chosen to be and the quartic couplings λ 2,6,8 for the Z 2 -odd scalar bosons.
For the scenario with Y D = 0, the singlet field η S is singly-charged (η ± S ≡ η S ), so that the charged components of the inert doublet field η ± can mix with η ± S . Similar to the above scenario, the mass eigenstates are defined through  with The mass matrix M 2 ± is expressed in the basis of (η ± , η ± S ) as The squared masses of the scalar fields are then given by Some of the parameters in the potential can be rewritten in terms of the physical parameters as , Therefore, the ten independent parameters in the scalar potential can be chosen as and the quartic couplings λ 2,6,8 for the inert scalar fields.
The parameters in the scalar potential are subject to the constraints of perturbative unitarity and vacuum stability. In order for our models to be perturbative, we require all the quartic couplings λ i in the potential to satisfy To impose the tree-level unitarity constraints, we consider all possible 2 → 2 elastic scatterings for the bosonic states in the high energy limit, and obtain thirteen independent eigenvalues of the s-wave amplitude matrix, expressed in terms of the scalar quartic couplings. By demanding the magnitude of each eigenvalue to be smaller than 8π [22], we find the following conditions for the quartic couplings 2 ; where a 1,2,3 are the eigenvalues for the following 3 × 3 matrix  with the coefficients (c 1 , c 2 , If we take λ 6,7,8 = 0, the above expressions are reduced to those in the two-Higgs doublet model (see, e.g., Ref. [24]).
To ensure the stability of the SM vacuum, besides requiring the quadratic terms µ 2 D and µ 2 S to be positive, we further require the potential to be bounded from below. The bounded-from-below conditions are given by [23] where For the convenience of discussions, we define the scalar trilinear coupling λ φ 1 φ 2 φ 3 to be the coefficient of the φ 1 φ 2 φ 3 term in the Lagrangian, where φ i are the physical scalar bosons in our model.
Before closing this section, we briefly comment on neutrino masses in our model. Under the charge assignments given in Table I, the structure of the dimension-5 operator is strongly In order to obtain nonzero values for all the elements of the 3 × 3 neutrino mass matrix for the observed mixing pattern, two additional Higgs doublet fields, denoted by Φ e and Φ µ , are required. Taking the U (1) charge for Φ e and Φ µ to be −q e and −q µ , respectively, we can write down all the dimension-5 effective Lagrangian as where Φ τ = Φ, and c ij and Λ are respectively dimensionless couplings and the cutoff scale.
Note that if we consider the case with one of the three Higgs doublets being absent, the neutrino mass matrix has the texture with three zeros; that is, one diagonal and two offdiagonal elements including their transposed elements are zero. It has been known that such textures cannot accommodate the current neutrino oscillation data [25]. Hence, at least three Higgs doublets are required. In the following discussions, we consider the model defined with just the Higgs doublet in Table I by assuming the Φ e and Φ µ fields to be completely decoupled.

III. MUON/ELECTRON MAGNETIC DIPOLE MOMENTS
The anomalous magnetic dipole moment of lepton is usually denoted by a ≡ (g −2) /2.
Currently, the differences between the experimental value a exp and the SM prediction a SM for = µ, e are given by presenting about 3.3σ [9] and 2.4σ [11] deviations, respectively. In our model, the new contribution to a , denoted by ∆a NP , mainly comes from the one-loop diagrams shown in Fig. 1, with Z 2 -odd particles running in the loop. These contributions are calculated to be where g ,k L,R denote the Yukawa couplings for theχ P L,R η k (χ P L,R η ± k ) vertices in the model with Y D = 1 (0). More explicitly, The loop functions are defined as follows: for the scenario of Y D = 1 at the 1σ (darker color) and 2σ (lighter color) levels.
where at any given x, In both Eqs. (30) and (31), the coefficient of Re(g ,k L g ,k * R ) can be much larger than that of |g ,k L | 2 + |g ,k R | 2 by a factor of M χ /m , and becomes the dominant factor for ∆a NP . We note that for a fixed value of M χ and the Yukawa couplings, a larger magnitude of the dominant term is obtained for a smaller mass of the scalar boson running in the loop. In addition, the contribution to the dominant term from the lighter scalar boson (η 0 1 or η ± 1 ) is opposite in sign to that from the heavier one (η 0 2 or η ± 2 ) due to the orthogonal rotation of the scalar fields, as seen in Eq. (32). Therefore, the sign of ∆a NP is determined by Re(g ,1 L g ,1 * R ). We thus take Re(g µ,1 L g µ,1 * R ) < 0 and Re(g e,1 L g e,1 * R ) > 0 in order to obtain ∆a NP µ > 0 and ∆a NP e < 0, as required by data.
This in turn can be realized by taking f µ L > 0, f µ R < 0, f e L,R > 0, and the mixing angle θ to be in the first quadrant. Note here that with a degenerate mass for η 1 and η 2 , ∆a NP would vanish due to the cancellation between the contributions of the two scalar bosons.
Therefore, a non-zero mass splitting between η 1 and η 2 is required. For simplicity, we take |f L | = |f R |(≡ f ) in the following analyses.
In Fig. 2, we show the regions in the plane of f and the mass M χ that can explain the corresponding (g − 2) anomalies in the scenario with Y D = 1. The left and right panels show the allowed regions for a mass difference ∆m η ≡ m η 2 − m η 1 of 100 GeV and 300 GeV, respectively. In this scenario, the lighter scalar η 0 1 can be the DM candidate and its mass m η 1 is fixed to be 80 GeV. In the next section, we will see that this choice of the DM mass is   Fig. 2, but in the scenario of Y D = 0. The mass of the lighter charged scalar η ± 1 is set to be 200 GeV.
compatible with both the observed relic density and the direct search experiments. It is clear that a smaller value of ∆m η results in a larger cancellation between the ∆a NP contributions from the two scalar bosons, thus pushing the required Yukawa couplings higher for the same M χ . Also, for a fixed M χ , the required value of f e is smaller than f µ by roughly a factor of 4. This can be understood in such a way that from Eq. (30) the ratio ∆a NP µ /∆a NP e is roughly given by m µ /m e × |f µ /f e | 2 200 × |f µ /f e | 2 if we take M χµ = M χe . Therefore, with the required ratio ∆a µ /∆a e by data to be about 3000, the Yukawa coupling for the muon needed to explain the data should indeed be about 4 times larger than that for the electron.
In Fig. 3, we show the results for Y D = 0. In this scenario, the lighter charged scalar boson η ± 1 would not be a DM candidate and its mass m η ± 1 would not be strongly constrained by the relic density and the direct search experiments. However, O(1) TeV of m η ± 1 requires a large Yukawa coupling f µ to explain the muon g − 2 anomaly, which leads to too small a relic density to explain the observed density of DM as we will see in the next section. We thus take m η ± 1 = 200 GeV as an successful example. In Fig. 3, we also observe a similar trend that for a fixed M χ , the required f e is smaller than f µ by roughly a factor of 4 and both are pushed higher for smaller ∆m η ± . Unlike the scenario of Y D = 1, the contours turn around at M χ ∼ 150 GeV in this scenario. This is because the dominant term in Eq. (31) reaches its maximum at M χ = m η ± k , so that the required value of f becomes smallest at M χ ∼ 150 GeV 3 . Note that this turning point is lower in the left plot because of the larger cancellations for the case with ∆m η ± = 100 GeV (left) than that with ∆m η ± = 300 GeV (right).
We note that, in both scenarios with Y D = 1 and 0, the charged Z 2 -odd particles can be pair produced at colliders and their leptonic decays are subject to constraints from the experimental searches at the LHC. These constraints will be discussed in Sec. IV B.
Lastly, we comment on the contributions from two-loop Barr-Zee type diagrams [26].
In our model, new contributions to the Barr-Zee type diagrams can enter via the Z 2 -odd particle loops in the effective hγγ, hZγ and W + W − γ vertices. The first two vertices, in particular, may give rise to sizable contributions to ∆a NP , if the scalar trilinear couplings are taken to be large. However, such large values are highly constrained by the Higgs data to be discussed in Sec. IV B. Together with the smallness of the Yukawa couplings for muon and electron, we find that contributions from these two types of diagrams are negligible.
The contributions from diagrams with the W + W − γ effective vertex have been examined in detail in Ref. [27]. It is shown that the contributions are at least two orders of magnitude smaller than the experimental measurements and can also be safely neglected.

IV. PHENOMENOLOGY
In this section, we discuss the phenomenological consequences of our models, focusing on the physics of DM and collider signatures of the new particles.

A. Dark Matter Phenomenology
As alluded to in Sec. II, the lightest neutral Z 2 -odd particle can be a DM candidate and corresponds to η 0 1 (η 0 H or χ ) in the scenario of Y D = 1 (Y D = 0). Current measurements of the cosmic microwave background radiation by the Planck satellite show the DM relic density to be [28] Ω DM h 2 = 0.120 ± 0.001, assuming the cold DM scenario. 3 For Y D = 1, the dominant term in Eq. (30) reaches its maximum at M χ ∼ 0.12m η1 . Thus, the turning behavior is not observed as we take η 0 1 to be the lightest particle. We first discuss the relic density of DM in the scenario of Y D = 1. The important DM annihilation processes are shown in Fig. 4. The amplitude of the s-channel Higgs-mediated process is proportional to the η 0 1 η 0 1 h coupling calculated as where the λ 3 and λ 7 parameters are chosen as independent parameters [see Eqs. (11) and (17)] in our analyses. Therefore, the λ η 0 1 η 0 1 h coupling can be taken to be any value as far as it satisfies the theoretical bounds discussed in Sec. II. This process can be particularly important when the DM mass is close to half of the Higgs boson mass due to the resonance effect. The amplitude of the t-channel process mediated by the heavier Z 2 -odd scalar bosons becomes important when the DM mass is larger than about 80 GeV because of the threshold of the weak gauge boson channels. The t-channel process mediated by the vector-like lepton χ is sensitive to the Yukawa couplings f L,R , while weakly depending on the mass of the lighter vector-like lepton. In addition to the processes shown in Fig. 4, we also take into account the contributions from DM co-annihilations with the heavier Z 2 -odd particles, i.e., η 0 A , η 0 2 , η ± and χ ± . For numerical calculations, we have implemented our model using FeynRules [29,30] and derived the relic density and direct search constraints using MadDM [31][32][33]. It is worth mentioning that in the Inert Doublet Model (IDM), another solution of the DM mass to satisfy the relic density may exist in a TeV region when the mass splitting among the Z 2 -odd scalar particles is small, typically less than 10 GeV [34]. In such a scenario, DM dominantly annihilates into a pair of weak gauge bosons whose annihilation cross section decreases by O(1/m 2 DM ), while the annihilation into the Higgs bosons is highly suppressed due to small Higgs-DM couplings. In our model, such a high mass solution cannot be realized, because the additional η 0 2 state cannot have the mass close to η 0 1 in order to explain the g − 2 anomaly as discussed in Sec. III. As a result, the (co)annihilation into a pair of the Higgs bosons is not suppressed at the high mass region. This situation can be clearly seen in the right panel of Fig. 6 in which we take ∆m = 30, 60, 120 GeV that can explain the g − 2 anomalies. Indeed, the predicted density is well below the observed value at the high mass region. In fact, we confirm that solutions do not appear even at a few hundred TeV of m η 1 .
In addition to the DM annihilation, the λ hη 0 1 η 0 1 coupling contributes to the scattering of DM with nuclei via the mediation of the Higgs boson, allowing our DM candidate to be probed by the direct search experiments. Fig. 7 shows the spin-independent DMnucleon scattering cross section and its upper limit at 90% confidence level obtained from the XENON1T experiment with a 1-tonne times one year exposure [35]. We find that λ hη 0 1 η 0 1 /v has to be smaller than 0.0026, 0.0034, and 0.0047 for the DM η 0 1 to has a mass around 50, 65 and 80 GeV, respectively, by which we can explain the observed relic density.
In conclusion, the mass of η 0 1 should be about 50, 65 or 80 GeV while having f 0.34 and λ hη 0 in order to satisfy both the relic density and the direct search experiment in the scenario with Y D = 1.
Next, we discuss the scenario with Y D = 0 assuming η 0 H to be the DM candidate. In this scenario, the properties of DM are quite similar to those of the scenario with Y D = 1 discussed above, where the annihilation processes can be obtained by replacing (η 0 1 ,η ± ,e/µ) with (η 0 H ,η ± 1,2 ,ν e /ν µ ) in Fig. 4. The η 0 H η 0 H h coupling is given as Again, this coupling can be taken to be any value due to the independent parameter processes χ χ → ν ν / + − mediated by a neutral or charged Z 2 -odd scalar boson. These processes alone, however, produce a cross section that is too small to account for the observed relic density. Thus, the scenario of having a fermionic DM in our model is ruled out.

B. Collider Phenomenology
We first discuss the constraints from direct searches for new particles at high-energy collider experiments. In our model, all the new particles are Z 2 -odd, and thus would only be produced in pairs at colliders. In addition, due to the new Yukawa interactions for the muon and the electron, their decays typically include a muon or an electron in association with missing energy carried away by the DM. Therefore, our model can be tested by looking for an excess of events with multiple charged leptons plus missing energy, which is identical to the signatures of slepton or chargino production in supersymmetric models.
We first focus on the pair production of the vector-like leptons χ ± at the LHC in the model with Y D = 1. The pair production occurs via the Drell-Yan process mediated by the photon and Z boson, so that its cross section is simply determined by the mass of χ . The Branching Ratio left panel of Fig. 8 shows the cross section of pp → γ * /Z * → χ + χ − with the collision energy of 13 TeV. The cross section is calculated at the leading order using MadGraph_aMC@NLO [36] with the parton distribution functions NNPDF23_lo_as_0130_qed [37]. It is seen that the cross section is about 900, 20 and 0.8 fb for M χ = 150, 300 and 600 GeV, respectively. On the other hand, the decays of χ ± strongly depend on the mass spectrum of the Z 2 -odd scalar bosons. For the case with (m η 0 1 , m η A , m η ± , m η 0 2 ) = (80, 200, 200, 380) GeV, the various decay branching ratios of χ ± are depicted in the right panel of Fig. 8. In this plot, we take θ = π/4 in which the branching ratios do not depend on f . We see that χ ± decay 100% into η 0 1 ± when M χ < 200 GeV because this is the only kinematically allowed channel. At higher masses, χ ± can also decay into η 0 2 ± , η 0 A ± and η ± ν . The heavier Z 2 -odd scalar bosons can further decay into the DM and a SM particle, i.e., η 0 2 → hη 0 1 , η 0 A → Zη 0 1 , and η ± → W ± η 0 1 . Therefore, when these channels are allowed, the final state of the χ ± decays can have 1 or 3 charged leptons. We note that the tri-lepton channel is highly suppressed by the small branching ratio of the leptonic decays of the Z boson or the Higgs boson.
In Fig. 9, we show the observed exclusion limit on the vector-like lepton masses M χ using the same set of parameters as in Fig. 8. The observed limit is derived based on the searches for events with exactly two or three electrons or muons and missing transverse momentum performed by the ATLAS experiment using the 36.1 fb −1 dataset of √ s = 13 TeV collisions [38]. We use MadGraph_aMC@NLO [36] to simulate the events and to compute the χ + χ − production cross section at the leading order. The events are further processed by Checkmate [39][40][41][42], which utilizes Pythia8 [43,44] for parton showering and hadronization and Delphes3 [45] for detector simulations and compares the number of events with the limit in a given signal region provided by the ATLAS experiment [46]. With our parameter choice, M χ 270 GeV is excluded. Note also that such lower bounds on the χ mass depend on the mass spectrum of the Z 2 -odd scalar bosons, and are usually lower than the bounds extracted in the literature (e.g., in Ref. [21]).  In Fig. 10, we summarize all the constraints discussed above in our model with Y D = 1.
The regions shaded by dark green and orange can explain, respectively, the electron and muon g − 2 within 1σ. The lower bound on M χ is derived from the observed direct search limit by the ATLAS collaboration (see Fig. 9), while the region shaded by brown cannot explain the DM relic density as the annihilation cross section of DM in this region is too large to reach the observed density (see Fig. 6).
We note that in addition to the pair production of χ ± , the inert scalar bosons can also be produced in pairs. When we consider the case where the vector-like lepton masses are larger than the masses of the inert scalar bosons, the signature of these scalar bosons become quite similar to that given in the IDM. As shown in Ref. [47], the upper limit on the cross section of multi-lepton final states given by the LHC Run-II data is typically one or more than one order of magnitude larger than that predicted in the IDM. Thus, we can safely avoid the bound from the direct searches for the inert scalar bosons at the LHC.
Let us briefly comment on the collider signatures in the model with Y D = 0. In this scenario, the vector-like lepton is electrically neutral, so that it is not produced in pair via the Drell-Yan process, but can be produced from decays of the inert scalar bosons, e.g., η ± 1,2 → ± χ 0 and η 0 H,A → ν χ 0 . The most promising process to test this scenario could then be a pair production of the charged inert scalar bosons pp → η ± i η ∓ j (i, j = 1, 2). However, we find that the production cross sections of η ± 1,2 are roughly one order of magnitude smaller than those of vector-like leptons shown in Fig. 8, so that such process is more weakly constrained by the current LHC data as compared with that in the model with Y D = 1.
Finally, we discuss an indirect test of our model by focusing on modifications in the Higgs boson couplings. Because of the Z 2 symmetry, the Higgs boson couplings do not change from their SM values at tree level. However, the loop-induced hγγ and hZγ couplings can be modified due to the new charged scalar boson loops, i.e., η ± (η ± 1 and η ± 2 ) in the model with Y D = 1 (Y D = 0). In order to discuss the modifications to the h → γγ and h → Zγ decays, we introduce the signal strength µ γγ and µ Zγ defined as follows: In our model, the production cross section of the Higgs boson should be the same as in the SM. Consequently, these signal strengths are simply given by the ratio of the branching ratio between our model and the SM. The decay rates of h → γγ and h → Zγ depend on the Higgs boson couplings to the charged scalar bosons, which are calculated as and for Y D = 0, The current global average of the Higgs diphoton signal strength is given by µ Exp γγ = 1.10 +0. 10 −0.09 [9], where the deviation of the central value from the SM expectation mainly originates from the CMS measurements [48]. On the other hand, the h → Zγ decay has not yet been observed, and the strongest limit is given by the ATLAS experiment [49], where the observed upper limit for the signal strength µ Zγ is 6.6 at 95% confidence level.
In Fig. 11 Fig. 11, it is clear that both scenarios of our model are able to accommodate the current experimental constraints from the h → γγ decay within a reasonably large range of parameter space without violating the perturbative unitarity and vacuum stability constraints.
As the decay rates of h → γγ and h → Zγ have different dependences on couplings, to see the correlation between µ γγ and µ Zγ would be useful in order to extract the structure of the model [50]. In Fig. 12, we show the correlation between µ Zγ and µ γγ for the scenario of Y D = 1 (left) and Y D = 0 (right). We only show the points which are allowed by the perturbative unitarity and vacuum stability bounds. For Y D = 1, we see that µ Zγ is strongly correlated with µ γγ . Within the 2σ region around the current measurements of µ Exp γγ , a signal strength for h → Zγ is predicted to be from 0.97 to 1.05. Such a prediction can be slightly modified by the choice of the mixing angle θ and the masses of the Z 2 -odd scalar bosons.
For Y D = 0, we observe no or little correlation between µ Zγ and µ γγ . This is because the contributions from the pure η ± 1 and η ± 2 loops are small in our particular choice of θ = π/4 due to smaller η + 1 η − 1 Z and η + 2 η − 2 Z couplings. On the other hand, the η ± 1 and η ± 2 mixed loop contribution, which appears in the h → Zγ decay but not the h → γγ decay, can be sizable.
The coupling λ hη ± 1 η ∓ 2 that contributes to this new diagram is given by With this additional mixed loop contribution, the model with Y D = 0 can predict µ Zγ = 1 even when µ γγ = 1. We note that our prediction on µ Zγ is sensitive to the choice of θ, because of the Zη ± i η ∓ j couplings. By scanning the mixing angle θ while imposing both theoretical and experimental constraints, we find that the model with Y D = 0 would predict an h → Zγ signal strength that is at most +10% larger than the SM value.

V. CONCLUSIONS
To explain the muon and electron g − 2 anomalies and the dark matter data, we have proposed a new model whose symmetry is enlarged to have a global U (1) and a discrete Z 2 symmetries and whose particle content is extended with two vector-like leptons and the inert scalar singlet and doublet fields. Depending upon the hypercharge assignment of the new fields, there are two different scenarios. Thanks to the new symmetries, we can safely avoid the lepton flavor-violating decays of charged leptons, while obtaining new contributions to the muon and electron g −2 with the desired signs and magnitudes for the data. In addition, the symmetries guarantee the stability of the DM candidate, which is the lightest neutral Z 2 -odd particle.
We have found that there are regions in the parameter space that can simultaneously h → Zγ signal strength that is at most +10% larger than the SM value.