Constraints on Light Leptophilic Dark Matter Mediators from Decay Experiments

We study the influence of leptophilic dark matter interactions on decays of muons and ground state mesons in existing experiments. We consider a secluded dark sector exclusively interacting with leptons via either a (leptophilic) scalar or vector mediator. These interactions will therefore influence leptonic decays and deform the energy spectra. We first study the Michel decay of muons, $\mu^+\rightarrow e^+\nu_e \bar{\nu}_\mu$, which allow us to constrain the parameter space reasonably well. Secondly, the rare $\pi^\pm$, $K^\pm$, $D^\pm$ and $D_s^\pm$ decays to $e\nu$ will be considered. Scalar mediators would remove the Standard Model helicity suppression, so that strong constraints can be derived. The resulting bounds on the couplings of the light mediators to electrons and muons still turn out to be somewhat weaker than those from searches at low--energy $e^+e^-$ colliders and the magnetic moment of the muon, respectively. Finally, we show that kaon and pion decays basically exclude a ''Co--SIMP'' scenario where a scalar dark matter particle has a dimension--5 coupling to electrons.


Introduction
Evidence for the existence of substantial non-baryonic mass in the universe, in addition to the baryonic contribution from the known Standard Model (SM) particles, has been piling up for decades [1]. Cosmological observations, from the CMB at the largest scales, the structure of galaxy clusters and gravitational lensing on intermediate scales, down to the rotation of single galaxies, all essentially only probe the gravitational interactions of this dark matter (DM), leaving the properties of the constituents of DM largely obscure. We do know that DM should be "cold", i.e. non-relativistic well before the CMB decoupled. Moreover, within the minimal cosmological framework, the overall DM density can be determined accurately, Ω DM h 2 = 0.120 ± 0.001 [2]; here Ω DM is the scaled DM mass density and h is the rescaled Hubble parameter.
In the absence of data pinning down the properties of DM, a wide range of models has been proposed. In spite of intensive efforts no clear signal for DM particles has yet been found in either direct or indirect detection experiments [1], leading to severe constraints on many models [3]. This is true in particular for models with weakly interacting massive particles (WIMPs), with masses very roughly at the weak scale. This has led to increased interest in sub-GeV masses [4], which for purely kinematic reasons are much less constrained.
Here we consider models that couple the potential dark sector to the Standard Model with light leptophilic mediators, i.e. mediators that couple directly only to leptons. One class of models assumes a gauged lepton-family number; after spontaneous symmetry breaking these models contain a massive vector mediator, sometimes called a "dark photon", that couples to some leptons of the Standard Model and in the dark sector to Dirac dark matter [5,6]. This lepton coupling introduces a kinetic mixing term with the ordinary photon, resulting in a small coupling between the dark photon and all electrically charged particles which can also be used to constrain the model. Another class of models assumes a scalar mediator, which again only couples directly to some or all charged leptons. Since this breaks the SU (2) gauge symmetry, models of this kind can at most be an effective theory. Finally, we consider the so-called "Co-SIMP" model [7] containing a light scalar DM particle with a non-renormalizable coupling to electrons.
We present a novel approach to constraining the parameter space using measurements of the decays of muons or ground-state flavored mesons into final states containing an electron. The spectrum of electrons produced in muon decays has been measured accurately; it agrees with SM predictions, which allows us to put upper bounds on the coupling of spin−1 mediators to electrons or muons. However, the resulting bounds turn out to be more than one order of magnitude weaker than the best constraint from e + e − colliders. New spin−0 particles coupling to electrons would remove the helicity suppression in charged meson decays into eν e final states; the resulting bounds on the renormalizable couplings of light scalar mediators are tighter than those from muon decay, but still somewhat weaker than those from e + e − colliders. However, bounds from pion and kaon decays suffice to exclude a thermal "Co-SIMP" for masses below 0.8 MeV; in the allowed mass range χ does not behave like a SIMP any more.
The remainder of this article is structured as follows: In section 2 we present the leptophilic models considered here. In sections 3 and 4 we describe the method of obtaining limits on the parameter space from the Michel decays of muons and pseudoscalar meson decays respectively. Section 5 presents our results and compares them to existing bounds. Section 6 finishes with some concluding remarks.

