2 CAPTURE OF DARK MATTER BY NEUTRON STARS

DM reaching terrestrial targets is expected to have low velocities $\beta _{\rm Earth}=\sqrt{\frac{2GM_{\rm Earth}}{R_{\rm Earth}}}$ = v _χ ~ 10⁻³ (taking c = ℏ = 1 units and being G the gravitational constant) due to the fact that the gravitational boost is small for the Earth and, accordingly, the Lorentz factor $\gamma _{\text{Earth}} =1/\sqrt{1-\beta _{\text{Earth}}^2} \sim 1$ . Instead, for more compact objects, i.e. with a larger ratio M/R, the boosting capability increases. For neutron stars (NSs) with masses M _NS ≃ 1.5 M_⊙ and radius R _NS ≃ 12 km, γ_NS ~ 1.26 or v _χ ~ 0.6, providing thus a mechanism to attain higher χ kinetic energies.

It is important to remark that DM could be in principle accreted (and retained) by the star not only during the collapsed stage but also during most of the previous stellar lifetime, being the capture rate, C _χ, provided the interaction strength with ordinary matter is not negligible. Different particle candidates display interaction couplings in a wide range (Bertone & Hooper Reference Bertone and Hooper2016). Regarding the progenitor stages, the progressively denser central core allows the rise of a finite spatial number density, n _χ(r), over time. Later, when a compact object is formed, the efficiency changes and the capture rate enhances, see for example, Gould (Reference Gould1987) and Kouvaris & Tinyakov (Reference Kouvaris and Tinyakov2011a). In this way, assumed an equation of state for regular SM matter in the interior of the NS, at a given galactic location, and with a corresponding ambient DM mass density ρ_χ ~ m _χ n _χ, following Kouvaris (Reference Kouvaris2008) and Guver et al. (Reference Guver, Erkoca, Reno and Sarcevic2014), the capture rate can be written for a weakly interacting DM particle as

where $\bar{v}$ is the average χ velocity in the existing distribution, with a local value ρ_{χ, 0} ≃ 0.3 $\rm \frac{GeV}{cm^3}$ and f _χ is the fraction of particles that undergo one or more scatterings while inside the star. It saturates to f _χ ~ 1 if σ_χN ≳ σ₀, otherwise $f_\chi \sim 0.45 \frac{\sigma _{\chi \text{N}}}{\sigma _0}$ , being $\sigma _0=\frac{m_n R_{\rm NS}^2}{M_{\rm NS}}\sim 10^{-45}$ cm² the geometrical cross-section. As thermalisation times for m _χ ~ 1 GeV particles are much smaller than dynamical cooling times, see Goldman et al. (1989), once inside, DM is believed to diffuse towards the denser central stellar regions according to the exponential law

with n _{0, χ} the central value and $\Phi (r)=\int _0^r \frac{GM(r^{\prime }) \text{d}r^{\prime }}{{r^{\prime }}^2}$ the gravitational potential. Finally, n _χ(r) = n _{0, χ} e ^{−(r/r _th)2} with a thermal radius $r_{\rm th}= \sqrt{\frac{9 k_B T}{8 \pi G \rho _n m_{\chi }}}$ , where ρ _n is the barionic density. Normalisation requires $\int _0^{R_{\text{NS}}} n_{\chi }(r) \text{d}V=N_\chi$ at a given time. N _χ is the number of DM particle population which resides inside. This number will depend, in general, on the capture, annihilation (or decay), and evaporation rates C _χ, C _a , C _s , respectively, see for example, Zentner (Reference Zentner2009). But, for DM particles with masses larger than ~ 2 keV, evaporation effects can be ignored as it is shown in Kouvaris & Tinyakov (Reference Kouvaris and Tinyakov2011b) and Krauss, Srednicki, & Wilczek (Reference Krauss, Srednicki and Wilczek1986). Then, as a function of time t, capture and annihilation processes compete to yield a population

where $\gamma (N_{\chi ,0})=\mathrm{tanh^{-1}} (\sqrt{C_a/C_\chi }N_{\chi ,0})$ and $\tau =1/\sqrt{C_\chi C_a}$ . For a typical progenitor of a compact star, one can estimate $N_{\chi ,0} \sim 10^{39} (\frac{m_\chi }{1\, \rm GeV})$ , see Kouvaris & Tinyakov (Reference Kouvaris and Tinyakov2010), assuming that local densities for average NS galactic distances peak around ~2 kpc where ρ_χ ~ 10²ρ_{χ, 0}. For times t ≳ τ ~ 10⁴ yr, the equilibrium sets $N_\chi (t)\sim \sqrt{C_\chi /C_a}$ .

In the case of asymmetric DM candidates, accretion of DM mass beyond a critical value, i.e. the Chandrasekhar mass, could destroy (collapse) the star over time as it has been studied in Kouvaris & Tinyakov (Reference Kouvaris and Tinyakov2011a), McDermott, Yu, & Zurek (Reference McDermott, Yu and Zurek2012), and Bramante, Fukushima, & Kumar (Reference Bramante, Fukushima and Kumar2013). Possible limiting values of N _χ arising from a fermionic nature provide a value of Chandrasekhar number N _Ch ~ ( M _Pl/m _χ)³ ~ 1.8 × 10⁵⁴ (1 TeV/m_χ)³ with M _Pl the Planck mass.

Another way quoted to limit the dark sector population and thus collapse the whole star stems from the fact that quark deconfinement could be triggered via DM seeding. Energy release from self-annihilating DM could induce spark formation energetic enough to nucleate stable bubbles of deconfined quark matter leading to a softening of the nucleon equation of state. This has been invoked as a mechanism that could drive NS to quark star (QS) conversion, see Pérez-García, Silk, & Stone (Reference Pérez-García, Silk and Stone2010), Pérez-García & Silk (Reference Pérez-García and Silk2012), and Pérez-García, Daigne, & Silk (Reference Pérez-García, Daigne and Silk2013). In addition, arguments along the same line state that unstable DM can also be constrained by structural stability of accreting objects (Pérez-García & Silk Reference Pérez-García and Silk2015).

