© ESO 2020

1. Introduction

Millisecond pulsars (MSPs) are rapidly rotating neutron stars (NSs) with relatively low magnetic field strengths. Many of the known MSPs are found in binary systems where the in-falling gas from the companion star spins up the central NS to extreme angular velocities (Radhakrishnan & Srinivasan 1982; Alpar et al. 1982). During the accretion phase, matter falls on the NS surface near the magnetic poles creating hotspots that emit polarized X-ray radiation. X-ray polarization from accretion-powered millisecond pulsars (AMPs) will be soon measured by upcoming space observatories like the Imaging X-ray Polarimeter Explorer (IXPE; Weisskopf et al. 2016) and the enhanced X-ray Timing and Polarimetry mission (eXTP; Zhang et al. 2019; Watts et al. 2019). The pulse profiles contain information about parameters of the NS, such as its mass and radius (e.g., Poutanen & Gierliński 2003; Watts et al. 2016; Salmi et al. 2018; Bogdanov et al. 2019). The relation between these parameters can give constraints on the equation of state of extremely dense matter constituting the inner parts of NSs (see, e.g., Lindblom 1992; Lattimer 2012).

The core theory on pulse profile modeling was presented in Pechenick et al. (1983), and was later extended for rapidly rotating spherical NSs by Poutanen & Gierliński (2003) and Poutanen & Beloborodov (2006). They used the Schwarzschild+Doppler approximation (S+D), where the effects of light bending and gravitational redshift are computed in Schwarzschild metric and special relativistic effects (aberration, Doppler effect) are included separately. In addition to the modulation of the X-ray flux, Viironen & Poutanen (2004) modeled the polarization degree (PD) and polarization position angle (PA) from an infinitely small hotspot (or two antipodal spots) on the surface of rapidly rotating spherical NS. Poutanen (2020a) presented a derivation of the formulae from Viironen & Poutanen (2004) using a common vector formalism and also studied the effects of rapid rotation on the polarization profile.

The effects of the deviation of the metric and stellar shape from spherical were considered by Braje et al. (2000), Cadeau et al. (2007), Bauböck et al. (2012), and recently by Vincent et al. (2018), Nättilä & Pihajoki (2018), and Pihajoki et al. (2018). For realistic NS equations of state and spins, the dominant effect comes from the oblateness of the NS surface. This prompted Cadeau et al. (2007) and Morsink et al. (2007) to suggest an oblate Schwarzschild (OS) approximation to describe the flux modulation from a rapidly rotating NS. This approximation accounts for the deviation of the NS shape from a sphere, but treats the surrounding space-time using the spherically symmetric Schwarzschild metric. The corrected OS formalism, also valid for finite spot size, was recently presented by Salmi et al. (2018), Bogdanov et al. (2019), and Suleimanov et al. (2020), but without considering polarization. Here we combine the techniques described in those papers to compute the pulse profiles with the methods suggested by Viironen & Poutanen (2004) and Poutanen (2020a) to compute polarization. With this combination we develop the complete theory needed to calculate the Stokes vector light curves from rapidly rotating oblate NSs in the OS approximation. We derive analytical expressions for apparent PA for a point source at the oblate surface of a rapidly rotating NS (which are also applied to calculate the PA from a finite-sized spot). Then we use these formulae to model pulse profiles and polarization and study the importance of the oblateness effects in such modeling.

The remainder of this paper is structured as follows. In Sect. 2, we present the theory that describes the transfer of polarized radiation from one or two hotspots at the surface of an oblate NS to the observer at infinity. In Sect. 3 we first confirm the accuracy of our model by comparing the polarization pulse profiles to those calculated with the general relativistic polarized radiative transfer code ARCMANCER (Pihajoki et al. 2018) using oblate star shape and higher order spacetime metric corrections. Then we use the theoretical model to fit synthetic PA data to estimate how well the NS parameters can be constrained and to see how the results differ if an incorrect shape of the star is assumed. We conclude in Sect. 4.

2. Theory

2.1. NS shape and mass

We consider a NS of mass M rotating at a high frequency ν_* > 300 Hz, at which the stellar shape flattens significantly (Cadeau et al. 2007; AlGendy & Morsink 2014). The shape of the star is given by the Schwarzschild coordinate R(θ), which is a function of the co-latitude. We use one of the recent rapidly rotating NS shape models presented by AlGendy & Morsink (2014). According to the model, the star shape is given by the approximation

$$\begin{array}{c}\hfill R(\theta )=R_{\mathrm{e}}(1+o_2(x,\overline{\mathrm{\Omega }}){cos}^2\theta ),\end{array}$$(1)

where the flattening parameter, which defines the difference between the equatorial radius R_e and the polar radius R_p, is given in the first-order approximation as

$$\begin{array}{c}\hfill o_2(x,\overline{\mathrm{\Omega }})=\frac{R_{\mathrm{p}}-R_{\mathrm{e}}}{R_{\mathrm{e}}}\approx {\overline{\mathrm{\Omega }}}^2(1.030x-0.788).\end{array}$$(2)

Here

$$\begin{array}{c}\hfill x=\frac{r_{\mathrm{s}}}{2R_{\mathrm{e}}}\end{array}$$(3)

is the dimensionless compactness parameter, r_s = 2GM/c² is the Schwarzschild radius of the NS, and

$$\begin{array}{c}\hfill \overline{\mathrm{\Omega }}=2\pi {\nu }_{\ast }{\left(\frac{R_{\mathrm{e}}^3}{\mathit{GM}}\right)}^{\frac{1}{2}}\end{array}$$(4)

is the dimensionless circular rotation frequency. This approximation gives the shape of the star with better than 1% accuracy when ${\overline{\mathrm{\Omega }}}^2<0.1$, which is fulfilled for most of the realistic equations of state for NSs with ν_* < 700 Hz and M > M_⊙.

When comparing results for oblate and spherical stars, we also account for the difference in mass and the equatorial radius of the same NS due to its rapid rotation. For this purpose, following Suleimanov et al. (2020), we relate the gravitational mass M and equatorial radius R_e of a rotating NS to the mass M₀ and radius R₀ of a spherical non-rotating NS of the same baryonic mass,

$$\begin{array}{c}\hfill M=M_0(a_0+\frac{a_1}{1.1-\overline{\nu }}+a_2{\overline{\nu }}^2)\end{array}$$(5)

and

$$\begin{array}{c}\hfill R_{\mathrm{e}}=R_0(0.9766+\frac{0.025}{1.07-\overline{\nu }}+0.07M_{1.4}^{3/2}{\overline{\nu }}^2),\end{array}$$(6)

where the coefficients are a₀ = 1 − a₁/1.1, $a_1=0.001M_{1.4}^{3/2}$, and a₂ = 10a₁. Here M_1.4 = M₀/1.4 M_⊙ and $\overline{\nu }={\nu }_{\ast }/{\nu }_{\mathrm{cr}}$ is the relative rotation frequency in units of the critical rotation frequency

$$\begin{array}{c}\hfill {\nu }_{\mathrm{cr}}=1278M_{1.4}^{1/2}{\left(\frac{10\mathrm{km}}{R_0}\right)}^{3/2}\mathrm{Hz}.\end{array}$$(7)

