© ESO 2020

1. Introduction

Thermonuclear explosions (flashes) of the freshly accreting matter on the surface of neutron stars (NSs) in low-mass X-ray binaries (LMXBs) are observed as type I X-ray bursts (Lewin et al. 1993; Strohmayer & Bildsten 2006). The most successful X-ray observatory for X-ray burst studies was Rossi X-ray Timing Explorer (RXTE). The large area and the high time resolution of its Proportional Counter Array (PCA; Jahoda et al. 2006) allowed investigation of the X-ray burst flux variability with unprecedented accuracy (Galloway et al. 2008). Thanks to RXTE/PCA, it was also possible to perform accurate measurements of X-ray burst spectral evolution. The burst spectra, usually well fitted with a blackbody, give a detailed view of the evolution of the main parameters, such as the blackbody temperature T_BB, its normalization K, and the blackbody flux F_BB (see e.g. Galloway et al. 2008). However, local spectra of hot NS atmospheres are not actually blackbodies, but are rather close to a diluted blackbody with the surface flux proportional to the Planck function wπB_E(T_c). The colour temperature of the spectrum T_c is typically higher than the effective temperature T_eff by the colour-correction factor f_c = T_c/T_eff >  1. As a result, the dilution factor has to be smaller than unity to conserve the bolometric flux. The observed blackbody temperature T_BB = T_c/(1 + z) is just a gravitationally redshifted colour temperature. The blackbody normalization K of the observed burst spectrum is proportional to the dilution factor w characterizing properties of radiation escaping from the NS surface.

The dilution and colour-correction factors depend on the NS luminosity relative to the Eddington luminosity, and change most significantly near the Eddington limit. Therefore, the most powerful X-ray bursts with luminosities exceeding the Eddington luminosity, so-called photospheric radius expansion (PRE) bursts, have to demonstrate a significant evolution of the blackbody normalization after the moment when the photosphere settles down at the NS surface (the touchdown point). Indeed, this kind of evolution has been observed in some X-ray bursts and the cooling tail method was designed to obtain the apparent NS size and the Eddington flux from the data (Suleimanov et al. 2011a; Poutanen et al. 2014). They were then converted to constraints on the NS basic parameters such as mass M and radius R. However, due to a divergence of the Jacobian the resulting distribution of M and R is biased (Özel & Psaltis 2015). This problem was avoided by fitting the evolution of the blackbody normalization during the cooling tail using M and R as parameters (Nättilä et al. 2016; Suleimanov et al. 2017a). A more detailed description, the history, and applications of the method are reviewed in Suleimanov et al. (2016) and Degenaar et al. (2018). Further development of the method is the direct fitting of the observed spectral sequence with the model spectra of hot NS atmospheres. It allowed us to obtain NS radii with unprecedented accuracy of 0.5 km (Nättilä et al. 2017). Such a narrow range of allowed NS radii already gives interesting constraints on the equation of state (EoS) of cold dense matter (see the recent review by Lattimer & Prakash 2016).

The studies quoted above assume that the NS is spherical and its emission is isotropic. However, the discovery by RXTE of burst oscillations (Strohmayer et al. 1996) and accreting millisecond pulsars (Wijnands & van der Klis 1998) imply that NSs in LMXBs are rapidly rotating at a rate in the range ν_* ≈ 200−600 Hz (Watts 2012; Patruno & Watts 2012; Campana et al. 2018). Rapid rotation distorts the NS shape and produces a potentially important systematic effect on the NS radii determination. Radiation from the approaching side of the NS becomes Doppler boosted shifting the emission peak to higher energies as, for example, was extensively discussed in the studies of pulse profiles from millisecond pulsars using the so-called Schwarzschild+Doppler (S+D) (Poutanen & Gierliński 2003; Poutanen & Beloborodov 2006) and oblate Schwarzschild (OS) approximations (Morsink et al. 2007; Bogdanov et al. 2019). Rotation also affects the value of the local gravitational acceleration and therefore the value of the flux in Eddington units and as a consequence the escaping spectrum and the colour correction.

Models for rotating NS have been constructed for a rather long time, with the current state being presented by Paschalidis & Stergioulas (2017). The accurate shape of a rotating NS depends on the assumed EoS, but for the astrophysical applications some approximate rotating NS models as suggested by Morsink et al. (2007) and AlGendy & Morsink (2014) are sufficient. The influence of the oblate form of the rotating NSs on the pulse profiles of millisecond pulsars was considered using accurate general relativistic ray-tracing in the curved space-time (see e.g. Cadeau et al. 2007; Nättilä & Pihajoki 2018). The effects of rapid rotation on the observed spectra and the line profiles were studied by Bauböck et al. (2013, 2015). In particular, Bauböck et al. (2015) computed the angle-averaged correction to the measured NS radius as a function of the rotation frequency assuming local isotropic blackbody emission of the same temperature. Vincent et al. (2018), on the other hand, used the local model atmosphere spectra from Majczyna et al. (2005), but assumed no latitudinal variation. We note that a modest change in the gravitational acceleration leads to large change in the colour-correction factor, especially close to the Eddington limit. Thus, we suggest that in order to compute realistic observed X-ray burst spectra, latitudinal variations of the local atmosphere spectra (e.g. in terms of variations of the colour-correction factor) should be accounted for.

In this paper we suggest a new method for computing spectra from rapidly rotating NSs. We base our approach on a slow-rotation approximation for the NS shape and space-time metric proposed by AlGendy & Morsink (2014). We introduce the modified oblate Schwarzschild (MOS) approximation where light bending is computed in Schwarzschild metric, but for calculations of the gravitational redshift and Doppler boost we approximately account for the frame dragging and quadrupole moment of a NS. In Sect. 2 we propose analytical formulae connecting the gravitational mass and equatorial radius of a rotating NS to the mass and radius of a non-rotating NS of the same baryonic mass. We then compute the spectra of the blackbody emitting, rapidly rotating NSs and obtain the colour-correction and dilution factors as a function of rotation rate. We extend the direct cooling tail method to rapidly rotating NSs in Sect. 3. The modified method is applied to the X-ray burst from LMXB SAX J1810.8−2609 in Sect. 4. We summarize the results in Sect. 5. Technical details are presented in the Appendices.

2. Observational appearance of rapidly rotating NSs

One of the goals of the current work is to develop the modified direct cooling tail method applicable to rapidly rotating NSs. The final result of the method would be an estimation of the equivalent radius of a non-rotating NS of a given mass using observed data on the X-ray bursts from the surface of a rapidly rotating NS. Two additional parameters appearing in this task are the rotational angular velocity Ω_* = 2πν_* of the NS and the inclination angle i between the rotation axis and the line of sight to the observer. The exact model of the rapidly rotating NS depends on the assumed EoS in its inner core, which is not well known. Moreover, constraining it from the observations is one of our goals. Fortunately, the shape of a rapidly rotating NS depends on the actual EoS only slightly, if the rotation is not very rapid, ν_* <  700 Hz. In this case, a slow-rotation approximation can be used, and the model of a rapidly rotating NS depends only on its gravitational mass M, equatorial radius R_e, and the angular velocity Ω_*, as discussed by Morsink et al. (2007) and AlGendy & Morsink (2014). The approach suggested in these works is used in this paper (see Appendix A for details). We describe now the effects arising due to rapid rotation.

2.1. Connection between basic parameters of rotating and non-rotating NSs

The first step in understanding the impact of rapid rotation is to compute the shape of a NS rotating with a given observed frequency ν_*. The most important input parameter is the equatorial radius R_e (see Appendix A). We would like to find a connection of the radius R and the mass M of a non-rotating NS with the equatorial radius R_e and gravitational mass M′ of a NS with the same baryonic mass , rotating at a rate ν_*. Cook et al. (1994) computed models of rapidly rotating NSs in general relativity for several EoS and a few values of NS masses, including 1.4 M_⊙ and the maximum possible NS mass for a given EoS. The ratio of R_e to the radius of the non-rotating configuration R computed from these models is shown in Fig. 1 as a function of the relative rotation frequency , where

(1)

Fig. 1. Dependences of the ratio Re/R on the relative rotational frequency for numerical NS models with various EoS (Cook et al. 1994) for the maximum possible mass (top panel) and the mass M = 1.4 M⊙ (bottom panel). Red solid curves show the fits with Eq. (2) for the masses M = 2.5 M⊙ (top panel) and 1.4 M⊙ (bottom panel). The fitting curve for M = 1.4 M⊙ is shown in the top panel as dashes. The relative errors of the fits are also shown in the narrow additional panels below the main panels.

