Heavy Dark Matter in White Dwarfs: Multiple-Scattering Capture and Thermalization
Nicole F. Bell Giorgio Busoni Sandra Robles Michael Virgato



KCL-PH-TH/2024-20

INT-PUB-24-016

FERMILAB-PUB-24-0196-T






a]Nicole F. Bell,
b,c]Giorgio Busoni,
d,e]Sandra Robles
a]and Michael Virgato
[a]ARC Centre of Excellence for Dark Matter Particle Physics, 

School of Physics, The University of Melbourne, Victoria 3010, Australia
[b]ARC Centre of Excellence for Dark Matter Particle Physics, 

Research School of Physics,
The Australian National University, Canberra ACT 2601, Australia
[c]
ARC Centre of Excellence for Dark Matter Particle Physics, 

Department of Physics, University of Adelaide, South Australia 5005, Australia 
[d]Theoretical Particle Physics and Cosmology Group, Department of Physics, King’s College London, Strand, London, WC2R 2LS, UK
[e]Particle Theory Department, Theory Division, Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA

n.bell@unimelb.edu.au
giorgio.busoni@anu.edu.au
srobles@fnal.gov
mvirgato@student.unimelb.edu.au






We present an improved treatment for the scattering of heavy dark matter from the ion constituents of a white dwarf. In the heavy dark matter regime, multiple collisions are required for the dark matter to become gravitationally captured.  Our treatment incorporates all relevant physical effects including the dark matter trajectories, nuclear form factors, and radial profiles for the white dwarf escape velocity and target number densities. Our capture rates differ by orders of magnitude from previous estimates, which have typically used approximations developed for dark matter scattering in the Earth.
We also compute the time for the dark matter to thermalize in the center of the white dwarf, including in-medium effects such as phonon emission and absorption from the ionic lattice in the case where the star has a crystallized core. We find much shorter thermalization timescales than previously estimated, especially if the white dwarf core has crystallized.
We illustrate the importance of our improved approach by determining the cross section required for accumulated asymmetric dark matter to self-gravitate.







