Photon-axion mixing in thermal emission of isolated neutron stars

Thermally emitting neutron stars represent a promising environment for probing the properties of axion-like particles. Due to the strong magnetic fields of these sources, surface photons may partially convert into such particles in the large magnetospheric region surrounding the stars, which will result in distinctive signatures in their spectra. However, the interaction depends on the polarization state of the radiation and is rather weak due to the low experimentally allowed values of the coupling constant $g_{\gamma a}$. In this work, we compute the degree of photon-axion transition in the case of 100% O-mode polarization and spectral energy distribution of an isotropic blackbody with uniform surface temperature. The stellar magnetic field is assumed to be dipolar. We show that with the maximum effect reached for the magnetic fields $\sim10^{13}$ - $10^{14}$ G (typical for X-ray dim isolated neutron stars) and $g_{\gamma a} = 2 \times 10^{-11}$ GeV$^{-1}$, the optical flux is reduced by 30 - 40%, while the high-energy part of the spectrum is not affected. The low-energy decrease exceeds 5% at $g_{\gamma a} \geq 2 \times 10^{-12}$ GeV$^{-1}$ and $m_a \leq 2\times 10^{-6}$ eV, which is below the present experimental and astrophysical limits on axion parameters. To obtain the actual observational constraints, rigorous treatment of the radiative surface layers is required.


Introduction
Pseudoscalar axion-like particles appear in various extensions of the Standard Model of particle physics. Proposed in the late 1970s to explain the absence of CP violation in strong interactions [1,2], they have transcended the original context and were employed to solve a variety of problems in physics and astrophysics [3,4,5,6]. In this work, we refer to all such particles simply as axions.
The defining property of axions is the conversion into photons and vice versa in the presence of an external magnetic field [7]. Many laboratory-based searches exploited this interaction in attempts to prove their existence, but so far no satisfactory observational evidence has been found, and the stringent bounds on the coupling constant over a wide axion mass range were obtained by the CERN Axion Solar Telescope (CAST), g γa ≤ 6.6 × 10 −11 GeV −1 (95% C.L.) [8]. For this reason, studying manifestations of axions in various astrophysical environments has become increasingly popular. A common strategy used in such works is to search for spectral features in sources of radiation observed through large-scale regions with relatively strong magnetic fields. These include active galactic nuclei such as NGC 1275 [9,10,11] (see also [12]), blazars and quasars [13,14], supernovae [15], etc. By comparing the detected spectra with the ones arising after the mixing was taken into account, several authors were able to impose astrophysical constraints on a range of axion parameters, exceeding the CAST limit.
More recently, considerable research has been devoted to studying the axion-photon interaction in strongly magnetized objects, namely magnetic white dwarfs (mWDs) and neutron stars. In the case of white dwarfs, photons emitted by the photosphere may undergo conversion as they propagate through the encompassing gaseous atmosphere [16]. If the plasma is sufficiently tenuous, this could lead to an enhancement of the observed degree of linear polarization. Since in most white dwarfs its value is very low, such interaction can be used to constrain a range of axion parameters. In particular, for mWD PG 1015+014 the induced polarization degree should exceed the observed value of 5% at g γa ≥ 2 × 10 −11 GeV −1 , which is below the CAST limit.
In the case of neutron stars, the presence of axions can manifest itself in a number of ways. First, their magnetospheres can be a site for the conversion of dark matter axions [17,18]. Also, axions could be produced thermally inside the core with energies from several keV to several hundred keV through the nucleon bremsstrahlung processes, and escape the interior due to the feeble interaction with matter. Their subsequent transition into photons in the magnetic field surrounding the star could lead to observable signatures in its hard X-ray spectrum [6,19,20]. Magnetars are natural targets for such investigation due to the extremely strong magnetic fields and high internal temperatures, and their study led to further astrophysical limits on axion parameters [21,22,23].
At the same time, manifestations of the reverse process have not yet been fully explored. Spectra of neutron stars commonly include a thermal component, and as the surface radiation propagates through the large magnetospheric region with a relatively strong magnetic field, a fraction of the initial flux may oscillate into axions. Although in most of the magnetosphere the probability of conversion is suppressed by the QED effects, previous studies have shown that the mixing can lead to distinctive spectral signatures [24,25]. However, the coupling values they considered are now excluded experimentally. The purpose of our work is to find out whether current constraints on axion parameters allow for observable effects, and to estimate their possible magnitude.
X-ray dim isolated neutron stars (XDINSs), also known as the "Magnificent Seven", are particularly suited for such investigation. These isolated, radio-silent sources exhibit purely thermal spectra with no high-energy power-law component, which is often found in other classes of neutron stars. Their X-ray luminosity of 10 31 − 10 32 erg s −1 is consistent with the surface blackbody emission with temperatures kT ∞ ∼ 40 − 100 eV, and the radiation radii of a few kilometers [26,27]. The spin periods P ∼ 3 − 12 s and period derivativesṖ ∼ 10 −14 − 10 −13 s s −1 translate into spin-down magnetic fields B ∼ 10 13 − 10 14 G.
In the strongly magnetized vacuum around a neutron star, the photons propagate in two normal modes, the extraordinary (X) and ordinary (O) mode, and only those polarized in the latter can convert into axions. Although the physical conditions of the outermost layers of XDINSs are still a debated issue, several models have been proposed to explain the observed properties of these sources [28], and the polarization patterns they predict are different. The conventional picture is that cooling neutron stars are covered with a gaseous atmosphere that reprocesses the emission coming from the outermost surface layers. As a result, the radiation acquires a net polarization predominantly in the X-mode, since the strong magnetic field substantially reduces its opacity while leaving that for the O-mode almost unchanged [29,30,31]. On the other hand, it has recently been appreciated that the stars with a low surface temperature ( 100 eV) and a strong magnetic field ( 10 13 G) may have a liquid or solid condensed surface due to a phase transition in the outermost layers [27,32,33,34]. Thus, in this case, the atmosphere may be absent. The emissivity of a naked condensed surface in both X-and O-mode is of the same order in a broad energy range [29,35].
Since the magnitude of photon-axion transition will depend on the composition of stellar surface layers, it is essential to estimate its possible magnitude at the currently allowed values of axion parameters. To this end, here we consider a simple picture in which the thermal emission is an isotropic blackbody with a uniform surface temperature, and is 100% polarized in the O-mode. These conditions provide the highest degree of conversion, and our results will show whether a more detailed study can lead to astrophysical constraints that exceed the present experimental limits.
Our paper is organized as follows. In Section 2, we summarize the main properties of the photon-axion interaction in a strong magnetic field. Then in Section 3, we calculate the resulting degree of conversion and determine the range of axion parameters m a , g γa at which the mixing has a significant effect. The observability of our results and the prospects for future work are discussed in Section 4.

