Signatures of photonaxion conversion in the thermal spectra and polarization of neutron stars
Abstract
Conversion of photons into axions under the presence of a strong magnetic field can dim the radiation from magnetized astrophysical objects. Here we perform a detailed calculation aimed at quantifying the signatures of photonaxion conversion in the spectra, light curves, and polarization of neutron stars (NSs). We take into account the energy and angledependence of the conversion probability and the surface thermal emission from NSs. The latter is computed from magnetized atmosphere models that include the effect of photon polarization mode conversion due to vacuum polarization. The resulting spectral models, inclusive of the generalrelativistic effects of gravitational redshift and light deflection, allow us to make realistic predictions for the effects of photon to axion conversion on observed NS spectra, light curves, and polarization signals. We identify unique signatures of the conversion, such as an increase of the effective area of a hot spot as it rotates away from the observer line of sight. For a star emitting from the entire surface, the conversion produces apparent radii that are either larger or smaller (depending on axion mass and coupling strength) than the limits set by NS equations of state. For an emission region that is observed phaseon, photonaxion conversion results in an inversion of the plane of polarization with respect to the noconversion case. While the quantitative details of the features that we identify depend on NS properties (magnetic field strength, temperature) and axion parameters, the spectral and polarization signatures induced by photonaxion conversion are distinctive enough to make NSs very interesting and promising probes of axion physics.
Subject headings:
stars: neutron — Xrays: stars — cosmology: miscellaneous1. Introduction
Axions with low mass (below 1 meV) are a direct prediction of the solution for the strong CP violation problem. Peccei & Quinn (1977) proposed axions as PseudoGoldstone bosons arising from the spontaneous breakdown of the U(1) symmetry. The PecceiQuinn axion occupies a narrow region in the axion parameter space of mass and coupling strength. However, axions (or pseudoscalar particles) can exist with more generic values for their mass and coupling strength, and are one type of dark matter candidate (e.g. Arvanitaki et al. 2009). This motivates efforts to constrain axion properties even if the parameters do not reach into the PecceiQuinn regime.
Several experimental avenues have been used to detect axions and constrain their properties. These fall into two categories: groundbased experiments and cosmologicalastrophysical efforts. Current limits on the axion mass () come from the latter approach. The cooling rate observed in supernova 1987A implies a direct axion mass upper limit: if axions are more massive, then cooling would be predominantly due to axions rather than neutrinos (Eidelman et al. 2004). The lower limit comes from cosmological considerations: if axions constitute dark matter, then they must be cold and not overclose the Universe (Preskill et al. 1983). Constraints on the coupling constant can be obtained with both astrophysical observations and groundbased experiments. The existence of horizontal branch stars with extended morphology, that is, stars with a range of surface temperatures, constrains the coupling constant of axionphoton interaction to be GeV (Raffelt 2008). Groundbased experiments also place similar constraints, and future experiments are anticipated that will obtain stronger limits (Andriamonje et al. 1983; De Panlis et al 1987; Wuensch et al. 1989; Asztalos et al. 2004, 2010).
Recent developments in understanding radiation dimming (Lai & Heyl 2006; Jimenez et al. 2011) and modification of the polarization pattern (Gill & Heyl 2011) have put further constraints on the plane. Limits on cosmological radiation dimming have also been used to place constraints on much lighter axion masses (Avgoustidis et al. 2010).
The potential use of neutron stars (NSs) as probes of axion parameters has been noted by a number of authors (e.g., Lai & Heyl 2006; Chelouche et al. 2009; Jimenez et al. 2011). In particular, occultation and eclipsing or transiting of a companion and dimming of the spectrum can potentially produce observable signatures. However, the phenomenon of occultation (when a background object passes through the influence of a NS magnetic field) was deemed to have too low a probability to be astrophysically relevant. In addition, there are no known binary systems involving a NS that are detached enough to yield a clean constraint (Jimenez et al. 2011). On the other hand, Lai & Heyl (2006) showed how the thermal Xray spectrum of a NS (modeled as simple blackbody emission) can be modified in a significant way by the presence of axions. Spectral features due to axionphoton oscillations have been discussed by Chelouche et al. (2009), and shown to be relevant for highly magnetized NSs in the submm wavelength range for an observationally interesting range of axion parameters.
In this paper, we focus on the soft Xray band in which NSs are routinely observed. The goal of our work is to perform a detailed calculation of the observable effects that the presence of axions would have on NS spectra, light curves, and polarization. Our analysis goes beyond previous computations in the literature for the subject by coupling the calculation of photonaxion conversion probabilities with detailed models of magnetized atmospheres and general relativistic models for emission from the NS surface. For the local surface brightness of the NS, we use accurate magnetized atmosphere models, which yield the angular and energydependent photon intensities in the two photon polarization modes. These models include the effect of mode conversion due to the vacuum polarization. The local intensity is then imported into a code for computing spectra and light curves from finite emission regions on the NS surface. The analysis incorporates the effects of gravitational redshift and gravitational light deflection and allows for arbitrary viewing geometries. The resulting spectra are then modified by the photonaxion conversion process, accounting for both the energy and angular dependencies of the conversion probability. Note that we solve the conversion between photon modes and the conversion between axions and photons independently, which is shown by Lai & Heyl (2006) to be a very good approximation for the regimes we explore in this paper. With our calculations, we are able to produce realistic models for the spectra, light curves, and polarization signal of NSs, and use these to identify specific features that carry the telltale signs of photonaxion conversion. We discuss these signs in the context of NS observations in §5.2.
Our paper is organized as follows: In §2, we discuss the formation of the NS thermal spectrum within a magnetized atmosphere, placing particular emphasis on the emergence of two photon polarization modes. In §3, we summarize the calculation of the photonaxion conversion probability for a magnetized object. The computation of the ‘observed’ thermal spectrum, light curves, and polarization with and without photonaxion conversion is described in §4, and the results are presented in §5. We summarize our findings in §6.
2. Photon propagation in the atmosphere of magnetized neutron stars
In this section, we describe the main properties of magnetized NS atmospheres and the processes that influence photon propagation in them. Axions do not directly couple to the processes described in this section. However, the details presented below are important for calculating the polarized radiation emerging from the NS atmosphere. As discussed in §3, axions couple only to photons that are linearly polarized in the plane formed by the directions of propagation and the magnetic field; the physics of photon propagation in the NS atmosphere determines the relative fraction of such photons.
Radiation emerging from the hot surfaces of magnetized NSs is the result of radiative transfer through a layer of ionized plasma. Due to the strong gravitational field of the NS, heavy elements settle quickly (Alcock & Illarionov 1980; Brown et al. 2002), and this atmospheric layer is composed of light elements, unless nuclear burning on the NS surface consumes the light elements (Chang et al. 2010, and references therein). Only recently have selfconsistent magnetic atmosphere models for partially ionized hydrogen and mid elements been constructed (see Mori & Ho 2007; Ho et al. 2008, and references therein for details; see also Suleimanov et al. 2009). For simplicity, we focus our calculations on fully ionized hydrogen atmospheres (see below). In this case, the scale height is
(1) 
where is the effective temperature of the star, is the gravitational acceleration, and is the angle that the ray makes with the surface normal.
In a NS atmosphere, radiation propagates in two distinct polarization modes: the extraordinary (X) mode and ordinary (O) mode, which are very nearly linearly polarized perpendicular and parallel, respectively, to the plane formed by the directions of propagation and the magnetic field. The X mode absorption opacity is reduced by a factor of relative to the O mode opacity, where the electron cyclotron energy is keV and . Thus the X mode radiation decouples from deeper, hotter layers in the NS atmosphere than the O mode radiation. For rays propagating at intermediate angles relative to the magnetic field, the net emission is significantly polarized and dominated by the X mode.
To produce surface emissivities for the X and O photon modes, we use the model of van Adelsberg & Lai (2006), which quantitatively incorporates the effects of vacuum polarization on the radiative transfer for fully ionized atmospheres with an external magnetic field parallel to the axis of symmetry. This model is reasonably accurate for temperatures K. At lower temperatures, partial ionization of atomic species becomes important, and partially ionized models are needed (see Ho et al. 2008, and references therein). Therefore, we use the fully ionized model only in the high temperature regime. It should be noted, however, that this is the effective temperature of the star, and that the measured value at the observer is lower due to gravitational redshift. Colder stars are less interesting for the purpose of this study: First, being much dimmer, their measured spectra are generally of lower statistical significance for detailed spectral studies. Second, and more importantly, since the fraction of Omode photons increases with temperature, spectra of colder stars are less sensitive to axion effects (see the discussion below).
For magnetic field strengths G, there is a significant contribution from vacuum polarization to the dielectric properties of the medium (Adler 1971; Tsai & Erber 1975). A ray traversing the density gradient of a magnetized NS atmosphere will eventually encounter a layer of the medium in which the plasma and vacuum contributions to the dielectric tensor are of the same magnitude, leading to a resonance. For the models studied in this paper, the resonance occurs at lower density than the X and O mode decoupling depths (where the optical depth ) for most photon energies and propagation angles. In addition, the integrated opacity across the vacuum resonance is negligible at these field strengths. A discussion of vacuum resonance effects on the mode opacities is given by Ho & Lai (2003). As shown by Lai & Ho (2002) and Lai & Ho (2003a), there is coherent mixing of the modes at the resonance, analogous to the MikheyevSmirnovWolfenstein (MSW) effect for neutrino oscillations (see, e.g., Bahcall 1989; Haxton 1995). Using the geometric optics approximation, and neglecting damping terms in the dielectric tensor (which only affect the width of the resonance), the amplitudes of the two modes, and , evolve according to
(2) 
where (e.g. Lai & Heyl 2006)
(3)  
(4)  
(5) 
In equations (3)(5), is the angle between the magnetic field and the photon direction (defined to be the direction), is the photon frequency, and , where and are the electron plasma and cyclotron frequencies, respectively. The field dependent functions and are (Potekhin et al. 2004)
(6) 
and
(7) 
where is the fine structure constant and is the magnetic field strength in units of the critical field G. The equations above assume that and . Since and at 1 keV, these constraints are easily satisfied throughout the soft Xray band in which NS thermal spectra typically peak.
Equation (2) determines the evolution of the photon mode amplitudes through the atmosphere (see Ho & Lai 2003; Lai & Ho 2003, for details). The amplitudes are related to the specific intensities of each mode by . As the photons encounter the resonance, the polarization mode ellipticities, , experience a discontinuity at the resonance. As discussed in Ho & Lai (2003), an alternative description of the polarization state uses ellipticities , which vary continuously through the vacuum resonance. If the variation of the atmosphere density profile is gradual enough, an adiabatic condition is satisfied (see below), and the mode properties are fixed along the curves described by . Before the resonance, a photon corresponds to an X mode polarized photon, while a photon corresponds to an O mode photon. After the resonance, the correspondence between and switches, so that under adiabatic conditions the character of each modes changes. This is analogous to the correspondence between flavor and mass eigenstates in the MSW mechanism for neutrino mixing. The adiabatic condition is set by evaluating the coefficient matrix in Eq. (2) at the resonance and setting the magnitude of the diagonal term equal to the offdiagonal term (Lai & Ho, 2002; 2003a). In the nonadiabatic regime, the diagonal terms have a much larger magnitude than the offdiagonal terms, leading to a set of decoupled equations for . In the adiabatic regime, the offdiagonal terms dominate over the diagonal ones, leading to a coupled set of equations for and mixing of the polarization states.
The conversion probability is derived from Eq. (2) by setting the radiation completely in one mode (e.g., , ) and evolving the equations to . If the atmosphere density profile is linear in z, then () is a linear function of z, and the conversion fraction is given by the well known LandauZener formula (see Lai & Ho, 2002; 2003a, and the references therein):
(8) 
The adiabatic energy is defined by
(9) 
where is the ion cyclotron energy, and is a slowly varying function of the magnetic field (see van Adelsberg & Lai, 2006, and the references therein). Because the vacuum resonance is very narrow in energy (and hence density, for a given photon energy and propagation angle), the linear approximation for the density profile is always satisfied for the atmosphere models considered in van Adelsberg & Lai (2006). In addition, the conversion occurs on a distance scale much less than the scale height of the atmosphere; thus the conversion process can be regarded as occurring exactly at the resonance density, and the asymptotic solution to Eq. (2), in the form of the conversion probability Eq. (8), can be used to treat mode conversion effects. This was shown to be an accurate approximation in van Adelsberg & Lai (2006), and we adopt it here in our calculations of mode conversion.
A photon with energy encounters the vacuum resonance at a density
(10) 
For a given energy, is the density at which the vacuum and plasma contributions are of equal magnitude (i.e., when ). If , then conversion between the modes at this depth is effective. If and are the values of the mode specific intensities at the resonance, then after conversion they become:
(11) 
This process is automatically incorporated into the models with G and G ( K) using the code of van Adelsberg & Lai (2006). For models with G, we use the code of Ho & Lai (2001). This code does not include partial mode conversion, but we have included the effect according to the following prescription: The code of Ho & Lai (2001) produces atmosphere profiles and at a set of discrete points in Thomson optical depth, , where . Since the conversion process occurs at lower densities than both the X and O mode decoupling depths, it does not affect the atmosphere structure and can be treated as a modification operating on the emergent specific intensities. For a given photon energy and propagation angle, we calculate the vacuum resonance density using eq. (10), and then locate a grid point , such that . We then compute a linear approximation to the resonance Thomson depth , using the formula:
(12) 
where , and . We relate the resonance temperature to the Thomson depth and density in a hydrostatic atmosphere using
(13) 
and calculate the atmosphere scale height at the vacuum resonance according to Eq. (1). Finally, the adiabatic energy and mode conversion probabilities are calculated according to Eqs. (8) and (9). The emergent intensities are then mixed using Eq. (2).
When the surface magnetic field exceeds the quantum critical field G, there are significant vacuum contributions to the dielectric properties of the medium. At the resonance, if the radiation energy is much greater than the adiabatic energy [Eq. (9)], there is a high probability of conversion between the X and O modes. For magnetic field strengths G, the vacuum resonance lies between the photospheres of the X and O modes (where the mode opacities ). The modes mix at the resonance with a fraction of Xmode photons converting into Omode photons and viceversa according to the probability in Eq. (8). The Omode photons, after converting into Xmode photons at the resonance, decouple from the atmosphere. The Xmode photons that convert to the Omode, subsequently interact strongly with the atmosphere, as they now experience the large opacity for O mode polarization. Any Xmode photons that do not convert encounter a large integrated opacity at the resonance (under typical conditions; see Ho & Lai, 2003). Thus, one of the net effects of vacuum polarization is to shift the location of the Xmode photosphere to the resonance density, which is at lower density and temperature than the original Xmode photosphere. As a result, the high energy spectrum is softer, and the relative fraction of X to O mode photons decreases. At magnetic field strengths G, radiation in both modes decouples from the NS atmosphere before encountering the vacuum resonance. In this case, mode conversion essentially switches the photospheres of the two modes, with the Omode photons now being produced in the deeper (hotter, denser) layers of the star. Therefore, if the emission is mostly polarized in the X mode, for , the plane of polarization switches to the O mode after resonant conversion (Lai & Ho 2003a, 2003b). The above considerations are particularly relevant for the subsequent mixing between Omode photons and axions, as emission which would otherwise have a small fraction of O mode radiation now predominantly consists of O mode photons, depending on the magnetic field strength and NS geometry (see van Adelsberg & Perna 2009, for more discussion).
Figure 1 shows the fractional contribution to the emergent flux of the X and O modes for three values of the field strength, G, G, and G. The effective temperature is K in all cases. Each panel shows the results for a separate emission angle relative to the magnetic field. In general, the X mode dominates the emergent flux at low energies, except around the ion cyclotron feature . According to the formula for [eq. (9)], there is an angledependent point at which mode conversion becomes efficient and the radiation is composed mostly of O mode photons. Since for a magnetic field perpendicular to the surface, at larger values of , mode conversion is efficient only at high energies. Note that as , the properties of the two modes become similar, and the difference in the opacities decreases.
Figure 2 shows the fractional contributions to the emergent flux of the X mode and O modes, for a fixed magnetic field strength, G, at three effective temperatures K, K, and K. The general trend is an increase in the O/X flux ratio as the temperature increases. This is straightforward to understand: since a higher temperature atmosphere produces more photons with higher , the average ratio decreases. Thus, the Xmode opacity becomes more similar to the Omode, leading to an increase in the O/X ratio.
3. Photonaxion conversion in the magnetic field of neutron stars
In the following we discuss the effect of photonaxion mixing as a photon beam propagates through the dipolar field of a magnetized NS. Firstly, we note that an important assumption of our analysis is that the conversions between the O and X modes and between the O mode and axion can be treated separately (each as a twostate mixing) rather than as a coupled threestate mixing. This assumption is valid as long as the vacuum resonance and the photonaxion resonance are well separated. This was shown to be the case by Lai & Heyl (2006) for highly magnetized NSs at soft Xray energies and axion masses eV. Thus, the assumption of decoupled mixing is satisfactory for our regime of interest, and it will be adopted in our analysis.
Our formalism for photonaxion conversion in a magnetized medium follows that of Raffelt & Stodolski (1988). To be consistent with the notation adopted in the previous section for the two photon polarization states, we denote the amplitudes of photons polarized parallel and perpendicular to the magnetic field as and , respectively. Under the assumption that the length scale over which the magnetic field varies is much larger than the photon and axion wavelengths, the evolution of the photon and axion fields in the direction, with frequency , is given by
(14) 
in units with , where is the axion field. The diagonal matrix elements are given by
(15) 
(16) 
where is the axion mass and is defined by Eq. (6). The offdiagonal component is given by
(17) 
where describes the coupling strength between the Omode photons and the axions, and is the distance from the stellar surface. Since mixing occurs only between the axion field and the Omode photon field, we neglect the evolution of the Xmode field in the discussion that follows.
The mixing strength can be measured in terms of the ratio between the offdiagonal term and the difference between the diagonal terms in Eq. (14); a useful parametrization can be made in terms of the mixing angle defined by
(18) 
We begin our analysis by considering a region of space over which the magnetic field can be approximated as homogeneous. Following Raffelt & Stodolsky (1988), the discussion is greatly simplified if we define phases relative to the unmixed component, and neglect a common phase. The solution can be found by first performing a matrix rotation to an eigenstate basis where the propagation matrix is diagonal, propagating the two eigenstates independently, and then rotating back to the photon–axion basis. This yields an evolution equation for the mixing components:
(19) 
where
(20) 
and
(21) 
(22) 
and the plus and minus sign indicates and , respectively.
The solution provided by Eq. (19) assumes that the field is homogeneous. The generalization for photon propagation within the dipolar field around the NS, with , can be made by assuming that the total evolution operator is the product of evolution operators within spatial slices of width wherein the field can be approximated as constant.
Within each shell, we use Eq. (19) to evolve the photon and axion amplitudes from and in shell to and , where is the subsequent shell. We tested the code for convergence and found that, for km, the amplitudes at infinity converge for the range of , , and of interest.
Fig. 3 shows the conversion probability for combinations of the axion parameters in the soft Xray band of interest for NS observations. The probability is plotted as a function of photon energy and propagation angle . There is a range in the parameter space of axion mass, coupling strength, and NS properties (such as surface magnetic field strength) where the photonaxion conversion probability is nonnegligible. This implies that photonaxion conversion can leave distinctive signatures in NS spectra, making them interesting sources for constraining axion physics. In the next section, we perform a detailed analysis of such signatures.
4. Thermal spectra, light curves, and polarization of neutron stars with photonaxion conversion
The observed spectrum of a NS is obtained by summing the local emissivities over the emission region. For simplicity, we consider an emitting region that has a constant temperature and magnetic field strength, with the field direction along the surface normal (e.g., a hot magnetic polar cap). Such a region allows us to directly examine the spectral dependencies of photonaxion conversion on these physical quantities. It is often the case that the emission is dominated by a region smaller than the entire star surface (i.e., a hot spot). The spectrum and lightcurves will then be functions of the viewing angle between the observer and spot axes. As we will show in § 5, phase dependent methods add considerable diagnostic power for detecting the presence of axions.
The calculation of phasedependent emission from an extended region on the NS surface follows the formalism developed by Pechenick et al (1983), with generalizations by Perna & Gotthelf (2008; see also Bernardini et al. 2011, for a similar geometry). We define the timedependent rotational phase as the azimuthal angle subtended by the magnetic dipole vector around the axis of rotation. It is related to the modulus of the NS angular velocity, , by . We choose the coordinate system so that the observer is located along the axis; the inclination angle of the rotation axis, , with respect to the line of sight is denoted by , while the angle between the magnetic dipole vector and the rotation angle is denoted by . Thus, the angle between and the line of sight is given by
(23) 
The emission region is assumed to be circular, with opening angle , and is centered around the magnetic dipole axis. Hence, also indicates the timedependent angular separation between the spot axis and the line of sight as the star rotates. We describe points on the NS surface by means of the polar angle , and the azimuthal angle, , in spherical polar coordinates (see Fig. 1 in Perna & Gotthelf 2008 for a visual representation of this geometry).
The emission region is restricted to for , while for it is identified by the conditions
(24) 
with
(25) 
and by the condition
(26) 
In the latter case, the outer boundary of the spot, , must be determined numerically from the expression
(27) 
Due to the intense gravitational field of the NS, photons emitted at the NS surface are substantially deflected as they travel to the observer. A photon emitted from colatitude will reach the observer (at infinity) if emitted at an angle with respect to the surface normal; the relation between the two angles is given by the raytracing function
(28) 
where , and are the NS radius and mass, respectively, and is the Schwarzchild radius.
The spectrum at the observer is obtained, as a function of the viewing angle, , by integrating the local emission over the observable surface; this procedure yields the flux seen by an observer at a distance . Accounting for the gravitational redshift of the radiation, this integral takes the form (Page 1995 and generalizations by Pavlov & Zavlin 2000, Heyl et al. 2003):
(29) 
in units of photons cm s keV. In the equation above, the index represents any of the Stokes parameter (U, V, Q, I), and is the energy as seen by the distant observer. The energy emitted at the stellar surface is given by , with
(30) 
The local specific intensity is equal to zero outside of the boundaries for and defined by equations (24)(26).
The observed flux depends on the geometric angles and through the angle in Eq. (23). In the following, we consider an orthogonal rotator, for which . This implies that . Given that the emission region is centered around the dipole axis, the angle will be identified with the angle between the photon propagation and magnetic field directions (see §3).
If photonaxion conversion occurs with probability , the phaseresolved spectrum is given by
(31) 
For a fixed energy (or integrated over an energy band), (or ) yields the light curve as a function of the viewing angle . The phase averaged spectrum is readily obtained from
(32) 
Finally, the linear polarization is computed as
(33) 
To simplify our notation, we will omit the subscript from the observed energy. In the discussion that follows, all quoted energies are redshifted energies for a NS of mass and radius km.
Before concluding this section, we need to remark that
Eq. (31) and Eq.(33) assume that the O and X
modes remain decoupled once they have emerged from the atmosphere of
the NS. The evolution of the modes in the changing dipole magnetic
field of the star was computed in papers by Heyl & Shaviv (2000,
2002), Heyl, Shaviv & Lloyd (2003), Lai & Ho (2003) and van
Adelsberg & Lai (2006). All of the above studies found that, in the
birefringent, magnetized vacuum near the star surface, the photon
modes are decoupled. In this region, the mode eigenvectors evolve
adiabatically along the changing direction of the magnetic field.
Further from the star, the mode evolution continues to be adiabatic,
up to the ”polarization limiting radius”, , which is defined
as the distance at which the polarization modes mix (and are frozen
thereafter). Heyl, Shaviv & Lloyd (2003) derive the following
expression for
(34) 
where and . For radii larger than , the mode eigenvectors are fixed, and the polarization no longer evolves (i.e., it is the observed polarization). The polarization modes do couple near the polarization limiting radius. The effect of this coupling was explicitly computed, for example, in van Adelsberg & Lai (2006). However, except for very rapidly rotating stars (i.e., with frequencies of tens of milliseconds) the coupling results in only a very small quantitative change. For example, for a star with period sec ( G, keV), the mode intensities and (, ) couple at the level of , and the amount of mixing further decreases at longer periods. Therefore, for the purpose of our study and the observational tests we propose (requiring NSs with relatively long periods in order to perform timeresolved spectroscopy), the formalism that we use provides a very good description of the observed flux and the linear polarization fraction computed here (see also Lai & Heyl 2006). We further note that, if the radius of the emission region is much smaller than the polarization limiting radius, then summing over the Stokes parameters at the detector is well approximated by summing over a mode (X or O) at the NS surface (Heyl et al. 2003). Since our emission region (a hot spot) is always much smaller than for the field strengths and energies under consideration, we use this approximation in our calculations.
5. Results: theoretical predictions and observational tests
5.1. Theoretical predictions
We present the results of our calculations of NS spectra, light curves, and polarization signals, with and without photonaxion conversion. We explore the dependencies on the stellar magnetic field strength and effective temperature, with particular emphasis on the magnetic field since it has a stronger influence than the temperature on the functional form of the modal energy dependence (see Figs. 1 and 2). We consider emission from a circular region, 10 deg in angular size, that is centered on the magnetic pole of the star. Such hot spot emission is common in NSs, and exhibits a strong dependence on phase. Given the significant angle dependence of the O and X mode intensities (see Figs. 1 and 2), the study of emission from a finite region allows us to extract more information on the effects that photonaxion conversion has on the observable properties of NSs.
As shown in Fig.3, the conversion probability is a strong function of the axion mass, , and the coupling constant, . For the thermal spectra of NSs to be affected by photonaxion conversion, the probability of this process must be nonnegligible in the soft Xray band, where the NS thermal spectrum peaks. Therefore, we perform a detailed analysis of the consequences of photonaxion conversion on NS spectral properties for a particular set of axion parameters; these parameters are chosen to yield a substantial conversion probability in the soft Xray band. We choose the values eV, and GeV; the study of this specific case will uncover the main signatures of axion conversion. Once these are identified, we explore a wider area of parameter space to identify the main region that can be probed by means of NS observations in the Xray band.
Fig. 4 shows phaseresolved spectra for three magnetic field strengths, , , and G. For each value of , we plot the spectra at two phases, coinciding with the angles between the line of sight and the magnetic field axis for the chosen geometry. We show angles deg and a typical large viewing angle deg. The former angle is essentially a ’phaseon’ spectrum (i.e. a spectrum observed in correspondence to the maximum of the pulsation); we do not plot the results for since the conversion probability in that case is zero. For each combination of viewing angle and magnetic field strength, the figure shows the fluxes from the O and Xmodes separately, as well as the total composite spectra. Results are plotted with and without photonaxion conversion.
The relative fraction of O and X mode photons, as a function of energy and angle, plays a fundamental role in determining the effects of axions on phasedependent spectra and light curves. Before discussing the spectral effects of axions in further detail, we comment on the energy and angular dependence of the two photon modes. Although only the Omode photons are affected by axion conversion, the properties of the Xmode photons also matter for our purposes because of vacuum polarizationinduced mode conversion. The Xmode has a complex beam pattern, consisting of the sum of a narrow “pencil” beam and a broad “fan” beam. The width of the pencil beam scales as , decreasing with magnetic field strength (Pavlov et al. 1994). Therefore, the modal intensities vary from being approximately equal to significantly different between widely separated emission angles.
For a phaseon spectrum, the intensities of the O and X modes peak at roughly similar energies. Thus, the presence of photonaxion conversion results in a dimming of the spectrum without significant spectral distortion. The situation is different at larger viewing angles, as evident in the deg panels of the Fig. 4. In these cases, the X mode dominates the spectrum at (redshifted) energies keV, while the O mode dominates at higher energies; thus, an inversion of the dominant spectral mode occurs. Suppression of the O mode by photonaxion conversion therefore produces not only a dimming of the spectrum but also a spectral distortion, suppressing the highenergy tail.
Since axion conversion occurs at larger distances from the star than mode conversion due to the vacuum resonance, only Omode photons are affected as they emerge from the star atmosphere. This is shown in Fig. 4, where the Xmode photons are untouched by photon conversion, while the Omode intensities are noticeably suppressed. The overall effect on the spectrum is that of flux suppression. However, since the O/X flux ratio, as well as the conversion probability, has an energy dependence which varies with viewing angle, the shape of the spectrum is also affected, causing a shift of the energy peak. While the shift is very small at phaseon viewing angles (since the Omodes and Xmodes have a similar energydependent intensity for deg), it is very pronounced at wider viewing angles, where the Omode largely dominates the intensity at high energies. In this case, enhanced suppression of high energy photons produces an effective shift of the peak of the spectrum towards lower energies.
Figures 5 and 6 quantify the phasedependent spectral effects discussed above. For the same values of magnetic fields and axion parameters considered in Fig. 4, Fig. 5 shows the ratio of the Xray flux in the 0.810 keV band photonaxion conversion to that without conversion. As expected, at all phases, the ratio is smaller than 1. This is because the contribution to the flux from the Omode photons is nonnegligible at all viewing angles with respect to the magnetic field axis. The precise value of the flux ratio, however, shows a significant dependence on the viewing angle, with a mimimum around deg. This is due to the Omode flux, which is a large fraction of the total flux at keV, for viewing angles in the deg range.
The left panels of Fig. 6 show the (redshifted) energy
where the spectrum peaks, , as a function of viewing
angle. We consider only angles deg since the flux from
larger angles, which is still visible due to relativistic
light deflection, is too low for phaseresolved spectroscopy.
We begin by examining the behavior of without photonaxion
conversion. At the magnetic field strengths for the models in this paper,
the spectral peak of the total (X and O mode) specific intensity follows
the beam pattern of the X mode: a narrow pencil beam at small angles
and a broad fan beam at large angles.
This is because conversion between the X and O modes occurs above the
photosphere of both modes and thus does not change the total intensity.
The observed variation of with viewing angle deviates somewhat from the
emission beam pattern due to the narrowness of the pencil beam (a few deg in width
around ), the gravitational deflection of photon trajectories,
and the finite size of the emitting region, which smears out features from nearby emitting
points
When axion conversion is taken into account for G, a common feature is a shift of the peak energy to lower values at all viewing angles except for those close to phaseon and edgeon. For deg, the intensity of the O mode photons peaks at similar or slightly lower energies than that of the Xmode photons (see Fig. 4); hence, suppression of O mode photons primarily affects the amplitude of the spectral intensity, as opposed to the location of the spectral peak. In fact, because the Omode intensity peaks at lower energies, especially at lower fields, at small angles there is actually a shift of the peak towards higher energies after photonaxion conversion.
Another notable signature in the variation of is a strong decrease (or flattening in the higher field case) for viewing angles deg. Around these angles, the contribution from the O mode photons is dominant at high energies, and the reduction of these photons through photonaxion conversion results in a significantly cooler spectrum. As the viewing angle increases further, the Xmode photons become increasingly dominant (see Fig. 1), and the reduction of high energy photon lessens. As a result, the effective temperature of the spectrum – as measured through – approaches the value it would have without photonaxion conversion.
To examine the dependence of our results on the effective temperature of the star, we consider the values , , and K (where is the unredshifted temperature at the stellar surface) for the magnetic field G. We find that the shift in due to axion conversion (compared to the case with equal temperature but no conversion) is more pronounced at higher temperatures. This is not surprising given that the relative fraction of O to Xmode photons is larger at higher temperatures.
The right panels of Fig. 6 show the phaseresolved effective area [i.e., ] that an observer measures when fitting spectra with a blackbody model. To eliminate the dependence of this diagnostic on the often poorly determined distance to the NS, we normalize to its value near phaseon. Since the photonaxion conversion probability is zero for precisely , we normalize to a small angle deg which represents a “phaseon” spectrum as an average of small angles around deg. When photonaxion conversion is included, the variation of the effective area with phase results from the interplay of two effects: a decrease in due to flux dimming and an increase in due to decreases in the peak energy and inferred temperature. Because of the dependence, the latter effect dominates the former. This is particularly the case for lower fields, which has the strongest shift with phase.
The variation of inferred effective area with phase is a particularly interesting and peculiar feature of photonaxion conversion. The specific details of the variation, however, depend on NS and axion parameters. In order to explore a wider region of parameter space, we compute the effective area for the phaseaveraged spectrum () with photon axion conversion, normalized to that without axion conversion. This represents a measure of the “correction” that photonaxion conversion makes to the inferred effective area (as determined through measurements of phaseaveraged flux and effective temperature). For our purposes, the phaseaveraged spectrum of a hot spot provides an adequate approximation for emission from the entire NS surface. The results of this calculation are displayed in Fig. 7 for a combination of magnetic fields and effective temperatures: G, K; G, K; and G, K. It is interesting to note that, in all cases, there are regions of parameter space in which the inferred effective areas are smaller/larger than the uncorrected areas. This is expected given the two counteracting effects described above. The region shown in Fig. 7 is one for which, at typical NS magnetic fields, there is a larger probability of photonaxion conversion in the soft Xray band (where NSs are routinely observed). The correction to the effective area becomes negligible for the largest values of and displayed in Fig. 7. This is expected since the photonaxion conversion probability tends to zero for these values at the keV energies of interest here. The results summarized by this figure have a strong potential to constrain axion signatures for NSs with good distance measurements (see the discussion in §5.2).
In Fig. 8, we examine the effects of photonaxion conversion on NS light curves at different energies, for the same set of axion parameters as in Fig. 4. The lightcurves are all normalized to the maximum value of the flux among the set of energies considered (0.5, 1.5, 2.5, and 4 keV). The main effect of photonaxion conversion is the relative suppression of the flux at higher energies with respect to that at lower energies.
Finally, we study the effect of photonaxion conversion on the magnitude of linear polarization [see Eq. (33)]. For consistency with Fig. 4, we consider the same set of parameters. The polarization signal with and without conversion is shown in Fig. 9. Since the main effect of axion conversion is suppression of Omode photons, the net result is a reduction of the magnitude of linear polarization. In addition, it is noteworthy that the presence of photonaxion conversion causes the plane of polarization to rotate from parallel (with respect to the plane formed by the propagation direction and magnetic field) to perpendicular for a phaseon spectrum (small , left panels of Fig. 9). This can be understood with reference to Fig. 4: for a NS observed phaseon with respect to the emission region, the intensity of the O mode photons slightly exceeds the Xmode photon intensity. When photonaxion conversion is considered, the O mode photon intensity is suppressed and becomes smaller than the Xmode intensity, which results in an inversion of the plane of polarization with respect to the case without conversion. Clearly, the extent to which this occurs depends on the specific axion parameter values, since the magnitude of the conversion probability is a strong function of and (see Fig. 3). For example, at G and K, the Xmode intensity of a phaseon spectrum is about 90% of the Omode intensity at keV. Therefore, if photonaxion conversion occurs with probability , it results in an inversion of the plane of polarization. For the combination of axion parameters adopted in Fig. 3, a probability of is achieved when GeV (for eV). When GeV, the probability is for the entire range of axion masses considered ( eV).
At larger viewing angles on the other hand (right panels of Fig. 9), it becomes evident one of the unique signatures of vacuum polarization on the polarization signal of NSs with magnetic fields G, that is the rotation at high energies of the plane of polarization from perpendicular to parallel (Lai & Ho 2003b). Photonaxion conversion does not erase this feature; however, as a result of the suppression of the Omode photons which dominate at higher energies, the shift of the plane occurs at higher energies than it would without conversion. For the axion parameters adopted in Fig. 9, the energy at which the shift occurs with the conversion differs by a factor of with respect to the case with no conversion, with the precise value depending on the magnetic field strength.
5.2. Observational tests
Our theoretical results provide several methods for constraining axion parameters from actual measurements of NS spectra. The foremost method is spectral analysis, often performed on phaseaveraged spectra when the resolution low. We showed in §5.1 that the relative suppression of the Omode photons with respect to the Xmode photons noticeably distorts the spectral shape with respect to the case without conversion. This distortion is strongest in the highenergy spectral tail, where the intensity of the Omode photons is larger for most viewing angles (see Fig. 4). Observed enhancements of high energy spectral tails are generally modeled with a power law (often of unknown origin; e.g. Manzali et al. 2007, Durant et al. 2011 for young and middleaged pulsars; Perna et al. 2001, Juett al. (2002), Rea et al. 2007 for magnetars). However, a deficit of high energy photons in a thermal fit could be a signature of photonaxion conversion, and the axion parameters could then be determined by means of detailed spectral fits that account for the effect. For this type of analysis, the distance to the NS is not important since the main discriminant is the spectral shape, in particular the deficit of highenergy photons.
For sources with good distance estimates
To illustrate how to use Fig. 7 with NS observations,
suppose is the maximum radius allowed by a particular EOS
(see, e.g., Lattimer & Prakash 2007, for a review), and
is the corresponding surface area. Then a NS whose inferred emission
area is would rule out the region of the
parameter space for which the effective area ’correction’ due to
photonaxion conversion is larger than 1. For G
and K (corresponding to a redshifted value K), the parameter space with eV and GeV (or GeV) would be ruled out. For NSs whose
radius determination relies on methods other than flux and peak energy
measurements
As an aside, we note that the the effects of axions on the inferred area discussed above suggest that the use of NS effective radii to rule out nuclear EOSs may not be entirely reliable. The existence of axions could severely bias such conclusions.
For sources with phaseresolved Xray spectra, observations can be fit with either blackbody or model atmosphere spectra to obtain peak energies and effective areas as a function of rotation phase or viewing angle . Figure 6 shows a clear signature of photonaxion conversion, i.e., a significant rise and fall of effective area at intermediate viewing angles, especially for . The important quantity for this diagnostic is the relative change of the effective area with phase; thus, this method does not require knowledge of the distance to the NS. It is also rather insensitive to the NS radius/mass (i.e., surface gravity), since variations in this quantity produce energy shifts that are the same at all phases, leaving the relative phasedependent signatures of the conversion unaltered. An exact estimate of the range of axion parameters that can be constrained by means of this method cannot be made in absolute terms, since it clearly depends on the quality of the data. However, from a theoretical point of view, photonaxion signatures are imprinted in the NS spectra as long as the conversion probability is nonnegligible in the typical Xray band of observations. Fig. 4 shows that this is the case for an axion coupling strength GeV and for a large range of axion masses, up to eV if GeV.
The theoretical framework and methods for extracting axion parameters
from observations of magnetic NSs is presented above. Detailed
analyses of specific sources require a model
appropriate for the specific viewing and emission geometry, magnetic
field strength, and surface temperature for each source;
a comparison with the data can
then be performed by folding the model through the response function of the
detector. Such an analysis will be the subject of
future work. Below, we briefly review NSs that serve as potential
candidates for obtaining constraints on axion physics. These sources
possess bright Xray thermal emission from the NS
surface. Among the observed NS population, such sources can be
grouped into several classes (which may not be mutually exclusive
groups; see, e.g., Kaspi 2010; Perna & Pons 2011; Pons &
Perna 2011). One class is the magnetars, which generally possess
superstrong magnetic fields
.
The isolated NSs are another interesting class. The seven confirmed isolated NSs have inferred from timing measurements (Kaplan & van Kerkwijk 2009, 2011), while features in their Xray spectrum suggest if the features are due to the proton cyclotron resonance (Haberl 2007). We note that the observed Xray spectra of isolated NSs are generally wellfit (except around spectral lines) by blackbodies, including Wienlike behavior at high energies (van Kerkwijk & Kaplan 2007). Various explanations have been proposed to explain the soft high energy tails as compared to the hard tails predicted by atmosphere models (e.g., Romani 1987), such as emission from a partially optically thin atmosphere (e.g., Motch et al. 2003) or the effect of vacuum polarization (see §2 and Ho & Lai 2003). Here we have shown that photonaxion conversion can suppress highenergy emission and produce softer tails (see Fig. 4). However, this cannot explain the spectra of isolated NSs. These NSs have surface temperatures K, and, as we discussed in Sec. 2, the conversion is less effective at low temperatures due to the low fraction of Omode photons. However, if soft tails (softer than the blackbody) were detected in NSs with relatively high surface temperature, then photonaxion conversion could indeed be a competing explanation.
A final class is the central compact object neutron stars (CCOs), which have , obtained from timing and spectral analyses (see, e.g., Halpern & Gotthelf 2010). At these low fields, the vacuum resonance is outside of the atmosphere, and the emerging intensity has a high fraction of Omode photons; their surface temperatures are constrained to be K, while hot spot temperatures are quite high [ K]. These characteristics, combined with the fact that the photonaxion conversion probability at these lower fields is still high for a wide range of axion parameters (Jimenez et al. 2011), make CCOs potentially good candidates. Current data on many of the NSs belonging to the classes described above are already of sufficient quality to allow axion constraints to be derived.
Finally, we note that Xray polarization measurements, like those that will be performed by the forthcoming Gravity and Extreme Magnetism Small Explorer (GEMS) satellite (Jahoda 2010), will increase the robustness of axion constraints, as well as provide an independent measure of the stellar magnetic field strength and geometry (van Adelsberg & Perna 2009).
6. Summary
Constraints on the axion mass and coupling strength by means of observations of magnetized objects have been discussed in the recent literature (e.g., Lai & Heyl 2006; Chelouche et al. 2009; Gill & Heyl 2011; Jimenez et al. 2011) as a means to complement and independently test constraints obtained from other methods. The studies above emphasized different observational aspects and tests. The modification to the Xray spectra of highly magnetized NSs induced by photonaxion conversion was discussed by Lai & Heyl (2006) using a blackbody spectrum to model the photon spectrum from NSs. Chelouche et al. (2009) discussed spectral features produced in the submm/IR wavelength regime for magnetized NSs, while Gill & Heyl (2011) considered limits that can be derived from polarization measurements of white dwarfs. Jimenez et al. (2011) focused their analysis on constraints that can be obtained through observations of eclipsing white dwarfs in binaries, but they also discussed qualitatively the potential of NSs. Since most of the photonaxion conversion happens at a distance of many stellar radii from the NS, they pointed out that there are two configurations that should be considered: light from a background object passing through the magnetic field of the NS (i.e., occultation) and a binary containing a NS where the companion transits close enough to the line of sight for its light to be influenced by the NS magnetic field. In these two cases the treatment is much simpler as the physics of photon propagation in the NS atmosphere is irrelevant. However, they determined that the occultation probability is too low to be astrophysically relevant and that there are no known binary systems involving a NS that are detached enough to yield a clean constraint.
Since NSs are routinely observed in the Xray domain, the goal of our paper has been an extensive exploration of the actual constraints on axion physics that can be obtained from observations of NS thermal spectra. We have also explored what can be learned from future Xray polarization measurements. For these purposes, we used detailed, magnetized atmosphere models, which properly account for the energy and angledependent emerging intensities of the two polarization modes, including mode conversion due to vacuum polarization effects. The emergent O mode intensity was coupled with the energy and angledependent photonaxion conversion probability, and this allowed us to make theoretical predictions for phaseresolved spectra with photon axionconversion in the relevant regime of magnetic fields and temperatures for NSs. Since the relative intensities of the O and X mode photons depend strongly on the angle at which they emerge with respect to the magnetic field direction, we considered emission from a region consisting of a hot spot with an axis coincident with the NS magnetic field axis. Thermal emission from most Xray bright NSs does indeed indicate that we are observing hot spots (this from both nonnegligible pulse fractions and small inferred emission areas). In order to extend axion studies by means of comparison with observations of any NS (especially those emitting from the entire surface, with nondipolar, complicated B topologies, and fast rotators), one would have to include, locally, magnetic fields nonnormal to the surface in the atmospheres, investigate how nondipolar fields influence both the emission spectra as well as the photonaxion conversion, and, for the polarization, explicitly integrate the Stokes parameters to . However, since the magnetic field structure over the entire surface of a NS is not fully known apriori, determination of the axion parameters would be more degenerate under these more general conditions.
Our analysis has identified some features in NS spectra that bear the telltale signs of photonaxion conversion:

The NS spectral shape is noticeably distorted compared to the case without photonaxion conversion, with suppression of the highenergy tail for viewing angles deg (with respect to the center of the emission region). Detailed spectral fits can yield axion parameters. This analysis does not require knowledge of the distance to the NS, as the axion signatures are imprinted in the spectral shape.

The spectral suppression of the Omode photons by photonaxion conversion dominates the high energy tail of the spectrum for a range of viewing angles. In addition, it results in a shift of the peak energy, , towards lower values for a range of rotation phases (i.e., viewing angles). As a result, the effective area of the spot (as inferred from measurements of the flux and peak spectral energy) is larger than its value without axion conversion. We demonstrated that a clear signature of photonaxion conversion is a significant rise and fall of effective area (normalized to the phaseon value) at intermediate viewing angles, especially for . Such a measurement requires high signaltonoise phaseresolved spectra, but does not require knowledge of the distance to the NS, since the axion signature appears in the relative change of the effective area with phase.

For a star emitting from the entire surface, photonaxion conversion can substantially affect the inferred NS emission area by , and the inferred radius by . These values are measured using the flux and peak energy of the thermal spectrum. If the distance is well constrained, and the NS radius is known through different methods, then the inferred emission area can be directly translated into a constraint in the plane. If the NS radius is not independently known, then an inferred emission area that exceeds the maximum value allowed by NS equations of state could indicate the presence of photonaxion conversion.

