Loopy Constraints on Leptophilic Dark Matter and Internal Bremsstrahlung

A sharp and spatially extended peak in an astrophysical gamma ray spectrum would provide very strong evidence for the existence of dark matter (DM), given that there are no known astrophysical processes that could mimic such a signal. From the particle physics perspective, perhaps the simplest explanation for a gamma ray peak is internal bremsstrahlung in DM annihilation through a charged t-channel mediator eta close in mass to the DM particle chi. Since DM annihilation to quarks is already tightly constrained in this scenario, we focus here on the leptophilic case. We compute the electromagnetic anapole and dipole moments that DM acquires at 1-loop, and we find an interesting enhancement of these moments if the DM particle and the mediator are close in mass. We constrain the DM anapole and dipole moments using direct detection data, and then translate these limits into bounds on the DM annihilation cross section. Our bounds are highly competitive with those from astrophysical gamma ray searches. In the second part of the paper, we derive complementary constraints on internal bremsstrahlung in DM annihilation using LEP mono-photon data, measurements of the anomalous magnetic moments of the electron and the muon, and searches for lepton flavor violation. We also comment on the impact of the internal bremsstrahlung scenario on the hyperfine splitting of true muonium.


I. INTRODUCTION
One of the cleanest signatures in indirect dark matter (DM) searches are peaks in the cosmic gamma ray spectrum from the Galactic Center or other regions of high DM density. On the one hand, there are no known astrophysical sources that could mimic such a signal. 1 On the other hand, gamma ray observatories are making tremendous progress in terms of statistics, resolution and control of systematic uncertainties.
From the particle physics point of view, peaks in the gamma ray spectrum can originate from DM annihilation or decay to two photons, a photon and a Z boson, or a photon and a Higgs boson. However, since DM is electrically neutral, these processes can only happen at the 1-loop level, making it likely that DM is first discovered in other annihilation or decay channels. There is, however, a class of models where the first experimental hint for DM is a gamma ray peak. Namely, this can happen in models where DM annihilates via a charged t-channel mediator, so that a photon can be emitted from the mediator, see Fig. 1. This process is called virtual internal bremsstrahlung (VIB) [2][3][4]. If the mediator mass m η and the DM mass m χ are close to each other, the resulting photon energy is strongly peaked (see Fig. 2) and can yield a line-like gamma ray signal if the width of the peak * jkopp@mpi-hd.mpg.de † lisa.michaels@mpi-hd.mpg.de ‡ juri.smirnov@mpi-hd.mpg.de 1 The authors of Ref. [1] show that a particular composition of a pulsar wind nebula could generate a peaked gamma ray signal, but an observation of a peak at the same energy in different regions of the galaxy would rule out this possibility.
is below the detector resolution.
In the present paper, we derive new constraints on leptophilic DM, and we translate these constraints into bounds on the cross section for internal bremsstrahlung. We also discuss the prospects for probing the parameter space of leptophilic DM even further with future experiments. We work in a simplified model which augments the Standard Model (SM) by a fermionic DM candidate χ and a charged scalar mediator η, with a coupling of the formχ η + h.c., where is a charged lepton field. This effective scenario can be realized in supersymmetry (SUSY) (see for instance [4]), where χ could be identified with the lightest neutralino, and η would be a slepton. It also applies to certain radiative neutrino mass models, whose direct detection phenomenology has been discussed in [46]. A simplified framework of the form used here has been employed, for instance, to explain an anomalous line-like feature at ∼ 135 GeV in the Fermi-LAT gamma ray data [4,47,48]. Even though the statistical significance of this feature is not yet convincing [48], and there are (inconclusive) indications that poorly understood systematic effects may play a role [49][50][51][52][53], it demonstrates the relevance of internal bremsstrahlung signatures as considered here if anomalous peaked features are found in future gamma ray observations. Our starting point is the observation that even in leptophilic models, loop processes endow the DM with nonzero electromagnetic moments, which in turn allow it to interact in direct detection experiments. If DM is a Majorana fermion, only an anapole moment is generated [54,55], while for Dirac fermions, also a magnetic dipole moment can exist. DM with anapole interactions has been studied previously in [56][57][58][59] using an effective field theory framework, and DM with magnetic dipole moments has been investigated in [46,57,[60][61][62][63][64]. The importance of loop processes even for hadrophilic DM has been studied in the context of LHC searches in [65].
Loop processes involving DM particles can also modify electromagnetic properties of leptons, in particular their anomalous magnetic moments and the energy levels of dilepton systems such as positronium and muonium. In the most general case, also lepton flavor violation could be induced by DM loops. Finally, if DM couples to electrons, it can be directly produced at LEP or at a future linear collider, allowing us to derive constraints from searches for mono-photons plus missing energy.
The paper is organized as follows. In Sec. II, we introduce the simplified model framework which we will use in the rest of the paper, we establish its connection to supersymmetric scenarios, and we review the expected indirect detection (internal bremsstrahlung) signals from DM annihilation in this model. In Sec. III, we compute the elec-tromagnetic form factors of DM and the resulting direct detection cross sections. We compare these to LUX [66] and XENON100 [67] data, and to the expected future sensitivity of XENON1T and LUX-ZEPLIN to derive constraints. We translate these constraints into limits on the intensity of possible internal bremsstrahlung signals. To illustrate the strength of direct detection limits, we show that for flavor-universal DM couplings to leptons, the explanation of the aforementioned 135 GeV feature in the Fermi-LAT data [4,47,48] in terms of internal bremsstrahlung is severely constrained. We then investigate in Sec. IV the complementary constraints from e + e − collider data, and in Sec. V the bounds from measurements of the anomalous magnetic moment of the muon and the electron, from searches for lepton flavor violation, and from possible future experiments on true muonium spectroscopy. We summarize our findings and conclude in Sec. VI.

II. INTERNAL BREMSSTRAHLUNG IN SIMPLIFIED MODELS
A. The simplest model The simplest theoretical models that feature internal bremsstrahlung in DM annihilation extend the Standard Model by a neutral DM candidate χ and a charged mediator η [4]. χ can be either a Majorana fermion (as in most supersymmetric theories) or a Dirac fermion (as for example in supersymmetric theories with preserved Rsymmetry [68,69]). As explained above, we are mostly interested in leptophilic models because DM couplings to quarks are already tightly constrained. In the simplest case, we thus start with the interaction Lagrangian where A µ denotes the photon field, χ is the fermionic DM candidate, e is the unit electric charge, is a SM lepton field, P R = (1 + γ 5 )/2 is the right-handed chiral projection operator, and y is the Yukawa coupling constant of the DM-lepton interaction. Unless indicated otherwise, we assume χ to be a Majorana fermion. Note that we have omitted couplings to left-handed leptons here which are more strongly constrained (though not ruled out) by collider searches and electroweak precision test [70]. We also do not consider the scalar potential for η since these terms are irrelevant to our discussion. The simplified model (1) has been introduced in [4], and it has been shown that the model could explain the 135 GeV feature in the Fermi-LAT data. The fit from [4] results in a preferred DM mass of m χ = 149 ± 4 (stat) +8 −15 (syst) GeV and an annihilation cross section σv rel χχ→ ¯ γ = (6.2 ± 1.5 +0.9 −1.4 ) · 10 −27 cm 3 s −1 . (Here, v rel is the relative velocity of the two annihilating DM particles, and the average · is taken over v rel .) The interactions in eq. (1) lead to annihilation of DM particles into pairs of SM leptons via t-channel exchange of the charged scalar η. This 2 → 2 process can be decomposed into an s-wave part and a p-wave part, the latter of which can usually be neglected because it is suppressed by the square of the small velocity v rel ∼ few × 100 km/s of DM particles in the Milky Way. The s-wave contribution is unsuppressed for Dirac DM, while for Majorana DM, it is helicity-suppressed by the small mass of the final state lepton [4]. This can be understood by noting that DM annihilation through the Yukawa interaction in Eq. (1) produces two leptons of the same chirality. For Majorana DM, however, Pauli blocking in the initial state requires the incoming DM particles to have opposite spin. Angular momentum conservation therefore requires a mass insertion on one of the final state lepton lines. Thus, for Majorana DM, higher order annihilation processes become important, in particular the 2 → 3 process χχ → ¯ γ, with two charged leptons and a photon in the final state (see Fig. 1). Since the photon carries away one unit of angular momentum, it can lift the helicity suppression, see for instance [71,72].
If the mediator mass m η and the DM mass m χ are nearly degenerate, the emission of an internal bremsstrahlung photon (first and fourth diagram in Fig. 1) is strongly peaked if the photon energy E γ gets close to m χ . The reason is that, in this case, one of the final state leptons is very soft, and one of the η propagators gets close to the mass shell. In other words, internal bremsstrahlung with m η m χ and E γ ∼ m χ can be viewed as DM annihilation into a lepton and a photon, with the emission of a soft lepton as a form of initial state radiation. While the spectral peak is thus due to internal bremsstrahlung only, it is important to take into account also the final state radiation diagrams to guarantee gauge invariance of the process. Note that, in contrast to gamma ray lines from DM annihilation to γγ, γZ or γH, the peaked signal from internal bremsstrahlung is not loop-suppressed, hence the cross section can be sizeable. The differential three-body cross section for χχ → ¯ γ in the case of Majorana DM has the following form [4] v rel dσ χχ→ ¯ γ dx with the electromagnetic fine structure constant α em , the number of final state lepton flavors N , and with the definitions x ≡ E γ /m χ and µ ≡ m 2 η /m 2 χ . In Eq. (2), we have neglected the lepton mass m and the DM velocity v rel . v rel dσ/dx is plotted in Fig. 2 for different values of µ. It is clear that, in order to have a distinct peak, a small degeneracy parameter µ 1.1 is necessary. Integrating over x, we immediately obtain also the full cross section [4] v rel σ χχ→ ¯ γ y 4 α em N 64π 2 m 2 Here, Li 2 is the dilogarithm function. The approximate expression for the relic density of Majorana DM in our toy model is [4] Ω χ h 2 0. 11 1 N 0.35 y 4 m χ 100 GeV for µ 1.2. For smaller µ, Ω χ h 2 is smaller by an O(1) factor due to coannihilations [4]. Eq. (4) implies that in the interesting parameter range 0.1 y 1, m χ 100 GeV, the model (1) naturally predicts a relic density comparable to the observed value 0.089 < Ω χ h 2 < 0.1227. Here, the quoted upper limit on Ω χ h 2 is taken from Planck [73], whereas for the lower limit, we conservatively use the WMAP value [74]. We thus account in a qualitative way for the uncertainty in Ω χ h 2 from the yet unresolved tension between different measurements of the Hubble constant H 0 = h · 100 km s −1 Mpc −1 .
For Dirac DM, v rel σ χχ→ ¯ γ is not a well-defined quantity in the limit m → 0 due to infrared divergences in the phase space region where the photon is soft or collinear with one of the leptons. For m = 0, we can evaluate v rel σ χχ→ ¯ γ numerically, see Fig. 2. 2 We find that the spectrum is entirely dominated by final state radiation and no internal bremsstrahlung peak is discernable at E γ ∼ m χ . This means in particular that no sharp spectral features are expected for Dirac DM. In the following, we will therefore use the two-body annihilation cross section v rel σ χχ→ ¯ = y 4 N 32πm 2 as a figure of merit for indirect detection of Dirac DM.

B. Extended models and connection to supersymmetry
A natural realization of scenario (1) is provided by the leptonic sector of supersymmetric extensions of the SM. There, the mediator η is the lightest slepton and the DM candidate χ is the lightest neutralino, which is given in terms of its bino (B), wino (W 3 ) and higgsino (H 0 Here, N ij are elements of the neutralino mixing matrix. The next-to-lightest slepton, as well as the squarks, are assumed to be much heavier than m η . In the MSSM, the Yukawa coupling y can be written in terms of the  unit electric charge e, the Weinberg angle θ W , and the neutralino mixing matrix element N 11 as [75] If, instead of Eq. (1), we were considering couplings to left handed leptons and their corresponding sleptons, the Yukawa coupling in the MSSM would be given by Since for conventional mechanisms of supersymmetry breaking, slepton masses of one chirality tend to be similar, we will also generalize (1) to include all three lepton flavors α and slepton flavors η α of one chirality, where α = e, µ, τ : Finally, we will also consider a more general model (which cannot be realized in the MSSM), in which couplings to both left-handed and right-handed fermions are included, and couplings are allowed to be flavor off-diagonal. The Lagrangian for this generalized toy model is Here, (y αj L/R ) are the Yukawa matrices, and η j are the mass eigenstates of the scalar mediators, of which an arbitrary nmber could exist. The index α runs over e, µ, τ , while j runs over all η j mass eigenstates.
Since our main motivation is the possibility of observing internal bremsstrahlung signals in future gamma ray observations, we will mostly focus on the case where the mass scale of the mediator(s), m η , is similar to the DM mass. The reason is that in this case the photon spectrum from internal bremsstrahlung is strongly peaked. Note that models with m η ∼ m χ are notoriously difficult to probe at colliders because the charged leptons produced in slepton decays are very soft. In the supersymmetric context, a model with nearly degenerate neutralino and slepton masses has been studied with a different goal in [76]. We will comment on this model also at the end of Sec. VI.

III. ELECTROMAGNETIC FORM FACTORS OF DARK MATTER AND DIRECT DETECTION CONSTRAINTS
We now establish the connection between indirect gamma ray signatures of DM in our toy model and direct laboratory searches on Earth. Connecting the final state fermion lines in the internal bremsstrahlung and final state radiation diagrams from Fig. 1, we obtain an effective vertex coupling the DM particle to the photon through loops of the form shown in Fig. 3. At dimension 5 and 6, the most general form of this effective interaction for a neutral fermion χ is [77]  where d M is the magnetic dipole moment, d E is the electric dipole moment, and A is the anapole moment. For Majorana DM, only the anapole term is nonzero [54,55], as can be seen by using the fact that a Majorana field is invariant under the charge conjugation operationĈ, i.e.ĈχĈ ≡ −iγ 2 χ * = χ. Applying this identity to the fermion fields in Eq. (9), it is straightforward to show that the magnetic and electric dipole terms vanish.

A. One Loop Contribution to the Electromagnetic Moments
We will now compute the loop induced electromagnetic interactions for the DM particles in our toy model Eq. (1).

Anapole moment for Majorana fermions
We begin by evaluating the diagrams in Fig. 3 to obtain the anapole form factor A in Eq. (9) for Majorana DM. For negligible 4-momentum transfer q we find Taking into account the behavior of the arctanh function when its argument approaches 1, it is easy to see that for 1 µ − 1 → 0, the anapole moment diverges logarithmically as . This behavior can be qualitatively understood by noting that, if m η m χ and q 2 0, all three propagators in the loops of Fig. 3 can be close to the mass shell simultaneously. In the limit µ − 1 1, on the other hand, the leading term in A is proportional to 1/ √ . Note that in this limit, the expression in Eq. (10) requires analytic continuation of the arctanh function into the complex plane. The dependence of A on the degeneracy parameter µ is shown in Fig. 4 (a) for y = 1 and m χ = 100 GeV.
If |q 2 | m 2 , a case that is relevant for instance in DM-nucleus scattering through loops containing electrons, the approximation q 2 → 0 underlying Eq. (10) is not applicable. In this case, it is instead convenient to set m = 0 and keep only to the leading term in ξ ≡ |q 2 |/m χ , which leads to At very small or ξ, one may wonder whether a calculation at fixed order in perturbation theory is still valid. However, in the case of interest to us, namely µ − 1 , the divergent logarithms in Eqs. (10) and (11) are at most of order 10 even for DM couplings to electrons.

Dipole moment for Dirac fermions
If χ is a Dirac fermion rather than a Majorana particle, only the two diagrams on the left in Fig. 3 exist. They generate an anapole moment A that is half as large as the one for Majorana DM, Eq. (10), and a magnetic dipole moment d M given by for q 2 → 0. The dipole moment will turn out to be numerically much more important than the anapole mo-ment in scattering processes involving Dirac DM. If m is We show results for DM couplings to electrons, muons, and tau leptons. Note that for couplings to electrons, the divergence in A is regularized by the momentum transfer q 2 rather than me because in typical DM-nucleus scattering processes, |q 2 | m 2 e . We have assumed y = 1 and mχ = 100 GeV. neglected compared to m χ , i.e. → 0, Eq. (12) simplifies to Note that, unlike the anapole moment A, the dipole moment d M is not divergent for → 0. For µ − 1 1, the leading term in d M is proportional to 1/ √ . The behavior of d M as a function of µ is shown in Fig. 4 (b).

B. Direct detection signals
In this section we will discuss the experimental limits on dark matter scattering through anapole and magnetic dipole interactions. This has been done previously at the effective field theory level for instance in Refs. [56,57,59,61,62,[78][79][80][81][82][83]. Here, we carry out a similar analysis using the latest LUX [66] and XENON100 [67] data, and we then translate the resulting constraints into new limits on the expected indirect detection signals in our toy model. Since the differential DM-nucleus scattering cross section dσ/dE r (where E r is the nuclear recoil energy) for anapole and dipole interactions differs from the conventional spin-independent or spin-dependent scenarios, we cannot directly use the published exclusion limits from LUX and XENON100, but instead have to fit the data at the event level. We do this by using a framework developed in Refs. [43,84,85], which we have extended by including LUX data and by implementing anapole and dipole interactions.
The differential cross section for DM-nucleus scattering through an anapole interaction is (cf. also [56,59]) while for dipole interactions we have [79,81,87,88] dσ dipole In both equations, the first line corresponds to scattering on the nuclear charge Z, while the second line describes scattering on the nuclear dipole moment d A . The nuclear mass is denoted by m N , and v is the velocity of the incoming DM particle. We have also included the nuclear charge form factor F Z (E r ) and the spin form factor F s (E r ). We parametrize F Z (E r ) as [89] F Z (E r ) = 3e −κ 2 s 2 /2 [sin(κr) − κr cos(κr)]/(κr) 3 , where κ = √ 2m N E r , s = 1 fm, r = √ R 2 − 5s 2 , R = 1.2A 1/3 fm (with the nuclear mass number A). For F s (E r ), we use [81] F s (E r ) = sin κR s /(qR s ) for κR s < 2.55 and κR s > 4.5, and F s (E r ) = 0.217 otherwise. Here, R s = A 1/3 . Note that nuclear dipole moments are subdominant in many target materials, including xenon, which we mostly focus on in this paper. The contribution from the nuclear dipole moment may be comparable to the contribution from the nuclear charge for instance in fluorine, sodium and iodine [88]. Note that Eq. (14) can be integrated over E r to yield a total cross section, while Eq. (15) has an infrared divergence, which makes the contribution from DM-dipole scattering would be correct only for a truly pointlike nucleus with magnetic dipole moment e/(2m N ). Here, instead, this spurious DM-dipole scattering term must be subtracted out and replaced by the correct term for scattering on the dipole moments of extended nuclei (second line of Eqs. (14) and (15)). total cross section for dipole interactions an ill-defined quantity.
The differential DM-nucleus scattering rate per unit target mass is given by where ρ 0 0.3 GeV/cm 3 is the local DM density, v min = m N E r /2/M χN is the minimal DM velocity required to yield a recoil energy E r , M χN = m χ m N /(m χ + m N ) is the reduced mass of the DM-nucleus system, and f ⊕ ( v) is the DM velocity distribution in the rest frame of the detector. We obtain f ⊕ ( v) by a Galilean transformation of the DM velocity distribution in the Milky Way rest frame, f MW . For the latter, in turn, we assume the conventional Maxwell-Boltzmann form with a smooth cutoff, , with velocity dispersion v 0 = 220 km/s and escape velocity v esc = 550 km/s. We expect the dependence of our results on this choice of velocity profile to be similar to what was found for DM scattering through contact interactions in the literature, see for instance [90][91][92][93].
In Fig. 5, we compare the differential reaction rates dR/dE r for anapole, dipole and spin-independent contact interactions, both with and without including nu-clear form factor and detector effects. For easier comparison, all rates are normalized to a total rate of 1 event above 10 keV per kg per day before taking into account nuclear form factor and detector effects. We see that anapole and contact interactions lead to similar event spectra, while dipole interactions are strongly enhanced at low energies due to the 1/E r dependence of the first term in Eq. (15). The nuclear form factor leads to a suppression of dR/dE r at higher energies. Note that at low energies, the scattering rate remains sizeable down to few keV even because such low energy events can occasionally produce a detectable number of photoelectrons due to Poisson statistics.
We conclude that with sufficient statistical power direct detection experiments could relatively easily distinguish dipole interactions from other interaction structures, while discriminating between anapole and contact interaction is challenging.
In the absence of a signal, we next derive limits on the anapole moment A, the dipole moment d M and the total DM-nucleon scattering cross section for contact interactions, σ χp .

C. Constraints from direct detection data
In Fig. 6 we show the constraints on the anapole and magnetic dipole moments of dark matter from 85.3 days of LUX data [66] and from 225 days of XENON100 data [67]. For the statistical analysis, we have used Yellin's maximum gap method [94]. The code employed to derive limits has been developed in [43,84,85]. Note that the qualitative shape of the exclusion curves is similar to the well-known exclusion limit for scattering through contact interactions. At low DM mass, the loss in sensitivity is slightly less steep for dipole interactions due to the enhancement of the scattering rate at low energies (see Fig. 5).
We now derive our main results by translating the LUX constraint on the anapole moment from Fig. 6 (a) into a constraint on the annihilation cross section σv rel χχ→ ¯ γ into two charged leptons plus an internal bremsstrahlung photon using Eqs. (10) and (3). Similarly, we convert the LUX limits on the dipole moment of Dirac DM from Fig. 6 (b) into bounds on the DM annihilation cross section into two charged leptons, σv rel χχ→ ¯ using Eqs. (12) and (5). Note that the total cross section for the 3-body final state ¯ γ is ill-defined in the Dirac case due to infrared divergences. Moreover, annihilation into ¯ γ is subdominant for Dirac DM. Our results are shown in Fig. 7  For Majorana DM, Fig. 7 (a) clearly reflects the increase in the anapole moment for small = m 2 /m 2 χ , which here translates into stronger limits on the model parameters and on σv rel χχ→ ¯ γ for coupling to electrons than for coupling to µ or τ . We also clearly see the effect of degenerate m η and m χ : for µ = m 2 η /m 2 χ close to unity, the anapole moment is significantly larger than for well separated m η and m χ (see Eq. (10) and Fig. 4). Comparing to the preferred parameter region from the gamma ray line search in [4], we find that this region is still marginally compatible with direct detection constraints if µ = 1.1. For µ = 1.01, it is disfavored at the 5σ confidence level if DM has couplings to electrons and at the 3σ confidence level if DM couples only to muons. Comparing to the cross section required for thermal relic DM (horizontal blue band in Fig. 7 (a)), we see that direct detection limits are just starting to probe this region. Note that our estimates for the thermal relic cross section are based on Eq. (4). They do not include the effect of co-annihilations [4], which would move the thermal relic cross section to smaller values. Note also that our perturbative calculations become inaccurate close to the gray regions in Fig. 7, inside of which y 2 is larger than 4π.
Comparing direct detection constraints to limits from gamma ray searches ( Fig. 7 (b)), we find that for flavoruniversal couplings and µ not too far from unity, direct searches are significantly more sensitive than continuum gamma ray searches in dwarf galaxies [4] and competitive with the bounds from gamma ray line searches [9]. (Note that in Refs. [4,9] these bounds are shown only for m χ 50 GeV, even though in principle, Fermi-LAT and H.E.S.S. are sensitive also to lower DM masses.) At m χ 10 GeV, direct detection limits are superseded by constraints from the anomalous magnetic moment g − 2 of the electron and the muon (see Sec. V A).
Looking into the future, Fig. 7 (c) illustrates that the sensitivity of direct detection experiments can be expected to improve by more than two orders of magnitude in the coming years thanks to the planned XENON1T [95] and LUX-ZEPLIN (LZ) [96] experiments. This will make direct DM searches highly sensitive to thermal relic DM. For XENON1T, we have assumed a total exposure of 2 200 kg yrs, while for LZ we use 10 000 kg yrs. In both cases, this corresponds to roughly 2 years of data taking. For comparison, we plot in Fig. 7 (c) also contours of constant σ ana (gray dotdashed curves), where σ ana is the direct detection cross section averaged over the DM velocity distribution: Note that direct detection limits on σ anapole χN are more than an order of magnitude weaker than direct detection limits on the cross section for DM-nucleon scattering through contact interactions. The reasons are the velocity dependence in σ anapole χN as well as the fact that anapole interactions are possible only with protons, whereas contact interactions occur also on neutrons.
For Dirac DM, Fig. 7 (d) shows that the qualitative picture is similar to Majorana DM, but the dependence on the lepton mass m is less strong. Comparing the direct detection limits to constraints from the Fermi-LAT analysis of gamma ray signals from dwarf galaxies [6], we find that for DM masses > 10 GeV, direct detection provides significantly stronger limits if m χ and m η are not too different. In this case, also thermal production (horizontal blue band in Fig. 7 (d)) is excluded for 10 GeV m χ few × 100 GeV.

IV. COLLIDER SEARCHES FOR LEPTOPHILIC DARK MATTER
A set of constraints on leptophilic DM complementary to the limits from direct detection can be obtained from collider data. Since tree level production of DM at hadron colliders [10][11][12][97][98][99][100] is impossible in the leptophilic case, the strongest constraints are expected to come from mono-photon events at LEP [101]. In the future, mono-photon searches at a linear collider may improve on these bounds [102].
Here, we apply the procedure described in [101] to our toy model, Eq. (1). We simulate the process e + e − → χχγ in CalcHEP 3.4 [103] including the effect of initial state radiation and beamstrahlung (with default parameters) on the beam energy. We analyze the simulated events in a modified version of MadAnalysis 1.1.2 (from the MadGraph 4 package) [104] that implements the efficiencies and resolutions of the DELPHI detector at LEP [105,106], see [101] for details. We have checked that our simulation reproduces the predicted γνν background from [106] to very good accuracy. To set limits, we add our signal prediction to the background prediction from [106], and compare to the DELPHI data from Fig. 1 of [106], which corresponds to an integrated luminosity of 650 pb −1 . Following [101] we use a simple χ 2 analysis to set limits on the Yukawa coupling y as a function of m χ and m η , and then convert these limits into constraints on σv χχ→ ¯ γ , which are shown in Fig. 7 (b) and (d). Systematic uncertainties are subdominant compared to statistical uncertainties in DELPHI and are therefore neglected in our analysis.
We also estimate the sensitivity of a future linear collider with a center of mass energy √ s = 500 GeV to leptophilic DM in our toy model. We simulate the signal and the dominant γνν background in CalcHEP 3.4 [103] while for the γγνν final state (with one photon escaping undetected) and for γe + e − events (with an undetected e + e − pair) we follow [102]: we qualitatively include the γγνν background by simply increasing the γνν background by 10%. For γe + e − events, we reweight the γνν spectrum by the energy dependent factor 0.825 [1 − E/(0.9 GeV)] 2 . Negative reweighting factors are excluded. The detector response of an ILC detector is modeled according to the information given in [102,107,108]. We assume an energy resolution of ∆E/E = 0.011 ⊕ 0.166/ E/GeV, where the notation ⊕ means that the different terms correspond to separate, statistically independent Gaussian distributions. We restrict our analysis to the photon energy range 10 GeV < E γ < 220 GeV (divided into 5 GeV bins) to remove events with on-shell Z production, and to the rapidity range |y| < 2.3. The detection efficiency is , and only to taus (thin dotted). The upper and lower boundaries of the colored bands correspond to µ ≡ m 2 η /m 2 χ = 1.1 and µ = 1.01, respectively. For illustration, we also show the approximate cross section required for a thermal relic (roughly estimated using Eq. (4)), and the tentative best fit region from Bringmann et al. [4]. The gray region corresponds to y 2 > 4π and thus cannot be reached in our toy model. In (b), we compare the LUX bounds on Majorana DM with flavor-universal couplings to limits from LEP mono-photon searches (see Sec. IV), g − 2 measurements (see Sec. V A), a Fermi-LAT search for continuum gamma rays from dwarf galaxies [4], and Fermi-LAT (solid) and H.E.S.S. (dotted) searches for gamma ray lines from the Galactic Center [9]. In (c), we project the future sensitivities of ton-scale direct detection experiments and of a future linear collider for Majorana DM with flavor-universal couplings and with µ = 1.01. For illustration, we have also drawn contours of constant velocity-averaged direct detection cross section σχp (see Eq. (17)). In (d), we summarize direct detection constraints induced by magnetic dipole interactions for Dirac DM with flavor-specific couplings, and we compare again to the thermal relic cross section, to LEP mono-photon limits, and to Fermi-LAT limits from dwarf galaxies [6]. Note that no sharp features are expected in the gamma ray spectrum from annihilation of Dirac DM.
given by 0.941 − 0.00129E γ /GeV. We derive limits using a simple χ 2 analysis, assuming an integrated luminosity of 50 fb −1 and neglecting systematic uncertainties. Our projected ILC limits are included in Fig. 7 (c) and (d).

V. CONSTRAINTS FROM PRECISION EXPERIMENTS
A. Lepton magnetic dipole moments The extension of the SM by a DM particle and a charged mediator in our toy model Eq. (1) leads to a new contribution to the anomalous magnetic moment (g − 2) of leptons via the vertex correction loop shown in Fig. 8. This has been used previously in [4,77] to constrain DM annihilation through charged mediators. In the case of complex Yukawa couplings, there can also be contributions to electric dipole moments, but we will not consider this possibility here. In the limit m m η , m χ , the anomalous magnetic moment of charged leptons is modified by [4] ∆a For DM couplings to electrons, we compare Eq. (18) to the difference between the SM prediction for a e and the experimentally measured value, a exp e − a SM e = (−1.06 ± 0.82) × 10 −12 [109] to derive the exclusion bound shown in Fig. 7 (b) for µ = 1.1 (lower edge of colored band) and for µ = 1.01 (upper edge of colored band).
For the g − 2 of the muon, the difference between the measured best fit value and the theoretical prediction is a exp µ − a SM µ = [2.87 ± 0.63 (exp.) ± 0.49 (theor.)] × 10 −9 [110]. We add the experimental and theoretical uncertainties in quadrature. To account for the significant discrepancy between theory and experiment, we artificially inflate the error by linearly adding an ad-hoc uncertainty given by the central value of the discrepancy, 2.87 × 10 −9 . Note that the discrepancy has a sign opposite to the one predicted by Eq. (18). The resulting constraint on σv rel χχ→ ¯ γ is shown in Fig. 7 (b).
We see that g − 2 constraints are competitive with direct and indirect searches only at DM masses < 10 GeV.

B. Positronium and muonium spectroscopy
Lepton-antilepton bound states such as positronium (e + e − ) and true muonium (µ + µ − ) are interesting laboratories for precision tests of QED because they can be studied accurately using spectroscopy, but are theoretically simpler than atoms. In particular, there are no nuclear effects that need to be taken into account. In our toy model for leptophilic DM, the box diagrams shown in Fig. 9 lead to an effective contact interaction of the form with This contact interaction contributes to the electrostatic potential between the + and − , thus modifying the hyperfine splitting E hfs between the energy of the orthostate (parallel spins, 3 S 1 ) and the para-state (antiparallel spins, 1 S 0 ). To obtain the new contribution ∆E hfs to E hfs , we first calculate the new term in the Hamilton operator of the system by plugging explicit expressions for the + and − wave functions into (19), integrating over d 3 x and adding a minus sign from the Legendre transform that converts the Lagrangian into the Hamiltonian as well as a factor 4 from the different ways in which the lepton fields can be contracted with the incoming and outgoing fermion states. The lepton wave functions are given by where ξ is a non-relativistic particle or antiparticle Dirac spinor normalized to unity. We find that the energy of the ortho-state remains unchanged while the energy of the para-state is increased. The splitting between the two states is thus reduced, with which is an O(10 −12 ) correction to E e + e − hfs = [203.3941 ± 0.0016 (stat) ± 0.0011 (syst)] × 10 9 Hz [111], well below the experimental precision and the precision of the SM prediction. The reason for the low sensitivity is that positronium is a relatively large system, whereas the contact interaction is effective only at very short distance. The same is true for e ± µ ∓ bound states.
More promising as a probe for contact interactions of the form of eq. (19), and of new physics in the lepton sector in general, seems to be "true muonium", i.e. a µ + µ − bound state. Even though true muonium has never been directly produced and studied in the laboratory, precision experiments seem feasible [112]. For true muonium, we have which is only an O(10 −7 ) correction to the leading term E µ + µ − hfs 4.23 × 10 7 MHz [113]. Using Eq. (4) and comparing to Eq. (24), we obtain that to exclude thermal relic dark matter with m χ = 130 GeV, µ = 1.1, E µ + µ − hfs needs to be measured with an accuracy of 0.2 MHz.

C. Lepton Flavor Violation
Even though in the simplest versions of our toy model motivated by supersymmetry, Eqs. (1) and (7), DM couplings to leptons are flavor diagonal, we now consider also the general Lagrangian Eq. (8) including flavor offdiagonal couplings. We derive constraints on these couplings from searches for the rare decays µ → eγ, τ → eγ and τ → µγ, which are mediated by the diagram shown in Fig. 8. Computing this diagram, we obtain for the decay rate where are Wilson coefficients in the effective Lagrangian and the loop functions J(µ), I(µ) are given by We have used the definition µ j ≡ m 2 ηj /m 2 χ , where m ηj are the masses of the charged mediators (see Eq. (8)).
We consider for illustrative purposes the special case where only three charged mediator η 1 , η 2 , η 3 exist, and where y L = 0. This can be realized in supersymmetry if all left-handed sleptons are too heavy to be phenomenologically relevant. We obtain in this special case for the branching ratios BR α → β γ Γ α → β γ /Γ SM (with the SM width Γ SM ) The expression for BR τ →eγ is identical to the one for BR τ →µγ , with the replacements y µj * R → y ej * R . With the current experimental limits BR µ→eγ < 5.7 × 10 −13 [114], BR τ →µγ < 4.4 × 10 −8 [115] and BR τ →eγ < 3.3 × 10 −8 [115], and using m χ = 100 GeV, we then obtain the following constraints on the elements of y R at µ = 1.1:

Process
Coupling Limit µ → eγ j (y µj R y ej * R ) 2 1/2 < 1.0 × 10 −4 τ → µγ j (y τ j R y µj * R ) 2 1/2 < 7.0 10 −2 τ → eγ j (y τ j R y ej * R ) 2 1/2 < 6.1 10 −2 We have seen in Eq. (4) that in our simplified model setup, at least one of the Yukawa couplings should be of order 0.1-1 to avoid DM overproduction. The above constraints show that flavor off-diagonal Yukawa couplings are therefore always subdominant. This justifies our neglecting them in the preceding sections. We have also studied the decay µ → 3e, which constrains a different combination of Yukawa couplings because it also receives contributions from box diagrams similar to Fig. 9. If we assume that flavor-diagonal Yukawa couplings are O(1), we obtain limits on the flavor off-diagonal couplings that are about a factor of 8 weaker than the limit from µ → eγ. To arrive at this estimate, we have used Ref. [116] to express BR(µ → 3e) in terms of the Wilson coefficients of the effective operators in Eqs. (28) and (19). We have then compared the predicted branching ratio to the current experimental limit from [110,117]. Note that planned searches for µ → 3e will improve the limit on BR(µ → 3e) by up to four orders of magnitude [118].

VI. CONCLUSIONS
In this paper, we have studied leptophilic dark matter models in which DM annihilation proceeds through a charged mediator and can therefore be accompanied by emission of a virtual internal bremsstrahlung photon. Such models are of great interest for indirect dark matter searches because internal bremsstrahlung can lead to spectral peaks in the gamma ray sky, a feature which is easily distinguishable from the large astrophysical gamma ray flux. Leptophilic DM models are also well motivated theoretically: they can be realized for instance in supersymmetric scenarios or radiative neutrino mass models, and in most cases, their parameter space is relatively unconstrained.
Here, we have established a connection between internal bremsstrahlung signals and loop-induced electromagnetic form factors of DM particles in leptophilic models. In particular, upon connecting the charged lepton lines in the internal bremsstrahlung diagrams in Fig. 1 to a loop, one immediately obtains the electromagnetic vertex corrections in Fig. 3. For Majorana DM, these lead to an anapole moment, while for Dirac DM, both anapole and magnetic dipole moments are generated, with the dipole moment being dominant in DM scattering processes. Interactions of the anapole and dipole moments with atomic nuclei then allow us to constrain the internal bremsstrahlung cross section using DM-nucleus scattering data from direct detection experiments. We have carried out this analysis for the most recent LUX and XENON100 data, and have found that direct detection constraints can be competitive with internal bremsstrahlung searches. This is true in particular if the mass splitting between the DM particle χ and the charged mediator η is very small-the case which is also most interesting for internal bremsstrahlung searches due to the peaked gamma ray spectrum.
If DM is a Majorana fermion that couples universally to all charged leptons, direct detection limits are of the same order as limits from gamma ray line searches, and better than continuum gamma ray constraints from dwarf galaxies (see Fig. 7 (b)). Specifically, for small mass splitting m 2 η /m 2 χ 1.1, LUX constrains the internal bremsstrahlung cross section σv rel χχ→ ¯ γ to be below few × 10 −28 cm 3 /s at m χ ∼ 20 GeV. At DM masses of order 100 GeV, which have been invoked previously to explain a bump in Fermi-LAT gamma ray data [4], LUX constraints imply that this interpretation is disfavored if DM couples to electrons or muons and if m η and m χ differ by few %. If the last condition is significantly violated, however, the expected bump in the gamma ray spectrum becomes relatively broad, making line searches less sensitive. If m η /m χ 1, also direct searches for the charged mediator η at colliders will impose impor-tant constraints, disfavoring m η few × 100 GeV [4,70]. These constraints are ineffective if m η ∼ m χ because the leptons from η decay will be very soft in this case and thus hard to detect.
We note an interesting connection between our results and the scenario studied by Konishi et al. [76] to solve the cosmological lithium-7 problem in the Constrained Minimal Supersymmetric Standard Model (CMSSM) with sleptons that are nearly mass degenerate with the lightest neutralino. For the preferred mass range from [76], 300 GeV m χ 500 GeV, this scenario would predict σv rel χχ→ ¯ γ ∼ 10 −28 cm 3 /s, well within the region testable by next generation direct detection experiments.
If DM is a Dirac fermion and the masses of χ and η are of the same order of magnitude, direct detection constraints disfavor thermal relic production of DM for m χ between 10-20 GeV and up to a few hundred GeV (see Fig. 7 (d)). For m χ > 20 GeV, direct detection limits are also significantly stronger than astrophysical limits from gamma ray line searches and from continuum gamma rays searches in dwarf galaxies.
In the future, we expect the XENON1T and LUX-ZEPLIN experiments to improve these direct detection limits by about two orders of magnitude. These experiments will thus test the thermal relic hypothesis for DM masses of order 10 GeV m χ few × 100 GeV. If a signal is detected, the spectrum of recoil events can be used to discriminate between anapole and dipole interaction and hence between Majorana and Dirac DM.
We have also studied constraints on our simplified model from low energy precision experiments. We confirm that bounds from the anomalous magnetic moment g − 2 of the electron and the muon are weaker than the direct detection constraints at m χ 10 GeV. Searches for the lepton flavor violating decays τ → µγ, τ → eγ, µ → eγ and µ → 3e are very powerful in setting bounds on DM annihilation into flavor violating final states. Finally, we have studied the possibility of obtaining constraints from a future measurement of the hyperfine splitting in true muonium (a µ + µ − bound state). We have found such a measurement to be challenging for heavy DM (m χ ∼ 100 GeV), where excluding thermal relic DM would require a measurement with a relative accuracy better than 10 −7 (see Eq. (24)). For lighter DM (m χ 10 GeV), however, requirements are weaker and an interesting measurement may be possible.
In summary, our results show that direct dark matter searches are powerful tools to search for leptophilic DM even though DM-nucleus scattering occurs only at the loop level in this case. They are complementary to, and sometimes significantly superior to, indirect searches and precision experiments. Particularly in a scenario where a peak is observed in the cosmic gamma ray spectrum, but no other indirect hints for DM are found, virtual internal bremsstrahlung in a leptophilic DM model provides an attractive explanation. Our results show how this scenario can be confirmed in direct detection experiments by looking for the electromagnetic moment interactions of DM with nuclei. This illustrates once again that the search for Dark Matter is an interdisciplinary task, and that only a combination of different search strategies can yield optimal results.