Exploring dark sector parameters in light of neutron star temperatures

Using neutron stars (NS) as a dark matter (DM) probe has gained broad attention recently, either from heating due to DM annihilation or its stability under the presence of DM. In this work, we investigate spin-$1/2$ fermionic DM $\chi$ charged under the $U(1)_{X}$ in the dark sector. The massive gauge boson $V$ of $U(1)_{X}$ gauge group can be produced in NS via DM annihilation. The produced gauge boson can decay into Standard Model (SM) particles before it exits the NS, despite its tiny couplings to SM particles. Thus, we perform a systematic study on $\chi\bar{\chi}\to2V\to4{\rm SM}$ as a new heating mechanism for NS in addition to $\chi\bar{\chi}\to2{\rm SM}$ and kinetic heating from DM-baryon scattering. The self-trapping due to $\chi V$ scattering is also considered. We assume the general framework that both kinetic and mass mixing terms between $V$ and SM gauge bosons are present. This allows both vector and axial-vector couplings between $V$ and SM fermions even for $m_V\ll m_Z$. Notably, the contribution from axial-vector coupling is not negligible when particles scatter relativistically. We point out that the above approaches to DM-induced NS heating are not yet adopted in recent analyses. Detectabilities of the aforementioned effects to the NS surface temperature by the future telescopes are discussed as well.

DM self-interaction can be understood phenomenologically as an exchange of the U (1) X gauge boson V between dark matter particles. Assuming DM χ is a spin-1/2 fermion carrying U (1) X dark charge g d , its interaction with V is given by the Lagrangian where the associated DM self-interaction cross section σ χχ can be calculated from L DM and constrained by Eq. (1). Here we consider the scenario of symmetric dark matter [53] where the numbers χ andχ are identical. It has been proposed that the vector boson V in the dark sector can mix with SM photons and Z bosons through kinetic [54][55][56][57][58][59] and mass mixing terms [60][61][62]; the latter generally arise from extended Higgs sectors. Ref. [61] provides an example by introducing two Higgs doublet Φ 1 , Φ 2 , with Φ 1 coupled to SM fermions, and one scalar Higgs singlet φ. Both Φ 2 and φ carry the dark charge g d and the mixing among Φ 1 , Φ 2 and φ are neglected for simplicity. Hence neither Φ 2 nor φ couples to SM fermions.
The vacuum expectation values of Higgs scalars give mass terms for Z, V , and their mixing, which, together with kinetic mixing terms, are given by where B µν ≡ ∂ µ B ν − ∂ ν B µ is the U (1) Y field strength in SM while ε γ and ε Z are the kinetic and V − Z mass mixing parameters respectively. It is important to note that we do not invoke spontaneous symmetry breaking for generating DM mass as can be seen from Eq. (2) [63]. Since none of the scalar fields mentioned above couple to χ, the only mediator between the dark and visible sector is the dark boson V . Explicitly speaking, the electromagnetic (EM) and neutral-current (NC) interactions between V and SM fermions f resulting from mixing terms in Eqs. (3) and (4) are given by where g 2 is the SU (2) L coupling and J EM µ and J NC µ are the SM electromagnetic and neutral currents, respectively. The coefficientε Z is a linear combination of two mixing parameters and it reduces to ε Z for m V m Z . Its general expression is given in Appendix A.
In this paper, we examine the effect of DM heating due to the above phenomenological setup for a nearby three giga-year-old (Gyr-old) and isolated NS. The associated temperature is around 100 K according to the standard cooling mechanism if there is no other heating source. Therefore, any temperature deviation from this benchmark value can be potentially due to DM annihilation in the star. DM annihilation channels in this regard include not only χχ → ff but also χχ → 2V , provided m V < m χ and the decay length of V → ff is smaller than the radius of the star. Since V couples to neutral scalar bosons to acquire its own mass, the annihilation process χχ → V * → V +s with s one of the neutral scalar bosons is also possible for m V + m s < 2m χ . Such a process yields comparable heating effect to that given by χχ → 2V but involves an additional mass parameter m s . For simplicity in our discussions, we shall not consider this kinematic region. We refer the readers to Refs. [64,65] for the phenomenology of such an annihilation channel.
Searching the nearby old and cold NS can improve our understanding about DM. The new dynamics emerging from the above phenomenological setup will be discussed in the following sections. For completeness, we also analyze the signal to noise ratio (SNR) in the James Webb Space Telescope (JWST) [66]. Future telescopes such as European Extremely Large Telescope (E-ELT) and Thirty-Meter Telescope (TMT) [68] will constrain DM properties with unprecedented sensitivities. In the following sections, we employ the NS mass M 0 = 1.4M and and the radius R 0 = 12 km. We also replace g d with α χ = g 2 d /4π and all equations are expressed in terms of natural units = c = k B = 1.

