Light dark states with electromagnetic form factors

New particles $\chi$ that are electrically neutral but couple to the electromagnetic current via higher-dimensional operators and that are sufficiently light, at or below the GeV-mass scale, can be produced in pairs in a number of dedicated high-intensity experiments. In this work we consider the production of $\chi$ through magnetic- and electric-dipole moments as well as through anapole moment and charge radius interactions in electron beams. We derive new constraints from BaBar, NA64 and mQ and forecast the future sensitivity on the existence of such states, from Belle-II, LDMX and BDX. We present for the first time a detailed treatment of the off-shell production of photons in electron beams with subsequent decay into a $\chi\bar\chi$ pair in a 2-to-4 process. These direct limits are then compared to the effects on SM precision observables, as well as to bounds from flavor physics and high energy colliders. Finally, we consider the possibility that $\chi$ is dark matter and study ensuing astrophysical and cosmological constraints. We find that a combination of all considered probes rule out $\chi$ particles with mass-dimension five and six photon interactions as dark matter when assuming a standard freeze-out abundance.


I. INTRODUCTION
Through a long history of astronomical and cosmological observations that date back into the first half of the previous century, it has now firmly been established that Newton's laws -when applied to the observed distribution of luminous matter -fail on galactic (kpc) length scales and beyond. Whereas there is an ongoing debate whether the gravitational force on small acceleration scale requires modification or if the Standard Model (SM) of particle physics is incomplete, predictions that are based on the existence of a new, neutral, and cold matter component stand unchallenged in explaining the precision observations of the cosmic microwave background (CMB). Indeed, the success of a model in which the Universe is filled with baryons, photons, neutrinos, dark matter (DM), and dark energy is so distinct that it has become a cornerstone of modern physics, termed the "standard cosmological model" or ΛCDM [1,2]. This new paradigm that purports the existence of a new form of matter, the dark matter, has become the center of much experimental and theoretical activity [3].
Even if gravity as it stands today will pass the test of time, new forces are likely at play in a dark sector of particles. While the electromagnetic interaction and the photon as its carrier stand out as the most important messengers in astronomy, dark matter -as the name suggests -must be to large degree electrically neutral and hence be non-luminous. However, even if DM carries no charge, it may still be coupled to the photon e.g. through various moments such as magnetic-or electric-dipole moments (MDM or EDM), through an * xiaoyong.chu@oeaw.ac.at † josef.pradler@oeaw.ac.at ‡ lukas.semmelrock@oeaw.ac.at anapole moment (AM) or through a charge radius interaction (CR). Indeed, "how dark" neutral DM needs to be in terms of its coupling to photons is a quantifiable question that has been addressed previously in [4][5][6]. In this work we will revisit those ideas in light of the much recent interest in sub-GeV dark sector searches, notably at the intensity frontier [7,8]. This interest not only spans DM detection, but more generally aims to probe light new particles and in this spirit we do not require that χ is DM per se but only long-lived on time-scales pertinent to terrestrial experiments.
In this paper we shall consider Dirac fermions χ that interact with an external electromagnetic field or current through MDM, EDM, AM and CR. In the non-relativistic limit, the corresponding respective Hamiltonian operators for the particle χ are H MDM = −µ χ (B·σ χ ), H EDM = −d χ (E·σ χ ), H AM = −a χ (J ·σ χ ) and H CR = −b χ (∇·E). We will refer to µ χ , d χ , and a χ as the magnetic dipole-, electric dipole-, and anapole moment, respectively, and to b χ as the charge radius. These operators differ most notably in the symmetry properties with respect to space reflection (P), charge conjugation (C), and time-reversal (T). Since the magnetic field B and the spin S χ = σ χ /2 are axial vectors, an MDM interaction is P-and T-even. Hence CP is a good symmetry, and, as is inherent in the field-theoretical description of particles, any state χ that has electromagnetic charge Q χ will automatically carry a magnetic moment µ = gQ χ /(2m χ ) S χ with g = 2 up to anomalous contributions; m χ is the mass of χ. In turn, any EDM interaction with the electric field E maximally violates P and T and is a new source of flavor-diagonal CP violation -one that has not been observed in any SM particle yet. Hence, searches for EDMs are considered essentially "background free" probes of new physics [9] and any detection of a dark particle through its EDM would be particularly striking. If χ has a non-vanishing anapole moment, it will interact with external electromagnetic currents J ; in that case χ may even be a Majorana parti-cle. Unlike MDM and EDM, AM does not correspond to a multipolar distribution of charge. It was first proposed by Zel'dovich [10] as a P-violating (but CP conserving) electromagnetic interaction, and has since been discovered in atomic nuclei [11]. Finally, the CR interaction appears at the second order in the expansion of a charge form factor, and as such it is a quantity e.g. well present in composite particles of the SM, and it is a searched for property of neutrinos [12]. Its transformation property with respect to discrete symmetries is one of a scalar.
The effective interactions of neutral χ with the photon field may arise in a variety of UV-complete descriptions. An imminent possibility is that χ is a composite particle so that a dipole moment arises through the particle's internal structure, such as, e.g., in technicolor theories [13][14][15]. Another possibility are frameworks with an extended set of particles such that electromagnetic moments are produced perturbatively, by the loops of charged states [16]. Here, a sizable MDM can be generated e.g. through axial or vector type Yukawa interactions y A,V with a scalar and a fermion at a common mass-scale M in the loop, so that parametrically µ χ ∼ Q|y A,V | 2 /M where Q is the electric charge of the mediators. The appearance of an EDM would then of course be tightly connected to CP violation in the new physics sector and d χ ∼ Q Im[y V y * A ]/M . Finally, AM and CR interactions may also arise radiatively; in the above example, a χ , b χ ∼ Q|y A,V | 2 /M 2 , possibly be enhanced through a DM-mediator mass-degeneracy [17]. Indeed, a significant amount of attention has been directed to the phenomenology of those interactions, see also [18][19][20][21][22] besides the above mentioned works and references therein. In this work, we shall remain agnostic about the origin of the MDM, EDM, AM and CR of χ. While any embedding into a concrete setup likely implies further constraints on the existence of the effective values of µ χ , d χ , a χ and b χ , this way we remain modelindependent.
Much of the effort on the intensity frontier goes into the search of light, sub-GeV mass dark sector states. New particles can be produced in electron or positron beams on fixed targets, at e + e − colliders, in hadron beams from protons on fixed targets or in proton-proton collisions at LHC. Here, much attention has been devoted to the onshell or the s-channel resonant production of particles that either decay to electrons or muons inside the sensitive detector volume or escape or decay invisibly. A prominent example is the search for dark photons [23][24][25][26] in fixed target experiments, e.g. at APEX [27] or NA48/2 [28] and at e + e − colliders, e.g. at BaBar [29] just to name a few; for a recent comprehensive compilation see [30].
In this work we study the pair-production of χχ through the MDM, EDM, AM and CR interaction with photons in experiments involving electrons as projectiles. For dark sector experiments, there are three principal detection methods: 1. search for missing momentum in e + e − colliders such as BaBar [31], Belle [32], Belle-II [33] and BESIII [34], 2. search for missing energy in e − fixed-target experiments such as NA64 [35] and the proposed LDMX detector [36], 1 3. direct search for χN or χe − scattering of χ particles produced in e − fixed-target experiments such as the previous mQ [40] detector and E137 at SLAC [41], as well as the proposed beam dump experiments at JLab (BDX) [42] or at the MESA accelerator [43].
We consider all three detection methods in turn, work out current constraints and estimate future sensitivity. As the production in fixed targets proceeds via the emission of off-shell photons, we also present a detailed exposition of the underlying 2-to-4 scattering processes, thereby filling a gap in the existing literature on BSM intensity frontier searches. We then contrast these direct limits with predictions and constraints from SM precision observables and from limits originating from flavor physics. The experimental searches for dark states on the intensity frontier have in no little part been motivated by the observation that light new physics can induce a shift in the anomalous magnetic moment of the muon, see e.g. [44][45][46][47] among other works, and hence explain a long-standing discrepancy between the SM prediction and experiment [48]. We investigate this possibility in the light of dark-sector electromagnetic form factors and finally, we return to the possibility that χ is cosmologically long-lived and may constitute dark matter. In this case, further restrictions on the parameter space of χ from cosmological and astrophysical observations take hold -such as energy injection to the SM sector, altering the CMB or influencing structure formation, as well as a local dark matter abundance that could lead to direct and indirect detection signals.
The paper is organized as follows. In Sec. II we present the relativistic versions of the interaction operators for MDM, EDM, AM and CR. In Sec. III we present our cross-section calculations on the production of χχ in e + e − collisions as well as e − -scattering on a nuclear target with details relegated to Appendix A. In Sec. IV we then specialize to the various experiments, work out constraints and sensitivity projections and discuss other complementary probes. In Sec. V we discuss indirect probes of χ-particles and in Sec. VI we comment on the possibility that χ is DM and highlight the additional restrictions on a successful model of DM. In Sec. VII we present our conclusions and provide an outlook to further going works.

