Constraints on Axions from a Relativistic Model of Spatially Extended Gamma-Ray Emission from Neutron Stars

Bijan Berenji

SLAC National Accelerator Laboratory, Menlo Park, CA, USA.

DOI: 10.4236/jmp.2024.1511082                

Abstract

Axions are hypothetical particles proposed to solve the strong CP problem in QCD and may constitute a significant fraction of the dark matter in the universe. Axions are expected to be produced in superfluid neutron stars and subsequently decay, producing gamma-rays detectable by the Fermi Large Area Telescope (Fermi-LAT). Considering that light QCD axions, as opposed to axions > 1 eV, may travel a long range before they decay into gamma rays, neutron stars may appear as a spatially extended source of gamma rays. We extend our previous search for gamma rays from axions, based on a point source model, to consider the neutron star as an extended source of gamma rays. The extended consideration of neutron stars leads to higher sensitivity to searches for axions, as it will be shown. We investigate the spatial emission of gamma rays using phenomenological models of neutron star axion emission. We present models including the fundamental astrophysics and relativistic, extended gamma-ray emission from axions around neutron stars. A Monte Carlo simulation of the LAT gives us an expectation for the extended angular profile and spectrum. For a source of $\simeq$ 100 pc, we predict a mean angular spread of $\simeq$ 2˚ with gamma-ray energies in the range 10 - 200 MeV, due to the cutoff of the spin-structure function $S_{\sigma }\left(\omega \right)$ . We demonstrate the feasibility of setting more stringent limits for axions in this mass range, excluding a range not probed by observations before. We consider projected sensitivities for mass limits on axions from RX J1856-3754, a neutron star at a distance of 130 pc. The limit based on 7.9 years of Fermi-LAT data is 3.9 meV for an inner temperature of the neutron star of 20 MeV.

Keywords

Axions, Particle Astrophysics, Dark Matter

1. Introduction

The QCD axion is being investigated for compelling theoretical reasons, and many promising methods have been investigated for its detection. The axion was postulated to solve the strong CP problem [1], by the mechanism of spontaneously broken U(1) PQ symmetry, and may constitute a significant fraction of the dark matter in the universe [2]. Axions can be studied by means of neutron stars, from which they are theorized to be produced by nucleon-nucleon bremsstrahlung. The coupling of axions to the electromagnetic field can also generate axions by the Primakoff effect [3].

Observations with the Fermi-LAT are crucial to an axion search or setting limits on axion parameters. Here, for the first time, we use Fermi-LAT observations of neutron stars with an extended source model to search for signatures of axions. The Fermi-LAT is an imaging, wide field-of-view, pair-conversion telescope that detects gamma rays with energies from 20 MeV to over 300 GeV [4]. This energy range includes the energies of photons from decaying axions, roughly 30 to 200 MeV in our model.

Extended gamma-ray sources have been extensively studied with the Fermi-LAT, e.g. Ref. [5], including pulsar wind nebulae and supernova remnants. A search for extended sources in the galactic plane, detecting 46 sources, has also been performed with the Fermi-LAT [6]. In addition, a search for extended high-latitude sources has found 24 sources that demonstrate extension [7]. Further, dark matter in galaxies may be modeled as extended sources of gamma rays, e.g. Ref. [8]. In addition, Andromeda has recently been observed as an extended source [9]. We may note that spatially extended emission from axions may occur in the vicinity of supernova remnants. Decays that occur at a distance from supernova remnants have been considered in Ref. [10]. In addition, there have been many investigations of photon-axion and photon-ALP conversion in extragalactic magnetic fields, with large distance scales, e.g. Ref. [11]. We consider the decays, but do not consider the oscillations, because the distances are not large [12]. Here, we consider variation on the point-source neutron star model considered previously, and consider extended emission due to axions decaying at a certain distance away from the source. As the distribution of gamma rays arising from axion decays falls off not as rapidly as a point source, according to this theory, we model the distribution of axions as a spatially extended source. According to the theory of convolution, the convolution of a flux with a delta function gives back the original flux function; however, the convolution of a flux with a distribution with a width of 2˚ - 3˚ is noticeably different, even if the Point-Spread Function (PSF) is larger than the width. If there is no signal detected, the limits may potentially be improved with respect to the point source analysis of neutron stars, due to the photons from axion decay potentially being spread out over a larger solid angle in the Regions of Interest (ROIs) corresponding to the neutron stars we wish to investigate.

There are theoretically and observationally justifiable reasons for investigating this model. The theoretical lower bound on axion decays from supernova energy-loss arguments has been placed at ~10 meV [13]. The 100 meV to 1 meV range has been mentioned as a promising region for future axion searches [14] [15]. The possibility of diffuse emission from axions produced by NN-bremsstrahlung in supernova cores has been theorized to yield axion mass limits in the meV range [14]. We provide a model that leads to more restrictive constraints on the axion mass when considering extended axion emission from neutron stars. There are recent theoretical constraints from neutron star cooling that predict axion masses in a close mass range [16]. It is possible to detect gamma rays in the energy range of 30 - 200 MeV with the Fermi-LAT [17]. Furthermore, the model projects a flux from 10 meV axions that can be measured by Fermi-LAT, as it will be shown.

Several investigators have recently studied axions via nucleon coupling, such as with axion decay, neutron star cooling, and solar axions. In Ref. [18], a sample of 17 isolated neutron stars was investigated with similar spectral models for axion decay as determined by Ref. [19], to determine a 95% C.L. upper limit on $m_a<9.6\times {10}^{-3}$ eV. From Ref. [20], an upper limit of 0.01 eV was placed from neutron star cooling of the neutron star at the center of the supernova remnant Cas A, relying upon a model of neutron superfluidity. In Ref. [21], the 8.41 keV line for ¹⁶⁹Tm resonant excitation was studied for solar axions produced via the Primakoff effect, yielding a 90% C.L. upper limit of 24 eV.

