Generalised form factor dark matter in the Sun

# Generalised form factor dark matter in the Sun

[    [    [
###### Abstract

We study the effects of energy transport in the Sun by asymmetric dark matter with momentum and velocity-dependent interactions, with an eye to solving the decade-old Solar Abundance Problem. We study effective theories where the dark matter-nucleon scattering cross-section goes as and with or , where is the dark matter-nucleon relative velocity and is the momentum exchanged in the collision. Such cross-sections can arise generically as leading terms from the most basic nonstandard DM-quark operators. We employ a high-precision solar simulation code to study the impact on solar neutrino rates, the sound speed profile, convective zone depth, surface helium abundance and small frequency separations. We find that the majority of models that improve agreement with the observed sound speed profile and depth of the convection zone also reduce neutrino fluxes beyond the level that can be reasonably accommodated by measurement and theory errors. However, a few specific points in parameter space yield a significant overall improvement. A 3–5 GeV DM particle with is particularly appealing, yielding more than a improvement with respect to standard solar models, while being allowed by direct detection and collider limits. We provide full analytical capture expressions for - and -dependent scattering, as well as complete likelihood tables for all models.

1]Aaron C. Vincent,

2]Aldo Serenelli

3]and Pat Scott

Prepared for submission to JCAP

Generalised form factor dark matter in the Sun

• Institute for Particle Physics Phenomenology (IPPP), Department of Physics, Durham University, Durham DH1 3LE, UK

• Institut de Ciències de l’Espai (ICE-CSIC/IEEC), Campus UAB, Carrer de Can Magrans s/n, 08193 Cerdanyola del Valls, Spain

• Department of Physics, Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK

## 1 Introduction

A non-relativistic relic dark matter particle, with a mass of a few GeV or more, is the leading candidate to explain astrophysical and cosmological phenomena ranging from cluster kinematics and galactic rotation curves, to gravitational lensing and the heights and positions of the acoustic peaks in the cosmic microwave background. Dark matter (DM) may have been produced in a similar way to Standard Model (SM) particles, either via chemical freeze-out (as in the weakly-interacting massive particle – WIMP – scenario) or via an initial asymmetry, analogous to baryogenesis (as in the asymmetric DM – ADM – scenario). If either of these scenarios is correct, it is possible that DM interacts weakly with SM particles. Such interactions would be seen most easily as a small elastic scattering cross-section between DM and quarks.

The search for DM-quark interactions has been the focus of terrestrial direct detection experiments such as DAMA [1], CoGeNT [2], CRESST-II [3] and CDMS II [4], who have all reported excess events above their expected backgrounds. On the other hand, XENON10 [5], XENON100 [6], COUPP [7], SIMPLE [8], LUX [9] and SuperCDMS [10] have all established strong limits on the DM-nucleon cross-section, seemingly in contradiction to the excesses observed by other experiments. These underground detectors are typically most sensitive to DM particles with masses of order 50 to 100 GeV, which lead to the largest recoil energies on the heavy nuclei used as targets. For the same reason, they are also best suited to probing fast-moving DM particles, resulting in a threshold velocity of a few tens of km s for an incoming DM particle to create a nuclear recoil event.

Collisions between DM and nuclei can also lead to capture and accumulation of DM in the solar core. For this to occur, collisions between DM and nuclei in the Sun must result in sufficient energy transfer for the DM velocity to be brought below the local escape velocity. The kinematic region probed by the Sun is quite different to the one probed by direct detection: optimal energy transfer, leading to optimal capture rates, occurs for DM particle masses closely matching the solar composition, i.e. a few GeV. Because DM gains speed as it falls into the solar potential well, the low-velocity tail of the local DM distribution is the dominant contributor to solar capture; the opposite is true for direct detection experiments. Direct detection and solar physics are therefore highly complementary laboratories for the study of DM-quark scattering.

### 1.1 Dark matter in the Sun

The effect of DM on stars has been the subject of investigation for some time (see [11, 12, 13] for reviews). There are two main scenarios where DM might create observable effects on the Sun: 1) the capture and annihilation of WIMP-like particles; and 2) the accumulation of an “asymmetric” species. The first case is characterised mainly by a search for high-energy neutrinos with [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31] and MeV-scale decay products [32, 33], using detectors such as IceCube and SuperKamiokande. Annihilation of DM in a star can also provide an energy source through the release of other SM particles [34, 35, 36, 37, 38, 39, 40, 41, 42, 11, 43, 44, 45, 46, 47], leading to changes in the core structure. The capture rate of DM in the Sun is however far too low for the energy released this way to have any significant effect on solar structure [11].

Our focus in this paper is therefore on the case where DM either cannot annihilate, or does so far more slowly than it is captured. In this case DM is effectively asymmetric (i.e. ADM, [48, 49]), and large quantities of DM may have built up inside the Sun over its lifetime. The weak interactions of DM with quarks give the DM particles relatively large inter-scattering distances, making them potentially significant conductors of energy. Just like the capture process, energy transfer is most efficient for lighter DM masses ( GeV), as momentum transfer is maximised when the masses of the colliding particles are equal, and the Sun is mostly H () and He (). It is interesting to note that this is roughly the mass range expected in the most common models of ADM [50], in order to explain the 1:5 cosmological relic baryon-to-DM density ratio. In contrast, earth-based direct detection experiments lose sensitivity in this range, as they make use of high-mass elements such as germanium () and xenon (), for which high incoming DM velocities are necessary to create a measurable recoil.

There is a window of elastic DM-nucleon scattering cross-sections [51, 52, 53, 54, 55] at cm for spin-independent (SI) couplings ( cm in the spin-dependent – SD – case) for which is large enough to allow sizeable capture of DM, but small enough that the mean interscattering distance is still large. A large inter-scattering distance allows heat to be redistributed away from the solar core by DM-nucleon scattering. This has the effect of reducing the temperature of the solar core , and increasing the central density and pressure. Changes of state variables near the solar centre affect the production of neutrinos from fusion processes, especially the flux of B neutrinos, which goes as with 20–25. The temperature and density changes in the core reduce the local sound speed, and force an overall mass redistribution that impacts the solar structure at other radii. This leads to modifications of the sound speed profile over the entire depth profile of the Sun, and shifts the height of the base of the convection zone. Both the sound speed profile and the depth of the convection zone have been independently measured using helioseismology.