The dimensionless speed (in units of speed of light c) of the NS surface at a given co-latitude as measured by a static observer is

$$\begin{array}{c}\hfill \beta (\theta )=\frac{v}{c}=\frac{2\pi R(\theta )sin\theta }{\sqrt{1-u(\theta )}}\frac{{\nu }_{\ast }}{c},\end{array}$$(8)

where we have corrected the observed pulsar frequency for gravitational redshift, which can be expressed as

$$\begin{array}{c}\hfill \frac{1}{\sqrt{1-u(\theta )}}={(1-\frac{r_{\mathrm{s}}}{R(\theta )})}^{-1/2}=1+z(\theta ).\end{array}$$(9)

The corresponding Lorentz factor at a given θ is

$$\begin{array}{c}\hfill \gamma (\theta )=\frac{1}{\sqrt{1-{\beta }^2(\theta )}}.\end{array}$$(10)

The angle η between the local normal to the NS surface and the radial direction depends only on latitude. The values of trigonometric functions for this angle can be represented as

$$\begin{array}{c}\hfill cos\eta =\frac{1}{\sqrt{1+f^2(\theta )}}\mathrm{and}sin\eta =\frac{f(\theta )}{\sqrt{1+f^2(\theta )}},\end{array}$$(11)

where the function f(θ) can be expressed in terms of function R(θ) and its derivative

$$\begin{array}{c}\hfill f(\theta )=\frac{1+z(\theta )}{R}\frac{\mathrm{d}R}{\mathrm{d}\theta }.\end{array}$$(12)

2.2. Vectors and angles

Let us next recall the vectors and the angles between them in the OS approximation. The geometry of the star is the same as in Salmi et al. (2018), and it is shown in Fig. 1. The axis of pulsar rotation coincides with the symmetry axis, and we choose it to be the applicate axis z in our fixed laboratory frame. The unit vector $\widehat{\mathit{k}}$ pointing to the direction of the observer from the star center lies in the xz-plane and makes angle i with the rotation axis:

$$\begin{array}{c}\hfill \widehat{\mathit{k}}=(sini,0,cosi).\end{array}$$(13)

Fig. 1. Geometry of the oblate star and the lab frame coordinate system xyz. A spot at the NS surface at co-latitude θ has a rotation phase ϕ from the observer vector $\widehat{\mathit{k}}$ around rotational axis z. Three unit vectors are defined for the spot: $\widehat{\mathit{n}}$ is normal to the surface, $\widehat{\mathit{r}}$ is the radius vector, and ${\widehat{\mathit{k}}}_0$ is the direction of emitted photon propagation.

The y-axis completes the right-handed Cartesian system.

We consider first a small hotspot being placed on the star surface and a constant co-latitude θ. For larger spots we integrate observed Stokes parameters over the spot area. The phase of the pulsar ϕ is the angle measured counterclockwise around the z-axis. We define the phase to be zero when the spot (or its center for large circular spots) passes through the xz-plane. When we consider more than one spot we identify the spot closest to the observer as the main emission region, which also defines the phase.

For every moment of time we define three local unit vectors on the surface related to the spot. They are the radius vector $\widehat{\mathit{r}}$ of the spot, the vector normal to the spot $\widehat{\mathit{n}}$, and the vector ${\widehat{\mathit{k}}}_0$ denoting the initial direction of light propagation emitted so that a photon eventually arrives at the observer. Correspondingly, we have three angles between these vectors, so that η is between $\widehat{\mathit{n}}$ and $\widehat{\mathit{r}}$, α is between $\widehat{\mathit{r}}$ and ${\widehat{\mathit{k}}}_0$, and σ is between $\widehat{\mathit{n}}$ and ${\widehat{\mathit{k}}}_0$. The applied coordinate system, unit vectors, and related angles are visualized in Fig. 1.

Each vector can be represented by a set of coordinates in the lab system. Thus, the radius-vector $\widehat{\mathit{r}}$ is

$$\begin{array}{c}\hfill \widehat{\mathit{r}}=(sin\theta cos\varphi ,sin\theta sin\varphi ,cos\theta ),\end{array}$$(14)

and the normal is

$$\begin{array}{c}\hfill \widehat{\mathit{n}}=(sin\lambda cos\varphi ,sin\lambda sin\varphi ,cos\lambda ),\end{array}$$(15)

where λ = θ − η is the angle between the normal vector and the spin axis. The vector ${\widehat{\mathit{k}}}_0$ lying between $\widehat{\mathit{r}}$ and $\widehat{\mathit{k}}$ is expressed as

$$\begin{array}{c}\hfill {\widehat{\mathit{k}}}_0=\frac{1}{sin\psi }(sin\alpha \widehat{\mathit{k}}+sin(\psi -\alpha )\widehat{\mathit{r}}),\end{array}$$(16)

where the angle ψ between $\widehat{\mathit{r}}$ and $\widehat{\mathit{k}}$ is determined by the scalar product of these vectors:

$$\begin{array}{c}\hfill cos\psi =\widehat{\mathit{r}}·\widehat{\mathit{k}}=cosicos\theta +sinisin\theta cos\varphi .\end{array}$$(17)

The relation between angle α and ψ depending on the radius of the emission point R(θ) is given by the light bending integral (Pechenick et al. 1983; Beloborodov 2002), which can be computed either exactly as described in Salmi et al. (2018), for example, or via the approximate analytical formulae in Poutanen (2020b). The angle σ can be expressed as (Morsink et al. 2007; Suleimanov et al. 2020)

$$\begin{array}{cc}& cos\sigma =\widehat{\mathit{n}}·{\widehat{\mathit{k}}}_0=cos\alpha cos\eta \hfill \\ \hfill & +\frac{sin\alpha }{sin\psi }sin\eta (cosisin\theta -sinicos\theta cos\varphi ).\hfill \end{array}$$(18)

Because the star is rotating, the spot is moving in the lab frame with dimensionless velocity $\mathit{\beta }=\widehat{\mathit{\beta }}\beta$, which we decompose into the product of the unit velocity vector

$$\begin{array}{c}\hfill \widehat{\mathit{\beta }}=(-sin\varphi ,cos\varphi ,0)\end{array}$$(19)

and the dimensionless velocity β given by Eq. (8).

2.3. Observed flux and polarization degree

The flux observed at energy E from a small spot of area dS′ (measured in the comoving frame) can be obtained from the following equation (Poutanen & Beloborodov 2006; Morsink et al. 2007; Salmi et al. 2018; Suleimanov et al. 2020):

$$\begin{array}{c}\hfill \mathrm{d}{F}_E={(1-u)}^{3/2}{\delta }^4{I}^{\prime }(E^{\prime },\mu )cos\sigma \mathcal{D}\frac{\mathrm{d}{S}^{\prime }}{D^2}.\end{array}$$(20)

Here u is the compactness at the emission point;

$$\begin{array}{c}\hfill \delta =\frac{1}{\gamma (1-\beta cos\xi )}\end{array}$$(21)

is the Doppler-factor, where ξ is the angle between the spot velocity vector and the photon direction,