is the maximum possible rotation frequency for a given non-rotating NS mass and radius (Haensel et al. 2009) and we defined M_1.4 = M/1.4 M_⊙. The tabulated values can be well represented by an approximate formula

(2)

which is shown in Fig. 1 by red curves. The relative fitting errors are also shown in the additional panels. The maximum errors of fitting using Eq. (2) are not larger than 5% for the most massive and rapidly rotating NSs. We note that we use the extreme EoSs A and L from Cook et al. (1994), which bracket the most probable modern EoS (see e.g. Nättilä et al. 2016), giving us a wider range of values to test our formulae.

The mass of the rotating configuration is also affected. It is obvious that the gravitational mass of the NS is not a simple sum of the baryon masses , but is reduced roughly as due to a deep gravitational well. Because the radius of the rotating NS is larger than that of a non-rotating NS, the gravitational correction is reduced and the gravitational mass increases. The ratio of the gravitational mass of a rotating star M′ to the mass of a non-rotating star M increases with the rotation rate as shown in Fig. 2 based on calculations of Cook et al. (1994). The tabulated values can be fitted by a relation (see red solid curves in Fig. 2)

(3)

Fig. 2. Dependences of ratios M′/M on the relative rotational frequency for numerical NS models with various EoS (Cook et al. 1994) for the maximum possible mass (top panel) and the mass M = 1.4 M⊙ (bottom panel). Red solid curves show the fits with Eq. (3) for the masses M = 2.5 M⊙ (top panel) and 1.4 M⊙ (bottom panel). The fitting curve for M = 1.4 M⊙ is shown in the top panel as dashes. The relative errors of the fits are also shown in the narrow additional panels below the main panels.

where the fitting coefficients are a₀ = 1 − a₁/1.1, , and a₂ = 10a₁. The accuracy of Eq. (3) is typically better than 1% (see additional panels in Fig. 2). We also note that relations (2) and (3) are valid for the NSs of fixed baryonic mass .

2.2. Apparent area of a rotating NS

The centrifugal force deforms the shape of a rapidly rotating NS, which becomes flattened at the poles and extended at the equator. This means that the projection area on the sky of a rapidly rotating NS S(ν_*, i) is different from the visible area of the non-rotating NS of the same baryonic mass S(ν_* = 0) = πR²(1 + z)², and it depends on the inclination angle i. The apparent area of the given family of models of the rapidly rotating NSs of the same baryonic mass is computed via Eq. (B.34) using the specific intensity and the energy ratio equal to one, and the approach described further in Appendix B. The relation between the equatorial radius R_e and the radius of the corresponding non-rotating NS is given by Eq. (2), while the masses are related via Eq. (3). The NS shape is described by approximate formulae from AlGendy & Morsink (2014), see Appendix A. The relativistic computations of the ratio S(ν_*, i)/S(ν_* = 0) for a few inclination angles i are presented in the top panel of Fig. 3. We note that the apparent area is larger for low inclination and the deformation is smaller for more compact stars. Because this ratio is greater than unity, in order to obtain the apparent area of a non-rotating NS, the observed apparent area has to be reduced by a factor depending on the assumed or measured ν_* and i.

Fig. 3. Top panel: ratio of the apparent area of a rapidly rotating NS to that of the non-rotating star (of the same baryonic mass) as a function of the rotational frequency for inclination angles i= 0°, 45°, and 90° is shown by black, red, and blue curves, respectively. The mass of a non-rotating star was assumed M = 1.5 M⊙ and two radii R = 14 (solid curves) and 12 km (dashed curves) are considered. Bottom panel: apparent area of a rapidly rotating NS for the same family of models as in the top panels divided by the apparent area of a non-rotating NS of the same gravitational mass M′ (from Eq. (3)) and the radius equal to equatorial radius Re given by Eq. (2). The red and the blue dashed curves overlap in both panels.

The bottom panel of Fig. 3 shows the ratio of the apparent area of a rotating NS to that of a spherical star of the same equatorial radius and same gravitational mass (and thus different baryonic masses). This ratio is obviously lower than unity. Bauböck et al. (2012, see their Figs. 5 and 6) also presented similar dependences. They used a rather different approach to describe rapidly rotating NSs, but their results for the case Kerr + Obl are similar to ours. For a better comparison we computed the ratio S(ν_*, i)/S(R_e(ν_*)) for the model with parameters M′ = 1.8 M_⊙, R_e = 15 km, ν_* = 500 Hz (corresponding to M ≈ 1.79 M_⊙ and R ≈ 14.13 km) and i = 90° and got the value of 0.963. We note that Bauböck et al. (2012) used another approximation for the NS shape from Morsink et al. (2007). Using the same approximation we obtained the value of 0.971, which is very close to the corresponding value of ≈0.975 in their Fig. 5.

2.3. Emission spectrum of a NS

In the calculations presented in this paper we consider a number of different cases for the local specific intensity. The simplest model for the surface emission is that of the blackbody emission, with the specific intensity measured in the frame associated with the NS surface being

(4)

where B_E′(T) is the Planck function of local temperature T at photon energy E′ and σ′ is the zenith angle measured from the local normal to the NS surface¹. The corresponding flux in this case is

(5)

and the bolometric flux corresponding to this case is given by the Stefan–Boltzmann law F_bol = σ_SBT⁴. Another case for the angular distribution that we would like to explore is that corresponding to the electron-scattering dominated atmosphere. We keep the Planck function to describe the energy dependence of the local specific intensity and approximate the dependence on the zenith angle by a linear function of cos σ′ (Chandrasekhar & Breen 1947; Sobolev 1949):

(6)

The chosen normalization of the angular factor gives the same emergent flux as a unit constant.

Furthermore, we consider the case of a diluted blackbody corresponding to the intensity

(7)

where w is the dilution factor, f_c is the colour-correction factor, T_eff is the effective temperature, and T_c = f_cT_eff is the colour temperature. The surface flux is now

(8)

Similarly to Eq. (6) we consider a diluted blackbody with a beamed emission pattern

(9)

The surface flux is also given by Eq. (8). Equation (9) provides a good approximation to the specific intensity from the NS atmosphere models computed by Suleimanov et al. (2012); Suleimanov et al. (2017a) for NS luminosities close enough to the Eddington value (L ≳ 0.1L_Edd). These accurate models for the intensity for various effective temperatures T_eff and gravities g_eff are also used in our calculations. They are also fitted by Eq. (9) with free w and f_c, which are tabulated on a grid of T_eff and g_eff and used in the calculations instead of accurate models for simplicity.

The flux from a bursting NS as observed at Earth can also be well represented by a diluted blackbody:

(10)

Here the normalization factor K is proportional to the dilution factor w, while the observed colour temperature T_BB is related to the colour-temperature f_cT_eff as measured at the NS surface.

2.4. Apparent spectrum of a blackbody emitting rotating NS

The apparent size of the NS on the sky does not fully determine how the spectral parameters change due to rotation. Let us assume that the NS surface emits locally a blackbody of constant local temperature T = T_c = T_eff with specific intensity given by Eq. (4). The observed spectrum for a rotating NS deviates from the blackbody for two reasons. Firstly, the surface gravitational redshift varies over the surface of a rotating NS. Secondly, Doppler boosting affects the amplitude and energy of the spectral peak. The Doppler boost is insignificant for pole-on sources and reaches its maximum level at i = 90°. The change in the spectrum thus depends not only on the NS rotation rate which affects its shape, but also on the inclination angle i.

Let us construct families of rotating NS models which have the same baryonic mass, but have different rotational frequencies. We first choose the mass M = 1.5 M_⊙ and consider two different radii R = 14 km and 12 km for a non-rotating star. Then for each spin, the values for the rotating NS mass M′ and equatorial radius R_e were computed using Eqs. (2) and (3). The observed spectrum is computed following the approach described in Appendix B. The resulting spectrum F_E in the frame of the observer at infinity was fitted in the observed energy range (0.1–10)kT with the diluted blackbody wπB_E(f_cT_∞) giving the values of f_c and w. Here T_∞ = T/(1 + z) is the gravitationally redshifted local temperature, where z is the gravitational redshift of the non-rotating NS, for which f_c = 1. The dilution factor w was then divided by the value obtained for a non-rotating NS.

We plot the corresponding colour-correction factor f_c, correction to the dilution factor w, and the correction to the bolometric flux with the solid curves in Fig. 4. The top panels demonstrate that the colour temperature decreases with the spin for a face-on observer (i = 0) because the NS oblateness increases at higher rotation rates resulting in a corresponding increase in the gravitational redshift of the polar region. On the other hand, the edge-on observer sees higher temperature comparing to a non-rotating star because of the Doppler effect. The effect is greater for larger NS radii.