3 LIGHT DARK MATTER SCATTERING AT HIGH DENSITIES

The DM capture rate Equation (1) mainly determines the magnitude of the novel effects due to the interplay of both matter sectors. For a given DM candidate, the interaction with the relevant degrees of freedom of the system considered will be crucial. This, in turn, comes phenomenologically described through couplings and mathematical operators in the underlaying theory. So far, no satisfactory answers exist coming from a unique model description, although there is a general consensus that the particle or particles in the dark sector should be electrically neutral and cold. The most popular extension to the SM, supersymmetry (SUSY), predicts that each particle has a partner of different spin but similar interactions. The lightest superpartner or LSP is stable in many cases and is often a neutralino, constituting an excellent DM candidate (Bertone & Hooper Reference Bertone and Hooper2016).

In a nuclear medium such as the core of a compact-sized (spherical) object of mass M and radius R, the (in) elastic scattering of χ out of nucleons, σ_χN , mainly determines the relative fraction of DM to ordinary matter. In general, most of the works calculate this quantity without considering in-medium effects, and obtain the mean free path of these DM particles, λ_χ, approximating as λ_χ ≃ 1/σ_χN n, where n is the ordinary nucleon number density, see for example, Spergel & Press (Reference Spergel and Press1985). This is usually considered as being sufficient to obtain knowledge about the internal processed happening inside the most efficient opaque environments. It is also important to know that the isospin charge is relevant as the coupling strengths to protons or neutrons may not be the same or even suffer from cancellations.

In the central core in an NS (with a content ≳ 90% neutrons), is described as having densities higher than several times nuclear saturation density n ₀ ≃ 0.17 fm⁻³. By neglecting medium corrections, one misses important features and qualitative insight into the system. Let us briefly comment on some of the possible issues. To begin with, Fermi-blocking, due to partial restriction of the outgoing nucleon phase space, can play a role diminishing the χN scattering cross-section. Finite temperature effects will allow the population of higher energy states in the nucleon sector with respect to the case in which no temperature effects are considered. It is important to note here that temperatures in the range T ≲ 50 MeV are usually achieved in the very early stages of proto-NS evolution, see for example, Page et al. (Reference Page, Lattimer, Prakash and Steiner2004). Later, after a primary neutrino cooling era, temperatures fall to the ~ keV range. This will effectively set at large times a T ≈ 0 configuration, as thermal energies are indeed much smaller than nucleon Fermi energies k _B T ≪ E _FN in the dense medium (k _B is the Boltzman constant). These in-medium effects are taken into account, for example, in Bertoni, Nelson, & Reddy (Reference Bertoni, Nelson and Reddy2013), Cermeño, Pérez-García, & Silk (Reference Cermeño, Pérez-García and Silk2016a), and Cermeño, Pérez-García, & Silk (Reference Cermeño, Pérez-García and Silk2016b).

As we have mentioned before, usually, the incoming DM particle is supposed to be thermalised in the galaxy with the Maxwellian mean velocities v _χ ~ 10⁻³. However, if the scenario considered is an accreting dense NS, the associated wavelength of the incoming DM particle decreases as $\lambda =\frac{2 \pi }{\sqrt{\gamma ^2 -1}\,m_\chi }$ . This expression sets, in practice, a measure of the validity of the quasielastic approximation, DM particles see nucleons as a whole, since matter is tested to sizes around λ ~ 1 fm, i.e. in the light DM mass range m _χ ≲ 5 GeV. Although further modelling would be required for the description of the inner hadron structure, for this mass range, a monopolar form factor in momentum space, $F(|{\vec{q}}|)=\frac{\Lambda ^2}{\Lambda ^2+q^2}$ with a cut-off parameter Λ = 1.5 GeV can be used, so that F(0) = 1, which can somewhat mitigate the short-range correlations arising in the calculation.

To check the importance of the modifications of the medium in the LDM scattering events, one has to calculate the differential cross-section per unit volume for the interaction between a fermionic particle of Dirac type, χ, and a nucleon field N (protons and neutrons). It is necessary to use the general expression, which can be found in Patrignani & Particle Data Group (Reference Patrignani2016),

where the phase space volume element is

$p^{\prime \mu }_N=(E^{\prime }_N, \vec{p^{\prime }_N})$ and $p^{\mu }_N=(E_N, \vec{p_N})$ are the four momenta for the outgoing and incoming nucleon, and $p^{\prime \mu }_\chi =(E^{\prime }_\chi , \vec{p^{\prime }_\chi })$ and $p^{\mu }_\chi =(E_\chi , \vec{p_\chi })$ those analogous for the DM particle, respectively. Momentum transfer is denoted by q ^μ = p′^μ _χ − p ^μ _χ. In this way, q ₀ = E′ _N − E_N = E _χ − E′_χ and $\vec{q}=\vec{p^{\prime }_N}-\vec{p_N}=\vec{p_\chi }-\vec{p^{\prime }_\chi }$ . The four-dimensional delta assures the energy and momentum conservation in the collision and $|\mathcal {\overline{M}_{\chi \rm N}}|^2$ is the squared matrix element for the interaction considered. The factor $\mathcal {F_{FB}}=f(E_N)(1-f(E^{\prime }_N))$ , with $f(E_i)=\frac{1}{1+e^{({E_i-\mu ^*_i})/{k_{B}T}}}$ i = p,n, and μ* _i the effective nucleon chemical potential for a particle with isospin of ith-type, see Reddy, Prakash, & Lattimer (Reference Reddy, Prakash and Lattimer1998), accounts for the Fermi blocking term. It takes into account the occupation of states and in this scenario affecting only to the nucleons. For the dark sector, it can be assumed that all outgoing DM particles states are in principle allowed and 1 − f _χ(E) ≈ 1 since the fraction of DM inside the star remains tiny at all times. The validity of this approximation is discussed in Cermeño et al. (Reference Cermeño, Pérez-García and Silk2016a). Let us mention here that effective values of nucleon mass and chemical potential define the quasi-particle nature of the nucleon in the medium and differ from the nude values by the presence of average meson fields (Serot & Walecka Reference Serot and Walecka1986).

