Dark Matter Thermalization in Neutron Stars
Abstract
We study how manybody effects alter the dark matter (DM) thermalization time inside neutron stars. We find that Pauli blocking, kinematic constraints, and superfluidity and superconductivity in the neutron star significantly affect the DM thermalization time, in general lengthening it. This could change the final DM mass and DMnucleon cross section constraints by considering black hole formation in neutron stars due to DM accretion. We consider the class of models in which DM is an asymmetric, complex scalar particle with a mass between 1 keV and 5 GeV which couples to regular matter via some heavy vector boson. Interestingly, we find that the discovery of asymmetric, bosonic DM could motivate the existence of exotic neutron star cores. We apply our results to the case of mixed sneutrino DM.
pacs:
95.35.+d, 95.30.Cq, 97.60.JdI Introduction
Over the past 80 years there has been increasing evidence–from galactic rotation curves, the Bullet cluster, the anisotropies in the cosmic microwave background, the distribution of galaxies, etc.–for the existence of cold, nonbaryonic dark matter (DM). While many constraints have been placed on DM, its precise nature remains unknown (see Bertone et al. (2005); Bergstrom (2000); Feng (2010) for a review). In particular, the mass of the DM particle is highly unconstrained. For elementary particle DM, we have that
(1) 
where the lower bound comes from bosonic DM being confined on galaxy scales ( kpc) and the upper bound is the Planck scale. Phase space density bounds fermionic DM to be heavier than several keV Tremaine and Gunn (1979). Other constraints on DM apply, even for bosons, if the galaxy core/cusp problem is explained by warm DM Boyanovsky et al. (2008). DM candidates with GeV must be either black holes, e.g. primordial black holes Ivanov et al. (1994); Carr et al. (2010); Capela et al. (2013); Griest et al. (2013), or extended objects like composite DM or solitons, e.g. Qballs Kusenko and Shaposhnikov (1998); Kusenko et al. (2005).
Currently, direct detection experiments most strongly constrain the DMnucleon cross section for DM with a mass 50 GeV. The best constraint comes from XENON100 which requires the DMnucleon cross section to be for 50 GeV Aprile et al. (2012). Lighter DM candidates with are less constrained by such experiments due to large backgrounds. For this paper we will focus on the lessconstrained parameter space of DM particles with mass between 1 keV and 5 GeV. In this region, astrophysical observations can be used to constrain DM mass and cross section parameter space which is difficult for direct detection experiments to probe.
Such astrophysical observations were initially used to constrain hypothetical weakly interacting particles by studying their interactions inside of planets and stars Press and Spergel (1985a, b); Goldman and Nussinov (1989), and in the past few years there has been a renewed interest in using these methods to place constraints on DM by considering neutron star accretion of DM particles Gould et al. (1990); Bertone and Fairbairn (2008); Kouvaris (2008); de Lavallaz and Fairbairn (2010); Kouvaris and Tinyakov (2010, 2011a, 2011b, 2012); McDermott et al. (2012); Kouvaris (2012); Guver et al. (2012); Bramante et al. (2013); Bell et al. (2013). Under certain conditions this accretion of DM results in a black hole that then destroys the neutron star. Observationally we know that the oldest neutrons stars are 10 billion years old Manchester et al. (2005), hence the DM parameter space which allows for neutron star destruction in less than 10 billion years is ruled out.
In what follows, we only consider bosonic DM since, due to the absence of Fermi pressure, it becomes gravitationally unstable for fewer accumulated particles than fermionic DM. We also only consider asymmetric DM (for recent reviews see Zurek (2013); Petraki and Volkas (2013)), in which there is an initial asymmetry between particles and antiparticles so that today only particles remain and DM annihilation in the neutron star can be ignored. Coannihilation of DM with other particles in the neutron star was considered in Bell et al. (2013) and will be neglected in our analysis. Using black hole formation in neutron stars due to the accretion of asymmetric, bosonic DM, references Kouvaris and Tinyakov (2011a, b, 2012); McDermott et al. (2012); de Lavallaz and Fairbairn (2010); Guver et al. (2012); Bramante et al. (2013); Bell et al. (2013) have put otherwise modelindependent constraints on DM. Here we show that an improved treatment of particle kinematics and manybody effects including Pauliblocking, superfluidity, and superconductivity can have a significant effect on the DM thermalization time, changing the final DM constraint. We also show that in some phases of high density matter, such as color superconducting quark matter, thermalization times can be surprisingly large.
This paper is organized as follows. In section II, we briefly review the nature of the DM constraint–the process of DM capture by a neutron star, black hole formation, and destruction of the neutron star. In section III, we discuss the effective theory for our generic DM model in light of this scenario. In section IV, we define the thermalization time and compare our results to those previously obtained. We conclude in section V and discuss mixed sneutrino DM in appendix A and B as an example of a DM model which can be constrained in this way.
Ii Dark Matter Capture and Black Hole Formation
Here we review this process as previously discussed by many others (e.g. Kouvaris and Tinyakov (2011a); McDermott et al. (2012)). There are number of steps involved. First, because of its gravitational interactions, DM is accreted by the neutron star. We can estimate the velocity of the incident DM particle at the surface of the neutron star by using classical energy conservation:
(2) 
where and is the particle’s velocity infinitely far from the neutron star. We will take . For a typical neutron star, kg and km. Using these standard values, we find
(3) 
This implies that the energy that a typical DM particle has at the surface of a neutron star is
(4) 
so we see that the incident DM energy is set by its mass and that typical DM particles are semirelativistic. These incident DM particles will scatter with quasiparticles inside the neutron star, lose energy, and become bound to the star.
Next DM thermalizes inside the neutron star. Since the incident DM particle is at most semirelativistic, and it must lose energy in order to be captured by the neutron star, it is safe to assume that the typical DM particle is nonrelativistic during the latter collisions that determine its thermalization time. As the DM thermalizes, it collects within a sphere of radius which satisfies
(5) 
where is the mass of the neutron star enclosed within a radius and is the temperature of the neutron star. We can estimate this by considering a neutron star with a constant core density GeV/fm and we find Kouvaris and Tinyakov (2010)
(6) 
This tiny sphere of DM at the center of the neutron star can begin to selfgravitate and collapse into a black hole. Gravitational collapse is accelerated if the captured DM forms a BoseEinstein condensate inside the star Kouvaris and Tinyakov (2011b); Jamison (2013). Once the black hole is formed, it must be massive enough to avoid evaporation due to Hawking radiation and then it may consume the neutron star. The observational signatures of a neutron star collapsing into a black hole is still an interesting, open question.
In previous works, Kouvaris and Tinyakov (2011a); McDermott et al. (2012), two calculations to constrain the DMneutron cross section as a function of DM mass are done: 1) the thermalization time calculation: years, in which is the time necessary for DM thermalization with the neutron star and 2) an accretion time calculation: years, in which is the time needed for the neutron star to accrete enough DM to form a black hole which will destroy the star, assuming thermalization occurs in a negligible amount of time. The second calculation sets the final constraint, and the first is used to find regions where the second constraint is not valid. In this paper we will consider only the first calculation and its application to the particular class of DM models to be discussed next.
Iii Dark Matter Model
We consider a model in which DM is a complex scalar particle which couples to regular matter by exchanging some heavy spin one boson. The effective Lagrangian for the interaction between DM and the fermions (nucleons, electrons, etc.) that are found in neutron stars is then given by
(7) 
where is the DM current, and are the vector and axialvector currents for the fermions, and is the coupling of to the mediator divided by the coupling of to the mediator. For simplicity we take to be the standard model value for fermions coupling to the Z boson. is the coupling constant after the heavy mediator has been integrated out. In general,
(8) 
where is the mass of the heavy mediator particle, is the coupling of the mediator to , and is the coupling of the mediator to .
Use of this effective theory is welljustified. In order for the effective theory to capture the relevant physics, one needs that the magnitude of the fourmomentum transfer squared, , is much less than in the DMfermion scattering processes. Based on the arguments in the previous section we know that the initial DM energy is at most and hence the maximum that the DM can give up is . Since the fermions inside the neutron star are highly degenerate, scattering events in which the DM gains energy and momentum from them are rare and have typical for the DM masses we are considering. Note that we will always take keV as it was shown in Kouvaris and Tinyakov (2011b) that for keV one must worry about captured DM escaping the star. Hence as long as , then and our effective theory is valid. For this reason we will take GeV, consistent with the assumption that the heavy vector mediator is either a standard model or pair, or a heavier, undiscovered particle.
Unless forbidden by some symmetry, this effective theory also includes DM selfinteractions. These have been shown to affect the critical number of DM particles needed for black hole formation Bramante et al. (2013); Bell et al. (2013), but are not relevant for thermalization time calculations and hence will not be discussed further here.
Iv Thermalization Time Calculation and Results
iv.1 The Thermalization Time
For scattering between DM and fermions, let be the initial DM fourmomentum, be the final DM fourmomentum, be the initial fermion fourmomentum, and be the final fermion fourmomentum. We define the thermalization time as the average time it takes for an incident DM particle to start having collisions in which the average energy transfer is less than the temperature of the neutron star, i.e. where is the zeroth component of the fourmomentum transfer . Note that we will assume that the DM particles are confined to the neutron star interior. While initially DM particles which have become bound to the neutron star may have an orbit that goes outside the star Kouvaris and Tinyakov (2011a), unless DM is extremely light ( keV), the DM particles are confined to the neutron star interior for later stages of cooling. Since we are considering the oldest, coldest neutron stars, we take eV.
To derive a formula for the thermalization time, we will make use of the DM scattering rate, . Using Fermi’s Golden Rule, the scattering rate for DM scattering with a medium of spin fermions is given by
(9)  
where is the amplitude for the process, is the fermion distibution function (FermiDirac for noninteracting fermions), and is the DM distribution function (BoseEinstein for noninteracting bosonic DM). We will neglect the Bose enhancement factor for the final DM state for simplicity since the distribution function for the DM particles is a complicated function of time due to the accumulation of DM–note that this means that (9) as is, is actually a lower bound on the scattering rate.
The tree level squared amplitude, , is averaged over initial and summed over final fermion spins and is given by
(10)  
(11) 
where we have used the notation with the mostly minus metric and is the fermion mass which could be the neutron or the electron mass, or respectively. Using finite temperature formalism, the scattering rate can also be expressed as Kapusta and Gale (2006); Doniach and Sondheimer (1974); Reddy et al. (1998)
(12) 
where contains the DM currents and is the fermion retarded polarization tensor. For noninteracting fermions, these are given by
(13a)  
(13b) 
where is the free fermion propagator at finite temperature and density. The form for this polarization tensor has been worked out in detail in Reddy et al. (1998) and Saito et al. (1989) and we use their results in our calculations. (If using the derivation in Reddy et al. (1998) note Reddy (2013).)
The polarization tensor, , characterizes the medium response to the DM probe. The fermion propagators contain the Pauli blocking factors (c.f. the factor of in (9)) which restrict the fermion phase space due to the Pauli exclusion principle, i.e. the incident fermion that interacts with the DM particle must come from the initial fermion distribution and the scattered fermion must occupy phase space that is not already filled by the initial fermion distribution. The polarization tensor also contains information about the inmedium fermionfermion interactions since is a fermion currentcurrent correlation function which includes a sum over all possible intermediate states.
Given an expression for the scattering rate ((9) or (12)), we can now define a discretized version of the thermalization time, , based on the physical reasoning that the average thermalization time is simply the sum of the average times for subsequent DM collisions until the average energy transfer per collision is less than the temperature of the neutron star. Thus we may write
(14) 
where is the initial DM energy, which we will always estimate to be (note that this assumes a decrease in initial DM velocity due to prior collisions necessary for DM capture) and for is the average final energy of a DM particle which had initial energy . The final energy is determined by calculating the scattering rate for a DM particle with initial energy weighted by the final DM energy, and dividing by the unweighted scattering rate for a DM particle with initial energy , i.e.
(15) 
The summation in (14) ends once . We expect that this generally results in . Expression (14) is used for all of our numerical work.
We also define an approximate, continuous version of the thermalization time as
(16) 
We use (16) for our analytic results, where instead of finding as described above, is fit to a value that approximates the numerical result well.
Now that the thermalization time is well defined, we can calculate the DM thermalization time inside a neutron star. This includes DM scattering with a liquid of neutrons, protons, and electrons, a neutron superfluid, and a proton superconductor. It also includes DM scattering with the matter in the core of the neutron star–possibly hyperons, pion or kaon condensates, quark gluon plasma, etc. For a review of the constituents of a neutron star see Glendenning (1997); Lattimer and Prakash (2004); Page and Reddy (2006); Page et al. (2011). The majority of the neutron star (roughly 85%) is made up of neutrons, so we will first consider DM thermalization by scattering with neutrons, in both the normal (Fermi gas phase) and the superfluid phase.
iv.2 Scattering with a Fermi Gas of Neutrons
From nucleonnucleon scattering data we know that the neutronneutron interaction can be either attractive or repulsive depending on the spin and spatial angular momentum of the neutrons and on the neutron density Takatsuka and Tamagaki (1993). At sufficiently low temperature, attractive interactions can lead to superfluidity and dramatically alter the lowlying excitation spectrum, and hence the DM scattering mechanism, and we discuss this in detail in the next section. Here, to calculate DMneutron scattering, we will ignore nuclear interactions and approximate the neutrons as a dense, noninteracting Fermi gas. This will provide a baseline result since we know from Fermiliquid theory that corrections due to strong interactions in the normal phase do not qualitatively change the nature of scattering or the kinematics Bardeen et al. (1957). From earlier work, relating to neutrino scattering in dense, normal neutron star matter Burrows et al. (2006), we expect that the DM scattering rates in the Fermi gas approximation are sufficient to provide a reliable order of magnitude estimate.
The fiducial calculation is done for neutrons at saturation density () which corresponds to a nonrelativistic neutron chemical potential of GeV. This implies that neutrons at saturation density are to a good approximation, nonrelativistic. Deep in the core, neutrons become mildly relativistic and but these relativistic corrections are modest. We calculate the thermalization time and then enforce years, which gives a constraint of the form as a function of . We then use this constrained coupling constant in the formula for the DMfermion cross section in the limit in which both the DM and fermion momenta tend to zero (a good approximation of what takes place in direct detection experiments):
(17) 
This gives the DMfermion cross section as a function of DM mass alone, with the constraint that DM thermalization takes longer than years.
For noninteracting neutrons it is simplest to use expression (9) for the scattering rate in the calculation of the thermalization time. Eqn. (9) was used for numerical calculations and an approximate analytic result was obtained as follows. For thermalization time scatterings it is a good approximation that both the neutrons and DM are nonrelativistic, so neglecting all momentum dependence in the amplitude in (10) and rewriting the scattering rate we find
(18) 
where is the fourmomentum transfer and . is the neutron response function, here given by
(19) 
where is the FermiDirac distribution function. Additionally in the limit of completely degenerate neutron matter (in reality so the neutrons really are quite degenerate) and for , we have Reddy et al. (1998):
(20) 
where , is the Heaviside step function, and is the neutron Fermi velocity.
Note that the step function is just enforcing nonrelativistic, low momentum transfer neutron kinematics, i.e. that . That this inequality holds can be seen simply from
(21) 
where is the angle between and . These neutron kinematics must be consistent with the same nonrelativistic, low momentum transfer DM kinematics () and since always by construction, the DM kinematics constrain the phase space more and the neutron step function in (20) can simply be set to 1. These kinematics are shown in Fig. 1.
Using (20) in (18), setting to zero as the thermalization time definition always has and completing the angular integrals gives
(22) 
Since the neutrons are approximated as completely degenerate, DM cannot lose energy to them, hence , and using , we find
(23) 
We can now use this to calculate the denominator in (16):
(24) 
Integration gives
(25) 
Using this in (16) we find
(26) 
Setting , using to match to numerical calculations, and enforcing years gives the final result for . Our numerical and analytic results are shown in Fig. 2 along with previous results for comparison. Note that the result shown from McDermott et al. (2012) is their full analytic result and not the approximation that they plot in their figures.
In order to compare with analytic expressions from previous works McDermott et al. (2012); Kouvaris and Tinyakov (2010), we neglect with respect to in (26) and insert (17) into the expression to obtain
(27) 
where . To get a feel for typical scales, this can be recast as
(28) 
which is generically longer than previous calculations by several orders of magnitude.
From Fig. 2 one can see that the results obtained here differ appreciably from those in previous works–in particular some regions of DM parameter space that were disallowed in McDermott et al. (2012) are allowed from this calculation due to an increase in thermalization times. This is because the proper inclusion of kinematics and Pauli blocking are essential for calculating the low energy and momentum transfer scattering processes that lead to thermalization. In the past Pauli blocking has been included only roughly, and kinematic effects have been neglected. To estimate the contribution of Pauli blocking and kinematic constraints to our calculations we define an effective suppression factor , given by
(29) 
where is the classical expression for the scattering rate, is the number density of the neutrons, is the cross section given in (17), and is the magnitude of the relative velocity between the incident DM and incident neutron. (see (9)) is the actual scattering rate which is an integrated version of . Using (17), (23) after integrating over , and using a thermal for the incident DM momentum, we find
(30) 
In Fig. 3 we compare our suppression factor to the Pauli blocking suppression factor used in McDermott et al. (2012).
Fig. 3 shows that the inclusion of Pauli blocking and kinematic effects in the properly integrated scattering rate makes a difference. These scattering kinematics and Pauli blocking are unique to DM scattering with a noninteracting, nonrelativistic Fermi gas and will change once interactions between the fermions are included, especially in the case of attractive interactions which can give rise to superfluidity or superconductivity. We discuss these effects next.
iv.3 Scattering with Superfluid Neutrons
From BCS theory we know that attractive interactions at the Fermi surface leads to the formation of Cooper pairs at low temperature and results in a phase transition to either a superfluid or superconducting state Bardeen et al. (1957). In this superfluid or superconducting state, there is a nonzero ground state expectation value (or condensate) of Cooper pairs which produces a gap in the fermion excitation spectrum and a Goldstone boson due to the spontaneous breaking of the symmetry associated with fermion number Bardeen et al. (1957); Weinberg (1986).
For neutrons in the core of the neutron star the dominant attractive interaction is in the pwave channel and is expected to lead to the formation of spintriplet Cooper pairs Page et al. (2011). Model calculations predict that the energy gap, , is roughly MeV, though this remains somewhat uncertain Schwenk and Friman (2004). The condensate of these pairs is expected to be spatially anisotropic and Goldstone bosons associated with the breaking of rotational invariance arise in addition to the Goldstone boson from the spontaneous breaking of fermion number Bedaque et al. (2003). Since in our model, DM couples only to the neutron density in the nonrelativistic limit (in (13a), ), the only relevant excitation at energies small compared to the gap is the Goldstone boson, or superfluid phonon, associated the breaking of the fermion number symmetry.
The superfluid phonon manifests as spikes in the densitydensity neutron response function (, c.f. (13b)) at , where is the speed of the superfluid phonon in the nuclear medium. Based on this, we can make an ansatz for the neutron response function:
(31) 
where and are normalization constants which can be fixed by enforcing the principle of detailed balance and the fsum rule Plischke and Bergersen (2006)
(32) 
This gives
(33) 
In our calculations we will only use the part of which allows the DM to lose energy. The part is also suppressed for thermalization scatterings with . We take
(34) 
Note that this response function only characterizes DM emission of a single phonon. Multiphonon processes are quite suppressed and will be discussed in section IV.5. This response function can be used in place of the neutron part of the scattering rate in (9) and then the DM thermalization time in a neutron superfluid can be computed. In doing so, we varied between and . The value is the leading order speed of the superfluid phonon Son and Wingate (2006). In general , where is the pressure in the nuclear medium and is the energy density; is expected to vary inside the neutron star due to phonon interactions and as a function of density.
We find in general that DM particles will only scatter with the neutron superfluid once or twice, leaving the DM with too much energy to be considered thermal. This result is due to the highly restricted kinematics of the neutron superfluid, see Fig. 4. Since the single phonon mode can only respond with , once the DM particle loses enough energy such that , neutron superfluid and DM kinematics are no longer compatible and no further scattering can occur.
Since DM cannot thermalize by single phonon emission in the neutron superfluid, the DM particle must scatter with something else inside the neutron star in order to thermalize. In addition to suppressed multiphonon scattering, the other standard options are scattering with protons or electrons. The protons are likely to be in a superconducting state with Cooperpaired protons and a massless, coupled proton and electron mode Baldo and Ducoin (2011). The paired protons have a gapped energy spectrum and hence do not contribute much to DM thermalization. The massless protonelectron mode has kinematics similar to that of the neutron superfluid phonon mode, so it also does not allow DM to thermalize. This only leaves the electrons to thermalize the DM.
iv.4 Scattering with a Fermi Gas of Electrons
Electrons in a neutron star have a vanishingly small critical temperature for pairing, so the lowenergy spectrum of particlehole excitations is ungapped and well described by that of a noninteracting Fermi gas. Thus electronDM scattering can be treated in the same way as the neutronDM scattering in section IV.2. Roughly of a neutron star is made up of electrons and for neutrons at saturation density, electrons have a chemical potential of GeV, indicating that the electrons are highly relativistic with the Fermi velocity . Since DM is nonrelativistic, its dominant coupling is to the electron density but the kinematics differs qualitatively from the neutron case because and DMelectron scattering is always kinematically allowed inside the neutron star.
The electron response function to leading order in the velocity of the DM particle is given by
(35)  
(36) 
where is the electron FermiDirac distribution function and is the angle between and . In Fig. 5 we show the numerical results obtained from setting years (using (12),(14), and (17)) for the low energy DMelectron cross section as a function of DM mass. The DMneutron cross section results are plotted for comparison.
Interestingly, if DM couples with equal strength to neutrons and electrons (i.e. if is fixed), then we find that thermalization times for DM scattering with electrons are roughly 50% of thermalization times for DM scattering with neutrons, so regardless of the presence of a superfluid, DMelectron scattering would be the most efficient process for DM thermalization.
iv.5 Scattering in Exotic Neutron Star Cores
So far we have considered DM thermalization with electrons and also neutrons, both in the normal phase and in the superfluid phase. However, the phase structure of matter in the neutron star core remains uncertain Page and Reddy (2006). In this section we study two specific phases of high density matter in order to explore their influence on DM thermalization. At asymptotically large densities where the strange quark mass can be considered small and perturbation theory is applicable, it is now well established on theoretical grounds that the ground state of quark matter is the color flavor locked (CFL) phase in which the approximate symmetry of QCD is spontaneously broken down to its vector subgroup due to the formation of a condensate of diquark pairs Alford et al. (1999, 2008). This is a color superconducting phase in which all nine (3 flavors 3 colors) light quarks form Cooper pairs and there is a gap in the particlehole excitation spectrum.
At densities of relevance to neutron stars, the strange quark mass is dynamically important, perturbation theory fails, and whether the CFL phase is present at these densities remains an open question. If it were present, the CFL phase could additionally contain a condensate of K mesons Bedaque and Schäfer (2002); Kaplan and Reddy (2002). Both the CFL and the CFLK phases are characterized by similar low energy properties at temperatures of relevance to old neutron stars. They are both devoid of electrons and the only relevant low energy degrees of freedom are the massless Goldstone bosons associated with the breaking of global symmetries in the ground state Alford et al. (1999, 2008). There is one massless phonon mode, sometimes called the boson, due to the breaking of the symmetry in the CFL phase and two massless phonon modes in the CFLK phase, one due to the breaking of (called ) and another due to the breaking of the hypercharge symmetry by the condensate (called ). The velocity of the mode in the relativistic limit is approximately given by and the velocity of the is where is proportional to the number density of the kaon condensate Kaplan and Reddy (2002).
Earlier we found that the gap in the nucleon spectrum due to pairing implied that DM thermalization would proceed via superfluid phonon emission processes as long as , where was the speed of the Goldstone mode in the nuclear medium. For , this process is kinematically forbidden and electron scattering dominates. In the CFL and CFLK phases electrons are absent and relevant DM thermalization processes can only involve the massless Goldstone bosons. As in the superfluid nuclear phase, thermalization in the CFL and CFLK phases proceeds by the phonon emission process shown by the second diagram in Fig. 6 (akin to the Cherenkov radiation of fast particles) as long as if DM couples to the baryon number current and if it also couples to the hypercharge current. When , DM thermalization cannot proceed by phonon emission and the dominant thermalization process is the two phonon process shown by the third diagram in Fig. 6. Here, the initial state phonon is thermal with energy , and the intermediate phonon is offshell.
For simplicity we will consider DM that couples only to the the baryon number. In this case, the low energy effective Lagrangian has the form where is half of the overall phase of the condensate that breaks the symmetry and where MeV is the quark chemical potential in the neutron star core Son and Stephanov (2000). Then the amplitude for the DM + phonon DM + phonon process can be approximated by
(37) 
where is the dimensionless constant that sets the strength of the leading order three phonon vertex in the low energy theory, is the energy of the intermediate phonon, is the magnitude of the momentum of the intermediate phonon, and and are the initial and final phonon energies respectively. We used an approximate form of the phonon propagator, , since the phonon must be offshell with . Using (37) the scattering rate can be calculated for arbitrary initial DM velocities and to first order in . Taking typical values for our parameters: , , , and , we find that the dimensional estimate (ignoring factors of and ) for the DMphonon scattering rate is given by
(38) 
Using , and in (38) for a single DMphonon scattering process, and estimating we find that the DM thermalization time is approximately
(39) 
indicating that the DM scattering rate is too low to allow for thermalization even for the oldest neutron stars with ages years if the core contains either the CFL or CFLK phase.
V Conclusions
We considered a relatively generic DM model in which the DM particle was a complex scalar that coupled to regular matter via some heavy vector boson, with the regular matter vector and axialvector couplings to the heavy mediator taken to be those to the Z boson. We then calculated DM thermalization times inside a neutron star for DM scattering with electrons, neutrons in the normal phase, neutrons in the superfluid phase, and color superconducting quarks.
We found several important results:

Including kinematics in the thermalization time calculation resulted in DM thermalization times that were qualitatively different from past results.

Previously neglected DMelectron scattering in ordinary neutron star cores is actually quite important. It is a more efficent DM thermalization mechanism than DMneutron scattering when the neutrons are in a nonsuperfluid state (for a fixed ) and it is the only relevant DM thermalization mechanism when the neutrons form a superfluid and the protons form a superconductor.

Exotic neutron star cores with color superconducting quark matter and no electrons give rise to very large thermalization times which protects neutron stars from their possible destruction as a result of DM accretion. Hence the discovery of asymmetric, bosonic DM could motivate the existence of exotic neutron star cores.
Acknowledgements.
This work was supported in part by the U.S. Department of Energy Grants DEFG0200ER41132 (SR and BB) and DEFG0296ER40956 (AEN and BB). BB would also like to acknowledge partial support from an ARCS Foundation Fellowship.Appendix A Nonthermal Dark Matter
Light DM candidates such as mixed sneutrinos can be produced nonthermally via the AffleckDine mechanism. The AffleckDine mechanism was originally associated with baryogensis Affleck and Dine (1985) but it can be applied to any scalar field that can have a large vev and whose interactions are negligible. In a cosmological context, the end result of this mechanism is to nonthermally produce a large number of nearly zeromomentum particles–exactly what is needed for cold DM.
To see this explicitly, consider the Lagrangian density for a complex scalar field with mass in curved space time with metric ,
(40) 
In the following we neglect the interaction terms, i.e. we take , and we assume that inflation smooths out any spatial dependence in our field so that . Working under the assumption that the universe is flat, isotropic, and homogeneous, we use the FriedmannRobertsonWalker metric with scale factor . Thus the classical equation of motion for our field is
(41) 
where dots indicate time derivatives and is the Hubble parameter.
This is a damped harmonic oscillator equation, so that for (which occurs early in the universe) is overdamped and hence is frozen at its initial vev, . For late times with , oscillates (only in time) with its natural frequency and we find that the energy density of the field, , is proportional to and scales like just like the energy density for highly nonrelativistic matter:
(42) 
For very nonrelativistic matter, , where is the number density of particles, thus we also have that so that for the proper choice of the initial vev , we can match ’s number density today to what we observe for DM, i.e.
(43) 
where is the time today, is the gravitational constant, is the Hubble parameter today, and Ade et al. (2013). We also note that the DM masses considered in this paper are large enough so that to match the energy density today, the vev is small enough such that DM is present (beginning roughly at the time that satisfies ) before matter domination, as required by observations of the CMB.
It is easy to explain why interaction terms for are negligible, as well as why a macroscopically large vev can form during the early universe if can be nonzero along a flat direction in the scalar potential. Such flat directions are common in supersymmetric theories; for example, it is possible to find flat directions with combinations of squarks, sleptons, and Higgs fields in the Minimal Supersymmetric Standard Model Gherghetta et al. (1996).
Appendix B Light Sneutrino Dark Matter
In this appendix we consider AffleckDine produced, supersymmetric DM. Specifically we take the DM candidate to be the lightest mass eigenstate of some linear combination of active sneutrinos, and an additional sterile sneutrino, . Such mixed sneutrino DM is discussed, for example, in ArkaniHamed et al. (2001); Hooper et al. (2005); Gopalakrishna et al. (2006); Belanger et al. (2010). This DM particle carries lepton number of +1 and based on the initial vev for the DM field, there can be an initial asymmetry between the number of DM particles and antiparticles, making DM annihilation negligible today.
Since the DM field, , is a superposition of and we can write
(44)  
(45) 
where is the heavier mass eigenstate. Note that for , is unconstrained by the invisible Z width Hebbeker (1999).
is a weak isosinglet with only gravitational interactions and is in a weak doublet. We take the dominant interactions of to be with the weak gauge bosons (i.e. its couplings via the superpotential are negligible and its coupling to is kinematically suppressed). We can find ’s weak interactions with gauge bosons by using the covariant derivative from the symmetry of the standard model, so that kinetic terms for and in the Lagrangian are given by
(46) 
where and . The relation of the and gauge bosons to the standard model photon , the , and the bosons is