Fig. 4. Top panels: colour correction to the observed blackbody temperature compared to a non-rotating NS for a NS of constant local temperature as a function of the rotational frequency for inclination angles i= 0°, 45°, and 90° is shown by the black, red, and blue curves, respectively. The left panels correspond to radius of a non-rotating NS R = 14 km and the right panels to R = 12 km. The mass of the non-rotating NS was assumed M = 1.5 M⊙ and both families of models have the same baryonic mass. The solid curves are for emission pattern corresponding to the isotropic intensity (4), while the dashed curves are for the electron-scattering dominated atmosphere with the intensity given by Eq. (6). The colour-correction factor is indistinguishable for the two considered emission patterns. Middle panels: correction factor to the blackbody normalization (dilution factor) as a function of the rotational frequency. Bottom panels: the correction to the bolometric flux. Same notations as in the top panels.

The middle panels of Fig. 4 show how the blackbody normalization is affected by rotation. We see that for low inclinations the normalization increases, which results from an increase in the apparent area as discussed in Sect. 2.2. However, for high inclinations the normalization decreases slightly, which is opposite to the behaviour of the apparent area. This is a natural consequence of a large increase in f_c because of the Doppler effect. However, the bolometric flux (shown in bottom panels of Fig. 4) increases a little more slowly than resulting in a reduction of w. The bolometric flux monotonically increases with the inclination being highest for an edge-on observer. It follows very closely the product , with the deviation from that dependence (i.e. bolometric correction to the blackbody) being negligible for i = 0° and reaching only 1–2% for i = 90° at ν_* = 700 Hz.

In addition to the blackbody, we also considered the case of an electron-scattering dominated atmosphere, where the intensity given by Eq. (6) is much stronger beamed along the surface normal. The results are shown with dashed curves in Fig. 4. Interestingly, the colour-correction factor turned out to be identical to that obtained for isotropic emission. The correction to the blackbody normalization, on the other hand, deviates farther from unity than for the isotropic emission (see middle panels of Fig. 4). Also in this case, the bolometric flux increases with the NS spin. For smaller NS radii, R = 12 km (right bottom panel), its dependence on the inclination angle is similar to the isotropic case, i.e. it monotonically increases with inclination. However, for larger radii, R = 14 km, it shows the opposite behaviour, with the flux being highest for a face-on observer. We note that the dashed curve below the blue solid curve in the left bottom panel of Fig. 4 is black and corresponds to i = 0°. The reason is a stronger radiation beaming along the normal and greater deformation of the NS surface. The increase in the apparent area for a face-on observer makes a stronger effect than the increase in the Doppler boost for an edge-on observer.

At this point we can try to compare our results with those of Bauböck et al. (2015). However, there is a large difference in our approaches that complicates direct comparison. Bauböck et al. (2015) computed the corrections to the colour temperature and bolometric flux too, but compared the results obtained for a rotating star with equatorial radius R_e to those for a non-rotating star of the same radius. In this case the overall gravitational redshift is similar for the two stars, and the colour temperature changes mostly due to the Doppler effect resulting in a 1% deviation. However, we fix the mass and radius of a non-rotating star and consider a family of models as a function of spin. Thus, at high spin rates, our R_e is significantly larger and polar radius is significantly smaller than the radius R of a non-rotating configuration. For low inclinations, this results in an effectively higher gravitational redshift leading to a decrease in the colour temperature (and f_c) by 1% compared to a non-rotating case. On the other hand, at high inclination, the effective gravitational redshift is much lower resulting in an observed temperature that is higher by 7% for a rotating star. We also note that corrections to the metric due to the frame dragging and quadrupole moment play no role for spins below 500 Hz.

As an illustration we also show in Fig. 5 the images of NSs rotating at 700 Hz (see also Vincent et al. 2018). The colour-coding and the white lines are the contours of constant observed bolometric intensity. It is equal to the product of the intensity as measured at the NS surface I(σ′) and the fourth power of the total redshift (i.e. the product of Doppler factor and the surface gravitational redshift, given by Eq. (B.35)). We consider the same two cases of an isotropically emitting star (see upper panels of Fig. 5) and an electron-scattering dominated atmosphere (lower panels of Fig. 5). The results are normalized by the factor σ_SBT⁴/π. The face-on image of an isotropically emitting star (left upper panel) shows the minimum at the pole because the redshift there is highest. The intensity is also low at the equator because of the transverse Doppler effect. At high inclinations, the contours of constant intensity nearly follow the lines of constant surface radial (projected) velocity which are nearly straight vertical lines. The highest and the lowest intensities are reached at the edges of the image, where the absolute values of the projected velocities are highest. The stronger beaming of radiation along the normal to the surface in the second case dominates over the reduction of intensity due to higher redshift, resulting in a peak of emission at the pole for a face-on star (left lower panel). At higher inclinations, the lines of constant intensity are no longer straight because of the strong influence of beaming. Along the vertical line the angle σ′ at which we see the surface element is different, resulting in different local intensity. Thus, the observed intensity peak is now not at the edge (where σ′=π/2), but is reached at the point where the angle of the line of sight to the local normal is smaller. We note here that our computation method is very efficient: computing one image of resolution 70 × 400 in latitude and azimuth on a standard laptop with IDL takes less than 0.3 s, which is about 3000 times faster than would be needed if an exact ray-tracing code is used (e.g. Vincent et al. 2018).

Fig. 5. Images of a NS rotating at 700 Hz. The mass and the radius of a non-rotating NS of the same baryonic mass are M = 1.5 M⊙ and R = 14 km. The white curves and the colours show the contours of constant bolometric intensity. Upper panels: blackbody case of local isotropic intensity. Lower panels: electron-scattering dominated atmosphere. Panels from the left to the right correspond to the inclinations i = 0°, 45°, and 90°. The lines of constant latitudes (every 10°) and longitudes (every 15°) are shown in black. The colour scheme is different for the two rows.

3. Modified direct cooling tail method

3.1. Non-rotating X-ray bursting neutron stars

The direct cooling method (Suleimanov et al. 2017a) assumes that the investigated NS with given mass M and radius R rotates slowly, its shape is spherical, and the surface has a uniform effective temperature T_eff. In this case a single model atmosphere with the given T_eff, the surface gravity g = GM(1 + z)/R² and the chemical composition is enough to describe the spectrum of the whole NS. The gravitational redshift factor is related to the NS Schwarzschild radius R_S = 2GM/c² as 1 + z = (1 − R_S/R)^−1/2. Therefore, we can connect the local atmosphere parameter T_eff with the whole observed NS luminosity L = 4πR² σ_SB , and use the relative NS luminosity ℓ = L/L_Edd instead of T_eff. Here

(11)

is the NS Eddington luminosity, κ_e = 0.2(1 + X) cm² g⁻¹ is the Thomson opacity, and X is the hydrogen mass fraction in the atmosphere.

It is well known (London et al. 1986; Lapidus et al. 1986) that the local spectra of hot NS model atmospheres are well described with a diluted blackbody (8), where the dilution and the colour-correction factors approximately follow the relation and f_c >  1. This explains why the observed spectra of X-ray bursting NSs are usually very well fitted with a blackbody given by Eq. (10) (see e.g. Galloway et al. 2008), where T_BB = f_cT_eff/(1 + z) is the observed colour temperature. As a result, the blackbody normalization parameter K used in the fits depends on the dilution factor w (see details in Suleimanov et al. 2011a)

(12)

where D is the distance to the source and

(13)

is a geometrical dilution factor proportional to the solid angle occupied by the NS on the sky.

The parameters w and f_c are found by fitting the emergent model spectra of hot NS atmospheres in the energy band of some X-ray instruments, like RXTE/PCA, blueshifted to the NS surface. The values of these parameters depend mainly on the relative luminosity ℓ, although they depend on the surface gravity g and the chemical composition as well. We have computed an extended grid of models (Suleimanov et al. 2011b, 2012; Nättilä et al. 2015) for a range of log g from 13.7 to 14.9 with the step 0.15 and for various atmosphere chemical compositions (pure hydrogen, pure helium, solar abundance, solar H/He mix with reduced and increased heavy element abundances, and pure iron). The model spectra were fitted with a blackbody resulting in dependences w − ℓ and f_c − ℓ.