II. ELECTROMAGNETIC FORM FACTOR INTERACTIONS
We consider the following interaction terms in a Lagrangian of a Dirac fermion χ that interacts with the photon gauge field A µ or its field strength tensor F µν , magnetic dipole (MDM): electric dipole (EDM): anapole moment (AM): − a χχ γ µ γ 5 χ∂ ν F µν , (1d) charge radius (CR): Here e > 0 is the electric charge, 2 µ χ and d χ are the dimensionful coefficients of the mass-dimension 5 dipole interactions, and a χ and b χ are the dimensionful coefficients of the mass-dimension 6 anapole moment and charge radius interaction; as usual, In what follows we shall measure the EDM and MDM moments in units of the Bohr magneton, µ B ≡ e/(2m e ) = 1.93 × 10 −11 e cm; m e is the mass of the electron. All coupling strengths are real.
The interactions in (1) are assembled in the matrix element of the effective electromagnetic current of χ, where p i,f and q = p i − p f are four-momenta. For a neutral particle χ the vertex function reads, 2 The sign convention in (1a) is such that χ (χ) carries a fraction of the electron (positron) charge.
We regard the various moments as being generated at an energy scale that is well above the typical center-ofmass (CM) energies of the high-intensity beams. The coefficients in (1) are hence the static limits of the form factors above, A millicharged particle χ will have an additional form factor Q(q 2 )γ µ with Q(0) = e. The next term in the q 2 -expansion of Q(q 2 ) is then the charge radius interaction and the equivalence with V (0) above can be seen when the Dirac equation is applied for on-shell χ, u(p f ) q 2 γ µ − q µ / q u(p i ) =ū(p f )q 2 γ µ u(p i ). As should be clear from above, an electric charge of χ is not mandatory for χ to possess CR interactions, and we shall consider both interactions separately below. We note in passing that the Dirac and anapole form factors V (q 2 ) and A(q 2 ) are related to the so-called vector and axial vector charge radii as r 2 V = −6V (0) and Finally, in addition to the operators listed above, at mass-dimension-7 there are the Rayleigh (or susceptibility) operators where scalarχχ or pseudoscalarχγ 5 χ bilinears multiply either the CP-even or CP-odd squared field strength tensors F µν F µν and F µνF µν . Such terms imply two-photon interactions with the χ-field and processes with a single photon in the final state proceed through a loop. In the processes arising from (1) that we shall consider in the paper, the latter interactions are then not only suppressed by the higher dimensionality of the coupling, but also by a loop-factor. As the constraints on dim-6 operators already turn out to be weak, we shall not consider Rayleigh operators any further in this work.

III. PAIR-PRODUCTION IN ELECTRON BEAMS
In this section we present the main expressions that are used to set bounds on new particles χ interacting through electromagnetic form factors. These involve the calculation of the χχ pair-production cross section, first, from e + e − annihilation accompanied by initial state radiation as shown in Fig. 1a and, second, from Bremsstrahlung in the scattering of electrons on nuclear targets as shown in Fig. 1b. They are the dominant production modes in BaBar or Belle II and in NA64, LDMX, mQ or BDX, respectively. Finally, since mQ and BDX search for new particles via their elastic scattering within the detector, Fig. 1c, we furthermore establish the nuclear and electron recoil cross sections in χ-scattering.
A. Pair-production at e + e − colliders In B-factories such as BaBar, Belle, and Belle II, fermion pairs can be produced in e + e − collisions accompanied by initial state radiation (ISR). Production in association with final state radiation (FSR) is suppressed with respect to ISR by 2 for Q and scales as µ 2 χ s/(4πα) ∼ (10 4 µ χ /µ B ) 2 for MDM and EDM. For the values of µ χ and d χ that are constrained through ISR by the various considered experiments, this gives a suppression of 10 −1 and smaller for FSR with MDM and EDM interactions. For AM and CR, the FSR diagrams vanish identically at tree level. Therefore, we are allowed to neglect FSR and only consider χχ production with associated ISR. Finally, mono-photon resonant production through a putative process Υ(nS) → γγ * → γχχ vanishes by Furry's theorem since Υ(nS) are C-odd states.
The ISR cross section approximately factorizes into the process without ISR, e + e − →χχ, times the improved Altarelli-Parisi radiator function that takes into account all soft and collinear radiative corrections up to order m 2 e /s [49,50], with the radiator function Here, x γ = E γ / √ s is the energy fraction carried away by the ISR, √ s 10 GeV and √ s χχ are the CM energies of the e + e − and χχ system, respectively, with s χχ = (1 − x γ )s; θ γ is the CM angle of the ISR photon. The χ-pair production cross section without ISR reads  (4) by integrating over cos θγ as a function of the ISR energy fraction, xγ = Eγ/ √ s, multiplied by xγ for Q (black), MDM/EDM (pink) and CR/AM (blue). All curves are normalized to the same total cross section.
Here, f (s χχ ) is the squared amplitude of γ * → χχ, integrated over the 2-body phase space of the fermion pair, as defined in (C7b). This factor will appear repeatedly and reads, Plugging (7) into (4), for m e → 0, we recover the expression in [51] for particles carrying Q and the expressions in [52] for MDM and EDM. In order to illustrate the difference between the various interactions, in Fig. 2 we show the differential cross sections with respect to the ISR energy for m χ = 0 (top) and m χ = 3 GeV (bottom). The different behaviors of the cross sections with respect to the fractional photon energy x γ suggest, that stronger limits on Q can be expected from high-energy mono-photon searches, while for AM and CR the low energy region appears more prospective in constraining these interactions (compare Fig. 4; to be discussed in more detail in Sec. IV). In the limit of m χ → 0 the models with the same dimensionality of the interaction operators show the same behavior, i.e. for µ χ = d χ (a χ = b χ ) the cross sections for MDM and EDM (for AM and CR) as a function of x γ are identical. Finally, we note in passing that in the presented numerical results we neglect the running of α since it is a few percent effect for energies up to 10 GeV.

B. Pair-production at fixed target experiments
In the fixed target experiments we consider in this work, a pair of fermions is produced via an off-shell photon emitted by an electron as it hits a nuclear target (e − N → e − N γ * → e − N χχ). The Feynman diagrams of the process are depicted in Fig. 1b. In principle, also inelastic scatterings to an excited nuclear state N * or to an all-inclusive final state X are possible. In the following we neglect such contributions, but present the general formalism to account for the latter in App. A and B.
Note that we only consider the emission of the photon from the electron leg. The emission from the nuclear leg is parametrically smaller by a factor of (Zm e /m N ) 2 ∼ 10 −6 for coherent photon emission and receives a further suppression from the nuclear form factor for photon momenta above a few 100 MeV. 3 Finally, in order to properly cover the full kinematic range that is accessible to χ andχ at a given CM energy, we compute the underlying 2-to-4 process exactly, without relying on approximations such as the Weizsäcker-Williams method. It is still possible, however, to relate the 2-to-4 cross section to a 2-to-3 process with the emission of a vector boson of virtual squared mass s χχ , in which the directions of χ-particles have been integrated out in f (s χχ ). 4 Equation (8) is exact, and can in principle be generalized for higher differential cross sections that are not sensitive to the direction of χ orχ. Our actual calculations require higher differential cross sections than (8). For NA64 and LDMX the full 4-body phase space is presented in App. C [Eq. C1] and for mQ and BDX where the momentum of χ enters in the elastic detection channel it is presented in App. D [Eq. D1].
In principle, also pions and η-mesons (among other mesons) can be produced in this scattering process, which can subsequently decay into χχ if kinematically allowed. However, while meson production dominates in fixed target experiments with a proton beam, it presents a subleading contribution for the electron beam experiments we consider. In addition, note that NA64 and LDMX should veto events that involve produced hadrons in the final state. For these reasons we neglect the contribution of hadronic production channels of χχ to the event rate.
In Fig. 3 we show the inclusive χχ-production cross sections as a function of m χ for several target materials and beam energies corresponding to the experiments discussed in Sec. IV. As expected from (7), the dimensionality of the operator determines the behavior of the cross section for m 2 χ s. Since models with lower dimensional operators show a flatter dependence on s χχ in (7), there is no visible difference in the inclusive cross section for Q at low masses between 4, 30 and 100 GeV beam energy in the top panel of Fig. 3, whereas for all other models, shown in the center and bottom panel, the cross section visibly rises with increasing beam energy. It is noteworthy that the cross section for an aluminium (Al) target is generally smaller than for tungsten (W) and lead (Pb) targets since the cross section scales with Z 2 and the charge of the Al nucleus is smaller by roughly a factor of 5-6. However, due to its smaller mass, the suppression at the kinematic edge is less pronounced in comparison to the heavier targets. The charge and mass differences between W and Pb are negligible in this case. Finally, while the cross section for millicharged particles rises with decreasing mass, it becomes independent of m χ for MDM and EDM, and even quicker so for AM and CR. When m χ increases towards the kinematic endpoint, the operators containing a factor of γ 5 yield velocity factors that result in a stronger suppression of the cross sections.

