Dark Matter and $(g-2)_{\mu,e}$ in radiative Dirac neutrino mass models

The origin of neutrino mass is a mystery, so is its nature, namely, whether neutrinos are Dirac or Majorana particles. On top of that, hints of large deviations of the muon and the electron anomalous magnetic moments (AMMs) are strong evidence for physics beyond the Standard Model. In this work, piecing these puzzles together, we propose a class of radiative Dirac neutrino mass models to reconcile $(g-2)_{\mu,e}$ anomalies with neutrino oscillation data. In this framework, a common set of new physics (NP) states run through the loops that generate non-zero neutrino mass and, due to chiral enhancement, provide substantial NP contributions to lepton AMMs. In addition, one of the three models studied in this work offers a Dark Matter candidate automatically stabilized by the residual symmetry, whose phenomenology is non-trivially connected to the other two puzzles mentioned above. Finally, our detailed numerical analysis reveals a successful resolution to these mysteries while being consistent with all colliders and cosmological constraints.

These discrepancies are large in magnitude, and the opposite sign between them is somewhat puzzling and hints towards physics beyond the SM (BSM). For attempts to solve these discrepancies simultaneously in BSM frameworks, see, e.g., Refs. .
The observation of neutrino oscillations [78][79][80][81][82][83][84] was the first conclusive evidence that the SM is incomplete and must be extended. Although the existence of non-zero neutrino masses 2 has been firmly established, the nature of neutrinos, viz. Dirac or Majorana is still unknown. As widely known, observation of neutrinoless double beta decay (see, e.g., Ref. [86]) would settle this issue and establish the Majorana nature of neutrinos; however, all experiments so far have null results. Similarly, despite the discovery that about eighty percent of the Universe's gravitating matter is non-luminous, we are yet to know anything about the nature of DM (see, e.g., [87]).
This work considers neutrinos as Dirac particles, and non-zero neutrino masses originate from quantum corrections. This framework proposes a simultaneous solution to the muon and the electron AMMs that is nontrivially linked to the neutrino mass generation mechanism. New physics (NP) contributions to the lepton g −2 and non-zero neutrino masses arise via one-loop correc- tions mediated by a common set of BSM particles. We present a class of radiative Dirac neutrino mass models that share these same features and study in detail a particular model (dubbed as Model-A) belonging to this class that also addresses the DM puzzle and offers rich collider phenomenology. Remarkably, the stability of the DM is guaranteed by the residual symmetry that emerges after U (1) B−L gauge symmetry is spontaneously broken. In Model-A, all these puzzles are deeply intertwined, which is illustrated in Fig. 1.
The paper is organized as follows: in Sec. II, we discuss the general framework, and in Sec. III, we provide details of the set of models under investigation. Next, we discuss the experimental constraints in Sec. IV and carry out a detailed DM phenomenology for Model-A in Sec. V. Finally, we conclude in Sec. VI.