Heavy Dark Matter in White Dwarfs: Multiple-Scattering Capture and Thermalization
    [
    April 24, 2024
=================================================================================








§ INTRODUCTION





White dwarfs (WDs) are the most abundant stellar remnants; in fact more than 90% of the stars in the Galaxy are WDs. Due to their high density, extreme conditions and the existence of observational data, they have recently been used to test and constrain dark matter (DM) models. This generally involves the accumulation of DM particles, which leads to observational signatures such as an increase in their luminosity <cit.> 
or DM-triggered supernova ignition/black hole formation <cit.>. The latter scenario specifically requires the capture of heavy DM of order 100 TeV, or above. For these masses, multiple collisions of the DM with the WD constituents are required for the DM to lose sufficient energy to become gravitationally bound to the star.

Analytical approaches in the literature that deal with capture of heavy DM via multiple scattering in stars <cit.> are based on Gould's seminal work for capture in the Earth <cit.>. It is worth noting that in the case of the Earth, the main contribution comes from the scattering of DM with iron nuclei. 
Gould's formalism was derived under the following assumptions: (i)  DM trajectories are unaffected by collisions, 
(ii) constant escape velocity, 
(iii) constant iron density, 
(iv) DM follows linear trajectories outside and inside the  Earth's core, thereby neglecting gravitational focusing. 
All of these assumptions approximately hold in the Earth's core, where the iron density profile is rather flat, and the  gravitational potential of the Earth is weak, meaning that the DM escape velocity is much smaller than the average DM halo speed. However, in the case of the Earth, a Monte Carlo simulation is in fact best suited to tackle multi-scattering capture of light DM <cit.>. 


For WDs, on the other hand, we have a very different kinematic regime. This enables an analytical approach to be attempted. For capture in WDs, the infalling DM particles are accelerated to velocities that are a sizeable fraction of the speed of light and hence orders of magnitude greater than the DM halo speed. In this regime, the assumptions (ii)–(iv) above do not hold. (Assumption (i) remains well justified.) 
In this paper, we relax all the non-valid simplifying assumptions, and derive a sound formalism for multiple scattering capture in WDs.  We show that the simplifying assumptions would lead to an overestimation of the DM capture rate, with the difference becoming more severe for heavier DM. 

Our formalism is based on a response function that properly encodes the relevant physical effects, including the cross section suppression due to the nuclear form factors; non-constant escape velocity; and gravitational focusing when calculating the probability for DM to undergo more than one collision and lose enough energy in the process to become bound to the WD. We shall begin by considering DM scattering from a single type of WD constituent. 
As we shall see, this method can be generalized to the more realistic case of WDs made of more than one ionic species.

Following gravitational capture, the DM particles will continue to scatter with the WD constituents until an equilibrium state is reached. The timescale for this process to occur has previously been estimated <cit.> by using the average DM energy loss in an energy regime that is suited for capture, i.e. for high energy transfers. We show that due to finite temperature effects the average energy loss has a different scaling with the DM kinetic energy and DM mass for high and low energy momentum transfers. The latter is more important for thermalization. In addition, if the WD is old enough that its core has begun to crystallize, then the ions will form a Coulomb lattice, in particular in the region near the center of the star. In this case, collective excitations of the ionic targets introduce an additional correction which we find to increase the thermalization time by less than an order of magnitude, much less than previously estimated.

Finally, to showcase the effect of using our revised calculations for capture and thermalization in the multiple-scattering regime, we consider the accumulation of heavy non-annihilating dark matter.  Such non-annihilating dark matter is expected in asymmetric dark matter models, where the DM abundance in the present-day Universe consists entirely of DM particles and no DM antiparticles (or vice-versa). 
We shall use our formalism to estimate the size of the DM-nucleon cross section required for non-annihilating DM accumulated in the WD core to self-gravitate. Since we are dealing with a DM population that grows over time, the time evolution of the WD core temperature and hence that of thermalization time is considered. Depending on the composition and compactness of the WD, we find at least one order of magnitude difference with previous results.  

This paper is organized as follows. In section <ref>, we briefly review the physics of WDs relevant to this work and how to model it in order to use observational data. In section <ref>, we derive the response function for DM collisions in the multiple-scattering regime, both for scattering from a single target species, and from multiple targets. The timescale for thermalization of heavy DM is discussed in section <ref>. The condition for self-gravitation of asymmetric DM is examined in section <ref> and our conclusions are given in section <ref>. 








§ WHITE DWARFS



White dwarfs are the most abundant stars in the Galaxy, being the end points of main sequence stars with mass ≲ 8-10   M_⊙, depending on metallicity. This compact stellar remnant is supported from gravitational collapse by electron degeneracy pressure. WDs are born at high temperatures, ∼ 10^8, and cool over billions of years to ∼ 10^5. 
Below we discuss how the internal structure is modeled from an appropriate equation of state, and summarize the observational data used in this work.




 §.§ Evolution and Internal Structure




WDs are born  as extremely hot, compact objects, whose fate is to cool slowly by radiating light. There is no longer nuclear burning in the core of these collapsed objects. The only energy source, responsible for the star's luminosity, is of thermal origin: the kinetic motion of the non-degenerate ions which are completely decoupled from the electrons. The degenerate electrons do not contribute to the cooling process. They determine the star's stellar structure and hence the WD mass radius relation. In addition, as excellent heat conductors, electrons are very efficient at thermalizing the WD core, and hence we can model this region of the WD as an isothermal sphere that accounts for ∼99% of the star's mass. 
As the WD evolves (cools), its ionic constituents undergo phase transitions from a hot gas to a liquid and finally to a crystal, i.e an ion lattice. A significant amount of energy is released in the latter first-order phase transition, which slows down the cooling process. The WD solidifies from the center to the surface. 
The remaining ∼1% of the WD mass is in the form of a heat-blanketing envelope, whose composition, mainly hydrogen and helium, plays a key role in the cooling process.


In ref. <cit.>, we modeled the  WD core as being composed of only one element, namely carbon, by solving the structure equations, the Tolman-Oppenheimer-Volkoff (TOV) equations <cit.> coupled to the 
zero-temperature Feynman-Metropolis-Teller (FMT) equation of state (EoS) <cit.>. 
In the present work, we adopt a more realistic approach to match observations. Specifically, we consider stratified WDs made of two elements, 50% each, namely carbon-oxygen (CO) for WDs of mass ≲1.05, and oxygen-neon (ONe) above this mass threshold. In addition, we use the finite temperature extension of the FMT EoS presented in ref. <cit.>. While the zero-temperature approach is adequate for WDs with masses ≳ 0.7 - 0.8 M_⊙, the predicted mass-radius relation in this approximation begins to deviate from observations for lower masses. These deviations are expected to stem from finite temperature effects <cit.>. 
In the relativistic FMT EoS at finite temperature, the thermal energy of the nuclei is considered when calculating the energy density and pressure within a Wigner-Seitz cell[A Wigner-Seitz cell encloses a cloud of electrons surrounding a finite size nucleus so that the system is neutral.]. The electron energy density and pressure depend on the degeneracy parameter and ultimately on the temperature of the system; for further details see ref. <cit.>. EoSs for C, O and Ne were obtained and later used together with the TOV equations to obtain WD radial profiles.
This includes the radial profiles of the escape velocity, , an example of which can be found in Fig. 1 of ref. <cit.>.








 §.§ Observations








Over the past few decades, our understanding of white dwarfs has greatly been improved with the aid of large-area surveys, in particular the Sloan Digital Sky Survey (SDSS) <cit.> which has produced the largest catalogue of spectroscopically confirmed WDs, see e.g. refs. <cit.>.
Spectroscopy provides measurements of the effective temperature  and surface gravity g of a WD, from which the WD mass  and radius  can be inferred.  
Spectroscopic data also give us information about the composition of the atmosphere of a white dwarf (spectral type). These outer layers are key to determining the core temperature , cooling age and other WD properties, with the help of model atmospheres. 



In Table <ref>, we show three old WDs  within 200 pc from the Sun, distances measured by Gaia EDR3 <cit.>. These WDs have spectral type DA, i.e. a hydrogen rich atmosphere. 
In order to model these stars, rather than simulating their full evolution, one can track how the core temperature and other relevant properties change over time by using evolutionary tracks based on model atmospheres such as those given in ref. <cit.> for CO WDs and ref. <cit.> for ONe WDs [Evolutionary sequences for CO WDs of various masses can be found at: <https://www.astro.umontreal.ca/ bergeron/CoolingModels/> and for ONe WDs at <http://evolgroup.fcaglp.unlp.edu.ar/TRACKS/ultramassive.html>.]. 
The most massive WDs are expected to be made of oxygen and neon, hence their effective surface temperature, mass, core temperature (), cooling age and age at the onset of core crystallization have been calculated using the evolutionary tracks for massive WDs of ref. <cit.>. The remaining less massive WDs, on the other hand, are expected to be a mixture of carbon and oxygen,  so that their corresponding inferred properties were estimated using the cooling models of ref. <cit.>. The evolution of the core temperature as a result of the cooling processes  for the WDs in Table <ref> is depicted in the left panel of Fig. <ref>. On the right panel of this figure, we show the radial density profiles of the WD core for the same stars, obtained  as outlined in section <ref>,  using the finite temperature FMT EoS for the core temperature  at the WD cooling age  listed in Table <ref>.
















§ MULTIPLE SCATTERING CAPTURE










In this section, we review the formalism for single-scattering capture of DM in white dwarfs. Building upon it, we derive a proper treatment for multiple-scattering capture, by including the radial dependence of the escape velocity and target number density, the effect of the nuclear form factors,  
as well as considering that DM trajectories follow the rules of gravity. 
Next, we generalize this formalism to the case of multiple targets. 








 §.§ Single Scattering








We begin by considering the situation where a single collision with the non-relativistic, non-degenerate ions in the WD core is enough for DM to become gravitationally bound to the star. 
In this case, the most general expression for the DM capture rate   C_1 in the optical thin limit is given by <cit.>

    C_1 =
    ρ_χ/m_χ∫_0^R_⋆ dr 4π r^2
         ∫_0^∞ du_χw(r)/u_χ(u_χ) Ω_T^-(w),

with

    Ω_T^-(w)    =   ∫_0^v_esc(r)dv   R_T^-(w→ v),
         
        R_T^-(w→ v)    =    ∫_0^∞ ds ∫_0^∞ dt  32 μ_+^4/√(π) k^3 n_T(r) dσ_Tχ/dcosθ_cmv t/w e^-k^2 u_T^2Θ(t + s - w) Θ(v - |t - s|),

where ρ_χ is the DM density, m_χ is the DM mass, n_T is the target number density,  and w(r)=√(u_χ^2+v_esc^2(r)) and v are the DM velocity before and after the collision, respectively. Kinematic quantities are defined in the center of mass (cm) frame, where θ_cm is the center of mass angle, σ_Tχ is the DM-ion cross section,  
s is the cm speed, t is the velocity of the DM in the cm frame prior to the interaction, k^2 = m_T/2,  m_T is the mass of the target and u_T is the speed of the target before the collision  in the cm frame: 

    u_T^2     =    2 μμ_+ t^2 + 2 μ_+ s^2 - μ w^2,
    μ   =   m_χ/m_T,    μ_± = μ± 1/2.




We assume that the DM  speed distribution  in Eq. <ref>, u_χ, follows a Maxwell-Boltzmann (MB) distribution, which reads <cit.> 

    (u_χ)du_χ=u_χ/v_dv_⋆√(3/2π)(exp[-3/2v_d^2(u_χ-v_⋆)^2]-exp[-3/2v_d^2(u_χ+v_⋆)^2])du_χ,

where  v_d is the velocity dispersion of the halo and  is the WD velocity in the Galactic rest frame. 
 In the T_⋆→0 limit the interaction rate Ω^-,  Eq. <ref>, becomes <cit.>

    Ω_T^-(w) = 4μ_+^2/μ wn_T(r)∫_w(r)|μ_-|/μ_+^(r)dv v  dσ_Tχ/dcosθ_cm (w,q_tr^2(w,v)),

where q_tr is the momentum transferred during the collision. 
Note that the differential cross section depends on the DM initial velocity w and the momentum transfer through the form factors F_i,j^(N,N').  We express the cross section as 

    dσ_Tχ/dcosθ_cm(w,q_tr^2)   =   1/32π(w,q_tr^2)/(m_χ+m_T)^2, 
    (w,q_tr^2)    =   m_T^2/m_N^2∑_i,j∑_N,N' C_i^N  C_j^N' F_i,j^(N,N')(w,q_tr^2),

where N={p,n}, the spin averaged squared matrix elements are given in the basis of  non-relativistic (NR) operators <cit.>, and the C_i are the coefficients that accompany the corresponding NR operators as given in Table <ref>.
The nuclear response functions
F_i,j^(N,N')(w,q_tr^2) are given in appendix C of ref. <cit.>. 
Note that we have assumed fermionic DM colliding with ion nuclei, with interactions described by the 6-dimensional Effective Field Theory (EFT) operators of Table <ref>. Note also that since WDs are composed mainly of ionic species with no nuclear spin, namely ^4He, ^12C, ^16O and ^20Ne, only scattering operators with spin-independent (SI) interactions are relevant.  For these elements in particular, only three response functions are non-zero, namely W_M, W_MΦ^”, W_Φ^”<cit.>.

At low energy, the cross section for scattering on target nuclei, σ_Tχ, can be related to the SI cross section on a single nucleon, σ_Nχ, using 

    dσ_Nχ/dcosθ_cm(w,q_tr^2)=1/32π1/(m_χ+m_N)^2∑_i,j C_i^N  C_j^NF̂_i,j^(N,N)(w,q_tr^2),

where F̂ is the nuclear response function for a nucleon N. Thus, in the low energy limit (w→0,q_tr^2→0) we have 

    dσ_Tχ/dcosθ_cm=(m_χ+m_N)^2/(m_χ+m_T)^2m_T^2/m_N^2∑_i,j∑_N,N' C_i^N  C_j^N' F_i,j^(N,N')(w→0,q_tr^2→0)/∑_i,j C_i^N  C_j^NF̂_i,j^(N,N)(w→0,q_tr^2→0)dσ_Nχ/dcosθ_cm.
 
For operators D1, D2, D5 and D6, for which only the W_M response function is present, this simplifies to 

    dσ_Tχ/dcosθ_cm=(m_χ+m_N)^2/(m_χ+m_T)^2m_T^2/m_N^2 A^2 dσ_Nχ/dcosθ_cm→  A^4 dσ_Nχ/dcosθ_cm,

where the final expression is obtained by taking the large m_χ limit. For D10, the precise coefficient that multiplies the DM-nucleon cross section in the RH expression of Eq. <ref> depends on the target, being γ A^4 with γ=1 for He, γ=1.56 for C, γ=1.02 for O and γ=1.22 for Ne. 





In most cases, for  WDs in the solar neighborhood, we have (r)≫ u_χ and hence w(r)∼(r). Using this approximation, we can rewrite Eqs. <ref> and <ref> in terms of the recoil energy of the target E_R 

    E_R = q_tr^2/2m_T,       E_R = E_R^max(1-cosθ_cm/2),

where

    E_R^max   =    2m_T (r)^2 (m_χ/m_T+m_χ)^2, 
    
        E_R^min   =   0.

Trading θ_cm and the DM final speed v for the recoil energy E_R, we obtain

    dσ_Tχ/dcosθ_cm   =   E_R^max/2dσ_Tχ/dE_R, 
    Ω_T^-(w) 
           =   4μ_+^2/μ wn_T(r) E_R^max(,m_χ,m_T)/2m_χ∫_E_R^min^E_R^maxdE_R dσ_Tχ/dE_R(,E_R).




In the original derivation of Eq <ref>, an average over the DM angular momentum J (angle of incidence) was performed <cit.>. Here,  we reintroduce this average explicitly since it will become relevant in the optically thick and multiple scattering regimes that depend on the specific trajectory and hence  angle of incidence. Note that this does not affect the result for single scattering in the optically thin limit. Reintroducing the average over J  and substituting the last expression into Eq. <ref>, leads to 

    C_1 = ρ_χ/m_χ∫_0^R_⋆ dr 4π r^2 n_T(r) σ_Tχ((r)) ^2(r) ∫_0^1 ydy/√(1-y^2)∫_0^∞ du_χ(u_χ)/u_χ ,

where y=J/J_max, and J_max is the maximum value of J at a given distance from the center of the star. Note that, thanks to the assumption u_χ≪, the integral over u
_χ can be  factored out  from the integral over the volume. 

When the DM-target cross section is larger than a threshold value, the capture rate approaches the geometric (optically thick) limit. In this regime, the flux of DM particles that traverses the star is considerably attenuated, hence the optically thin approximation is no longer appropriate. 
To account for this,  an optical factor η(r) is introduced in Eq. <ref>, which removes captured DM particles from the incoming flux. This η(r) factor,  defined in terms of the optical depth for DM-target  interactions, depends on the trajectory followed by the DM particle until it is captured <cit.>. 
There exist two equally probable paths  a DM particle can follow to reach a given point r within the WD. Therefore, there are two corresponding optical depths τ^-_χ 
(shortest path) and τ^+_χ (largest path) given by <cit.> 

    τ^-_χ(r,y)
            =∫_r^R_⋆dx/√(1-y^2J_max(r)^2/J_max(x)^2)Ω^-(w(x))/v_esc(x)√(1-^2(x)),
    τ^+_χ(r,y)
           =∫^r_min_r + ∫^_r_mindx/√(1-y^2J_max(r)^2/J_max(x)^2)Ω^-(w(x))/v_esc(x)√(1-^2(x)) = 2τ^-_χ(r_min,y)-τ^-_χ(r,y),

where r_min is the position of the perihelion. 
For details on the derivation of these expressions, see Appendix C.1 of ref. <cit.>. 
Then, the optical factor averaged over the two possible trajectories and the DM angular momentum is 

    η(r)=1/2∫_0^1y dy/√(1-y^2)(e^-τ_-(r,y)+e^-τ_+(r,y)).

Finally introducing this factor in Eq. <ref>, the full expression for the capture rate in the single scattering regime reads 

    C_1 = ρ_χ/m_χ∫_0^∞ du_χ(u_χ) /u_χ∫_0^R_⋆ dr 4π r^2 n_T(r) ^2(r) σ_Tχ((r)) η(r).

For large DM masses, however, this expression will overestimate the true capture rate. 
 This is because it assumes that a single collision is enough to capture a DM particle, consistent with the use of E_R^min=0. In the next section, we go beyond this limitation by developing a formalism to handle capture in the multiple-scattering regime. 



 §.§ Response Function for Multiple Scattering





When the DM is heavy enough, Eq. <ref> (which assumes a capture probability ∼1) would overestimate the correct result. In this regime, more than one collision is necessary for the DM to be captured.  
To account for this, we construct a function that incorporates both of the non-linear effects into the capture rate, namely the star opacity and the multiple-interactions with the ionic target species. 


First, we define a probability density function for the energy lost by a DM particle in a collision as 

    f(E_R)= 1/σ_Tχdσ_Tχ/dE_R(E_R).

Then, the probability that the DM loses an amount of energy of at least δ E =m_χ u_χ^2/2 after a  single collision is given by

    F_1(δ E)=∫_δ E^∞dE_R f(E_R),

where we have extended the upper integration limit from E_R^max to infinity. This is valid for white dwarfs, as the escape velocity is very large. 
In the same way, the probability of losing the same amount of energy after N scatterings is 

    F_N(δ E)=∫_0^δ EdE_R  F_N-1(δ E-E_R)f(E_R).


For simplicity, we assume that the DM-target cross section is  well approximated by 

    dσ_Tχ/dcosθ_cm∝ e^-E_R/E_0,
 
where E_R is the recoil energy and E_0 depends on the specific nuclear target. 
That is, we assume exponential nuclear form factors similar to the Helm approximation.   
This leads to 

    f(E_R)   =   Θ(E_R)/E_0e^-E_R/E_0,
     F_1(δ E)    =    e^-δ E/E_0.

Defining the dimensionless quantity

    δ = δ E/E_0=m_χ u_χ^2/2 E_0,

and taking the  Laplace transform of the F functions written in terms of δ, we find 

    F̃_1(s)=1/1+s,     F̃_N(s)=1/(1+s)^N,

where the last expression corresponds to 

    F_N(δ) = e^-δδ^N-1/N-1!.


To account for multiple collisions, we define the probability that the DM undergoes N scatterings with nuclei. For DM passing through the same path length τ_χ, this is given by the Poisson distribution <cit.>
 
    p_N(τ_χ) = e^-τ_χτ_χ^N/N!.

Substituting F_1(δ) and p_0 into Eq. <ref>, we obtain for single scattering
 
    C_1 = ρ_χ/m_χ∫_0^R_⋆ dr 4π r^2 n_T(r) σ_Tχ((r)) ^2(r)∫_0^1 ydy/√(1-y^2)∫_0^∞ du_χ(u_χ) /u_χp_0(τ_χ) F_1(δ).

This is equivalent to Eq. <ref>, which accounts for the star opacity, once we average over the two equally probable optical depths, τ_χ^+ and τ^-_χ. Note that F_1 indicates the probability for DM to be captured after a single collision.
We can now generalize the above expression to the case of N scatterings as
 
    C_N = ρ_χ/m_χ∫_0^R_⋆ dr 4π r^2 n_T(r) σ_Tχ((r)) ^2(r)∫_0^1 ydy/√(1-y^2)∫_0^∞ du_χ(u_χ)/u_χp_N-1 (τ_χ)  F_N(δ).

The total capture rate is given by the sum over all N collisions,

    C=∑_N C_N.



Next, instead of first evaluating the integrals in Eq. <ref> and then summing over N, we sum the series first by introducing the response function, G(τ_χ,δ)

    G(τ_χ,δ)    ≡   ∑_N=1^∞ p_N-1(τ_χ)  F_N(δ) = ∑_N=1^∞e^-τ_χτ_χ^N-1/(N-1)!e^-δδ^N-1/(N-1)!
       =    e^-τ_χ-δI_0(2√(τ_χδ)),

where I_0 is the  modified Bessel function of the first kind I_n for n=0. This response function encodes the probability to lose an amount of energy of at least δ E through multiple collisions. Therefore, G(τ_χ,δ)dτ_χ is the probability that DM collides with an ionic target and loses at least an energy δ E in the scattering process, after travelling a differential length δτ_χ in the medium, from a region of the stellar interior with optical depth τ_χ. 

To account for the two possible optical depths, given by Eqs. <ref> and <ref>, we average over them and over the DM angular momentum, as in section <ref>. This leads to

    G̃(r,δ) = ∫_0^1 dy y dy/√(1-y^2)G(τ_χ^-(r,y),δ)+G(τ_χ^+(r,y),δ)/2.

Note that the response function depends on the radial coordinate through the optical depth and on the DM velocity at infinity, u_χ. Finally, the total capture rate that accounts for the stellar structure, star opacity and multiple scattering reads

    C = ρ_χ/m_χ∫_0^R_⋆ dr 4π r^2 n_T(r) ^2(r) σ_Tχ((r))  ∫_0^∞ du_χ(u_χ) /u_χG̃(r,m_χ u_χ^2/2E_0),

where the average over the DM angular momentum has been absorbed in the definition of G̃(r,m_χ u_χ^2/2E_0).

Recall that we assumed exponential nuclear form factors to derive Eq. <ref>. 
As in ref. <cit.>, 
we now introduce the nuclear response functions given in appendix C of ref. <cit.> in Eq. <ref>. In principle, we can proceed in a similar fashion to the exponential form factor case and derive their corresponding F_N functions, as well as the response function for multiple scattering, which is more involved and is given in terms of generalized hypergeometrical functions. A way to circumvent this is to approximate the form factor with an exponential. 
 For instance, if we consider WDs made of carbon, and interactions described by the scalar-scalar operator (D1), we have 

    dσ_TχdE_R   ∝    e^-E_R/E_1(1-x_0E_R/E_1)^2,

with x_0∼0.22 and E_1∼1.25, where these factors come from the nuclear response function for ^12C. 
This function is well approximated by e^-E_
R/E_0
where  E_0 is given by 

    E_0 = lim_m_χ→∞∫_0^∞ dE_R dσ_Tχ/dE_R E_R∫_0^∞ dE_R dσ_Tχ/dE_R.

This  value of E_0 is then used in the response function G̃. We have checked that this approximation holds for the spin independent operators D1 and D5 for carbon, oxygen and neon targets. 






In Fig. <ref>, we compare our approach (magenta) 
with previous calculations in the literature that made use of simplifying assumptions, namely those of ref. <cit.> (light blue) and ref. <cit.> (orange). We consider scalar-scalar interactions of 10^3 and 10^6 mass DM, for a 1 WD made of carbon,  located in the solar neighborhood, i.e., we assume ρ_χ=0.4/^3,  =220/s and v_d=270/s. Radial profiles for the carbon number density and escape velocity for this WD were obtained in ref. <cit.>, 
Our results are given in terms of the maximal capture rate that can be achieved, the so-called geometric limit, which is given by
  
    C_geom   =   π R_⋆^2 ρ_χ/m_χ∫_0^∞w^2()/u_χ(u_χ)du_χ,
       =   π R_⋆^2 ρ_χ/3v_⋆ m_χ[(3 v_esc^2()+3 v_⋆^2+v_d^2) (√(3/2)v_⋆/v_d) +√(6/π) v_d e^-3^2/2 v_d^2].



As can be seen from Fig. <ref>, approaches that do not incorporate the radial dependence of the escape velocity, the ionic target number density and form factors, as well as a DM-target relative velocity distribution, overestimate the capture rate by at least one order of magnitude (orange lines) for DM-nucleon cross sections  σ_pχ≲10^-42^2. We note that ref. <cit.> assumed that the nuclear form factor saturates to ∼0.3, while ref. <cit.> considered 0.5, which we find is the main cause of the large discrepancy. Recall that we have retained the full radial dependence of the form factors and recoil energies, rather than averaging them over the NS radius. We find that F_i,j^(N,N') saturates to  ∼0.034,  for the scalar-scalar operator and recoil energies lower than  order MeV. 
Note also that the approach of ref. <cit.> does not reach the geometric limit for large cross sections; this occurs because they 
truncate the sum over C_N at a certain maximum value of N,  N_max, 
which neglects important terms when the cross section is sufficiently large. 
There is another substantial difference when comparing our results with those calculated using the prescription of ref. <cit.> (light blue lines), which neglects the motion of the WD. That approach overestimates both the capture rate and the threshold cross section at which the capture rate saturates the geometric limit.


Regarding the remaining operators in Table <ref>, operators D2 and D10 are momentum suppressed in the non-relativistic limit, meaning that their matrix elements are proportional to the squared momentum transfer q_tr^2. Therefore, the approximation in Eq. <ref> does not hold. Following the procedure outlined above for non-suppressed interactions, we find that the response function takes the form of a linear combination of hypergeometric functions. For D6, whose cross section also depends on the DM velocity, we can still use the response function in Eq. <ref>, but in this case the value of E_0 will depend on u_χ. 









 §.§ Multiple scattering with multiple targets






The response function obtained in the previous section is valid for white dwarfs made of a single ionic species. However, the core of real WDs is composed mainly of two elements plus small traces of other elements. 
In this section, we  generalize our previous result to the case of scattering on two ionic species. 
Each target i will have a different optical depth, τ^i_χ; the total optical depth τ_χ is simply the sum of the contributions from each of the target species. 



First, as in the previous section, we consider the probability for DM to interact with a target i and lose energy of at least δ E, while travelling a length dτ_χ^i, starting from a layer in the WD with optical depth τ_χ^i. This 
is given by the differential element G(τ_χ^i,δ_i)dτ_χ^i, where 

    δ_i =δ E/E_0^i,

and the energy scale E_0^i depends on the target i and is calculated using Eq. <ref>. Thus,  the probability to interact and lose the same amount of energy when DM travels a path-length  τ_χ^i is simply  the integral of the differential element over the trajectory, i.e. 

    G(τ_χ^i,δ_i) = ∫_0^τ_χ^i dτ  G(τ,δ_i).

Note that this probability satisfies all the expected properties:

    lim_τ_χ^i→ 0 G(τ_χ^i,δ_i)    =    0,
    lim_τ_χ^i→∞ G(τ_χ^i,δ_i)    =   1, 
    lim_δ_i→ 0 G(τ_χ^i,δ_i)    =    1-e^-τ_χ^i  .

The first property applies when the DM particle travels through a region that does not contain the target i. The second limit tells us that by taking a sufficiently large optical depth, it is possible to make the probability to interact and lose an arbitrary amount of energy as close to 1 as desired. Finally, the third property gives us the probability that DM  interacts with the target i; this expression takes the expected form, given that e^-τ_χ^i is the probability of no interaction with the medium.

Note that G(τ_χ^i,δ_i), can be rewritten as

    G(τ_χ^i,δ_i) = ∫_δ_i^∞ dδ  k(τ_χ^i,δ),

where k(τ_χ^i,δ) is a probability density function, obtained  by differentiating Eq. <ref>

    ∂/∂δ_i G(τ_χ^i,δ_i) = - k(τ_χ^i,δ_i).





Next, we introduce a second target species. In the presence of these two ionic targets, the cumulative probability of DM to lose an energy δ E after travelling an optical depth τ_χ^i in the target i and τ_χ^j in  the second target j  is found to be 

    G_2,ij(δ E) = ∫_0^δ E/E_0^j dz   G(τ_χ^i,δ E - z E_0^j/E_0^i) [-∂/∂ z G(τ_χ^j,z)].
 
Since the previous expression is a cumulative probability, 
to calculate the capture rate when the last scattering is over the target i, we differentiate Eq. <ref> with respect to τ_χ^i to obtain[Note that, technically, we should differentiate Eq. <ref> with respect to τ_χ^j. However, once we add  G(τ_χ^i,δ_i)e^-τ_χ^j (which is already differentiated with respect to τ_χ^i) to Eq. <ref> and arrive at Eq. <ref>, we realize that the integral of the latter is symmetric in the indices i,j. Thus,  one can differentiate with respect to either τ_χ^i or τ_χ^j, so we chose the former for convenience.] 

    G_2,ij(δ E) = ∫_0^δ E/E_0^j dz   G(τ_χ^i,δ E - z E_0^j/E_0^i) [-∂/∂ z G(τ_χ^j,z)].

This is the probability for DM to be captured between τ_χ^i and τ_χ^i+dτ_χ^i. 



We are now ready to obtain the capture rate for scattering with two targets. In the previous expression, we have considered that the DM energy is reduced below the threshold for capture to occur after subsequent scatterings with two targets i and j; however, capture can happen after collisions with only one ionic species. This is accounted for in the following expression 

    G_12,ij = G(τ_χ^i,δ_i)e^-τ_χ^j + G_2,ij,

where the first term represents capture after interactions only with the element i. 
As in section <ref>, we average Eq. <ref> over the 2 possible trajectories τ^-_χ and τ^+_χ and the DM angular momentum to obtain the multiple target response function G̃_12,ij. 
Thus, the capture rate for the case where scattering from species i dominates is given by

    C_i = ρ_χ/m_χ∫_0^R_⋆ dr 4π r^2 ^2(r) ∫_0^∞ du_χ(u_χ) /u_χ∑_i n_i(r)σ_iχ(r)G̃_12,ij(r,m_χ u_χ^2/2 E_0^i).

The contribution of the second target, C_j, is obtained by interchanging i with j in Eq. <ref>. Adding the two contributions results in a total capture rate of 

    C_tot = C_i+ C_j.

In principle, this procedure can be generalized to any number of targets.

 






In Fig. <ref>, we compare our approach  with that of ref. <cit.> which addresses multi-target multiple-scattering capture in stars, neglecting the stellar structure. We consider a 0.94 WD, which is composed mainly of carbon and oxygen. The radial profiles of the constituent number densities were obtained by running a simulation of the evolution of a 6 main sequence star with solar metallicity, using 
the public stellar evolution code Modules for Experiments in Stellar Astrophysics () <cit.> package, version 21.12.1. 
As in the previous section, we assume the D1 EFT operator, since ref. <cit.> only deals with the case of constant DM-nucleon cross section. As we can see, the difference between the two approaches is not as striking as that observed in Fig. <ref>, with capture rates from ref. <cit.> (light blue lines) being at most a factor ∼2 larger than our results (magenta lines) for m_χ=10^6 (right panel). This is due to the fact ref. <cit.> adopted the parametrization of the Helm form factor given in ref. <cit.>, and averaged the form factor and scattering angles for every DM mass considered, which yields more realistic results than those of refs. <cit.> and <cit.>. Even so, the discrepancy with our results increases with the DM mass. 








In Fig. <ref>, we show the contribution of each
WD constituent to the total capture rate for the WD in Fig. <ref>. We immediately notice that at large cross sections, in the region where the capture rate saturates the geometric limit,  scattering off carbon dominates the capture rate for both DM masses considered. This is due to the fact that in the outer layers of the WD core there are significantly more carbon than oxygen targets, and in this regime capture occurs closer to the surface. For smaller cross sections,  we observe that 
collisions with oxygen nuclei contribute a significant fraction of the total capture rate, even surpassing the carbon contribution for m_χ=10^3.  This occurs because, for smaller cross sections, DM can be captured deeper inside the core, where the oxygen capture rate exceeds its carbon counterpart, thus resulting in a more sizeable oxygen contribution. It is also worth noting that oxygen and carbon form factors do not suppress the capture rate to the same degree; the  oxygen form factor gives a larger suppression than that of carbon at high recoil energies, and less at low recoil energies.  








§ THERMALIZATION






Once captured, the DM will continue to scatter within the WD, gradually losing energy until it thermalizes in the innermost regions of the star. 
In section <ref>, we outline the calculation of the thermalization timescale for the DM-ion interactions listed in Table <ref>. As WDs cool, their ions begin to crystallize, forming an ion lattice. The effect of the lattice structure on the thermalization time is discussed in section <ref>.






 §.§ Thermalization time



The thermalization process can be broadly divided into two stages. The first stage is where the orbit of the captured DM is not wholly contained within the star. The second stage begins once the orbit is fully contained within the star and ends when the DM has reached thermal equilibrium at the WD center. The first stage typically proceeds significantly faster than the second, making up less than 1% of the total thermalization time <cit.>. Therefore, we focus on the second stage to determine the thermalization time.  
During this last stage, there are two distinctive kinematic regimes depending on whether the DM or target velocity dominates the relative velocity of the interactions. 
The core of a WD is expected to be quite homogeneous in terms of density (see right panel of Fig. <ref>) and temperature. As most of the time required to thermalize is confined to the second kinematic regime when DM reaches small energies, we can consider interaction rates in the central region of the WD. 



 To estimate the thermalization time, we use the temperature dependent differential interaction rate Eq. <ref>, which accounts for the thermal motion of the nuclei. 
 Summing over the average time between interactions until the energy transferred by the DM reaches the core temperature of the WD <cit.>, ⟨Δ K_χ⟩ <, we obtain  

    = ∑_n=0^N 1/Ω^-(w(x_n)),

where x_n is the ratio of the DM kinetic energy K_χ to the WD core temperature  after n scatterings, 
and the sum ends for N such that x_N<1. Each subsequent energy transfer is computed by first calculating the average energy lost in each collision using

    ⟨Δ x ⟩   =   E(x)/Ω^-(w(x)),
    
        E(x)    =   ∫_0^x dy (x-y) R^-(x→ y),
 
where the normalized differential interaction rate (Eq. <ref>) is given in terms of  x=K_χ/T_⋆ and represents the probability distribution function for scattering from an initial energy x T_⋆ to a final energy y T_⋆. 
From now on, and in Appendix <ref>, we will indicate with Ω^-(x) the interaction rate as a function of the dimensionless variable x, and similarly, we will indicate with R^-(x→ y) the differential scattering rate from an initial energy x to a final energy y.



While the summation in Eq. <ref> can formally be evaluated analytically for some interactions, in general it must be computed numerically. 
When doing so,  numerical instabilities  in the computation of the thermalization time arise from the small nuclear velocities and large DM masses. A sounder approach is then to 
analytically calculate the dominant terms of the interaction rate at high and low (close to thermalization) energy transfers, and use these results to approximate the effects of the thermal motion of the targets. 


To obtain analytic expressions for the thermalization time, we can instead approximate the sum in Eq. <ref>  to an integral <cit.> 

    dx/dt   =    E(x), 
       =   ∫_1^∞dx/E(x) ,

where we have set the initial energy to infinity for simplicity, as the initial energy of the dark matter will be significantly larger than when it is thermalized. 
The validity of this approximation improves for smaller fractional energy losses as an increasing number of collisions are required to thermalize the DM. Specifically, at low energies, we have 
⟨Δ x⟩≈ O(1)×√(x/μ), while at high energy we have ⟨Δ x ⟩≈ 2 x/μ;  both approximations are more accurate for large DM masses. 
Note also that we neglect the up-scattering rate. We expect up-scattering to be negligible at large energies and only start to play a role close to thermalization, affecting the result by an 𝒪(1) factor. Hence, all the results presented in this section should be taken as order of magnitude estimates.

To further simplify the analysis, we neglect the effect of the nuclear form factors on the thermalization time. We discuss the effect of including them in Appendix <ref>. In summary, the exponential suppression of the form factors is relevant for energy losses greater than ⟨Δ K_χ⟩∼ 2 K_χ/μ≳ E_0≈𝒪(). As the thermalization time is dominated by the interactions at lower energies, 
it is not impacted by the form factors.




Eq. <ref> can be evaluated analytically for two distinct scenarios: when K_χ≫ T_⋆ (high energy) and when K_χ∼ T_⋆ (low energy), which correspond to x≫1  and x∼ 1, respectively (see Appendix <ref>). 
We set E(x) equal to the sum of the high-energy and low-energy contributions to obtain an approximate expression for the energy transfer E(x).  
For differential cross sections proportional to powers of the DM-ion relative velocity  v_rel^2m, i.e., cross sections of the form 

    dσ_T χ/dcosθ_cm = σ_T/2 v_rel^2m,

where σ_T is a proportionality constant, 
we obtain 

    E(x) ∼  n_T^c σ_T  v_T^2m+1√(x/μ)[2( x/μ)^m+1+1/√(π)Γ(m+3/2)   ],

where v_T^2=3/m_T.  
Similarly, for differential cross sections  proportional to powers of the momentum transfer q_tr^2m, 

    dσ_Tχ/dcosθ_cm = σ_T(m+1)/2^m+1q_tr^2m/q_0^2m,

where σ and q_0 are constants, 
we find

    E(x) ∼  2 n_T^c σ_T  v_T (m+1/m+2)
        (2 v_T^2 m_T^2/q_0^2)^m√(x/μ)[2( x/μ)^m+1+1/√(π)Γ(m+3/2)   ].



  We can now calculate the corresponding thermalization times using Eq. <ref>; see Appendix <ref>. 
 It is worth noting that the high energy term does not affect this timescale; it is only necessary to ensure that the integral converges at large energies. 
The resulting thermalization time for differential cross sections proportional to powers of the relative speed v_rel^2m is given by

    ∝μ/n_T^cσ_T v_T 1/v_T^2m,

where n_T^c=n_T(0) is the ion number density at the WD center. For interactions with cross sections proportional to powers of the momentum transfer q_tr^2m, we find that the  thermalization time reads

    ∝μ/n_T^cσ_T v_T(q_0^2/2 v_T^2 m_T^2)^m.


The precise value of the  𝒪(1) proportionality constant  in  Eqs. <ref> and <ref> is determined by the coefficients of the expansion of E(x) in x/μ at all orders, whereas we have only included the highest and lowest order. The exact thermalization time also depends on the up-scattering rate, which we expect to have an impact near the point of thermalization (x∼ 1). 
It is worth noting that neglecting the low energy term contribution in E(x),  would have resulted in an incorrect scaling of the thermalization time with respect to μ at large mass, with a scaling of ∝μ^3/2+m instead of the correct scaling of ∝μ. 



To derive  for the operators in Table <ref>, we use the proportionality constants σ_T and q_0 in Eqs. <ref> and <ref>, their corresponding thermalization times given by Eqs. <ref> and <ref>, Eq. <ref> to relate the DM-nucleus to the DM-nucleon cross section, and the expressions for the differential DM-nucleon cross sections for each operator in the non-relativistic limit given in Appendix A of ref. <cit.>. 
We obtain

    1   ∼   π  m_χ m_N^2 Λ^4/n_T^c (c_N^S)^2√(1/3m_T^9T_⋆), 
    2   ∼   2π m_χ^3 m_N^2 Λ^4/3n_T^c (c_N^S)^2√(1/3m_T^11 T_⋆^3),
    5   ∼   π m_χ m_N^2Λ^4/n_T^c (c_N^V)^2√(1/3m_T^9T_⋆),
    6   ∼   π m_χ m_N^2 Λ^4/2n_T^c (c_N^V)^2√(1/3m_T^7 T_⋆^3),
    10   ∼   π m_χ m_N^2Λ^4/12γ n_T^c (c_N^T)^2√(1/3m_T^7 T_⋆^3),

where γ in Eq. <ref> is a constant that depends on the specific target (see section <ref>). 




 §.§ Effect of the lattice structure





If crystallization has already started  at the very central region of the WD then, because the momentum transfer becomes increasingly smaller during thermalization, the lattice structure of the WD becomes relevant. In this regime, we can no longer assume DM scattering off a free gas of nuclei, and instead must consider the Coulomb lattice of ions as a collective medium. 
This amounts to considering  phonon excitations within the lattice. To account for these effects, we include the dynamic structure function of the lattice S(q_tr) in  the differential cross section 

    dσ_Tχ/dcosθ_cm→ S(q_tr)dσ_Tχ/dcosθ_cm.

This structure function for a Coulomb lattice is well known <cit.>, and consists of contributions from elastic Bragg scattering and inelastic phonon scattering/absorption. However, only the latter contribute to this current application. For momentum transfers greater than the Brillouin zone, q_B = (6π^2 n_T)^1/3, the structure factor can be approximated as

    S(q_tr) ≈ 1 - e^-2 W(q_tr),

where W(q_tr) is the Debye-Waller factor, which for body centered cubic lattices is given by 

    W(q_tr)  = ⟨ r_T^2⟩ q_tr^2/6,

where ⟨ r_T^2 ⟩ is the mean squared separation of ions in the Coulomb lattice. An analytical fit to W(q_tr) with 1% error is given in  ref. <cit.> as

    W(q_tr)    ≃q_tr^2/2m_Tω_p (1.4 e^-9.1 /ω_p+13 /ω_p) = q_tr^2/2q_sup^2,
    
        q_sup^2     = m_T ω_p/1.4 e^-9.1 /ω_p+13 /ω_p, 
    ω_p     = √(4π Z^2 e^2 n_T/m_T),

where ω_p is the ion plasma frequency. 

If  q_tr^2≪ q_sup^2, one can power expand the exponential in the structure factor,  and therefore the interaction rate effectively becomes suppressed by an additional power of q_tr^2. Following a similar treatment to the case of cross sections proportional to q_tr^2, we find that Eqs. <ref> and <ref> now receive an additional suppression factor[Note that we neglect the additional 𝒪(1) factor 2(m+1)/(m+2).] 

    q_sup^2/6 m_T .

This amounts to increasing the thermalization time by an O(1) factor that depends on the equilibrium temperature when accounting for phonon emission. 
Note that this approximation is valid at O(6m_T T/q_sup^2) since we have performed an expansion of the Debye-Waller structure factor. When 6m_T T/q_sup^2>1, which occurs for some period of time after the onset of crystallisation in the lighter WDs in Table <ref>, a sound approximation for the thermalization time is obtained by multiplying  Eqs. <ref> and <ref>  by
max[1,q_sup^2/(6 m_T )].



Introducing the factor in Eq. <ref>  into Eqs. <ref>-<ref>, we find the following thermalization times for the SI EFT operators

    1   ∼   π  m_χ m_N^2 Λ^4 q_sup^2/6n_T^c (c_N^S)^2√(1/3m_T^11T_⋆^3), 
    2   ∼   π m_χ^3 m_N^2 Λ^4 q_sup^2/9n_T^c (c_N^S)^2√(1/3m_T^13 T_⋆^5),
    5   ∼   π m_χ m_N^2Λ^4 q_sup^2/6n_T^c (c_N^V)^2√(1/3m_T^11T_⋆^3),
    6   ∼   π m_χ m_N^2 Λ^4 q_sup^2/12n_T^c (c_N^V)^2√(1/3m_T^9 T_⋆^5),
    10   ∼   π m_χ m_N^2Λ^4 q_sup^2/72γ n_T^c (c_N^T)^2√(1/3m_T^9 T_⋆^5).









Fig. <ref> illustrates the variation of the thermalization timescale across the evolution (cooling time t_ cool) of the 1.11 WD SDSS J232257.27+252807.42, using the temperature evolution in the left-hand side of Fig. <ref>. We consider the scalar-scalar operator (D1), DM of mass m_χ=10^6 and a DM-proton cross section of σ_pχ=10^-41^2. Before the core of the WD begins  solidifying (at t_ cool < 2.73 for this particular WD; see Table <ref>) we use Eq. <ref>, and after the onset of crystallisation 
we use Eq. <ref>. Note the sudden increase in the thermalization time (magenta line) once the crystallisation  front starts moving from the WD center onwards. 

For comparison, we also show in Fig. <ref>  results using other prescriptions in the literature. The prescription  for the second stage of thermalization from  ref. <cit.>, which  accounts for the effect of the lattice structure so that it holds after the onset of crystallisation, is shown in light blue. This result is orders of magnitude greater than ours. This is mainly due to the fact that this prescription is based on the  expression for the average energy loss at high energy, which, as mentioned in section <ref>, scales as  ⟨Δ K_χ⟩∼ 2 K_χ/μ. Conversely, we find that at the energy relevant for thermalization and, in particular, for phonon emission and absorption (i.e., the low energy transfer regime in which the thermal motion of the targets – which is included in the scattering rate – cannot be neglected) the average energy loss is ⟨Δ K_χ⟩≈ O(1)×√(K_χ/μ). 
Results from the prescriptions in  refs. <cit.> and <cit.>, which are quite similar and do not consider the lattice structure factor, are depicted in orange and green, respectively.  
Since these two  studies do not account for phonons, comparison with our results should be fair for t_ cool < 2.73 (right hand side of the dashed blue line). 
The main difference between these two calculations is a factor A^2 missing in the definition of the DM-target cross section in ref. <cit.>, which is also one of the main sources of discrepancy with our findings. 
The other main  difference with our results is that both refs. <cit.> and <cit.> use the DM average energy loss at high energy, thereby neglecting finite temperature effects.







§ SELF-GRAVITATION OF ASYMMETRIC DM WITHIN WDS





In this section, to showcase the importance of using a proper formalism for capture and thermalization of heavy DM in WDs, we consider the case of asymmetric DM that accumulates in the WD core over time by scattering with the ionic target species. Further collisions, following capture, allow the DM to thermalize to the WD core temperature. In this way, DM settles in the WD interior,  forming a DM sphere close to the center of the star. The general form of the potential energy of this DM cloud, including the WD gravitational potential and self-gravitation, is given by <cit.>

    U_N_χ = -4π G m_χ N_χ(t)∫_0^∞ M(r)n_χ(r) r dr - 4π G m_χ N_χ(t) ∫_0^∞ M_χ(r) n_χ(r) r dr,

where N_χ is the number of accumulated DM particles. The quantity M(r) is the WD mass enclosed in a radius r, while M_χ(r) is the mass of DM particles enclosed in the same radius, given by

    M_χ(r)=4π m_χ N_χ(t)∫_0^r n_χ(r') r'^2 dr'.

The number density of DM particles, n_χ, follows a Maxwell-Boltzmann distribution 

    n_χ(r) = exp[-m_χϕ(r)/]∫_0^∞ 4π r^2 exp[-m_χϕ(r)/]≃1/π^3/2 r_χ^3exp[-r^2/r_χ^2],

where ϕ(r)=-∫_r^∞ G M(r')/r'^2 dr' is the gravitational potential and 

    r_χ=√(3 /2π Gρ_c m_χ)

is the scale radius of the DM sphere, which can be obtained using the virial theorem, in the absence of the self-gravitation term.  
In the right-hand side of Eq. <ref>, we have assumed a constant density in the innermost regions of the WD core <cit.>.  
Using this expression in Eq. <ref>, 
the mean potential energy per DM particle from Eq. <ref> reads 

    U =  -2 π G ρ_c m_χ r_χ^2- G m_χ^2 N_χ(t)/√(2π) r_χ,

where ρ_c=ρ(0) is the central density of the WD. 




The DM core will grow in mass as the accretion process continues, because asymmetric DM does not self-annihilate. 
When considering the full expression for the potential energy per DM particle (Eq. <ref>) in the virial theorem, we find that the necessary condition for collapse is [This condition is obtained from the largest real solution of the cubic equation that results from the virial theorem.]  

    N_χ(t) ≥4√(2)π^3/2r_χ^3ρ_c/3√(3) m_χ = N_crit.

The number of DM particles that have thermalized within the WD center is given by <cit.>

    dN_χdt =  C(t) ( 1 + ddt)^-1,         t≥(t) 
    
        0,         t < (t).

Integrating this expression, we can calculate the number of  DM particles accumulated within the WD throughout its lifetime. Comparing this result with the collapse condition of Eq. <ref> allows us to obtain the cutoff scale Λ required for collapse and, hence, the corresponding DM-nucleon cross section. 

Note that, in principle, the capture rate C depends on time through the WD temperature (see Fig. <ref>). We have checked that for the DM mass range considered here, C(t) can be taken as constant and equal to the value in the →0 approximation derived in section <ref>. The thermalization time is also a function of the cooling time through the WD core temperature, as shown in Fig. <ref>. This time dependence of  will be accounted for in our results.






In Fig. <ref>, we show the DM-proton cross section required to reach N_crit as a function of the DM mass, for two of the local, old WDs in Table <ref>, assuming DM-nucleon scalar-scalar interactions. We arrive at these conclusions by integrating Eq. <ref>, using the radial profiles for the target number densities obtained in section <ref> to compute the capture rates given in section <ref>, together with the thermalization time in Eqs. <ref> or <ref> using the WD cooling curves given in the left panel of Fig. <ref>. 




We note that the time derivative of the thermalization timescale in Eq. <ref>, which is proportional to m_χ,  starts to have an impact for m_χ≳10^12 (0.62 WD) and m_χ≳10^11 (1.11 WD). 
As we move to heavier DM masses, we reach a DM mass above which the thermalization is greater than the age of the WD. This is shown by the shaded grey area in Fig. <ref>, and it is responsible for the changes of slope in σ_pχ. 
We compare these results with the prescription given in ref. <cit.> (light blue lines) which neglects the requirement that the DM particles captured in a  time interval dt must thermalize before being added to the DM cloud that has sunk to the center of the WD. The main source of discrepancy with our results, responsible for cross sections that are smaller by ∼1 (top panel) and ∼3 (bottom panel) orders of magnitude, stems from the computation of the capture rate. This highlights the importance of including the relevant physics of our more complete treatment.  Specifically, including the radial profiles of the escape velocity and target number densities, as well as the nuclear form factors, in the multiple scattering  capture rate, together with a proper estimation of the thermalization time. It is the latter which imposes an upper bound on the mass of non-annihilating DM that can reach self-gravitation in a WD. 














§ CONCLUSIONS









In this paper we have revisited the calculation of the  capture rate of heavy dark matter (DM) in white dwarfs (WDs). In this mass regime, more than one collision is required for the DM to become gravitationally bound to the star. We have extended our formalism for DM capture from the single-scattering to the multiple-scattering regime, by introducing a response function that encodes the cumulative probability for DM to lose an amount of energy of at least δ E through multiple collisions.
Our treatment incorporates gravitational focusing, nuclear form factors, the variation of the escape velocity along the WD interior and the DM-target relative velocity. We have shown the inclusion of these effects, which are often neglected in the literature, alters the capture rate by orders of magnitude, with the corrections being more critical for very heavy DM.
The response function method is able to handle DM-nucleon interactions that are momentum or velocity suppressed, and can be extended to the case of DM-capture via collisions with multiple targets. 

Following capture, the DM will continue to scatter in the star, progressively losing energy until it settles at the center of the star. Eventually, it will reach thermal equilibrium. 
We have estimated the time required for this steady state to be reached, for both a non-crystallized and a crystallized core. In the case of the latter, in-medium effects such as phonon emission and absorption delay thermalization, because the final stages of that process are characterized by low-momentum-transfer DM-target interactions. However, we find that this delay amounts to less than an order of magnitude increase in the thermalization time, which is much smaller than previous estimates. 


Finally, to highlight the importance of correctly calculating the multi-scattering capture rate and the thermalization time, we have applied our approach to the case of non-annihilating DM accumulated in a WD core throughout its lifetime. In doing so, we have found orders of magnitude corrections to the DM-nucleon cross sections for which the accumulated DM achieves self-gravitation.



§ ACKNOWLEDGEMENTS

We thank Maura Ramirez-Quezada for detailed discussions on an earlier version of this manuscript. 
NFB and GB were supported by the Australian Research Council through the ARC Centre of Excellence for Dark Matter Particle Physics, CE200100008.
SR acknowledges partial support from the UK STFC grant ST/T000759/1 and from the Fermi National Accelerator Laboratory (Fermilab), a U.S.
Department of Energy, Office of Science, HEP User Facility. 
MV  was supported by an Australian Government Research Training Program Scholarship. 
SR and MV  thank the  Institute for Nuclear Theory at the University of Washington for its hospitality and the Department of Energy for partial support during the completion of this work.







§ ANALYTIC DERIVATION OF THE THERMALIZATION TIME



In this section, we derive analytic expressions for the thermalization time for heavy DM. 
In section <ref>, we calculate the interaction rate for finite temperature in the high and low energy transfer regimes. In section <ref>, we estimate the corresponding average energy transfer, and the derivation of the thermalization time can be found in section <ref>. We discuss the effect considering nuclear form factors in the above mentioned calculation in section <ref>. 






 §.§ Interaction rate for thermalization





To avoid numerical precision issues, 
it is convenient to obtain an expansion of the interaction rate for low values of T_⋆. 
To this end, we revisit the expression for the differential interaction rate for finite temperature Eq. <ref>, 

    R_T^-(w→ v)    =   ∫_0^∞ ds  ∫_0^∞ dt    F(s,t)4μ_+^2/μn_T(r)v/wdσ_Tχ/dcosθ_cm(s.t,w,v) Θ(v - |t - s|), 
    
       F(s,t)    =   8μ_+^2/√(π)k^3tμ  e^-k^2u_T^2Θ(t+s-w).
 
 Next, we define  the following functions

    δ_EXP(x,x_0,c)    =     c   e^-c(x-x_0)Θ(x-x_0), 
    δ_G(x,x_0,c)    =   c/√(π)  e^-c^2(x-x_0)^2,
 
where x, x_0, and c are generic variables. In the limit  c→∞, these functions tend to delta functions, i.e.

    lim_c→∞∫_-∞^∞ dx δ_EXP(x,x_0,c) f(x)    → f(x_0),
    lim_c→∞∫_-∞^∞ dx δ_G(x,x_0,c) f(x)    → f(x_0),

where f is a generic function. 
Using the functions in Eqs. <ref> and <ref>, we rewrite F(s,t) 

    F(s,t) ds dt    =   δ_EXP(t^2,(w-s)^2,2μμ_+k^2) dt^2  δ_G(s,μ w/2μ_+,2μ_+k)ds.




To express the differential interaction rate in terms of the initial and final kinetic energy of the DM, we make the following substitutions

    w   =√(x/μ)v_T,     
        v=√(y/μ)v_T,
    
        s   =a/√(1+μ)v_T,     
        t=b/√(1+μ)v_T,

where v_T^2=3/m_T, x is the ratio of the DM initial kinetic  energy to T_⋆, which is O(1) close to thermalization, and y represents the final DM kinetic energy ratio. This leads to 

    R_T^-(x→ y) =   ∫_-∞^∞ da  ∫_-∞^∞   db    
       F(a,b) 2μ_+^2/μn_T(r)/√(x)v_T/√(μ) Θ(√(y)-√(μ)|a-b|/√(1+μ)) 
    
           ×   dσ_Tχ/dcosθ_cm(a/√(1+μ)v_T,b/√(1+μ)v_T,√(x)/√(μ)v_T,√(y)/√(μ)v_T),

where

    F(a,b) da db = δ_EXP[b^2,(√(1+μ)/√(μ)√(x)-a)^2,μ]db^2 δ_G(a,√(μ x)/√(1+μ),√(1+μ)) da.





In the large DM mass regime, the third parameter in  δ_G and δ_EXP in the previous expression is  large since μ≫1, 
hence  we can expand these functions around the following points

    a   =   √(μ x)/√(1+μ), 
    
        b    =   √(1+μ)/√(μ)√(x)-a = √(x)/√(μ(1+μ)).

Integrating over the delta functions, the resulting expression for the  interaction rate is 

    R_T^-(x→ y) =  2μ_+^2/μ^3/2n_T(r) v_T/√(x)dσ_Tχ/dcosθ_cm(√(x)/√(μ)v_T,√(y)/√(μ)v_T) Θ(y-μ_-^2/μ_+^2x).




We consider  differential cross sections proportional to powers of the DM-ion relative  velocity  v_rel^2m=w^2m and to powers of the momentum transfer q_tr^2m, namely

    dσ_Tχ/dcosθ_cm = σ_T2 w^2m = σ_T2(x v_T^2/μ)^m,     dσ_Tχ∝ v_rel^2m
    σ_T(m+1)2^m+1 q_0^2m( m_χ^2 w^2-v^2μ)^m =σ_T(m+1)2^m+1(m_T^2 v_T^2q_0^2)^m (x-y)^m,     dσ_Tχ∝ q_tr^2m

where σ_T and q_0 are normalization constants. 

Plugging the appropriate differential cross section in Eq. <ref> and integrating over y, we recover the interaction rate at high energy.
Thus, for dσ_Tχ∝ w^2n 

    Ω_T^-(x) = ∫_0^xdy   R_T^-(x→ y) =  n_T(r)σ_T w^2n+1,

and for  dσ_Tχ∝ q_tr^2m, we obtain

    Ω_T^-(x)  = n_T(r)σ_T w [2 μ^2 w^2  m_T^2/(μ+1)^2 q_0^2]^m≃  n_T(r)σ_T w (2  w^2  m_T^2/q_0^2)^m.


To obtain the previous results, we have made two assumptions. First, by approximating the  integral over dt to a delta function integral, we are assuming that the DM-target relative speed is equal, in magnitude, to the difference between the center of mass velocity and the DM speed. This approximation only holds when one of the two particles has negligible velocity compared to the other. 
Second,  by approximating  the integral over s to a theta function, we are effectively assuming that for μ≫1 the relative velocity is equal to the DM speed w, thus a negligible target speed is assumed. 
To calculate the interaction rate in the regime where the DM velocity is negligible compared to the target speed, we only make the first assumption and approximate δ_EXP to a delta function. Thus, 
we return to Eq. <ref>, and  only integrate over b,  so that the interaction rate becomes

    Ω_T^-(x)    =   2μ_+^2/μ^3/2n_T(r) v_T/√(x)∫_-∞^∞ da  ∫_0^x dy  δ_G(a,√(μ x)/√(1+μ),√(1+μ))
     Θ(√(y)-|2a√(μ)/√(1+μ)-√(x)|) 
    
           ×   dσ_Tχ/dcosθ_cm(a/√(1+μ)v_T,√(x(1+μ))-a√(μ)/√(μ(1+μ))v_T,√(x)/√(μ)v_T,√(y)/√(μ)v_T),

where we have swapped the integration order. The theta function sets a lower limit on the integration intervals for y and a, which leads to

    Ω_T^-(x)    =   2μ_+^2/μ^3/2n_T(r) v_T/√(x)∫_0^√(x(1+μ)/μ) da   δ_G(a,√(μ x)/√(1+μ),√(1+μ)) 
          ×∫_y_min(a)^x dy  dσ_Tχdcosθ_cm(a   v_T/√(1+μ),√(x)√(μ)v_T-a   v_T√(1+μ),√(x)√(μ)v_T,√(y)/√(μ)v_T)

where

    y_min(a) = (√(x)-2a√(μ)/√(1+μ))^2.
 
Next, we rescale a→ (z + √(μ x))/√(1 + μ) for δ_G to become a pure Gaussian.


As in the high energy transfer regime, we consider cross sections proportional to powers of the DM-target relative velocity which in the low energy regime read 

    dσ_Tχ/dcosθ_≃σ_T/2u_T^2m = σ_T/2v_T^2m z^2m,

and cross sections that depend on the momentum transfer, 

    dσ_Tχ/dcosθ_cm   =   σ_T(m+1)/2^m+1 q_0^2m (2 m_χ^2 t^2)^m [   1-(s^2+t^2-v^2)(s^2+t^2-w^2)/4 s^2 t^2]^m
       =   σ_T(m+1)/2^m+1(m_T v_T/q_0)^2m( √(μ x) - zμ/√(μ x)+z)^m (x-y)^m.


To estimate the interaction rate in this regime, we keep the dominant terms in μ,  for the case of cross sections proportional to v_rel^2m, this leads to

    Ω_T^-(x)    ∼     -2μ_+^2/μ^3/2n_T(r) v_T/√(x)∫_-√(μ x)^√(x/μ) dz δ_G(z,0,1) 4z √(x)/√(μ)σ_T/2(v_T^2z^2)^m
       ∼   1/√(π) n_T(r) σ_T v_T^2m+1∫_0^∞ dz   z^2m+1 e^-z^2  = m!/2√(π) n_T(r) σ_T v_T^2m+1.

Similarly, for cross sections that depend on the momentum transfer  q_tr^2m, we find 

    Ω_T^-(x)    ∼     -2μ_+^2/μ^3/2n_T(r) v_T/√(x)∫_-√(μ x)^√(x/μ) dz δ_G(z,0,1) σ(m_T/q_0)^2m2^m+1z^2m+1v_T^2m√(x)/√(μ)
       ∼   m!/2√(π) n_T(r)σ v_T (2m_T^2 v_T^2/q_0)^m.




 §.§ Average energy transfers




From the previous section,  the energy transfer in units of T_⋆ is given by x-y. Hence, the average energy lost per collision is given by 

    ⟨Δ x ⟩ = E(x)/Ω^-(w(x))=  ∫_0^xdy   R_T^-(x→ y)(x-y)/∫_0^xdy   R_T^-(x→ y).

As we shall see in the next section,  the computation of the thermalization time only requires evaluating the energy transfer E(x). This can be done in the same manner as for the interaction rate in section <ref>. 
In the high energy transfer regime x≫1,  we obtain 

    E(x) ≃ 2n_T(r)σ_T  (x /μ)^m+3/2 v_T^2m+1,       dσ_Tχ∝ v_rel^2m
    4(m+1)/m+2 n_T(r)σ_T  v_T  (x /μ)^m+3/2(2 m_T^2 v_T^2 /q_0^2)^m,       dσ_Tχ∝ q_tr^2m




In the low energy regime, i.e. x∼1 we find 
for cross sections proportional to powers of the DM-ion relative velocity v_rel^2m

    E(x) ∼Γ(m+3/2)n_T(r)σ_T /√(π)√(x/μ) v_T^2m+1,

while for differential cross sections proportional to q_tr^2m we have 

    E(x) ∼Γ(m+3/2)2(m+1)/m+2n_T(r)σ_T v_T/√(π)√(x/μ)(2m_T^2 v_T^2/q_0^2)^m.






 §.§ Thermalization times





Having  calculated the expression for the energy transfer in the high-energy and low-energy regimes, we can now evaluate  Eq. <ref>. 
The high energy term of E(x) primarily serves as a regulator, allowing us to set the initial temperature to infinity for convenience, which has a negligible impact on the overall result. For cross sections proportional to v_rel^2m, we obtain 

    E(x) ∼ n_T^c σ_T  v_T^2m+1√(x/μ)[2( x/μ)^m+1+1/√(π)Γ(3/2+m)   ].

Plugging this result in Eq. <ref>,  we calculate the thermalizaton time 

    =   ∫_1^∞dx/n_T^c σ_T  v_T^2m+1√(x/μ)[2( x/μ)^m+1+1/√(π)Γ(3/2+m)   ]
       ∼    C_m μ/n_T^cσ_T v_T 1/v_T^2m,

where C_m is a constant of order 𝒪(1) for m=0,1,2 

    C_m = π/2(m+1)(π/2m+2) [2√(π)/Γ(m+3/2)]^2m+1/2m+2.




For differential cross sections proportional to powers of the momentum transfer q_tr^2m, we obtain a similar result 

    ∼    C_m μ/n_T^cσ_T v_T(q_0^2/2 v_T^2 m_T^2)^m,
    
        C_m    =   π(m+2)/4(m+1)^2( π/2m+2) [2√(π)/Γ(m+3/2)]^2m+1/2m+2.



Note that the exact value of the C_m coefficients actually depends on all terms of the series expansion, while here, we have calculated only the dominant terms at high and low energy. Thus, we will consider them to be 𝒪(1) numbers. Here, we have neglected up-scattering, which we expect it to have also an 𝒪(1) effect. 




 §.§ Effect of the form factors




We can determine the effect of the nuclear form factors on the thermalization time by calculating the new value of E_FF(x) when accounting for them, and comparing it to the expression without them, i.e. 

    E_FF(x)/E(x)   =   ∫_0^x dy (x-y)R_T^-(x→ y) exp[-(x-y)/E_0]/∫_0^x dy (x-y)R_T^-(x→ y)
       ∼    1+𝒪(T_⋆/E_0),

where E_0 is at most  𝒪(), their exact value depending on the type of the interaction and the ion target, much greater than the WD core temperature (see Fig. <ref>). 


As stated in the previous section, the thermalization timescale only depends on the low energy behaviour of the energy transfer, so that when calculating thermalization times, it does not matter that the form factor suppresses the energy transfer at high energies, as this regime does not play any role close to thermalization.



















Bibliography 

 