The direct cooling tail method is based on the results of these computations. It is possible to demonstrate (Suleimanov et al. 2017a) that the observed dependence K − F_BB during the thermonuclear burst cooling tail has to be fitted with the model dependence computed for each pair of parameters (M, R) and the actual chemical composition of the atmosphere². Formally, there are two fitting parameters Ω and the observed Eddington flux

(14)

but actually for every pair (M, R) they both depend on the distance to the source D alone because the z value is known. The fitting procedure provides a χ² map, which can be used to estimate the most probable values for the NS mass and radius (see examples in Suleimanov et al. 2017b,a).

3.2. Spectra of bursting rapidly rotating NSs

Once we understand the observational appearance of the blackbody emitting NS as discussed in Sect. 2, we can discuss the emission of a bursting rotating NS. One of the effects that needs to be accounted for is the change in the effective gravity along the latitude. The effective gravity can be presented approximately as , where R(θ) is the NS radius at co-latitude θ. The approximate expressions of the effective surface gravity using relativistic computations were suggested by AlGendy & Morsink (2014) (see also Appendix A). It is clear that the effective surface gravity at the poles is stronger than at the equator. It means that the value of the local Eddington flux depends on latitude. This also implies that it is not possible to unambiguously determine the touch-down point in PRE bursts because it is not known a priori how the actual radiation flux is distributed over the NS surface. We can only make some suggestions. We consider two cases. The first option is that the bolometric flux in units of the local Eddington flux F_Edd = c g_eff/κ_e is constant over the NS surface (i.e. F(θ)/F_Edd(θ) = ℓ). In this constant relative flux (CRF) case the photosphere of the extended envelope touches the NS surface simultaneously and we observe this moment as a touch-down point in the X-ray burst light curve. The second case assumes constant absolute flux (CAF), F(θ) = const, at every co-latitude equal to the flux of a non-rotating NS with a given relative luminosity ℓ. In this case, the photosphere touches the NS surface at the poles first and we can associate the touch-down point with this moment (see Fig. 6). The real latitudinal flux distribution is not known, and we can only guess that it lies between the two limiting cases considered above. In this work we consider both possibilities as equiprobable and estimate the effect of rapid rotation on the observed spectra and derived NS parameters for both cases. We note here that constant surface flux does not mean same spectrum of the emitted radiation, because the real atmosphere spectrum depends not only on the effective temperature, but also on local gravity (i.e. the colour-correction factor varies along the latitude). Thus, even this problem cannot be reduced directly to the case of a blackbody emitting NS considered in Sect. 2.4.

Fig. 6. Qualitative difference between the shapes of the extended photospheres (dashed contours) at the touch-down point for two different hypothesis on the latitudinal distribution of the radiation flux over the surface of a rapidly rotating NS, CRF case (top panel) and CAF case (bottom panel).

Now let us discuss the appearance of the bursting rotating NS to a distant observer at inclination i for a given relative luminosity ℓ. The details of the method to compute the observed spectrum are presented in Appendix B. The local specific intensity is either taken from the atmosphere models computed by Suleimanov et al. (2012); Suleimanov et al. (2017a) for the given T_eff(θ) and g_eff(θ) or approximated by a diluted blackbody given by Eq. (9); the dilution factor w and the colour-correction factor f_c are found by interpolation in the existing grid of the precomputed values (see Suleimanov et al. 2017a). The accuracy of this approximate method is demonstrated in Fig. 7 where the results of the accurate computations taking real atmosphere models (solid curves) are compared with the approximate ones (dashed curves) for two relative luminosities (ℓ= 0.9 and 0.1), two inclination angles, i = 0° and 90°, and assuming the CRF case. The errors in the bolometric fluxes are typically below 0.5% and 1.5% for the high- and low-luminosity models, respectively. They are not significant for the following discussion.

Fig. 7. Spectra of the rotating NS with the parameters M′ = 1.519 M⊙, Re = 15.48 km, and ν* = 700 Hz computed for two inclination angles i= 0° (blue solid curve) and i= 90° (red solid curves), and for two relative luminosities (same over the surface) ℓ= 0.9 (top panel) and ℓ= 0.1 (bottom panel). The solid curves represent spectra computed taking accurate local spectra for atmosphere models from Suleimanov et al. (2012). The dashed curves correspond to the diluted blackbody approximation of the local spectra. Solar chemical composition of the atmosphere is assumed.

The final stage is the approximation of the observed spectrum from a rapidly rotating NS by a diluted blackbody w′πB_E in the observed energy range 3–20 keV, where T_eff, ∞ is the redshifted effective temperature of the corresponding non-rotating NS with the relative luminosity , which is an independent input parameter of the problem. The observed bolometric luminosity of the rotating NS model L(ν_*, i) can also be computed, together with the relative luminosity ℓ′  =  L(ν_*, i)/L_Edd(0), where L_Edd(0) is the observed Eddington luminosity of the corresponding non-rotating NS of mass M. Repeating this procedure at a grid of we obtain the model dependences w′−ℓ′ and , which are shown in Figs. 8 and 9. The dependence w′−ℓ′ has to be compared with an observed dependence K − F_BB for some appropriate X-ray bursting NS. We chose the mass and the radius of a non-rotating NS to provide the surface gravity log g= 14.15. Then we took the rotation frequency to be equal 700 Hz, and computed the curves w′−ℓ′ and for three inclination angles, 0°, 45°, and 90° assuming solar chemical compositions with heavy element abundances reduced 100 times. The curves were computed for both limiting cases of the flux distribution over the surface, CRF (Fig. 8) and CAF (Fig. 9). In the second case, the maximum luminosities that can be reached correspond to the local Eddington limit at the NS equator.

Fig. 8. Dependences of and w′ on the relative luminosity for a rapidly rotating NS for three inclination angles. The corresponding dependences for a non-rotating NS with the surface gravity log g0= 14.15 are shown with the dashed curves. Here chemical composition is solar with the heavy element abundances reduced by 100 times. The constant relative flux distribution over the NS surface is assumed (CRF case).

Fig. 9. Same as Fig. 8, but computed with constant absolute flux distribution over the NS surface (CAF case).

The qualitative behaviour of the colour-correction factor and the dilution factor w′ is similar to those shown in Fig. 4. The colour correction shows a greater increase for highly inclined NSs (i = 90°), whereas the greater increase in the dilution factor occurs for the face-on NS (i = 0°). The results for the CAF case and the relatively low luminosity (L/L_Edd ∼ 0.1) are very close quantitatively to the results shown in Fig. 4 (R = 14 km, the electron scattering case). In this case, the local diluted blackbodies at different latitudes are very close to each other. At the higher relative luminosities the local blackbodies differ more as the relative flux F(θ)/F_Edd(θ) near the equator approaches the local Eddington limit.

3.3. Influence of NS rotation on determination of basic parameters

Using computed dependences w′−ℓ′ we can estimate the effect of rapid rotation on determination of the NS radius. Let us assume that the computed dependence w − ℓ for a non-rotating NS (see the dashed curves in the bottom panels of Figs. 8 and 9) closely fits the observed dependence K − F_BB for some X-ray burst occurring at the surface of a rapidly rotating NS with ν_* = 700 Hz. Let us also suggest that the fitting parameters (see Sect. 3.1) are F_Edd, ∞= 6.76×10⁻⁸ erg s⁻¹ cm⁻² and Ω= 1261 (km/10 kpc)². They can be combined to obtain a distance-independent quantity, the observed Eddington temperature, which corresponds to the dashed curve at the M − R plane in the bottom panels of Figs. 10 and 11.

Fig. 10. Top panel: results of fitting the w − ℓ dependence with w′−ℓ′ dependences computed for different inclination angles i and assuming constant relative flux over the NS surface (CRF case, see Fig. 8). The corresponding correction factors are (AF, AR) = (1.053, 0.856) for i = 0°, (1.1, 0.972) for 45°, and (1.15, 1.12) for 90°. Bottom panel: corresponding solutions at the (M, R) plane obtained assuming a non-rotating NS (dashed curve) and using additional correction factors AF and AR for different inclination angles.

Fig. 11. Same as Fig. 10, but for constant absolute flux over the NS surface (CAF case, see Fig. 9). The corresponding correction factors are (AF, AR) = (1.12, 0.85) for i = 0°, (1.16, 0.97) for 45°, and (1.2, 1.12) for 90°.