(47) 
where and , where is the magnitude of the electron charge and is the weak mixing angle. Replacing and with and , we find that the interaction Lagrangian for is given by
(48)  
If the the fourmomentum transfer squared in a DM interaction is less than a few GeV, then the first term in the above interaction Lagrangian is dominant. We may also integrate out the in the remaining term, giving us the generic effective Lagrangian form that we had in (7). In this case, the effective coupling constant is given by
(49) 
where for a neutron and for an electron and all the results of the previous sections can be applied.
References
 G. Bertone, D. Hooper, and J. Silk, Physics Reports 405, 279 (2005), ISSN 03701573, URL http://www.sciencedirect.com/science/article/pii/S03701573040%03515.
 L. Bergstrom, Rept.Prog.Phys. 63, 793 (2000), eprint hepph/0002126.
 J. L. Feng, Ann.Rev.Astron.Astrophys. 48, 495 (2010), eprint 1003.0904.
 S. Tremaine and J. E. Gunn, Phys. Rev. Lett. 42, 407 (1979), URL http://link.aps.org/doi/10.1103/PhysRevLett.42.407.
 D. Boyanovsky, H. de Vega, and N. Sanchez, Phys.Rev. D77, 043518 (2008), eprint 0710.5180.
 P. Ivanov, P. Naselsky, and I. Novikov, Phys.Rev. D50, 7173 (1994).
 B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys.Rev. D81, 104019 (2010), eprint 0912.5297.
 F. Capela, M. Pshirkov, and P. Tinyakov, Phys.Rev. D87, 123524 (2013), eprint 1301.4984.
 K. Griest, A. M. Cieplak, and M. J. Lehner (2013), eprint 1307.5798.
 A. Kusenko and M. E. Shaposhnikov, Phys.Lett. B418, 46 (1998), eprint hepph/9709492.
 A. Kusenko, L. C. Loveridge, and M. Shaposhnikov, JCAP 0508, 011 (2005), eprint astroph/0507225.
 E. Aprile et al. (XENON100 Collaboration), Phys. Rev. Lett. 109, 181301 (2012), URL http://link.aps.org/doi/10.1103/PhysRevLett.109.181301.
 W. H. Press and D. N. Spergel, Astrophys. J. 294, 663 (1985a).
 W. H. Press and D. N. Spergel, Astrophys. J. 296, 679 (1985b).
 I. Goldman and S. Nussinov, Phys. Rev. D 40, 3221 (1989), URL http://link.aps.org/doi/10.1103/PhysRevD.40.3221.
 A. Gould, B. T. Draine, R. W. Romani, and S. Nussinov, Phys.Lett. B238, 337 (1990).
 C. Kouvaris and P. Tinyakov, Physical Review D 83, 083512 (2011a).
 C. Kouvaris and P. Tinyakov, Phys.Rev.Lett. 107, 091301 (2011b), eprint 1104.0382.
 C. Kouvaris and P. Tinyakov, arXiv preprint arXiv:1212.4075 (2012).
 S. D. McDermott, H.B. Yu, and K. M. Zurek, Physical Review D 85, 023519 (2012).
 A. de Lavallaz and M. Fairbairn, Phys.Rev. D81, 123521 (2010), eprint 1004.0629.
 C. Kouvaris, Phys.Rev. D77, 023006 (2008), eprint 0708.2362.
 C. Kouvaris and P. Tinyakov, Phys.Rev. D82, 063531 (2010), eprint 1004.0586.
 G. Bertone and M. Fairbairn, Phys.Rev. D77, 043515 (2008), eprint 0709.1485.
 C. Kouvaris, Phys.Rev.Lett. 108, 191301 (2012), eprint 1111.4364.
 J. Bramante, K. Fukushima, and J. Kumar, Phys.Rev. D87, 055012 (2013), eprint 1301.0036.
 N. F. Bell, A. Melatos, and K. Petraki, Phys.Rev. D87, 123507 (2013), eprint 1301.6811.
 T. Guver, A. E. Erkoca, M. H. Reno, and I. Sarcevic (2012), eprint 1201.2400.
 R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, Astron.J. 129, 1993 (2005), eprint astroph/0412641.
 K. Petraki and R. R. Volkas, Int.J.Mod.Phys. A28, 1330028 (2013), eprint 1305.4939.
 K. M. Zurek (2013), eprint 1308.0338.
 A. O. Jamison (2013), eprint 1304.3773.
 J. I. Kapusta and C. Gale, FiniteTemperature Field Theory Principles and Applications (Cambridge University Press, Cambridge, 2006).
 S. Doniach and E. H. Sondheimer, Green’s Functions for Solid State Physicists (W. A. Benjamin, Inc., Massachusetts, 1974).
 S. Reddy, M. Prakash, and J. M. Lattimer, Phys.Rev. D58, 013009 (1998), eprint astroph/9710115.
 K. Saito, T. Maruyama, and K. Soutome, Phys.Rev. C40, 407 (1989).
 S. Reddy, private communication (2013), in reference Reddy et al. (1998), some of the equations had typos. In particular, eqn. (69) should read , eqn. (78) should read , and eqn. (79) should read .
 N. K. Glendenning, Compact Stars (SpringerVerlag, New York, 1997).
 J. Lattimer and M. Prakash, Science 304, 536 (2004), eprint astroph/0405262.
 D. Page and S. Reddy, Ann.Rev.Nucl.Part.Sci. 56, 327 (2006), eprint astroph/0608360.
 D. Page, M. Prakash, J. M. Lattimer, and A. W. Steiner (2011), eprint 1110.5116.
 T. Takatsuka and R. Tamagaki, Prog.Theor.Phys.Suppl. 112, 27 (1993).
 J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175 (1957), URL http://link.aps.org/doi/10.1103/PhysRev.108.1175.
 A. Burrows, S. Reddy, and T. A. Thompson, Nucl.Phys. A777, 356 (2006), eprint astroph/0404432.
 Many pulsars have been observed in globular clusters Camilo and Rasio (2005) which can have DM densities as large as GeV/cm (or perhaps larger) as predicted for the globular cluster M4 McCullough and Fairbairn (2010) however the presence of such large DM densities in globular clusters is uncertain–see also Conroy et al. (2011).
 S. Weinberg, Prog.Theor.Phys.Suppl. 86, 43 (1986).
 A. Schwenk and B. Friman, Phys.Rev.Lett. 92, 082501 (2004), eprint nuclth/0307089.
 P. F. Bedaque, G. Rupak, and M. J. Savage, Phys.Rev. C68, 065802 (2003), eprint nuclth/0305032.
 M. Plischke and B. Bergersen, Equilibrium Statistical Physics (World Scientific Publishing Co. Pte. Ltd., Singapore, 2006).
 D. T. Son and M. Wingate, Annals of Physics 321, 197 (2006), eprint arXiv:condmat/0509786.
 M. Baldo and C. Ducoin, Phys.Rev. C84, 035806 (2011), eprint 1105.1311.
 M. G. Alford, K. Rajagopal, and F. Wilczek, Nucl.Phys. B537, 443 (1999), eprint hepph/9804403.
 M. G. Alford, A. Schmitt, K. Rajagopal, and T. Schäfer, Rev.Mod.Phys. 80, 1455 (2008), eprint 0709.4635.
 P. F. Bedaque and T. Schäfer, Nucl.Phys. A697, 802 (2002), eprint hepph/0105150.
 D. B. Kaplan and S. Reddy, Phys.Rev. D65, 054042 (2002), eprint hepph/0107265.
 D. Son and M. A. Stephanov, Phys.Rev. D61, 074012 (2000), eprint hepph/9910491.
 I. Affleck and M. Dine, Nuclear Physics B 249, 361 (1985), ISSN 05503213, URL http://www.sciencedirect.com/science/article/pii/055032138590%0215.
 P. Ade et al. (Planck Collaboration) (2013), eprint 1303.5076.
 T. Gherghetta, C. F. Kolda, and S. P. Martin, Nucl.Phys. B468, 37 (1996), eprint hepph/9510370.
 N. ArkaniHamed, L. J. Hall, H. Murayama, D. TuckerSmith, and N. Weiner, Phys.Rev. D64, 115011 (2001), eprint hepph/0006312.
 D. Hooper, J. MarchRussell, and S. M. West, Phys. Lett. B605, 228 (2005), eprint hepph/0410114.
 S. Gopalakrishna, A. de Gouvea, and W. Porod, JCAP 0605, 005 (2006), eprint hepph/0602027.
 G. Belanger, M. Kakizaki, E. Park, S. Kraml, and A. Pukhov, JCAP 1011, 017 (2010), eprint 1008.0580.
 T. Hebbeker, Phys.Lett. B470, 259 (1999), eprint hepph/9910326.
 F. Camilo and F. A. Rasio, ASP Conf.Ser. (2005), eprint astroph/0501226.
 M. McCullough and M. Fairbairn, Phys.Rev. D81, 083520 (2010), eprint 1001.2737.
 C. Conroy, A. Loeb, and D. Spergel, Astrophys.J. 741, 72 (2011), eprint 1010.5783.