Probing dark gauge boson with observations from neutron stars

We present an investigation on the production of light dark gauge bosons by the nucleon bremsstrahlung processes in the core of neutron stars. The dark vector is assumed to be a $U(1)_{B-L}$ gauge boson with a mass much below keV. We calculate the emission rate of the dark vector produced by the nucleon bremsstrahlung in the degenerate nuclear matter. In addition, we take into account the photon-dark vector conversion for the photon luminosity observed at infinity. Combining with the observation of J1856 surface luminosity, we find that a recently discovered excess of J1856 hard x-ray emission in the 2-8 keV energy range by XMM-Newton and Chandra x-ray telescopes could be consistently explained by a dark vector with gauge coupling $e^{\prime}=5.56\times 10^{-15}$, mixing angle $\varepsilon=1.29\times 10^{-9}$, and mass $m_{\gamma^{\prime}}\lesssim 10^{-5}$ eV. We also show that the mixing angle $\varepsilon>7.97 \times 10^{-9}$ for $m_{\gamma^{\prime}} \lesssim 3 \times 10^{-5}~\rm eV$ and the gauge coupling $e^{\prime}>4.13 \times 10^{-13}$ for $m_{\gamma^{\prime}} \lesssim 1~\rm keV$ have been excluded at 95% confidence level by the J1856 surface luminosity observation. Our best-fit dark vector model satisfies the current limits on hard x-ray intensities from the Swift and INTEGRAL hard x-ray surveys. Future hard x-ray experiments such as NuSTAR may give a further test on our model.