When discussing the effect of density (and T) in the diffusion of DM in NSs, the quantity we should analyse is the mean free path λ_χ to be obtained from the integrated cross-section per unit volume $\frac{\sigma (E_\chi )}{V}$ , as $\lambda _\chi =(\frac{\sigma (E_\chi )}{V})^{-1}$ .

In a previous step, one obtains the differential value as

where only the lowest order terms are retained. Since v ² _χ ~ v ² _N ≪ 1 given $\frac{|\vec{p_i}|}{E_i}=v_i$ , i = χ, N and nucleon Fermi velocities v _N ~ v _FN = p _FN/E _FN. On the other hand, expressing the energy conservation as a function of the dispersion angle for the outgoing DM particle θ, one obtains the minimun value of $\vec{p_N}$

and

These expressions, Equations (7) and (8), are discussed step by step in Cermeño et al. (Reference Cermeño, Pérez-García and Silk2016a).

We have to remark that this scenario is restricted to temperatures and densities typical for the thermodynamical evolution of the stellar core region, that is T ≲ 50 MeV and n ≃ (1 − 3)n ₀. Let us note that if finite temperature is considered, a detailed balance factors must be added to the medium response to weak probes (Horowitz & Pérez-García Reference Horowitz and Pérez-García2003; Horowitz & Wehrberger Reference Horowitz and Wehrberger1991) under the form

This factor provides the relation between the dynamical nuclear structure factor for positive and negative energy transfers q ₀ as the thermodynamic environment can donate energy to the outgoing particle.

Now, to obtain the total cross-section per unit volume and the inverse of it, i.e. the mean free path, the integrations over all possible outgoing energy transfer values and solid angle have to be performed. In this way,

When particularising the calculation for a spin-independent (scalar-vector) interaction model, the interaction Lagrangian is $\mathcal {L_I}=\sum _{N=n,p}g_{s,N}\chi \overline{\chi } N \overline{N}+g_{v,N}\chi \gamma ^{\mu }\overline{\chi } N\gamma _{\mu } \overline{N}$ , where g _{s, N} (g _{v, N} ) are the effective scalar (vector) couplings of the DM particle to the N field.

The effect of the Pauli blocking can be seen, for example, in Figure 1 of Cermeño et al. (Reference Cermeño, Pérez-García and Silk2016a) where the differential cross-section per unit volume as a function of the energy transfer q ₀ for different values of the transferred momentum $|\vec{q}|$ is shown for a pure neutron system. We use n = n ₀, setting T = 0 and m _χ = 0.5 GeV. The limiting upper value of the energy transfer is E _χ − m _χ ≈ 130 MeV. This is because, in general, −∞ < q ₀ < E _χ − m _χ, since m _χ < E′_χ < ∞. Instead, at T = 0, the energy transfer cannot be negative, so that q ₀ > 0 and 0 < q ₀ < E _χ − m _χ. The triangular shape is due to the Heaviside Fermi distribution at T = 0. Beyond q ₀ values limited by real-valued angles in Equation (7), the scattered states are not allowed since it is kinematically impossible to scatter a nucleon due to lack of empty states. To understand the value of the maxima of this plot, we have to remind ourselves that at T = 0 our allowed states will be given by the factor Θ(E_F − E_N )Θ(1 − f_N (E_F − E′ _N )), and, on the other hand, there is a minimum allowed value for $\vec{p_N}$ given by the kinematics of the interaction, $|\vec{p}_{-}|$ . Then, it is easy to verify that maximum values are reached when $|\vec{p}_{-}|=|\tilde{p_N}|$ being $|\tilde{p_N}|$ the value of $|\vec{p_N}|$ when $|\vec{p^{\prime }}_N|=|\vec{p}_F|$ .

Figure 1. Differential cross-section per unit volume as a function of the energy transfer q ₀ for values of $|\vec{q}|=20, 41, 207$ , and 290 MeV. The DM particle mass is m _χ = 0.5 GeV and T = 0 at n = n ₀. From Cermeño et al. (Reference Cermeño, Pérez-García and Silk2016a).

Finite temperature effects can be observed in Figure 2 where the detailed balance factors have been included and the transferred momentum $|\vec{q}|=20$ MeV, chemical potential μ = E _FN, DM particle mass m _χ = 0.5 GeV and nucleon density n = n ₀ are fixed.

Figure 2. Differential cross-section per unit volume as a function of the energy transfer q ₀ at T = 0, 5, 10 MeV for a nucleon density n = n ₀. We set $|\vec{q}|=20$ MeV and m _χ = 0.5 GeV. From Cermeño et al. (Reference Cermeño, Pérez-García and Silk2016a).