Models
The strongest bound on many DM models comes from "direct" search experiments, which look for elastic scattering of ambient DM particles off the nuclei in a detector [1]. Leptophilic dark matter models, where the dark matter particles primarily couple to the Standard Model leptons, either directly or via another "mediator" particle, avoid most of these bounds. In such models the DM particles can interact with nucleons only via loop diagrams.
One way of incorporating these ideas is a hidden sector that contains only singlets under the Standard Model gauge group. However, the simplest (thermal) DM production mechanism requires some coupling to SM particles. To this end one may introduce additional fields which mediate interactions between both sectors. These are then fittingly called portals.

Vector Mediator
The extension of the Standard Model by a new vector boson is well motivated both from a bottom up as well as from a top down perspective, e.g. from grand unified theories [8]. Here we consider scenarios where the gauge group is extended by another U (1) D gauge group which is spontaneously broken such that the associated particles become massive. The new vector boson A is often called a dark photon. By assumption the Dark Matter particles χ are charged under U (1) D . In this article we are concerned with the production of a light A which decays invisibly. The exact nature of χ therefore is not relevant for us. Assuming it to be a Dirac fermion for simplicity's sake, the resulting Lagrangian can the be written as Here F µν and F µν are the field strength tensors of QED and of U (1) D , respectively. The new coupling constants g D and e l are free parameters of the model, and the term ∝ describes kinetic mixing between the new gauge boson and the photon. In order to avoid anomalies, one may chose to gauge any combination X = yB − x i L i of baryon number B and lepton family number L i , with constraint 3y = x e + x µ + x τ [9]. Popular choices include gauged B − L or L µ − L τ [10,11]. While the former is highly constrained by direct collider searches and other 5th force experiments, the latter still exhibits rather weak constraints with a region that is even favored by the g µ − 2 anomaly [12][13][14]. The term proportional to has been added to work with the most general renormalizable gauge invariant Lagrangian. Even though this term might be absent at tree level, as is the case for some GUT theories, it can be generated by loop contributions such as the diagram shown in fig. 1, where the fermion running in the loop is charged both under the Standard Model and dark U (1)-group. In the case at hand this leads to = l ee l 12π 2 ln m 2 l µ 2 . (2.2) In anomaly-free theories the dependence on the renormalization scale µ cancels in the sum. The term ∝ in eq.(2.1) leads to additional effective interactions of the form where J µ em = ψ Q ψψ γ µ ψ is the electromagnetic current; hence every electrically charged particle interacts with the dark photon as it is now millicharged under U (1) D [15]. This can be used to put strong constraints on the parameter space.
If m A < 2m χ the new gauge boson will mostly decay to the leptons to which it couples directly: where d l = 1 for charged leptons while d l = 1/2 for left-handed neutrinos. If m A > 2m l ± one may search for visible A → l + l − decays. Here we are instead interested in scenarios with mostly invisible A decays, either because m A < 2m l ± for the relevant charged lepton, or because m A > 2m χ and g D e l . A → χχ decays are also described by eq.(2.4), with the obvious replacements e l → g D , m l → m χ and d l → 1.
Within a given cosmological scenario the χ relic density imposes one constraint on the parameters of the model [4,14]; even the case e l = 0, in which case A couples to SM particles only via eq.(2.3), can lead to the correct relic density in minimal cosmology [16].
Here we implicitly assume that this constraint is used to determine the DM mass m χ , which allows us to vary the mass and couplings of A freely.
Additional motivation for direct interactions with the muon specifically come from the anomalous magnetic moment of the muon, a µ = (g −2) µ /2. The Standard Model prediction [17] differs from the experimental result [18]: The additional 1−loop contribution from a vector boson coupling to muons is [19,20]: (2.6)

