Stellar Basins of Gravitationally Bound Particles

A new physical phenomenon is identified: volumetric stellar emission into gravitationally bound orbits of weakly coupled particles such as axions, moduli, hidden photons, and neutrinos. While only a tiny fraction of the instantaneous luminosity of a star (the vast majority of the emission is into relativistic modes), the continual injection of these particles into a small part of phase space causes them to accumulate over astrophysically long time scales, forming what I call a"stellar basin", in analogy with the geologic kind. The energy density of the Solar basin will surpass that of the relativistic Solar flux at Earth's location after only a million years, for any sufficiently long-lived particle produced through an emission process whose matrix elements are unsuppressed at low momentum. This observation has immediate and striking consequences for direct detection experiments---including new limits on axion parameter space independent of dark matter assumptions---and may also increase the prospects for indirect detection of weakly interacting particles around compact stars.


INTRODUCTION
Stars are poor photon emitters. Their photon opacity is so high that the only effective radiating component is a thin shell near the stellar surface. Stellar energy losses per unit volume are thus suppressed by the surface-tovolume ratio. They are further diminished by thermal self-insulation: photon luminosity scales as the fourth power of the surface temperature, which is typically several orders of magnitude lower than the core temperature (a factor of about 2000 in the Sun).
Stars can thus serve as sensitive "astrophysical laboratories" of weakly coupled particles [1][2][3], whose contributions to the overall luminosity-via volumetric emission-and thermal transport-via long mean free paths-can be disproportionately large. Indeed, neutrino emission is the main energy loss mechanism for the first 10 5 y after the birth of a neutron star [4], despite neutrinos' tiny coupling. Even the Sun has a fractional luminosity of order 10 −9 due to thermal neutrino-pair production [5], in addition to the 3% from fusion neutrinos.
A vast number of experimental efforts (notable examples include [43][44][45][46][47][48][49][50][51][52]) are ongoing to detect (in a terrestrial laboratory setting) the relativistic Solar flux and a potential Galactic DM population of these particles. One leading experiment, sensitive to both types, recently reported [53] a statistically significant excess of electronrecoil events with energies just above a keV, intriguingly near the Sun's core temperature T . If the excess were in fact due to new fundamental physics beyond the Standard Model, many simple models appear in conflict with the aforementioned stellar cooling constraints, discounting these exotic hypotheses in favor of more mundane explanations such as unforeseen radioactive backgrounds.
In this work, I investigate a hitherto unknown effect in astroparticle physics that seems benign at first, but gives rise to dramatic observable consequences in laboratory experiments-including the limits and putative excess of Ref. [53]-and astrophysical observations. I stress that no new particle physics model is introduced; rather, I posit the occurrence of a new phenomenon that is generic within a large class of well-motivated models. The main point of this paper is simple: stars are able to emit massive particles into bound orbits. The energy loss rate for this bound emission component is typically very small (see Eq. 5). However, this peculiar ensemble of particles populates parts of phase space that may survive for millions to billions of years around isolated stars, even in the inner Solar System. Over time, the density of this "stellar basin" can exceed the density of the relativistic stellar flux, including within the Solar System (cfr. Eq. 8).
In what follows, I describe how stellar basins of lowmass, weakly interacting particles can form and evolve. I discuss general aspects of soft emission near a particle mass threshold, before delving into a case study of Solar axions coupled to electrons. Based on these results, I present new limits on axion parameter space, and suggest the possible re-interpretation of the excess of Ref. [53] as due to a Solar axion basin. I conclude with potential avenues for future work. A generic spectrum dQ/dω, the differential energy loss rate per unit volume per energy ω of the emitted particle, shown in black. If the particle has mass m, then stellar emission just above threshold, specifically m < ω < m(1 + v 2 esc /2) with vesc the escape velocity, consists of nonrelativistic modes gravitationally bound to the star (blue region). While only a small fraction of the per-volume energy loss rate Q, given by the integral under the curve and most of which escapes to infinity (yellow region), the bound-orbit emission can accumulate over long times to form a "stellar basin". The vertical axis has arbitrary units, and the escape velocity is greatly exaggerated for illustrative purposes. A massless spectrum is shown in dashed gray. The energy is in units of temperature T , with characteristic Boltzmann suppression at ω T .