Theoretical framework
In this section we outline the physical basis for our calculations. Throughout the paper we use Lorentz-Heaviside units with = c = 1 and α = e 2 /4π ≈ 1/137.

Polarization of radiation
In the presence of an external magnetic field, the vacuum around the star behaves as a birefringent medium, in which the photons are polarized in two normal modes: the ordinary mode (O-mode), with the electric field oscillating in the plane of the propagation vector k and the local magnetic field B, and the extraordinary mode (X-mode), with the electric field oscillating perpendicularly to both k and B. According to the QED, the dielectric and magnetic permeability tensors of the vacuum are modified by virtual electron-positron pairs, which affects the polarization properties of radiation primarily in two ways.
First, the photons emitted at the surface maintain their initial polarization state up to a characteristic distance [36,37,38], the polarization radius r p , which depends on the photon energy and magnetic field strength. In this region, the wave electric field can instantly adapt its direction to that of B. Around r p , the coupling weakens, and the photon electric field is no longer able to properly follow the variation of the local magnetic field, until for r ≫ r p , its direction freezes. If the radiation propagates radially in a dipolar field, its initial polarization state is preserved because the direction of transverse part of B is always constant. This does not hold for the photons, trajectory of which has been modified by the gravitational effects [39]. In general, to account for the QED-induced depolarization at r > r p , it is necessary to solve the system determining the evolution of the polarization state [38,40]. However, the magnitude of this effect is relatively small, because the dipolar magnetic field does not significantly change its direction outside the polarization radius. According to our calculations, for the energies considered in this paper, the angle between the wave electric field and the transverse part of B for both modes changes only slightly from its initial value in the region that will be of interest to us. For this reason, we can assume that the radiation always retains its initial mode without losing significant physical accuracy.
Second, the refractive indices for each mode are different and can deviate from their vacuum value of unity. For the photon energies well below the electron rest mass, the Euler-Heisenberg effective Lagrangian can be applied [41,42], which for arbitrary values of B = |B| leads to where Θ is the angle between k and B, b is the magnetic field strength in units of the critical field B crit = m 2 e /e = 4.4 × 10 13 G.
The magnetized plasma surrounding the star modifies the dielectric tensor similarly to the vacuum polarization. In the following, we assume that it has the Goldreich-Julian charge den- where Ω is the stellar angular velocity [43]. Note that ρ GJ estimates only the net charge density, and the total charge density may be larger due to the presence of electronpositron pairs [44]. However, since there is no reliable probe of the plasma density in the magnetospheres of neutron stars, we will only consider the Goldreich-Julian value as a conservative lower limit.