$$\begin{array}{c}\hfill cos\xi =\widehat{\mathit{\beta }}·{\widehat{\mathit{k}}}_0=-\frac{sin\alpha }{sin\psi }sinisin\varphi ;\end{array}$$(22)

μ denotes the cosine of the inclination σ′ of the line of sight to the normal in the co-rotating reference frame expressed as (Poutanen & Gierliński 2003; Poutanen & Beloborodov 2006)

$$\begin{array}{c}\hfill \mu =cos{\sigma }^{\prime }=\delta cos\sigma ;\end{array}$$(23)

and I′(E′,μ) is the spectral intensity measured in the comoving frame at the zenith angle σ′. Furthermore,

$$\begin{array}{c}\hfill \mathcal{D}=\frac{1}{1-u}\frac{\mathrm{d}cos\alpha }{\mathrm{d}cos\psi }\end{array}$$(24)

is the lensing factor and D is the distance to the star (we use D = 1 since it does not affect the PA profile). We note that D ≫ R to ensure the validity of the presented equations because we implicitly assume that the observer’s image plane is at infinity.

The focus of this paper is the effects of NS oblateness on the PA, which means that the choice of the spectral model is not very important. We approximate the spectrum of the emitted photons in the local corotating frame at energy E′ with the Planck function B_E′(T_bb). However, the dependence of intensity on cosine of the zenith angle μ = cosσ′ is taken from that expected for X-ray bursting millisecond pulsars (Chandrasekhar & Breen 1947; Chandrasekhar 1960; Sobolev 1949, 1963; Suleimanov et al. 2012):

$$\begin{array}{c}\hfill I(E^{\prime },\mu )=B_{E^{\prime }}(T_{\mathrm{bb}})(0.421+0.868\mu ).\end{array}$$(25)

Because the PD does not change along photon trajectory, the observed PD from a small spot is fully determined by the PD at the zenith angle σ′ as measured in the comoving frame. For the simulations below we assume that the angular dependency of the PD corresponds to that of the electron-scattering dominating atmosphere, which we approximate as (Viironen & Poutanen 2004)

$$\begin{array}{c}\hfill P(\mu )=-\frac{1-\mu }{1+3.582\mu }11.71\%.\end{array}$$(26)

The negative value for the PD means that the dominant direction of the oscillation of the electric vector is perpendicular to the plane formed by the local normal and the photon momentum.

For a rapidly rotating NS, time delays of photon propagation around the star have to be accounted for. The observed phase (i.e., the phase when photon arrives to the observer) differs from the rotation (emission) phase ϕ,

$$\begin{array}{c}\hfill \phi =\varphi +2\pi {\nu }_{\ast }\mathrm{\Delta }t,\end{array}$$(27)

where Δt is the time delay calculated as explained in Salmi et al. (2018), among others. However, unlike there, we use the radius at the center of the main spot R(θ) as a reference radius, and a reference emission angle is chosen to correspond to the emitted photon direction from that point at ϕ = 0. By this definition, the time delay Δt is zero for a single small spot at ϕ = 0, and thus same for both oblate and spherical stars at that phase.

2.4. Polarization angle

Let us now present derivation of the PA, which differs from that presented for spherical stars in the previous papers (Viironen & Poutanen 2004; Poutanen 2020a). There are several steps that need to be made in order to transform the Stokes vector from the co-rotating frame of the star to the sky plane. Let us first define the main polarization basis, in the sky plane, formed by $\widehat{\mathit{z}}$ and $\widehat{\mathit{k}}$,

$$\begin{array}{c}\hfill {\widehat{\mathit{e}}}_1^{\mathrm{m}}=\frac{\widehat{\mathit{z}}-cosi\widehat{\mathit{k}}}{sini},{\widehat{\mathit{e}}}_2^{\mathrm{m}}=\frac{\widehat{\mathit{k}}\times \widehat{\mathit{z}}}{sini}.\end{array}$$(28)

The vector ${\widehat{\mathit{e}}}_1^{\mathrm{m}}$ is colinear to the projection of the star’s rotation axis at the sky plane and the second vector ${\widehat{\mathit{e}}}_2^{\mathrm{m}}$ is perpendicular to the first one. This basis is fixed in the lab frame because the pulsar rotation axis (along the $\widehat{\mathit{\omega }}$) is very stable and the line of sight (along the $\widehat{\mathit{k}}$) does not change during the observation. The polarization (pseudo-)vector is observed in this basis. The PA defines the direction of the polarization vector in a basis, and it is customary to measure it in the counterclockwise direction from the first basis vector ${\widehat{\mathit{e}}}_1^{\mathrm{m}}$.

Next let us consider the basis formed by $\widehat{\mathit{r}}$ and $\widehat{\mathit{k}}$:

$$\begin{array}{c}\hfill {\widehat{\mathit{e}}}_1^{\psi }=\frac{\widehat{\mathit{r}}-cos\psi \widehat{\mathit{k}}}{sin\psi },{\widehat{\mathit{e}}}_2^{\psi }=\frac{\widehat{\mathit{k}}\times \widehat{\mathit{r}}}{sin\psi }.\end{array}$$(29)

This basis is defined in the same plane as the previous one, as it is also perpendicular to $\widehat{\mathit{k}}$. However, this basis is rotated relative to the main one. The angle between these two bases we denote by χ₀, and the corresponding trigonometric functions are

$$\begin{array}{c}\hfill cos{\chi }_0={\widehat{\mathit{e}}}_1^{\mathrm{m}}·{\widehat{\mathit{e}}}_1^{\psi }=\frac{sinicos\theta -cosisin\theta cos\varphi }{sin\psi }\end{array}$$(30)

and

$$\begin{array}{c}\hfill sin{\chi }_0=-{\widehat{\mathit{e}}}_1^{\mathrm{m}}·{\widehat{\mathit{e}}}_2^{\psi }={\widehat{\mathit{e}}}_2^{\mathrm{m}}·{\widehat{\mathit{e}}}_1^{\psi }=-\frac{sin\theta sin\varphi }{sin\psi }.\end{array}$$(31)

Thus, we get the standard expression for the PA of the rotating vector model (RVM; Radhakrishnan & Cooke 1969):

$$\begin{array}{c}\hfill tan{\chi }_0=\frac{-sin\theta sin\varphi }{sinicos\theta -cosisin\theta cos\varphi }.\end{array}$$(32)

For the spherical star and slow rotation this is sufficient; for the oblate NS shape and rapid rotation additional corrections are needed.

To take into consideration the surface curvature we define another basis formed by the radius vector $\widehat{\mathit{r}}$ and the direction of the light propagation ${\widehat{\mathit{k}}}_0$ near the NS surface:

$$\begin{array}{c}\hfill {\widehat{\mathit{e}}}_1^{\alpha }=\frac{\widehat{\mathit{r}}-cos\alpha {\widehat{\mathit{k}}}_0}{sin\alpha },{\widehat{\mathit{e}}}_2^{\alpha }=\frac{{\widehat{\mathit{k}}}_0\times \widehat{\mathit{r}}}{sin\alpha }={\widehat{\mathit{e}}}_2^{\psi }.\end{array}$$(33)