In the absence of photonaxion conversion, radiation from a hot spot observed phaseon at energies keV is linearly polarized in the plane formed by and the direction of photon propagation. Conversion rotates the plane of polarization and leads to radiation polarized perpendicular to the plane of and the direction of photon propagation, for a range of axion parameters.
We concluded §5.2 with a discussion of appropriate sources to study within the theoretical framework developed here for probing axions. While this paper has outlined the main elements for connecting theoretical ideas of photonaxion conversion to actual observations that are made of NSs, studying each suitable source will require a specific suite of models, tailored to the particular object (e.g., its magnetic field, surface temperature, emission and viewing geometry). The theoretical models are then convolved with the detector response, and a process of minimization identifies the most likely values of the parameters. Detailed analysis of the most promising NS candidates will be the subject of future work.
Footnotes
 affiliation: JILA and Department of Astrophysical and Planetary Science, University of Colorado at Boulder, 440 UCB, Boulder, CO, 80304, USA
 affiliation: School of Mathematics, University of Southampton, Southampton, SO17 1BJ, UK
 affiliation: ICREA and ICC, University of Barcelona (IEECUB), Spain
 affiliation: Center for Relativistic Astrophysics and School of Physics Georgia Institute of Technology, Atlanta, Georgia 30332, USA
 affiliation: ICREA and ICC, University of Barcelona (IEECUB), Spain
 A similar expression for the polarization limiting radius is derived in van Adelsberg & Lai (2006): . The scaling of in this equation differs slightly from that in Eq. (34) due to different adiabatic conditions applied by the two sets of authors. Heyl et al. (2003) constrain the change in the relative difference in the mode indices of refraction compared to the magnitude of that difference [see their Eq. (9)], while van Adelsberg & Lai (2006) consider the distance at which the mode evolution equations have equal diagonal and offdiagonal terms [see their Eq. (62)]. This results in slightly different magnitudes and scaling relations; for slowly rotating NSs, the precise numerical value of does not change the value of the observables, and both approaches result in the same physics.
 To be able to clearly observe the very narrow pencil beam, both the size of the emission region and the viewing angle would need to be smaller than (or comparable to) the width of the beam.
 Distances to isolated NSs are usually accurate to, at best, , while the distances to NSs in globular clusters are known to much higher precision.
 Examples of such methods of NS radius measurement (or more generally of the radiustomass ratio) include spectral lines (e.g. Cackett et al. 2009), thermal Xray pulse profiles (e.g. Bogdanov et al. 2008; note that at low energies keV, these pulse profiles are not very sensitive to the presence of photon axion conversion, as shown in Fig. 8), gravitational wave emission (Lenzi et al. 2009), quasiperiodic oscillations in accreting NSs (Miller et al. 1998), and neutrino emission from protostars (Lattimer & Prakash 2007). Masses can be independently inferred from radial velocity studies for NSs in binary systems (e.g. van Kerkwijk et al. 2011) and pulsar timing (a comprehensive review of all these methods can be found in Lattimer & Prakash 2007).
 Note however that the bursting behaviour used to characterize magnetars has also been seen in a lower magnetic field object (see Rea et al 2010).
