Faint Light from Dark Matter: Classifying and Constraining Dark Matter-Photon Effective Operators

Even if Dark Matter (DM) is neutral under electromagnetism, it can still interact with the Standard Model (SM) via photon exchange from higher-dimensional operators. Here we classify the general effective operators coupling DM to photons, distinguishing between Dirac/Majorana fermion and complex/real scalar DM. We provide model-independent constraints on these operators from direct and indirect detection. We also constrain various DM-lepton operators, which induce DM-photon interactions via RG running or which typically arise in sensible UV-completions. This provides a simple way to quickly assess constraints on any DM model that interacts mainly via photon exchange or couples to SM leptons.


Introduction
The origin and nature of Dark Matter (DM) remains a key puzzle of modern physics [1][2][3] and a varied program is underway, searching for particle DM through direct detection [4], indirect detection [5] and at colliders [6]. Though referred to as 'Dark', DM may still have couplings to the Standard Model photon. This coupling may be at tree-level, in which case the DM is electrically charged [7], a possibility which has received much attention as a possible explanation of the EDGES 21cm anomaly [8][9][10][11]. Such a renormalisable coupling to photons would be valid up to an arbitrarily high scale, though the DM in this case must be 'milli-charged' [12][13][14][15][16], with a charge substantially smaller than that of the electron [17,18]. Alternatively, DM may couple to the photon through higherdimensional, non-renormalizable operators, which arise at low energy after integrating out heavy Beyond-the-Standard-Model (BSM) degrees of freedom.
Such higher-dimensional DM-photon interactions allow us to constrain the scale of New Physics beyond the Standard Model, which is the aim of the current work. There are a wide range of these interactions and it is not always clear which constraints should be the most relevant for a given interaction or which interaction is most constraining for a given UV model. Our philosophy in this work will therefore be to classify these DM-photon operators and present updated constraints on the ones most relevant for phenomenology. We classify these operators as follows: in terms of their dimension, anticipating that the lowest dimensional operators will typically dominate; in terms of the nature of the DM candidate (Majorana or Dirac fermion, as well as real or complex scalar); and in terms of the CP parity of the interactions. Such a classification allows one to determine which operators (and therefore which bounds) are the most relevant for a given UV model, without detailed calculation. This classification will be a useful tool in constraining New Physics beyond the Standard Model through DM-photon couplings.
Many DM-photon operators have been discussed previously in the literature. These include electric and magnetic dipoles (e.g. [19][20][21][22][23]), anapole (e.g. [19,24,25]) and Rayleigh interactions (e.g. [26][27][28][29]). However, it is useful to update and summarize the current bounds from direct detection (DD) and indirect detection (ID). Moreover, we consider simple UV completions in which DM couples to a heavy particle and a lepton. We provide completely general expressions for the Wilson coefficients, and use the bounds on the effective operators to constrain the parameters of this model, which are relevant in any scenario in which DM couples to charged particles. These so-called 'lepton portal' models have been studied previously [30]; we update the bounds using the latest results from direct, indirect and collider DM searches. We also aim to carefully make the connection between such simple UV models and DM-photon interactions in effective field theories. In some cases, the strongest constraints come from DM-photon operators which are expected to be sub-dominant based on naïve dimension-counting. In other cases, the dominant constraints do not come from DM-photon interactions themselves but instead from other operators appearing in sensible UV completions.
The paper is structured as follows. We begin in Sec. 2 with a classification of nonrenormalisable DM-photon interactions; the operators we focus on are summarized in Ta-ble 1. In Sec. 3, we present bounds on the coefficients of these operators from direct and indirect detection experiments. We also present bounds on a number of 'UV-related' operators, which may appear in typical UV-complete models. In Sec. 4, we present two simple UV models which give rise to DM-photon interactions at low energy. We also map these UV models onto the non-renormalisable operators presented in Sec. 2, allowing us to compare experimental bounds on such UV models coming from a full calculation and from an effective field theory approach. Finally, in Sec. 5, we summarise our results and present our conclusions.

Classifying DM-photon Interactions
Our aim is to classify and study models of Dark Matter (DM) in which possible signals in direct and indirect detection arise through photon exchange. We begin in Sec. 2.1 by listing the DM-photon interactions which produce signals at the energy scales relevant for direct and indirect searches. These are described by an effective field theory (EFT) that is invariant only under the electromagnetic gauge symmetry, dubbed EMSM χ [31]. In general a given DM-photon interaction in EMSM χ can be induced in two ways, either by a treelevel matching or by renormalisation group (RG) running. In Sec. 2.2, we will therefore also identify all operators in EMSM χ that i) radiatively induce DM-photon interactions and ii) are constrained dominantly by photon exchange at low energies. Finally in Sec. 2.4, we identify a number of additional operators which do not lead to DM-photon interactions but which are expected to appear in any UV-complete model due to SM gauge invariance. In some cases, these 'UV-related' operators may lead to stronger bounds that those from photon exchange.

Photon Operators
Here, we list all possible operators that couple DM to photons up to dimension 7. We organise the operators according to the nature of the DM particle. We first discuss "Majorana fermion" operators (χ), which are allowed in the case of Majorana DM. We then discuss "Dirac fermion" operators (ψ), which are additional operators which vanish for a Majorana fermion. Similarly, we distinguish between "real scalar" operators (R) and "complex scalar" operators (S), which vanish for a real scalar. All operators can only be induced at loop level, so we include suitable loop factors of e/(16π 2 ) or e 2 /(16π 2 ) in the Lagrangians. The UV scale Λ denotes the scale of generic heavy physics that has been integrated out, above which the EFT is expected to break down.
Majorana Fermion: For Majorana DM, only operators with scalar, pseudo-scalar and axial-vector bilinears do not vanish identically. The Lagrangian that describes the interac-tions with photons is given by whereF µν = 1 2 µνρσ F ρσ . The leading operators we consider here are labelled with underbraces. The term in the first line is the Anapole moment, which is the only allowed dimension-6 operator and describes the lowest order interactions of Majorana DM with photons. There are also four dimension-7 Rayleigh operators, which can be divided into CP-even (second line) and CP-odd operators (third line).
Dirac Fermion: For Dirac DM, alo tensor and vector bilinears are allowed. In addition to the interactions described in Eq. (2.1) one can write three other leading order operators. The Lagrangian is where σ µν = i 2 [γ µ , γ ν ]. The first and the second terms are the dimension-5 magnetic and electric dipoles respectively, while the third one is the dimension-6 fermion charge radius operator. The factor 2 appearing in the first line has been chosen so that the Feynman rules for Majorana and Dirac fermions are the same.
Real Scalar: For real scalars, the lowest order interactions with photons are described by two dimension-6 Rayleigh operators. The Lagrangian is given by Complex Scalar: For complex scalars, in addition to the dimension-6 Rayleigh interactions in Eq. (2.3), one can directly couple a complex scalar current to the electromagnetic field strength tensor, which is the scalar charge radius. The Lagrangian reads As for fermion DM, the factor of 2 multiplying L RealScalar ensures the same Feynman rules for complex and real scalars.