This basis describes the plane perpendicular to the initial photon direction close to the NS surface. Because in the Schwarzschild metric the light trajectories are planar, vectors ${\widehat{\mathit{e}}}_2^{\alpha }$ and ${\widehat{\mathit{e}}}_2^{\psi }$ coincide. The PA measured in this basis near the surface is the same as that measured by the instrument in the basis denoted by upper index ψ in Eq. (29) because the polarization vector is parallel-transported along photon trajectory. In addition, we introduce the basis associated with the local normal vector $\widehat{\mathit{n}}$:

$$\begin{array}{c}\hfill {\widehat{\mathit{e}}}_1^{\sigma }=\frac{\widehat{\mathit{n}}-cos\sigma {\widehat{\mathit{k}}}_0}{sin\sigma },{\widehat{\mathit{e}}}_2^{\sigma }=\frac{{\widehat{\mathit{k}}}_0\times \widehat{\mathit{n}}}{sin\sigma }.\end{array}$$(34)

This basis lies in the same plane as the previous one, perpendicular to ${\widehat{\mathit{k}}}_0$. The angle between them, χ₁, can be computed from

$$\begin{array}{c}\hfill cos{\chi }_1={\widehat{\mathit{e}}}_1^{\sigma }·{\widehat{\mathit{e}}}_1^{\alpha }=\frac{cos\eta -cos\alpha cos\sigma }{sin\alpha sin\sigma }\end{array}$$(35)

and

$$\begin{array}{cc}& sin{\chi }_1={\widehat{\mathit{e}}}_2^{\alpha }·{\widehat{\mathit{e}}}_1^{\sigma }=\frac{\widehat{\mathit{n}}·({\widehat{\mathit{k}}}_0\times \widehat{\mathit{r}})}{sin\alpha sin\sigma }=\frac{{\widehat{\mathit{k}}}_0·(\widehat{\mathit{r}}\times \widehat{\mathit{n}})}{sin\alpha sin\sigma }\hfill \\ \hfill & =\frac{sin\alpha }{sin\psi }\frac{\widehat{\mathit{k}}·(\widehat{\mathit{r}}\times \widehat{\mathit{n}})}{sin\alpha sin\sigma }=\frac{sin\eta sinisin\varphi }{sin\psi sin\sigma }.\hfill \end{array}$$(36)

Thus, we get

$$\begin{array}{c}\hfill tan{\chi }_1=-\frac{sin\eta cos\xi }{cos\eta -cos\alpha cos\sigma }.\end{array}$$(37)

If the star is spherical (i.e., η = 0), angle χ₁ is exactly zero.

The last step is to take into account the positional vector rotation due to the relativistic motion of the NS surface. The relativistic correction of the polarization vector we denote by χ′. It can be shown that Lorentz transformation of the polarization basis vector ${\widehat{\mathit{e}}}_1^{\sigma }$ is (see, e.g., Nagirner & Poutanen 1993)

$$\begin{array}{c}\hfill {\widehat{\mathit{e}}}_1^{\prime }=\frac{\widehat{\mathit{n}}-\mu (\delta {\widehat{\mathit{k}}}_0-\gamma \mathit{\beta })}{\sqrt{1-{\mu }^2}}.\end{array}$$(38)

The expressions for the trigonometric functions of relativistic rotation angle of the polarization plane are

$$\begin{array}{cc}& cos{\chi }^{\prime }={\widehat{\mathit{e}}}_1^{\prime }·{\widehat{\mathit{e}}}_1^{\sigma }=\frac{1-{cos}^2\sigma -\delta \gamma {cos}^2\sigma \beta cos\xi }{sin\sigma \sqrt{1-{\mu }^2}}\hfill \\ \hfill & =\frac{{sin}^2\sigma -\beta cos\xi }{sin\sigma \sqrt{1-{\mu }^2}(1-\beta cos\xi )}\hfill \end{array}$$(39)

and

$$\begin{array}{cc}\hfill sin{\chi }^{\prime }& ={\widehat{\mathit{e}}}_1^{\prime }·{\widehat{\mathit{e}}}_2^{\sigma }=\frac{\mu \gamma \beta }{sin\sigma \sqrt{1-{\mu }^2}}\widehat{\mathit{\beta }}·({\widehat{\mathit{k}}}_0\times \widehat{\mathit{n}})\hfill \\ \hfill & =\frac{\beta cos\sigma }{sin\sigma \sqrt{1-{\mu }^2}(1-\beta cos\xi )}{\widehat{\mathit{k}}}_0·(\widehat{\mathit{n}}\times \widehat{\mathit{\beta }}).\hfill \end{array}$$(40)

The vector product $\widehat{\mathit{n}}\times \widehat{\mathit{\beta }}$ is essentially the meridional vector pointed towards the north pole:

$$\begin{array}{c}\hfill \widehat{\mathit{m}}=\widehat{\mathit{n}}\times \widehat{\mathit{\beta }}=(-cos\lambda cos\varphi ,-cos\lambda sin\varphi ,sin\lambda ).\end{array}$$(41)

Then using Eq. (16) we can express the scalar triple product as

$$\begin{array}{c}\hfill {\widehat{\mathit{k}}}_0·(\widehat{\mathit{n}}\times \widehat{\mathit{\beta }})={\widehat{\mathit{k}}}_0·\widehat{\mathit{m}}=\frac{1}{sin\psi }(sin\alpha cos\zeta -sin(\psi -\alpha )sin\eta ),\end{array}$$(42)

where ζ is the angle the vector $\widehat{\mathit{k}}$ makes with the meridian

$$\begin{array}{c}\hfill cos\zeta =\widehat{\mathit{m}}·\widehat{\mathit{k}}=cosisin\lambda -sinicos\lambda cos\varphi .\end{array}$$(43)

Thus, we get

$$\begin{array}{c}\hfill sin{\chi }^{\prime }=\beta cos\sigma \frac{cos\zeta sin\alpha -sin(\psi -\alpha )sin\eta }{sin\psi sin\sigma \sqrt{1-{\mu }^2}(1-\beta cos\xi )}.\end{array}$$(44)

Combing Eqs. (22), (39), and (44) we get

$$\begin{array}{c}\hfill tan{\chi }^{\prime }=\beta cos\sigma \frac{cos\zeta -sin\eta \frac{sin(\psi -\alpha )}{sin\alpha }}{\frac{sin\psi }{sin\alpha }{sin}^2\sigma +\beta sinisin\varphi }.\end{array}$$(45)

For the case of a spherical star we can set η = 0 and recover

$$\begin{array}{c}\hfill tan{\chi }_{\mathrm{sph}}^{\prime }=\beta cos\alpha \frac{cosisin\theta -sinicos\theta cos\varphi }{sin\psi sin\alpha +\beta sinisin\varphi },\end{array}$$(46)

which coincides with Eq. (29) from Viironen & Poutanen (2004) and derived by Poutanen (2020a). An alternative derivation of Eq. (45) is presented in Appendix A.

The total PA for each spot is then obtained by summing the angles between the intermediate bases since they are constant for parallel transport along the photon trajectory:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{obl}}={\chi }_0+{\chi }_1+{\chi }^{\prime }.\end{array}$$(47)