Scalar mediator
Another possibility is that a scalar particle mediates interactions between the dark matter and the SM particles. Here we consider a real scalar field φ with mass m φ that is a singlet under the Standard Model gauge group, and again a dark matter Dirac fermion χ. The Lagrangian is: This Lagrangian respects QED gauge invariance, but the interactions of the scalar with the leptons break electroweak gauge invariance explicitly. These couplings might originate from gauge invariant (but non-renormalizable) dimension 5 operators [21]: where Φ is the Standard Model Higgs field. Once Φ obtains a vacuum expectation value v, the interactions in (2.8) lead to the scalar couplings in eq.(2.7), with Another possibility is to have the light scalar φ mix with the neutral component of a (second) scalar doublet, which can in principle have renormalizable O(1) couplings to leptons. 1 Most lepton-specific scalar mediator models considered in the literature assume the effective scalar couplings c l to be proportional to the charged lepton masses, in which case the e l follow the lepton mass hierarchy.
A nonvanishing e l leads at 1−loop to an effective coupling of the scalar to two photons, through the diagram shown in figure 2. If m φ < 2m l , this decay can be used to search for φ in the diphoton invariant mass distribution. The decay width to photons is 2 [24] and and the loop function F 1/2 reads For g D > e l and 2m χ < m φ the scalar mediator decays mostly invisibly to the dark sector. The corresponding decay width is The main avenue for collider experiments to constrain this model is then missing energy searches. On the other hand, for 2m χ > m φ > 2m l the mediator decays mostly visibly to leptons it directly couples to, with decay width We'll be interested in scenarios where |e l | > ∼ 10 −3 for l = e or µ. The decay width (2.13) then corresponds to a lifetime τ φ < ∼ 1.5 · 10 −14 (1 MeV)/m φ seconds; even accounting for a Lorentz boost, φ → l + l − decays will then usually be "prompt" if the corresponding decay is kinematically allowed.

Co-SIMP
Finally, we consider the Co-SIMP mechanism proposed by Smirnov et al. [7]. A real scalar particle χ with strong self-interactions is assumed as dark matter. An interaction with the Standard Model of the form χχe → χe is introduced in order to dissipate entropy from the dark sector whilst a Z 3 symmetry stabilizes χ. We consider an electrophilic version of this model, in which case the relevant interaction is described by the effective operator Once again this does not respect the electroweak gauge symmetry; as before gauge invariance can be restored by replacing the dim−5 operator of eq.(2.16) by a dim−6 operator ∝ēΦeχ 3 + h.c.. Observation of cosmological structures imply m χ > ∼ 5 keV, while m χ < ∼ m e is required in order to avoid a WIMP-like freeze-out [7]. Note that this model does not contain additional free parameters that allow to tune the relic density independent of the laboratory limits which are the main topic of this work. In order to compute the cosmologically preferred value of Λ we therefore solve numerically the Boltzmann equation describing the freeze-out of χ [7]: Here x = m χ /T , T being the temperature, s is the total entropy density, and Y ψ = n ψ /s where n ψ is the number density of ψ particles, with ψ ∈ {χ, e} in our case; the superscript eq denotes the equilibrium value of the corresponding quantity. The thermally averaged cross section [25] is also obtained numerically from the integral Here g i denotes the number of internal degrees of freedom of particle i (1 for χ and 4 for e), f i = 1/(exp E i /T ± 1) is the phase space distribution function for the i−th particle in the initial state, and |M| 2 is the averaged squared matrix element for the relevant process χ + χ + e → χ + e.