Basin Formation
Consider a volume-emission process of a weakly interacting particle with mass m, with dQ/dω the differential energy loss rate per unit volume per energy ω = m 2 + k 2 of the emitted particle. The next two sections will explain how to calculate dQ/dω in general and in a specific model, respectively; for now, assume it exists with some spectrum like the one shown in Fig. 1 for example. The total luminosity for this process is given by a volume integral over the stellar interior For any physical process of interest, the vast majority of the luminosity is radiated away to infinity, leading to an energy density in unbound particles that falls off as the inverse square radius R outside the stellar surface at R * . However, a small fraction of the luminosity is emitted just above threshold (ω m, |k| m), into bound orbits. Expanding in small kinetic energy per unit massω k ≡ (ω − m)/m k 2 /2m 2 , one finds in general that where n p is a positive integer, the sum is over different emission processes p, and energy dependence is extracted out such thatQ p only depends on the radius R within the stellar interior (assumed to be spherically symmetric). Particles emitted into bound (radial) orbits are those with negative asymptotic energy per unit massẼ ≡ ω k +Φ(R) where Φ(R) is the gravitational potential. The probability density P(R,Ẽ) for a particle of normalized energyẼ to be at radius R is proportional to the inverse of its velocity v(R,Ẽ) = 2[Ẽ − Φ(R)]: The normalization constant C(Ẽ) is fixed by requiring d 3 R P(R,Ẽ) = 1; the quoted expression is valid for a gravitational potential Φ(R) −G N M * /R, and is a good approximation for orbits with maximum radius far outside the radius R * of the gravitating body.
The bound energy density at radius R ≥ R * grows at a rate: For R R * = max{R }, one can approximateẼ − Φ(R ) −Φ(R ), and C(Ẽ) as in Eq. 3, yielding the simplification: In most cases of interest, the integral above evaluates to a result of order p L p min{m 3 /T (0) 3 , 1}|Φ(0)| np/2 with T (0) and Φ(0) the core temperature and gravitational potential of the star, and L p the luminosity for each process (L ≡ p L p ). The crux of stellar basins is that the weakly coupled bound particles can accumulate for astrophysically long times, compensating for their lower production rate relative to the unbound flux from Eq. 1. Denote by τ the effective 1/e lifetime of a bound orbit at radius R. This lifetime is set by the most efficient of the three dominant processes for depleting the basin, namely (gravitational) ejection with lifetime τ eject , (effective) absorption with lifetime τ abs , and radiative decay with lifetime τ rad . More precisely, one has that τ −1 = τ −1 eject + τ −1 abs + τ −1 rad . The ejection lifetime in the Solar System at R = 1 AU is at least τ eject (AU) 10 7 y for example (see next section), often much faster than decay or absorption channels, as discussed later around Eqs. 16 and 18. After some time > τ , the basin will be filled and reach a (quasi-)steady-state density ρ b =ρ b τ . At this saturation point, the ratio of bound-to-unbound energy densities is: The second line contains the parametric estimate for the integral mentioned below Eq. 5, took only one dominant emission process with exponent n p , and is written in terms of the escape velocity v esc (R) ≡ −2Φ(R).
The salient feature of stellar basins-the large accumulation time τ -is responsible for the huge enhancement in the first factor in Eq. 7, which can easily outweigh the remaining small factors due to phase-space suppression. Near the surface of an isolated neutron star (R * ∼ 10 km), this factor can in principle be as large as its age over its light-crossing time, τ /R ∼ 10 21 at τ ∼ 10 9 y. Compact remnants like neutron stars or white dwarfs furthermore have large escape velocities, mitigating the phase-space suppression.
Even at Earth's location in the Solar System, the bound energy density in the basin may exceed that of the unbound Solar flux, for processes that are not suppressed near threshold (i.e. n p = 1) and for masses m not far below the core temperature T (0). Numerically, the escape velocities are v esc (AU) ≈ 1.40 × 10 −4 and v esc (0) ≈ 4.61 × 10 −3 , while τ /AU ≈ 6.3 × 10 11 (τ /10 7 y). Eq. 7 therefore indicates that the bound population can start dominating the unbound Solar flux of particles with mass m ∼ T (0) ∼ keV after only a million years: up to O(1) prefactors. Finally, the relative enhancement of the bound state population at Earth's surface from particle production in the core of Earth itself is even more striking. With surface and core escape velocities In absolute terms, ρ b in Earth's basin will turn out to be small, though it can equal or exceed that of the Solar basin at small masses m T ⊕ (0), depending on the ratio of gravitational ejection timescales.

