Effects of Compton scattering on the neutron star radius constraints in rotation-powered millisecond pulsars

TL;DR

This study evaluates the impact of using the Thomson approximation versus exact Compton scattering in neutron star atmosphere models, finding minimal bias in radius constraints but some bias in hot spot size and temperature estimates.

Contribution

It compares the effects of Compton and Thomson scattering models on neutron star parameter constraints, demonstrating the adequacy of the simpler Thomson model for radius measurements.

Findings

Radius constraints are not significantly affected by the scattering model.

Biases in hot spot size and temperature estimates occur when using the Thomson approximation.

The Thomson model remains adequate for radius measurements with NICER data.

Abstract

The aim of this work is to study the possible effects and biases on the radius constraints for rotation-powered millisecond pulsars when using Thomson approximation to describe electron scattering in the atmosphere models, instead of using exact formulation for Compton scattering. We compare the differences between the two models in the energy spectrum and angular distribution of the emitted radiation. We also analyse a self-generated synthetic phase-resolved energy spectrum, based on Compton atmosphere and the most X-ray luminous rotation-powered millisecond pulsars observed by the Neutron star Interior Composition ExploreR (NICER). We derive constraints for the neutron star parameters using both the Compton and Thomson models. The results show that the method works by reproducing the correct parameters with the Compton model. However, biases are found in size and the temperature of the emitting hot spot, when using the Thomson model. The constraints on the radius are still not significantly changed, and therefore the Thomson model seems to be adequate if we are interested only in the radius measurements using NICER.

Tables

Table 1. Table 1: Most probable values and 68% and 95% credible limits for Compton and Thomson models applied to synthetic data.

Quantity	95% lower limit	68% lower limit	Most probable value	68% upper limit	95% upper limit
Compton fit to Compton model
$R_{\mathrm{eq}}$ (km)	$11.99$	$12.00$	$12.01$	$12.02$	$12.03$
$\theta$ ($\mathrm{deg}$)	$35.6$	$35.8$	$35.9$	$36.0$	$36.2$
$\rho$ ($\mathrm{deg}$)	$4.99$	$5.00$	$5.00$	$5.01$	$5.01$
$T$ (MK)	$3.121$	$3.125$	$3.128$	$3.131$	$3.135$
$\log \left({\sigma }_i\right)$	$0.868$	$0.869$	$0.869$	$0.870$	$0.872$
McPHAC fit to Compton model
$R_{\mathrm{eq}}$ (km)	$11.99$	$12.00$	$12.01$	$12.02$	$12.03$
$\theta$ ($\mathrm{deg}$)	$35.5$	$35.6$	$35.8$	$35.9$	$36.0$
$\rho$ ($\mathrm{deg}$)	$4.95$	$4.96$	$4.96$	$4.96$	$4.97$
$T$ (MK)	$3.136$	$3.139$	$3.142$	$3.145$	$3.148$
$\log \left({\sigma }_i\right)$	$0.980$	$0.993$	$1.007$	$1.020$	$1.032$

Full text

1 Introduction

The equation of state of cold matter beyond nuclear densities can be constrained using astronomical observations of masses and radii of neutron stars (NSs) (Steiner et al., 2010; Lattimer, 2012; Nättilä et al., 2016; Özel & Freire, 2016; Suleimanov et al., 2016; Watts et al., 2016; Degenaar & Suleimanov, 2018; Watts et al., 2019). In case of rapidly rotating NSs having radiating ‘hotspots’ around their magnetic poles, we can model the observed pulses using general relativity and obtain constraints for their mass and the radius (Poutanen & Gierliński, 2003; Miller & Lamb, 2015). However, detailed modelling requires knowledge of the spectral energy distribution and of the angular emission pattern of radiation emitted by the hotspots. The radiation escaping the hotspots is affected by energy-dependent absorption as well as by the anisotropic and energy-dependent scattering of photons by electrons in the atmospheres of NSs.

There have been several studies aiming to constrain NS masses and radii using pulse profiles of accreting millisecond pulsars (AMPs), in which the matter from a low-mass companion star accretes onto the magnetic poles of the NS (see e.g. Poutanen & Gierliński, 2003; Leahy et al., 2008; Poutanen, 2008; Morsink & Leahy, 2011; Salmi et al., 2018). However, these approaches suffer from a relatively high number of unknown NS parameters and from the uncertainties in the atmospheric structure, and therefore also in the angular and energy distribution of the emitted radiation.