The relation between $m_a$ and the axion coupling $f_a$ is given by:

$m_a=\frac{0.6\text{meV}}{f_a/\left({10}^{10}\text{GeV}\right)}$ (1)

Our search for axions from neutron stars depends on the axion-coupling to quarks via NN-bremsstrahlung, where the derivative couples to the axion field in the Lagrangian as:

$ℒ\subset \frac{1}{f_a}{g}_{ann}\left({\partial }_{\mu }a\right)\overline{N}{\gamma }^{\mu }{\gamma }_{5}N.$ (2)

The axion-nucleon coupling may be parametrized in terms of $m_a$ :

$g_{ann}={10}^{-8}\left(\frac{m_a}{1\text{eV}}\right).$ (3)

We consider KSVZ axions [22] [23], to be distinguished from the DFSZ axion model [24] [25]. According to the KSVZ “hadronic” axion model, the heavy quarks are electrically neutral and carry PQ charges. On the other hand, in the DFSZ model, there are at least two Higgs doublets and ordinary quarks have PQ charges [26]. The axion field should be a Bose-Einstein Condensate (BEC) [27] and is expected to be responsible for the nucleon Electron Dipole Moment (EDM).

In this article, we investigate the spatial emission of gamma rays using phenomenological models in order to determine the projected sensitivities of Fermi-LAT observations from photon fluxes of neutron stars. We present the fundamental astrophysical model, the model for extended gamma-ray emission from axions around neutron stars, the Monte Carlo simulation model. We demonstrate the feasibility of setting more stringent limits for QCD axions than previous literature values [28], which could potentially exclude a range not probed by observations before.

There have been studies in the literature with RX J1856-3754 as a superfluid neutron star [29]. The superfluidity gives the surface temperature as 100 eV, which, allowing for a modest temperature gradient, would give a core temperature of 10 MeV. In addition, the phase diagram of superfluidity in the neutron star gives temperatures of order $\mathcal{O}$ 10 MeV.

The following is the organization of this article. In Section 2, we present the underlying astrophysical model. In Section 3, we discuss the extended emission of gamma rays due to axions from neutron stars. In Section 4, we discuss projected limits from neutron star RX J1856-3754. In Section 5, we discuss the relevance of the limits in the astrophysical context as well as to other astrophysical limits on the axion mass.

2. Astrophysical Model

Axions may be produced in neutron stars by the NN-bremsstrahlung reaction $nn\to nna$ , where $n$ is a neutron [30]. The axions produced in this manner would be relativistic (see below). For a physical description of this process, we follow the phenomenology of Hanhart, Philips, and Reddy [31], who model this process as a nucleon-nucleon scattering process. Furthermore, the quarks $|n\rangle =|u\rangle |d\rangle |d\rangle$ are free, inside a neutron, according to the principle of asymptotic freedom. We developed an astrophysical model to derive an energy flux from axions emitted from neutron stars, which subsequently decay to photons in Ref. [19]. In deriving the differential photon flux (Φ), we consider the differential emissivity with respect to axion energy. In the case of radiative decay of axions $a\to 2\gamma$ , we determine the photon energy from the axion mass, the relativistic boost $\gamma$ , and the angle of photon direction with respect to the axion direction, ${\theta }_a$ . In addition, we consider a neutron star of volume $V_{NS}$ as a uniform density sphere with a radius of 10 km, a timescale for axion emission $\text{Δ}t$ as described below, a neutron star at a distance $d$ , and the axion decay width ${\text{Γ}}_{a\gamma \gamma }$ . We consider ${\text{Γ}}_{a\gamma \gamma }$ as given by [32]:

${\text{Γ}}_{a\gamma \gamma }=1.1\times {10}^{-24}{\text{s}}^{-1}{\left(\frac{m_a}{1\text{eV}}\right)}^5.$ (4)

The energy flux is related to the axion emissivity from nucleon-nucleon bremsstrahlung, as well as the timescale of axion emission from the nuclear medium, which both depend on the axion mass:

$E\frac{\text{d}\Phi }{\text{d}E}=2\frac{\text{d}{ϵ}_a}{\text{d}\omega }\delta \left(E-\omega /2\right)\frac{V_{NS}\text{Δ}t{\text{Γ}}_{a\gamma \gamma }}{4\pi d^2},$ (5)

where $\omega$ is the energy of emitted axions, ${ϵ}_a$ is the axion emissivity of the neutron star matter. This phenomenologically accounts for the nucleon-nucleon bremsstrahlung process as a nucleon-nucleon scattering. $V_{NS}$ is the volume of the neutron star, and $d$ is the distance to it. We model the timescale of axion emission from a neutron star as the mean time $\text{Δ}t$ between successive axion emissions in the nuclear medium. It is shown in Ref. [19] that $\text{Δ}t\simeq 23.2$ s. Upon simplification of Equation (5), we obtain the following equation for the spectral energy distribution:

$\begin{array}{c}E\frac{\text{d}\Phi }{\text{d}E}=1.8\times {10}^{-2}{\left(\frac{m_a}{\text{eV}}\right)}^5\left(\frac{\text{Δ}t}{23.2\text{s}}\right){\left(\frac{100\text{pc}}{d}\right)}^2\\ \times {\left(\frac{2E}{100\text{MeV}}\right)}^4\left(\frac{S_{\sigma }\left(2E\right)}{{10}^7{\text{MeV}}^2}\right){\text{cm}}^{-2}\cdot {\text{s}}^{-1}.\end{array}$ (6)

where $S_{\sigma }\left(\omega \right)$ is the spin-structure function, which accounts for the energy and momentum transfer and includes the spins of the nucleons. The function $g\left(\omega \right)={\omega }^4{S}_{\sigma }\left(\omega \right)$ , which we use in our extended analysis, is shown in Figure 1. For the purpose of this investigation, this function has been fit to an analytic functional form, according to a log-likelihood minimization using the MINUIT optimizer, given by:

$g\left(\omega \right)=\left(\alpha +\beta \omega \right)\gamma \exp \left(-\frac{1}{2}{\left(\left(\omega -\delta \right)/ϵ\right)}^2\right)\exp \left(\eta +\theta \omega \right)$ (7)

In Table 1, we present the values of the parameters of $g\left(\omega \right)$ .

Figure 1. The function $g\left(\omega \right)={\omega }^4{S}_{\sigma }\left(\omega \right)$ , whose parameters have been fit to an analytic functional form given by Equation (7). Dots refer to Monte Carlo simulated points, while the solid line refers to fitted model described by Equation (7).

We may consider the effect on axion mass limits due to variations in the model parameters. In the model of neutron stars that we are considering [33], we may consider $T=20$ MeV, and $\mu /T=10$ [33]. In the model, the neutron star matter has superfluidity, and its phase diagram is described by QCD. Since the observed surface temperature of RX J1856-3754 is 0.28 keV [34], our assumption of the core temperature of $T=20$ MeV with a modest temperature gradient can be justified. The source RX J1856-3754 is an isolated neutron star that is detected in X-rays but not in γ-rays [34]. Furthermore, superfluidity has been demonstrated to reduce the late-time cooling of neutron stars, up to an age of 10⁹ yr., by the mechanism of frictional heating, as well as suppressing the neutrino emission energy loss mechanism [35].

Table 1. Best-fit parameters of the function $g\left(\omega \right)$ as parametrized by Equation (7).

Parameter	Best-fit value
$\alpha$	−1.29 × 1011
$\beta$	−0.364
$\gamma$	−1.37 × 104
$\delta$	31.3
$ϵ$	37.0
$\eta$	22.6
$\theta$	3.49 × 10−3

3. Extended Emission of Axions from Neutron Stars

Axions decay with finite width ${\Gamma }_{a\gamma \gamma }=1.1\times {10}^{-24}{\text{s}}^{-1}{\left(\frac{m_a}{1\text{eV}}\right)}^5$ ; the probability that axions decay from the point of emission from the neutron star increases with distance from the source. In other words, the survival probability decreases with angular distance from the neutron star. Thus, the gamma rays arising from the axion decay would render the neutron star as an extended source in gamma rays. The differential survival probability $\text{d}P$ is related to the probability $P$ as follows:

$\text{d}P=-P\frac{{\Gamma }_{a\gamma \gamma }}{\beta \gamma c}\text{d}r.$ (8)

In the preceding equation, we divide $\text{d}r$ by $\beta c$ , to obtain the time to traverse a radial distance $\text{d}r$ . We also divide the decay rate ${\text{Γ}}_{a\gamma \gamma }$ by $\gamma$ and $\beta$ , the commonly used relativistic parameters, to account for time dilation.

$P_0=\frac{{\Gamma }_{a\gamma \gamma }}{\beta \gamma c}$ (9)

The energy conservation condition, that the sum of kinetic energy plus potential energy of a radiated axion be equal to the energy radiated by the axion from the neutron star, can be expressed as:

$K+U=-E_{rad}$ (10)

The preceding equation leads us to the following expression in terms of energy densities:

$\frac{1}{2}\rho {\beta }^2-\frac{G\rho M_{NS}}{c^{2}r}=-\frac{{\displaystyle \int \text{d}t}{\displaystyle \int \text{d}\Omega }{\displaystyle {\int }_0^{R_{NS}}\text{d}r}{r}^2{ϵ}_a}{4/3\pi r^3{c}^2},$ (11)

where $\rho$ is the mass density of axions, $M_{NS}$ is the mass of the neutron star, and $G$ is the gravitational constant. In the limit $G{M}_{NS}/\left(r{c}^2\right)\gg {\beta }^2$ , we obtain a distribution of axions $\rho ~1/r^2$ . This is assuming a sufficient emissivity of axions from the neutron stars. Thus, we may convolve the function $P\left(r\right)$ with the function $1/r^2$ , which describes the spatial density distribution of axions, to obtain $f\left(r;\gamma \right)$ . This is justified because the probability of being found at a distance $r$ ($~1/r^2$ ) is mutually exclusive of the probability of survival at a distance $r$ , i.e. the joint probability distribution is the convolution of these two functions.

The energy and spatial dependence of the flux may be factorized:

$\frac{\text{d}\Phi }{\text{d}E\text{d}\Omega }=\frac{1}{2\pi }\frac{\text{d}{P}_r\left(\theta ;m_a,\omega \right)}{\text{d}\cos \theta }\frac{\text{d}{\Phi }_0}{\text{d}E}.$ (12)

Figure 2. The geometry of the axion decays into photon. The NS-LAT line defines the focal plane. The axion is emitted on a radial path with a colatitude of ${\theta }^{\prime }$ . The decay photon is emitted at an angle ${\theta }_a$ , and the $\theta$ is a measure of the extension of the source.

The geometry of the decays with respect to the neutron star and the LAT is shown in Figure 2. The distribution of the opening angles of the photons, ${\theta }_a$ , with respect to the axion momentum direction, is given by [36]:

$P\left({\theta }_a\right)=\frac{1}{4\beta \gamma }\frac{\cos \frac{{\theta }_a}{2}}{{\sin }^2\frac{{\theta }_a}{2}}\frac{1}{\sqrt{{\gamma }^2{\sin }^2\frac{{\theta }_a}{2}-1}}.$ (13)

This distribution is strongly peaked in the forward (${\theta }_a=0$ ) direction, and has a characteristic width $1/\gamma$ . We may determine the spatial distributions of γ-rays using the following procedure. We sample a $\gamma$ parameter from the distribution $g\left(\omega \right)={\omega }^4{S}_{\sigma }\left({\omega }^4\right)$ . The radial coordinate of decay is sampled randomly from $f\left(r;\gamma \right)$ . We can determine $\theta$ simply from geometrical considerations according to:

$\theta =\text{arcsin}\left(\frac{r}{d}\sin \left(\pi -{\theta }_a\right)\right).$ (14)

Geometrically, the condition for acceptance of the photon event is:

$\left|\theta -\left(\pi /2-{\theta }_a-{\theta }_p\right)\right|<0.4\text{rad}$ (15)

This corresponds to a condition on the polar angle ${\theta }_p\lesssim {23.0}^{\circ }$ , where ${\theta }_p$ is the angle between the normal vector $\widehat{n}$ to the top of the LAT and the vector of the photon momentum ${\widehat{p}}_{\gamma }$ .

This condition on the polar angle is derived as follows: define a colatitude angle to ${\theta }^{\prime }$ , ${\theta }_1$ . From Figure 2, it follows that:

${\theta }_1={\theta }_a-\theta ,$ (16)

and

${\theta }^{\prime }=\pi /2-{\theta }_a+\theta ,$ (17)

Considering that the field of view (FOV := $\text{Δ}{\Omega }_{\text{LAT}}$ ) of the LAT is 2.4 sr [4], where:

$\text{Δ}{\Omega }_{\text{LAT}}={\displaystyle {\int }_{\text{LAT}}\text{d}\varphi }{\displaystyle {\int }_{\text{LAT}}\text{d}\left(\text{cos}\theta \right)}$ (18)

$\text{Δ}{\Omega }_{\text{LAT}}=2.4=2\pi {\displaystyle \int \text{d}\left(\text{cos}\theta \right)}$ (19)

Thus, we derive the condition of the angular acceptance of the photon events in the LAT as:

$\Delta \cos \theta \simeq \mathrm{0.4.}$ (20)

Figure 3. Angular profile of 1 meV and 10 meV axions: theoretical distribution (blue) and theoretical distribution convolved with the Pass8 LAT PSF at 60 MeV (green).

In Figure 3, we plot the angular distributions for $P\left(\theta \right)$ and $P_r\left(\theta \right)$ , which were convolved numerically, for $m_a=1$ meV. The Monte Carlo simulation was carried out with 10⁹ events. The distribution of $\theta$ was determined from sampling ${\theta }_a$ according to Equation (14). From this, $\theta$ was determined from Equation (15). The theoretical distribution derived from Monte Carlo simulation can be convolved with the Point Spread Function (PSF) of the LAT at 60 MeV, which is approximately 6˚ for events that convert in front of the tracker. The angular spread of the theoretical and convolved probability distributions plotted in Figure 3 may be parametrized by the following quantity:

$\sqrt{\langle {\theta }^2\rangle }=\sqrt{\frac{{\displaystyle {\sum }_i{P}_i{\theta }_i^2}}{{\displaystyle {\sum }_i{P}_i}}},$ (21)

where $P_i$ is the probability per $i^{th}$ bin, and ${\theta }_i^2$ is the squared value of $\theta$ in that bin. We assume, of course, that the distribution is assumed symmetric about 0.

From inspection of Figure 3, for 1 meV axions, $\sqrt{\langle {\theta }^2\rangle }={1.8}^{\circ }$ , and for 10 meV axions, $\sqrt{\langle {\theta }^2\rangle }={2.46}^{\circ }$ . Observe that the distribution for the 1 meV axions is diminished, but non-negative, between 2˚ and 4˚, due in part to the convolution kernel of the LAT PSF.

We may determine the γ-ray energy from the following equation:

$E_{\gamma }=\frac{1}{2}m\gamma \left(1+\beta \cos \left({\theta }_a\right)\right)$ (22)

We determine the spectral energy distribution from modifying Equation (6), by considering instead of the distance $d$ of the neutron star, the distance $r^{\prime }$ from the LAT at which the decay vertex $a\to 2\gamma$ occurred, which is given by:

$r^{\prime }=\sqrt{r^2+d^2-2rd\sin {\theta }^{\prime }}$ (23)

Thus, we obtain:

$\begin{array}{c}E\frac{\text{d}\Phi }{\text{d}E}=1.8\times {10}^{-2}{\left(\frac{m_a}{\text{eV}}\right)}^5\left(\frac{\text{Δ}t}{23.2\text{s}}\right){\left(\frac{100\text{pc}}{r^{\prime }}\right)}^2\\ \times {\left(\frac{2E}{100\text{MeV}}\right)}^4\left(\frac{S_{\sigma }\left(2E\right)}{{10}^7{\text{MeV}}^2}\right){\text{cm}}^{-2}\cdot {\text{s}}^{-1}.\end{array}$ (24)

We may note the dependence on the fifth power of the axion mass, as was derived in Ref. [19]. The angular probability distribution $P\left(\theta \right)$ falls off rapidly with increasing angle. The smaller the axion mass $m_a$ , the narrower the distribution $P\left(\theta \right)$ . It may be observed that the larger the γ-ray energy, the narrower the angular distribution, as shown in Figure 3. Our limits will provide larger values of $f_a$ for smaller values of $m_a$ that we choose.

4. Projected Limits from the Extended Model for RX J1856-3754

4.1. Simulation Experiment

We attempt a simulation experiment in order to test the feasibility of determining a signal from a given simulated flux. One simulation model was considered: the extended model presented earlier. This simulated model was generated using gtobssim from the ScienceTools, as described in Ref. [37]. This was done using the energy-dependent spatial templates described earlier in this paper in Section 2, while the spectral model was generated using the function in Equation (24). The same instrument response functions were used for the data analysis from the experimental observations; in the second case, the extended model for axions was considered.