At temperatures T > 0, the negative energy transfer states get increasingly populated and the sharp nucleon distribution is smoothed. As q ₀ → 0, the inverse detailed balance factor S ⁻¹(q ₀, T) → 0. The corresponding divergence will, however, be integrable in order to obtain a finite integrated cross-section. This, can be seen, for example, in Figure 3, where the mean free path for the χ particle as a function of kinetic energy K _χ = E _χ − m _χ for temperatures T = 0 (solid line) T = 10 MeV (dashed line) and T = 30 MeV (dotted line) is plotted. m _χ = 1 GeV and n = n ₀ are considered. The notation used states that E _χ ≡ ω. Note that for T = 0 in the limit of vanishing kinetic energy for the incoming DM particle E _χ → m _χ, the phase space available for the outgoing particles vanishes as the energy transferred q ₀ → 0 due to filled population levels, therefore providing a null cross-section, while this is not true in the finite T case as more channels are available.

Figure 3. DM particle mean free path as a function of kinetic energy for m _χ = 1 GeV at n = n ₀ for T = 0, 10, 30 MeV. From Cermeño et al. (Reference Cermeño, Pérez-García and Silk2016a).

In Figure 4, we show the variation of the DM particle mean free path as a function of density (in units of n ₀) for two values of temperature and fixing m _χ = 1 GeV considering effective nucleon masses for the T = 0 case while not for the finite temperature case. We can see competitive effects from the temperature and effective mass. A steady decrease is obtained in case the naked nucleon mass is considered. Temperature effects, which are relevant in the early stages of dense star evolution, tend to increase the opacity of nucleon matter to prevent DM nearly free streaming. The cross-section is greatly affected by the finite density of matter, namely by the effect of a smaller effective nucleon mass m* _N < m_N . Temperature effects are important although to a lesser extent relative to density ones.

Figure 4. DM particle mean free path as a function of density (in units of n ₀) for two values of temperature, T = 0 and T = 10 MeV. Effective (naked) nucleon mass has been used in the zero (finite) T calculation. From Cermeño et al. (Reference Cermeño, Pérez-García and Silk2016a).

It is important to point out that in this LDM mass range, scattering is diffusive to very good approximation as λ_χ/R _NS ≪ 1. On the other hand, for the choice of couplings strengths consistent with collider constraints as Khachatryan et al. (Reference Khachatryan2015) and Aad et al. (Reference Aad2015), once fixed K energy, an SM neutrino displays typical mean free path somewhat smaller, see Reddy et al. (Reference Reddy, Prakash and Lattimer1998), Horowitz & Wehrberger (Reference Horowitz and Wehrberger1991), and Horowitz & Pérez-García (Reference Horowitz and Pérez-García2003).

4 LIGHT DARK MATTER SCATTERING AT LOW DENSITIES

At low densities, ≲ 10¹³ g cm − 3 inside the star crust matter arranges in a series of irregular shapes or structures called generically pasta phases (Horowitz, Pérez-García, & Piekarewicz Reference Horowitz, Pérez-García and Piekarewicz2004). In particular, the outer crusts in NSs are formed by periodically arranged nuclei with typical densities ranging from ρ ≃ 2 × 10⁶–4 × 10¹¹ g cm³. In the single-nucleus description, see Rüster, Hempel, & Schaffner-Bielich (Reference Rüster, Hempel and Schaffner-Bielich2006), a series of nuclei with increasing baryonic number, A, from Fe to Kr form a lattice before neutrons start to leak out of nuclei. At these high densities, electrons form a degenerate Fermi sea. Around ρ₀ ≃ n ₀ m _N ≃ 2.4 × 10¹⁴ g cm³, and beyond that, matter can be considered a homogeneous system.

When the phase space is restricted to fermionic DM particles with masses m _χ ≲ 100 MeV, the relevant cross-section of the interaction between DM particles and the outer stellar layers is that scattering off nuclei. Taking this into account, the possibility of production of quantised lattice vibrations (phonons) in the periodically arranged structures in the NS outer crust by LDM scattering, arises. These extra phonons can impact subsequent quantities of interest, such as the ion thermal conductivity, that are relevant for computing the cooling behaviour of NSs.

Phonons are quantised vibrational modes characterised by a momentum $\vec{k}$ and polarisation vector $\vec{\epsilon }_{\lambda }$ appearing in a nuclear periodic system, as it can be seen in Ziman (Reference Ziman1960). They can have a number of different sources. They can be excited due to non-zero temperature T in the medium. The Debye temperature allows us to evaluate the importance of the ion motion quantisation. For a bcc lattice, see Carr (Reference Carr1961), T _D ≃ 0.45T _p, being $T_{\text{p}}=\omega _{\text{p}} /k_{\rm B}=\sqrt{\frac{4 \pi n_AZ^{2} e^{2}}{k_{\rm B}^{2} m_{A}}}$ the plasma temperature associated to a medium of ions with number density n _A , baryonic number A, electric charge Ze and mass m _A . At low temperatures T < T _D, the quantisation becomes increasingly important and the thermal phonons produced are typically acoustic modes, following a linear dispersion relation $\omega _{{k},\lambda }=c_{l, \lambda }|{\vec{k}}|$ , where $c_l=\frac{\omega _{\text{p}}/3}{ (6\pi ^{2} n_A)^{1/3}}$ is the sound speed. In addition, phonon production can be caused by an external scattering agent, for example, SM neutrinos. In this respect, weak probes such as cosmological neutrinos with densities n _ν ~ 116 cm⁻³ per flavour have been shown to provide small phonon production rates in a crystal target, see Ferreras & Wasserman (Reference Ferreras and Wasserman1995). Due to the tiny mass of the neutrino, the experimental signature of this effect seems however hard to confirm.

The single phonon excitation rate (per unit volume) is obtained as Cermeño et al. (Reference Cermeño, Pérez-García and Silk2016b)

