Exploring New Physics with O(keV) Electron Recoils in Direct Detection Experiments

Motivated by the recent XENON1T results, we explore various new physics models that can be discovered through searches for electron recoils in O(keV)-threshold direct-detection experiments. First, we consider the absorption of axion-like particles, dark photons, and scalars, either as dark matter relics or being produced directly in the Sun. In the latter case, we find that keV mass bosons produced in the Sun provide an adequate fit to the data but are excluded by stellar cooling constraints. We address this tension by introducing a novel Chameleon-like axion model, which can explain the excess while evading the stellar bounds. We find that absorption of bosonic dark matter provides a viable explanation for the excess only if the dark matter is a dark photon or an axion. In the latter case, photophobic axion couplings are necessary to avoid X-ray constraints. Second, we analyze models of dark matter-electron scattering to determine which models might explain the excess. Standard scattering of dark matter with electrons is generically in conflict with data from lower-threshold experiments. Momentum-dependent interactions with a heavy mediator can fit the data with a GeV dark matter mass but are in tension with collider constraints. Next, we consider dark matter consisting of two (or more) states that have a small mass splitting. The exothermic (down)scattering of the heavier state to the lighter state can fit the data for keV mass splittings. Finally, we consider a subcomponent of dark matter that is accelerated by scattering off cosmic rays, finding that dark matter interacting though an O(100 keV)-mass mediator can fit the data. The cross sections required in this scenario are, however, generically challenged by complementary probes of the light mediator. Throughout our study, we implement an unbinned Monte Carlo analysis and use an improved energy reconstruction of the XENON1T events. ar X iv :2 00 6. 14 52 1v 1 [ he pph ] 2 5 Ju n 20 20


Introduction
The quest to identify the particle nature of dark matter (DM) by detecting DM in terrestrial experiments has been ongoing for more than three decades. Despite numerous searches at direct-detection, indirect-detection, and collider experiments, no convincing signal for DM has been found to date. Given the profound implications for our understanding of the DM particle's properties if we were to find it in the laboratory, any claim for a possible DM signal in one of these experiments deserves to be studied carefully.
The XENON1T collaboration has recently observed an unexplained excess of electronic recoil events with an energy of O(keV) [1]. While the most likely explanation is a neglected background source or a statistical fluctuation, the possibility that the excess could be the first sign of new physics (not necessarily even a sign of DM) is intriguing. The excess of events does not appear in the traditional search for nuclear recoils from elastic DMnucleus scattering. Rather, it appears as an excess in a search for electron recoils (ER). The XENON1T search has an exposure of 2.36 × 10 5 kg − day in the 1 − 30 keV energy range. The background rate is reported to be 76 ± 2 events/tonne year keV implying a total of ∼ 1476 background events. An excess of 53 events has been observed at the 1 − 7 keV low energy region (corresponding to roughly a 3σ excess), with the excess mainly located in the 2-keV and 3-keV energy bins.
In this paper, we explore several possibilities for the origin of this signal. We will focus mostly on the possibility that the origin is attributable to DM, but will also consider bosonic particles (pseudo-scalar, scalar and vector) produced in the Sun, which do not necessarily have to be a DM component. We discuss in the context of the XENON1T excess several models previously considered in the literature: O(keV) bosonic DM that is absorbed by an electron in the xenon atom [2][3][4][5][6][7][8][9], bosonic DM that is emitted from the Sun [10][11][12][13][14], and DM scattering off electrons in xenon [15][16][17][18][19]. For the absorption of scalar particles, we explore systematically the implications of the XENON1T hint on realistic models of bosonic DM. We show that only the dark photon or a "photophobic" axion-like particle can fit the XENON1T hint. We also discuss new model scenarios that have not been previously considered in the literature: "exothermic" DM scattering off electrons (for previous work focused on nuclear scattering see [20,21]), a pseudo-scalar with a density-dependent potential to avoid stellar cooling bounds while explaining the XENON1T data (for a similar effect see e.g. [22][23][24][25][26][27][28][29]), and cosmic-ray accelerated DM that here interacts with electrons through an intermediatemass mediator (for previous work focused on heavy mediators or light mediators interacting with nuclei see [30][31][32][33]; see also [34]). These models deserve further study in future dedicated papers, but we provide their salient features focusing on the XENON1T excess.
This paper is organized as follows. In §2, we describe the requirements that new physics needs to satisfy in order to explain the XENON1T excess, and also detail the models that we will discuss. In §3, we describe important features of the XENON1T data, our method for reconstructing the energy, and our statistical analysis. In §4, we focus on the absorption of bosonic particles that are either (non-relativistic) DM particles in our halo or emitted from the Sun. §5 investigates how a density-dependent potential can be used to circumvent the stellar cooling bound. In §6, we discuss DM-electron scattering, reviewing the "standard" case and then focusing on multi-component DM with small mass splittings. We will see that "exothermic" DM scattering off electrons has a rich phenomenology. §7 considers a subdominant DM component that is accelerated by scattering off cosmic rays.

Models and Summary
The XENON1T excess motivates us to consider various known as well as novel new physics scenarios, focused mostly, but not solely, on DM models that can be discovered via a highthreshold ( keV) ER searches. We first summarize the relevant features of the excess and then identify possible mechanisms that may explain it.
The following considerations are important when considering a prospective new physics signal: • The excess events have an energy of 2-3 keV. The measured spectrum suggests that a potential signal should contribute to more than a single bin. However, the finite energy resolution of the experiment makes this (statistically weak) observation less sharp, allowing for rather narrow spectra to provide a reasonable fit. The 1 keV bin and the bins with energies ≥4 keV are approximately consistent with background expectations. See Fig. 4 left.
• Numerous direct-detection experiments place stringent constraints on any accompanying nuclear recoil signal (for a recent compilation of low-mass DM limits on nuclear interactions see [46]). Models that predict such a signal must evade them.
Prospective models that could produce the observed excess and satisfy its features can be separated into models that predict an absorption signal ( §4 and §5) and those that predict a scattering signal ( §6 and §7). We consider several scenarios: 1. Absorption. We will consider the case that an electron absorbs a bosonic particle: pseudo-scalar (axion), a scalar, or a vector. The boson may be either non-relativistic or relativistic. The former may occur if the particle constitutes a component of the DM; in this case the ER spectrum is peaked at the mass of the DM, and can fit the data only due to the experiment's finite energy resolution. We find that a vector and a pseudo-scalar can explain the XENON1T excess, while a scalar is in conflict with stellar cooling constraints. Next, light bosons may be produced in the Sun, which has a temperature of around 2 keV. A non-zero mass around 2.6 keV could also cut the solar emission kinematically, providing the best fit of the data. However, for bosons produced in the Sun strong constraints arise from stellar cooling, strongly disfavoring the couplings needed to explain the XENON1T excess for the vanilla axion and dark photon models.

2.
Chameleons. The stellar cooling constraints on light bosons may be evaded if the couplings of SM particles to the corresponding bosons are screened inside high-density or high-temperature stellar objects. Such chameleon-like particles have a rich phenomenology and can revive the Solar explanation of the XENON1T hint.
3. DM scattering. The DM-electron scattering rate depends on the momentumtransfer-dependent atomic form-factor. This steeply-falling function is highly suppressed for momenta q 1/a 0 = α EM m e , where a 0 is the Bohr radius, α EM is the fine structure constant and m e is the electron's mass. As a consequence, DM scattering through a light mediator or a velocity-independent heavy mediator predict a steeply rising spectrum at sub-keV energies and are thus disfavored. 4. Velocity-suppressed DM scattering. Models that exhibit velocity-or momentumdependent heavy-particle-mediated DM-electron scattering are allowed by experimental data at lower energies and provide an adequate fit to the XENON1T excess. However, such models are likely in tension with collider bounds on new particles that generate this operator [59]. 5. Exothermic DM. An unsuppressed high-energy spectrum from DM-electron scattering may stem from an exothermic scattering of DM off electrons, the result of DM consisting of two or more states whose masses are slightly split by an amount denoted as δ. The atomic form factor, together with the scattering kinematics, imply a rather narrow electron recoil spectrum that is peaked near δ for a wide DM mass range and can explain the XENON1T excess for δ ∼ O(keV). The spectrum can be broadened if the DM-electron interaction increases with increasing momentum transfer q, if there are three or more DM states whose mass is split by different amount of O(keV), or if the DM mass is well below the GeV-scale. We will discuss the phenomenology of the exothermic scenario, leaving to future work a detailed investigation in concrete models of the relic abundance of the heavier state(s) and the constraints. 6. Accelerated DM. A small subcomponent of DM may be accelerated through its interactions in the Sun [41,42] or with Cosmic Rays (CRs) [30][31][32]. While we find that the component accelerated from the Sun cannot explain the XENON1T excess without being in conflict with lower-threshold direct-detection searches, we find that CR scattering of DM with non-trivial momentum-dependent form factor can address the XENON1T excess while evading other direct-detection constraints. In the scenario we consider here, direct constraints on the mediator exclude robustly this explanation.
In Fig. 1, we summarize the goodness-of-fit of the various absorption scenarios discussed above to the XENON1T measurement. We see that the bosonic DM scenarios (red curve) can fit the data well with a predicted mass of m X = 2.6 keV and coupling to electrons. Among these, the scalar DM case is excluded by current direct detection constraints from XENON1T and PandaX [43,60] while the dark photon and the axion are good explanation of the XENON1T excess. In the latter case the anomalous axion coupling to photon should be set to zero to avoid X-rays constraints.
In all the solar cases, the adddition of a non zero mass ameliorates the fit by cutting off the spectrum kinematically, in better agreement with the 1 keV bin being consistent with the background prediction. For pure electron coupling the axion explanation is disfavored compared to the scalar or the dark photon. The reason is that the peaks in the spectrum of the axion ABC production [49] are not observed in the data. If the solar production happens through Primakoff process [57] the scalar provide a very good fit of the data while the axion explanation is disfavored. As we will discuss, the reason can be traced back to  the different energy dependence of the axion and scalar absorption rate in xenon. The scalar rate grows fast at low energies a good fit can be obtained by cutting the raise with a mass of m φ = 1.9 keV, hence generating a bump between 2 and 3 keV. Viceversa, the axion absorption rate is suppressed at low energies and the resulting spectrum is too flat at energies above 3 keV to provide a good fit of the signal, independently on the mass of the axion. Consistently, in the massive axion fit tends to prefer a massless axion.
In Fig. 2, we summarize the goodness-of-fit of the different scattering scenarios presented above. First, we notice that elastic scattering cannot explain the XENON1T hint for F (q) ∝ q n form factors with n ≤ 0, since the electron recoil spectrum rises at low energy, in tension with complementary direct-detection experiments at lower energy thresholds. However, for F (q) ∝ q n with n > 0 the spectrum falls fast enough towards lower energies and provides an adequate fit to the XENON1T data.
Second, we show that exothermic scattering can fit well the data when the heavy and light DM states are split in mass by a few keV. Once the splitting is marginalized to the best fit value, the p-value is essentially independent of the DM mass as long as it is heavier than the splitting itself. The spectrum is peaked near 2 keV and fits well the data, without being trivially excluded by complementary direct detection experiments. In concrete models, the rich phenomenology of these DM scenario could provide other handles of testing them at beam dump experiments or in nuclear recoil.
Third, we discuss accelerated DM by scattering with cosmic rays. In such a case, the challenge is again to find a scenario where the accelerated spectrum falls sufficiently Figure 2. Summary of the scattering scenarios considered here, with their p-value as a function of the mass. We show in green momentum suppressed DM-electron scattering with form factor F (q) ∝ q discussed in Sec. 6.1, in blue exothermic scattering with light mediator F (q) ∝ 1/q 2 and momentum suppressed form factors F (q) ∝ q discussed Sec. 6.2. Here, the splitting between the heavy state and the light state in the dark sector is marginalized to minimize the p-value. In dark red we show accelerated DM by cosmic rays scattering for a fixed mediator mass m X = 100 keV and a fixed ratio between the mediator and the DM mass m X /m DM = 1/15. This is discussed in Sec. 7.
rapidly at energies lower than 2 keV. We achieve this by considering axial-scalar interactions between the accelerated DM and the SM, mediated by a light new mediator with mass around 100 keV. Just as other models of accelerated DM, this scenario is likely to be challenged by other observation probes. We leave a more in depth study of this scenario for future work. We now present our data analysis framework, before discussing each of these model scenarios in detail.

