Scalar Dark Matter Candidates -- Revisited

We revisit the possibility of light scalar dark matter, in the MeV to GeV mass bracket and coupled to electrons through fermion or vector mediators, in light of significant experimental and observational advances that probe new physics below the GeV-scale. We establish new limits from electron colliders and fixed-target beams, and derive the strength of loop-induced processes that are probed by precision physics, among other laboratory probes. In addition, we compute the cooling bound from SN1987A, consider self-scattering, structure formation, and cosmological constraints as well as the limits from dark matter-electron scattering in direct detection experiments. We then show that the combination of constraints largely excludes the possibility that the galactic annihilation of these particles may explain the long-standing INTEGRAL excess of 511 keV photons as observed in the galactic bulge. As caveat to these conclusions we identify the resonant annihilation regime where the vector mediator goes nearly on-shell.


I. INTRODUCTION
The history of scalar dark matter (DM) in the MeV-GeV mass bracket is a long one. It has its roots when main-stream literature was primarily focusing on electroweak-scale new (supersymmetric) physics. Highenergy colliders had long explored the GeV-scale, and naive cosmological considerations suggested that thermal DM needed to have a mass of several GeV for the least [1,2]. However, light-scalar DM [3,4] turned out to be a perfect possibility-and still is. For example, it can couple to the Standard Model (SM) fermions either by a Yukawa-type interaction of some heavy fermions F or by a new gauge interaction mediated by a new vector particle Z [4].
Cosmologically, such DM retained its right of existence by achieving a sufficient annihilation cross section through an equally light Z (which has an even longer history [5,6]), or by a possible near independence of DM mass in the annihilation cross section when F is involved [4]. Experimentally, a light Z was viable, because the neutral current phenomenology remained largely unaffected, either from a suppression with center-of-mass energy arising from momentum-resolved diagrams in processes at high-energy, or, at low-energy, by ensuring that the effective strength of the interaction is smaller than the weak interactions of the SM. In turn, a model with F -mediation was even simpler to retain as the mass-scale of these particles can be in the TeV.
This, at the time seeming "niche physics" quickly gained momentum. When the SPI spectrometer on board of the INTEGRAL satellite confirmed a strong flux of 511 keV photons at the level of almost 10 −3 /cm 2 /s [7,8] (see also [9,10]) coming from the galactic bulge, a MeVscale DM origin was suggested on the basis of its spatial * celine.boehm@sydney.edu.au † xiaoyong.chu@oeaw.ac.at ‡ jui-lin.kuo@oeaw.ac.at § josef.pradler@oeaw.ac.at morphology and its general compatibility with the relic density requirement [11] while at the same time obeying soft gamma-ray constraints [3,12,13]. Concretely, the signal, especially its high bulge-to-disk ratio, is unexpected from known astrophysics [9,14] and calls for a new production mechanism of low-energy positrons. This can be achieved through DM annihilation into e + e − pairs. Decomposing the annihilation cross section in terms of the relative velocity as σ ann v = a + bv 2 and assuming a NFW dark matter halo profile, the observations suggest that the best fit values for the a or b parameters are [15], with a strong preference for a constant cross section (avalue) [16] albeit large uncertainties and an additional dependence on the cuspiness of the inner DM halo profile [17]. With these numbers in mind, the p-wave is roughly commensurate with the value required for a successful thermal relic σ ann v ∼ few × 10 −26 cm 3 s −1 (where v ∼ 0.3 at freeze-out). The second piece of early impetus for such models was their connection to low-energy precision physics, in particular to the electron and muon anomalous magnetic moments (g − 2) e,µ . There is a curious coincidence between the DM-viable parameters of the models, in part suggested by (1), and the lepton-mass dependent shifts to (g − 2) e(µ) on the order of 10 −11 (10 −9 )-essentially at the level of their observed magnitudes [11]. 1 This is a nice example of how laboratory probes of SM quantities inform us on the astrophysical and cosmological viability of new physics and vice versa [15].
A lot of work has been done since the introduction of the aforementioned scalar DM models. Regarding the INTEGRAL interpretation, the DM mass is now strongly constrained, e.g. from annihilation in flight, and the DMmass now generally needs to be below tens of MeV, and can even be fully excluded in specific models and/or under certain assumptions of the state of the early Universe; for more details see, e.g., [3,12,19,20] and references therein. In turn, the measurement [21] of the anomalous magnetic moment of the muon now stands in (3)(4)σ tension with the SM predicted one [22]. In fact, the imminent experimental update for (g−2) µ is much in the limelight today, especially after its connection with GeV-scale dark sector physics became more broadly appreciated following [23]. Finally, a beyond the Standard Model (BSM) sector that contains light dark states has found further motivations, such as the astrophysical core/cusp problem [24][25][26], various galactic cosmic ray excesses [27][28][29][30], among others. These motivations paired with the to-date absence of new physics at the electroweak (EW) scale have acted as a great innovation driver for devising laboratory and observational tests for sub-GeV dark sector physics, see, e.g. [31][32][33] and references therein.
In light of the significant amount of activities in the past two decades that has gone into the exploration of the MeV-GeV mass range and the large amount of results, it seems timely to revisit and confront the originally proposed models of sub-GeV scalar DM [4] to this wealth of data. The objective of this paper is to study the production of light scalar DM pairs at the intensity-frontier, consider the models' influence on SM precision tests, study the sensitivity to the new physics in high-energy colliders, at direct detection experiments and through astrophysical and cosmological observations. This will provide a more comprehensive assessment as whether light dark matter particles could explain INTEGRAL or the anomalous magnetic moment of the muon in this setup and also a summary of the constraints on new physics in this mass range. A summary of results for an exemplary DM mass of 10 MeV is shown in Figs. 1 and 2.
The paper is organized as follows: in Sec. II we introduce the models together with the parameter regions of central interest. In Sec. III the bounds from the intensityfrontier experiments, from precision observables and from LEP are derived. While such bounds constrain parts of the parameter regions of interests, complementary limits arise from cosmological and astrophysical observations, discussed in Sec. IV. Section V is devoted to a study of the low mediator mass region in the Z model. Conclusions are presented in Sec. VI. Several appendices provide additional details on our calculations.