In case of rotation-powered millisecond pulsars (RMPs), more independent information of the model parameters (e.g. mass and inclination) is often attained from radio data and the existing NS atmospheric models without effects of accretion may be used. In many RMPs, the bulk of X-ray radiation is thermal emission coming from the polar caps that are heated by a return flow of relativistic electrons and positrons in the open field line region (see e.g. Harding & Muslimov, 2002; Bogdanov, 2018). Nevertheless, few RMPs exhibit nearly pure, non-thermal emission generated most probably by synchrotron emission from pulsar magnetospheres (Zavlin, 2007). We focus on the thermally emitting RMPs, where the composition of the atmosphere is more confidently known than in AMPs (the RMP atmosphere likely consists of pure hydrogen instead of a mixture with heavier elements), and the temperature is low enough that the electron scattering presumably can be described using Thomson scattering approximation. The angular and energy distribution of the escaping photons can be described by using, for example, a plane-parallel atmosphere model in local thermodynamic equilibrium.

This type of model for RMPs that assumes Thomson scattering has previously been implemented in the McGill Planar Hydrogen Atmosphere Code (McPHAC), as described by Haakonsen et al. (2012) (see also e.g. Zavlin et al., 1996; Heinke et al., 2006). This code was also used by Miller (2016) to simulate the data for RMP PSR J1614$-$2230 that can be provided by the Neutron star Interior Composition ExploreR (NICER), and to study the constraints on the NS mass and radius that can be obtained with those data. The question we ask in this paper is how an approximate treatment of Compton scattering affects the radiation spectra escaping from NS atmosphere and how this in turn affects the constraints on NS mass and radius from the NICER data. We note that exact treatment of Compton scattering is very important when considering NS atmospheres heated by accretion (Suleimanov et al., 2018), or by magnetospheric return currents as, for example, was recently discussed by Bauböck et al. (2019) using a very simplified atmosphere model (see also e.g. Zel’dovich & Shakura, 1969; Alme & Wilson, 1973; Ruderman & Sutherland, 1975; Zampieri et al., 1995; González-Caniulef et al., 2019). Modelling the heated RMP atmospheres is, however, beyond the scope of this work and will be discussed elsewhere.

The remainder of this paper is structured as follows. In Sect. 2, we discuss the methods, including modelling the NS atmosphere, ray-tracing, and our method to create and analyse synthetic data. In Sect. 3, we first compare our spectral results to those computed with McPHAC, and then obtain NS parameter constraints fitting the data, which are created with the full Compton model, with both the full Compton and approximate Thomson scattering models. We conclude in Sect. 4.

2 Methods

We first constructed a model for NS atmosphere consisting of pure hydrogen. This is justified by the fact that, without the effects of continuing accretion, gravitational stratification leaves only the lightest elements in the atmospheric layers, which determine the properties of the escaping radiation. We computed the atmosphere model and the angular distribution of the specific intensity of the escaping radiation using three different approaches. In the first one we used our code (Suleimanov et al., 2012), which treats Compton scattering using the exact relativistic Klein-Nishina cross-section and redistribution function derived and presented in details by Aharonian & Atoyan (1981), Prasad et al. (1986), Nagirner & Poutanen (1994), Poutanen & Svensson (1996), and Poutanen & Vurm (2010). As a second model, we used the same code, except where Compton scattering is treated in the Thomson limit. This simplifies and accelerates the calculations dramatically. The third model was constructed using McPHAC code, which also treats Compton scattering in the Thomson approximation (we used their anisotropic version of the model).

The parameters of the model are the effective temperature $T_{\mathrm{eff}}$ (which we will call just $T$ for brevity) and the surface gravity $g$. The solution of the equations that describe the NS atmosphere (see e.g. Suleimanov et al., 2012) provides us with the intensity of the escaping radiation. We tabulated these intensities over a grid of 360 photon energies, equally spaced in $\log E$ (keV) from $-$3.4 to 1.3, and 7 points in the cosine of the zenith angle $\mu$ (in the interval between zero and one using Gaussian nodes) for 11 values of temperature $T$ (K) spaced equally in $\log T$ from 5.5 to 6.6, and ten values of surface gravity $g$ (cm s*-2*) spaced equally in $\log g$ from 13.7 to 14.6.

