Superdense beaming of axion dark matter in the vicinity of the light cylinder of pulsars

In this article we treat the non-adiabatic photon-to-axion resonant conversion of curvature radiation, synchrotron emission and inverse Compton scattering dominating the spectral density function of pulsars. First, we introduce emission models and benchmark observational data.vWe adopt a state-of-the-art density profile that relieves tension with the quantum electrodynamics vacuum polarization effect in highly magnetic stars, leading to efficient mixing. Then, we estimate the dark matter flux induced by photon-axion oscillation across the light cylinder of the neutron star. We find that pulsars might produce axion overdensities many orders of magnitude over the occupation number of dark matter in the Galactic halo within a broad parameter space. We point out possible new methods for axion detection derived from these results and other future lines of work.


Introduction
Axions [1,2] are hypothetical pseudo-scalar bosons theorized as a consequence of the dynamic solution to the problem of charge and parity (CP) in the strong interaction [3]. Furthermore, axion is a well-grounded candidate for cold dark matter (cDM) [4][5][6][7]. Axion and axion-like particles (ALPs) mix with photons by action of an external magnetic field. The axion-photon coupling Lagrangian density is being F µν the field strength tensor andF its dual, φ the axion field, m φ the pseudo-scalar mass and g φγ the axion-photon coupling constant, j is density of current and the standard photon field is A µ . Classically, the interaction term in Eq. 1.1 simplifies to g φγ φ E · B, with E the photon field and B the external magnetic field. The quantum chromodynamics (QCD) axion presents a light mass which scales inversely to a typical energy scale, f φ . This parameter determines the axion-photon coupling strength. Differently, ALPs present mass-coupling independent parameters, creating a large parameter space which has been only superficially explored -cf. [4].
In the present work we investigate the non-adiabatic resonant conversion of curvature radiation, synchrotron emission and inverse Compton scattering photons dominating the spectral density function of isolated neutron stars. We apply this concept to spectral functions from astronomical observations in order to estimate the overdensity of axion dark matter that might be beamed from somewhere around the light cylinder of rotation-powered pulsars and B-powered magnetars.
The manuscript is structured as follows. First, we briefly explain the fundamentals of neutron stars relevant for this work in section 2. The axion-photon mixing equations are introduced in section 3 in order to obtain an analytical expression which allows us to estimate the photon-to-axion conversion probability in the surroundings of the star, and we apply this formalism to realistic spectral density functions. The results are presented in section 4, and discussed in section 5.

Rotation-powered pulsars and magnetars
Pulsating stars (pulsars) are rotating neutron stars with regular periods. 'Canonical' neutron stars would present a radius about R ∼10 km, a mass of the order of M ∼1.4 M . In a pseudo-aligned rotator, the field becomes stationary, or B z B 0 R r 3 , and the dielectric tensor does not vary significantly during a period. Here, B 0 is the field strength at the surface of the NS, of the order of 10 8−10 T, or 10 10−11 T in the case of magnetars [29][30][31]; r is the radial coordinate. The angular velocity is Ω = 2π/P , being P a rotation period from tens to milliseconds. The surface temperature is around 10 6 K. Their thin atmosphere is composed of a 'cold' plasma -mainly H and He - [32]. A neutron star is illustrated in Fig. 1, where some aspects that are relevant to this article are highlighted.  Even if discovered half a century ago [33], pulsar emission mechanisms remain full of incertitude [34]. The standard picture is based on the creation of pair cascades and particle migration trough magnetic forces, accelerating the charges to relativistic velocities. Radiative processes are dominated by curvature radiation, synchrotron emission and inverse Compton scattering. However, acceleration and emission regions, or the specific altitude -sometimes referred to as h 0 [35] -at which the pair cascade and charge acceleration start, is a prominent point of discrepancy. The classic picture considers radio emission from pulsars is originated at 'low altitudes' over the polar cap [36][37][38]. The high-energy emission is regularly associated to a wider interval of higher altitudes. Nevertheless, the precise concept of 'low altitude' and 'high altitude' is diffuse, and different models suggest that high-energy emission originates somewhere inside the light cylinder (LC), around the LC, etc. -cf. [8,34] and references therein.
The light cylinder is filled with plasma at a corotation density. Goldreich & Julian (GJ) calculated the minimal corotational charge densityn GJ -in a neutron star under a number of idealizations, including a perfect conductivity and homogeneity [39]. Frequently, modern approaches are over-dense compared to the GJ model, and a pair multiplicity factor κ = n e /n GJ arises with typical values 10 4−7 [40][41][42][43]. In that manner, the charge density in our model of star reads while the characteristic frequency of the plasma is Observational campaigns are a reliable alternative in order to provide luminosity curves and spectral densities in the absence of a well-established theoretical model. In Fig. 2 we present the reference spectral density functions for this article. The radio-frequency (RF) spectra is modeled using Eq. A.1 [44][45][46][47][48], or data from [49][50][51] for magnetar PSR J1550-5418 -an object of interest as its rotation and magnetic axis appear to be lined upand magnetar PSR J1809-1943 [52]. The X-ray spectra, adapted from [53,54], corresponds with magnetar PSR J0146+61. The Geminga -rotation-powered -pulsar PSR J0633+17 spectral function is adapted from [54,55].
The polarization state of synchrotron emission is governed by the magnetic vector. Even though that the polarization of the photon beam from a neutron star can be affected by Faraday rotation and other spurious effects through line-of-sight to observer, recent research could shed some light about its properties. Interestingly, the persistent signal received from neutron stars presents a strong linear polarization in the direction of the rotation axis, with Stokes V 0 [8,34,56,57]. This parallel component of the photon field refers to the ordinary mode -O-mode -that we confine our concern throughout this work.