Photon-axion mixing
The oscillations of a photon with energy ω and an axion in the external magnetic field are described by the following interaction: where F µν is the electromagnetic stress tensor,F µν is its dual, a is the axion field, m a is the axion mass, and g γa is the coupling constant [7]. By adding the above term to the Euler-Heisenberg effective Lagrangian and linearizing the wave equations, we obtain the following system: where A and a denote the amplitudes of the ordinary photon state and the axion, respectively, and the z− axis is along k. Note that only the radiation polarized in the O-mode can convert into axions, and the X-mode flux always remains intact. The components of the mixing matrix are: where ω 2 p = 4παn e /m e is the plasma frequency, n e = ρ GJ /e is the electron density [45].
While the converted fraction of the initial photon flux can only be obtained by numerically integrating Equations 5, the strength of the mixing can be estimated as follows. The external magnetic field varies along the photon trajectory over a lengthscale l B = B/|k · ∇B| = r/3, where a hat denotes a unit vector. If we take the value of B to be constant over l B , the probability of a photon-axion transition after traveling this distance is where ∆ 2 osc = ∆ + ∆ p − ∆ a 2 + 4∆ 2 M , so that the oscillation length is l osc = 2π/∆ osc . Strong mixing occurs when the conversion probability is close to unity, which is satisfied if

Numerical simulations
In this section, we calculate the observational signatures of photon-axion mixing in the O-mode flux of a thermally emitting neutron star, and determine the range of axion parameters m a , g γa at which the effect is significant.

Computational model
Our treatment of the surface radiation relies on a number of simplifications which are discussed below. We stress that our primary goal is to assess whether the mixing can have a noticeable effect on spectra of XDINSs at the currently unconstrained values of m a , g γa and whether there is room for future investigations, rather than to derive theoretical predictions to be compared with the observational data. For this reason, we avoid dealing with the complex structure of the outermost stellar layers and consider the most favorable conditions for the detection of photon-axion conversion.
The stellar magnetic field is assumed to be an aligned dipole, components of which in polar coordinates (r, θ) take the form where B p is the surface polar field strength. The functions f dip , g dip account for the relativistic corrections [46,47]. The stellar radius, mass, and period are R s = 12 km, M s = 1.4M ⊙ , and P s = 7 s, respectively. Initially, the radiation is emitted from the cooling stellar surface with an isotropic blackbody distribution and uniform temperature kT ∞ = 50 eV. The seed photons are 100% polarized in the O-mode, which provides the maximum degree of conversion into axions. As the radiation moves away from the emission point, we numerically integrate Equations 5 along its trajectory, accounting for the gravitational redshift. The observer is located at infinity, with the line of sight inclined at an angle χ with respect to the spin axis, and collects the contribution of each part of the surface which is into view due to the effects of ray bending [39].