II. SETUP
Neutrinos being Dirac in nature requires the presence of right-handed partners ν Ri (i = 1 − 3), which automatically allows for tree-level neutrino mass via the term −L y ⊃ Y ν L H * ν R when the SM Higgs, H acquires a vacuum expectation value (VEV); here L is the SM lepton doublet, and is the Levi-Civita tensor. This, however, demands y ν ∼ O(10 −11 ) to be consistent with experimental data, which is seemingly unnatural [88] since the Yukawa couplings of the charged fermions in the SM are typically in the range 10 −6 − 1. On the contrary, it is aesthetically attractive to generate Dirac neutrino mass radiatively that would naturally require the corresponding Yukawa couplings typically in the range 10 −3 − 1. Symmetry arguments can naturally forbid the aforementioned tree-level term to achieve this. In this work, we accomplish this by extending the SM gauge symmetry by U (1) B−L [89,90]; in the literature, various types of symmetries are imposed to realize radiative Dirac mass, see, e.g., Refs. .
Gauge anomaly cancellation conditions 3 then require the right-handed neutrinos to carry charges which are either ν R1,2,3 = {−1, −1, −1} or ν R1,2,3 = {5, −4, −4} [137][138][139]. We choose the latter charge assignment since the former allows the unwanted tree-level term in the Lagrangian. To spontaneously break U (1) B−L , we employ a SM singlet scalar that carries three units of B − L charge: σ ∼ (1, 1, 0, 3) 4 . Then non-zero mass for the neutrinos appears through loop diagrams when both the electroweak (EW) and U (1) B−L symmetries are broken (in our setup, the only two fields that acquire VEVs are H  and σ). These loop diagrams originate from ultra-violate (UV) completion of the following unique dimension five operator: where i = 1 − 3 and j = 2, 3. In the following, we very briefly summarize how to construct UV-complete oneloop models utilizing this d = 5 operator; for details, we refer the reader directly to Ref. [126] (we adopt the nomenclature used therein). One-loop Dirac neutrino mass models can be constructed out of two independent topologies: T1-i and T1ii that are shown in Fig. 2. Depending on the Lorentz structures (i.e, fermion-fermion-scalar or scalar-scalarscalar interaction) associated with these vertices, three different diagrams (T1-i-1, T1-i-2, and T1-i-3) can be drawn for topology T1-i. On the other hand, for T1-ii, there is a unique diagram labeled as T1-ii-1. Moreover, by interchanging the external scalar legs of some of these diagrams, in total, eight minimal models can be fabricated [126].
In this work, we focus 5 on the diagram T1-ii-1 and propose explicit models in light of the muon and the electron g − 2. In particular, we formulate three minimal models (labeled as Model-A, Model-B, Model-C), each of which, in addition to correctly reproducing neutrino oscillation data, addresses both (g − 2) µ,e . Among these three models, DM can also be incorporated within Model-A. Moreover, this model shows profound correlations among the neutrino mass, DM, and (g − 2) µ,e as well as offers rich collider phenomenology. Stunningly, no ad hoc symmetry needs to be imposed by hand to realize this dark matter; instead, its stability is assured by a leftover discrete symmetry resulting from the breaking of U (1) B−L gauge symmetry.
Since topology T1-i in Fig. 2 has a 4-particle vertex, there is a unique choice of attaching external H and σ lines to it to form the d = 5 operator of Eq. (4). Completion of the neutrino mass diagram then requires (i) a vector-like 6 Dirac fermion ψ (of three generations) and 5 Conclusions obtained in our work are very general and applicable to most of the models fabricated from topology T1-i. 6 Due to vector-like nature, these fermions do not alter the anomaly cancellation conditions. (ii) two distinct BSM scalars Φ 1,2 ; see the top diagram in Fig. 3. With only these new states, corrections to muon and electron AMMs are too small to be consistent with experimental findings. However, large NP contributions to lepton g − 2 can naturally arise within this setup via chirally enhanced terms that are proportional to vectorlike fermion mass by introducing (iii) the third scalar Φ 3 ; see the bottom diagram in Fig. 3. In Fig. 3, blue lines correspond to the iso-doublet flow drawn to indicate the most economical way to build this class of models. As shown in Fig. 3, this choice requires a single BSM scalar field to be iso-doublet (and no BSM iso-doublet fermion is required). This way, the least number of new degrees of freedom is introduced in a given model belonging to this class (this is our minimality criterion). Minimal models of this class then contain four BSM fields that propagate inside the loops and contribute to neutrino mass as well as lepton g − 2, and have the following quantum numbers: where, Y and β are the hypercharge and the B − L charge, respectively, carried by the vector-like fermion.
It is important to note that: (a) if Y = −1 then β = 2 is not allowed. In this case, a cubic term of the form V ⊃ H † Φ 1 σ is allowed, which would lead to an induced VEV for Φ 1 resulting in a tree-level Dirac mass for the neutrinos via L y ⊃ L Φ * 1 ν R . Furthermore, (b) if Y = 0 then β = −1 is not allowed. In this scenario, a combination of three terms L y ⊃ L H * ψ R , ψ L σν R , and m ψ ψ L ψ R in the Lagrangian would generate neutrino mass via treelevel Dirac seesaw (once VEVs of H and σ are inserted).