The observed spectra depend on the NS mass, equatorial radius $R_{\mathrm{eq}},$ and spin (which determine gravitational acceleration $g$ as a function of co-latitude), and on the properties of the emitting spot, that is the local temperature $T$, the angular radius $\rho$, and the centroid (magnetic) co-latitude $\theta$. The spectra also depend on the observer inclination (i.e. the angle between the line-of-sight and the NS rotation axis) and the distance to the source. To compute the observed phase-resolved spectra and pulse profiles, we used ‘oblate Schwarzschild’ approximation (see e.g. Poutanen & Beloborodov, 2006; Morsink et al., 2007; Miller & Lamb, 2015; Salmi et al., 2018), taking into account the deformed shape of the star together with the special and general relativistic corrections to the photon trajectories and angles. For calculations of the total observed flux, integration over the spot surface is needed. However, in order to speed up the computations, the surface gravity $g$ for the atmospheric model was assumed to be constant within the spot (using the correct value for the spot centre). Thus, for each model we first needed only one piece-wise, two-dimensional, linear interpolation from the set of pre-calculated spectral tables to obtain a single two-dimensional array of intensities as functions of $E$ and $\mu$ only (corresponding to a given temperature $T$ and surface gravity $\log g$). In our examples this is justified, because we considered only relatively small spots, where the changes in the NS radius within the spot are small. Then, separately for each position within the spot, we again made a piece-wise, two-dimensional, linear interpolation to obtain intensities corresponding to a required photon emission zenith angle and energy.

The synthetic data were created keeping in mind the most promising NICER targets. Instead of PSR J1614$-$2230 used by Miller (2016) as an example case, we focus on PSR J0437$-$4715 (the closest known RMP), or a similar pulsar with an expected high count rate (needed in order to observe any possible differences in the parameter constraints from the two spectral models). This pulsar has a complicated pulse profile presumably produced by two small high-temperature spots surrounded by a cooler annular region, and also an additional power-law component (Bogdanov, 2013). However, since we aim only to compare the Thomson and Compton models, and are not necessarily interested in modelling this particular pulsar, we ignored these complications, and assumed two spots with constant temperature and pure thermal spectrum. PSR J0437$-$4715 was mainly used to obtain typical values for the parameters of the synthetic data. The model parameters are the following: spot temperature $T\approx$ 3.133 MK (0.27 keV), spot angular radius $\rho=5.0^{\circ}$, spot co-latitude $\theta=36^{\circ}$, equatorial radius of the star $R_{\mathrm{eq}}=12$ km, and an arbitrary phase shift. These parameters were treated as free when fitting the data. Other model parameters were the NS mass $M=1.76\leavevmode\nobreak\ {M}_{\sun}$ , NS spin frequency $\nu=173.6$ Hz, the distance to the star $D=156.3$ pc, the inclination $i=42.4^{\circ}$, and neutral hydrogen column density for interstellar absorption $N_{\mathrm{H}}=7\times 10^{19}\mathrm{cm}^{-2}$ (see e.g. Bogdanov, 2013; Deller et al., 2008; Verbiest et al., 2008). They were regarded as fixed because they are or can be determined from other (radio) observations with relatively good accuracy. We assumed that the observation of the source is long enough to accumulate the total number of observed counts of $4\times 10^{7}$.

The fitting procedure of the data is mostly the same as presented in Salmi et al. (2018). We used Bayesian analysis and an affine invariant ensemble sampler (Goodman & Weare, 2010) to obtain posterior probability distributions for the free model parameters. The only exception is the phase shift, for which the maximum likelihood solution in each fit was found. Additionally, the intrinsic scatter of the model was set as a free parameter $\log\sigma_{\mathrm{i}}$. This is a measure of the systematic errors from the choice of the model (see e.g. Salmi et al., 2018). We assumed the prior probability distributions to be uniform in all of the parameters. The limits of the priors were set to (11 km, 13 km) in $R_{\mathrm{eq}}$, (0∘, 90∘) in $\theta$, (1∘, 40∘) in $\rho$, (0.928 MK, 4.062 MK) in $T$, and (0.868, 5.212) in $\log\sigma_{\mathrm{i}}$. The synthetic pulse-profile data were binned into 16 phase bins and NICER energy channels located between 0.3 and 10 keV. In addition, we required each modelled energy-phase bin to have more than 20 observed counts.