Gravitational Ejection
A critical parameter determining the maximum energy density to which the stellar basin can fill up is the lifetime τ (R). Contributions to the inverse lifetime from the radiative decay rate τ −1 rad (Eq. 18) and from the effective absorption rate τ −1 abs (Eq. 16) in the stellar interior are relatively straightforward to calculate, given a specific particle physics model. However, in an important case of interest, namely that of ultra-weakly coupled particles in the Solar System, the gravitational ejection timescale τ eject (R) is the principal limiting factor.
The main obstacle towards an accurate determination of τ eject (R) is due to the well-known chaotic nature of orbits in the Solar System, with a Lyapunov time of approximately 5 million years at R = 1 AU [54]. Taking this to be the timescale at which the entire phase space volume is explored efficiently, a conservative estimate would be that orbits migrate to an ejection-inducing orbital resonance with a characteristic timescale of two Lyapunov times: Even with such an overly efficient ejection rate, the parametric estimate of Eq. 8 predicts a local fractional overdensity of the stellar basin over the relativistic flux density by up to one order of magnitude. At first glance, the estimate of Eq. 9 appears to be confirmed by numerical simulations of forward evolution of Near-Earth Objects (NEOs) [55,56], with Ref. [57] finding a median lifetime of 10 million years for 117 NEOs with perihelia q < 1.3 AU and high-quality initial orbital elements. However, the vast majority of NEOs ended their lives due to Sun-grazing orbits (56%) or planetary impacts (18%), both of which leave weakly coupled particles unharmed; 16% of NEOs survived the total integration time of 60 million years, while only 10% were ejected from the Solar System [57]. A fiducial 1/e ejection timescale would be based on the NEOs unaffected by impacts, a restriction that could admittedly introduce bias.
There are four reasons that suggest even the fiducial value of Eq. 10 is an underestimate. First, the stellar basin is populated with vastly different orbits than those of asteroids or captured comet remnants. The main contribution toρ b (AU) comes from directionally isotropic, radial orbits with aphelia very near 1 AU that are later processed-by secular perturbations from the planetsinto high-eccentricity e ≈ 1 elliptical orbits with semimajor axes a ≈ 0.5 AU at all inclinations i. In contrast, NEOs are preferentially located near the ecliptic plane i ≈ 0, and near unstable parts of phase space, since they were resonantly driven away from the asteroid belt or (to a lesser extent) captured by close encounters in the inner Solar System [58]. Second, the known NEO population is known to exhibit strong observational biases against (separately) high-e, high-i, and low-a orbits, all of which are associated with higher dynamical lifetimes [57][58][59]. Third, the mortality rate of NEOs does not follow a simple exponential decay [57], indicating that various parts of phase space may be protected (e.g. be locked into metastable Lidov-Kozai resonances [60]) and have longer ejection timescales. Fourth, since the basin energy density injection rate scales asρ b ∝ 1/R 4 , the total energy density in logarithmic radial shells (e.g. between R and 2R) scales as 1/R. So even if there were a depletion of ρ b (R), orbits interior to R may be responsible for compensating replenishment as they migrate outward.
Indeed, Refs. [61,62] studied analytically the dynamics of diffusion throughout all Earth-orbit-crossing phase space (but without account of orbital resonances), in the context of DM capture into the Solar System (see also Refs. [63,64]). The dominant Solar basin injection for that part of phase space is those orbits which reach aphelion at or just above 1 AU, and which therefore have nearly vanishing speed in the heliocentric frame when they cross Earth's orbital radius, or a velocity of v ⊕ ∼ 30 km/s (opposite Earth's motion) in the geocentric frame. Ref. [62] found the combined action of Venus and Earth leading to a time τ diff ∼ (M /M ⊕ ) 2 P ⊕ (with P ⊕ ≡ 1 y) to diffuse out of this region, i.e. longer than the age of the Solar System, leaving this part of phase space "unfilled" from Galactic DM capture (inner solid contour of Fig. 3 in Ref. [62]). By the timereversed argument, it would also not be depleted over that time. As these orbits do not cross Jupiter's, energy pumping due to the Jovian disturbing potential also does not occur before third order in secular perturbation theory away from orbital resonances [65], so at least as slow as τ pump ∼ P 4 J /P 3 ⊕ 3 × 10 9 y, where P J ≡ (a 3 J /G N M J ) 1/2 with M J and a J the Jovian mass and semi-major axis. Preliminary simulations [66] confirm that orbits in Ref. [62]'s "unfilled" region typically do remain there for 4.6 × 10 9 y, even given secular and resonant perturbations.
I therefore reckon that the actual gravitational ejection time is somewhere between the fiducial value of Eq. 10 based on forward simulations of (biased) NEO orbits and the maximally optimistic estimate τ eject (AU) ∼ 4.5 × 10 9 y (optimistic).
Nevertheless, this work prudently uses Eq. 9 for the purposes of recasting limits. The widely varying benchmarks of Eqs. 9, 10, and 11 are clearly not satisfactory, and analytic arguments are unlikely to settle this issue. Dedicated simulations are sorely needed to establish not just the typical ejection time, but also to characterize the statistical behavior of ρ b as well as its temporal intermittency and modulation due to Earth's weakly eccentric orbit.

