Current status and muon $g-2$ explanation of lepton portal dark matter

In this paper, we summarize phenomenology in lepton portal dark matter (DM) models, where DM couples to leptons and extra leptons/sleptons. There are several possible setups: complex/real scalar DM and Dirac/Majorana fermion DM. In addition, there are choices for the lepton chirality that couples to DM. We discuss the prediction of each model and compare it with the latest experimental constraints from the DM, the LHC, and the flavor experiments. We also propose a simple setup to achieve the discrepancy in the anomalous magnetic moment of muon.


Introduction
Dark matter (DM) is one of the biggest mysteries in our universe. While evidences of the existence have been built up since the first suggestion in 30's, these are obtained via gravitational force and the fundamental nature of DM, e.g. spin and mass, is known very little. This has motivated theorists to propose a lot of possibilities for particle DM, and encouraged many attempts to reveal the physical characters.
As a good DM candidate, a massive elementary particle that weakly interacts with particles in the Standard Model (SM) has been discussed. If the size of the interaction is compatible with the weak coupling strength and DM mass ranges from a few GeV to a few TeV, the thermal relic density of DM is surprisingly in agreement with the observed value. This fascinating candidate is called weakly interacting massive particle (WIMP) DM, and significant attention has been paid to this kind of DM. In fact, there are many possible extensions of the SM predicting WIMP DM candidates, e.g. wino, bino and higgsino in supersymmetric models [1,2], Kaluza-Klein particles in extra dimension models [3,4], inert scalars in extended Higgs models [5][6][7][8][9], and so on. In these extended models, DM is often accompanied with extra charged or colored particles, which mediates interaction between DM and SM. Such extra particles provide good opportunities to test DM models at the LHC experiment.
Many of those extended models are originally motivated by theoretical problems of the SM, such as the gauge hierarchy problem, rather than by explaining DM itself. However, compatibility of the models with theoretical problems and experimental results considerably restrict the parameter space where the DM problem is resolved. In contrast, once going away from the original motivation and working only on the DM problem, one can take a more effective way using simplified models, in which an electromagnetically (EM) neutral stable particle is simply added as a DM candidate, and interaction of DM to the SM sector is introduced in an ad hoc manner. An additional symmetry may also be imposed in another way to guarantee the DM stability. Given a specific setup in the framework of simplified models, one can study DM physics, such as thermal relic density and DM signatures, with a limited number of parameters, and examine the viability of the setup. The studies based on simplified models can cover many theoretically motivated models.
To examine simplified models, it is helpful to classify DM candidates according to spin and interaction with SM particles #1 . As an example, let us assume that DM is not charged under the electroweak (EW) symmetry and a (complex) scalar field. The simplest possible interaction between this DM candidate and the SM sector is via a scalar quartic coupling, λ X |X| 2 |H| 2 , where X is the DM field and H the SM doublet Higgs field. The quartic coupling is responsible for DM thermal production and, once it is fixed by the observed value, one can unambiguously predict signatures of this DM candidate at high-energy collider, direct and indirect detection experiments. This setup is known as a Higgs portal DM model and has been studied well [10][11][12][13][14][15]. In another example, the scalar DM can be coupled directly to SM quarks/leptons by newly introducing vectorlike #1 In some cases, it is also important whether DM is self-conjugate or not, as discussed later. charged (or colored) fermions. Yukawa couplings λ i f XF f i involving the DM X, the vectorlike fermions F and quarks/leptons f i (i: flavor) are crucial to DM physics. This model makes different DM signatures from that of the Higgs portal model. The vectorlike fermions can be produced at high energy collider, and can be tested directly. Besides, the non-trivial flavor structure of the Yukawa couplings λ i f can induce flavor changing processes, which may bring other constraints into play. These models are called fermion portal models [16,17] and what we study in this paper.
In the latter example, DM annihilation is induced by t-channel exchange of new particle. We can consider similar setups: DM can be either of scalar or fermion field; it may be self-conjugate or not; SM fermions coupled to DM are either of quarks or leptons, which are further SU (2) L doublets or singlets. A certain new fermion or sfermion (scalar) field is introduced, so that the portal Yukawa coupling is allowed and the t-channel annihilations are turned on. There are many works on fermion portal models, but most of them focused on only limited setups or particular phenomena. More importantly, all these studies were done several years ago. In the present paper, therefore, we examine all possible setups in fermion portal models by taking into account the latest data from LHC searches, flavor physics and DM physics, including higher-order corrections to DM searches that have not been studied well. Then, we update the current status of each setup. In particular, we concentrate on namely lepton portal DM models [18,19], given no positive results of direct searches at the LHC and direct detection experiments. For quark portal models, see e.g. Refs. [16,17,[20][21][22][23][24][25][26][27][28][29] and references therein.
In the present paper, we also study impact of lepton portal DM models on precision observables. There is an explicit correlation between these observables and DM thermal relic density, since both of them are induced via DM-lepton Yukawa couplings. Once DM density is fixed, some observables are predicted. Of these, we especially discuss the anomalous magnetic moment (g − 2) for muon. There is a longstanding discrepancy in the muon g − 2 [30]: where a exp µ and a SM µ are the experimental result and the SM prediction of a µ = (g − 2)/2 for muon. #2 It is well-known that the size of the discrepancy is the same order as that of the EW corrections. On the other hand, in the thermal freeze-out scenario, an EW scale mass and coupling are required to produce the correct DM abundance. This coincidence suggests that the DM problem and the discrepancy in ∆a µ can be resolved in this class of models. This issue has been addressed in the literature working on lepton portal DM models [33][34][35][36]. It was observed, however, that due to strong chirality suppressions, ∆a µ cannot be accommodated in the minimal models, where either an SU (2) L singlet or doublet leptonic mediator is introduced. We point out in this paper that if we modify the minimal models and introduce both mediators, the suppression is removed, and the #2 Recently, QED NLO corrections to a pion form factor have been calculated in Ref. [31], where a possibility of the corrections accommodating the discrepancy is examined. However, the authors have shown that the corrections are too small to diminish the existing discrepancy in the SM. For BSM interpretation of the discrepancy, see e.g. Ref. [32].  discrepancy in ∆a µ can be resolved. This possibility is studied in the scalar DM model in Refs. [35,37]. In this paper, we examine both scalar and fermion DM models, and improve the analysis of the DM direct detection taking into account the loop correction.
To our knowledge, we identify for the first time which setup and which parameter space can achieve the discrepancy in the fermion DM models. This paper is organized as follows. In Sec. 2, we classify lepton portal DM models, based on the spin of DM, gauge invariance and renormalizability. In our study, DM is assumed to be either of scalar or fermion. The EW and color neutrality of DM is also required there. Then, we study the phenomenology in each model. In Sec. 3, the LHC constraint on each model is reviewed. In Sec. 4, we study the DM physics, i.e. relic density, direct detection and indirect detection, within the minimal models, and summarize the updated status of each model. Contributions to ∆a µ in the minimal models are also calculated there, and we reconfirm difficulty in inducing a large ∆a µ . In Sec. 5, we examine extended models that is specialized in explaining the discrepancy in ∆a µ . Sec. 6 is devoted to summary. In Appendix, the detail of our calculation concerned with DM physics and renormalization group (RG) are shown.

Lepton portal DM models
We study DM that is a singlet under the SM gauge group and has spin 0 or 1/2. The scalar DM is denoted by X and the fermion DM is denoted by χ. The DM field is also classified into whether it is self-conjugate or not. For the self-conjugate scalar (fermion) DM, X † = X (χ c = χ). The self-conjugation nature is important in DM physics. We further introduce a leptonic mediator whose spin is 1/2 or 0 for the scalar or fermion DM, respectively. The classification of DM is shown in Table 1. In the following, we shall introduce scalar and fermion DM models separately.

Scalar DM models
We shall set up scalar DM models in which DM is an EM and color neutral real or complex scalar. To couple DM to SM leptons at renormalizable level, extra leptons are introduced, which are (2, −1/2) or (1, −1) under the EW gauge symmetry SU (2) L × U (1) Y . Then, the DM particle has Yukawa couplings with these leptons. The extra leptons should be vectorlike to make the models anomaly-free. The matter content in the scalar DM models is shown in Table 2. In order to stabilize DM, a parity symmetry Z 2 or global U (1) symmetry is imposed on the models, under which DM and the extra leptons non-trivially transform while all the SM fields are trivial. This symmetry distinguishes the extra leptons L L (E R ) from the corresponding SM ones i L (e i R ). Note that the Higgs portal coupling |X| 2 |H| 2 is always allowed in the scalar models, and it could have some impacts on DM annihilation and direct detection. In this study, to highlight phenomenology that lepton portal couplings induce, we assume that the Higgs portal couplings is negligibly small, and do not consider any combinational effects with it.
We introduce relevant part of Lagrangian in the scalar DM models. The mass terms for the vectorlike leptons are given by The Yukawa interactions involving the vectorlike leptons are given by where i L and e i R are the SU (2) L doublet and singlet leptons in the SM. The index i = 1, 2, 3 runs over the SM three generations. Here,H ≡ iσ 2 H † . The first two terms are the portal couplings of DM to the SM leptons. The latter two terms in Eq. (3) generate a mass mixing between the vectorlike charged leptons after the EW symmetry breaking. The mass matrix for them is given by where v H ≡ H 0 is the Higgs vacuum expectation value (VEV). The primed field E L (E R ) is the charged component in the doublet vectorlike lepton L L (L R ). We define the mass eigenstates as where c X , s X (X = L, R) satisfy c 2 X + s 2 X = 1. The left-and right-handed fermions can be combined to Dirac fermions as E a ≡ (E La , E Ra ), where a = 1, 2. Their masses are denoted by m E 1 and m E 2 . In the mass base, the Yukawa couplings involving the DM and vectorlike leptons E 1,2 are given by where e i is a Dirac fermion for a SM charged lepton in the i-th generation. The model has chance to resolve the discrepancy in the muon anomalous magnetic moment. The new physics contribution to ∆a µ in the scalar DM models is given by [38][39][40] When m X ∼ O (100 GeV), the current discrepancy can be explained only if chirality-flip effect proportional to the vectorlike lepton mass m E 1,2 is sizable. This arises only in the model with a non-vanishing mass mixing, i.e. s L = 0 or s R = 0. A discrepancy in the electron anomalous magnetic moment ∆a e is also reported [41] recently, although it is less significant than that of the muon. We expect that ∆a e can be explained in the lepton portal models instead of ∆a µ . However, these cannot be explained simultaneously, since the lepton flavor violating (LFV) decay µ → eγ is induced just as the portal Yukawa coupling to the electron is turned on #3 . We do not pursue this possibility in this paper. The simultaneous explanation for both anomalies are studied in Refs. [45][46][47][48][49][50][51][52][53][54].
Throughout this paper, the "minimal" setups refer to the models in which either doublet or singlet leptonic mediator exists, and the other is decoupled. There is no mass mixing between them. In the model Lagrangian introduced above, the minimal setups are therefore obtained by simply neglecting either L or E. As alluded above, the non-minimal setup that involves both the vectorlike leptons in the model is interesting because it can accommodate the discrepancy in the muon g − 2.

Fermion DM models
The fermion DM can be Majorana or Dirac. In the fermion DM models, mediators are complex scalars and are (2, −1/2) or (1, −1) under the EW symmetry. The SU (2) L doublet (singlet) mediator, named slepton, is denoted by L ( E). The matter content #3 Similar conclusion is obtained in models which ∆a e,µ is explained by the loop corrections involving vectorlike leptons and Z boson [42][43][44].  Table 3. As in the scalar DM models, the minimal setups are obtained by neglecting either of the sleptons. The slepton masses and interaction terms are given by Since the quartic coupling in Eq. (9) is irrelevant to DM physics, we simply neglect these terms. The Yukawa couplings involving the DM are given by For simplicity, we use the common notation, λ i L and λ i R , for the portal Yukawa couplings in both scalar and fermion DM models. As in the scalar DM models, we assume that one of the portal Yukawa couplings is sizable in the minimal setup. Models with the both couplings are discussed in Sec. 5.
If both L and E exist in the model, the trilinear terms induce the mixing between these states. The mass matrix is given by where c θ , s θ satisfy c 2 θ + s 2 θ = 1. Their masses are denoted as m E 1 and m E 2 . The Yukawa couplings involving the DM and sleptons are given by The new physics contribution to ∆a µ in the fermion DM model is given by [38][39][40] The loop functions are defined as As in the scalar DM models, both singlet and doublet sleptons are necessary to explain ∆a µ so that the chirality-flip effect proportional to m χ is substantial.
Hereafter, when we analytically study generic features of the models, we will not assume any specific structure of the portal Yukawa couplings. The couplings are supposed to have arbitrary values and, therefore, entirely flavor violating there. On the other hand, if DM has sizable couplings to two or more generations at the same time, such setups will be excluded in reality due to large LFV processes induced. Therefore, when we numerically analyze the models to calculate physical quantities, derive constraints and show some plots, we will identify a structure of the portal couplings such that DM couples exclusively to one of the three generations. Then, strong constraints from LFV processes are satisfied. In the analysis later, we will mainly focus on the namely muon-philic case, i.e. λ µ L,R λ e L,R , λ τ L,R , motivated by the discrepancy in ∆a µ .

The LHC limits on extra lepton and slepton
We discuss experimental limits from the LHC in this section. In the lepton portal models, the lightest mediator decays to a SM lepton and a DM particle. The pair production of the vectorlike leptons E 1 and sleptons E 1 is respectively. Both processes give two SM leptons and a large missing energy which has been studied in slepton searches at the LHC [55-60]. We will study a case that the (s)lepton dominantly decays to a muon and a DM particle. The LHC limits in the electron-philic case will be similar to the ones in the muon-philic case, while the limits will be much weaker in the tau-philic case.
We estimate the experimental limits on sleptons and vectorlike leptons from the results in Ref. [57]. To illustrate our method, we study a limit on slepton first. The limit is estimated from the ATLAS analysis for a doublet smuon µ L pair production. We define an efficiency SR in each signal region (SR) as where σ ref prod is a cross section of µ L µ L production and σ exp SR is an experimental upper bound on the effective cross section in the signal region. Here, µ decays to a muon and a DM particle exclusively. The efficiency SR is a probability that pair-produced smuons pass the kinematic-cut of signal region.
We estimate the efficiency SR from the experimental limits as a function of mass difference between slepton and DM, ∆m ≡ m E 1 − m χ . On the experimental bound in (m E 1 , m χ ) plane, the efficiency is given by Here, the efficiency is assumed to be independent of slepton mass approximately. The production cross section is, on the other hand, determined by slepton mass. When a mass difference is larger than a limit on slepton mass with a massless DM, denoted by m max E , the efficiency is set at SR (m max E ). By using this efficiency function of ∆m, a point ( (19) in any signal region. In the analysis of Ref.
[57], there are two types of signal regions. One is designed for a mass spectrum with large mass difference and requires no jets in an event. The other is designed for a mass spectrum with small mass difference and requires a jet in an event. We refer to the signal regions without (with) jet for ∆m > 200 (≤ 200) GeV.
For vectorlike lepton, we estimate a limit by analogy with slepton search. We estimate the efficiency in Eq. (18) from the limit on the degenerate slepton scenario, where both left-and right-handed selectron and smuon have a common mass. Our estimations could be improved, using a limit which gives more similar shape as that for the vectorlike lepton. Figure 1 shows production cross sections of a pair production of sleptons (left) and vectorlike leptons (right) with √ s = 13 TeV. For the production cross sections of sleptons, we refer to the result of LHC SUSY Cross Section Working Group [61][62][63][64][65]. The mixed slepton is defined as the lighter slepton when s L = s R = s θ = 1/ √ 2. Production cross sections for vectorlike leptons are calculated by MadGraph5 2 6 5 [66] based on an UFO [67] model file generated with FeynRules 2 3 32 [66,68]. As in the vectorlike lepton case, the mixed extra lepton is defined as the lighter one when s L = s R = s θ = 1/ √ 2. Figure 2 shows the estimated 95% C.L. upper limits on the mediator lepton masses and DM masses. The red, blue and green lines show the limit for the doublet, singlet and maximally mixed sleptons (vectorlike leptons), respectively. The doublet leptons are more constrained than the other cases because of the larger production cross sections. The limit on singlet slepton shows good agreement with the experimental limit in Ref. [57].
The analysis in Ref.
[57] excludes parameter space where the mass difference is larger than about 100 GeV. More degenerate region would be excluded by more dedicated searches exploiting an initial state radiation and soft leptons [56,69]. The limits, however, exist only in restricted parameter space where the mass difference is about 10 GeV and the tightest limit is about 250 GeV for the degenerate four sleptons scenario. The limits from these searches are not shown in the following analyses, but we shall note that a shallow parameter space with ∆m 10 GeV would be constrained by these searches. Let us comment on the heavier mediator leptons. The heavier states E 2 and E 2 can decay to a SM lepton and DM particle as the lightest ones. In addition, these may decay to a lighter mediator lepton and a SM boson. For instance, E 2 → E 1 Z → µZX is possible and can give clean signals with three charged leptons and large missing energy per one vectorlike lepton E 2 . The branching fractions to a SM boson can be comparable with or even dominate over that to a lepton and a DM particle depending on the couplings and mass spectrum. This is an interesting possibility to discover the lepton portal models but this is beyond the scope of this paper. In the following, we only show the limits from a pair production of the lightest mediator leptons in Fig. 2.

Dark matter physics in the minimal models
In this section, we study DM-related physics and discuss constraints on the models with the minimal matter contents, where either SU (2) L doublet or singlet mediator field exists. There are four types of DM (real/complex/Dirac/Majorana), and the mediator field is (i) an SU (2) L doublet L/ L or (ii) an SU (2) L singlet E/ E. The next-to-minimal models with both mediator fields will be discussed in the next section.

Annihilation cross sections
Pair annihilation is a basic property of particle DM. It governs the DM abundance based on the thermal freeze-out mechanism. If annihilation occurs in halo, it also contributes . to cosmic ray flux which can be probed by indirect searches at telescopes. Since DM particles are non-relativistic at freeze-out and also in halo, it is useful to expand the pair annihilation cross section in terms of relative velocity v of DM particles: where a A , b A and c A are dubbed as partial s-wave, p-wave and d-wave contributions, respectively. Here, the subscript A represents a final state of the annihilation process.
In the lepton portal models, the Yukawa couplings in Eqs. (3) and (10) induce various annihilation processes. In this section, we summarize the important features of the DM annihilation to i j , i j V and V V , where i is a SM charged lepton or neutrino and V ( ) is a SM gauge boson. Here, i = 1, 2, 3 is the flavor index. Note that the neutrinos have only left-handed component. The sample diagrams are shown in Fig. 3. The full analytical formulas are shown in Appendix A. Figure 4 shows ratios of thermal averaged cross sections, σv V and σv V V , to that of the tree-level two-body annihilations σv , with m DM = 500 GeV, the SU (2) L complex DM doublet mediator doublet mediators and a muon-philic coupling (λ µ L = 0, λ e L = λ τ L = 0). Here, σv (V ) = σv µμ(V ) + σv νµνµ(V ) and σv ¯ W = σv µνµW + + σv μνµW − . The velocity suppressed processes are evaluated at the freeze-out temperature: v 2 0.24 and v 4 0.1. The plots are independent of the value of the coupling λ µ L , since it cancels out among the numerator and denominator. In the real, complex and Majorana cases, these are also independent of the choice of the lepton flavor as far as we assume only one of the couplings to be non-vanishing. This is not the case for the Dirac DM, however. There is a slight flavor dependence in σv ¯ γ due to a collinear divergence in the limit m → 0. This will be discussed later more concretely. Phenomenologically important effects in these processes are commented in the following.
-DM pair annihilation into i¯ j where m DM denotes the DM mass: m DM is identical to m X in the scalar DM model and m χ in the fermion DM model. r and i are defined as r ≡ m L /m X and i = m 2 i /m 2 X in the scalar DM model, and r ≡ m L /m χ and i = m 2 i /m 2 χ in the fermion DM model, with m i being mass of the lepton i . Here the sub-leading order in i are neglected. In the models except for the Dirac DM model, the s-wave contribution is helicity suppressed ( i 1). In the real DM, the p-wave contribution b i j is also helicity suppressed and the leading contribution is the d-wave, so that the annihilation cross section is very suppressed and other processes discussed below are relatively important. The expressions for the case (ii) is obtained by replacing λ i L → λ i * R , m L → m E and m L → m E . The full analytical expressions of the expansion coefficients are in Appendix A.
-DM pair annihilation into i¯ j V Figure 3 (middle) shows a diagram for the DM annihilation into a pair of leptons, accompanied with a gauge boson V , where V = γ, Z, W . In all types of DM, these processes have s-wave contributions without the helicity suppressions, while these are suppressed by a gauge coupling strength and three-body phase space by a factor ∼ α/π. Using the parametrization of Eq. (20), the relative importance at freeze-out is evaluated by where n means the dominant partial wave of σv ¯ in each model. As seen in Fig. 4, this ratio is O(1-0.1) in the real DM model (n = 4) depending on r, while no more than 0.1 in the other three models. We will therefore include these processes in calculating DM thermal abundance only in the real DM model. The concrete expressions of the cross sections as well as their squared amplitudes are listed in Appendix A.
In general, the three-body processes are superposition of the final state radiation (FSR) from on-shell leptons, and an emission from the off-shell intermediate state, the namely virtual internal bremsstrahlung (VIB). The differential cross section of the FSR is related to the two-body cross section [70], independently of the types of DM. Here, x is defined as x = 2E γ / √ s with the photon energy, E γ . If (σv) is helicity suppressed, the FSR contribution is also suppressed and negligible. The real, complex and Majorana DM models meet this condition. In these models, the three-body processes are dominated by the VIB, and exhibit a sharp spectrum of emitted vector bosons V around E V = m X(χ) if the DM and mediator masses are nearly degenerate [71][72][73][74][75][76]. Figure 5 shows a photon spectrum from the three-body annihilation χχ → i j γ in the Majorana DM model. The spectrum in the scalar DM models is the same as this. An emitted photon in the VIB process fairly looks like a monochromatic spectrum within detector resolution. This sharp spectrum will be a distinctive signal of these types of DM models as will be discussed in Sec. 4.3. In the Dirac model, (σv) ¯ is not helicity suppressed, then the FSR entirely dominates the three-body process. It follows from Eq. (23) that in the massless lepton limit, the FSR cross section is not well-defined, so that we have to keep the lepton mass finite to evaluate it. This leads to the slight flavor dependence of σv ¯ γ / σv ¯ , as mentioned above. Moreover, when we integrate over x, we encounter a divergence at the infrared edge x = 0. This can be eliminated by taking radiative corrections to the two-body processes into account. In Fig. 4 (bottom-left), to evaluate the size of the three-body process, we just regularized the infrared divergence by integrating over the range 0.1 < x < 1. This corresponds to introducing a sharp infrared cutoff at E γ = 0.1m DM . We have confirmed that the numerical result obtained in that way is smaller than the so-called double logarithm approximation, (σv) ¯ (α/π) log(4m 2 DM /m 2 ) log(m 2 DM /Λ 2 IR ), by a factor ∼ 1.5.

-DM pair annihilation into V V
The DM can pair-annihilate into two gauge bosons via loop diagrams shown in the Fig. 3 (right). The possible final states are V V = γγ, γZ, ZZ, W + W − . We do not consider the annihilation into hZ, hh, since these are further suppressed by the small Higgs Yukawa couplings of the charged leptons. These processes cannot be leading contributions at the freeze-out in any case due to the suppression via the gauge couplings and loop factor.
Nevertheless, XX → γγ, γZ will be significant at indirect detection of the real scalar DM. For large r = m E 1 /m X , the cross section is scaled as Hence, the loop annihilation XX → γγ, γZ can be a sizable fraction of the total annihilation cross section for large r. For example, (σv) γγ /(σv) ∼ 0.01 is realized for r = 3. 1.× 10 -26 This small fraction is irrelevant to the DM abundance, but it can be crucial in gammaray searches since produced photons in XX → γγ, γZ are monotonic and can be strongly constrained.

Relic density
In the thermal freeze-out scenario, it is assumed that DM abundance is produced in the thermal bath. The produced number density is determined by the Boltzmann equation, with H the Hubble rate and n eq DM the equilibrium density of DM. Here, σv eff is the effective annihilation cross section of DM, and is expressed in terms of DM pair annihilation and coannihilation. As a concrete example, in the real scalar DM model, it is given by where ∆m = m L − m X . σv = σv + σv V + σv V V represents thermal average of DM pair annihilation cross section mentioned in the previous section, while σ XL(XL) v and σ LL v those of coannihilation cross section whose initial state is XL (XL) and LL, respectively. Since the DM number density is frozen at T f m DM /20, the effective cross section σv eff at that temperature is crucial. It follows from Eq.(26) that coannihilation processes are important when the exponential suppressions are not strong. Naive estimate suggests e is not too small, i.e. ∆m/m X O(0.1) is a necessary condition that coannihilation is effective. Indeed, as we will see in Sec. 4.5, coannihilation processes help to deplete DM thermal abundance down to the observed value with a smaller Yukawa coupling in such a mass region. In our analysis, to take all coannihilation processes into account, we employ micromegas 4.3.5 [77] to numerically solve the Boltzmann equation. In the real scalar DM, we also consider the higher order processes, such as three-body process and loop processes, by including only their s-wave contributions. In the other types of DM, we ignore those contributions, since these are no more than 10 % of the leading process as seen in Fig. 4.

Indirect detection
Indirect dark matter searches look for cosmic ray fluxes, such as gamma ray, anti-proton, positron and neutrino, originated from DM annihilation on top of astrophysical backgrounds. These are good complementary tools of direct detections to probe DM, although these suffer from large systematic uncertainties of astrophysical contributions. In this paper, we will not perform any new data analysis, and will simply rescale the sensitivity curves derived in the literature. Here, we shall discuss the constraints from the indirect searches on the lepton portal models. As discussed below, the meaningful constraints are obtained from the DM annihilation into or γ. In these processes, the constraint on the complex scalar DM is very similar to that on the Majorana DM, because the squared amplitudes are the same in these two types #4 . Thus, we only mention real scalar, Majorana and Dirac fermion DM, to avoid repeating the same comments.
The most prominent target in cosmic ray searches for DM annihilation is a spectral feature, such as gamma line or VIB photon. Such a sharp spectrum can be well disentangled from uncertain astrophysical backgrounds, since attributing it to one astrophysical process is difficult in general. The search for a spectral feature is often a unique way to observe DM signatures from the sky, particularly when the tree-level two-body annihilation is suppressed.
The spectral features in the lepton portal models have been examined in the literature. It has been found that the sensitivity can be improved when we exploit the photon spectrum instead of the continuum photon flux. In the Majorana DM, for example, dedicated searches for the spectral features set orders of magnitude stronger upper bounds on σv ¯ γ than the bounds from the continuum gamma-ray observation of dwarf spheroidal (dSph) galaxies [71]. The study of the spectral feature is then extended to the pair annihilation into γγ, showing the upper limits on a combined annihilation cross section: σv ¯ γ + 2 σv γγ 10 −26 -10 −27 cm 3 /s for DM mass ranging from 40 GeV to 10 TeV [72]. Note that EW gauge invariance requires the existence of weak VIB processes, such as ¯ Z, that exhibit a spectral feature of an anti-proton flux from decays and hadronization of the weak bosons, similarly to that of photon flux. The impacts of the weak VIB emission on the PAMELLA anti-proton search are studied in Majorana DM in [78,79].
#4 See Appendix A for the explicit forms.
Unfortunately, these cosmic-ray signatures of the Majorana DM cannot be observed even at future telescopes unless a boost factor is larger than O (10).
The search for the spectral features can impose good complementary bounds on the real scalar DM. As will be discussed in Sec. 4.4, this type of DM is almost free from the direct detection bound. On the other hand, the fraction of annihilation into ¯ γ and γγ in the total annihilation is sizable because of the strong suppression in the two-body process XX → ¯ . This suggests that more fluxes are predicted than in the Majorana DM, opening up a possibility to discover DM signatures in gamma-ray spectrum. It is shown in Ref. [76] that the future GAMMA-400 [80] and CTA [81] experiments can probe the real scalar DM where m X ≥ 100 GeV and m E 1 /m X ≥ 1.2. The sensitivity prospect will be shown in Fig. 7 later.
For the Dirac DM, the spectral feature will not be applicable, since the gamma-ray flux is dominated by the smooth FSR and a continuum secondary gamma-ray. The latter originates from the decay and fragmentation of the SM particles produced by DM annihilation. Gamma-ray observation of dwarf spheroidal (dSph) galaxies has therefore the best sensitivity. The strongest limit is set on the tau-philic DM. In that case, the lower limit on σv ττ lies below the canonical thermal relic cross section σ 0 ∼ 3×10 −26 [cm 3 /s] for m χ 100 GeV. In the muon-philic or electron-philic DM models with σv µµ,ee = σ 0 , the lower bound on the DM mass is around 10 GeV [82]. These bounds can be simply applied to the Dirac DM, since the unsuppressed two-body process χχ → ¯ is mainly responsible for the thermal freeze-out. In the other DM types, the limit from the observation of dSph galaxies is much weaker than those from the gamma line limit [72].
We briefly comment on other possible constraints. Of particular importance is constraint from positron flux observations. In [83,84], upper limits on annihilation cross section into leptonic final states are derived based on the AMS-02 data of the positron fraction [85]. The most stringent limit is posed on e + e − channel, since the positron spectrum is so sharp even after propagation in galactic space that the spectral search is applicable. The resulting 95% C.L. lower limit on DM mass is 100 GeV, when the thermal relic annihilation cross section into e + e − is assumed [84]. The constraint on µ + µ − and τ + τ − channels (or associated VIB processes) is much weaker, since the positron spectrum is broader and it becomes difficult to disentangle the DM signature from smooth astrophysical background. For further detail of the analysis, see Refs. [83,84]. In the lepton portal models, the positron bound is relevant only to the Dirac DM and gives the limit on m X 100 GeV for e + e − channel. Note that the other three types predict a sharp positron spectrum in eēγ process similarly to the VIB photon, nevertheless the cross section is too small to bring the positron constraints into play.
The lepton portal DM models predict neutrino flux as well. If the leptonic mediator is a weak doublet, the cross section can be sizable since tree-level annihilation into νν is possible. The produced flux can be detected at neutrino telescopes. So far, the observations at neutrino telescopes have found no significant excess of neutrinos over the background. This is interpreted as an upper bound on annihilation cross section. For instance, the ANTARES neutrino telescope has searched for self-annihilation of DM in the center of the Milky Way, and reported bounds on the five representative annihilation channels [86,87]. Of these, τ + τ − , µ + µ − , and νν are relevant to us. Using the latest 11 years data, and assuming the NFW halo profile and 100 % branching ratio, the upper limit on σv νν is 10 −23 -10 −24 [cm 3 /s] for DM mass ranging from 50 GeV to 100 TeV [87]. The limits on τ + τ − and µ + µ − channels are weaker than the former under the same assumption. The searches at IceCube neutrino observatory set similar upper bounds on them [88]. In our models, σv νν is less than 10 −26 [cm 3 /s] in every setup, so that the constraint from the neutrino flux has no impact on the models.

Direct detection
In the lepton portal DM models, there is no tree-level scattering between the DM and a nucleus, so that the leading contribution arises at loop levels. As is well known, primary contribution is photon exchanging in all types of the DM models. The diagrams are shown in Fig. 6. The relevant DM-photon effective interactions are dependent on the DM types and masses. There are also similar Z exchanging contributions in which photons in the leading diagrams are replaced by Z bosons. These are, however, suppressed by lepton masses and therefore sub-leading. In the following, we will summarize typical features of the leading process in each type of DM. See Appendix A for the full analytical expressions for the photon and Z contributions.
In the complex scalar DM model, the photon penguin diagram shown in Fig. 6 (left) is the leading contribution. The induced DM-nucleon effective operator is given by where N = p, n. Here, we define In the limit m L m i , m X , the coefficient C V,N is dominantly given by the photon-penguin diagram, C γ V,N , in the case (i). The expression for the case (ii) is obtained by formally replacing λ i L → λ i * R and m L → m E . When the vectorlike lepton is heavy, the cross section can be considerably enhanced by the logarithmic term, leading to strong limits from direct detection experiments.
In the real scalar DM model, there is no penguin-type contribution. The leading contribution arises at two-loop level via two photon-exchanging [89]. The diagram is shown in Fig. 6 (right). The DM-nucleon scattering cross section is so suppressed that there are no significant constraints from the direct detection. When DM couples to the electron, a DM-electron scattering is induced at tree level. This scattering is limited in light mass region typically m X 1 GeV. In this paper, however, we focus only on m X 100 GeV, since such a light DM scenario would be strictly constrained by, e.g. the LEP experiment in our model. We note that the recent study of limits on DM-electron scattering includes Refs. [90][91][92][93].
In the fermion DM models, there are multi-pole interactions, in addition to contact-type DM-nucleon effective interaction, The differential cross section for elastic DM scattering off a nucleus is expressed in terms of these Wilson coefficients: where we assumed the Dirac DM. Here, µ red = m χ m N /(m χ + m N ) is the reduced mass of DM and a nucleus, Z an atomic number of a nucleus, and m N , J A and µ A nuclear mass, spin and magnetic moment, respectively. In our models, C S,N and C V,N are induced by Z and Higgs penguin diagrams. We use the Helm form factor normalized with F (0) = 1 for F (E R ), and use a spin form factor with thin-shell approximation derived in [94] for F spin (E R ). The first two lines are the spin-independent contributions and the third line is the spin-dependent one.
For the spin-independent one, there are non-contact contributions which appear with 1/E R and lead to infrared enhancement, while these are absent in the spin-dependent part. It should be noted that due to the non-contact contributions, we cannot simply refer to the exclusion curve reported in the experimental papers, in which the contact type interaction is assumed. The different dependencies on E R and v from the contact ones should be taken into account, if the rate is affected by the non-contact contributions. The method to translate the null results at direct detection experiments into limits on the parameter space in our model is explained in Appendix B.
In the Dirac DM model, interactions via the charge radius b χ , the magnetic dipole µ χ , and electric dipole d χ are particularly important, since these are not suppressed by the DM velocity v ∼ 10 −3 and nuclear recoil energy E R = O(10 keV). The anapole a χ is suppressed by v or E R , thus it has a negligible effect on the Dirac case. In the minimal setup, the electric dipole d χ is also vanishing. It is CP-violating, while the photon penguin contribution in Fig. 6 (left) is always proportional to |λ i L | 2 in the case (i), or |λ i R | 2 in the case (ii). Thus, CP-violating contribution is not generated in the minimal setup. If there are both double and singlet mediators, a non-vanishing d χ proportional to Im(λ i L λ i R ) appears. The asymptotic behaviors of b χ and µ χ , as m L → ∞, are given by in the case (i). The expression of the case (ii) is obtained by λ i L → λ i * R and m L → m E . There is a logarithmic enhancement in b χ as in the complex DM, while such an enhancement is absent in µ χ . As pointed out in Ref. [21], the charge radius operator gives dominant contribution to the cross section for m χ 1 TeV, while the magnetic dipole operator becomes dominant one for larger masses. Since the former is a dimension six operator and the latter is a dimension five, their asymptotic behaviors in large m χ are scaled as 1/m 2 χ and 1/m χ , respectively. As a consequence, the magnetic dipole interaction remains more relevant than the charge radius operator, as m χ → ∞.
For the Majorana DM, only the anapole moment a χ is non-vanishing in Eq. (30) due to the Majorana condition, χ c = χ. The differential cross section is obtained by setting b χ = µ χ = d χ = 0 and replacing a χ → 2a χ in Eq. (32). The reason for the latter replacement is that there are twice as many possible contributing diagrams as the Dirac DM. The Z and Higgs penguin contributions are also non-vanishing, but these are suppressed by lepton masses and thus negligible. The direct detection bounds on the Majorana case can be obtained by the same method as in the Dirac DM, although these are essentially very weak due to the suppressions by the DM velocity v and nuclear recoil energy E R .

Current status
In this section, we summarize current limits on the minimal models. Free parameters in the models are the DM mass, the mediator mass and the portal Yukawa couplings λ i L(R) . In this paper, we assume that the Yukawa coupling to one generation of SM leptons is sizable and those to the other generations are negligible, to suppress the LFV processes. In this section, we show plots in the muon-philic case, λ µ L(R) λ e,τ L(R) . The qualitative behavior will be similar to the ones in the electron-and tau-philic cases. For simplicity, the Yukawa couplings to the muon λ µ L,R will be simply denoted by λ L,R in the following. Figures 7 and 8 show the limits from the direct and indirect detection as well as new physics contribution to the muon anomalous magnetic moment ∆a µ and a scale of Landau pole in the four types of DM models. The coupling λ L (or λ R ) is fixed to explain the central value of the observed DM abundance: Ω CDM h 2 = 0.1186 ± 0.0020 [30]. We define a scale of Landau pole Λ where λ L(R) (Λ) = √ 4π. The beta functions for the Yukawa coupling constants are listed in Appendix C. Similar study can be found in Ref. [76] for the real DM, in Ref. [29] for the complex DM, in Refs. [20,33] for the Majorana DM and in Ref. [21] for the Dirac DM.
In these figures, the latest result of the XENON1T experiment [95] excludes the red region. In the gray region on the upper corner, the DM thermal abundance requires a non-perturbatively large coupling λ L(R) > √ 4π below TeV scale. In the green region, coannihilation is too efficient to explain the thermal abundance observed by the Planck collaboration [30]. The orange contours stand for the Landau pole scales and the purple contours stand for ∆a µ . The brown regions are excluded by the slepton searches at the LHC, where the limits are projected from Fig. 2.
The shaded blue regions in the bottom panels are excluded by the current gamma line observations at the Fermi-LAT [96] and the HESS [97]. The cyan lines show the future sensitivity at the GAMAM-400 [80] and the CTA [81]. To obtain these limits, we refer to the 95% C.L. upper limit on the combined cross section σv ¯ γ + 2 σv γγ shown in Fig. 5 of [72]. The referred limit is obtained in the Majorana DM model with r = 1.1, so that strictly speaking, it cannot simply be applied to the real DM case. As shown in Appendix A, however, the photon spectrum of XX → γ in the real DM model is the same as that of χχ → γ in the Majorana DM. In addition, the three-body processes dominate over the two-body annihilation to γγ for r 3 as shown in Fig. 4. Thus, we apply the upper limit on the Majorana DM to the real DM, assuming that the difference between two limits is marginal as far as the models are perturbative. In the following, we discuss more details of the current limits type-by-type.

-Complex scalar DM
Two top panels in Fig. 7 show the results of the complex DM. The two-body annihilation XX † → ¯ is the leading one at freeze-out. O(1) Yukawa coupling is required to explain the DM abundance, since s-wave contribution is helicity suppressed and the p-wave is dominant. The other processes, such as the VIB process XX † → ¯ V , are sub-leading and less than 10% of the total rate, as mentioned above. The large Yukawa coupling lowers the Landau pole scale. For instance, the DM mass should be smaller than 1 TeV so that the Yukawa couplings are perturbative up to the GUT scale around 10 16 GeV. In compressed regions with r ∼ 1, the coannihilation further reduces the thermal abundance and the smaller Yukawa coupling λ L,R is enough to explain the observed value.
The direct detection at the XENON1T [95] has already excluded wide parameter space. In particular, it fully covers the theoretically allowed parameter space in the minimal setup  In the gray region on the upper corner, the DM thermal abundance requires a non-perturbatively large coupling λ L(R) > √ 4π below TeV scale. In the green region, coannihilation is too efficient to explain the thermal abundance observed at the Planck [30]. The orange contours stand for the Landau pole scales and the purple contours stand for ∆a µ . The brown regions are excluded by the slepton searches at the LHC, where the limits are projected from Fig. 2. The shaded blue regions in the bottom panels are excluded by the current gamma line observations at the Fermi-LAT [96] and the HESS [97]. The cyan lines show the future sensitivity at the GAMAM-400 (dotted) [80] and the CTA (dashed) [81].
with the singlet vectorlike lepton. It should be noted that the compressed mass region (m F ≈ m X ) looks already excluded by the XENON experiment in the figure. As pointed out in Ref. [29], allowing O(1%) fine-tuning between m X and m F , one should be able to find a narrow allowed region in this mass regime. In this paper, we do not focus on such a fine-tuned case. For readers who are interested in the fine-tuned case, see e.g. Ref. [29]. In the SU (2) L doublet vectorlike lepton case, we can find the allowed region where 1 TeV m X 3 TeV. It will be covered by the future XENON1T experiment, whose projected limit is shown by a dashed red line, assuming that the sensitivity is 4.5 times better than the current one. The direct detection limit on the tau-philic case will be slightly weaker than that shown in Fig. 7, because of the smaller logarithmic factor ln m τ /m F than in the muon-philic case. The indirect detection is not sensitive to this type of DM as discussed in Sec. 4.3.
The discrepancy of the muon anomalous magnetic moment is hardly explained in this case. As mentioned in Sec. 2.1, the sizable ∆a µ can be obtained only if it is enhanced by the vectorlike lepton mass due to the chirality flip. This is a common feature in the minimal lepton portal models, irrespectively of DM type. This fact motivates us to consider models with both singlet and doublet mediators.

-Real scalar DM
The bottom panels in Fig. 7 show the results of the real DM model. Both the sand pwave contributions in the two-body annihilation XX → ¯ are helicity suppressed and the d-wave is the dominant. The VIB processes XX → V and the loop processes V V have also too small cross sections to be dominant. As a result, the annihilation rates at freezeout temperature are so small that the observed DM abundance can only be explained when m L(E) ≤ 3m X while keeping the perturbative coupling. In wide parameter space, the DM abundance is correctly produced with help of coannihilation.
At the DM direct detection, there is no contribution from the photon penguin diagram, unlike the complex DM. The leading DM-nucleon scattering is induced at 2-loop via diphoton exchange [89]. As a result, the direct detection gives no constraints. On the other hand, the real DM can be probed by using cosmic ray fluxes at the indirect detection experiments. As discussed in Sec. 4.3, the sharp spectral feature of gamma ray is the promising signal. The blue regions where r 1.1 are excluded by the Fermi-LAT [96] and HESS [97]. A combination of future observations at GAMMA-400 [80] and CTA [98] would test the SU (2) L singlet vectorlike lepton case when r 1.1. The sensitivity to the SU (2) L doublet vectorlike lepton case is slightly weaker, because σv γ + 2 σv γγ is smaller than the singlet case.

-Dirac fermion DM
The results of the Dirac fermion DM case are shown on the top panels in Fig. 8. The prime difference from the other types is that the partial s-wave in χχ → ¯ annihilation is not helicity suppressed. The smaller Yukawa coupling is predicted to account for the DM abundance, and thus the Landau pole scale is higher than the other cases.  The constraints from the DM direct detection are much weaker than in the complex DM because of the smaller Yukawa coupling. The current limits on the DM mass from the XENON1T are 300 (220) GeV for the singlet (doublet) slepton case. The XENONnT experiment [99] will probe the region below the red dashed lines and will cover most of the parameter space with perturbative Yukawa coupling. To refer to the XENONnT sensitivity, we assumed it is 50 times better than the current limit.
There are no limits from the indirect searches in the muon-philic case as shown in Fig. 8. The limit from the observations of the dSph galaxies at Fermi-LAT [82] is about 10 GeV in the muon-philic case, while it reaches about 100 GeV in the tau-philic case. In the electron-philic case, the stringent limit about 100 GeV will be set by the positron flux searching for annihilation into e + e − .

-Majorana fermion DM
The pair annihilation of the Majorana DM is similar to that of the complex DM. The pwave in χχ → ¯ dominates the freeze-out processes. The VIB and loop annihilations are subdominant and are no more than 10 % contributions. Hence, sizable Yukawa couplings are required to achieve correct DM density. This causes lower Landau-pole scales in analogy with the complex DM.
The constraint from the DM direct detection is, on the other hand, so weak that there are no exclusion lines in the figures. At 1-loop level, the non-vanishing contributions to DM-nuclei scattering are via the anapole interaction, Z-penguin and Higgs exchanging. The anapole interaction induces the spin-independent scattering. This gives the leading contribution although it is suppressed by DM velocity. The Z-penguin induces contacttype interactions, (χγ µ γ 5 χ)(N γ µ N ) and (χγ µ γ 5 χ)(N γ µ γ 5 N ). The former contributes to the spin-independent scattering and the latter to the spin-dependent one. These interactions, however, give only sub-leading effects due to the suppression by the lepton mass. The former contribution is further suppressed by the DM velocity. The Higgs exchanging contribution is also suppressed by the leptons mass. Altogether, the current sensitivity of direct detection cannot probe the Majorana DM. For a reference, we show the red dotted lines that correspond to the neutrino floor, which is assumed to have 10 times better sensitivity than the XENONnT reach. The regions between the dotted lines are within the neutrino floor sensitivity.
Moreover, the indirect detection hardly constrains the parameter space since the χχ → ¯ process is suppressed by small DM velocity v ∼ 10 −3 in our Galaxy. It might be possible that the photon sharp spectral features from the VIB and loop processes are probed at the gamma ray telescopes if a boost factor is O (10) or larger [72]. We conclude that the Majorana DM is most invisible from the DM searches among the minimal lepton portal models.

Lepton portal models for ∆a µ
We have seen that it is hard to explain the discrepancy of the muon g − 2 in the minimal lepton portal models that have either SU (2) L doublet or singlet mediator. In this section, we study the extended model with both the doublet and singlet mediators, and examine correlation between DM physics and ∆a µ . We focus on the real scalar and Majorana fermion DM, since otherwise most of parameter space will be excluded by the DM direct detection.

Real scalar DM
The presence of the vectorlike lepton mixing influences the DM annihilation. In particular, the s-wave contribution to XX → e iēj appears as This is not suppressed by the light SM lepton masses and dominates the annihilation cross section. Note that the annihilation into a neutrino pair is not changed because of the absence of singlet vectorlike neutrino in our model. Let us discuss a correlation between ∆a µ and the DM abundance. First of all, we assume for simplicity that the DM abundance is determined solely by Eq.(35), i.e. the s-wave contribution of XX → µμ, and then estimate the induced ∆a µ . This assumption implies that two portal couplings are not hierarchical (λ L λ R ), the doublet-singlet mixing is not tiny, and any coannihilation process is not effective. Hence, we may regard r = m E 1 /m X as large in the estimate here. In this case, based on Eq. (35), the total annihilation cross section is approximately evaluated as where σ 0 3 × 10 −26 [cm 3 /s] is the canonical value to explain the observed DM density. Inserting this relation into Eq. (7), we obtain Thus, ∆a µ ∼ 5.0 × 10 −8 is predicted, which is clearly too large. This conclusion is independent of DM and mediator masses, the portal Yukawa couplings and mixing angle, as far as the partial s-wave of XX → µμ is responsible for DM production. Of course, the above estimate is not exact and there is a numerical error in fact, since our treatment is too rough. In Eqs. (36) and (37), we ignored the finite r correction and did not take into account a loop function of ∆a µ etc.. Those corrections are, however, not so large, so that the consequence that the prediction of ∆a µ is too large is quite robust. Therefore, to generate the favorable value of ∆a µ , we have to relax some of the above assumptions and prevent the s-wave contribution Eq.(35) from dominating the DM annihilation. One simple way to do that is to employ large coannihilation. It requires a considerable finetuning between m X and m E , while it allows for small pair-annihilation contributions in return, so that the Yukawa couplings become small, leading to a suppression of ∆a µ . Another way is that we consider a hierarchical coupling or small vectorlike lepton mixings. It is obvious from Eq.(35) that the s-wave contribution is proportional to (λ L λ R ) 2 . If we take λ L = ελ R (ε 1), the contribution is suppressed by a factor of ε 2 . Therefore, as ε decreases, the relative importance of the partial d-wave grows, since it has the unsuppressed contribution ∝ λ 4 R . The relative growth of the d-wave contribution breaks the correlation in Eq. (37), and enables us to explain the ∆a µ discrepancy. This is the case with small vectorlike lepton mixings. The s-wave part is also suppressed for the small mixings, while the d-wave is not. Then, we easily reach the same conclusion through a parallel discussion. In this paper, we will focus on the case with a hierarchical coupling, assuming an unsuppressed vectorlike mixing.
To see quantitative details, we show ∆a µ in Fig. 9 as a function of the DM mass m X (left) and the mass ratio r = m E 1 /m X (right). In the analysis, we assume for simplicity that m L = m E and κ =κ which lead to the maximal mixings s L = s R = 1/ √ 2. The resulting condition s L = s R is also in favor of constraints from the EW precision observables (EWPOs). Furthermore, the mass difference between the vectorlike leptons are set to m E 2 − m E 1 = 2κv H = 100 GeV. The relative size of the Yukawa couplings is fixed to λ L = 0.01λ R , and the absolute size of the couplings is determined via the observed DM abundance. ∆a µ is explained within 1σ (2σ) uncertainties on the (light) purple band. We see that ∆a µ is successfully explained together with the DM density.
The behavior of ∆a µ is understood as follows. With the relation λ L = 0.01λ R , the DM pair annihilation is essentially dominated by the d-wave, which is scaled as (σv) µμ ∝ λ 4 R /(m 2 X r 8 ). Then, as far as coannihilation is irrelevant to the DM production, we find the scaling of ∆a µ , where we assumed that λ R is determined via the DM abundance. It follows from the equation that ∆a µ decrease (increase) as the DM mass m X increase (decrease) in this regime. In Fig. 9 (left), we can observe such a behavior in fact. Once coannihilation operates and becomes superior to the pair annihilation, smaller Yukawa couplings are predicted to explain the DM abundance and thus ∆a µ becomes small. Since coannihilation is significant as r → 1, ∆a µ is smaller as r is closer to unity. We see this effect via the coannihilation in the solid (r = 1.01) and dashed (r = 1.1) lines in Fig. 9 (left). As the DM mass decreases, ∆a µ increases until m X = 700 (200) GeV for r = 1.01 (1.1), and then ∆a µ starts decreasing due to the coannihilation dominance. This behavior cannot be found when r = 2, because the coannihilation is not effective when r 1.2. We see the effects of the coannihilation more explicitly in Fig. 9 (right). With DM mass fixed, a small ∆a µ is induced when r ≈ 1, while ∆a µ increases monotonically as r goes away from unity. Remarkably, as r → ∞, ∆a µ looks approaching the asymptotic value Eq.(37), independently of DM mass. To understand that, we shall take a closer look at the DM pair annihilation. Keeping only the relevant contributions, we find the cross section, in the limit of r 1 and ε 1, Despite of the small ε, the s-wave contribution can overcome the d-wave one, if the following condition is satisfied: Since we fix the vectorlike mass splitting to m E 2 − m E 1 = 100 GeV and v 4 ∼ 0.1 at freezeout, the right-hand side is about 10 −4 when m X = 100 GeV and r = 3. Thus, in the case with λ L = 0.01λ R in Fig. 9 (right), the above condition can be satisfied when r 3.
Once the s-wave contribution dominates the DM production, we can repeat the previous discussion to derive Eq.(37), and then arrive at the asymptotic value ∆a µ ∼ 5.0 × 10 −8 for large r. Figure 10 shows the parameter space, where ∆a µ is explained consistently with the DM physics, on the (m X , r − 1) plane (top) and the (m X , λ µ L /λ µ R ) plane (bottom). The size of the Yukawa couplings is fixed to explain the DM abundance. The meanings of green and gray regions are the same in Figs. 7 and 8. The purple lines show the induced values of ∆a µ in unit of 10 −9 . We also highlight the regions with the (light) purple band, where ∆a µ is explained within 1σ (2σ). Since ∆a µ can easily be reduced by relaxing our assumption for the maximal mixing, we conclude that the discrepancy ∆a µ can be   resolved in the real DM model. Note, however, O (10%) mild tunings may be required to realize the DM abundance.
We comment on the constraint from the EW precision observables (EWPOs). In the scalar DM case, O(1) Yukawa coupling is required to achieve the correct relict density since the annihilation cross section is strongly suppressed. We study the EWPOs concerned with the T parameter and the partial Z boson decay width to two µ. We estimate the bound on the EWPOs based on the study in Refs. [30,35]. We conclude that the most stringent constraint comes from the partial decay width of Z boson to 2 µ. The deviation of the decay width of Z → µµ is much less than O(0.1)% that is below the limit [30], in the case with s L = s R = 1/ √ 2. If s L is not the same as s R , the deviation is enhanced and may be tested by the future experiment [35]. In the same analogy with the scalar DM model, the constraint from the EWPOs in the Majorana DM case is not so strong.

Majorana Fermion DM
As in the real DM model, the annihilation χχ → e i e j has the s-wave contribution, that is not suppressed by the lepton masses. Let us study a correlation between DM and ∆a µ in the Majorana DM model. In this case, the DM pair annihilation χχ → µµ will dominantly contribute to the DM production as far as the DM is enough lighter than the sleptons and the coannihilation is negligible. We neglect these effects for simplicity. In this case, the annihilation cross section is given by with r ≡ m E 1 /m χ . Here, the second term is the p-wave contribution from the chiralityconserving interaction. If the s-wave part is dominant, ∆a µ is estimated as where we used the fact that the cross section is approximately equal to σ 0 . Thus, if the above assumption is valid, ∆a µ ∼ 10 −7 is predicted as in the real DM model. However, the Majorana case is not so simple as the real scalar case. Even if λ L = λ R and the maximal mixing s θ = 1/ √ 2, the p-wave can have comparable contribution with the s-wave, in spite of a mild velocity suppression v 2 ∼ 0.24. The cross section in this case is expressed by Majorana DM, s θ = 1/ 2 , λ L μ = 0.1λ R μ Figure 11: The same plots for the Majorana DM as Fig. 9.
As r increases, the s-wave contribution is decaying more rapidly than the p-wave one. This indicates that the latter can still be leading for large r. For example, when we consider 100 GeV DM and the slepton mass difference m 2 = (100 GeV) 2 , two contributions have the similar size at r 1.1. With such a parameter set, the p-wave will therefore be leading contribution where r 1.1. Note that r = 1.1 is not large, so that coannihilation is operative to some extent in this regime. It may be illuminating to derive an asymptotic behavior in the Majorana DM case as well as the real scalar one. In general, both the sand p-wave contributions are comparably important, but in the limit of r → ∞ or m χ → ∞, we find that the p-wave one is dominant. In this limit, ∆a µ is evaluated by The expression in the Majorana DM case is more complicated than the one in the real DM case. ∆a µ depends on m χ , r, λ L /λ R , the slepton mixing and the slepton mass difference. Figure 11 shows ∆a µ as a function of the DM mass m χ and the ratio r ≡ m E 1 /m χ . As in the real DM case, we consider m 2 L = m 2 E , so that the slepton mixing is maximized, )/2. The ratio of the two portal couplings is λ L = 0.1λ R , and the absolute size is determined to explain the DM density. We see that ∆a µ > O (10 −9 ) can be realized where the DM is lighter than a few TeV.
The dependence on the DM mass and the ratio r is slightly different from the one in the real DM case. Concerning the DM mass dependence, as m χ increases, ∆a µ decreases more rapidly than in the real DM. This is understood by comparing Eqs. (45) and (38). ∆a µ ∝ 1/m 2 χ in the Majorana case, while ∆a µ ∝ 1/m χ in the real case. ∆a µ increases with decreasing m χ , but when coannihilation becomes active, ∆a µ starts decreasing. This behavior is similar to the real DM case. Further, we can see in Fig.11 (right) that the dependence on r is also different. Since ∆a µ is scaling as 1/r 2 for large r, ∆a µ decreases Δa μ = 10 -10 with increasing r in the Majorana case. That is in contrast to the real case where ∆a µ is nearly independent of r in the large r limit as shown in Eq. (37). With decreasing r, ∆a µ is increasing, but when r gets close to unity, coannihilation becomes effective, which makes the portal Yukawa couplings small. Then, Eq. (45) is no longer valid. Altogether, ∆a µ is maximized at r − 1 =0.1-1 as shown in Fig. 11. Figure 12 is plotted in the same manner as Fig. 10. We choose the same parameter set as used in Fig. 11. Incompatibility of mass fine-tuning with coupling hierarchy is maintained in the Majorana case as well, but the restriction on the mass degeneracy is a bit loose. More concretely, when we take λ L /λ R = 0.1 (the top-left of Fig. 12), a parameter space, where (m E 1 − m χ )/m χ = O(1), is still allowed. As a consequence, we can accommodate both DM and (g − 2) µ explanation without terrible fine-tuning in mass and hierarchical coupling, but this time we have to cost a very large Yukawa coupling close to √ 4π, which in turn causes a low Landau pole scale.

Summary
In this paper, we have examined the lepton portal DM models, in which scalar or fermion DM couples to SM leptons via Yukawa interaction involving DM and leptons. New EW charged fields, namely vectorlike leptons or sleptons, are introduced, so that the renormalizable Yukawa interactions are allowed by the EW symmetry. Depending on the spin of DM, the new EW charged fields are different. We have classified our DM models as shown in Table 1 and discussed the phenomenology in each model. The DM relic abundance can be achieved thermally in all DM models. One important issue is how we can evade strong bounds from DM direct detection experiments, without spoiling the thermal production of DM. We have found that the complex DM model has been almost excluded by the latest XENON1T result. In the Dirac DM model, s-wave contribution is dominant in the annihilation, so the bound from the DM direct detection is relatively weak, and several hundreds GeV DM is still allowed. The future sensitivity of the XENONnT/LZ experiments will probe the remaining parameter space in the Dirac model. The real DM and the Majorana DM, on the other hand, could evade the direct detection bound, although small mass difference between the DM and the extra EW charged particle is required to evade too low Landau pole scale. One interesting point is that we can test the real DM model at the indirect detection.
We have also investigated a possibility that the lepton portal DM models can explain the discrepancy in the muon anomalous magnetic moment, together with the DM thermal production. We have found that large enough ∆a µ cannot be induced in the minimal models, where either singlet or doublet vectorlike lepton (or slepton) exists. This consequence is common to the four DM types. However, when we consider the extended model with both singlet and doublet, ∆a µ can easily be accommodated. In Sec. 5, we have demonstrated this possibility in both real and Majorana DM models, and clarified explicit correlations between the annihilation cross section of DM and ∆a µ , as shown in Eqs. (37) and (45). If there is DM behind the discrepancy, there are both SU (2) L doublet and singlet fields in addition to DM, and their couplings with DM are large. As far as detectability is concerned, it may be difficult to prove the explanations in both cases, depending on the parameter space. If DM is lighter than a few TeV and the mass difference is larger than 10 %, the indirect searches for DM probably conclude the real DM scenario. The LHC experiments possibly prove both real and Majorana DM scenarios, if DM mass is several hundreds of GeV. In other parameter region, we may be able to test our scenarios in flavor physics, turning on the other couplings between DM and other leptons.
In our study, we have focused on heavy DM region, m DM ≥ 100 GeV. There would be allowed regions, even if DM is lighter than EW gauge bosons, since the direct detection experiments become insensitive to light mass DM. In such a light DM case, however, it may be difficult to achieve the relic abundance of DM thermally, due in part to lower limits on masses of the EW charged new particles from the collider experiments. The resulting large mass splitting between DM and the new charged particles gives rise to the inactive coannihilation mechanism. Without the coannihilation, the DM production relies on the DM pair annihilation, whose cross section scales as σv ∝ λ 4 m 2 DM /m 4 E . It is easy to see that DM cannot be so light when m E 100 GeV while keeping perturbativity. This suggests that thermal production of DM will impose a lower limit on DM mass. In addition, as DM gets lighter, cosmological and astrophysical searches, such as CMB observations and indirect detection, will grow in importance. The bounds depend on detail of DM models and need a further dedicated analysis, that is beyond the scope of this paper.

A Analytic expressions
In this Appendix, we summarize analytic expressions and intermediate results, which have been omitted in the main text. We focus on the minimal models, and consider only the case that DM couples to the left-handed leptons i , introducing the vectorlike doublet L ( L), whose neutral and charge component are assumed to be degenerate m E = m N (m E = m N ). However, the results can easily be translated into the case of right-handed leptons e i R , by replacing λ i L → λ * i R and taking the weak charges of leptons into account. For abbreviation, we use µ ≡ m 2 L /m 2 X and i ≡ m 2 i /m 2 X . In a part of calculation of 1-loop annihilations into V V , we exploit Package-X 2.0 [100].

A.1 Complex scalar DM
For scalar DM, the relevant interaction is given by where i = e i , ν i L denote charged leptons or neutrinos.
A.1.1 Annihilation For the velocity expansion, we find where we only keep the leading order terms in i and j .
2. XX † → i¯ j V The differential cross section is expressed by where x ≡ 2E V / √ s and y ≡ 2E f / √ s. The cross section is obtained by performing x and y integrals over The calculation has been done in the s-wave limit and we neglect lepton masses.
Squared amplitudes: where we mean i¯ j W = e iνj W + +ē i ν j W − and Differential cross sections: where The results are consistent with [73,74,76].
Cross sections: where Li 2 (z) = − 1 0 dt log(1 − zt)/t is the dilogarithm function. The cross sections for the Z and W emissions are obtained by numerical integrations.
where C 0 is the Passarino-Veltman function defined by The above results are consistent with Ref. [76].
The loop function A S ZZ is more complicated. To the best of our knowledge, the explicit expression is given for the first time, which is described by where In the limit of m Z → 0, the latter function is not contributing B ≈ O(m 2 Z ), while the former function is approximate to which is equivalent to the loop function appearing in the photon contribution Eq. (60), as it should be. For W W , the amplitude is obtained by the replacement

A.1.2 Direct detection
Charge radius operator: The photon penguin diagram induces the so-called DM charge radius operator, which in turn provides the DM coupling to the quark vector current, C V,q = −eQ q b X , through the equation of motion for photon. The penguin contribution involving L and i is given by Loop functions arê with ∆ ≡ µ 2 + ( − 1) 2 − 2µ( + 1).

Z-penguin:
The Z penguin diagram also induces the DM coupling to the quark vector current, The contribution is expressed by with g V,q = (T 3 ) q − Q q s 2 W and g A, = (T 3 ) , and the loop function iŝ It should be noted that, if DM couples to the singlet leptons e i R , the sign of the Z penguin contribution is flipped: with g A,e i = −1/2.

A.2 Real scalar DM
The notation of Yukawa interaction is the same as the complex case, but DM is selfconjugate in this case (X † = X).

XX → i¯ j
Similarly, for the following velocity expansion, the coefficients are 2. XX → i¯ j V and V V Cross sections for these processes are just four times larger than the corresponding cross sections in the complex DM. The spectrum of vector boson in XX → i¯ j V and the dependence on µ and ξ V in both processes are completely same.

A.3 Dirac DM
Relevant interaction is given by All results given here are consistent with Ref. [21].

A.3.2 Direct detection
Effective DM interactions to photon are described by The contributions involving L and i are written as follows: Charge radius operator: Magnetic dipole operator: Anapole operator: a χ (µ, ) = 1 12 Z-penguin: The couplings to vector and axial vector currents, are induced. The Z penguin contribution involving L and i is and Like the real DM case, if DM couples to the singlet leptons e i R , the sign of the contribution to C V,q should be flipped, while the contribution to C A,q is unchanged: Higgs exchanging: withĉ H (r, ) = ∆ + µ A.4 Majorana DM The notation of Yukawa interaction is the same as the Dirac case.
A.4.1 Annihilation 2. χχ → i¯ j V Squared amplitudes: where Differential cross sections: where where A F γZ , A F ZZ and A F W W are and Yukawa coupling. Then, we can translate the null results at direct detection into limits on our models by comparing our calculation with experimental results at scattering event level. To pull the number of exclusion events out of the public limit (the solid black line in Fig. 5 of [95]), we required that contact-type scattering reproduces the exclusion curve. The number of exclusion events extracted in this way is N exc. 5.6, 12, 14 for m DM = 30 GeV, 100 GeV, 1 TeV, respectively. This is the way we have drawn the red lines in Figs. 7 and 8 in the text. We have confirmed that our approach produces consistent results with [21] which studied the same type of DM candidate before #5 . We would like to note that even if we take statistical methods more carefully, the results are not affected so much. Taking the statistical fluctuations into account, the number of DM-nuclei scattering events is expressed by with˜ (E R ) an efficiency function #6 , and ν(E R ) the expected number of photoelectrons (PEs) for a given recoil energy E R , for which, respectively, we can take the green line in Fig. 1 of [95] and the black dotted line in Fig. 13 of [99]. The PMT efficiency can be taken to σ PMT = 0.4 [113,114], for example. With this expression and using S min 1 = 3 PEs and S max 1 = 70 PEs, we would be able to calculate the expected number of scattering events in the parameter space considered. We have checked that, for some reference points, this statistical treatment gives practically the same results as those of our simpler analysis.

C Renormalization Group Equations
#5 The definition of our signal region differs by a factor 2 from the definition in Ref. [21], in which the authors defined their signal region as the lower half of the mean nuclear recoil band (e.g. the red solid line in Fig. 4 of [112]). Thus, our bound is aggressive by the factor. Once we use their definition instead of ours, we can find practically the same exclusion line as what they gave in [21]. #6 This is a different function from Eq. (111).
The RGEs of the Yukawa couplings are given by where Here, Y u and Y d are the up and down quark Yukawa matrices in the SM, respectively. We show the full 1-loop RGEs for completeness. We neglect the Yukawa couplings except λ 2 L and λ 2 R in our numerical evaluation. Note that sizable κ and/orκ will make the Landau pole scale lower.

C.1.2 Real Scalar DM
The interaction is where X is now a real scalar field. Compared with the complex scalar DM case, there are Yukawa coupling renormalization fro λ L , λ R , and the other parts are same as the complex case.

C.2 Fermion DM
The Yukawa interactions are where the SM Yukawa matrix Y e is neglected.

C.2.2 Majorana Fermion
The RGEs are given by [55] ATLAS Collaboration, G. Aad et al., Search for direct stau production in events with two hadronic τ -leptons in √ s = 13 TeV pp collisions with the ATLAS detector, arXiv:1911.06660.
[56] ATLAS Collaboration, G. Aad et al., Searches for electroweak production of supersymmetric particles with compressed mass spectra in √ s = 13 TeV pp collisions with the ATLAS detector, arXiv:1911.12606.
[57] ATLAS Collaboration, G. Aad et al., Search for electroweak production of charginos and sleptons decaying into final states with two leptons and missing transverse momentum in √ s = 13 TeV pp collisions using the ATLAS detector, arXiv:1908.08215.
[58] CMS Collaboration, C. Collaboration, Search for supersymmetry in events with tau leptons and missing transverse momentum in proton-proton collisions at sqrt(s)=13 TeV, .