Thermal conduction by dark matter with velocity and momentum-dependent cross-sections

# Thermal conduction by dark matter with velocity and momentum-dependent cross-sections

## Abstract

We use the formalism of Gould and Raffelt GouldRaffelt90a () to compute the dimensionless thermal conduction coefficients for scattering of dark matter particles with standard model nucleons via cross-sections that depend on the relative velocity or momentum exchanged between particles. Motivated by models invoked to reconcile various recent results in direct detection, we explicitly compute the conduction coefficients and for cross-sections that go as , , , , and , where is the relative DM-nucleus velocity and is the momentum transferred in the collision. We find that a dependence can significantly enhance energy transport from the inner solar core to the outer core. The same can true for any -dependent coupling, if the dark matter mass lies within some specific range for each coupling. This effect can complement direct searches for dark matter; combining these results with state-of-the-art solar simulations should greatly increase sensitivity to certain DM models. It also seems possible that the so-called Solar Abundance Problem could be resolved by enhanced energy transport in the solar core due to such velocity- or momentum-dependent scatterings.

## I Introduction

If cosmological dark matter (DM) is similar in character to a WIMP (weakly-interacting massive particle), it may possess a weak-scale cross-section for interaction with ordinary nucleons. In this case, collisions between DM particles in the halo of the Milky Way and nuclei in the Sun will slow some of the DM enough to gravitationally bind it to the Sun Press85 (); Griest87 (); Gould87b (). If enough is captured and remains in the Sun, DM may impact the structure and evolution of the Sun itself Steigman78 (); Spergel85 (); Faulkner85 ().

The effects of captured DM particles on the Sun and other stars have been a subject of investigation for many decades (see e.g. Scott09 (); Turck12 (); Zurek13 () for reviews). Much work has concentrated on the detectable signals of annihilating DM in stellar cores, such as GeV-scale neutrinos escaping from the Sun Krauss86 (); Gaisser86 (); Griest87 (); Gandhi:1993ce (); Bottino:1994xp (); Bergstrom98b (); Barger02 (); Desai04 (); Desai08 (); IceCube09 (); IceCube09_KK (); IC40DM (); SuperK11 (); IC22Methods (); Silverwood12 (); IC79 () and the effects of core heating by annihilation products SalatiSilk89 (); BouquetSalati89a (); Moskalenko07 (); Spolyar08 (); Bertone07 (); Fairbairn08 (); Scott08a (); Iocco08a (); Iocco08b (); Scott09 (); Casanellas09 (); Ripamonti10 (); Zackrisson10a (); Zackrisson10b (); Scott11 (), or the MeV-scale neutrinos from these decay products Rott13 (); Bernal:2012qh (). The weakly-interacting nature of WIMP-like DM1 can also make it a medium for energy transport in the Sun or other stars Steigman78 (); Spergel85 (); Faulkner85 (), the implications of which have been studied extensively Gilliland86 (); Renzini87 (); Spergel88 (); Faulkner88 (); Bouquet89 (); BouquetSalati89b (); Salati90 (); Dearborn90b (); Dearborn90 (); GH90 (); CDalsgaard92 (); Faulkner93 (); Iocco12 (). Even when that energy transport is highly localised (i.e. conductive), WIMPs may cool the core by absorbing heat, which can then be deposited at larger radii. These changes can be potentially constrained with the temperature-sensitive B neutrino flux, or complementarily, with precision helio/asteroseismological measurements of the sound speed, small frequency separations and gravity modes Lopes02a (); Lopes02b (); Bottino02 (); Frandsen10 (); Taoso10 (); Cumberbatch10 (); Lopes10 (); Lopes:2012 (); Lopes:2013 (); Casanellas:2013 (); Lopes:2014 ().

So far all studies of energy transport by DM particles have assumed that the cross-section for scattering between WIMPs and nuclei is independent of their relative velocity and the momentum exchanged in the collision. This is also a common assumption in analyses of direct searches for dark matter. Among these, DAMA Bernabei08 (), CoGeNT CoGeNTAnnMod11 (), CRESST-II CRESST11 () and CDMS II Agnese:2013rvf () have reported excess events above their expected backgrounds, consistent with recoils caused by WIMP-like dark matter. On the other hand, XENON10 Angle:2011th (), XENON100 XENON2013 (), COUPP COUPP12 (), SIMPLE SIMPLE11 (), Super-Kamiokande SuperK11 () and IceCube IC79 (), under the assumption of constant scattering cross-sections and a standard halo model, would together exclude such “detections” from being observations of WIMPs. However, there are assumptions built into these exclusions. Some recent work (e.g. Bozorgnia:2013hsa (); Mao:2013nda ()) has focused on the astrophysical uncertainties. A greater deal of interest has developed in scattering cross-sections with a non-trivial dependence upon the relative velocity or momentum transfer of the collision, as it appears that such velocity– or momentum-dependent scattering could reconcile the various direct detection results Masso09 (); Feldstein:2009tr (); Chang:2009yt (); Feldstein10 (); Chang10 (); An10 (); Chang:2010en (); Barger11 (); Fitzpatrick10 (); Frandsen:2013cna ().2