where $f_{\chi }(\vec{p_\chi })=\frac{n_\chi \mu }{4 \pi m_\chi ^3 K_2(\mu )} e^{-\mu \sqrt{1+\frac{|\vec{p_\chi }|^2}{m_\chi ^2}}}$ is the Maxwell–Juttner distribution (Juttner Reference Juttner1911) function for relativistic incoming DM particles. With $\mu =\frac{m_\chi }{k_{B\rm }T}\approx 6.7$ for $\sqrt{<v^2>}\sim 0.6$ , see Hakim (Reference Hakim2011) and Cercignani & Kremer (Reference Cercignani and Kremer2002), and K ₂(μ) the modified Bessel function of second kind. The value of a takes into account the length of the interaction when a DM particle approaches a nucleus in the periodic lattice. Using the Born approximation for the cross-section (off protons and neutrons) in the centre of mass frame, it can be written as

with $\int _{-1}^{1} 2\pi \text{d}( \cos \; \theta _{\chi }) |\mathcal {\overline{M}_{\chi \rm N}}|^2 \equiv |\tilde{\mathcal {{M}_{\rm N}}}|$ and the averaged matrix element considering a prescribed model interaction. In this expression, Pauli blocking effects are not considered but, as an approximation, density dependence can be retained using a parameterisation of the nuclear Fermi momentum $|\vec{p}_{\rm FN}| \sim (3 \pi ^2 n_0 Y_N)^{1/3}$ and the nuclear fractions Y _p = Z/A, Y _n = (A − Z)/A. It is important to say that the allowed values for the phonon momentum are restricted for the Born approximation, $|\vec{k}|\ll \frac{1}{a}$ , as well as for kinematical restrictions when energy conservation is imposed, $0\le |\vec{k}|\le 2 \left( \frac{c_lE_\chi }{(c_l^2-1)}+\frac{|\vec{p_\chi }|}{|c_l^2-1|}\right)$ .

The single phonon excitation rate (per unit volume) R ⁽⁰⁾ _k due to the accretion of DM particles is almost constant with $|\vec{k}|$ , whereas for neutrinos $R^0_\nu (|\vec{k}|)=R_{\nu 0} e^{ -b|\vec{k}|}$ with b a constant which depends on neutrino mass. This means that, for a fixed value of $|\vec{k}|$ , contribution from DM particles to the phonon excitation rate will be much higher than neutrino contribution.

Now, it will be important to discuss the astrophysical impact of these phonons. As we have mentioned before, phonon production can be crucial for determining further transport properties, in particular, thermal conductivity in an ion-electron system such as that in the outer NS crust. As an important contribution to the total ion conductivity, κ _i , partial ion conductivities due to ion–ion, κ_ii ≡ κ_ph, and ion-electron collisions, κ_ie, must be added, see for example, Chugunov & Haensel (Reference Chugunov and Haensel2007), under the prescription $\kappa ^{-1}_{\text{i}}=\kappa ^{-1}_{\text{ii}}+\kappa ^{-1}_{\text{ie}}$ . Standard mechanisms to produce lattice vibrations include thermal excitations, as analysed in detail in previous works (Negele & Vautherin Reference Negele and Vautherin1973; Potekhin, Chabrier, & Yakovlev Reference Potekhin, Chabrier and Yakovlev1997; Baiko et al. Reference Baiko, Kaminker, Potekhin and Yakovlev1998). In an NS, the outer crust can be modelled under the one-component-plasma description. This low-density solid phase can be classified according to the Coulomb coupling parameter Γ = Z ² e ²/ak _B T, where a = (4πn _A /3)^1/3 is the ion sphere radius. It is already known that typically for Γ ⩾ Γ_m ≃ 175, or below melting temperature T < T _m, single-ion systems crystallise, see Nagara, Nagata, & Nakamura (Reference Nagara, Nagata and Nakamura1987).

There are a number of processes that can affect thermal conductivity in the medium. The so-called U-processes are responsible for modifying the electron conductivity such that for high temperatures, T > T _U, electrons move almost freely, more details in Ziman (Reference Ziman1960). Assuming a bcc lattice, T _U ≃ 0.07T _D. Thus, in the scenario depicted here, the temperature range must be T _U < T < T _D < T _m for each density considered. According to kinetic theory, the thermal conductivity can be written in the form

as it can be seen in Ziman (Reference Ziman1960), where $C_A=9 \left(\frac{T}{T_{\text{D}}}\right)^3 \int _0^{T_{\text{D}/T}} \frac{x^4 e^x \text{d}x}{(e^x -1)^2}$ is the phonon (dimensionless) heat capacity (per ion) and L _ph is an effective phonon mean free path that includes all scattering processes considered: U-processes and impurity (I) scattering processes (both dissipative) and the phonon normal (N) scattering which are non dissipative $L^{-1}_{\text{ph}}=L^{-1}_{\text{U}}+L^{-1}_{\text{I}}+L^{-1}_{\text{N}}$ . Typically, the thermal conductivity is related to the thermal phonon number at temperature T, L _ph ~ 1/N _{0, kλ}, where N _{0, kλ} = (e ^{ω _kλ/k _B T} − 1)⁻¹. The contribution from DM can be obtained by the net number of phonons that results from the competition between thermal phonons and scattering excitation and stimulated emission, as it is shown in Ferreras & Wasserman (Reference Ferreras and Wasserman1995), in a 4-volume δVδt. Using the averaged rate per unit volume and weighting with the incoming distribution providing the frequencies of different values of momenta, the total number of phonons can be written as

see Cermeño et al. (Reference Cermeño, Pérez-García and Silk2016b). Where ${\tilde{R}^{(0)}_k}$ is the single phonon excitation rate for each particular momentum value (not averaged over incoming χ momenta). Since the source (NS) is in relative motion to the LDM flux, there is a Doppler shift characterised by the source velocity v ≡ v _NS ~ 10⁻² i.e. galactic NS drift velocity. As we have mentioned before, the distribution of NSs in our galaxy peak at distances ≲ few kpc ((Lorimer et al. Reference Lorimer2006)) where the DM density is enhanced with respect to the solar neighbourhood value.