XENON1T
In this section, we review the relevant aspects of the XENON1T experimental apparatus and the electron recoil analysis, with a focus on describing our treatment of the energy reconstruction and statistical analysis that is used throughout this work.

Energy Reconstruction Method
The experiment utilizes a dual-phase xenon Time Projection Chamber [43,[61][62][63][64][65][66][67][68][69][70], to search for weakly interacting particles. When one of the xenon atoms in the Liquid Xenon (LXe) phase recoils or is ionized due to a collision, photons are emitted and detected by photomultiplier tubes (PMTs). This signal is called the prompt scintillation signal (S1). In addition to the photons emitted close to the interaction point, ionized electrons drift inside the detector due to an external electric field. When the electrons reach the Gaseous Xenon layer (GXe) at the top of the detector, they are extracted across the liquid-gas interface, collide with xenon atoms, and produce a proportional scintillation light, known as the S2 signal, which is also measured by the PMTs.
The ratio of S2/S1 provides a handle that enables one to differentiate between Nuclear Recoil (NR) and ER events. Further information about a given event can be inferred by its location inside the PMTs, the time difference between the arrival of the S1 and S2 signals, and the S1 and S2 signal shapes. This complementary information is taken into account in the analysis by the XENON1T collaboration, however, it is not publicly available. When the XENON1T collaboration reports their data, they use the corrected S1 (cS1) and corrected S2 (cS2), which takes into account this additional information.
In their analysis of the Science Run 1 (SR1) data, the XENON1T collaboration provides a scatter plot of (cS1, cS2 b ) (the 'b' subscript signifies that only the PMTs at the bottom of the detector were used for the S2 reconstruction). Rather than using their reconstructed keV-binned energy spectrum, we will use the data from this scatter plot to reconstruct the energies for each event. We do this, since the keV-binned data results in a loss of information, as the XENON1T detector resolution is as low as ∼ 0.3 keV [70] at their analysis threshold ∼ keV. In order to reconstruct the energies, we use the procedures laid out by the XENON1T collaboration in [68] (with additional data taken from [71,72]), which allows us to simulate the detector response and the effects of reconstructing the signal. We use a Monte Carlo (MC) simulation to determine how an ER with a given energy is distributed on the (cS1,cS2 b ) plane, and use a maximum likelihood estimator to find the energy of the event. Below, we refer to this way of reconstructing the energy as "our method", even though it is based on information provided in previous XENON1T papers; we do so to differentiate it from the way the energy was reconstructed by the XENON1T collaboration in their ER analysis paper [1], where they simply use where g 1 = 0.142 and g 2 = 11.4 are the probabilities for one photon to be detected as a photo-electron in the PMT and the charge amplification factor, respectively, and the mean energy to produce a detectable quanta is W = 13.8 eV.
In Fig. 3, we reproduce the (cS1,cS2 b ) scatter plot for events tagged in [1] to have an energy below 9 keV (above this energy, the resolution is ∼keV so the binning leads to only marginal information loss; moreover, the excess is concentrated below this energy, so we will not be concerned with events at higher energies). The color of the points signifies the difference between the energy reconstructed by our method and the simplified formula used in [1], Eq. (3.1). The colored points on the plot are in units of the energy resolution, calculated with our method (see below). Due to the finite size of our MC sample, large numerical errors may occur in rare cases where the calculated likelihood for the reconstructed energy of a given event is small (≤ 5% C.L.) for all energies. To avoid such errors, in those cases we use the simplified reconstruction method, Eq. (3.1). Observed events by the XENON1T collaboration [1] in the (cS1,cS2 b ) plane, for events they tagged as having an energy ≤ 9 keV. The colors of the points correspond to the difference in the reconstructed energy (in units of the energy resolution) between our energy reconstruction calculation and the simplified equation used by XENON1T (Eq. (3.1)). The energy resolution is estimated using our energy reconstruction calculation. The gray dashed lines show constant energy lines using the simplified energy reconstruction given in Eq. (3.1). In blue is the expectation value for the energy interval [keV, 9keV]. The black points are more than 2σ away from this expectation value, and we did not sample the parameter space finely enough with our MC to reliably reconstruct their energy; for these points, we assume simply that their energy is given by Eq. (3.1).
We show in Fig. 4 (left) our calculation of the keV-wide binned energy spectrum and compare it with the XENON1T spectrum. The two spectra are nearly identical. This provides confidence in our energy reconstruction method, and allows us to use the full unbinned energy information for our new physics analyses below. We also include in this plot the background model from [1].
For the formula Eq. (3.1) to be the best estimator for the energy, the variables cS1 and cS2 b should be anti-correlated. While this has been validated for high energies, preliminary measurements appear to suggest that there is only a weak anti-correlation for low energies. In particular, this can be seen from measurements of the 37 Ar line at 2.83 keV presented in [71]. Our MC simulation of 2.8 keV events (assuming a uniform distribution in z) agrees well with the contours in the (cS1, cS2 b ) plane for the 37 Ar data found in [71]: our agreement is better than 3% for the central value, and we can find even better agreement if we change the simulation parameters slightly from those given by the XENON1T collaboration in [68] within their error margins. Our simulation also agrees well with the observed weak correlation between cS1 and cS2 b . This provides further confidence in our energy reconstruction method, especially at the O(keV) energies relevant for the excess events. Fig. 4 (right) shows the energy resolution estimated from our MC (black points), a fit to these MC data (red line), and the energy resolution estimated in [1] (blue line). Our energy resolution is slightly better than that used in [1]. We leave a more detailed investigation of this difference to future work.  Figure 4. Left: Comparison of the naive spectrum reconstructed, and the one by the MC, for energies below 9 keV. In blue is the background model from [1]. While the biggest disagreement, at the lowest-energy bin, seems to be a ≥ 1σ disagreement, we note that this is misleading, as many of the points were reconstructed on the edge of the bin causing small differences to be magnified. An overall good agreement between our MC method and the one used by the XENON1T collaboration is observed, enabling us to use the full unbinned energy information throughout this paper. As an aside, we also note that our binned energy spectrum does not have the same monotonically decreasing spectrum in the ∼5 − 10 keV energy bins. Right: Energy resolution estimated from our MC (black points), a fit to these MC data (σ E [keV] = −0.21 + 0.39 √ E, red line), and the energy resolution estimated by the XENON1T collaboration in [1]

Statistical Method
For our analyses, we use a likelihood ratio test, with unbinned likelihoods. For each signal model, s, that depends on parameters θ s , we find the likelihood of the signal+background hypothesis for the data as a function of the model parameters,

2)
where E i are the reconstructed energies, n is the number of observed events, dN b /dE (dN s /dE) is the background spectrum (signal spectrum), and µ b = dN b /dE (µ s = dN s /dE) are the total expected background (signal) events. We maximize the likelihood to find the best fit points. In order to estimate the significance and quality of our fits, we assume the asymptotic formulas found in [73]; we therefore assume that the negative of twice the log-likelihood-ratio of the signal+background ratio hypothesis compared to the background-only hypothesis is distributed according to a χ 2 distribution, with the number of degrees of freedom set equal to the number of model parameters.
In each of the following sections, we describe how to derive the spectrum of events. The measured spectrum will be modified by detector response effects. In particular, for a given theoretically predicted signal, we modify the spectrum by the effective exposure, E(ω), of the xenon detector. We then smear the resulting spectrum by a gaussian with the resolution presented by the red line in Fig. 4 (right), and then calculate the likelihood.
Two important effects that the XENON1T collaboration treats in their analysis and that we, however, will ignore is an unconstrained contribution to the background from tritium decays and the look-elsewhere effect. Including the tritium background when reporting the significance of a potential DM signal is technically correct, but we will report significances without this contribution. The look elsewhere effect is important for determining the global significance of a particular model to explain the excess; however, we will only report the local significance, which is sufficient for the purposes of this paper.