3 Results

3.1 Spectral properties

We began our calculations by checking that our code gives similar results to McPHAC when we used Thomson scattering instead of Compton. This is shown in Figs. 1 and 2, with the former showing the emergent spectrum and the latter the angular dependencies of the emitted radiation. Figure 1 also shows the results computed with the full Compton model. The comparison between the angular dependencies given by that model and those of McPHAC are shown in Fig. 3. The parameters of the model, temperature and surface gravity, were chosen to be $T=$ 3.1623 MK and $\log g=14.3856$, which are reasonable for RMPs.

From the aforementioned figures, we see that the calculations with McPHAC agree with the Thomson version of our code within a few per cent, albeit displaying a small systematic discrepancy that increases with energy, and is probably connected to the increasing error at high zenith angles (i.e. small $\mu$) seen in Fig. 2. The largest difference is about 3%. In addition, the effective temperature produced by McPHAC code is slightly higher, meaning the energy conservation is not extremely accurate. However, this should only have a minor effect to the fitted effective temperature. A much larger difference is seen between our Compton and the Thomson models, which also becomes more significant at higher energies (above 3 keV for the chosen temperature) and small $\mu$. This difference in spectrum is similar to that presented in Suleimanov & Werner (2007).

Taking into account the energy response matrix of the NICER instrument, we also show the modelled phase-averaged count spectra in Fig. 4 with $T=$ 2.0 MK and in Fig. 5 with $T=$ 3.1 MK (other parameters being the same as in Sect. 2). The data produced with the Compton model in Fig. 5 also represent our synthetic data in the following sections. From the figures, we see that the discrepancy between the models at the highest energies can be partly hidden because of only a few detected counts, and therefore large statistical errors. We also assumed the calibration error of the instrument to be 1%. In any case, a clearly observable difference above 3 keV remains when $T=$ 3.1 MK.

3.2 Parameter constraints with the correct model

We applied the method described in Sect. 2 with full treatment of Compton scattering in the atmosphere, to fit the synthetic data created using the same model. Due to the ignorance of heating by magnetospheric return currents, our data do not resemble what is expected in real sources. We still confirm the robustness of our method by getting no strong biases in the constraints for radius and other parameters. The fitted pulse profiles are shown in Fig. 6 (integrated to three energy bins), and the posterior probability distributions are shown in Fig.  7. The credible limits of all parameters are also listed in Table 1. The best-fit solution presented in Fig. 6 has $\chi^{2}/\mathrm{d.o.f.}=6711/(6736-6)\approx 1.00$ (for six free parameters including the phase shift), when ignoring the calibration error, which is used only for fitting purposes but not actually present in the synthetic data.

As expected, the model accurately describes the synthetic data, as seen in the posterior probability distribution for intrinsic scatter $\sigma_{\mathrm{i}}$. The mean $\log\sigma_{\mathrm{i}}<1$ of the posterior translates to an error of less than ten counts in each phase-energy bin. This effectively means zero intrinsic scatter, as it is significantly smaller than the Poisson noise of the data, which is $56$ counts on average in each fitted phase-energy bin. For the radius, we find the 68% (95%) limits and the most probable value as $R_{\mathrm{eq}}=12.01^{+0.01\leavevmode\nobreak\ (0.02)}_{-0.01\leavevmode\nobreak\ (0.02)}$ km. We note that this, and the other limits presented here and in the following section, are considerably tighter than what is expected, if comparing, for example, to the approximation in Eq. (5) by Psaltis et al. (2014), which has been used to predict 5% accuracy for the NICER targets. With our model parameters and the amount of detected counts ($4\times 10^{7}$), we should have about 1% accuracy. Our even tighter limits could be due to the anisotropic effects (ignored in the aforementioned equation), which can strongly increase the second harmonic of the pulse profile signal (Poutanen & Beloborodov, 2006), and thus decrease the uncertainty in its measurement (as we regard the atmospheric effects to be known). In any case, this is not critical given that we are only interested in the differences between the two spectral models. Similarly, tight constraints are found for other parameters so that the correct point remains inside their 68% limits, except for the temperature where the correct point is slightly offset towards smaller values, but is still inside the 95% limits.