In Figure 5, the phonon thermal conductivity as a function of barionic density is shown for two different typical temperatures for the base of the crust, using m _χ = 100 MeV. Solid lines are the standard thermal result with no DM. Dash-dotted and dashed lines correspond to n _χ/n _{0, χ} = 10, 100, respectively. It is clear to see that at the largest LDM local densities considered, there is an enhancement over the thermal result well inside the outer crust. This corresponds to the site where the DM-induced effects have the most influence ((Pérez-Azorín, Miralles, & Pons Reference Pérez-Azorín, Miralles and Pons2006)) as this is the most massive part of the outer crust. Although, below these densities, there is a negligible change. At lower T, the effect of a perturbation over the thermal phonon population is more important. Enhanced (decreased) conductivities at moderate LDM densities are due to a net reduction (increase) of the number of phonons in the lattice as a result of cancellation of modes.

Figure 5. Phonon thermal conductivity as a function of density (in units of 10¹⁰ g cm³) for temperatures T = 510⁷ K (blue), 510⁸ K (red), and m _χ = 100 MeV. Dash-dotted and dashed lines depict the impact of a LDM density n _χ/n _0χ = 10, 100. Solid lines are the standard thermal result with no DM for each case. From Cermeño et al. (Reference Cermeño, Pérez-García and Silk2016b).

The same type of plot can be shown in Figure 6 in a more realistic stellar context. Here, the phonon thermal conductivity as a function of baryonic density at T = 10⁸ K and m _χ = 65 MeV can be seen. Solid, dot-dashed, and dashed lines correspond to cases with no DM, n _χ/n _{0, χ} = 10, 100, respectively. Electron thermal conductivity is also shown for magnetised realistic scenarios in the perpendicular direction to a magnetic field B of strength B = 10¹⁴ G (dotted) and B = 10¹⁵ G (doble dotted), these NSs with so high magnetic fields are called magnetars and they are being studied by many groups, see for example, a review from Turolla, Zane, & Watts (Reference Turolla, Zane and Watts2015). Ions are mostly not affected by the presence of a magnetic field for the range of temperatures and magnetic fields considered, this can be checked in Chugunov & Haensel (Reference Chugunov and Haensel2007). The parallel direction electronic contribution is not depicted here since it is typically much larger k _e∥ ~ 10¹⁷ − 10¹⁹ erg cm⁻¹ s⁻¹ K⁻¹. Values of electronic contribution have been obtained from Viganò (Reference Viganò2013). Since the global conductivity is κ = κ_e + κ_ph, this result is expected to contribute to the reduction of the difference in heat conduction in both directions and thus to the isotropisation of the NS surface temperature pattern, which is studied for example in Perna et al. (Reference Perna, Viganò, Pons and Rea2013) and Viganò et al. (Reference Viganò, Perna, Rea and Pons2014), as it would be smoothly driven towards more isothermal profiles for latitudes among pole and equator. It is already known, see for example, Kaminker et al. (Reference Kaminker, Yakovlev, Potekhin, Shibazaki, Shternin and Gnedin2007), that the outer crust plays an important role in regulating the relation among temperature in the base of it and the surface. The detailed calculation remains, however, for future works.

Figure 6. Phonon thermal conductivity as a function of density (in units of 10¹⁰ g cm³) at T = 10⁸ K and m _χ = 65 MeV. Solid, dot-dashed, and dashed lines correspond to cases with no DM, n _χ/n _{0, χ} = 10, 100. Perpendicular electron thermal conductivity is also shown for B = 10¹⁴, 10¹⁵ G. From Cermeño et al. (Reference Cermeño, Pérez-García and Silk2016b).

5 IN-MEDIUM ENERGY RELEASE FROM ANNIHILATING DARK MATTER

In previous sections, we consider a scattering interaction between DM and ordinary matter. Apart from the spin-independent interaction model considered before (corresponding to effective scalar and vector couplings), more effective operators for fermionic DM particles can be studied. For example, the so-called Coy DM model proposed by Boehm et al. (Reference Boehm, Dolan, McCabe, Spannowsky and Wallace2014) which belongs to the family of simplified models, see for instance, Dolan et al. (Reference Dolan, Kahlhoefer, McCabe and Schmidt-Hoberg2015), Bauer, Haisch, & Kahlhoefer (Reference Bauer, Haisch and Kahlhoefer2017), and Baek, Ko, & Li (Reference Baek, Ko and Li2017), includes two new particles to the SM, i.e. a Dirac fermion DM candidate, χ which interacts with a boson a, with mass m _a , through the pseudoscalar coupling of strength g _χ as proposed in Arina, Del Nobile, & Panci (Reference Arina, Del Nobile and Panci2015). The interaction lagrangian of the model reads

A final state of fermion–antifermion is formed through the boson mediator. The vertex involves a pseudoscalar coupling with the SM fermions f with strength g ₀ g _f , where g _f denotes a multiplicative factor common to all couplings of a with SM fermions. The most common type of couplings are the flavour-universal coupling g _f = 1, independent of the fermion type, and Higgs-like coupling, proportional to the fermion mass $g_f=\frac{m_f}{174 \; \rm GeV}$ . However, there are other schemes where a couples either to quarks or leptons exclusively, and with a flavour structure (Dolan et al. Reference Dolan, Kahlhoefer, McCabe and Schmidt-Hoberg2015).