III. MODELS
This section discusses three different versions of models belonging to the class introduced in the previous section. For simplicity, we restrict ourselves to the case of |Y | ≤ 1. We label these models as Model-A (Y = 0), Model-B (Y = −1), and Model-C (Y = +1) for which the full quantum numbers of NP states are specified in Table I. In the following text, we provide all the necessary details of these models.
Yukawa interactions:-In three of these models, the new Yukawa part of the Lagrangian takes the following general form: Here Y 1,2,3 are in general 3 × 3 arbitrary matrices, and we define their entries by, Without loss of generality, we choose to work in a basis where the vector-like fermion mass matrix is diagonal, From Fig. 3, it can be seen that Y 1 and Y 2 are responsible for neutrino mass generation, whereas, Y 1 and Y 3 provide NP contributions to the lepton g − 2 that are chirally enhanced. For sizable Yukawa couplings, lepton flavor violating (LFV) processes provide stringent constraints on the off-diagonal couplings of these matrices and force them to take almost diagonal form. To be consistent with the experimental data of (g − 2) e,µ , entries of Y 1 and Y 3 coupling matrices are required to be substantial; hence to suppress LVF, we adopt diagonal textures for these two matrices. On the other hand, entries of Y 2 are required to be somewhat smaller to incorporate correct neutrino mass scale. Hence, for the rest of the analysis, the Yukawa coupling matrices are chosen to have the following form: Note that, due to different B − L charge assignments of the right-handed neutrinos, the first column of Y 2 is zero. For the simplicity of our work, we treat all couplings to be real. Scalar interactions:-Owing to the B − L charge assignments, the scalar potential of this theory takes a simple form. Instead of writing the entire potential, we only provide the relevant interactions required to generate neutrino mass as well as lepton AMMs, here, is the Levi-Civita tensor. Since the SM Higgs doublet transforms trivially under B − L, it does not mix with the BSM scalars.
As a result of these spontaneous symmetry breakings, BSM neutral (charged) scalars originating from Φ 1,2,3 mix as can be seen from Eq. (13). Then the corresponding mass-squared matrices can be written as, where we have defined the following quantities: Model-A: (no doubly charged scalar) Model-B: (no doubly charged scalar) Model-C: (no neutral scalar) Furthermore, we diagonalize these matrices as, where we use the notation: x = {0, +, ++} and M 1 > M 2 . We denote these two mass eigenstates S x 1,2 by (i) Explicitly, the flavor and the mass eigenstates are related via the following identities: It is important to note that due to the simplified form of the scalar potential, there is no mass splitting between Re[φ 0 ] and Im[φ 0 ] components; this is why, neutral scalar cannot serve as a viable DM candidate in Model-A and -B (Model-C does not contain any neutral scalar within Φ 1,2,3 ). Consequently, the only model that provides a DM candidate is Model-A (see Sec. V for details), which is a Dirac fermion DM (Model-B and Model-C do not contain electrically neutral BSM fermion).
Neutrino mass:-The leading contributions to neutrino masses in this theory appear at the one-loop order, as shown in Fig. 3 (Feynman diagram on the top). In this Feynman diagram, BSM neutral (singly-charged) scalars and fermions run through the loop in Model-A (Model-B and Model-C). It is straightforward to compute the (29) and x = 0, +, + for Model-A, -B, -C, respectively. Since neutrinos are Dirac particles, it is simple to solve for the Yukawa couplings x ij ∈ Y 2 in terms of neutrino observables that are known quantities and in terms of couplings y i (to be determined from lepton g − 2) as follows: here and U is the left-rotation matrix that diagonalizes the neutrino mass matrix (recall that Dirac neutrino mass matrix is not symmetric), i.e., The solution given in Eq. (30) corresponds to normal mass ordering for neutrinos. Analogously, the solution for inverted ordering can be trivially constructed by relabelling the B −L charges of the right-handed neutrinos, which would correspond to the third column being zero in Eq. (12) instead of the first column. Note that, due to the non-universal charge assignments of the right-handed neutrinos, one of them carrying five units of B − L charge remains massless (as well as the lightest SM neutrino). However, within this framework, non-zero m 1 is generated via dimension-7 operator 7 of the form Lepton magnetic dipole moment:-The NP contributions to lepton AMMs in this theory appear at the one-loop order, as shown in Fig. 3 (Feynman diagram on the bottom). These contributions are typically large due Here, for illustration, we have only plotted the dominant chiraly enhanced contribution; however, for numerical analysis, we have consider the full expressions.
to vector-like fermion mass insertion in the loop. The outgoing photon in this Feynman diagram is emitted either by an internal scalar or fermion, or by both scalar and fermion, depending on the model. In Model-A, -B, and -C, scalars (fermions) running in the loop are singly-charged (neutral), neutral (singly-charged), and doubly-charged (singly-charged) states, respectively. It is straightforward to evaluate the contribution arising from BSM states, which yields, where summation over the BSM scalars and fermions must be understood. Furthermore, we have defined, with Q ψ = 0, −1, +1 for Model-A, -B, -C, respectively, and, Finally, the re-defined Yukawa couplings (Y * L , Y R ) appearing in Eq. (33) are given by, Model-B : Model-C : In Fig. 4, we present beyond SM contributions to the muon and the electron anomalous magnetic moments arising in three versions of models under consideration. It is clear from this plot that the required large contributions as observed in the experiments can be naturally provided within this framework without requiring large Yukawa couplings. Furthermore, opposite signs of the muon and the electron AMMs can be incorporated via an appropriate choice of the signs of the associated Yukawa couplings.