Non-adiabatic photon-axion resonant oscillation in neutron stars
In this section we focus on O-mode waves propagating in the 'weak dispersion' limit with a refractive index close to unity, or equivalently ω∼|k|. In a strongly magnetized 'cold' plasma, the conductivity is non-zero only along the direction of the magnetic field [58] -appendix B -. The photon-axion oscillation takes place 'on-the-spot', or in a narrow projection along r c and hence the gravitational gradient is negligible. In that manner, the wavelength becomes stationary and arbitrarily shifted to the red. The mixing equations are derived from Eq. D.2b [12] i d dz g φγ ] 1/2 and the mixing angle tan(2θ m ) = 2∆ g φγ /(∆ φ − ∆ p − ∆ ). The term ∆ vanishes in the weak dispersion limit as the refractive index is n k/ω = 1 + ∆ /ω. Quantum electrodynamics (QED) modifications to Maxwell's equations due to vacuum polarization dissipate in dense sectors, where vacuum effects are negligible compared to plasma effects and, therefore, the vacuum refractive index -∆ vacvanishes from the dispersion relations of Eq. 3.1 in certain circumstances -see appendix C.
Landau-Zener-Stueckelberg-Majorana formalism [59][60][61][62] provides an analytic solution to the mixing equations of a two-state system with an energy gap between the states which depends linearly on time. 1 The transition probability between the states -ζ, η -reads P η↔ζ = 1 − e −2πΓ , with Γ the adiabaticity parameter of the system. Now, in order to writedown a compact analytic expression from Eq. 3.1, we consider pure non-adiabatic conversion so that the characteristic parameter results Γ<< 1 at resonance [15,21,25,26]. A secondorder series truncation yields P γ φ 2π|Γ|, with The expression in Eq. 3.2 unveils interesting physics. The conversion probability scales with the frequency of the photon from the NS, ω. Contrarily, the axion mass range for which resonant conversion takes place is set by the density profile through the magnetosphere, since ω p = m φ at resonance. That is conceptually similar to the case of haloscopes, where the scanning frequency is only given by the resonant frequency of the cavity [64]; or light-shiningthrough-walls (LSW) experiments, where the mass of the axion emerging in the axion-photon vertex is determined by the resonant frequency of the cavity independently of the nominal frequency of the laser pump [65][66][67]. Consequently, the conversion probability depends only on two parameters for the stated model of star and axion, i.e., z c and ω. Moreover, given the ω factor in Eq. 3.2, non-adiabatic conversion is inconsistent when entering the γ-ray energy frame. That is the reason why a γ-ray analysis is out of the scope of the article. Finally, the production of axions with lower coupling factors is enhanced at high frequency.
In order to render feasible an analytic approach, we consider a stationary magnetic field and a strongly magnetized plasma [15, 16, 18-21, 25, 26]. On the other hand, QED corrections suppress the cross-section by orders of magnitude when vacuum polarization effects are significant [12]. Consequently, our non-adiabatic resonant model is ruled out in several sectors. First, O-modes are evanescent below the cut-off frequency of the plasma [68,69], when ω p > ω. Second, adiabatic regions with Γ 1 appear near the stellar surface for Xrays. Third, for sectors with ω ω B , with ω Bc = eB(z c )/m e , standard simplifications in the dielectric tensor do not hold -see appendix B -. Finally, the parameter space at which QED corrections govern the dispersion relations in Eq. 3.1, or equivalently ω 2 p < Q QEDconsult appendix C.
The time-scale of the photon-axion system at resonance, τ s ∼| ω − ω p | −1 , is typically much shorter than the pulsation time-scale of the star -from the tens ns up to few ms in the radio-frequency regime [70]; significantly larger in the case of X-rays [30] -. Then, τ s << t pulse and the conversion takes place in a pseudo-stationary regime instead of transitory, where the model could loss generality.