It has been shown that high-precision solar evolution models including capture, transport and (minimal) annihilation of DM can be built to satisfy the observed solar age, radius, and luminosity [54, 53, 55]. At the same time, it appears possible for the inclusion of DM in such models to affect the (less well constrained) B flux in an observable way, and even improve agreement with the observed sound speed profile and the depth of the convection zone. Neither of these latter two observables are well reproduced by standard solar models computed with the latest surface compositions [56, 57, 58, 59, 60, 61]. This issue, known as the “solar composition” or “solar abundance” problem, has been brought about by the 20–25% reduction in the measured solar metallicity in recent years [62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73], and is one of our motivations for this paper.

Although several studies have indicated that ADM can alleviate some of this tension, the cross-sections required to do so are typically far higher than allowed by limits from direct detection. Here we investigate whether broader consideration of the kinematic structure of the DM-nucleon vertex might provide a way around this. In the process, we provide first rigorous limits from solar physics on such DM models, which we refer to as ‘generalised form factor dark matter’. In a separate paper [74] we discussed a specific realisation of momentum-dependent dark matter that leads to a 6 improvement over the Standard Solar Model (SSM). We revisit this model in Sec. 6.1.

### 1.2 Generalised form factor dark matter

The kinematic differences between direct detection and the Sun become even more marked if the DM-nucleon interaction is not assumed to be independent of the DM-nucleon relative velocity, or the momentum transferred in the collision. There is indeed no guarantee that the standard SI and SD operators correctly represent the DM-quark interaction. In particle physics, the interaction cross-section generally depends on the centre of mass energy and the transferred momentum, parameterised using the Lorentz-invariant Mandelstam variables , and . In the non-relativistic limit, these become the centre of mass momentum, proportional to the relative velocity , and the transferred momentum . As these are small quantities, the constant term usually dominates in a series expansion of the cross-section. However, many models with non-trivial dependencies on and exist, typically motivated by theoretical arguments or attempts to reconcile experimental results.

To make quantitative predictions, a specific form of must be chosen. Because we wish to remain as general as possible, we choose to focus on couplings of the form and , with . The and forms are respectively called -wave and -wave interactions, and correspond to the cases where the initial state particles possess 1 and 2 units of relative angular momentum, respectively. These are always present, but normally only dominate when all lower-order terms in the scattering matrix element – including the constant (-wave) term – are suppressed due to cancellations. Cross-sections depending on can arise, for example, from a non-zero particle radius (the analogue of a nuclear form factor), from parity-violating couplings like , and , or from a small anapole or dipole interaction between the dark and visible sectors [75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87]. We refer to the class of models where DM-nucleon scattering cross-sections depend on some combination of and/or as ‘generalised form factor DM’ because it generalises the effects of form factors to arbitrary powers of and .

Concretely, we focus on the couplings

 σ=σ0(qq0)2n (1.1a) σ=σ0(vrelv0)2n. (1.1b)

These can lead to either spin-dependent (SD) or spin-independent (SI) interactions, depending on the axial structure of the DM-nucleon interaction vertex. The normalisation must be defined with respect to some reference velocity or momentum . We will choose  km s, the typical halo DM velocity, and MeV, corresponding to a nuclear recoil energy of around 10 keV in an underground direct detection experiment.

The DM-nucleus cross-section is related to the above DM-nucleon cross-sections via:

 σN,i=m2nuc(mχ+mp)2m2p(mχ+mnuc)2[σSIA2i+σSD4(Ji+1)3Ji|⟨Sp,i⟩+⟨Sn,i⟩|2], (1.2)

where and are respectively the atomic number and total angular momentum of nuclear species ; and are the spin expectation values of its proton and neutron systems. Given the dependence of Eq. 1.2 we will find that, in spite of the Sun’s small metallicity, spin-independent DM can have a significantly larger effect than a spin-dependent DM candidate, which couples mostly to hydrogen.111See also caveats to this treatment in Sec. 2.2.1. Beyond a few studies, e.g. [88], this fact has not been emphasised very much in the literature. The full effect of a momentum-dependent cross-section will furthermore depend crucially on the composition of heavier elements.

The full impact of momentum and velocity-dependent DM on solar observables has not been studied before. The authors of Ref. [30] computed the effect of such couplings on the capture rate of spin-dependent DM, and computed neutrino fluxes from an annihilating species. However, they did not include a treatment of heavier elements, nor of the crucial energy transport by DM once captured. In Refs. [89, 90], the authors computed capture and transport rates for ADM models with long-range interactions. To account for the impacts of the non-trivial scaling of these cross-sections with momentum and velocity, they employed effective cross-section scaling factors to decouple the cross-sections entering the capture and transport calculations, and avoid modifying the standard velocity-and-momentum-independent treatment of capture and energy transport. As we show later, it happens that this rescaling can indeed be done without any loss of generality for capture in the velocity-dependent case, but it is not possible in the momentum-dependent case. It is also not possible to account for the effects on energy transport of either a velocity or momentum dependence in this manner; rather, a full recalculation of the transport coefficients must be performed [91].

The structure of this paper is as follows: in Section 2, we review the capture equations for DM in the Sun, and present the necessary modifications for velocity and momentum-dependent scattering of DM with nucleons. In Section 3 we review the theory of conductive heat transport by DM developed in Refs. [92, 91], along with its application to solar modelling. Section 4 describes the DarkStec computer code that we have developed for simulating the effects of generalised form factor DM on the Sun. We present results in Section 5, and discuss their implications with regards to the Solar Abundance Problem, current experimental limits and previous work on the topic in Section 6. We summarise in Section 7.

## 2 Capture of dark matter by the Sun

