Solar axions cannot explain the XENON1T excess

We argue that the interpretation in terms of solar axions of the recent XENON1T excess is not tenable when confronted with astrophysical observations of stellar evolution. We discuss the reasons why the emission of a flux of solar axions sufficiently intense to explain the anomalous data would radically alter the distribution of certain type of stars in the color-magnitude diagram in first place, and would also clash with a certain number of other astrophysical observables. Quantitatively, the significance of the discrepancy ranges from $3.3\sigma$ for the rate of period change of pulsating White Dwarfs, and exceedes $19\sigma$ for the $R$-parameter, which measures the ratio between the numbers of Horizontal Branch and Red Giant Branch stars in globular clusters.

Introduction. The XENON1T collaboration [1] has reported an excess in low-energy electronic recoil data below 7 keV, rising towards low energies with a prominent peak around 2-3 keV. The collaboration cautions that the excess could be explained by β decays from a trace amount of tritium, representing an unaccounted source of background, but they also explore the possibility that the signal is due to different types of new physics. The most intriguing interpretation, which also provides the best fit to the data, is given in terms of solar axions, favoured over the background-only hypothesis at the 3.5σ level.
There are basically three production mechanisms that contribute to the total solar axion flux: i) Atomic recombination and deexcitation, Bremsstrahlung, and Compton (ABC) interactions [2] that are controlled by the axion-electron coupling g ae , ii) Primakoff conversion of photons into axions [3] induced by the axion-photon coupling g aγ , iii) axion emission in the M1 nuclear transition of 57 Fe [4] that produces mono-energetic 14.4 keV axions, and is controlled by and effective axion-nucleon coupling g eff an . Since this last process cannot play any role in producing events below 10 keV, we will not include in our analysis astrophysical observables sensitive to g eff an . Conversely, axions produced through i) and ii) feature a wide spectrum peaking around a few keV. The production rates are independent of the axion mass for m a < ∼ 100 eV. As regards detection, electron recoils occur via the axioelectric effect which is controlled by g ae . Because of this, and because the location of the peak around 2-3 keV corresponds roughly to the maximum of the axion energy spectrum for the ABC processes, the Primakoff and 57 Fe components are both allowed to be absent as long as there is a nonzero ABC component. This selects g ae as the crucial coupling to attempt to explain the data in terms of the QCD axion [5][6][7][8]. Taken at face value the strength of the XENON1T excess requires g ae > ∼ 10 −12 , corresponding to an axion decay constant f a < ∼ 10 8 GeV, and in turn to an axion mass m a > ∼ 0.06 eV. However, astrophysical considerations indicate that such a large value of g ae is not tenable. While we will eventually quantify in terms of standard deviations the discrepancy between the solar axion interpretation of the XENON1T excess and astrophysics observations, it should be kept in mind that stellar evolution would be so drastically affected by the exceedingly large energy losses via axion emission, that the nice agreement between stellar evolution codes and the observed color-magnitude diagram (CMD) for clusters of stars would be destroyed. The strategy that we will follow consists in assuming that g ae and g aγ lie in the 90% C.L. regions resulting from the XENON1T fit [1]. We will then estimate the effects of extra energy losses on a set of astrophysical observables related to Red Giants Branch (RGB) and Horizontal Branch (HB) stars, and to White Dwarfs (WDs), that are particularly sensitive to axion emission via ABC and Primakoff processes.
Astrophysical observables and axion couplings. The axion interactions with photons and electrons read where the couplings can be related to model-dependent dimensionless coefficients as g aγ = α 2π Caγ fa and g ae = C ae me fa . In benchmark axion models C aγ and C ae are typically of O(1), although a strong enhancement/suppression is possible in specific cases [9][10][11][12][13]. In the following, we will adopt the notation g γ10 ≡ g aγ × 10 10 GeV and g e13 ≡ g ae × 10 13 (g e12 ≡ g ae × 10 12 ).
Astrophysical considerations have been systematically used to place severe bounds on light and weakly interacting particles, such as neutrinos and axions [14]. Noticeably, a set of anomalous observations have recently led to speculations that new physics is at play [13,15,16], and the axion case appears especially compelling [17] since it fits particularly well these observations [18]. The observables that are most effective to constrain g ae and g aγ and on which we focus our analysis are described below.
• Tip of RGB stars in globular cluster. We denote by M I,TRGB the luminosity of the tip of the RGB in glob-ular clusters. RG stars are characterized by a He core and a burning H shell. During the RGB evolution, the ashes from the burning shell increase the He core mass, while the star luminosity (determined by equilibrium at the surface of the He core between thermal pressure supporting the non-degenerate envelope against the gravity pull from the core) keeps growing. The process continues until the core reaches sufficiently large temperatures and densities (T ∼ 10 8 K, ρ = 10 6 g cm −3 ) to ignite He, an event known as the He-flash. At this stage the star has reached the maximum luminosity, that is the RGB tip, after which it shrinks and moves to the HB. If an additional core-cooling mechanism is at play, He ignition is delayed, the core would accrete a larger mass, and the star would reach higher luminosities. Therefore, measurements of M I,TRGB allow to test the rate of cooling during the RGB phase. The method is particularly effective for constraining g ae since in RG cores axions can be efficiently produced via electron bremsstrahlung. The most recent analyses [19][20][21] have derived comparable constraints.
Here we adopt the result of the analysis of the M5 cluster (NGC 5904) in Ref. [19] which provides the most conservative bound M I,TRGB = −4.17 ± 0.13 mag. In terms of g ae this observable can be written as [19] M theo which results from an analytic fit to ten evolutionary track points reaching close to the RGB tip obtained from numerical sinulations, and corresponding to different values of g ae up to g e13 = 9 [19,22]. The associated theoretical uncertainty is σ 2 = 0.039 2 + (0.046 + 0.012g e13 ) 2 .
• R-parameter. After He ignition the RG core expands and its density decreases by two orders of magnitude. The star migrates to the HB and remains supported by He burning in a non-degenerate core. The ratio R = N HB /N RGB between the number of stars in in globular clusters in the HB and in the upper portion of the RGB directly measures the duration of He burning in the HB phase. The value R = 1.39 ± 0.03 was obtained in Ref. [23] from the analysis of 39 clusters.
The duration of the HB phase can be affected by g aerelated processes both directly and indirectly. If g ae is sufficiently large, axion emission would directly produce extra cooling of the He core. The star self-regulates by slightly contracting, the core temperature increases speeding up the He burning rate. Once the He fuel is exhausted, the star turns into a WD. The indirect effect is related to the growth of the degenerate He core during the RGB phase previously discussed. HB stars would unavoidably inherit a more massive core from the parent RGs, resulting in am increased He burning rate to contrast the larger gravitational pull, and shortening further the duration of the HB phase. Note that the indirect effect of g ae is so important that for g e13 ∼ 15 it would suffice to depopulate almost completely the HB in the CMD (R ≈ 0). Cooling of HB stars can also proceed via the Primakoff effect ∝ g 2 aγ , which is particularly efficient at the typical temperatures and densities of HB cores (T ∼ 10 8 K, ρ = 10 4 g cm −3 ), For sufficiently large values of g aγ , R can still decrease well below the observed values even when g ae ≈ 0. Hence an accurate determination of this observable allows to probe the axion coupling to both photons and electrons. In terms of g ae and g aγ the R-parameter can be written as [18,24] R theo = 7.33Y − 0.095 21.86 + 21.08g γ10 δM c = 0.024 where δM c is the change in the He-core mass, and Y 0.255 ± 0.002 is the primordial He abundance. Similarly to Eq. (2) this expression is derived from an analytic fit to evolutionary track points calculated with stellar evolutionary codes modified to account for axion emission [18,24], and thus it cannot be reliably extrapolated to values of g ae and g aγ much larger than those corresponding to the last point fitted. For definiteness, Eq. (2) and (for g γ10 < ∼ 1) Eq. (3) can be trusted up to g e13 ∼ 9. Thus, we will not input into these expressions the much larger XENON1T values g e13 ∼ 30. Rather, very conservatively, we will limit ourselves to estimate the tension between the observed values of M I,TRGB and R, and the values resulting from Eqs. (2) and (3) when evaluated at g e13 ∼ 9 (g aγ ≈ 0). As regards values of g aγ too large to be used in Eq. (3), they can be directly constrained from the lifetime of HB stars which, in the presence of extra cooling, scales as ∼ L 0 /(L 0 + L a ) with L 0 (L a ) the standard (the axion) core luminosity [14]. Hence for g γ10 > ∼ 1, rather than Eq. (3), we will instead use with a = (6.26Y − 0.12) and b = 0.41 [23]. Note that Eq. (4) neglects both direct and indirect effects of g ae on HB and on RGB stars, and hence it would also yield conservative limits.
• White Dwarf luminosity function. The third observable we consider is the distribution of WDs as a function of their luminosity. The WD luminosity function (WDLF) has been used for decades as a tool to measure the WD cooling efficiency, since the WD number distribution depends on the speed of their evolution. Observations of the WDLF thus allow to place strong bounds on new exotic cooling processes, including axion emission (see Ref. [25] for a recent review). WDs are compact objects whose hydrostatic equilibrium is supported by electron degeneracy pressure, hence axion emission from WDs would dominantly depend on g ae . Here we will use the bound g WDLF e13 ≤ 2.8 obtained in Ref. [26].
• Rate of period change of WD variables. WD variables (WDV) are WDs whose luminosity varies periodically, with a period Π ranging from a few to several minutes, depending on the particular star. Because the oscillation period depends on the luminosity, an observed secular changeΠ of the period tracks the rate of decrease of the star luminosity. To a very good approximationΠ/Π is directly proportional to the cooling rateṪ /T , hence an accurate measurement ofΠ allows to set bounds on possible sources of extra cooling (see Ref. [27] for a comprehensive review). Here we consider four WDVs: G117-B15A [28], R548 [29], L19-2 [30] (for whichΠ has been measured for two different pulsation modes), and PG1351+489 [31]. We list in Table I the corresponding measured values oḟ Π/Π. Theoretically, the rate of change in the WD pulsating period as a function of g e13 can be parametrized as [17]:Π theo WDi = a i + b i g 2 e13 , where a i and b i are parameters specific for each WD.
XENON1T vs. Astrophysics. Fig. 1 shows contours of the axion energy-loss rates per unit mass in a temperature vs. density plane, for a pure He plasma. Contour iso-lines for energy loss due to Compton (dashed blue) and Bremsstrahlung (solid red) processes, which are controlled by g ae , are also shown. For reference, we have fixed the axion-electron coupling to g e13 = 4.3, which corresponds to the RGB bound from M5 [19]. Energy loss rates for different values of the couplings can be obtained recalling that the rates scale as g 2 ae . The labelled disks in the figure show the position of the RGB tip and of a typical HB star (of mass 0.8M ) and a range of WDs with luminosities varying roughly from 5 × 10 −4 to 5 × 10 −1 L (dashed gray rectangle). The blue disk indicates the temperature/density of a typical WD variable [27]. The location of the Sun with respect to the temperature/density is marked with a yellow disk on top of the broken gray line which locates Main Sequence (MS) stars of different masses. Note that since MS stars, including the Sun, are supported by H burning cores, their position with respect to the energy loss iso-lines for the He plasma is approximate, and slightly shifted towards larger rates. The picture shows clearly that the Sun is a relatively faint axion emitter with respect to other stellar objects, so that values g ae > ∼ 10 −12 as required to account for the XENON1T excess would unavoidably turn other stars into bright 'axion lighthouses'. The RGB would extend to higher luminosities than the ones observed, and the decreased duration of the He burning phase would depopulate the HB, to the point that for smaller clusters with relatively few stars, already for g e13 ∼ 15 we would expect R ≈ 0. In short, regardless of other details, a value g e13 ∼ 30 would definitely destroy the agreement between stellar evolution models and the observed CMD.
Quantifying the tension. The projections of the XENON1T 90% C.L. best fit region onto the (g ae , g aγ ), Contours of the axion energy-loss rates per unit mass, in erg g −1 s −1 , for a pure He plasma. We highlight different stellar systems as well as the two main processes involving the axion coupling to electrons, fixed to gae = 4.3 × 10 −13 .
(g ae , g eff an ) and (g ae g aγ , g ae g eff an ) planes are given in Fig. 8 of Ref. [1]. Since only g ae and g aγ can be responsible for the anomalous XENON1T data below 7 keV, we focus on the best fit region for these two couplings, which is reproduced in Fig. 2 by the blue band. In the figure we also show the 2σ limits on g ae , g aγ obtained from the single astrophysical observables, as well as the result of a global fit to the full set of stellar cooling data. The curve depicting the CAST [32] limit in the (g ae , g aγ ) plane in Ref. [1] was taken from Ref. [33]. We update this bound with the most recent CAST results from Ref. [34] which, in the g ae 0 limit, and for m a < ∼ 20 meV (m a < ∼ 0.7 eV), correspond to g aγ < 0.66 (2.0) × 10 −10 GeV −1 . These limits are represented in Fig. 2 by the two green lines, in which we have folded in the effects of a non-zero g ae that would increase the production of solar axions and strengthen the bounds. The vertical dashed line is LUX limit [35], while The grey horizontal line at g γ10 = 4.1 corresponds to the limit from a global fit to solar data, which includes the measured flux of 8 B and 7 Be neutrinos as well as additional data inferred from helioseismology observations [36]. This is about a factor of two stronger than the bound labeled "solar ν" in the upper panel of Fig. 8 in Ref. [1] which is taken from Ref. [37]. To assess quantitatively the discrepancy between the values of g ae and g aγ needed to reproduce the XENON1T excess we proceed as follows: we first extract the allowed ranges from the 90% C.L. region of Ref. [1] falling in the region allowed by solar data (cf. the blue area in Fig. 2). This region can be parametrized by means of an effective FIG. 2. XENON1T 90% C.L. fit (blue region). 3σ exclusion limit from solar data (grey hatched region). 2σ LUX limit (grey dashed line) and CAST limits for ma < 20 meV and ma < 0.7 eV (green lines). Individual 2σ limits from R-parameter, TRGB, WDLF, WDVs (grey lines) and global bounds from stellar cooling: 1σ (dark red) to 4σ (light red).  coupling [13] The 90% C.L (68% C.L.) region of XENON1T is then well represented by the range g e12 ∈ [2.6, 3.7] (g e12 ∈ [2.8, 3.5]). Varying g ae and g aγ subject to the condition that g e12 remains within this range, we derive the range of expected values for the astrophysical observables implied by the XENON1T data, and we confront them with the measured values. Our results are collected in Table I where we also give, for each observable, the tension be-tween the values expected from the fit to XENON1T data and astrophysical observations. It can be seen that the large values of g ae required to fit the XENON1T excess are in strong conflict with all the astrophysical observables. The discrepancy is at the level of ∼ 3.4σ for the WDVs cooling rates (last five rows in the Table), and reaches ∼ 6σ for the WDLF in the third row. As regards the first two rows, we recall that the values of R theo and of M theo I,TRGB are obtained respectively from Eq. (3) and Eq. (2) by setting g e13 = 9, rather than by inserting the much larger values g e13 ∼ 30 needed to account for the XENON1T data. As explained above this is because for such large values of g ae we cannot trust the theoretical parametrizations in Eqs. (2) and (3). In the Table we have then marked with a the corresponding tensions, since sizeably larger tensions can be expected for values of the observables in agreement with the XENON1T solar axion fit. To give an example, if we assume that Eq. (3) can be trusted up to g e13 ≈ 15, we would obtain R ≈ 0, corresponding to a complete depopulation of the HB, and 46σ away from observations.
Conclusions. In this work, we have explained why astrophysical observations firmly exclude that solar axions could account for the XENON1T excess. Other explanations based on solar production of new light particles or on modifications of neutrino properties (such as a neutrino magnetic moment) are also prone to severe astrophysical constraints, and as long as the corresponding new physics processes would also occur in RG, HB and WD stellar cores, they can likewise be excluded. However, other mechanisms involving either absorption or scattering of new particles of non-solar origin off target electrons [38][39][40][41][42], although less compelling than the QCD axion, might still provide viable explanations for the XENON1T data.
Note added. After completing this letter, Refs. [43,44] appeared claiming that besides the axio-electric effect, also the inverse Primakoff process can contribute to the detection of solar axions by XENON1T. This would relax the best fit region towards lower values of g ae at the cost of increasing g aγ which, however, can still remain below the solar data limit. Although the tension with astrophysical bounds would get relaxed, using the results of Ref. [43,44] we have verified that the discrepancy with the R-parameter remains quite large, at least at the level of 8σ.