For various injected values of the axion mass $m_a$ into a Monte Carlo simulation, the Test Statistic (TS) and the test statistic for extension (TS_ext) have been tabulated. The experimental values, for $m_a\simeq 1$ meV, the TS for a point source is 13.56, and the TS_ext = 27.11. From simulations of the ROI, for a 1 meV axion, the TS = 0.015, and TS_ext = 0.031. These values are computed from the likelihood function $ℒ$ , where $ℒ$ is the likelihood function. We refer here and henceforth to the log-likelihood test statistic $\text{TS}\propto \log ℒ$ . The customary interpretation of TS in sigma ($\sigma$ ) is $\sigma \simeq \sqrt{TS}$ . Thus, we establish that the source may be extended . However, the full simulation of the ROI corresponding to J0108 using gtobssim doesn’t match, probably because the fluxes of the point sources are not optimized for this ROI. The value of the Test Statistic (TS) for extension is for an axion mass of 0.34 meV and flux of 1.84 × 10⁷ MeV∙cm⁻²∙s⁻¹, where ${ℒ}_0$ is the value of the likelihood function, and ${ℒ}_{ext}$ is the value of the likelihood function for the extended model. This value of TS_ext is marginally significant, and signals that if extension of such a source exists, then it would be feasible to quantify this numerically. The TS for detection of an actual source would be 89 for a putative axion mass of 1 meV, which corresponds to a $\sigma =9.4$ . The flux determined from the simulated extended model is 330 × 10⁻⁶ cm⁻²∙s⁻¹∙MeV⁻¹, over a range of gamma-ray energies. This energy flux compares well with the expected SED in Figure 6.

4.2. Simulation of Spectral Model and Spatial Template

We chose a near galactic neutron star in order to consider the optimum sensitivity possible with this model.

RX J1856-3754 has a surface temperature of 40 eV and an age of 3.8 × 10⁶ yr [38]. Assuming a modest temperature gradient for this medium-aged pulsar, we may obtain a core temperature of 10 MeV, which may be a reasonable depending on the inputs to this model. This star may also be a quark star, which may help account for the high core temperature [38] [39].

The neutron star RX J1856-3754 was chosen because it is one of the nearest neutron stars at a distance of 160 pc [40], and it lies greater than 15˚ degrees above the galactic plane. This is justified on the grounds that the possible background contamination to the putative signal, the galactic diffuse emission, is greater near the galactic plane. From the point of view of future Fermi-LAT observations, nearby gamma-ray point sources are less than 1.5˚ degrees away from RX J1856-3754 [41] in celestial coordinates.

The spectral model corresponding to Equation (24) is plotted in Figure 4. Gamma-rays with energy 60 - 200 MeV are produced from this model. The spectral model peaks around 75 MeV, and is convolved with the point-spread function of energy.

Observe that the spatial distribution around the images in Figure 5 has energy-dependent radii. It is all the more important to consider this when optimizing an astrophysical analysis.

Figure 4. PSF-convolved spectrum for extended emission of 10 meV axions from neutron star RX J1856-3754, see text for details.

Figure 5. Spatial distribution maps of gamma-rays around neutron star J0108-1431. E = 36.27, 53.00, 77.46, 113.20, 165.44. The units of the colorbar are density (normalized to 1) on a logarithmic scale. The pixels are 0.4˚ × 0.4˚.

Figure 6. The expected spectral energy distribution of 10 meV axions from the spatially extended model compared to the point source sensitivity of the LAT. The region is determined by 10% hardening or softening of the spectrum curve. The LAT curve is the 10-year sensitivity using Fermi-LAT Pass 8 instrument response functions, with a region corresponding to 95% C.L.

In Figure 6, we plot the expected sensitivity of the model. With the projected 10-year sensitivity of Fermi-LAT using Pass 8 instrument response functions, it may be possible in principle to observe or set upper limits on the axion mass below 10 meV.

In Table 2, the upper limit on flux for comparable sources, is shown. From this value, the upper limits on flux and axion-nucleon-nucleon coupling, are shown.

Table 2. For neutron star RX J1856-3754, the flux upper limit and the upper limit on the axion mass.

Quantity	95% C.L.
$\text{Φ}$	<1 × 10−12 cm−2∙s−1
$m_a$	< 3.9 meV
$f_a$	>6 × 109 GeV
$g_{ann}$	<3.9 × 10−11

5. Discussion

The model-dependent observational limits derived here for the hypothesis of extension of $m_a<3.9$ meV at the 95% C.L. are a substantial improvement upon the previous Fermi-LAT point-source limit with neutron stars [19]. The point-source limits for this source under the same data analysis conditions as for the extended source model would not be changed much from the previous point-source limits, as the energies where the spectrum has a considerable contribution (i.e. above 60 MeV) were considered in the previous analysis. Furthermore, the data-dependent limits also fall below the projected limits derived here. From consideration of Figure 7, the limits from the Fermi-LAT using 7.9 years of observations are improved by a factor of nearly 100 with respect to the point source limits. Also, these limits fall in a range that has not yet been excluded by previous observations, and represent a substantial improvement over the SN 1987A-derived upper limit of $m_a<16$ meV [30]. If the hypothesis of extension is valid, however low the flux may be, then we report a highly significant detection. This should be explored in a possible follow-up study. The signal for RX J1856-3754 could be contaminated with known nearby point sources, according to the 3FGL Fermi-LAT catalog. It should be noted that the possible detection of a signal depends on the extended model energy-dependent spatial templates shown and derived earlier in this study, and that the point source study of a putative signal is probably too naïve for searching for axions. Needless to say, future studies should focus on neutron stars with farther distances and that are also younger in age.

Fundamentally, we would expect to obtain better limits in this analysis using an extended source model than modeling neutron stars as axion point sources. Namely, Fermi-LAT would subtend a larger solid angle for an extended source than a point source. If no signal is observed with a larger solid angle, the upper limit on the mass $m_a$ would be lower.