II. DM CAPTURE AND NS TEMPERATURE
When a NS swipes through space, the DM particles in the halo can scatter with the baryons and leptons inside the star. Once DM loses an appreciable fraction of kinetic energy, it will be gravitationally captured by the NS. This capture process has been investigated extensively with contributions from neutrons, protons and leptons as well as relativistic corrections included in Refs. [69,70]. In this paper, only neutron contribution to the capture rate C c is considered. Contributions from other particle species are ignored due to their small yields. The DM number N χ in the star satisfies the differential equation while the anti-DM number Nχ evolves according to Here C a is the DM annihilation rate. Both coefficients C c and C a are well studied and the expressions can be found in Refs. [23,69,70] and references therein. We do not reproduce here. Thus, the exact solutions to Eqs. (6) and (7) are obtained where τ eq = 1/ √ C c C a is the equilibrium timescale. Once t > τ eq , dN χ /dt = 0, and N χ (t > τ eq ) = C c /C a according to Eq. (6). The total annihilation rate at this stage only depends on the capture rate since Γ a = C a N χ Nχ = C c . 1 Note that C c depends on σ χn and σ χn ≤ 1 The C a depends on the annihilation cross section σv explicitly but its effect only appears in the total annihilation rate Γ A through tanh(t/τ eq ) as τ 2 eq ∝ 1/ σv . No matter what the value of σv is, one always has tanh(t/τ eq ) ≈ 1 as long as the NS age is greater than the equilibrium time scale. In this case the total annihilation rate Γ A no longer depends on σv , but rather it is solely determined by σ geom χn ≈ 10 −44 cm 2 where σ geom χn is the geometric cross section. In principle, the maximum capture rate is determined by C c (σ geom χn ) = C geom c . Besides, when DM falls into the NS surface, it is accelerated up to 0.3c − 0.5c. The nonrelativistic (NR) limit for calculating σ χn is not applicable. Furthermore one has to consider contributions from axial-vector coupling due to V −Z mass mixing given by L mass . We have thoroughly included these effects. A brief discussion on how to compute σ χn in terms of relativistic kinematics is given in Appendix A.
NS is known to suffer from eternal cooling due to neutrino and photon emissions. Without extra energy injection, the NS temperature drops until it releases all its heat. However, if SM particles are produced due to DM annihilation in the star, these particles can become a heat source and potentially prevent the star from inevitable cooling. Therefore, the evolution of NS interior temperature T b is governed by the equation where ν ≈ 2.1 × 10 4 erg cm −3 s −1 (T b /10 7 K) 8 is the neutrino emissivity, γ ≈ 1.8 × 2 is the photon emissivity, χ is the DM emissivity that is responsible for the heating from DM annihilation, and c V the NS heat capacity [18]. Additionally, the surface temperature T s observed by a distant observer is related to T b by T s ≈ 8.7 × 10 5 K (g s /10 14 cm s −1 ) 1/4 (T b /10 8 K) 0.55 where g s = GM/R 2 ≈ 1.85 × 10 14 cm s −2 accounts for the redshift correction from the surface gravity of the star. It is also pointed out that when T b < 3700 K, there is no distinction between T b and T s [23].
During each annihilation, a pair of DMs release 2m χ of energy in a form of SM particles or dark bosons depending on which channels are kinematically allowed. The total energy release rate by DM is E χ = 2m χ Γ a i b i where b i is the branching ratio of a specific channel, e.g., e ± , µ ± , τ ± or qq, and i b i ≤ 1. Neutrino pair νν is also part of the annihilation channel in the presence of V −Z mass mixing, but it cannot contribute to the heating. In addition to the annihilation, DM also loses kinetic energy E k to the star through the capturing process.
This has been realized as the kinetic heating [21] with the rate the capture rate C c . We have carefully examined that even σv is an order of magnitude smaller than the thermal relic one given in Ref. [71], the NS can still attain equilibrium at t ∼ Gyr and have no sensitivity on σv anymore. The typical timescale for τ eq is about hundreds to thousands of years for ( σv , σ χn ) = (6 × 10 −26 cm 3 s −1 , 10 −47 cm 2 ) for m χ around MeV to TeV. where γ = 1/ √ 1 − v 2 is the Lorentz factor. 2 Thus, DM emissivity χ is given by where V is the NS volume.