### 2.1 Standard (velocity and momentum independent) treatment

The population of DM particles in the Sun is follows the differential equation

 dNχ(t)dt=C⊙(t)−A(t)−E(t), (2.1)

where is the capture rate, is rate at which annihilations occur and represents evaporation. Unless DM is strongly self-interacting (not the case we consider here, but discussed in Ref. [93]), does not depend on the DM population already captured by the Sun. is the rate of annihilation, and is proportional to the square of the DM population. Here we consider the case where DM is fully asymmetric, so , although we comment briefly in Section 4.1 on the implications of allowing DM to self-annihilate.

Evaporation occurs when a DM particle gains enough energy from a scattering event to overcome the Sun’s gravitational potential and escape, so is linear in the DM population (being simply the product of the single-particle evaporation probability and the number of candidates for evaporation). Evaporation requires a significant gain in momentum, as the typical velocity of a thermalised DM particle is  km s, whereas the escape velocity in the solar core approaches 1400 km s. This means that evaporation is only significant if DM is similar in mass to the nucleus with which it scatters in the Sun. In practice this means that DM about the mass of helium (4 GeV) and lighter is typically most prone to evaporation, but other masses closely matched to significant elements in the Sun (C, N, O, Fe , 14, 16, 56 GeV) can also in principle be affected [94]. The evaporation rate depends on the interaction cross-section, mean free path and thermal regime (LTE vs Knudsen); extending the standard analyses to generalised form factor DM is therefore non-trivial. In practice the rate of evaporation is extremely low in almost all cases where the nuclear scattering cross-section is allowed by direct detection; for very specific analyses it should be taken into account, but for the purposes of this paper we assume that the evaporation rate is zero. We intend to return to this point in detail in a future paper.

We now turn to the capture rate as it is implemented in our simulations. As we are using the DarkStars code [95], we closely follow Refs. [11, 96]. We take the local distribution function of DM to be Maxwell-Boltzmann, with dispersion  km s. In the frame of the Sun, moving at  km s relative to the Galactic rest frame,

 f⊙(u)=(32)3/24ρχu2π1/2mχu30exp(−3(u2⊙+u2)2u20)sinh(3uu⊙/u20)3uu⊙/u20. (2.2)

As it falls into the gravitational potential well of the Sun, a DM particle acquires a velocity . It becomes gravitationally captured by the Sun if it loses enough kinetic energy in a scattering event for to fall below the local escape velocity . For this to occur, the fractional energy lost by the DM particle , corresponding to nuclear recoil energy , must be in the interval

 u2w2≤2ERmχw2≤μμ2+, (2.3)

where and . The local capture rate of particles with velocity is the sum over nuclear species at radius , of the rate of scattering in the interval of Eq. 2.3, i.e.

 (2.4)

where is the nuclear form factor. For hydrogen this is a constant, whereas for heavier nuclei we use the usual Helm form factor

 |Fi(ER)|2=exp(−EREi). (2.5)

Here is a constant quantity for each nuclear species , given by

 Ei=5.8407×10−2mN,i(0.91m1/3N,i+0.3)2GeV. (2.6)

Integrating over the phase space, the total capture rate of DM in the Sun is then

 C⊙(t)=4π∫R⊙0r2∫∞0f⊙(u)uwΩ(w)dudr. (2.7)

For a constant cross-section and a Maxwell-Boltzmann velocity distribution, this can be solved analytically [11, 96]. However, in the following we generalise the above equations to allow momentum and velocity-dependent and .222We will use as a shorthand for either or , as the spin and kinematic dependence can be factorised. For an explicit treatment, see [97]. In this case, numerical integration becomes necessary.

### 2.2 Velocity and momentum-dependent treatment

We begin with a momentum-dependent cross-section of the form Eq. 1.1a, with . To include such a cross-section in Eq. 2.4, the constant cross-section must be replaced with Eq. 1.1a and the dependence on the nuclear recoil energy moved inside the form-factor integral. This explicitly illustrates the equivalence of a momentum-dependent cross-section to a change in form factor, and we refer to the corresponding integral and associated multiplicative factors as the ‘generalised form factor integral’ (). For hydrogen, , so the change required is

 GFFIn=0=∫|F(ER)|2dER → GFFIn≠0,H=(pq0)2nmχw22μn∫μ/μ2+u2/w2ΔndΔ, (2.8)

where we have expressed the form factor integral in terms of instead of for compactness. We do this by noting that the transferred momentum is , so the fractional energy change can be written , where is DM particle’s incoming momentum. Performing the integral in Eq. 2.8 yields

 (2.9)

For heavier elements, the integrand includes the Helm form factor (Eq. 2.5), so

 GFFIn≠0,i≠H=(pq0)2nEi(Bμ)n[Γ(1+n,Bu2w2)−Γ(1+n,Bμμ2+)], (2.10)

where , and is the (upper) incomplete gamma function. To gain a more intuitive understanding of Eqs. 2.9 and 2.10 we show in Fig. 1 the overall enhancement or suppression to the capture rate with respect to the constant case,

 Fcapture=∑ifiA2iGFFIn≠0,i∑ifiA2iGFFIn=0,i, (2.11)

for three values of . In this example we just take and use the present-day AGSS09ph Standard Solar Model333Publicly available at http://www.ice.csic.es/personal/aldos/Solar_Models.html. [70, 61]. Here represents the fractional composition of each species with atomic number . For the spin-dependent cross-sections, the sum over the species only includes hydrogen. Of course the actual capture rate must be accurately integrated over the entire star, and we compute it precisely with Eq. 2.7 in our simulation. In each case, the overall behaviour roughly follows the dependence, as heavier DM will have gained more momentum as it falls into the solar gravitational well.

In the velocity-dependent case (Eq. 1.1b), where , the modification to the capture rate is much simpler. The partial capture rate is simply transformed by the replacement

 σ→σ0(wv0)2n. (2.12)

