Updated constraints on velocity and momentum-dependent asymmetric dark matter

We present updated constraints on dark matter models with momentum-dependent or velocity-dependent interactions with nuclei, based on direct detection and solar physics. We improve our previous treatment of energy transport in the solar interior by dark matter scattering, leading to significant changes in fits to many observables. Based on solar physics alone, DM with a spin-independent $q^{4}$ coupling provides the best fit to data, and a statistically satisfactory solution to the solar abundance problem. Once direct detection limits are accounted for however, the best solution is spin-dependent $v^2$ scattering with a reference cross-section of 10$^{-35}$ cm$^2$ (at a reference velocity of $v_0=220$ km s$^{-1}$), and a dark matter mass of about 5 GeV.


Introduction
Despite tremendous success in predicting neutrino fluxes and describing the bulk structure of the Sun, the Standard Solar Model (SSM) still fails to reproduce key observables relating to helioseismology. After the downward revision of solar photospheric abundances over 10 years ago [1][2][3][4][5][6][7][8][9][10][11][12] it has become clear that the predicted sound speed profile, position of the base of the convection zone and surface helium abundance are all several standard deviations away from the values obtained by direct helioseismological inversion [13][14][15][16][17][18]. Despite over a decade of effort, no adequate solution has been found to reconcile solar modelling with these precision observables [19][20][21][22], prompting a search for solutions beyond the Standard Model.
As the Sun travels through the dark matter halo of the Milky Way, it inevitably interacts with the dark matter population 1 . If a DM particle scatters with a nucleus in the Sun to a velocity less than the local escape velocity, it becomes gravitationally bound. From there, it will quickly settle into an equilibrium orbit near the solar centre, governed by the thermodynamics of the weakly-interacting DM gas in the steep gravitational potential of the Sun. The effects of dark matter capture in stars has been studied in depth since the 1980s . In the absence of self-annihilation (i.e. asymmetric dark matter, ADM [62,63]), a large population of dark matter can accumulate. This can act as a heat conductor, acquiring kinetic energy from the hot core, and releasing it via interactions with nuclei in the cooler outer regions. This can lead to a slightly shallower temperature gradient with height in the star. Despite the small population of DM (at most ∼1 particle per 10 10 baryons), these small adjustments in the thermal gradient can have measurable effects on our own Sun. These include the solar structure itself -including the sound speed c s (r) and convective zone radius r CZ -and on neutrino fluxes from fusion processes, due to their strong dependence on the core temperature.
Precision solar models that include capture and heat transport from standard spindependent or spin-independent ADM can be built to satisfy the solar radius, age and luminosity [64][65][66]. Whilst these models can bring some helioseismological observables into better agreement with data than the SSM, this comes at the cost of an elastic scattering cross-section that is several orders of magnitude higher than allowed by direct detection experiments. These models also lead to drastic underproduction of solar neutrinos.
However, the kinetic regime probed by solar capture is very different from earth-based direct detection experiments: the Sun preferentially captures slower-moving particles from the DM halo, as these are more likely to scatter to sub-escape velocities, whereas direct detection relies on large enough velocities to produce observable recoils in a target material. This observation is partly what led to efforts to examine DM models with less trivial interactions with the Standard Model [67][68][69][70][71]. Such interactions are generic in particle physics: scattering mediated by a massive particle generally gives rise to non-relativistic elastic scattering cross section that is proportional to some positive power of the momentum transfer q or relative velocity v (see e.g. [72]). Conversely, new long-range forces can yield cross sections proportional to negative powers of these quantities (e.g. [73]).
Recently, we found [70,71] that a light asymmetric particle with a momentum-dependent interaction with quarks could be captured in large enough quantities in the Sun and conduct heat in such a way that helioseismic observables could be brought into excellent agreement with data. Although this led to a reduction in the predicted neutrino fluxes, the significant improvement in sound speed and convective zone depth was sufficient to produce a marked improvement over the SSM of 6 standard deviations and a potential solution to the solar composition problem. New searches by direct detection experiments at low scattering threshold have since ruled out the best fit from these studies. CRESST-II [74] analysed the specific model highlighted in [70,71] excluding it to high significance. The even more sensitive CDMSlite analysis [75] released by the SuperCDMS collaboration around the same period confirms this result.
We have furthermore revised the formalism of [67], and found that the conductive luminosity had been over (under) estimated in models with a cross section proportional to a positive (negative) power of q or v. We have corrected this error, finding that the best fit point in each model has moved in the parameter space. The best fits that we find are as good as those we found in [70,71], though the required cross sections are higher, in some cases by as much as an order of magnitude. Given the strength of the direct detection bounds that we derive in this work, this is of little consequence as the previously-favoured parameter space is now entirely ruled out.
Our goals in this work are therefore 1) to update our solar simulations, correcting the error in luminosity, and 2) to confront these results with constraints from direct detection, to determine whether such models can indeed lead to a solution to the solar composition problem. We find that the parameter space that remains, after taking into account the recent CDMSlite results, is highly restricted. Although the models that obey these bounds do not fit the solar data as well as some that are excluded by direct detection, they are nonetheless capable of producing very significant improvements over the SSM (see also [69] for similar results). This paper is structured as follows. In Sec. 2 we outline the class of models that we are considering, and present the relevant equations for capture by the Sun and heat transport within it. In Sec. 3 we briefly discuss constraints from direct detection on light DM particles able to affect solar structure. We present the results of our simulations in Sec. 4, and conclude in Sec. 5.

