Shedding light on dark matter with recent muon $(g-2)$ and Higgs exotic decay measurements

Recently, we have witnessed two hints of physics beyond the standard model: A 3.3$\sigma$ local excess ($M_{A_0} = 52$ GeV) in the search for $H_0\to A_0 A_0\to b\bar{b}\mu^{+}\mu^{-}$ and a 4.2$\sigma$ deviation from the SM prediction in the $(g-2)_\mu$ measurement. The first excess was found by the ATLAS Collaboration using 139 fb$^{-1}$ data at $\sqrt{s}=13$ TeV. The second deviation is a combination of the results from the Brookhaven E821 and the recently reported Fermilab E989 experiment. We attempt to explain these deviations in terms of a renormalizable simplified dark matter model. Inspired by the null signal result from dark matter (DM) direct detection, we interpret the possible new particle, $A_0$, as a pseudoscalar mediator connecting DM and the standard model. On the other hand, a new vector-like muon lepton can explain these two excesses at the same time while contributiong to the DM phenomenology.


I. INTRODUCTION
The existence of dark matter (DM) is now well established by old and new astrophysical and cosmological evidence. Conversely, its particle properties remain unclear, in particular, the way to incorporate DM and its interactions into the standard model (SM) of particle physics. In the case of DM-quark interaction, current DM direct detection experiments, such as XENON1T, have not observed any DM-nuclei scattering evidence [1]. Furthermore, the Muon g − 2 Collaboration at Fermilab has recently reported an eye-catching measurement of the anomalous magnetic dipole moment of µ ± with a 3.3σ deviation from the SM prediction achieving a combined experimental average of ∆a µ = (2.51 ± 0.59) × 10 −9 [2]. This confirms a long standing tension between SM and experimental data that was previously reported by the E821 experiment at Brookhaven National Laboratory, inspiring several beyond the SM (BSM) extensions (see Ref. [3] for a review). The new average is consistent with a 4.2σ deviation from the SM prediction strongly motivating new proposals . Furthermore, the ATLAS collaboration has reported a search for the Higgs boson exotic decay channel (1.7σ) is reported at M A 0 = 52 GeV with the branching fraction given by BR(H 0 → A 0 A 0 → bbµ + µ − ) ∼ 3.5×10 −4 . Therefore, the null observation from XENON1T and the two possible BSM observations (∆a µ and H 0 → A 0 A 0 → bbµ + µ − ) can be interpreted as hints that the new spin-0 particle A 0 can serve as a mediator to connect DM with SM, especially for its non-trivial coupling with muons.
One can simply consider the additional spin-0 particle A 0 coupling to SM fermion pairs via the mixing with SM-like Higgs boson only or like the scalar/pseudoscalar in the conventional two Higgs doublet models (Type-I, II, X, and Y 2HDM [54], 2HDM+S [55]) with M A 0 = 52 GeV. However, their predicted signal strength for BR(H 0 → A 0 A 0 → bbµ + µ − ) cannot reach the measured 3.5×10 −4 when including all other experimental constraints [56,57]. We go beyond the models mentioned above and involve a new vector-like muon lepton (VLML) not only to enhance BR(A 0 → µ + µ − ) but also to contribute to ∆a µ [60], thus, explaining both excesses. Motivated by these observations and from theoretical considerations, we propose a renormalizable simplified DM model based on extending the SM with three SM singlet fields: a Dirac DM, a VLML, and a pseudoscalar mediator. One important advantage of adopting a pseudoscalar mediator here is that the elastic cross section of DM-nuclei collision is suppressed by the small recoil energy making easy to escape the current XENON1T stringent limit.
The paper is organized as follows. First, we briefly discuss the model setup in Section II.
Next, we consider relevant constraints for this model used in our likelihood functions in Sec. III. In Sec. IV, we present our numerical analysis and the 2σ allowed regions. Finally, we conclude in Sec. V.