C. Detection via elastic scattering on nuclei
In mQ and BDX the produced fermions can directly be detected via elastic χ-nucleus and/or χ-electron scattering. In this case, the signal of χ particles would be a number of recoil events in the detector beyond the background, instead of missing momentum/energy. A presentation of the expressions of the recoil cross section dσ/dE R in the lab frame as a function of the recoil energy E R and valid for relativistic χ-energies E χ is relegated to App. E. The cross sections are functions of the nuclear mass (magnetic moment) m N (µ N ) and the nuclear spin I N . In the small-velocity limit, i.e., for |p χ | = m χ v and E χ m χ + m χ v 2 /2 with v 1, the expressions in (E2) are in agreement with the formulae found in Refs. [54,55]. Since the expressions are exact with respect to phase space, they apply to both χN and χe − scatterings.

IV. INTENSITY FRONTIER SEARCHES WITH ELECTRON BEAMS
In this section we place constraints on the discussed fermion models with electromagnetic form factors from BaBar, NA64 and mQ and present projections for Belle II, LDMX and BDX. All limits in the presented plots are obtained at 90% confidence level (C.L.). The limits and projections for the models with dim-5 (dim-6) operators are shown in the summary plots in Fig. 7 (Fig. 8). 5

A. BaBar and Belle 2
BaBar at SLAC [31] and Belle-II at SuperKEKB [33] are experiments at asymmetric e + e − colliders. Longlived states χ that are produced through their (feeble) electromagnetic form factor interaction can escape from the detectors. The primary signature is hence the missing energy in χχ production in association with a single photon (mono-photon missing energy search). The main background processes are SM final states in which electrons or photons leave the detector undetected, i.e., e + e − → γ / γ, which produces a peak at s χχ = 0 and e + e − → γ/ e + / e − as well as e + e − → γ / γ / γ which constitute a continuum background that rises with s χχ . The irreducible background e + e − → γνν is suppressed by m −4 Z and can be neglected.
For deriving constraints from BaBar, we use the data of the analysis of mono-photon events performed by the collaboration in a search for invisible decays of a light scalar at the Υ(3S) resonance [57,58]. 6 The CM energy is thus 10.35 GeV. Two search regions with different trigger and cut efficiencies, a high photon energy region with 3.2 GeV < E γ < 5.5 GeV and low energy region with 2.2 GeV < E γ < 3.7 GeV in the CM frame, were presented. The high energy search uses the full data set of 28 fb −1 integrated luminosity and the low energy search a subsample of 19 fb −1 . The geometric cuts of the search regions are −0.31 < cos θ γ < 0.6 (high energy) and −0.46 < cos θ γ < 0.46 (low energy), where θ γ is the photon angle in the CM frame with respect to the beam line as in (4).
With the respective angular cuts applied, the number of signal events is given by where eff is the total trigger and cut efficiency, L is the integrated luminosity. The boundaries of s χχ are given by the bins in Fig. 4, and the production cross section dσ e + e − →χχγ is given in (4). Following Ref. [59], in our analysis we apply a non-geometric cut eff of 30% (55%) in the high (low) energy region for BaBar and place conservative constraints on the models in (1) by requiring that the expected number of signal events does not exceed the observed number N (i) obs at 90% C.L. in any bin i. Concretely, we require N obs is given by the statistical and systematical uncertainties of N (i) obs added in quadrature. 7 For the Belle II projections, we follow Ref. [59] and scale up the Babar background from Ref. [57] (shown by the gray histogram in Fig. 4) to an integrated luminosity of 50 ab −1 with the CM energy of 10.57 GeV, 6 In more recent [58], Boosted Decision Trees (BDT), trained on simulated massive dark photon signal events, were used for event selection. Our cases of main interest have different event shapes as shown in Fig. 4, and hence differ significantly from the training set used in [58]. For this reason, we use the smaller data set [57] with cut-based event selection. 7 As is evident from Fig. 4, a signal will simultaneously be present at similar strength in many bins. It is hence safe to neglect the look-elsewhere effect. employ a constant efficiency cut of 50% in both search regions and take identical geometric cuts for Belle II and BaBar. While the exact background subtraction in the mono-photon search of Belle II needs to be provided by event generators [60], here we adopt the two projectionmethods of Ref. [59]: first is a "statistics limited" projection, where we assume a full subtraction of the background yielding a limit that is only governed by statistical fluctuations. It comprises the best case scenario with the strongest ensuing limits on the new couplings. Second is a "systematics limited" projection, where we assume a 10% systematic uncertainty for the γ / γ peak and a 20% uncertainty for the continuum backgrounds and exclude parameter regions where the signal exceeds the the systematic plus statistical uncertainties at 90% C.L.. With good control over the backgrounds, the actual reach of Belle II should lie between those two projections.
In Fig. 4 we show exemplary event shapes as a function of s χχ at BaBar for Q with = 0.2 and for AM with a χ = 4×10 −3 GeV −2 . As can be seen, the limit on millicharged particles is driven by the high mono-photon energy search region (top panel of Fig. 4) while for particles carrying an anapole moment (and in analogy the other operators) the limits are set in the low mono-photon energy search region as well as from the bins with lower photon energy in the high energy region. The ensuing constraints, as well as projected ones from Belle II, are shown as the parameter space of χ-mass and couplings in Fig. 7 for MDM and EDM and in Fig. 8 for CR and AM models.
For the Q model, we find that in the kinematically unsuppressed region, m χ 1GeV, the BaBar data excludes ≥ 6.4 × 10 −2 at 90% C.L. Moreover, the reach of the Belle II experiment is ∼ 2 × 10 −2 and ∼ 4 × 10 −3 for the systematics and statistics limited scenarios, respectively.

B. NA64 and LDMX
The NA64 experiment at the CERN SPS [35] and LDMX at SLAC [36] are current and proposed fixed target experiments using electron beams, where a χχ pair can be produced via off-shell photon emission. An active veto system together with cuts on the search region makes both experiments essentially backgroundfree. Therefore, in the following we place limits on the operators in (1) by assuming zero background in the regions of interest.
In the NA64 experiment, an electron beam with E 0 = 100 GeV hits an electromagnetic calorimeter (ECAL) which is a matrix of 6 × 6 modules consisting of lead (Pb) and scandium (Sc) plates, each module having a size of approximately 40 radiation lengths (X 0 ), with X 0 = 0.56 cm for Pb. The radiation length of Sc is about one order of magnitude larger, so the scattering of electrons with Sc is sub-leading, and will not be considered here. In the initial part of the ECAL there is a pre-shower (PS) detector of 4 X 0 in length with the same setup, functioning as target and detector simultaneously. Missing energy signals are required to have a starting point of the EM shower localized within the first few radiation lengths of the detector [61]. In our analysis, we conservatively take a thin target limit, with L target = X 0 , meaning that the fermion pair is essentially produced within the first radiation length of the target material. After scattering, the final state electrons are detected in the main ECAL with a search region of 0.3 GeV < E 4 < 50 GeV, where the lower boundary is due to the electron identification threshold in the PS detector and the upper boundary requires a missing energy larger than 0.5E 0 [35]. We adopt a polar angular coverage of θ 4 < 0.23 rad so that the final state electron should pass through the whole ECAL. The total cross section, however, is strongly forward-peaked and therefore the contribution of large angle scattering is negligible, meaning that a variation of θ max 4 does not affect our limits.
The number of signal events in the thin target limit and with the geometric and angular cuts as specified above is given by where dσ prod is the χ-pair production cross section, e − N → e − N χχ, given by (C6), N EOT = 4.3 × 10 10 is the number of electrons on target (EOT), ρ target is the target density and m N the target mass. The detector efficiency for NA64 only marginally depends on the electron energy, and we take it to be constant with eff = 0.5. The constraints from NA64 on the interactions in (1) are then derived by requiring that parameters that generate more than 2.3 events in the detector are excluded.
The LDMX proposal consists of a silicon tracker system, an ECAL and a hadron calorimeter (HCAL). In our analysis we use the benchmark points of the proposed LDMX phase I with 4 × 10 14 EOT at E 0 = 4 GeV, and phase II with 3.2 × 10 15 EOT at E 0 = 8 GeV [36]. LDMX will use a tungsten (W) target (X 0 = 0.35 cm) with a thickness of 0.1 X 0 in phase I and an aluminium (Al) target (X 0 = 8.9 cm) with 0.4 X 0 in phase II. The first four track sensors with an area of 4 cm×10 cm are located 7.5 − 52.5 mm behind the target along the beam line. We take a polar angular coverage in the forward direction of θ 4 < π/4 and a search region for the final state electron energy of 50 MeV < E 4 < 0.3 E 0 . While the actual signal efficiency varies mildly with the dark state mass and track requirements [36], a constant eff = 0.5 is assumed here. In the end, the number of signal events is calculated according to (10), and the projected constraints from LDMX are derived analogously, requesting at least 2.3 events in the detector.
In Fig. 5 we show the angular (top panel) and energy distributions (bottom panel) of the final state elec-trons for LDMX phase I obtained from (10). The event rates are rather similar for all models when normalized to the same total cross section; we normalize it to σ = 10 −36 cm 2 yielding approximately 1 predicted signal event in LDMX phase I. The signal is strongly forward peaked for s χχ m 2 χ , but in the transverse direction and at the kinematic edge, the curves start to deviate slightly. We also show the geometric acceptance region of LDMX in the top panel of Fig. 5 and the energy cut in the bottom panel, indicating that LDMX is searching in a region in the (θ 4 , E 4 )-plane where the event shapes of all models are similar.
Both the current limits from NA64 and projections for the LDMX phases are shown in Fig. 7 for MDM and EDM and in Fig. 8 for CR and AM models. Whereas the NA64 limit is weaker than the limit derived from BaBar data, LDMX will significantly improve the sensitivity to µ χ and d χ but will be on par with BaBar for a χ and b χ . We note in passing that for the Q model, our projections for LDMX agree with Ref. [36] and for NA64 are slightly weaker, by about a factor 2, than the limits presented in Ref. [62]. This is most likely due to a combination of a different assumption on detection efficiencies and the neglect of an energy threshold in the final state electron in [62].