Dark matter in the Sun
The method of obtaining the capture rate and subsequent energy transport due to momentum or velocity-dependent DM in the Sun is presented in detail in Refs. [67,71]. Here, we show the main results, along with a crucial correction to the transported luminosity formula, Eq. 2.12. Reflecting the nature of most concrete models so far proposed in the theory literature, we focus on models with isoscalar (identical proton and neutron) DM-nucleon couplings with the Standard Model. We parameterise the resulting differential cross-sections as Each of these can be spin-independent (SI), coupling coherently to all nucleons; or spindependent (SD), coupling only to the spin of unpaired nucleons. This "form factor" approach [76], encompasses both non-relativistic effective operators that can come from effective pointlike interactions [77][78][79], as well as certain classes of long-range forces (e.g. [73]). The cross section for scattering with a nucleus is then

Capture in the Sun
In the absence of annihilation or evaporation, the population of DM in the Sun is simply given byṄ χ (t) = C (t): Here, f (u) is the DM speed distribution in the Sun's frame, and Ω(w) encodes the kinematics of scattering below the local escape velocity v esc (r). w(r) ≡ u 2 + v 2 esc is the local DM velocity inside the star's gravitational potential. Expressing the DM-to-nucleon mass ratio as µ i ≡ m χ /m N i , this quantity is: where GF F I is the generalized form factor integral, and n i (r) is the number density of each nuclear species in the Sun. In the case of a constant cross section, the GF F I is simply an integral over the nuclear form factor F i (E R ): where the integral limits are set by the kinematics required to downscatter the DM velocity enough for it to be captured. For velocity-dependent scattering, the required modification of the capture rate is straightforward: the integrand in Eq. 2.3 is simply multiplied by an overall factor of [w(r)/v 0 ] 2n .
When the cross section is momentum-dependent, the factors of (q/q 0 ) 2n must be included inside the GF F I. For scattering with hydrogen (assuming a point-like particle, i.e.
where p = m χ w. For heavier nuclei, a parameterisation of the form factor must be chosen. Significant work over the past few years has gone into computing accurate form factors that model the nuclear response function using effective field theory; these yield interactions close to those in Eqs. 2.1 and Ref. [80] has computed these response functions for capture in the Sun. The standard Helm [81] parameterisation is where E i is a constant quantity for each nuclear species i, given by GeV. (2.8) For spin-independent capture and transport, we use the following 15 elements: H, He, C, N, O, Ne, Mg, Na, Al, Si, S, Ar, Ca, Fe and Ni. Though heavier elements are suppressed in abundance, the A 2 coherence factor in Eq. 2.2 can enhance the capture rate off these elements by as much as four orders of magnitude.
We have furthermore checked in Ref. [71] that for most models, the deviations with respect to Helm are minimal, as low-momentum scattering tends to dominate the capture rate. There is one exception: when computing spin-dependent capture rates, we neglect elements heavier than hydrogen. These do not benefit from the A 2 coherent enhancement present for spin-independent interactions, and their very low abundance further suppresses their contribution. However, in the q 4 SD case, which can be matched to an O 6 ≡ ( S χ · q/m N )( S N · q/m N ) coupling, the total capture rate in an evolved solar model is actually dominated by scattering with nitrogen. This does not turn out to be very helpful, as we will find that the cross-section required to produce an effect in the Sun with this operator is ruled out by five orders of magnitude. We thus use the Helm form factor for every case, and the generalised form factor integral for momentum-dependent interactions becomes: where B ≡ 1 2 m χ w 2 /E i , and Γ(m, x) is the (upper) incomplete gamma function. The Sun cannot capture more dark matter than the so-called geometric limit, in which all of the DM that intercepts the solar disk is captured. The total capture rate is therefore the smallest of Eq. 2.3 and where u 0 = 270 km s −1 is the dispersion of the Maxwell-Boltzmann DM velocity distribution, u = 220 km s −1 is the velocity of the Sun relative to the DM halo, and ρ χ = 0.38 GeV cm −3 . is the local density of dark matter. We show the cross section that saturates this bound in Fig. 1, for each of the models we consider in this paper. Strictly speaking, the transition from Eq. 2.3 to 2.10 should be smooth, rather than a sharp cutoff as we impose here. To accurately model this transition, the full optical depth must be including in Eq. 2.3. This leads to a suppression of the capture rate as σ 0 approaches the saturation cross section, with the ultimate effect of slightly suppressing the effect of DM close to this threshold. As this computation is relatively involved, we address it in an upcoming work [82]. Spin-independent Spin-dependent Figure 1. Cross section at which the capture rate computed in DarkStec saturates the geometric limit 2.10. Left: spin-independent interactions, coupling to every element; right: spin-dependent interactions with hydrogen only. Despite their small abundances, heavier elements contribute significantly to spin-independent capture, thanks in large part to the coherent σ N ∝ A 2 enhancement.

Conductive energy transport by dark matter
In the local thermal equilibrium (LTE) regime, where the mean DM interscattering distance l χ is much smaller than the scale of the DM distribution r χ , the equilibrium distribution of DM particles in the gravitational potential φ(r) of the Sun is given by [55,67,83] n χ,LTE (r) = n χ,LTE (0) T (r) where r = 0 represents the centre of the Sun. The conductive luminosity is: where we have removed the erroneous factor ζ 2n which was present in [67,71]. Below, we will find that this has the effect of displacing the best fit regions in parameter space, to higher (lower) values of σ 0 for positive (negative) values of n. v T (r) ≡ 2k B T (r)/m χ is related to the typical thermal velocity [83]; v 0 and q 0 are respectively the reference velocity and momentum defined in Eq. 2.1. The effect of our correction to Eq. 2.12 is to suppress transport by models with positive powers of q and v, and enhance it for negative powers. The rate of energy transported per unit mass of stellar material is: This quantity is usually expressed in units of erg g −1 s −1 . The molecular diffusion coefficient α(r, µ) and conduction coefficient κ(r, µ) are computed by perturbatively solving a Boltzmann collision equation in the diffuse gas limit. This method was introduced by Gould and Raffelt [83] and generalised in Ref. [67]; as before, we use the tabulated coefficients given in the latter reference.
In the non-LTE, "Knudsen" regime K ≡ l χ /r χ 1, the DM distribution becomes isothermal: Gould and Raffelt [83] found that the transition to the LTE regime could be well described by the interpolating functions [55,84]: and Energy transport has the opposite behaviour as a function of the interaction cross section in the LTE and Knudsen regimes. In the LTE regime (large σ 0 ), a higher scattering rate suppresses the typical distance travelled, confining DM and suppressing energy transport. Conversely, in the Knudsen (small σ 0 ) regime, the inter-scattering distance is already large, and energy transfer is instead dominated by the collisional efficiency: luminosity thus increases with increasing cross section. Heat conduction is therefore maximised at the boundary between these regimes, when inter-scattering distances are large, but interaction rates are still large enough to allow efficient energy transfer. We show this behaviour in Fig. 2. It follows that the largest effect of momentum or velocity-dependent dark matter on solar observables will occur when the cross-section is closely matched to the Knudsen peak (as long as a sufficient amount of DM can be captured with this cross section).
In Fig. 3 we show the equivalent for the impact of DM mass, illustrating the strong enhancement of energy transport at low masses due to mass-matching between DM and hydrogen or helium. This corresponds to transitions in α and κ at m χ ∼ m N for different interactions [67], and can also be seen in the fact that transport is maximised at lower masses for SD scattering, reflecting the fact that H plays a role due to its spin, but He does not. The combination of this effect with the transition from LTE to non-local transport can be seen in Fig. 4. Whilst this figure shows contours for momentum and velocity-independent scattering only, other interactions lead to a very similar pattern, just scaled or shifted in the x, y and/or z directions.
Finally, as the DM mass is lowered below 5 GeV, it may pick up sufficient speed from thermal collisions to escape the sun. This can lead to significant loss through evaporation.
However, the exact amount of evaporation depends sensitively on the mass and cross section. Since this dependence is not not known, we also explore this region of the parameter space.
In an upcoming work [82] we aim to compute these precise evaporation masses, for every model explored here.

Direct detection bounds
Recently, the SuperCDMS collaboration has presented their lowest-threshold analysis ever, [75]. This 70 kg day run was sensitive to recoil energies above 56 eVee (or around 300 eV nuclear recoil energies), allowing them to set exclusion bounds on DM masses extending below 2 GeV. Though CRESST-II [74] have presented analyses of specific models in this work, we use CDMSlite data because it is both more sensitive and easier to recast into new limits. Since CDMS uses a germanium crystal target, it presents the further advantage of being sensitive to spin-dependent interactions: 7.6% of natural Ge is in the form of 73 Ge (J = 9/2), allowing spin-dependent bounds to be derived with the same data set. In contrast, the CaWO 4 target used by CRESST-II presents no appreciable target mass with nonzero spin.
Assuming different SI and SD couplings, but equal couplings to protons and neutrons for each type of interaction (i.e. isoscalar couplings), the differential nuclear recoil rate is where µ p ∼ µ n is the DM-nucleon reduced mass and v min ≡ m N E R /2µ 2 N is the minimum DM velocity required to produce a recoil with energy E R . In this expression, µ N is the DMnucleus reduced mass. Using the binned data, efficiency curve and ionisation yield presented in Ref. [75] and a simple Poisson analysis, along with the same halo parameters as in our solar analysis, we obtain the exclusion curves presented in Fig. 5. For SI limits (left panel) we use the same form factor as for solar capture. For the spin-dependent limits (right panel), we assume scattering only on 7.6% of the target, and use the spin expectations and isoscalar structure function S 00 from Ref. [85], such that Although we show limits for isoscalar couplings only, it is worth noting that most of the unpaired spin contribution from 73 Ge comes from neutrons, so limits on ADM that couples mainly to protons would be significantly weaker. The SI limits in Fig. 5 are slightly more constraining than the limits presented by CRESST-II [74], for the spin-independent q 2 and q 4 cases. Spin-independent Spin-independent   Spin-independent

-28
Spin-dependent Figure 5. CDMSlite upper limits on σ 0 for the models that we consider, from the event rates and instrument details given in Ref. [75]. To compare with our Solar results, we use the same values of v 0 = 220 km/s and q 0 = 40 MeV. We also show results obtained by the CRESSST-II low-threshold analysis taken from [74], which are slightly less constraining than CDMSlite for DM masses above 1.4 GeV.

Solar bounds
To produce our updated constraints, we employ the DarkStec code, developed for and described in Ref. [71]. DarkStec is a combination of the legendary GARSTEC solar evolution code [86,87], and the lightweight capture and transport routines from DarkStars [52,55,88]. These were modified to include the full capture and transport technology described above, for momentum and velocity dependent dark matter. To obtain our results, we perform a scan over a grid of masses m χ = {3, 5, 10, 15, 20, 25} GeV and one cross section per decade from σ 0 = 10 −40 to 10 −30 cm 2 , totalling over 900 simulations, or approximately 2.5 CPU years. A full description of the observables that we use is given in Ref. [71]. Here we summarise the salient points. Solar neutrino fluxes and small frequency separations are highly sensitive to the impacts of DM on the core of the Sun, and typically provide the strongest constraints on DM models. The depth of the convection zone is sensitive to impacts of DM on the temperature gradient closer to the surface, and generally depends more subtly on the overall diffusivity of DM in the Sun. The surface helium abundance and radial profile of the sound speed profile probe the entire radiative region, but contribute comparatively little additional constraining power over other observables.
Plots of the neutrino fluxes, depth of the convection zone, sound-speed and smallseparation likelihoods vs m χ and σ 0 are presented in Appendix A. We also provide the full results of all our simulations, including each of these observables, their likelihoods and the combined likelihood presented in Sec. 4.1, as machine-readable supplementary material.
Solar neutrino fluxes: Minute changes in the core temperature can have substantial effects on the solar neutrino fluxes, due to their steep temperature-dependence. The main fusion process in the Sun, p + p → 2 H + e + + ν e , is quite insensitive to changes due to DM energy transport, as the luminosity condition enforces a fairly constant reaction rate. However, the subdominant flux of neutrinos from the decay of 8 B produced in the pp and pep chains is an especially good probe of the core temperature, as are the two lines at E ν = 862 keV and 384 keV from the decay of 7 Be. While these are measured at 3% and 5% accuracy, respectively, the uncertainties in their production rates are closer to 14% and 7%. We add these errors in quadrature when assessing the goodness of fit. As the SSM is in very good agreement with the measured neutrino fluxes, this serves as a constraint on the total effect DM can have in the solar core. The changes in neutrino fluxes in different models are shown in Figs. 10-13. When uncertainties are compared on the same footing, our constant, SD results are in agreement with [64].
Depth of the convection zone: The radius at which hydrostatic equilibrium breaks down, radiative pressure is no longer balanced by gravity, and convective energy transport begins, depends crucially on the radial temperature and pressure gradients. The height of the convection zone can be inferred from helioseismological measurements, R CZ, = 0.713 ± 0.001 R , while the SSM predicts a much higher value: R CZ,SSM = 0.722 ± 0.004 R . As DM deposits energy at larger radii, the absolute temperature gradient increases slightly, thus lowering the location of the base of the convection zone. This is shown in Figs 14 and 15.
Surface helium abundance: As the initial mass fractions of hydrogen, helium and metals cannot be directly measured, they (along with the mixing length parameter α M LT ) are allowed to vary as DarkStec attempts to find a solution for the input parameters. The present-day surface helium abundance is thus predicted by these parameters, when combined with the physics entering into the SSM. The measured value is Y s = 0.2485 ± 0.0034, while the SSM yields Y s = 0.2356 ± 0.0035. The surface helium abundance essentially reflects the initial helium abundance, which is strongly constrained by the requirement that the model must reproduce the observed solar luminosity. Any change in initial helium abundance must come with a corresponding change in the initial hydrogen fraction, which modifies the resulting solar luminosity unless the core temperature is also drastically altered. In extreme cases, DM can have enough of an effect on the core temperature to change the resulting surface helium abundance expressed by a model -so we do include this observable in our fits -but typically, the effects on the core temperature are not drastic enough for this observable to be significantly affected.
Sound speed: The radial sound speed profile can be inferred by inverting the oscillation frequencies of the different angular modes projected onto the solar surface. Fig. 6 shows the sound speed profile obtained with our SSM (blue solid line), when compared with the errors from helioseismological measurement and inversions (green regions) and modelling (blue regions). The change in the temperature, pressure and density due to the redistribution of heat by DM can be quite large -as long as the cross section is large enough to saturate the capture rate and hit the Knudsen transition. Fig. 6 also shows the effect of three models which are allowed by DD constraints, along with one which is not (SI, q 4 , solid magenta line), but which illustrates the strong effect DM can have. We also show in Figs. 16 and 17 the goodness of fit of each of our models to the sound speed profile from helioseismology observations with SoHO and BiSON [89], using a chi-squared with 5 equally-spaced points between R = 0.1R and 0.67R , where helioseismic errors are minimal. Given the large modelling errors, however, and the fact that c s is highly correlated with the more precise frequency separation ratios, we do not include the sound speed in our total χ 2 .
Frequency separation ratios: It is possible to remove a large fraction of the systematic errors that enter sound speed inversions near the core by directly using specific combination of frequencies. Two most commonly used quantities are the so-called small frequency separation ratios [90]: where ∆ l (n) ≡ ν n,l − ν n−1,l and These have the additional advantage of being insensitive to details in the structure of the external layers of the Sun, which are not properly modeled in SSMs. The 1/r dependence ensures sensitivity to the properties of the solar core, and the modeling and observational uncertainties on r 02 and r 13 are much smaller than for Solar neutrinos. These frequency separation ratios are shown in Fig. 7 for the SSM (red), and several models including energy transport by ADM.

Combined constraints
We define a combined chi-squared, including all of the quantities listed above except the sound speed profile: where φ ν B and φ ν Be correspond to the boron-8 and beryllium-7 neutrino fluxes, r CZ is the depth of the convection zone, Y S is the surface helium abundance, and r ij are the small frequency separation ratios. The sound speed profiles are strongly correlated with the small separations, and less precise, so we do not use them. The subscript "obs" indicates the observed value. The errors σ i include both modelling and observational error, added in quadrature.
The results of these combined constraints are given in Fig. 8 for spin-independent couplings, and Fig. 9 in the spin-dependent case. The thick magenta lines in these figures are the constraints that we obtained above from the latest CDMSlite data. There is very little parameter space in which DM is unconstrained by direct detection, yet has a strong enough effect on the Sun to change the fit.
The best-fit points for each model are given in Table 1. In all cases, the SI direct detection limits are too strong to allow any of our best-fit points. In only two spin-dependent cases (v 2 and v 4 ) are the best-fit points below the CDMSLite limits. In one other case (SD q −2 ), there is a region that is almost as good as the best fit; we list this point instead. These allowed models are labelled with a dagger ( †) in Table 1.
Despite the low p-values, each of the allowed models still represents a remarkable improvement over the SSM. This improvement is driven by the much better fits to the small Modelling error Helioseismology error No DM SI, const., m = 15 GeV, < 0 = 10 !37 cm 2 SD, q 2 , m = 3, < 0 = 10 !39 cm 2 SD, v 2 , m = 5, < 0 = 10 !35 cm 2 SI, q 4 , m = 3, < 0 = 10 !32 cm 2 Figure 6. Deviation of radial sound speed profile from values inferred from combined SoHO and BiSON helioseismology data [89] using our SSM. We show the best fit for a constant cross section in black, and the two best fits that are also allowed by CDMSlite data, the spin-dependent v 2 and q 2 models (indicated by daggers in Table 1). We finally show the best overall fit, a spin-independent q 4 DM model. However, this model is ruled out by direct detection data by many orders of magnitude. The dark and light blue bands represent the 1 and 2σ errors on modelling, estimated from error propagation of uncertainties in SSM inputs; the green bands represent the 1 and 2σ errors on helioseismological inversions. frequency separations, assisted by a drop in the base of the convection zone. The former observation tells us that the main improvement comes from better modelling of the region near the solar core. However, the increased tension with measured neutrino fluxes suggests that the decrease in core temperature predicted by models of ADM in the Sun may be too strong. A full DM solution to the solar abundance problem must therefore strike a delicate balance between softer thermal gradients and a lower core temperature. We end this section by noting that we also provide as supplementary material online a complete table of the boron and beryllium neutrino fluxes, surface helium, radius of the convection zone and small frequency separations for all of our models, in addition to the partial and total chi-squared values defined in Eq. 4.3. This table can be used for quick lookup and interpolation, e.g. for global fits of DM models.  Figure 7. Small frequency separations r 02 (left) and r 13 (right) as defined in Eq. 4.1, for various dark matter models. DM models correspond to the best-fit parameters for two couplings returning the best overall p-values (see Table 1), without being excluded by CDMSlite, along with the best constant case (black), and the overall best fit case (pink). The latter two are however excluded by direct detection.
Predictions are compared with helioseismological observations from the BiSON experiment [91]. 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. Though the v 2 model with m χ = 5 GeV yields substantial improvement in r 13 for all frequencies, it does less well at low frequencies ( 2500 Hz) with respect to the measured r 02 .

Conclusions
We have presented new limits on momentum and velocity-dependent dark matter from direct detection, and an updated study of their effects on state of the art simulations of the Sun. A correction of the LTE conduction formalism (Eq. 2.12), has shifted the best fits in the parameter space, and the very strong limits from the low-threshold analysis of CDMSlite have ruled out all but a small part of the parameter space of models that can improve the SSM. Our solver could not find solutions for certain values in the parameter space, ostensibly because they led to changes in the solar structure that were too large with respect to the SSM (which it was originally designed to compute). This is not to say that these models are ruled out, however. We note that in some cases (notably spin-dependent v −2 ), the non-convergence region lies on top of the region that is probed by direct detection. A better exploration of this parameter space could be very interesting indeed.
Of the models that we have found that provide a good fit, only some can be obtained from the current paradigm of simplified models of dark matter. For instance, a spin-dependent v 4 cross section does not occur in such setups. A fermion DM candidate interacting via a vector mediator can give rise to an SD, v 2 cross section, however. This requires a purely vector coupling to the DM, and an axial coupling at the quark level. However, this setup also gives a similar contribution proportional to q 2 (see e.g. [93]), which would be three orders of magnitude larger using our normalisations of q 0 and v 0 , and thus excluded by experiment. This is not to say that this model cannot be obtained, only that it is not straightforward using current tools. A q −2 cross section, on the other hand, is closer to a long-range force, which is again not covered in the simplified model approach. A model that can yield a similar type of behaviour would be an electromagnetic dipole or anapole interaction. These models actually couple with mixed powers of v and q, and thus require explicit recalculation of the molecular diffusion and conduction coefficients α and κ, and is the subject of a dedicated study [94].
Finally, we note that the effects of evaporation have not been included here, as the full evaporation formalism must be generalised to the form factor case from scratch. Because the relevant DM masses are so close to the evaporation limit for constant cross-sections (m evap ∼ 4 GeV), this may have a an important impact on the remaining parameter space. Full evaporation rates for momentum-and velocity-dependent interactions will be computed in detail in an upcoming publication [82].
It is clear from the fits we present here that the addition of heat transport by dark matter in the Sun can improve the SSM significantly. The improvement is well beyond what one would naively expect from adding just two degrees of freedom. Even if DM is not the solution, the fact that the true solution can be accurately mimicked by DM would seem to indicate something fundamental about the physical processes behind the solution. Given the types of DM interactions that work, it is clear that if DM is to explain the solar composition problem, then a strong model-building effort is required to connect our results with a more robust model of dark matter that can be embedded in a UV-complete theory.

A Contour plots of solar observables
In this appendix we show contour plots for individual quantities impacted by the presence of dark matter in the Sun. In order, we show the boron-8 and beryllium-7 neutrino fluxes, the depth of the convection zone and the goodness of fit (χ 2 ) computed for the sound speed profile and for the small frequency separations r 02 and r 13 . Figure 10. The ratio of the predicted 8 B neutrino flux to the measured value φ ν B,obs = 5.00 × 10 6 cm −2 s −1 [92], for each type of spin-independent dark matter coupling defined in Eq. 2.1. In every case the white and green 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. Be,obs = 4.82 × 10 9 cm −2 s −1 [92], for each type of spin-independent dark matter coupling defined in Eq. 2.1. In every case the white and green 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. Regions masked in light blue correspond to parameter combinations where models did not converge. Figure 14. Ratio between the modelled and measured location of the bottom of the convection zone r CZ , for spin-independent couplings. Darker regions represent a better fit to the observed value than the SSM. The white and green lines represent the contours at which the predicted value falls within 1σ and 2σ of the measured value, respectively. The theoretical uncertainty on r CZ (0.004 R ) is much larger than the experimental error (0.001 R ), so the former dominates when we add them in quadrature. Regions masked in light blue correspond to parameter combinations where models did not converge.