The focus on constant cross-sections in stellar dark matter studies is therefore not because these are the only interesting DM models, but simply because the required theoretical background for including the effects of velocity or momentum-dependent couplings to nucleons in stellar structure simulations does not exist in the literature. This is in contrast to direct detection, where such effects can be easily computed (e.g. Cirelli13 ()). Here we provide such calculations for a range of non-constant cross-sections. This work will allow the full transport theory to be incorporated into a solar evolution code that also treats the capture and annihilation of DM, in order to derive self-consistent limits on DM models with non-standard cross-sections.

Our other motivation is the well-known Solar Abundance Problem: recent photospheric analyses APForbidO (); CtoO (); AspIV (); AspVI (); AGS05 (); ScottVII (); Melendez08 (); Scott09Ni (); AGSS () indicate that the solar metallicity is 10–20% lower than previously thought, putting predictions of solar evolution models computed with the improved photospheric abundances at stark odds with sound-speed profiles inferred from helioseismology Bahcall:2004yr (); Basu:2004zg (); Bahcall06 (); Yang07 (); Basu08 (); Serenelli:2009yc (). A deal of effort has been devoted to finding solutions to this quandary Bahcall05 (); Badnell05 (); Arnett05 (); Charbonnel05 (); Guzik05 (); Castro07 (); Christensen09 (); Guzik10 (); Serenelli11 (); Frandsen10 (); Taoso10 (); Cumberbatch10 (); Vincent12 (), but so far no proposal has really proven viable. Upcoming solar neutrino experiments such as SNO+, Hyper-Kamiokande and LENA should allow the solar core metallicity and temperature to be measured even more accurately than can be done with helioseismology alone Haxton13 (); Lopes13b (); Smirnov13 (), hopefully shedding some light on the discrepancy. One might postulate that the disagreement is in fact caused by enhanced energy transport in the solar core mediated by non-standard WIMP-nucleon couplings; the calculations we carry out here allow this proposition to be tested.

This paper is structured as follows: first, in Section II we provide some background on WIMP-nucleon couplings and how one goes about calculating energy transport by WIMPs, introducing the thermal conduction coefficients and . In Section III we briefly summarise the relevant equations from Ref. GouldRaffelt90a () that allow calculation of and . In Section IV we modify this treatment to account for non-trivial velocity and momentum structure in the scattering cross-section, and show the effect on and . In Section V we provide some examples of the effect on energy transport in the Sun, then finish with some discussion and concluding remarks in Section VI.

## Ii Preliminary theory

The quantity that describes microscopic interactions between WIMPs and Standard Model nuclei is the total cross-section for WIMP-nucleon scattering. Determining the impacts of energy conduction by a population of weakly-interacting particles in a star means working out their influence on bulk stellar properties like the radial temperature, density and pressure profiles. Initial solutions were obtained by way of simple analytic approximations Faulkner85 (); Spergel85 (); Gilliland86 (), but this quickly gave way to the established approach of solving the Boltzmann collision equation (BCE) inside the spherically symmetric potential well of a star Nauenberg87 (); GouldRaffelt90a (); GouldRaffelt90b (). The formalism and a general solution for standard WIMP-quark couplings were described by Gould and Raffelt GouldRaffelt90a (). There are three quantities that, together, are sufficient to describe the transport of energy by a diffuse, weakly-interacting massive particle in a dense medium such as a star:

• , defined as the inverse of the mean number of WIMP-nucleon interactions per unit length; this gives the typical distance travelled by a WIMP between scatterings;

• , related to the thermal diffusion coefficient; this parameterises the efficiency of a species’ diffusion inside the potential well;

• , the dimensionless thermal conduction coefficient; this describes the efficiency of energy transfer from layer to layer in the star.