C. mQ and BDX
In the mQ experiment at SLAC [40,53], N EOT = 8.4 × 10 18 electrons with E 0 = 29.5 GeV energy were incident on a W target of thickness 6 X 0 . The produced dark particles are then searched for in their scattering off electrons and nuclei in the main detector, approximately 110 m down the beam line. The detector, made of plastic scintillators, has a polar angular coverage of θ χ < 2 mrad, and a total path length of L det = 131 cm. The experiment found 207 recoil events above the estimated background, but below the estimated background uncertainty σ bkg = N bkg 382, within the signal time window. While millicharged particles would mostly produce single photo-electron (PE) events, such assumption is not justified for the recoil spectra of the models considered in this work. Therefore, similarly to [63], we focus on nuclear/electron recoil events with a deposited energy larger than 0.1 MeV. Following [53], 90% C.L. exclusion limits can then be obtained by requiring the number of signal events induced by χ, N sig < 207 + 1.28σ bkg 697. Here we will not consider a further possible reduction of single-PE backgrounds [63]; they may strengthen our limits by a factor of 2-3.
The proposed BDX experiment at Jefferson Lab [42] will use an electron beam of E 0 = 11 GeV incident on an Al target that comprises 80 layers with thickness of 1-2 cm each. The χ particles produced in the bump pass through thick shields and reach an ECAL 20 m downstream. The ECAL, composed of CSI(Tl) crystals, will have an area of 50×55 cm 2 and a depth of L det = 260 cm, Right: Corresponding electron and proton recoil spectra dNχ/dER as obtained from (12) for the same parameters.
implying a polar angle coverage of θ χ 12.5 mrad. The BDX collaboration has estimated that for N EOT = 10 22 the number of background events with deposited energy above 500 MeV is approximately 10, coming mostly from high energy neutrinos that scatter off electrons [64]. In the following, we conservatively take the uncertainty in the background estimation σ bkg = N bkg = 10, and assume a benchmark value for the projected number of observed events N obs = 15. It allows to set 90% C.L. limits on models studied here by requiring N sig ≤ 18.
In both experiments, the energy-differential expected number of particles χ that reach the detector is given by where the production cross section dσ prod given by (D9) has been integrated over the angular coverage of the detector, and is the integrated energy distribution of electrons during its propagation in the target [25]. For the mQ experiment, t 0 = L target /X 0 = 6 (X 0 = 0.35 cm for W), and for BDX we take t 0 = 15 as fiducial value (X 0 = 8.9 cm for Al). Finally, m N and ρ target are the target mass and mass density, respectively. The energy spectra (11) are shown in Fig. 6 for several benchmark parameter values. One can observe that most of the χ particles entering the detector carry kinetic energies in the GeV-ballpark. The total number of signal events in the downstream detectors that are produced in the scattering of χ on target i is then given by where (i) eff is the detection efficiency of an event with recoil energy E R . N (i) T is the target number density. For mQ, the signal is dominated by the scattering on e − and by the (coherent) scattering on C atoms. Hence i = e − , C and we compute N (i) T according to the chemical composition of the detector. For BDX, the detector is made of CsI, but given a threshold in deposited energy of E th R = 500 MeV, we consider i = p, e − as targets. The corresponding detection cross sections dσ (i) det /dE R are given by the elastic recoil cross sections found in App. E; the form-factors F E,M that are to be used for the respective targets i are provided in App. B. The energy of χ particles, E χ , is limited by the electron beam energy, E max χ , as given in (D10). In turn, E χ determines the maximal recoil energy of target i, E max R , see (E3). In the right panel of Fig. 6 we show the theoretical recoil spectra dN (i) /dE R that are obtained from (12) by dropping the integration over E R and before accounting for detection efficiencies, i.e. by setting eff = 1 for the same values of parameters as in the left panel of Fig. 6. For m χ 1 GeV the scattering on electron is most important for setting any limit on the interactions. Once m χ 1 GeV the maximum electron recoil energy E max R falls below the detector threshold E th R and any further sensitivity relies on recoil of protons (and, in principle, on neutrons). For completeness, we also note that for Fig. 6 we have neglected the momentum distribution of protons inside the target nucleus. The Fermi momentum of nu-   cleons is p F ∼ 250 MeV and it will hence affect the shape of the proton recoil spectrum at energies E R < 1 GeV, but be negligible for the bulk of the shown spectra.
For the value of the detection efficiency eff , we use 0.5 for mQ [53]. Based on the simulations in [42], a numerical function of E R has been adopted for electron recoils in BDX, giving eff ∼ 5% for E R = 0.5 GeV and eff ∼ 1% for E R = 1.5 GeV. However, there is currently no data available for proton recoils in BDX with E R ≥ 0.5 GeV, so we derive the BDX projections based on electron recoils only.
By requiring N sig , calculated from (12), to be below the values given above, the 90% C.L. (projected) limits on the electromagnetic factors of χ have been derived and are shown in Figs. 7 and 8. For m χ ≥ 0.5 GeV and E 0 = 11 GeV, no incoming χ can deposit energy more than 0.5 GeV by its scattering off an electron. So that the high-m χ end in the parameter space of Figs. 7-8 will not be constrained by electron recoils in BDX. In general, BDX will have stronger sensitivity than mQ to the models studied here. In contrast, the mQ experiment remains more sensitive to millicharged particles due to its low threshold in recoil energy.
The summary plots, Figs. 7 and 8, show that while the projections of BDX and LDMX for dim-5 interactions reach beyond Belle II for small m χ , for dim-6 interactions they flatten out below ∼ 0.1 GeV without surpassing Belle II in sensitivity. The difference of slopes of the BDX and LDMX curves in Figs. 7 and 8 is understood as follows: first, note that the χ-pair production rate receives significant contributions from the e-N forward scattering part which in turn is regulated by m χ ; the lower boundary of s χχ is 4m 2 χ . Second, as seen in Fig. 3, the m χ -dependence of the production rate becomes weaker for higher mass-dimension interactions, since the s χχ factors in the emission pieces, (7b-7e), tend to reduce the contribution of the forward scattering part. Taken together, this explains the difference in slope. It is also for these reasons that the millicharged case (dim-4) [53] has an even bigger slope than the dim-5 ones.

V. SM PRECISION OBSERVABLES AND HIGH ENERGY SEARCHES
Focusing on a mass of χ below several GeV, SM precision observables and LEP measurements in general play an important role in constraining the existence of the various electromagnetic interactions. In the following we update bounds on EDM and MDM SM precision tests first obtained in [5], going beyond parametric estimates where possible, and provide new ones on CR and AM. In addition we compute the constraints from flavor physics, as well as from missing energy searches at LEP.