Absorption
We consider first models of bosonic DM, confronting them with the XENON1T measurement. Three cases are considered: pseudo-scalar (axion), scalar, and vector bosons. For each we explore the non-relativistic case, in which the boson constitutes the DM, and the relativistic case, for which the boson is produced in the Sun.
In the case of bosonic DM, the rate of events in the XENON1T detector per unit energy is dR DM where the energy, ω, is kinematically constrained to equal the DM mass (ignoring small non-relativistic corrections of order the DM energy), and is then smeared to account for the detector resolution as described in §3. The total number of events in the relevant XENON1T energy window is then obtained by convolving the above rate with the effective exposure E(ω) reproduced in §3 and integrating over energy. The DM flux, Φ DM , is the same for all bosons, and depends only on the DM relic density, ρ χ , and the mass of the light boson (4. 2) The absorption cross section, σ I abs , depends, however, on the interaction of a given light boson I with the bounded electrons in the liquid xenon. Here we consider three cases: For light bosons produced in the Sun, the differential event rate per unit energy can be written as dR Sun where the differential solar flux dΦ Sun I /dω depends on the production mechanisms of the light bosons I inside the Sun's environment and needs to be treated case by case. Below we discuss solar axions in Sec. 4 Figure 5. Solar flux spectra for the axion production from ABC processes [49] (blue) and Primakoff production [57] (yellow), for scalars from bremsstrahlung [11,14] (green), and for dark photons (red). On the left, we assume a massless boson, while on the right, the kinematical threshold due to a finite boson mass, m I = 2.6 keV is shown with solid lines. To highlight the spectral features, the plots are normalized such that the total integrated flux in the energy window of interest for XENON1T , ω ∈ [1, 10] keV equals 1 event/cm 2 /s.

Axion-Like Particles
We consider an axion-like particle (ALP) of arbitrary mass m a that couples to photons and electrons, The ALP can be absorbed inside the detector material leading to an ioniziation signal. The cross section for this so called axio-electric (AE) effect [4,5,74,75] can be written as [8] σ where ω a = m 2 a + k 2 a is the energy of the ALP and v a is its velocity. We take the photoelectric cross section, σ PE (ω a ), from [76], which agrees reasonably well with experimental data above 30 eV. The above formula is approximate, and chosen to correctly reproduce the results obtained in the non-relativistic limit, v a 1, and in the relativistic limit, v a → 1. In what follows, we will derive the XENON1T best-fit regions for g aee as a function of the ALP mass. Theoretically, however, g aee is often related to g aγγ and for the ALP DM case, X-rays measurements can then be used to exclude part of the parameter space. It is therefore interesting to understand the theoretical relation between the two couplings, which will allow us to identify viable ALP models. As we shall see below, three conclusions can be drawn: 1. Fitting the data with QCD axion DM requires a high degree of fine tuning of its ultraviolet (UV) couplings to electrons and the UV anomaly with respect to electromagnetism.
2. More general ALP DM requires suppressed couplings to photons in the UV, which typically implies a non-anomalous global symmetry with respect to QED.
3. Standard solar ALPs are excluded by stellar constraints, motivating chameleon-like ALPs to be discussed in Sec. 5.
To understand these statements, let us briefly discuss the origin for g aee and g aγγ . The parametrization of Eq. (4.4) can be mapped to concrete models where the pseudo-Nambu-Goldstone boson (pNGb) of a spontaneously broken global symmetry couples to the photons and electrons. An arbitrarily small mass m a can be introduced as a soft breaking of the pNGb shift symmetry. More explicitly, we can write where f a is the ALP decay constant and E eff parametrizes the effective coupling to photons, which is related to the UV parameters through Here C UV is the UV coupling of the axion to electrons, while E UV is the UV anomaly with respect to electromagnetism, which is model dependent. A(x) parametrizes the electron loop function, for m a m e . This feature can be traced back to the fact that in the presence of a purely derivative coupling to electrons, only the effective operator ∂ 2 aFF is generated below the electron threshold [77]. If E UV is non-zero, the electron coupling is modified by the running contribution induced by the photon coupling. At low energies, one finds For the QCD axion, the coupling to the gluon field strength gives further contributions to the effective photon and electron couplings generated by the mixing of the axion with the QCD mesons below the confinement scale [78], The strong X-ray limits on E eff together with the O(1) contribution from QCD explains why QCD axion DM must be tuned to address the anomaly. Various axion models have been studied, where the different hierarchies between the electron and photon couplings are realized: Figure 6. Left: Allowed parameter space for ALP dark matter in the (m a , g aee ) plane. The ALP decay constant is plotted on the right y-axis. The red star is our best fit point in Eq. (4.10) and and the dark red regions are the 1σ and 2σ regions. In gray, we show the bounds from star cooling of red giants [54,55] and white dwarfs [51], and the current direct detection constraints from Xenon1T and PandaX [43,60]. The gray dotted contours show the X-rays constraints for different values of E UV , the shaded region on the top right is excluded by X-rays even for E U V = 0. Dashed brown contours show the initial misalignment necessary to get the right DM relic abundance. Right: Signal shape for the best fit point in Eq. (4.10). The black dots are the XENON1T data, the gray shaded region is the expected background, the blue line is the signal shape after smearing, and the blue shaded region is the resulting signal plus background distribution.
• KSVZ models, where C UV = 0, and the electron coupling is only generated from the photon coupling via the running [81,82].
• Photophobic models where E UV = 0 and the electron coupling dominates the phenomenology. See [83] for a general discussion of photophobic ALPs and the Majoron [84,85] as a particularly motivated example of this coupling structure.
Of the above, and in the absence of tuning, only the Photophobic ALPs can fit the XENON1T hint without being excluded, if they are DM.

ALP Dark Matter
If the ALP is DM, the axio-electric effect should be treated in the non-relativistic limit with E m a , and thus the energy absorbed by the bounded electron in the detector is equal to the axion mass. Consequently, in order to explain the XENON1T signal the ALP masses must be around m a ∼ 1 keV. The predicted spectrum is a narrow peak around the ALP mass, with the observed signal spreading into several bins from detector resolution effects that smear the predicted signal, as shown in Fig. 6 (right). The 1σ and 2σ bands of our likelihood fit is shown in red in Fig. 6 (left), where the best fit point is m a = 2.5 keV , g aee = 2.5 × 10 −14 , 2log(L(S + B)/L(B)) = 15.7 . (4.10) The number of signal events is given by, where we used that σ PE = 1133 cm 2 /gram and the effective XENON1T exposure, E(E).
The predicted coupling to electrons fixes the decay constant to be f a /C eff 10 10 GeV shown on the right y-axis. We further show constraints from white dwarfs [51] (dark blue) and red giant [54,55] (light blue) cooling as well as terrestrial limits from PandaX [43,60] (light green) and the XENON1T S2-only analysis [43] (darker green).
If the coupling to photons is non-vanishing, as predicted in any untuned UV complete scenario, the ALP DM with the desired range of masses and decay constants is severely challenged by its large decay rate into di-photons, Imposing that the ALP is stable on timescales of our Universe gives which gives a strong upper bound on the coupling to photons in order for our best fit point to be stable. Even stronger constraints on the diphoton width come from observations of the cosmic X-ray background (CXB) [86]. The best fit ALP is predicted to produce monochromatic photon-lines at frequency ν a = 3.1 × 10 17 Hz m a 2.6 keV . (4.14) A very conservative bound can be extracted by requiring the intensity of the photon line to be less than the measured CXB background at that frequency, which is ν a I νa (2.3 ± 0.2) × 10 −11 W m −2 rad −1 . Using this procedure, we find This bound is very similar to the one obtained in [87] and could be improved by looking at individual sources and performing background subtraction. Future X-ray missions like Athena [88], as well as new techniques as line intensity mapping [89,90], can then be an important to further constrain this model. On the left of Fig. 6, we illustrate this limit with dotted gray lines and gray region, for different values of E UV . Interestingly, this bound excludes a portion of the parameter space even if E U V = 0, due to the irreducible one-loop contribution to the photon coupling in Eq. (4.7). We see that a very small E UV value is needed to explain the XENON1T anomaly, disfavoring most existing ALP models, and in particular the QCD axion, and hinting towards photophobic ALPs such as the Majoron. We note that for the Majoron, the XENON1T signal is correlated with possible future signals in µ + → e + + a that could be seen at future high intensity muon facilities like Mu3e Figure 7. Predicted spectrum for the solar production of a photophobic axion with a best-fit value of m a = 1.3 keV (left) and a photophilic axion with a massless axion as the best fit model (right). The dashed and solid lines show the signal spectrum before and after detector smearing effects, respectively. The measured XENON1T data is shown as black dots, while the gray-shaded region is the expected binned background, and the blue-shaded region is the predicted binned signal.
(see [58] for further details). Depending on the final seesaw scale one could further explore the parameter space looking at µ → eγ at MEGII [85].
It is also interesting to ask what are the conditions for an ALP DM addressing the anomaly to have the observed DM relic abundance. If one considers a generic axion-like particle with a non-dynamical mass m a , the correct relic abundance can be generated in the region of interest via the misalignment mechanism [87,[91][92][93]] On the left of Fig. 6, we show in dotted brown lines two O(1) values for the misalignment angle, θ 0 , for which the observed DM relic abundance is obtained. We conclude that no tuning is needed to address the ALP DM relic density in the region of interest.