3.3 Parameter constraints with the incorrect model

We also applied the method described in Sect. 2, using NS atmosphere model McPHAC to fit the synthetic data that were created using the full Compton scattering model. The fitted pulse profiles, integrated to three energy bins, are shown in Fig. 8 for illustration. We see that the fits are worse at the highest energies due to a large difference in the spectral shapes. The posterior probability distributions are shown in Fig. 9, and the credible limits are listed in Table 1.

We find that the constraints for radius are still not biased, but very close to those obtained in the previous section, since $R_{\mathrm{eq}}=12.01^{+0.01\leavevmode\nobreak\ (0.02)}_{-0.01\leavevmode\nobreak\ (0.02)}$ km. However, the credible limits for the temperature, and especially for the size of the spot, are clearly different: the temperature is higher while the spot size is smaller than the correct values. Neither of them agrees with the $95$% limits. The best-fit solution presented in Fig. 8 has $\chi^{2}/\mathrm{d.o.f.}=11142/6730\approx 1.66$. In addition, according to the notably higher intrinsic scatter $\sigma_{\mathrm{i}}$, the model does not describe the synthetic data as well as the correct model that includes Compton scattering (although $\sigma_{\mathrm{i}}$ is still effectively very small compared to the Poisson noise).

We also calculated the results assuming two other NS masses. With otherwise a similar setup as discussed above, the synthetic data were created with masses $1.4{M}_{\sun}$ and $2.0{M}_{\sun}$ and then fitted with both Thomson and Compton models. We find no major difference compared to the already presented results using an NS mass of $1.76{M}_{\sun}$. Although, in the case of the Thomson model, biases in temperature $T$ and spot size $\rho$ are found to depend on the NS mass. The bias is always significant and tends towards the same direction, but it is higher with higher masses. Further, the slope between radius $R_{\mathrm{eq}}$ and magnetic co-latitude $\theta$ in the two-dimensional posterior probability histogram is different for every mass because the star is more oblate with higher masses. In all cases the input radius is still obtained at least within the 95% limits.

4 Conclusions

We studied the possible outcomes of using Thomson scattering approximation in the atmosphere calculation instead of a full Compton scattering model when trying to constrain NS parameters from RMP pulse profile observations of NICER. Our spectral comparisons showed that the difference in the observed spectrum may not be detected, due to the low count rate at highest energies, if the temperature of the emitting hotspot is $T=2$ MK. However, in case of $T\approx 3$ MK, a significant discrepancy can be observed.

We simulated and fitted synthetic data, based on the Compton atmosphere with $T=3.1$ MK and the NICER target PSR J0437$-$4715, which is expected to give some of the most constraining limits to NS radius. Fitting with the same Compton model, we obtained very tight limits for the NS parameters without strong biases, demonstrating the robustness of our method. Likewise, fitting with the Thomson model resulted in very similar constraints on the radius. However, the obtained size and the temperature of the hotspot were significantly different. The exact credible limits should not be taken too seriously, as we have exaggerated the predicted count rate in order to emphasize the differences between the two spectral models.

According to our results, Compton scattering seems to be unimportant in obtaining accurate radius constraints for RMPs, at least in the case of a similar model and comparable data to that used here. However, for the interpretation of the data from a mission that is more sensitive at high energies and observes more counts at the energies around and above 3 keV, the effects of Compton scattering would need to be taken into account in a precise manner. These effects will be even more important for atmospheres heated in the surface layers by bombarding particles, which were not considered in this paper but are expected to be present in real sources.

Acknowledgments

This research was supported by the University of Turku Graduate School in Physical and Chemical Sciences (TS), by the Ministry of Science and Higher Education of the Russian Federation grant 14.W03.31.0021 (JP, VFS), the Deutsche Forschungsgemeinschaft (DFG) grant WE 1312/51-1, the travel grant of the German Academic Exchange Service (DAAD, project 57405000), and the Academy of Finland travel grant 317552 (JP). The computer resources of the Finnish IT Center for Science (CSC) and the Finnish Grid and Cloud Infrastructure project are acknowledged.