A. Running of α
In the electroweak theory, Fermi's constant G F and the masses of W and Z bosons are related and receive calculable quantum corrections, summarized in the ∆r parameter [65]. The SM prediction is ∆r SM = 0.03672 ± 0.00019 [66] and the observed value can be obtained from where For the actual number we have used the global averages of M W , M Z , and A 0 and their errors quoted in [66]. The range in which new physics can contribute is hence limited at 95% C.L. as 8 −0.0038 < ∆r new < 0.00018.
Electromagnetic interactions of χ will contribute to ∆r through the running of α(Q 2 ) they induce. Explicitly, is the polarization function (Fig. 9a) at photon momentum q 2 = −Q 2 . Details on the calculation of the photon vacuum polarization are found in App. F. Equation (14) then implies for MDM and EDM independent of m χ , saturating the upper limit in (14) and improving the previously obtained limits in [5] by half an order of magnitude. In turn, for AM and CR we obtain saturating the lower limit in (14). The latter interactions hence contribute with opposite sign. The above limits are at or even stronger than the projected sensitivity from the various experiments considered above. As a quantum effect, they are, however, susceptible to a model-dependence that a direct observation of a particle is not. If the theory contains additional states at a mass scale below M Z or if χ carries more than one EM form factor, cancellations in ∆r new can occur. It is for this reason that we do not explicitly show those limits in Figs. 7 and 8. Finally, for completeness we also record a limit obtained from (14) on the positive contribution from Q, 0.1, weakening for increasing χ-mass to 0.3 at m χ = 10 GeV.
Here we consider the effect of the various interactions on the anomalous magnetic moment of the muon a µ = 1 2 (g − 2) µ . There is a tantalizing long-standing discrepancy between the measured value [67] and prediction in the SM (see [48] for a summary), indicating a 3−4σ tension. Hence, any positive contribution at the level of 3 × 10 −9 may solve the (g − 2) puzzle, and much anticipated upcoming experiments [68,69] will further weigh in on the anomaly in the near future.
In the current context, the lowest-order contribution to g − 2 enters through the vacuum polarization diagram, Fig. 9b where χ is in the loop of the internal photon line. In QED, the 2-loop leptonic contribution is of course well known [70,71], readily evaluated by noting the correction from the vacuum polarization to the photon propagator, and one may in principle proceed in analogy using (F3), appropriately renormalized. A quicker way is to estimate the correction ∆a µ from χ by using the dispersion relation and optical theorem from the total cross section of χ-pair creation, σ e + e − →χχ , similar to what is done with hadronic contributions in SM, where K(s)/(πα) is the contribution to the muon anomalous magnetic moment from a photon of mass √ s, 9 .
The cross section σ e + e − →χχ is given in (6) with s = s χχ . The integral (18) has to be cut off such that the validity of the effective theory is respected. This yields a reasonable estimate for dim-5 operators, and we show the resulting 90% bands in Figs. 7. Since K(s) ∼ 1/s, for dim-6 operators the energy dependence makes the same integral linearly sensitive on the chosen cut-off so that we do not show any region in Fig. 8. However, (18) is sufficient to convince oneself that both AM an CR only yield interesting contributions in an otherwise deeply excluded region; we leave a systematic study of this to a dedicated work.

C. Induced electric dipole moments
A sensitive and "background-free" probe of new physics are electric dipole moments of SM particles. Most recently, the ACME collaboration improved the limit on the electron EDM d e by an order of magnitude [73], Given the strong limit in (20), for simplicity, we restrict our attention to the electron EDM; for a general discussion of neutron and atomic EDMs cf. [9] and references therein. Electron electric dipole moments may be generated by the loops of χ-particles. Since the interaction is T-and P-odd it can only be induced in a loop diagram with an odd number of insertions of EDM vertices of χ, as it is the only T-odd interaction among the ones we consider.
Using dim-5 MDM and EDM operators, Ref. [5] then argued that such dipole moments appear at the 3-loop level with four photons attached to the χ-loop, as lower order diagrams vanish by virtue of Furry's theorem. The interaction d e hence appears through a combination d χ µ 3 χ and d 3 χ µ χ . Furry's theorem relies on C-invariance of the considered interactions and moving to dim-6 operators, we observe that the AM interaction is C-odd. Thereby, a χloop with 3 photons that are attached through EDM, AM and one out of the set of interactions { Q, MDM, CR} as shown in Fig. 9c does not vanish a priori. 10 However, when attached to the electron line, this does not imply a 2-loop contribution to d e as the resulting operator would be P-odd; we also verify the vanishing 2-loop contribution through explicit calculation using a projector onto d e [75].
Hence, we conclude that d e is only induced at the 3loop level, an example of which is shown in Fig. 9d. A parametric estimate of its size would then be, where the factor of m e originates from the necessary helicity-flip on the electron line to which 3 photons are attached (e 3 -factor). The factor (mass) −1 is built from the other interactions present in the diagram, saturated by powers of m χ or, possibly, m e to obtain the correct dimension, such as 3 e 3 /m χ , µ 3 χ m 2 χ , d 2 χ µ χ m 2 χ , b χ d 2 χ m 3 χ and so on. We have omitted any loop factors (16π) −2n with n = 3 as we do not know if other large combinatorial factors are present. 11 Saturating the present limit (20) we infer that the suppression mass scale in (21) is at the level of ∼ 10 13 GeV (d χ /µ B ) or larger. For example, if the electron EDM then solely appears through the combination of dim-5 operators, such estimate points to a related inverse effective scale (10 −5 − 10 −6 ) µ B GeV/m χ , implying an important constraint on the presence on joint EDM and MDM interactions. A systematic study of this goes beyond the scope of this work, but similar limits on other combination of operators are readily evaluated from the above arguments.

D. Rare meson decays
Measurements and searches of the invisible decay width of various mesons place a constraint on light dark 10 A way to see it is by following a proof for a generalized version of Furry's theorem [74]. Also note that the photons attached via AM and CR need to be virtual. 11 de induced by a possible EDM of the SM τ lepton, dτ , at (finite) 3-loop level has explicitly been calculated in [76] with the result de = a(me/mτ )dτ (α/π) 3 , where a turned out to be a number close to unity. For MDM/EDM interactions of χ, Ref. [5] saturates the dimensions by powers of me and includes a doublelogarithmic enhancement log 2 (mχ/me).
particles [77]. Strong constraints are in particular obtained from data on B + and K + decays. The decays K + → π + χχ and B + → K + χχ are closely related to the semi-leptonic decay with a charged lepton pair in the final state, K + → π + l + l − and B + → K + l + l − , as both are accompanied by the emission of γ * in the flavor changing s → d and b → s transitions. Hence, by evaluating explicitly the decay-details are found in App. G-we find for the MDM and EDM interactions and With the experimental value Br(K → πe + e − ) = (3.00 ± 0.09) × 10 −7 [66] together with the strong constraint on the branching ratio K + → π + + inv 4 × 10 −10 [78,79] which is applicable in the pion momentum ranges 211 MeV < p π < 229 MeV and 140 MeV < p π < 199 MeV allows to set a constraint in the χ-mass range m χ < 58 MeV and 76 MeV < m χ < 130 MeV. We obtain, approximately, and show the numerically obtained curves in Figs. 7 and 8. Equation (24) also shows that only MDM and EDM interactions are reasonably well constrained by rare Kdecays.
To estimate a constraint on the electromagnetic interactions from a limit on Br(B + → K + + inv) < 7 × 10 −5 [80] one can proceed in a similar fashion and relate the branching ratio Br(B + → K + χχ) to Br(B + → K + µ + µ − ) = (4.41 ± 0.23) × 10 −7 [66]. We obtain, for the MDM and EDM interactions and From those numbers we limit the interactions away from the kinematic threshold as, where the latter limit on CR/AM again implies a UVscale that would be well below the electroweak scale. The numerically obtained curves are again included in in Figs. 7 and 8.

E. Invisible Z-width
A model-dependent constraint arises if χ has fundamental interactions with the hypercharge through similar operators as in (1), but with F µν replaced by the hypercharge field strength F Y µν and A µ replaced by the associated gauge boson B µ . The Z-decay into χχ is possible for m χ < M Z /2 and the partial width is given by where the functions f (M 2 Z ) for the various operators are found in (7). The experimental measurement of the invisible with Γ(Z → inv) exp = 499.0 ± 1.5 MeV together with the SM prediction Γ(Z → inv) SM = 501.44 ± 0.04 MeV [66] limits any additional contribution to Γ(Z → inv) new < 0.56 MeV at 95% C.L. In the kinematically unsuppressed region this implies,

F. Missing Energy at high-energy colliders
At high-energy colliders, sub-GeV particles χ can be produced kinematically unsuppressed, yielding limits on the couplings that are independent of m χ . Among them, the strongest bounds on interactions studied here come from the LEP collider. Ref. [52] has studied the L3 data with an integrated luminosity of 619/pb at CM energies √ s = 188.6 − 209.2 GeV [81]. Using the mono-photon channel, it was found that MDM and EDM models are constrained to The results can also be generalized to AM and CR interactions. Here we only adopt the high-energy single photon selection, which requires one photon with transverse momentum larger than 0.02 √ s. Within the polar angular range 14 • ≤ θ γ ≤ 166 • , in total 1898 events were reported while the SM expectations are 1905.1 events, mostly from e + e − → ννγ(γ). Following [52], we use the the eight data subsets with different kinematic regions of √ s studied in the single photon search in [81] and place 90% C.L. bounds on the couplings from the data subset that leads to the best constraints using the CLs method. We take a flat selection efficiency of 71% and a systematic background uncertainty of 1.1%, and obtain |a χ | or |b χ | < 1.5 × 10 −5 GeV −2 (33) for the AM and CR models. Taking into account the low-energy selection would not significantly change these results. The limits are also shown in Figs. 10 and 11. While the bounds above are stronger than those derived from LHC data [82,83], future experiments, such as ILC and HL-LHC, may improve the sensitivity on these couplings by one to two orders of magnitude [84][85][86].

