# Enhanced neutrino emissivities in pseudoscalar-mediated
Dark Matter annihilation in neutron stars

M. Cermeño1<EMAIL_ADDRESS>M. A. Pérez-García1<EMAIL_ADDRESS>and R. A. Lineros2<EMAIL_ADDRESS>1 Department of Fundamental Physics,
University of Salamanca, Plaza de la Merced s/n, E-37008 Salamanca, Spain
2 Instituto de Física Corpuscular – CSIC/U. Valencia,
Parc Científic, calle Catedrático José Beltrán 2, E-46980 Paterna, Spain.

(Received January 1, 2019; Revised January 7, 2019)

###### Abstract

We calculate neutrino emissivities from self-annihilating dark matter ($\chi$)
in the dense and hot stellar interior of a (proto)neutron star. Using a model
where dark matter interacts with nucleons in the stellar core through a
pseudoscalar boson ($a$) we find that the neutrino production rates from the
dominant reaction channels $\chi\chi\rightarrow\nu\bar{\nu}$ or
$\chi\chi\rightarrow aa$, with subsequent decay of the mediator
$a\rightarrow\nu\bar{\nu}$, could locally match and even surpass those of the
standard neutrinos from the modified nuclear URCA processes at early ages. We
find that the emitting region can be localized in a tiny fraction of the star
(less than a few percent of the core volume) and the process can last its
entire lifetime for some cases under study. We discuss the possible
consequences of our results for stellar cooling in light of existing dark
matter constraints.

dark matter, neutron star, neutrinos, emissivity

##  1 Introduction

Dark Matter (DM) is an essential ingredient of the standard cosmological
model. We now know it constitutes nearly 85% of the Universe matter density.
However, despite the tremendous amount of progress that has been made in the
search for this missing type of matter, both on the theoretical and
experimental fronts, its true nature remains an open question. The Standard
Model (SM) of particle physics alone cannot explain the nature of DM,
suggesting that it must be extended. Many theoretical model proposals have
arisen aiming to explain the existing phenomenology (Bertone & Hooper, 2016).
The possible interplay between ordinary and DM could reveal interesting novel
features, thus serving as a smoking gun evidence for the existence of a dark
sector. As examples of the previous, one could cite a possible contribution to
the reionization of the Universe and the increase of gas temperature prior to
the reionization epoch, leaving a potentially detectable imprint on the
cosmological 21-cm signal, as studied in Chuzhoy (2007) and Mapelli et al.
(2006). More in particular, the production of SM neutrinos from annihilation
of proposed dark candidates, generically $\chi$, with energy $E_{\nu}\lesssim
m_{\chi}$ is of paramount importance for the description of internal dynamics
and energetic balance in stellar scenarios. In order to be specific, for
example, one can consider the solar context. A DM particle will be
gravitationally captured by the Sun if, in scattering against solar nuclei, it
falls below the local escape velocity. This accumulation mechanism can lead to
a local stellar DM density higher than that of the galactic halo where it
resides, potentially providing us an opportune region in which to search for
visible signatures (Kouvaris, 2007; Vincent et al., 2015; Rott et al., 2015).
Additionally, low density environments in solar-type stars could yield
interesting features in the neutrino channel (Aartsen et al., 2012; Palomares
& Pascoli, 2008).

On the other hand, in denser stellar environments such as those leading to the
formation of neutron stars (NSs), neutrinos are vastly produced as they are
very efficient at releasing the excess of gravitational energy when a compact
stellar object is formed from a more massive progenitor. It has been now more
than 30 years since the supernova SN1987A event allowed us to glimpse the
complex behaviour of the neutrino internal dynamics and obtain confirmation of
the existence of a preliminary neutrino trapping phase followed by a
transparency era (Gusakov et al., 2004) from the neutrino telescopes on Earth
(Yuksel & Beacom, 2007). In addition, X-ray satellite measurements have also
provided indications of the cooling sequence, for a catalog of isolated
cooling NSs, see (Vigano et al., 2013) and (Yakolev & Levenfish, 1995).
Although a global understanding of the extracted temperatures for these
objects is still missing, the so-called minimal cooling mechanism has been
successful at reproducing the trends of observed cooling curves (Page et al.,
2004). When solving for the internal temperature profile $T(r,t)$ as a
function of stellar radius $r$ and time $t$, one of the key ingredients that
can dictate the energetic balance is the local energy emissivity,
$Q_{E}=\frac{dE}{dVdt}$, i.e. the energy produced per unit volume per unit
time, through a prescribed particle physics reaction.

In this work we will be interested in obtaining astrophysical neutrino
emissivities related to novel reaction channels involving DM undergoing self-
annihilation processes inside the star. In particular we focus on models in
which dark matter particles communicate with the visible sector through a
pseudoscalar mediator. These have been quoted to be well-motivated both from
theoretical and from phenomenological grounds. Some of these models belong to
a set of the so-called simplified type including Boehm et al. (2014), Wild
(2016), Bauer et al. (2017) and Baek et al. (2017). As mentioned, they extend
the SM by (at least) two particles, a DM candidate as well as a state that
mediates the DM interactions with the visible sector, and are able to capture,
with a minimal set of assumptions, some important features of more
ultraviolet-complete (UV) theories while providing a (semi-)consistent
framework in order to analyse the experimental results (Banerjee et al.,
2017).

In this setting we will be interested in the dominant neutrino production
processes, i.e. the s-wave process $\chi\chi\rightarrow\nu\bar{\nu}$ and the
p-wave process $\chi\chi\rightarrow aa$, with subsequent decay
$a\rightarrow\nu\bar{\nu}$. Although the previous reactions constitute the
main neutrino emission channels in this model setting, additional reactions
like e.g. radiative $a$-emission or $\chi\chi\rightarrow aaa$ could also
happen but we will not consider them here as they are subdominant. As we will
show, the two main reactions could provide a contribution to the standard
astrophysical neutrino emissivities in NS environments of sizable magnitude at
early times. The observability of such -indirect- effects caused by DM seems
nowadays difficult as it could be critically relying on the finest
capabilities of current and future X-ray and gamma satellites (NICER, eXTP,
LOFT, ATHENA, CHANDRA).

##  2 Neutrino emissivities from DM annihilation

In this work we are interested in calculating neutrino emissivities from DM
self-annihilation in dense and hot stellar interiors, i.e. that of a (proto)
NS. In order to carry out our calculation we choose a model where DM particles
interact with SM particles through a pseudoscalar mediator. This kind of
models are well-motivated, both from the theoretical and phenomenological
point of view. With direct detection bounds being typically subleading in such
scenarios, the main constraints arise from collider searches (meson bounds)
and from indirect detection experiments (Banerjee et al., 2017). Examples
recently used along this line include the coy dark matter model of Boehm et
al. (2014) and others (Wild, 2016; Bauer et al., 2017; Baek et al., 2017).
Although popular, we must stress that these simplified models have some
limitations, regarding construction itself, and when confronted with bounds
from collider searches (Goncalves et al., 2017; Dolan et al., 2015).

We now introduce our concrete model realization. We consider a model where the
SM field content is extended by a Dirac fermion, $\chi$, with mass $m_{\chi}$,
which plays the role of a dark matter candidate, and a pseudoscalar field,
$a$, with mass $m_{a}$, which mediates the interaction of ordinary and dark
sectors. The interaction lagrangian of the model reads

$\mathcal{L_{I}}=-i\frac{g_{\chi}}{\sqrt{2}}a\bar{\chi}\gamma_{5}\chi-
ig_{0}\frac{g_{f}}{\sqrt{2}}a\bar{f}\gamma_{5}f\,,$ (1)

where $g_{\chi}$ is the DM-mediator coupling, $g_{f}$ corresponds to the
couplings to the SM fermions, $f$, and $g_{0}$ is an overall scaling factor.
From the usual schemes used for matter couplings when introducing Beyond-
Standard-Model motivated physics we will restrict for simplicity to the so-
called flavour-universal, which sets $g_{f}=1$ for all SM fermions. Let us
recall, however, that there are other schemes where $a$ couples either to
quarks or leptons exclusively, and with a flavour structure which will be
treated elsewhere.