Similarly, for the spherical case this reads

$$\begin{array}{c}\hfill {\chi }_{\mathrm{sph}}={\chi }_0+{\chi }_{\mathrm{sph}}^{\prime }.\end{array}$$(48)

2.5. Observed Stokes parameters

Let us finally describe how to obtain the final Stokes parameters and the total PA. Because we do not consider circular polarization, the Stokes vector that describes observed radiation from an infinitely small spot of area dS′ contains three components (dF_I, dF_Q, dF_U). Without losing the generality, we assume that the PA equals zero in the polarization basis defined in the co-moving frame of the spot because we expect azimuthal symmetry of the radiation in that frame. For the case of the electron-scattering dominated atmosphere the dominant direction of oscillations of the electric vector is perpendicular to the normal; we can still take PA = 0°, but assign a negative value for the PD P. Because the PD does not change along photon trajectory, the observed Stokes vector is fully defined by the flux dF_I, PD P, and PA χ:

$$\begin{array}{c}\hfill \left(\begin{array}{c}\mathrm{d}{F}_{\mathrm{I}}\\ \mathrm{d}{F}_{\mathrm{Q}}\\ \mathrm{d}{F}_{\mathrm{U}}\end{array}\right)=\mathrm{d}{F}_{\mathrm{I}}\left(\begin{array}{c}1\\ Pcos2\chi \\ Psin2\chi \end{array}\right).\end{array}$$(49)

When the spots have a finite size, we need to integrate over the spot area in the comoving frame expressed as (Nättilä & Pihajoki 2018; Salmi et al. 2018; Bogdanov et al. 2019; Suleimanov et al. 2020)

$$\begin{array}{c}\hfill \mathrm{d}{S}^{\prime }=\gamma R^2(\theta )\mathrm{d}cos\theta \mathrm{d}{\varphi }^{\prime }/cos\eta .\end{array}$$(50)

For each element of the spot we need to compute a separate light curve as a function of observed phase for the flux, PD, and PA, accounting for different time delays, and to sum them up. The total Stokes vector (denoted by the index “tot”) is obtained by integrating the Stokes vectors given by Eq. (49) over the NS surface. The corresponding total PD is then

$$\begin{array}{c}\hfill P^{\mathrm{tot}}=\frac{\sqrt{{(F_{\mathrm{Q}}^{\mathrm{tot}})}^2+{(F_{\mathrm{U}}^{\mathrm{tot}})}^2}}{F_{\mathrm{I}}^{\mathrm{tot}}},\end{array}$$(51)

where $F_{\mathrm{I}}^{\mathrm{tot}}$ represents the total flux, and the PA of the total Stokes vector χ^tot can be computed from

$$\begin{array}{c}\hfill cos2{\chi }^{\mathrm{tot}}=\frac{F_{\mathrm{Q}}^{\mathrm{tot}}}{F_{\mathrm{I}}^{\mathrm{tot}}},sin2{\chi }^{\mathrm{tot}}=\frac{F_{\mathrm{U}}^{\mathrm{tot}}}{F_{\mathrm{I}}^{\mathrm{tot}}}.\end{array}$$(52)

3. Results

3.1. Comparing PA for oblate and spherical stars

We begin by comparing the analytic expressions for the PA in Eqs. (47) and (48) in Figs. 2–4, for different co-latitudes, inclinations, and spin frequencies, and show them as functions of the rotation phase ϕ. The resulting pulse profiles look very similar if we present them as functions of the observed phase φ. For a spherical star we took the radius R₀ = 12 km and mass M₀ = 1.4 M_⊙. For an oblate star we computed the gravitating mass and equatorial radius using Eqs. (5) and (6).

Fig. 2. Top panel: polarization angles computed for the spherical (dashed) and the oblate (solid) star as a function of the rotation phase ϕ for different spot co-latitudes θ = 15° (red), 30° (green), 45° (blue), and 60° (purple). Second panel: relativistic rotation for the spherical (dashed) and the oblate (solid) stars. Third panel: component χ1 for the oblate star. Bottom panel: difference between the PAs for the oblate and spherical stars. Parameters of the non-rotating NS are R0 = 12 km and M0 = 1.4 M⊙, while in the oblate case of a NS rotating at a rate of ν* = 600 Hz, the mass and radius are computed according to Eqs. (6) and (5). The inclination is i = 40°.

Fig. 3. Same as Fig. 2 (only the top and bottom panels), but for different inclinations i = 15° (red), 30° (green), 45° (blue), and 60° (purple). Here the spot co-latitude is θ = 35°.

Fig. 4. Same as Fig. 3, but for different frequencies ν* = 1 (red), 200 (green), 400 (blue), and 600 Hz (purple). Here the spot co-latitude is θ = 40° and inclination is i = 50°.

Fig. 5. Comparison of Stokes parameters and PA profiles for our analytical OS approximation against the numerical OS approximation calculated with ARCMANCER. Three upper panels: flux (blue curve), two Stokes parameters (FQ in red and FU in orange), and residuals in all three Stokes parameters. All the corresponding results obtained using ARCMANCER are shown with black dashed curves that overlap almost exactly the colored lines. Two lower panels: PA profile and residuals between the results of the two codes. Our results are shown by a blue curve and the results of ARCMANCER with an almost exactly overlapping black curve. The profiles are computed for one spot at a NS with ν* = 600 Hz, Re = 12 km, M = 1.4 M⊙, θ = 60°, i = 40°, ρ = 10°, Tbb = 1 keV, and E = 5 keV.

From Figs. 2 and 3 we see that deviation between the spherical and oblate PA profiles (for ν_* = 600 Hz) depends strongly on the combination of the angles i and θ and on the phase ϕ. The difference is the smallest (less than about 10 degrees) when either θ ≪ i or θ ≫ i. However, for θ ≈ i we see a large deviation in the vicinity of ϕ ≈ 0, where the total angle χ changes the sign. In these situations the dashed and solid lines (for spherical and oblate stars, respectively) have opposite monotonicities. The explanation is that for the oblate case the condition of λ = θ − η < i can be satisfied even though θ > i for a spherical star. For example, at a spin rate of ν_* = 600 Hz, the angle η can be as large as 7° at θ = 45°. Nonetheless, we note that the PD is also the smallest at ϕ ≈ 0, since the direction towards the observer is closest to the local normal to the surface. The small difference between angles θ and i leads to the closest alignment between the normal and the trajectory of the light on these phases, and thus the high difference in the PA may still be difficult to observe.

The precise rotation phase ϕ, where χ changes sign, is slightly smaller than zero because χ′ reaches the extremum at an earlier phase where the denominator of Eq. (45) reaches the minimum. For the oblate case that phase is generally smaller than the corresponding phase for the spherical star, partly because χ′ is generally larger than ${\chi }_{\mathrm{sph}}^{\prime }$ close to ϕ = 0, as can be seen from Eqs. (45) and (46) and the second panel of Fig. 2. The angle χ₁ is zero at ϕ = 0. Around phases ϕ ≈ ±π/2 the correction from χ₁ dominates the differences in the total χ. Thus, for instance, the PA for the oblate NS is closer to zero in comparison to its spherical counterpart for those cases when the PA does not go the full circle (i.e., when i > λ). This is a direct consequence of the oblate shape of the star: the angle χ₁ defined by Eq. (37) has a different sign than χ₀ in Eq. (31). In addition, as can be seen from Fig. 4, the difference in the total PA χ gets larger for faster rotating NSs, and it can reach ten degrees at about 400 Hz.

