Probing dark gauge boson with observations from neutron stars

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 Planck2016AA . The most popular class of DM candidates is the weakly interacting massive particles (WIMPs) Lee1977PRL ; Hut1977PLB 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 XENON1T2018 and indirect Ackermann2015PRL 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 Holdom1986PLB ; Chiang2020PRD , which can naturally arise from the breaking of grand unified theory (GUT) groups Dienes1997NPB ; Ibe2019JHEP or the string theory Dienes1997NPB ; Lukas2000JHEP ; Blumenhagen2005JHEP ; Abel2008JHEP ; Burgess2008JHEP ; Goodsell2009JHEP ; Burrage2009JCAP ; Cicoli2011JHEP . 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 Fabbrichesi2020 .

Supernovae serve as a powerful factory for the emission of dark photons with masses $\lesssim 100$ MeV  Bjorken1009PRD ; Dent2012 ; Kazanas2014NPB ; Rrapaj2016PRC ; Chang2017JHEP ; Hardy2017JHEP ; DeRocco2019JHEP ; Hong2020 . 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. Chang2017JHEP ; Hardy2017JHEP ; DeRocco2019JHEP , 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 $\lesssim$ keV. For instance, Refs. An2013PLB ; An2013PRL show that the measured luminosity of the Sun, $L_{\odot}=3.83\times 10^{26}$ W, has constrained the dark photon parameter space to the range of $\varepsilon m_{\gamma^{\prime}}<4\times 10^{-12}~{}{\rm eV}$ An2013PRL .

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 Iwamoto1984PRL ; Morris1986PRD ; Brinkmann1988PRD ; Iwamoto2001PRD ; Lai2006PRD ; Fortin2018JHEP ; Sedrakian2016PRD ; Sedrakian2019PRD ; Carenza2019JCAP ; Harris2020JCAP ; Buschmann2021PRL ; Leinson2019JCAP ; Leinson2021JCAP ; Fortin2021 . The core of a NS is composed of highly degenerate nuclear matter that has an average density of $\sim 3\times 10^{14}$ $\rm g/cm^{3}$, which corresponds to an average distance of $\sim 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 Chang2017JHEP ; Hardy2017JHEP ; DeRocco2019JHEP , 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 Chang2017JHEP ; Hardy2017JHEP . 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 Davidson1979PRD ; Marshak1980PLB ; Mohapatra1980PRL ; Davidson1987PRL ; Montero2009PLB ; Ma2014PLB . 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^{\prime}$ 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 Dessert2020APJ 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. Buschmann2021PRL 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 Page2016 , 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 Fortin2019JCAP . We present an analytical formula of the dark vector-photon conversion probability for dark vector with mass $m_{\gamma^{\prime}}\lesssim(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 Buchner2021 . 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 photon-dark vector conversion. In Sec. V, we describe the NS simulation based on the modified NSCool code in detail. In Sec. VI, we compute the hard x-ray spectrum from the dark vector-photon conversions. The statistical analysis of data based on the Bayesian inference is carried out in Sec. VII, and the 95% upper limits on the dark vector model are derived in Sec. VIII. Finally in Sec. IX, we summarize our findings. The plasma effects on dark vector’s production are discussed in Appendix A. The dark vector-photon conversion probability is derived in Appendix B. In Appendix C, we show the analytical results for the dark vector-photon conversion probability. The detailed calculations of the dark vector emission rate are given in Appendix D. Discussions on the one-pion approximation are presented in Appendix E.

We are interested in the extension of the SM with an additional $U(1)$ gauge symmetry, having a new gauge boson $A_{\mu}^{\rm D}$ called the dark vector. The relevant effective Lagrangian at low scales is written as

where we assume the dark vector has a Stueckelberg mass $m_{\gamma^{\prime}}$, $F_{\mu\nu}=\partial_{\mu}A_{\nu}^{\rm SM}-\partial_{\nu}A_{\mu}^{\rm SM}$ and $F_{\mu\nu}^{\prime}=\partial_{\mu}A_{\nu}^{\rm D}-\partial_{\nu}A_{\mu}^{\rm D}$ are, respectively, the field strengths of the SM photon and dark gauge boson, and $\varepsilon$ represents the mixing angle between the SM photon and dark vector. The parameter $e^{\prime}$ denotes the coupling between $A_{\mu}^{\rm D}$ and the current $J^{\mu}$, 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 Davidson1979PRD ; Marshak1980PLB ; Mohapatra1980PRL ; Davidson1987PRL ; Montero2009PLB ; Ma2014PLB . The kinetic terms of gauge boson in the Eq. (1) can be diagonalized by rotating the gauge fields as

where $\left(A_{\mu}^{\rm D},~{}A_{\mu}^{\rm SM}\right)^{T}$ and $\left(A_{\mu}^{\prime},~{}A_{\mu}\right)^{T}$ are the so-called interaction eigenstates and mass (propagating) eigenstates, respectively. For a massive dark vector, the rotation angle $\theta$ is locked at zero Fabbrichesi2020 . 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 $\varepsilon 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 $\omega_{p}^{2}=4\pi\alpha_{\mathrm{EM}}\sum_{i}n_{i}/E_{F,i}$, where the sum goes over the charged particles with a number density $n_{i}$ and a Fermi energy $E_{F,i}^{2}=m_{i}^{2}+\left(3\pi^{2}n_{i}\right)^{2/3}$. The nonzero 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_{\mu}^{\prime}$ in the medium is given by Hong2020

where $f$ denotes the SM fermions with an electric charge $q_{f}$ and a dark gauge quantum number $q_{f}^{\prime}$ in units of $e$ and $e^{\prime}$, respectively. The proton and neutron have the electric charges $q_{p}=1$ and $q_{n}=0$, and the dark gauge quantum numbers $q_{p}^{\prime}=q_{n}^{\prime}=1$. The EM and dark gauge quantum numbers for the electron are $q_{e}=q_{e}^{\prime}=-1$. The visible photon mass is encoded in the polarization tensor $\pi_{T,L}$ given in Appendix A.

From Eq. (3), the effective couplings of electron, proton, and neutron to $A_{\mu}^{\prime}$ in a dense medium (i.e., with $\pi_{T,L}\gg T^{2},~{}m_{\gamma^{\prime}}^{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_{\mu}^{\prime}$ 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 Potekhin2020MNRAS 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. Buschmann2021PRL 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_{\gamma^{\prime}}\lesssim(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_{\gamma^{\prime}}$ cannot take the vanishing limit and the numerical value for $\varepsilon$ 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^{\prime}$ 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.

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.

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 Chiu1964PRL ; Bahcall1965PLBa ; Bahcall1965PLBb ; Friman1979APJ ; Maxwell1987APJ includes several neutrino emission processes (see Yakovlev2001PR 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_{\rm D}$, that is significantly larger than the canonical NS mass of $1.4~{}M_{\odot}$ Lattimer1991PRL ; Beloborodov2016APJ . For the NS with mass less than $M_{\rm 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 Friman1979APJ

where $\rho$ is the NS mass density profile, $\rho_{0}=2.8\times 10^{14}\mathrm{~{}g}\cdot\mathrm{cm}^{-3}$ is the nuclear saturation density, $T_{c}$ is the NS core temperature, and the suppression factor $R_{\rm M}\leq 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. Yakovlev1995AA by using the same formalism as that in Ref Friman1979APJ . The expression for the neutrino emissivity in the proton branch can be approximated by the following rescaling relation Yakovlev2001PR

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 $\Theta_{F}=1$ for $p_{{F},n}<3p_{{F},p}+p_{{F},e}$; otherwise, $\Theta_{F}=0$. As an example, take $p_{{F,e}}=p_{{F},p}=p_{{F,n}}/2$, we have $Q_{\nu}^{p}\simeq 0.5Q_{\nu}^{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. Yakovlev2001PR ). Summing up all of the processes, the total neutrino emission rate can be estimated as $Q_{\nu}\simeq RQ_{\nu}^{n}$, with a rescaling factor $R\simeq 1.5$ Yakovlev2001PR . 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_{\gamma^{\prime}}<L_{\nu}$ Raffelt1996 . 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 Buschmann2021PRL . In Fig. 5, we depict the relation between core temperature and surface temperature at infinity using NSCool code Page2016 (with the default NS model). We observe that the relation is independent of the new gauge coupling $e^{\prime}$ (with $\varepsilon=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}\simeq(10^{9}\mathrm{~{}K})\left(t/\mathrm{yr}\right)^{-1/6}$ Beloborodov2016APJ . 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\gtrsim 100$ kyr, the photon emission will exceed the neutrino one Raffelt1996 , 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 $\sigma$ is the Stefan-Boltzmann constant, $R_{s}$ is stellar radius, and $T_{s}$ is stellar surface temperature. Without photon-dark vector conversions, the observed surface photon luminosity at infinity is accordingly redshifted as

where ${\phi_{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 $\sim 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_{\gamma\to\gamma^{\prime}}(\omega)=P_{\gamma^{\prime}\to\gamma}(\omega)$ 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 $\omega$, the differential surface luminosity of the photon is given by

We will show in Appendix C that the $\gamma-\gamma^{\prime}$ conversion probability can be written as $P_{\gamma\to\gamma^{\prime}}(\omega)=P_{0}\omega^{2}$ for dark vector in the mass range $m_{\gamma^{\prime}}\lesssim(\omega/R_{s})^{1/2}\sim 3\times 10^{-5}$ eV, with $\omega\sim T_{s}\sim 50$ eV. With these, we finally obtain the observed surface photon luminosity at infinity by taking into account the $\gamma-\gamma^{\prime}$ conversion

where $A_{s}=4\pi R_{s}^{2}$. 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 $\gamma-\gamma^{\prime}$ conversion once we know the values of $L_{\gamma}$ 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 $\omega\gtrsim m_{\gamma^{\prime}}^{2}R_{s}$. However, we have integrated over $\omega$ in the range of $(0,\infty)$ for $P_{\gamma\to\gamma^{\prime}}dL_{\gamma}/d\omega$. It is easy for the reader to confirm that there is nearly no difference between the integration over $\omega$ in the ranges of $(0,\infty)$ and $(0.9T_{s},\infty)$, which means that the small values of $\omega\lesssim 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_{\gamma^{\prime}}\lesssim 3\times 10^{-5}$, with $\omega\sim T_{s}\sim 50$ eV.

In Fig. 6, we show the variation of the surface luminosity observed at infinity $L_{\gamma}^{\infty}$ with the mixing angle $\varepsilon$. We assume the surface luminosity $L_{\gamma}=8.4\times 10^{31}$ erg/s and temperature $T_{s}=47.1$ eV for the NS. As shown in this figure, the $\gamma-\gamma^{\prime}$ conversion strongly reduces the observed surface luminosity when the mixing angle reaches values of $\sim 10^{-8}$.

In this work, we employ the NSCool code Page2016 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, nucleon-nucleon 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_{\odot}$ (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\pm 13$ pc away from us and has an age around $(4.2\pm 0.8)\times 10^{5}$ yr Kerkwijk2008PAJ ; Walter2010APJ (a summary can be found in Refs. webNScool ; Potekhin2020MNRAS ). Recent studies Potekhin2020MNRAS 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. Buschmann2021PRL . 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 $e^{\prime}=10^{-13},~{}5\times 10^{-13}$, and $10^{-12}$, and the results are represented by the solid, dashed, and dotted curves, respectively. The blue, dark, and yellow curves represent the neutrino, dark vector, and photon surface luminosities, respectively. The red point represents the observation of the J1856 surface luminosity. It is shown by the simulations that the photon luminosity is the dominant component for J1856. We find that the model with a new gauge coupling $e^{\prime}\lesssim 10^{-13}$ predicts the same age-luminosity profile as the standard NS cooling model since in this case the cooling by the dark vector emission is negligible. There is a significant deviation in the photon luminosity from the observation when $e^{\prime}$ increases to $5\times 10^{-13}$ and, therefore, $e^{\prime}$ should be constrained by the observation. We will show the constraints on the model by fitting to the observations below.

Recent analyses of the data from the XMM-Newton and Chandra x-ray telescopes Dessert2020APJ 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. Dessert2020APJ that the NS J1856 hard x-ray spectrum has a $\sim 5\sigma$ 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 Buschmann2021PRL . Observations by the ROSAT All Sky Survey Haberl2007ASS have shown that all the M7 members have soft spectra that are well described by near-thermal distributions with surface temperatures $\sim 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. Buschmann2021PRL 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 $\sim 10^{13}$ G. Therefore the NSs provide an excellent place to test the $\gamma^{\prime}-\gamma$ oscillations. Having calculated both the differential luminosity of dark vectors (12) and the dark vector-to-photon conversion probability $P_{\gamma^{\prime}\to\gamma}$ 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 $\theta$. In this work we adopt the $\theta$-averaged conversion probability $\bar{P}_{\gamma^{\prime}\to\gamma}=\frac{1}{2\pi}\int_{0}^{2\pi}d\theta P_{\gamma^{\prime}\to\gamma}(\theta)$ as an approximation for our calculations. In the limit of small dark vector mass, $m_{\gamma^{\prime}}\lesssim(\omega/R_{s})^{1/2}\sim 10^{-4}$ eV for energy $\omega\sim 1$ keV and NS radius $R_{s}\simeq 11$ km, $\bar{P}_{\gamma^{\prime}\to\gamma}$ becomes independent of $m_{\gamma^{\prime}}$ (see Eq. (56) in Appendix C). Note that in order to compare with the observed hard x-ray spectrum, the differential flux $F_{\gamma}(\omega)$ is calculated at energy $\omega=E_{\rm obs}/e^{\phi_{s}}$, where $E_{\rm obs}$ is the observed photon energy. Then the differential flux observed at infinity is given by $F_{\gamma}^{\infty}(E_{\rm obs})=e^{3\phi_{s}}F_{\gamma}(\omega)$. Compared with the photon differential luminosity $dL_{\gamma}/d\omega$ whose maximum value is obtained at energy $\omega\simeq 2.82T_{s}\sim 141$ eV, with a surface temperature $T_{s}\sim 50$ eV, the hard x-ray spectrum $F_{\gamma}(\omega)$ from $\gamma^{\prime}\to\gamma$ conversion has its maximum value at a much higher energy $\omega\simeq 3.31T_{c}\sim 6.6$ keV, with a core temperature $T_{s}\sim 2$ keV (the surface-core temperature relation can be found in Fig. 5). We thus expect the flux from $\gamma^{\prime}-\gamma$ 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.

Our starting point for the statistical analysis of the J1856 data in the context of the dark vector model is the Bayesian inference Bayesian1992 . For the model consisting of a set of parameters $\vec{\theta}=\left\{\theta^{1},\theta^{2},\cdots,\theta^{n}\right\}$, the Bayes theorem is given by

where $\mathcal{P}(\rm{data})$ is the data probability, $\mathcal{P}(\vec{\theta})$ is the prior probability that indicates the degree of belief one has before observing the data, and $\mathcal{P}(\vec{\theta}\mid{\rm 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 $\mathcal{P}({\rm data}\mid\vec{\theta})=\mathcal{L}(\vec{\theta})$ links the posterior probability to the likelihood of the data, where the likelihood function $\mathcal{L}(\vec{\theta})$ takes the form of

where

where $\lambda_{i}^{\exp}$ denotes the experimental value with an uncertainty in the measurement $\sigma_{i}$, $m$ represents the number of data points, and $\lambda_{i}^{\text{the }}$ 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 Dessert2020APJ , $T_{s}^{\infty}$ Buschmann2021PRL , $L_{\gamma}^{\infty}$ Ho2007MNRAS ; Mignani2013MNRAS ; Yoneyama2017PASJ , and the distance measurements Kerkwijk2008PAJ ; Walter2010APJ . Noting that $L_{\gamma}^{\infty}$ inferred from the observations is in a narrow band $(0.05-0.08)\times 10^{33}$ erg/s, we take $L_{\gamma}^{\infty}=0.065\times 10^{33}$ erg/s as the central value and assume a Gaussian distribution for the luminosity. The parameter set is $\vec{\theta}=\{d,e^{\prime},\varepsilon\}$, where the distance $d$ enters the spectrum function. We fix the mass of dark vector $m_{\gamma^{\prime}}=10^{-6}$ eV to establish the approximation of conversion probability, i.e., Eq. (56). Note that our analysis is independent of $m_{\gamma^{\prime}}$ as long as $m_{\gamma^{\prime}}\lesssim\left(T_{s}/R_{s}\right)^{1/2}\sim 3\times 10^{-5}~{}\rm{eV}$, with surface temperature $T_{s}\sim 50$ eV.

The performance of the Bayesian statistical analysis is carried out using the UltraNest package Buchner2021 , which implements a nested sampling Monte Carlo technique Skilling2004 . 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^{\prime}$ and $\varepsilon$. 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 $\chi^{2}_{\rm min}/{\rm 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 $\chi^{2}_{0}/{\rm 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 1.

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\sigma$ 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.

The likelihood ratio test Rolke2005 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^{\prime}$ and determine the “nuisance parameter” $d$ by minimizing the chi-square. Upper limits at the 95% confidence level on $\varepsilon$ are derived by increasing the chi-square from its minimum value of the model until it changes by 2.71, $\chi_{\rm upp}^{2}(\varepsilon)=\chi_{\rm min}^{2}(\bar{\varepsilon})+2.71$, with $\bar{\varepsilon}$ denoting the parameter that minimizes the chi-square at fixed $e^{\prime}$. The black curve in Fig. 11 represents the constraints on $e^{\prime}-\varepsilon$ 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\sigma$ and $2\sigma$ confidence intervals. The conclusions inferred from this figure are summarized as follows:

• •

  There is an upper limit on the mixing angle $\varepsilon<7.97\times 10^{-9}$, which is independent of the gauge coupling $e^{\prime}$ and dark vector mass $m_{\gamma^{\prime}}$ (as long as $m_{\gamma^{\prime}}\lesssim 3\times 10^{-5}$ eV). This is because the observed surface luminosity would be strongly suppressed by the $\gamma-\gamma^{\prime}$ conversion when $\varepsilon\gtrsim 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^{\prime}<4.13\times 10^{-13}$, which is independent of the mixing angle $\varepsilon$ and dark vector mass $m_{\gamma^{\prime}}$ (as long as $m_{\gamma^{\prime}}\lesssim 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^{\prime}\gtrsim 5\times 10^{-13}$, as shown in Fig. 7.

• •

  For the gauge coupling in the range of $6\times 10^{-15}\lesssim e^{\prime}\lesssim 10^{-13}$, the x-ray spectrum data dominates the contributions to chi-square. We have shown that the model with parameters around the mean values $e^{\prime}=5.56\times 10^{-15}$ and $\varepsilon=1.29\times 10^{-9}$ can explain the x-ray spectrum excess in the energy band $2-8$ keV. Since the x-ray spectrum from the new physics is proportional to $(e^{\prime}\varepsilon)^{2}$, constraints on $\varepsilon$ become stronger with the increase of $e^{\prime}$.

In Fig. 12, we compare our constraints with those limits from cosmological, astrophysical, as well as terrestrial observations. In the left plot of the figure, we summarize the constraints on the mixing angle $\varepsilon$ in the low mass range $\lesssim 10^{-4}$ eV. The oscillation between the ordinary photon and the massive dark photon $\gamma\to\gamma^{\prime}$ induces deviations on the black body spectrum in the cosmic microwave background, which has been constrained by the COBE/FIRAS experiment Fixsen1996APJ ; Caputo2020PRL ; Aramburo2020JCAP ; McDermott2020PRD (light-green region). The constraints from detecting modifications of Coulomb force by the new gauge force in atomic and nuclear experiments are depicted by the orange region Jaeckel2010PRD . Using the phenomenon of light shining through a wall for dark photons, the grey region has been bounded by the experiment CROWS Betz2013PRD at CERN. We observe that our constraint (light-blue region) on $\varepsilon$ 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^{\prime}$. 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 Sushkov2011PRL ; Murata2015CQG ; Chen2016PRL . Furthermore, light boson states emitted from the Sun may be probed in the DM direct detection experiments An2020PRD ; XENON1T2020 ; PandaX2021 . We compare our NS cooling constraint from J1856 (blue region) with the constraints given by the Sun (red) Hardy2017JHEP , HB (orange) Hardy2017JHEP , XENON1T (green) An2020PRD , and the current fifth force experiments (light-grey) Sushkov2011PRL ; Murata2015CQG ; Chen2016PRL . We observe that the fifth force experiments give the most stringent constraints for $2\times 10^{-5}~{}{\rm eV}\lesssim m_{\gamma^{\prime}}\lesssim 3$ eV. 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 $\lesssim 1$ eV. In summary, we conclude that the $2\sigma$ parameter space with $m_{\gamma^{\prime}}\lesssim 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 2 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|\leq 17.5^{\circ}$ survey Oh2018APJ and the 14-year INTEGRAL galactic plane survey Krivonos2019MSAI . We also adopt the projected sensitivity (provided in the supplementary materials of Ref. Buschmann2021PRL ) 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.

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_{\gamma^{\prime}}^{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}\sim 50$ eV. On the other hand, the peak of the spectrum from $\gamma^{\prime}\to\gamma$ conversion is obtained at a higher energy that is determined by the core temperature $T_{c}\sim 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_{\odot}$. 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 $\gamma^{\prime}\to\gamma$ 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^{\prime}=5.56\times 10^{-15}$ and $\varepsilon=1.29\times 10^{-9}$. Furthermore, the $2\sigma$ parameter space with $m_{\gamma^{\prime}}\lesssim 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 $\gamma\to\gamma^{\prime}$ conversion, there exists an upper limit on the mixing angle, $\varepsilon<7.97\times 10^{-9}$, for $m_{\gamma^{\prime}}\lesssim(T_{s}/R_{s})^{1/2}\sim 3\times 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^{\prime}$, and thus can be applied to the dark photon model.

• •

  The emission of the dark vectors accelerates the NS cooling and reduces the surface luminosity of J1856, which leads to an upper limit on the gauge coupling, $e^{\prime}<4.13\times 10^{-13}$, at 95% confidence level for $m_{\gamma^{\prime}}\lesssim 1$ keV.

Our best-fit dark vector model predicts much lower hard x-ray intensities than the current limits from the Swift and INTEGRAL hard x-ray surveys. Future hard x-ray observations of J1856 by NuSTAR in particular may constrain or provide additional evidence for the best-fit dark vector from this work. NSs are promising targets for testing the weakly-interacting light dark vector particle. The constraints on the dark vector model can be further improved with more hard x-ray observations of the NSs from future telescopes.