III. DECAYS OF DARK BOSON
Here we discuss the case of V produced by DM annihilation. V is usually produced in DM rich environment. If V can subsequently scatter off the surrounding DM multiple times, it could lose energy and be self-trapped. It then decays promptly as shown in Fig. 1b.
However, such self-trapping effect is in general inefficient since the χV scattering length χV is much larger than the thermal radius r th . Hence the scattering rate is suppressed and irrelevant to the heating. We show detailed discussions in Appendix C. Another trapping is due to the scattering between V and neutrons. On the other hand the relevant cross section is further suppressed by the factorε 4 Z and the scattering length is expected to be much larger than the NS radius. It is safe to omit this effect in our calculation as well.
However, V can decay into other SM particles before it propagates to the surface as long as the decay length dec is shorter than R 0 (see Fig. 1a). The decay length is given by Since V is produced on shell, we do not consider V decaying back to χ due to m V < m χ . The probability for V to convert into SM particles after a propagation distance r is We took r = R 0 in the calculation. However, if neutrino is the decay product, it cannot be considered as the heating source and must be subtracted. By examining the numerical results for F , we found that V can decay before it exits the star in most of our parameter space of interest. This implies that χχ → 2V also plays an important role in NS heating (see Appendix C for details). Generally speaking, NS contain muons and electrons that are degenerate. To enable the decay V → ¯ with corresponding to either the electron or the muon, V should not only be heavier than 2m , it also has to be energetic enough so that the kinetic energy of exceeds the chemical potential of for preventing the Pauli blocking effect. This condition has been implemented in our study.
Given the information in this section, we summarize that even when χχ → 2V dominates the annihilation channel for m V < m χ , the heating effect is still efficient due to V decays.
However, the self-trapping is generally unimportant due to χV r th in this paper.

