Direct detection of dark energy: the XENON1T excess and future prospects

We explore the prospects for direct detection of dark energy by current and upcoming terrestrial dark matter direct detection experiments. If dark energy is driven by a new light degree of freedom coupled to matter and photons then dark energy quanta are predicted to be produced in the Sun. These quanta free-stream towards Earth where they can interact with Standard Model particles in the detection chambers of direct detection experiments, presenting the possibility that these experiments could be used to test dark energy. Screening mechanisms, which suppress fifth forces associated with new light particles, and are a necessary feature of many dark energy models, prevent production processes from occurring in the core of the Sun, and similarly, in the cores of red giant, horizontal branch, and white dwarf stars. Instead, the coupling of dark energy to photons leads to production in the strong magnetic field of the solar tachocline via a mechanism analogous to the Primakoff process. This then allows for detectable signals on Earth while evading the strong constraints that would typically result from stellar probes of new light particles. As an example, we examine whether the electron recoil excess recently reported by the XENON1T collaboration can be explained by chameleon-screened dark energy, and find that such a model is preferred over the background-only hypothesis at the $2.0\sigma$ level, in a large range of parameter space not excluded by stellar (or other) probes. This raises the tantalizing possibility that XENON1T may have achieved the first direct detection of dark energy. Finally, we study the prospects for confirming this scenario using planned future detectors such as XENONnT, PandaX-4T, and LUX-ZEPLIN.