Results
We have established in previous sections that, once the axion and NS models are assigned, the non-adiabatic resonant conversion probability in Eq. 3.2 depends only on the frequency of the photon emitted by the NSω -and the distance from its surfacez c -. The latter parameter determines the charge density via Eqs. 2.1 and 2.2 -. The resulting spectral flux density of the axion beam produced in the axion-to-photon conversion reads Finally, we define the dark matter (DM) spectral flux density where ρ ⊕ is the occupation number of dark matter around the Solar System position, about ρ ⊕ ∼ 300 MeVcm −3 , and ∆ν φ ∼ 10 −6 ν is the line-width of the axion signal for the virial velocity v φ ∼ 10 −3 c.
In Fig. 4 we report 1D calculations of the axion flux density induced by non-adiabatic resonant axion-to-photon conversion -Φ φ -in units of local dark matter flux density -Φ DM -through the light cylinder of neutron stars for benchmark parameters. Several inferences can be drawn from the analysis of the different panels. For a fixed pair multiplicity factor, we consider both integrated radio-wave spectral density functions and persistent X-ray flux densities from Fig. 2. The characteristic frequency of the plasma scales with the root of the charge density -see Eq. 2.2 -. Consequently, near the stellar surface, for κ 10 6 cutoff frequencies of the order of a few hundred GHz arise; and the corresponding parameter space is neglected for frequencies of a few tens of GHz as ω < ω p and hence O-modes are not propagated while the resonant conversion is disrupted -lower panels -. On the other hand, the ratio of plasma effects to vacuum effects on the mixing relations depends on both plasma density and photon energy -see appendix C for details on the QED suppression of axion-photon cross-sections in magnetic stars -. From Eq. C.1, it is found that axion-photon conversion is strongly suppressed below κ 10 4−5 for keV photons since ω 2 p << Q QED and vacuum birefringence governs the dispersion relations -upper panels -; or when κ < 10 7 for photons with higher frequencies of several keV -bottom-right panel.