Solar ALPs
ALPs can also be produced in the Sun through processes involving the electron and photon couplings of Eq. (4.4). Here we study solar production, not making any assumptions on the ALPs relic density. We consider both the relativistic case, m a T , for which the energy absorbed by the bounded electrons is independent of the ALP's mass, as well as the nonrelativistic case, m a T , in which the spectrum is significantly modified, improving the fit to the XENON1T data. See Fig. 5 for the different spectra with a massless (relativistic) and massive (non-relativistic) ALPs.
Two production mechanisms are of interest: (i) the "ABC" processes: atomic recombination and de-excitation, bremsstrahlung, and Compton scattering, all depending on the value of g aee [49]. (ii) The Primakoff process [57], which is the conversion of photons into axions in the electromagnetic fields of the electrons and ions making up the solar plasma. This is the dominant production mechanism in the energy range relevant for XENON1T, which depends on g aγγ .
We discuss photophobic ALPs where both production and absorption are controlled by g aee , so that the total signal rate scales as Here Φ ABC, Xe is the integrated ABC flux in the energy window that is relevant for the XENON1T experiment, calculated for m a = 1.3 keV. We also consider photophilic ALP models, where the ALP coupling to photons contributes substantially to the production, while the ALP coupling to the electrons controls the absorption rate. Here one finds (4.18) where Φ P, Xe is the integrated Primakoff flux, once again, in the energy window that is relevant for the XENON1T experiment, and with a massless ALP. We find that for the energy range of interest (ω keV), the ABC productions are subdominant for g aee /g aγγ 16 MeV, or equivalently, E eff 27 C eff . This is satisfied in many standard QCD axion models (see also [94] for an explicit model where E eff takes on very large values).
The best fit points in these scenarios are where the spectrum for the two cases is shown on the top of Fig. 7. The peaked structure of these signals is due to the convolution of the solar fluxes with the detector smearing and efficiency, suggesting that in principle, one may be able to differentiate between the two solar production mechanisms with more data. We note that the solar production of a massless photophobic axion does not reproduce well the spectral shape of the data. The reason can be traced back to Fig. 5 where one can clearly see that the ABC production does not shut off fast enough below 2.5 keV, leaving an excess signal in the lowest energy bin. On the other hand, the massless photophilic model provides the best-fit model. The parameter space for the solar production of the photophilic and photophobic axions is shown in Fig. 8. On the left plot, we show in red the photophilic 1σ and 2σ best-fit regions in the g aγγ g aee versus m a plane. As mentioned above, the best-fit value lies outside the plot at m a = 0. Stellar cooling constraints [10,51,54,55] are shown in blue. For the same model with the best-fit value m a = 0, the middle plot shows constraints in the g aγγ − g aee plane. The 1σ and 2σ best fit regions are shown in red, white dwarf (WD), horizontal branch (HB), and sun cooling limits are marked with dashed-blue lines, and limits form Figure 8. Left: 1σ and 2σ best-fit regions (red) for the photophilic solar axion scenario. The best fit point corresponds to the m a = 0 case and lies outside the plot. White dwarfs (WD) [51], red giants (RG) [54,55], horizontal branch (HB) [95], and Sun [10] stellar cooling constraints are shown by the blue-shaded regions, and the CAST limit [96] is shown in orange. Middle: 1σ and 2σ best-fit regions as in the left plot (red) for the best-fit massless photophilic solar axion model, but here shown in the g aγγ − g aee plane. Stellar cooling constraints are indicated with dashed blue lines, while limits from CAST [96] are shaded in orange; arrows point to regions that are allowed. The theoretical axion model lines are shown in the bottom-right part of the plot. Right: Same as the left plot, but for the photophobic solar axion model. The red star indicates the best fit point in this case.
CAST are shown in orange. We also show the predicted model lines for the DFSZ and KSVZ axion models. Finally, on the right plot we show the 1σ and 2σ best-fit regions for the photophobic case in red and the stellar constraints in blue. We conclude that for all cases, the solar axion explanation to the XENON1T anomaly is in severe tension with stellar cooling constraints. In Sec. 5, we discuss briefly a possible mechanism to circumvent these bounds.

The Scalar
Consider now a scalar, φ, that couples to photons and electrons L scalar = g φγγ 4 φF µν F µν + g φee φēe . (4.21) The cross section for scalar-electric (SE) effect can be written in terms of the photoelectric one as [14,97] where ω φ = m 2 φ + k 2 φ is the energy of the scalar φ, v φ its velocity, and σ PE (ω φ ) is again the photoelectric cross section already used in Eq. (4.5). Notice that in the case of scalar DM, the expression above leads to a suppression of the absorption rate of v 2 DM 10 −6 . The parametrization of Eq. (4.21) can be mapped to concrete models. Two particularly motivated scenarios are (i) a light SM singlet mixing with the SM Higgs doublet, and (ii) the dilaton from a spontaneously broken conformal-invariance. Below we briefly review these models, pointing to the distinct nature of their photon and electron couplings. Figure 9. Left: Allowed parameter space for scalar dark matter in the (m a , g aee ) plane. The red star is our best fit point in Eq. (4.10) and and the dark red regions are the 1σ and 2 sigma regions around it. In blue we show the bounds from star cooling of red giants [54,55] and white dwarfs [51] and in green the present direct detection constraints from Xenon1T and PandaX [43,60]. Right: Signal shape for the best fit point in Eq. (4.24). black dots are the Xenon1T data. The gray shaded region is the expected background, the blue line is the signal shape and the blue shaded region is the resulting signal plus background distribution.
A singlet obtaining a VEV would generically mix with the Higgs through the quartic λ φH φ 2 H † H. The mixing can be written in terms of the ratio of the Higgs and the singlet VEVs, sin θ = v/f , and the final couplings of the singlet to photons and electrons are generated once the mixing is resolved [98][99][100] (4.23) Here κ SM φγγ 11/3 + O(m φ /m e ) 2 is the asymptotic value of the SM loop functions from W 's and Standard Model fermions for m φ m e , and we fixed the coupling of the Higgs to electrons to be y e = m e /v, ignoring possible deviations from its predicted SM value. In this simple framework, the ratio between the photon and the electron coupling is fixed to g φγγ m e /g φee = 4.2 · 10 −3 , and a large coupling to nucleons is also generated from the couplings of the Higgs to gluons.
Conversely, if the scalar in Eq. (4.21) is a dilaton, its coupling to the SM are more model dependent and controlled by the infrared (IR) trace anomaly contributions induced by direct UV couplings between the CFT and the SM. In this framework, the dilaton mixing with the Higgs can be arbitrarily suppressed [101], and the prediction of Eq. (4.23) are changed. In particular, for a dilaton, one can entertain the possibility of a loop-suppressed photon coupling, which decouples as m φ /m e . Thus, analogously to the ALP case, we consider two possibilities: • The Higgs-mixing scenario, where the ratio of the relative strength of photon and electron couplings is fixed.  Figure 10. Predicted spectrum for the solar production of a photophobic scalar with a best-fit value of m φ = 2.3 keV (left) and photophilic scalar with a best fit value m φ = 1.7 keV (right). The dashed and solid lines show the signal spectrum before and after detector smearing effects respective. The measured XENON1T data is shown as black dots while the gray-shaded region is the expected binned background and blue-shaded region is the predicted binned signal.
• The photophobic dilaton scenario, where g φγγ is suppressed as m φ /m e and the electron coupling dominates the phenomenology.

Scalar Dark Matter
As in the ALP case, the absorption spectrum of the scalar is sharply peaked around its mass, as can be seen in the spectrum plotted for the best-fit scalar DM model on the right of Fig. 9, with values, The number of signal events is given by The predicted coupling to electrons corresponds to a mixing angle with the Higgs of order sin θ 2.1×10 −7 . On the left of Fig. 9, we show in red the 1σ and 2σ best-fit regions for the scalar DM case in the g φee -m φ plane. On the right y-axis, we map g φee to the mixing angle for the doublet-singlet model, Eq. (4.23). Regions excluded by RG cooling constraints [54,55] are shown in light blue, while the exclusion regions due to the XENON1T S2-only analysis [43] and PandaX-II analysis [60] are shown in dark and light green, respectively. As one can see from Fig. 9 the scalar DM cannot explain the XENON1T excess because of the large suppression of its absorption rate compared to the ALP case.

Solar scalar
Much like ALPs, light scalars can be produced in the Sun, whether or not they constitute DM. For a photophobic scalar, the production in the Sun is dominated by electron-nucleus scalar-bremsstrahlung N + e → N + e + φ. The rate can be obtained through the rescaling of the regular photon-bremsstrahlung by the ratio of the matrix elements squared. Doing so we find The above agrees numerically with the one given in [49]. Similarly, a photophilic scalar is produced via the Primakoff process, with a rate similar to that of the ALP. The predicted fluxes are shown in Fig. 5. We fit both the photophilic and photophobic scalar to the XENON1T data. We find the best-fit points for which we show with dashed and solid blue lines the predicted spectrum before and after smearing respectively in Fig. 10. As before, the gray region shows the expected binned background while the blue fillings show the binned contribution of the signals. The XENON1T data are shown in black. In Fig. 11, we show in red the 1σ and 2σ best-fit regions of the solar production of a photophilic (left) and photophobic (right) scalars. In both cases, only a massive scalar can explain the XENON1T anomaly. This is in contrast to the photophilic ALP case for which the massless ALP provided the best fit. The reason for this can be traced back to the rapidly falling absorption rate at high energies, Eq. (4.22). This implies a soft spectrum, which must be cut off at production through kinematic effects from a massive particle. For the photophilic case, combined stellar cooling constraints are shown in blue while those are shown separately for the photophobic scalar.