II. REPRESENTATIVE MODELS
In this paper we shall focus on a complex scalar DM candidate. The Galactic 511 keV gamma ray line can then be explained by either t-channel or s-channel an-nihilation processes [4]. The former process necessarily involves an electrically charged particle, taken as a fermion below. Without loss of generality, the s-channel case assumes the presence of an intermediate gauge boson, which we shall take as leptophilic.
A. Heavy fermion mediator F In the first model that we consider, the scalar DM particle, denoted by φ, couples to the Standard Model (SM) via heavy fermionic mediators. For the sake of generality we take φ to be complex, but mention applicable formulae for real φ along the way. Concretely, φ and its antiparticle φ * may couple to the SM charged and neutral leptons l = (l − L , l − R ) T and ν l through a Yukawa-like interaction with the introduction of new electrically charged and neutral fermions F ± and F 0 , arranged as part of Here, (ν l , l − L ) and l − R are the SU(2) L doublets and singlets of lepton flavor l = e, µ, τ ; P L = (1−γ 5 )/2 and P R = (1+ γ 5 )/2 are the projection operators. We take all couplings to be real. In the presence of right-handed neutrinos ν R , additional interactions become possible, For the purpose of this paper, we shall not consider the latter option in any detail, but mention applicable results in passing.
There are a number of options related to (2), see e.g. [35][36][37][38][39][40][41]. In what follows, we usually drop the superscript on c l L,R for the coupling to electrons and electronneutrinos as we consider them as always present, c L,R ≡ c e L,R . Non-zero couplings to the second and third generations are a priori not the main focus of the paper, but they lead to further interesting consequences. Among them is a contribution to the anomalous magnetic moment of the muon, discussed below. If there is a single generation of heavy fermions F -which is the way how the Lagrangian is written-one may additionally induce lepton-flavor violating processes between the electron sector and muon or tau sector for c µ = 0 or c τ = 0, respectively (see below). At the expense of considering three generations of heavy fermions, F l , the flavor symmetry can be restored. Finally, we note that there is also a global dark U(1)-or Z 2 -symmetry in (2) between φ and F ; the former (latter) applies for F 0 being Dirac (Majorana).
the effective UV-scale Λ F = (c L c R /m F ) −1 . Before constraining the model, we infer the normalization points for the couplings from two particularly important predictions: the contribution to the anomalous magnetic lepton moment, as well as the DM annihilation cross section corresponding to the INTEGRAL signal and the thermal relic.

Anomalous magnetic moment
Under the assumption that m F m φ ≥ m l and that all c l -couplings are real, the one-loop contribution to the leptonic anomalous magnetic moment, shown in the left panel of Fig. 3, is given by in agreement with previous calculations [4,43]; note that a l ≡ (g l − 2)/2. Therefore, to address the long-standing muonic g − 2 anomaly [21] a exp µ − a SM µ = (290 ± 90) × 10 −11 , the corresponding favoured region is c µ F ≡ |c µ L c µ R | ∼ (5.5 − 7.6) × 10 −2 with m F = 1 TeV. The full expression without assuming the mass hierarchy is given in App. A. In anticipation of the constraints to be derived below, we point out that the contribution (4) to the electron anomalous magnetic moment will be of central importance when assessing the viability of explaining various anomalies.

DM annihilation
In the model with heavy fermionic mediators F ± and F 0 , the non-relativistic DM annihilation cross section into e + e − via F ± exchange or into Dirac electron neutrinos via F 0 exchange with the participation of a (kine- Belle II, stat.

FIG. 2.
Summary of constraints obtained in this paper for the vector-mediated model as a function of the effective scale Λ Z = ( √ g φ g l /m Z ) −1 for a fixed DM mass of m φ = 10 MeV in heavy mediator limit, similar to Fig. 1. A star indicates, that the bound only applies under certain conditions. The new/additionally shown bounds here are from the Z-invisible width (for m Z = 10 GeV only) and from parity violation using E158 under the assumption g φ = gL in the section " Precision Tests/LEP". Section "Astro/Cosmo/Direct Detection" now shows the annihilation constraint from Voyager 1 data. The region of interest for (g − 2)e, µ, where the bound for (g − 2)e is based on [34], also requires further assumptions of couplings g φ = gL = ±gR. The INTEGRAL interpretation is excluded. matically unsuppressed) light right-handed state given in Eq. (3),ν e ν R orν R ν e , reads where v M = 2(1 − 4m 2 φ /s) 1/2 is the Møller velocity. The s-wave component agrees with the one in Eq. (1) in [4], while the p-wave component is different, due to the fact that we expand in the Lorentz invariant product σ ann,F v M rather than σ ann,F v rel ; see App. B 1 for the full expressions, as well as those for real scalar DM. Above we have omitted terms that are suppressed by (m l, φ /m F ) 4 as well as higher-order terms. For the special case c L c R = 0 and for m l → 0 the above cross section vanishes, and the process becomes d-wave dominated [44][45][46], scaling as v 4 rel m 6 φ /m 8 F . Given a TeV-scale F and m l, φ well below GeV-scale, the latter terms do not contribute to the annihilation cross section in any appreciable way. Finally, for real scalar φ, a factor of four should be multiplied to the expression in Eq. (6) as both t-and u-channel processes contribute.
The annihilation to a pair of left-handed neutrinosν l ν l , mediated by F 0 , is either suppressed by neutrino mass or 1/m 4 F or v 2 /m 4 F , for Dirac neutrinos; see Eqs. (B1) and (B2). However, if neutrinos are Majorana fermions, one may additionally annihilate to ν l ν l or, equivalently, ν lνl with an s-wave cross section similar to (6) [47]. We comment on this possibility when considering cosmological constraints.
While we do not presume any production mechanism of the observed DM relic abundance, we will show the required parameters for thermal freeze-out below. Here, the DM abundance is Ω φ h 2 = 0.1198 [48], where Ω φ is the density parameter of φ and h is the Hubble constant in units of 100 km/s/Mpc. The observed relic density is achieved with c 2 [4]. The parameter regions that yield the required annihilation cross section within uncertainties for both the thermal freeze-out and the INTEGRAL 511 keV line [20] are shown in Fig. 6 for m F m φ .
Although the 511 keV line prefers a DM mass below tens of MeV, we scan over the entire MeV-GeV mass range. Meanwhile, we take the latest analysis of the INTEGRAL X-ray data [49], which suggests that DM The couplings g L,R and g φ are understood as a product of gauge coupling g and charge assignments q L,R and q φ so that g L,R = gq L,R and g φ = gq φ , respectively. Again, there are many options available with (7). They generally differ by the Z mass m Z , by their chiral couplings, by the absence or presence of family universality and/or kinetic mixing, by their (extended) Higgs sector, by potential additional fields that are required to cancel associated gauge anomalies in the UV and so forth; see e.g. [50] and references therein.
Here, we are primarily focused on the phenomenology associated with the Z coupling to electrons, and shall take g L and g R as flavor blind for when muons are involved. The special cases g l ≡ g L = g R and g L = −g R correspond to a pure vector and axial-vector interactions, respectively. 2 For the purpose of illustration, we consider m Z ≥ 10 GeV in most of our discussions. As will be shown, only a Z below the EW scale is of relevance for the INTEGRAL signal, so appears on-shell at highenergy colliders. As a result, although Z is generally off-shell for the low-energy phenomenology, and bounds derived below can be represented using √ g φ g l /m Z , results from LEP need to be treated with caution. For the latter, we provide bounds both on √ g φ g l /m Z in the heavy mediator limit, and on √ g φ g l for m Z m Z . The possibility of a Z below 10 GeV will be discussed separately in Sec. V.
Similarly as above, for the case m Z m l , the oneloop contribution to (g − 2) l is given by, in agreement with [54,55] if a pure vector coupling g L = g R ≡ g l is assumed. The full expression of Eq. (8) is found in App. A. The associated diagram of interest is shown in the right panel of Fig. 3, and the (g − 2) µ favoured parameter space is g l /m Z ∼ (4.6-6.4) TeV −1 . The constraint from (g − 2) e will be evaluated in Sec. III C. For a flavor-blind g l assumed here, the combination of several experiments excludes the possibility that this simple model explains the muonic g − 2 anomaly. This conclusion holds irrespective of if Z decays dominantly into SM leptons or into DM particles, as the leading constraint comes from the measurements of electron-neutrino scattering [56]. 3