It is possible to fit the dependence w − ℓ with the dependences w′−ℓ′ computed for the rapidly rotating NS (see upper panels in Figs. 10 and 11) and to obtain the correction factors A_F and A_R to the fitting parameters F_Edd, ∞ and Ω. This means that the new parameters, which should actually be used to determine M and R of the non-rotating NS, are and Ω′=A_RΩ. The new parameter is larger than the old one (A_F >  1) for all inclination angles because the critical Eddington limit is reached at a rotating NS at a lower luminosity. The new observed solid angle Ω′ maybe smaller or larger than Ω depending on the inclination: A_R >  1 for large inclination angles and A_R <  1 for a face-on rotating NS (see Figs. 8 and 9, bottom panels). The value of A_R depends on two factors. It depends first on the apparent area, which is always larger for a rotating NS (see Fig. 3), which would make A_R smaller than unity. In addition, the dilution factor depends on the colour temperature of the spectrum of a rotating NS. The higher the ratio of the colour to the local effective temperatures, the smaller the dilution factor w′. For highly inclined rapidly rotating NSs, the colour temperature is higher due to the Doppler boosting (see Fig. 4), and its influence is more important than the increase in the apparent area resulting in A_R >  1.

As a result the solutions at the M − R plane for the different inclination angles are shifted (see Figs. 10 and 11, bottom panels). The interpretation of these shifts is as follows. If we use the non-rotating NS model curve w − ℓ, we obtain an incorrect solution presented by the dashed curve instead of the correct one shown by the solid curves for different inclination angles. For example, if we treat a face-on rapidly rotating (with ν_* = 700 Hz) NS as a non-rotating NS, we obtain that the NS radius R is 3–3.5 km larger than that obtained using the method that accounts for rotation. For an edge-on system, the radius is nearly the same within 0.5–1 km. We note here once more that the equatorial radius of the rapidly rotating NS is larger than the radius of the non-rotating NS with the same baryonic mass. For instance, a non-rotating NS with R = 12 km and M = 1.5 M_⊙, will have the equatorial radius R_e ≈ 13 km if it rotates at ν_* = 700 Hz. Because there is an observational bias which restricts the inclination for real observed NSs in LMXBs to be i <  70°, treating NSs of unknown inclination as non-rotating will always result in an overestimate of the radius. We finally note that the two models for the flux distribution over the NS surfaces give similar radius values, within 0.5 km (compare Figs. 10 and 11), but the CAF model systematically gives a smaller radius.

3.4. Essentials of the method

Now we describe the method for obtaining NS parameters from the data, using the models for rotating NSs. We need to have a reasonable estimate of the rotation frequency ν_*, the inclination, and the chemical composition of the atmosphere. Let us choose a pair of mass and radius (M, R) of the corresponding non-rotating NS, and compute the model curve w′−ℓ′ using the method described above. Then we fit the observed K − F_BB dependence with this model curve using the same fitting parameters F_Edd, ∞ and Ω as for a non-rotating star. The corrections for the rapid rotation are accounted for by differences between the model curve w′−ℓ′ and the corresponding model curve w − ℓ for a non-rotating NS. Both fitting parameters depend on the distance to the source D, which is actually the only fitting parameter for the given pair (M, R).

We repeat the fitting procedure at a grid (M, R) obtaining the χ² map which determines the most probable values of M and R of a non-rotating NS and corresponding confidence intervals. This map thus depends on the assumed values of the inclination angle i, the rotational frequency ν_*, and the chemical composition of the atmosphere.

4. Application of modified cooling tail method to X-ray bursts from rapidly rotating NSs

SAX J1810.8−2609 was used as a key source for NS radius determination in two works (Nättilä et al. 2016; Suleimanov et al. 2017a), where we found the radius to be in the range 11.0–12.5 km assuming M = 1.5 M_⊙ and that the NS rotates slowly. Unfortunately, there is only one X-ray burst that could be used for the cooling tail method. Recently, Bilous et al. (2018) using archival RXTE observations discovered in that source burst oscillation at frequency ν_*= 531.8 Hz. Such a high rotation rate may affect the radius determination significantly, and we can use it as a test bed for our new modified cooling tail method.

The burst that we use for our study has Obs ID 93044-02-04-00, and it was observed on 2007 August 13. The details about the burst, the data reduction, and blackbody fitting can be found in Nättilä et al. (2016). The results of the application of the modified cooling tail method are shown in Fig. 12. We consider three different cases. First of all, in order to test the new method and compare it with the previous results, we assumed that the NS is non-rotating, ν_*= 0 Hz, but used the new method modified for rapid rotation. The resulting χ² map, presented in Fig. 12 with solid black curves, is identical to the map produced with the direct cooling tail method for a non-rotating NS shown in Fig. 5 in Suleimanov et al. (2017a).

Fig. 12. Confidence regions at the M–R plane for SAX J1810.8−2609 obtained from the χ2 map using the cooling tail method modified for rapid rotation. The solid, dashed, and dotted curves indicate the 68%, 90%, and 99% regions, respectively. The red contours correspond to the inclination i = 60°, while the blue contours are for the face-on NS. The black contours give the constraints assuming a non-rotating NS. (a) CRF case, i.e. constant relative flux. (b) CAF case, i.e. constant absolute flux. The triple dot-dashed curves are the contours of the constant distance of 4.0 and 4.5 kpc (distance is the only free parameter for a given pair of M and R).

Then we computed χ² maps taking the actually observed rotational frequency and assuming two inclination angles, i= 0° and 60°. We considered two limiting cases for the flux distribution over the NS surface as described in Sect. 2.4: the constant relative flux (CRF, Fig. 12, left panel) and the constant absolute flux (CAF, Fig. 12, right panel). In both cases the assumption about the face-on inclination gives the NS radius that is approximately 1 km smaller than that assuming a non-rotating NS. The best estimate of the NS radius is 11.2 ± 0.5 km at M = 1.5 M_⊙, and the corresponding equatorial radius R_e is 11.6 ± 0.6 km. Taking a higher inclination of 60° leads to the χ² maps closer to the map assuming no rotation and the NS radius of 11.8 ± 0.5 km at M = 1.5 M_⊙. In this case the corresponding equatorial radius, 12.3 ± 0.6 km, is close to the radius of the NS determined under assumption that it is non-rotating. We note that the CRF case gives smaller χ² values than the CAF model, and also a smaller χ² is achieved at higher inclination (i = 60°). However, the statistical significance of the improvement is below 3σ.

We recall that the only free parameter in the direct cooling tail method (in addition to the mass and radius) is the distance, which therefore can be estimated. The solutions that correspond to the distances of 4.0 and 4.5 kpc are shown by the triple dot-dashed lines in Fig. 12. There is only one estimation for the distance to SAX J1810.8−2609 (Natalucci et al. 2000), but it cannot be considered as an independent one, as they used the maximum fluxes of type I X-ray bursts assuming that they correspond to the Eddington flux for a NS with M = 1.4 M_⊙. Nevertheless, their estimation of 4.9 ± 0.3 kpc is close to ours.

Another X-ray bursting source where rapid rotation should be accounted for is 4U 1608−52, which rotates at ν_*= 620 Hz (Muno et al. 2002). The cooling tail method applied to this source assuming no rotation gave the NS radius in the range 13–16 km depending on the selected data points for M = 1.5 M_⊙ (Poutanen et al. 2014). We found that the rapid rotation with ν_*= 700 Hz can increase the visible NS radius by 2–3 km for face-on systems and by 1 km for edge-on systems. Thus, accounting for the rotation, the NS radius in 4U 1608−52 will be in range 10–13 km for low inclination and 12–15 km for an inclination around 70°. This is consistent with the current best limits on the NS radii 10.5–13.5 km coming from other sources and observables (Steiner et al. 2013; Nättilä et al. 2016, 2017; Suleimanov et al. 2017b; Abbott et al. 2018).

5. Summary

We have presented a new version of the cooling tail method, modified for rapid rotation, for determining NS parameters from X-ray bursts. The aim of the method is to find the most probable values for the mass and radius of a non-rotating NS, whose rotating model describes the observed spectral evolution of X-ray bursts taking place on the surface of a rapidly rotating NS. The rotation frequency ν_*, the inclination angle of the NS rotation axis to the line of sight i, and the atmospheric chemical composition are the input parameters of the model.

First, we developed a model that transforms a non-rotating NS with given M and R to the approximate model of a rapidly rotating NS with the same baryonic mass of the increased equatorial radius R_e and the slightly increased gravitational mass M′. We derived the approximate relations between R and R_e, and between M and M′ using the models of rapidly rotating NSs computed by Cook et al. (1994) for various EoS.