Typically, in these models DM phenomenology is controlled by four parameters,
$m_{\chi}$, $m_{a}$, $g_{\chi}$, and $g_{0}g_{f}$. In the range
$m_{\chi}<m_{\rm Higgs}$ and $m_{a}<m_{\chi}$, the relevant annihilation
processes into two-body final states (Abdullah et al., 2014; Arina et al.,
2015) are s-wave $\chi{\chi}\rightarrow f\bar{f}$ and p-wave
$\chi{\chi}\rightarrow aa$. As a remark it is worth mentioning that, as
presented, the most straightforward UV-completion of this setup would be in
the framework of the two Higgs doublet model or models involving even more
extended scalar sectors. However, one should keep in mind that additional
interactions with extra scalars arise at tree-level and that can introduce
important phenomenological model-dependent features (Haber & O’neil, 2011).

Despite the limitations of simplified models, in our particular realization it
is reasonable to expect that the very light mediators will not distort the
relic density predictions due to the presence of additional annihilation
channels involving these extra scalars as discussed in Banerjee et al. (2017).

DM abundance in our universe is likely to be fixed by the thermal freeze-out
phenomenon: DM particles, initially present in our universe in thermal
equilibrium abundance, annihilate with one another until chemical equilibrium
is lost due to the expansion of the universe. The present-day relic density of
these particles is predictable and it has been measured by Planck (Ade et al.,
2015) to be $\Omega_{\rm CDM}h^{2}=0.1198\pm 0.0015$.

Due to the pseudoscalar portal considered here, this model provides spin
dependent interactions with nucleons (N) at tree level. In this way the
$\chi$-N interaction considered in direct searches is suppressed because it is
momentum dependent, see Freytsis & Ligeti (2011), Cheng & Chiang (2012) and
Gresham & Zurek (2014) for details. Instead, the spin independent cross
section is not present at tree level but the effective interaction at one-loop
can be constructed (Ipek et al., 2014). Estimations of both cross sections in
vacuum are given in Freytsis & Ligeti (2011). Both features regarding the
behaviour of the cross section impact the capability of the star to capture DM
during the stage of progenitor and in the collapsed configuration, although
they can compensate each other in the star lifetime in order to have a finite
meaningful amount of DM populating the object (Kouvaris & Tinyakov, 2011).

Usual model analysis considers sets of parameters with a variety of bounds at
different level of significance. Here, in order to be definite, we will
restrict our analysis to three different sets of flavour-universal parameters
that are not in conflict with existing phenomenology to describe light DM
($m_{\chi}\lesssim 30$ GeV) interactions with ordinary matter. We consider
constraints from direct detection experiments (Bertone et al., 2005),
cosmological bounds (Zel’dovich et al., 1965) and collider bounds (Dolan et
al., 2015). The masses and couplings used in this work appear in Table 1.

Table 1: Parameters used in this work as appearing in the interaction lagrangian in Eq. (1).  Model aa We use flavour-universal $g_{f}=1$. |  $m_{\chi}$ [GeV] | $m_{a}$ [GeV] | $g_{\chi}$ | $g_{0}$
---|---|---|---|---
A | 0.1 | 0.05 | $7.5\times 10^{-3}$ | $7.5\times 10^{-3}$
B | 1 | 0.05 | $1.2\times 10^{-1}$ | $2\times 10^{-3}$
C | 30 | 1 | $6\times 10^{-1}$ | $5\times 10^{-5}$

Model sets A and B are mainly determined by DM relic abundance since the dark
candidate mass is in the region where direct detection experiments are less
restrictive (Ipek et al., 2014). The couplings in set C are chiefly
constrained by LUX results (Akerib et al., 2017) in spin independent and spin
dependent cross sections and, in addition, they respect restrictive rare meson
decays (Dolan et al., 2015) as well. In the beforementioned cases we estimate
the parameters using MicroOmegas (Belanger et al., 2010) and direct detection
cross sections at one-loop level (Freytsis & Ligeti, 2011; Ipek et al., 2014).

According to the current stage of exploration of the phase space of masses and
cross sections for DM candidates interacting with nucleons in ordinary matter,
dense compact stars are believed to be suitable places to find this kind of
matter. NSs are believed to be efficient DM accretors (Gould, 1987). One of
the key quantities that can dictate their internal stellar energetic balance
is the local energy emissivity, $Q_{E}=\frac{dE}{dVdt}$ (energy produced per
unit volume per unit time, through a prescribed particle physics reaction). In
this work we will be interested in the annihilation reactions of DM into two-
body fermionic states ($f$), $\chi{\chi}\rightarrow f\bar{f}$ and two
pseudoscalar boson states $\chi{\chi}\rightarrow aa$ with subsequent decay
$a\rightarrow f\bar{f}$. Furthermore, we will discuss possible astrophysical
consequences particularizing to the $f=\nu$ neutrino channel.

Formally, the expression for $Q_{E}$ generically denotes the energy emission
rate per stellar volume arising from fermionic or pseudoscalar pair production
and can be written as (Esposito et al., 2002)

$Q_{E}=4\int
d\Phi(E_{1}+E_{2})\,|\overline{\mathcal{M}}|^{2}\,f(f_{1},f_{2},f_{3},f_{4}),$
(2)

with

$d\Phi=\frac{d^{3}\vec{p_{1}}}{2(2\pi)^{3}E_{1}}\frac{d^{3}\vec{p_{2}}}{2(2\pi)^{3}E_{2}}\frac{d^{3}\vec{p_{3}}}{2(2\pi)^{3}E_{3}}\frac{d^{3}\vec{p_{4}}}{2(2\pi)^{3}E_{4}}\,(2\pi)^{4}\delta(p_{1}+p_{2}-p_{3}-p_{4}),$
(3)

the 4-body ($12\rightarrow 34$) phase space element and
$|\mathcal{\overline{M}}|^{2}$, the spin-averaged squared matrix element of
the reaction considered. The additional factor $f(f_{1},f_{2},f_{3},f_{4})$
accounts for the global phase space blocking factor due to the initial and
final particle distribution functions, $f_{i},\,i=1,...,4$ we will discuss
below. $\delta(x)$ is the 4-dimensional delta function. We will denote
$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. The detailed associated Feynman diagrams
are shown in Fig.(1). Let us note that besides the quoted annihilation
processes we consider, there may be additional pseudoscalar boson production
s-wave $\chi{\chi}\rightarrow aaa$ (Abdullah et al., 2014), initial/final
state radiation and internal bremsstrahlung processes $\chi{\chi}\rightarrow
f\bar{f}a$ or $\chi{\chi}a\rightarrow f\bar{f}$ (Bell et al., 2017). However,
since the cross sections for these processes are proportional to
$g_{\chi}^{2}g_{f}^{4}$ and $g_{\chi}^{4}g_{f}^{2}$, respectively, they are
subdominant in the case of a Dirac fermion DM candidate (Bringmann et al.,
2017; Ibarra et al., 2013). Similarly, radiative $a$-production can arise from
the SM particles interaction inside the star, but this process is found to be
only relevant in the case of very light mediators ($\lesssim$ eV) like axions
or Majorons (Sedrakian, 2016; Farzan, 2003).

Figure 1: Feynman diagrams for DM annihilation reactions
$\chi{\chi}\rightarrow f\bar{f}$ and $\chi{\chi}\rightarrow aa$ considered in
this work.

Specifically, for the case of annihilation into fermionic pairs (left diagram
in Fig.(1)) we label the emissivity as $Q_{E}^{f\bar{f}}$. It includes the
expression for the spin-averaged squared matrix element as