Bounds from µ − Decays
Having introduced the models we will consider, we turn to a discussion of the decays which we use to derive bounds on the parameters of these models. We begin with a discussion of the muon decay spectrum. Since this has widely been used as a high precision test of the electroweak theory [26], it might provide a good chance to constrain our models. The double differential width for µ − → e − ν µνe decays at rest can be written as [1]: Here W eµ = (m 2 µ + m 2 e )/2m µ is the maximum electron energy (neglecting possible neutrino masses), x = E e /W eµ is the rescaled electron energy, x 0 = m e /W eµ is its minimum value, P µ is the degree of muon polarization and θ is the angle between the polarization vector of the muon and the outgoing electron. The functions for the isotropic part F IS (x) and the anisotropic part F AS (x) are given by [1]: Here the Michel parameters ρ, η, ξ and δ have Standard Model values 3/4, 0, 1 and 3/4 respectively. Measurements of these parameters have been used to put constraints on the effective parameters of additional four fermions interactions.
Here we use precise measurements of the muon decay spectrum to put constraints on four-body decays µ − → e − ν µνe X where X is a light leptophilic mediator. The contributing Feynman diagrams for the case of a vector mediator coupling to electron number are shown in figure 4. We assume that X is either long lived or decays invisibly, so that the four-body final state has the same basic signature as the three-body final state. However, the observable electron spectrum of this four-body mode differs from the spectrum predicted by the SM. The obvious Standard Model background is the radiative muon decay, µ − → e − γν µνe , where the photon escapes detection but nevertheless carries some energy, thereby also altering the electron spectrum. Figure 4. Diagrams contributing to µ − → e − A ν µνe decays for the case that A couples to electron number.
Since our signal involves a four particle final state, this case cannot be mapped directly onto eqs.(3.1) and (3.2). However, when the dark mediator remains invisible the observable final state has the same topology as in the SM, i.e. the combination will be observed. In order to derive estimates of possible constraints we studied the sensitivity of fitting procedures similar to those employed in the experimental determination of the spectral parameters.
The four body matrix elements were derived with help of the Mathematica package FeynCalc [27] followed by the numerical integration over all kinematic parameters except x and cos θ. The resulting two dimensional distribution is then added to the pure Standard Model spectrum. The Michel parameters that best describe this new distribution were extracted from eqs.(3.1) and (3.2) by means of a χ 2 fit, using the same binning as in ref.
[26]. More exactly, we minimized Here Γ i is the muon decay width in the i−th bin; Γ i (Michel) is computed from eqs. (3.5) We assume muon polarization P µ = 1, as predicted by the SM for the relevant case of muons produced in meson decays.
A technical subtlety arises because experiments do not fit to the whole spectrum: cuts have to be applied in order to cover detector inefficiencies and blind spots. We modeled the effects of these cuts by restricting our kinematical fit to the fiducial region covered by the TWIST detector [26]: approximately x ∈ [0.45, 0.98] and |cos θ| ∈ [0.54, 0.96]. We found that these cuts affect the sensitivity limit on the coupling only by an O(1) factor, for the mediator masses considered here. Our sensitivity limits should not be confused with experimental bounds; we did not use real data, nor did we include QED corrections when modeling the SM prediction for the decay spectrum. Our procedure should nevertheless give a reasonable estimate of the sensitivity of the measurements of muon decays to the new mediators.
We finally note that the measurement of the muon lifetime cannot be used directly to constrain our models. In the SM this measurement is used to determine the experimental value of G F , which is a free parameter of the theory. A deviation from the SM could therefore only be detected by comparing this measurement with a second, independent determination of G F . Assuming unitarity of the quark mixing (CKM) matrix, the experimental "CKM unitarity test" can be recast as a measurement of G F -with, however, much poorer precision [1]. Moreover, the emission of collinear mediators can give rise to ln(m µ /m X ) enhanced terms in the decay distribution, which cancel in the total muon decay width by the Kinoshita-Lee-Nauenberg (KLN) theorem [28,29] once loop diagrams are included. We therefore expect that a comparison of different measurements of G F has much poorer sensitivity to the light mediators we consider than the measurement of the muon decay spectrum discussed above.

Bounds from Leptonic Decays of Charged Pseudoscalar Mesons
At tree level the total width for the decay of a charged pseudoscalar meson P ± to a lepton pair is given by [30] Here G F is the Fermi constant, m P and m l are the meson and lepton mass, respectively, F P is the P decay constant, and V q 1 q 2 is a CKM matrix element, P + being a (q 1q2 ) bound state. Owing to the V − A nature of charged current weak interactions, both the charged lepton and the neutrino "like" to be left-handed, which however is forbidden by angular momentum conservation. This leads to the well known helicity suppression represented by the factor m 2 l in eq.(4.1). Evidently this factor strongly suppresses the decay to an electron and neutrino. This is often exploited for tests of lepton universality in the ratio of decays to electrons and muons. The SM prediction for the ratio of decay widths is where δR QED describes the effect of QED corrections, including real photon emission. It is important to note that the emission of a spin−1 boson does not change the helicity structure of the amplitude, and therefore does not lift the helicity suppression. In contrast, if a scalar or pseudoscalar particle couples to the electron, the helicity suppression is removed, which can enhance the electronic decay mode significantly. Since measurements are in agreement with the SM prediction (4.2), this lifting of the helicity suppression can be used to derive bounds on the couplings of new light spin−0 particles. In the limit of vanishing electron and neutrino masses, the total width for P → eν e φ decays can easily be computed analytically: 3 This expression manifestly avoids the m 2 e suppression. In the Co-SIMP model one instead has to emit three scalar χ particles, leading to a considerably more complicated phase space integral; however, since the new vertex again violates chirality, also in this case the helicity suppression is lifted.
Whenever the new scalars remain invisible the event will have the same topology as a rare decay into e + ν e , albeit with a softer electron energy spectrum. This change of the electron spectrum can reduce the sensitivity due to kinematic cuts employed by the experiments. Bounds are then extracted by saturating the maximal allowed difference ∆ P between the theory prediction R SM p and the experimental result R exp p , i.e. ∆ P ≥ Γ(e + ν e + X) Γ(µν µ (γ)) .