SOFT EMISSION NEAR MASS THRESHOLD
For a general emission process of a weakly coupled boson with four-momentum k = (ω, k), the net energy loss rate per unit volume is: where |M| 2 is the spin-averaged square matrix element, P in ≡ (E in , P in ) = in i p i is the total four-momentum of incoming particles, dP in ≡ in i g i d 3 p i /(2E i ) are the usual momentum integral measures (similarly, P out and dP out for outgoing particles), with factors g i accounting for internal degrees of freedom such as spins. Shorthand d ≡ d/(2π) andδ() ≡ 2πδ() is used. Generalization to multi-particle emission, e.g. pair production of neutrinos, is straightforward. For some models of weakly coupled particles (such as kinetically mixed photons), special care must be taken to include medium-dependent dispersion relations and interactions, complications ignored here.
The phase-space distribution functions-f (k) of the weakly coupled particle and f i (p i ) of the in-and outgoing particles-are collected into the function: The three terms correspond to spontaneous emission, stimulated emission, and absorption, respectively.
with η i = ±1 if particle i is a boson/fermion to include Bose enhancement or Pauli blocking. All distribution functions are normalized such that n i = g i d 3 p i f i (p i ) gives the local number density. Those of the in-and outgoing particles are assumed to be thermally equilibrated at the same temperature T : with (relativistic) chemical potentials µ i . At early times, when the stellar basin is nearly empty f 1, only spontaneous emission is operable-the limit taken in the case study of the next section. At nonnegligible occupation numbers f 1, stimulated emission and absorption should be considered-in equilibrium (i.e. using Eq. 14), they famously cancel each other in the limit of zero recoil [67,68]. However, at finite ω, absorption always "wins" eventually: where Eq. 14 and conservation of both energy and chemical potential were used to get to the second line. Therefore, any (bound) mode k eventually saturates to a level f (k) = 1/(e ω/T − 1) where F = 0 and detailed balance is reached. (Of course, unbound modes with |k|/m > v esc continue to escape to infinity and never reach this level.) The saturation condition implicitly defines the absorption timescale roughly the time it takes for the stellar core to reach densities of order the second fraction. The effective absorption time for locations far outside the star is parametrically longer for stars with planets, which will secularly perturb the emitted particles' trajectors into nonstar-traversing orbits, and thus diffuse into other parts of phase space with lower occupation numbers. I treated the case for single boson emission; occupation levels for fermions would saturate earlier, at f = 1/(e ω/T + 1). I am not aware of a regime in dense astrophysical media where Bose-enhanced exponential growth can occur (like astrophysical masers in dilute molecular gas clouds), as F st > F abs necessitates population inversion, and thus large departures from thermal equilibrium. Finally, processes with constant matrix elements in the soft limit, i.e. |M| 2 const. as k → 0, are precisely those with n p = 1 that efficiently fill stellar basins. In the soft limit, one can ignore k in the three-momentum delta function and replace ω with m everywhere in Eq. 12. This allows the phase space integral to be factored out d 3 k (m 3 / √ 2π 2 ) dω kω 1/2 k , realizing the form of Eq. 2 with n p = 1.