Given a form of the WIMP-nucleon cross-section as a function of the microscopic kinetic variables (, etc.), the coefficients and are quantities that are averaged over the kinetic gas distributions, thus depending only on , the ratio between the WIMP mass and the mass of the nucleus with which it scatters. They can therefore be computed once and for all for each type of WIMP-nucleon coupling.

In the following sections, we compute the thermal conduction coefficients and of WIMPs with non-standard couplings to the standard model, going beyond the case computed explicitly by Ref. GouldRaffelt90a (). Because these computations must be done on a case by case basis, here we focus on the specific velocity- and momentum-dependent cross-sections and , with .3 Here is the WIMP-nucleus relative velocity, and is the momentum transferred during a collision. Our choice of these types of couplings is motivated in part by a number of recent direct-detection experiments Bernabei08 (); CoGeNTAnnMod11 (); CRESST11 (); Agnese:2013rvf (), but they can easily occur theoretically Fitzpatrick13 (); Kumar13 (). Mixed couplings (whether linear combinations of these couplings or cross-terms between them) also each require explicit calculation, and cannot be treated using combinations of the results we present here; we leave these for future work.

The dominant operators that arise in the most basic models of interactions between WIMPs and standard model particles (namely quarks, denoted ) are and . These respectively lead to the regular spin-independent and spin-dependent cross-sections, which have no velocity or momentum structure beyond standard kinematic factors, and have been the focus of exclusion efforts in direct experimental searches to date.

There is, however, no guarantee that the dominant coupling between the dark sector and the SM is as simple as or . Possibilities include a finite particle radius (the dark sector analogue of the nuclear form factor), a vector coupling to quarks, parity-violating interactions such as or or even a small dipole or anapole interaction with the standard model Pospelov00 (); Sigurdson04 (); Chang:2009yt (); Feldstein:2009tr (); Feldstein10 (); Chang:2010en (); Fitzpatrick10 (); Barger11 (); Frandsen:2013cna (); DelNobile:2013cva (). If these terms dominate, then the interaction cross-section will be proportional to some power of the momentum transfer :

 σ=σ0(qq0)2n, (1)

where is the cross-section at some normalisation momentum . Models charged under multiple gauge fields can display interactions that go as a sum of different powers of Feldstein:2009tr (). Even for scattering via standard or couplings, scattering between WIMPs and nuclei can be strongly momentum-dependent due to the nuclear form factor. This is especially important at high momentum transfer, so has a large impact on rates at which WIMPs can be scattered to below the local escape velocity and captured by a star. In this paper we neglect the influence of nuclear form factors, mainly because the majority of WIMP scattering in stellar interiors is with hydrogen and helium, at much lower momentum transfer than successful capture events. In the cases where momentum is exchanged with heavier nuclei (mainly C, N and O), the suppression due to a Helm form factor is only a few percent at relevant values of – much smaller, for example, than uncertainties on the the local interstellar WIMP density or velocity distribution. The form factors should in principle be taken into account if the impact of a given WIMP on a particular star is to be calculated in its fullest, goriest detail. Operators also exist which produce interactions that depend on the relative velocity between the WIMP and the nucleus:

 σ=σ0(vrelv0)2n. (2)

The collisional cross-section is proportional to the squared matrix element amplitude , which is typically a function of the centre-of-mass energy. This dependence can be expanded as Taylor series of , where the only sizeable contribution in the non-relativistic limit is the constant term. However, in some models the form of causes a partial or exact cancellation of the constant term, leaving . This is called p-wave suppression. The case in which the second term is also suppressed – leading to , is referred to as d-wave. On the other hand, DM models which scatter through the exchange of a massive force carrier exhibit a Sommerfeld-like enhancement. These can display “resonant” behaviour, which leads to a scattering cross-section proportional to Sommerfeld (); Hisano05 (); AHDM ().

Given a certain WIMP mass, a good choice for the reference velocity in a star is the temperature-dependent velocity

 vT(mχ,r)=[2kBT(r)/mχ]1/2. (3)

The typical thermal WIMP velocity can be obtained by multiplying by a further factor of GouldRaffelt90a (). Here refers to the stellar temperature profile as function of radius . In this paper we adopt  km s in the Sun, and provide a simple means for rescaling our results to other values of . Our chosen value is comparable to the typical relative velocities in direct detection experiments, where  km s is often employed. Similarly for , where our default is  MeV: this is about the typical momentum transferred both in thermal collisions in the Sun and direct detection. This value corresponds to nuclear recoil energies of  keV in direct detection experiments. As for , our results can be easily rescaled to any other ; the process is described in the passage following (74).