IV. IMPLICATION OF DM ON NS TEMPERATURE
In this section, we describe how NS surface temperature T s is affected by the DM annihilation. If χ is negligible, the standard cooling mechanism gives T s ≈ 100 K for a 3-Gyr-old NS. But when χ is large enough to counterbalance γ,ν , T s could remain at a relatively higher temperature. We present the numerical results of T s for both α χ = 1 and 0.01 in Figs. 2 and 3 respectively. The adjacent DM density around NS is assumed to be the same as that of it requires multiple scatterings to capture the DM and implies that the single-scattering capture is inefficient. Thus, we restrict our discussion below to the TeV DM where NS has better sensitivity and can be complementary to current DM direct searches.
From left to right, we have m V /m χ = 10 (heavy mediator), 1 (equal mass) and 0.1 (light mediator). T s is indicated by the color bar placed on the right and the lowest temperature is 100 K. Without annihilation, e.g. no anti-DM exists, solely kinetic heating can raise T s up to 1750 K. If DM annihilation is included, T s can maximally reach to 3100 K.
Various constraints are also plotted, including XENON1T [8], XENON LDM (low mass DM) based on the ionization [10], and of Migdal [9] effects, SIDM [48][49][50][51][52], SN1987A [78] and beam dump experiments [75][76][77]. The parameter curve rendering DM annihilation cross section at the thermal relic value [71] in the early Universe is plotted in green on each figure for comparison. 3 We adopted the method given in Ref. [80] for computing the Sommerfeld enhancement factor. The DM relative velocity in the early Universe is taken to be c/3. See Appendix B for details. Here we present the thermal relic cross section as a reference point and refer the readers to Refs. [53,58,81] for detailed discussions. In addition, although the captured DMs can have relatively large Sommerfeld enhancement due to low velocities, 4 the enhanced σv only shortens the equilibrium timescale τ eq . When t τ eq , the total 3 Since we have taken DM as Dirac fermions in the symmetric scenario [53], the thermal relic cross section is therefore two times larger than that in the Majorana DM case. Note that the thermal relic σv is in general m χ dependent and slightly deviates from the canonical value 6 × 10 26 cm 3 s −1 [71]. We have incorporated this for determining the green curve on each figure. 4 Assuming DMs are thermalized with the NS core where T χ = T b . Thus the mean velocity is about annihilation rate only depends on the capture rate with Γ A = C c . The NS is generally insensitive to the Sommerfeld enhancement as long as the DM is in equilibrium. when m V = m Z with m Z the SM Z boson mass. In fact the value forε Z at this point is Thus, the DM-neutron scattering cross section σ χn depends onε Z and is proportional to in the NR limit. (See Eq. (A12) for reference.) 5 The last term shows the suppression factor due to Pauli blocking where ξ ∼ q/µ F with q being the momentum transfer during the scattering and µ F the neutron chemical potential.
In the equilibrium epoch, t τ eq , the total annihilation rate Γ A = C c ∝ σ χn . 6 Whenε Z is at the resonant point, σ χn is enhanced drastically by the factor m 2 Z /Γ 2 Z so does the DM heating resulted from DM emissivity χ . This accounts for the dip at m V = m Z in each figure.
On the other hand, DM heating for m χ in the sub-GeV region is much stronger. It can be understood that, as m χ m n , q ∝ m χ while m V /m χ is held fixed, we have σ χn ∝ε 2 Z /m χ according to Eq. (12). Hence a smaller m χ leads to a larger σ χn as well as a more effective DM heating. However, the effect of DM heating will not grow indefinitely withε Z as σ χn ≤ Nonetheless, in the presence of V −Z mass mixing, neutrinos are also part of the annihilation products and take a significant branching ratio in the DM annihilation in such a mass region. Neutrino cannot contribute to the heating-this explains why T s is much colder when m χ O(170) MeV. The DM heating in this region is mainly due to kinetic heating.
As σ χn = σ geom χn , the resulted T s is around 1750 K from pure kinetic heating. Various η values in Fig. 2 characterize the contributions from ε γ,Z toε Z . 7 Both η = 1 and 0 are similar because even when ε γ = 0, its effect onε Z is suppressed by m 2 V /m 2 Z as 5 In the numerical calculation, we used the general expression for σ χn , Eq. (A7), and the derivation is given in the same appendix. Nonetheless, Eq. (12), or Eq. (A12), is simpler and suitable for our discussions in the main text. 6 We found that the equilibrium condition holds in most of the parameter space in this work. However, in the calculation we adopted Γ A = C a N χ Nχ with N χ (Nχ) given by Eq. (8), instead of simply assuming Γ A = C c . 7 Since neutron charge is neutral, Q = 0, the effect of kinetic mixing Qε γ in Eq. (A4a) has zero contribution seen from Eq. (A5). For η = 1, the kinetic mixing can contribute comparably to the V − Z mass mixing unless m V > m Z . This can be clearly seen in Fig. 2 that the difference between η = 1 and 0 is apparent only in m V > m Z region, which is the region to the right of the dip.
To the left of the dip, the contribution from ε γ toε Z for η = 1 is negligible.
For η = ∞, ε Z vanishes so that the only contribution toε Z comes from ε γ . As discussed earlier, the effect of kinetic mixing term is suppressed by m 2 V /m 2 Z and thus σ χn ∝ε 2 Z ∝ ε 2 γ m 4 V /m 4 Z . The associated DM heating is, in general, much weaker than the cases with η = 1 and 0. However, the advantage of η = ∞ is that no neutrinos can be produced by the DM annihilation due to the absence of V − Z mass mixing. The energy released from DM annihilation can be fully deposited into the NS. This accounts for the higher T s than when η = 1 and 0 in terms of the same σ χn , but the difference is not apparent. Numerical calculation shows it is around ten to O(100) K.