The integral Eq. 2.7 can then be evaluated to obtain the modified capture rate. We later show with our solar simulation code that the enhancement due to velocity-dependence agrees very well with what is obtained using the standard () treatment and an average of the cross-section throughout the volume of the star

 σ(vrel)≃⟨σ(w0)⟩=1M⊙4π∫R⊙0σ(w0)ρ⊙(r)r2dr, (2.13)

where is the typical velocity acquired by an in-falling WIMP as it reaches radius . This indicates that a cross-section proportional to or – typically thought of as a suppression – actually yields a considerable enhancement. We note that has been erroneously approximated to rather than in the literature [90], although does at least provide the correct enhancement to the capture rate to within a factor of a few.

#### 2.2.1 Form Factors

We finally comment on the use of the Helm form factor 2.5 in our capture equations. Very recent work [97] has shown that a more accurate computation of the nuclear response functions to DM-nucleus scattering can have a substantial impact on the capture rate. These response functions amount to corrections to the Helm form factor by terms proportional to powers of . These will dominate at large momentum transfers. We have checked that in most cases considered here – and in all cases where the effect of DM is large enough to modify the solar observables – this quantity is much smaller than one. This is because the low DM masses that yield the largest conduction effects are mainly sensitive to the leading, constant term in the NR expression. However, we caution the reader that accurate modification of the capture rate for DM with larger masses than those considered here, should make use of the full nuclear response functions detailed in Ref. [97].

We have also used the standard approach for spin-dependent capture in the Sun, wherein it is assumed that hydrogen is the dominant contributor to the capture rate. It was shown in [97] that, while this is strictly true for velocity-dependent cross sections (non-relativistic operator ) and holds to a few percent accuracy for a constant cross section (), momentum-dependent interactions can lead to a different behaviour. In the case (), nitrogen dominates capture by an order of magnitude or more for all DM masses. The net effect is that the spin-dependent capture rates computed here for have been underestimated with respect to the case when the full nuclear response functions are included.

### 2.3 The geometric limit

In all cases, the total effective cross-section of the Sun to collisions with DM particles cannot exceed the geometric “cutoff”

 σmax=πR2⊙(t), (2.14)

which corresponds to the case where the Sun is optically thick to DM. This places a fundamental limit on the capture rate

 Cmax(t) = πR2⊙(t)∫∞0f⊙(u)uw2(u,R⊙)du = 13πρχmχR2⊙(t)⎛⎜⎝e−32u2⊙u20√6πu0+6GNM⊙+R⊙(u20+3u2⊙)R⊙u⊙Erf[√32u⊙u0]⎞⎟⎠.

The actual capture rate to be used must therefore be the lesser of Eqs. 2.3 and 2.7. Assuming a steady radius and a local DM density of GeV cm, the maximum amount of DM that can accumulate in the Sun by its current age is therefore

 Nχ,max≃1.5×1047(GeVmχ)≃1.3×10−11(GeVmχ)Nb, (2.16)

where is the total number of baryons in the Sun.

## 3 Conductive energy transport by dark matter

If enough DM is captured by the Sun, its large typical inter-scattering distance means that it is a more efficient carrier of heat over long distances than ordinary baryonic material, so it can act as an additional mechanism for heat transport alongside photons. The formalism we use to compute the effect of microscopic energy transport by DM conduction was developed by Gould and Raffelt [92]. By solving a perturbative expansion of the Boltzmann collision equation (BCE) in the Sun’s gravitational potential, they showed that the thermal conduction by a weakly-interacting species can be expressed in terms of two quantities: a dimensionless molecular diffusivity , and thermal conductivity , where is again the ratio between the DM and nucleon masses. If more than one nuclear species is present, and get replaced with effective values, which are weighted by the number densities of each species in the plasma mixture at each height in the Sun.

In the local thermal equilibrium (LTE) regime, where is much smaller than both the inverse of the local temperature gradient and the DM scale height , the values of and are found by solving the first order expansion of the BCE. This is done in terms of the quantity , via the formal inversion of the Boltzmann collision operator , where represents the net change in the DM phase space distribution due to collisions with nuclei. The collision operator is defined via phase space integrals over collision rates, meaning that the dependence of the collisional cross-section on and must be explicitly included in the calculation. In Ref. [92] Gould and Raffelt computed and tabulated and for a constant scattering cross-section (). In Ref. [91], we extended the formalism to include velocity and momentum-dependent cross-sections (), showing that the resulting changes in , and can have potentially large effects on heat transfer in the Sun.

The equilibrium distribution of DM particles in the gravitational potential of the sun is given by [92, 11, 91]

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

where represents the centre of the Sun. The conductive luminosity is:

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

where the factor accounts for a velocity-dependent or momentum-dependent cross-section, and is related to the typical thermal velocity [92]; and are respectively the reference velocity and momentum defined in Eq. 1.1. The rate of energy transported per unit mass of stellar material is:

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

This quantity is usually expressed in units of ergs g s.

If the condition is violated, LTE is no longer valid and the system exists in the Knudsen regime of non-local transport, where the DM distribution is essentially isothermal. It was shown in Refs. [92, 98] by Monte Carlo simulation that Eqs. 3.1 and 3.2 can be corrected to properly account for a large Knudsen number . Here we have formally defined

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

Following [52, 11], we define:

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

where and are empirical values taken from the numerical results of [98]. The DM distribution is then a combination of the isothermal and LTE distributions:

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

where

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

Finally, the Knudsen-corrected luminosity is:

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

This should then be used in place of in Eq. 3.3 to compute the energy injected or removed at each radius by DM-nucleon collisions. In Fig. 1 we illustrate the effective enhancement or suppression of transport in a toy solar model due to a momentum-dependent cross-section, plotting

 Ftransport≡∫|ϵ(r,n≠0)|r2dr∫|ϵ(r,n=0)|r2dr, (3.10)