$|\mathcal{\overline{M}}_{f\bar{f}}|^{2}=\frac{g_{\chi}^{2}g_{f}^{2}}{4}\frac{s^{2}}{(s-m_{a}^{2})^{2}+E_{|\vec{q}|}^{2}\Gamma^{2}},$
(4)

where $q^{2}=s=(p_{1}+p_{2})^{2}=(p_{3}+p_{4})^{2}$ is the Mandelstam variable
and $E_{|\vec{q}|}=\sqrt{|\vec{q}|^{2}+m_{a}^{2}}$. In this case

$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})),$
(5)

and $f_{\chi},f_{f}$ are the local stellar distribution functions for DM and
fermionic particles, respectively, containing density and temperature
dependence we will discuss further below. $\Gamma$ is the pseudoscalar
particle decay width in the local medium through the reaction $a\rightarrow
f\bar{f}$. It is obtained using the optical theorem as

$\Gamma=\frac{1}{E_{|\vec{q}|}}{\rm Im}\Pi(\vec{q}),$ (6)

where $\Pi(\vec{q})$ is the pseudoscalar polarization insertion given by

$\Pi(\vec{q})=\frac{ig_{f}^{2}}{2}\int\frac{d^{4}k}{(2\pi)^{4}}tr[\gamma_{5}G^{0}(k)\gamma_{5}G^{0}(k+q)],$
(7)

and the corresponding cut of the associated tadpole diagram involves the
fermion propagator $G^{0}(k)$ including a vacuum and matter contribution
(Chin, 1977; Matsui & Serot, 1982).

Using the Eq.(2) and Eq.(4) we can obtain an expresion for the emissivity
produced by the annihilation of DM particles into ${f\bar{f}}$,
$Q_{E}^{f\bar{f}}$. Let us first deal with the integration over $\vec{p_{4}}$
so that

$\int\frac{d^{3}\vec{p_{3}}}{2E_{3}(2\pi)^{3}}\int\frac{d^{3}\vec{p_{4}}}{2E_{4}(2\pi)^{3}}(2\pi)^{4}\delta(p_{1}+p_{2}-p_{3}-p_{4})=\int\frac{2\pi|\vec{p_{4}}|^{2}d|\vec{p_{4}}|}{4E_{3}E_{4}(2\pi)^{2}}d(\cos\,\theta)\delta(q_{0}-E_{3}-E_{4})$
(8)

where $\theta$ is the angle between $\vec{p_{4}}$ and $\vec{q}$, and
$q_{0}=E_{1}+E_{2}$ for this annihilation channel. Besides, we can express the
energy delta function as

$\delta(q_{0}-E_{3}-E_{4})=\frac{\sqrt{m_{f}^{2}+|\vec{q}|^{2}+|\vec{p_{4}}|^{2}-2|\vec{p_{4}}||\vec{q}|\cos\,\theta}}{|\vec{p_{4}}||\vec{q}|}\delta(\cos\,\theta-\cos\,\theta_{0})\Theta(q^{2}\geq
4m_{f}^{2}),$ (9)

where

$\cos\;\theta_{0}=\frac{1}{2|\vec{p_{4}}||\vec{q}|}\left(|\vec{q}|^{2}-q_{0}^{2}+2q_{0}E_{4}\right),$
(10)

and $q^{2}=q_{0}^{2}-|\vec{q}|^{2}$. $\Theta(x)$ is the Heaviside function.
Eq.(9) has been obtained using

$\delta[f(x)]=\sum_{i}\frac{\delta(x-x_{0i})}{|f^{\prime}(x)|_{x_{0}i}|},$
(11)

with $x_{0i}$ the zeros of $f(x)$. Now, we use that
$|\vec{p_{4}}|d|\vec{p_{4}}|=E_{4}dE_{4}$. Imposing $\cos^{2}\theta_{0}\leq
1$, we obtain limits for the integration over $E_{4}$

$E_{4}\pm=\frac{1}{2}\left(q_{0}\pm|\vec{q}|\sqrt{1-\frac{4m_{f}^{2}}{q^{2}}}\right).$
(12)

In the same way, we use $|\vec{p_{1}}|d|\vec{p_{1}}|=E_{1}dE_{1}$ and
$|\vec{p_{2}}|d|\vec{p_{2}}|=E_{2}dE_{2}$. After that, Eq.(2) takes the form

$\displaystyle Q_{E}^{f\bar{f}}$ $\displaystyle=$
$\displaystyle\frac{1}{2(2\pi)^{5}}\int_{m_{\chi}}^{\infty}dE_{1}\sqrt{E_{1}^{2}-m_{\chi}^{2}}\int_{m_{\chi}}^{\infty}dE_{2}\sqrt{E_{2}^{2}-m_{\chi}^{2}}\int_{-1}^{1}d(\cos\,\phi)q_{0}\Theta(q^{2}\geq
4m_{f}^{2})$ (13)
$\displaystyle\times\int_{E_{4}-}^{E_{4}+}dE_{4}\int_{-1}^{1}d\cos\theta\delta(\cos\,\theta-\cos\,\theta_{0})f(f_{1},f_{2},f_{3},f_{4})\frac{|\mathcal{\overline{M}}_{f\bar{f}}|^{2}}{|\vec{q}|}$

where $\phi$ is the angle between $\vec{p_{1}}$ and $\vec{p_{2}}$,
$E_{3}=E_{1}+E_{2}-E_{4},$

$|\vec{q}|=\sqrt{|\vec{p_{1}}|^{2}+|\vec{p_{2}}|^{2}+2|\vec{p_{1}}||\vec{p_{2}}|\cos\,\phi},$
(14)

and $|\vec{p_{1}}|=\sqrt{E_{1}^{2}-m_{\chi}^{2}}$,
$|\vec{p_{2}}|=\sqrt{E_{2}^{2}-m_{\chi}^{2}}$.

Instead, for the DM annihilation into two pseudoscalars (middle and right
diagrams in Fig.(1)) now the emissivity is labeled as $Q_{E}^{a{a}}$. In this
case, when calculating the spin averaged matrix element one must note that

$|\mathcal{\overline{M}}_{aa}|^{2}=\frac{1}{2}\frac{1}{2}\sum_{s,s^{\prime}}|\mathcal{M}_{aa}|^{2},$
(15)

where $s$ and $s^{\prime}$ are the spin states of the dark matter particle.
The squared matrix element finally reads,

$\displaystyle|\mathcal{\overline{M}}_{aa}|^{2}$ $\displaystyle=$
$\displaystyle\frac{-g_{\chi}^{4}}{2}\\{\frac{(t-m_{a}^{2})^{2}-m_{\chi}^{2}(m_{\chi}^{2}+2m_{a}^{2})}{(t-m_{\chi}^{2})^{2}}+\frac{(u-m_{a}^{2})^{2}-m_{\chi}^{2}(m_{\chi}^{2}+2m_{a}^{2})}{(u-m_{\chi}^{2})^{2}}+2\frac{2m_{\chi}^{2}-s}{u-m_{\chi}^{2}}$
(16) $\displaystyle+$
$\displaystyle\frac{(s-2m_{\chi}^{2})(2m_{a}^{2}-s)+2m_{\chi}^{2}(m_{\chi}^{2}+2m_{a}^{2}-2s)-2(t-m_{a}^{2})^{2}}{(t-m_{\chi}^{2})(u-m_{\chi}^{2})}\\},$

where $s=k^{2}=(p_{1}+p_{2})^{2}$, $t=(p_{1}-p_{3})^{2}=(p_{4}-p_{2})^{2}$ and
$u=2m_{\chi}^{2}+2m_{a}^{2}-s-t$.

As we can see, now the matrix element not only depends on $s$ but also on $t$
and $u$. Dealing with the integration to obtain the emissivity it is
convenient to write these variables as

$s=2m_{\chi}^{2}+2E_{1}E_{2}-2|\vec{p_{1}}||\vec{p_{2}}|\cos\,\theta_{12},$
(17)

and

