Issue 
A&A
Volume 678, October 2023



Article Number  A99  
Number of page(s)  9  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202346833  
Published online  12 October 2023 
Polarized radiation from an accretion shock in accreting millisecond pulsars using exact Compton scattering formalism
^{1}
Tuorla Observatory, Department of Physics and Astronomy, University of Turku, 20014 Turku, Finland
email: anna.a.bobrikova@utu.fi
^{2}
Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
Received:
8
May
2023
Accepted:
11
August
2023
Pulse profiles of accreting millisecond pulsars can be used to determine neutron star (NS) parameters, such as their masses and radii, and therefore provide constraints on the equation of state of cold dense matter. Information obtained by the Imaging Xray Polarimetry Explorer (IXPE) can be used to decipher pulsar inclination and magnetic obliquity, providing ever tighter constraints on other parameters. In this paper, we develop a new emission model for accretionpowered millisecond pulsars based on thermal Comptonization in an accretion shock above the NS surface. The shock structure was approximated by an isothermal planeparallel slab and the Stokes parameters of the emergent radiation were computed as a function of the zenith angle and energy for different values of the electron temperature, the Thomson optical depth of the slab, and the temperature of the seed blackbody photons. We show that our Compton scattering model leads to a significantly lower polarization degree of the emitted radiation compared to the previously used Thomson scattering model. We computed a large grid of shock models, which can be combined with pulse profile modeling techniques both with and without polarization included. In this work, we used the relativistic rotating vector model for the oblate NS in order to produce the observed Stokes parameters as a function of the pulsar phase. Furthermore, we simulated the data to be produced by IXPE and obtained constraints on model parameters using nested sampling. The developed methods can also be used in the analysis of the data from future satellites, such as the enhanced Xray Timing and Polarimetry mission.
Key words: methods: numerical / polarization / stars: neutron / techniques: polarimetric / Xrays: binaries
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Rapidly rotating accreting neutron stars (NSs) can be used to study the properties of the extremely dense matter inside their cores. Some of these are accretionpowered millisecond pulsars (AMPs), which show Xray pulsations at the stellar spin frequency. Their emission profiles can be modeled, accounting for light bending and other relativistic effects, and from these profiles, the information about the mass and radius of the NS can be extracted (see e.g., Pechenick et al. 1983; Miller & Lamb 1998; Poutanen & Gierliński 2003; Poutanen & Beloborodov 2006; Morsink et al. 2007; Lo et al. 2013; Watts et al. 2016; Salmi et al. 2018; Bogdanov et al. 2019). The mass and radius constraints can be translated to the equation of state (EOS) constraints for matter inside the NS core (e.g., see Lindblom 1992; Lattimer 2012; Hebeler et al. 2013; Baym et al. 2018).
This technique was recently applied to the case of rotationpowered millisecond pulsars using data from the Neutron Star Interior Composition Explorer (NICER; Miller et al. 2019, 2021; Riley et al. 2019, 2021; Salmi et al. 2022). In these sources, the hot regions on the NS surface are heated by a magnetospheric return current and the emission escaping the NS surface can be computed using selfconsistent atmosphere models (see e.g., Zavlin et al. 1996; Ho & Lai 2001; Heinke et al. 2006; Ho & Heinke 2009; Haakonsen et al. 2012; Salmi et al. 2020). For AMPs, these models are not adequate because the emitted radiation is Comptonized in an accretion “shock” formed above the hot region. In reality, a standard gasdynamical shock might not exist at all, but the kinetic energy of incoming particles is released in the surface layer through Coulomb collisions and plasma instabilities (Zel’dovich & Shakura 1969; Suleimanov et al. 2018; Salmi et al. 2020). Comptonization of soft photons generated within the layer as well as by the underlying NS surface is then responsible for the powerlaw spectra extending up to 100 keV as observed from AMPs (Poutanen & Gierliński 2003; Gierliński & Poutanen 2005; Falanga et al. 2005a,b, 2007, 2011, 2012). Electron scattering also causes the radiation to be significantly polarized. Variations of polarization with pulsar phase can be used to determine geometrical parameters such as inclination angle and magnetic obliquity (Viironen & Poutanen 2004), which in their turn allow us to improve constraints on the NS mass and radius.
Previously, polarized emission models for AMPs used the Thomson scattering approximation in an optically thin NS atmosphere (see Sunyaev & Titarchuk 1985; Viironen & Poutanen 2004; Salmi et al. 2021). However, the polarization degree due to Compton scattering on hot electrons strongly depends on the electron temperature (Nagirner & Poutanen 1994; Poutanen 1994); for example, for 100 keV electrons, polarization is less than 50% of that for Thomson scattering. On the other hand, the models for nonpolarized radiation – in the context of pulse profile modeling – typically employ approximate formulas for the anisotropy of the radiation and empirical models for the Comptonized spectra (see e.g., Poutanen & Gierliński 2003; Leahy et al. 2008; Steiner et al. 2009; Salmi et al. 2018). Selfconsistent accretionheated atmosphere models have also been developed (see e.g., Zampieri et al. 1995; Deufel et al. 2001; Suleimanov et al. 2018), although these are relatively computationally expensive and have a large number of free parameters.
In this work, we applied the formalism for Compton scattering in a hot slab (Nagirner & Poutanen 1993; Poutanen & Svensson 1996) to compute the Stokes parameters of the emergent radiation as a function of energy and emission angle using only three model parameters: the electron temperature and optical depth of the hot slab on top of the NS surface, and the temperature of the seed blackbody photons coming from the NS. By employing pulse profile modeling for polarized radiation from rapidly rotating oblate NSs (Poutanen 2020; Loktev et al. 2020), we also simulated data for the Imaging Xray Polarimeter Explorer (IXPE; Weisskopf et al. 2022) and updated the NS geometry parameter constraints previously predicted in Salmi et al. (2021) using the Thomson scattering model. These methods can also be applied when analyzing the observations from IXPE^{1} or from future Xray polarimetric missions such as the enhanced Xray Timing and Polarimetry mission (eXTP; Zhang et al. 2019; Watts et al. 2019). The atmosphere lookup tables produced here can also be combined with any AMP pulse profile modeling codes^{2}.
The remainder of this paper is structured as follows. In Sect. 2, we present the theory that describes the formation of polarized radiation in a Comptonizing slab above the NS surface. We then apply this theory to obtain the intensity and polarization of the radiation escaping from the slab. We then describe the method to obtain pulse profiles using our new radiation model. We compare our resulting spectra and pulse profiles to those obtained using previous models. We also discuss pulse profiles obtained for different NS parameters. In Sect. 3, we generate synthetic data and apply our fitting routine to determine the NS parameters from the data. We discuss the applications of our model and predictions for the upcoming observations in Sect. 4. In Sect. 5, we conclude and summarize our findings.
2. Emission model
2.1. Radiative transfer equation
We consider a simple NS atmosphere model in which the atmosphere is a plane parallel slab consisting of electrons and lying above an optically thick source that radiates as a blackbody. The hot slab has the Thomson optical depth τ_{T} = σ_{T}n_{e}H, where H is the vertical height, σ_{T} is the Thomson cross section, and n_{e} is the electron concentration. The electron gas in the atmosphere is considered to be isotropic and isothermal with the electron temperature T_{e}. The momentum distribution of electrons is given by relativistic Maxwellian distribution characterized by the dimensionless temperature Θ_{e} = kT_{e}/m_{e}c^{2}:
where γ is the electron Lorentz factor and K_{2} is the modified Bessel function of the second kind. This distribution is normalized to unity:
where is the electron dimensionless momentum. The radiation field can be described by the Stokes vector I(τ, x, μ), where μ is the cosine of the zenith angle, an angle between the normal to the slab and the direction of photon propagation, and x = E/m_{e}c^{2} is the photon energy measured in the units of the electron rest mass.
The Stokes vector usually contains four components I, Q, U, and V, but because of the azimuthal symmetry and the absence of sources of circular polarization, the latter two are equal to zero. We therefore use just two Stokes parameters I(τ, x, μ) = (I, Q)^{T}, where superscript T denotes the transposed vector.
The photon distribution at the bottom of the slab is considered to be Planckian of the temperature T_{bb}. The incident radiation (for zeroth scattering order, n = 0) at the bottom of the slab is given by the Stokes vector
where Θ_{bb} = kT_{bb}/m_{e}c^{2}. In order to describe the propagation of polarized radiation through the hot electron slab, we solve the radiative transfer equation (RTE; Nagirner & Poutanen 1994; Poutanen & Svensson 1996):
where dτ = σ_{T}n_{e}dz is the Thomson optical depth, S is the source function (also a Stokes vector), and σ(x) is the dimensionless Compton scattering crosssection (in units of the Thomson crosssection σ_{T}). We solve the RTE using the iterative scattering method of Poutanen & Svensson (1996), where the intensity is represented as a series expansion in scattering orders. The Stokes vector of unscattered radiation at all optical depths is
From the known Stokes vector at the nth scattering order, we find the source function for the following scattering order as
where is the 2 × 2 azimuth averaged redistribution matrix describing Compton scattering by isotropic electrons (see Appendix A.1 in Poutanen & Svensson 1996). If the source function is known, the Stokes vector is found via the formal solution of the RTE:
Iterations proceed until the required accuracy of the total Stokes vector I = ∑_{n}I_{n} is achieved so that the maximal contribution of the next scattering is less than 1% of the total spectrum in all energies and angles.
Finally, we can compute the polarization degree (PD) of the emergent radiation as
The positive values of P imply that the polarization vector lies in the meridional plane defined by the slab normal and the photon momentum, while for the negative values, the polarization vector is perpendicular to that plane, as in the case of optically thick electron scattering atmosphere (Chandrasekhar 1960; Sobolev 1963).
2.2. Spectropolarimetric properties of accretion shocks
We can now apply the derived formalism to predict the beaming patterns, energy spectra, and polarization of the accretion shock at the NS surface. We take as a fiducial set of parameters: τ_{T} = 1, T_{e} = 50 keV, T_{bb} = 1 keV. Figures 1 and 2 show the spectrum and the PD as functions of the cosine of the zenith angle and energy, respectively, for this set of parameters. The parameters of the source are the same as those used in the Thomson model of Salmi et al. (2021; see their Figs. 1 and 2), and so these figures can be used to illustrate the impact of using an exact Compton scattering description instead of the Thomson scattering model in the slab of hot electrons.
Fig. 1. Intensity of the emergent radiation (left) and the PD (right) as functions of the cosine of the zenith angle. The atmosphere parameters are τ_{T} = 1.0, electron temperature T_{e} = 50 keV, and seed photon temperature T_{bb} = 1 keV. Black, blue, green, orange, and red solid lines show the model for photon energies of 2, 5, 8, 12, and 18 keV, respectively. The angular dependence of the intensity is normalized so that . 
Fig. 2. Intensity of emergent radiation (left) and the PD (right) as functions of photon energy. The same atmosphere parameters are used as in Fig. 1. Black, blue, green, orange, red, and magenta solid lines show the model for μ = 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0 respectively. 
When comparing our results with those of Salmi et al. (2021), we note that the PD differs significantly while the spectrum remains mainly unchanged. This effect is most obviously seen in the angular dependence of the PD. The absolute values of the PD are reduced by a factor of two for the photon energies of 2, 5, and 8 keV, and increased at higher energies. For Compton scattering, the PD changes the sign already at the photon energy of 8 keV and is positive at higher values, while for Thomson scattering the sign changes at photon energies above 18 keV. Last but not least, the maximum in the absolute value of the PD for the photon energies higher than 8 keV has shifted from μ = 0 to higher values of the cosine of the zenith angle. Keeping in mind the energy range of the IXPE satellite (2–8 keV), we can say that compared to the previously developed model for AMPs that used Thomson scattering, the Compton scattering model predicts a lower PD. This result is in agreement with our expectations, as for hot electrons the relativistic aberration becomes important (Poutanen 1994). In the case of cold electrons, the PD of the Thomson singly scattered radiation is highly dependent on the scattering angle, reaching 100% at 90°. For relativistic electrons (but still in the Thomson scattering regime), the same happens in the electron rest frame; however, because of the relativistic motion of the electron, the photon in the external frame is preferentially scattered in the direction of the electron motion, resulting in zero final PD independent of the scattering angle (Bonometto et al. 1970; Nagirner & Poutanen 1993). In our case of 50 keV electrons, the electrons are not yet relativistic and the maximum PD is reduced to 60% from the Thomson case. We can predict here that with the higher electron temperature, the PD will decrease further. This will be examined in the following section.
2.3. Grid of models
The model we consider for the accretion shock has three parameters: the optical depth of the Comptonizing plasma τ_{T}, the dimensionless electron temperature Θ_{e}, and the dimensionless characteristic photon temperature Θ_{bb} (associated with T_{e} and T_{bb}, respectively). To determine these parameters from the upcoming IXPE data, we need a dense grid of models. For the optical depth of the shock τ_{T}, we chose the lower and higher value of 0.5 and 3.5, respectively, and varied the parameter in steps of 0.1. For the electron temperature, we varied Θ_{e} from 0.04 to 0.2 (i.e., T_{e} from ∼20 to 100 keV) with a step of 0.004 (i.e., ∼2 keV). For the seed photon temperature, the dimensionless Θ_{bb} was varied in the interval (1–3)×10^{−3} (i.e., T_{bb} = 0.5–1.5 keV) in steps of 2 × 10^{−4} (i.e., ∼0.1 keV). Altogether, we computed 13 981 models.
The impact of the NS atmosphere parameters on the PD can be seen in Fig. 3. For the electron temperature, we see the expected decrease in the PD with T_{e} as electrons become more relativistic. For the two other parameters, the main feature is that with increasing the parameter value, the change of sign of the PD happens at higher energies. Within the mentioned intervals of parameter values, the impact is most significant with increasing the optical depth of the atmosphere and least significant with changes to the seed photon energy. For the IXPE observational range of 2–8 keV, this means that, in general, PD increases with the increase in these parameters.
Fig. 3. Angular dependency of the PD for different atmosphere parameters. For each plot, we vary only one parameter (marked on the plot) while other parameters are from the fiducial set (τ_{T} = 1, T_{e} = 50 keV, T_{bb} = 1 keV). Cyan, blue, green, orange, and red solid lines show the PD on photon energies of 2, 5, 8, 12, and 18 keV, respectively. 
In order to use the calculated Stokes parameters more effectively, we wrote a 3D interpolation routine over this grid of parameters, so that for any arbitrary combination of these parameters within the intervals mentioned, we can obtain the intensity and PD of the emission as functions of energy and zenith angle. We tested that the results of this interpolation match the direct calculation results for any combination of the parameters within the grid to within 8% (within 1% in the IXPE energy range).
2.4. Pulse profiles and polarization
We calculated the pulse profiles from the known beaming pattern and PD exactly as in Salmi et al. (2021). After obtaining the beaming pattern and PD of the emission on top of the NS atmosphere, we transport them from the hot emitting regions of the star (called hot spots) to the observer following the formalism described in Morsink et al. (2007), AlGendy & Morsink (2014), Poutanen (2020), and Loktev et al. (2020). We consider an oblate, rapidly rotating star in the Schwarzschild metric (hence, the effects of rotation on spacetime are neglected). Gravitational redshift and the lightbending effect, as well as the effect of the fast rotation of the star, are accounted for in the polarization angle (PA) transfer. The only difference in approach is the reference radius used in timedelay calculation. While Loktev et al. (2020) use the exact radius of the star at the center of the spot, we follow the path described in Poutanen & Beloborodov (2006) and applied in Salmi et al. (2021) and use the equatorial radius as a reference radius.
We again return to the results from Salmi et al. (2021) and compare the pulse profiles coming from the model with Compton scattering in the atmosphere with previously obtained ones using Thomson scattering. Figure 4 presents this comparison for the case of two antipodal spots, both combined and individually from each spot, for the NS with atmospheric parameters from the fiducial set (inclination i = 60° and magnetic colatitude θ = 20° as in Salmi et al. 2021). We see that the flux does not change significantly, which is in agreement with the beaming patterns not changing either. The PD is now decreasing from 2 to 5 keV and then slowly increasing towards 8 keV. From inspection of the corresponding lines in Fig. 1, and accounting for the gravitational redshift, it is clear that the PD is following the same pattern for the corresponding energies of the emitted photons. The ‘jump’ by 90° in the PA for the 8 keV photons comes from the PD changing sign.
Fig. 4. Pulse profiles of the observed flux, PD, and PA for two antipodal spots shown for three different energies (2, 5, 8 keV). The solid black curves correspond to the total flux (PD, PA), the dashed blue lines correspond to the contribution from the primary spot, and the dotted red lines are for the secondary spot. The dashdotted gray lines correspond to the Thomson model from Salmi et al. (2021; combined from two spots). The NS parameters in both models are those shown in Table 1 of Salmi et al. (2021). 
The main conclusion here is that we expect a significant drop in the observed PD compared to that predicted by Salmi et al. (2021). The impact of this reduction in PD on the possibility of observing the AMPs is further discussed in Sects. 3 and 4.
3. Synthetic data and their analysis
3.1. Generating the synthetic data
The developed model can be further used to analyze the data. Polarimetric measurements of the AMPs are yet to be carried out. In the absence of observations, we used the IXPEOBSSIM framework (Baldini et al. 2022; version 16) with the PCUBE algorithm to generate several samples of synthetic data imitating the data to be obtained from the IXPE observations of an AMP. We substituted the source with the model described above and obtained the event list for a specified observing time of 600 ks for a source with a flux of 100 mCrab. In this part of our study, we do not investigate the energy dependency of PD, and collect all the data in one energy bin. These data sets are in exactly the same form as publicly available IXPE data products, so we will be able to apply the developed routines to the real data as soon as they become available.
As mentioned above, applying the Compton scattering formalism reduces the expected PD as compared to the previous studies where Thomson scattering in the NS atmosphere was considered. We synthesized the data for the same parameter set as used in Fig. 4 (see Table 1, Model 1) and compared the PD to the socalled minimal detectable polarization (MDP) at 99% significance (see Weisskopf et al. 2010 for a definition of the MDP and Kislat et al. 2015 for further details) values for the corresponding phase bins. Figure 5a illustrates that the obtained values of the PD are below the MDP. According to Weisskopf et al. (2010), this means that we cannot have a proper measurement of the PD in this case. From this point, there are three main ways to improve the quality of the potential measurement: (i) we can increase the number of counts to reduce the MDP values, (ii) study different energy bands to see if the PD values are higher there, and/or (iii) choose a more promising source, which in terms of simulated data means selecting different source parameters.
Fig. 5. Simulated PD and normalized Stokes q and u profiles for the NS parameter set 1 (panels a and b are for the 2–8 keV and 2–5 keV energy range, respectively) and set 4 (panel c) from Table 1. The exposure time is 600 ks in all three cases. The blue lines show the theoretical model, the purple dots are the synthesized observed data, and the red dots are the MDP values for the corresponding phase bins. 
Parameter sets for all the computed models and the most probable values obtained in the fitting.
Increasing the number of counts by prolonging the exposure time is impossible with AMPs, as they are transient sources observable only during the outbursts, which typically last for a few weeks (see e.g., Patruno & Watts 2021). Also, observations lasting longer than 600 ks are unrealistic for scheduling reasons. The two other options are explored below.
Studying different energy bands is relevant as, according to our predictions from Sect. 2.1, the PD changes sign within the 2–8 energy band. This means that averaging over all energies decreases the total PD value. Figure 5b illustrates the results of analyzing the 2–5 keV energy band instead of the 2–8 keV represented in Fig. 5a. The blue line corresponding to the model prediction is now significantly higher and closer to the MDP values. However, the improvement is slight, and we note that the PD values are still below the MDP values in this case.
In our search for a more optimistic (from the point of view of potential observations) yet still valid set of both geometrical and atmosphere parameters, we aimed to reproduce the known spectra of the AMPs with the parameters that, coming from our previous investigations, are expected to give a higher PD of the emission. From Fig. 3, we see that the optical depth affects the PD strongest from all the atmosphere parameters, and so we aimed to increase it but preserve the overall spectrum behavior presented in Poutanen & Gierliński (2003). We also attempted to increase the electron temperature to match the predictions in Poutanen & Gierliński (2003), but the scenario proposed there gives PD values below the MDP values, and so we chose not to continue with it. We developed two possible adjustments to the fiducial scenario (Model 1): we can choose the geometry such that the angle between the inclination and colatitude of the spot is almost 90° (as in Model 2) or use ‘more promising’ atmosphere parameters (as in Model 3). We can also combine these two (Model 4). The parameter values for these four models are given in the left part of Table 1 in the columns labeled correspondingly 1–4. Figure 5c shows that for Model 4, the PD now exceeds the MDP values. We also checked this for Models 2 and 3, and for Model 2 the PD is slightly higher than the MDP values, while for Model 3 the PD is similar to MDP. In all the models, we only calculated the emission coming from one spot on the surface of the NS.
In addition, we attempted to improve the quality of the results for Model 1 illustrated in Fig. 5a,b by reducing the number of phase bins, as the PD varies less in these cases compared to other models. We find that for Model 1, case a, only reducing the number of bins to one gives us the PD value that barely exceeds the MDP, meaning that we can detect something from such a source, but we would lose the phase variability completely. For Model 1, case b, we can obtain a detection with three phase bins, but only for the bin that includes the peak (in Fig. 5b the corresponding phase shift is 0.5), and again, the PD barely exceeds the MDP. We conclude that reducing the number of bins can support the observations with barely measurable PD, but its impact is not sufficient for us to rely upon it in cases where the source is otherwise unobservable with IXPE.
Nevertheless, the most recent development of the IXPE data analysis shows that for multiplepoint observations (i.e., observations split in different phase bins), information about polarization can be extracted even from the data points where PD values do not exceed MDP (GonzálezCaniulef et al. 2023; Suleimanov et al. 2023). As IXPE measures the Stokes parameters q and u rather than PD, we can fit these parameters directly and obtain some constraints on the geometry of the source. This approach is further explored in the following section.
3.2. Fitting the data
After synthesizing the data, we applied the MULTINEST (Feroz et al. 2009) multimodal nested sampling algorithm (specifically, the PYMULTINEST package, see Buchner et al. 2014) to obtain new constraints on both the geometrical and atmosphere parameters of the source. As in Salmi et al. (2021), we fit the geometrical parameters, such as inclination, colatitude, phase shift, and position angle, but unlike there, we also perform a simultaneous fit over the parameters of the atmosphere, such as characteristic photon temperature, electron temperature, and optical depth of the thermalization layer. Table 1 presents the most probable values and 68% credible intervals for these parameters. Figure 6 shows the results of the fitting presented for the fiducial set of parameters (Model 1) and for a ‘more promising’ scenario (Model 4). In both cases, the full energy band of IXPE (2–8 keV) is considered, and the exposure time is 600 ks.
Fig. 6. Posterior probability distributions for NS parameters for Model 1a (left panels) and Model 4 (right panels) from Table 1 using one data realization. The twodimensional contours correspond to 68% and 95% credible regions with the corresponding regions colored black and gray. The histograms show the normalized onedimensional distribution for a given parameter derived from the posterior samples. The mean value and the 68% and 95% credible levels for the parameters are shown by the vertical dashed black, red, and orange lines, respectively. Blue lines indicate the true values of the parameters that were used for the data simulations. 
Figure 6 clearly illustrates the improvement of the fitting results if comparing the data synthesized from the Model 1type of source (left) and the data coming from the Model 4type of source (right). As we use the uniform prior for all our fittings, in the former case, the constraints originate from our direct fitting of the q and u Stokes parameters. This illustrates the possibility to extract some information even when the MDP values are not reached. In the latter case, we can clearly see the significant constraints coming from the PD values exceeding the MDP by approximately a factor of two.
We can also conclude that the constraints on the atmospheric parameters are notably less restrictive than the ones on the geometrical parameters of the NS. The main reason for obtaining such broad constraints on atmospheric parameters lies in the energy range of the IXPE satellite, as we are unable to detect the whole spectrum and fit it correctly. A solution to this issue could be to use the observational results from several missions simultaneously; for instance, fitting the NuSTAR data in its higher energy range to determine the electron temperature (and consequently the optical depth through the spectral slope). Then, using the values obtained, we can fit the IXPE data with the fixed (or constrained) atmospheric parameters. However, we see that even with large uncertainties on these parameters, we are able to derive constraints on the geometry of the source from the IXPE data only.
Table 1 also shows that the constraints from the Model 3 data are less restrictive than the ones from Model 2, and so we can conclude that the atmospheric parameters are less important in our search for the ‘promising’ source than the geometric parameters of the NS. We also attempted to improve the constraints by reducing the number of simultaneously fitted parameters (e.g., we fixed the atmospheric parameters and only fitted the geometrical parameters). Decreasing the number of fitted parameters reduces the calculation time, but the credible intervals for the fitted parameters remain the same.
In the fittings presented here, we fixed the mass and radius of the NS. This is done for the simplicity of the calculations, as they are quite timeconsuming. Moreover, we only fit the normalized Stokes q and u parameters, and they are not particularly sensitive to the mass or radius. We could also fit the Stokes I parameter (i.e., flux) as well and add mass and radius as free parameters. Adding two more parameters would increase the computational time significantly, but we do not expect any significant impact of this on the credible intervals of other parameters.
4. Discussion
Our fiducial set of NS parameters is based on the SAX J1808.4−3658, and Fig. 4 shows that, for the parameters selected, our predictions for the PD of the SAX J1808.4−3658like source are low. However, with the uncertainties we still have on both the geometry and the atmospheric parameters of the source, there is still a possibility to observe a higher PD. Furthermore, with about 20 AMPs discovered, we can put forward several conclusions about the perspective of observing measurable PD from this class of sources.
It is known that the product of the electron temperature and optical depth for the AMPs is almost invariant (Poutanen 2006). Nevertheless, we have learned that all the atmospheric parameters have their own contribution to the PD values; see Fig. 3. Optical depth is the most important parameter here; we are searching for a source with an atmosphere that is sufficiently optically thick for many Compton scattering events to occur, yet not so thick that the pulsation pattern is excessively influenced by these scatterings. In this sense, objects such as IGR J17591−2342 and IGR J17511−3057 are promising. The former has an estimated optical depth of 1.59–2.3 and a pulsed fraction of the emission of 10%–17% (Kuiper et al. 2020; Manca et al. 2023). The latter has a lower optical depth of 1.34, but higher electron and seed photons temperatures (T_{e} = 51 keV, T_{bb} = 1.36 keV) and a pulsed fraction of 14% (Papitto et al. 2010). Constraints on the geometrical parameters of the AMPs are among the goals of the polarimetric observations. So far our main source of information on the possible geometry is the amplitude of the pulsation, which depends approximately on the product of the sines of the inclination and the magnetic colatitude (Poutanen & Beloborodov 2006). Moreover, while for the inclination we usually have some constraints, little is known about the geometrical configuration of the magnetic field of the AMP. From the emission patterns shown in Fig. 1, the highest PD is received when the line of sight and the normal to the emitting area are orthogonal. One of the most interesting objects to study is Swift J1749.4−2807, for which the inclination is wellconstrained, 74.4° < i < 77.3°, and the colatitude is estimated to be θ ≈ 50° (Altamirano et al. 2010). In this case, we can expect high polarization and we can test whether the inclination obtained from the polarimetry actually matches the orbital one.
We performed all our fittings with 600 ks of exposure time. It is typical for AMPs to go into outburst for a couple of weeks, and so in order to collect enough photons we need either a rather soft spectrum of emission or a very bright source. We can exclude hard sources such as IGR J18245−2452 with a photon index of Γ ≈ 1.4 (Papitto et al. 2013).
Figure 4 illustrates an additional feature of the PD behavior we expect to observe: the observed PD is different between cases where the secondary spot is invisible (e.g., the very beginning of the outburst) and where it is visible (at the later stages of the outburst). In order to get additional information about the geometry of the source, it is crucial to observe the outburst from the earliest to the latest possible stages.
We purposely focused on phaseresolved analysis and did not investigate the possibilities that energyresolved analysis could open up. IXPE has an energy resolution of 0.57 keV (at 2 keV; Ramsey et al. 2021), and so with sufficient statistics it is possible to study the behavior of PD with energy and obtain more information about the geometry and atmospheric properties of the source.
The model presented in this paper can be improved further. For instance, the isothermal (static) slab approximation of the accretion shock can be replaced with a more physical model where the dynamics of the accreting gas is computed selfconsistently with the electron temperature. Also, the model for seed soft photons can be modified to account more accurately for the reprocessing of the hard Xrays at the NS surface, accounting for the polarization properties of this radiation. We intend to work on these developments in the future.
5. Summary
We present a new model to compute the flux and polarization of the emission coming from the accretionpowered millisecond pulsars. We used the formalism for Compton scattering in a planeparallel isothermal hot slab to describe the propagation of photons in the accretion shock above the NS surface. Figures 1 and 2 show the result of such calculations for the fiducial set of parameters taken from Salmi et al. (2021). Comparing our results with the previous studies using Thomson scattering formalism shows that, at least for the same NS parameters, our model predicts lower PD values.
We calculated the Stokes parameters of the radiation as a function of two variables, the zenith angle of the emitted photon and its energy, for almost 14 000 combinations of three parameters: electron temperature, the temperature of the seed blackbody photons, and the Thomson optical depth of the slab. The tables of the spectral energy distribution of the Stokes parameters for all the parameters can now be used for further analysis of the accretionpowered millisecond pulsars, both in studying the pulse profiles and the polarimetric properties of the emerged radiation. Figure 3 illustrates one of the ways to use this new data set. We can understand the impact of the NS atmosphere parameters on the angular dependency of the PD from this particular study.
We then used the relativistic rotating vector model for the oblate NS to produce observed Stokes parameters as a function of the pulsar phase. Figure 4 shows the observed fluxes coming from the two antipodal spots of an AMP and the combined total flux as calculated using our model and a total flux calculated using a previously developed Thomson model. A significant drop in the PD between these two models is again present here.
We used the developed model of the AMP emission to generate the synthetic data sets imitating the one that the IXPE satellite could produce after observing an AMP. We note that the NS with the parameters from our fiducial set would emit the light with PD values below the MDP values of IXPE (which still does not exclude the possibility to extract information about the geometry of the source); however, for a different set of NS parameters, the PD is high enough to be detectable by IXPE, as shown in Fig. 5. We then applied the multimodal nested sampling technique to fit the simulated data with our model and obtain constraints on both the geometrical and atmospheric parameters of the NS. The results of the fitting are presented in Fig. 6. The main conclusion of this study is that we can obtain reasonably good constraints on the geometrical parameters from the IXPE data for some of the AMPs; however, the parameters of the atmosphere of the NS need to be constrained from the supporting observations by other Xray missions such as NuSTAR.
Our findings contribute to the advancement of understanding AMPs and have significant implications for future research in this field. We are waiting for the AMPs to be observed by IXPE so that we may apply the developed formalism to analyze the resulting data. We also hope to see the eXTP mission observing the AMPs, as it would solve several of the issues discussed in this article: the larger effective area will lower the MDP, and so the fainter sources will become observable; the broader energy range will allow us to develop better constraints on the atmospheric parameters; and finally, the wide field instrument of eXTP will be able to catch the transients, such as AMPs, and observe the outbursts from their very beginning.
The tables are available at https://github.com/AnnaBobrikova/ComptonSlabTables
Acknowledgments
This research was supported by the Academy of Finland grant 333112 and by the Finnish Cultural Foundation grant 00220175 (AB). TS acknowledges support from ERC Consolidator Grant (CoG) No. 865768 AEONS (PI: A. Watts). The computer resources of the Finnish IT Center for Science (CSC) and the FGCI project (Finland) are acknowledged. We thank Anna Watts, Bas Dorsman, Alessandro Di Marco, and John Rankin for discussions and testing our models.
References
 AlGendy, M., & Morsink, S. M. 2014, ApJ, 791, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Altamirano, D., Cavecchi, Y., Patruno, A., et al. 2010, ApJ, 727, L18 [Google Scholar]
 Baldini, L., Bucciantini, N., Lalla, N. D., et al. 2022, SoftwareX, 19, 101194 [NASA ADS] [CrossRef] [Google Scholar]
 Baym, G., Hatsuda, T., Kojo, T., et al. 2018, Rep. Prog. Phys., 81, 056902 [NASA ADS] [CrossRef] [Google Scholar]
 Bogdanov, S., Lamb, F. K., Mahmoodifar, S., et al. 2019, ApJ, 887, L26 [NASA ADS] [CrossRef] [Google Scholar]
 Bonometto, S., Cazzola, P., & Saggion, A. 1970, A&A, 7, 292 [NASA ADS] [Google Scholar]
 Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chandrasekhar, S. 1960, Radiative Transfer (New York: Dover) [Google Scholar]
 Deufel, B., Dullemond, C. P., & Spruit, H. C. 2001, A&A, 377, 955 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Falanga, M., BonnetBidaud, J. M., Poutanen, J., et al. 2005a, A&A, 436, 647 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Falanga, M., Kuiper, L., Poutanen, J., et al. 2005b, A&A, 444, 15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Falanga, M., Poutanen, J., Bonning, E. W., et al. 2007, A&A, 464, 1069 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Falanga, M., Kuiper, L., Poutanen, J., et al. 2011, A&A, 529, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Falanga, M., Kuiper, L., Poutanen, J., et al. 2012, A&A, 545, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601 [NASA ADS] [CrossRef] [Google Scholar]
 Gierliński, M., & Poutanen, J. 2005, MNRAS, 359, 1261 [Google Scholar]
 GonzálezCaniulef, D., Caiazzo, I., & Heyl, J. 2023, MNRAS, 519, 5902 [CrossRef] [Google Scholar]
 Haakonsen, C. B., Turner, M. L., Tacik, N. A., & Rutledge, R. E. 2012, ApJ, 749, 52 [CrossRef] [Google Scholar]
 Hebeler, K., Lattimer, J. M., Pethick, C. J., & Schwenk, A. 2013, ApJ, 773, 11 [CrossRef] [Google Scholar]
 Heinke, C. O., Rybicki, G. B., Narayan, R., & Grindlay, J. E. 2006, ApJ, 644, 1090 [NASA ADS] [CrossRef] [Google Scholar]
 Ho, W. C. G., & Heinke, C. O. 2009, Nature, 462, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Ho, W. C. G., & Lai, D. 2001, MNRAS, 327, 1081 [Google Scholar]
 Kislat, F., Clark, B., Beilicke, M., & Krawczynski, H. 2015, Astropart. Phys., 68, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Kuiper, L., Tsygankov, S. S., Falanga, M., et al. 2020, A&A, 641, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lattimer, J. M. 2012, Annu. Rev. Nucl. Part. Sci., 62, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Leahy, D. A., Morsink, S. M., & Cadeau, C. 2008, ApJ, 672, 1119 [Google Scholar]
 Lindblom, L. 1992, ApJ, 398, 569 [CrossRef] [Google Scholar]
 Lo, K. H., Miller, M. C., Bhattacharyya, S., & Lamb, F. K. 2013, ApJ, 776, 19 [Google Scholar]
 Loktev, V., Salmi, T., Nättilä, J., & Poutanen, J. 2020, A&A, 643, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Manca, A., Gambino, A. F., Sanna, A., et al. 2023, MNRAS, 519, 2309 [Google Scholar]
 Miller, M. C., & Lamb, F. K. 1998, ApJ, 499, L37 [CrossRef] [Google Scholar]
 Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24 [Google Scholar]
 Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, ApJ, 918, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Morsink, S. M., Leahy, D. A., Cadeau, C., & Braga, J. 2007, ApJ, 663, 1244 [NASA ADS] [CrossRef] [Google Scholar]
 Nagirner, D. I., & Poutanen, J. 1993, A&A, 275, 325 [NASA ADS] [Google Scholar]
 Nagirner, D. I., & Poutanen, J. 1994, Astrophys. Space. Phys. Rev., 9, 1 [Google Scholar]
 Papitto, A., Riggio, A., Di Salvo, T., et al. 2010, MNRAS, 407, 2575 [Google Scholar]
 Papitto, A., Ferrigno, C., Bozzo, E., et al. 2013, Nature, 501, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Patruno, A., & Watts, A. L. 2021, Astrophys. Space Sci. Lib., 461, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Pechenick, K. R., Ftaclas, C., & Cohen, J. M. 1983, ApJ, 274, 846 [Google Scholar]
 Poutanen, J. 1994, ApJS, 92, 607 [NASA ADS] [CrossRef] [Google Scholar]
 Poutanen, J. 2006, AdSpR, 38, 2697 [NASA ADS] [Google Scholar]
 Poutanen, J. 2020, A&A, 641, A166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poutanen, J., & Beloborodov, A. M. 2006, MNRAS, 373, 836 [NASA ADS] [CrossRef] [Google Scholar]
 Poutanen, J., & Gierliński, M. 2003, MNRAS, 343, 1301 [NASA ADS] [CrossRef] [Google Scholar]
 Poutanen, J., & Svensson, R. 1996, ApJ, 470, 249 [CrossRef] [Google Scholar]
 Ramsey, B. D., Attina, P., Baldini, L., et al. 2021, Proc. SPIE, 11821, 118210M [NASA ADS] [Google Scholar]
 Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21 [Google Scholar]
 Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJ, 918, L27 [CrossRef] [Google Scholar]
 Salmi, T., Nättilä, J., & Poutanen, J. 2018, A&A, 618, A161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Salmi, T., Suleimanov, V. F., Nättilä, J., & Poutanen, J. 2020, A&A, 641, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Salmi, T., Loktev, V., Korsman, K., et al. 2021, A&A, 646, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Salmi, T., Vinciguerra, S., Choudhury, D., et al. 2022, ApJ, 941, 150 [CrossRef] [Google Scholar]
 Sobolev, V. V. 1963, A Treatise on Radiative Transfer (Princeton: Van Nostrand) [Google Scholar]
 Steiner, J. F., Narayan, R., McClintock, J. E., & Ebisawa, K. 2009, PASP, 121, 1279 [Google Scholar]
 Suleimanov, V. F., Poutanen, J., & Werner, K. 2018, A&A, 619, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suleimanov, V. F., Forsblom, S. V., Tsygankov, S. S., et al. 2023, A&A, in press, https://doi.org/10.1051/00046361/202346994 [Google Scholar]
 Sunyaev, R. A., & Titarchuk, L. G. 1985, A&A, 143, 374 [NASA ADS] [Google Scholar]
 Viironen, K., & Poutanen, J. 2004, A&A, 426, 985 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Watts, A. L., Andersson, N., Chakrabarty, D., et al. 2016, Rev. Mod. Phys., 88, 021001 [NASA ADS] [CrossRef] [Google Scholar]
 Watts, A. L., Yu, W., Poutanen, J., et al. 2019, Sci. China Phys. Mech. Astron., 62, 29503 [Google Scholar]
 Weisskopf, M. C., Elsner, R. F., & O’Dell, S. L. 2010, Proc. SPIE, 7732, 77320E [NASA ADS] [CrossRef] [Google Scholar]
 Weisskopf, M. C., Soffitta, P., Baldini, L., et al. 2022, J. Astron. Telesc. Instrum. Syst., 8, 026002 [NASA ADS] [CrossRef] [Google Scholar]
 Zampieri, L., Turolla, R., Zane, S., & Treves, A. 1995, ApJ, 439, 849 [Google Scholar]
 Zavlin, V. E., Pavlov, G. G., & Shibanov, Y. A. 1996, A&A, 315, 141 [NASA ADS] [Google Scholar]
 Zel’dovich, Y. B., & Shakura, N. I. 1969, Soviet Ast., 13, 175 [Google Scholar]
 Zhang, S., Santangelo, A., Feroci, M., et al. 2019, Sci. China Phys. Mech. Astron., 62, 29502 [Google Scholar]