SOLAR AXION BASIN
As a case study of interest, consider a pseudoscalar particle a (hereafter "axion") coupled derivatively to the axial electron current with dimensionless coupling g aee : If the axion mass m is below twice the electron mass m e , its main radiative decay channel is to two photons (minimally induced from a diagram with an electron loop and the aψ e ψ e vertex from Eq. 17), with a radiative rate assuming m 2m e [69]. The numerical estimate is evaluated near parameters of interest, at which the axion can be considered cosmologically stable, and the Xray photon rate is beyond current observational capabilities [70,71]. The decay to two photons could be even further suppressed (at the expense of tuning), but could also be much higher if the Peccei-Quinn symmetry is anomalous and the axion has a direct coupling to photons [70]. Observable implications of rapid decays of stellar basin axions to photons are left to future work.
In the remainder of this section, I compute the contributions to the local basin density ρ b from bremsstrahlung in the Sun and Earth's core, which are the primary sources for basin axions with masses below 5 keV. Photoproduction (the Compton process) in the Sun is the main production channel at higher masses. (I have not yet included other atomic processes, such as those of Ref. [72], within this framework.) I first calculateQ p for each process (massage Eq. 12 into the form of Eq. 2), and subsequently evaluate the volume integrals of Eq. 4 over the Standard Solar Model of Ref. [73] and the Preliminary Reference Earth Model of Ref. [74].
The results in terms of energy density at Earth's surface are displayed in Fig. 2 for a reference coupling of g aee = 3 × 10 −13 . Evidently, the (coupling-independent) parametric estimate of ρ b /ρ ∞ in Eq. 8 holds, validating the conclusion that the local Solar axion density is dominated by the Solar basin for keV-scale axion masses, at least by one order (possibly four orders) of magnitude. The Solar basin densities remain low enough that one can ignore axion re-absorption at Earth's location (cfr. Eq. 16). The reader not interested in the details of the computation may skip to the next section, where it is shown that the basin energy densities of Fig. 2 are potentially detectable in current direct detection experiments.