B. Case for m V < m χ
For the light mediator case, the channel χχ → 2V dominates over χχ → ff due to spectral flux density with T s at a given frequency ν is given by [21] f where R 0 is the NS radius, γ ≈ 1.35 is the relativistic factor of DM on the NS surface, d is the distance between the NS and the Earth, and a = ω/k B T s . It is easy to estimate that f ν peaks at a ≈ 3 and the peak value is clearly proportional to T 3 s . Taking T s = 2000 K and d = 10 pc as an example, f ν peaks at ν −1 = 2 µm with the peak value of 0.84 nJy. The SNR for JWST-like telescope scales as √ t exp for a given ν where t exp is the exposure time. This scaling stems from the fact that, for the exposure time t exp , the signal photon number in the frequency range [ν, ν + dν] is proportional to f ν · dν × t exp . On the other hand, the statistical fluctuation of the photon number, arising mainly from in-field and scattered zodiacal light, scattered thermal emission from the telescope, and the scattered starlight, scales as √ t exp for a given frequency ν. This results into the above mentioned scaling for SNR. It is important to note that SNR is a complicated function of ν because the above mentioned backgrounds are also wavelength dependent [66].
To reach the peak of f ν for T s = 2000 K, i.e., 0.84 nJy with SNR = 10, one requires 1.4 × 10 6 s of exposure time. For SNR = 2, which is the criterion for our presentation below, the required exposure time is 5.6 × 10 4 s.
In Fig. 4, we plot the t exp for obtaining SNR = 2 over d − T s plane. The region enclosed by the red line represents t exp < 10 5 s. There are multiple filters available for NIRI with ν −1 centered at various different values [67]. We select the filter with ν −1 most suitably matching the corresponding blackbody wavelength at T s . For instance, F200W filter is used for T s close to 2000 K, while F277W filter is adopted for T s around 1500 K. The switching of filters when appropriate is reflected in the zigzag behavior of the red sensitivity curve in Fig. 4. In principle, as σ χn ∼ σ geom χn , kinetic heating can maximally warm the NS up to 1750 K without DM annihilation. For NS that is located within 10 pc, JWST can achieve SNR = 2 with t exp ≤ 10 5 s for T s ≥ 1750 K