3.2. Comparison with polarized light tracing

Let us now compare our results of polarized light curves to those computed with the general relativistic polarized radiative transfer code ARCMANCER (Pihajoki et al. 2018). In that code, the PA profile is calculated numerically along the geodesics using a user-specified metric. The method relies on numerically propagating a given set of coordinate vectors along the geodesic starting from the observer’s image plane until the photon trajectory intersects with the NS surface. This comparison was used to verify that our analytical formulae include all the necessary special and general relativistic sources of polarization plane rotation.

We found good agreement between the results of the codes when using the OS approximation for both codes as shown in Fig. 5. The PA profiles and the Stokes parameters were calculated for one hotspot with ν_* = 600 Hz, R_e = 12 km, M = 1.4 M_⊙, θ = 60°, i = 40°, angular radius of the spot ρ = 10°, temperature of the spot T_bb = 1 keV, and photons observed at E = 5 keV energy. We observe a difference in the Stokes parameters of about 1% and less than $0\stackrel{°}{.}2$ difference in the PA. However, we note that the polarization profiles of ARCMANCER are more sensitive to the distance of the image plane than in the case of flux (comparisons between fluxes were made in Nättilä & Pihajoki 2018). Therefore, to avoid systematic biases in the residuals, we had to increase the image plane distance from 100, which was suitable for light curve comparison, to more than 1000 r_s.

We also compared the results for the PA profiles computed in the OS approximation, spherical S+D approximation, and in the Butterworth & Ipser (1976) metric, around NS as described in AlGendy & Morsink (2014, AGM), where the shape of the star and the space-time itself are no longer spherical (see Fig. 6). We found that our analytical OS model matches the AGM model computed with ARCMANCER much more closely than the spherical S+D model. Most of the polarization discrepancies are removed by correcting the shape of the star, and only small deviations due to the space-time quadrupole deformations are left, as noted previously by Pihajoki et al. (2018).

Fig. 6. Stokes parameters and PA profiles for our analytical OS model (blue curve), the exact result of ARCMANCER with AGM model (dashed black curve), and the S+D model from Viironen & Poutanen (2004; green curve). The panels, from top to bottom, show the flux normalized to the maximum, the normalized Stokes parameters FQ and FU, the PA, and the difference in PA (χAGM − χOS in red and χAGM − χS + D in orange), respectively. The blue and black curves almost completely overlap. The parameters used to calculate the profiles are the same as in Fig. 5, except for a spherical star we use the mass and radius of the corresponding non-rotating NS of the same baryonic mass: M0 = 1.394 M⊙ and R0 = 11.43 km (see Sect. 2.1).

3.3. Determining geometrical parameters from polarization data

New constraints on the geometrical parameters (and thus also on the mass and radius) of the accreting millisecond pulsars are expected to be obtained with the upcoming polarimeters and X-ray missions like IXPE and eXTP. Therefore, we made simple estimates to predict how accurately the main geometrical parameters, namely the observer inclination i and the spot co-latitude θ, can be measured.

Using our model for polarized radiation from two antipodal hotspots (presented in Sect. 2), we produced synthetic PA profile data, assuming a fixed uncertainty of the PA in all the phase bins. In addition, we considered profiles only at 5 keV energy because the PA profile is almost independent of energy. The PA profiles were fitted using an affine invariant ensemble sampler with EMCEE (Foreman-Mackey et al. 2013) package implemented in PYTHON. We kept the mass M, the equatorial radius R_e, the observer inclination i, the spot co-latitude θ, and the angular radius of the spot ρ as free parameters, and assumed either a constant 15° or a 2° standard deviation in the measured PA. The latter error is achievable, for instance, with eXTP in the energy range 2–10 keV, assuming brightness of 0.1 Crab, 4% PD, 1 Ms exposure, and a negligible amount of background counts. The 15° error would correspond to a 1 ks exposure but otherwise the same assumptions. In the case of IXPE the errors are expected to be the same order of magnitude.

Instead of fitting PA profiles, it is preferable to use Stokes parameters because they are additive and their distribution is very well described by a bivariate normal distribution, even though they are neither the directly observed quantities from X-ray polarimeters (see Kislat et al. 2015). For simplicity, here we assumed that the errors in the measured PA would also be normally distributed, which is not true in general. In actual observations the error in PA will also be much higher at those pulsar phases when the polarization degree is closer or below the detection limit (which happens when the flux is at the maximum). A more accurate study of simulated observations, using Stokes parameters and the energy response of the IXPE detector, will be made in a forthcoming publication (Salmi et al. 2020).

The synthetic data were created assuming a rapidly rotating star (ν_* = 600 Hz to emphasize the difference between an oblate and a spherical star) with M = 1.4 M_⊙, R_e = 12 km, i = 40°, θ = 60°, and ρ = 10°. When fitting the data, we also assumed a very conservative time resolution and used only ten phase bins. We treated the phase shift as a nuisance parameter. Instead of marginalized likelihoods over all phase shifts, we used a bisection method to maximize the likelihood, which is much faster, while the difference between these two methods is negligible (Salmi et al. 2018).

The results of these fits are shown in Figs. 7 and 8. The credible intervals are also presented in Table 1. As expected, we found a significant improvement in the parameter constraints assuming a smaller error in PA (i.e., longer exposure time of the detector). The 68% credible interval for the inclination angle i changes from 23°– 50° to 38°– 43° and for spot co-latitude from 42°– 70° to 57°– 64°. In the case of short exposure, the mass and radius are not very well constrained. However, in the case of a long exposure we also get limits for radius $R_{\mathrm{e}}=12.5_{-1.3(2.4)}^{+1.5(3.0)}$ km and even for mass $M=1.{39}_{-0.22(0.36)}^{+0.22(0.43)}{M}_{\odot }$ with 68% (95%) limits. This is mainly due to the influence of mass and radius on the moment when the second spot becomes visible and starts to contribute to the PA profile. The size of the spot ρ is also not very well constrained, although with a long exposure we get an upper limit of 22° with 95% certainty. These limits should be taken as estimates because the model describes the synthetic data perfectly, which is not the case with real observations. We also note that much tighter constraints for the mass and radius are expected when combining polarization measurements with the modeling of the pulse profiles, i.e., flux (which is not considered here). The most essential benefit of polarization data is expected to be in breaking the degeneracy in i and θ (Viironen & Poutanen 2004).

Fig. 7. Posterior distributions for an oblate fit to an oblate star assuming 15° standard deviation in the measured PA. In the two-dimensional posterior distributions the 68% (in black) and the 95% (in gray) highest posterior density credible regions are shown. In the one-dimensional posterior distributions the dashed red contour shows a 68% and the dashed orange contour a 95% highest posterior density credible interval. The blue lines show the input values.

