Generalised form factor dark matter in the Sun

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 $v_{\rm rel}^{2n}$ and $q^{2n}$ with $n = -1, 0, 1 $ or $2$, where $v_{\rm rel}$ is the dark matter-nucleon relative velocity and $q$ 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 $\sigma_{SI} \propto q^2$ is particularly appealing, yielding more than a $6\sigma$ improvement with respect to standard solar models, while being allowed by direct detection and collider limits. We provide full analytical capture expressions for $q$- and $v_{\rm rel}$-dependent scattering, as well as complete likelihood tables for all models.


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 freezeout (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 m χ = 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 −1 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.
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), 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 (1 − 10 GeV), as momentum transfer is maximised when the masses of the colliding particles are equal, and the Sun is mostly H (A = 1) and He (A = 4). It is interesting to note that this is roughly the mass range expected in the most common models of ADM [48], 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 (A ∼ 72) and xenon (A ∼ 131), for which high incoming DM velocities are necessary to create a measurable recoil.
There is a window of elastic DM-nucleon scattering cross-sections [49][50][51][52][53] at σ ∼ 10 −36 cm 2 for spin-independent (SI) couplings (10 −34 cm 2 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 T c , 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 8 B neutrinos, which goes as T β c 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 [51][52][53]. At the same time, it appears possible for the inclusion of DM in such models to affect the (less well constrained) 8 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 [54][55][56][57][58][59]. 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 [60][61][62][63][64][65][66][67][68][69][70][71], 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 [72] 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.

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 s, t and u. In the non-relativistic limit, these become the centre of mass momentum, proportional to the relative velocity v rel , and the transferred momentum q = ∆p. 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 v rel and q exist, typically motivated by theoretical arguments or attempts to reconcile experimental results.
To make quantitative predictions, a specific form of σ(v rel , q) must be chosen. Because we wish to remain as general as possible, we choose to focus on couplings of the form σ ∝ v 2n rel and σ ∝ q 2n , with n = {−1, 1, 2}. The v 2 rel and v 4 rel forms are respectively called p-wave and d-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 (s-wave) term -are suppressed due to cancellations. Cross-sections depending on q can arise, for example, from a non-zero particle radius (the analogue of a nuclear form factor), from parity-violating couplings likeχγ 5 χQQ,χγ µ γ 5 χQγ µ Q andχσ µν γ 5 χQσ µν Q, or from a small anapole or dipole interaction between the dark and visible sectors [73][74][75][76][77][78][79][80][81][82][83][84][85]. We refer to the class of models where DM-nucleon scattering cross-sections depend on some combination of q and/or v rel as 'generalised form factor DM' because it generalises the effects of form factors to arbitrary powers of q and v rel .
Concretely, we focus on the couplings 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 σ 0 must be defined with respect to some reference velocity v 0 or momentum q 0 . We will choose v 0 = 220 km s −1 , the typical halo DM velocity, and q 0 = 40 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: where A i and J i are respectively the atomic number and total angular momentum of nuclear species i; S p,i and S n,i are the spin expectation values of its proton and neutron systems. Given the A 2 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. 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. [86,87], the authors computed capture and transport rates for ADM models with long-range interactions. To account for the impacts of the nontrivial 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 [88].
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 momentumdependent scattering of DM with nucleons. In Section 3 we review the theory of conductive heat transport by DM developed in Refs. [88,89], 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

Standard (velocity and momentum independent) treatment
The population of DM particles in the Sun N χ (t) is follows the differential equation where C (t) is the capture rate, A(t) is rate at which annihilations occur and E(t) represents evaporation. Unless DM is strongly self-interacting (not the case we consider here, but discussed in Ref. [90]), C (t) does not depend on the DM population already captured by the Sun. A(t) 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 A(t) = 0, 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 E(t) 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 ∼100 km s −1 , whereas the escape velocity in the solar core approaches 1400 km s −1 . 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 =⇒ m χ ∼ 12, 14, 16, 56 GeV) can also in principle be affected [91]. 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 [92], we closely follow Refs. [11,93]. We take the local distribution function f (u) of DM to be Maxwell-Boltzmann, with dispersion u 0 = 270 km s −1 . In the frame of the Sun, moving at u = 220 km s −1 relative to the Galactic rest frame, (2.2) As it falls into the gravitational potential well of the Sun, a DM particle acquires a velocity w = u 2 + v 2 esc (r, t). It becomes gravitationally captured by the Sun if it loses enough kinetic energy in a scattering event for w to fall below the local escape velocity v esc (r, t). For this to occur, the fractional energy lost by the DM particle ∆ = 2E R /m χ w 2 , corresponding to nuclear recoil energy E R , must be in the interval where µ ≡ m χ /m nuc and µ ± ≡ (µ ± 1)/2. The local capture rate of particles with velocity w is the sum over nuclear species at radius r, of the rate of scattering in the interval of Eq.
where |F i (E R )| 2 is the nuclear form factor. For hydrogen this is a constant, whereas for heavier nuclei we use the usual Helm form factor Here E i is a constant quantity for each nuclear species i, given by GeV. (2.6) Integrating over the phase space, the total capture rate of DM in the Sun is then For a constant cross-section and a Maxwell-Boltzmann velocity distribution, this can be solved analytically [11,93]. However, in the following we generalise the above equations to allow momentum and velocity-dependent σ SI and σ SD . 1 In this case, numerical integration becomes necessary.