An existing DM distribution inside the star capable of self-annihilating can be a source of energy emission for the dynamical cooling behaviour, once the star is formed. In the range m _χ < m _Higgs and m _a < m _χ, the relevant annihilation processes are $\chi {\chi } \rightarrow f\bar{f}$ and χχ → aa, see Arina et al. (Reference Arina, Del Nobile and Panci2015).

In this context, these annihilation reactions, $\chi {\chi } \rightarrow f\bar{f}$ and χχ → aa, with subsequent decay $a \rightarrow f\bar{f}$ can be considered to explore possible astrophysical consequences in the neutrino fermionic channel as dense stars are efficient DM accretors. One of the key quantities that can govern the internal stellar energetic balance is the local energy emissivity, $Q_{\text{E}}=\frac{\text{d}E}{\text{d}V\text{d}t}$ , which can be understood as the energy emitted per unit volume per unit time, through a prescribed particle physics reaction.

Considering the above quoted reactions which involve annihilating DM accreted by a NS, with baryonic density ρ _b ~ 2ρ₀, from an existing galactic distribution. Q _E denotes the energy emission rate per volume due to fermionic or pseudoscalar pair emission in the final states. In general, its expression can be written under the form (Esposito et al. (Reference Esposito, Mangano, Miele, Picardi and Pisanti2002))

with $\text{d}\Phi =\frac{\text{d}^3p_1}{2(2 \pi )^{3}E_1} \frac{\text{d}^3p_2}{2(2 \pi )^{3}E_2} \frac{\text{d}^3p_3}{2(2 \pi )^{3}E_3} \frac{\text{d}^3p_4}{2(2 \pi )^{3}E_4}\, (2 \pi )^4 \delta ^4 (p_1+p_2- p_3-p_4)$ the 4-body (12 → 34) phase space element and $|\mathcal {\overline{M}}|^2$ the spin-averaged squared matrix element of the reaction considered. The additional factor f(f ₁, f ₂, f ₃, f ₄) accounts for the initial and final particle distribution functions contribution which will be different for each annihilation channel. The notation used is the following $p_1=(E_1, \vec{p_1})$ , $p_2=(E_2, \vec{p_2})$ as the incoming 4-momenta while $p_3=(E_3, \vec{p_3})$ , $p_4=(E_4, \vec{p_4})$ are the outgoing 4-momenta, respectively. Feynman diagrams and matrix elements can be seen in Cermeño, Pérez-García, & Lineros (Reference Cermeño, Pérez-García and Lineros2017).

For the case of annihilation into fermionic pairs $f(f_1,f_2,f_3,f_4) =f_{\chi }(E_1)f_{\bar{\chi }}(E_2)(1-f_{f}(E_3))(1-f_{\bar{f}} (E_4))$ and f _χ, f _f are the local stellar distribution functions for DM and fermionic particles, respectively, containing density and temperature dependence. Whereas for the annihilation into pseudoscalars, the factor $f(f_1,f_2,f_3, f_4)=f_\chi (E_1)f_{\bar{\chi }}(E_2)f_a (E_3)\\f_a(E_4)$ . Given the fact that inside the star, the small DM fraction provides a Fermi energy E _{F, χ} ≪ k _B T, in both reactions, its distribution function can be approximated by a classical Maxwell–Boltzmann type $f_\chi =f^{\text{MB}}_\chi (|\vec{p_i}|)=\left( \frac{1}{2\pi m_\chi k_B T}\right)^{\frac{3}{2}}n_\chi (r) e^{\frac{-|\vec{p_i}|^2}{2m_\chi T}}$ . We can see that the phase space factor f(f ₁, f ₂, f ₃, f ₄) in Equation (16) introduces DM density and T dependence into the calculation as a thermalised DM distribution exists inside the NS core. For the outgoing fermions, the medium density effects have to be taken into account too through the phase space blocking factors as we have reviewed in Section 3, whereas for the outgoing pseudoscalars (bosons) no blocking has to be considered.

In order to analyse an interesting astrophysical consequence, Cermeño et al. (Reference Cermeño, Pérez-García and Lineros2017) restrict the analysis to three different sets of flavour-universal (g _f = 1) parameters that are not in conflict with existing phenomenology of direct detection experiments (Bertone & Hooper Reference Bertone and Hooper2016) nor cosmological bounds (Zeldovich Reference Zeldovich1965; Zeldovich, Okun, & Pikelner Reference Zeldovich, Okun and Pikelner1965). The masses and couplings considered can be seen in Table 1. Models A and B are mainly determined by DM relic abundance due to the fact that the DM mass is in the region where direct detection experiments are less restrictive but can be more constrained with Kaon decays. On the contrary, the couplings in set C are constrained by LUX results, Akerib et al. (Reference Akerib2017), in spin-independent and spin-dependent cross-sections. For more details, see Cermeño et al. (Reference Cermeño, Pérez-García and Lineros2017). On the other hand, they fixed the final fermion states to neutrinos. These weakly interacting SM fermions play a key role in releasing energy from NSs. At the end of the life of a very massive star, the central core collapses gravitationally and it is well known that its gravitational binding energy is emitted in neutrinos (and antineutrinos) of the three families. A very efficient cooling scenario emerges in the first ~10⁵ yr. Standard processes such as the URCA cooling or the modified URCA (MURCA) cooling, see Friman & Maxwell (Reference Friman and Maxwell1979) and Yakovlev & Levenfish (Reference Yakovlev and Levenfish1995), can release neutrinos with associated emissivities $Q^{\rm URCA}_{\text{E}}\sim 10^{27} R (\frac{k_{\rm B}T}{0.1 \, \rm MeV})^6$ erg cm⁻³ s⁻¹ and $Q^{\rm MURCA}_{\text{E}}\sim 10^{21} R (\frac{k_{\rm B}T}{0.1 \, \rm MeV})^8$ erg cm⁻³ s⁻¹, respectively. R is a control function of order unity. These standard processes release energy from the baryonic system, having a very different effect in the energetic balance when compared to average energy injection from stellar DM annihilation processes (Kouvaris Reference Kouvaris2008). It is important to remark that in this scenario, it can be assumed that neutrinos do not get trapped after being produced and therefore their Fermi-Dirac distributions fulfill f _ν ~ 0.