II. RENORMALIZABLE SIMPLIFIED DARK MATTER MODEL
In this section, we show our model configuration. We consider a SM singlet Dirac fermion χ as a DM candidate. A new pseudoscalar mediator A is also introduced to explain the null signal result in DM direct detection and possible excess in Higgs exotic decay. Additionally, we introduce a VLML, ψ ± , that will contribute to the (g − 2) µ excess. Thus, the renormalizable Lagrangian for this simplified DM model can be written as where H is the SM scalar SU(2) L doublet and L µ is the second generation left-handed lepton doublet. Note that the dimension-3 terms with µ A and µ A break the parity [61]. The tadpole for A is removed and A = 0 is assumed in this model.
After electroweak symmetry breaking (EWSB), the pseudoscalar A and the SM Higgs boson h mix with each other via the µ A term. The relation between the mass eigenstates, H 0 and A 0 , and the interaction states, h and A, is where the mixing angle, α, is given by with v ∼ 246 GeV the vacuum expectation value. We assign H 0 as the SM-like Higgs boson with M H 0 = 125 GeV and A 0 with M A 0 = 52 GeV to explain the ATLAS Higgs boson exotic decay excess [53]. In the next section, we will see that the LHC Higgs boson measurements can put a strong upper limit on sin 2α in this model. Similarly, the VLML and muon will mix together after EWSB [60]. According to Z → l + l − precision measurements [62], an upper limit for the mixing between the left-handed muon and VLML is set at where M ψ ∼ v implies κ O(10 −2 ). In contrast, κ doesn't suffer from this constraint and its value can be larger.
The Lagrangian to describe the interactions between the SM sector and DM sector via H 0 and A 0 portal can be written as where s α = sin α, c α = cos α and in the first-order approximation of κ, and δ V = 1(2) for V = Z(W ± ). Note that, excepting muon pair, there are only axial (scalar) couplings for H 0 and A 0 with DM pair (SM fermion pairs). The term with g A is important to enhance BR(A 0 → µ + µ − ). Additionally, the CP-violation effect in the muon Yukawa interactions is a feature of this model. For the first order of κ approximation, the interactions of VLML can be expanded as where Note that the {H 0 , A 0 }µψ terms in Eq. (7) can contribute to ∆a µ via one-loop Feynman diagrams [60].
From Eq. (1) we can read off ten undetermined parameters in this model: We can further fix µ A = 5.0 GeV, λ A = 1.0 and y = 0.01 for this analysis because their changes are neither relevant nor significant for this work. Considering that the LHC Higgs boson measurements constrain sin 2 α < 0.12 at 95% C.L. [63,64], we assign the upper bound sin α < 0.3 in our parameter scan 1 . As required by Eq. (4), κ cannot be large and one has to introduce a larger κ to explain the ∆a µ deviation. Furthermore, a very massive ψ ± could suppress the contribution to ∆a µ when we require perturbativity of κ with κ ≤ √ 4π. Finally, for simplicity, we assume M A 0 /2 < M χ < M ψ such that A 0 → χχ and the annihilation χχ → ψ + ψ − are kinematically forbidden. Taking all these facts into consideration, the scan range for the non-fixed parameters is given by the followings bounds where the star ( * ) indicates that the parameter is scanned logarithmically in base 10.

III. EXPERIMENTAL CONSTRAINTS
Recently, the Muon g − 2 Collaboration at Fermilab reported a new measurement of the anomalous magnetic dipole moment of the muon achieving the new combined result [2] ∆a µ = (2.51 ± 0.59) × 10 −9 consistent with a 4.2σ deviation from the SM. Note that the SM uncertainties have been taken into account in the error bar 2 . The dominant BSM contributions to (g − 2) µ in this model comes from one-loop Feynman diagrams with ψ ± , A 0 , and H 0 as [60] where f (t) = (2t 3 + 3t 2 − 6t 2 lnt − 6t + 1)/(t − 1) 4 . Because of the constraint in Eq. (4) and large CP-violation effect of muon EDM, we can safely assume κ κ .
Furthermore, the ATLAS collaboration has detected a local excess in the Higgs decay with a local 3.3σ significance. By taking a conservative approach, we assume that the central value is located at 3.5 × 10 −4 and the likelihood is a Gaussian distribution with error bar equal to 3.5 × 10 −4 divided by 3.3 at M A 0 = 52 GeV. In the on-shell limit, we can calculate individually and then multiply them We can divide the major constraints used for this study in five categories: (1) the LHC Higgs boson measurements, (2) the LEP and LHC A 0 searches, (3) the DM phenomenology, (4) the ATLAS multi-lepton search and (5) electron, muon electric dipole moments (EDMs).