At the next step we computed the emergent spectra of the obtained rotating NS models for different relative luminosities ℓ of the non-rotating NS with the given M and R. As the effective temperature distribution over the rotating NS surface at any burst moment is not known a priori, we considered two limiting cases: CRF (i.e. the same ratio ℓ of the bolometric flux to the local Eddington flux) and CAF (i.e. the same absolute bolometric flux at each latitude, the same as for the non-rotating NS with a given ℓ). We note that the local Eddington flux depends on the effective gravity which includes the centrifugal force. We showed that approximating NS atmosphere spectra by a diluted blackbody gives a good description of the observed spectra from a rotating NS. Therefore, for further modelling we used a diluted blackbody approximation for the local spectra. The parameters of these spectra, the dilution factor w and the colour-correction factor f_c, are found by interpolating in the precomputed set of models at a grid of the effective temperatures, effective surface gravities, and chemical compositions of the atmosphere from Suleimanov et al. (2012); Suleimanov et al. (2017a).

We used the universal (independent of EoS) slow-rotation approximations for description of the shape of the rotating NS and the effective surface gravity distribution suggested by Morsink et al. (2007) and AlGendy & Morsink (2014). To compute the observed spectra and integrate over the visible surface of a rapidly rotating NS, we used a simplified approach computing approximately the total redshift factor using a metric of the rotating NS in a slow-rotation approximation and making corrections for the oblate shape (see Morsink et al. 2007; Nättilä & Pihajoki 2018; Salmi et al. 2018). The light bending angle and the lensing factor were computed using accurate analytical formulae from Poutanen (2019), which allowed us to speed up calculations by orders of magnitude.

The computed spectra are then fitted with the diluted blackbody model. The fitting parameters, the dilution factor w′, and the colour-correction factor , are found from comparison with the fiducial spectrum, which was the blackbody spectrum with the effective temperature of a non-rotating NS with the same relative luminosity ℓ. We thus obtained the model dependence w′−L(ν_*, i)/L_Edd(0), where L_Edd(0) is the Eddington luminosity of a non-rotating NS. This model curve can then be compared to the observed dependence K − F_BB using the distance to the source as the only fitting parameter, in addition to M and R of a non-rotating NS. As a result we get a χ²-map at the M − R plane, which allows us to determine the most probable values and the confidence region for M and R of the corresponding non-rotating NSs. These constraints depend on the assumed inclination i, rotation rate ν_*, and chemical composition.

We applied our method to an X-ray burst observed from SAX J1810.8−2609 and investigated earlier by Nättilä et al. (2016) and Suleimanov et al. (2017a), who assumed no NS rotation. We revised the NS radius determination using our modified cooling tail method accounting for rapid rotation of this NS with ν_*= 531.8 Hz (Bilous et al. 2018). We found that the NS radius is smaller by ≈1 km (11.2 ± 0.5 km instead of 12.2 ± 0.5 km assuming no rotation) for the face-on inclination and M = 1.5 M_⊙. The corresponding equatorial radius R_e of this rapidly rotating NS, 11.6 ± 0.6 km, is still slightly smaller than the radius of the NS obtained assuming no rotation. At high inclination the effect of rapid rotation is small.

Acknowledgments

This research has been supported by the grant 14.W03.31.0021 of the Ministry of Science and Higher Education of the Russian Federation. V. F. S. also thanks Deutsche Forschungsgemeinschaft (DFG) for financial support (grant WE 1312/51-1). We thank the German Academic Exchange Service (DAAD, project 57525212), the Academy of Finland (projects 317552 and 331951) and the Magnus Ehrnrooth foundation for travel grants.

References

Appendix A: Shape and gravity of a rotating neutron star

Let us consider a rotating NS with the gravitational mass M′, equatorial radius R_e, and observed angular velocity Ω_* = 2πν_*. For completeness, we present here the whole collection of formulae from AlGendy & Morsink (2014) that we use to define the NS shape and gravity. They demonstrate that basic properties of the rotating NS model can be represented using two dimensionless parameters

(A.1)

The metric of a stationary axisymmetric rotating NS is usually written in the form (Butterworth & Ipser 1976)

(A.2)

where metric coefficients depend on coordinates and co-latitude θ. The circumferential radius is related to by Friedman et al. (1986)

(A.3)

The metric coefficients for a non-rotating stationary NS are known as the isotropic Schwarzschild metric coefficients

(A.4)

where

(A.5)

A well-known relation exists between the circumferential r and isotropic Schwarzschild radial coordinate :

(A.6)

The metric coefficients in Eq. (A.2) can be expanded in powers of up to order 2 (AlGendy & Morsink 2014)

(A.7)

where P₂ is the Legendre polynomial of order 2, and b and q are the dimensionless coefficients

(A.8)

(A.9)

For a given NS described by parameters x and , for given circumferential radius r and angle θ, the metric coefficients are computed using Eqs. (A.3), (A.4), (A.7), and (A.5) via iterations. The angular velocity of the local zero angular momentum observer with respect to an observer at rest at infinity ϖ can also be expanded with the same accuracy

(A.10)

where 𝒥= IΩ∗ is the NS angular momentum, is the moment of inertia and

(A.11)

Eq. (A.10) can be written in dimensionless form as

(A.12)

Following AlGendy & Morsink (2014) we use the approximation for the circumferential radius of a rapidly rotating NS

(A.13)

and for the effective gravity in a slow-rotation approximation

(A.14)

where

(A.15)

is the surface gravity for a spherical NS with mass M′ and radius R_e, and

(A.16)

We note that the approximate fitting formulae for the effective surface gravity (A.14) and the rapidly rotating NS shape (A.13) are not self-consistent as they are independent approximations to the accurately computed values. Slightly different approximations for the gravity and for the shape are given by Eq. (50) in AlGendy & Morsink (2014) and Eqs. (8)–(10) in Morsink et al. (2007), respectively. A comparison of these approximations to those given by Eqs. (A.14) and (A.13) for two rotating NSs is shown in Fig. A.1. The approximation to the NS shape derived by Morsink et al. (2007) (dashed curves) gives a less oblate NS shape. The difference in the effective gravity approximations is insignificant for slowly rotating NSs. We also compared the shape with the models published by Cook et al. (1994). For example, for M′ = 1.4103 M_⊙, R_e = 11.71 km, and Ω = 5033.6 s⁻¹ (ν_* = 801.12 Hz), Cook et al. (1994, see their Table 15) give the eccentricity ϵ = 0.488, while Eq. (A.13) gives ϵ = 0.495, which is about 2% accurate. We conclude that we can use approximations (A.13) and (A.14) for slowly rotating NSs with .

Fig. A.1. Top panel: shape of the rotating NS computed with two different approximations. The solid curves correspond to Eq. (A.13), while the dashed curves to Eqs. (8)–(10) in Morsink et al. (2007). Bottom panel: effective surface gravity distributions. The solid curves correspond to Eq. (A.14), while the dashed curves to Eq. (50) in AlGendy & Morsink (2014). The red curves correspond to a relatively slow rotating NS with (ν* = 842 Hz, M′ = 1.8 M⊙, Re = 13.6 km). The blue curves are for a faster rotating NS with (ν* = 1084 Hz, M′ = 1.99 M⊙, Re = 15.1 km).

Appendix B: Radiation observed from a rotating neutron star

Accurate computations of the observed radiation from a rapidly rotating NS are rather complicated. The main reason is a complex form of geodesics in the vicinity of the star, which requires using the full ray tracing approach (see e.g. Bauböck et al. 2012; Nättilä & Pihajoki 2018). This is a time-consuming approach that is not appropriate for the extensive computations presented in this work. A simplified approach was proposed by Morsink et al. (2007). They considered a local Schwarzschild metric at every surface point of an oblate NS. They computed gravitational redshift and the light bending effect as for a non-rotating NS adding the Doppler boosting factor that accounts for rotation (oblate Schwarzschild or OS approach). This approach gives acceptable results for slowly rotating NSs with ν_* <  700 Hz. However, we believe that it is possible to account for the main effects introduced by the metric of the rapidly rotating NS without significant complication of the OS approach cited above. Nättilä & Pihajoki (2018) showed that it is important to have an accurate treatment of the redshift effect (i.e. the ratio of the observed to the emitted photon energy E/E′), which is affected both by the non-zero quadrupole momentum and by the frame dragging impact on the Doppler boosting factor. Thus, we introduce here a new modified OS (MOS) approximation, using the metric of a rapidly rotating NS instead of the local Schwarzschild metric for the photon energy change, but keeping simplified Schwarzschild description for the light bending, which is the main simplification comparing to the exact treatment.

