1 Introduction↩︎

The Standard Model (SM) of particle physics has achieved remarkable experimental success over the past several decades, yet several unresolved issues remain, among which the strong CP problem in quantum chromodynamics (QCD) is particularly compelling[1]–[3]. This problem arises from the theoretically allowed but experimentally unobserved CP-violating term in the QCD Lagrangian: \[\mathcal{L}_{\text{CP-vio}}=\frac{\alpha_{s}}{4 \pi} \theta \operatorname{Tr}\left(G_{\mu v} \tilde{G}^{\mu v}\right)\] where \(G_{\mu\nu}\) is the gluon field strength tensor and \(\tilde{G}^{\mu\nu}\) its dual. To resolve this, the Peccei-Quinn mechanism was proposed, giving rise to the axion—a pseudoscalar Nambu-Goldstone boson. Extensions of the SM predict more general axion-like particles (ALPs) that couple to photons via \[\mathscr{L}_{a \gamma \gamma} = -\frac{1}{4} g_{a \gamma \gamma} a F_{\mu \nu} \tilde{F}^{\mu \nu} = g_{a \gamma \gamma} a \boldsymbol{{E}} \cdot \boldsymbol{B},\] where \(a\) is the axion field, \(g_{a\gamma\gamma}\) is the coupling constant, \(F\) and \(\tilde{F}\) are the electromagnetic field-strength tensor and its dual, E and B are the electric field and magnetic field. For an electromagnetic wave, the interaction specifically requires the electric field component parallel to the external magnetic field (\(E_\parallel B_\perp\)). This property enables photon-to-ALPs conversion in strong magnetic fields; the resulting photon-ALP conversion probability in a magnetic field is determined by the mixing of the photon and ALP states as they propagate through the magnetosphere [4], [5].

ALPs are probed by experiments targeting either non-relativistic or relativistic regimes. Non-relativistic ALPs, with low velocities are primary dark matter candidates, detected by haloscopes like Axion Dark Matter eXperiment (ADMX) and Taiwan Axion Search Experiment with a Haloscope (TASEH) [6]–[10] in resonant cavity as well as constrained by cosmological evolution like Big Bang Nucleosynthesis (BBN) [11] and cosmic background [12] and so on.Relativistic ALPs are explored by helioscopes like CERN Axion Solar Telescope (CAST) [13], by light-shining-through-walls (LSW) experiments like Any Light Particle Search (ALPS)[14], and by pulsarscope studies [15].

Pulsars are uniquely powerful source for axion-like particle research due to their extreme magnetic fields. Young rotation-powered pulsars typically have surface dipole magnetic fields of \(10^{11}\)–\(10^{13}\) G, far stronger than any achievable in terrestrial experiments and significantly higher than those of millisecond pulsars (MLPs). Such intense fields greatly enhance the probability of photon–ALP conversion, which may result to fluctuation features in their X-ray spectra.

In this work we probe possible fluctuation in pulsar spectra arising from photon-ALPs conversion in their magnetospheres. Focusing on the soft X-Ray band accessible to NASA’s Neutron Star Interior Composition Explorer (NICER, 0.3-12 keV), we exclude the non-relativistic dark-matter ALP regime and instead explore the relativistic regime. This enables us to improve the constraints on coupling constant \({g}_{a \gamma \gamma}\) within a refined parameter space compared to previous results[16]–[18].

2 Formula↩︎

The evolution of the coupled photon-ALP system as it traverses the pulsar magnetosphere can be described by a Schrödinger-like equation:\[-i\frac{d}{dz}\Psi(z) = \mathcal{M}\Psi(z), \quad \Psi(z) = \begin{pmatrix} A_{\perp}(z) \\ a(z) \end{pmatrix}\] where \(A_{\perp}(z)\) and \(a(z)\) represent the photon and ALP field amplitudes with path length \(z\), respectively. The mixing matrix \(\mathcal{M}\) accounts for the interaction between the states: \[\mathcal{M} = = \omega\mathbb{1} +\begin{pmatrix} \Delta_{\gamma}(z) & \Delta_{a\gamma}(z) \\ \Delta_{a\gamma}(z) & \Delta_{a}(z) \end{pmatrix}\] In this framework, \(\omega\) is the oscillation frequency in vacuum, \(\Delta_{\gamma}\) and \(\Delta_{a}\) are the refractive indices for photons and ALPs, and \(\Delta_{a\gamma}(z) = \frac{1}{2}g_{a\gamma\gamma}B_{\perp}(z)\) represents the interaction strength mediated by the transverse magnetic field.