I. INTRODUCTION
One of the significant puzzles in particle physics that demand new theory beyond the Standard Model (SM) is the nature of dark matter (DM), whose existence has been confirmed from various cosmological and astrophysical observations [1]. The most popular class of DM candidates is the weakly interacting massive particles (WIMPs) [2,3] with a mass at around the electroweak scale and a coupling strength similar to the weak coupling. However, the WIMP hypothesis has suffered stringent constraints due to the null results from both direct [4] and indirect [5] DM detection experiments. As an alternative candidate of DM, a hidden sector consisting of very low-mass particles has drawn more attention in recent years.
One of the simplest extensions of the SM is to include an additional U (1) gauge group [6,7], which can naturally arise from the breaking of grand unified theory (GUT) groups [8,9] or the string theory [9][10][11][12][13][14][15][16]. The gauge boson associated with the new U (1) group, dubbed the dark photon (which mixes with the SM photon), may play the role of a mediator between the SM and dark sectors or as a dark matter candidate [17].
Supernovae serve as a powerful factory for the emission of dark photons with masses 100 MeV [18][19][20][21][22][23][24][25]. With a sufficiently weak interaction, dark photons generated within the core of a supernova can escape the progenitor star before their decays, and will contribute to the stellar energy transport. The supernova cooling would be accelerated if the dark photons could transport energy out of the core more efficiently than the standard supernovae cooling via the emission of neutrinos. The bounds on dark photons by using the measurements of SN 1987a have recently been updated in Refs. [22][23][24], with the consideration of plasma effects for the production of dark photons. The Sun, horizontal branch (HB) stars, and red giants, on the other hand, can serve as important laboratories for dark photons in the mass range of keV. For instance, Refs. [26,27] show that the measured luminosity of the Sun, L = 3.83 × 10 26 W, has constrained the dark photon parameter space to the range of εm γ < 4 × 10 −12 eV [27].
In this work, we will focus on the production of dark gauge bosons in the core of neutron stars (NSs), which are the remnants of the supernova explosion and have long been recognized as excellent laboratories to search for axionlike particles [28][29][30][31][32][33][34][35][36][37][38][39][40][41]. The core of a NS is composed of highly degenerate nuclear matter that has an average density of ∼ 3 × 10 14 g/cm 3 , which corresponds to an average distance of ∼ 1 fm between nucleons. Dark gauge bosons can be produced through bremsstrahlung of proton and neutron scattering in the NS cores.
In the previous literature [22][23][24], the nondegenerate limit has been assumed to obtain the dark vector emission rate as well as the constraints on the dark vector from the supernova, where the plasma is thought to be diluted and nonrelativistic. However, the nuclear matter in the NS core is in fact highly compressed, and therefore, in this work, we will present the calculation of the production of the dark vector in the NS core in the degenerate limit.
The in-medium effects can lead to a suppression on the dark vector production rate if the dark vector which couples coherently to the charged SM plasma has a mass much below the stellar plasma frequency [22,23]. In additional to the electromagnetic (EM) current through mixing, the new gauge boson could couple to the B − L current in the well-known theory of U (1) B−L extension of the SM [42][43][44][45][46][47]. For such a theory, the new vector boson will dominantly be produced by the neutron bremsstrahlung in the NS core. This is because, for one thing, neutron is the dominant component of the NS, and for another, the vector boson emitted by the neutron in the initial or final state is not suppressed by the plasma effects. In this work, we will assume that the dark gauge boson is the B − L gauge boson with a mass much below keV, the core temperature of the NS, and the gauge coupling e is sufficiently small.
Recently, a significant excess of hard x-ray emission in the 2 − 8 keV energy range was found by analyzing the data from the XMM-Newton and Chandra x-ray telescopes [48] observing on the nearby magnificent seven (M7) x-ray dim isolated NSs. The excess was interpreted as the emission of axions by the NS in Ref. [38] and was found to be consistent with the current constraints. In this work, we explore the interpretation of a dark vector scenario for the J1856 hard x-ray excess. We calculate the dark vector emission rate in the NS core and simulate the evolution of the NS based on the modified NSCool code [50], including the dark vector emissivity. Due to the strong magnetic field of the magnetospheres surrounding the NSs, the dark vectors produced in the NS cores may convert into x-rays when they propagate outwards [49]. We present an analytical formula of the dark vectorphoton conversion probability for dark vector with mass m γ (T s /R s ) 1/2 , where T s and R s are the surface temperature and radius of the NS, respectively. We also calculate the surface luminosity observed at infinity by taking into account the photon-dark vector conversion.
Based on the Bayesian inference, the analysis of J1856 hard x-ray data, as well as surface luminosity observation are carried out by the UltraNest package [51]. We find that the dark vector scenario is favored by the hard x-ray excess, and we also derive 95% upper limits on the model parameters from the data.
The structure of this paper is arranged as follows. Section II defines the dark gauge boson model considered in this work, and sets up the required notations. In Sec. III, we calculate the production of dark vectors from nucleon bremsstrahlung, including the squared matrix element and emission rate. In Sec. IV, we briefly review the neutrino and photon luminosities for the NS cooling. We also calculate the photon surface luminosity observed at infinity that includes where we assume the dark vector has a Stueckelberg mass m γ , are, respectively, the field strengths of the SM photon and dark gauge boson, and ε represents the mixing angle between the SM photon and dark vector. The parameter e denotes the coupling between A D µ and the current J µ , which in this work is assumed to be the SM B−L current arising from a U (1) B−L symmetry, as widely investigated in the literature [42][43][44][45][46][47]. The kinetic terms of gauge boson in the Eq. (1) can be diagonalized by rotating the gauge fields as where A D µ , A SM µ T and A µ , A µ T are the so-called interaction eigenstates and mass (propagating) eigenstates, respectively. For a massive dark vector, the rotation angle θ is locked at zero [17]. We see that the mixing leads to the interactions between the dark vector and the SM charged particles, such as electron, proton and charged pions, with a coupling strength εe. Furthermore, due to the mixing, the dark vector can convert to a photon during their propagation.
In the plasma of supernovae, the visible photon develops a nonzero mass that is determined by the plasma frequency ω 2 p = 4πα EM i n i /E F,i , where the sum goes over the charged particles with a number density n i and a Fermi energy plasma mass of the photon can lead to inequivalent effective couplings to the charged particles in the plasma and in the vacuum (the Lagrangian parameter). The effective coupling between the SM fermions and A µ in the medium is given by [25] e f eff = e q e q f − q f q e + (εe − q e e ) q f where f denotes the SM fermions with an electric charge q f and a dark gauge quantum number q f in units of e and e , respectively. The proton and neutron have the electric charges q p = 1 and q n = 0, and the dark gauge quantum numbers q p = q n = 1. The EM and dark gauge quantum numbers for the electron are q e = q e = −1. The visible photon mass is encoded in the polarization tensor π T,L given in Appendix A.
From Eq. (3), the effective couplings of electron, proton, and neutron to A µ in a dense medium (i.e., with π T,L T 2 , m 2 γ ) are explicitly given by We observe that the dark vector couplings to the electrically charged fermions are suppressed by both the mixing term in the Lagrangian and medium-induced effects, while the interaction of A µ to the neutrons is not suppressed by the plasma effects. As a result, there exists a suppression from the plasma effect in the emission of dark vectors from the electron and the proton currents. On the other hand, the production of dark vectors via the neutron current does not have such a suppression, and the emission rate from neutrons in the medium can be approximated by that in the vacuum.
For a NS of interest here, because of the large number of neutrons in a NS core, the dark vectors can be copiously produced via the bremsstrahlung process involving the neutron current. Furthermore, since we are concerned with the case of low dark vector masses, the contributions from the electron and proton bremsstrahlung processes in the crust can and will be neglected. Note that the nucleon superfluidity in the NS core is still under debate. At temperatures below the critical temperature of Cooper pair formation, the nucleon bremsstrahlung rate may be suppressed by nucleon superfluidity. However, recent studies [52] on NS cooling indicate that neutron superfluidity in the middle-aged NS core is weak and the critical temperatures are possibly too low to be relevant for our analyses.
Thus, following Ref. [38] the nucleon superfluidity effect will be ignored in our work.
There exists oscillations between dark vector and photon due to the mixing term in our model. In the presence of an inhomogeneous external field, the approximation for the dark vector-photon conversion probability in the weak-mixing limit is given by The integral starts from the stellar surface r 0 = R s since the x-ray photons produced inside the NS would be absorbed. The parameters and detailed calculations for Eq. (5) are presented in Appendix B. For m γ (T s /R s ) 1/2 , we find an analytical formula of Eq. (5), which can be found in Appendix C. We will see that the NS is an excellent test bed for the dark vector hypothesis since the conversion probability can be largely enhanced under a strong magnetic field. For the approximation to make sense, m γ cannot take the vanishing limit and the numerical value for ε should be sufficiently small to satisfy the weak-mixing condition.
A few remarks are in order. The traditional U (1) B−L model is gauge-anomaly free if right-handed neutrinos are introduced. While anomaly free at higher energies, the effective Lagrangian (1) is anomalous due to the mixing term at low energies. And it is the phenomenology of this low-energy effective Lagrangian (1) that we want to discuss in this work.
We also note that the constraints on e that we obtain from NS cooling in this work do not depend on the mixing term and can be applied to the traditional U (1) B−L model.

III. DARK VECTOR PRODUCTION
In this section, we calculate the emission rate of dark vectors from the neutron bremsstrahlung in the NS core. Supplements for the calculations are given in Appendix D.