Summary and conclusions
In this article we have treated the non-adiabatic photon-to-axion resonant conversion in the state-of-the-art picture of neutron stars. Relying on observational data in both the radiowave domain and the X-ray range, we have applied the model, of remarkable simplicity, to rotation-powered pulsars and magnetars. Within a realistic parameter space, results in Fig.4 suggest that persistent flux densities up to 19 orders of magnitude over the dark matter flow in the Galactic halo, around 14-15 Wm −2 , might originate somewhere in the polar caps of those astrophysical objects, region from which a dense beam of axions could be produced by photon-to-axion conversion. If we take into account that the relativistic axion velocity and the velocity dispersion of the surrounding dark matter differ by three orders of magnitude, this equates to dark matter overdensities of up to 16 orders of magnitude located somewhere in the vicinity of the star's light cylinder, in a large frequency-mass parameter space.
The conversion probabilities obtained from Eq. 3.2 would be inaccurate in sectors in tension with the following limits: ω p > ω, Γ 1 and ω ω B -i.e., sectors where O-modes are evanescent, adiabatic sectors or distant regions where the field becomes relatively weaker, respectively -. Consequently, we ruled out the model in these regions. In spite that the oscillation amplitudes tend to be lower beyond these limits and therefore our interest is reduced accordingly, methods that can be used to extend the model to those sectors -shaded in Fig. 4 -have been suggested, e.g., multi-dimensional numerical calculations can extend the ω << ω B sector by enabling a 'weak' field and also the propagation of oblique modes extending the ω p > ω space, ray-tracing technique, etc. [16,20,73,74]. On the other hand, the claim that Quantum Electrodynamics (QED) modifications to Maxwell's equations due to vacuum polarization suppresses axion-(X-ray)photon conversion in highly magnetic neutron stars has been entrenched for many years [12]. As a result of this, the possibility of relativistic axions with keV energies originating from photon-axion conversion in magnetic stars, or X-rays, bombarding the Earth causing observable effects has been neglected for decades. However, we have found, due to a higher pair density consistent with the current picture of neutron stars, a broad ω 2 p > Q QED sector at which the tension with QED vanishes in a more realistic neutron star model beyond the anachronistic Goldreich&Julian's density profile adopted firstly in [12] and later in more contemporary works. The numerical approach of [75] would allow the simulations to be extended to the parameter space where vacuum polarization effects are relevant, i.e., ω 2 p < Q QED , where the conversion probability is, however, suppressed by many orders of magnitude. Finally, some estimates might be in tension with the altitude from the stellar surface at which the radiation starts -i.e., regions with z c h 0 -. We note two points. First, different parametric configurations result in a higher/lower resonance altitude. Second, we have discussed in section 2 that this parameter presents a significant uncertainty and can be ad hoc at times. Aiming at keeping our work as model-independent as possible, it was not ingrained more deeply throughout the manuscript.
Neutron stars have the stronger magnetic fields known in the universe, about B 10 11 T at surface and axion convert efficiently to photon, and vice-versa, in a strong magnetic field. As ambient axions fall on the star, a narrow line, with ∆ν < < ν, appears at a frequency corresponding to the axion mass. This signal has been searched for by different radio telescopes, without success. Unfortunately, the authors of the corresponding works adopted a star model with a Goldreich & Julian (GJ) density profile to compare their observational data against a simulation that predicts the shape, frequency and intensity of the emission line originating from axion-to-photon conversion -see the most recent and restrictive results in [76][77][78] -. Crucially, astronomical observations of neutron stars and emission models allow one to estimate the observer-frame pair density, n e , reliably and, by this method, the Goldreich & Julian density profile adopted by those authors, dated in the late 1960s, has been shown to be under-estimated by many orders of magnitude in more recent works [43]. In other words, typical densities in neutron star magnetospheres concordant with the currently available data establish plasma densities above κ 10 4 times denser than the Goldreich & Julian model, with κ = n e /n GJ the so-called 'pair multiplicity factor'. This has important consequences. First, it is well-known that O-mode electromagnetic waves do not propagate in a plasma below a cut-off frequency, ω p , which scales with the plasma density profile with distance, n e (r), in the form ω p ∝ n 1 /2 e → ω p ∝ κ 1 /2 -see Eq. 2.2 -; hence suppressing the conversion of axions with a lower mass near the star surface, where the field would be stronger enhancing the conversion probability, P, as P φ γ ∝B 2 , since the conversion takes place at m φ ∼ ω p resonance. Moreover, such a large pair density implies that the exclusion bounds resulting from those previous works, set at m φ ∼ 15−30 µeV, are overestimated by several orders and, in fact, fall below the CAST exclusion sector and vanish [79], since the resonance of axions of dozens of µeV of mass is only possible far away from the star surface, where the field is much weaker, as B∝r −3 , and the conversion probability is reduced accordingly. Therefore, a significant blueshift of the spectral line emerging from axion-to-photon conversion is expected with respect to previous estimates, from frequencies of several GHz to frequencies of several hundred GHz, or m φ 100 µeV. In consequence, caution should be exercised with earlier simulations that adopt the Goldreich & Julian density profile in predicting radio signals originating from the resonant conversion of axions into O-mode photons in magnetospheres; while higher-frequency telescopes must be used for scanning, possibly leading to new exclusion limits for axion at high-frequency.
The direct detection of neutron star axions at the Earth position by dark matter halohelio-tele-scopes [64,80,81] is an interesting question. Relativistic aberration transforms the regular semi-isotropic dipole pattern of radio-emission into a narrower beam of light with a high directionality in the observer's line-of-sight. The emerging axion beam is shaped accordingly because of momentum conservation, ideally. It also suffers from the inverse square law with distance. Geometric dilution would truncate prospects for the direct detection of the stellar axion flow even for features as giant radio pulses (GRPs) [82][83][84] with a typical duration of the order of few-tens ns with a scatter time-scale of the order of the ms and a transient flux up to 10 4−5 times over the unpulsed component observed by radio-telescopes. The high-energy frame -more specifically soft X-ray -should be examined in more detail. Dark matter astroparticles released from neutron stars could be analyzed in units of the axion rate received from the Sun [64]. The solar axion spectra presents a distribution peaking around 2-3 keV with a flux density about 10 9 cm −2 s −1 for g φγ ∼ 10 −11 GeV −1 [13,80]. In this band, the quiescent fluence of magnetars can be O(10 −1 ) cm −2 s −1 -see Fig. 2 -while the photon-to-axion conversion amplitude can be >10 −1 for a broad parameter space. Thus, solar axion flux would be 10 12 times denser than axion magnetar flux at the Earth position due to geometric dilution. Nonetheless, intermediate bursts (IBs) featured by magnetars, with luminosities up to 10 43 erg s −1 and a characteristic time of tens of seconds, and mostly giant flares (GFs), which present a spike lasting ∼0.1-1 s that can reach powers up to 10 47 erg s −1 and a long tail -few hundred seconds -have been observed at ω>2 keV [30,31,[85][86][87][88][89]. They represent a contrast to their quiescent luminosities in the range 2-10 keV, with a break about 10 33 erg s −1 for -the less bright -'transient' magnetars, and 10 33 erg s −1 forthe more potent -'persistent' sources. In spite of we emphasize that this result must be viewed with caution, 2 numbers could add up so that the direct detection of magnetar axions triggered by X-ray IBs/GFs cannot be discarded [90].
self-absorption spectral function of an optically thick pulsar presents a power-law with a very steep slope The spectral indexes are related by 2α+1 = δ. For simplicity, we have approximated ω b ∼ ω τ , corresponding ω τ to the frequency which maximizes the optical path. This assumption would lead to deviations within the precision that is desired for the present work.