VI. χ DARK MATTER
If χ is long-lived on cosmological time-scales, it may be DM. Additional laboratory, astrophysical, and cosmological limits then apply, which we discuss in this section.

A. Direct detection
Our primary interest is in the GeV and sub-GeV mass scale in m χ . The elastic scattering cross sections on nuclei are given by the non-relativistic expansion of (E2) in the incident energy of χ, E χ m χ + m χ v 2 /2. In the case of the cross sections for EDM or MDM, the interaction of the dipole with the charge of the nucleus can be significantly enhanced by a small relative velocity. Direct detection experiments have put very stringent constraints on such models for m χ above several GeV [87][88][89].
The sub-GeV region is better probed through DMelectron scattering as it allows to extend the conventional direct detection down to masses of ∼ 10 MeV [90,91], below which the halo DM kinetic energy falls below atomic ionization thresholds. First limits on even lower mass DM have been achieved by the use of semiconductor targets [92,93], by utilizing either a solar reflected velocity component [94] or a cosmic-ray accelerated component [95].
Limits on DM electron scattering are conventionally quoted in terms of a reference cross section [90] where |M χe (q)| 2 is the scattering amplitude on a free electron, evaluated at a typical atomic squared momentum transfer q 2 = α 2 m 2 e ; we list various |M χe (q)| 2 in App. E.
The q-dependence of the actual ionization cross sections is moved into a DM-form factor and exclusion limits on σ e have been derived for a constant form factor [94,96] and for one proportional to 1/q [91]. For MDM and CR, which correspond to constant form factors, we combine the bound from [96] based on the XENON10 data [97] and the bound in the low mass region by considering solar reflection [94] based on the XENON1T data [98]. For EDM, corresponding to a 1/q form factor, the XENON10 bound derived in [91] has been adopted here. These limits are shown in Figs. 10 and 11. For AM, the χ-electron scattering is velocity suppressed, and we estimate a relatively weak limit in the range a χ 0.1 − 10 GeV −2 ; we hence omit it in the figures below. It should be noted that for the operators in (1) it is also possible that they induce transitions between mass-split (Majorana) states χ and χ * . Once ∆m ≡ |m χ * − m χ | ≥ 100 keV, such DM candidates can easily avoid any constraints from current direct detection experiments, see e.g. [54,99].

B. Dark radiation (Neff )
It is well known that (sub-)MeV mass particles, when they are thermally populated in the early Universe, will face severe limits from cosmology. χ particles with mass m χ few × MeV and couplings µ χ , d χ ≥ 10 −9 µ B or a χ , b χ ≥ 5×10 −5 GeV −2 reach chemical equilibrium with the SM thermal bath once the early Universe's temperature exceeded T γ ∼ m χ . For example, the annihilation of χ into electron and photon pairs-if it continues after neutrino decoupling-heats the photon bath relative to the one of neutrinos, and hence modifies the effective number of neutrino species, N eff , which in SM is given by N eff = 3.045 [100,101]. A general lower limit m χ 7 − 10 MeV has been derived for a Dirac fermion χ [102] based on the measurements by the Planck satellite, yielding N eff = 3.15 ± 0.23 [2]; see also [103][104][105] for earlier work. Since the bound is subject to a modeldependence when annihilation into other dark states (including neutrinos) is open, we do not show this bound explicitly in the figures below.

C. Kinetic decoupling
DM is required to kinetically decouple from the thermal bath after the photon temperature falls below 0.5 keV in order to avoid the over-damping of large-scale structures (LSS) [106,107]. The kinetic decoupling temperature can be obtained by considering the transport scattering cross section between χ and target i = e, p, where θ is the CM scattering angle, and the cross section have been factored as a product of a constant σ χi 0 and the n-th power of the relative velocity v. The expressions for the various operators in leading order of v are listed in App. E. We find that n = −2 for EDM, n = 0 for MDM and CR, and n = 2 for AM. Then, assuming T χ = T and demanding that the specific heating rate of DM particles derived in [108], be smaller than the Hubble rate, H(T ), from T ∼ 0.5 keV to recombination, we obtain upper limits on the couplings. For dim-5 operators, they are shown as dark shaded regions labeled "LSS overdamping" in Fig. 10.
As can be seen, the scaling of the cross section with v −2  makes the EDM case more stringently constrained. For fermions with dim-6 operators the bounds are too weak to be visible in Fig. 11. Our results on MDM and EDM agree with earlier studies [108,109].
The scattering processes tend to cool the gas in reverse. One can estimate the corresponding cooling rate using (36) by replacing n i with f i n χ , where f i is the fraction of electron/proton particles involved in the interac-tion and n χ gives the DM number density [110,111]. Therefore, cosmological 21 cm observations, such as recently reported by the EDGES experiment [112], may also constrain the couplings concerned here. Nevertheless, due to the minute free-electron fraction at redshift z ∼ 20, f e,p ∼ 10 −4 , the bounds will be much weaker than the ones from considering over-damping. Even for the EDM case with n = −2, in order to have the cooling rate comparable to the corresponding Hubble rate, one needs d χ ∼ µ B for m χ = 10 MeV. 12 D. Self-scattering of χ At low-redshift, a strong self-interaction among light DM particles modifies the shape and density profile of DM halos, as well as the kinematics of colliding clusters. A fair amount of attention has been devoted to the study of millicharged DM, which exhibits a strong velocity dependence of the self-scattering cross section. For the models studied here, however, the corresponding interactions are not enhanced at low relative velocities, which can be seen explicitly the expression of the cross sections listed in App. H. Hence, we may directly use the constraint on the self-scattering cross section obtained from from Bullet cluster observations, σ SI /m χ 1.25 cm 2 /g [114,115], to derive the bounds on MDM, EDM, CR and AM. Again, while the limits on CR and AM interactions are too weak to be shown, in Fig. 10 the ensuing constraints on MDM and EDM are shown by the shaded region labeled "σ SI /m χ ". To be general, we do not take into account the scattering between χ andχ particles, which, in any case, does not change our conclusion.

E. Supernova cooling
Light χ particles with m χ 400 MeV can be produced in pairs inside supernovae (SN), and, if such particles escape from the core, it increases the SN cooling rate. In turn, if the coupling is only strong enough, the particles may reach thermal equilibrium, be trapped, and an energy loss argument does not immediately apply.
We estimate the region in parameter space where χparticles affect the SN cooling as follows: if χ is thermalized within the SN core, a radius r d can be estimated as the radius where the blackbody luminosity of χ equals the luminosity in neutrinos L ν = 3 × 10 52 erg/s [116]. This is a maximum permissible luminosity if χ were to carry away this energy freely ("Raffelt criterion"). However, χ will typically be further deflected in elastic scatterings for r > r d , engaging in a random walk before it reaches a free streaming radius that we take as r inf = 50 km. In order to estimate the efficiency of the trapping, as a rough criterion, we impose as an upper limit above which energy-loss fails to be efficient; a similar argument has been used in [117] which we largely follow. We evaluate (37) by taking the SN numerical model of [118] with T χ = T (r), and protons, with mass density ρ p (r), are assumed to be at rest as we focus on m χ m p . The value of r d is then estimated from the SN model for each value of m χ , and it ranges from 22 km for m χ ∼MeV to 11.5 km for m χ ∼ 0.4 GeV.
On the flip side, the lower boundary in couplings for energy loss to be efficient can be estimated by calculating the energy production rate of χ particles by electronpositron annihilation within the SN core, following [84,119,120]. 13 The corresponding momentum distributions of electron and positron, f e − and f e + , as functions of the temperature and (opposite) chemical potentials for each r, are also given by the same SN 13 For interactions studied here, the contribution from plasmon decay is either comparable or sub-leading, and neglected for simplicity.
model [118]. Setting the core size r core = 10 km, we derive the lower boundary by requiring Both the lower and upper bounds discussed above are presented in Figs. 10 and 11. While these results are indicative of the relevant region in parameter space, a more detailed analysis may moderately but not qualitatively alter the results (e.g. [119,120]); we defer a more precise study of SN cooling to a dedicated future work.
In passing, we note that, besides SN, DM particles traversing the solar system may also be captured by the sun, affecting the stellar evolution. This has been studied in [55,121] for m χ above 4 GeV. For lighter χ particles as we are concerned with in this work, any effect will be strongly suppressed due to immediate DM evaporation; see also [94].