Let us start with discussing the emission of a rapidly rotating spherical NS in Schwarzschild metric. Then we consider the correction coming from the shape of the NS. Finally, the correction to the redshift due to the non-zero quadrupole momentum and by the frame dragging is discussed.

B.1. Rapidly rotating spherical star

The flux observed by an observer can be computed in two different ways. We can integrate over the spherical surface in the non-rotating frame, or make integration over the surface in the co-rotating frame of the NS. Both approaches should produce the same result.

Let us start from the first approach. We introduce a spherical coordinate system with the polar axis along the NS rotation axis. An observer is situated at distance D in the direction given by the unit vector , with i being the inclination (see Fig. B.1). A surface element in the static frame can be defined in terms of co-latitude θ and azimuthal angle ϕ with the unit vector of its position being . Its area is dS = R²d cos θ dϕ. The solid angle of this element as observed at distance D is

(B.1)

Fig. B.1. Geometry of the problem.

where the impact parameter and the polar angle φ are defined using polar coordinates on the image plane centred to the observed NS centre. The impact parameter is related to the emission angle α, the angle between the radial direction , and the unit vector along photon momentum close to the NS surface as (see e.g. Beloborodov 2002)

(B.2)

where u = R_S/R is the NS compactness with R_S = 2GM/c² being the Schwarzschild radius. Substituting Eq. (B.2) to Eq. (B.1) and given that the element area can also be written as dS = R²d cos ψ dφ, with ψ being the angle between the local radius-vector at the NS surface and the line-of-sight in flat space-time,

(B.3)

we get (Beloborodov 2002; Poutanen & Beloborodov 2006)

(B.4)

with the lensing factor

(B.5)

The relation between α and ψ is given by the integral (Misner et al. 1973; Pechenick et al. 1983; Beloborodov 2002)

(B.6)

It can be computed using an accurate numerical procedure (e.g. as is given in the Appendix of Salmi et al. 2018), but here we use for simplicity an approximate formula suggested by Poutanen (2019)

(B.7)

where y = 1 − cos ψ, and

(B.8)

This formula is an improvement over the analytical approximation suggested by Beloborodov (2002) and is 0.05% accurate over the whole range of angles of interest and for any realistic NS compactness. Using this relation we obtain the expression for the lensing factor

(B.9)

The unit vector along photon momentum close to NS surface can be found from

(B.10)

The observed flux at energy E from the element will be

(B.11)

The observed specific intensity can be related to the intensity measured in the co-rotating frame as

(B.12)

where the ratio of energies combines gravitational redshift and the Doppler effect

(B.13)

and in principle we can consider the possibility that the local intensity also depends on co-latitude θ and the zenith angle α′ measured in the corotating frame of the element. The Doppler factor

(B.14)

depends on the NS velocity at this latitude relative to the external non-rotating frame,

(B.15)

and the Lorentz factor is

(B.16)

The angle ξ between the photon momentum and the spot velocity vector β = β(− sin ϕ, cos ϕ,0) in the external static frame can be expressed as (Poutanen & Gierliński 2003; Poutanen & Beloborodov 2006)

(B.17)

The angle α′ in Eq. (B.12) the photon momentum makes to the radial direction as measured in the co-rotating frame of the spot is related to the similar angle measured in the static frame as

(B.18)

We finally obtain the total observed flux by integrating Eq. (B.11) over the whole visible stellar surface

(B.19)

Now let us consider an alternative derivation when we integrate over the NS surface in the co-rotating frame. Let us define an element on the NS surface at co-latitude θ and azimuthal angle ϕ′ in the frame corotating with the NS, with the extent given by dθ and dϕ′. The area of this element is dS′=γR²d cos θ dϕ′, where the Lorentz factor γ appears because the area expressed in angular coordinates is measured by co-moving observers (Nättilä & Pihajoki 2018; Lo et al. 2018; Bogdanov et al. 2019). The solid angle dS′ occupies on the observer’s sky is given by Eq. (B.4) with the only difference that cos α dS should be substituted by cos α′ dS′, so that we obtain

(B.20)

The observed flux is thus (Poutanen & Beloborodov 2006)

(B.21)

We note here that it is important to distinguish between the emission phase ϕ (i.e. position of the element when photons were emitted as measured by external static observer) and the arrival (observed) phase ϕ_obs when the flux from this element is actually observed. Their differentials are simply related by the standard time contraction formula as (Rybicki & Lightman 1979)

(B.22)

Now averaging the flux given by Eq. (B.21) over the observed phase and changing the integration variable to the emitted phase using Eq. (B.22) we obtain

(B.23)

We now notice that the factor γ(1 − β cos ξ) has cancelled out with one of δ. Further integration over ϕ′ reduces to just removing the factor 1/2π because the flux does not depend on the choice of ϕ′. Thus, for the total flux observed from a star we arrive at Eq. (B.19).

B.2. Rapidly rotating oblate star

For the oblate star with the radius depending on co-latitude as R(θ), we need to make small modifications to the formulae presented in the previous section (Morsink et al. 2007; Salmi et al. 2018). Now the local normal does not coincide with the radial direction, and the angle between them η (see Fig. B.1) can be expressed as

(B.24)

where

(B.25)

(B.26)

and the derivative of the circumferential radius could be found from Eq. (A.13):

(B.27)

The area of the surface element (as measured in the static frame) is

(B.28)

where the factor cos η accounts for deviation of the stellar surface from the spherical surface.

The local atmosphere of a rotating oblate NS is axisymmetric relative to the local normal. Therefore, we have to transform the angle α measured relative a local radius-vector to the angle σ measured relative to a local normal. This angle can be found from the spherical law of cosines

(B.29)

where

(B.30)

and thus

(B.31)

Similarly to the spherical case correcting for relativistic aberration, we obtain the emission zenith angle in the frame comoving with the surface

(B.32)

Thus, we arrive at a modified expression for the solid angle

(B.33)

If we know the local specific emergent radiation intensity I_E′(σ′,θ) at the surface of a rotating NS as a function of energy E′, co-latitude θ, and the angle σ′ between the emergent ray and the surface local normal, we can compute the total observed spectrum F_E from the whole surface

(B.34)

We note here in the formulae for the light bending angle (B.7) and the lensing factor (B.9), both u and R are now latitude-dependent. In the calculations presented in this paper we considered different cases for the local specific intensity as described in Sect. 2.3. The total photon energy change now can be expressed as

(B.35)

The three terms in this expression correspond to the Doppler factor δ, a pure gravitational redshift e^ν (which accounts for quadrupole moment, see Eq. (A.7)), and a phenomenological factor with

(B.36)

being the frame dragging dimensionless velocity. The expression for ϖ is given by Eq. (A.10). In the formulae for the Lorentz (B.16) and Doppler factors (B.14), instead of surface velocity given by Eq. (B.15) we use the velocity relative to a zero angular momentum observer:

(B.37)

These ad hoc formulae nevertheless give accurate values of the total energy redshift 1 + z ≡ E′/E at the pole and for the maximum and minimum values at the equator for an observer at i = 90°. For example, for a NS model published in Table 15 of Cook et al. (1994) with M′ = 1.4103 M_⊙, R_e = 11.71 km, and Ω = 5033.6 s⁻¹ (i.e. ν_* = 801.12 Hz), our redshift at the pole z_p = 0.2796 is lower by just 0.1% than the value z_p = 0.280 given there. Our maximum and minimum redshifts at the equator (for i = 90°), z_b = 0.6094 and z_f = −0.0244, are very close to those computed by Cook et al. (1994), z_b = 0.606 and z_f = −0.025.

We also computed the observed profiles of a narrow emission line emitted by a rapidly rotating NS. In particular, we reproduced Fig. 5 from Nättilä & Pihajoki (2018), where M′ = 1.4 M_⊙, R_e = 10 km, ν_* = 700 Hz, and i = 20° were assumed. They computed the line profile using four different set-ups: the exact solution, the exact solution with the quadrupole moment artificially increased by four times, an oblate NS in local Schwarzschild metric, and a spherical rotating NS (R = R_e) in Schwarzschild metric with Doppler boosting taken into account. As intrinsic surface emission we took a Gaussian line at unit energy with the standard deviation σ = 0.002. Our results presented in Fig. B.2 demonstrate that our approach allows us to reproduce the main features of the correct solution. The line profiles computed at different inclination angles are shown in Fig. B.3 for M′ = 1.5 M_⊙, R_e = 14 km, and ν_* = 600 Hz. The results closely reproduce the main features of the line profiles shown in Fig. 6 of Nättilä & Pihajoki (2018).