The Dark Photon
As the final absorption scenario of this section, let us consider the dark photon A , a massive gauge boson of a broken (dark) gauge group U (1) . The dark photon may couple to ordinary matter via its kinetic mixing with the visible photon [102]. Much as in the previous sections, we consider the absorption of a dark photon DM and the production of a dark photon in the Sun as explanations to the XENON1T anomaly.
The relevant interactions are where m A is the mass of the dark photon, F µν and F µν are the photon and dark photon field strength respectively, and is the kinetic-mixing parameter. After the kinetic terms are Figure 11. Left: 1σ and 2σ best-fit regions (red) for the photophilic solar scalar (left) and photophobic solar scalar (right) scenarios. White dwarfs (WD) [51], red giants (RG) [54,55], and horizontal branch (HB) [97] stellar cooling constraints are shown by the blue-shaded regions. The red stars indicate the best fit points in both cases. In contrast to a photophilic ALP, the scalar must be massive in order to explain the data, due to its sharply rising absorption rate at low energy.
diagonalized, the dark photon couples to the electron vector current with a coupling strength e, and the dark-photon absorption cross-section can be related to the SM photelectric cross section by a simple rescaling, σ DP (E) = · σ PE (E) . Inside a medium, the propagation of electromagnetic fields is determined by the polarization tensor Π µν = e 2 J µ EM , J ν EM , which can be decomposed into longitudinal and transverse components as, where L,T are the polarization vectors. In general, in-medium effects should be accounted for in order to correctly compute the dark photon absorption rate. We implement these effects following the discussion in [11,103]. For dark photon DM with mass near 1 keV, we find that the absorption is dominated by the transverse modes, and the inclusion of the longitudinal ones modifies the rate by less than 10%.
For m A larger than the typical solar plasma frequency ω pl 0.3 keV, the production of dark photons in the Sun is dominated by the transverse modes at energies ω ∼ keV. In such a case the flux at the Earth is found to be [11] where the interaction rate Γ T is dominated by free-free absorption and Compton scattering. At lower masses, the behavior of the flux from the Sun depends crucially on the nature of  Figure 12. Left: The 1σ and 2σ best-fit regions (red) for dark photon DM with a Stuckelberg mass. Light and darker blue represent the RG and HB cooling limits, respectively, and the light green region is excluded due to the XENON1T S2-only analysis [43]. Right: An example of the predicted spectrum for dark photon DM using the best-fit value m A = 2.5 keV. The dashed and solid lines show the signal spectrum before and after detector smearing effects, respectively. The measured XENON1T data is shown as black dots, while the gray-shaded region is the expected binned background. The blue-shaded region is the predicted binned signal.
the dark photon mass [103]. For a non-dynamical Stuckelberg mass, the dark and visible sectors decouple in the m A → 0 limit for an on-shell A . As a consequence, the rate of production/absorption of the transverse modes falls off as (m A /T ) 4 , where T is the Sun's temperature. Adding a Stuckelberg mass to the dark photon will then cut off the solar flux around ω m A as shown in Fig. 5. Conversely, if the dark photon's mass is generated through the VEV of a dark Higgs, then the ratio between the dark photon mass and the dark Higgs mass is controlled by the ratio of the Higgs quartic and the dark gauge coupling m h /m A ∼ √ λ/e . For m h ∼ m A the production/absorption of a dynamical dark photon therefore goes predominantly through the radial component in the m A → 0 limit [4]. This case shares many features with the absorption scenarios discussed so far, and we will not discuss it here for the sake of brevity.

Dark Photon Dark Matter
If the dark photon plays the role of DM, its predicted absorption spectrum in the XENON1T detector is very similar to the other bosonic DM cases discussed in the previous subsections. On the right panel of Fig. 12, we show an example for the best-fit model, Dashed and solid lines represent the unsmeared and smeared spectrum, respectively. We see that, as with the axion and scalar, the spectral shape is peaked around the dark photon mass, and detector resolution allow for a reasonable fit to data. As in previous plots, the expected binned background is shown in the figure in gray, while the binned signal is shown in blue. The XENON1T data is presented with black dots.
On the left plot of Fig. 12, we show in red the 1σ and 2σ best-fit regions for dark photon DM with a Stuckelberg mass. In light and darker blue, we show the RG and HB cooling limits respectively and in light green, the constraint from the XENON1T S2-only analysis [43]. We learn that the explanation of the XENON1T anomaly with dark photon DM is viable.
Finally, two remarks are in order. First, a major advantage of dark photon DM compared to the ALP and scalar cases is that the decay rate of a keV dark photon into SM particles is extremely suppressed [4]. The only decay channel allowed kinematically is A → 3γ, which is induced by dimension eight operators generated at one loop from the electron coupling. The width of this process is suppressed by ∼ α 5 2 (m A /m e ) 8 , and the dark photon explanation to the XENON1T anomaly is safely outside any bound from decaying DM. Second, the misalignment mechanism, which comfortably explains the scalarand axion-DM relic densities, fails to generate the observed dark photon abundance unless non-minimal coupling to gravity are taken to be unnaturally large [87]. Also the contribution from inflationary fluctuations cannot explain the DM relic abundance in this mass range [104]. Non-minimal mechanisms of dark photon productions have been suggested recently [105][106][107][108]. In particular, the mechanism in [105] can accommodate the correct DM abundance for a keV dark photon by postulating a coupling φF F between the inflaton, φ and the dark photon.

Solar Dark Photon
For a Stuckelberg dark photon produced in the Sun, the best fit point is As for the case of the scalar, the presence of a mass cuts off the low-energy flux to reduce the signal yield in the lower XENON1T bins. The unsmeared and smeared spectrum, together with the binned background, signal, and data is shown in Fig. 13 (right). In the left plot, we show the best fit region for the model, together with the HB and RG stellar cooling bounds. We learn that as for the scalar and ALP, the best-fit regime is robustly excluded by the astrophysical bounds.

Chameleon-like ALPs: Circumventing the Stellar Cooling Bounds
As discussed in the previous section, particles produced in the Sun are excluded as an explanation for the XENON1T anomaly due to stringent stellar cooling constraints. These constraints arise from the energy loss induced by the emission of light bosons in the star environment. In principle, the constraints can be evaded if the properties of these particles depend on the environment, thereby allowing for a suppressed production in stars. Such   Table 1. Summary of the bounds on the electron coupling g aee from star cooling with the rough value of the density at the core.
Chameleon-like particles have been studied extensively in a broader context, for example in order to evade fifth-force constraints or play the role of dark energy (see e.g. [22,109]), but also for the particular case of ALPs [23][24][25][26][27][28][29].
Here we focus on the specific case of chameleon-like ALPs (cALPs). While most previous work has focused on suppressing the axion-photon couplings in stars, we choose to study the suppression of the axion-electron coupling, which is sufficient to open up the parameter space for solar ALP models that predict either only the latter or both couplings (see Fig. 8). Below we entertain a simple novel model of this kind, leaving a more general framework as well as possible generalizations for future work.
For the axion-electron coupling, g aee , four stellar cooling bounds may need to be addressed: RG, WD, HB stars, and Sun cooling. The resulting bounds on the ALP electron couplings are summarized in Table 1. Among the four, the solar cooling bound is the least constraining and does not exclude the ALP explanation of the XENON1T anomaly (see for instance Fig. 8). The HB bound is in marginal tension with the XENON1T explanation if one accounts for the potentially large systematical uncertainties.
For this reason, we focus here mostly on evading the RG and WD bounds. The energy losses in RG and WD are dominated by the production of light bosons in the highly degenerate core, where the central density is of order ρ WD,RG ∼ MeV 4 , roughly four orders of magnitude larger than the core density of the Sun (see Table 1). Therefore, a model that suppresses production only in high density stars while keeping it unaltered in low density ones may evade RG and WD constraints and, at the same time, leave the ALP production in the Sun unchanged. To illustrate this point, we now discuss a simple model for which production in high-density objects is suppressed. A more thorough study of the constraints, as well as a UV-completion of this model, is left for future work.
We note that an alternative possibility, which we do not pursue here, and which was discussed in [25], is to suppress solar production in order to relax stellar cooling bounds with respect to direct detection. Indeed, for a given suppression factor, S 1, in the solar production of ALPs, the solar flux scales as Sg 2 aee , while the the solar detection rate scales as Sg 4 aee . Increasing g aee , while keeping the detection rate fixed, implies a relative suppression in the solar cooling bound, which scales as S 1/2 . This observation could play an important role in generalizing cALPs to the case of light scalars and dark photons.
Consider a complex Standard Model (SM) singlet, S, charged under a Peccei-Quinn (PQ) symmetry [110] and a real SM-singlet X. The two fields are odd under the same Z 2 , and X couples to density. Below a given cutoff scale, M , we assume that the following Z 2 -invariant interactions are generated The interaction term with the electrons can be induced in a Froggatt-Nielsen construction [111], where the SM electrons carry charges under the same U (1) PQ that rotates the complex singlet S. Ensuring that under that symmetry [e L ] + [e R ] + 1 = 0, allows the operator above while forbidding unwanted others (we normalize the singlet charge to be [S]=1). The cut-off scale in such a construction would correspond to the scale of the vector like-fermions required to generate this interaction [112,113]. For simplicity, we consider the theory below the Higgs mass scale, ignoring further complications that might arise above it. The potential V (S) is such that S develops a VEV, S = 1 √ 2 (f a + s)e ia/fa , where s is the massive singlet with mass m s = √ λ s f a and a is the ALP, which is massless up to the addition of operators breaking the U (1) P Q explicitly. For m S m X , we can neglect the s dynamics and write the effective coupling of the ALP to the electrons where ρ is the matter density and Θ(x) = 0 if x < 0 and 1 otherwise. The second term in Eq. (5.1) expresses nothing more than the idea discussed in [114]: at low densities, X has a allowed by star cooling bounds (we fix λ S = 1) once the ALP coupling to electrons g aee is fixed at its best fit point. Right: Parameter space of the chameleon ALP produced in the Sun. The star cooling bounds from HB, WD and HB stars summarized in Table 1  negative mass, obtaining a VEV. Conversely, at high densities, its squared mass is positive, and the Z 2 symmetry is restored. As shown in Eq. (5.2), for ρ M 2 m 2 X one finds X = 0, and the coupling of the ALP, a, to electrons vanishes, shutting down its production in stars.
Several conditions limit the parameter space of the example above: • First, in accordance with the discussion above if we want to avoid WD, RG, or HB constraints while keeping the Sun flux unsuppressed. The allowed parameter space in the (m X , M ) plane is shown in the white band of Fig. 14. • Second, the quartic λ SX X 2 |S| 2 was omitted from Eq. (5.1) even though it is allowed by all symmetries. When S obtains a VEV, such a quartic induces a new mass term for X that could destroy the density-dependent VEV of X. To avoid this, we require λ SX f 2 a m 2 X . Independently of its bare value, this quartic will be generated at one • Finally, we need to avoid the phenomenological constraints on X. In the limit m s m X the coupling of the chameleon field X to electrons g Xee = S m e c ee /M 2 is enhanced compared to the one of the ALP and is bounded from below by A conservative bound on the parameter space can be obtained by requiring g min Xee to satisfy the stellar cooling constraints [97]. Setting g aee to the XENON1T best fit and setting λ s = 1, we get the maximal value of λ max X allowed by stellar cooling constraints. In the mass range 10 −4 keV m X 10 keV the RG bounds are the most stringent, and we find, λ X λ max X ≡ 7 × 10 −8 2.6 × 10 −12 g aee 2 g Xee 6.7 × 10 −16 The above reveals a hierarchy between the quartic of the PQ-breaking field λ S , and that of the chameleon, λ X , needed in order to make this model phenomenologically viable. This hierarchy might be difficult to realize quantum mechanically. For instance, two loop contributions to the singlet and chameleon quartics induced by their electron couplings, will act to make them both of the same order. Higher chameleon masses weaken the phenomenological bounds, allowing for a milder hierarchy between the couplings, but at the price of lowering the cut-off scale M as in dictated by Eq. (5.3).
In summary, cALPs could avoid stellar cooling bounds. As shown in Fig. 14 right, the stellar cooling from dense stars can be circumvented if a new light scalar X controls the coupling of the ALP to matter. If chameleon-like scalar X lies in the mass vs cut-off range shown in Fig. 14 left, its potential is modified by density dependent effects. In the simplest construction, the chameleon-like scalar can be light and the cut-off of the theory can be arranged to be sufficiently high if a hierarchy between the quartic of the PQ radial mode and the quartic of the chameleon is arranged as shown in Fig 14 left.