F. Dark matter annihilation and indirect search
Here we consider the case that χ is a symmetric DM candidate, so that χχ annihilation into light charged SM fermions and photons is operative both at CMB decoupling and at present day. The associated release of electromagnetic energy modifies the recombination history at high redshift and generates an excess in cosmic rays at low redshift.
For sub-GeV DM, the most stringent bounds typically come from CMB observations [122]. For MDM and CR DM particles, the annihilation to e + e − is s-wave, resulting in stringent limits on the relevant couplings. In contrast, for EDM and AM, where the annihilation cross sections to e + e − is p-wave and hence velocity suppressed, the leading bounds come from Voyager 1 data [123,124], and are relatively weaker. Moreover, CMB observations also strongly constrain the other annihilation channel, χχ → γγ. In contrast, gamma-ray telescope searches lead to bounds that are one to two orders of magnitude weaker in the sub-GeV region (see e.g. [125]), and are thus not shown in the figures. This channel is operative for MDM and EDM; for AM and CR, tree-level annihilation into physical photons vanishes identically. Here we will not further consider annihilation at loop-level (see e.g. [126]).
The bounds are summarized in Figs. 10 and 11. They show that MDM and CR, as symmetric DM candidates, are strongly constrained by the current observational data. Moreover, on-going and future experiments [112,127] have the potential to further improve the constraints on light DM annihilations soon.

VII. CONCLUSIONS AND OUTLOOK
In this paper, we consider the phenomenology of a new long-lived Dirac fermion χ in the MeV-GeV mass range and which interacts with the photon through mass dimension 5 and 6 operators MDM, EDM and CR, AM, respectively, but is otherwise neutral. In a bottom-up approach, the existence of such states χ is constrained from intensity frontier and high energy collider experiments, from SM precision observables, from flavor physics, and, once the lifetime exceeds roughly one second, from cosmology and astrophysics. We have considered all these possibilities in turn and derived constraints on the dimensionful couplings. When not present in the current literature, we also derive analogous constraints on a potential millicharge of χ.
At intensity frontier experiments where these fermions can be produced in pairs, we focus on electron beams and show that the most stringent constraints on these models are currently set by mono-photon searches at BaBar setting upper bounds on electric and magnetic dipole moments of ∼ 4 × 10 −5 µ B and on anapole moments and charge radius interactions of ∼ 2 × 10 −3 GeV −2 . Furthermore, we find that ongoing and future experiments, such as Belle II, LDMX and BDX, may extend the sensitivity by more than one order of magnitude. The projections for e + e − -colliders discussed in this work are basically independent of the fermion mass m χ across the whole MeV-GeV mass range, while the constraints from fixed target experiments exhibit a mass dependence that makes the interactions better testable in the low MeVmass range. These results need to be compared to monophoton searches at high energy colliders. We find that LEP is superior in terms of present constraining power, but Belle-II and LDMX can improve the reach to smaller couplings in the dim-5 case. In the dim-6 case, however, due to the higher dimensionality of the coupling, we find that LEP remains the best probe for AM and CR interactions.
Next, we use a number of SM precision observables as an indirect test to constrain the parameter space further, albeit in a possibly more model dependent way than a direct observation can offer. For example, corrections to the vacuum polarization of the photon from the new interactions induce a running of the fine structure constant α, and we place a constraint on the respective new physics models that is even stronger than the bounds from accelerator experiments. However, we find that some of the operators contribute with a relative sign and cancellations are in principle possible. We also estimate the contribution to (g−2) µ of the muon, concluding that the new interactions cannot contribute at the level of 3 × 10 −9 -reflecting the current tension in this observable -without being in conflict with other direct searches, most notably with mono-photon limits derived from BaBar-data. In turn, relating K ± → π ± χχ and B ± → K ± χχ meson decays to analogous channels with missing energy we find that flavor observables yield only relatively weaker limits and are superseded by other intensity frontier searches. On the other hand, if χ also couples to the hypercharge field tensor with similar strength as to F µν , the LEP observation on the invisible width of the Z-boson strongly constrains the whole parameter space which will be reachable by any of the future experiments discussed here. Clearly, the existence of such constraint hinges on the concrete UV-realization that induces the studied dim-5 and dim-6 operators.
Finally, we have addressed the question on the possibility of χ being DM. If χ is to provide 100% of the DM abundance and has a mass below the GeV-scale, direct detection constraints derived from DM-electron scattering are stronger than the limits derived from accelerators by 2 to 3 orders of magnitude for MDM and EDM. In contrast, we find that direct detection experiments have presently little sensitivity to AM and CR interactions when compared to other direct tests. In all cases, a combination of cosmological, astrophysical, direct detection and accelerator constraints rule out fermions with electromagnetic form factors making up 100% of the DM assuming a standard freeze-out scenario. If χ, however, only makes up a small fraction, 0.1%, of the DM abundance, one can evade the astrophysical constraints and will only be left with the accelerator and the supernova cooling constraints, which are independent of the DM relic density. For couplings smaller than µ χ , d χ ∼ 10 −9 µ B or a χ , b χ ∼ 5 × 10 −5 GeV −2 , DM particles of MeV mass do not reach thermal equilibrium with the SM anymore, allowing for alternative production mechanisms in the early Universe and therefore in this region of parameter space fermions with EM form factors could still be DM without violating any constraints. Note that ∆N eff bounds also do not apply in this region of parameter space.
This work constitutes a systematic study of sub-GeV mass dark states carrying electromagnetic form factors. Yet, many further avenues exist, which we are going to address in upcoming works: First, we have focused on electron beams employed in intensity frontier searches. Once we consider protons as projectiles, a whole number of relevant experiments for similar analyses become available. Most notable are the proton beam dump facilities, such as E613 at Fermilab [128], SHIP at CERN's SPS [129] and the proposed MilliQan experiment at the LHC [130], as well as the neutrino detectors like MiniBooNE [131], DUNE [132], LSND [133] and CEνNS [134]. Studies of dark sectors along similar lines have already been performed in [56,[135][136][137]. In order to systematically study χ-pair production from hadron beams, simulations are needed to generate event spectra because hadron physics cannot be neglected anymore. For example, in this set of experiments, contributions of the on-shell production and subsequent decay of mesons π 0 , η → γχχ and J/ψ, Υ → χχ will be non-negligible for the production of a χχ pair when kinematically allowed. Working along the lines of our presented appendices for the 2-to-4 processes, a systematic write-up of the production process and derivation of the relevant kinematic quantities that feed into the computation of observables is in preparation, filling a gap in the existing literature.
Second, in this work we have considered a Dirac fermion χ for concreteness. A complex, but electrically neutral scalar particle φ can have dim-6 CR interactions. In addition, at the same mass dimension, a vector particle may carry electric and magnetic quadrupole moments. Finally, at mass-dimension 7, the Rayleigh or susceptibility operators alter the phenomenology: the elementary interaction involves two photons, so that either an additional loop-suppressing factor is present in signals with one external photon, or new signal topologies with two external photons need to be considered. We leave a systematic study of those possibilities for a future paper.
Third, in this paper we have chosen a bottom-up approach, leaving the UV physics that generates the dim-5 and dim-6 operators unspecified. Such UV completion will result in additional constraints, most prominently from high energy colliders, in particular when the value of the effective dimensionful coupling points towards generating physics that must be at or below the TeV-scale. Such program for electroweak-scale particles has recently been started in [22]. Given the relatively weak intensity frontier limits on dim-6 couplings, our considered parameter space is likely further constrained by the presence of additional particles carrying electromagnetic charge and/or hypercharge.
Fourth, a more systematic study of cosmological and astrophysical constraints can be performed. For example, the limits on supernova cooling are derived on simplifying assumptions, and a more detailed study incorporating finite temperature effects and a better simulation of the efficiency of trapping may alter the inferred regions in parameter space, albeit only moderately. Finally, we have shown that a χ particle that freezes-out at the appropriate relic DM density is already ruled out. On the other hand, once the couplings diminish, χ may freeze-in, and a detailed computation of its abundance as a function of mass and interaction strength as well as a discussion on its detectability is still outstanding.
Dark sector particles with mass below the GeV-scale that are perfectly electrically neutral may still interact with the photon through a number of higher-dimensional operators. The question how strong such coupling can be and "how dark is dark" is an experimental one. In this paper we have answered it utilizing currently available experimental observables and provided forecasts to clarify the detectability with the future experimental program that is coming online or is considered.

ACKNOWLEDGMENTS
We thank M. Battaglieri, A. Hoang, Z. Ligeti, M. Passera, A. Ritz for helpful discussions. The authors are supported by the New Frontiers program of the Austrian Academy of Sciences. LS is supported by the Austrian Science Fund FWF under the Doctoral Program W1252-N27 Particles and Interactions. We acknowledge the use of computer packages for algebraic [138,139] and numeric [140] evaluations.