VI. SUMMARY AND OUTLOOK
In this work we have investigated the new dynamics arising from the kinetic mixing and V − Z mass mixing between the dark gauge boson V of the broken U (1) X symmetry and neutral gauge bosons in SM. In particular, V − Z mass mixing induces a resonance at m V ≈ m Z , which can be seen from the pole ofε Z at m V = m Z . The axial-vector part of the coupling between V and SM fermions has been included in our calculations. As χχ → 2V dominates the annihilation channel for m V < m χ , V can decay into a pair of SM fermions before it exits NS and induces NS heating in addition to χχ → ff . Although this contribution appears naturally in the dark boson model considered here, it is usually not included in the model-independent analysis, such as the one performed in Ref. [23]. We also demonstrated numerically that NS can provide constraints on sub-GeV DM with feeble coupling to SM particles complementary to the current direct search. The detectability with reasonable t exp in JWST telescopes is discussed. Similar conclusion can be drawn for the future JWST-like telescopes.
We note that this work only considers χn scattering in the capture rate. This explains why NS is not sensitive to the dark sector when ε Z = 0 (η = ∞). Neutrons interact with DM only through NC interaction governed byε Z . Once ε Z = 0, NC interaction becomes suppressed since ε γ inε Z is oppressed by m 2 V /m 2 Z . However, NS also consists of protons, although the fraction of them is rather small. When protons are included, charged current interaction will be involved for the capture of DM and NS remains sensitive to the dark sector even for ε Z = 0. In general, NS sensitivity will be improved by including proton contributions. We leave this for future studies. gives rise to relations between fields in the gauge basis and those in mass eigenstate basis.
However, since we are only interested in interactions up to O(ε γ ) or O(ε Z ), we do not need to perform the diagonalization but rather treating the mixing terms ε γ B µν V µν /(2 cos θ W ) and ε Z m 2 Z Z µ V µ as perturbations. These two mixing terms generate the following two-point functions at the tree level where k is the four-momentum of V entering into kinetic mixing or V − Z mixing vertex. Hence the electromagnetic coupling of V to SM fermions results from multiplying the two-point function iΠ µν V γ , the photon propagator iD γ αµ (k), and the electromagnetic coupling ieA α J α EM , as shown in Fig. 6. This multiplication leads to Similarly, NC coupling of V to SM fermions is given by multiplying the two-point function iΠ µν V Z , the Z boson propagator iD Z αµ (k), and the NC coupling igZ α J α NC / cos θ W . This gives rise to . (A3) Here we have used the physical conditions k 2 = m 2 V and k µ µ V = 0. We have also chosen unitary gauge for the Z-boson propagator. Therefore, the interaction vertex between dark bosons and neutrons in Fig. 5 have the following Lorentz structure ieψ n γ µ (a f +b f γ 5 )ψ n with and Γ Z is the Z boson decay width, Q and I 3 are the electric charge and the weak isospin respectively. In Table I, we list Q and I 3 for various particles. The values for neutrons can be obtained by summing the corresponding quantum numbers of three quarks udd in the low energy limit.
Mixing parameters ε γ andε Z are responsible for electromagnetic and NC interactions, respectively. Electromagnetic interaction does not contribute to σ χn since Q = 0 for neutrons.  On the other handε Z has a feeble dependence on ε γ with a suppression factor m 2 V /m 2 Z when m V m Z . This explains why σ χn is still nonzero when ε Z = 0 (η = ∞).
The spin-averaged χn scattering amplitude is given by where s, t, and u are the Mandelstam variables. DM scatters with neutrons when its velocity boosted to 0.3c − 0.6c by the NS gravity. It must be treated relativistically. However, neutrons can be treated as at rest since the chemical potential is O(200) MeV in the star.
Therefore, from the method in Ref. [82], we are able to write down the DM-neutron scattering cross section as is the Källén function, and where m 1 = m 3 = m χ and m 2 = m 4 = m n for the χn scattering. In Eq. (A10), the energy E 1 = γm 1 is the total energy carried by particle one, which is DM.

Pauli blocking in the χn scattering
Note that if the momentum transfer √ −t in Eq. (A7) is smaller than the Fermi momentum, the suppression by Pauli blocking takes effect. We include this in the numerical  [69]. Our result agrees with Ref. [69] in the three benchmark scenarios that |M χn | 2 are constant, t-dependent and t 2 -dependent.

Axial-vector contribution in the NR limit
If χ can be treated nonrelativistically as well, we have s = m 2 χ + m 2 n + 2m χ m n , u = m 2 χ + m 2 n − 2m χ m n and t = 0. Therefore the amplitude and the cross section become, which are independent of b f where it determines the strength of axial-vector coupling.