A. The LHC Higgs boson measurements
For the LHC Higgs boson measurements, one can further classify them as
[57]. The partial decay width of H 0 → A 0 A 0 can be written as where If m χ < M H 0 /2, the Higgs boson can also decay to a pair of DM as Therefore, we require Note that, in this study, we fix of sin α and λ HA only.
The decay width Γ(H 0 → µ + µ − ) can be modified by changing the µ-ψ mixing. We require the measured signal strength of H 0 → µ + µ − , relative to the SM prediction, to be 1.19 +0.44 −0.42 (stat) +0.15 −0.14 (syst) as reported by the CMS collaboration [65]. In this work, we focus on the result from CMS because of the rather larger error bar in the one from the ATLAS collaboration [66].
The VLML ψ ± can also contribute H 0 → γγ and further reduce the SM prediction.
The partial decay width of H 0 → γγ in our model is where where α em is the fine-structure constant and N c = 3(1) for quarks (charged leptons).

B. The LEP and LHC A 0 searches
The pseudoscalar A 0 with a mass of 52 GeV can be explored at the LEP and LHC experiments. The decay modes of A 0 for these searches are A 0 → bb, τ + τ − , µ + µ − . For the search at LEP, the most stringent limit is sin 2 α × BR(A 0 → bb) 3.5 × 10 −2 from the e + e − → ZA 0 production [68]. On the other hand, the multilepton final states from the LHC . We will see in Sec. IV that the allowed range for sin α remains smaller than these A 0 limits from LEP and LHC.

C. The DM phenomenology
The DM phenomenology can be classified in three parts: • DM relic density Since we fixed the mass of A 0 to be M A 0 = 52 GeV and the annihilation channel χχ → ψ + ψ − is kinematically forbidden, one would expect that the Higgs resonance regions χχ → f f can play an important role to enhance the annihilation cross section in the early universe. However, the Higgs exchange is suppressed by sin α when compared with A 0 exchange and the dominant channel contributing to the relic density is χχ → µ ± ψ ∓ via A 0 exchange when it is open. The most favorable size of g χ κ in Eq. (7) can be determined using Planck's relic density measurement Ω χ h 2 = 0.12 ± 0.001 [71]. For example, a value of κ 2.46 with g χ 0.13 can fulfill the Planck constraint for m χ around 200 GeV.
• DM direct detection In Eq. (5), DM interacts with quarks via H 0 /A 0 exchange resulting in a suppressed treelevel amplitude for DM-nucleon elastic scattering due to small momentum transfer.
With a simple estimation based on Refs. [74,75], we find that the loop contribution with the condition g χ sin α 0.05 is still below the neutrino floor. Our upper limit g χ sin α 0.6 could bring the loop contribution above the neutrino floor but is still well below the current XENON1T limit [1]. Note that the complete two-loop DM-nucleon elastic scattering can be found in Ref. [76]. Our model corresponds to the special case with CP phases φ χ = π/2 and φ SM = 0 in Eq. (2.1) of the aforementioned reference.
They have shown that the full two-loop calculation can lead to a smaller cross section than previous approaches [74,75]. As we will see in Sec. IV, our actual result may be below g χ sin α = 0.05. Hence, we simply ignore DM direct detection constraints in this study.
• DM indirect detection The dominant channel of DM annihilation in the present universe within our explored can successively decay to µ ∓ bb. The final state µ + µ − bb can produce a soft photon or electron spectrum. Therefore, one may indirectly detect DM annihilation by using dSphs gamma ray data from Fermi [77] and electron-positron data from AMS02 [85].
In addition, the produced photons and relativistic electrons may ionize and heat the universe gas at the recombination epoch. Thus, one may constrain DM annihilation cross section σv by using Cosmic Microwave Background (CMB) anisotropy data from Planck [71]. In this study, because the cosmic ray propagation is rather uncertain and CMB limit is weaker than the one from Fermi dSphs gamma ray data-in our DM mass range of interest-we will only include the Fermi dSphs data in our analysis.