$t=m_{\chi}^{2}+m_{a}^{2}-2E_{1}E_{3}+2|\vec{p_{1}}||\vec{p_{3}}|\cos\,\theta_{13},$
(18)

being $\theta_{ij}$ the angle between $\vec{p_{i}}$ and $\vec{p_{j}}$.

In the same way that we did for the annihilation into fermions we can write

$\int\frac{d^{3}\vec{p_{3}}}{2E_{3}(2\pi)^{3}}\int\frac{d^{3}\vec{p_{4}}}{2E_{4}(2\pi)^{3}}(2\pi)^{4}\delta(p_{1}+p_{2}-p_{3}-p_{4})=\int\int\int\frac{|\vec{p_{3}}|^{2}d|\vec{p_{3}}|d\phi_{3}}{4E_{3}E_{4}(2\pi)^{2}}d(\cos\,\theta_{3})\delta(k_{0}-E_{3}-E_{4}),$
(19)

where we are denoting the four momentum $k=(k_{0},\vec{k})$ and $\theta_{i}$
as the angle between $\vec{p_{i}}$ and $\vec{k}$. As obtained in Eq.(9), we
find

$\delta(k_{0}-E_{3}-E_{4})=\frac{E_{4}}{|\vec{p_{3}}||\vec{k}|}\delta(\cos\,\theta_{3}-\cos\,\theta_{3,0}),$
(20)

being

$\cos\,\theta_{3,0}=\frac{1}{2|\vec{p_{3}}||\vec{k}|}(|\vec{k}|^{2}+2k_{0}E_{3}-k_{0}^{2}).$
(21)

Now, we can write the emissivity into two pseudoscalars as

$\displaystyle Q_{E}^{a{a}}$ $\displaystyle=$
$\displaystyle\frac{1}{2(2\pi)^{6}}\int_{0}^{\infty}|\vec{k}|d|\vec{k}|\int_{m_{\chi}}^{\infty}\sqrt{E_{2}^{2}-m_{\chi}^{2}}dE_{2}\int_{-1}^{1}d\cos\,\theta_{2}\frac{E_{1}+E_{2}}{E_{1}}\Theta(k^{2}-4m_{a}^{2})$
(22)
$\displaystyle\times\int_{E_{3}-}^{E_{3}+}dE_{3}\int_{-1}^{1}d\cos\,\theta_{3}\delta(\cos\,\theta_{3}-\cos\,\theta_{3,0})\int_{0}^{2\pi}d\phi_{3}f(f_{1},f_{2},f_{3},f_{4})|\mathcal{\overline{M}}_{aa}|^{2}$

where
$|\vec{p_{1}}|=\sqrt{|\vec{k}|^{2}+|\vec{p_{2}}|^{2}-2|\vec{k}||\vec{p_{2}}|\cos\,\theta_{2}}$,
$t=m_{\chi}^{2}+m_{a}^{2}-2E_{1}E_{3}+2|\vec{p_{3}}|(|\vec{k}|\cos\,\theta_{3}-|\vec{p_{2}}|\cos\,\theta_{23})$
and
$\cos\,\theta_{23}=\cos\,\phi_{3}\,\sin\,\theta_{2}\,\sin\,\theta_{3}+\cos\,\theta_{2}\cos\,\theta_{3}$.
We have also used the trigonometric relation
$|\vec{p_{1}}|\cos\,\theta_{13}=|\vec{k}|\cos\,\theta_{3}-|\vec{p_{2}}|\cos\,\theta_{23}$.
The limits for the outgoing energy in the integral are

$E_{3}\pm=\frac{1}{2}\left(k_{0}\pm|\vec{k}|\sqrt{1-\frac{4m_{a}^{2}}{k^{2}}}\right).$
(23)

In the case of annihilation into pseudoscalars the phase space factor reads

$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}).$
(24)

In this case one should also take into account the further decay of each
pseudoscalar into fermionic pairs and the availability of kinematical phase
space through and additional Pauli blocking factor. Although not explicit,
there is also a further local dependence on the DM density in the distribution
function that will be discussed later in the manuscript.

###  2.1 Dense and hot stellar scenario

In order to explain the physical relevance of the quantities under scrutiny
obtained in the previous section, at this point we will particularize to that
of a dense and hot stellar scenario. We will focus on a (proto) NS. Briefly, a
NS is mostly constituted by nucleons forming a central core at a density in
excess of nuclear saturation density, $\rho_{0}\simeq 2.4\times 10^{14}\,\rm
g/cm^{3}$. An average NS has a radius $R\lesssim 12$ km and mass
$M\sim(1-1.5)M_{\odot}$ (mostly in its core) being thus a star with large
compactness ratio $\sim M/R$. For the sake of our discussion we will consider
a typical baryonic core density value $\rho_{b}=\rho_{N}\sim 2\rho_{0}$.
Regarding internal temperature and composition, NSs are born as hot lepton-
rich objects with temperatures $T\sim 20$ MeV evolving into cold $T\sim 10$
keV neutron-rich ones, after a deleptonization era. Assuming dark and ordinary
matter have coupling strengths at the level of current experimental search
bounds, NSs are believed to be capable of accreting (and retaining) DM
particles whose masses are larger than a few GeV from an existing galactic
distribution.

Accretion of a dark component will proceed not only during the collapsed stage
but also during most of the previous progenitor stellar lifetime at different
epoch-dependent capture rates, $C_{\chi}$. First, in the progenitor stages,
the progressively denser nuclear ash central core is effectively opaque to DM
and allows building up an internal finite DM number density over time,
$n_{\chi}(r)$, being $r$ the radial stellar coordinate. Briefly, the
progenitor with a mass $\sim(10-15)M_{\odot}$ is able to fuse lighter elements
into heavier ones and thus its composition changes through the burning ages.
Hydrogen first, and later the He, C, O, and rest of heavier elements up to Si
proceed through the burning stages. Spin-dependent (mostly from H) as well as
spin independent $\chi$-N cross sections allow the gravitational capture of DM
population inside the star. Coherence effects may play a role for slowly
moving, low ${m_{\chi}}$ incoming DM particles scattering nuclei off when
their associated de Broglie wavelength is comparable to the nuclear size, and
in this case the spin-independent cross section bears a multiplicative factor
$\simeq A^{2}$ where $A$ is the baryonic number. Since the later burning
stages proceed rapidly, the ${\rm He-C-O}$ stage gives the main contribution
to the DM capture in the progenitor. As the thermalization times during this
set of stages can follow the internal dynamics the collapsed star will have as
a result a non-zero, mostly inherited, initial DM population.

Most in detail, the DM particle population number inside the star, $N_{\chi}$,
will not only depend on the capture rate $C_{\chi}$ (Gould, 1987) but also on
the self-annihilation rate, $C_{a}$. Note that in the range of masses in the
parameter sets we consider, evaporation effects (Krauss et al, 1986) as well
as decay (Perez-Garcia & Silk, 2015) do not substantially modify the DM
population as the kinetic to gravitational potential energy ratio remains
small.

Then the DM particle number, $N_{\chi}$, can be obtained a function of time
$t$ by solving the differential equation

$\frac{dN_{\chi}}{dt}=C_{\chi}-C_{a}N^{2}_{\chi},$ (25)

considering the two competing processes, capture and annihilation (Kouvaris &
Tinyakov, 2010)

$N_{\chi}(t)=\sqrt{\frac{C_{\chi}}{C_{a}}}\rm\,tanh\left[\frac{t}{\tau}+\gamma(N_{\chi,0})\right],$
(26)

where

$\gamma(N_{\chi,0})={\rm
tanh}^{-1}\left(\sqrt{\frac{C_{\chi}}{C_{a}}}N_{\chi,0}\right)$ (27)

and

$\tau^{-1}=\sqrt{C_{\chi}C_{a}}.$ (28)

At $t=0$, when the protoNS is born, a typical progenitor may have already
provided an initial population