Dark Matter-Electron Scattering
If DM interacts with electrons, it can scatter off the electrons in the target material and produce an electron recoil signal [15]. Due to the distinctive kinematics of this process, the electron recoil signal for "standard" DM-electron scattering peaks at recoil energies well below the keV energies needed to explain the XENON1T data; this standard process is thus in conflict with lower threshold direct-detection searches. However, we will investigate here several novel signals to check whether they could explain the XENON1T data: exothermic scattering off electrons as well as DM-electron interactions that increase as a function of the momentum transfer (up to some cutoff scale). We will find that exothermic scattering off electrons work well, and momentum-dependent interactions also provide a potential explanation of the XENON1T excess.

Standard DM-Electron Scattering
We begin by reviewing the standard DM-electron scattering kinematics and formalism discussed in [15,17], before discussing momentum-dependent and exothermic interactions. Consider a DM particle with mass m χ and initial velocity v which scatters off a bound electron transferring a momentum q. Energy conservation of the DM-atom system gives, where ∆E e is the energy transferred to the electron and m N is the mass of the nucleus. This can be written as As the electron initially is in a bound state, it can have arbitrary momentum and hence the momentum transfer q could take any value. The maximum energy that can be deposited is then found by maximizing the above equation with respect to q, and we get For m χ m N , µ χN m χ , and almost the entire kinetic energy of the incoming DM particle can be transferred to the electron. Since the typical DM halo velocity is v ∼ 10 −3 , a DM particle with mass of a few GeV can in principle then produce a O(keV) electron recoil. However, for DM with masses above the MeV scale, the typical momentum-transfer scale is set by the electron's momentum, given by q typ ∼ Z eff αm e ∼ Z eff × 4 keV, where Z eff is the effective charge seen by the electron. From Eq. (6.2) (neglecting the second term, which is usually small), ∆E e ∼ 10 −3 q typ ∼ Z eff × few eV. While higher momentum transfers are possible, they are dramatically suppressed, since it is unlikely for the electron to have a momentum that is much higher than the typical momentum.
We can see this behavior in more detail by calculating the atomic form factor f 1→2 (q), which captures the transition from state 1 to state 2, where ψ 1 (x) (ψ 2 (x)) is the initial bound-state (final-state) electron wavefunction. There are various methods to calculate the wave functions. For simplicity, we here follow [15], taking the initial bound-state wave functions from [115] and the outgoing wave functions to be plane waves. We here also do not subtract the identity operator from the operator e iq·x ; our form factors will therefore not be correctly behaved for q O(keV), since we are using plane waves as the outgoing wave functions, which are not orthogonal with the bound state wave functions. However, this does not affect our results below, since DM-electron scattering does not typically sample the atomic form factor in the momentum region q O(keV). In order to calculate the scattering rates, we will multiply these form factors by a Fermi function (see below). In future work, it will be interesting to calculate these form factors using more accurate wave functions, especially for the outgoing electron. We will comment on including relativistic corrections for high q and high ∆E e below. We plot |f 1→2 (q)| 2 in Fig. 15 (without the Fermi function) for different initial electron shells {n,l} of the xenon atom and for two different outgoing energies of the electron E ≡ ∆E e − E nl , where E nl are the binding energies of the respective shells. We see that the form factor drops sharply for q αm e for the 5p shell, but also that the form factor peaks at larger q for deeper shells (which have a higher velocity). For every shell, the peak also shifts to higher q for higher ∆E e . In order for an electron in any shell to give ∆E e 1 keV, we need q MeV (see Eq. (6.2)). This is possible, but highly suppressed.
We can now write the cross section for the scattering rate as [15,17] where |F DM (q)| 2 is the DM-electron interaction form factor and σ e is the reference DM-electron cross section defined as where |M free (αm e )| 2 is the absolute value squared of the matrix element describing the elastic scattering between DM and a free electron. We also include a Fermi factor,  {5p, 5s, 4d, 4p, 4s, 3d, 3p, 3s, 2p, 2s, 1s} [116,117]. The differential scattering rate will then be given by, where we sum over all the occupied initial shells {n, l} with respective binding energies E nl . The η(v min ) is defined by, where v min is given by, v min = ∆E e q + q 2m χ (6.11) and (normalized as d 3 v g χ (v) = 1) where v χ is the DM velocity in the Earth frame, and v E is the Earth's velocity in the galactic rest frame. We take a peak velocity of v 0 = 220 km/s, an average Earth velocity of v E = 240 km/s, and a galactic escape velocity of v esc = 544 km/s. We set ρ χ = 0.4 GeV/cm 3 . The DM form factor depends on the precise DM-electron interaction, but we will consider F DM = 1 "heavy" mediator (6.13) light" mediator (6.14) F DM = q αm e q-dependent "heavy" mediator , (6.15) where "heavy" and "light" refer to the mass of the mediator, which is respectively above or below the typical momentum transfer. We have also included a momentum-dependent Figure 16. Electron recoil spectra for standard DM-electron scattering for m χ = 10 GeV and σ e = 10 −40 cm 2 for three different form factors. Solid lines show results calculated with plane waves as outgoing electron wave functions while the dotted lines include relativistic effects using results from [118]. See text for details.
form factor from an interaction that grows at higher q; this is predicted, for example, with a pseudo-scalar-scalar operator, iχγ 5 χēe mediated by a heavy scalar field. The resulting differential rates for m χ = 10 GeV are shown in Fig. 16 for σ e =10 −40 cm 2 .
For standard DM-electron scattering, relativistic corrections are important for high ∆E e and therefore for high q. We calculate the rates using the available atomic form factors with relativistic corrections found in [118]; these form factors are given for q ≥ 100 keV, and we show the spectra in Fig. 16 with dotted lines. We see that they predict a larger signal rate in the region relevant for explaining the XENON1T excess than those predicted with nonrelativistic form factors (solid lines). Nevertheless, we see that the spectra for F DM ∝ 1 and especially F DM ∝ 1/q 2 are unable to explain the XENON1T signal without being in dramatic conflict with lower-threshold direct-detection searches from, e.g., XENON1T (S2-only analysis) [43] (for heavy mediators) and SENSEI [44] (light mediators), due to the steeply rising form-factor for lower q. The q-dependent heavy mediator, with F DM = q αme , however, does provide an adequate fit to the XENON1T excess. The best-fit point is given by We show in Fig. 17 (left) the 1σ and 2σ regions in the σ e versus m χ parameter space The red line is our best fit values in Eq. (6.16) and the dark red regions are the 1σ and 2σ regions. We show the spectrum of the best-fit point in Fig. 17 (right). We note that the q-dependent form-factor requires a rather low cutoff to produce the needed cross sections (see right y-axis of Fig. 17 (left)), which is likely in tension with collider bounds [59]. In the above left figure we show on the right y-axis the values for the cutoff divided by the mediator's coupling to DM and electrons. Figure 17. Left: Allowed parameter space for DM-electron scattering through a q-dependent heavy mediator, with F DM = q αme . The red line is our best fit values in Eq. (6.16) and the dark red regions are the 1σ and 2σ regions. On the right y-axis we show the required mediator mass divided by its couplings to DM and electrons, needed in order to obtain the corresponding cross section. Right: Signal shape for the best fit point in Eq. (6.16). The black dots are the XENON1T data, the gray shaded region is the expected background, the blue line is the signal shape after smearing, and the blue shaded region is the resulting signal plus background distribution.

��� ��� ��� ��� ��� ��� ���
We include a rough estimate of the signal yield of the S2-only analysis [43]. We consider two bins: (0.2, 0.5) keV and (0.5, 1) keV. We avoided considering the S2-only analysis beyond 1 keV to ensure that the dataset is independent from the one used to fit the signal. We impose a conservative bound by requiring a signal yield of less than 22 events in the (0.2, 0.5) keV bin, as well as less than 5 events in the (0.5, 1) keV bin. To evade the rising spectrum at low, sub-keV energy, we now turn our attention to exothermic DM, which can explain the XENON1T signal, and also has several interesting features that deserve more in-depth study.