Velocity and momentum-dependent treatment
We begin with a momentum-dependent cross-section of the form Eq. 1.1a, with σ ∝ q 2n . 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' (GF F I). For hydrogen, |F (E R )| 2 = µ 2 + /µ, so the change required is where we have expressed the form factor integral in terms of ∆ instead of E R for compactness. We do this by noting that the transferred momentum is q = √ 2m nuc E R , so the fractional energy change can be written ∆ = µq 2 /p 2 , where p = m χ w is DM particle's incoming momentum. Performing the integral in Eq. 2.8 yields (2.9) 1 We will use σ as a shorthand for either σSI or σSD, as the spin and kinematic dependence can be factorised. For an explicit treatment, see [94].    [59]), due to a momentum-dependent cross-section, as computed with Eq. 2.11. The black line at F = 1 corresponds to σ = const. Solid lines are spin-independent (SI) and dashed lines are spin-dependent (SD), coupling only to hydrogen. Right: effective enhancement/suppression of energy transport as defined in Eq. 3.10 , in the LTE regime, with σ 0 = 10 −35 cm 2 . We note that away from LTE, the behaviours reverse, and a q −2 cross-section actually yields an enhancement with respect to the constant case. This can be seen in Fig. 2. The full effect of a q-dependent cross-section is then the combined effect of the left and right-hand panels.
For heavier elements, the integrand includes the Helm form factor (Eq. 2.5), so (2.10) where B ≡ 1 2 m χ w 2 /E i , and Γ(m, x) 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, for three values of n = 0. In this example we just take w = w(r = 0) and use the present-day AGSS09ph Standard Solar Model 2 [59,68]. Here f i represents the fractional composition of each species with atomic number A i . For the spin-dependent cross-sections, the sum over the species i 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 (p/q 0 ) 2n 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 σ ∝ v 2n rel , the modification to the capture rate is much simpler. The partial capture rate Ω(w) ∝ σ is simply transformed by the replacement 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 (n = 0) treatment and an average of the cross-section throughout the volume of the star where w 0 ≡ v 2 0 + v 2 esc (r) is the typical velocity acquired by an in-falling WIMP as it reaches radius r. This indicates that a cross-section proportional to v 2 or v 4 -typically thought of as a suppression -actually yields a considerable enhancement. We note that σ(v rel ) has been erroneously approximated to σ( w 0 ) rather than σ(w 0 ) in the literature [87], although [v esc (r = 0)/v 0 ] 2n does at least provide the correct enhancement to the capture rate to within a factor of a few.

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 = πR 2 (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 The actual capture rate to be used must therefore be the lesser of Eqs. 2.15 and 2.7. Assuming a steady radius and a local DM density of ρ χ 0.38 GeV cm −3 , the maximum amount of DM that can accumulate in the Sun by its current age is therefore 16) where N b M /m p is the total number of baryons in the Sun.