All Tables
Parameter sets for all the computed models and the most probable values obtained in the fitting.
All Figures
Fig. 1. Intensity of the emergent radiation (left) and the PD (right) as functions of the cosine of the zenith angle. The atmosphere parameters are τ_{T} = 1.0, electron temperature T_{e} = 50 keV, and seed photon temperature T_{bb} = 1 keV. Black, blue, green, orange, and red solid lines show the model for photon energies of 2, 5, 8, 12, and 18 keV, respectively. The angular dependence of the intensity is normalized so that . 

In the text 
Fig. 2. Intensity of emergent radiation (left) and the PD (right) as functions of photon energy. The same atmosphere parameters are used as in Fig. 1. Black, blue, green, orange, red, and magenta solid lines show the model for μ = 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0 respectively. 

In the text 
Fig. 3. Angular dependency of the PD for different atmosphere parameters. For each plot, we vary only one parameter (marked on the plot) while other parameters are from the fiducial set (τ_{T} = 1, T_{e} = 50 keV, T_{bb} = 1 keV). Cyan, blue, green, orange, and red solid lines show the PD on photon energies of 2, 5, 8, 12, and 18 keV, respectively. 

In the text 
Fig. 4. Pulse profiles of the observed flux, PD, and PA for two antipodal spots shown for three different energies (2, 5, 8 keV). The solid black curves correspond to the total flux (PD, PA), the dashed blue lines correspond to the contribution from the primary spot, and the dotted red lines are for the secondary spot. The dashdotted gray lines correspond to the Thomson model from Salmi et al. (2021; combined from two spots). The NS parameters in both models are those shown in Table 1 of Salmi et al. (2021). 

In the text 
Fig. 5. Simulated PD and normalized Stokes q and u profiles for the NS parameter set 1 (panels a and b are for the 2–8 keV and 2–5 keV energy range, respectively) and set 4 (panel c) from Table 1. The exposure time is 600 ks in all three cases. The blue lines show the theoretical model, the purple dots are the synthesized observed data, and the red dots are the MDP values for the corresponding phase bins. 

In the text 
Fig. 6. Posterior probability distributions for NS parameters for Model 1a (left panels) and Model 4 (right panels) from Table 1 using one data realization. The twodimensional contours correspond to 68% and 95% credible regions with the corresponding regions colored black and gray. The histograms show the normalized onedimensional distribution for a given parameter derived from the posterior samples. The mean value and the 68% and 95% credible levels for the parameters are shown by the vertical dashed black, red, and orange lines, respectively. Blue lines indicate the true values of the parameters that were used for the data simulations. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.