In this paper we will not focus on specific particle physics models or the WIMP-nucleon couplings they give rise to. Rather, we will take the general forms (1) and (2), and compute the effect on thermal conduction by an additional heavy species inside a star.

## Iii Thermal Conduction in a Star

### iii.1 Theory

The theoretical framework of energy transport by WIMPs was constructed in detail by Gould and Raffelt GouldRaffelt90a (). Transport by a weakly-coupled species can occur in two distinct regimes. If the typical inter-scattering distance of a particle is much larger than the typical geometric length scale (which in this case is the length scale of the WIMP distribution in a star, not the stellar radius), one reaches the Knudsen limit of non-local transport. In the opposite regime, the Knudsen number is , and WIMPs scatter many times per scale height as they travel through the medium. Here we are interested in the latter regime, where local thermal equilibrium (LTE) holds.

In GouldRaffelt90a (), Gould and Raffelt analysed the conduction of energy by a collection of massive particles in the LTE regime. They modelled energy transport using the Boltzmann collision equation:

 DF=l−1χCF. (4)

here is the phase space density for particle velocity and position at time . The typical inter-scattering distance for a WIMP is the inverse of the mean number of interactions it undergoes per unit length

 lχ(u,r)=[∑ini(r)⟨σi(u)⟩]−1, (5)

were the index denotes the individual nuclear species present in the stellar gas, and the angular brackets represent thermal averaging. is simply the number density of species at location . The relative weighting of different WIMP velocities in the thermal average depends on the local WIMP velocity distribution, which is itself a function of the temperature with height in the star ; the thermal averaging is therefore itself understood to be dependent on . The differential operator is , where is the gravitational acceleration.

To proceed, we make three key approximations:

1. The dilute gas approximation: WIMP-WIMP interactions can be neglected because they are much less frequent than WIMP-nucleon interactions.4

2. Local isotropy: this allows us to drop any dependence in the collision operator .

3. The conduction approximation: energy transport is conductive if is smaller than the other two length scales in the problem. That is,

• , i.e. local thermal equilibrium.

• , i.e. the typical inter-scattering distance is much smaller than the length scale over which the temperature changes. This means that , allowing a perturbative expansion in powers of .

Actually, approximation 3 is only strictly necessary for the computation of and rather than their application, as correction factors exist for making the conduction treatment approximately applicable even in the non-local regime when (see e.g. Bottino02 (); Scott09 ()). In practice both conditions in approximation 3. break down for cross-sections  cm and smaller – we will discuss the correction factors explicitly in Section V. The collision operator has a straightforward intuitive definition: is a function of the WIMP velocity , position and time , and represents the local rate of change of the WIMP phase space distribution. It can be formally written in terms of the components and ,

 CF = ∫d3vC′(u,v,r,t)F(v,r,t),with (6) C′ ≡ Cin(u,v,r,t)−Cout(v,r,t)δ3(v−u). (7)

is the rate at which particles with velocity are scattered to velocity . is the rate at which particles with velocity are scattered to any other velocity. Performing the integral of in (6) over incoming velocities gives the total rate of scattering to velocity from any velocity, whereas the integral over gives the total rate of scattering from velocity to any other velocity. The difference in these two rates is therefore the net rate of change in the WIMP phase space distribution.

The BCE was solved by GouldRaffelt90a () perturbatively, with the expansion

 F=F0+εF1+... (8)

where is an expansion parameter (as explained above),5 and is a regular Boltzmann distribution. Such a distribution is in thermal equilibrium by definition, so . Using this expansion to first order, dropping the term because it is much smaller than , the first order BCE is

 DF0=εlχCF1. (9)

Assuming that the thermal timescale in a star is sufficiently long compared to the timescale for conductive energy transport by multiple WIMP scatterings, the solution of the BCE will be time-independent (stationary), such that the time-derivative in can be ignored. In this case, the left-hand side of (9) is

 DF0(u,r)=εlχ[α(r)+mχu22T(r)∇lnT(r)|∇lnT(r)|]uF0(u,r), (10)