Appendix A: Matrix elements for emission
In this appendix we provide the detailed calculation of the pair production of two long-lived particles χ through their electromagnetic form factors. From the interaction operators (1) we derive the Feynman rules iΓ µ (q) with incoming photon four momentum q and Hermiticity of the electromagnetic current implies Γ µ (q) = Γ µ (−q) whereΓ = γ 0 Γ † γ 0 is the Dirac adjoint. We wish to compute the scattering cross section of e − on a nucleus N with emission of a χχ pair and into a potentially inclusive hadronic n-particle final state X n , where p 3 = n i=1 p 3,i is the total inclusive four momentum of the recoiling target, potentially fragmented into n states. We introduce the momentum transfer variables such that q = q 1 + q 2 which follows from overall energymomentum conservation. Since q, q 1 , q 2 are linearly dependent, we express the scattering by the six independent momenta p 1 , p 2 , q, q 1 , q 2 , p χ and their scalar products. The matrix elements for the processes depicted in Fig. 1b can split into the emission piece of the χχ-pair, M ν χ =ū(p χ )Γ ν (q)v(pχ), and the scattering piece on the nucleus, M N , where p 3 |J ρ (0)|p 1 = J ρ is the hadronic matrix element of the electromagnetic current. The overall matrix element for the emission is hence, where we have dropped all gauge-dependent pieces as is manifest for when q gets dotted into the vertices (A1).
In addition, the q 1 gauge-dependent part of the photon propagator in (A5) drops out when being dotted into the hadronic matrix element, q 1,ν p 3 |J ν (0)|p 1 = 0, so that we omitted those terms above as well.
The differential cross section for the process (A2) then reads, where v 1,2 are the velocities of the incoming particles and the total phase space Here p 5,6 = pχ ,χ . One can further introduce to (A7) an integral with respect to s X = m 2 X = p 2 3 to obtain dΦ = 2ds X dΦ 4 where dΦ 4 is the 4-body phase space of p j (j = 3, 4, 5, 6) -which will explicitly be solved analytically in Appendices C and D -and where the second line will be absorbed by the hadronic tensor W ρσ defined below.
In the lab-frame where |v 1 | = 0, the cross section can be written as where χ µν is the χ emission piece. The electron scattering, averaged (summed) over the initial (final) spins is described by the terms L ρσ,µν ab = 1 2 Tr ( / p 4 + m e )γ µ ( / p 4 + / q + m e ) The emission of the χχ pair is captured in the final state spin-summed matrix element χ µν given by,

Appendix B: Hadronic tensor and form factors
The response of the nuclear target is described by a hadronic tensor in (A9). In its general form, W σρ (−q 1 ) is given by, where the δ-function is from overall energy-momentum conservation and in this context guarantees that the fragments have the overall momentum p 3 . The most general (parity conserving) form for (B1) is where the form factors W 1 and W 2 are functions of q 2 1 and s X . In the elastic limit s X = m 2 N and q 2 1 = 2(p 1 · q 1 ) and W 1,2 are only functions of q 2 1 . For spin-1/2 fermions, they are related to the usual electric and magnetic form factors F E and F M via whereas for a scalar target they read, When considering the scattering on protons p or neutrons n, the form factors can be taken in the dipole form [141], F n E (t) = − µ n t (4m 2 n + 5.6t) and F p,n M = µ p,n F p,n E where t = −q 2 1 and µ p,n are the magnetic moments of proton and electron in units of the nuclear magneton µ N = e/2m p , respectively; with µ p = 2.79 and µ n = −1.91.
For nuclear targets of charge Z, the electric form factor is given by where A is the mass number and a(Z) = 111Z 1/3 /m e and d(A) = 0.164 GeV 2 A −2/3 [142,143]. The above expression includes a factor that accounts for the screening of Z by electrons such that F E (t) = 0 for t → 0. Finally, in our numeric evaluations we neglect the magnetic form factor for Z 1, as it is generally subleading. We note in passing, that for heavy states χ, inelastic scattering may start to play a significant role [25].
Appendix C: 4-body phase space for LDMX and NA64 NA64 and LDMX measure the energy and momentum of the final state electron. Therefore we want to show the most important steps in our calculation of the covariant phase space for the elastic scattering process in Fig. 1b with X n = N .
While a 4-body phase space has 12 degrees of freedom, 4 of them can be reduced by the momentum-energy conservation, and another one, the rotation along the axis of p 1 + p 2 , is redundant. That is, only 7 of them are necessary to express the squared amplitude, and we obtain here: The phase space can be described with the 5 independent Lorentz invariant integration variables s 3χχ = (p 3 + p χ + pχ) 2 , s χχ = q 2 = (p χ + pχ) 2 and u 2q = p 2 · q as well as q 2 1 = (p 1 − p 3 ) 2 and q 2 2 = (p 2 − p 4 ) 2 with all momenta defined as in App. A. The other 2 degrees of freedom come from Ω Rχχ χ , the solid angle of p χ in the CM frame of the χχ pair (which will be integrated out). φ R3χχ 3 is the azimuthal angle of p 3 in the frame where p 3 +p χ +pχ = 0, λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + ac + bc) is the triangle function. The Jacobian of the transformation from E 4 and cos θ 4 to the invariant variables s 3χχ and q 2 2 is given by The integration boundaries for the invariant masses s χχ and s 3χχ are The integration boundaries of the t-channel variables q 2 1 and q 2 2 are .
At last, the angular variable u 2q is given by from which one can easily obtain the boundaries of u 2q with φ R3χχ 3 ∈ {0, π}, as well as the expression of ∂φ R3χχ 3 ∂u2q . In (C5), G n is the Gram determinant of dimension n and ∆ n the respective Cayley determinant (or symmetric Gram determinant) [144]. The double differential cross section in the lab frame obviously follows from (C1) as For LDMX and NA64 where we are not interested in the direction of the χ pair, we can integrate out the part Ω Rχχ χ immediately, using the quantity defined below with f (s χχ ) found in the main text (7). This is because the rest of the squared amplitude is independent of Ω Rχχ χ .
Appendix D: 4-body phase space for mQ and BDX Since mQ and BDX ought to detect the produced dark fermions directly via elastic scattering on nuclei and electrons, we are interested in the direction of one of the fermions -denoted as θ χ below -instead of the direction of the final state electron.
The integration ranges of s 34χ are given by and that of t 2χ are And the upper and lower bounds of t 13 are given by At last, the integration ranges of φ 3 and φ 4 are both from 0 to 2π. Note that one Lorentz invariant needs to be given by the equation det[M] = 0, where the (i, j) entry of the 5 × 5 matrix M is the scalar product of p i and p j . This equality is due to the requirement that in fourdimensional space-time, any five 4-vectors cannot be linearly independent. We use this requirement to obtain the value of (p χ · pχ) from the seven Lorentz invariants above. Moreover, there are two solutions of p χ · pχ from det[M] = 0. This can be understood in the frame where p 1 + q 23χ = 0, as shown in Fig. 12. The two constraints derived from t 14 and p 3 · p 4 , illustrated by the circles at the bottom of the two cones, cannot uniquely fix p 4 (and thus pχ). Nevertheless, the degeneracy can be broken by fixing the rotation direction of φ 4 , as can be seen from Fig. 12. Here we take one solution of p χ ·pχ for φ 4 ∈ [0, π) and the other for φ 4 ∈ [π, 2π). Other Lorentz invariants are not affected by φ 4 ↔ 2π − φ 4 .
In the lab frame, we have where in the limit of m e → 0 the available range of E χ is and the corresponding range of cos θ χ is It is straightforward to also take into account the contribution of inelastic scattering between electron beam and the target by setting m X , or equivalently s X , variable, and by introducing the inelastic form factor of the target. Nevertheless, such a contribution is in general sub-leading, so following [42,53] we neglect it here and take m X = m N to obtain our numerical results .

Appendix E: Elastic scattering cross sections
Here, we present the full relativistic form of the differential recoil cross section of a nucleus with charge Z, electric and magnetic form factors F E and F M , mass m N and spin I N when scattering on a fermion χ with mass m χ and energy E χ , where A = E 2 χ − m 2 χ (2m N + E R ) and the functions g E (E R ) and g M (E R ) for all models are given by Q: EDM: AM: where t = 2m N E R ; the form factors F E,M (t) are given in App. A. The expressions above also apply to χe − scattering with F E,M = 1 and Z = −1 and for which the pre-factors of the electric and magnetic pieces agree by setting I N = 1/2 and µ N = µ B . For millicharged particle scattering in the ultra relativistic limit, i.e. for m χ → 0, we recover the Rosenbluth formula. Finally, the maximal recoil energy induced by χ carrying energy E χ is given by where m i = m N or m e . Equivalently, the minimal value of E χ to deposit energy E R in the target at rest is E min For completeness, we also list the squared scattering amplitudes |M χe (q)| 2 on free electrons, averaged over initial states and summed over final states, as they go into the definition of a reference cross sectionσ e (34) for estimating the direct detection limits. The leading terms are 14

EDM:
|M χe (q)| 2 = 64πd 2 χ αm 2 e m 2 AM: |M χe (q)| 2 = 16πa 2 χ αm 2 χ q 2 , CR: |M χe (q)| 2 = 64πb 2 χ αm 2 e m 2 χ . (E5d) Finally, in the limit of small velocities the transport cross sections σ χi T as they are used in (35) read, Here the momentum transfer of interest is, on average, at the same order as the temperature T . Hence we use F E = 1 and F M = −2.79 for the proton, neglecting the dependence on t, as we are interested in T = O(keV) in this context.