B Dielectric tensor of a cold plasma in a strong homogeneous magnetic field
Longitudinal modes in a plasma are compressive waves. In the equation of motion, their respective restoring force is described by a ∇p term. For considering a cold homogeneous plasma, we neglect the ∇p pressure. Withm≡ẑ the coordinate along the magnetic field B, x andŷ directions across the magnetic field, the effective dielectric tensor of a cold plasma presents the form [58] ε xx = ε yy = 1 + where ω p is defined by Eq. 2.2,ω = ω − k z v , γ = (1 − v 2 /c 2 ) −1 /2 and ω B = e|B|/m e with |B|= B 0 , e the elementary charge and m e the electron mass.
The electrodynamic properties of a pulsar magnetosphere can be described in the strong field limit B 0 → ∞, where ω B >> ω and ω B >> ω p . In this approximation, easily fulfilled in the vicinity of the star, where the field is stronger, or for sufficiently low frequency, the dielectric permittivity becomes ε xx = ε yy = 1, ε xy = ε yx = ε yz = ε zy = 0, and This ensures that the conductivity is non-zero only along the direction of the magnetic field.

C The role of the density profile in the weight of QED corrections
The claim that QED modifications to Maxwell's equations due to vacuum polarization suppresses axion-X-ray conversion in highly magnetic stars has been entrenched for decades [12]. However, the tension with QED is relaxed in the modern neutron star model beyond the Goldreich & Julian density profile adopted throughout this manuscript [39,43].
Corrections due to QED effects arise from the Euler-Heisenberg Lagrangian L QED = α 2 90m 4 e [(F µν F µν ) 2 + 7 4 (F µνF µν ) 2 ] in strong external fields. Thus, the leading order terms of the refractive indices for plasma and vacuum are ∆ p ∼ ω 2 p/2ω and ∆ vac ∼ 7α 90π ω( B(r) /B crit ) 2 , respectively; being α = e 2 /4π the fine-structure constant and B crit = m 2 e/e the 'critical' field. Therefore, the ratio of plasma effects to vacuum effects -∆p /∆vac -is read [17]  where we have introduced the plasma frequency profile from Eq. 2.2 and defined Q QED = 2ω∆ vac . From the analysis of Eq. C.1 it follows that no observable vacuum birefringence effects are expected for low-frequency photons, as ω 2 p >> Q QED and plasma effects govern the dispersion relations in Eq. 3.1; while in the case of high-energy photons the weight of vacuum polarization effects gradually vanishes at relatively large distances from the stellar surface for pair multiplicity factors above κ 10 5 , consistent with the current picture of neutron stars.

D Axion electrodynamics in neutron stars
A modification of Maxwell's equations arises from a light, pseudo-stable axion [92] ∇ · (E + g φγ φB) = ρ , where the monopole density ρ m and the monopole current J m were replaced by zero. From the linearisation of the equations of motion from Eq. 1.1 for a stationary and free of electric charge neutron star background in the aligned rotator approximation, the following system of equations emerges [21] + m 2 φ φ = g φγ E · B 0 , (D.2a) where it was introduced A ← A 0 + A in fields and densities while J = σ · E, being σ the conductivity tensor. Note the Klein-Gordon formalism in Eq. D.2a.
E Gradient of the plasma frequency in the magnetosphere of neutron stars In the aligned rotator limit -m ẑ -and on-the-spot resonance approximation -where the magnetic field presents a negligible gradient -, from Eqs. 2.1 and 2.2 we get In that form, it follows that ∂ ∂z ω p = − 3 /2 z −5/2 ξ. At resonanceω p = m φ -it is straightforward to obtain where we have used m φ = z −3/2 c ξ from Eq. E.1.