To quantify the photon-to-ALP conversion, we consider the probability \(\mathcal{P}(\gamma \rightarrow a)\) of a photon emitted from the stellar surface (\(z=0\)) converting into an ALP as it propagates to infinity. Following standard magnetospheric emission models, we assume the X-ray photons originate at or near the surface, where the dense stellar interior is opaque to X-ray propagation. The momentum transfer between the photon and ALP is given by \(q = \frac{m_a^2 - \omega_{pl}^2}{2E}\) , where the \(\omega_{pl}\) is the plasma frequency. In the small-mixing limit(\(q \approx 0\)), where the conversion efficiency is low, the probability can be expressed in an integral form:

\[\begin{align} \nonumber \mathcal{P}(\gamma\rightarrow a) = & \left| \int_{0}^{\infty} \Delta_{a\gamma}(z) e^{iqz} dz \right|^{2} \\ \nonumber & = \left| \int_{0}^{\infty} \frac{1}{2}g_{a\gamma\gamma} B_{\perp}(z) e^{iqz} dz \right|^{2} \\ & \approx \left| \int_{0}^{\infty} \frac{1}{2}g_{a\gamma\gamma} B_{\perp}(z) \right|^{2} \text{if } q\approx 0\;. \end{align}\]

For a standard dipole distribution, the transverse magnetic field component is \(B_{\perp}(r) = B(r) \sin \theta\), where \(\theta\) is the angle between the magnetic axis and the line of sight. The field scales as \(1/r^3\). By substituting \(r = R_{\text{NS}} + z\), where \(R_{\text{NS}}\) is the neutron star radius, we anchor the coordinate system to the emission site at the surface (\(z=0\)). Under this approximation, the integral simplifies to a clean analytical form:

\[\begin{align} &\mathcal{P}(\gamma\rightarrow a) = \left| \frac{1}{2}g_{a\gamma\gamma} \sin\theta B_{0} R_{NS}^{3} \int_{0}^{\infty} \frac{1}{(R_{NS}+z)^{3}} dz \right|^{2} \nonumber \\ & = \left| \frac{1}{4}g_{a\gamma\gamma} \sin\theta B_{0}R_{NS} \right|^{2} \;. \end{align}\]

Beyond the dipole assumption, we can expand this to a general multipole form. Since there is no free current outside the neutron star, the magnetic potential \(\Phi_M\) satisfies Laplace’s equation, \(\nabla^2\Phi_M = 0\). In spherical coordinates, requiring the field to vanish at infinity implies a radial scaling of \(B \propto 1/r^{l+2}\) for a multipole of order \(l\): \[B_{\perp}(r,\theta) = B_{0} \left( \frac{R_{\text{NS}}}{r} \right)^{l+2} \sin\theta\] Using this generalized field in the conversion integral yields the final probability for any multipole order \(l\): \[\mathcal{P}(\gamma\rightarrow a) = \left| \frac{1}{2(l+1)} g_{a\gamma\gamma} \sin\theta B_{0} R_{\text{NS}} \right|^{2}\;.\] This generalized result allows us to evaluate the sensitivity of our constraints to different magnetospheric models, with \(l=1\) representing a standard dipole field and \(l=2\) representing the quadrupole field.

3 Analysis↩︎

We analyzed the soft X-ray spectra of three rotation-powered pulsars—PSR J2229+6114, PSR J1849-0001, and PSR B0531+21—using data from NASA’s NICER mission [19]. All spectra (see the Table 1 in Appendix) in this work were processed using the HEASoft version 6.35.2 and NICERDAS version 11 pipeline, with event filtering performed via NICERL2 and spectral extraction through NICERL3-spect. For our background subtraction, we utilized the "3c50" empirical model[20]. All of the spectra show in Fig.  1, the yellow banding region means the fit range we used. Following the magnetospheric emission models for young pulsars [21], [22], we assume the X-ray photons in the 0.8–10 keV range originate from the polar cap or the near-surface vacuum gaps. Because these emission heights are negligible compared to the conversion path length, and given the star’s opacity to X-rays, we define the starting point of our integral at the neutron star surface.