$N_{\chi,0}=1.5\times
10^{39}\left(\frac{\rho_{\chi}}{\rho^{ambient}_{\chi,0}}\right)\left(\frac{1\;\rm
GeV}{m_{\chi}}\right)\left(\frac{\sigma_{s}}{10^{-43}\;\rm cm^{2}}\right),$
(29)

where $\sigma_{s}\equiv\sigma_{\chi-N}$ is the $\chi-N$ scattering cross
section. As this quantity is currently unknown, only experimental constrains
exist for it. In the range of DM masses used in this work
$\sigma_{s}\in[10^{-46}-10^{-33}]\,\rm cm^{-2}$ (B. J. Kavanagh, 2017). Eq.
(29) assumes that the majority of the NS population can be found at
galactocentric distances of a few kpc where $\rho_{\chi}\sim
10^{2}\rho^{ambient}_{\chi,0}$. We use $\rho^{ambient}_{\chi,0}\simeq 0.3$
$\rm\frac{GeV}{cm^{3}}$ as the solar-circle DM density value.

Let us mention that both capture and annihilation rates, will be intimately
determined by the parameters of the model at hand, i.e.
$m_{\chi},m_{a},g_{0},g_{\chi}$ (we set $g_{f}=1$). In particular, the DM
capture rate on the progenitor depends on the scattering cross section on
nuclei (nucleons) that is proportional to the product of the couplings $(\sim
g^{2}_{0}g^{2}_{\chi})$ and the annihilation cross section proportional to the
sum of $(\sim g^{2}_{0}g^{2}_{\chi})$ and $(\sim g^{4}_{\chi})$ terms for the
two reactions considered $\chi{\chi}\rightarrow f\bar{f}$ and
$\chi{\chi}\rightarrow aa$, respectively (Buckley et al, 2015).

For the three models considered in this work appearing in Table 1, the average
progenitor capture rate allows a non-vanishing initial DM population,
$N_{\chi,0}$, since the annihilation rate, proportional to $n^{2}_{\chi}(r),$
is negligibly small at that stage. Later, in the NS collapsed state and at a
given galactic location with a corresponding ambient DM density $\rho_{\chi}$,
the capture rate, $C_{\chi}$, it is approximated up to factors of order unity
by the expression (Gould, 1987; de Lavallaz et al., 2010).

$C_{\chi}\simeq 1.8\times 10^{25}\left(\frac{1\,\rm
GeV}{m_{\chi}}\right)\left(\frac{\rho_{\chi}}{\rho^{ambient}_{\chi,0}}\right)f_{\chi,N}\,\,\rm
s^{-1}.$ (30)

A few remarks are due regarding this expression. $f_{\chi,N}$ denotes a
phenomenological factor dealing with the opacity of stellar matter.
$f_{\chi,N}$ depends on the ratio of the leading contribution of $\chi-N$
scattering cross section $\sigma_{s}$ to the minimum geometrical cross section
of a NS made of nucleons of mass $m_{N}$ and defined as

$\sigma_{0}=\frac{m_{N}R^{2}}{M}\sim 10^{-45}\rm cm^{2}.$ (31)

Thus this factor saturates to unity, $f_{\chi,N}\sim 1$, if
${\sigma_{s}}\gtrsim{\sigma_{0}}$. Otherwise,
$f_{\chi,N}\sim\frac{\sigma_{s}}{\sigma_{0}}$. Using Boehm et al. (2014) and
appendix D in Dolan et al. (2015) we consider the scattering cross section (at
one-loop) in the appropriate kinematical limit in our compact star so that for
the parameters used in this work $f_{\chi}$ is effectively in the saturated
regime.

The expressions in the literature for DM capture rates in dense objects are
based on interactions with quarks (nucleons) that in practice happen via a
contact term (Freytsis & Ligeti, 2011), possibly including form factors. Note
that a more ellaborate treatment would involve the calculation of the non-
relativistic limit of the (full) series of operators included in the
lagrangian under study. Such a detailed analysis is beyond the scope of this
work and remains to be done. The usual phenomenological treatment, through the
$f_{\chi,N}$ factor, makes use of a lower bound to the global cross section
with all relevant contributions in this realization picture. It is important
to emphasize that the strength of the computed emissivity will depend on the
number of dark matter particles remaining inside in the star at any given
time.

As thermalization times for DM particles in the light mass range we consider
are consistently smaller than dynamical cooling times (Goldman & Nussinov,
1989), inside the star the DM particle number density takes the form

$n_{\chi}(r)=n_{0,\,\chi}e^{-\frac{m_{\chi}}{k_{B}T}\Phi(r)},$ (32)

with $n_{0,\,\chi}$ the central value, $T$ the NS temperature and $k_{B}$ the
Boltzmann constant. The gravitational potential is given by

$\Phi(r)=\int_{0}^{r}\frac{GM(r^{\prime})dr^{\prime}}{{r^{\prime}}^{2}},$ (33)

where $M(r^{\prime})$ is the NS mass inside a spherical volume of radius
$r^{\prime}$. So that assuming an approximately constant density core,

$n_{\chi}(r)=n_{0,\,\chi}e^{-(r/r_{\rm th})^{2}},$ (34)

with a thermal radius

$r_{\rm th}=\sqrt{\frac{3k_{B}T}{2\pi G\rho_{N}m_{\chi}}}.$ (35)

Normalization requires $\int_{0}^{R}n_{\chi}(r)dV=N_{\chi}$ at a given time,
as reflected by Eq.(26). Note that potential limiting values of $N_{\chi}$ may
arise from the fact that a fermionic $\chi$ would involve the existence of a
Chandrasekhar critical mass for collapse (MacDermott et al., 2012). This
possibility is safely not fulfilled as long as $N_{\chi}(t)<N_{\rm Ch}$, where
$N_{\rm Ch}\sim(\,M_{\rm Pl}/m_{\chi})^{3}\sim 1.8\times
10^{57}\,(1{\rm\,GeV}/m_{\chi})^{3}$ with $M_{\rm Pl}$ the Planck mass.

Let us now comment on the fact that inside the star the tiny DM fraction can
be described by a distribution function of a classical Maxwell-Boltzmann type

$f_{\chi}=f^{MB}_{\chi}(|\vec{p_{i}}|,r)=\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}},\,i=1,2.$
(36)

and in the non-relativistic scheme
$E_{i}=\frac{\vec{p_{i}}^{2}}{2m_{\chi}},\,\rm i=1,2$. The annihilation rate,
$C_{a}$, depends on the thermally averaged annihilation cross section inside
the star, $\langle{\sigma_{a}v}\rangle$, for the two reactions considered in
this work, see Fig.(1). Therefore, the stellar $\chi$-distribution contained
in the thermal volume region $\sim r^{3}_{\rm th}$ determines the annihilation
rate $C_{a}\sim\langle{\sigma_{a}v}\rangle/V_{\rm th}$ (Goncalves et al.,
2017; Arina et al., 2015). Note that the presence of the phase space factor
$f(f_{1},f_{2},f_{3},f_{4})$ in Eq. (2) will introduce further DM density and
$T$ dependence into the vacuum standard calculation as a thermalized DM
distribution exists inside the NS core. As for the outgoing fermions, the
medium density effects will generally arise from the phase space blocking
factors and collective effects (Cermeño et. al, 2016). In case of neutrinos we
assume $f_{\nu}\sim 0$, although in cases where a trapped fraction $Y_{\nu}>0$
exists it would further decrease the response.

Figure 2: Energy emissivity from DM self-annihilation channel
$\chi\chi\rightarrow\nu\nu$ as a function of temperature. Standard emission
refers to MURCA processes. $N_{\chi}=N_{0,\chi}$ is assumed. See text for
details.

##  3 Results