A. Dark vector emission rate
The dark vector energy emission rate per unit volume for the process N 1 + N 2 → N 3 + N 4 + γ (N = n, p) in a strongly degenerate nuclear matter is given by [30,37] where the symmetry factor S = 1/4 for the identical particles in the initial and final states and S = 1 for mixed processes, the squared matrix element |M| 2 = spins |M | 2 sums over initial and final spins, ω is the energy of emitted dark vector, p k is the momentum associated with the nucleon N k (k = 1, 2, 3, 4), and p γ is the momentum of the dark vector. In the presence of a nuclear mean field U N , the energy dispersion relation of the nucleon is modified to be where m * is the effective nuclear mass and E * = p 2 + m 2 * is the effective energy. Note that the nuclear energies that enter the energy δ-function and the denominator on the right-hand side of Eq.(6) are the effective energies [37]. The fermion phase-space distribution function is given by where µ k is the nuclear chemical potential of N k . With Eq.
The occupation numbers f 1,2 and the Pauli blocking factors (1 − f 3,4 ) are for the incoming and outgoing nucleons. The dark vector are assumed to escape freely so that a Bose stimulation factor as well as its absorption are neglected.
In the nonrelativistic (NR) limit, the nucleon effective mass m * is much larger than all other energy scales such as the temperature or Fermi energies, so that the energy taken by the nucleon is E * m * + E kin , where E kin = p 2 /2m * is the nuclear kinetic energy. In a bremsstrahlung process, the radiation typically carries the momentum |p γ | E kin |p| and, therefore, the outgoing nucleons carry essentially all of the momentum and the radiation momentum can be neglected. With these approximations, we show in the next section the scattering matrix element of the nucleon bremsstrahlung, and in Appendix D we calculate the emission rate (6) in the degenerate limit.

B. Squared matrix element
As shown above, the effective couplings of A µ to electron and proton in the medium are suppressed in the limit of light m γ . The emission of dark vector from the neutronneutron bremsstrahlung is the primary production mechanism in the inner core of a NS.
The nucleon-nucleon interaction in the nucleon-nucleon bremsstrahlung was treated by the one-pion-exchange (OPE) approximation in Ref. [53] and was used in the calculation of the axion emission rate from the nucleon-nucleon bremsstrahlung [30]. For our present work, it is sufficiently reliable to calculate the dark vector emission rate using the OPE approximation (see also Ref. [19]). Conservatively, we also include a factor of 1/4 [54] in the emission rate to account for the overestimation in the OPE approximation at lower temperatures [55].
Additional remarks on the OPE approximation are provided in Appendix E.
We show the Feynman diagrams for the process n + n → n + n + γ (where γ ≡ A ) in where m π is the pion mass, the momenta k = P 2 − P 4 , l = P 2 − P 3 , with P k being the four-momenta associated with n k (k = 2, 3, 4), and k = p 2 − p 4 and l = p 2 − p 3 . In the NR limit, the nucleon mass is much larger than the temperature, the dark vector mass, and the transferred momenta, i.e., m 2 * k 2 , l 2 . Furthermore, the four-momenta transferred among the nucleons are much larger than the momentum carried by the dark vector, i.e., k 2 , l 2 p 2 γ , ω 2 , m 2 γ . For the model considered in this work, the coefficients are given by where f nn = f and f 1.05 is the pion-neutron coupling, the parameter g α = g β = e ≡ g n is the coupling between A and neutron. We see that only the coefficient C k is nonzero. The Feynman diagrams for the process n + p → n + p + γ are depicted in Fig. 2. The squared matrix element for this process is given by where C k = C l = f 4 np g 2 n , with f np = √ 2f as required by isospin invariance. Due to the plasma effects, we have neglected the diagrams where the dark vector is emitted from an electron or proton.

C. Dark vector luminosity
Given the emission rate, to compute the luminosity of dark vectors from the NS core we need to know the NS equation of state (EOS), which determines the matter distribution in the NS, the chemical potential of the matter, and the profiles of neutron and proton Fermi momenta. In our analysis, we employ the Akmal-Pandharipande-Ravenhall (APR) EOS [56] to model the uniform nuclear matter in the NS core and assume a NS mass of 1.4 M , which corresponds to a NS core radius of 11.0 km. It is expected that the nucleon bremsstrahlung in the NS core dominates the production of dark vectors. For the weak coupling theory, the mean free path of dark vector in the medium is much larger than the size of the NS. So in the calculation we have neglected the dark vector reabsorption effect.
The differential dark vector luminosity is given by integrating the differential emissivity over the volume of the NS dL γ dω = 4π Integrating the differential luminosity over ω, we obtain the total dark vector luminosity.
With the APR EOS for the NS and a light dark vector with mass m γ keV, the dark vector luminosity from the NS core can be estimated as where T c is the core temperature of the NS. The emission of dark vectors can result in NS cooling. Since the standard NS cooling scenario in terms of neutrino and photon emissions fits rather well to the observation of the NSs, the stellar cooling can be used to constrain the dark vector model.
In Fig. 3, we plot the dark vector "energy emission rate" Q γ /I, where I is given by   As shown in the left plot of Fig. 4, the differential luminosity of dark vectors is nearly a constant for ω T c and cuts off at ω ∼ T c due to the factor e −ω/Tc in the differential luminosity. In the right plot of the figure, we observe that both the constant regime and the cutoff point of dL γ /dω increase with the NS core temperature.