To identify potential spectral fluctuations indicative of photon-ALP conversion, each pulsar’s continuum was modeled with an absorbed power-law function: \(F(E) = AE^{-\Gamma}\). A sliding-window fitting approach was applied to detect localized deviations from this smooth baseline. We used a window width of 4 bins for energies below 2 keV and increased this to 10 bins for higher energies to account for decreasing photon statistics. Due to the overlapping nature of the windows, each energy bin is included in multiple fits and thus associated with more than one set of best-fit parameters. To determine the most representative local fit, we adopted the minimum \(\chi^2\) criterion, selecting for each energy bin the fit corresponding to the lowest reduced \(\chi^2\) value among all overlapping windows, all calculation are calculated at the 95% confidence level. The statistical significance of deviations between the data and our baseline model was quantified using the "pull" (\(z_i\)): \[z_{i}=\frac{y_{\text{data,}i}-y_{\text{model,}i}}{\sigma_{i}s_{i}}, \quad s_{i}=\sqrt{\frac{\chi^{2}}{d.o.f}}\] where \(s_i\) is the scaling factor from the local fit. From this, we calculated a two-tailed Gaussian probability \[p_{\text{two-tailed,}i} = 2[1-\Phi(|z_i|)]\] to determine the likelihood that a deviation arose from random noise, where \(\Phi\) is the Cumulative Distribution Function, which is used to represent random fluctuation would push the measurement at least this far away from prediction. Finally, we defined an empirical conversion probability,\(P_{\text{conversion,}i} = 1 - p_{\text{two-tailed,}i}\), to represent the chance that an observed fluctuation is not consistent with random noise and could potentially be interpreted as a photon-axion conversion signature.

4 Result↩︎

For a typical pulsar radius between 10 to 13 km, we take \(R_{\text{NS}} = 13\) km when estimating the conversion scale. To refine our constraints, we incorporated pulsar-specific geometry by assuming tilt angles (\(\theta\)) based on existing multi-wavelength research. We adopted \(\theta \approx 60^\circ\) for PSR B0531+21 based on its torus-jet morphology, \(\theta \approx 45^\circ\) for PSR J2229+6114 from pulse profile fitting, and a conservative \(\theta \approx 10^\circ\) for the nearly aligned PSR J1849-0001.[23]–[26]

With different angle assumption mentioned above, the coupling constant distribution show in Fig. [fig:proce951], and the combined upper bound is also showed in Fig. [fig:proce952]. A central component of our analysis is the intense magnetic environment of these targets; the surface magnetic fields (\(B_0\)) were estimated using the standard magnetic dipole braking model[27], yielding values of \(2.03 \times 10^{12}\) G for PSR J2229+6114 , \(7.5 \times 10^{11}\) G for PSR J1849\(-\)​0001, and \(1.2 \times 10^{12}\) G for PSR B0531+21.[28]–[30]

Figure [fig:result951] shows the derived constraints for each pulsar and compares them with non-dark matter axion constraints near X-ray range and this result might provide new extension of the existing bounds. Figure [fig:result952] shows the constraints compared with dark matter axion constraints from the astrophysical bounds.

5 Conclusion↩︎

In summary, We analyzed NICER observations of three rotation-powered pulsars, PSR J2229+6114, PSR J1849\(-\)​0001 and PSR B0531+21, to search for spectral signatures of photon–ALP conversion. Using a sliding-window power-law analysis, we derived constraints of \(g_{a\gamma\gamma} \sim 10^{-10}\) \(–\) \(10^{-13}\,\mathrm{GeV^{-1}}\), by assuming the angle of pulsar tilt and the magnetic field distribution. The PSR J2229+6114 and B0531+21 providing the strongest limit based on the assumptions. Though the study gives us the rough story about the coupling constant, the further study toward these pulsar is still necessary. These findings highlight the potential of X-ray observations of strongly magnetized pulsars as sensitive probes of axion-like particles in the relativistic regime.

6 Acknowledgment↩︎

We thank the Institute of Physics, Academia Sinica, and the Department of Physics, National Cheng Kung University, for their support. We also acknowledge funding from the National Science and Technology Council (NSTC), Taiwan.

7 Summary of Observations↩︎

References↩︎