D. The ATLAS multi-lepton search
The VLML ψ ± with 150 ≤ M ψ ≤ 450 GeV in our model is in a mass range that can be explored at the LHC via single and pair productions [72]. The major single production of ψ ± is via gg → H * 0 /A * 0 → µ ± ψ ∓ process. Since the allowed value of κ is tiny, the single production from pp → W ± * → νψ ± and pp → Z * → µ ± ψ ∓ are highly suppressed regarding Eq. (8). Alternatively, the major pair production of ψ ± is via pp → γ * , Z * → ψ + ψ − process.
We have checked that the pair production from gg → H * 0 /A * 0 → ψ + ψ − is much smaller than the Drell-Yan type process with our parameter settings. The production cross sections of Drell-Yan type process are only dependent on M ψ , but the single production process is a function of M ψ , sin α, κ and κ . Therefore, we fix sin α = 0.1, κ = 5 × 10 −2 and κ = 2.0 but vary M ψ from 150-450 GeV to show both single and pair productions at √ s = 13 TeV in the left panel of Fig. 1. Since the cross sections of pair production are much larger than single production ones, we will focus on the constraint of ATLAS multi-lepton search [73] on cross section (pb) Excluded region from recasting of SR5L Excluded region from recasting of SR5L plane.
the pair production channel in this analysis. We remark that ψ ± → µ ± A 0 is the dominant decay mode (BR(ψ ± → µ ± A 0 ) 99%) for our parameter settings.
We closely follow the simulation and analysis in ATLAS multi-lepton search [73] based on a data sample with L = 139 fb −1 at √ s = 13 TeV. The package MadAnalysis 5 [78,79] is used to recast the signal region SR5L for the following process First, the signal processes in Eq. (20) with up to two extra partons were generated from the leading order matrix elements by using Madgraph5 aMC@NLO [92]. We apply Mangano's prescription (MLM) [80,81] for jet-parton matching with M ψ /4 as the matching scale 3 . The Pythia8 [82] is used for parton showering and merging as well as hadronization.
We have modified the ATLAS template in Delphes3 [83] to be consistent with the setup of Ref.
[73] for fast detector simulation. The event selections in this analysis are summarized below: • First, in order to suppress background muons from semileptonic decays of c-and bhadrons, any muon within ∆R = 0.4 of a jet is removed.
• One of the following triggers with efficiency ≥ 90% is used: (1) P T (µ 1 ) > 27 GeV, (2) 3 The K-factor (QCD corrections) for Drell-Yan type process ranges between 1.2 and 1.5. Since its effect is mild, we will not include this K-factor in our recasting. P T (µ 1,2 ) > 15 GeV and (3) P T (µ 1 ) > 23 GeV, P T (µ 2 ) > 9 GeV where the subscript of muon is in the order of P T .
• Finally, if two muons are found in ∆R = 0.6 and P T (µ) < 30 GeV for one of them, both muons are discarded. This selection is used to suppress leptons from a decay chain with multiple heavy flavour quarks backgrounds.
The model-independent limit is reported as σ 95 obs = 0.129 fb as calculated at 95% C.L. from the signal region SR5L 4 . We apply this limit to constrain BR(ψ ± → µ ± A 0 ) × BR(A 0 → µ + µ − ) for varying M ψ as shown in the right panel of Fig. 1. It is expected to find that BR(A 0 → µ + µ − ) already suffers from strong constraints such that A 0 cannot be muonphilic.

E. The EDM of electron and muon
On the EDM side, the latest electron EDM measurement from the ACME collaboration is reported as |d E e | < 1.1 × 10 −29 ecm at 90% C.L. [86]. In the case of the muon EDM, the measurement from the Muon g − 2 Collaboration is |d E µ | < 1.9 × 10 −19 ecm at 95% C.L. [87]. The detailed formulas for both electron and muon EDMs in this model are given in Appendix A. We adopt both electron and muon EDMs constraints in our analysis. Note that the constraints from tau and neutron EDMs from the Belle collaboration [88] and Ref. [89], respectively, are much weaker than electron and muon EDM measurements in this model and can be safely ignored.