(4.4)
Here X stands for either a single spin−0 mediator φ or for the three Co-SIMP scalars χ, and is an acceptance correction factor due to the softer electron spectrum. For the most precise measurement of R π , by the PIENU collaboration [34], this factor is computed as follows. This experiment analyses decays of a stopped π + beam. In order to discriminate between direct π + → e + decays and the dominant background from π + → µ + → e + , an energy cut E e ≥ 52 MeV was used in the experimental definition of π + → e + ν e decays, which have a nominal positron energy of E e = 69.8 MeV. If we want to apply this analysis to our π + → e + νφ decay, the same cut on the positron energy should be applied. The resulting acceptance correction is shown in figure 5.  . Acceptance correction for π + → e + ν e φ and π + → e + ν e χχχ decays due to the cut on E e by the PIENU experiment.
For the Co-SIMP model the acceptance is approximately constant in the allowed mass range, with ≈ 11%. This limitation does not apply to K ± , D ± and D ± s decays, which are studied in flight, so that even the two-body decay mode has a broad energy spectrum. The most precise measurement of R K comes from the NA62 collaboration [35]. Here a Kaon beam decays in flight inside the detector. Since the accepted range of electron energies has a width of several 10 GeV we simply assume that all decay modes have the same acceptance.
The decay of D ± or D ± s mesons to eν has not yet been observed, but 90%c.l. upper bounds on the corresponding branching fractions have been set by the CLEO and Belle Collaboration respectively [36,37]. We use these to derive limits on the couplings of our light spin−0 mediators.