DM annihilation
For the annihilation cross section via a s-channel Z , the s-wave component vanishes as scalars have no spin, and the p-wave component reads It agrees with Eq. (3) in [4] when taking v rel 2v φ , where v φ is the DM velocity. Since the cross section only varies by about a factor of two when either g R = 0 or g L = 0, we do not distinguish the left-and right-chiral couplings any further for annihilation, and simple take g L = g R ≡ g l in the remainder. For real scalar DM, the annihilation would be extremely suppressed since the Z does not couple to a pair of real scalars at tree level.
Taking the DM annihilation φφ * → Z * → l − l + with cross section as above, for m Z m φ > m l , the parameter region of interest for INTEGRAL is shown by the red and blue bands in Fig. 7 for NFW and Einasto profiles, respectively. Finally, we note that the observed DM relic abundance is achieved when [4] Pair production of φ in electron-positron annihilation in association with initial state radiation. Photon emission from the intermediate charged F is suppressed and hence neglected.
Pair production of φ in electron-beam fixed-target experiments. We consider φ emission from both initial and final state electrons (but not from the heavy F particle). Note that a global dark symmetry in Eq. (2) forbids the diagram with φ and φ * interchanged for the left process; see main text.
in the limit of m Z m φ . Obviously, m Z around or above the EW scale puts us into the non-perturbative regimes and is not of interest for us. Depending on the Z decay width, resonant annihilation at the point m Z 2m φ introduces additional velocity dependence in the annihilation. A detailed investigation of the resonant point, such as performed in [28,58,59], is beyond the scope of this paper. Due to these reasons, we focus on m Z ≥ 2.1m φ .

III. LABORATORY CONSTRAINTS
The dark sector particles may be produced in the laboratory, especially at electron-beam facilities, through electron-positron annihilation in colliders ( Fig. 4) or electron-nuclei bremsstrahlung in fixed-target experiments (Fig. 5). Moreover, they can also appear virtually through loops, affecting EW precision measurements. Such considerations thus put upper bounds on the coupling of SM particles to the dark sector.
We briefly introduce the experimental data of interest and our methods to derive the related constraints below, and refer to [60] and our appendices for further details of relevant cross sections.

A. Electron-beam facilities
We first consider intensity frontier experiments, including low-energy electron-positron colliders and electronbeam fixed-target experiments. For the values of m F and m Z concerned above, we can only produce φ via off-shell mediators in these experiments.
Following our previous work [60], we derive the expected number of signal events and constraints from current experiments such as BaBar [61], NA64 [62,63] and mQ [64,65], as well as projected sensitivities for future ones, including Belle II [66,67], LDMX [68] and BDX [69]. Depending on the observable signatures, these experiments can be put in three categories described below.
The first category is to look for large missing transverse momentum/energy, accompanied by a mono-photon signal, in low-energy electron-positron colliders, such as BaBar and Belle II. The expected number of signal events in each energy bin reads where eff is the efficiency, L is the integrated luminosity, θ max,min γ are the cuts on the photon angle in centre-ofmass (CM) frame w.r.t. the beam axis, is the energy fraction carried away by the photon with s and s φφ = (1 − x γ )s being the CM energy square of e − e + and φ-pair, respectively. The differential cross section is found in App. B 2. For BaBar, we take the data of the analysis of mono-photon events in a search for invisible decays of a light scalar at the Υ(3S) resonance [61]. The CM energy is 10.35 GeV, with two search regions of 3.2 GeV ≤ E γ ≤ 5.5 GeV and 2.2 GeV ≤ E γ ≤ 3.7 GeV. 4 For Belle II, we follow [70] and derive the projection by scaling up the BaBar background to an integrated luminosity of 50 ab −1 with a similar CM energy. The second category describes the missing energy search in electron-beam fixed-target experiments, such as NA64 and future LDMX. The expected number of signal single-electron events is given by × cos θ max state electron and θ 3 is the scattering angle w.r.t. the beam axis of the final state electron in the lab frame, with its detection efficiency given by eff (E 3 ). The differential cross section is derived in App. B 3. The background in such experiments is usually negligible after imposing stringent selection criteria. The NA64 experiment uses an electron beam with E beam = 100 GeV and has collected data of N EOT = 4.3 × 10 10 . We select events with only a final state electron, with its energy between [0.3, 50] GeV and θ 3 ≤ 0.23 rad. For the proposed LDMX experiment, we use the benchmark values of phase I with N EOT = 4 × 10 14 at E beam = 4 GeV and phase II with N EOT = 3.2 × 10 15 at E beam = 8 GeV.
The energy and geometry cuts on final state electrons are 50 MeV < E 3 < 0.5E beam and θ 3 < π/4. A constant eff = 0.5 is taken for both experiments. The last category includes mQ and BDX, which are electron-beam fixed-target experiments designed to directly observe φ-e (or φ-nucleon) recoil events in a downstream detector. The expected number of electron recoil events is given by where n e is the electron number density in the detector, L det is the detector depth. The threshold recoil energy E th R depends on the experiment and E max with the exact differential recoil cross section given in App. B 4. The production spectrum of φ is computed by in which the factor 2 accounts for the production of the φ-pair, θ φ is the scattering angle w.r.t. the beam axis of the produced φ in the lab frame with boundaries given by the geometry of the downstream detector and I(E) is the integrated energy distribution of electrons during their propagation in the target [60]. The differential cross section for the φ energy distribution is listed in App. B 3.
For the mQ experiment, the incoming electron with energy E beam = 29.5 GeV impinges upon a tungsten target with N EOT = 8.4 × 10 18 . The collaboration has reported 207 recoil events above the background, which is below the uncertainty of the latter σ bkg = 382 within the signal time window. We derive the upper bounds on the dark sector couplings from events with electron recoil energy E R ≥ 0.1 MeV. For BDX, electrons with E beam = 11 GeV are incident on an aluminium target which comprises 80 layers with thickness of 1-2 cm each.
The BDX collaboration estimated that for N EOT = 10 22 the number of background events with E R ≥ 0.35 GeV is about 4.7 [71]. Again, we only consider electron recoil events, with a detection efficiency of 100% for mQ [65], and 20% for BDX [71].