for the three cases of momentum-dependent cross-section, with cm. Once again, this is illustrated using a present-day SSM with the AGSS09ph abundances [70, 61]. Note that this cross-section leads mainly to transport near the LTE regime. As is decreased further towards the Knudsen regime, the behaviour reverses, as the enhancement provided by the longer inter-scattering distance is overcome by the Knudsen suppression. We illustrate the Knudsen behaviour in Fig. 2 by plotting the total energy transport for different types of generalised form factor DM as a function of the cross-section , for a constant DM-to-baryon ratio . This shows how the peak energy transport varies in each model. We note that this peak occurs for a much smaller cross-section in the spin-independent case, due to the reduced mean free path caused by scattering with helium and metals rather than just hydrogen.

The full effect of ADM on energy transport is then a combination of the effects illustrated in the left and right panels of Fig. 1, keeping in mind the degree of non-locality indicated by Fig. 2.

## 4 The DarkStec solar dark matter code

In order to accurately model the effects of generalised form factor DM on solar observables, we implemented full velocity- and momentum-dependent DM capture and energy transport in the high-precision solar evolution code GARSTEC [99, 100]. We took DM routines from the public dark stellar evolution code DarkStars [95], producing a hybrid code DarkStec.

GARSTEC is the descendant of the legendary Kippenhahn code. Numerical aspects and physics inputs are described in detail in Ref. [99] and the modified version of GARSTEC used for this work is the same described in Ref. [100]. Here we just give a summary of the most relevant physical inputs. It includes the nuclear energy generation routine exportenergy.f444Publicly available at http://www.sns.ias.edu/~jnb., updated with the astrophysical factors recommended in the Solar Fusion II [101] compilation. It makes use of the Opacity Project radiative opacities [102], complemented at low temperatures with those from Ref. [103]. The equation of state is the 2005 update of OPAL [104]. Microscopic diffusion of elements, including gravitational settling, thermal and concentration mixing, is treated according to Ref. [105].

The calibration of a solar model generally implies adjusting a number of free parameters in the model to match an equal number of observables. In the present case, the observables are the present-day ( Gyr) solar luminosity , solar radius and the metal-to-hydrogen mass fraction . The latter is a critical quantity, as it determines the composition, namely the metallicity, of the calibrated model. In this work, we adopt the photospheric solar abundances from Ref. [70], for which . The free parameters in the model are the mixing length parameter and initial helium and metal mass fractions and respectively. The latter two suffice to determine the initial abundances of all elements in the model because the relative metal abundances are taken from Ref. [70] with acting as the normalisation factor, and by definition ( being the initial hydrogen abundance).

In practice, the solar model calibration starts with a homogeneously-mixed 1 M pre-main sequence model that is evolved assuming no mass loss until it reaches , the current solar system age. At that age, model predictions are compared with the observables, and a Newton-Raphson scheme is implemented to iteratively find the solution. This is generally achieved to 1 part in 10 within two to three iterations for standard solar models (i.e. with no DM). More iterations are necessary as the effects of DM become more important. In the most extreme cases, no physical solutions are found (e.g. resulting in negative ). Note that each iteration requires four evolutionary calculations to evaluate the partial derivatives needed for the Newton-Raphson scheme. Each evolutionary calculation requires 700–800 timesteps. At each timestep, the solar structure is discretized in about 2000 shells. These requirements for the integration of solar models, both in spatial and time resolution, guarantee a numerical precision better than 1% in all model predictions [58].

DarkStars [95] is a Fortran95 package that implements capture, annihilation and energy transport by regular () WIMP dark matter in a general stellar evolution code, as described in Refs. [11, 40, 39, 106]. The capture routines were originally adapted from DarkSUSY [107]. The underlying evolutionary code [108] is a Fortran90 rewrite of the the venerable Fortran77 Cambridge STARS package [109, 110, 111]. These codes use the relaxation method to solve the coupled 1D ordinary differential equations of stellar structure over an adaptive grid. DarkStars is the state of the art in dark stellar evolution for the case, as it features the full capture calculation (Eq. 2.7) for SI and SD scattering on the 22 most important elements, various options for the DM velocity distribution (including user-defined distributions), and proper Gould-Raffelt treatment of conductive energy transport (Eq. 3.2), including self-consistent density profiles and the Knudsen-dependent interpolation between the LTE and isothermal (non-local) regimes (Eqs. 3.7, 3.9).

For DarkStec, we adapted the DarkStars capture routines to implement capture of generalised form factor DM (Eq. 2.92.10, 2.12) rather than just the case. We used the conductive transport routines from DarkStars essentially unaltered, except that we included the additional factor of in Eq. 3.2 and utilised the and tables that we computed earlier for [91].

At each regular GARSTEC timestep, DarkStec computes the total DM capture rate by solving Eq.  2.7, assuming a local halo DM density of 0.38 GeV/cm. This input rate is then used to update the DM population in the Sun following Eq. 2.1. DarkStec uses the numerical version of the capture routines from DarkStars to evaluate the modified capture equation Eq. 2.7. It computes the thermal diffusivity and conductivity coefficients and at each height in the star by interpolating in the tables of Ref. [91], as these quantities depend on the specific mixture of plasma species with which the DM particles interact. The code then computes the DM density using and , which it uses to determine energy transport. DarkStec then interpolates the resulting values of (Eq. 3.3) to the grid used by GARSTEC, where they are treated as an additional energy source at each height in the star.

We considered seven ADM models with SI interactions with nucleons, and seven with SD interactions: the constant cross-section case, and with . For each coupling, we simulated one solar model for each point on a grid of and . We computed a subset of models using a stringent convergence criterion of one part change in for the the solar luminosity, radius and surface metallicity (), with 4 kyr and 10 Myr minimum and maximum time steps. Because it was extremely computationally expensive to do every simulation including the full treatment of capture and transport at this accuracy, we carried out all other simulations with a more relaxed convergence criterion of a part in , with minimum and maximum time steps of 40 kyr and 40 Myr. We then corrected these lower-accuracy results to consistently reproduce the observables and chi-squared values of the higher-accuracy models, using the systematic differences we saw between models computed both ways. In total, our calculations took CPU years.

### 4.1 Annihilation