IV. EXPERIMENTAL CONSTRAINTS
This section briefly describes the phenomenological implications of the proposed models and the current experimental bounds on the BSM states, along with future collider prospects.
LHC bounds on scalars and fermions:-Model-A contains a DM candidate (see Sec. V for details) via which it can be tested in colliders. Specifically, the singly charged scalars can be efficiently pair-produced at the LHC through the s-channel γ/Z exchange. Once pairproduced, each of them decays into a DM and a SM lepton, i.e., pp → + − + / E T . This process mimics the standard slepton searches carried out by ATLAS as well as CMS collaborations [140][141][142] and non-observation of any such processes lead to a bound of m S ± i ≥ 450 GeV [141]. On the contrary, Model-B/C does not contain a DM candidate. Consequently, collider probes of these models are distinct from Model-A. Following Model-A, we assume that BSM scalars are heavier than BSM fermionic states in Model-B/C. Then pair-produced singly (singly and doubly) charged scalars in Model-B (Model-C) decay into a pair of SM lepton (νν or + − depending on singly or doubly charged scalar) and a pair of BSM singly charged fermions (ψ + ψ − ). In fact, the singly charged fermions also get pair-produced through the s-channel γ/Z exchange, which provides the relevant bounds for these models. Note, however, that for a general charge assignment with an arbitrary value of β, a renormalizable coupling responsible for the decay of these fermions may not be present; hence ψ ± are expected to be long-lived.
To make them decay, we fix the B − L charge such that β = −1 for Model-B/C, therefore a mixing term of the form L ⊃ m ψ L R is allowed. Its contribution to SM lepton masses can be fully neglected if 1, where we define m ≡ v H / √ 2. Through this mixing, the vectorlike leptons will promptly decay if 2 × 10 −7 [143]. For a quasi-stable vector-like lepton, assuming chargino like interactions, ATLAS search [144] provides a bound of m ψ ± ≥ 750 GeV [143]. On the other hand, for the prompt decay scenario, each of the pair-produced vectorlike lepton decays into ψ → h , ψ → Z , and ψ → W ν , for which currently there is no collider bound. However, HL-LHC (14 TeV with integrated luminosity of 3ab −1 ) will probe these processes and should be able to exclude up to about m ψ ± ≥ 460 GeV for the first two generations [145] and m ψ ± ≥ 600 GeV for tau-like ψ [143].
Electroweak precision constraints:-Neutrino mass and BSM contributions to the lepton anomalous magnetic moments heavily depend on the splitting between the neutral (and charged) BSM scalars; hence electroweak precision measurements provide stringent constraints on the model parameters. The so-called Tparameter is the most crucial among these oblique parameters, which we compute following Refs. [146][147][148][149] and allow it to vary within the range given by [150], ∆T = 0.03 ± 0.12.
LEP bounds on vector boson:-Since we study gauged B − L extension of the SM, this theory contains a heavy gauge boson Z , which is significantly constrained from the non-observation of any direct signature at LEP and LHC. This gauge boson couples to quarks as well as leptons, thus Z can be produced and searched for at the LEP via s-channel ee → Z → f f processes. To assure LEP-II bound [151], we impose M Z g > 6.94TeV at 95% C.L., which is derived from effective four-lepton operator [152] and valid for M Z 200 GeV. LHC bounds on vector boson:-Furthermore, at hadron colliders, Z can be efficiently produced in the s-channel due to its couplings to quarks, which would show up as resonances in the invariant mass distribution of the decay products. Searching for a massive resonance at the LHC decaying into di-lepton/di-jet imposes stringent limits on the respective production crosssection. The current data from 13 TeV LHC search for a heavy resonance decaying into two leptons (assuming a 100% branching ratio) via pp → Z → + − provides the tightest constraint of M Z > 4.9 TeV. This bound is somewhat relaxed for other branching ratios, which is depicted in Fig. 8 using the current data (ATLAS [153] and CMS [154]) and future projections (HL-LHC [155] and FCC-hh [156]) in the coupling versus mass plane [157].
Cosmological constraints on vector boson:-Since neutrinos are Dirac particles in our scenario, the existence of right-handed neutrinos ν R is implied. Since these right-handed neutrinos are mass degenerate with left-handed neutrinos ν L , they could contribute to the effective number of relativistic degrees of freedom N ef f in the early Universe. In the case of purely SM interactions, the contributions are completely negligible. However, in our model, ν R have gauge interactions with Z , through which they can be in thermal equilibrium with the SM plasma in the early Universe via s-channel f f ↔ ν R ν R processes that increases N ef f . Cosmological data, however, requires that ν R s decouple from the SM plasma much earlier than the ν L . To be specific, Planck 2018 data [158,159] requires that ν R s must decouple at temperatures greater than T > 0.6 GeV [160].
The best current measurement implies N ef f = 2.99 ± 0.17 [158,159], which is in complete agreement with the SM prediction N SM ef f = 3.045 [161][162][163]. Then a 2σ C.L. corresponds to ∆N ef f < 0.285, which we adopt in our analysis. Future sensitives of CMB-S4 experiments will reach a precision of ∆N ef f ∼ 0.03 [164,165] that can probe NP scale up to about 50 TeV [160,166]. For heavy Z , i.e., M Z > 20 GeV, assuming that ν R decouples before T ∼ 0.6 GeV, one obtains M Z /g > 11.4 TeV, for the case of the standard B − L theory [166] (this limit becomes much stronger in the lower mass range, see, e.g., Ref. [160]). We re-scale this bound using the unconventional charge assignment of the right-handed neutrinos in our theory, and the corresponding Planck 2018 bound is presented in Fig. 9.