Results
Here we show our estimated sensitivity to the new coupling as function of the mass of the postulated new boson, derived from muon and charged meson decays. We also show existing limits found in the literature, where we focus on the strongest bounds for the case at hand.
The relevant current bounds come from: • BaBar: Dedicated search for the dark photon in e + e − → γA ; A → invisible with CM energies near the Υ resonances. 90% c.l. limits were derived on the dark photon coupling constant ε 2 for m A ≤ 8 GeV [38]. This can directly be applied to our model with a vector mediator, with the replacement e e = eε [6].
• CCFR: Measurement of neutrino trident production events νN → νN µ − µ + using a muon-neutrino beam with average energy E ν = 160 GeV. The observed N CCFR = 37.0 ± 12.4 events agree with the Standard Model prediction N SM = 45.3 ± 2.3 [39]. This is used to set limits on additional contributions from a vector mediator coupling to the muon-neutrino and to the muon [9].
• Belle II: Search for e + e − → µ + µ − Z with beam energies of 4 and 7 GeV where the Z is radiated off of one of the muons and decays invisibly with m Z < 6 GeV [40]. The limits can then be recast to muonphilic scalar mediators [41].
• PIENU: Search for the three body decay π + → l + νX where l is an electron or muon and X an invisible neutral boson. Pions were stopped in a detector and the spectrum of the charged lepton was measured. A search for the smooth signal spectrum was then carried out below the energy of the two body decay. This sets limits on the branching ratio Γ (π + → e + νX) /Γ (π + → µ + ν) for X−masses in the range 0 < m X < 120 MeV [42]. . The blue curves indicate our estimated 90% c.l. sensitivity limit on the coupling constant between electron or muon and the spin−1 mediator A from existing measurements of muon decay. The grey shaded region in the left frame is excluded by the BaBar search for e + e − → γA (A → invisible) [38], rescaled to e e = εe, while the red shaded region in the right frame is excluded by an analysis [9] of CCFR data on νN → νN µ − µ + . The green band in the right frame indicates parameters that bring the experimental and theoretical values of (g − 2) µ within 2σ.
Our estimated sensitivity of existing muon decay data to the couplings of a new vector mediator are shown in Fig. 6. Unfortunately these estimated sensitivities are considerably weaker than the best existing bounds; the electron coupling is constrained by BaBar, while the muon coupling is constrained by the CCFR trident data. We find a considerably better sensitivity to the electron coupling, since near-collinear A emission off the electron, which is enhanced by a large logarithm, reduces the energy of the electron and thus leads to an observable effect. In contrast, near-collinear A emission off the muon neutrino leaves the electron energy essentially unchanged. Below a mediator mass of O(10 MeV) the bound on the Michel parameter δ determines the sensitivity limit. For larger mediator masses the sensitivity limit is set by the ξ parameter. Of course, the sensitivity is worse for larger mediator masses due to the closing phase space. In the left frame, the dashed black, and solid yellow, red and orange lines indicate limits from additional contributions to P + → e + ν e decays (dashed black: P = π; yellow: P = K; orange: P = D; red: P = D s ). The dot-dashed black line represents the limits derived from the result of the search of PIENU for π + → eνX; the drop of sensitivity at m φ ∼ 55 MeV is due to the signal being similar to π + → µ + → e + [42]. The region shaded in gray in the left frame is excluded by an older BaBar limit recast to scalar mediators [43]. In the right frame the orange shaded region is excluded by a recast of a Belle II search for e + e − → µ + µ − A , (A → invisible) [41], and the green band indicates parameters that bring the experimental and theoretical values of (g − 2) µ within 2σ.
Our projected bounds on the couplings of the scalar mediator φ are depicted in Fig. 7. Muon decay is much less sensitive to these couplings than to those of a spin−1 mediator shown in Fig. 6, since collinear emission of a soft spin−0 boson off a fermion is suppressed. This also explains why the bounds on e e and e µ from muon decay are now quite similar.
However, bounds on P + → e + ν e decays, where P + is a pseudoscalar ground state meson, lead to quite stringent bounds on new scalar couplings of the electron. For m φ ≤ 400 MeV the strongest bound originates from Kaon decays; decays of D + and D s mesons are considerably weaker, but extend to larger mediator masses. However, even the constraint from K + decays is somewhat weaker than that from an older BaBar search for e + e − → φγ with invisible φ. We do not show bounds from charged meson decays on the muon coupling. Since the helicity suppression of P + → µ + ν µ is much weaker than for the electron mode, the resulting bounds on e µ are considerably less stringent than those on e e , and are thus not competitive with the constraint from the measurement of g µ − 2. Finally, constraints on the leptonic Co-SIMP scenario are summarized in Fig. 8. The blue curve indicates the value of Λ required to obtain the correct relic density via χχe → χe reactions, which change the number of χ particles without changing the number of SM particles; this is the essence of the Co-SIMP scenario of ref. [7]. The constraint from K + → e + ν e decays excludes this value of Λ unless m χ > 0.8 MeV. Note that for m χ > 2m e /3 another reaction becomes possible, χχχ → e + e − ; however, this is not a SIMP scenario any more. Meson decays therefore exclude the whole region of parameter space where the correct χ relic density is determined by the Co-SIMP reaction. 4

Conclusion
In this paper we investigated the influence of leptophilic dark matter models on two well measured decays of Standard Model particles: deviations from the muon decay spectrum due to additional undetected particles and the removal of the helicity suppression in the leptonic decays of charged ground state pseudoscalar mesons. The measured muon decay spectrum proved to be less sensitive to the new couplings than published searches at e + e − colliders and from neutrino trident production. Limits on the electrophilic scalar coupling derived from meson decays are an order of magnitude stronger than those from muon decays but still weaker than existing limits by a factor of ∼ 2. Hence a moderate precision improvement in leptonic kaon decays could probe new regions of parameter space. On the other hand, the measurement of the muon decay spectrum would have to become a lot more accurate to yield competitive limits, which currently looks unlikely. The muonphilic scalar limit at m φ < 1 MeV is competitive with the recast search from Belle II at e µ ∼ 5 · 10 −2 , but lies well above the upper bound from the measurement of the magnetic dipole moment of the muon.
For the purely electrophilic Co-SIMP we showed that the model is ruled out for m χ < 800 keV by the kaon decay branching ratio because at allowed couplings χ is overproduced by freeze-out. Above this mass the DM particle χ doesn't behave like a SIMP any more, since its relic density is greatly affected by annihilation into SM particles (namely e + e − pairs). Our analysis therefore excludes the electrophilic Co-SIMP scenario.