IV. RESULTS
In this section we will present our numerical results. In Eq. (10), we show our prior ranges for the model parameters while in Sec. III all the observables that are used in our analysis are discussed. The corresponding model file required for DM phenomenology is generated with FeynRules [90]. The DM relic density, annihilation cross section, and Higgs decay width at tree-level are computed by using micrOMEGAs [91]. We perform our global scan with the Markov Chain Monte Carlo package emcee [93].
To determine the allowed parameter space at 1σ and 2σ we analyse the data obtained via emcee using the profile likelihood (PL) method. This approach allows us to remove the unwanted parameters as nuisance parameters by maximizing the likelihood over them.
After profiling the unwanted parameters we are left with the likelihood of the parameters we are interested in. The likelihood for n parameters of interest, L(x 1 , x 2 , . . . , x n ), can be integrated as where C is the smallest n-volume with probability , x k are placeholder parameters and 1/N is a normalization factor with N the result of integrating L inside C → ∞.
Our likelihood function can be modeled as a pure Gaussian distribution using the total χ 2 total in the form L = exp(−∆χ 2 total /2), ∆χ 2 total = χ 2 total − min(χ 2 total ). (22) In this study χ 2 total is made with the combination of the constraints mentioned in Sec. III and the Appendix A: In the following subsections, we will discuss our result based on the 1σ and 2σ allowed regions of the previous BNL (g − 2) µ data (orange unfilled contours) and the combined results from BNL and the recent E989 (blue filled contours). However, all other constraints given in Sec. III and Appendix A are applied. The VLML mass M ψ vs. κ . The inner and outer contours are 1σ (solid) and 2σ (dashed) confidence limits, respectively. The orange unfilled and blue filled contours correspond to the previous BNL (g − 2) µ data and the combined one from BNL and the recent E989, respectively.
In addition to (g − 2) µ , all the other constraints are considered.
A. The impact from (g − 2) µ results on κ and M ψ From Eq. (12), after fixing M H 0 and M A 0 , it is easy to note that the dependence of ∆a µ is reduced to the parameters κ , M ψ , and sin α. Furthermore, when compared with M ψ and κ , smaller values of the parameter sin α do not have an important effect on the value of ∆a µ due to the dominant term with cos α. In Fig. 2, we present the correlation between M ψ and κ . For heavier ψ ± , a larger κ is required to satisfy the measured value of ∆a µ . It is clear that, when the new combined (g − 2) µ data is applied, the 2σ shrinks significantly while the 1σ region changes slightly. The new combined (g − 2) µ data constrains these parameters to the limits κ > 1.8 and M ψ < 315 GeV with an almost linearly correlation between them.
In Sec IV D we will describe how these facts can help in searches for the VLML ψ ± at the LHC. Colors and constraints are as described in Fig. 2. In Fig. 3, we can clearly see that the difference between the blue and orange contours is small. As mentioned in the previous section, ∆a µ is not sensitive to sin α, λ HA , and κ.
On the other hand, the Higgs boson phenomenology is rather interesting for these three parameters. First of all, the Higgs boson exotic and invisible decay widths are functions of both sin α and λ HA as can be seen in Eq. (13) and (15). From the left panel of Fig. 3 we can see the limits sin α 0.15 and λ HA 0.01 . Trying to be consistent with the ATLAS data for Higgs decaying to 2b2µ, we find that κ is proportional to sin α as presented in the right panel of Fig. 3.
Considering that the decay width of A 0 to any fermion except the muon changes only by a constant mass factor and that all BRs add up to 1, we expect an universal relationship between BR(A 0 → µ + µ − ) and BR(A 0 → f f ) to be linear. The reason is that if BR(A 0 → bb) increases then BR(A 0 → µ + µ − ) decreases rapidly, before BR(A 0 → bb) reaches 1. One can see, from Eq. (6), that enhancing BR(A 0 → µ + µ − ) requires larger κ, while we also need to increase sin α in order to satisfy ATLAS 2b2µ measurement, which depends on In Fig. 4 bb) × BR(A 0 → µ + µ − )) plane. When the product of both BRs is in the range 10 −3 -4 × 10 −2 , we find the limit 0.85 < BR(A 0 → bb) < 0.89. Finally, κ and sin α are also constrained by other searches such as LHC Higgs decay (H 0 → µ + µ − and H 0 → γγ), and EDMs constraints (both electron and muon). Therefore, their values are restricted to be small as can be seen in Fig. 3.
C. The impact from DM measurements on g χ , M χ , M ψ , and sin α In general, the relevant model parameters for DM phenomenology are g χ , M χ , κ , M ψ and sin α. Since the allowed range of κ is already restricted to 1.8 < κ √ 4π, here we focus on the analysis of g χ , M χ , M ψ and sin α. In the left panel of Fig. 5, we can clearly observe that both ranges of M χ and M ψ shrink if the new measurement of (g −2) µ is applied.
Note that only ψ ± enters to the (g − 2) µ loop calculation, but DM mass M χ and M ψ can be further restricted by Planck relic density constraint. In this region (larger M χ ), the dominant channel of DM annihilation is χχ → µ ± ψ ∓ .
For the Higgs resonance region, where the χχ → f f is relevant, the presence of terms with g χ sin α will bring some regions with large g χ ∼ O(10 −1 ) into the 2σ and 1σ regions as can be seen in the right panel of Fig. 5 and the g χ column of Fig. 7.
In Fig. 6 can be explored through mono-X processes via the off-shell exchange of H 0 , A 0 [94]. In this model, the cross sections for mono-X processes are proportional to sin 2α, making this exploration channel a challenging one at the LHC. We then turn to the possible new spin-0 particle A 0 . Since A 0 → bb is the dominant decay mode, the search for H 0 → A 0 A 0 → bbbb, as shown in Refs. [95,96], is crucial to confirm or rule out the excess presented in Ref. [53]. Other option is to use pp → ZA 0 → (l + l − )(bb) to confirm the existence of A 0 with M A 0 = 52 GeV. The production cross section is about sin 2 α × 7.67 pb. The well-known jet substructure techniques of Ref. [97] can be applied to this search for A 0 → bb.
Finally, even if the model is already constrained by the search of the pair production of VLML ψ ± with multi-lepton signature as presented in the right panel of Fig. 1, we can still explore multi-b jets processes at the LHC in the near future. The signature for the single production of ψ ± is 2µ2b and the possible SM backgrounds are tt, tb, bbZ. Thanks to the larger cross sections of the ψ ± pair production, one can explore this channel by two signatures: 2b4µ and 4b2µ. The possible SM backgrounds for the former one are ttZ and The simplified DM models are common approach for DM phenomenological studies. DM and at least one mediator are two indispensable ingredients inside these models. This opens up the possibility of discovering a mediator before finding the actual DM, thus helping us narrow down the regions worth exploring and the possible interactions between DM and the SM. Motivated by a local 3.3σ deviation at M A 0 = 52 GeV in the H 0 → A 0 A 0 → bbµ + µ − search at ATLAS Run 2, we proposed a model that this spin-0 particle A 0 is a pseudoscalar mediator. Moreover, recently the Muon g − 2 Collaboration at Fermilab hinted at BSM physics with a reported 4.2σ deviation from the SM in the combined (g − 2) µ measurements.
We found that a new vector-like muon lepton (VLML) can explain both ATLAS Higgs boson exotic decay excess and (g − 2) µ . In this renormalizable DM model, we involve a Dirac fermion χ and a pseudoscalar A, both SM singlets, plus an extra VLML ψ ± .
We comprehensively constrain the parameter space from this model using LHC Higgs boson data, DM measurements and electron and muon EDMs. We found that due to the close relationship between κ , M ψ ± and (g − 2) µ these parameters are bounded as κ > 1.8 and M ψ < 315 GeV at 2σ. Both sin α and λ HA are strongly affected by limits on Higgs boson exotic and invisible decays resulting in the upper bounds sin α 0.15 and λ HA 0.01. To have a prediction consistent with BR(H 0 → A 0 A 0 → bbµ + µ − ) ∼ 3.5 × 10 −4 from ATLAS data, we found that κ is positively correlated with sin α as shown in Fig. 3. We determined the limit to which we can enhance BR(A 0 → µ + µ − ) with larger κ considering that we need to increase sin α as well to satisfy 10 Note that κ and sin α are also constrained by H 0 → µ + µ − , H 0 → γγ and electron and muon EDMs measurements, therefore, their values are quite restricted. In addition to the Higgs mediated resonance regions χχ → f f , the major DM annihilation channels in our model are χχ → µ ± ψ ∓ via H 0 /A 0 exchanges and χχ → H 0 A 0 . Thanks to the pseudoscalar mediator, the constraints from DM direct detection can be safely ignored.
In summary, the renormalizable simplified DM model presented here simultaneously explains the ATLAS Higgs boson exotic decay excess and the recently reported (g − 2) µ result.
We summarize the 1, 2 and 3σ allowed regions of all the model parameters in the triangle plot Fig. 7. Moreover, we have proposed ways to further confirm the existence of A 0 with a mass M A 0 = 52 GeV and searches for VLML ψ ± at the LHC. Additionally, DM annihilation to 2µ2b is an interesting signature for indirect detection. Here, the first muon is primary produced but the second muon together with a pair of b-quarks come from VLML decay.
The raised either electron or gamma ray spectra can be very different with the conventional DM annihilation scenario whose two final state particles carry the same energy. We will return to this in a future work.
In this appendix, we develop the contributions to electron and muon EDMs in our model. The effective Lagrangian for the lepton l EDM can be written as From the interactions in Eq. (5), we determine the explicit contributions from two-loop Barr-Zee type diagrams [100][101][102] for both electron and muon EDMs and display the resulting expressions in what follows. The electron EDM can be calculated from the two-loop Barr-Zee type diagram in the middle of Fig. 8

Two-loop Barr-Zee EDMs
with τ xy = M 2 x /M 2 y . The notation H 0 /A 0 denotes summation of the contributions from H 0 and A 0 with the upper(lower) overall sign for H 0 (A 0 ).
Similarly, the muon EDM from all the two-loop Barr-Zee type diagrams in Fig. 8 can be written as