Exothermic Dark Matter and Electron Recoils
DM could consist of two or more states with similar masses, see e.g. [20,21,[119][120][121][122][123][124]. We consider two states, χ 1 and χ 2 , with masses m χ 1 and m χ 2 = m χ 1 + δ, respectively, with |δ| m χ 1 and |δ| m χ 2 . For example, χ 1 and χ 2 could be two Majorana fermions that originated from a Dirac fermion that is charged under a new U (1) gauge symmetry; if there are mass terms for the Dirac fermion that break the U (1) symmetry, it is possible to split them into the two Majorana fermions, with the gauge boson coupling off-diagonally to χ 1 and χ 2 . Similarly, one can consider two real scalars that originated from a complex scalar. If χ 1 is the incoming state, then χ 1 can convert to χ 2 when scattering off ordinary matter, allowing for both the "inelastic" DM scenario [119] (δ > 0) and the "exothermic" DM scenario (δ < 0); the latter was discussed in the context of nuclear interactions in direct-detection experiments in [20,21].
The relic abundance of the two states depends on the precise model parameters. For δ sufficiently small (typically 2m e ), the lifetime of the heavier state for decays via the (off-shell) mediator into the lighter state plus two neutrinos, or for decays into the lighter state plus three photons, is easily much longer than the age of the Universe [122,123]. However, the fractional abundance of the heavier state after freeze-out in the early Universe will depend sensitively on the precise DM-mediator interaction strength and the DM and mediator masses [122,123]. For sub-GeV DM, the abundance of the heavier state will typically be small. However, even a small fractional abundance of the heavier state can leave dramatic signals in direct-detection experiments, since, as we will see, the mass splitting δ can be entirely converted into kinetic energy of the electron when scattering off an electron in a target material. The exothermic scenario allows all relic particles in the halo to scatter, while the inelastic up-scatter of the lighter to the heavier state will be highly suppressed for δ 10's of eV. Inelastic and exothermic DM scattering off electrons has not been discussed previously in the literature. We focus here on exothermic scattering, since it is able to explain the XENON1T excess. We provide here only a few of the salient features of exothermic scattering, focusing on the phenomenology and leaving to future work a more detailed study of both the exothermic and inelastic DM scenarios in the context of electron recoils.
We assume that the incoming DM transfers momentum q to the target. In contrast to Eq. (6.2), the energy-conservation equation now reads where ∆E e is again the energy transferred to the electron. Assuming a small mass-splitting compared to the mass scale of the DM i.e. |δ| m χ 1 ∼ m χ 2 , we can simplify this as, In contrast to the "standard" DM-electron scattering discussed above (δ = 0) (and in contrast also with exothermic nuclear scattering, see below), ∆E e can be well above the "typical" energy transfers of ∆E e ∼ 10 −3 q typ ∼ Z eff × few eV applicable for δ = 0. In particular, for δ ∼ O(−keV), the electron recoil spectrum will be peaked at O(keV), and can explain the XENON1T excess. The differential scattering rate is given by where the minimum velocity to scatter is given by Figure 18. Differential recoil spectra for "exothermic" DM, for two DM states that are split in mass by δ=−2.5 keV (left) and δ=−4 keV (right) for two DM masses, m χ =1 GeV (solid) and m χ = 1 MeV (dashed). We consider three DM form factors, F DM =1 (blue), F DM = (αm e /q) 2 (orange), and F DM = (q/αm e ) (green).
As there is an upper bound of v max =v esc + v E on the DM halo velocity, we get upper and lower bounds on the allowed values of q for a given m χ and a fixed ∆E e , Consider first DM masses of O(GeV). From Eq. (6.22), we see that the value of q max is near m χ v max ∼ O(MeV), which is much higher than q typ . Thus, the recoil spectrum in ∆E e depends on the behavior of q min as a function of ∆E e . From Eq. (6.21), we see that, for δ = 0 and ∆E e ∼ O(keV), q min q typ and therefore the integral over q misses the peak of the form-factor, leading to strongly suppressed scattering rates (as discussed above). However, for δ ∼ O(−keV), we see that q min = 0 when ∆E e = |δ|. For ∆E e smaller or larger than |δ|, q min increases and the available phase space decreases again, thus giving a suppression in the rate. Hence, for m χ ∼ O(GeV) and δ ∼ O(−keV), we get a sharp peak in the spectrum at ∆E e ∼ |δ|. Moreover, we see that the value of the atomic form factor actually increases towards larger ∆E e , and its peak shifts towards larger q (see Fig. 15). This means that after integrating over q, the total scattering rate is actually larger for δ < 0 than for δ = 0 (for fixed σ e ).
In Fig. 18, the solid lines show spectra for m χ = 1 GeV for δ = −2.5 keV (left) and δ = −4 keV (right) for f χ 1 σ e = 10 −40 cm 2 and for three different form factors. Here the fractional abundance of χ 1 is f χ 1 = nχ 1 nχ 1 +nχ 2 , with n χ 1 (n χ 2 ) being the number density of χ 1 (χ 2 ). We see that the spectrum is sharply peaked at δ and is reminiscent of a DM absorption signal, which provides an adequate fit to the XENON1T excess. Figure 19. The 1σ and 2σ best-fit regions that explain the XENON1T excess for exothermic DM and a heavy mediator (F DM = 1) in the f χ1 σ e versus δ plane for m χ = 0.8 MeV (f χ1 = nχ 1 nχ 1 +nχ 2 ) (left) , and in the f χ1 σ e versus m χ plane for δ = 4.2 keV (middle). In the right plot, we show an example of the predicted spectrum for the best-fit value from Eq (6.23). The dashed and solid lines show the signal spectrum before and after detector smearing effects, respectively. The measured XENON1T data is shown as black dots while the gray-shaded and blue-shaded regions are the expected binned background and signal respectively. We next consider m χ ∼ O(MeV). In this case, if δ ∼ O(−keV), q min = 0 at ∆E e = |δ|. However, q max is now only a few keV, and we see from Fig. 15 that the q typ for ∆E e ∼ O(keV) is higher. Note that for such small m χ , there is barely any kinetic energy in the DM to give recoil energies larger than |δ|, which is already around keV. So the spectrum sharply cuts off at ∆E e ∼ |δ|. For ∆E e < |δ|, both q min and q max increase, and at some ∆E e below |δ|, the allowed values of q cross the peak of the form factor. Hence, we see a peak in the spectrum for ∆E e slightly below |δ|. Moreover, the spectrum is not as sharply peaked as it is for heavier DM, since for heavier DM the allowed values of the momentum transfer are always q ∼ O(keV)∼ q typ . We see this in Fig. 18, where the dashed lines show spectra for m χ = 1 MeV for δ = −2.5 keV (left) and δ = −4 keV (right), both for f χ 1 σ e = 10 −40 cm 2 and for three different form factors. The spectrum has a wide peak, wider than for heavier exothermic DM and wider than a DM absorption signal, which (for the larger value of δ) provides a very good fit to the XENON1T data.
In Fig. 19, we show the 1σ and 2σ best-fit regions that explain the XENON1T excess for a heavy mediator (F DM = 1) in the δ-m χ 1 plane (left) and the σ e -m χ 1 plane (middle). We see that the best-fit point is given by (6.23) In the right plot, we show how the signal at the best-fit point compares with the XENON1T data and background model. In Fig. 20, we show the corresponding plots for a light mediator (F DM ∝ 1/q 2 ). Here the best-fit point is given by Figure 20. The 1σ and 2σ best-fit regions that explain the XENON1T excess for exothermic DM and a light mediator (F DM = (αm e /q) 2 )) in the f χ1 σ e versus δ plane for m χ = 3 MeV (f χ1 = nχ 1 nχ 1 +nχ 2 ) (left), and in the f χ1 σ e versus m χ plane for δ = 2.7 keV (middle). In the right plot, we show an example of the predicted spectrum for the best-fit value from Eq (6.24). The dashed and solid lines show the signal spectrum before and after detector smearing effects, respectively. The measured XENON1T data is shown as black dots while the gray-shaded and blue-shaded regions are the expected binned background and signal respectively.
Finally, in Fig. 21, we show the corresponding plots for a q-dependent heavy mediator (F DM ∝ q); here the best-fit point is given by (6.25) We see that exothermic DM can explain well the observed XENON1T ER spectrum. We now make a few comments: • The inclusion of relativistic corrections when calculating the atomic form factors is not essential for exothermic scattering, since q is not forced to be large to obtain a large ∆E e and the form factors typically peak at values of q below which relativistic corrections become important. We therefore neglect relativistic corrections in our calculations.
• While m χ ∼ 1 GeV provides an adequate fit to the XENON1T excess, one can obtain an even better fit for heavy DM by imagining that DM consists of three or more states. For example, for three states χ 1 , χ 2 , and χ 3 , with mass splitting δ 21 ≡ m χ 2 − m χ 1 , δ 31 ≡ m χ 3 −m χ 1 , and δ 32 ≡ m χ 3 −m χ 2 , with δ 21 , δ 31 , δ 32 all negative, the electron recoil spectrum would show up to three peaks. Of course, the actual size of the various peaks will depend sensitively on the relic abundances of the three DM states, and hence depend sensitively on the model parameters.
• If the DM couples also to nuclei (for example, if the mediator is a dark photon), DM could scatter exothermically off nuclei. We can contrast the kinematics for exothermic DM scattering off electrons with the kinematics for exothermic DM scattering off nuclei. For exothermic scattering off nuclei, the mean recoil energy is E R ∼ |δ|µ χ 1 ,N m N , where m N is the mass of the nucleus and µ χ 1 ,N is the reduced mass of χ 1 and the nucleus; the spread in energy around the mean recoil energy is given by [20,21]. Figure 21. The 1σ and 2σ best-fit regions that explain the XENON1T excess for exothermic DM and a momentum-dependent heavy mediator (F DM = (q/αm e )) in the f χ1 σ e versus δ plane for m χ = 4.5 MeV (f χ1 = nχ 1 nχ 1 +nχ 2 ) (left), and in the f χ1 σ e versus m χ plane for δ = 3 keV (middle). In the right plot, we show an example of the predicted spectrum for the best-fit value from Eq (6.25). The dashed and solid lines show the signal spectrum before and after detector smearing effects, respectively. The measured XENON1T data is shown as black dots while the gray-shaded and blue-shaded regions are the expected binned background and signal respectively.
For χ 1 scattering off a xenon atom, with m χ 1 ∼ 1 GeV and δ ∼ 1 keV, E R ∼ 8 eV, while the typical spread in energy around the mean recoil energy for the same parameters and a DM velocity of v ∼ 10 −3 is ∆E R ∼ 21 eV. This is below the XENON1T and many other experimental thresholds. The lowest-threshold search thus far for low-mass DM scattering off nuclei was achieved by CRESST-III; they obtained a threshold of 19.7 eV in [125] and a threshold of 30.1 eV in [126]. Therefore, there will be a constraint on DM with masses near the GeV scale that scatters exothermically off nuclei from CRESST-III; this bound will disappear for DM masses below ∼500 MeV.
• It is possible to obtain electron recoils from the Migdal effect when DM scatters exothermically off nuclei; this could lead to additional constraints, which requires a careful study that we leave to future work.
• As mentioned above, the fractional abundance of the heavier state after freeze-out in the early Universe will depend sensitively on the precise DM-mediator interaction strength and the DM and mediator masses. Moreover, in a concrete model there will typically also be other constraints from searches at beam dumps, fixed-target experiments, and colliders. We leave a detailed investigation of concrete models to future work.