Running to the Nuclear Scale
The Lagrangians in the previous section are defined at the UV scale Λ, which denotes the scale of generic heavy physics that has been integrated out. From this scale we should run the theory down to the relevant IR scale µ IR (e.g. the nuclear energy scale µ n ∼ 1 GeV for direct detection 1 ). Since all relevant photon operators can only arise at loop level at the energy scale Λ, we should also take into account the RG contribution of possible treelevel operators that can mix into the photon operators. Focusing on operators that are constrained dominantly by photon exchange, we restrict to tree-level operators that do not involve quarks and Higgs fields. Operators with quark fields, even if not valence quarks, will always induce DM-gluon couplings at 1-loop [34]. On the other hand operators with Higgs fields will give rise to DM couplings to valence quarks already at tree-level due to Z-boson exchange [31]. Thus, the only operators that are relevant for our purposes are tree-level operators that involve vector currents of SM leptons. We therefore consider also the following additional operators: Note that to avoid strong constraints from lepton flavor-violation, we will assume a coupling only to a single DM lepton ( = e, µ, τ ) at a time.
These operators induce fermion charge radius, anapole and scalar charge radius operators through RG evolution, according to ). The precise value is calculated using the runDM code [35], which takes into account the running of couplings between Λ and the nucleon level [31][32][33]. We have checked that there is a very good numerical agreement between runDM and the analytic expression for K(Λ) given above 2 .
Since the operators in the UV can only arise at loop-level themselves, the logarithmically enhanced contribution from the tree-level operators with leptons typically dominates. In the next section we will demonstrate that this enhancement of the radiatively-induced operators can be significant, particularly in the case of Majorana fermion DM, and lead to stronger direct detection bounds than on the fermion charge radius, anapole and scalar charge radius operators themselves.

Summary
We collect the relevant operators we have discussed so far in Table 1. We arrange the operators according to the type of DM particle and, for a given DM particle, we distinguish between those operators which are CP-even (top row) and CP-odd (bottom row). We also order the operators according to their dimension (Dim. 5, Dim. 6 or Dim. 7). We focus on operators that lead to the strongest constraints, contributing to signals in either direct detection or indirect detection. Therefore, when presenting bounds in Sec. 3, we will neglect operators that have the same CP quantum numbers as the dominant operators but give weaker constraints. That is, for two operators with the same dimension and CP properties, we present bounds only on the more strongly constrained of the two; the subdominant operator we denote with round brackets in Table 1. In particular, this concerns the fermion charge radius, the anapole and the scalar charge radius operators, which at low energy are dominantly generated from the tree-level operators via RG evolution, denoted by −→ in the table. In addition, as we discuss in more detail later, some of the Rayleigh operators give sub-leading constraints (and in some cases give no constraints at all).

UV-related Operators
Up to now we have taken into account only operators that are dominantly constrained by photon exchange at low energies. However, in a UV-complete model SM gauge invariance implies that these operators are necessarily accompanied by additional operators that are potentially more strongly constrained than the operators in Table 1. In this section we briefly comment on this issue and discuss the possible operators arising in UV-complete 2 Running effects also play a role in the calculation of the direct detection signal for Rayleigh interactions (∝ F µν Fµν , F µνF µν ), as we discuss in more detail in Appendix A.2.

CP
Dim. 5 Dim. 6 Dim. 7 models respecting eletroweak gauge invariance, which we therefore refer to as "UV-related" operators. First of all, due to electroweak symmetry it is clear that all photon operators will be accompanied by operators involving the Z-boson. However, such operators typically give effects in DD and ID which are subleading compared to the corresponding photon operators, since they are suppressed by the Z-boson mass appearing the propagator 3 . We will therefore not explicitly consider the bounds on these kinds operators in the following.
The other class of operators that typically arise in UV models are tree-level operators that couple the DM field to lepton currents. For example, the operator O ψψ = ψγ µ ψ γ µ generically comes together with the operator ψγ µ ψ γ µ γ 5 , since only left-handed and righthanded leptons have definite SM quantum numbers. In order to be complete, we will 3 In direct detection this is a general statement, while there may be exceptions for indirect detection.
In particular, in cases where the DM-photon coupling gives a gamma-ray line signature, the corresponding coupling with the Z-boson will produce a gamma-ray continuum extending down to energies below the DM mass. In such cases, it may be possible to constrain heavier DM (at a given gamma-ray energy) in the case of DM-Z interactions. See also Ref. [26]. consider all possible tree-level operators with a similar structure to the ones in Eq. (2.5), i.e. for fermion DM all dimension six operators with two DM fermions and two leptons, and for scalar DM all dimension six operators with two DM scalars, two leptons and one derivative. According to this approach we define Using the EOM for the lepton, / D = −im , the scalar operators involving derivatives can be re-written as (2.8) Finally we add the four-fermion operators with scalar currents, which must involve additional fermion masses if the UV completion does not contain additional sources of chiral symmetry breaking 4 . Under this assumption these operators arise only at dimension-8, and are given by Note that all of the operators in this section are CP-even, apart from those which contain one pure axial bilinear (e.g. ψψ · γ 5 ).
In the next section we will show that the bounds on these additional UV-related operators are rarely stronger than the ones in Table 1 when considering DD (where their contributions are loop-suppressed). In ID, however, they may be relevant as they induce tree-level annihilation of the DM into SM leptons.