Fig. B.2. Line profiles from a NS with M′ = 1.4 M⊙, Re = 10 km rotating at 700 Hz, and seen by a distant observer at inclination i = 20°. The line profiles for a spherical NS with R = Re and for an oblate NS shape both with Schwarzschild exterior metric are shown by the black and blue curves, respectively. Accounting for frame dragging and quadrupole moment gives the profile shown by the red curve. The line obtained with a quadrupole moment artificially increased by a factor of 4 above that given by Eq. (A.9) is shown by the dashed red curve. Compare to Fig. 5 in Nättilä & Pihajoki (2018).

Fig. B.3. Line profiles observed at different inclination angles, i = 90° (black), 60° (orange), 40° (violet), 20° (magenta), 10° (blue), and 5° (red). The adopted NS parameters are M′ = 1.5 M⊙, Re = 14 km, ν* = 600 Hz. The energy range is shifted by the redshift factor . Compare to Fig. 6 in Nättilä & Pihajoki (2018).

All Figures

Fig. 1. Dependences of the ratio Re/R on the relative rotational frequency for numerical NS models with various EoS (Cook et al. 1994) for the maximum possible mass (top panel) and the mass M = 1.4 M⊙ (bottom panel). Red solid curves show the fits with Eq. (2) for the masses M = 2.5 M⊙ (top panel) and 1.4 M⊙ (bottom panel). The fitting curve for M = 1.4 M⊙ is shown in the top panel as dashes. The relative errors of the fits are also shown in the narrow additional panels below the main panels.

In the text

Fig. 2. Dependences of ratios M′/M on the relative rotational frequency for numerical NS models with various EoS (Cook et al. 1994) for the maximum possible mass (top panel) and the mass M = 1.4 M⊙ (bottom panel). Red solid curves show the fits with Eq. (3) for the masses M = 2.5 M⊙ (top panel) and 1.4 M⊙ (bottom panel). The fitting curve for M = 1.4 M⊙ is shown in the top panel as dashes. The relative errors of the fits are also shown in the narrow additional panels below the main panels.

In the text

Fig. 3. Top panel: ratio of the apparent area of a rapidly rotating NS to that of the non-rotating star (of the same baryonic mass) as a function of the rotational frequency for inclination angles i= 0°, 45°, and 90° is shown by black, red, and blue curves, respectively. The mass of a non-rotating star was assumed M = 1.5 M⊙ and two radii R = 14 (solid curves) and 12 km (dashed curves) are considered. Bottom panel: apparent area of a rapidly rotating NS for the same family of models as in the top panels divided by the apparent area of a non-rotating NS of the same gravitational mass M′ (from Eq. (3)) and the radius equal to equatorial radius Re given by Eq. (2). The red and the blue dashed curves overlap in both panels.

In the text

Fig. 4. Top panels: colour correction to the observed blackbody temperature compared to a non-rotating NS for a NS of constant local temperature as a function of the rotational frequency for inclination angles i= 0°, 45°, and 90° is shown by the black, red, and blue curves, respectively. The left panels correspond to radius of a non-rotating NS R = 14 km and the right panels to R = 12 km. The mass of the non-rotating NS was assumed M = 1.5 M⊙ and both families of models have the same baryonic mass. The solid curves are for emission pattern corresponding to the isotropic intensity (4), while the dashed curves are for the electron-scattering dominated atmosphere with the intensity given by Eq. (6). The colour-correction factor is indistinguishable for the two considered emission patterns. Middle panels: correction factor to the blackbody normalization (dilution factor) as a function of the rotational frequency. Bottom panels: the correction to the bolometric flux. Same notations as in the top panels.

In the text

Fig. 5. Images of a NS rotating at 700 Hz. The mass and the radius of a non-rotating NS of the same baryonic mass are M = 1.5 M⊙ and R = 14 km. The white curves and the colours show the contours of constant bolometric intensity. Upper panels: blackbody case of local isotropic intensity. Lower panels: electron-scattering dominated atmosphere. Panels from the left to the right correspond to the inclinations i = 0°, 45°, and 90°. The lines of constant latitudes (every 10°) and longitudes (every 15°) are shown in black. The colour scheme is different for the two rows.

In the text

	Fig. 6. Qualitative difference between the shapes of the extended photospheres (dashed contours) at the touch-down point for two different hypothesis on the latitudinal distribution of the radiation flux over the surface of a rapidly rotating NS, CRF case (top panel) and CAF case (bottom panel).
In the text

Fig. 7. Spectra of the rotating NS with the parameters M′ = 1.519 M⊙, Re = 15.48 km, and ν* = 700 Hz computed for two inclination angles i= 0° (blue solid curve) and i= 90° (red solid curves), and for two relative luminosities (same over the surface) ℓ= 0.9 (top panel) and ℓ= 0.1 (bottom panel). The solid curves represent spectra computed taking accurate local spectra for atmosphere models from Suleimanov et al. (2012). The dashed curves correspond to the diluted blackbody approximation of the local spectra. Solar chemical composition of the atmosphere is assumed.

In the text

Fig. 8. Dependences of and w′ on the relative luminosity for a rapidly rotating NS for three inclination angles. The corresponding dependences for a non-rotating NS with the surface gravity log g0= 14.15 are shown with the dashed curves. Here chemical composition is solar with the heavy element abundances reduced by 100 times. The constant relative flux distribution over the NS surface is assumed (CRF case).

In the text

	Fig. 9. Same as Fig. 8, but computed with constant absolute flux distribution over the NS surface (CAF case).
In the text

Fig. 10. Top panel: results of fitting the w − ℓ dependence with w′−ℓ′ dependences computed for different inclination angles i and assuming constant relative flux over the NS surface (CRF case, see Fig. 8). The corresponding correction factors are (AF, AR) = (1.053, 0.856) for i = 0°, (1.1, 0.972) for 45°, and (1.15, 1.12) for 90°. Bottom panel: corresponding solutions at the (M, R) plane obtained assuming a non-rotating NS (dashed curve) and using additional correction factors AF and AR for different inclination angles.

In the text

	Fig. 11. Same as Fig. 10, but for constant absolute flux over the NS surface (CAF case, see Fig. 9). The corresponding correction factors are (AF, AR) = (1.12, 0.85) for i = 0°, (1.16, 0.97) for 45°, and (1.2, 1.12) for 90°.
In the text

Fig. 12. Confidence regions at the M–R plane for SAX J1810.8−2609 obtained from the χ2 map using the cooling tail method modified for rapid rotation. The solid, dashed, and dotted curves indicate the 68%, 90%, and 99% regions, respectively. The red contours correspond to the inclination i = 60°, while the blue contours are for the face-on NS. The black contours give the constraints assuming a non-rotating NS. (a) CRF case, i.e. constant relative flux. (b) CAF case, i.e. constant absolute flux. The triple dot-dashed curves are the contours of the constant distance of 4.0 and 4.5 kpc (distance is the only free parameter for a given pair of M and R).

In the text

Fig. A.1. Top panel: shape of the rotating NS computed with two different approximations. The solid curves correspond to Eq. (A.13), while the dashed curves to Eqs. (8)–(10) in Morsink et al. (2007). Bottom panel: effective surface gravity distributions. The solid curves correspond to Eq. (A.14), while the dashed curves to Eq. (50) in AlGendy & Morsink (2014). The red curves correspond to a relatively slow rotating NS with (ν* = 842 Hz, M′ = 1.8 M⊙, Re = 13.6 km). The blue curves are for a faster rotating NS with (ν* = 1084 Hz, M′ = 1.99 M⊙, Re = 15.1 km).

In the text

	Fig. B.1. Geometry of the problem.
In the text

Fig. B.2. Line profiles from a NS with M′ = 1.4 M⊙, Re = 10 km rotating at 700 Hz, and seen by a distant observer at inclination i = 20°. The line profiles for a spherical NS with R = Re and for an oblate NS shape both with Schwarzschild exterior metric are shown by the black and blue curves, respectively. Accounting for frame dragging and quadrupole moment gives the profile shown by the red curve. The line obtained with a quadrupole moment artificially increased by a factor of 4 above that given by Eq. (A.9) is shown by the dashed red curve. Compare to Fig. 5 in Nättilä & Pihajoki (2018).

In the text

	Fig. B.3. Line profiles observed at different inclination angles, i = 90° (black), 60° (orange), 40° (violet), 20° (magenta), 10° (blue), and 5° (red). The adopted NS parameters are M′ = 1.5 M⊙, Re = 14 km, ν* = 600 Hz. The energy range is shifted by the redshift factor . Compare to Fig. 6 in Nättilä & Pihajoki (2018).
In the text