It deserves to be mentioned that assuming somewhat higher or lower temperatures could alter this limit somewhat. The simulation-derived limit of 10 meV is not at the minimum range possible, according to this model, but suggests an upper bound on values that could be derived in principle.

Figure 7. Comparison of exclusion ranges compared with the possible range of masses presented in this article. Exclusion regions for axions: the Fermi-LAT point source limits (light red) from neutron stars, projected limits from the Fermi-LAT using the spatially extended model (dark red) of the neutron star RX J1856-3754, compared with previous astrophysical limits.

Interestingly, the limit constrains the allowed parameter space for axions as cold dark matter. These projected limits are better than the bounds of $m_a<20$ meV reported by CAST [42]. While the ADMX projection excludes a smaller range of masses, it is probing DFSZ axions, not KSVZ axions. This model sets more restrictive limits than this and other neutron star cooling observations [16]. Neutron star cooling by axions is a quite distinct process from the process of emitted axions decaying outside the neutron star, for which current limits have been reported as $m_a<60$ meV. In future observations, the projected limits could potentially be improved by statistically combining limits from multiple neutron stars, as shown in Ref. [19]. Although there may be some uncertainty over the precise temperature of the neutron star, which we assume as $T=20$ MeV, this applies generally to neutron stars with hadronic physics. In extended models of neutron stars, which contain free quarks in a QCD phase [43], it is generally accepted that there is a range of temperatures that are higher, generally between 10 MeV to 60 MeV [33]. In this article, we set limits on $f_a~{10}^{10}$ GeV. Axions with masses ~1 meV may be a source of dark matter, although it cannot comprise the majority of the dark matter because the abundance is too low [44]. This bound is much stronger than the weak upper bound of $f_a\lesssim {10}^{12}$ GeV from cosmological arguments [45]. In future work, a consideration of the possible background signal from the pulsar in addition to the axion signal could result in a better limit on the axion mass. In addition, simulation studies could enhance the selection of neutron star targets that would yield the best limits. If Equation (1) is relaxed, so that both $m_a$ and $f_a$ are free parameters, then ALPs may also be considered in a generalization of the model presented here. Qualitatively, increasing (decreasing) the assumed $T$ would tend to shift the spectrum towards higher (lower) energies. The model flux depends on ${\omega }^4{S}_{\sigma }\left(\omega ;\mu ,T\right)$ , which increases with $T$ , but the timescale depends on ${\left({\displaystyle \int \text{d}\omega }\omega S_{\sigma }\left(\omega ;\mu ,T\right)\right)}^{-1}$ , which decreases with $T$ . Thus, a simple calculation finds that the limits on $m_a$ would be smaller for $T=50$ MeV, and larger for $T=10$ MeV. The order of magnitude of the limits would still be the same for these changes in temperature. If the assumed $T$ was higher, the Fermi-LAT could be used to set more stringent limits. For $T=10$ MeV, we estimate the limits would be a factor of 4 less restrictive. Increasing the degeneracy parameter $\mu$ would tend to decrease the amplitude of the spin-structure function. At $\mu /T=11$ , the limits would be larger, and at $\mu /T=9$ , the limits would be smaller. Changing the $k$ parameter would not affect the limits substantially.

Similar limits have been achieved by other recent axion searches. Supernova constraints for SN 1987A yield a limit of $m_a\lesssim 16$ meV [30]. From vacuum magnetic birefringence, a limit of $m_a<5\times {10}^{-4}$ eV has been reached [46]. From magnetically induced dichroism, $m_a=1\text{-}1.5$ meV has been reached [47]. From globular cluster studies, a limit of $m_a{\sin }^2\beta <15$ meV has been attained [48]. From the white dwarf luminosity function, an inclusion range of $2.5\lesssim m_a{\sin }^2\beta \lesssim 7.5$ meV has been achieved [49] [50]. From the supernova remnant Cas A, the axion mass was set at $m_a=2.3\pm 0.4$ meV/C_n [51]. These limits are complementary to the results shown in our work, and we hope that better limits can be obtained in future work. The work here assumes superfluidity in the neutron star, which does apply. This allows for higher temperatures in the phase diagram than is the case with conventional models for neutron stars.

Acknowledgements