References
 Adler, S. L. 1971, Ann. Phys., 67, 599
 Alcock, C. & Illarionov, A. 1980, ApJ, 235, 534
 Andriamonje, S. et al. 2007, J. Cosm. Astropart. Phys., 4, 10
 Arvanitaki, A., Dimopoulos, S., Dubovsky, S., Kaloper, N., MarchRussell, J. 2010, Phys. Rev., D, 81, 123530.
 Asztalos, S. J. et al., 2004, Phys. Rev. D, 69, 011101
 Asztalos, S. J. et al., 2010, Phys. Rev. Lett. 104, 041301
 Avgoustidis, A. Burrage, C. Redondo, J., Verde, L., Jimenez, R. 2010, JCAP, 1010, 24, preprint arXiv:1004.2053
 Bahcall, J. 1989, Neutrino Astrophysics, (Cambridge: Cambridge Univ. Press)
 Bernardini, F., Perna, R., Gotthelf, E., Israel, G. Rea, N., Stella, L. 2011, MNRAS, 418, 638
 Bogdanov, S., Grindlay, J. E., Rybicki, G. B. 2008, ApJ, 689, 407
 Brown, E. F., Bildsten, L., & Chang, P. 2002, ApJ, 574, 920
 Cackett, E. M., Miller, J. M., Bhattacharyya, S., Grindlay, J. E., Homan, J., van der Klis, M., Miller, M. C., Strohmayer, T. E., Wijnands, R. 2008, ApJ, 674, 415
 Chang, P., Bildsten, L., & Arras, P. 2010, ApJ, 723, 719
 Chelouche, D., Rabadan, R., Pavlov, S. S., Castejon, F. 2009, ApJ, 180, 1
 De Panlis, S. et al. 1987, Phys. Rev. Lett. 59, 839
 Durant, M., Kargaltsev, O., Pavlov, G. G., 2011, arXiv eprint (arXiv:1109.1984)
 Eidelman, S. et al. 2004, Phys. Lett. B592, 1
 Gill, R., & Heyl, J. S. 2011, arXiv:1105.2083
 Haberl, F. 2007, Astrophys. and Space Science, 308, 181
 Halpern, J. P., Gotthelf, E. V. 2010, ApJ, 709, 436
 Haxton, W. C. 1995, ARA&A, 33, 459
 Heyl, J. S., Shaviv, N. J. 2000, MNRAS, 311, 555
 Heyl, J. S., Shaviv, N. J. 2002, Phys. Rev.D, 66, 023002
 Heyl, J. S., Shaviv, N. J., Lloyd, D. 2003, MNRAS, 342, 134
 Ho, W. C. G. & Lai, D. 2001, MNRAS, 327, 1081
 Ho, W. C. G. & Lai, D. 2003, MNRAS, 588, 962
 Ho, W. C. G., Potekhin, A. Y., Chabrier, G. 2008, ApJS, 178, 102
 Jahoda, K. 2010, Proc. SPIE, 7732, 77320W
 Juett, A. M., Marshall, H. L., Chakrabarty, D., Schulz, N. S 2002, ApJ, 568L, 31
 Jimenez, R., PeñaGaray, C., & Verde, L. 2011, Physics Letters B, 703, 232
 Kaplan, D. L., van Kerkwijk, M. H. 2009, ApJ, 705, 798
 Kaplan, D. L., van Kerkwijk, M. H. 2011, ApJ, 740, 30
 Kaspi, V. M. 2010, Publ. of the Nat. Acad. of Science, 107, 7147
 Lai, D. & Heyl, J. 2006, Phys. Rev. D, 2006, 74, 13003
 Lai, D. & Ho, W. C. G. 2002, ApJ, 566, 373
 Lai, D. & Ho, W. C. G. 2003a, ApJ, 588, 962
 Lai, D. & Ho, W. C. G. 2003b, Phys. Rev. Lett., 91, 071101
 Lattimer, J. M., & Prakash, M. 2007, Phys. Rep., 442, 1
 Lenzi, C. H., Malheiro, M., Marinho, R. M., Marranghello, G. F., Providencia, C. 2009, Journ. of Phys. 154, 1
 Manzali, A., De Luca, A., Caraveo, P. A. 2007, ApJ, 669, 570
 Mereghetti, S. 2007, MmSAI, 78, 644
 Miller, M. C., Lamb, F.K., Psaltis, D. 1998, ApJ, 508, 791
 Mori, K. & Ho, W. C. G. 2007, MNRAS, 377, 905
 Motch, C., Zavlin, V. E., Haberl, F. 2003, A&A, 408, 323
 Page, D. 1995, ApJ, 442, 273
 Pavlov, G. G., Sanwal, D., Garmire, G. P., Zavlin, V. E. 2002, Neutron Stars in Supernova Remnants, ASP Conference Series, Vol. 271, held in Boston, MA, USA, 1417 August 2001. Ed. Patrick O. Slane and Bryan M. Gaensler. San Francisco: ASP, p.247
 Pavlov, G. G., Shibanov, Yu. A., Ventura, J., Zavlin, V. E. 1994, A&A, 289, 837
 Pavlov, G. G. & Zavlin, V. E. 2000, ApJ, 529, 1011
 Peccei, R. D., Quinn, H. R. 1977, Phys. Rev. Lett., 38, 1440
 Pechenick, K. R., Ftaclas, C., Cohen, J. M. 1983, ApJ, 274, 846
 Perna, R., Heyl, J. S., Hernquist, L. E., Juett, A. M., Chakrabarty, D. 2001, ApJ, 557, 18
 Perna, R. & Gotthelf, E. V. 2008, ApJ, 681, 522
 Perna, R. & Pons, J. A. 2011, ApJ, 727, L51
 Pons, J. A. & Perna, R. 2011, ApJ, 741, 123
 Potekhin, A. Y., Lai, D., Chabrier, G., Ho, Wynn C. G. 2004, ApJ, 612, 1034
 Preskill, J. Wise, M. B., Wilczek, F. 1983, Phys. Lett. B, 120, 127
 Raffelt, G. G. 2008, Lect. Notes Phys., 741, 51
 Raffelt, G. & Stodolski, L. 1988, Phys Rev. D, 37, 1237
 Rea, N. et al. 2007, MNRAS, 381, 293
 Rea, N. et al. 2010, Science, 330, 944
 Romani, R. W. 1987, ApJ, 313, 718
 Suleimanov, V., Potekhin, A. Y., & Werner, K. 2009, A&A, 500, 891
 Tsai, W. Y. & Erber, T. 1975, Phys. Rev. D, 12, 1132
 van Adelsberg, M. & Lai, D. 2006, MNRAS, 373, 1495
 van Adelsberg, M. & Perna, R. 2009, MNRAS, 399, 1523
 van Kerkwijk, M. H. & Kaplan, D. L. 2007, Ap&SS, 308, 191
 van Kerkwijk, M. H., Breton, R. P., Kulkarni, S. R. 2011, ApJ, 728, 95
 Wuensch, W. U. et al., 1989, Phys. Rev., D, 40, 3153
 Zavlin, V. E. 2009, in Neutron Stars and Pulsars, Astrophysics and Space Science Library, Vol. 357. Springer Berlin Heidelberg, p. 181