DM models with momentum- and velocity-dependent nuclear scattering need not necessarily be entirely asymmetric. To investigate the implications of energy injection from annihilation in such models, we also implemented annihilation in DarkStec, following [11] and DarkStars. In general however, allowing annihilation simply weakens the limits that we obtain from solar physics on generalised form factor DM, and does not improve the overall fit to solar data in any significant way. We therefore do not show these results, nor discuss annihilation beyond this subsection. It is worth noting that although we assume zero annihilation cross-section for all the results we show, some small (sub-thermal) annihilation cross-section is certainly still allowed in each model, and would make no impact on the solar observables.

## 5 Limits on generalised form factor dark matter from solar physics

In this section we present a systematic overview of the results obtained from our simulations. In each subsection, we review the effects of velocity and momentum-dependent asymmetric dark matter on a specific observable. These are: the boron-8 and beryllium-7 neutrino fluxes, the depth of the convection zone , the sound speed profile , the small frequency separations and the surface helium abundance . Agreement between the predicted and observed values of the neutrino fluxes is typically unchanged or worsened by thermal transport by DM. The same is true for . In contrast, the predictions of , the sound speed profile and the small frequency separations are often closer to the observed values when conductive energy transport by DM is included. In Sec. 5.7 we construct a combined likelihood, which encompasses the overall improvement or degradation in the fit to solar data for each cross-section form and combination of and . We also present a few interesting benchmark cases in which the agreement is significantly improved.

For every form of the cross-section, we ran simulations over a grid of DM masses and cross-sections: 5, 10, 15, 20 and 25 GeV, and each decade in between and cm. The colour scales in the figures of this section are interpolations from this grid. For some specific cases, where low-mass points gave a good overall improvement over the Standard Solar Model, we extended the mass axis down to 1 GeV. We once again caution that although evaporation is expected to have an effect near these small masses, a full kinematic analysis would be necessary to determine its exact and dependence for each model; the importance of evaporation for these models remains essentially unexplored. Given the importance of kinematic matching with individual nuclei, we caution against placing too much store in quick estimates of this effect.

We show the impacts of generalised form factor DM on capture rates and their saturation in Fig. 3, and on observables in Figs. 4 to 17.

Simulations that led to a reduction in of approximately or more did not converge, because the effects on the overall solar evolution are too large for the numerical solver to handle. In many cases, the excess energy transport actually led to a density inversion in the core. In some of these models, a solution probably exists, and could be found with a better solver. In others, the Sun probably cannot exist as a stable body. We simply mask the the non-converged region in blue in all our plots.

### 5.1 Saturation of capture

A velocity or momentum dependence in the nuclear scattering cross-section has a striking impact on the rate of DM capture by the Sun. The minimal required to render the Sun opaque to dark matter – and thus to saturate the capture rate Eq. 2.3 – can be lowered by as much as two orders of magnitude for a cross-section, and approximately one order of magnitude in the and cases. On the other hand, negative powers of and positive powers of a spin-dependent -coupling yield a large suppression in the capture rate, seen as a much larger required cross-section to achieve saturation. We illustrate the required value of in Fig. 3, based on the output of our simulations. Interestingly, a spin-independent or cross-section can yield an enhancement or suppression of the capture rate depending on the DM mass; this is a consequence of the behaviour illustrated in Fig. 1, and simply reflects the fact that -dependent scattering is most efficient for large momentum transfers, which is much easier when the DM mass is closely matched with the masses of heavier elements.

### 5.2 Solar neutrino fluxes

In general, energy transport by DM removes energy from the solar core and reduces its temperature, causing a reduction in neutrino production rates. Given its strong temperature-dependence, the B neutrino production rate is the first place to look for changes in solar observables.

In Figs. 4 and 5, we show the effect on the B neutrino flux for spin-independent and spin-dependent dark matter with velocity and momentum-dependent couplings. Fluxes are normalized to the observed value, cms; lighter colouring represents a reduced flux. Although the measurement error on the B neutrino flux is only 3% [112], the overall uncertainty from modelling is around 14%. We add the absolute uncertainties in quadrature to obtain the 1 region. White lines therefore represent the isocontour where the neutrino flux falls below 85% of the measured value; black lines are the 2 (71%) contours.

Although the flux of Be neutrinos is not as temperature-sensitive as that of B neutrinos, they can also be used as a weaker, independent, probe of the solar core temperature. In Figs. 6 and 7 we show the ratio of the predicted Be neutrino fluxes to the observed value cms. Here again we plot and contours using the theoretical and observational uncertainties, which are respectively 7% and 5% for Be neutrinos.

As expected, the effect of ADM on neutrino fluxes is larger for low DM masses, mainly due to the increased efficiency of momentum transfer between DM and H/He as the DM takes on a similar mass to these two dominant nuclei. Ignoring evaporation, transport becomes even more efficient as mass is decreased to 1 GeV for SD dark matter; in the SI case, the largest effect is seen at 3 to 4 GeV, around the average nucleus mass in the Sun.

The dependence on the cross-section is also as expected: below a minimum value of , not enough DM can be captured to significantly affect the B production rate. For constant cross-sections, this is around cm in the SI case, and cm for SD dark matter. This highlights the importance of heavier elements, to which only the SI DM couples in our simulations: the large fraction of helium, combined with the dependence of the DM-nucleus cross-section (Eq. 1.2) means that the behaviour of DM inside a star is highly sensitive to elements heavier than hydrogen.

The upper edges of the regions shown in Figs. 47 highlight the transport behaviour: if the cross-section falls below a certain value, it allows the DM to more efficiently penetrate the dense plasma, carrying energy from hot to cooler regions. The location of the maximum reflects a combination of being near the transport efficiency peak (where the transition from the LTE to the Knudsen regime occurs) while simultaneously maintaining a large capture rate through a large enough cross-section.

The impacts of velocity and momentum-dependent scattering of dark matter can be understood in terms of this balance between capture and transport. Positive powers of lead to an enhancement in the capture rate, moving the window in which DM transport has a significant effect on neutrino fluxes to lower cross-sections. However, the important suppression in the overall transport rate, as seen in Fig. 5 of [91], yields an overall suppression in the effect relative to the constant case. The opposite behaviour is seen for , although the enhanced transport is not sufficient to compensate for the suppression in capture.