Bounds on EFT Operators
We now present bounds from direct and indirect detection on the DM-photon operators described in the previous section. We also briefly discuss bounds on 4-fermion effective operators from Bhabha scattering at LEP-2.
Direct Detection: For direct detection, we show current bounds from the full LUX exposure ("LUX (2016)") [36], the one tonne-year XENON1T exposure ("XENON1T (2018)") [37] and the projected sensitivity of the LZ experiment ("LZ-projected") [38,39]. Similar constraints for high mass DM come from PANDAX-II [40]. Here, we will focus on DM at the GeV-scale and above, so we do not consider a number of constraints which are relevant only at low mass [41][42][43]. Details of direct detection cross sections and bounds are given in Appendix A. In all cases, we calculate an approximate single bin Poisson limit using the approach of Ref. [44] to translate the bounds between different interactions. Details of reproducing the bounds from LUX can be found in Appendix B.2 of Ref. [45] and details of calculating the LZ projections are given in Appendix D of Ref. [33]. The corresponding calculation for XENON1T is given here in Appendix A.3.
Indirect Detection: Bounds from indirect detection are presented for the FERMI search for spectral lines using 5.8 years of data and an isothermal DM density profile ("FERMI 5.8 yrs (2013)") [46] and the H.E.S.S.-I search for high-energy gamma ray lines using 112 hours of data and an Einasto DM profile 5 ("H.E.S.S.-I 112h (2013)") [48]. For operators which allow DM annihilation directly into charged leptons (such as O ψψ ), we also show conservative bounds from FERMI observations of dwarf Spheroidal galaxies ("FERMI dSphs (2016)") [49]. In this case, FERMI searched for a continuum spectrum of gamma rays (rather than a line) and we differentiate between limits for different final states: e, µ or τ 6 . When the DM couples to a single F µν , we set limits by summing over all fermion final states, weighting by the photon-mediated cross section, using the publicly available likelihood tools (in the Supplemental Material to Ref. [49], available at https://wwwglast.stanford.edu/pub data/1048/). We also show conservative limits from the FERMI Galactic Halo search using an isothermal DM density profile ("FERMI GH (2013)") [51]. While this constraint is typically slightly weaker than the FERMI dSphs search for s-wave annihilation, in some cases we study operators which lead to p-wave annihilation, σv ∝ v 2 . Given the larger velocities in the Galactic Halo compared to those expected in dwarfs, the Fermi GH search may come to dominate in these cases. Note that limits for e ± final states are not reported in Ref. [51]. However, the µ ± can be taken as a conservative estimate of the e ± bound; in the latter case the gamma-ray flux is expected to also include a stronger inverse compton component (see e.g. [52,53]). Apart from the bounds obtained by the H.E.S.S.-I search for high-energy gamma ray lines, we do not consider any other limits 5 For conservative choices of the DM density profile (e.g. an isothermal sphere), the limits from H.E.S.S. become very weak (see for example the right panel of Fig. 5 in Ref. [47]). 6 Loop-induced annihilation into quarks is discussed in Ref. [50]. We have checked explicitly that these bounds are subdominant, apart from when annihilation occurs very close to the Z-pole (mDM ∼ mZ /2) in which case loop-induced annihilation may be competitive with tree-level annihilation to charged leptons.
from Cherenkov telescope arrays because they strongly depend on assumptions about the DM density profile. Annihilation cross sections for the operators we consider in this work are given in Appendix B.
Colliders: The only relevant constraint from colliders comes from measurements of Bhabha scattering at LEP-2, giving bounds on 4-fermion operators involving electrons where v = 246 GeV. These operators can in turn be induced at one-loop with 2 insertions of DM-DM-e-e operators, giving at leading order in m 2 DM /Λ 2 where we have used the UV scale Λ to cut off (quadratically divergent) loops in Euclidean momentum space. As discussed in the previous section, more realistic UV models involve chiral combinations of operators, i.e. have correlated operator coefficients These induce chiral operators involving only electrons with coefficients given by, analogously to Eq. (3.3), For the low-energy bounds on c V V , c AA , c LL , c RR we use Refs. [54,55], which have reanalyzed the LEP data in Ref. [56] to derive the constraints (assuming the presence of one operator at a time) The resulting upper bounds on the tree-level DM-electron operator coefficients are shown in Table 2. These numbers should be taken with a grain of salt, since they arise from  quadratically divergent loops, which are sensitive to the specific UV completion. Nevertheless the above procedure gives indicative bounds, which will reproduce the results of the particular UV completions that we consider in the next section. Note that these bounds are quite stringent, as a result of the preference of LEP-2 data for positive coefficients of 4-electron operators that already disfavor the SM, while the new loop contributions are always negative (in particular for LL operators, where at 95% CL the coefficient is still positive, so we decide to set a bound on 99% CL). We do not show these bounds explicitly in the rest of Sec. 3 (though we will in Sec. 4). However, we remind the reader that these rather stringent bounds apply for all DM masses in the scenario where DM couplings to electrons at tree-level.
EFT Validity: Note that in principle, we should be careful about the validity of the EFT formalism, especially in the case of indirect detection. In this case, the typical energy scale √ s ∼ 2m DM can easily be comparable to the EFT cut-off Λ. When √ s > Λ, the EFT is no longer valid and we must worry about the internal structure of the EFT operators defined in Sec. 2.1 (which would lead us to considering more UV complete models, as discussed in Sec. 4). However, we do not mark this invalid region in the following figures, because they will depend on the couplings involved. In the EFT formalism, we can only constrain the combination Λ/C (1/n) . If the coupling C is small (as for weakly-coupled theories), then the constrained value of Λ is reduced and the breakdown of the EFT sets in at smaller DM masses. Instead, for strongly coupled theories, the EFT may be valid up to much higher DM masses. With this caveat, we now present constraints on the DM-photon effective operators.

Dirac Fermion Dark Matter
In Fig. 1, we show constraints on the 4 most relevant high-scale operators which give rise to low energy interactions between photons and Dirac fermion Dark Matter, ψ.
In the top left panel of Fig. 1, we show constraints on the magnetic dipole operator, O ψψF . This operator gives rise to a coherently enhanced direct detection cross section through the coupling of the DM magnetic dipole to the nuclear charge, in addition to a spin dependent dipole-dipole interaction. These interactions have both a long-range contribution and a contact contribution and lead to constraints on the EFT cut-off at the level of Λ/C ψψF 8000 GeV for m DM ∼ 50 GeV. The cross section for Dark Matter annihilation into 2 photons scales as Λ −4 , meaning that indirect detection bounds are suppressed relative to those from direct detection (where σ ∼ Λ −2 ) 7 .
In the top right panel of Fig. 1, we show constraints on the electric dipole operator, O ψ5ψF . This too gives rise to coherently-enhanced spin-independent scattering with the nucleus, though in this case the interaction is purely long range and therefore enhanced at low recoil energies. This results in strong constraints from direct detection, Λ/C ψ5ψF 7 × 10 6 GeV. Constraints on direct annihilation into gamma rays are the same as in the case of the magnetic dipole operator and therefore substantially weaker than direct detection constraints. Annihilation into fermions is p-wave suppressed, so that constraints from dwarf Spheroidals are negligible.
In the bottom left panel of Fig. 1, we show constraints on the 4-fermion interaction O ψψ , which induces the fermion charge radius operator O ψψ∂F at low energy. This gives rise to the standard (coherently-enhanced) spin-independent interaction, though coupling only to the proton content of the nucleus. Despite being a dimension-6 operator, the 4fermion interaction gives constraints which are not much weaker than the magnetic dipole, Λ/C 1/2 ψψ 3000 GeV. As discussed in Sec. 2.2, this is because O ψψ∂F is induced through RG evolution, giving an enhancement of ∼ log(Λ/µ) relative to the magnetic dipole (induced only at the loop level at the scale Λ). Finally, the dominant annihilation channel of O ψψ is directly into charged leptons, leading to competitive constraints from FERMI searches for diffuse gamma rays in both dwarf Spheroidal galaxies and the Milky Way Galactic halo.
In the bottom right panel of Fig. 1, we show constraints on the 4-fermion interaction, O ψ5ψ . In this case, RG evolution induces the anapole operator O ψ5ψ∂F . The resulting interaction with nuclei is velocity suppressed and therefore gives much weaker constraints that the other operators we have considered in this section, Λ/C 1/2 ψ5ψ 50 GeV. As in the case of O ψψ , the dominant constraint in indirect detection comes from annihilation directly into charged leptons, though in this case the cross section is p-wave suppressed, leading again to rather weak constraints. Constraints from the FERMI Galactic halo search dominate over that in dwarf Spheroidals by a factor of a few, owing to the larger DM velocities expected in the Milky Way halo.

Majorana Fermion Dark Matter
In Fig. 2, we show constraints on the 4 most relevant high-scale operators which give rise to low energy interactions between photons and Majorana fermion Dark Matter, χ. In contrast with the Dirac fermion case, the magnetic and electric dipole operators are forbidden, so we now include the Rayleigh interactions (∝ F µν F µν , F µνF µν ). We also include the anapole 10 GeV compared to Λ/C 1/2 χ5χ 50 GeV (corresponding to a difference in the rate of a factor of ∼ 600). This justifies our omission in the previous section of the fermion charge radius operator itself, in favour of the 4fermion operator. The two induce the same low-scale effective operator between DM and the photon, but the latter receives a logarithmic enhancement from RG evolution.
In the bottom left panel of Fig. 2, we show constraints on the Rayleigh operator O χχF F . As detailed in Appendix A (and originally in Refs. [26,28]), this interaction leads to coherent nuclear scattering, with a cross section scaling as Z 4 /Λ 6 for a nucleus of charge Z. However, the Rayleigh operator is dimension-7 and the scaling with Λ dominates (along with the loop-suppression prefactor), resulting in constraints at the level of Λ/C 1/3 χχF F 4 GeV. The cross section for annihilation into 2 photons also scales as Λ −6 , as well as being p-wave suppressed. Even so, above DM masses around 100 GeV, indirect detection constraints begin to dominate over direct searches, as the annihilation cross section is enhanced as m 4 ψ . In the bottom right panel of Fig. 2, we show constraints on the Rayleigh operator O χ5χFF . In this case, the direct detection signal is negligible [28] as the cross section is momentum-suppressed, while the annihilation cross section is s-wave, leading to much stronger constraints than in the case of O χχF F . The cut-off scale is constrained to be Λ/C 1/3 χ5χFF 100 GeV at m DM = 100 GeV, again with the cross section growing rapidly with m ψ . Note that constraints on the Rayleigh operator O χ5χF F are identical to those on O χ5χFF (shown in the bottom right panel of Fig. 2); the two share the same momentumsuppressed direct detection cross section and s-wave annihilation cross section. As we note in Sec. A.2, direct and indirect constraints on the remaining Majorana Rayleigh operator O χχFF are negligible.
We note here that the Rayleigh operators O χ5χF F and O χ5χFF are more strongly constrained than the (induced) anapole operator, due to the s-wave annihilation cross section and scaling with the DM mass. Of course, the relative importance of the different operators depends on the exact UV completion, which we explore later in Sec. 4.

Scalar Dark Matter
In Fig. 3, we show constraints on the most revelant high-scale operators which give rise to low energy interactions between photons and complex scalar Dark Matter S (top row) or real scalar Dark Matter R (bottom row).
In the top left panel of Fig. 3, we show constraints on the scalar charge radius operator, O ∂S∂SF . In direct detection, this gives rise to a short range coherently-enhanced coupling to the protons in the nucleus (up to O(1) factors, the cross section is the same as for the fermion charge radius). The resulting constraint is Λ/C 1/2 ∂S∂SF 400 GeV for complex scalar DM with a mass around 50 GeV. For indirect detection, the charge radius operator couples only to electromagnetic currents, so annihilations into photons is not possible at tree-level. This means that searches for gamma ray lines are not constraining.
In the top right panel of Fig. 3, we show constraints on the complex scalar coupling to Standard Model leptons, O S∂S , which induces the scalar charge radius operator radiatively. As we saw in previous sections (and in particular in the top row of Fig. 2), the logarithmic enhancement to the scalar charge radius operator coming from the RG evolution dominates over the contribution from the tree-level O ∂S∂SF defined at the high scale. The resulting constraint is Λ/C S∂S 2000 GeV, stronger than the constraint on O ∂S∂SF . Indirect detection constraints come from annihilation directly into charged leptons, although the cross section is p-wave suppressed, leading to much weaker constraints than those coming from direct detection.
In the bottom panels of Fig. 3, we show constraints on the dimension-6 real scalar Rayleigh operators, O RRF F and O RRFF . The resulting constraints follow a similar pattern as for the Majorana fermion Rayleigh operators (bottom row of Fig. 2). While the scalar Rayleigh operators are dimension-6 (rather than dimension-7 as in the case of the Majorana fermion), direct detection constraints on O RRBB are still very weak. They also weaken rapidly with increasing m DM , as the cross section scales as σ RRF F ∼ m −2 DM Λ −4 . For both operators O RRF F and O RRFF , the annihilation cross section is s-wave and scales as m 2 DM Λ −4 , meaning that constraints from indirect detection are similar for the two operators and dominate over direct detection in both cases.

UV-related Operators
Finally we discuss the bounds on the UV-related tree-level operators given in Eqs. (2.7) and (2.9). In particular, we highlight cases where bounds on UV-related operators are stronger than bounds on DM-photon operators with the same CP properties. In such cases, realistic UV models may be more strongly constrained than would be expected from considering only the DM-photon interactions which they give rise to.
In Fig. 4, we show bounds on dimension-7 (top row) and dimension-8 (bottom row) UV-related operators for Dirac Fermion DM ψ. In the top left panel, we show constraints on the operator O ψψ 5 . These come predominantly from annihilation into charged leptons at tree-level, with constraints from dwarf Spheroidals and the Galactic halo being comparable (Λ/C 1/2 ψψ 5 200 − 1000 GeV). Direct detection constraints arise only at loop-level, driven by the charged lepton Yukawa [50]. For coupling to e and µ, DD constraints are negligible, while for coupling to the τ constraints may extend up to O(50 GeV). We see from Fig. 4 that this operator O ψψ 5 is the most strongly constrained of the Dirac fermion UV-related operators. Even so, these constraints are much weaker than the corresponding constraints from the magnetic dipole operator O ψψF (CP-even) or the electric dipole operator O ψ5ψF (CP-odd), presented in the top row of Fig. 1. For Dirac DM then, we expect that the DM-photon EFT should capture the most relevant constraints, even including additional operators which may appear in the UV.
In Fig. 5, we show constraints on the UV-related operators for Majorana Fermion DM χ. Considering the operator O χ5χ 5 (top), we see that DD constraints are again loopsuppressed and non-zero only in the case of coupling to τ leptons. In ID, the annihilation cross section has both an s-wave contribution (proportional to m 2 ) and a p-wave contribution (proportional to m 2 χ ). For large DM masses, the FERMI constraints are similar to those from O χ5χ (shown in Fig. 2), arising from the p-wave contribution (with the limits being almost universal for e, µ and τ for m χ > 2 TeV). Instead, for lighter DM, the s-wave contribution becomes dominant. Constraints for different charged leptons then diverge, substantially strengthening in the case of coupling to the heavier τ . These ID constraints may be competitive with the DD constraints arising from the operator O χ5χ (which induces the anapole operator at low energy). In contrast to the Dirac case, the strongly-constrained electric and magnetic dipoles operators are forbidden for Majorana DM so constraints from the UV-related O χ5χ 5 may be as relevant as the strongest constraints from DM-photon interactions.
In the bottom row of Fig. 5, we show constraints on the dimension-8 UV-related operators. These couple DM (pseudo-)scalar currents to lepton (pseudo-)scalar currents and do not give rise to DD signals. For operators coupling to the scalar DM current χχ (bottom left panel), the annihilation cross section is p-wave suppressed, leading to rather weak constraints, at the level of Λ/C 1/4 χχ (5) −8 10 GeV, depending on the lepton. Instead, annihilation through the pseudo-scalar DM current χγ 5 χ (bottom right panel) is s-wave, leading to constraints on Λ up to 100 GeV for couplings to the τ . Again, such constraints may be competitive with constraints from DD and ID shown in Fig. 2, coming from the operator O χ5χ (CP-even) or the Rayleigh operators O χ5χFF (CP-even) and O χ5χF F (CPodd). This suggests that the UV-related operators may be particularly relevant in the case of Majorana DM (though typically only for coupling to the τ ).
In Fig. 6, we show constraints on UV-related operators for complex (top row) and real (bottom) scalar DM. For the operator O S∂S 5 (top left), the dominant constraints come from direct detection (which is loop-suppressed and only relevant for coupling to the τ , as in the case of fermionic DM); in this case the annihilation cross section is p-wave suppressed. For the remaining operators, the bounds arise from s-wave annihilation into leptons at tree-level, with negligible DD constraints. For the case of complex scalar DM, these limits are always weaker than limits from the (induced) scalar charge radius (top row of Fig. 3): the UV-related operators can safely be ignored and the EFT describing DM-photon interactions should capture all the relevant effects. For the real scalar, these bounds may be competitive with constraints from the Rayleigh operators O RRF F and O RRFF (bottom row of Fig. 3), at low DM mass and for coupling to the τ only. The annihilation cross section for the dimension-6 scalar operators scales as m 2 , meaning that constraints for coupling to e and µ are much weaker and can typically be neglected.

Example UV Completions
In this section we consider two simple examples of UV models that induce DM-photons interactions at one-loop by coupling DM to SM leptons. This mainly serves as a illustration of the usefulness of the EFT analysis for deriving constraints from DD and ID, which we compare to bounds obtained from a fixed-order calculation. This allows us to discuss the region of UV parameter space where the EFT is valid and identify the regime of EFT breakdown, in particular in DD. We complete this section by considering also collider constraints, where Bhabha scattering at LEP turns out to provide an important bound on the UV completion with electrons.

Fermion DM: DD and ID Constraints
We add to the SM a singlet fermion χ that is either Dirac or Majorana and a complex scalar S with hypercharge −1. The relevant part of the Lagrangian is given by where i denotes a SM lepton. In order to avoid constraints from lepton flavor-violating processes, we only consider the case of a single non-vanishing λ i . We can then always make λ i real through a phase redefinition of S, so this theory conserves CP. We take M DM M S , so χ is a DM candidate that is stable because of a Z 2 symmetry.
We now match this theory to the relevant (dominant) effective operators in Table 1, by integrating out the heavy scalar and expanding in powers of external momenta or equivalently M DM . As a result of CP invariance, we have C ψ5ψF = 0, while for the non-vanishing operators the parametric dependence on λ i , M DM , M S can be inferred from dimensional analysis, and only the numerical factors have to be calculated. Below we show only the results for the most relevant operator coefficients 8 . For the Dirac case, one obtains while in the Majorana case one finds Using these coefficients, one can derive the allowed regions in the UV parameter space M DM − M S /λ i , just as we derived bounds on the individual effective operators from Sec. 3. These are shown as solid lines in the upper (Dirac) and lower (Majorana) panels of Fig. 7, for all three SM leptons. It is instructive to compare these bounds, which have been obtained from matching to the EFT at leading order in M DM /M S and m i /M S , to the constraints obtained from a fixed order calculation in the full UV theory. This gives the result for the Dirac fermion DM-quark amplitude where p (k) is the momentum of the incoming DM fermion χ (incoming quark q = u, d), p (k ) is the momentum of the outgoing DM fermion χ (outgoing quark q) and q = p − p = k − k. Furthermore Q q denotes the electric charges of the SM quarks, and the functions C α (λ i , m i , q, M DM , M S ) can be found in Appendix C. The case of Majorana DM is given by the Dirac DM result with the replacements C σ → 0, C γ → 0, C γ5 → 2 C γ5 . From these expressions one can recover the EFT coefficients by taking the limit of small lepton masses, small DM masses, and small momentum transfer, m 2 i , M 2 DM , q 2 M 2 S . In this limit some diagrams acquire an IR divergence, which is cut-off by m 2 i or q 2 , depending on which is larger. For m 2 i > q 2 one finds at leading order in M 2 DM /M 2  Figure 7: Constraints on the UV models for Dirac (upper panels) and Majorana DM (lower panels), for couplings to e, µ and τ . In grey is shown the excluded region of parameter space where M DM > m S , which depends on λ i due to the choice of the y-axis. Regions below the brown lines are excluded by DD constraints from XENON1T, while the blue (magenta) lines denote the bounds from ID, coming from FERMI observations of dwarf Spheroidal galaxies (Galactic Halo). Constraints shown as solid lines arise from considering just the effective operators, and only depend on the combination m S /λ i . The full fixed order calculation also depends on λ i and m S separately, and we show the resulting constraints for λ i = 1/3 (dot-dashed), λ i = 1 (dotted) and λ i = 3 (dashed).
The resulting bounds on the UV parameter space from DD (XENON1T) are shown in the upper (Dirac) and lower (Majorana) panels of Fig. 7, for all three SM leptons. In contrast to the EFT calculation, these bounds depend not only on the combination λ 2 i /M 2 S but also on the ratio M 2 DM /M 2 S . We have therefore chosen to present the bounds for three cases of the UV coupling λ i = {1/3, 1, 3}, which are denoted by brown dot-dashed, dotted and dashed lines, respectively.
Before we discuss the results, we also need to consider the bound from ID due to the direct annihiliation into leptons in the UV theory. The annihilation cross-sections into leptons i are given by: for the case of Dirac/Majorana fermion, respectively, with v denoting the relative velocity and r F = M DM /M S < 1. This implies a dependence on three independent parameters M DM , M S and λ i , therefore we again present the bounds from ID (FERMI) in Fig. 7 for three cases of the UV coupling λ i = 1/3, 1, 3, which are denoted by blue dashed-dotted, dotted and dashed lines, respectively. Besides the bounds obtained from the EFT and the full calculation for both DD (brown) and ID (blue and magenta), we have indicated in Fig. 7 the (grey) regions of the parameter space that are excluded for stable DM, M DM > M S (which depend on λ i due to the choice of our y-axis). It is clear that close to the border of this region the EFT starts to break down, which is illustrated by the discrepancy between the solid and the various dashed lines. In fact, close to the physical boundary the bound in the full theory starts to grow, which is due to the fact that the amplitude in Eq. (4.4) acquires an infrared divergence in this limit that is regulated by the charged lepton mass. Indeed in this region the DD bound is dominated by the dipole coefficient C σ in Eq. (4.4), which for The solid curves in the EFT calculation are therefore an excellent approximation of the full calculation, and allow us to infer an approximate bound on the UV model without doing the full calculation. It turns out that a tree-level calculation is enough to estimate the correct bound (using the results in Sec. 3) in almost the entire parameter space, as the constraints from DD and ID are largely dominated by the tree-level operators C ψψ (Dirac) and C χ5χ (Majorana). In the case of ID, this is true in all points of the parameter space, since they are dominated by the tree-level annihilation into leptons (and annihilation into photons at one-loop gives subleading bounds). For Majorana DM, ID gives the most stringent bound only for coupling with the τ , where the cross section is s-wave and proportional to m 2 τ . In this case, DD bounds are very weak, because of an approximate cancellation between the contributions of the anapole and the four-fermion operator (via RG running), see Eq. (4.5). For Dirac DM, it is DD that sets the strongest bounds (except for very light DM masses), which at DM masses below roughly 300 GeV arise from the tree-level operator C ψψ . Above this mass they are surpassed by the constraints on the magnetic dipole C ψψF , since this coefficient grows with M DM , cf. Eq. (4.2). This behaviour is somewhat surprising; typically the dipole operator is assumed to give the strongest constraints, while in this case we must include the higher-dimensional C ψψ to fully capture the behaviour of the DM-photon interactions.

Scalar DM: DD and ID Constraints
We now consider a second UV model, in this case for scalar DM. The logic follows closely that of the Fermion DM Lepton portal discussed in Sec. 4.1. However, we include the full discussion again for completeness.
We add to the SM a singlet scalar s that is either complex or real, and a Dirac fermion F with hypercharge 1. The relevant part of the Lagrangian is given by where i denotes a SM lepton. In order to avoid constraints from lepton flavor-violating processes, we only consider the case of a single non-vanishing λ i . We take M DM M F , so that s is a DM candidate stabilized by a Z 2 symmetry.
We now match this theory to the dominant effective operators in Table 1, integrating out the heavy fermion and expanding in powers of external momenta or equivalently M DM . As a result of CP invariance, we have C SSFF = C SS 5∂ = 0, while for the non-vanishing operators the parametric dependence on λ i , M DM , M F can be immediately inferred from dimensional analysis, and only the numerical factors have to be calculated. For the complex case, one obtains while in the real scalar case one obtains Using these coefficients, one can derive the allowed regions in the UV parameter space M DM − M F /λ i , just as we derived bounds on the individual effective operators from Sec. 3. These are shown as solid lines in the upper (complex scalar) and lower (real scalar) panels of Fig. 8. It is instructive to compare these bounds, which are obtained from matching to the EFT at leading order in M DM /M F and m i /M F , to the constraints obtained from an fixed order calculation in the full UV theory. The result for the DM-quark amplitude is (for a complex scalar) where p (k) is the momentum of the incoming DM scalar s (incoming quark q = u, d), p (k ) is the momentum of the outgoing DM scalar (outgoing quark q) and q = p − p = k − k. Furthermore Q q denotes the electric charges of the SM quarks, and the function C S (λ i , m i , q, M DM , M F ) can be found in Appendix C. In the case of real scalar DM, the above amplitude vanishes, C S → 0. From this expression one can recover the EFT coefficient C ∂S∂SF by taking the limit of small leptons masses, small DM masses, and small momentum transfer, m 2 i , M 2 DM , q 2 M 2 F . In this limit some diagrams acquire an IR divergence, which is cut-off by m 2 i or q 2 , depending on which is larger. For m 2 i > q 2 one finds at leading order in M 2 DM /M 2 Figure 8: Same as Fig. 7, but for complex scalar (upper panel) and real scalar DM (lower panel). We also show the limit from gamma ray line searches with FERMI (grey diamonds) and H.E.S.S. (red circles). Note that the bounds from line searches do not strictly apply whenever the continuum limits dominate.
Again this expression yield the EFT result in Eq. (4.2), apart from the logarithm that is reproduced in the EFT by running down the tree-level operator, recovering Eq. (2.6). The resulting bounds on the UV parameter space from DD (XENON1T) are shown in the left (complex scalar) and right (real scalar) panel of Fig. 8. In contrast to the EFT calculation the bounds do not depend only on the combination λ 2 i /M 2 F but also on the ratio M 2 DM /M 2 F , therefore we have chosen to present the bounds for three cases of the UV coupling λ i = {1/3, 1, 3}, which are denoted by brown dot-dashed, dotted and dashed lines, respectively.
Finally, we also need to consider the bound from ID due to direct annihiliation into leptons in the UV theory. The annihilation cross-sections into leptons i are given by: for the case of real/complex scalar, where v is the relative velocity and r = M DM /M F < 1.
Besides the bounds obtained from the EFT and the full calculation for both DD (brown) and ID (blue), we have indicated in Fig. 8 the regions of the parameter space that are excluded for stable DM, M DM > M F (which depend on λ i due to the choice of our y-axis). It is clear that close to the border of this region the EFT starts to break down, which is illustrated by the discrepancy between the solid and the various dashed lines. In fact, close to the physical boundary the bound in the full theory starts to grow, which is due to the fact that the amplitude in Eq. (4.12) acquires an infrared divergence in this limit that is regulated by the charged lepton mass. Indeed in this region, where As in the fermion DM case, the solid curves in the EFT calculation are an excellent approximation of the full calculation. In the complex scalar case, the bounds are dominated by DD, with ID playing no role for any DM mass. In contrast to the fermion DM case, however, these DD bounds come predominantly from a single tree-level operator O S∂S as expected from the EFT analysis of Sec. 3.3. For real scalar DM, the Rayleigh operator O RRF F appears with a coefficient of 1/3, suppressing the (already weak) DD bounds and making them negligible. For coupling to e and µ, the Rayleigh operator still gives rise to substantial annihilation into gamma ray lines, so line searches dominate. Instead, for coupling to the τ , the most relevant bounds at low DM mass come from FERMI diffuse searches (blue and magenta) in which case the ss → cross section scales with m 2 /M 2 F . The H.E.S.S. line search may still be relevant for strongly coupled theories at high DM mass, though we note that the H.E.S.S. analysis makes optimistic assumptions about the DM density profile in the Milky Way, as discussed at the start of Sec. 3.

Fermion and Scalar DM: Collider Constraints
Finally we discuss also the constraints from colliders on the UV models in the previous sections. We consider constraints from direct searches at the LHC and LEP and indirect constraints from Bhabha scattering at LEP-2. The bounds are summarized in Fig. 9 and discussed in the following. Note that in other UV scenarios, such as Z -mediated UV completions, the constraints from collider searches can be significantly stronger [33,50].

Direct Searches at LHC and LEP
As discussed in e.g. Ref. [30], the dominant collider constraints on these UV completions come from scenarios where the charged mediator is pair produced, each sub-sequently decaying into a charged lepton and a DM particle. The signature is therefore + − plus missing energy.
Fermion DM: In the case of fermionic DM (scalar mediators), this signature has been considered in searches for sleptons in the MSSM at LHC energies of 8 TeV [64][65][66][67] and 13 TeV [68][69][70]. We use the 95% CL bounds from the lower right panel of Fig. 5 (Fig. 6) of Ref. [69] for selectrons (smuons) decaying with a 100% branching fraction into right-handed electrons (muons) and a neutralino. We show the resulting bounds on M S as green regions in the upper panel of Fig. 9 (note that they do not depend on the coupling λ i as long as the decay is prompt). In the case of staus, LHC searches, e.g. in Ref. [71], are not yet sufficiently sensitive to derive any constraints.
Direct searches for final state leptons with missing energy have also been carried out at LEP, yielding constraints on heavy charged scalars of about 80-90 GeV. For the 95% CL limits we use constraints on sleptons provided by the PDG in Ref. [72]. For different values of the slepton-DM mass difference ∆m, these bounds are mẽ > 73 GeV for any ∆m, and mμ > 88 GeV, mτ > 79 GeV for ∆m > 15 GeV, and shown as blues regions in the upper panel of Fig. 9.
Scalar DM: In the case of scalar DM (fermionic mediators), one can recast the LHC slepton searches. In the case of couplings to right-handed muons, this has been done in e.g. Ref. [73] (see Fig. 5 therein). The case of couplings to right-handed electrons should be similar (as one can see from the slepton case), so for simplicity we use the bounds from the muon case also for electrons, shown as green regions in the lower panel of Fig. 9. For fermions coupling to right-handed taus, one can recast the stau searches at CMS [71], which yield upper bounds on the stau pair production cross section. These analyses can be potentially relevant, because the production cross section of a heavy fermion is larger than that of heavy scalars (with the same mass and gauge quantum numbers). For the precise values of the cross section we use Fig. 2.1 in Ref. [74], which is about a factor 10 larger than stau pair production, and allows us to derive bounds on this scenario, shown in green in the lower right panel of Fig. 9.
In this case LEP provides bounds for any mass splitting and all three flavors. We use the pure higgsino bound yielding m F > 92.4 GeV at 95% CL [75], shown as blue regions in the lower panel of Fig. 9.

LEP Indirect Searches
As in section 3, we use the EFT analysis of Refs. [54,55] to constrain the UV model with couplings to electrons. Integrating out the heavy fermion and scalar and one-loop, one can match to the low-energy 4-electron operator with a coefficient given in terms of UV parameters as Here f D,M,S denote the loop functions for the case of Dirac/Majorana fermion and complex scalar DM while for real scalar DM c RR = 0. The EFT coefficient c RR is constrained by Bhabha scattering at LEP-2 according to and we use this bound at 95% CL to derive bounds on the UV parameters, shown as red contours in Fig. 9 for different values of λ e = {1, 3}, distinguishing Dirac (dark red) and Majorana DM (light red). It is instructive to compare this result with the EFT estimate in Section 3, see Eq. (3.7). Using the tree-level matching conditions in Eqs. (4.2) and (4.10), one can check that the EFT result exactly reproduces the results of the full calculation in the EFT limit M 2 F,S → 0 for all four cases, upon the identification Λ = M S,F in the EFT.

Comparison with DD and ID Constraints
In this section we compare the bounds from colliders in Fig. 9 to those from direct and indirect DM searches in Figs. 7 and 8.
Dirac DM: For Dirac fermions, bounds from DD prevail for λ i = 1 for all SM leptons and for M DM > 10 GeV (see upper panels in Fig. 7). This holds even more for λ µ,τ > 1 in the cases of µ and τ , while for e the LEP-2 indirect searches give stronger constraints than DD for λ e 4. Instead for λ i < 1 colliders typically provide stronger constraints than DD. For e, µ LHC limits prevail for λ e,µ 0.2, while LEP starts to become more relevant than DD for all leptons at around λ i 0.1. Note that close to the degenerate limit DD bounds always dominate, since in this case collider searches loose sensitivity while DD cross sections are enhanced, cf. Eq. (4.8).
Majorana DM: For Majorana fermions, bounds from colliders prevail for λ i < 1 for all SM leptons (see lower panels in Fig. 7). ID bounds are only important for taus with λ τ ≥ 1, while DD bounds only play a role close to the degenerate limit (as discussed above). Note in particular that in the electron case, Bhabha scattering provides stringent bounds (in the TeV regime) for λ e > 1.
Complex Scalar DM: The complex scalar case is very similar to the Dirac fermion case, apart from stronger limits from colliders. In particular, bounds from DD prevail for λ i = 1 for all SM leptons and for M DM > 10 GeV (see upper panels in Fig. 8). This holds even more for λ µ,τ > 1 in the cases of µ and τ , while for e the LEP-2 indirect searches give stronger constraints than DD for λ e 4. Instead for λ i < 1 colliders typically provide stronger constraints than DD. LHC limits prevail for λ e,µ 0.5 and λ τ 0.4 , while LEP starts to become more relevant than DD for all leptons at around λ i 0.1. Close to the degenerate limit DD always dominate, since in this case collider searches loose sensitivity while DD cross sections are enhanced, cf. Eq. (4.8) (although in contrast to fermion DM LEP completely closes the kinematic gap).
Real Scalar DM: The real scalar case is very similar to the Majorana fermion case, apart from stronger limits from colliders (especially the presence of LHC bounds in the τ case), and no bounds from DD and Bhabha scattering, In particular, bounds from colliders prevail for λ i < 1 for all SM leptons (see lower panel in Fig. 8). ID bounds are only important for taus with λ τ 2. In the degenerate limit, LHC (but not LEP) bounds can be avoided, and for strong coupling λ e,µ ≈ 5 the constraints from annihilation into gamma ray lines can play an important role.

Summary and Conclusions
In this work, we have studied the scenario where DM couples to quarks only via photon exchange. To this end, we have considered the most general effective interactions of fermion and scalar Dark Matter with the photon, distinguishing between Dirac/Majorana and complex/real scalar DM, see Eqs. (2.1)-(2.4). We have systematically classified these operators according to their CP properties and dimension, as we summarize in Table 1. Since DM is neutral under electromagnetism, these operators can only arise at loop level at the cutoff scale. This implies that for direct detection (DD) and indirect detection (ID), there are other relevant operators, which arise at tree-level at the cutoff scale and mix into DM-photon operators via RG evolution. According to our philosophy, these can only be four-fermion operators coupling DM to vector currents of SM leptons, which we therefore include in the classification and list in Eq. (2.5). Moreover, such operators are usually accompanied by other DM-lepton currents due to SM gauge invariance. We also include these operators in the classification, referring to them as "UV-related operators", see Eqs. (2.7)-(2.9).
We have then studied the model-independent constraints from direct and indirect searches in Section 3, considering one operator at a time. The results are shown for Dirac DM in Fig. 1, for Majorana DM in Fig. 2 and for complex and real scalar DM in Fig. 3. For the operators coupling directly to photons, like the dipole operators, we update existing bounds in the literature, in particular showing bounds from the latest Xenon1T (2018) results. We find that tree-level operators are strongly constrained by both DD and ID, depending on the Lorentz structure. In the case of DD the resulting bounds provide stringent constraints, due to large RG mixing from the SM vector current. For example, the 4-fermion operator O ψψ = ψγ µ ψ · γ µ mixes into the charge radius operator O ψψ∂F = ψγ µ ψ · ∂ ν F µν and bounds from the former interaction are stronger than those from the latter. In the Dirac DM case, these bounds (on dimension-6 operators) can even compete with those on the electromagnetic dipole (of dimension-5), which has the same CP quantum numbers. We also note the bounds on the UV-related operators can be important in the case of Majorana or real scalar DM coupling to tau leptons (see Figs. 5 and 6), which in general are more weakly constrained. Finally, we also derive bounds on four-fermion operators involving DM and electrons from Bhabha scattering at LEP-2, which can be quite stringent due to the fixed sign of this new contribution relative to the SM.
In Section 4 we have studied two examples of explicit UV completions, in which DM couples to SM leptons and a new heavy fermion or scalar. The main purpose of this section is to demonstrate the usefulness of the EFT analysis, which allows us to gain further insight into the resulting constraints from DD and ID, shown in Figs. 7 and 8. Here we compare the bounds obtained in the EFT prediction to those from the full fixed-order calculation, finding excellent agreement except in the degenerate limit for DD. In particular, in most of the parameter space the dominant bounds come from 4-fermion operators (as expected from the results of Section 3), which can be obtained by a simple tree-level calculation. Finally, we have discussed the constraints from colliders (see Fig. 9) on these UV completions, which can be more important that DD and ID for Majorana fermion and real scalar DM. In particular, we also studied new constraints from LEP-2 on the electron portal DM scenario, which give important constraints for sizable couplings.
In conclusion, we have provided model-independent constraints on a wide class of effective operators coupling DM to photons or leptons. This compendium will allow the reader to quickly assess the constraints on any DM model that interacts with quarks mainly via photon exchange. We have illustrated this approach in two explicit UV models. Updating the bounds on these UV models, we find that DM-photon interactions are typically constrained by multiple complementary probes. The EFT analysis provides a simple way of taking into account all relevant constraints and shedding light on DM-photon interactions without exhaustive calculations.

Acknowledgments
We thank Eugenio Del Nobile and Andrew Cheek for helpful discussions concerning the anapole interaction, Adam Falkowski for providing details on Refs. [54,55], and Lorenzo Calibbi and Emmanuel Stamou for useful discussions on collider constraints and technical aspects of the loop calculations. BJK acknowledges funding from the Netherlands Organization for Scientific Research (NWO) through the VIDI research program "Probing the Genesis of Dark Matter" (680-47-532).

A Direct Detection
In this appendix we collect the relevant expressions that have been used to derive the bounds from direct detection. We begin by reviewing the non-relativistic effective field theory (NREFT) formalism, which gives the differential scattering cross section in terms of non-relativistic operator coefficients. Then these coefficients are matched to the relativistic operators below the weak scale, which allows us to provide the contribution of a single operator to the differential scattering cross section. Finally we summarize the approximate bounds from XENON1T.

A.1 NREFT Formalism
Here we list the NREFT operators [76][77][78][79][80][81] we consider in this work, with normalisations chosen to match that of Ref. [44] (this is the same normalisation as in the public tool NRopsDD tools, which we have used to obtain the bounds on relativistic operators in Section 3). In the context of DM-photon interactions, the most relevant operators 10 describing DM interactions with nucleons N = p, n are: (A.1) Here q denotes the momentum transfer, v ⊥ the transverse WIMP-nucleon velocity, S DM the DM spin and S N the nucleon spin. The transverse velocity is given by is the reduced mass of the WIMP-nucleon system and m N is the nucleon mass. Note that O NR 1 is the standard spin-independent (SI) operator, which is coherently enhanced and therefore typically dominates. The spin-dependent (SD) operator corresponds to O NR 4 . The DM-nucleon matrix element can then be expressed as a sum over operators, where the coefficients c N =p,n i can be calculated from a given relativistic Lagrangian (see e.g. Ref. [81]). The DM-nucleus matrix element is then obtained by summing over nucleons in a given target nucleus T . The differential scattering cross section can then be written in terms of the appropriate nuclear response functions F , which encode the nuclear structure: The nuclear recoil energy E R and recoil momentum are related by q 2 = 2m T E R . The nuclear response functions F i,j can be written in terms of nuclear form factors F X , where X = M, Σ , Σ , ∆,Φ , Φ , Σ ∆ correspond to sums over different nucleon properties in the target nucleus T . We employ the following relations (suppressing the N = p, n indices and q 2 dependence in all but the first expression): Here, C(j DM ) = 4j DM (j DM + 1)/3, where j DM is the DM spin. The transverse velocity appearing here is the WIMP-nucleus velocity: where µ DM,T = m DM m T /(m DM + m T ) is the reduced mass of the WIMP-(target nucleus) system with m T denoting the target nucleus mass. When the operator O NR i appears with additional powers of q n , the corresponding response function is q 2n F i,i .
The response F M is the standard SI form factor which is coherently enhanced and therefore scales as A 2 (for equal couplings to protons and neutrons). The standard SD form factor is a combination of F Σ and F Σ . Nuclear form factors have been calculated in Refs. [76,77,82]. For Xenon, we use the form factors listed in Appendix A.3 of Ref. [76], apart from the spin-dependent form factors, which we take from Ref. [83] which includes two-body currents.
Finally, the differential recoil rate with a given target T is obtained from the differential cross section as: where ρ χ is the DM density, f (v) the DM velocity distribution and v min = m T E R /(2µ 2 χT ). Further details can be found in e.g. Ref. [84].

A.2 Non-relativistic Operator Matching and Cross Sections
Full details of the matching of the photon interactions with quark-level and nucleon-level interactions can be found in a number of references [44,78,85,86]. Here we highlight a subtlety in the matching which is not always considered.
We can use the equations of motions in order to rewrite operators involving ∂ ν F µν in terms of currents of Standard Model fermions f . We are considering interactions with nucleons, so we restrict ourself to interactions with quark currents q: Note that this expression describes the interaction with free fermions, while the nucleons we are interested in consist of bounds states of quarks. The expectation value of the quark vector current inside the nucleon can be written as (see e.g. Eq. (22) of Ref. [86]): where u N are nucleon spinors. The two form factors F (q,N ) 1 (q 2 ) ('Dirac') and F (q,N ) 2 (q 2 ) ('Pauli') encode the internal nucleon structure (as a function of the momentum transfer q 2 ). The first describes the contribution of quarks q to the charge of nucleon N while the second describes the contribution to the nucleon magnetic moment. Galactic Dark Matter is not energetic enough to probe the internal nucleon structure, so we can safely set q 2 → 0, in which case we obtain (see Appendix A of Ref. Now, using Eq. (A.9), we see that the quark electromagnetic current embedded in the nucleon is equivalent to: where a p = 1.793 and a n = −1.913 are calculated from the Pauli form factors given above (these are in fact the anomalous magnetic moments of the nucleons). Typically when coupling to the DM vector current, the first term in Eq. (A.11) leads to the standard spin-independent interaction, in which case the second term is sub-dominant and may be neglected [44]. However, in general (and in particular in the case of the anapole interaction O χ5χ∂F ), both terms must be included. The relativistic interactions at the DM-nucleon level can then be matched onto the NREFT operators described in Sec. A.1 using standard 'dictionaries' (e.g. Ref. [81]). Below we list the non-relativistic matrix elements for DM scattering with nucleons N = p, n (as in Eq. (A.3)) which are induced by the EMSM χ operators after this matching has been performed: 11 The corresponding expression for neutrons is obtained by isospin symmetry: n ↔ p, u ↔ d, s ↔ s.
Here Q p = 1, Q n = 0 are the nucleon electric charges and g p = 5.59, g n = −3.83 are the nucleon g-factors [44]. The factor C denotes the dimensionless numerical coefficient of each of the operators appearing in the relativistic Lagrangian. The precise definitions of the operators and their numerical coefficients (loop factors and O(1) factors) are given in Sec. 2.1. As an example, in the case of the magnetic dipole interaction, we have C ≡ e C ψψF /(32π 2 ). The operators in Sec. 2.2 may induce some of the interactions appearing in Eq. (A.12), with the appropriate coefficients given in Eq. (2.6). Using Eq. (A.4), a single EMSM χ operator leads to a contribution to the differential cross-section for a target nucleus T given by: where α = e 2 /(4π). Following Ref. [26], the differential cross section for the Majorana and Dirac Rayleigh operators O χχF F and O ψψF F is given by: where the nuclear coherence scale is Q 0 = √ 6R N and the nuclear radius is given by [1] R N = (0.3 + 0.91A 1/3 ) fm .
(A. 15) We normalize the form factor F associated with the Rayleigh operator such that F(0) = 1. An explicit expression is given in Appendix A of Ref. [28] (which corrects a number of mistakes in the derivation of Ref. [26]). As detailed in Ref. [28,87], the Rayleigh operator mixes into the operators χχqq and χχG µν G µν , which in turn lead to standard SI scattering with nucleons. These different operators may interfere and we include this effect in our analysis 12 . However, the impact on the limits turns out to be small. The energy scale Λ probed by direct detection experiments is rather low, so the mixing effects are not substantial.
Reference [88] pointed out that for Rayleigh scattering, the effects of 2-nucleon scattering should not be ignored in calculating the direct detection rate. While O(1) corrections are expected, the exact form of the proton-proton form factor (which accounts for protonproton correlations in the nucleus) is not known, so we neglect these effects here.
The remaining fermion Rayleigh operators give a negligible contribution to direct detection cross sections. The operator O χχFF gives a vanishing DM-nucleon matrix element, while O χ5χF F and O χχFF give rise to momentum-suppressed scattering [28]. Given the already weak constraints on Λ which arise from the unsuppressed operator O χχF F , we therefore neglect direct detection limits coming from the remaining operators. Note that the operators O χ5χF F and O χχFF are still constrained by indirect detection, where they give rise to an s-wave annihilation cross section.
For scalar DM, the Rayleigh operators O RRF F and O SSF F lead to a similar cross section as in the fermion case: The other scalar Rayleigh operators O RRFF and O SSFF gives a vanishing contribution to the direct detection cross section (due to the anti-symmetry properties ofF µν ) [28].

A.3 Approximate XENON1T Bounds
The XENON1T [37,89] detector is a liquid Xenon time projection chamber, operating with a total Xenon mass of 3.2t and a fiducial target mass of 1.30t. Table I of Ref. [37] reports the number of observed events and the number of expected background events for a number of different cuts. We use here the 0.9t 'reference mass' and events only in the 'reference' signal region, between the median and −2σ quantile in (cS2 b , cS1) space. We therefore calculate the number of expected WIMP signal events for a total exposure of 278.8 days × 0.9t, corrected by an overall factor of 0.475 to account for the nuclear recoil acceptance of the 'reference' region in (cS2 b , cS1). The nuclear recoil detection efficiency is taken from Fig. 1 of Ref. [37]. Approximate limit (0.9t) Figure 10: Approximate XENON1T bound (dashed blue) calculated from a single-bin Poisson upper limit. We show for comparison the median expected sensitivity (dashed black) and observed limit (solid black) of the XENON1T 2018 exposure [37]. Shaded bands show the 1σ and 2σ bands on the sensitivity. Our limit matches the official limit within a factor of a few for all DM masses above 10 GeV.
For this exposure, the number of observed events is N obs = 2, while the best fit number of expected background events is N BG = 1.62. We use this to set a single-bin Poisson upper limit on the number of signal events N 90% sig . Neglecting background uncertainties, this limit is calculated by solving [90]: k N obs +1 P (k|N BG + N 90% sig ) = 90% , (A. 17) where P (k|N ) is the Poisson probability of observing k events, given N expected events. This limit on the number of signal events is then converted into a limit on the relevant operator coupling. We show in Fig. 10 the resulting approximate limit (dashed blue) for the standard spin-independent interaction, as well as the median expected sensitivity and observed limit as reported by the XENON1T collaboration (dashed and solid black, respectively). Our approximate limit is within a factor of a few of the official XENON1T limit for DM masses heavier than 10 GeV and is consistent with other approximate limits presented in the literature [91]. The XENON1T 2018 exposure observed an upward fluctuation at high recoil energies, leading the observed limit to be weaker than the expected sensitivity at high DM mass. Our limit matches the official limit closely at high DM mass but deviates for low masses, as we do not include any recoil energy information about the individual events. A more detailed treatment would require more information about the background and signal distributions in (cS2 b , cS1) and is beyond the scope of this work.