V. DARK MATTER
In this theory, the Diracness of neutrinos is protected by the U (1) B−L symmetry. Remarkably, this same symmetry is also responsible for stabilizing the DM candidate. The spontaneous symmetry breaking of U (1) B−L by the VEV of σ leaves a residual discrete symmetry Z D Fields (such that ω D = 1) and stabilizes the DM. In particular, stability is guaranteed if D > 2 as well as if D is even [112]. In such a case, all SM fields transform as even powers of ω, and the lightest particle transforming as an odd power of ω will be automatically stable. For concreteness, for Model-A, we fix β = −1/2 in our analysis, and the corresponding charges of all fields under the residual symmetry are presented in Table II. As can be seen from this Table, all particles in the dark sector (ψ, φ, S 0 , η + ) carry charges that are odd powers of ω. The physical states in the dark sector of the Model-A consists of three Dirac fermions: ψ i , (i = 1 − 3), two singly-charged scalars, η i and two complex neutral scalars, S i with i = 1, 2. The DM candidate is the lightest Dirac fermion, which we choose to be ψ 2 ≡ χ. A typical DM annihilation channel in our model is demonstrated in Fig. 5. Before delving into the thermal freezeout scenario of the DM, we delineate the relevant parameter space and constraints set by the DM direct detection experiments.