Conductive energy transport by dark matter
If enough DM is captured by the Sun, its large typical inter-scattering distance l χ 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 [89]. 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 µ ≡ m χ /m nuc 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 l χ (r) is much smaller than both the inverse of the local temperature gradient |∇ ln T (r)| and the DM scale height r χ , the values of α and κ are found by solving the first order expansion of the BCE. This is done in terms of the quantity ε ≡ l χ (r)|∇ ln T (r)|, via the formal inversion of the Boltzmann collision operator C(u, r), where C(u, r)F (u, r) represents the net change in the DM phase space distribution F (u, r) 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 v rel and q must be explicitly included in the calculation. In Ref. [89] Gould and Raffelt computed and tabulated α(µ) and κ(µ) for a constant scattering cross-section (n = 0). In Ref. [88], we extended the formalism to include velocity and momentumdependent cross-sections (n = 0), showing that the resulting changes in α, κ and l χ can have potentially large effects on heat transfer in the Sun.
The equilibrium distribution of DM particles in the gravitational potential φ(r) of the sun is given by [11,88,89] where r = 0 represents the centre of the Sun. The conductive luminosity is: where the factor ζ 2n accounts for a velocity-dependent cross-section, and v T (r) ≡ 2k B T (r)/m χ is related to the typical thermal velocity [89]; v 0 and q 0 are respectively the reference velocity and momentum defined in Eq. 1.1. The rate of energy transported per unit mass of stellar material is: This quantity is usually expressed in units of ergs g −1 s −1 .
If the condition l χ r χ 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. [89,95] by Monte Carlo simulation that Eqs. 3.1 and 3.2 can be corrected to properly account for a large Knudsen number K ≡ l χ /r χ . Here we have formally defined Following [11,50], we define:   Figure 2. Illustration of the transition from the local thermal equilibrium (LTE) to the Knudsen (non-local, isothermal) regime of energy transport by DM scattering, for momentum-dependent (left) and velocity-dependent scattering (right). Leftwards of the peak in each curve corresponds to the Knudsen regime, whereas rightwards is the LTE regime. The total energy transport is plotted for a fixed DM mass (m χ = 10 GeV) and number ratio of DM to baryons (n χ /n b = 10 −15 ). Solid lines are spin-independent couplings, whereas dashed lines represent the spin-dependent case where the DM scatters only on hydrogen. This has the effect of increasing the mean free path, leading to a transition to the Knudsen regime at a much higher value of σ 0 .
where K 0 = 0.4 and τ = 0.5 are empirical values taken from the numerical results of [95]. The DM distribution is then a combination of the isothermal and LTE distributions: where Finally, the Knudsen-corrected luminosity is: This should then be used in place of L χ,LTE 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 for the three cases of momentum-dependent cross-section, with σ 0 = 10 −35 cm 2 . Once again, this is illustrated using a present-day SSM with the AGSS09ph abundances [59,68]. Note that this cross-section leads mainly to transport near the LTE regime. As σ 0 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 σ 0 , for a constant DM-to-baryon ratio n χ /n b = 10 −15 . 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.

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 [96,97]. We took DM routines from the public dark stellar evolution code DarkStars [92], 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. [96] and the modified version of GARSTEC used for this work is the same described in Ref. [97]. Here we just give a summary of the most relevant physical inputs. It includes the nuclear energy generation routine exportenergy.f 3 , updated with the astrophysical factors recommended in the Solar Fusion II [98] compilation. It makes use of the Opacity Project radiative opacities [99], complemented at low temperatures with those from Ref. [100]. The equation of state is the 2005 update of OPAL [101]. Microscopic diffusion of elements, including gravitational settling, thermal and concentration mixing, is treated according to Ref. [102].
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 (τ = 4.57 Gyr) solar luminosity L , solar radius R and the metalto-hydrogen mass fraction (Z/X) . 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. [68], for which (Z/X) = 0.0180. The free parameters in the model are the mixing length parameter α MLT and initial helium and metal mass fractions Y ini and Z ini 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. [68] with Z ini acting as the normalisation factor, and X ini + Y ini + Z ini = 1 by definition (X ini being the initial hydrogen abundance).
In practice, the solar model calibration starts with a homogeneously-mixed 1 M premain 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 5 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 Y ini ). 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 [56].
DarkStars [92] is a Fortran95 package that implements capture, annihilation and energy transport by regular (n = 0) WIMP dark matter in a general stellar evolution code, as described in Refs. [11,39,40,103]. The capture routines were originally adapted from DarkSUSY [104]. The underlying evolutionary code [105] is a Fortran90 rewrite of the the venerable Fortran77 Cambridge STARS package [106][107][108]. 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 n = 0 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.9-2.10, 2.12) rather than just the n = 0 case. We used the conductive transport routines from DarkStars essentially unaltered, except that we included the additional factor of ζ 2n in Eq. 3.2 and utilised the α and κ tables that we computed earlier for n = 0 [88].
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 3 . 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 α(r, t) and κ(r, t) at each height in the star by interpolating in the tables of Ref. [88], as these quantities depend on the specific mixture of plasma species with which the DM particles interact. The code then computes the DM density n χ (r, t) using α(r, t) and κ(r, t), which it uses to determine energy transport. DarkStec then interpolates the resulting values of χ (r) (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, σ ∝ q 2n and σ ∝ v 2n rel with n = {−1, 1, 2}. For each coupling, we simulated one solar model for each point on a grid of m χ and σ 0 . We computed a subset of models using a stringent convergence criterion of one part change in 10 5 for the the solar luminosity, radius and surface metallicity (Z/X), 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 10 3 , 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 ∼ 1.5 CPU years.

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.

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 r CZ , the sound speed profile c s (r), the small frequency separations and the surface helium abundance Y S . 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 Y S . In contrast, the predictions of r CZ , 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 crosssection form and combination of σ 0 and m χ . 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: m χ = 5, 10, 15, 20 and 25 GeV, and each decade in σ 0 between 10 −40 and 10 −30 cm 2 . 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 m χ and σ 0 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 φ ν B of approximately ∼60% 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.

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 σ 0 required to render the Sun opaque to dark matter -and thus to saturate the capture rate Eq. 2.15 -can be lowered by as much as two orders of magnitude for a v 4 rel cross-section, and approximately one order of magnitude in the q −2 and v 2 rel cases. On the other hand, negative powers of v rel and positive powers of a spin-dependent q-coupling yield a large suppression in the capture rate, seen as a   Figure 3. Cross-section normalisation σ 0 at which dark matter capture saturates the geometric limit of Eq. 2.15 for momentum-dependent (left) and velocity-dependent (right) cross-sections. Spinindependent couplings are shown with solid lines, and spin-dependent with dashed lines. The SI q −2 , v 2 and v 4 cases are not shown, as they saturate the capture rate below σ 0 = 10 −40 cm 2 (the smallest cross-section that we simulated). much larger required cross-section to achieve saturation. We illustrate the required value of σ 0 in Fig. 3, based on the output of our simulations. Interestingly, a spin-independent q 2 or q 4 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 q-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.

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 temperaturedependence, the 8 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 8 B neutrino flux for spin-independent and spin-dependent dark matter with velocity and momentum-dependent couplings. Fluxes are normalized to the observed value, φ ν B,obs = 5.0 × 10 6 cm −2 s −1 ; lighter colouring represents a reduced flux. Although the measurement error on the 8 B neutrino flux is only 3% [109], 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 1σ isocontour where the neutrino flux falls below ∼ 85% of the measured value; black lines are the 2σ (71%) contours.
Although the flux of 7 Be neutrinos is not as temperature-sensitive as that of 8 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 7 Be neutrino fluxes to the observed value φ ν Be,obs = 4.82 × 10 9 cm −2 s −1 . Here again we plot 1σ and 2σ contours using the theoretical and observational uncertainties, which are respectively 7% and 5% for 7 Be neutrinos.  Figure 4. The ratio of the predicted 8 B neutrino flux to the measured value φ ν B,obs = 5×10 6 cm −2 s −1 , for each type of spin-independent dark matter coupling defined in Eq. 1.1. In every case the white and black lines show the isocontours where the flux is respectively 1 and 2σ lower than the observed values, based on observational (3%) and modelling (14%) errors, added in quadrature. The cross-sections are normalized such that σ = σ 0 (v/v 0 ) 2n or σ = σ 0 (q/q 0 ) 2n , with v 0 = 220 km s −1 and q 0 = 40 MeV. Simulations carried out in the masked regions did not converge, due to the significant heat conduction by the DM particles, leading in extreme cases to density inversions in the core.   Figure 6. The ratio of the predicted 7 Be neutrino flux to the measured value φ ν Be,obs = 4.82 × 10 9 cm −2 s −1 , for each type of spin-independent dark matter coupling defined in Eq. 1.1. In every case the white and black lines show the isocontours where the flux is respectively 1 and 2σ lower than the observed values, based on observational (5%) and modelling (7%) errors, added in quadrature. The cross-sections are normalized such that σ = σ 0 (v/v 0 ) 2n or σ = σ 0 (q/q 0 ) 2n , with v 0 = 220 km s −1 and q 0 = 40 MeV.  Figure 7. As per Fig. 6, but for spin-dependent couplings.
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 m χ 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 σ 0 , not enough DM can be captured to significantly affect the 8 B production rate. For constant cross-sections, this is around 10 −38 cm 2 in the SI case, and 10 −36 cm 2 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 A 2 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. 4-7 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 v rel 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 [88], yields an overall suppression in the effect relative to the constant case. The opposite behaviour is seen for σ ∝ v −2 rel , although the enhanced transport is not sufficient to compensate for the ∼10 −2 suppression in capture.
Again, the effect of a momentum-dependent cross-section is more subtle. Although the boosted capture rate for a q −2 cross-section moves the region of interest to lower crosssections, 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 q 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 q −2 plots that they can also inhibit energy transport, by reducing the mean free path for conduction.

Depth of the convection zone
The boundary between the convective and radiative zones r CZ is the location at which the radiative and adiabatic temperature gradients are equal. Above this location, the temperature gradient becomes super-adiabatic, hydrostatic equilibrium breaks down and convection sets in. The amount of energy transported by DM is negligible at heights r ∼ r CZ . However, the changes in the density and temperature gradient near the core lead to a small but significant increase in the radiative temperature gradient at much larger radii, causing it to exceed the adiabatic gradient at a slightly lower depth and shift the lower boundary of the convection zone downwards.
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, r CZ, = 0.713 R .  Although the error on the helioseismological inference is just 0.001 R , the theoretical error on the modelled r CZ, is much larger: σ CZ, ,th = 0.004 R . Adding these errors in quadrature, we plot in white and black the contours containing the regions where r CZ, is within 1σ and 2σ, 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 (r CZ, = 0.722 R ) and the depth inferred from helioseismology amounts to almost a 3σ discrepancy. Except for a small region at high q 2 cross-section (σ 0 = 10 −33 cm 2 ), spin-dependent DM struggles to bring r CZ to within much better than 2σ of the inferred value. However, energy conduction from SI interactions does substantially better: v ±2 rel and q 2 models can produce good agreement with the observed value of r CZ, , leading in some regions to agreement at better than 1σ.

Surface helium abundance
The surface helium abundance Y s is another observable that has fallen into disagreement with observations since the revision of the solar composition. Our SSM prediction is Y s = 0.2356, whereas the observed value is 0.2485 ±0.0034. With theoretical errors of ±0.0035 taken into account, this amounts to a discrepancy of 3σ. The addition of dark matter does very little to change Y s , and for most models the discrepancy remains at the 3σ level. For some very few cases (not shown), it actually worsens the discrepancy by up to an additional ∼2σ. This only occurs for models where the fit to sound speed observables is also substantially worsened, notably for large SD q 2 and q 4 cross-sections. The reduction in Y s is caused by a corresponding increase in X i , the initial hydrogen fraction. The increase in X i 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 L .
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.

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 where the sum goes over 5 equally-spaced radial points between r = 0.1 R and 0.67 R . We do not include radii below 0.1 R because the reconstructed sound speed in this region is highly uncertain due to the low number of low-degree (low-l) p-modes reaching the innermost radii of the solar core [110]. 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 [111]. The errors on inversions were taken from [112]. The values of this χ 2 , 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.

Frequency separation ratios
The inverted sound speed profile c s (r) obtained from helioseismological measurements is not very accurate near the solar core because not enough low-l 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 [113]. This is because for radial order n 1 and low angular degree l, 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 [114,115]. In particular, two very useful quantities are the frequency separation ratios σSI ∝ const.  Figure 11. Effective sound-speed χ 2 cs defined in Eq. 5.1, showing the improvement in goodness-of-fit of the radial sound speed profile in solar models with different types of spin-independent ADM. The χ 2 value of the best fit point is indicated by a green line on the colour bar and is shown as a green star in each figure. ∆χ 2 = 2.3 and 6.18 contours, corresponding to 1σ and 2σ deviations from the best fit, are respectively shown in white and black.  Figure 12. Same as Fig. 11, but for spin-dependent couplings.  Figure 13. Small frequency separations r 02 (top) and r 13 (bottom) as defined in Eq. 5.2, for the best fit constant (left) and generalised form-factor dark matter models (right). The latter correspond to the best-fit models of the three couplings returning the best overall p-values (see Tab. 1). Predictions are compared with helioseismological observations from the BiSON experiment [114]. Inner black error bars correspond to observational error, whereas outer (green) bars also include modelling error. Below each figure we show the residuals with respect to BiSON data, in units of the total error. The m χ = 3 GeV, q 2 case with σ 0 = 10 −37 cm 2 yields the best improvement, bringing the largest discrepancy from nearly 4σ to little more than a standard deviation. and ∆ l (n) ≡ ν n,l − ν n−1,l . The sound speed gradient is weighted by 1/r in the integral above, so r 02 (n) and r 13 (n) 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 r 02 (n) and r 13 (n) measured by the BiSON experiment [115]. 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 q 2 scattering with m χ = 3 GeV, σ 0 = 10 −37 cm 2 . In this model, the largest discrepancy in r 02 and r 13 falls to barely 1σ.
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 χ 2 r 02 + χ 2 r 13 , where .
Formally, there is a correlation between values of r 02 (n) and r 13 (n) 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 χ 2 r 02 + χ 2 r 13 . 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 q 2 SI coupling, which also clearly produces the best agreement with the observed frequency separations.

Combined limits
It is clear from the results we have presented in Secs. 5.2-5.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 s 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: where φ ν B,obs = 5.00 × 10 6 cm −2 s −1 , φ ν Be,obs = 4.82 × 10 9 cm −2 s −1 , r CZ,obs = 0.713 R and Y S,obs = 0.2485. 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 r , +2 , 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.  Figure 14. Combined likelihood χ 2 r02 + χ 2 r13 for the small frequency separations defined in Eq. 5.2, for different models of ADM with spin-independent scattering cross sections. The green line on the colour bar shows the χ 2 value of the best fit, indicated by a green star in each panel, whereas white, black and cyan lines show the preferred regions at 1, 2 and 4 σ confidence levels, respectively. At low masses, many DM models can yield improvements in the overall chi-squared with respect to the Standard Solar Model, but the q 2 model provides a clear best fit.    Figure 17. Same as Fig. 16, but for spin-dependent couplings. Table 1. Standard Solar Model (SSM) and best fit (b.f.) values for each of the models we consider, along with observable quantities. The DM mass and cross-sections are in GeV and cm 2 , respectively. The boron-8 neutrino flux φ ν B is in units of 10 −6 cm −2 s −1 and the beryllium-7 neutrino flux φ ν Be is expressed in 10 −9 cm −2 s −1 . The full chi-squared is defined in Eq. 5.5 and includes the neutrino fluxes, surface helium abundance Y S , depth of the convection zone and small frequency separations. 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 χ 2 and implied p-value. The best-fit regular SI and SD cases shown in Table 1 differ from the best-fit cases in Ref. [72] because here we actually choose them on the basis of the full likelihood Eq. 5.5. In Ref. [72] 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 crosssections, including the frequency separations shifts the best-fit masses and cross-sections. This is to be compared to the q 2 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.

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 q. This is confirmed to very high significance by the small frequency separations, seen in Fig. 14. At low masses m χ ∼ 3 − 5 GeV, an asymmetric DM candidate with σ = σ 0 (q/40MeV) 2 and σ 0 = 10 −37 cm 2 leads to a substantial improvement with respect to the SSM. The best fit occurs at 3 GeV, and leads to a convective zone depth r CZ = 0.718R , only 1.2σ away from the inferred value from helioseismology. The 8 B and 7 Be fluxes are respectively 3.78 × 10 6 cm −2 s −1 and 4.29 × 10 9 cm −2 s −1 , both within 2σ of the measured values. Surface helium remains low, at Y S = 0.233. The large improvement to the fit to c s (r) 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. [72].
In the next section, we will discuss constraints from other sectors, and show that the best-fit q 2 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 q-dependent models, as well as constant and velocity-dependent models we examined, nearly all led to poor overall fits: the effect on φ ν and Y S 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 p-values greater than 10 −3 . As pointed out in Ref. [87], a v −2 rel SD cross-section does indeed lead to an improved fit to the sound speed profile and r CZ 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 3σ reduction in the 8 B neutrino fluxes, as well as significant reduction in 7 Be neutrinos. In fact, we find that a v 2 rel SD cross-section yields a much more plausible overall fit to observables than v −2 rel SD scattering, resulting in a p-value of 3.1 ×10 −4 as compared to 4.2 ×10 −10 (vs. 0.85 for the q 2 SI model). It is also important to note that the region where v −2 rel 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. [52], 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.

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 q and v rel -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 (q n scattering, m χ ∼ 3-5 GeV, σ SI ∼ 10 −37 -10 −38 cm 2 ) are some of the few that survive the direct search limits. For velocity-dependent and constant couplings, the curves lie almost entirely below the plotted regions, indicating that in these cases, the entire section of parameter space that turns out to have an impact on solar physics is disfavoured.
It is worth remembering, however, that limits from direct detection and solar physics probe very different parts of the dark matter halo velocity distribution: solar capture preferentially occurs from the low-velocity end of the distribution, whereas direct detection probes the high-velocity tail. The compatibility of solar and direct searches is therefore highly sensitive to the chosen velocity distribution, and might be modified substantially by choosing something other than the standard Maxwell-Boltzmann distribution.
Limits on both momentum-dependent and velocity-dependent DM-quark couplings also exist from colliders [116][117][118][119], although in this case a specific effective operator or UVcomplete model must be assumed in order to calculate definite event rates. In the classification scheme of [120]( [116]), operators F2(D2) and F10(D10) have leading contributions to the nuclear recoil rate from q 2 SI terms, operator F6(D6) from v 2 rel SI terms, operators F3(D3), F6(D6) and S2(C2/R2) from q 2 SD terms, operator F4(D4) from q 4 SD terms, and operator S4(C4) from v 2 rel SD terms. For the most part, ATLAS and CMS have not focused on such operators due to the weakness of the limits on them from direct detection. Limits from Tevatron data have been calculated by CDF [121] on the combination of F4(D4) and the 'regular' constant SI scalar couplingχχqq (F1/D1), but so far only a few phenomenological groups [116,119] have placed limits on the other operators using up-to-date collider data.
From the point of view of solar physics, given the preference we see for q 2 SI scattering in ameliorating the Solar Abundance Problem in Fig. 16, operators F2(D2) and F10(D10) would seem to be the most interesting. Existing limits on F10(D10) from ATLAS data [119] are quite stringent at m χ below a few hundred GeV, excluding mass scales of below a TeV for the physics associated with the particle mediating the interaction. Such limits would likely make achieving our preferred reference cross-section rather difficult. Limits on F2(D2) are much weaker though, with only mass scales below 42 GeV excluded at m χ < 50 GeV (assuming that the contact operator approximation even remains valid in this regime). It therefore seems feasible that a DM particle interacting with quarks by this effective coupling (χγ 5 χqq) might explain the recent anomalies seen in some direct detection experiments, and help rectify the known discrepancies between theory and data in solar physics.

Comparison with previous work
The effect of DM heat conduction in the Sun, especially by asymmetric dark matter, has been studied in several works. Most recently, the main qualitative features of such models were pointed out by [122] while accurate computations using solar models were done in Refs. [51][52][53], with somewhat contradictory conclusions. In this section we will endeavour to settle the differences. Among these studies, only [52] used the proper thermal conduction mechanism of Gould and Raffelt [89]. Our energy transport and 8 B neutrino flux computations are in good agreement with that work. However, Ref [52] did not include the (small) effect of molecular settling, nor did they account for the uncertainty in the initial helium mass fraction Y , metallicity Z or mixing length parameter α, all of which are free parameters of solar models. Allowing Y , Z and α to vary has the effect of partially compensating for the increased luminosity from DM transport, thus reducing the overall effect of DM and leading us to weaker overall constraints.
In contrast, Refs. [51,53] used the approach developed by Spergel and Press (SP, [123]). However, comparison with Monte Carlo simulations found this approach to be flawed [95]. This is due to the violation of all three assumptions that go into the SP approach: 1) a Maxwellian velocity distribution; 2) spatial uniformity of the WIMP velocity distribution; and 3) local isotropy. To further quantify the errors introduced by the use of the SP treatment, we have implemented the SP transport mechanism into our code, following Ref. [49] for a constant, spin-dependent (SD) cross-section. We plot the total energy transport for a given DM population in the Sun as a function of the spin-dependent cross-section in the lefthand panel of Fig. 18. Solid lines represent the accurate Gould and Raffelt (GR) approach (which we use in this paper), whereas dashed lines are the SP calculation. As in Fig. 2, the peak energy transport occurs at the transition between local (LTE) and non-local (Knudsen) transport. Although both approaches roughly yield the same location for this peak, SP slightly underestimates transport in the LTE regime, whereas it significantly overpredicts transport in the non-local regime. In the large cross-section LTE regime, the SP case predicts DM luminosities that are 1.5 (for m χ = 5 GeV) to 4 (m χ = 80 GeV) times too small; in the Knudsen limit, this becomes an overestimate by a factor of 3 (m χ = 80 GeV) to 9 (m χ = 5 GeV). It is worth noting that the largest discrepancies happen very near the location of interest, where low masses and cross-sections potentially lead to effects in the Sun and simultaneously evade constraints from direct detection. In the right-hand panel of Fig. 18 we show the resulting discrepancy between predicted 8 B neutrino fluxes with SP and GR transport, as computed with our solar model. The red region is where SP overpredicts an effect, whereas blue corresponds to the LTE regime where SP underestimates the energy conduction.
Two other recent papers [86,87] have also investigated the impacts of specific dark matter models with non-constant cross-sections on solar physics. In both these papers, the authors employed existing capture and transport schemes (SP in [86] and GR in [87]) without properly accounting for the dependence of the capture and transport processes on the kinematic structure of the cross-section. Ref. [86] examined a magnetic dipole dark matter model with a cross-section that scales as v 2 rel q −2 relative to the standard constant case. 4 They then applied the 'momentum transfer cross-section', which amounts to multiplying the cross-section by (1 − cos θ) in order to cancel the divergence induced by forward scattering when integrating over all scattering angles to obtain a total cross section. This process is only useful for obtaining a rough scattering rate, as it simply transmutes a momentumdependence into an equivalent v rel -dependence (see e.g. calculations of κ for q −2 couplings in [88] for details). The correct treatment is of course to include the full momentum-and velocity-dependent differential cross-section in the integral over the DM velocity distribution function. In the case of Ref. [86], this treatment simply converted a v 2 rel q −2 cross-section to a v 2 rel v −2 rel = constant cross-section, effectively ignoring the momentum and velocity-dependence of the model. Similarly in Ref. [87], where the application of the momentum transfer crosssection treatment converted what was actually a v 2 rel q −4 cross-section in the long-range limit into a v 2 rel v −4 rel = v −2 rel cross-section, relative to the standard case. The authors then took a thermal average over the velocity distribution and applied it as a constant scaling to the standard capture and transport treatments. As we have shown previously [88], this is valid  [89,95], overpredicts the transported energy by as much as a factor of 9 in the non-local (Knudsen) regime, while underestimating it in the LTE case. Right: comparison between predicted 8 B fluxes once these transport mechanisms are implemented into the full DarkStec solar capture code. In the red region, using SP yields a larger reduction in the neutrino flux than GR.
for the capture rate (or would be, if the cross-section were truly dependent only on v rel and not q), but qualitatively incorrect when dealing with conductive energy transport.
In light of these shortcomings in the way capture and transport were modelled, and in the theoretical treatment of the scattering cross-sections involved, it is difficult to estimate how the results of Refs. [86,87] would change were the kinematics of the cross-sections treated correctly. Calculating the transfer coefficients as per Ref. [88] and carrying out the full phenomenological analysis as we have here, for v 2 rel q −2 and v 2 rel q −4 couplings, would therefore be an interesting topic for future study.

Conclusions
In this paper, we have explored the prospects for solving the Solar Abundance Problem using generalised form factor dark matter, confirming that it can indeed provide a solution to this long-standing issue [72]. In the process, we derived the full capture rates in the Sun for dark matter with momentum-or velocity-dependent interactions with nucleons, and used solar physics to place limits on those couplings. This is the first extensive, physically-consistent analysis of the effects on the Sun of momentum-or velocity-dependent dark matter, taking into account the full capture kinematics, along with the self-consistent thermal conduction framework developed in Ref. [88]. We incorporated this framework into a new state-of-the art dark solar evolution code, DarkStec.
Our best-fit model, asymmetric DM of m χ ∼ 3 GeV and σ 0 = 10 −37 cm 2 interacting with nucleons proportionally to q 2 , shows excellent promise for solving the Solar Abundance Problem. Our model provides a remarkable fit, and a full investigation of the effects of evaporation, which has not been done before for non-constant cross-sections, should serve as a cross-check on the consistency of such a scenario. The prospects for detecting such a particle at the next round of direct detection experiments, and in the next run of the LHC, also appear very promising.