Bremsstrahlung
For axion bremsstrahlung production e − + Z j → e − + Z j + a from collisions of electrons with nuclear ions of charges Z j and masses m j , and corresponding momenta p 1 + p 2 = p 3 + p 4 + k, the squared matrix element in the soft limit is: Screening effects are accounted for by the Debye-Hückel scale κ 2 s (4πα/T )( j n j Z 2 j + n e ) ≡ 4πα(n N + n e )/T , proportional to the number density of different ions (n j ) and electrons (n e ). The three-momentum transfer is defined as q ≡ p 2 − p 4 . The ions have low occupation numbers f 2 , f 4 1, and may be considered infinitely heavy relative to electrons so that nuclear recoil can be neglected. That leaves a trivial p 2 integral, p 3 = p 1 − q fixed by the delta function (ignoring axion momenta in the soft limit k → 0, ω → m), and the integral over p 4 can be replaced by one over q. Plugging everything into ) are shown for three benchmark basin lifetimes (Eqs. 9, 10, 11; respectively thin, thick, dashed). Black curves are totals of these subprocesses, including a component ρ B,⊕ b at low energies from bremsstrahlung in the Earth's core (for τ = 1, 5 × 10 9 y). The blue shaded region spans the range of these total predictions; its conservative lower boundary is used to set conservative upper limits on couplings in Fig. 3. Relativistic Solar axion flux energy densities (in the massless limit, no attempt was made to compute finite-mass suppression) are shown as dot-dashed lines. The Compton (ρ C ∞ ) and bremsstrahlung (ρ B ∞ ) values agree with those of Ref. [72], whose luminosity ratios are used to get ρ ABC ∞ that also includes atomic recombination subprocesses. For axion masses m 0.5 keV (and perhaps as low as tens of eV), the local Solar basin energy density exceeds that of the relativistic Solar flux, and may reach up to 1% of the cosmic DM energy density ρDM (gray line) at this reference coupling.
Eq. 12 yields: In the first line, the electrons are taken to be nonrelativistic. The second line makes use of the following definitions: ω p ≡ p 2 1 /2m e as the kinetic energy of the incoming electron, ≡ m/ω p as the ratio of axion mass over available electron kinetic energy, and ξ ≡ κ 2 s /2m e ω p as a dimensionless screening measure.
Nondegenerate Medium.-Electrons in the Solar plasma are nondegenerate, as n e (m e T ) 3/2 , so one can approximate f 3 0 (no Pauli blocking) and Eq. 14 as a Maxwell-Boltzmann distribution: With these approximations and after converting d 3 k (m 3 / √ 2π 2 ) dω kω 1/2 k as at the end of the previous section, comparison of Eq. 20 with Eq. 12 yields: The screening measure ξ is quite small in practice. I find that the integral in Eq. 22 is approximated by the empirical formula 4.7 exp{−1.38s 1/2 − 0.704s − 0.024s 3/2 } for s ≡ m/T 20 to better than ±5% accuracy throughout the most strongly emitting regions of the Solar interior. Axion production from electron-electron bremsstrahlung gives another additive contribution, with the same result provided one make the replacement n N ↔ n e / √ 2, similar as in Ref. [75]. The results from the volume integral of Eq. 5 with substitution of Eq. 22 over the Solar Model of Ref. [73] are plotted in Fig. 2 as the gold curves.
For reference, the total energy loss per unit volume in the massless axion limit is with screening corrections ignored [1]. The luminosity found by integrating over the Solar volume leads to a relativistic flux energy density ρ B ∞ plotted as the gold dashed curve in Fig. 2 (cfr. Eq. 1 and integral above).
Degenerate Medium.-Earth's core also sources a local "planetary basin" of axions via electron-ion bremsstrahlung. It consists primarily of metallic iron, wherein the electron gas is strongly degenerate, with a Fermi energy E F = p 2 F /2m e = µ−m e T . The electron density is n e = p 3 F /(3π 2 ). At standard conditions, iron's mass density is ρ N = m N n N = 7.874 g cm −3 , with Fermi energy E F ≈ 11.1 eV, Fermi momentum p F ≈ 3.37 keV, and electron density n e ≈ 0.411 keV −3 . From charge neutrality, I derive an effective ionic charge Z N n e /n N ≈ 0.636 (the core electrons of iron shield the rest of the nuclear charge). I adopt a Thomas-Fermi approximation for the screening scale: κ s (e/π) √ m e p F ≈ 4.00 keV. One complication is that Earth's core is considerably denser than iron at ambient pressures; I deal with this by keeping Z N fixed and scaling up the density and Fermi levels by appropriate factors. I ignore nickel and other elements, but they would be straightforward to include by summing over N = Fe, Ni, . . . . I take a fiducial temperature for both the inner and outer core of T ⊕ ≈ 5000 K, and the mass densities from Ref. [74]. For a strongly degenerate medium, the result of Eq. 20 simplifies considerably because m/E F may be taken to be much less than unity, as high values would be Boltzmann suppressed. I find the closed-form expression: Electron-electron bremsstrahlung production of axions can be ignored entirely, as it is suppressed by the Pauli blocking factor d 3 p 2 2f 2 (1 − f 4 )/n e 3E F T /p 2 F 1. The associated basin energy density at R = R ⊕ ≈ 6371 km is plotted as green curves in Fig. 2.

Compton Process
At high axion masses significantly above the temperature (and plasma frequency) of the Sun, the Compton process γ + e − → a + e − is the dominant axion production channel, with an energy production per unit of 6D phase space volume of: with the abbreviation P ≡ p γ + p 1 − k − p 2 . The polarization-averaged square matrix element is |M C | 2 | e 2 g 2 aee m 2 /m 2 e in the soft limit. The electrons are taken to be nonrelativistic (E 1,2 m e ), nondegenerate (f 2 0) bystanders that absorb momentum but have negligible energy recoil. Then the p 1 integral over 2f 1 gives a factor of n e , the p 2 integral is trivial, and the p γ integral over f γ and δ(E γ − m) can be done analytically, taking into account the dispersion relation E 2 γ = ω 2 pl + p 2 γ at plasma frequency ω 2 pl = 4παn e /m e . I find: kinematically allowed only when m > ω pl , which is approximately 0.3 keV in the Solar core. The resulting Solar basin energy densities from this process are depicted in Fig. 2 by red curves. For reference, the total energy loss rate per unit volume from the Compton process is as calculated for massless axions in Ref. [75]. Its Solar volume integral leads to a relativistic flux energy density shown as the red dashed curve in Fig. 2.