DM Direct Detection and Parameter Space
DM Direct Detection:-The fermionic DM, being charged under U (1) B−L , will scatter with the nucleon via the exchange of Z at the tree-level, and in the limit of zero momentum transfer one can deduce the following effective interaction term in the Lagrangian between the DM and the nucleon, where, N is the nucleon. The spin independent scattering cross-section associated with this effective interaction is given as, where, the reduced mass is µ r = m N MDM m N +MDM , and the crosssection is insensitive to the DM mass. Consequently, for a fixed DM mass, the limits on the σ SI set by the currently operating DM direct detection experiments will constrain the region of two dimensional M Z − g parameter plane as shown in Fig. 8.
DM Parameter Space:-We want to correlate the constraints on the electron and muon magnetic dipole moments with the parameter space associated with the DM relic density. Hence, based on the choice of the Yukawa matrices given in Eq. (12), and the scalar masses and mixing angles given in Eq. (20), we choose the parameter space for the DM relic density analysis which is spanned by Besides, the Yukawa couplings Y 2 being connected with the neutrino mass generation are much smaller than Y 1 and Y 3 , and therefore don't play any significant role in the DM relic density analysis.

DM Relic Density
We consider the standard thermal freeze-out to achieve the correct DM relic density, Ωh 2 = 0.12 ± 0.001 (68% C.L.) observed by Planck [159]. The dominant (co)annihilation processes which control the DM relic density for our considered DM parameter space, where Y 1 and Y 3 Yukawa matrices being diagonal, are enumerated in the following.
• χ χ → qq, l + l − , ν l ν l i.e. to the SM quark (q), charged lepton (l) and neutrino (ν l ) pairs via the exchange of Z at the s-channel.
When the mass splitting between the DM candidate and another dark sector particle carrying same dark charge is small, the coannihilation channels also become important. In our case, the dominant coannihilation channels involving the DM candidate χ and other dark sector particles: ψ 1 , η + 1,2 and S 1,2 are, • χ ψ 1 → µ − e + , ν µ ν e , χ η − 1,2 , χ S * 1,2 → X X where X, X are the SM particles, and the conjugated channels.
where i, j = 1, 2, and the relevant conjugated channels.
As the DM relic density involves multiple (co)annihilation channels, which also depend on the multi-dimensional parameter space, we consider two scenarios to capture the dynamics in a simplified manner. First, we consider the mass of the Z quite heavy compared to DM and other BSM particles so that the contribution of the Z mediated channels in the thermal freeze-out become negligible. We call it the Z decoupled scenario, and in this case, the thermal freeze-out of the DM is controlled by the BSM Yukawa sector of the Model-A. In the second scenario, dubbed here as BSM Yukawa + Z gauge boson scenario, we take into account both the simultaneous contributions of the BSM Yukawa sector and the Z gauge bosons to determine the DM relic density. Z decoupled scenario:-Although we consider the Z contribution to be negligible in determining the DM relic density for this scenario as mentioned above, the relevant parameter space is still large enough to disentangle the effects of the BSM Yukawa couplings, Y 1 and Y 3 , the mass-splittings between the DM and the charged scalars η + 1,2 and neutral scalars S 1,2 , and the scalar mixing angles θ + and θ 0 on the DM relic density concretely just by scanning over the parameters randomly. Therefore, we consider a few benchmark points of the parameter space and discuss the impact of the variations of the parameters on the DM relic density.
First we consider a benchmark point where we fixed the following parameters as follows, and vary only the Yukawa couplings, Y 1ee ∈ (0.01, 1) and Y 1µµ ∈ −(0.1, 1) randomly to select the points which simultaneously satisfy the constraints on the electron and muon magnetic moments ∆a e and ∆a µ and Tparameter ∆T . Afterwards, we calculate the relic density of the DM F 2 for each of these selected points using micrOMEGAsv5.2 [167] with the model files generated by FeynRules [168].  fig.) and on the mass splittings between the DM and charged and neutral scalars for a fixed value of Y3 µµ (right fig.). Here we vary Y1 µµ ∈ −(0.1, 1). For simplicity, in this presentation, we only consider the parameter points satisfying the constraints on the muon magnetic moment ∆aµ and the T parameter ∆T simultaneously. Again the black horizontal line represents the correct DM relic density. The respective collider bounds are not shown in this plots, see text for details.
Subsequently, for the parameter set given in Eq. (47), the correct DM relic density is obtained for M DM = 310 GeV which is excluded by the LHC limits, and for 1.4 − 1.7 TeV mass range as seen in Fig. 6. To illustrate it, first, we write down the representative annihilation crosssection at low-velocity approximation that is associated with the DM, χ annihilating into the SM leptons (charged and neutrinos) via the exchange of charged (η 1,2 ) and neutral (S 1,2 ) scalars at t-channel, where the terms having small masses of the final-state leptons are neglected, andỸ 1a andỸ 3a are the Yukawa couplings of electron or muon sector (here denoted by a) redefined by absorbing the appropriate scalar mixing angles associated with charged or neutral scalar (here expressed withS i ). As we choose Y 3µµ = 1 for our parameter set in Eq. (47), when the DM mass is in 100−680 GeV range, the term likeỸ 4 3µµ M 2 DM in the numerator of Eq. (48) dominates the cross-section, and therefore we see the narrow band for that mass range. In contrast, when M DM > 680 GeV, the DM mass is large enough to make the terms likeỸ 4 1µµ M 2 DM in the numerator of Eq. 48 also comparable to the terms with Y 3µµ , and thus we see the band at larger masses in the DM relic density plot as |Y 1µµ | takes value in the range (0.1, 1). Besides, as the mass-splittings are small (150 GeV) when the DM mass is close to O(TeV), the coannihilation channels also contribute significantly to the DM relic density, so the simplified Eq. 48 is not applicable for higher DM mass ranges with small mass-splittings.
In Fig. 7 (left), we can see that if Y 3µµ is set to 0.5 instead of 1, the annihilation cross-section decreases for a wide range of the DM mass, and for this reason, the DM relic density remains overabundant for DM mass close to TeV scale as opposed to the case with Y 3µµ = 1 for which it remains underabundant for similar DM mass range. In addition, when we increase the mass-splittings between the DM and the charged and neutral scalars, for lower mass range, the denominator in Eq. 48 becomes more significant, and the annihilation cross-section decreases, which results in the overabundant DM. Besides, larger mass-splittings suppress the coannihilation processes during the thermal freeze-out resulting in the overabundance of the DM for masses in the TeV range.
Furthermore, to capture the impact of Z on the DM relic density, first, we consider four benchmark points presented in Table. III which simultaneously satisfy the constraints on the electron and muon magnetic moments ∆a e and ∆a µ and T-parameter ∆T in the Z decoupled scenario.
For benchmark point I and II of Table III with M DM = 500 GeV, the thermal freeze-out of the DM is dominated Figure 8. Correlation between the mass of the Z gauge boson M Z and the DM relic density Ωh 2 for four benchmark points given in Table III by the processes, χχ → ν µ ν µ , µ + µ − . On the other hand, for benchmark point III and IV with M DM = 3000 GeV, the thermal freeze-out processes are dominated by the coannihilation, for example, by BSM Yukawa sector + Z gauge boson:-Now we calculate the DM relic density for the each benchmark points in Table III by varying the mass of the Z gauge boson M Z ∈ (0.5, 20) TeV and the gauge coupling g ∈ (0.02, 1).
From Fig. 8, we can see that for a fixed M Z , the increase of the gauge coupling g decreases the DM relic abundance determined in the Z decoupled scenario. Therefore, even if the BSM Yukawa sector sets an overabundant DM relic density, one can achieve its observed value by adjusting the Z mass and coupling g as seen in Fig. 8 (upper and lower right figures). Moreover, we can see that the lower the value of the gauge coupling g , the lower the minimum value of the M Z where the Z 's contribution to the DM relic abundance becomes negli-gible. Besides, we see the dip in the DM relic abundance due to the resonant enhancement of DM annihilation processes involving the exchange of Z at the s-channel when M Z ∼ 2M DM .
In Fig. 9, we present the correlation between M Z and g for our four benchmark points with all relevant constraints set by the collider searches and DM direct detection experiments, the observed DM relic density, and cosmological limit on the extra radiation. By correlating the Figs. 8 (upper left and right) and Fig. 9 (left), for  fig.), the cyan and blue points correspond to (M Z , g ) pairs that satisfy the constraints on DM relic density for benchmark I and benchmark II points in the Z decoupled scenario, respectively, whereas for mDM = 3 TeV (right fig.), the yellow-green and brown points represent those with same attribute for benchmark III and benchmark IV points, respectively. The shaded regions are ruled out by current collider bounds, cosmological constraints, and from dark matter direct detection. Future colliders will probe regions inside green and purple lines (not shaded).
M DM = 500 GeV, we see that the current limits provided by ATLAS and CMS have already ruled out the (M Z , g ) values for which the Z gauge boson give significant contributions to the DM relic abundance. Again the close inspection of Figs. 8 (lower left and right) and Fig. 9 (right) for m DM = 3 TeV indicates that apart from the small parameter space close to the resonance point around M Z ∼ 6 TeV, almost all the (M Z , g ) values are ruled out by LEP, LHC, DM direct detection (e.g. PandaX-4T [169] which has set the most stringent limit on the spin independent DM-nucleon cross-section for 10 GeV -10 TeV DM mass range), and Planck constraint on N eff . Furthermore, rest of the parameter space where Z contributions could play an important role in determining DM relic abundance will be fully probed by the HL-LHC and future collider like FCC-hh.
Since the DM particle is inducing the anomalous magnetic dipole moment via mass-insertion, the chirality flip enhancement is strongly correlated to the mass generation mechanism of the associated lepton and causes related loop contributions to its mass. Hence, if M DM O(3) TeV, typically, fine-turning of a large degree is required to adjust the lepton mass correctly. For details, see, e.g., [170] and references therein. This is why in this work, we explored DM phenomenology in the range 100 GeV < ∼ M DM < ∼ 3 TeV.