Fig. 8. Posterior distributions for an oblate fit to an oblate star assuming 2° standard deviation in the measured PA. The colors, contours, and other symbols are the same as in Fig. 7.

3.4. Effects of oblateness on the polarization

To test the effect of oblateness on the PA, we also fitted the synthetic PA profile (for oblate star) using an S+D model for a spherical star. The free parameters of our model were the same as in Sect. 3.3. The synthetic PA profile is shown in Fig. 9 with a blue curve, and the corresponding profile with a spherical star is shown with a green curve (the former being the same profile used in Sect. 3.3). As in Fig. 6, we see again evident differences between the two models.

Fig. 9. Flux (top panel), PD (middle panel), and PA (bottom panel) as a function of observed phase for an oblate (in blue) and spherical (in green) NS with identical model parameters at 5 keV (given in Sect. 3.3). For the spherical star corresponding non-rotating mass and radius are used: M0 = 1.394 M⊙ and R0 = 11.43 km.

The results of the fitting are shown in Fig. 10 (spherical fit to an oblate star) for a 1σ error of 2° in the measured PA, and they can be compared with those obtained using the correct model with the same exposure time in Fig. 8. The credible intervals are also found from Table 1. The results show no bias when the correct model is used, but with an incorrect model for the star shape we see an obvious bias in the spot co-latitude θ giving too small values for the angle. The correct value is outside of the measured 95% credible interval. This is likely explained by the co-latitude appearing in terms of λ = θ − η in PA Eq. (45). Fitting a spherical model to an oblate star tends to find a correct λ rather than correct θ. It shows that the effects of oblateness are important to include when calculating the PA profiles, at least for stars rotating as fast as ν_* = 600 Hz.

Fig. 10. Posterior distributions for a spherical fit to an oblate star assuming 2° standard deviation in the measured PA. The colors, contours, and other symbols are the same as in Fig. 7.

In addition, we see a small but noticeable shift in the measured radius R_e towards higher values (even though R_e should be smaller for a corresponding non-rotating star as seen in Sect. 2.1). However, the correct value is still inside the 68% credible interval. A minor dependence of the PA profile on R_e is not surprising because the radius affects the time interval when the second spot is visible. The visibility condition is also different for oblate and spherical stars, and thus it can be responsible for biasing the constraints in the radius.

4. Conclusions

We developed the formalism to compute the light curves of Stokes parameters observed from a rapidly rotating NS in the OS approximation. We derived analytical formulae describing the rotation of the polarization plane as a result of special and general relativistic effects accounting for the oblate shape of a NS. The expression for the PA is summarized in Eqs. (32), (37), (45), and (47). In the case of a spherical star, our results simplify to those presented in Viironen & Poutanen (2004) and Poutanen (2020a). We also verified the accuracy and correctness of our formalism by comparing the resulting polarization pulse profiles to numerical calculations computed with the general relativistic polarized ray-tracing code ARCMANCER.

In addition, we estimated that the accuracy of a few degrees in the measured observer inclination and spot co-latitude can be obtained from the upcoming X-ray polarization measurements, assuming a correct theoretical model and using data with high precision. We caution that the resulting system parameters can be seriously biased if a spherical star shape is erroneously used to compute the profiles emitted by an oblate NS.

We also note that we only considered the simple case with two antipodal hotspots assuming a dipolar magnetic field configuration. The presented model, however, can also be applied to more complicated spot configurations, such as recently inferred in MSPs based on the observations of PSR J0030+0451 by the Neutron star Interior Composition ExploreR (Bilous et al. 2019; Miller et al. 2019; Riley et al. 2019). In that case, the polarized pulse profiles would show a more complicated structure. Because the non-antipodal spots have largely different PAs, the observed total PA would show abrupt breaks when any of the spots eclipses. Future polarimetric observations will thus be able to test the magnetic field geometry in MSPs.

Acknowledgments

This research was supported by the Magnus Ehrnrooth Foundation, the University of Turku Graduate School in Physical and Chemical Sciences (T. S.), the Academy of Finland grants 322779 (J. P.) and 333112 (V. L., J. P.), and the Russian Science Foundation grant 20-12-00364 (V. L., J. P.). The computer resources of the Finnish IT Center for Science (CSC) and the FGCI project (Finland) are acknowledged.

References

Appendix A: Alternative derivation for the general relativistic correction of PA

The general relativistic correction of PA for the Schwarzschild metric can also be derived using the common vector formalism presented in Poutanen (2020a). Here we reproduce the necessary equations from their Sect. 3.2.2, but considering an oblate NS instead of a spherical one. The common vector perpendicular to both photon momentum and the spot velocity vector in both frames is

$$\begin{array}{cc}& \widehat{\mathit{M}}=\frac{{\widehat{\mathit{k}}}_0\times \widehat{\mathit{\beta }}}{|{\widehat{\mathit{k}}}_0\times \widehat{\mathit{\beta }}|}=\frac{{\widehat{\mathit{k}}}_0^{\prime }\times \widehat{\mathit{\beta }}}{|{\widehat{\mathit{k}}}_0^{\prime }\times \widehat{\mathit{\beta }}|}\hfill \\ \hfill & =\frac{1}{sin\xi sin\psi }\left(\begin{array}{c}-cos\varphi [cosisin\alpha +cos\theta sin(\psi -\alpha )]\\ -sin\varphi [cosisin\alpha +cos\theta sin(\psi -\alpha )]\\ sinicos\varphi sin\alpha +sin\theta sin(\psi -\alpha )\end{array}\right).\hfill \end{array}$$(A.1)

The unit vector of the photon momentum in the spot comoving frame comes from the Lorentz transformation:

$$\begin{array}{c}\hfill {\widehat{\mathit{k}}}_0^{\prime }=\delta [{\widehat{\mathit{k}}}_0-\gamma \beta \widehat{\mathit{\beta }}+(\gamma -1)\widehat{\mathit{\beta }}(\widehat{\mathit{\beta }}·{\widehat{\mathit{k}}}_0)].\end{array}$$(A.2)

These expressions do not depend on the NS shape. The polarization basis related to photon momentum in the comoving frame $\widehat{\mathit{k}}{\prime }_0$ and the surface normal is (for an oblate star)

$$\begin{array}{c}\hfill {\widehat{\mathit{e}}}_1^{\prime 0}=\frac{\widehat{\mathit{n}}-cos{\sigma }^{\prime }{\widehat{\mathit{k}}}_0^{\prime }}{sin{\sigma }^{\prime }},{\widehat{\mathit{e}}}_2^{\prime 0}=\frac{{\widehat{\mathit{k}}}_0^{\prime }\times \widehat{\mathit{n}}}{sin{\sigma }^{\prime }}.\end{array}$$(A.3)

The angle between vector $\widehat{\mathit{M}}$ and the basis vector ${\widehat{\mathit{e}}}_1^{\prime 0}$ is

$$\begin{array}{cc}& cos{\chi }_M^{\prime }=\widehat{\mathit{M}}·{\widehat{\mathit{e}}}_1^{\prime 0}=\frac{C}{sin\xi \sqrt{1-{\delta }^2{cos}^2}\sigma },\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill & sin{\chi }_M^{\prime }=-\widehat{\mathit{M}}·{\widehat{\mathit{e}}}_2^{\prime 0}=\frac{cos\sigma }{sin\xi \sqrt{1-{\delta }^2{cos}^2\sigma }}\frac{cos\xi -\beta }{1-\beta cos\xi },\hfill \end{array}$$(A.5)