DIRECT DETECTION SIGNALS
Given that energy densities of the Solar axion basin at Earth's location can be competitive or exceed those of the relativistic (unbound) Solar axion flux, searches for Solar axions must be recast to take into account this "new" population component. To first order, the stellar basin resembles a DM halo with a peculiar density profile centered on the star that is continually replenished. At Earth's location within the Solar System, the velocity dispersion within the basin would be significantly smaller-bounded by [v esc (AU)] 2 ∼ 10 −8 as opposed to [v galactic DM ] 2 ∼ 10 −6 . For all practical purposes, any one axion absorption event from the Solar basin would thus be indistinguishable from an axion DM absorption event.
Dark-matter absorption line searches such as those of Refs. [49,53] that set limits g DM aee on the electron-axion coupling (depicted as thin lines in Fig. 3) may thus be recast as coupling limits g basin aee that do not make use of the assumption that the axions make up the cosmic DM. Suppose the DM event rate is C(g DM aee ) 2 ρ DM for some proportionality constant C. Then the event rate for Solar basin absorption is C(g basin aee ) 4ρ b g=1 τ with the same constant C. Equating the two produces the recasting map: The blue region is conclusively excluded by the recasting procedure of Eq. 28 with the conservative estimate of the gravitational ejection time of Eq. 9, independent of assumptions about dark matter. This limit is a new "Solar basin" interpretation of dark-matter axion absorption line searches by the XENON1T [49,53], LUX [45], and PandaX-II [46] collaborations (upper limits, assuming the axion makes up all DM, depicted as thin dotted lines). The region above the thick blue curve is also disfavored if the fiducial estimate of Eq. 10 is adopted. Both XENON1T's S2-only [49] and S1+S2 [53] searches could have plausibly detected an axion in the Solar basin near or above the dashed blue curve (with the optimistic estimate of Eq. 11), and below anomalous cooling constraints of white dwarfs (WD) and red giants (RG) (orange and gray regions) extrapolated to high axion masses. The tentative excess of Ref. [53] is in gross violation of those cooling constraints when attributed to absorption of a relativistic Solar axion flux (red region), but could be mutually compatible with absorption of Solar basin axions in the m ∈ [1.5, 3] keV range (inside the green rectangle), as well as with hints of anomalous WD cooling (pink region).
Above,ρ b g=1 is the energy injection rate at g aee = 1 (merely a mathematical trick: at such high couplings the Sun would become opaque to axions). The resulting Solar axion basin limit in Fig. 3 is depicted by the light blue region (for the conservative estimate of τ in Eq. 9); the lower blue lines should be viewed as delineating parameter space that is disfavored ("fiducial", with Eq. 10) and where signals are conceivable ("optimistic", with Eq. 11). Fig. 3 displays cooling constraints from WDs [25] and RGs [22] as orange and gray regions respectively. They are crudely extrapolated to finite axion masses by a factor proportional to the inverse square root of ∞ m dω ω 2 √ ω 2 − m 2 /(e ω/T − 1) as in Refs. [76,77]. This extrapolation might be overly stringent, as the temperature variation among white dwarfs in particular will skew cooling analyses based on luminosity functions. Colder WDs are more strongly affected by the finite axion mass than hotter ones, and because they are more numerous, these are also the better constrained observationally. There may be some indications of excess WD cooling, both from luminosity distribution [25,26] and pulsation period drift [28][29][30] analyses, that could potentially be explained by an electron-axion coupling; this parameter space is plotted as the pink region below the constraints. On the other hand, Ref. [78] has contested these cooling hints, leveraging a potentially more sensitive population of hot WDs in a single globular cluster, and claims a limit of g aee 8.3 × 10 −14 , just above the lower boundary of the pink region in Fig. 3.
The low-energy excess of electron-recoil events recently reported by the XENON1T experiment [53] is not compatible with stellar cooling constraints when interpreted as a signal from the relativistic Solar axion fluxsketched as the red "flux excess" region in Fig. 3 for m 1.5 keV and 2.6 × 10 −12 g aee 3.7 × 10 −12 . Note that the conservative Solar axion basin limit (recast from the S2-only DM analysis from the same collaboration) further excludes m 0.7 keV within this space. If the optimistic estimate for the gravitational ejection time scale in Eq. 11 is indeed correct, then the excess of events could instead be due to an absorption line from the Solar basin if the axion mass were roughly between 1.5-3 keV, a region highlighted as the "basin excess" green box in Fig. 3. The global significance for the axion DM line search was reported [53] to be lower primarily due to the look-elsewhere-effect within the search region m ∈ [1, 210] keV; Fig. 3 shows that the region of interest for a Solar basin axion search should instead be m 20 keV, with a correspondingly lower trials factor and higher global significance. With a generous stretch of the imagination, these axion parameters might also be responsible for the WD cooling hint, though the analysis of Ref. [78] disfavors this scenario.
I conclude this section with the main takeaway: axioelectric direct detection experiments are much more sensitive to keV-scale Solar axions than previously anticipated. If simulations confirm the optimistic estimate for the gravitational ejection time of Eq. 11, then Fig. 3 shows that the line search of Ref. [53] already has the best discovery potential to g aee for m ∈ [6, 10] keV, surpassing the strongest claimed WD and RG cooling bounds [22,25,78]. DISCUSSION I have identified and described a previously overlooked process in astroparticle physics. Amusingly, it is not just feeble interaction strengths, which allow for volumetric stellar emission, that make stars interesting sources of weakly coupled particles. A small but non-zero mass opens up the possibility of bound-orbit emission that fills up a "stellar basin" over time, a spectacular channel not available to massless photons. I have performed a case study for CP-odd pseudoscalars coupled to electrons in some detail, and investigated the consequences for direct detection experiments.
It would be interesting to study stellar basins in the context of other particles, couplings, and experiments, as well as for stars other than the Sun. Other direct detection experiments-such as CAST [43], which probes the axion-photon coupling-might also be affected by the Solar basin. Preliminary estimates indicate that axions would quickly saturate to maximum occupation numbers (see Eq. 16) around isolated neutron stars from nuclear bremsstrahlung production, even with incredibly tiny couplings to nuclear matter. This "neutron star axion basin" could potentially boost indirect detection signatures such as those of e.g. Ref. [79,80]. While the effects on initial cooling rates are negligible in practice, the basin could potentially act as a reservoir that slows the cooling process, as the basin is re-absorbed by the star when it cools down. If part of the basin is protected on a long-lived orbit outside the stellar interior, secular perturbations may kick it back inside (or into other surrounding objects) at a later time. The Sun also fills up a "neutrino basin", but at energy densities far below that of even the cosmic neutrino background, rendering detection all but impossible.
Finally, the most pressing loose end in the present analysis is the poorly known gravitational ejection time scale. Simulations of the orbital dynamics (beyond several Lyapunov times) are needed to characterize the expected Solar basin density and its fluctuations, which would be crucial input to direct detection experiments. Solar basin signals (with ρ b ∝ 1/R 4 ) will exhibit stronger annual modulation than Solar flux signals (with ρ ∞ ∝ 1/R 2 ) due to Earth's eccentric orbit, which is also out of phase with the modulation expected from a DM signal. This temporal variation combined with the ultra-narrow kinetic energy spread would be a smoking gun of direct detection signals from the Solar basin.
KVT is especially indebted to Timothy Wiser for sharing insights regarding gravitational ejection dynamics and phase space diffusion, and thanks Asimina Arvanitaki, Masha Baryakhtar, Nathaniel Craig, Isabel Garcia Garcia, Junwu Huang, Andrew Jayich, Ania Jayich, Amalia Madden, Cristina Mondino, Aaron Pierce, Oren Slone, Anna-Maria Taki, and Neal Weiner for discussions and comments. This research is funded by the Gordon and Betty Moore Foundation through Grant GBMF7392, and supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.