Again, the effect of a momentum-dependent cross-section is more subtle. Although the boosted capture rate for a cross-section moves the region of interest to lower cross-sections, closer to what is allowed by underground experiments, the overall suppression in transport means that there is very little overall effect on solar observables. Positive powers of fare much better, with an enhancement both in capture and transport.

Comparison of Figs. 4 and 5 highlights the importance of heavier elements in such processes. In most cases, they lead to a significant enhancement of the capture and transport rate; however, we note in comparing the plots that they can also inhibit energy transport, by reducing the mean free path for conduction.

### 5.3 Depth of the convection zone

In Figs. 8 and 9, we show the ratio of the location of the lower boundary of the convection zone predicted by each model to the value inferred from helioseismology, .

Although the error on the helioseismological inference is just 0.001 , the theoretical error on the modelled is much larger: . Adding these errors in quadrature, we plot in white and black the contours containing the regions where is within 1 and , respectively, of the inferred value.

As can be seen in Figs. 8 and 9, the difference between the depth of the convection zone in the Standard Solar Model () and the depth inferred from helioseismology amounts to almost a 3 discrepancy. Except for a small region at high cross-section ( cm), spin-dependent DM struggles to bring to within much better than 2 of the inferred value. However, energy conduction from SI interactions does substantially better: and models can produce good agreement with the observed value of , leading in some regions to agreement at better than 1.

### 5.4 Surface helium abundance

The surface helium abundance is another observable that has fallen into disagreement with observations since the revision of the solar composition. Our SSM prediction is , whereas the observed value is 0.2485 . With theoretical errors of taken into account, this amounts to a discrepancy of . The addition of dark matter does very little to change , and for most models the discrepancy remains at the level. For some very few cases (not shown), it actually worsens the discrepancy by up to an additional . This only occurs for models where the fit to sound speed observables is also substantially worsened, notably for large SD and cross-sections. The reduction in is caused by a corresponding increase in , the initial hydrogen fraction. The increase in is demanded by the reduction of the core temperature, which requires a greater amount of hydrogen to be present in the core in order to maintain the same nuclear burning rate as in the SSM, and to thereby match the observed solar luminosity .

Because it changes little, we do not show the contour plots for the surface He abundance. For completeness though, we include it anyway in our full likelihood computation in Sec. 5.7.

### 5.5 Sound speed profile

In Fig. 10 we show the deviation of the modelled radial sound speed of some models that produced the best overall fits with respect to the measured values from helioseismology. To illustrate the agreement of each model with the observed sound speed profile, we construct an effective chi-squared measure

 χ2cs=∑ri[cs,model(ri)−cs,hel.(ri)]2σ2cs,hel.(ri) (5.1)

where the sum goes over 5 equally-spaced radial points between and . We do not include radii below because the reconstructed sound speed in this region is highly uncertain due to the low number of low-degree (low-) p-modes reaching the innermost radii of the solar core [113]. The upper limit of the range is the location of the largest discrepancy in the sound speed of the Standard Solar Model, at the base of the convection zone; above this point essentially all models agree very well with the observed sound speed because the temperature gradient is adiabatic and therefore does not depend on the detailed composition of the Sun. We added errors from modelling and inversions for each point in quadrature. We obtained modelling errors by using models for which one input parameter was varied at a time to obtain partial derivatives at each radial point, and then combined quadratically given the uncertainty in the AGSS09 abundances and errors in each parameter quoted in [114]. The errors on inversions were taken from [115]. The values of this , meant to show the relative improvement to the sound speed modelling, are shown in Figs. 11 and 12.

From these figures, along with the sound speed profiles illustrated in Fig. 10, it is clear that the addition of momentum or velocity-dependent dark matter to the solar model does indeed help alleviate the discrepancy between modelling and observation in some specific cases. We will return to these cases in Sec. 6.1.

### 5.6 Frequency separation ratios

The inverted sound speed profile obtained from helioseismological measurements is not very accurate near the solar core because not enough low- p-modes are available for a precise and accurate inversion. Instead, information about the solar core can be gained by using the so-called frequency separation ratios. A large advantage of using these ratios is that, unlike individual frequencies, they are not affected by the detailed structure of the solar surface, which is poorly described by solar models [116]. This is because for radial order and low angular degree , the surface effects are functions of the eigenfrequency, so they cancel out when considering frequency differences between modes of similar frequencies. In addition, by taking ratios of appropriate frequency differences, the solar core structure becomes the dominant effect in the observed signal [117, 118]. In particular, two very useful quantities are the frequency separation ratios

 r02(n)=d02(n)Δ1(n),r13(n)=d13(n)Δ0(n+1), (5.2)

where

 dl,l+2(n)≡νn,l−νn−1,l+2≃−(4l+6)Δl(n)4π2νn,l∫R⊙0dcsdrdrr. (5.3)

and . The sound speed gradient is weighted by in the integral above, so and are most sensitive to changes there.

We show a set of these ratios in Fig. 13 for the same examples shown in Fig. 10. These are compared with the values of and measured by the BiSON experiment [118]. Model uncertainties are computed following the same procedure as for the sound speed. As can be seen in the residual subplots, the SSM overestimates each of these ratios by as much as 4, once modelling errors are included. Thermal transport by ADM can smooth the sound speed gradient, yielding smaller differences in Eq. 5.3. As in the case of the sound speed profile, the best improvement is provided by scattering with GeV, cm. In this model, the largest discrepancy in and falls to barely .

In Figs. 14 and 15, we show the overall ability of different models to fit the observed frequency separations. We quantify this with the combined chi-squared , where

 χ2rℓ,ℓ+2=∑n[rℓ,ℓ+2,th.(n)−rℓ,ℓ+2,obs.(n)]2σ2obs.(n)+σ2th.(n). (5.4)