IV. NEUTRINO AND PHOTON LUMINOSITIES
One of the most fascinating stars known in the Universe is the NS, whose birth in supernova explosions is accompanied by the most powerful neutrino outburst. The standard NS cooling scenario based upon the pioneering works [53,[57][58][59][60] includes several neutrino emission processes (see [61] for a review). The direct Urca cooling process consists of beta decay and electron capture that can induce a huge sink of energy in the NS core. However, a sufficiently small separation between the Fermi levels of protons and neutrons is required to activate this mode, which generally requires a minimum mass of the NS, denoted by M D , that is significantly larger than the canonical NS mass of 1.4 M [62,63]. For the NS with mass less than M D , the modified Urca cooling mode can take place everywhere and dominate the NS cooling process. The modified Urca process is similar to the direct Urca, but involves an additional nucleon spectator, The modified Urca by the neutron branch is shown by Eq. (14). Its neutrino emission rate in the NS core is given by [53] Q n where ρ is the NS mass density profile, ρ 0 = 2.8 × 10 14 g · cm −3 is the nuclear saturation density, T c is the NS core temperature, and the suppression factor R M ≤ 1 appears with the onset of superfluidity. Cooper pair cooling could become dominant if the superfluidity occurs. However, as argued above, the critical temperature for the Cooper pair formation is possibly too low to be relevant for this work. We thus do not take the superfluidity into account.
The proton branch, Eq. (15), was first analyzed in Ref. [64] by using the same formalism as that in Ref [53]. The expression for the neutrino emissivity in the proton branch can be approximated by the following rescaling relation [61] Q p where m p * and m n * are the effective masses of proton and neutron, respectively, and p F,e , p F,p , and p F,n are the Fermi momenta of electron, proton, and neutron, respectively. Note that Θ F = 1 for p F,n < 3p F,p + p F,e ; otherwise, Θ F = 0. As an example, take p F,e = p F,p = p F,n /2, we have Q p ν 0.5Q n ν , i.e., the proton branch is nearly as efficient as the neutron branch. In addition to the Urca processes, the neutrino's emission processes can also be induced by the bremsstrahlung in nucleon-nucleon collisions, which, however, are subdominant for the NS cooling (a summary of these processes can be found in Figs 11 and 12 of Ref. [61]). Summing up all of the processes, the total neutrino emission rate can be estimated as Q ν RQ n ν , with a rescaling factor R 1.5 [61]. With the neutrino emission rate, we can then determine the neutrino luminosity by the integration where r 0 = 11 km is the NS core radius. For the APR EOS, the neutrino luminosity is then given by Previous literature has obtained the constraints on new physics models, e.g., the axion coupling, by requiring that the emissivity of the new particle should not exceed that of the neutrino, i.e., L γ < L ν [65]. Otherwise, the evolution of the NS will be strongly altered by the new physics, which is in conflict with the observations. As shown above, in addition to the NS density profile, the luminosities of dark vector and neutrino strongly depend on the core temperature, which, however cannot be determined by the EOS or by direct measurements. The core temperature may be inferred from the NS surface temperature observations but suffers from large uncertainties [38]. In Fig. 5, we depict the relation between core temperature and surface temperature at infinity using NSCool code [50] (with the default NS model). We observe that the relation is independent of the new gauge coupling e (with ε = 0). However, as we will show below, the luminosity at infinity cannot be solely determined by the luminosity at NS surface if there are conversions between photon and dark vector during their propagation toward the observer. In this case, there is no direct relation between core temperature and surface temperature at infinity. Furthermore, in the standard NS evolution, the additional heating effects can be neglected for the middle-aged NS, and the evolution of the NS core temperature with the time is expected to follow a power law T c (10 9 K) (t/yr) −1/6 [63]. Such a relation may also break down in new physics if the cooling of the NS is strongly affected by the emissions of new weakly-interacting particles.
In this work, we will perform numerical simulations of the NS cooling based on the modified NSCool code that includes additional energy loss via the dark vector emission. In this way, we determine the core temperature and the surface luminosity for the NS with a given age.
This will be described in detail below.
In the standard NS cooling model, the emissions of neutrinos and photons lead to the NS cooling. In the early stage, the neutrino emission is the dominant mode for the NS cooling.
As the stellar temperature drops, the NS cooling by photon emissions becomes significant since the neutrino luminosity decreases much faster than that of photon. For NSs with age t 100 kyr, the photon emission will exceed the neutrino one [65], and thus dominates the NS cooling process.
We now turn to the photon emission of the NS. The NS cooling due to the emission of photons is directly measured by the NS surface photon luminosity, which is given by where σ is the Stefan-Boltzmann constant, R s is stellar radius, and T s is stellar surface tem-perature. Without photon-dark vector conversions, the observed surface photon luminosity at infinity is accordingly redshifted as where φ s is the gravitational potential at the stellar surface, and where G is the gravitational constant, and M is the stellar mass.
However, the observed surface photon luminosity will be changed if the photons and dark vectors can convert into each other during their propagation to the Earth. For the NS with age ∼ 10 6 yr that is concerned in this work, we will show that the surface luminosity of the NS will be dominated by the photon emissions. For the parameter space we are interested in, the dark vectors' luminosity is always a subdominant component and, therefore, their contributions to the observed photon luminosity will be neglected. In this case, the surface photon luminosity observed at infinity is given by where P γ→γ (ω) = P γ →γ (ω) is the photon-dark vector conversion probability, which will be given in Appendix B. Since the Stefan-Boltzmann law (20) can be obtained by integrating the Planck's law over ω, the differential surface luminosity of the photon is given by We will show in Appendix C that the γ − γ conversion probability can be written as P γ→γ (ω) = P 0 ω 2 for dark vector in the mass range m γ (ω/R s ) 1/2 ∼ 3 × 10 −5 eV, with ω ∼ T s ∼ 50 eV. With these, we finally obtain the observed surface photon luminosity at infinity by taking into account the γ − γ conversion where A s = 4πR 2 s . The surface temperature at infinity is then given by Using these results, we can determine the observed surface luminosity and temperature by including the redshift by the gravitational potential as well as the γ − γ conversion once we know the values of L γ and T s at the NS surface.
The validation of Eqs. (25) and (26) should be further clarified. One may be concerned that the approximation of the probability requires ω m 2 γ R s . However, we have integrated over ω in the range of (0, ∞) for P γ→γ dL γ /dω. It is easy for the reader to confirm that there is nearly no difference between the integration over ω in the ranges of (0, ∞) and (0.9T s , ∞), which means that the small values of ω 0.9T s make negligible contributions to the integration. We thus conclude that Eqs. (25) and (26) are valid for our calculations provided that m γ 3 × 10 −5 , with ω ∼ T s ∼ 50 eV.
In Fig. 6, we show the variation of the surface luminosity observed at infinity L ∞ γ with the mixing angle ε. We assume the surface luminosity L γ = 8.4 × 10 31 erg/s and temperature T s = 47.1 eV for the NS. As shown in this figure, the γ − γ conversion strongly reduces the observed surface luminosity when the mixing angle reaches values of ∼ 10 −8 .

V. NS COOLING SIMULATION
In this work, we employ the NSCool code [50] for the numerical simulation of the thermal evolution of the star. The current version of the NSCool code incorporates all the corresponding neutrino cooling reactions, including direct and modified Urca processes, nucleonnucleon bremsstrahlung, as well as the thermal Cooper pair breaking and formation (PBF).
We modify the code to include the emissivity of the dark vector so that it can take part in the NS cooling along with the neutrino emission processes. For the core of the NS, we adopt the APR EOS with a mass 1.4 M (Prof_APR_Cat_1.4.dat). For the crust EOS, we use the default profile implemented by the file Crust_EOS_Cat_HZD-NV.dat.
We focus on the observations from one of the M7 members, RX J1856.5-3754 (J1856 for short), which locates at a distance 123 ± 13 pc away from us and has an age around (4.2±0.8)×10 5 yr [66,67] (a summary can be found in Refs. [52,68]). Recent studies [52] on NS cooling indicate that the critical temperature is very low and the neutron superfluidity is very weak for the middle-aged NS. We thus turn off the PBF processes in the simulations with NSCool, as is done in Ref. [38]. We find that the results from the standard NS cooling without including the PBF processes can account for the surface luminosity observations on J1856 quite well. However, the simulation including the PBF processes predicts a much lower surface luminosity than the observation.
In Fig. 7, we depict the neutrino (blue), dark vector (dark), as well as photon (yellow) surface luminosities as a function of the NS age. In the simulations, we assume  where L ∞ s is the observed surface luminosity of J1856 at infinity, and e 2φs is the redshift factor.

VI. X-RAY FROM DARK VECTOR-PHOTON CONVERSION
Recent analyses of the data from the XMM-Newton and Chandra x-ray telescopes [48] show a significant excess of hard x-ray emissions, in the energy range of 2 − 8 keV, from the nearby M7 x-ray dim isolated NSs. In particular, it has been shown in Ref. [48] that the NS J1856 hard x-ray spectrum has a ∼ 5σ excess, which is the most significant hard x-ray excess among the M7 members. The fact that the hard x-ray excess is observed in some NSs but not others depends on the experimental measurements as well as the NS properties [38].
Observations by the ROSAT All Sky Survey [69] have shown that all the M7 members have soft spectra that are well described by near-thermal distributions with surface temperatures ∼ 50 − 100 eV and, therefore, they are expected to produce negligible hard x-ray flux. The explanation of the hard x-ray excess in the context of an axion model has been explored in Ref. [38] and is found to be consistent with the current constraints. In this work, we explore the interpretation for the excess in the dark vector scenario.
The conversion probability between dark vector and photon is proportional to the square of the magnetic field strength. Observations show that all the M7 NSs have strong magnetic fields with a characteristic dipole magnetic field strength around ∼ 10 13 G. Therefore the NSs provide an excellent place to test the γ − γ oscillations. Having calculated both the differential luminosity of dark vectors (12) and the dark vector-to-photon conversion probability P γ →γ in Appendix B, the differential flux of hard x-ray photons produced from dark vector-photon conversion is given by where d is the distance between the NS and Earth. The probability depends on the dipole magnetic field orientation θ. In this work we adopt the θ-averaged conversion probabilitȳ P γ →γ = 1 2π 2π 0 dθP γ →γ (θ) as an approximation for our calculations. In the limit of small dark vector mass, m γ (ω/R s ) 1/2 ∼ 10 −4 eV for energy ω ∼ 1 keV and NS radius R s 11 km,P γ →γ becomes independent of m γ (see Eq. (C3) in Appendix C). Note that in order to compare with the observed hard x-ray spectrum, the differential flux F γ (ω) is calculated at energy ω = E obs /e φs , where E obs is the observed photon energy. Then the differential flux observed at infinity is given by F ∞ γ (E obs ) = e 3φs F γ (ω). Compared with the photon differential luminosity dL γ /dω whose maximum value is obtained at energy ω 2.82T s ∼ 141 eV, with a surface temperature T s ∼ 50 eV, the hard x-ray spectrum F γ (ω) from γ → γ conversion has its maximum value at a much higher energy ω 3.31T c ∼ 6.6 keV, with a core temperature T s ∼ 2 keV (the surface-core temperature relation can be found in Fig. 5).
We thus expect the flux from γ − γ conversion can make contributions to the hard x-ray excess in the energy range 2 − 8 keV and at the same time reproduce the correct spectrum shape.

VII. STATISTICAL ANALYSIS
Our starting point for the statistical analysis of the J1856 data in the context of the dark vector model is the Bayesian inference [70]. For the model consisting of a set of parameters θ = {θ 1 , θ 2 , · · · , θ n }, the Bayes theorem is given by where P(data) is the data probability, P( θ) is the prior probability that indicates the degree of belief one has before observing the data, and P( θ | data) is the conditional probability-density function that is a posterior probability representing the change in the degree of belief one can have after giving the measurement data. Note that P(data | θ) = L( θ) links the posterior probability to the likelihood of the data, where the likelihood function L( θ) takes the form of where where λ exp i denotes the experimental value with an uncertainty in the measurement σ i , m represents the number of data points, and λ the i denotes the predicted value for a given model.
For the dataset, we take into account the J1856 x-ray spectrum observations in 2 − 8 keV [48], T ∞ s [38], L ∞ γ [71][72][73], and the distance measurements [66,67]. Noting that L ∞ γ inferred from the observations is in a narrow band (0.05 − 0.08) × 10 33 erg/s, we take L ∞ γ = 0.065 × 10 33 erg/s as the central value and assume a Gaussian distribution for the luminosity. The parameter set is θ = {d, e , ε}, where the distance d enters the spectrum function. We fix the mass of dark vector m γ = 10 −6 eV to establish the approximation of conversion probability, i.e., Eq. (C3). Note that our analysis is independent of m γ as long as m γ (T s /R s ) 1/2 ∼ 3 × 10 −5 eV, with surface temperature T s ∼ 50 eV.
The performance of the Bayesian statistical analysis is carried out using the UltraNest package [51], which implements a nested sampling Monte Carlo technique [74]. It computes the (log-)likelihood as well as the marginal likelihood ("evidence") Z to perform model comparison. Meanwhile, the posterior probability distributions on model parameters are constructed to describe the parameter constraints of the data. We assume a uniform prior distribution for the distance and a log-uniform prior distribution for e and ε. We choose the number of live points to be 20000, and the total number of the samples called in the fitting is 225802.
In each running, we simulate the NS cooling based on the modified NSCool code that includes the dark vector emissivity. We then determine the surface luminosity as well as the surface temperature at the age of J1856. We obtain a minimum value (corresponding to the maximum value of the likelihood) of the reduced chi-square χ 2 min /DOF = 1.02/3 in the fitting to the data. On the other hand, when the emissivity of dark vector is turned off in the fitting, we obtain a value of χ 2 0 /DOF = 9.38/6 with the NS standard cooling model. We thus conclude that we have made the fitting much better by taking into account the emission of dark vectors. Therefore, the data supports the existing of a dark vector with couplings summarized in Table I.
The corner plot in Fig. 8 shows the posterior distributions of the parameters. In Fig. 9, we show the predictions of the model for the J1856 x-ray spectrum. The black curve denotes the prediction of the model with the median values, while the black band represents the 1σ confidence interval. Figure 10 depicts the evolution of the surface luminosity and temperature observed at infinity with the NS age. In this plot, we take the parameters to their best-fit values. As indicated in these figures, the x-ray spectrum as well as the luminosity observations can be well interpreted in the B − L dark vector model. We note that all of the PBF processes have been turned off in our simulations. Including these processes would lead to much lower surface luminosities than the observations.

VIII. MODEL CONSTRAINTS
The likelihood ratio test [75] is used to determine the limit on and the significance of a possible dark vector contribution to the J1856 observations. To do this, we fix e and determine the "nuisance parameter" d by minimizing the chi-square. Upper limits at the 95% confidence level on ε are derived by increasing the chi-square from its minimum value of the model until it changes by 2.71, χ 2 upp (ε) = χ 2 min (ε) + 2.71, withε denoting the parameter that minimizes the chi-square at fixed e . The black curve in Fig. 11 represents the constraints on e − ε and the grey region is excluded by the constraint. The green and yellow regions represent the parameter spaces that are favored by the observations at 1σ and 2σ confidence intervals. The conclusions inferred from this figure are summarized as follows: • There is an upper limit on the mixing angle ε < 7.97 × 10 −9 , which is independent of the gauge coupling e and dark vector mass m γ (as long as m γ 3 × 10 −5 eV). This is because the observed surface luminosity would be strongly suppressed by the γ − γ conversion when ε 10 −8 , as illustrated in Fig. 6. This constraint can also be applied to the dark photon model.
• There is an upper limit on the gauge coupling e < 4.13 × 10 −13 , which is independent of the mixing angle ε and dark vector mass m γ (as long as m γ 1 keV). This is because the surface luminosity would be strongly reduced due to the cooling of the NS by the dark vector emissivity with e 5 × 10 −13 , as shown in Fig. 7.   In Fig. 12, we compare our constraints with those limits from cosmological, astrophysical,  [80]. Using the phenomenon of light shining through a wall for dark photons, the grey region has been bounded by the experiment CROWS [81] at CERN. We observe that our constraint (light-blue region) on ε from J1856 surface luminosity observation is much stronger than the above constraints by about an order of magnitude.
In the right plot of Fig. 12, we summarize the constraints on the gauge coupling e .
The exchange of light bosons, such as gauge bosons, scalar axions, and dilatons among others generates a Yukawa potential, which is tested by the fifth force experiments [82][83][84].
The colored regions are excluded by the constraints.
Constraints from J1856 surface luminosity observation are much stronger than those given by the Sun luminosity observations and the DM direct detections, for the dark vector with mass 1 eV. In summary, we conclude that the 2σ parameter space with m γ 10 −5 eV, for explaining the hard x-ray excess, remains available when various observation constraints on the dark vector model are taken into account.
In Table II we summarize the current limits on hard x-ray intensity for NSs near the galactic plane (J1856, J0806, J0720, and J2143) from the 105-month Swift Burst Alert Telescope all-sky hard x-ray |b| ≤ 17.5 • survey [88] and the 14-year INTEGRAL galactic plane survey [89]. We also adopt the projected sensitivity (provided in the supplementary materials of Ref. [38]) at 95% confidence level for a 400 ks NuSTAR observation of J1856 in two energy bands. The third column provides the predicted intensities, assuming the median values of the parameters for the dark vector model. We observe that the predicted intensities are far below the current limits from Swift and INTEGRAL, and future observations on M7   by the NuSTAR experiment may be useful to test the dark vector model.

IX. SUMMARY AND CONCLUSIONS
The NS is recognized as one of the most excellent astrophysical laboratories for searching for new light particles that couple weakly to SM particles. In this work, we have presented the emission of light dark vectors from the nucleon bremsstrahlung processes in the core of the NS. The dark vector is assumed to be a B − L gauge boson and to have a mass much below about keV, the core temperature of the NS. The dominant production mode of dark vector in the NS core is the neutron bremsstrahlung since the production of dark vector by the charged particles in the NS core is suppressed by a factor of m 2 γ from the in-medium effect. Since the plasma is thought to be dilute and nonrelativistic in the Sun and supernova, the nondegenerate limit was employed to determine the dark vector emissivity and obtain constraints in previous literature. In the current work, we present the calculation of the emission rate of dark vectors from the nucleon bremsstrahlung processes in the degenerate limit, which is the case for the strongly-compressed circumstances in the NS core. In addition, we also calculate the photon luminosity observed at infinity by taking into account the photon-dark vector conversion during their propagation.
In this work, we attempt to interpret the J1856 hard x-ray spectrum excess in terms of the dark vector model while taking into account the J1856 surface luminosity and temperature observations. The thermal photon spectrum makes negligible contribution to the hard x-ray observations since the energy of the thermal spectrum peak is determined by the surface temperature T s ∼ 50 eV. On the other hand, the peak of the spectrum from γ → γ conversion is obtained at a higher energy that is determined by the core temperature T c ∼ 2 keV. Therefore, we expect that the dark vector emission can lead to the hard x-ray excess.
The evolution of the NS with time would be altered if its energy loss was dominated by the production of dark vectors rather than the standard neutrino emission. We perform numerical simulations of the NS cooling based on the modified NSCool code that includes additional energy loss via the dark vector emission, assuming the APR EOS for the NS core with a canonical mass of 1.4 M . In this way, we determine the surface luminosity, as well as the surface temperature, for the NS with a given age. Then we calculate the hard x-ray spectrum from γ → γ conversion and consider the inverse conversion for the surface luminosity observed at infinity. We carry out the Bayesian statistical analysis of the J1856 data using the UltraNest package and employ the likelihood ratio test to construct 95% confidence level upper limits on the parameters of the dark vector model. Our findings are summarized as follows: • The fit to the J1856 hard x-ray excess data favors the dark vector model with e = 5.56 × 10 −15 and ε = 1.29 × 10 −9 . Furthermore, the 2σ parameter space with m γ 10 −5 eV for the interpretation of the hard x-ray excess does not conflict with any of the currently observed constraints.
• Due to the γ → γ conversion, there exists an upper limit on the mixing angle, ε < 7.97 × 10 −9 , for m γ (T s /R s ) 1/2 ∼ 3 × 10 −5 eV from the J1856 surface luminosity observation. This constraint is independent of the production of the dark vectors, and therefore, is independent of the gauge coupling e , and thus can be applied to the dark photon model. This Appendix briefly reviews the results of photon self-energies in plasmas given in Refs. [90,91]. To leading order in the EM coupling constant α, the EM polarization tensor Π µν that determines the effects of a plasma on the propagation of photons is given by [90] Π µν (K) = 16πα where f e and fē are the Fermi distribution function for electron and positron, E = p 2 + m 2 e , K µ = (ω, k), P µ = (E, p), K 2 = ω 2 − k 2 , and P · K = Eω − p · k. The integration over the angular parts can be performed by ignoring the K 4 term. Taking into account the degenerate where v * denotes the typical electron velocity in the plasma, ω p is the plasma frequency, which is dominated by the electrons with E 2 F,e = m 2 e + 3π 2 n e 2/3 .
Under the high density circumstance in the NS core, the electron Fermi momentum p F,e = (3π 2 n e ) 1/3 ∼ O(100) MeV m e , T . Finally, the function, In the Coulomb gauge, the transverse and longitudinal components of the effective propagator for the EM field are given by With the explicit expression of the photon propagator, one can find that the emission of dark vectors is given by the vacuum matrix element for the emission of massive photons and multiplied by the effective coupling given by Eq. (3) [25], i.e., This equation shows that the dark vector produced by the EM charged current is suppressed by the dark vector mass. While for neutral currents, the effective coupling e n ef f = e , thus the dark vector produced from this process in medium is the same as that in the vacuum and the production is not suppressed [25].

Appendix B: Conversion probability
The dark vectors emitted from the neutron bremsstrahlung processes may be converted into x-ray photons as they propagate outwards through the magnetosphere around the magnetized NS. The stellar magnetic field is assumed to be dipolar where B 0 is the magnetic field at the surface of the NS, and r 0 is the NS radius. The evolution of the photon and the dark vector in the presence of an external magnetic field can be described in terms of the following system of first-order differential equations [49,92] where the term ∆ γ = − m 2 γ 2ω (ω being the energy of the fields) is due to the finite dark vector mass. The condition m γ ω is required to satisfy the weak-dispersion limit [32]. These equations can describe the propagations of both the perpendicular modes and the parallel modes. Strong-field QED effects in vacuum give rise to the term [93] where θ is the angle between the direction of propagation and the magnetic field, and q is a dimensionless function of b = B(r)/B c (where the critical QED field strength B c ≡ m 2 e c 3 /(e ) = 4.414 × 10 13 G) given by [92,94] where ⊥ and denote the perpendicular and parallel modes of the dark vectors, respectively.
These formulae have the correct b 1 and b 1 limits. Since the observer is far away from the source, we take the latter limit, in which we haveq ⊥/ 1. Note that these results are only valid for photon energies below the electron mass m e ∼ 500 keV, which is applicable for dark vector photon with energies in the hard x-ray frequency regime.
We follow Ref. [92] to resolve differential equations (B2) in the weak-mixing limit. Equation (B2) can be rewritten as a "Schrödinger equation": where A = (A, A ) and the "Hamiltonian" matrices are The Schrödinger equation in the interaction representation is given by where The "transformation operator" is defined as where r 0 is the initial position. The solution at order n + 1 can be obtained order by order: with the zero-order solution A 0 int (r) = A(r 0 ). For the first-order solution, we have The dark vector-photon conversion probability at the first-order (in the interaction representation) is Obviously, the conversion probability in the Schrödinger representation equals to Eq. (B15) with transformations (B9).

Appendix C: Analytical results for conversion probability
In this section, we show some analytical results for the dark vector-photon conversion probability to directly see how the probability is enhanced by a strong magnetic field. Since both modes of dark vector obey the same equations of motion (B2), we will focus on the parallel mode below. The dark vector-photon conversion probability in the weak-mixing limit is given by where y = 7α 45π B 0 Bc 2 ω sin 2 θ. We plot P γ →γ as a function of m γ and θ in Fig. 13. In both plots, we take B 0 = 10 14 G and ε = 10 −12 . As shown in the left plot of Fig. 13, the probability For small values of m γ , the oscillation term ∆ γ r can be neglected and the conversion probability is approximated as Averaged over θ, the probability can be parametrized as (C3) Note that we have used the θ-averaged probability in our calculations for the luminosity and spectrum. In the presence of an inhomogeneous external field, the probability is proportional to ω 2 for the dark vector-to-photon conversion [49], but is inversely proportional to the frequency in the zero background field limit [95][96][97]. Furthermore, when the external field is removed, the conversion probability approaches zero in the limit of low dark vector mass since in this case the probability is proportional to m 2 γ [95][96][97]. However, for the case of an inhomogeneous external field, the dark vector-to-photon conversion probability dose not depend on m γ in the low dark vector mass limit [49], which also appears in the axion-photon conversion [38].
Appendix D: Dark vector emission in strong degenerate plasma For the nucleon-nucleon bremsstrahlung emission of dark vectors in the NS core, we can safely assume a degenerate limit because of T c O(MeV) in the core of the NS. Let us first consider the process n + p → n + p + γ , the dark vector energy emission rate is given by the formula (6). Multiply the equation by one in the form [37] 1 = ∞ 0 dp 1 dp 2 dp 3 dp 4 δ ( where p k ≡ |p k | and p F n and p F p are the Fermi momenta for neutron and proton, respectively, and p k dp k = E * k dE * k has been used. The squared matrix can be expanded as where ξ = 2 for matrix (9) and (11). The energy emission rate (6) can be written as where the energy integral is In the strong degeneracy limit µ j /T → ∞, which is a good approximation for the processes taking place in the NS core, the energy integral is given by [31] The angular integral is given by The dark vector momentum p A has been neglected in the momentum-conserving δ-function.
Since the squared matrix element is in general a function of the momentum transfer k and l, we can convert the integral to k and l by inserting the unity [41] 1 = d 3 kd 3 lδ 3 (k − p 2 + p 4 ) δ 3 (l − p 2 + p 3 ) (D7) into the right-hand side of Eq. (D6). We can eliminate p 2 , p 3 , and p 4 one by one with the aid of the three 3-momentum-conserving δ-functions. Then where we have relabeled p 1 as p. In order to evaluate the integration, we choose the spherical coordinates to expand the three vectors as p = p(0, 0, 1), k = k(sin θ k , 0, cos θ k ), l = l(sin θ l cos φ l , sin θ l sin φ l , cos θ l ). (D9) Namely, p lies along the z axis, k lies in the x − z plane, and l points to an arbitrarily direction. Then we have |p + k| 2 = k 2 + 2kp cos θ k + p 2 , |p + l| 2 = l 2 + 2lp cos θ l + p 2 , |p + k + l| 2 = p 2 F p + 2kl sin θ k sin θ l cos φ l + 2kl cos θ k cos θ l .
We see that the term proportional to k · l does not contribute. This is the consequence of strong degeneracy limit.
The δ-functions in Eq. (D8) can be written as with the definitions x k = cos θ k , x l = cos θ l , and x φ = cos φ l , and . (D17) Using the above relations, the integration (D8) can be simplified as A = 32π 2 p 2 F p p 2 F n dk dl la(k, l) with the constraints on the parameters −1 ≤ x k , x l , x φ ≤ 1 and k, l > 0.
Similarly, for the process n + n → n + n + A the angular integral is given by (D19)

Appendix E: Remarks on OPE approximation
In order to realistically determine the production rate of new bosons (e.g. axion or dark vector) from nucleon-nucleon bremsstrahlung processes, the simplified treatment based on one-pion exchange and the use of the Born approximation for the nucleon-nucleon interaction was originally introduced in Refs. [28,30,53,98,99]. It has been realized in the literature [21,36,55,100,101] that such a simplified method ignores some relevant effects such as the multiple nucleon scatterings, and leads to a larger emission rate and thus a stronger constraint on the coupling.
One way to go beyond the OPE approximation has been performed in Refs. [55,100] by including a nonperturbative treatment of the two-nucleon scattering in the soft-radiation approximation. It was found that for the range of conditions encountered in a NS, this treatment results in an approximate modification of the axion emissivity by a factor of 1/4 when compared with the OPE results [54]. Using the soft-radiation approximation, Ref. [21] calculated the dark vectors emissivity from a supernova and found a factor of 10 decrease in the emission rate. Reference [36] refined the calculation based on the OPE approximation by systematically taking into account the effects that had been neglected previously, including the contribution of the two-pions exchange, effective in-medium nucleon masses and multiple nucleon scatterings. They found that the axion emissivity from a supernova was reduced by ∼ 10 with respect to the basic OPE calculation. From these results we observe that the reduction of emissivity depends on the condition of matter. For the highly compressed matter in a NS, the emissivity is reduced by a factor of 4 with respect to the OPE approximation, while the factor can be up to 10 for the dilute plasma in a supernova.
Further improvement of the soft-radiation approximation can be made by including the many-body effects, such as the Landau-Pomeranchuk-Migdal (LPM) suppression [102,103] owing to multiple scatterings in the medium. It has been shown in Ref. [104] that the LPM suppression of soft radiations increases with temperature and becomes relevant for T ∼ 5 MeV, above which the energy of the emitted boson is less than the nucleon decay width. Thus, the LPM effect is of particular importance for radiation in a supernova, but it becomes insignificant for the isolated NSs with a core temperature T ∼ 10 keV.