B. High-energy colliders
High-energy colliders may produce any of the dark sector particles studied here, leading to missing energy signatures. In the F model, a TeV-mass charged fermion F remains largely unconstrained by current bounds from LEP or LHC data, while the missing energy search in LEP [72] is able to constrain the overall coupling as which can be improved by investigating DM production via Drell-Yan processes with high-luminosity LHC [39], as well with ILC [73]. For the Z model, the LEP bound varies depending on the Z mass. For a heavy Z above the LEP energy scale, we obtain a bound from missing energy events induced by DM pair production as in agreement with previous results [74,75], which is stronger than the reach of low-energy beam experiments. Although this is shown in Fig. 7, it does not apply to m Z below the LEP energy scale, where a more proper LEP bound may come from missing energy induced by on-shell Z production, requiring g l 0.01 [76]. Its combination with the perturbative condition being comparable to the BaBar bound for m Z = 10 GeV. Naively speaking, these two LEP bounds, valid for different parameter regions of m Z , converge at m Z ∼ m Z . Projected sensitivities on a leptophilic Z portal have also been derived for future colliders, see [53,[77][78][79][80].

C. Precision observables induced by loops
The measurement of the fine structure constant, α, has been improved significantly with Cs atom interferometers [34]. Taking as input α ≡ α (Cs) , the SM prediction of the electron anomalous magnetic moment a Belle II, stat.
Belle II, syst. with a nominal 3σ requirement. This discrepancy could also be rephrased in a tension between α extracted from Cs and from a e measurements using the standard model prediction, α(a SM e ), for which α (Cs) − α(a SM e ) < 0. Both models-through their contribution to a e -then imply an inferred shift in the value of α. One should obtain the same constraints from both.
In the F model, positive (negative) c L c R yields a positive (negative) contribution; cf. (4) or the full expression in App. A. As shown in Fig. 6, either sign then puts a strong constraint on the model with a F mediator. For the positive product c L c R > 0 we obtain c L c R /m F ≤ 6.2 × 10 −5 TeV −1 and for the negative product c L c R < 0 we obtain |c L c R |/m F ≤ 6.1 × 10 −4 TeV −1 . For the Z model, where ∆a (Z ) e scales as 1/m 2 Z , the resulting 3 σ limit based on (19) requires g l /m Z 9.5 TeV −1 for g l = g R = g L , and the limit becomes even weaker if g R = −g L . It is always surpassed by the LEP bound above, for both m Z m Z and m Z ≤ m Z , and is hence not included in Fig. 7; see Fig. 2. 5 We may exercise some caution in applying (19), as it takes a positive half-σ shift to rule out any model by increasing the 2.5σ tension to 3σ. Here we stress that both the F -and Z -mediated models allow for both signs in 5 Another observable is the running of the fine structure constant, given by the photon vacuum polarization induced by the charged F -loop, Π(−M 2 Z ) − Π(0). This number needs to be below 0.00018 [82], requiring m F 80 GeV. The formula for Π(p 2 ) is given in Eq. (C8) with g l replaced by e. A dark U(1) gauge boson Z does not contribute to the running at one-loop. their contributions. Therefore, going in the other direction, one may first bring both measurements into reconciliation and in a further consequence, allow for a particularly large shift before the lower boundary in (19) is reached. In this sense, (19) entails both, an aggressive and conservative limit. In Fig. 6 we show both ensuing constraints as well as the region of interest for alleviating the tension.
The invisible decay Z → φφ * induced by the 1-loop diagram containing F or Z will alter the decay width of Z, see Fig. 9. Such additional contribution is bounded by experiments [83] to satisfy Explicit calculation of the relevant loop diagrams, detailed in App. C, reveals that the ensuing constraints (c F /m F < 26.6 TeV −1 and √ g φ g l < 0.35 for m Z = 10 GeV) are weaker than those above from general missing energy searches. We hence do not show this constraint in the figures.

D. Further precision tests
Before closing this section, we also mention that further limits arise when the models introduce new sources of parity or flavor violation.
The effect of parity violation can be parametrized as a deviation from the SM-predicted value of the weak an-gle θ W . 6 The E158 experiment at SLAC used polarized electron beams of E e 46-48 GeV, and measured the Møller scattering asymmetry with one polarized electron at low momentum-transfer Q 2 0.026 GeV 2 [87,88]. The subtracted value of 4 sin 2 θ W − 1 from the data is −0.0369 (52), slightly higher than the SM prediction, −0.0435(9) at low energy. For the Z model with m Z Q 2 , the resulting bound scales as g R(L) /m Z 0.38 (0.58) TeV −1 at g L(R) = 0 and almost vanishes at g L = ±g R . Moreover, at g L = g R cos 2θ W /(cos 2θ W − 1), the contribution of Z can be absorbed by re-scaling the Fermi constant. Although for m Z ≥ 10 GeV concerned in this section, it is at most comparable with the LEP bound above, future experiments, such as MOLLER [89] and P2 at MESA [90], have the potential to improve the limit by more than one order of magnitude [91,92]. In contrast, the F model does not induce additional electron scattering processes at tree level, and thus can hardly be constrained by such experiments.
Moreover, if only one generation of F is present, nonzero couplings to the muon and tau sector induce lepton flavor violation, similar to flavored DM [35]. For example, one may have the decay µ → eγ by closing the φ-loop, effectively via the magnetic dipole interaction [93]. The current strongest limit, from the MEG experiment [94], requires Br µ→eγ 4.2 × 10 −13 . This in For purely chiral interaction with c l L = 0, the bound is relaxed to |c µ R c e R |/m F 0.05 TeV −1 as an additional spin-flip is needed [95]; the same bound applies if one switches the chiralty subscripts.
Another example is the decay µ → eφφ at a rate compared to the SM mode, The first scaling applies for (c µ L c e R ) 2 + (c µ R c e L ) 2 = 0; the second scaling applies for when the same combination of couplings vanishes. Hence, because of the coupling structure and because of kinematic effects, one generally expects distortions of the electron spectrum. The latter is an important test for the V − A nature of weak interactions and has been mapped out well in the coefficients describing it [83,96]. The sensitivity is, however, likely superseded by the radiative decay above, and only applies to m φ ≤ m µ /2; a detailed study of it is beyond the scope of this paper. 6 Here we note that there is no sensitivity from the feats that detected atomic parity violation [84], because of the leptophilic nature of couplings involved; see [85,86] for when a Z coupling to quarks is present.
Finally, we mention that for a Z boson that couples to quarks/leptons with appreciable strength, precision observables were also investigated in [97,98]; note however that stringent constraints from dilepton resonance searches derived in the latter work are avoided, as in our setup Z dominantly decays into a φ-pair.

IV. COSMOLOGICAL AND ASTROPHYSICAL CONSTRAINTS
As the scalar φ is assumed to be the dominant DM component, the models are also constrained by astrophysical and cosmological observables, as well as from DM direct detection experiments. These constraints are discussed in the following.

A. BBN/CMB ∆Neff bounds
Here we take into account the BBN/CMB bounds on N eff from early Universe observations, while at the same time remaining agnostic about the state of the Universe for T MeV. Since m φ ∼ O(MeV), during Big Bang Nucleosythesis (BBN) φ can still be relativistic and contribute to the radiation density, summarized in the parameter N eff . Recall that we always set m Z ≥ 2.1m φ , so Z only plays a sub-leading role in the radiation density budget, even though it has three degrees of freedom. Currently, two relativistic degrees of freedom, like from a thermal complex scalar, are still considered to be marginally allowed by BBN on ∆N eff [99][100][101].
In contrast, the Planck measurement of the cosmic microwave background (CMB) spectrum requires that N eff = 2.99 ± 0.33 at the last scattering surface [48]. This limits the residual DM annihilation after neutrinodecoupling that injects energy either into the visible or into the neutrino sector [102,103]. In the F model, φ pairs annihilate into electrons. Under the assumption of a sudden neutrino decoupling at 1.41 MeV [104], we obtain a lower bound from N eff as m φ ≥ 5.1 MeV for a complex scalar, consistent with previous results [20,105]. 7 However, the CMB bounds from N eff become much weaker if the scalar DM annihilates into both electrons and neutrinos, which happens in the F model with Majorana neutrinos, as well as in the Z model. The underlying reason is that both, the photon-and neutrino-fluid is being heated so that the ensuing offset in the ratio of their respective temperatures is milder. In a flavor-blind set-up assuming DM annihilates to electrons and each species of SM neutrinos equally, we then estimate that the Planck bound on N eff only requires m φ 2.0 MeV. The latter possibility was not considered in [20,105].

B. Direct detection
As we focus on sub-GeV mass φ particles, we consider φ-e scattering signal as it has a lower threshold on DM mass in direct detection experiments. Exclusion limits are customarily presented in terms of a reference scattering cross section [106], where |M φ-e (q)| 2 q 2 =α 2 m 2 e is the squared matrix element of φ scattering on a free electron, summed over final state spins and averaged over initial state spin, evaluated at a typical atomic momentum transfer q = αm e . To order O(v 2 rel ) it is given by for the two representative models. Note that bounds on σ e have been obtained for the present case of constant DM form factors, most recently in SENSEI [107], which also summarizes previous bounds from XENON10 [108,109] and XENON1T [110], as well as from considering a solar-reflected DM flux [111]. The corresponding constraints, combining the results of experiments mentioned above, are shown in Figs. 6 and 7.

C. Indirect search
To explain the INTEGRAL signal, φ has to be a symmetric DM candidate, implying φ annihilation into SM leptons both at epochs of BBN and CMB, as well as at low redshift. Among them, bounds from BBN observables and DM annihilating to neutrinos [105,112] are very weak, and are not further considered. Since in the considered models φ does not annihilate into photons at tree-level (except when accompanied by final state radiation), we focus on the channel φφ * → e + e − .
For the F -mediated case, in which both s-wave and p-wave annihilation are present (6), it turns out that the constraint from CMB [113] is in general stronger than that from Voyager 1 data [114]. In the Z -mediated case, since the leading contribution of φφ * → e + e − is p-wave (9), the annihilation at CMB epoch is velocity suppressed and the bounds from present-day data such as from Voyager 1 [115,116] is more stringent, disfavoring DM masses above O(30) MeV to explain the INTEGRAL 511 keV line. This will be further improved by about one order of magnitude in future experiments, such as e-ASTROGAM [117,118] and AMEGO [119]. The CMB constraint for the F case and the Voyager 1 constraint for the Z case are shown in Figs. 6 and 7, respectively.

D. Structure formation
To avoid the collisional damping of DM primordial fluctuations [120,121], DM has to kinetically decouple from the observable sector in the early Universe. In the considered models, DM couples to electrons and neutrinos with similar strength. Since the number density of electrons is much lower than that of background neutrinos once T m e , the scattering on neutrinos hence governs the ensuing constraint. Here we take the bounds derived in [122,123] for both energy-independent and energy-dependent DM-neutrino scattering cross sections. Concretely, we require that for the F model and for the Z model The requirement consequently leads to c 2 F /m F 0.25 (m φ /MeV) 1/2 TeV −1 , as well as √ g φ g l /m Z 2.17× 10 4 (m φ /MeV) 1/4 TeV −1 . Both bounds are weaker than those obtained above, and are not shown in the figures.

E. DM self-scattering
If φ constitutes DM, its self-interaction may change the shape and density profile of DM halos, and the kinematics of colliding clusters. Such self-interaction is apparently very weak in the heavy F model. In the Z model, the DM particle φ can efficiently selfscatter via Z exchange. The self-scattering cross section averaged over φφ → φφ, φφ * → φφ * and φ * φ * → φ * φ * reads [124] where velocity-suppressed terms have been neglected. 8 However, the current bound, σ SI /m φ ≤ 0.5 cm 2 /g from cluster observations [125][126][127][128][129], is also not able to provide any meaningful bounds on the Z model with a 10 GeV m Z .  (18) applies instead. We do not show a band for (g − 2)µ, which would need an assumption on g φ /g l , since it is already excluded elsewhere (see main text and Fig. 2).

F. Anomalous Supernovae cooling
An important constraint arises from the anomalous energy loss via φ production in hot stars, especially inside supernovae (SN), as we consider m φ = O(MeV-GeV) which has overlap with the SN core temperature. To avoid the suppression of neutrino emission from the SN core after explosion, we impose the so-called "Raffelt criterion", which states that the energy loss via dark particle production has to be smaller than the luminosity in neutrinos, L ν = 3 × 10 52 erg/s [130]. 9 Here we follow the method in [60,133,134], and adopt the SN1987A numerical model of [135] with a total size r SN = 35 km, to derive the bounds on the leptophilic DM models above.
The dominant φ production channel is pair creation from electron-positron annihilation. As our mediator particles, F or Z , are much heavier than the core temperature of SN, we can safely neglect thermal corrections. Quantitatively, the lower boundaries of the exclusion regions are derived by requiring 9 The bounds from SN1987A are derived from the cooling of the proto-neutron star; doubts exist if SN1987A was a neutrinodriven explosion [131] in which case the limits become invalidated. Such speculation could be resolved once the remnant of SN1987A is firmly observed [132].
where r c is the core size of SN1987A, taken as 15 km here and f e − ,e + are Fermi-Dirac distributions for electron and positron.
On the flip side, if the coupling between φ and SM particles inside the SN core is so strong that the φ becomes trapped inside the core, the energy loss via φ emission diminishes and again drops below the neutrino luminosity. 10 To estimate the corresponding upper boundaries of couplings, we first define a radius r d , where a thermalized blackbody luminosity of φ equals L ν [60,133]. Consequently, if φ gets further deflected by scattering off electrons for r > r d , the effective dark luminosity drops below L ν . In practice, we impose Here, the mean-free-path λ φ-e at each radius is calculated by averaging the momentum distributions of both φ and electrons where the Pauli blocking factor induced by final state electrons has been taken into account. The relevant differential cross sections are listed in App. B 4.
The resulting SN1987A exclusion regions, combining both lower and upper boundaries, are given in the right panels of Figs. 6 and 7. Our lower boundaries agree well with previous results [137,138]. Regarding the upper boundaries, we find that Pauli blocking plays an important role in suppressing φ-electron scattering. Meanwhile, although there is little Pauli blocking in φ-nucleon scattering, φ only couples to quarks at loop level, yielding a suppression by another factor m 2 φ /m 2 F for the F model and α/π for the Z model. It hence turns out that φ-electron scattering dominates the capture in the parameter regions studied here.

V. Z MASS BELOW 10 GEV
Now we consider m Z 10 GeV, for which Z may appear on-shell in intensity-frontier experiments, CM energy permitting. In this case, the effective operator approach with its scale given by √ g φ g l /m Z does not apply any longer. Nevertheless, the DM interpretation of the INTEGRAL signal only relies on the product g φ g l , instead of on each of the two couplings individually, up to a minor effect of Z decay width. Therefore, to obtain the most conservative constraints on g l , we choose to maximize the dark coupling, g φ , using the cluster bounds [126,128,129] above We also impose a perturbativity bound, requiring that α D = g 2 φ /(4π) should not exceed 10. There are then three parameters left. To explain the INTEGRAL signal, Eq. (1) allows to solve for the value of m Z for each choice of (m φ , g l ). Now we directly adopt the existing constraints (e.g. NA64 and BaBar) and projected sensitivities (e.g. LDMX and Belle II) on an invisibly decaying Z from [139]. For direct detection, we take the bounds with electron recoils summarized in [107], and future projections of SENSEI at SNOLAB in [140]. Interestingly, for m φ below 100 MeV, this may be further improved by orders of magnitude using potential signals from neutron star heating [141]. Indirect bounds on p-wave DM annihilation are relatively weaker [105], and are not considered further.
The results are summarized in the left panel of Fig. 8. Below the dash-dotted line, a self-interaction cross section saturating (31) is achieved with α D < 10. Above the line, we set α D = 10 and the bound (31) is not saturated.
Assuming smaller values of α D necessarily requires larger values of g l while their product has to be large enough to explain the INTEGRAL signal. Therefore, this choice of α D gives the most conservative bounds, and choosing smaller α D would exclude more parameter space in Fig. 8. The green and red lines in the figure are the contours of m Z and m Z /m φ , respectively. While m Z = 3m φ has often been chosen in the literature, we also show the contour of m Z = 2.1m φ , below which the average DM kinetic energy may overcome the mass barrier to produce intermediate Z on-shell at freeze-out or after. As is shown in the figure, for each value of m Z the mass ratio already becomes important in determining the annihilation cross section at m Z = 3m φ .
Our results suggest that most of the parameter region with non-resonant DM annihilation has been excluded, while the LDMX experiment is projected to probe the whole region of m Z 2.1 m φ . This is different from the heavy m Z case above, where Belle II always has a better sensitivity in probing the DM particle. In presence of resonant annihilation, the required values of g l can be much smaller than what fixed-target experiments will reach, thus indirect DM searches, such as e-ASTROGAM [118], are expected to better probe this parameter region.
Similarly, instead of calibrating on INTEGRAL, we may alternatively fix the annihilation cross section by the thermal freeze-out requirement (right panel of Fig. 8).
In this case, a limited parameter region is still allowed, but will be probed by future fixed-target experiments, as well as future direct detection experiments using electron recoils. Our result is in agreement with previous studies which generally fix m Z ≡ 3m φ , and α D = 0.1 [139] (or 0.5 [142][143][144]), but have not considered dark matter selfinteraction. For both panels of Fig. 8, the excluded region from Z searches would reduce by choosing larger values of α D . By introducing the upper bound on α D from Eq. (31), the exclusion in Fig. 8 is hence robust in both cases, i.e. either when imposing the DM interpretation of INTEGRAL (left) or when requiring successful thermal freeze-out (right).

VI. CONCLUSION
In this work we consider the possibility that DM is a complex scalar particle φ with a mass below the GeVscale. The particle is assumed to couple to SM leptons, either via a charged fermion F or via a vector boson Z . These models fare among the simplest UV-complete extensions to the SM, and have been contemplated as sub-GeV DM candidates well before the field exploded with activity in this mass bracket. Among other reasons, they draw their attention from the fact that φ annihilation today might explain the galactic INTEGRAL excess and/or bring into reconciliation the prediction and observation of the anomalous magnetic moment of the muon.
Given the tremendous recent activity devoted to the search of light new physics, it is only timely to revisit these models of scalar DM in light of much new data. These particles can be probed in the laboratory such as in electron-beam experiments, and by astrophysical and cosmological observations. We collate the latest observational and experimental data and subject the model to all relevant bounds and provide forecasts on the sensitivity of proposed future experiments.
Respecting the bounds on charged particles from high energy colliders LEP and LHC, we consider F to be at or above the EW scale. The combination m F /c 2 F is inherent to most observables and can be interpreted as the effective UV-scale Λ F for that model. We calculate the production of φ-pairs, mediated by the exchange of off-shell F , in the fixed target experiments NA64, mQ, LDMX, and BDX as well in e + e − colliders BaBar and Belle II. When the production is kinematically unsuppressed, the best bound is Λ F 250 GeV by BaBar, currently surpassed by LEP with Λ F 1 TeV. LDMX-II can improve on this number to 5 TeV. Turning to the Z model, we consider both heavy and light vector mediators. If Z remains off-shell in all experiments, we may take the combination m Z / √ g φ g l as the effective UV scale Λ Z . In this case, BaBar points to Λ Z 35 GeV to be improved by Belle-II to Λ Z 170 GeV at best, weaker than the current LEP bound of Λ Z 346 GeV.
These direct limits are then compared to loop-induced precision observables, concretely, to g − 2, to the invisible width of the Z and to Z-boson oblique corrections. We explicitly revisit all those calculations, confirming previously presented scaling relations in the limit m φ,l /m F 1 or m φ,l /m Z 1, and, as an added value, provide the full expressions of the loop integrals. We find that for the F model, the improved limit obtained from g − 2 of the electron surpasses all direct observables, with Λ F 10 4 TeV, while for the Z model, they do not play a role in the phenomenology. We also complement those constraints with limits that arise from the freedom in the chiral structure of the models, using the parity asymmetry in polarized electron scattering. Finally, we derive limits from lepton flavor violation that are dependent on the concrete UV-content of the models.
Turning to astrophysical constraints, we derive the anomalous energy loss induced by φ-pair production in the assumed proto-neutron star of SN1987A. This adds strong and complementary new limits on the parameter space for m φ 100 MeV down to Λ F 10 5 TeV and Λ Z 3 TeV. We furthermore consider constraints from direct detection, structure formation, CMB energy injection, and DM-self scattering. Here, the CMB puts stringent constraints on the s-wave annihilation mediated by F . In turn, for the p-wave annihilation mediated by Z the bounds are sub-leading. For those reasons, a thermal freeze-out in the F -mediated model is firmly excluded, whereas the Z model remains little constrained from energy injection.
Regarding the DM interpretation of the INTEGRAL 511 keV line, we show that it is excluded in both the F -mediated model as well as in the Z -mediated model with m Z ≥ 10 GeV. In the model with charged F , the crucial constraints come from the (g − 2) e data, from the CMB, and from SN1987A. For the Z model, intensityfrontier experiments and direct detection via electron recoils play the major role. However, a caveat exists: if the annihilation is resonant, m Z 2m φ , the INTEGRAL signal may still be explained in conjunction with a light m Z ≤ 10 GeV while at the same time being experimentally allowed.
In a final part, we then entertain the possibility of a Zmass below 10 GeV. The phenomenology then changes, as the vector may go on-shell in various considered searches. To derive conservative constraints on the Z coupling to the observable sector, g l , we then saturate the dark coupling g φ , taking its maximally allowed value from DM self-scattering. Requiring a successful explanation of INTEGRAL by fixing m Z constrains the values of the remaining free parameters g l and m φ , driving m φ to the resonant annihilation regime; the combination of LDMX, Belle II, and direct detection via electron recoils will be able to completely cover the parameter region of m Z ≥ 2.1 m φ in the near future, i.e. the range where Z remains off-shell.
As an outlook, we comment on the resonant region which is not studied here. For m Z 2m φ the annihilation cross section is greatly enhanced and the required value on g l coming from the annihilation cross section diminishes. This hampers the direct experimental sensitivity considered in this work. In turn, however, it opens the possibility of using displaced vertex searches in fixedtarget experiments, depending on the decay mode of Z . Dialing down the Z -mass further, m Z < 2m φ the annihilation via φφ * → Z ( * ) Z ( * ) → 2e + 2e − will eventually come to dominate. As the process is not velocity suppressed, we then re-enter the regime where stringent CMB bounds apply. We leave those aspects to dedicated future work.
Acknowledgments XC and JP are supported by the New Frontiers program of the Austrian Academy of Sciences. JLK is supported by the Austrian Science Fund FWF under the Doctoral Program W1252-N27 Particles and Interactions. We thank the Erwin Schrödinger International Institute for Mathematics and Physics where this work was started for their hospitality. We acknowl-edge the use of computer packages for algebraic calculations [145,146].
Appendix A: Full expression for (g − 2)l From explicit calculation for the (g − 2) l contribution in the representative models as shown in Fig. 3 we obtain, Taking only the leading order terms and setting g L = g R ≡ g l , we recover Eq. (4) and Eq. (8) of the main text.
Appendix B: Formulae for cross sections

φ-pair annihilation
Here we give the non-relativistic expansion of the DM annihilation cross section via a t-channel fermion F , to second order of the relative velocity, assuming the hierarchy m F m φ > m l , Here we have neglected terms of O(m 5 l,φ /m 5 F ). By taking only the leading order in O(m l,φ /m F ) and replacing v rel with 2v φ , we retrieve the s-wave expression of Eq. (1) in [4], while our p-wave result differs. The difference arises from the general mismatch between v rel and v M . We prefer to use the non-relativistic expansion of the Lorentz invariant quantity σ ann v M [147]; our scattering amplitude agrees with the one given in the appendix of [4]. For the special case c L c R = 0 and for m l = 0, the s-wave component in (B1) vanishes and the process becomes p-wave dominated, scaling as v 2 rel m 2 φ /m 4 F . For real scalar DM, φ = φ * , the annihilation process via a u-channel fermion F needs to be added, and the cross section in the same approximation becomes 2. e − e + annihilation with initial state radiation First, we detail the annihilation cross section associated with initial state radiation (ISR), corresponding to Fig. 4. Following [60], the differential cross section with ISR is formulated as the cross section without ISR times the improved Altarelli-Parisi radiator function. The an-nihilation cross sections, without ISR and with s denoting the squared CM energy, to order O(m 3 e,φ /m 3 F ) and after an average over the initial state spins has been performed, reads For the fixed-target experiments, we consider the production of φ-pairs via electron-nucleus bremsstrahlung depicted in Fig. 5. Here, we need to compute the 2-to-4 cross section σ 2→4 for the process e − (p 1 ) 2 ≡ t 2 such that q = q 1 + q 2 given that total momentum is conserved. The differential cross section is then written as where g 1 and g 2 are the spin degrees of freedom of electron and nucleus, |M| 2 is the amplitude square, and dΦ is the total phase space. By introducing an integral w.r.t. s X ≡ m 2 X = p 2 4 , in the lab frame Eq. (B5) becomes where L µν,ρσ contains the leptonic average over g 1 , φ ρσ includes the φ-emission piece together with the heavy mediator propagator, W µν is the hadronic tensor with its concrete form given in [60] and dΦ 4 is the 4-body phase space which is also analytically computed in [60].
The leptonic tensor from the two diagrams and their interference can be expressed as where, for completeness, we spell out the individual traces in the following. For the F -mediated model they read, where indices ρ, σ can obviously be abandoned. For the vector-mediated model we obtain for the traces, L µν,ρσ ab,Z = 1 g 1 Tr ( / p 3 + m e )γ ρ (g L P L + g R P R )( / p 3 + / q + m e )γ µ ( / p 1 + m e )(g L P R + g R P L )γ σ ( / p 1 − / q + m e )γ ν + Tr ( / p 3 + m e )γ µ ( / p 1 − / q + m e )γ ρ (g L P L + g R P R )( / p 1 + m e )γ ν ( / p 3 + / q + m e )(g L P R + g R P L )γ σ . (B9) In the case that φ-emission proceeds by F mediation, we may simplify the calculation by utilizing the fact that m F = 1 TeV is much larger than any momentum transfer considered here and approximate the the F propagator as where q is the four-momentum flowing in the propagator. We hence write the llφφ interaction as an effective operator suppressed by m F . The φ-emission piece is thus φ (F ) ρσ = 1, and the integral over the product of Lorentzinvariant phase space of φ and φ * denoted by dΠ φ,φ * , is given by Note that the Lorentz indices ρ, σ do not have physical meaning in the F -mediated model, and are introduced to unify the notation of the two models.
The φ-emission piece φ (Z ) ρσ in the Z case reads Terms containing p φ and p φ * can be integrated over the φ-pair phase space using a modified version of Lenard's formula [134] for massive final states, yielding Here, the function f (s φφ ) following the convention used in [60,134,148] and for the Z -mediated model reads When considering the experiments NA64 and LDMX where final state electrons are measured, we utilize the following differential cross section e ds X ds φφ dt 2 dp 1q ∂φ R4q e ds X ds φφ dt 2 dp 1q ∂φ R4q where p 1q ≡ p 1 · q, the angle φ R4q 4 is the angle between ( q 1 , p 1 ) plane and ( q 1 , q) plane in the frame that p 4 + q = 0, s 4 ≡ (p 4 + q) 2 = (p 1 + p 2 − p 3 ) 2 , λ(a, b, c) = a 2 + b 2 + c 2 − 2ab − 2bc − 2ac is the Källén function, and E 1 = E beam in the lab frame. The integration boundaries for s φφ , t 2 and p 1q , and the Jacobian |∂φ R4q 4 /∂p 1q | are found in [60,134]. In the lab frame (| p 2 | = 0), the Lorentzinvariant variables s 4 and t 1 can be expressed in terms of E 1 , E 3 and cos θ 3 as For the experiments mQ and BDX we need the spectrum and distribution of the produced φ particles. Therefore, we use where s 36 ≡ (p 3 + p φ * ) 2 , φ * 436 is the angle between the planes of ( p 1 − p φ , p 1 ) and ( p 1 − p φ , p 4 ) in the frame that p 3 + p 4 + p φ * = 0 ranging from 0 to 2π, t 23 ≡ (p 2 − p 3 ) 2 , φ * 36 is the angle between the planes of ( p 1 − p 4 − p φ , p 4 ) and ( p 1 − p 4 − p φ , p 3 ) in the frame that p 3 + p φ * = 0 ranging from 0 to 2π, s 436 ≡ (p 3 + p 4 + p φ * ) 2 = (p 1 + p 2 − p φ ) 2 , and t 15 ≡ (p 1 −p φ ) 2 . The integration boundaries of these variables are given in [60]. Note that In the lab frame, one obtains the relations In the derivation of the squared amplitude, one needs to define an additional Lorentz-invariant variable p 65 ≡ p φ * · p φ so that every scalar product of momenta can be written in terms of a combination of Lorentz-invariant variables [60]. Using the fact that in 4-dimensional spacetime, any five 4-vectors cannot be linearly independent, we express p 65 as a function of other Lorentz-invariant variables by solving det(M ) = 0, where the (i, j) entry of the 5 × 5 matrix M is p i · p j . There are two solutions of p 65 corresponding to φ * 36 ∈ [0, π) and φ * 36 ∈ [π, 2π). Other Lorentz-invariant variables are not affected by φ * 36 → 2π − φ * 36 .

φ-e scattering
To leading order, the differential recoil cross section dσ φ-e /dE R as a function of the incoming φ-energy, E φ , and the recoil energy of the electron, E R , for the F (sand u-channel scattering) and Z (t-channel scattering) cases read as follows, where we use g L = g R ≡ g l for the sake of presenting a more compact formula. The equations above are used γ, Z(p) in computing the electron scattering signal in mQ and BDX.
To properly account for the Pauli blocking factor in the computation of the upper boundaries of the SN1987A exclusion region, we also provide the differential cross sections in terms of Mandelstam variable t, taking again the limit of heavy mediators, The presence of heavy new charged fermions F ± , induces effective couplings between φ and SM gauge bosons through a set of triangle loops containing F ± and SM charged leptons, demonstrated in Fig. 9. This loopinduced charge radius interaction can e.g. be probed in direct detection experiments [39,149].
Here we detail the calculation of the coupling of the φpair to an off-shell photon via the aforementioned triangle loops. The amplitude that needs to be dotted into the off-shell photon reads where iM l 1F,µ and iM l 2F,µ correspond to the diagrams containing one and two F ± in the loop, respectively. Note that we take one charged SM lepton as an exam-ple, and denote its mass by m l . Here left-and righthanded SM leptons contribute equally, so we further set c L = c R ≡ c F . After using Feynman parametrization and dimensional regularization to perform the loop integral, we find that the divergences in M l 1F,µ and M l 2F,µ mutually cancel. The remaining finite terms read with p µ being the four-momentum of the photon, p 2,µ the four-momentum of one of the outgoing φ particles, and ∆ 1F =−m 2 l (z−1)+z m 2 F +m 2 φ (z−1) +p 2 x(x+z−1), ∆ 2F =−m 2 F (z−1)+z m 2 l +m 2 φ (z−1) +p 2 x(x+z−1).
Note that if m l = m F , the total amplitude vanishes, as a manifestation of Furry's theorem. One can also directly write the combination of the two diagrams as an effective charge radius operator if we integrate out the heavy F , where F µν is the field strength tensor of SM photon. We have checked that our loop calculation correctly matches onto this effective operator, with the Wilson coefficient in numerical agreement with previous results [39,149]. Note that a sum over the contributions from all SM leptons needs to be performed in the actual evaluation. The above calculation can be generalized to infer the additional contribution to the invisible Z decay width by replacing the relevant couplings in the above amplitudes by the weak charges where s W = sin θ W , c W = cos θ W with θ W being the weak angle, and g V being the usual vector coupling of weak current; axial vector currents do not contribute when c L = c R . Note that one also needs to include the contribution from diagrams containing neutral lep-tons as shown in the bottom panel of Fig. 9. We have checked that the resulting bound on c F is rather weak, and is thus not shown in the constraint plots.
In the Z case, the SM photon or Z boson gains an effective coupling to φ via mixing with Z , originating from a SM lepton loop, shown in the bottom panel of Fig. 9. To estimate the mixing, we need to compute the self-energy iΠ µν given by the usual Lorentz structure, iΠ µν = i p 2 g µν − p µ p ν Π(p 2 ) .
The divergence in Π(p 2 ) can be cancelled by a counterterm, that is, after specifying a renormalization condition, the effective mixing can be evaluated. The mixing between Z and Z is computed in the same way but with Eq. (C6). The ensuing constraint on √ g φ g l from the Z invisible decay is fully covered by the LEP bound.