In this section we explain our results regarding the emissivities in the NS
astrophysical scenario particularizing to the case in which the final state
fermions produced in the reactions depicted in Fig.(1) are neutrino pairs.
Neutrinos are weakly interacting SM fermions known to play a key role in the
internal energy dynamics of a massive stellar progenitor undergoing
gravitational collapse. In such an event most of the gravitational binding
energy is emitted into neutrinos (and antineutrinos) of the three families. A
very efficient cooling scenario emerges in the first $\sim 10^{5}$ yr.
Standard processes such as those present in the URCA or the modified URCA
(MURCA) cooling (Friman & Maxwell, 1979; Yakovlev et al., 2005) among others
can release neutrinos with associated emissivities $Q^{\rm URCA}_{E}\sim
10^{27}\mathcal{R}(\frac{T}{0.1\,\rm MeV})^{6}$ $\rm erg\,cm^{-3}\,s^{-1}$ and
$Q^{\rm MURCA}_{E}\sim 10^{21}\mathcal{R}(\frac{T}{0.1\,\rm MeV})^{8}$ $\rm
erg\,cm^{-3}\,s^{-1}$, respectively. Typical energetic scales can be obtained
from the conversion factor $1$ MeV $\sim 10^{10}$ K. $\mathcal{R}$ is a
reduction function of order unity describing the superfluid effects in the
neutron and proton branches of those reactions (Yakolev & Levenfish, 1995). We
must keep in mind the fact that these neutrinos effectively cool off the star
as they leave, having scattered a few times with ordinary nucleon matter
(Horowitz & Perez-Garcia, 2003; Perez-Garcia, 2010) after a first rapid
trapping stage. In this way processes with neutrino production in reactions
involving nucleon components effectively release energy from the baryonic
system as the associated neutrino mean free path is relatively long
$\lambda\sim 28\,\rm cm\,(100\,\rm MeV/E_{\nu})^{2}$. In analogy with what
happens at the standard neutrino trapping stage in very young stars, when
neutrinos have energies of dozens of MeV, it is expected that energetic
neutrinos produced in reactions of DM annihilation could have mean free paths
very small, even at low stellar temperatures (in evolved stars) so they may
not escape so easily the dense medium (Kouvaris, 2007).

Figure 3: Energy emissivity from DM self-annihilation channel
$\chi\chi\rightarrow aa$ with subsequent decay $a\rightarrow\nu\nu$ as a
function of temperature. Standard emission refers to MURCA processes.
$N_{\chi}=N_{0,\chi}$ is assumed. See text for details.

In Fig.(2) the logarithm (base 10) of the energy emissivity for the process
$\chi\chi\rightarrow\nu\bar{\nu}$ is shown as a function of temperature for
the three sets of DM parameters A, B, C in Table 1, with dashed, dotted and
dash-dotted lines, respectively. Baryonic density is fixed at $2\rho_{0}$. We
also fix (for reference) the assumed number of DM particles at initial instant
$N_{\chi}=N_{0,\chi}$. Note, however, that $N_{\chi}$ is time (and model)
dependent as it rapidly decreases when the self-annihilation sets in. We give
below a suitable fit where the actual T and $N_{\chi}$ dependence is
reflected. The baryonic density dependence is, however, weak. The standard
physics cooling is depicted here by the MURCA emissivity (solid line) for the
sake of comparison. Although the latter is not the only process that could
possibly contribute to the effective standard cooling, all other processes
capable of are considered weaker at the temperature and density conditions
considered in this work. Thus we take the standard processes to be represented
by an upper limit to the currently used emissivities chosen as the MURCA
processes. Note we do not consider exotic meson codensates nor URCA
emissivities either, since stellar central densities required are usually
higher than the one taken as reference here in order to provide the $Y_{p}\sim
11-15\%$ proton fraction (Lattimer et al., 1991) to sustain the fast reaction.
Further, we assume that a possible neutrino trapping phase will not be
significant, particularly at low T. However, at high T it may produce a
further reduction that has to be accounted for through a Pauli blocking
factor. It must be included to take into account the time dependent non-
vanishing leptonic fraction of order $Y_{L}\sim 0.1$ (Pons et al., 1999) until
transparency sets in at $\sim 20$ s.

The trends depicted in Fig.(2) with temperature and DM particle population
dependence for models A, B and C and for both reaction channels can be fit as

$Q_{E}(T,N_{\chi})=Q_{0}\left(\frac{N_{\chi}}{N_{0,\chi}}\right)^{2}\left(\frac{T}{1\,\rm
MeV}\right)^{-3}.$ (37)

In Table 2 we give values for parameters $Q_{0}$ and $N_{0,\chi}$ for the
reactions and models considered. We can see that around $T\sim 0.1$ MeV
standard emissivities $\rm log_{10}\,(Q_{E})\sim 21$ are as powerful as those
from the DM annihilation processes in the thermal volume region $V_{\rm
th}\sim r^{3}_{\rm th}$. However, since DM population is a decreasing function
of time, one can expect that there will be a minimal number of $N_{\chi}$
population to beat the MURCA processes for $T<0.1$ MeV. We find
$Q_{E}(T,N_{\chi})>Q_{\rm MURCA}$ for $T\in[0.01,0.1]$ $\rm MeV$ for
$N_{\chi}/N_{0,\chi}\gtrsim 10^{-5},3.6\times 10^{-6},5.6\times 10^{-6}$ for
A, B and C models, respectively. One should note that the $N_{\chi}$ self-
consistently depends on temperature and how it dynamically changes with time.
A fully detailed cooling simulation would yield the temporal sequence to
determine the complete behaviour. As this is not the goal here, we give
instead an estimate on the time duration of the dominance of the DM
annihilation channel, i.e. where it could beat the local MURCA processes from
Eq.(26). We obtain $t\lesssim 50$ s for models A and B while this condition is
true at all times for model C. We note that all timescales are thus
overlapping the standard transparency window for SM neutrinos.

In Fig.(3) the logarithm (base 10) is shown for the reaction
$\chi\chi\rightarrow aa$ with subsequent decay $a\rightarrow\nu\nu$. The
number of DM particles is also fixed $N_{\chi}=N_{0,\chi}$ for reference using
the same argument as with the $\chi\chi\rightarrow\nu\nu$ reaction in Fig.(2).
In this case, for a population $N_{0,\chi}$, the neutrino emissivity is
largely enhanced with respect to the direct production of neutrinos
$\chi\chi\rightarrow\nu\nu$. For model C (dash-dotted line) $Q_{E}$ matches
and surpasses the standard MURCA emission below $T\sim 0.5$ MeV while for
models B (dotted line) and A (dashed line) that happens for $T\sim 0.3$ MeV
and $T\sim 0.1$ MeV, respectively. Values for the phenomenological fit in this
channel are also provided in Table 2.

If we now consider the running character of DM population number as in the
previous case we find $Q_{E}(T,N_{\chi})>Q_{\rm MURCA}$ in the interval
$T\in[0.01,0.1]$ $\rm MeV$ for $N_{\chi}/N_{0,\chi}\gtrsim 3\times 10^{-5}$
and model A, while for model B $Q(T,N_{\chi})>Q_{\rm MURCA}$ in
$T\in[0.01,0.3]$ $\rm MeV$ when $N_{\chi}/N_{0,\chi}\gtrsim 5\times 10^{-8}$.
Finally $Q(T,N_{\chi})>Q_{\rm MURCA}$ in $T\in[0.01,0.5]$ $\rm MeV$ for
$N_{\chi}/N_{0,\chi}\gtrsim 3\times 10^{-10}$ for model C. This last ratio for
model C is achieved during the entire lifetime of the star while not for the
other cases. For the bi-pseudoscalar production reaction and at the
thermodynamical conditions in the scenario considered we obtain that typical
values obtained for the $a-$decay length are of order of $\sim 10^{2}$ fm,
making it a negligible contribution to the neutrino transport as their decay
length is so tiny compared to stellar size.