Results
To get a clear picture of the conversion process, let us first study the mixing in an individual photon ray. The simplest case is the radiation emitted in radial direction (colatitude of the emission point is the same as that of the observer), so that its trajectory is not modified by gravity. We take the coupling constant g γa = 2 × 10 −11 GeV −1 (the highest value allowed by the study of mWD PG 1015+014 [16]), axion mass m a = 10 −8 eV, and photon energy ω ∞ = 1 eV. The polar magnetic field Figure 1: Evolution of the photon and axion square amplitudes in the case of radiation emitted radially, as a function of the distance r from the stellar surface. The parameters of the system are ω ∞ = 1 eV, g γa = 2 × 10 −11 GeV −1 , m a = 10 −8 eV. The magnetic field geometry is an aligned dipole with B p = 10 13 G and χ = 60 • . The probability of conversion beyond r = 500 km is vanishingly small, this region is not displayed. strength is B p = 10 13 G (representative of XDINSs [27]) and the observer inclination angle is χ = 60 • . The resulting plot of the photon and axion square amplitudes is shown in Figure 1.
As we can see, significant conversion takes place at a distance of several stellar radii away from the surface. The system does not enter the strong-oscillation regime (in this case, |A | 2 and |a| 2 would have the form of sinusoidal curves) because Conditions 8 are never met simultaneously. In particular, near the emission point, where l osc is small compared to l B , ∆ is much larger than 2∆ M , and therefore the mixing is suppressed. As the radiation moves away from the surface of the star, ∆ eventually becomes comparable with 2∆ M (∆ a is negligible and |∆ p | ∆ ), which is due to their different dependence on the radial coordinate: ∆ ∼ 1/r 6 , while ∆ M ∼ 1/r 3 . However, in this region, l osc is already 5-10 times greater than l B , which results in the conversion probability being less than unity. Still, these conditions are sufficient for the mixing to have a noticeable effect. At relatively large distances, the conversion probability becomes negligible.
We then perform the calculations for a broad energy range (from 1 eV to 10 keV), integrating the intensity over the visible part of the stellar surface. The result is conveniently expressed in terms of the differential modification factor Λ, which is defined as the fraction of the initial photon flux with a given energy, reaching the observer without conversion. In particular, if Λ = 1, the flux will not be modified, and if Λ = 0, all photons will turn into axions and thus no signal will be detected. The data for g γa = 2 × 10 −11 GeV −1 , m a = 10 −8 eV, and B p = 10 13 G, χ = 60 • are presented in Figure 2, as a function of the redshifted photon energy.
Our results show that the mixing is most significant at the optical energies, where the O-mode flux is almost halved. At the same time, it has little effect on the X-ray emission, which is due to the fact that the QED contribution ∆ increases linearly with increasing ω, while ∆ M does not depend on it. At higher photon energies, ∆ and ∆ M become of the same order in the re- gion where the oscillation length is much larger than l B , which results in a considerable decrease in the conversion probability. The shapes of the photon and axion square intensities are similar to those shown in Figure 1, except for the oscillations having smaller amplitudes. It is noteworthy that the converted fraction of the flux varies by no more than a few percent between different ray trajectories, which due to the fact that the system remains in the weak-oscillation regime. For the same reason, its value does not change noticeably for other observer inclination angles.
We then repeat our calculations for different values of m a , g γa to obtain the range at which the mixing has a significant effect on the optical radiation. Figure 3 illustrates the modification factor at ω ∞ = 1 eV for B p = 10 13 G and a given choice of axion parameters. The optical photon flux is reduced by more than 5% at g γa ≥ 6 × 10 −12 GeV −1 and m a < 2 × 10 −6 eV. This result can be explained using Conditions 8: for lower values of the coupling constant, ∆ M does not get close enough to ∆ in the region with a relatively small oscillation length; for higher axion masses, ∆ a becomes of the same order as both ∆ M and ∆ , suppressing the conversion probability. The plasma contribution ∆ p becomes important at lower values of g γa , where ∆ ≫ ∆ M almost everywhere in the magnetosphere. At some point, the terms ∆ p and ∆ become of the same order and cancel each other because of their opposite signs. Even though the oscillation length in this region is much larger than l B , this results in a slight increase in the conversion probability. We also perform the same calculations for a greater polar field strength B p = 10 14 G, which can also be present in XDINSs [27]. As shown in Figure 4, in this case the converted fraction reaches 5% at even lower values of the coupling constant, g γa ≥ 2 × 10 −12 GeV −1 and m a ≤ 2 × 10 −6 eV.