where is a separation constant. Assuming the stellar temperature gradient to be spherically symmetric, both and are pure dipoles. We use the standard notation for spherical harmonics: , where corresponds to the degree (related to total angular momentum) and is the order (angular momentum solely in the direction). We can therefore express the BCE entirely in terms of monopole () and dipole ( components of the various quantities; this is the Rosseland approximation.

It is convenient to transform our working variables to dimensionless quantities. We define the WIMP-to-nucleus mass ratio

 μ≡mχ/mnuc, (11)

and divide velocities by the local temperature-dependent velocity (3)

 x ≡ v/vT, (12) y ≡ u/vT, (13) z ≡ vnuc/vT, (14)

with denoting the velocity of a nucleus.

We also define a dimensionless (angle-dependent, differential) cross-section , and a total dimensionless (angle-independent, non-differential) cross-section . The total dimensionless cross-section is simply integrated over the centre-of-mass scattering angle :

 ^σtot(vrel)≡∫1−1dcosθCM^σ(vrel,q), (15)

remembering that is itself a function of both and . The original dimensionless differential cross-section () is defined by the requirement that the total dimensionless cross-section be unity at , i.e.

 ^σtot(vT)=1. (16)

Similarly, we define the normalised phase space distribution functions

 fj,mν(x,r)dx≡1nχ,LTE(r)Fj,mν(v,r)dv, (17)

where is the WIMP number density in the LTE (conductive) approximation, is the expansion order in , is the degree and is the order of the spherical harmonic expansion.

The first order equation (9) can then be expressed as the combination of

 (αmy−δm0y3)f0,00(y,r)=∫dxC(y,x,r)f1,m1(x,r) (18)

and the stationarity condition, expressed as the absence of a net WIMP flux across any given surface inside the star:

 ∫dxxf1,m1(x,r)=0. (19)

Here are thermal diffusivity coefficients corresponding to the dipole coefficients in the spherical harmonic expansion of . In the Dirac notation of GouldRaffelt90a (), where

 |f⟩ ≡ f(y), (20) Q|f⟩ ≡ ∫dxQ(y,x)f(x), (21) ⟨g|f⟩ ≡ ∫dxg(x)f(x), (22) ⟨g|Q|f⟩ ≡ ∫dydxg(y)Q(y,x)f(x), (23)

we see that (18) and (19) become quite compact:

 αm(r)|yf0,00⟩−δm0|y3f0,00⟩ = C|f1,m1⟩, (24) ⟨y|f1,m1⟩ = 0. (25)

In terms of the inverse of the collisional operator , we can write the solution to the first-order BCE as

 |f1,m1⟩=αmC−1|yf0,00⟩−δm0C−1|y3f0,00⟩. (26)

The diffusivity coefficients are then fully specified by multiplying (26) by :

 α±1 = 0, (27) α0 = ⟨y|C−1|y3f0,00⟩⟨y|C−1|yf0,00⟩. (28)

This can be used to obtain the stationary WIMP density profile GouldRaffelt90a ()

 nχ,LTE(r) = nχ,LTE(0)[T(r)T(0)]3/2 × exp⎡⎢ ⎢⎣−∫r0dr′kBα(r′)dT(r′)dr′+mχdϕ(r′)dr′kBT(r′)⎤⎥ ⎥⎦,

where refers to the gravitational potential at height in the star. Together with the thermal conductivity

 κ=√23⟨y3|f1,01⟩, (30)

one can then use (III.1) to finally obtain the luminosity carried by WIMP scattering,

 Lχ,LTE(r) = 4πr2ζ2n(r)κ(r)nχ,LTE(r)lχ(r) (31) × [kBT(r)mχ]1/2kBdT(r)dr,

and the corresponding energy injection rate per unit mass of stellar material:

 ϵχ,LTE(r)=14πr2ρ(r)dLχ,LTE(r)dr. (32)

Here is the stellar density, for velocity-dependent scattering and for momentum-dependent scattering. The factors of in (31) connect the velocity (momentum) scale at which the reference cross-section is defined, to the thermal scale () at which the dimensionless conductivity is computed. Note that as per Ref. Scott09 (), our sign convention is such that positive and refer to energy injection rather than evacuation, which differs from Ref. GouldRaffelt90a ().

The problem of calculating the rate of energy transport by WIMPs is now one of explicitly calculating and for the WIMP-nucleus mass ratio and cross-section of interest. The only position dependence in the dimensionless BCE arises through the temperature; the transformation to dimensionless variables , etc., means that the collision operator , as well as and , contain no -dependence at all. This allows and to be expressed solely in terms of the mass ratio . In a real star, more than one species of gas is present however, in different ratios at different . Thus, to compute (III.1) and (31) and must be computed for the specific mixture of gases at each height Scott09 ():

 α(r,t)=∑iσini(r,t)∑jσjnj(r,t)αi(μi) (33)

and

 κ(r,t)=[lχ(r,t)∑i{κi(μi)li(r,t)}−1]−1. (34)

Here is the number density of species , while and are respectively the scattering cross-section and the typical interscattering distance of WIMPs with species ; this is in contrast with the complete sum in (5). Henceforth we denote and .

Gould and Raffelt GouldRaffelt90a () calculated these quantities for a range of mass ratios, but only the case of velocity- and momentum-independent cross-sections; here we extend the treatment to more general cross-sections.

### iii.2 Calculations of α and κ

Because we plan on evaluating the collision operator numerically, it makes sense to express the functions as one-dimensional “vectors” and the operators as two-dimensional matrices. This allows direct calculation of a well-defined inverse operator , simply by inverting the corresponding matrix for . Although the integrals of the form are from zero to infinity, the exponential behaviour of means that the contribution to each integral is less than a fraction of a percent. We therefore truncate the vectors and matrices at .

is relatively straightforward to construct:

 Cout(x,r)=∫d3z|x−z|^σtot(vT|x−z|)Fnuc(z), (35)

where is the velocity distribution of the nuclei,

 Fnuc(z)=(πμ)−3/2e−|z|2/μ. (36)

An expression for is not as obvious. We transform to centre of mass (CM) coordinates, where is the velocity of the CM frame, and and are respectively the incoming and outgoing WIMP velocities in that frame.6 These quantities are illustrated in Figure 1.

Ref. GouldRaffelt90a () then finds

 Cjin(y,x,r) = (1+μ)4yx∫∞0da∫∞0dbFnuc(z)2πb⟨Pj^σ⟩ (37) × Θ(y−|a−b|)Θ(a+b−y) × Θ(x−|a−b|)Θ(a+b−x).

is the angle-averaged product of the differential cross-section with the -th Legendre polynomial, as expanded around transverse scattering angles in the lab frame. The angle averaging is performed in the azimuthal direction around the axis in the CM frame. The integral over the zenith angle is performed implicitly when integrating over and . Note that in this angle-average, is a function of , whereas is a function of and . The two angles are related as

 cosθCM = A+Bcosϕab′, (38) cosθlab = G+Bb2xycosϕab′, (39)

with

 A ≡ cosθab′cosθab = (x2−a2−b2)(y2−a2−b2)4a2b2, B ≡ sinθab′sinθab, (42) G ≡ (x2+a2−b2)(y2+a2−b2)4a2xy. (43)

As we are only solving the BCE to first order, the only relevant term in (37) is . Noting that , we then have

 ⟨Pj^σ⟩=⟨P1^σ⟩ = 12π∫2π0dϕab′(G+Bb2xycosϕab′) (44) × ^σ[(1+μ)bvT,A+Bcosϕab′].

## Iv conduction coefficients for non-standard WIMPs

We write the cross-sections of interest to us explicitly in terms of the dimensionless variables. We will find will usually have a tractable analytic form, but that the evaluation of the kinematic integrals for (37) must be done numerically, which we do using the CUBPACK multi-dimensional adaptive cubature integration package CUBPACK (). In every case, we produce an linearly-spaced matrix for , which we explicitly invert to find , and . We find to be sufficient: both and vary by less than one part in when going from to , across the full range of values we consider. For the standard case, our method of explicit matrix inversion yields identical results to the values of and presented in Table I of Ref. GouldRaffelt90a (), where a slightly different iterative method was used. These values are shown in the figures in this section and are given explicitly in Table 3 at the end of the text. Our values of agree exactly within given precision with Gould and Raffet for small , and to within a part in as . Our results agree within at small , up to as approaches 100. These differences are small enough not to have any appreciable impact on modelled energy transport. We have furthermore cross-checked our results for the collision operator and as well as the first-order solution presented in the figures of Ref. GouldRaffelt90a ().

### iv.1 Velocity-dependent scattering

The relative velocity is the difference between the incoming WIMP speed and the nucleus speed ; expressed in terms of these velocities, the dimensionless differential cross-section is just

 ^σv2n=12|x−z|2n. (45)

is then modified by extra powers of the relative velocity. The angular and radial integrals can be performed analytically, and the results are given in Table 1. is also modified in a simple way. The angular dependence drops out in the centre of mass frame where we compute and , so that

 ^σv2n=12(1+μ)2nb2n. (46)

These are related to (2) via:

 σ=2σ0ζ−2n^σv2n, (47)

where is defined below (32). There is one final complication in the computation of (37) for the case: the expression (46) diverges for , complicating the numerical integration of , even though the integral itself is finite. We address this by imposing a cutoff velocity such that:

 ^σv−2→1+ω22(|x−z|2+ω2)=1+ω22(1+μ)2b2+2ω2, (48)

where the factor in the numerator ensures that . Both and converge to finite values as . This convergence is illustrated in Figure 4. As is varied from to , varies by less than one part in , whereas changes by less than at low , up to for . The thermally-averaged cross-section, necessary for computing a particle’s typical inter-scattering distance , can potentially depend on non-zero . A large value is reasonable, for example, in the case of resonant Sommerfeld-enhanced scattering, where the cross-section can reach a saturation value for , below which becomes constant.

Values of for are illustrated on the left-hand side of Figure 2. For positive powers of , this results in enhanced scattering of particles with high incoming velocity ; conversely this enhancement is present a low incoming velocities for .

In Figure 3, the thermal conduction coefficients and are shown for scattering (blue dashed lines), and for (cyan dotted lines). For the case, and converge to well-defined curve as ; these are the values we show in the bottom panels of Figure 3 (dark red dashed lines), and later use to compute . Numerical values are provided in Table 3.

The effect of positive powers of is to suppress at low values of , and to suppress for large values of . is a measure of how well WIMPs can diffuse outward in the potential well of a star. This suppression ultimately means that the distribution will be more compact than in the standard case. The effect of a suppression, meanwhile, implies that the WIMP population will be less efficient at heat transport via scattering with light nuclei such as hydrogen. In the lower panel of Figure 3 we show that the opposite behaviour is true for : a “fluffier” WIMP distribution around the stellar core, with enhanced heat transport via collisions.

In the case of a large cutoff velocity , the thermal conduction coefficients lie in between the and cases. We illustrate this in Figure 4, for several values of the dimensionless cutoff velocity . Sommerfeld-enhanced DM models, therefore, would benefit from some of the enhanced conduction of the case.

### iv.2 Momentum-dependent scattering

The momentum transferred during a collision, is equal to the total momentum gained or lost by the WIMP:

 q2 = m2χ|u−v|2 (49) = q20ζ−2|x−y|2 = 2b2ζ−2q20(1−cosθCM).

The latter expression allows us to write the cross-section in terms of the centre-of-mass velocity of the incoming WIMP, , and the CM scattering angle , as the scattering event does not change the WIMP’s speed in the CM frame, only its direction.

Then

 ^σq2n=2−nHn|x−y|2n=Hnb2n(1−cosθCM)n, (50)

where the normalisation factor required to ensure that is

 Hn=(1+μ)2n∫1−1(1−cosθ′)ndcosθ′. (51)

The structure of (50) shows us that

 ^σq2n=2^σv2n(1−cosθCM)n∫1−1(1−cosθ′)ndcosθ′, (52)

so

 ^σtot,q2n=^σtot,v2n=(1+μ)2nb2n, (53)

as we would expect, given that is independent of . This also means that . In terms of the dimensionful cross-section (1),

 σ=2nHnσ0ζ−2n^σq2n. (54)

The component requires an explicit evaluation of :

 ⟨P1^σ0⟩ = 12G, (55) ⟨P1^σq2⟩ = 12b2(1+μ)2[G(1−A)−b2B22xy], (56) ⟨P1^σq4⟩ = 38b4(1+μ)4[{G(1−A)−b2B2xy}(1−A) (57) + GB22], ⟨P1^σq−2⟩ = 1ln2−lnξb−2(1+μ)−2 (58) × [(1−A+ξ)b2+Gxyxy√(1−A+ξ)2−B2−b2xy],

where , and are defined in (III.243). The parameter comes from the replacement . We make this substitution because for , the factor in (50) causes the cross-section to diverge for forward scattering at small angles; likewise for the integral (15). For the computation of in (28), the powers of allow the divergence to cancel. However, , which governs momentum transfer, formally diverges. Fortunately, the divergence at corresponds to forward scattering, in which no momentum is actually transferred. We regulate this divergence with the “momentum transfer cross-section,” which comes from plasma physics Krstic () but has also been applied to parameterise transport by WIMPs Tulin:2013teo ():

 ^σtot,T(vrel) ≡∫1−1dcosθCM^σT(vrel,q), (59) ^σT(vrel,q) ≡(1−cosθCM)^σ(vrel,q), (60)

where is as per (50), but with the replacement in the denominator of (51). This gives

 ⟨P1^σT,q−2⟩=12Gb2(1+μ)2, (61)

leading to a finite, well-defined value of .

The behaviour of as a function of and for momentum-dependent cross-sections is shown in the right-hand panels of figure 2. The two upper right-hand panels show an enhancement in the scattering rate when the momentum transfer is large, over both the constant and velocity-dependent cases, which favour collisions closer to the line. The lower-right panel shows the case computed with (58) and , illustrating the divergent behaviour in the forward-scattering rate, and a general suppression of collisions which lead to appreciable momentum exchange.

The values of and are shown in dot-dashed orange for the case and in dot-dashed green for in Figure 3. The case is also shown (dashed purple lines). The behaviour of and for a cross-section is qualitatively similar to the case. Indeed, is identical for the and cases, thanks to the momentum-transfer cross-section, which essentially deletes the angular dependence in the case and makes identical for momentum and velocity-dependent cross-sections when . In the next section we will show that, once this behaviour is combined with the different typical inter-scattering distances in models with momentum-dependent and velocity-dependent scattering, the conduction of energy can be greatly modified with respect to the fiducial constant cross-section case.

## V Effect on energy transport in the Sun

We now illustrate the effect of a non-standard WIMP-nucleon cross-section on the conduction of energy within the Sun. Our discussion applies equally well to any star, but modelling and measurements of the Sun and its properties are far more precise than for other stars. We use the temperature, density and elemental composition profiles computed with the AGSS09ph solar model7 AGSS (); Serenelli:2009yc ().

Scattering cross-sections  cm allowed by bounds from direct detection experiments can easily lead to energy transport outside the LTE regime, i.e. non-local transport. To account for this, based on the results of Gould & Raffelt GouldRaffelt90a (); GouldRaffelt90b (), Refs. Bottino02 (); Scott09 () introduce the quantities:

 h(r) = (r−rχrχ)3+1, (62) f(K) = 11+(KK0)1/τ, (63)

where is the Knudsen number, and . The WIMP scale height as a function of the central temperature and density is

 rχ=(3kBTc2πGρcmχ)1/2. (64)

The WIMP distribution becomes a combination of the isothermal and LTE distributions:

 nχ(r)=f(K)nχ,LTE+[1−f(K)]nχ,iso, (65)

where

 nχ,iso(r,t)=N(t)e−r2r2χπ3/2r3χ. (66)

The total luminosity can be finally written:

 Lχ,total(r,t)=f(K)h(r,t)Lχ,LTE(r,t). (67)

The expressions (62, 63) are fits to the results of direct Monte Carlo solutions to the Boltzmann Equation for in (1, 2). In principle similar solutions should be obtained and fitted for , but carrying out such Monte Carlo simulations is well beyond the scope of this paper. We expect that the form of the suppression functions will be broadly similar to (62, 63) for all , however. Conductive transport is generically expected to be most efficient for some close to but less than 1, because the larger the inter-scattering distance of a WIMP relative to the scale height , the further it will have travelled (in an effective sense) in the star, and therefore the greater impact it will have on the local luminosity function, when it deposits its energy – but beyond , the conduction approximation breaks down, and transport must be less efficient. Transport is also enhanced in the surface regions of a star due to the enhancement of (III.1) at large compared to (66), and suppressed in the stellar core generically, as the radial height is actually less than , causing the impacts of conductive transport to be predominantly felt ‘downstream’ (i.e. at radii larger than ; see GouldRaffelt90a () for details).

For the calculation of WIMP conductive luminosities, we require the typical inter-scattering distance (5), which depends on the local velocity-averaged scattering cross-section. Taking a Maxwell-Boltzmann distribution for the velocities of the WIMPs and nuclei once more, the distribution of their relative velocities, in normalised units, is

 Fχ−n(x−z)=[π(1+μ)]−3/2e−|x−z|21+μ. (68)

This then gives:

 ⟨σv2⟩ =σ0ζ−2⟨v2relv2T⟩ =32σ0ζ−2(1+μ), (69) ⟨σv4⟩ =σ0ζ−4⟨v4relv4T⟩ =154σ0ζ−4(1+μ)2, (70) ⟨σv−2⟩ =σ0ζ2⟨v2Tv2rel⟩ =2σ0ζ2(1+μ)−1, (71)

where . The averaged -dependent cross-sections follow from these, as