VI. CONCLUSION
The Fermilab's Muon g − 2 experiment has recently confirmed the longstanding tension of the muon AMM. Furthermore, the recent precise measurement of the electron AMM at the Berkeley Lab shows deviations from the theoretical prediction. These two anomalies together strongly hint towards physics beyond the Standard Model. Besides, the origin of neutrino mass remains a mystery even after the groundbreaking discovery of neutrino oscillation about twenty-five years ago. Moreover, even though we know DM exists, we do not yet know what it is at a fundamental level.
This work proposes a class of radiative Dirac neutrino mass models where neutrino mass arises at a one-loop level. Furthermore, NP states that participate in neutrino mass generation also run through the loops and significantly contribute to (g − 2) µ,e . These large contributions arise due to chirality enhancements required to simultaneously explain the (g − 2) µ and (g − 2) e data. For completeness, we have studied three benchmark models, one of which (Model-A) offers a Dark Matter candidate whose stability is naturally protected without imposing additional symmetries by hand. For Model-A, we have performed a detailed numerical analysis to investigate the correlations among the common model parameters accommodating neutrino oscillation data, the muon, the electron g-2, and dark matter relic abundance. Parameters that generate non-zero neutrino mass also play a non-trivial role in explaining the muon and the electron g − 2 simultaneously; furthermore, these same parameters take part in dark matter annihilations to the SM particles in reproducing the correct relic abundance. This model is subject to numerous constraints arising from colliders and cosmology, through which it can be probed in the current and upcoming experiments. A detailed study of lepton flavor violation and electric dipole moments of the electron and the muon and their possible links to the puzzles resolved in this work is left for future work.