Discussion and conclusions
In this paper we studied the photon-axion conversion in the magnetosphere of a thermally emitting neutron star. Assuming the dipolar magnetic field and the spectral energy distribution of Figure 3: Contour plot of the O-mode modification factor Λ for ω ∞ = 1 eV, as a function of the axion parameters m a , g γa . The magnetic field geometry is B p = 10 13 G, χ = 60 • . We also include the values of g γa allowed by CAST [8]. an isotropic blackbody with a uniform surface temperature, we have shown that the mixing can result in a significant decrease in the optical O-mode flux at the allowed values of m a , g γa . Further investigation may lead to astrophysical constraints on axion parameters that exceed the present experimental and astrophysical limits.
We found that the photon-axion system exhibits only weak oscillations at the currently allowed values of the coupling constant. The reason for the decrease in the conversion probability is either the dominance of the QED contribution ∆ , or the fact that the oscillation length is much larger than the scale length of the magnetic field. This is in agreement with the study of [25], who focused on higher photon energies and currently excluded values of g γa . According to their results, weak oscillations are present in the X-ray range at g γa = 10 −10 − 10 −9 GeV −1 and the converted fraction depends inversely on ω ∞ , similarly to the case shown in Figure 2. Only at a relatively high value of g γa = 10 −8 GeV −1 , the conversion probability becomes equal to unity in a large area of the magnetosphere, and the system enters the strong-oscillation mode. As a result, Λ depends on the photon energy in a more complex way and, in particular, can increase with increasing ω ∞ . Such relation should not be confusing, since it is simply a consequence of the strong-oscillation mode, which is substantially different from the weak-mixing regime discussed in our work [7,24]. Overall, greater values of the coupling constant could be more easily constrained by looking into the X-ray energies where the emission of XDINSs typically peaks, while at lower values of g γa one has to consider the mixing in the optical band.
Our results exceed the limits obtained by [16] in their study of mWD PG 1015+014. At the same time, observations of the magnetic white dwarfs with higher surface field strengths could potentially exclude a wider region of the axion parameter space, which will be surpassed by the limits shown in Figures 3 and  4 only at high axion masses. However, these results should be interpreted with caution, as they are highly dependent on the assumptions about the atmospheric plasma composition. In particular, the authors of [16] consider the case of a relatively sparse hydrogen plasma with a surface density of ρ 0 = 10 −10 g cm −3 and a barometric profile. They also note that there is no clear agreement on the surface plasma density of magnetic white dwarfs, and its actual value is expected to lie between 10 −11 and 10 −6 g cm −3 . In this regard, the study of XDINSs may provide a more robust method of imposing constraints on axion parameters, which, however, will depend on another kind of assumptions.
A major challenge in assessing the effect of mixing on spectra of XDINSs is the uncertainty in the structure of their radiative surface layers. A well-known feature of the "Magnificent seven" is that their optical and ultraviolet counterparts exceed the extrapolation of X-ray blackbody at low energies by a factor 5 -50 [48], and numerous authors have made extensive use of atmospheric and condensed surface models in attempts to explain its origin. A strongly magnetized atmosphere has been suggested for RX J1605.3+3249 [49]. The observed properties of RX J0720.4-3125 were attributed to the emission of a bare condensed surface [50], and a combination of the latter and a thin, magnetic, partially ionized hydrogen atmosphere was used to fit the spectral shape of J1856.5-3754 [51] and RX J1308.6+2127 [52].
The polarization pattern induced by both models can be conveniently expressed in terms of the intrinsic polarization degree (that at the source): where F X,O is the monochromatic, phase-averaged flux in each mode, obtained by integrating the intensity over the visible part of the stellar surface [29]. By definition, Π EM L = 1 for the radiation 100% initially polarized in the X-mode, Π EM L = −1 for that in the O-mode, and Π EM L = 0 for the unpolarized emission. According to the results of [29], in the case of a fully ionized hydrogen atmosphere, Π EM L ≈ 0.99 in the X-ray band and Π EM L ≈ 0.87 in the optical band for all possible viewing geometries. The radiation of a condensed surface in the X-ray band is almost unpolarized, |Π EM L | 0.07, while in the optical band Π EM L can reach −0.5 for the most favorable viewing geometries, with the emission being mostly polarized in the O-mode. Since the polarization state of thermal radiation directly affects the magnitude of photon-axion conversion, the effect will be observable only if the flux contains a considerable fraction of O-mode photons. From this point of view, the emission of a condensed surface is a promising place to search for its manifestations, but the result is likely to be lower than that obtained in Section 3, since we considered the liming case of 100% O-mode polarization. The presence of an atmosphere may make the interaction less noticeable because its radiation is predominantly polarized in the X-mode.
Rigorous treatment of the outermost stellar layers provides several ways of imposing constraints on a range of axion parameters. One approach is to simply measure the change in the low-energy spectra, as shown in Figure 3. Although the optical counterparts of XDINSs are rather faint, observational data have been presented for all sources [48]. However, this method will not be sensitive enough if the converted fraction is relatively small, which is expected at low values of the coupling constant g γa . In such cases, it may be more efficient to compute the polarization observables, linear polarization fraction Π L and polarization angle χ p [38], which accurately describe the polarization state of the detected radiation and can reflect a slight decrease in the O-mode flux. Since various types of stellar surface composition are expected to have different polarization patterns [29], constraints can be based on a comparison of Π L and χ p , predicted by the models which do and do not take the mixing into account. Additionally, one can make use of the observational polarimetric data, as described in [53] for the case of RX J1856.5-3754.
Another potentially promising approach is to consider the photon-axion interaction in magnetars. Their magnetic fields are stronger than those of XDINSs B p ≈ 10 15 G , and are believed to include higher multipoles and global or local shears (the "twisted magnetosphere" model, see [54]), which in general can increase the mixing strength. In addition, their lowenergy radiation may include a considerable fraction of the Omode photons [55,56]. Such investigation, however, is ham-pered by our present poor knowledge of both the composition of surface layers and the magnetospheric structure of magnetars, as well as the faintness of the optical counterparts of their thermal emission [57]. At the same time, considering other classes of neutron stars, such as radio pulsars, is unlikely to produce greater constraints on axion parameters. Due to their lower magnetic field strengths, the converted fraction will likely be smaller than in the case of XDINSs.