In this context, to compare the result with the standard physics cooling, in the left panel of Figure 7, the energy emissivity for the process χχ → νν as a function of temperature for three sets of DM parameters in Table 1 is shown. The MURCA energy emissivity (solid line) is plotted too for the sake of comparison. Although this latter is not the only process contributing to the effective cooling, it sets an upper limit to standard emissivities in the scenario depicted due to the fact that to sustain the URCA fast reaction central densities must be higher to provide the ~11% protons fraction required (Lattimer et al. Reference Lattimer, Pethick, Prakash and Haensel1991). It is easy to see that around T ~ 0.1 MeV local central emissivities $\rm log_{10}\, ({\it Q}_E)\sim 21\text{--}22$ are balanced by the heating processes. It is important to remind ourselves that T ~ 0.01 MeV for a NS of ~10⁵ yr. But the most important result is obtained for the reaction χχ → aa with subsequent decay a → νν, Figure 7, right panel. In this case, the neutrino emissivity is largely enhanced with respect to the direct production of neutrinos χχ → νν. For model C (somewhat heavier DM particles), Q _E matches and surpasses the standard MURCA emission at T ≲ 1 MeV while model B has analogous behaviour for T ≲ 0.3 MeV. Note that both models provide larger emissivities in this channel than the direct neutrino pair emission in the analysed T range in spite of being in a different range of masses. Model A provides similar results for both emission reactions, however. This effect of enhanced emissivities can have an impact on internal temperatures and matter phases depending critically on them, as can be seen in Page et al. (Reference Page, Lattimer, Prakash, Steiner, Bennemann and Ketterson2014). In this sense, works as Page et al. (Reference Page, Prakash, Lattimer and Steiner2011) and Shternin et al. (Reference Shternin, Yakovlev, Heinke, Ho and Patnaude2011) quote that the rapid cooling of the Cas A may be an indication of the existence of global neutron and proton superfluidity in the core.

Figure 7. Energy emissivity from channels χχ → νν (left) and via pseudoscalar mediators χχ → aa (right) with subsequent decay a → νν as a function of temperature. Standard emission refers to MURCA processes. From Cermeño et al. (Reference Cermeño, Pérez-García and Lineros2017).

As most of the fraction of DM particles will be in the central region of the NS, characterised by the length r _th, the impact of the radial extent of an emitting inner region through these new processes can be correlated to the ratio $\xi =(\sqrt{2} r_{\rm th}/R_b)$ , which make reference to the fraction where 95% of the approximately Gaussian distribution of DM particles can be found, versus R _b ~ 9 km, radial extent of the boundary, or limit of the core-to-crust region. Since the crust region has a tiny mass, in this scenario, it is not necessary to consider the refinement which it is taken into account in Cermeño et al. (Reference Cermeño, Pérez-García and Silk2016b). For the models analysed in Cermeño et al. (Reference Cermeño, Pérez-García and Lineros2017), this fraction ranges from ξ ∈ [0.03, 0.42] at T ~ 1 MeV to ξ ∈ [0.01, 0.14] at T ~ 0.1 MeV.

Focusing in a fixed stellar evolutionary time, t, to analyse the temperature evolution, the radial (internal) heat equation r < R _b can be written as

where κ is the thermal conductivity given by a contribution from baryons and another from electrons Shternin & Yakovlev (Reference Shternin and Yakovlev2007) κ = κ _b + κ _e and Q _ν is the standard MURCA emissivity, $Q_\nu \sim Q_{\text{E}}^{\rm MURCA}$ . The DM heating term, H(r), has the radial dependence induced by the DM density so that H(r) = H ₀ e ^{−(r/r ₀)2} with $H_0=Q_{\text{E}}^{aa}n^2_{0,\chi }$ and $r_0=r_{\rm th}(T)/\sqrt{2}$ . Although a detailed solution of the full evolution equation remains for future works that should consider the dynamical part, in Cermeño et al. (Reference Cermeño, Pérez-García and Lineros2017) an study of the radial T-profile at different fixed evolutionary times is done solving Equation (17) and assuming a flat initial profile to see how it is distorted from the presence of internal cooling and heating terms. They find that, for example, fixing a value $T(r)\equiv \bar{T}$ for r ∈ [0, R _b ] the central value T ₀ is related to it through the approximate solution $T(r)=T_0+{\alpha } \frac{r^2}{6}-\beta \frac{r_0^2}{2} [1-\sqrt{\pi }\frac{r_0}{2} \frac{\rm erf(r/r_0)}{r}]$ , where erf (x) is the error function and ${\alpha }=Q_\nu (\bar{T})/\kappa ({\bar{T}})$ and $\beta =H_0({\bar{T}})/\kappa ({\bar{T}})$ are coefficients that can be approximated at $\bar{T}$ . Typical values can be seen in Cermeño et al. (Reference Cermeño, Pérez-García and Lineros2017). The important thing here is that if we fixed, for example, $\bar{T}=10^9 K\sim 0.1\ \rm MeV$ , one can verify that the net effect of models B and C is heating the central volume while model A provides net cooling. However, it is important to remark that the star will self-adjust in a consisting way that has to be dynamically determined to see at what extent the subsequent thermal evolution is affected in the stellar volume. This could be studied through a full dynamical simulation to see how the cooling mechanisms adjust inside the star as a temporal sequence.