I. INTRODUCTION
More than two decades after the discovery that the expansion of the Universe is accelerating [1,2], the nature of the dark energy (DE) component driving this phenomenon and making up ∼ 70% of the energy budget of the Universe remains a mystery [3][4][5][6][7]. What is perhaps the simplest theoretical DE candidate, a cosmological constant (CC) resulting from the collective zero-point energy of quantum fields, suffers from a severe disagreement between its theoretical value suggested from quantum field theory considerations, and the tiny value inferred from cosmological observations [8][9][10]. This staggering discrepancy goes under the name of cosmological constant problem [9,10]. Several DE models have been proposed as alternatives to the CC, within most of which the cosmic acceleration phenomenon is ascribed to additional (typically scalar) degrees of freedom, either in the form of new fundamental particles and fields, or modifications to General Relativity (GR). Competing models predict a variety of novel cosmological phenomena which will be actively sought for by the next generations of cosmological surveys, including but not limited to the Simons Observatory [11,12], CMB-S4 [13,14], Euclid [15,16], DESI [17], the Vera C. Rubin Observatory Legacy Survey of Space and Time [18], and the Nancy Grace Roman Space Telescope [19].
In all but the simplest models, the scalar field leads to a modification of gravity where the associated new light scalar degree of freedom couples to matter [20][21][22]. Such theories predict the existence of fifth forces that are ostensibly excluded by solar system tests of GR [23,24]. To circumvent this problem, realistic models inevitably need to include some form of screening mechanism [6,[25][26][27][28]. Screening mechanisms dynamically suppress fifth forces in the solar system without the need to tune the model parameters, resorting to environmental effects in the case arXiv:2103.15834v3 [hep-ph] 23 Aug 2021 of several (but not all) of these mechanisms. This then allows for strong deviations from GR on cosmological scales, which could modify the growth of structure. The ability to account simultaneously for cosmic acceleration and satisfy solar system tests of GR have made modified gravity (MG) theories with screening mechanisms leading science targets for current and upcoming cosmological surveys [29]. In addition to their indirect effects on the cosmological background expansion and structure formation [30,31], 1 DE and MG theories with screening mechanisms have the attractive feature that they are also amenable to direct detection of their associated effects. The distinctive fifth forces they predict can be searched for in laboratory experiments [29,37,38] and astrophysical environments [27,28,[39][40][41][42][43][44][45][46][47][48][49][50][51], and their observation would unequivocally confirm the hypothesis that DE is linked to a modification of GR.
The purpose of this paper is to broaden the scope of new physics accessible to terrestrial dark matter (DM) direct detection experiments by exploring their potential to detect DE quanta produced in the Sun, and to open up a new frontier for the direct detection of dark energy. Screening mechanisms suppress the production of DE quanta in the core of the Sun due to the large scalar mass or weak coupling to matter, but they can be produced in regions of strongly magnetised plasmas through couplings of the scalar to photons. In particular, they can be produced in the solar tachocline, 2 in a manner analogous to the Primakoff process for axions [53,54]. This possibility opens up yet another exciting avenue for directly detecting DE fluctuations: DE particles produced in the solar tachocline could be observed in dark matter (DM) direct detection experiments. To date, searches for DE-like particles in direct detection experiments have mostly focused on axion-specific experiments such as the CERN Axion Solar Telescope (CAST) [55][56][57][58] and the Axion Dark Matter Experiment (ADMX) [59]. Other detection techniques currently being considered are at different stages of realisation, or have started to gather data [60][61][62][63][64][65][66]. For additional details, see Refs. [67][68][69][70].
All well-known screening mechanisms can be classified into chameleon [71], symmetron [72], , K-mouflage [74], and Vainshtein [75] screened theories. In order to encapsulate the diversity exhibited by screening mechanisms, we will use an effective theory for the fluctuations of the scalar about the background profile due to the environment that includes operators relevant for each mechanism, including generic couplings to the standard model (SM). Using this, we 1 See e.g. Refs. [32][33][34][35][36] for recent works studying the effect of direct couplings of DE to baryons on both local and cosmological observables. 2 The tachocline is the turbulent shear layer located at the base of the solar convection zone, marking the transition between the radiative interior and the differentially rotating convective zone [52]. The radial position of the tachocline is approximately 0.7R , with R the solar radius.
will lay out the formalism for computing the expected signal in the associated detectors following from the scattering/absorption of DE scalars produced in the solar tachocline by nuclei or electrons (depending on the specific detector details). We will assume that the DE scalar possesses a coupling to photons, allowing for its production from photons in the tachocline. We will also use the environmental dependence of the couplings in the effective theory for the scalar fluctuations. We note that screening mechanisms of the chameleon, symmetron and Damour-Polyakov types show a clear dependence of these coupling on the local density, i.e. on the environment. On the other hand, screening mechanisms of the K-mouflage and Vainshtein types depend on more global features of the matter distribution, such as the total mass of the stars, with little or no dependence on the local matter density [76]. In order to provide a well-studied and well-constrained example, we will later specialise to the chameleon mechanism. Theories that screen using the chameleon mechanism [29,71,77,78] do so by increasing the mass of the DE scalar in dense environments, making their fifth force too short-ranged to be relevant in the solar system or in terrestrial fifth force searches. As a case study, we will apply our methodology to the XENON1T DM direct detection experiment [79], which recently reported a ∼ 3.3σ excess in the electron recoil data above their expected background, in the energy range 1 − 7 keV [80]. The XENON1T collaboration finds that a fit to the signal which includes a solar axion component is preferred over the background-only hypothesis at a significance of 3.4σ. Alternative explanations proposed by the collaboration include an additional tritium background and a neutrino magnetic moment, both of which are preferred with a significance of 3.2σ. Unfortunately, for the solar axion (and, to a minor but still important extent, for the neutrino magnetic moment), the parameters that best fit the XENON1T signal are strongly excluded by astrophysical observations of stellar evolution [81], particularly by constraints from the cooling of horizontal branch stars in globular clusters, the cooling of white dwarfs, and the tip of the red giant branch I-band magnitude of globular clusters and galaxies. The difficulties faced by these three scenarios have spurred the proposal of various new physics models that might also be able to explain the XENON1T excess, with varying degree of motivation or plausibility: for a list of works in this direction, we refer the reader to e.g. Refs. .
Applying our formalism to XENON1T, we find that chameleon-screened DE provides a fit to the signal of only slightly lower quality than the aforementioned scenarios, and is preferred over the background-only hypothesis at a significance of 2.0σ. The stellar bounds that are debilitating for the solar axion interpretation of the XENON1T signal are avoided by chameleon-screened DE theories due to the associated fifth force screening, as well as by the dependence of the mass of the scalar field on the local density of matter [53,54]. In particular, the astro-physical objects from which the strongest axion bounds derive -red giant, horizontal branch, and white dwarf stars -have cores whose density is larger than that of the solar tachocline by a factor of > 10 6 . As a consequence, chameleon-screened DE scalars are too heavy to be thermally produced within these objects, and cannot act as a new source of energy loss. This raises the tantalizing possibility that XENON1T, originally devised to detect DM, may instead have detected dark energy, or in any case a propagating scalar degree of freedom of a theory of modified gravity.
The possibility that chameleon-screened scalars might indeed be at the origin of the XENON1T signal is independently testable by a number of upcoming lowthreshold DM direct detection experiments. These experiments are all sensitive to recoil energies of O(100 eV) or lower, and are expected to detect the chameleoninduced signal at much higher statistical significance than XENON1T. Examples of experiments in this class, some of which make use of cryogenic semiconductor detectors, include SuperCDMS [136], CDEX-10 [137], DAMIC [138], CRESST-III [139], PICO [140], LUX-ZEPLIN [141], EDELWEISS [142], SENSEI [143], PandaX-II [144], PandaX-4T [145], and of course XENONnT [146] among others. We will demonstrate that these experimental setups, often discussed in the context of searches for light (sub-GeV or even MeV) DM scattering off nuclei or electrons, are equally well-placed to detect the signal of solar chameleons. 3 The rest of this paper is organized as follows. Our effective theory for screened dark energy in the solar system is presented and discussed in Sec. II. The specific chameleon model we specialise to in this work is presented in Sec. III, where we calculate its flux from production in the solar tachocline, and detection cross-section in DM direct detection experiments. In Sec. IV we discuss the data analysis method we use to confront our model against the XENON1T signal. In Sec. V we present the results following from this analysis. We critically discuss these results and prospects for probing chameleons in current and future DM direct detection experiments in Sec. VI. Finally, we draw our concluding remarks in Sec. VII. In addition, our paper contains three technical appendices. Appendix A is devoted to a more detailed discussion of the chameleon mechanism. Appendix B revises the production of chameleon-screened DE scalars in the Sun, as well as the resulting flux and energy spectrum of solar chameleons on Earth. Appendix C presents details concerning the computation of the cross-section for what we refer to as the "chameleo-electric effect", a process which is analogous to the photoelectric and axioelectric effects for photons and axions respectively, and which is relevant for the computation of the expected 3 In the following we will use the word chameleon for all chameleonscreened scalars. When the original chameleon scalar with an inverse power potential is meant, this will be explicitly specified.
signal in the XENON1T detector. All the codes associated with this work are made publicly available online at github.com/lucavisinelli/XENONCHAM.

II. DARK ENERGY EFFECTIVE THEORY IN THE SOLAR SYSTEM
Here, we work under the assumption that dark energy arises as a result of the cosmological dynamics of a single scalar field, which we denote by ϕ. This scalar does not have to be a fundamental field. Instead, it could arise as a low energy degree of freedom emerging from more involved dynamics at higher energy. For example, the scalar could be the Stückelberg field of the broken time-diffeomorphism symmetry of the cosmological background [147]. Our further assumption is that this scalar is involved in generating the acceleration of the Universe and at the same time couples to gravity and/or the SM in a manner that is ghost-and pathology-free, as explicitly realised e.g. in models of the Horndeski [148], beyond Horndeski (GLPV) [149,150], or Degenerate Higher-Order Scalar-Tensor (DHOST) classes of scalar-tensor theories [151,152]. These theories are among the leading candidate DE theories accompanying a modification of gravity (see Refs. [3, 6, 20-22, 31, 153] for more general reviews concerning MG theories and cosmological tests thereof).
Multi-messenger astronomy involving a gravitational wave (GW) and an associated optical counterpart can be used to constrain the properties of theories of modified gravity. To date, there has been one such confirmed event. In 2017, the LIGO/Virgo collaboration observed the gravitational wave event GW170817 [154], resulting from a binary neutron star merger. Simultaneously, the Fermi Gamma-ray Burst Monitor and the INTEGRAL Anti-Coincidence Shield spectrometer detected the short gamma-ray burst GRB170817A [155,156], which was identified as being the electromagnetic counterpart to GW170817. The joint GW170817/GRB170817A detection restricts the speed of GWs to differ from the speed of light by no more than one part in 10 15 , setting strong constraints on theories of the Horndeski [157][158][159][160][161][162][163], beyond Horndeski [164,165], and DHOST types [166][167][168][169][170], although these constraints are subject to important caveats [171,172]. Other important constraints on such theories arise from potential instabilities in GW backgrounds [173,174], astrophysical bounds [27,28,[175][176][177], and cosmological constraints [178][179][180][181]. Nevertheless, large regions of parameter space remain observationally and theoretically viable [182][183][184][185][186]. Since we will work at the level of the effective theory in the solar system, our formalism is insensitive to the details of the fundamental theory responsible for driving cosmic acceleration. We will thus assume that the aforementioned bounds are satisfied. The chameleon-screened dark energy theories that we study in this paper predict that the speed of gravitational waves is luminal so the bounds are automatically satisfied without the need to tune parameters. Attempts to embed our more general solar system effective theory into a covariant model which can be extended to cosmological scales should ensure that these bounds are satisfied. We defer this study to future work.
We begin by expanding the field ϕ( x, t) around some background value ϕ 0 ( x, t) as ϕ( x, t) = ϕ 0 ( x, t) + φ( x, t), where φ( x, t) is a spacetime-dependent perturbation. Depending on the nature of the fundamental theory, ϕ 0 ( x, t) could either be the field in the cosmological background, the field sourced by the Milky Way, or even different in the Earth and the Sun. The latter scenario is realised by the chameleon screening mechanism. The relevant part where X = −g µν ∂ µ φ∂ ν φ and F µν is the photon fieldstrength tensor. The sum in Eq. (3) includes photons, but since their energy-momentum tensor is traceless they do not contribute to the first two terms. A direct coupling to photons through their field-strength tensor can arise through quantum anomalies [187], and we have therefore included this coupling with strength β γ . This term is of critical importance since it allows for the production of scalars in the solar tachocline. The couplings to matter arise from the Jordan frame metric, i.e. the metric coupling SM particles to the scalar when the graviton is canonically normalised in the underlying covariant theory: From a quantum field theory perspective, it is unlikely that the couplings to each matter species are universal, hence our choice to treat the couplings β i , c i , and d i as being species-specific. The term multiplying g µν is referred to as the conformal factor and the term multiplying ∂ µ φ∂ ν φ is referred to as the disformal factor. We therefore refer to the β i s and c i s as conformal couplings, and d i s as disformal couplings (c i can also be referred to as kinetic-conformal coupling). The effective theory in Eq. (3) includes environmental variations of the various coupling constants via their dependence on the background field ϕ 0 , and therefore on the local matter density when the energy-momentum tensor of matter is dominated by non-relativistic species, such as in the late-time Universe after matter-radiation equality, or in astrophysical situations through the virialised matter density. In particular, we have allowed the kinetic matrix Z µν (ϕ 0 ) to be non-diagonal and to depend on the background field. This structure typically emerges from ghost-free higher-derivative couplings of the field to curvature tensors, and gives rise to the Vainshtein [75] and K-mouflage [74] screening mechanisms. We have also allowed the mass to be background field-dependent. This gives rise to the chameleon mechanism [71,77]. Note that m(ϕ 0 ) is the mass of fluctuations about ϕ 0 . The mass of ϕ in the cosmological background is instead O(H 0 ), so that ϕ can act as a DE scalar. Finally, we have allowed the coupling constants β i , c i , and d i to be fielddependent too. This field-dependence is utilised by the symmetron [72] and dilaton [188] mechanisms.
In this work, we shall focus on theories that utilise the chameleon mechanism, and therefore set Z µν (ϕ 0 ) = η µν , while taking β i , c i , and d i to be background fieldindependent. We also introduce the energy scales M i ≡ M/d The mechanism by which DE is produced in the Sun and either scatters or is absorbed in DM direct detection experiments is the following. The Sun is necessarily screened, implying that the mass of the scalar in the Sun is m = m(ϕ ) > 10 3 H 0 [189,190]. In practice, the high density of the Sun with respect to the cosmological background (ρ /ρ c ∼ 10 29 ) implies that m is in fact much heavier than this deep inside the Sun. The exact value is model-dependent, as it depends on the exact density dependence of the scalar's mass. The high mass prevents the efficient production of chameleons via Compton or bremsstrahlung processes in the Sun's core via a Boltzmann suppression. However, the direct coupling to photons allows for production in the magnetic field of the tachocline via a mechanism analogous to the Primakoff process for axions [191]: we review the chameleon production process in Appendix B. The relevant operators for this production process are: where M γ = M/d 1/4 γ is the energy scale related to the disformal coupling to photons. As chameleons are produced in the tachocline and not in the core, we need to impose that the mass of the chameleon in the core be larger than the local temperature. This can be achieved using the density-dependence of the mass, as the ratio between the densities in the core and in the tachocline is around two orders of magnitude.
Once produced, solar chameleons free-stream out of the Sun, with a fraction passing through the Earth, and an even smaller fraction through the chambers of DM direct detection experiments. In these chambers the mass is mostly determined by the detector's geometry [71,77] and is typically small (m DC m e , where m e is the mass of the electron). Chameleons can therefore be treated as massless particles in the chambers of DM direct detection experiments for all intents and purposes. Chameleons passing through the detector chamber can scatter off or be absorbed by the particles utilised for the detection, via couplings of the form These couplings give rise to what we refer to as the "chameleo-electric effect", and correspondingly to electron recoils in the O(keV) range. Similar interactions with neutrons and protons will instead give rise to atomic recoils, which we will not explicitly consider in this paper.
In the subsequent analysis, we will neglect the coupling c i , which controls the strength of the kinetic-conformal coupling XT i . We find that this coupling has an effect similar to that of the conformal coupling b i at the level of detection. Moreover, from a statistical perspective, it does not lead to a substantial improvement in the fit to the XENON1T signal (which would otherwise warrant its inclusion as a free parameter), as the disformal detection channel dominates over the conformal one(s), for reasons we will discuss in Sec. IV below. The choice of setting c i = 0 for the purposes of this work should be viewed as a simplifying working assumption, which sets a minimal phenomenologically working model: (re)-including c i would not change our results, nor the goodness of the fit.

A. Theoretical considerations
As it stands, our effective theory is still too general to calculate the production and detection processes because we need to determine the free functions m 2 (φ 0 ) = m 2 ( x) = m 2 (ρ), β i (ρ), and d i (ρ). There are two possibilities for fixing the spatial-dependence of the scalar's mass. The first is to parameterise our ignorance by assuming a functional form for its density and constraining the parameters associated with this choice. The second is to provide an explicit model and thereby to calculate the spatial-dependence explicitly. We opt for the second choice for three reasons. First, it is not guaranteed that an arbitrary fitting function will reproduce the dynamics of any fundamental theory, so it is not clear what information is lost by making such a choice. Second, chameleon models are well-constrained by laboratory and astrophysical probes, so choosing a well-studied model allows us to explore complementarity with these bounds, and to determine the feasibility of our scenario by identifying the existence of regions of parameter space where our model can simultaneously satisfy these bounds and successfully explain the XENON1T signal. Finally, chameleons have β i and d i constant, so we can exemplify our scenario using a simple minimal model. Chameleon theories are subject to a no-go theorem [190] that excludes the possibility of self-acceleration defined strictly as acceleration in the Jordan frame but not the Einstein frame in the complete absence of any cosmological constant. Dark energy scenarios where the acceleration is driven by the scalar potential, i.e. a quintessence-like explanation, are not excluded by this theorem, but require a tuning of an overall additive constant. Of course, such a tuning is also present in selfacceleration scenarios as this additive constant is arbitrarily set to zero. The generalized couplings (disformal and kinetic-conformal, XT ) that we consider here are additional potential caveats to the no-go theorem since they were not considered when deriving it [192]. Other caveats are discussed in Ref. [190].
As discussed in more detail in Appendix A, the densitydependence of the chameleon's mass arises because its dynamics are governed by an effective potential where ρ is the density of the matter species coupled to the chameleon, β m is the strength of such coupling, and V (φ) is the bare potential, which would govern the dynamics if the field were not coupled to matter. A common model for the bare potential is the power-law chameleon [71,77,78] leading to a density-dependent mass at the minimum φ min of the effective potential where Λ is an energy scale and n is a power-law index (V (φ) ∝ φ −n ). Note that both n > 0 (inverse powerlaw) and n < 0 (power-law) can lead to the chameleon behaviour provided n = −1, −2, or an odd negative integer. We also assume that β m φ/M Pl 1, in order for the excursion of the chameleon field not to exceed M Pl /β m . Note that the swampland conjectures (see Ref. [193]) lead to a lower bound on β m cV /ρ, where c is a constant of order unity [194].

B. Production in the Sun
Chameleons can be resonantly produced in a dense magnetised plasma when the chameleon mass matches the plasma frequency of the environment. This process, governed by the chameleon-photon coupling in Eq. (6), occurs in the Sun within a narrow shell whose location depends on the chameleon rest mass [54,195]. Chameleon production can also occur through non-resonant processes, which occur in all magnetised regions inside the Sun. Here we adopt the non-resonant production mechanism and consider a magnetic field profile B = B(r), where r is the radial coordinate. 4 For the solar model, we adopt the profiles described in Ref. [198], which has also been used to derive the formula for the Primakoff flux used by the XENON1T collaboration, see Ref. [199]. We also note that there is some disagreement in the field between different solar models, see e.g. Refs. [200][201][202][203][204][205].
The resulting differential flux per unit energy of solar chameleons on Earth, resulting from isotropic production in the Sun, is given by where d = 1 A.U. is the Earth-Sun distance, and R t ∼ 0.7R is the tachocline radial coordinate. The flux of chameleons produced in the Sun, dΦ/dω, is calculated in Appendix B, see Eq. (B13). In principle, one could also consider production from toroidal magnetic modes deeper within the Sun. However, we note that there are significant uncertainties associated to the strength and profile of these modes [206][207][208]. The equipartition value for the large-scale toroidal magnetic field due to shearing is ∼ 1 T [206][207][208]. As we later assume a strength B t = 30 T for the tachocline magnetic field, we expect the associated contribution to the chameleon flux relative to the contribution we considered to be suppressed by a factor (1/30) 2 ≈ O(10 −3 ), and to appear at higher energies than those of interest for the XENON1T excess. While the strength of these modes may be amplified locally by up to O(10 2 ) their equipartition values in the convection zone to form sunspots [206], the significant uncertainties at play prevent us from fully quantifying the impact of these processes on our results. We thus conservatively choose to neglect the effect of toroidal magnetic modes on chameleon production deeper within the Sun, noting that these could lead to subdominant features in the event rate at higher energies, but deferring a full study to follow-up work.
We have not considered other more "standard" production mechanisms, which in the parameter space spanned are subdominant. For example, chameleons could be produced through so-called ABC reactions (atomic recombination and deexcitations, bremsstrahlung, and Compton), in a similar fashion to axions [209][210][211]. However, by virtue of the chameleon mass being densitydependent, we can always find large regions of parameter space where these processes are kinematically disfavoured (in particular by adjusting the energy scale Λ). As discussed at the end of Appendix A, this can be achieved by requiring that m 2 eff (ρ core ) Pl the plasma frequency squared given by Eq. (B6), ρ core 150 g/cm 3 the Sun's core density, and T core 1.5 × 10 7 K the Sun's core temperature. Typically we expect that the mass of the chameleons scale like m(ρ) ρ α where α = 3/4 for chameleons with n = 1. This implies that the mass deep inside the core increases typically by two orders of magnitude compared to the mass in the tachocline. If this condition is satisfied, production of solar chameleons in the deeper regions of the Sun is kinematically forbidden, with the overall flux being dominated by Primakoff-like production in the tachocline.
From Eq. (9), we see that the previous condition translates into constraints on β m , Λ, and n. Fixing β m and n > 0, we find that this condition sets an upper limit on the allowed value of Λ. Later in our analysis we will consider β m 10 2 and n = 1, since we find that the event rate in XENON1T is mostly sensitive to β γ and M e , and only weakly sensitive to β m , n, Λ, and M γ . For this choice of β m and n, we find that production of solar chameleons in the Sun's core is kinematically forbidden as long as Λ 1 µeV. We remark again that fixing Λ to other values has no appreciable effect on the XENON1T event rate, and hence on our results.

C. Chameleon detection
To reach the XENON1T detector and leave detectable imprints, solar chameleons with energies ω O(keV) need to propagate through various dense media unscathed. In general, chameleons with incoming energy ω will be able to traverse a dense barrier of energy ρ provided ω m(ρ). Let us focus on the benchmark point we discussed above and which we consider throughout the paper, where β m 10 2 , Λ 1 µeV, and n = 1, so that m ∝ ρ 3/4 . We thus need to ensure that chameleons make it through the densest material along their path. The highest density involved in the problem is that of lead, which the XENON1T detector is partially made of, and whose density is ρ Pb ∼ 10 g/cm 3 . For the above parameters, we find m(ρ Pb ) ∼ 0.6 keV, meaning that chameleons with energies ω 0.6 keV are able to reach the XENON1T detector. This is sufficient to ensure that chameleons are able to fit the XENON1T excess, which occurs at energies higher than this cut-off. Moreover, as rocks and the tachocline have densities respectively one and two orders of magnitude lower than that of lead, this also ensures that chameleons are able to escape the tachocline and propagate through the rock which makes up Gran Sasso (being mostly made of limestone, with density ρ ∼ 3 g/cm 3 ).
In the XENON1T detector, solar chameleons can be absorbed by electrons via the chameleo-electric effect. This is the chameleon analogue of the photoelectric and axio-electric effects for photons and axions respectively. The cross-section for the chameleo-electric effect in DM direct detection experiments, σ φe , is computed in Appendix C, see in particular Eq. (C30). The resulting differential event rate per unit production energy ω for chameleon absorption by electrons in the XENON1T detector is given by: where Θ(ω) is the energy resolution of the detector and (ω) is the XENON1T detection efficiency, given in Fig. 2 of Ref. [80]. 5 The "raw" differential event rate per unit production energy of chameleons in the XENON1T detector, i.e. not taking into account energy resolution and detection efficiency effects, is given by: where the expression for the flux at Earth is given in Eq. (10), and where the number of atoms per ton of xenon is N Xe = 4.6 × 10 27 ton −1 , so that the differential event rate is expressed in units of ton −1 yr −1 keV −1 . We have appended the subscript th to the differential event rate per production energy in Eq. (11) to stress that this is a theoretical event rate, which depends on the underlying chameleon parameters through the dependence of dΦ Earth /dω and σ φe in Eq. (12) on these parameters. Comparing the theoretical event rate against the event rate measured by XENON1T will allow us to set constraints on the underlying chameleon parameters. While we focus on the XENON1T apparatus, we stress that the results in Appendix B and Appendix C, from which we derive Eqs. (11)- (12), are more broadly applicable. In particular, they can be applied to future DM direct detection experiments, which we discuss in more detail in Sec. VI B. Note, moreover, that because we used an effective action, the expression for the cross-section we compute, as well as its derivation, are general results that can be applied to any mass and coupling and for different detector setups. Our goals are now to explore whether solar chameleons are able to account for at least part of the observed low-recoil excess observed in XENON1T and, if so, to determine benchmark examples of solar chameleon parameters which provide an adequate fit to the XENON1T signal.
As discussed in Sec. II, we do not include the kineticconformal c i coupling (i.e. the term proportional to XT i ), since we find a posteriori that including this operator does not lead to a substantial improvement in fit to the XENON1T signal, which would instead be required to warrant its inclusion as a free parameter. If we were to include this term, a fit to XENON1T data would strongly prefer setting the associated coupling to zero, leading to the one-parameter extension not being preferred from a statistical point of view. This may not be the case for other DE models but it is for the specific case of chameleon DE.
The physical reason why this coupling worsens the fit to the XENON1T signal is that the associated crosssection does not scale fast enough with energy ω. This leads to the corresponding peak in the resulting event rate being below the ≈ 2 keV required to explain the XENON1T excess. On the other hand, the disformal coupling leads to additional powers of ω in the associated cross-section, thereby moving the peak to the correct position to explain the XENON1T excess.

IV. ANALYSIS
In principle, the parameter space describing production in the Sun and subsequent detection in the XENON1T detector of solar chameleons is sixdimensional, and spanned by the following parameters: the coupling to matter (in this case electrons) β e ≡ β m , the scale governing the disformal coupling to electrons (hereafter "electron disformal scale") M e = M/d 1/4 e , the coupling to photons β γ , the scale governing the disformal coupling to photons (hereafter "photon disformal scale") M γ = M/d 1/4 γ , and finally the energy scale Λ and powerlaw index n describing the chameleon self-interaction potential as given in Eq. (A12).
To simplify our analysis, we set Λ = 1 µeV and n = 1. As explained earlier, requiring Λ O(µeV) ensures that the chameleon's effective mass is sufficiently high in the core of the Sun so that production of chameleons through the usual Compton and bremsstrahlung mechanisms is kinematically suppressed. On the other hand, n = 1 corresponds to the best-studied chameleon model, with bounds typically only being reported for this specific choice [37]. These extensive studies have excluded a large region of parameter space [29], and it is thus of interest to explore whether XENON1T is able to probe part of the remaining parameter space of this model. In any case, we have explicitly verified that fixing Λ and n to other values affects chameleon production and the resulting event rate well below the %-level, and thus has no appreciable effects on our results.
These choices leave us with 4 parameters: β e , M e , β γ , and M γ . However, we anticipate that XENON1T will be mostly sensitive to β γ and M e , and very weakly sensitive to β e and M γ , for the following reasons. Firstly, we expect the disformal detection channel to dominate over the conformal one, as the former scales with a higher power of energy than the latter, see Eq. (C30). This feature moves the peak in the event rate towards higher energies, better fitting the XENON1T excess. Therefore, detection in XENON1T is mostly controlled by the electron disformal scale M e rather than the coupling β e , with the associated cross-section scaling as 1/M 8 e , see Eq. (C30). Second, if we require that M γ O(keV) so that Primakoff production in horizontal branch stars does not dominate over neutrino losses [212], production in the Sun will predominantly proceed through the conformal channel. 6 As a result, production will mostly be controlled by the photon coupling β γ rather than the photon disformal scale M γ . In particular, the associated production flux scales as β 2 γ , see Eq. (B13). Finally, the expected event rate in XENON1T, which is the only quantity we can directly compare to observations, depends only on the product of the production flux and the detection cross-section, as can be clearly seen in Eq. (12). This product scales as We therefore expect that the XENON1T measurements will predominantly constrain the following parameter combination, which we denote by β eff , and refer to as the effective coupling: We can view β eff as being the chameleon equivalent of the product g aγ g ae for the case of solar axions produced via the Primakoff effect in the Sun and detected via the axioelectric effect in the XENON1T detector. In particular, the expected event rate in the XENON1T detector is proportional to β 2 eff , as we derive in Eq. (C31) in the appropriate limit discussed in the Appendix.
We now proceed to analyze the XENON1T measurements to verify whether solar chameleons can fit these measurements, and whether our previous expectations are met. We perform a Bayesian statistical analysis to constrain the four chameleon parameters, which we collectively denote by θ ≡ {β γ , M γ , β e , M e }. The likelihood L(θ|d) to observe the data d given a certain set of model parameters θ is where the χ 2 function entering the likelihood is given by with the sum being performed over the energy bins ω i at which XENON1T measure their event rates. In Eq. (15), (dR/dω) th denotes the theoretical event rate given by Eq. (11), (dR/dω) meas denotes the rate as measured by XENON1T (black re-binned datapoints of Fig. 4 in Ref. [80]), and B 0 denotes the background model (red curve of Fig. 4 in Ref. [80]), which is itself a function of energy. The XENON1T background model is described in more detail in Sec. IIIB of Ref. [80] (see in particular their Table 1 and Fig. 3), and includes contributions from ten different components, ranging from 214 Pb to solar neutrinos. We refer the reader to Ref. [80] for more a detailed discussion of B 0 .
For simplicity, we only consider the first 16 bins, in the recoil energy range 1.5 keV ω R 16.5 keV. We do not consider bins beyond the 16 th for two reasons: 1. the theoretical solar chameleon event rate drops quickly beyond the third bin, partly due to the limited width of the differential chameleon flux (see Fig. 2 below), and to the effects of energy resolution and detector efficiency, entering in Eq. (11) through Θ(ω) and (ω); 2. the measured rate in the bins beyond the third is highly consistent with the background model B 0 , and therefore does not call for new physics explanations: the only exceptions are a few anomalous bins (the 17 th , the 20 th , the 24 th , the 26 th , and the 29 th bins respectively), which none of the proposed theoretical models (including the solar axions, neutrino magnetic moment, and tritium explanations invoked by the XENON1T collaboration) have been able to explain. 7 Of these 16 bins, the second and third deviate the most from the XENON1T background model B 0 , and therefore contribute the most to the excess. We impose flat priors on log 10 β γ ∈ [0; 11], log 10 (M γ /keV) ∈ [0; 25], log 10 β e ∈ [1; 2], and log 10 (M e /keV) ∈ [0; 25]. For β γ , the upper prior edge is motivated by the latest results from the Kinetic Weakly Interacting Slim Particles (KWISP) detector on the CAST axion search experiment at CERN, which finds β γ < 10 11 [58], whereas the lower prior edge is arbitrary (we have checked that extending it to lower values does not qualitatively affect our final results). Due to the weak sensitivity of the XENON1T measurements to β e , we have chosen a narrow prior for this parameter. We have explicitly verified that extending the prior further has no effects on our results. Similarly, fixing β e (for instance to β e = 10 2 ) would also have no effect on our results, see later discussion below Eq. (16).
We allow the photon and electron disformal scales to span the range between the keV scale and the Planck scale. It is worth noting that limits on the disformal scale M obtained from collider searches, indicating M O(100 GeV) as for instance in Refs. [213][214][215][216] (including works from one of us), only apply to the chameleon-quark disformal coupling scale, and not to the scale governing the coupling to photons and electrons. The strongest bound on the scalar-photon disformal coupling comes from demanding that Primakoff production of scalars in horizontal branch stars does not dominate over neutrino losses, and the strongest bound on the scalar-electron disformal coupling similarly derives from demanding that losses from Compton and bremsstrahlung production do not significantly alter the properties of these objects. In both cases, the bounds impose M e , M γ O(0.1 GeV) as derived by one of us in Ref. [212]. Note, however, that these bounds only apply in the limit where the scalar's mass can be neglected. This is not the case in our scenario since we impose m eff (ϕ ) > T , where T is the core temperature of the Sun. The bounds derived in Ref. [212] therefore do not directly apply to our case. Moreover, as we will discuss in Sec. VI A, we expect production within these stellar objects to be kinematically suppressed for the benchmark point in parameter space we considered. We consequently conservatively choose to allow M e and M γ to be as low as O(keV), but not any lower. As discussed at the end of Appendix B, for M e , M γ O(keV), the back-reaction effect of the disformal coupling on the scalar field profile in the Sun can become non-negligible, resulting in the break-down of the resulting production flux computation. The O(keV) scale which determines whether or not this effect is negligible is set by the maximal temperature reached within the Sun.
ground model B 0 . The solar axion explanation of the low-energy excess improves the fit to the 14 th bin through the contribution of 57 Fe axions, see Fig. 7b in Ref. [80].
To sample the posterior distribution of the chameleon parameters we use Markov Chain Monte Carlo (MCMC) methods. We make use of the cosmological sampler Montepython3.3 [217], configured to act as a generic sampler. The convergence of the generated MCMC chains is monitored through the Gelman-Rubin parameter R−1 [218], and we require R−1 < 0.01 for the chains to be considered converged.
Finally, we also quantify the significance of the preference (if any) for the solar chameleon model over the background-only model B 0 . We do so by adopting the same test statistic q(s) used by the XENON1T collaboration, with s symbolically denoting the signal parameters, see Eq. (17) in Ref. [80]. This is essentially a profile log-likelihood test statistic. The statistical significance of the preference (if any) for the solar chameleon model is then determined by q(0), i.e. comparing the difference in goodness-of-fit of the best-fit solar chameleon model relative to a B 0 -only fit.

V. RESULTS
We now analyze the XENON1T data using the methodology, priors, and likelihood described in Section IV. We perform an MCMC run on the fourdimensional parameter space spanned by β γ , β e , M γ , and M e . This MCMC run confirms our earlier expectation that we are only sensitive to the parameter combination of β γ and M e given by β eff in Eq. (13), while not being sensitive to β e and M γ . We shall discuss the obtained constraints on β eff later on.
A very important result of our analysis is that we are able to identify regions/benchmark points in parameter space which provide a good fit to the XENON1T signal (to be quantified shortly). One such benchmark example is presented in Fig. 1, where the blue curve is obtained by fixing the chameleon parameters to β e = 10 2 , M e = 10 3.6 keV, β γ = 10 10 , M γ = 1000 TeV, Λ = 1 µeV, and n = 1. The black data points denote the XENON1T measurements, and the grey curve is the XENON1T background model B 0 (the measurements and background are taken from Ref. [80]). Overall, the resulting fit to the XENON1T signal is good, with a bestfit χ 2 min 13.2 for 16 datapoints. Moreover, we find an improvement in fit of χ 2 min − χ 2 B0 −4.0 with respect to a B 0 -only fit, with the latter delivering χ 2 B0 17.2. Under this signal model, and using the previously discussed q(s) profile log-likelihood test statistic, a B 0 -only fit to the signal is rejected at a significance of 2.0σ.
The detection rate in Eq. (12) depends on the spectrum of chameleons on Earth resulting from Primakofflike production in the tachocline, as given by Eq. (10). In Fig. 2 we show this solar chameleon flux on Earth using the same benchmark choice of chameleon parameters as in Fig. 1. Note that the production flux remains high even at energies below 1 keV, an aspect which will have important consequences for our subsequent discussion. Although the fit in Fig. 1 is visually adequate, some features require further investigation. In fact, focusing on the second and third bins, i.e. the two bins where the measured event rate deviates the most from B 0 , the quality of the fit is only slightly worse that of the solar axion, neutrino magnetic moment, and tritium explanations invoked by XENON1T (see Figs. 7a-7d in Ref. [80]). For comparison, the predicted signal resulting from the bestfit solar axion model is given by the red curve in Fig. 1. However, as discussed in Ref. [80], the background model is rejected at a significance of more than 3σ within all these signal models, much stronger than our 2.0σ.
The paramount difference between the solar axion and solar chameleon models is the available production channels. The effects of this are evident by comparing the blue and red curves in Fig. 1 (for solar chameleons and solar axions respectively). As discussed in the Introduction, screening prevents the production of chameleons in the core through bremsstrahlung, Primakoff, and Compton effects and proceeds via the Primakoff effect [54] in the tachocline. Axions, on the other hand, are not affected by screening, and are produced in the core of the Sun  through different mechanisms: these include ABC reactions [209][210][211], the Primakoff effect [219], and the 57 Fe transition line [220]. These differences appear clearly in Fig. 1, where the solar chameleon model (blue line) shows one peak at ω ≈ 2 keV, while the solar axion model (red line) shows three distinct peaks at three different energies corresponding to the different production mechanisms. These features result in a better fit of the solar axion model to the third, fourth, fifth, sixth, and 14 th bins, improving the overall fit and significance at which the model is preferred over the background. In particular, ABC reactions are responsible for the considerably improved fit to the bins from the third to the sixth. These differences between the solar chameleon and solar axion scenarios can be distinguished in future DM direct detection experiments by their different spectra. We further discuss the prospects for detection of solar chameleons in future DM direct detection experiments in Section VI B. Although the measured recoil rate in the first bin is perfectly in line with the background model, the chameleon model overshoots this first point by just over one standard deviation. The reason is that the solar chameleon differential flux on Earth remains appreciable at low energies O(keV), see Fig. 2. This results in an increase in the event rate over B 0 at lower energies, leading to a poorer fit to the first bin, which is perfectly in line with B 0 and would not call for any additional contributions to the fit.
The reason why the signal beyond the third bin is completely dominated by the background model B 0 (which at that point is in very good agreement with the XENON1T measurements) is that the production flux, after peaking at an energy slightly below ω = 1 keV, quickly drops for higher energies (see Fig. 2). In other words, in the sum in the numerator of Eq. (15), one has (dR/dω) th + B 0 B 0 (dR/dω) meas for i ≥ 4, and therefore (dR/dω) th + B 0 − (dR/dω) meas 1, resulting in small contributions to the χ 2 for i ≥ 4 in the numerator of Eq. (15). For this reason, computing the χ 2 over all 16 bins might be a misleading goodness-of-fit metric. For the same choice of benchmark parameters as mentioned above, and in Fig. 1, the contribution of the first three bins to the total χ 2 is 6.5, which better quantifies the imperfect fit to the first three bins visible in Fig. 1.
Finally, let us discuss how XENON1T constrains the effective coupling β eff , given by Eq. (13). We treat β eff as a derived parameter whose posterior distribution we infer from our MCMC samples of the four fundamental parameters. The normalized posterior distribution for log 10 β eff is given in Fig. 3, which shows that XENON1T is indeed able to set meaningful constraints on this parameter.  (13,16). A value of log 10 β eff −4.5 is required to provide a good fit to the XENON1T signal, of quality identical to that shown in Fig. 1.
The shape of the log 10 β eff posterior, with a tail as log 10 β eff → 0 and a plateau for large negative values of log 10 β eff , is worth explaining further. Moving log 10 β eff → 0 means that either β γ is being raised or M e is being lowered. In other words, either or both solar chameleon production and detection are being enhanced. Enhancing the signal sufficiently will make the total event rate too large compared to the XENON1T measurements, and hence increasingly unlikely, leading to the tail in the log 10 β eff posterior as log 10 β eff → 0.
On the other hand, moving log 10 β eff towards large negative values means that either or both solar chameleon production and detection are being suppressed. Within this regime, the resulting event rate would be too low compared to the XENON1T detector background, and therefore the total signal is dominated by B 0 . This im-plies that (dR/dω) th +B 0 B 0 for all i in the numerator of Eq. (15), regardless of the choice of model parameters. This behaviour leads to an extended "plateau" in parameter space where the likelihood is completely flat with χ 2 17.2. Along the plateau, the quality of the resulting fit to the XENON1T signal is not only identical for any choice of parameters, but also identical to the quality of a B 0 -only fit. The goodness-of-fit along the plateau is worse by only ∆χ 2 = +4.0 with respect to the solar chameleon best-fit (for which χ 2 min 13.2, see Fig. 1). This explains the shallow plateau in the log 10 β eff posterior for large negative values of log 10 β eff .
We find that the best fit to the XENON1T signal occurs for a value of β eff 10 −4.5 . This implies Because of the shape of the log 10 β eff posterior shown in Fig. 3 (the tail can extend indefinitely to large negative values of log 10 β eff ), we do not quote a confidence interval on log 10 β eff . Rather, we use Eq. (16) as indicative of what combinations of β γ and M e lead to a good fit to the XENON1T signal, of quality identical to that shown in Fig. 1. We have also verified our earlier expectation that β e and M γ play a negligible role in our analysis. We have fixed these parameters to β e = 10 2 and M γ = 1000 TeV respectively, and repeated the analysis. Doing so, we find essentially the same posterior for log 10 β eff as shown in Fig. 3, which was instead previously obtained by varying all four parameters. Demanding that solar chameleons fit the XENON1T signal, and combining the relation in Eq. (16) and the upper limit of β γ < 10 11 from CAST [58] sets an upper limit on M e 10 MeV. In other words, if M e 10 MeV, for any allowed value of β γ the event rate will be completely dominated by B 0 , and we will find ourselves along the plateau for large negative values of log 10 β eff in Fig. 3. This upper limit on M e is ostensibly in contradiction with the lower limit one obtains by demanding that horizontal branch stars are not affected by the disformal coupling, which sets M e 0.1 GeV, as found by one of us [212]. However, we note that the limit obtained in Ref. [212] is not applicable to our case, as it was obtained assuming that the scalar is massless. This is clearly not the case in our scenario, given the constraints we have imposed on the chameleon mass within the Sun, which in turn suppresses production, making it so that the bounds obtained in Ref. [212] do not apply. This is generically true for all of the relevant limits in Ref. [212], which were all derived assuming a massless scalar.
The upper limit M e 10 MeV is nominally in tension with LEP constraints, which require M e 3 GeV as derived by one of us in Ref. [212]. However, care must be taken with this bound, as it again was derived assuming a massless chameleon. A proper re-evaluation of the bounds of Ref. [212] would require a dedicated analysis, for instance determining the field profile in the LEP pipe simultaneously including conformal and disformal cou-plings, a calculation which is well beyond the scope of this paper. Therefore, in continuing our exciting program for the direct detection of dark energy quanta, re-evaluating the LEP bounds is of paramount importance, and will be the subject of follow-up work. 8 Let us summarize the main findings of this Section: 1. The expected event rate in the XENON1T detector is sensitive to the combination of the photon coupling β γ and the electron disformal scale M e given by the effective coupling β eff in Eq. (13).
2. On the other hand, XENON1T has only very weak sensitivity to β e , M γ , Λ, and n.

3.
A value of log 10 β eff −4.5 is required to fit the XENON1T signal well. This leads to an improvement in χ 2 of 4.0 compared to a B 0 -only fit (B 0 excluded at 2.0σ), and a quality of fit as shown in Fig. 1. In no region of parameter space can the quality of the fit to the XENON1T signal be better than that shown in Fig. 1. As CAST requires β γ < 10 11 , demanding that solar chameleons explain the XENON1T signal and therefore log 10 β eff −4.5 implies M e 10 MeV. 4. Given the previous points 1. and 2., we can identify various combinations of the chameleon parameters which lead to a fit of identical quality to that shown in Fig. 1. There is therefore a large window of parameter space which can account for part of the XENON1T excess, while remaining consistent with laboratory and astrophysical tests.

5.
With respect to the solar axion model invoked by XENON1T, the statistical significance of the preference for the solar chameleon model is lower because of the poorer fits to the first bin, as well as to a few bins at higher energies (due to the larger available number of production channels for solar axions). However, we stress that the solar chameleon model is not excluded by other bounds, unlike the solar axion model.
Our overall conclusion is that solar chameleons are able to provide an adequate fit to the XENON1T signal. This raises the tantalizing possibility that XENON1T, originally constructed to detect dark matter, may have achieved the first direct detection of dark energy quanta. In principle, one may also consider a hybrid chameleonaxion scenario where both particles are present and contribute to the XENON1T signal, which thus results from an incoherent sum of solar-produced axions and chameleons. This could be beneficial for the solar axion model, as it might be able to alleviate the tension in the axion-electron coupling between the XENON1T results and astrophysical constraints [81]. Within this scenario, the lower-energy excess would be fitted by solar chameleons, plus a smaller contribution from solar axions, allowing for a lower value of g ae . On the other hand, the higher-energy end would be fit by the solar axion via its couplings to photons and nucleons. We defer the study of this interesting hybrid possibility to followup work.

VI. DISCUSSION
In this section, we discuss other experimental bounds on our scenario, stellar bounds in particular, and the prospects for detecting dark energy in current and planned dark matter direct detection experiments.

A. Stellar cooling constraints
The stellar bounds on the axion-electron and axionphoton coupling are debilitating for the solar axion interpretation of the XENON1T excess [81,222]. The situation with chameleons is different. The paramount difference between the two models is the environmentdependence of the chameleon's mass. This ensures that chameleons are not produced in the Sun's core since Compton and bremsstrahlung processes are kinematically suppressed, leading to a severe Boltzmann suppression. Instead, chameleons are produced in the strong magnetic field of the solar tachocline. Similarly, the cores of red giant, horizontal branch (HB), and white dwarf stars are denser than the Sun's by a significant factor, implying an even stronger suppression in these objects. It is possible that some of these objects may have strong magnetic fields [223], but without dedicated individual observations and detailed stellar modelling, it is not possible to derive quantitative constraints on chameleons.
Additionally, in Tab. I we report typical core densities and temperatures for these objects, alongside the chameleon mass for the benchmark parameter space point we have considered throughout the paper. As we see, production of chameleons within these objects is strongly kinematically suppressed, even more so than within the Sun, implying that stellar cooling constraints are evaded. Therefore, the bounds obtained by one of us in Ref. [212], derived assuming a massless chameleon, may be safely evaded.
Finally, we note that the literature is rich with astrophysical (stellar and galactic) bounds on chameleons, see e.g. Refs. [27,28,[224][225][226][227][228][229]. These refer to searches for the effects of fifth forces rather than chameleon particle production. Chameleon models predicting fifth forces in astrophysical objects occupy a different region of param-  I. Typical core densities and temperatures for stellar objects of interest: the Sun, white dwarfs, red giants, and horizontal branch stars. The final column reports the chameleon mass within the core of these objects for the benchmark parameter space point we have considered throughout the paper, where βe = 10 2 , Λ = 1 µeV, and n = 1. It is clear that production of stellar chameleons is strongly kinematically suppressed within these objects, as mcore Tcore.
eter space than those that give rise to chameleon particle production in the solar tachocline considered in this work. The underlying reason for this is that fifth forces are only relevant in astrophysical bodies of radius R if m(ϕ 0 )R ∼ 1, where ϕ 0 is the background field in that body. The chameleon theories accessible to the direct detection experiments that we have discussed in this work have m(ϕ )R 1, implying a fifth force range too small to affect stellar structure. For this reason, it is generally the case that astrophysical fifth force searches do not constrain our proposed scenario. These theories may be subject to the bounds arising from laboratory tests [29,37], although such bounds are highly modeldependent. The specific model studied in this work is able to simultaneously satisfy all experimental bounds and account for part of the XENON1T signal.

B. Other dark matter direct detection experiments
Recent DM direct detection searches prior to XENON1T did not report any excess over the expected background. For example, the PandaX-II experiment with an exposure of ≈ 27 ton × day placed an upper limit on the axion-electron coupling g ae 4.35 × 10 −12 for an axion mass m a 1 keV [230]. 9 A competitive limit has also been placed by the LUX-ZEPLIN collaboration with an exposure of 11.2 ton × day, which lead to the result g ae 3.5 × 10 −12 at 90% confidence level [232]. These exposures are all significantly lower than the XENON1T exposure of 0.65 ton × yr, which could explain why these experiments did not observe any low-energy excess.
Various upcoming experiments plan to search for the signal reported by the XENON1T collaboration, and will 9 A more recent analysis with an exposure of 100.7 ton×day found a similar result, gae 4.6 × 10 −12 [231]. and electron recoil background in units of (ton × yr × keV) −1 for recoil energies 10 keV, expected for the upcoming XENONnT [146], PandaX-4T [145], and LUX-ZEPLIN [141] experiments, which will be able to confirm or disprove the XENON1T excess. The last column reports the number of excess events that are expected per year in each detector, in the energy range (1 − 30) keV. either confirm or disprove it. These experiments include XENONnT (the planned upgrade to XENON1T) [146], as well as PandaX-4T [145] and LUX-ZEPLIN [141], all of which use a dual-phase xenon time projection chamber. In Table II we report the expected exposure in units of ton × yr and electron recoil background in units of (ton × yr × keV) −1 for each of these experiments.
We focus on the benchmark point in parameter space considered in Sec. V, as well as Figs. 1 and 2. Adopting these values, the expected excess number of events per year due to a hypothetical signal is about 20 for XENON1T, 180 for XENONnT, 130 for PandaX-4T, and 250 for LUX-ZEPLIN. Note that we have not considered the effects of background noise or energy resolution in obtaining these estimates of excess number of events since the energy resolution specifications for these future experiments are currently unavailable.
As is clear from Table II, all of these next-generation xenon-based experiments will have a background level B 0 a factor of ≈ 5-6 lower than current levels, while the exposure will increase by over an order of magnitude. This combination of lower background and increased exposure results leads to the extremely high number of expected events. By virtue of this, future experiments will be able to confirm or disprove our hypothesis that solar chameleons are the origin of the XENON1T signal with extremely high statistical significance. More generally, it will be possible to test at high significance whether the XENON1T excess is due to a statistical fluke, a background contaminant, or new physics such as the scenario considered here.

VII. CONCLUSIONS
Most of our knowledge about dark energy (DE) arises from cosmological measurements which are mainly sensitive to its gravitational effects. Yet searching for non-gravitational signatures of DE by directly detecting DE quanta would be an extremely important step towards understanding the physics powering cosmic acceleration. In this paper, our aim has been to broaden the scope of new physics accessible to terrestrial dark matter (DM) direct detection experiments by investigating the intriguing possibility that these instruments may be able to detect DE quanta via their couplings to matter. Specifically, we have envisaged a scenario wherein DE particles produced in the strong magnetic fields of the solar tachocline travel to Earth and are absorbed by electrons or nuclei in terrestrial DM detectors. In this paper, focusing on DE scalars including screening mechanisms of the chameleon type, we have laid out the formalism for computing the expected signal from such a process, demonstrating that it can lead to measurable recoils in the O(keV) range, well within the sensitivity of current and upcoming DM direct detection experiments.
We have applied our results to the XENON1T experiment, which recently reported a ≈ 3.3σ excess in their electron recoil data at recoil energies of ≈ 1-2 keV. We have shown that solar chameleons can explain the XENON1T excess (see blue curve in Fig. 1), and are preferred over the background-only hypothesis at a significance of ≈ 2.0σ. Our results have been obtained using the code which we make publicly available at github.com/lucavisinelli/XENONCHAM.
Compared to the much discussed solar axion interpretation of the XENON1T signal, the statistical preference for solar chameleons is lower, mostly due to the reduced number of available production channels in the Sun (chameleons are only produced through Primakofflike processes in the tachocline). However, the stellar cooling constraints which are debilitating for the solar axion model do not apply to solar chameleons due to the environment-dependence of the chameleon mass, which within the dense environments of red giants and white dwarfs results in Compton and bremsstrahlung production processes being suppressed.
We have also studied prospects for testing this explanation in future DM direct detection experiments. If solar chameleons are indeed at the origin of the XENON1T excess then this will be confirmed at very high significance in upcoming experiments such as XENONnT, PandaX-4T, and LUX-ZEPLIN. What is perhaps more important is that future low-threshold DM direct detection experiments will be well-suited to detect the signatures of DE particles produced within the Sun.
There are several avenues for future research in this direction. While our study has focused on chameleonscreened scalars, we stress that the effective theory approach we have adopted is quite general, and our formalism can therefore be applied (with appropriate modifications) to several other physical scenarios of interest. More generally, we hope that this paper will stimulate further research aimed towards enabling direct detection of dark energy, searching for non-gravitational signatures of dark energy, and unraveling the physics of cosmic acceleration in terrestrial laboratories. Chameleon-screened theories provide an explicit example where the mapping between the effective theory in the solar system given in the main text and a more complete theory can be calculated [189] (see Refs. [29,233] for more details about screening). In this case, the scalar φ is canonically normalised with a scalar potential V (φ). The interaction with matter is obtained via the Jordan frame metric g J µν = A 2 (φ)g µν . We treat the disformal and derivative interactions as negligible perturbations compared to the dominant effect due to the conformal rescaling given by A(φ). This is valid provided the suppression scales of the disformal and derivative interactions are large enough. In the Einstein frame, where the graviton and scalar are canonically normalized but the coupling to matter is non-minimal, the dynamics of φ depend on the effective potential where ρ is the conserved matter density in the Einstein frame, related to the density in the Einstein frame as ρ E = Aρ. This relationship between ρ E and ρ follows from the non-conservation of the Einstein-frame energy momentum tensor which, for pressureless matter, gives [234]ρ which expresses the local conservation of the matter density ρ. The Jordan frame matter density ρ J = A −4 ρ E is conserved in the Jordan frame where dt J = Adt E and the local Hubble rate is given by The Jordan frame matter density is deemed to be the "physical" density as, in the local Jordan frame where the Jordan metric is nearly Minkowskian, the Lagrangian of the standard model reduces to the usual one. As long as A ≈ 1, the difference between ρ and ρ J is negligible. Chameleon-screened models have a minimum of the effective potential φ(ρ) which depends on the conserved matter density ρ in the Einstein frame. The minimum equation reads which can be used to obtained a parametric description of the value φ(ρ). Taking the derivative of the minimum equation with respect to ρ leads to where we have defined the effective mass It is convenient to introduce the effective coupling and to integrate (A6) This provides a one-to-one relationship between the density of matter and the value of the minimum φ(ρ). Moreover from (A6) we get allowing one to obtain A(ρ) as a decreasing function of ρ, i.e. A −1 is an increasing function of ρ with a positive derivative. Finally the minimum equation gives from which we can find V (ρ). Hence eliminating ρ between A(ρ), V (ρ), and φ(ρ), one can reconstruct A(φ) and V (φ). The potential V (φ) is defined up to an integration constant, i.e. the screening properties of the models do not depend on an additive cosmological constant. This constant has to be tuned to generate the appropriate acceleration of the Universe. As the screening properties are independent of this choice, this has no effects on the results obtained in this paper. As a specific example, let us consider the inverse power law potential.
where we have included the constant V 0 which needs to be adjusted to fit the current dark energy value. As we have seen, the screening properties only depend on the the inverse power law part. We note that models could arise from the strong dynamics of a confining supersymmetric dark sector at higher energy [235]. The configuration of the field at which the potential is minimized φ(ρ), and the chameleon rest mass squared obtained from the curvature of the effective potential, m 2 φ , are given by This is the mass m φ that we use as a template in the main text [see Eq. (9)].
Finally this reconstruction procedure of V (φ) and A(φ) from m φ (ρ) and β(ρ) allows one to design models where the production of scalars in the tachocline is favoured compared to very deep inside the Sun or in other even denser stars. In the Sun, all that is required is that in the core where the density ρ core ≈ 150 g cm −3 is large compared to the one in the tachocline ρ tach ≈ 1 g cm −3 , the production of scalars is kinematically forbidden. For chameleons produced in matter, e.g. by the Primakoff process deep in the electric field of a nucleus, their effective mass is modified by the presence of the surrounding plasma and their mixing with photons as m 2 eff = m 2 φ (ρ core ) − ω 2 Pl where the plasma frequency is defined in Eq. (B6). For large enough densities, the mass m φ (ρ core ) is generically larger than the plasma frequency which scales as ρ 1/2 . 10 As a result, only scalars of momenta k 2 T 2 core − m 2 eff can be produced. When m eff T core , production is highly suppressed. This also applies to the Compton and bremsstrahlung processes involving the direct coupling between scalars and electrons where the same kinematical obstruction is at play. This mechanism was first proposed in Ref. [236].

Appendix B: Production of chameleons in the Sun
Given the chameleon-photon disformal coupling in Eq. (6), the probability of a photon in the uniform magnetic field B converting into a chameleon after a distance is given by [195] where the coefficients are In the expressions above, ω is the energy of the produced chameleon, the dimensionless parameter b = B 2 t /M 4 γ is the ratio of the magnetic field in the solar tachocline B t to the UV-cutoff scale of the effective theory, and the plasma frequency is given in terms of the electron number density n e and the electron mass m e as The quantity B z is the z-component of the magnetic field which we fix by assuming an isotropic magnetic field distribution as B 2 z = B 2 /3. The thickness of the tachocline is much larger than the main free path of photons in the region, λ ≈ 0.3 cm [237,238], so that photon propagation in this region proceeds through a random walk process, which can be described as a Poisson diffusion process with mean free path λ. For a typical distance between two scatterings, the total number of scatterings per unit time is ∼ c/ . For a given length path , the differential probability of conversion in the solar interior is [54] where ls = ct 3 × 10 10 cm in the tachocline, i.e. approximately one light-second. Here,t is the typical time such that the photon flux at the tachocline n γ,t =vn t wheren t is the photon number density at the tachocline andv = (cλ/t) 1/2 is the typical radial velocity of photons due to their Brownian motion. The differential flux of chameleons per unit energy emitted by the Sun is where n γ (R) is the photon flux profile and the photon spectrum p γ (R) depends on the temperature profile of the plasma T = T (R) as We model the magnetic field profile inside the Sun as a thin shell around the solar tachocline, where the magnetic field is taken to be constant with a value B t = 30 T. The thin shell around the tachocline has radius R t = 0.7R , where R is the radius of the Sun, and the thickness ∆R = 0.01R . The integrand in Eq. (B8) at the tachocline is then where n γ,t = n γ (R t ) ≈ 10 21 cm −2 s −1 and p γ,t is the expression in Eq. (B9) evaluated at the tachocline temperature T ≈ 0.2 keV, and we have used Eq. (B7) which expresses the differential probability of conversion. Inserting Eq. (B1) into Eq. (B10) we obtain (B12) For most of the region of the parameter space we explore, the relation ω λ holds, for which the integral in Eq. (B12) is I( ω /λ) ≈ √ π. For the region of parameters allowed we have ∆ B ∆ pl . We also assume M γ O(keV), which corresponds to b 1. In this limit, the expression for the solar chameleon flux in Eq. (B11) reduces to dΦ dω = p γ,t n γ,t ∆R λ (B13) In the limit where m 2 φ ω 2 pl , Eq. (B13) reduces to dΦ dω = p γ,t n γ,t ∆R λ The computation we have just outlined is valid as long as the effect of the disformal coupling on the scalar field profile in the Sun is negligible. In other words, that the back-reaction effect due to the disformal coupling can be neglected. The disformal coupling leads to a contribution to the kinetic term of the scalar field proportional to P/M 4 i [239][240][241][242][243], with P ∝ T 4 the pressure of the solar photon gas, at a temperature T (where, in our units, the proportionality factor is smaller than unity). Therefore, requiring that the disformal coupling does not modify the field-profile in the Sun is tantamount to requiring that M i T core , and therefore M i O(keV). The condition in the tachocline is weaker as its temperature is lower than the core temperature by an order of magnitude.
In the analysis of Sec. IV, we have imposed priors which ensure that M γ , M e O(keV), to satisfy the previous bound. In addition, we note that for the benchmark point of parameter space we have discussed throughout the paper, and in particular in Figs. 1 and 2, we have set M e = 10 3.6 keV and M γ = 1000 TeV, such that the bound is well satisfied.

Appendix C: Chameleo-electric cross-section
In this Appendix we derive an expression for the detection cross-section of chameleons from the analogue of the photoelectric effect. The cross-section for the "chameleoelectric" effect receives contributions from each of the three terms in the last line of the Lagrangian in Eq. (5). We first consider the disformal coupling between φ and the electron, where M e = M/d 1/4 e is the energy scale related to the disformal coupling with electrons. The stress-energy tensor associated with an electron four-spinor ψ is where γ µ are the Dirac matrices and a bracket denotes a symmetrisation over the four-indices µ, ν. In the following, we adopt Feynman slash notation A = γ µ A µ for a four-vector A µ , and we defineψ = ψ † γ 0 .
We decompose the electron free field as ψ = s d 3 p 2E p u s (p)b s e −ip·x + v s (p)c s e ip·x , (C3) where the subscript s labels the spinor component and u s (p) and v s (p) are Dirac spinors following the normalisation condition s u sūs = p+m e and s v svs = p−m e . The operator b s (p) and its adjoint satisfy the anticommutation relation and similarly for the operator c s (p), where curly brackets denote the anti-commutation of the two operators. We have defined d 3 p = d 3 p/(2π) 3 . The non-relativistic electron bound state is where the Dirac spinor χ s (p) has the antiparticle entries equal to zero and ϕ(p) is a non-relativistic wave function in the momentum representation. We consider the ground state of a bound electron, where Z is the atomic number (Z = 131 for xenon) and a 0 is the Bohr radius. In momentum space, we obtain ϕ(p) = d 3 re −ip·r ϕ(r) = 8 √ π (p 2 + (Z/a 0 ) 2 ) The scattering vertex from the disformal coupling in Eq. (C1) is sketched in the Feynman diagram in Fig. 4 and amplitude given by where we have introduced the vector A µ = k µ (p + p ) ν k ν + k µ (p + p ) ν k ν .
The square of the amplitude summed over the spins of the final states and averaged over the spins of the initial states is with the cross-section Here, δ(x) = 2πδ(x) and we decomposed the four-vectors as follows. We consider the incoming and outgoing chameleon four-vectors k µ = (ω, k) and (k ) µ = (ω , k ). The electron is described by the systems p µ = (E, p) and (p ) µ = (E , p ), where E = m e − E b and E b is the binding energy of the atomic electron. In the following, we make use of the magnitudes k = |k|, k = |k |, p = |p|, and p = |p |. We also introduce the angle θ between k and p, the angle θ between k and k , and the azimuthal angle τ between the projections of k and p on the plane orthogonal to k. The orientations of these vectors are sketched in Fig. 5.
To proceed with the computation, we consider the conservation of the four-vector on shell E + k 2 + m 2 φ = (p ) 2 +m 2 e + (k ) 2 +m 2 φ ,(C14) (p ) 2 −(k ) 2 −p 2 = k 2 −2pk x+2pky−2kk y , (C15) where y = cos θ, y = cos θ , x = yy + sin θ sin θ cos τ . In the detector, the mass of the chameleon is expected to be set by a resonance condition involving the size of the cavity R [71,244,245], and thus to be of the order of m φ ∼ 1/R ≈ 10 −7 eV, which is much smaller than other energies in the system. For this reason, we neglect m φ in the rest of the computation. Combining Eqs. (C14) and (C15) we obtain where "≈" indicates the limit m e k |E b | m φ . We define the product α 1 = (p + p ) ν k ν = (E + E )ω − (p + p ) · k = = (E +E)ω−(2py+k−k y ) k , where in the last line the spatial part is p = p + k − k . Similarly, we define We then have A µ = α 1 k µ +α 2 k µ . In the limit considered, E ≈ E ≈ m e , so we obtain In the definitions of α 1 and α 2 , the temporal part of the four-product is the dominant one in the limit considered. The time component of the vector A µ is where in the last step we used Eq. (C19). The square of the vector A µ is A 2 = A µ A µ = (α 1 k µ + α 2 k µ ) α 1 k µ + α 2 k µ = = m 2 φ α 2 1 + α 2 2 + 2α 1 α 2 (ωω − kk y ) ≈ ≈ 2α 1 α 2 (ωω − kk y ) , where the last approximation assumes a massless chameleon at detection, k µ k µ = 0 or |k| = ω. Since ω ≈ ω, in the limit considered we have We evaluate the product where in the last step E ≈ E ≈ m e |p|, |k|, |k | . Using the approximation where the angular integral is trivial as there are no angles appearing in Eq. (C26). In the last step, we have normalised the wave function according to Eq. (C8). A second contribution to the cross section comes from the conformal term in Eq. (5) for which the absorption cross section depends on the photo-electric cross section σ photo in the limit ω m φ as [210,246,247] σ φe,conf = β 2 e ω 2 2παM 2 Pl σ photo . (C29) We have taken the energy-dependent photoelectric cross section from Ref. [248]. Note, that we have not considered the production/detection from the XT coupling and therefore set c e = 0. In terms of the parameters used in the MCMC analysis, Eqs. (C27) and (C29) combine to give the following expression for the cross-section σ φe = σ φe,dis + σ φe,conf = m 2 e ω 4 8π 2 M 8 e + β 2 e ω 2 2παM 2 Pl σ photo .
(C30) Although in this work we have focused on xenon-based detectors such as XENON1T, the results of the computation can be applied more broadly to any material.