Table 2: Parameters obtained for the fit in Eq. (37).  Channel | Model | $\rm log_{10}\,Q_{0}$ $[\rm erg\,cm^{-3}\,s^{-1}]$ | $N_{0,\chi}$
---|---|---|---
$\chi\chi\rightarrow\nu\nu$ | A | 17.3 | $4.1\times 10^{44}$
$\chi\chi\rightarrow\nu\nu$ | B | 18 | $2.4\times 10^{44}$
$\chi\chi\rightarrow\nu\nu$ | C | 17.6 | $5\times 10^{38}$
$\chi\chi\rightarrow aa\tablenotemark{*}$ | A | 18 | $4.1\times 10^{44}$
$\chi\chi\rightarrow aa$ | B | 22.5 | $2.4\times 10^{44}$
$\chi\chi\rightarrow aa$ | C | 27 | $5\times 10^{38}$

At this point it is worth noting that tighter restrictions in the validity of
the coupling of DM to u-d-s quarks coming from including complementary
experimental bounds, e.g. rare meson decays, could somewhat reduce the
validity of models A and B. In addition, from isospin considerations, the
$\chi$ coupling to neutrons and protons in the NS will also affect somewhat
the results as the proton-to-neutron ratio inside the NS core is smaller than
unity ($\sim 1/9$). As the ratio of coupling fulfills $|g_{p}/g_{n}|>1$ when
considering flavour universal fermion couplings (Dolan et al., 2015) we expect
that the computed emissivities could be increased by a factor $\sim 10-100$.

We have considered the local emissivities of the novel reaction involving
self-annihilating matter in the thermal volume. We must emphasize that since
there is no uniform distribution of DM inside the star, neutrinos being
produced from nuclear reactions in the majority of the core volume will wash
out the dark contribution early, when temperature is high enough. However,
later, as temperature decreases there is an effective competition of the very
efficient dark central engine (located in a few percent of the core volume)
and the colder core emission.

The radial extent of the DM annihilating inner region is correlated to the
ratio of the thermal radius to the NS radius. It is defined as
$\xi=\sqrt{2}r_{\rm th}/R$ and indicates the radial fraction where DM
particles can be found. Since the crust region has a tiny mass we will not
consider this refinement here (Cermeño et al., 2016). For the parametrizations
A, B, C analyzed in this work this ratio takes values e.g. $\xi\in[0.03,0.42]$
at $T\sim 1$ MeV, $\xi\in[0.007,0.11]$ at $T\sim 0.1$ MeV and
$\xi\in[0.003,0.04]$ at $T\sim 0.01$ MeV. The volume where the dark emitting
region resides shrinks as $\sim\sqrt{T}$.

As thoroughly studied, enhanced emissivities in the medium can have an impact
on internal temperatures, temporal cooling sequence and (un)gapped matter
phases (Page et al., 2013, 2006). In this regard, recent works (Page et al.,
2011; Shternin et al., 2011) 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. In addition, current observations of thermal relaxation of NS crusts
indicates that even a small stellar volume fraction where fast neutrino
emission reactions can take place would provide distinctive features. More in
detail, it has been shown (Brown et al., 2018) that even if there is a
relatively small local volume where distortion from the standard energetic
mechanisms is taking place inside the star, a fast neutrino reaction (Direct
URCA) in a volume of $\sim 1\%$ could explain the neutrino luminosities in the
cooling curve of some particular objects like MXB 1659-29.

In the case presented in this work the long-term dark engine reaction
$\chi\chi\rightarrow aa$ could provide emissivities $Q_{E}^{aa}\sim 4\times
10^{22}$ $\rm erg\,cm^{-3}\,s^{-1}$ which result higher than that estimated
for the $Q_{MX1659-29}\sim 1.7\times 10^{21}$ $\rm erg\,cm^{-3}\,s^{-1}$for
$T\sim 10^{8}$ K. The stellar volume affected in the annihilating DM mechanism
is, nevertheless, much smaller for this range of temperatures but still
providing the same powerful emission. Besides the process discussed in this
work, other ones such as rotochemical heating (Fernández & Reisenegger, 2005)
or hot blobs located at different depths in the crust in young NS (Kaminker et
al., 2014) have been also treated in the literature adding more sources of
energetic variability based on SM matter.

It remains for further work to explore more exhaustively the precise relation
between the model parameters $m_{\chi},m_{a},g_{0},g_{\chi}$ (and $g_{f}$) and
the energetic efficiency of the emission. This will impact the duration of the
dominance of such emission from DM annihilation over standard processes and
thus its potential observability. We believe that the qualitatively different
picture arising from the DM self-annihilation process inside NS may be
worthwhile to explore.

##  4  Conclusions

We have calculated the energy emissivity of self-annihilating dark matter from
an existing stellar distribution into final state SM fermions. A pseudoscalar-
mediated DM interaction with the ordinary nucleon matter has been used. Later,
as an astrophysically relevant case we have particularized to neutrinos as
final states, and we have considered in detail those produced from s-wave
channels $\chi\chi\rightarrow\nu\nu$ or via pseudoscalar mediators p-wave
$\chi\chi\rightarrow aa$, and subsequent decay $a\rightarrow\nu\nu$. In the
inner stellar regions the radiation engine can encompass about $\lesssim 7\%$
of the total stellar volume for $T\lesssim 10^{10}$ K and the energy
emissivity can be enhanced orders of magnitude compared to the MURCA standard
neutrino processes for parameter sets respecting constraints of direct
detection limits, cosmological bounds or even tighter rare meson decay bounds.
We have provided a phenomenological fit of emissivities including dependence
of temperature and DM particle number. Taking as reference the usual standard
temporal sequence of NS cooling behaviour we expect that, for the models
analyzed in this work, model C (with $m_{\chi}=30$ GeV, $m_{a}=1$ GeV) could
be effectively active during the whole life of the star. Although a detailed
solution of the full evolution equation is out of the scope of this work it is
reasonable to foresee that the contribution of this new dark mechanism to the
set of already known standard cooling reactions will drive the star into
internal dynamical self-adjustment that is likely to emerge with a distinctive
temperature sequence whose observability remains to be properly analyzed in
future works.

We thank useful discussions with J. Silk, J. Edsjö, H. Grigorian and C.
Albertus. This work has been supported by PHAROS Cost action, MINECO
Consolider-Ingenio Multidark CSD2009-00064, Junta de Castilla y León SA083P17
and FIS2015-65140-P projects. M. Cermeño is supported by a fellowship from the
University of Salamanca.