where

$$\begin{array}{c}\hfill C=-\frac{sin\alpha }{sin\psi }cos\zeta +\frac{sin(\psi -\alpha )}{sin\psi }sin\eta .\end{array}$$(A.6)

Therefore, we get

$$\begin{array}{c}\hfill tan{\chi }_M^{\prime }=\frac{cos\sigma }{C}\frac{cos\xi -\beta }{1-\beta cos\xi }.\end{array}$$(A.7)

The polarization basis in the laboratory frame is

$$\begin{array}{c}\hfill {\widehat{\mathit{e}}}_1^0=\frac{\widehat{\mathit{n}}-cos\sigma {\widehat{\mathit{k}}}_0}{sin\sigma },{\widehat{\mathit{e}}}_2^0=\frac{{\widehat{\mathit{k}}}_0\times \widehat{\mathit{n}}}{sin\sigma }.\end{array}$$(A.8)

Now we can evaluate the angle χ_M between the common vector $\widehat{\mathit{M}}$ and the projection of the stellar normal on the sky as viewed in that frame:

$$\begin{array}{cc}& cos{\chi }_M={\widehat{\mathit{e}}}_1^0·\widehat{\mathit{M}}=\frac{C}{sin\sigma sin\xi },\hfill \end{array}$$(A.9)

$$\begin{array}{cc}\hfill & sin{\chi }_M={\widehat{\mathit{e}}}_2^0·\widehat{\mathit{M}}=-\frac{cos\sigma cos\xi }{sin\sigma sin\xi },\hfill \end{array}$$(A.10)

and

$$\begin{array}{c}\hfill tan{\chi }_M=-\frac{cos\sigma cos\xi }{C}.\end{array}$$(A.11)

Thus, the total general relativistic correction angle ${\chi }^{\mathrm{\prime }}={\chi }_M^{\mathrm{\prime }}+{\chi }_M$ can be found as

$$\begin{array}{cc}& tan{\chi }^{\prime }=\frac{tan{\chi }_M^{\prime }+tan{\chi }_M}{1-tan{\chi }_M^{\prime }tan{\chi }_M}\hfill \\ \hfill & =-\frac{\beta cos\sigma C{sin}^2\xi }{C^2+{cos}^2\sigma {cos}^2\xi -\beta cos\xi (C^2+{cos}^2\sigma )}.\hfill \end{array}$$(A.12)

It is easy to show that C² + cos²σ = sin²ξ, so we finally get

$$\begin{array}{c}\hfill tan{\chi }^{\prime }=-\frac{\beta cos\sigma C}{{sin}^2\sigma -\beta cos\xi },\end{array}$$(A.13)

which is equivalent to Eq. (45).

All Tables

All Figures

Fig. 1. Geometry of the oblate star and the lab frame coordinate system xyz. A spot at the NS surface at co-latitude θ has a rotation phase ϕ from the observer vector $\widehat{\mathit{k}}$ around rotational axis z. Three unit vectors are defined for the spot: $\widehat{\mathit{n}}$ is normal to the surface, $\widehat{\mathit{r}}$ is the radius vector, and ${\widehat{\mathit{k}}}_0$ is the direction of emitted photon propagation.

In the text

Fig. 2. Top panel: polarization angles computed for the spherical (dashed) and the oblate (solid) star as a function of the rotation phase ϕ for different spot co-latitudes θ = 15° (red), 30° (green), 45° (blue), and 60° (purple). Second panel: relativistic rotation for the spherical (dashed) and the oblate (solid) stars. Third panel: component χ1 for the oblate star. Bottom panel: difference between the PAs for the oblate and spherical stars. Parameters of the non-rotating NS are R0 = 12 km and M0 = 1.4 M⊙, while in the oblate case of a NS rotating at a rate of ν* = 600 Hz, the mass and radius are computed according to Eqs. (6) and (5). The inclination is i = 40°.

In the text

	Fig. 3. Same as Fig. 2 (only the top and bottom panels), but for different inclinations i = 15° (red), 30° (green), 45° (blue), and 60° (purple). Here the spot co-latitude is θ = 35°.
In the text

	Fig. 4. Same as Fig. 3, but for different frequencies ν* = 1 (red), 200 (green), 400 (blue), and 600 Hz (purple). Here the spot co-latitude is θ = 40° and inclination is i = 50°.
In the text

Fig. 5. Comparison of Stokes parameters and PA profiles for our analytical OS approximation against the numerical OS approximation calculated with ARCMANCER. Three upper panels: flux (blue curve), two Stokes parameters (FQ in red and FU in orange), and residuals in all three Stokes parameters. All the corresponding results obtained using ARCMANCER are shown with black dashed curves that overlap almost exactly the colored lines. Two lower panels: PA profile and residuals between the results of the two codes. Our results are shown by a blue curve and the results of ARCMANCER with an almost exactly overlapping black curve. The profiles are computed for one spot at a NS with ν* = 600 Hz, Re = 12 km, M = 1.4 M⊙, θ = 60°, i = 40°, ρ = 10°, Tbb = 1 keV, and E = 5 keV.

In the text

Fig. 6. Stokes parameters and PA profiles for our analytical OS model (blue curve), the exact result of ARCMANCER with AGM model (dashed black curve), and the S+D model from Viironen & Poutanen (2004; green curve). The panels, from top to bottom, show the flux normalized to the maximum, the normalized Stokes parameters FQ and FU, the PA, and the difference in PA (χAGM − χOS in red and χAGM − χS + D in orange), respectively. The blue and black curves almost completely overlap. The parameters used to calculate the profiles are the same as in Fig. 5, except for a spherical star we use the mass and radius of the corresponding non-rotating NS of the same baryonic mass: M0 = 1.394 M⊙ and R0 = 11.43 km (see Sect. 2.1).

In the text

Fig. 7. Posterior distributions for an oblate fit to an oblate star assuming 15° standard deviation in the measured PA. In the two-dimensional posterior distributions the 68% (in black) and the 95% (in gray) highest posterior density credible regions are shown. In the one-dimensional posterior distributions the dashed red contour shows a 68% and the dashed orange contour a 95% highest posterior density credible interval. The blue lines show the input values.

In the text

	Fig. 8. Posterior distributions for an oblate fit to an oblate star assuming 2° standard deviation in the measured PA. The colors, contours, and other symbols are the same as in Fig. 7.
In the text

	Fig. 9. Flux (top panel), PD (middle panel), and PA (bottom panel) as a function of observed phase for an oblate (in blue) and spherical (in green) NS with identical model parameters at 5 keV (given in Sect. 3.3). For the spherical star corresponding non-rotating mass and radius are used: M0 = 1.394 M⊙ and R0 = 11.43 km.
In the text

	Fig. 10. Posterior distributions for a spherical fit to an oblate star assuming 2° standard deviation in the measured PA. The colors, contours, and other symbols are the same as in Fig. 7.
In the text