Formally, there is a correlation between values of and with common eigenfrequencies. In practice, however, the correlation is minimal (the Pearson correlation coefficient is always lower than 0.1), so we assume they are independent when computing . Figs. 14 and 15 show that all SI and SD couplings do exhibit areas of parameter space where the small frequency separations can be brought into better agreement with observations than in the SSM. In many cases however, these are not the same parameter combinations preferred by other observables. The notable exception to this is the SI coupling, which also clearly produces the best agreement with the observed frequency separations.

### 5.7 Combined limits

It is clear from the results we have presented in Secs. 5.25.6 that the disagreement between model predictions and helioseismology can be ameliorated, at the cost of a slightly worse agreement with neutrino fluxes. Surface helium Y is also affected to a lesser extent. It is therefore instructive to construct an overall likelihood function based on the fit of these data by our models. We define the combined chi-squared as:

 χ2 = (1−ϕνB,obs/ϕνB)2σ2B+(1−ϕνBe,obs/ϕνBe)2σ2Be (5.5) +(rCZ−rCZ,obs)2σ2CZ+(YS−YS,obs)2σ2YS +χ2r02+χ2r13.

where , , and . The uncertainties include both the observational and modelling errors, added in quadrature. These are given individually in the last two lines of Table 1. We include the small frequency separations , but not the sound speed profile, as the latter is less precise and is correlated with the former.

We show the combined limits from Eq. 5.5 in Figs. 16 and 17, comparing with limits from direct detection where such results exist (i.e. in the SI case; [30]). We discuss these results, and the comparison with other limits, in more detail in the following section.

In Table 1, we give a break-down of the values of each observable at the best-fit parameters of each model, along with the best-fit and implied -value. The best-fit regular SI and SD cases shown in Table 1 differ from the best-fit cases in Ref. [74] because here we actually choose them on the basis of the full likelihood Eq. 5.5. In Ref. [74] we only computed the small-frequency separations for models that provided the best fits to all other observables. Because the fits to different observables are inconsistent in the models with constant cross-sections, including the frequency separations shifts the best-fit masses and cross-sections. This is to be compared to the SI case, where adding in the small frequencies merely makes the best-fit point even better.

We also provide as supplementary material online a complete table of the boron and beryllium neutrino fluxes, surface helium, convective zone radius and small frequency separations for all of our models, in addition to the partial and total chi-squared values defined in Eq. 5.5. This table can be used for quick lookup and interpolation for, e.g. global fits of DM models.

## 6 Discussion

### 6.1 Solving the Solar Abundance Problem

Our results of Section 5 indicate that the impacts of different types of generalised form factor DM on solar observables is quite varied, and not always straightforward. What is clear, however, from Figs. 16 and 17 and Table 1, is that some models do indeed lead to a much better overall fit than the Standard Solar Model alone. Although there is always some reduction in the goodness of fit to observed neutrino fluxes, this can be accommodated in many cases by the theoretical errors.

In Fig. 10 we showed the reduction of the discrepancy between simulated and measured sound speed profiles of the best-fit models of the couplings that give the best overall improvement compared to the Standard Solar Model. On the left, we showed the best models for constant spin-dependent and spin-independent cross-sections, and on the right we showed the best-fit models of the three best generalised form factor couplings, according to the likelihood defined in Sec. 5.7. We find that the best improvement comes from a scattering cross-section that is proportional to the square of the transferred momentum . This is confirmed to very high significance by the small frequency separations, seen in Fig. 14. At low masses GeV, an asymmetric DM candidate with and cm leads to a substantial improvement with respect to the SSM. The best fit occurs at 3 GeV, and leads to a convective zone depth , only 1.2 away from the inferred value from helioseismology. The B and Be fluxes are respectively  cms and  cms, both within of the measured values. Surface helium remains low, at . The large improvement to the fit to can be seen in the right-hand panel of Fig. 10, and to the small frequency separations probing the solar core in Fig. 13. This specific case yields an overall 6 improvement compared to the Standard Solar Model, as we discussed explicitly in Ref. [74].

In the next section, we will discuss constraints from other sectors, and show that the best-fit model does indeed constitute a viable candidate for reconciling solar modelling, helioseismology and spectroscopy. It is worth noting, however, that the masked regions of non-convergence in our contour plots may hide parameter combinations that yield even better fits than this; careful work on the numerical solver is required to explore this possibility.

Other -dependent models, as well as constant and velocity-dependent models we examined, nearly all led to poor overall fits: the effect on and is always too large, in spite of the improved fit to helioseismological measurements. Every other ADM model we examined led to some significant improvement over the SSM, but failed to provide -values greater than 10. As pointed out in Ref. [90], a SD cross-section does indeed lead to an improved fit to the sound speed profile and for a small region in mass and cross-section space (see Figs. 9 and 12). Unfortunately, this case is disfavoured by our overall likelihood, mainly due to excessive reduction in the core temperature, which is reflected in a greater-than reduction in the B neutrino fluxes, as well as significant reduction in Be neutrinos. In fact, we find that a SD cross-section yields a much more plausible overall fit to observables than SD scattering, resulting in a -value of 3.1 as compared to 4.2 (vs. 0.85 for the SI model). It is also important to note that the region where SD improves helioseismological agreement is very close to the boundary beyond which our models failed to converge. This issue, also seen by the authors of Ref. [54], could indeed be hiding a better overall fit that is simply beyond the reach of our solver. Given the great reduction in neutrino fluxes that we observe in the bordering area however, we remain somewhat sceptical of this possibility.

### 6.2 Comparison with collider and direct limits

The dashed curves in Fig. 16 show independent limits from direct detection on SI scattering, as derived in Ref. [30]. To the best of our knowledge, similarly model-independent limits on and -dependent SD couplings from direct detection are not available in the literature. For momentum-dependent couplings, the direct limits suggest that all interesting masses (from the point of view of solar physics) above 5 GeV are in tension with direct searches. Curiously though, the areas that lead to the best overall improvements in solar observables ( scattering, –5 GeV,