Accelerated Dark Matter
A fraction of DM could be accelerated to high velocities, producing an energetic DM flux that impinges on the Earth [30-32, 41, 42]. Such an accelerated component may then be detected with experiments such as XENON1T, allowing for sensitivity to very light DM, which otherwise cannot be probed without sub-keV threshold experiments. Specifically, two distinct mechanisms have been suggested. In the first, DM interacts with the solar Figure 22. Accelerated dark matter flux due to interactions with cosmic ray electrons. The flux is shown for two different DM form factors: |F DM (q)| 2 ∝ q 2 /(q 2 + m 2 φ ) 2 (left) and |F DM (q)| 2 ∝ 1/(q 2 + m 2 φ ) 2 (right). The three different solid colored lines show the flux for varying values of the mediator mass, m φ =1 eV (blue), m φ =100 keV (orange), and m φ =100 GeV (green). In these plots the DM mass is set to 1 MeV and the DM-electron cross-section is taken to beσ e = 10 −30 cm 2 . The vertical dashed lines indicate the energy thresholds for the Super-K and XENON1T experiments.
interior to produce a significantly harder spectrum [41,42]. However, for standard DM with F DM = 1, the resulting flux ends at around 2 keV, thereby naively disfavoring a simple fit to the XENON1T data. A second energetic DM flux is generated through interactions with cosmic rays (CRs) [30][31][32]. This was used to derive world-leading limits on DM-electron couplings for DM in the eV to few keV mass range using the Super-K experiment [127].
Naively, DM acceleration from CRs cannot address the XENON1T anomaly either, for the following reason. For the previously studied DM-electron interactions with trivial form factor (F DM = 1), the predicted accelerated DM spectrum does not vary by more than an order of magnitude between keV and 100 MeV. However, the Super-K analysis uses 176 ktyears of data (with a ∼ 100 MeV threshold), about five orders of magnitude larger than the 0.65 tonne-year exposure available in XENON1T. Consequently, any signal at XENON1T would be naively excluded by Super-K.
The above argument does not hold for a DM interacting with electrons via a light mediator. Indeed, in such a case, both the produced flux and scattering rate predict a steeply falling spectrum towards higher energies, thereby easily compensating for the relative low exposure of XENON1T with its significantly lower threshold. However, experiments with lower thresholds may then be more constraining. We now study this possibility in detail.
We consider DM that interacts solely with electrons via a light mediator. For our purpose, in order to fit the XENON1T anomaly, the mediator mass must be lighter than a few MeV or else the benefit of having a low-threshold experiment in comparison to the Super-K experiment is lost. Furthermore, the mediator must also be heavier than roughly 1 keV in order to evade the S2-only analysis of XENON1T [43]. For such masses, stellar cooling Figure 23. Electron recoil spectra from cosmic ray accelerated DM flux for |F DM (q)| 2 ∝ q 2 /(q 2 + m 2 φ ) 2 (left) and |F DM (q)| 2 ∝ 1/(q 2 + m 2 φ ) 2 (right). Three different values of the mediator mass, m φ =1 eV (blue), m φ =100 keV (orange) and m φ = 100 GeV (green) are shown. The DM is fixed at MeV and the DM-electron cross section is taken to be 10 −30 cm 2 . As discussed in the text, only the |F DM (q)| 2 ∝ q 2 /(q 2 + m 2 φ ) 2 with an intermediate mediator mass can viably address the XENON1T data.
constraints are also irrelevant. A complementing study of this scenario with significantly lighter mediator masses is upcoming [34].
The DM flux obtained from interactions with CRs is given by, Here σ χe is the DM-electron (momentum-dependent) cross-section, ρ χ is the DM density profile, taken to be an NFW profile [128] with a scale radius of r s = 20 kpc and a local density of ρ = 0.4 GeV/cm 3 , and dΦ e /dE e is the CR-electron flux. In order to derive limits from the low threshold (∼150 eV) S2-only XENON1T analysis, one crucially needs to know the CR flux down to O(keV) energies. However, measurements only provide the spectrum down to MeV energies [129] and therefore an extrapolation must be used. Strictly speaking, this implies that systematic uncertainties hinder the possibility of using the S2-analysis to exclude the CR-accelerated DM solution. In what follows, we thus simply assume that the CR flux drops to zero below MeV energies.  Figure 24. Left: The 1σ and 2σ best-fit regions (red) for cosmic-ray accelerated DM that interacts with electrons via a form factor |F DM (q)| 2 ∝ q 2 /(q 2 + m 2 φ ) 2 . Exclusion regions due to the XENON1T S2-only analysis are shown by the light green-shaded region, while the limits from the Super-K data [32] are shown in darker shaded green. SN cooling strongly disfavors much of the parameter space due to constraints on the light mediator that couples to electrons. To illustrate this, in the middle plot we show the best fit region for the mediator-electron coupling y e as a function of the mediator mass, m φ . The blue region is excluded by HB cooling, while in the orange region the mediator thermalizes before electron decoupling, and may therefore suffer from (model-dependent) limits on N eff . The magenta shaded region is excluded by the (g − 2) measurement of the electron. Right: The signal spectral shape for the best fit point for this model. The black dots are the XENON1T data, the gray shaded region is the expected background, the blue solid line is the signal shape after detector smearing, and the blue dotted line is the signal before smearing. The blue shaded region is the resulting signal plus background distribution.
The model discussed here has three independent parameters: the DM mass, m χ , the mediator mass, m φ , and the DM-electron cross-section,σ e . In Fig. 22, we show the predicted accelerated DM flux for three different values of m φ and for two different form factors, fixing m χ = 1 MeV andσ e = 10 −40 cm 2 . We see that the expected flux is indeed significantly lower at the Super-K threshold for lighter mediators than for heavier mediators. We calculate next the expected electron spectra induced by the accelerated DM flux in xenon for the same two different form factors, and for the same values of m φ , with m χ = 1 MeV. The spectrum is shown in Fig. 23. For the |F DM (q)| 2 ∝ 1/(q 2 + m 2 φ ) 2 the predicted spectrum does not flatten out at low energies, and thus we find it to be excluded by the XENON1T S2-only analysis [43]. Conversely, for |F DM (q)| 2 ∝ q 2 /(q 2 + m 2 φ ) 2 , we see that we get a peak resembling the XENON1T excess for m φ ∼ 100 keV. This is because for high energies, the spectrum is suppressed due to the atomic form factor, while towards low energies (below the XENON1T excess), it falls due the q 2 suppression in the DM form factor. We therefore proceed with studying this case, fixing the mediator mass to be 100 keV as a benchmark value, and scanning over the DM mass and cross section.
In Fig. 24 (right), we show the expected signal for the best-fit value of the CRaccelerated DM with a mediator mass m φ = 100 keV. The dashed and solid blue lines show the unsmeared and detector smeared spectra. The gray region is the expected binned background while the blue shaded regions show the contribution of the binned signal. In the left plot, we show the corresponding 1σ and 2σ best-fit regions in red together with limits from the XENON1T S2-only analysis [43] and Super-K experiment [32] (shades of green). However, a vanilla mediator coupled to electrons with m φ = 100 keV, is strongly constrained by SN cooling [58]. In fact the relevant parameter space is subject to several constraints. We illustrate this in the middle of Fig. 24 where we plot the mediator's coupling to electrons, y e as a function of the mediator's mass. The blue region is excluded by HB while in the orange region the mediator thermalizes before electron decoupling and is thereby subject to (model-dependent) constraints from measurements of N eff . The magenta-shaded region is excluded by measurements of the g − 2 of the electron. Finally, we point out that in similar fashion to the Chameleon-like ALP of Sec. 5, here too it may be possible to evade stellar cooling constraints in a more sophisticated scenario. We postpone such a study for future work.
Before closing, a remark is in order. In deriving the limits above we used the nonrelativistic form factors for the DM-electron interactions. Corrections that arise in the relativistic limit can be found in [118]. To understand why it is justified to neglect relativistic corrections, we note that a DM with mass around 1 MeV must be accelerated to velocities above ∼ 0.03c. Using Eq. (6.11) and since E e keV, one finds q 30 keV. Since the atomic and DM form factors are both dominated at low q, this justifies neglecting the relativistic corrections, which become important only at significantly higher values of q.