BB acknowledges support from California State University, Los Angeles Department of Physics and Astronomy in the College of Natural and Social Sciences. BB acknowledges support from a generous grant from the Orthopedic and Arthritic Center, Oxnard, California.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Peccei, R.D. and Quinn, H.R. (1977) CP Conservation in the Presence of Pseudoparticles. Physical Review Letters, 38, 1440-1443. https://doi.org/10.1103/physrevlett.38.1440
[2]	Sikivie, P. (2011) The Emerging Case for Axion Dark Matter. Physics Letters B, 695, 22-25. https://doi.org/10.1016/j.physletb.2010.11.027
[3]	Skobelev, V.V. (2000) Primakoff Effect: Synchrotron and Coulomb Mechanisms of Axion Emission. Physics of Atomic Nuclei, 63, 1963-1968. https://doi.org/10.1134/1.1335094
[4]	Atwood, W.B., et al. (2009) The Large Area Telescope on the Fermi Gamma-Ray Space Tele-Scope Mission. The Astrophysical Journal, 697, Article 1071. http://stacks.iop.org/0004-637X/697/i=2/a=1071
[5]	Lande, J., Ackermann, M., Allafort, A., Ballet, J., Bechtol, K., Burnett, T.H., et al. (2012) Search for Spatially Extended FERMI Large Area Telescope Sources Using Two Years of Data. The Astrophysical Journal, 756, Article 5. https://doi.org/10.1088/0004-637x/756/1/5
[6]	Ackermann, M., et al. (2017) Search for Extended Sources in the Galactic Plane Using Six Years of Fermi-Large Area Telescope Pass 8 Data above 10 GeV. The Astrophysical Journal, 843, 139-163.
[7]	Ackermann, M., et al. (2018) The Search for Spatial Extension in High-Latitude Sources Detected by the Fermi Large Area Telescope. The Astrophysical Journal Supplement Series, 237, 32-68.
[8]	Abazajian, K.N. and Kaplinghat, M. (2012) Detection of a Gamma-Ray Source in the Galactic Center Consistent with Extended Emission from Dark Matter Annihilation and Concentrated Astrophysical Emission. Physical Review D, 86, Article ID: 083511. https://doi.org/10.1103/physrevd.86.083511
[9]	Ackermann, M., Ajello, M., Albert, A., Baldini, L., Ballet, J., Barbiellini, G., et al. (2017) Observations of M31 and M33 with the Fermi Large Area Telescope: A Galactic Center Excess in Andromeda? The Astrophysical Journal, 836, Article 208. https://doi.org/10.3847/1538-4357/aa5c3d
[10]	Giannotti, M., Duffy, L.D. and Nita, R. (2011) New Constraints for Heavy Axion-Like Particles from Supernovae. Journal of Cosmology and Astroparticle Physics, 2011, Article 15. https://doi.org/10.1088/1475-7516/2011/01/015
[11]	Meyer, M., Horns, D. and Raue, M. (2013) First Lower Limits on the Photon-Axion-Like Particle Coupling from Very High Energy γ-Ray Observations. Physical Review D, 87, Article ID: 035027. https://doi.org/10.1103/physrevd.87.035027
[12]	Raffelt, G. and Stodolsky, L. (1988) Mixing of the Photon with Low-Mass Particles. Physical Review D, 37, 1237-1249. https://doi.org/10.1103/physrevd.37.1237
[13]	Raffelt, G.G. (1996) Stars as Laboratories for Fundamental Physics: The Astrophysics of Neutrinos, Axions, and Other Weakly Interacting Particles. University of Chicago Press.
[14]	Raffelt, G.G., Redondo, J. and Maira, N.V. (2011) The Mev Mass Frontier of Axion Physics. Physical Review D, 84, Article ID: 103008. https://doi.org/10.1103/physrevd.84.103008
[15]	Redondo, J., Raffelt, G. and Maira, N.V. (2012) Journey at the Axion Mev Mass Frontier. Journal of Physics: Conference Series, 375, Article ID: 022004. https://doi.org/10.1088/1742-6596/375/1/022004
[16]	Sedrakian, A. (2016) Axion Cooling of Neutron Stars. Physical Review D, 93, Article ID: 065044. https://doi.org/10.1103/physrevd.93.065044
[17]	Atwood, W., et al. (2012) Pass 8: Toward the Full Realization of the Fermi-LAT Scientific Potential. arXiv: astro-ph/1303.3514.
[18]	Lloyd, S.J., Chadwick, P.M. and Brown, A.M. (2019) Constraining the Axion Mass through Gamma-Ray Observations of Pulsars. Physical Review D, 100, Article ID: 063005. https://doi.org/10.1103/physrevd.100.063005
[19]	Berenji, B., Gaskins, J. and Meyer, M. (2016) Constraints on Axions and Axionlike Particles from Fermi Large Area Telescope Observations of Neutron Stars. Physical Review D, 93, Article ID: 045019. https://doi.org/10.1103/physrevd.93.045019
[20]	Hamaguchi, K., Nagata, N., Yanagi, K. and Zheng, J. (2018) Limit on the Axion Decay Constant from the Cooling Neutron Star in Cassiopeia A. Physical Review D, 98, Article ID: 103015. https://doi.org/10.1103/physrevd.98.103015
[21]	Abdelhameed, A.H., Bakhlanov, S.V., Bauer, P., Bento, A., Bertoldo, E., Canonica, L., et al. (2020) New Limits on the Resonant Absorption of Solar Axions Obtained with a 169Tm-Containing Cryogenic Detector. The European Physical Journal C, 80, Article No. 376. https://doi.org/10.1140/epjc/s10052-020-7943-5
[22]	Kim, J.E. (1979) Weak-Interaction Singlet and Strongcpinvariance. Physical Review Letters, 43, 103-107. https://doi.org/10.1103/physrevlett.43.103
[23]	Shifman, M.A., Vainshtein, A.I. and Zakharov, V.I. (1980) Can Confinement Ensure Natural CP Invariance of Strong Interactions? Nuclear Physics B, 166, 493-506. https://doi.org/10.1016/0550-3213(80)90209-6
[24]	Dine, M., Fischler, W. and Srednicki, M. (1981) A Simple Solution to the Strong CP Problem with a Harmless Axion. Physics Letters B, 104, 199-202. https://doi.org/10.1016/0370-2693(81)90590-6
[25]	Zhitnitsky, A. (1980) On Possible Suppression of the Axion Hadron Interactions. Soviet Journal of Nuclear Physics, 31, 260.
[26]	Tanabashi, M., et al. (2018) The Review of Particle Physics. Physical Review D, 98, Article ID: 030001.
[27]	Ratra, B. (1991) Expressions for Linearized Perturbations in a Massive-Scalar-Field-Dominated Cosmological Model. Physical Review D, 44, 352-364. https://doi.org/10.1103/physrevd.44.352
[28]	Sánchez-Conde, M.A., Paneque, D., Bloom, E., Prada, F. and Domínguez, A. (2009) Hints of the Existence of Axionlike Particles from the γ-Ray Spectra of Cosmological Sources. Physical Review D, 79, Article ID: 123511. https://doi.org/10.1103/physrevd.79.123511
[29]	Blaschke, D., Klähn, T. and Sandin, F. (2007) Equation of State at High Densities and Modern Compact Star Observations. Journal of Physics G: Nuclear and Particle Physics, 35, Article ID: 014051. https://doi.org/10.1088/0954-3899/35/1/014051
[30]	Raffelt, G.G. (2008) Astrophysical Axion Bounds. In: Kuster, M., Raffelt, G. and Beltrán, B., Eds., Axions, Springer, 51-71. https://doi.org/10.1007/978-3-540-73518-2_3
[31]	Hanhart, C., Phillips, D.R. and Reddy, S. (2001) Neutrino and Axion Emissivities of Neutron Stars from Nucleon-Nucleon Scattering Data. Physics Letters B, 499, 9-15. https://doi.org/10.1016/s0370-2693(00)01382-4
[32]	Raffelt, G.G. (1990) Astrophysical Methods to Constrain Axions and Other Novel Particle Phenomena. Physics Reports, 198, 1-113. https://doi.org/10.1016/0370-1573(90)90054-6
[33]	Rüster, S.B., Werth, V., Buballa, M., Shovkovy, I.A. and Rischke, D.H. (2005) Phase Diagram of Neutral Quark Matter: Self-Consistent Treatment of Quark Masses. Physical Review D, 72, Article ID: 034004. https://doi.org/10.1103/physrevd.72.034004
[34]	Pavlov, G.G., Kargaltsev, O., Wong, J.A. and Garmire, G.P. (2009) Detection of X-Ray Emission from the Very Old Pulsar J0108-1431. The Astrophysical Journal, 691, 458-464. https://doi.org/10.1088/0004-637x/691/1/458
[35]	van Riper, K.A., Link, B. and Epstein, R.I. (1995) Frictional Heating and Neutron Star Thermal Evolution. The Astrophysical Journal, 448, Article 294. https://doi.org/10.1086/175961
[36]	McDonald, K.T. (2019) Neutral Pion Decay. http://www.physics.princeton.edu/mcdon-ald/examples/piondecay.pdf
[37]	(2015) FSSC: Fermi Data: Online Documentation: Science Tools: Cicerone. https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/
[38]	Drake, J.J., Marshall, H.L., Dreizler, S., Freeman, P.E., Fruscione, A., Juda, M., et al. (2002) Is RX J1856.5-3754 a Quark Star? The Astrophysical Journal, 572, 996-1001. https://doi.org/10.1086/340368
[39]	Trümper, J.E., Burwitz, V., Haberl, F. and Zavlin, V.E. (2004) The Puzzles of RX J1856.5-3754: Neutron Star or Quark Star? Nuclear Physics B—Proceedings Supplements, 132, 560-565. https://doi.org/10.1016/j.nuclphysbps.2004.04.094
[40]	Posselt, B., Arumugasamy, P., Pavlov, G.G., Manchester, R.N., Shannon, R.M. and Kargaltsev, O. (2012) XMM-Newton Observation of the Very Old Pulsar J0108-1431. The Astrophysical Journal, 761, Article 117. https://doi.org/10.1088/0004-637x/761/2/117
[41]	Acero, F., Ackermann, M., Ajello, M., Albert, A., Atwood, W., Axelsson, M., Baldini, L., Ballet, J., Barbiellini, G., Bastieri, D., et al. (2015) Fermi Large Area Telescope Third Source Catalog. The Astrophysical Journal Supplement Series, 218, 23-64.
[42]	Arik, M., Aune, S., Barth, K., Belov, A., Bräuninger, H., Bremer, J., et al. (2015) New Solar Axion Search Using the CERN Axion Solar Telescope with 4He Filling. Physical Review D, 92, Article ID: 021101. https://doi.org/10.1103/physrevd.92.021101
[43]	Alford, M.G. (2009) Quark Matter in Neutron Stars. Nuclear Physics A, 830, 385c-392c. https://doi.org/10.1016/j.nuclphysa.2009.09.034
[44]	Graham, P.W. and Rajendran, S. (2013) New Observables for Direct Detection of Axion Dark Matter. Physical Review D, 88, Article ID: 035023. https://doi.org/10.1103/physrevd.88.035023
[45]	Dine, M. and Fischler, W. (1983) The Not-So-Harmless Axion. Physics Letters B, 120, 137-141. https://doi.org/10.1016/0370-2693(83)90639-1
[46]	Maiani, L., Petronzio, R. and Zavattini, E. (1986) Effects of Nearly Massless, Spin-Zero Particles on Light Propagation in a Magnetic Field. Physics Letters B, 175, 359-363. https://doi.org/10.1016/0370-2693(86)90869-5
[47]	Zavattini, E., Zavattini, G., Ruoso, G., Polacco, E., Milotti, E., Karuza, M., et al. (2006) Experimental Observation of Optical Rotation Generated in Vacuum by a Magnetic Field. Physical Review Letters, 96, Article ID: 110406. https://doi.org/10.1103/physrevlett.96.110406
[48]	Viaux, N., Catelan, M., Stetson, P.B., Raffelt, G.G., Redondo, J., Valcarce, A.A.R., et al. (2013) Neutrino and Axion Bounds from the Globular Cluster M5 (NGC 5904). Physical Review Letters, 111, Article ID: 231301. https://doi.org/10.1103/physrevlett.111.231301
[49]	Bertolami, M.M.M., Melendez, B.E., Althaus, L.G. and Isern, J. (2014) Revisiting the Axion Bounds from the Galactic White Dwarf Luminosity Function. Journal of Cosmology and Astroparticle Physics, No. 2014, Article 69. https://doi.org/10.1088/1475-7516/2014/10/069
[50]	Isern, J., García-Berro, E., Torres, S. and Catalán, S. (2008) Axions and the Cooling of White Dwarf Stars. The Astrophysical Journal, 682, L109-L112. https://doi.org/10.1086/591042
[51]	Leinson, L.B. (2014) Axion Mass Limit from Observations of the Neutron Star in Cassiopeia A. Journal of Cosmology and Astroparticle Physics, No. 2014, Article 31. https://doi.org/10.1088/1475-7516/2014/08/031