## References

  * Bertone & Hooper (2016) Bertone, G., Hooper, D. , arXiv:1605.04909v2 [astro-ph.CO]
  * Chuzhoy (2007) Chuzhoy, L., arXiv:0710.1856 [astro-ph]
  * Mapelli et al. (2006) Mapelli, M., Ferrara, A. and Pierpaoli, E., Mon. Not. Roy. Astron. Soc. 369 (2006) 1719
  * Kouvaris (2007) Kouvaris,C., Phys. Rev. D 77 (2008) 023006
  * Rott et al. (2015) Rott, C. et al., JCAP11(2015) 039
  * Vincent et al. (2015) Vincent, A. C., Scott, P. and Serenelli, A., Phys. Rev. Lett. 114 (2015) 081302
  * Aartsen et al. (2012) Aartsen, M. G. et al., IceCube Collaboration, arXiv:1612.05949 [astro-ph.HE].
  * Palomares & Pascoli (2008) Palomares-Ruiz, S., Pascoli, S., Phys. Rev. D 77 (2008) 025025
  * Gusakov et al. (2004) Gusakov, M.E., Kaminker, A.D., Yakovlev, D.G., and Gnedin, O.Y. , Astron. Astrophys 423 (2004) 1063.
  * Yuksel & Beacom (2007) Yuksel, H. and Beacom, J. F., Phys. Rev. D 76 (2007) 083007
  * Vigano et al. (2013) D. Vigano et al., MNRAS 434, (2013) 123
  * Yakolev & Levenfish (1995) Yakovlev, D. G. and Levenfish, K. P., A&A 297 (1995) 717
  * Page et al. (2004) Page, D. , Lattimer, J. M., Prakash, M. and Steiner, A. W. , ApJ 155 (2004) 623
  * Boehm et al. (2014) Boehm, C., Dolan, M. J., McCabe, C., Spannowsky, M. and Wallace, C. J, JCAP 1405 (2014) 009.
  * Wild (2016) Wild, S., Ph. D. Thesis, Munich, Tech. U. (2016)
  * Bauer et al. (2017) Bauer, M., Haisch, U. and Kahlhoefer, F., arXiv: 1701.07427
  * Baek et al. (2017) Baek, S., Ko, P. and Li, J., arXiv: 1701.04131
  * Banerjee et al. (2017) Banerjee, S., Barducci, D., Bélanger, G., Fuks, B., Goudelis, A., and Zaldivar, B, JHEP 1707 (2017) 080
  * Arina et al. (2015) Arina, C. , Del Nobile, E. and Panci, P., Phys. Rev. Lett. 114 (2015) 011301
  * Abdullah et al. (2014) Abdullah, M., et al., Phys. Rev. D 90 (2014) 035004, arXiv:1404.6528 [hep-ph].
  * Haber & O’neil (2011) Haber, H. E. and O’Neil, D., Phys. Rev. D83 (2011) 055017
  * Ade et al. (2015) P. A. R. Ade et al. [Planck Collaboration], A&A 594, A13 (2016), arXiv:1502.01589 [astro-ph.CO].
  * Freytsis & Ligeti (2011) Freytsis, M. and Ligeti, Z., Phys. Rev. D. 83 (2011) 115009
  * Cheng & Chiang (2012) Cheng, H. and Chiang, C., JHEP 07 (2012) 009
  * Gresham & Zurek (2014) Gresham, M. and Zurek, K. M. , Phys. Rev. D 89 12 (2014) 123521
  * Ipek et al. (2014) Ipek, S., McKeen, D., and Nelson, A. E., Phys. Rev. D 90 (2014) 055021
  * Bertone et al. (2005) Bertone, G., Hooper, D. and Silk, J., Phys. Rep. 405 (2005) 279
  * Zel’dovich et al. (1965) Zel’dovich, Y. B. , Zh. Eksp. Teor. Fiz. 48 (1965), 986; Zel’dovich, Y. B. , Okun, L. B. and Pikelner, S. B. , Usp. Fiz. Nauk. 84 (1965), 113
  * Goncalves et al. (2017) Goncalves, D., Machado, P. A. N. and No, J. M., Phys. Rev. D 95 (2017) 055027
  * Dolan et al. (2015) Dolan, M. J., Kahlhoefer, F., McCabe, C. and Schmidt-Hoberg, K., JHEP 1503 (2015) 171
  * Akerib et al. (2017) Akerib, D. S. et al., Phys. Rev. Lett. 118 (2017) 021303
  * Belanger et al. (2010) Belanger, G., Boudjema, F., Pukhov, A. and Semenov, A., arXiv:1005.4133v1
  * Gould (1987) Gould, A., ApJ 321 (1987) 571
  * Esposito et al. (2002) Esposito, S., Mangano, M., Miele, G., Picardi, I. and Pisanti, O., Mod. Phys. Lett. A 17 (2002) 491
  * Bell et al. (2017) Bell, N.F. , Cai, Y., Dent,J. B., Leane, R. K. and Weiler, T. J., Phys. Rev. D 96 (2017) 023011
  * Bringmann et al. (2017) Bringmann, T., Calore, F., Galea, A., Garny, M., JHEP 1709 (2017) 041
  * Ibarra et al. (2013) Ibarra, A., Totzauer, M. and Wild, S., JCAP12 (2013) 043
  * Sedrakian (2016) Sedrakian, A., Phys. Rev. D 93 (2016) 065044
  * Farzan (2003) Farzan, Y., Phys.Rev. D 67 (2003) 073015
  * Chin (1977) Chin, S. A. , Annals Phys. 108 (1977) 301
  * Matsui & Serot (1982) Matsui, T. and Serot, B. D., Annals Phys. 144 (1982) 107
  * Krauss et al (1986) Krauss, L. M., Srednicki, M. and Wilczek, F., Phys. Rev. D 33 (1986) 2079
  * B. J. Kavanagh (2017) B. J. Kavanagh, arXiv:1712.04901v3 [hep-ph]
  * Perez-Garcia & Silk (2015) Perez-Garcia, M. A. and Silk, J., Physics Letters B 744 (2015) 13
  * Kouvaris & Tinyakov (2010) Kouvaris, C. and Tinyakov, P., Phys. Rev. D 82 (2010) 063531
  * Kouvaris & Tinyakov (2011) Kouvaris, C. and Tinyakov, P., Phys. Rev. D 83 (2011) 083512
  * Buckley et al (2015) Buckley, M. R., Feld, D. and Gonçalves, D., Phys. Rev. D 91 (2015) 015017
  * de Lavallaz et al. (2010) de Lavallaz, A. and Fairbairn, M., Phys.Rev. D 81 (2010) 123521
  * Goldman & Nussinov (1989) Goldman, I. and Nussinov, S., Phys. Rev. D 40 (1989) 3221
  * MacDermott et al. (2012) MacDermott, S. D., Yu, H.B. and Zurek, K. M., Phys. Rev. D 85 (2012) 023519
  * Cermeño et. al (2016) Cermeño, M., Perez-Garcia, M. A. and Silk, J., Phys. Rev. D 94 (2016) 023509
  * Friman & Maxwell (1979) Friman, B. L. and Maxwell, O. V., ApJ 232 (1979) 541
  * Horowitz & Perez-Garcia (2003) Horowitz, C. J. and Perez-Garcia, M. A., Phys. Rev. C 68 (2003) 025803
  * Perez-Garcia (2010) Perez-Garcia, M. A., Eur. Phys. J. A 44 (2010) 77
  * Cermeño et al. (2017) Cermeño, M., Pérez-García, M. A. and Silk, J., Publications of the Astronomical Society of Australia, 34 (2017)
  * Yakovlev et al. (2005) Yakovlev, D. G. et al., Nucl. Phys. A 752 (2005) 590
  * Lattimer et al. (1991) Lattimer, J. M., Pethick, C. J., Prakash, M. and Haensel, P., Phys. Rev. Lett. 66 (1991) 2701
  * Pons et al. (1999) Pons,J., Reddy, J., Prakash, M., Lattimer, J.M., Miralles, J.A., Astrophys. J. 513 (1999) 780
  * Cermeño et al. (2016) Cermeño, M., Perez-Garcia, M. A. and Silk, J., Phys. Rev. D 94 (2016) 063001
  * Page et al. (2013) Page, D., Lattimer, J. M. , Prakash, M. and Steiner, A. W., arXiv:1302.6626v2 [astro-ph.HE]
  * Page et al. (2011) Page, D., Prakash, M., Lattimer, J. M. and Steiner, A. W. , Phys. Rev. Lett.106, 081101 (2011)
  * Shternin et al. (2011) Shternin, P. S., Yakovlev, D. G. , Heinke, C. O., Ho, W. C. G. and Patnaude, D. J., Mon. Not. Roy. Astron. Soc. 412 (2011) L108
  * Brown et al. (2018) E. F. Brown et al., Phys. Rev. Lett. 120, 182701 (2018)
  * Page et al. (2006) Page, D., Reddy, S., Annu. Rev. Nucl. Part. Sci. 56 (2006) 327
  * Shternin & Yakovlev (2007) Shternin, P. S. and Yakovlev, D. G., arXiv: 0705.1963 [astro-ph]
  * Fernández & Reisenegger (2005) Fernández, R. and Reisenegger, A., ApJ, 625 (2005) 291
  * Kaminker et al. (2014) Kaminker, A. D., Kaurov, A. A., Potekhin, A. Y. and Yakovlev, D. G., MNRAS 442, (2014) 3484