Appendix B: DM annihilation
We can divide the DM annihilation into two categories, which are m V ≥ m χ and m V < m χ , respectively. For the prior case, DM can only annihilate into SM particles as shown in Fig. 7a. For the later one, as long as g d eε γ,Z , the dominant annihilation products are two dark bosons V as shown in Fig. 7b. The amplitude for χχ → ff is given by where Γ V is the V decay width. Assuming DM is at rest in the star, the amplitude can be simplified into The partial decay widths of V are given by for V → ff and for V → χχ. Note that we have omitted the Heaviside theta function θ(m V − 2m χ,f ) in the above expressions but it is always implemented when we perform the calculation to ensure the energy conservation. Besides, when m χ > m V , the channel χχ → 2V is allowed and the amplitude is In the NR limit, Thus, the general expression for annihilation cross section is obtained by using the Fermi golden rule, where m f is the final state particle mass and µ f F the chemical potential of fermion f in the star. There is no chemical potential for dark boson V . Therefore, we arrive at for χχ → ff and for χχ → 2V . The total annihilation cross section is the sum of both We note that the second term contributes when m V < m χ .

Appendix C: Dark boson in the star
Dark bosons can be produced from DM annihilation once m V < m χ . This channel is thought to have feeble effect on the heating since V interacts weakly with the NS medium and escapes without any trace. However, we found that, depending on the strength of ε γ,Z , V can decay into SM particles before it reaches the surface of the star. In the case that the decay length dec is much smaller than the radius of the star, the total energy released from the annihilation can be fully deposited to the star. See Fig. 1a. We also examine the case where V is produced in the DM rich region in the center of the star. V could undergo multiple scattering with the surrounding DM and self-trapped until it decays. See Fig. 1b.
This is another way to extract energy from V . We discuss both effects in the following.

Decay length
The dark boson decay length with time dilation effect is given by where v = 1 − m 2 V /m 2 χ is the V velocity and τ dec = Γ −1 V the V lifetime at rest. Let us assume that V is produced in the center of the star and its propagation distance is R 0 . Fig. 8 presents F defined in Eq. (11), i.e., the fraction of V converting into SM particles after traveling a distance r = R 0 , as functions of ε Z,γ and m χ for m V = 0.1m χ . We have subtracted neutrino contributions from F since they cannot generate heat. Since the branching ratio of V decays to neutrinos is nonzero in the case of V − Z mass mixing, F is generally smaller than 1 for η = 1. For η = ∞, no neutrinos can be produced, thus F can reach unity. In these figures, the chemical potential for electron µ e F is about O(170) MeV. For a dark boson at rest with m V ≤ µ e F , V → e + e − can be Pauli blocked even for m V ≥ 2m e . On the other hand, if V is highly boosted as a result of heavy DM annihilation, V → e + e − is not Pauli blocked as long as m χ ≥ µ e F . Therefore, to enable V → ff decays, two conditions are required. The first is m χ ≥ µ f F and the second is m V ≥ 2m f .

Dark boson-DM interaction length
Feynman diagrams contributing to χV scattering are shown in Fig. 9 and the amplitude is given by − m 4 χ (3s 2 + 14st + 3t 2 ) + m 2 χ (s 3 + 7s 2 t + 7st 2 + t 3 ) To compute the scattering cross section σ χV , it is fair to assume DM at rest. However, V is produced with relativistic velocity since m χ > m V . We follow the procedure given in Eqs. (A7)-(A10) and set m 1 = m 3 = m V and m 2 = m 4 = m χ . Thus, Note that χV scattering is not subject to Pauli blocking since DM does not become degenerate in the presence of annihilation.
The χV scattering length χV is given by with n χ ≡ N χ /V χ the average DM number density. The volume characterizing DM in NS is is the thermal radius. If χV r th , V can scatter with surrounding DM multiple times and gradually lose its kinetic energy. However, our numerical result shows that χV r th in all of our interested parameter space. In Fig. 10, we take α χ = 1 and T χ = 1000 K. The choice α χ < 1 makes χV even longer due to a weaker χV interaction. For η = 1, the region for χV /r th < 1 happens when m χ 300 MeV. However, even V can be self-trapped, it hardly decays into particles other than neutrinos because the allowed channels, eg. e ± and µ ± , are Pauli blocked. For η = ∞, only a very small parameter space leads to χV /r th < 1.
Therefore, we conclude that the self-trapping of V is insignificant, hence only V decays contribute to the energy injection.