© ESO 2018

1 Introduction

Accurate modeling of the emission from and around compact objects is a complicated combination of radiative processes and relativity. Not only is the object curving the space-time around it and hence affecting the trajectory of photons, but it can also affect the apparent observed radiation as the emitting surface can move with relativistic velocities. Existing convenient and modern frameworks for the emission around rotating (typically Kerr) black holes include GEOKERR (Dexter & Agol 2009), GYOTO (Vincent et al. 2011), GRAY (Chan et al. 2013), PANDURATA (Schnittman & Krolik 2013), ASTRORAY (Shcherbakov & McKinney 2013), HEROIC (Narayan et al. 2016), ODYSSEY (Pu et al. 2016), and GRTRANS (Dexter 2016), to name a few. Here we instead focus on the emission from rotating neutron stars, for which the Kerr metric is not a good appro- ximation if the star is rotating rapidly. By introducing BENDER¹, we aim to provide a similar publicly available platform for ray tracing problems focused on spinning compact objects.

Following the path of photons in a space-time of a rotating neutron star is a challenging task, both theoretically and numerically. Current and future observations, on the other hand, demand highly accurate models. For example, computing accurate pulse profiles of hot spots on spinning neutron stars has recently been intensively investigated, motivated by many upcoming or planned new space-borne X-ray observatories like ESA’s XIPE (Soffitta et al. 2013), eXTP of the China National Space Administration (CNSA; Zhang et al. 2016), and the already deployed Astrosat (Agrawal 2006) of the Indian Space Research Organization (ISRO) and NASA’s NICER (Gendreau et al. 2012). The expectation is that we may be able to better constrain the unknown neutron star (NS) equation of state (EoS) with accurate pulse profile observations, using the information encoded in the radiation (see, e.g., Lo et al. 2013).

Previous studies of emission from NSs are mainly formulated in a way that uses a non-rotating curved space-time metric for the bending of the photon paths, with special relativistic corrections to the rotational effects added separately (see, e.g., Pechenick et al. 1983; Page 1995; Miller & Lamb 1998; Weinberg et al. 2001; Poutanen & Gierliński 2003; Poutanen & Beloborodov 2006; Lamb et al. 2009a,b; Lo et al. 2013; Miller & Lamb 2015, but also Braje et al. 2000 for an alternative treatment). Space-times in these studies were also typically described by the spherically symmetric Schwarzschild metric.

As it has turned out, however, fast rotation can be a serious complication when considering the observed emission. First of all, because a finite pressure supports the rotating star, the star is squeezed into an oblate spheroid, and the oblateness increases with increasing rotation rate (Cook et al. 1994; Morsink & Stella 1999; Morsink et al. 2007; Bauböck et al. 2013a; AlGendy & Morsink 2014). This bulge in the equator will then distort the gravitational field outside the star. In Newtonian theory, the next-order correction to a non-spherical object (with azimuthal rotational symmetry and reflection symmetry along the equatorial plane) is defined by the quadrupole moment (see, e.g., Laarakkers & Poisson 1999; but also Pappas & Apostolatos 2012). Introduction of rapid rotation will then not only make the star oblate, but will also make the exterior space-time latitude dependent. Pioneering work in computing pulse profiles of such objects was made by Cadeau et al. (2005, 2007). Recently, a general ray tracing formulation in Hartle-Thorne metric-variant was given in Psaltis & Johannsen (2012) and Bauböck et al. (2012). This formulation was used to compute pulse profiles in Psaltis & Özel (2014). Here we seek to provide a similar, but open and publicly available code for solving similar types of problems. Moreover, we focus on building a connection between theprevious special relativistic formulations where different rotational effects are added separately by hand, and the full rotating general relativistic formulations where all of these effects naturally emerge from the theory.

The main focus of our framework is on the X-ray emission from accreting millisecond pulsars (AMPs; Wijnands & van der Klis 1998; Patruno & Watts 2012) and nuclear-powered millisecond pulsars (Watts 2012). We stress, however, that the whole framework presented in this paper is general enough to be applied to any problem of radiation originating from the vicinity of rotating compact objects. The radiation from AMPs emerges from hot spots on the surface of a rapidly rotating neutron star. The spots are heated by the infalling accreted material, which is being channeled to the magnetic poles by the neutron star’s magnetic field. The magnetic axis does not need to coincide with the rotational axis of the star, and hence pulsations can be observed from the spots that are rotating around the star. In the case of nuclear-powered millisecond pulsars, quasi-coherent oscillations are observed during a thermonuclear type I X-ray burst. The mechanism producing the pulses is, however, very similar to the case of AMPs, as an asymmetric bright patch in the burning surface layer is the origin of the observed pulsation. Accretion can also spin up these objects into extreme rotational velocities: spin frequencies of up to 620 Hz have been verified (4U 1608−52; Muno et al. 2002), whereas even a typical source has a spin around 400−500 Hz (Watts 2012; Papitto et al. 2014). Hence, if accurate emission is to be studied from these sources, one has to take the oblate shape and (in some cases) the second-order corrections to the space-time into account.

The paper is structured as follows. In Sect. 2 we introduce the framework of formulae and the theoretical background needed to compute the emission. We also describe the numerical methods used to solve the system of equations and present the publicly available code BENDER, which implements this framework. Next, we apply the code to various physical problems in Sect. 3. We compare our computations with results in the literature, when possible, to verify our calculations. Finally, in Sect. 4 we summarize our work.

2 Theory

2.1 Space-time metric

In our following derivations we use geometric units where G = c = 1 for the gravitational constant G and the speed of light c. We also assume the metric signature of (−, +, +, +) following the Misner et al. (1973) sign convention. Additionally, for the actual numerical calculations in the code, we set GM∕c² = 1, hence describing lengths in units of gravitational radius, where M is the mass of the compact object.

The exterior space-time of a static, non-rotating, spherically symmetric mass is described by the well-known Schwarzschild metric $$d{s}^2=-(1-u)d{t}^2+{(1-u)}^{-1}d{r}^2+r^2(d{\theta }^2+{\sin }^2\theta d{\varphi }^2),$$(1)

where r is the radial coordinate defined so that the area of a sphere at coordinate time t is 4πr², and we set u ≡ 2M∕r.

This metric is equivalent to an alternative solution known as isotropic Schwarzschild metric (see, e.g., Misner et al. 1973) $$d{s}^2=-{\left(\frac{1-\frac{\overline{u}}{2}}{1+\frac{\overline{u}}{2}}\right)}^{2}d{t}^2+{\left(1+\frac{\overline{u}}{2}\right)}^4(d{\overline{r}}^2+{\overline{r}}^2(d{\theta }^2+{\sin }^2\theta d{\varphi }^2)),$$(2)

where $\overline{r}$ is the so-called isotropic radial coordinate, and we set $\overline{u}\equiv M/\overline{r}$. This kind of isotropic metric has the useful feature that surfaces of constant time are conformally flat, and hence the angles are represented without distortion. However, this also means that angular isotropic coordinates do not faithfully represent the distances within the spheres, nor does the radial coordinate correspond directly to the radial distance. From here on, we mark all variables related to the isotropic radial coordinate with a bar on top.

We consider a rotating compact object. To describe our system, we need a dimensionless angular velocity $$\widehat{\Omega }=\Omega {\left(\frac{R_{\text{e}}{}^3}{M}\right)}^{1/2},$$(3)

where Ω is the angular velocity of a sphere with an equatorial radius R_e and a mass M scaled with the Newtonian mass shedding (Kepler) limit ${(M/R_{\text{e}}{}^3)}^{1/2}$ (see Friedman & Stergioulas 2013). Here R_e is described using the usual Schwarzschild radial coordinate, and it corresponds to the equatorial radius of the star for which 2πR_e gives the proper length of the circumference in the rotational equator as measured in the local static frame. The asymptotically flat metric near a stationary axisymmetrically rotating object in isotropic form is (Bardeen & Wagoner 1971) $$\begin{array}{ll}d{s}^2\hfill & =-e^{2\overline{\nu }}d{t}^2+{\overline{r}}^2{\sin }^2\theta {\overline{B}}^2{e}^{-2\overline{\nu }}{(d\varphi -\overline{\omega }dt)}^2\hfill \\ \hfill & +e^{2(\overline{\zeta }-\overline{\nu })}(d{\overline{r}}^2+{\overline{r}}^{2}d{\theta }^2),\hfill \end{array}$$(4)

where $\overline{\omega }$ is the angular velocity of the local inertial frame, and the functions $\overline{\nu }$, $\overline{B}$ and $\overline{\zeta }$ in the metric coefficients can be expanded in the powers of $\widehat{\Omega }$ and ū (Butterworth & Ipser 1976). Here $e^{-\overline{\nu }}$ is the time-dilation factor relating the proper time of the local observer to the time at infinity. Physical interpretation of $\overline{B}$ follows from the fact that the proper circumference of a circle around the axis of symmetry is $2\pi (e^{-\overline{\nu }}\overline{B}\overline{r}\sin \theta )$. Similarly, the interpretation of $\overline{\zeta }$ follows from the fact that $e^{\overline{\zeta }-\overline{\nu }}$ acts as a conformal (angle preserving) factor of the space-time. We also note that the time and space coordinates are connected in the isotropic metric via the $\overline{\nu }$-term that also enters both the radial and angular terms. The zeroth-order terms ($\widehat{\Omega }=0$) of the series expansions are the familiar Schwarzschild metric coefficients expressed in isotropic coordinates (see Table 1).

The first-order expansion in rotation (${\widehat{\Omega }}^1$) is qualitatively related to Kerr metric. In this case, we introduce an angular velocity term of the local inertial frame $\overline{\omega }$ that accounts for the frame-dragging effects. It can be defined as $$\overline{\omega }=\frac{2j}{M}({\overline{u}}^3-3{\overline{u}}^4),$$(5)

where the dimensionless quantity $j=\mathcal{J}/M^2$ and $\mathcal{J}=I\Omega$ is the star’s angular momentum with moment of inertia I(M, Ω).

The second-order expansion (${\widehat{\Omega }}^2$) corresponds to a similar approximation as the Hartle-Thorne slow-rotation space-time (Hartle & Thorne 1968), which introduces two quadrupole moments into the metric. These second-order multipole moments can be defined via the dimensionless quantities q and β, the dimensionless moments of energy density and pressure, respectively. These are, however, dependent on the selection of a coordinate system. The coordinate invariant quadrupole moment is a combination of these two quantities and is given in Pappas & Apostolatos (2012; see Eq. (11) therein; see also AlGendy & Morsink 2014 and Eq. (18)) as $$q_{\text{inv}}=q+\frac{4}{3}\beta .$$(6)

This detail should be taken into account when comparing the strength of the quadrupole deviations between different metric descriptions.

Yagi & Yunes (2013) showed that when the NS mass, radius, and spin are known, the star’s structure (and hence the surrounding metric) is almost fully characterized by these three quantities alone. In other words, regardless of the unknown microphysics of the underlying matter, the NS parameters (such as moment of inertia and coordinate-invariant quadrupole moment) are connected by what is called (approximative) universal relations. Bauböck et al. (2013a) defined empirical relations for these parameters in the Hartle-Thorne metric based on their computations of rotating NSs with various different EoS. AlGendy & Morsink (2014) later refined these relations for the metric representation (Eq. (4)) by Butterworth & Ipser (1976). In practice, we can then parameterize the previously presented quantities with great accuracy by using only the dimensionless angular velocity $\widehat{\Omega }$ and the compactness parameter x. We note that these two parameters are defined in terms of the equatorial circumferential radius R_e defined in the usual Schwarzschild coordinate system. To the lowest order, these parameterizations are (AlGendy & Morsink 2014) $$\begin{array}{l}q=-0.11\frac{{\widehat{\Omega }}^2}{x^2},\\ \beta =0.4454{\widehat{\Omega }}^{2}x,\end{array}$$

and $$I=\sqrt{x}(1.136-2.53x+5.6x^2)M{R}_{\text{e}}{}^2.$$(9)

We also note that both q and β are $\mathcal{O}({\widehat{\Omega }}^2)$, whereas j is $\mathcal{O}(\widehat{\Omega })$.

It is possible to transform between the Schwarzschild coordinate radius r and the isotropic $\overline{r}$ coordinates using the relation (Friedman et al. 1986) $$r=\overline{B}{e}^{-\overline{\nu }}\overline{r}.$$(10)

The relation between the differentials of the two radial coordinates is $$dr=e^{\overline{\zeta }}d\overline{r},$$(11)

which can be computed using the series representation of Butterworth & Ipser (1976).

Since the series expansions of the metric coefficients are expressed in terms of the isotropic radial coordinate $\overline{r}$, we favor this notation in our derivation. However, in some cases, we simplify the equations into a more intuitive form using the Schwarzschild radial r coordinate.

2.2 Oblate shape of the neutron star

Because of the rotation and finite pressure supporting the NS, it is not a perfect sphere when it is rotating. However, it retains axisymmetry and can be approximated with an oblate spheroid. Similarly to the Eqs. (7)–(9), AlGendy & Morsink (2014) constructed an approximate formulation for the shape of thesurface of a rotating neutron star. It is given by expressing the radius as a function of colatitude θ as $$\begin{array}{ll}R(\theta )\hfill & =R_{\text{e}}\left(1-\frac{R_{\text{e}}-R_{\text{p}}}{R_{\text{e}}}{\cos }^2\theta \right)\hfill \\ \hfill & =R_{\text{e}}[1-{\widehat{\Omega }}^2(0.788-1.03x){\cos }^2\theta ],\hfill \end{array}$$(12)

where R(π∕2) = R_e is the radius of the star in its rotational equator, and R_p is its radius as measured along the rotation axis.

The elemental surface area for a spheroid is given as (using the usual Schwarzschild coordinates) $$dS(\theta )=R^2(\theta )\sin \theta \sqrt{1+f{(\theta )}^2}d\theta d\varphi ,$$(13)

where $$f(\theta )=\frac{1}{\overline{R}(\theta )}\frac{d\overline{R}(\theta )}{d\theta }=\overline{B}{e}^{-\overline{\zeta }-\overline{\nu }}\frac{1}{R(\theta )}\frac{dR(\theta )}{d\theta },$$(14)

and $\overline{B}{e}^{-\overline{\zeta }-\overline{\nu }}\approx {\left(1-\frac{2M}{r}\right)}^{-1/2}+\mathcal{O}({\widehat{\Omega }}^2)$. The angle γ, defined as the angle between the radial unit vector r and the surface normal n is given by $$\cos \gamma ={\left(1+f{(\theta )}^2\right)}^{-1/2}.$$(15)

Then the normal to the surface can be defined using the radial vector r and the tangential vector θ as $$n=\cos \gamma r+\sin \gamma \theta .$$(16)

See Fig. 1 for a clarification of the angles.

Fig. 1 Geometry of the system. Here we show the underlying spherical coordinate system with ϕ and θ coordinates along with the Cartesian (x, y, z) coordinate system. In addition, the observer’s $\widehat{x}-\widehat{y}$ image plane is shown as viewed from an inclination angle i. The star is taken to rotate rapidly around the y -axis, which leads to an oblate (squeezed) shape for the emitting surface, and the radial vector r and the surface normal n therefore start to differ from each other.

2.3 Geodesic motion using the Hamilton-Jacobi equation

In general, the motion of light rays in curved space-time is governed by the second-order geodesic equation. In this section we present an equivalent theoretical formalism based on Hamilton-Jacobi (sometimes also known as super-Hamiltonian) description (Misner et al. 1973; Chandrasekhar 1998). The advantage of this alternative representation is its physical intuitiveness as the formalism relies heavily on identifying and using the constants of motion of the problem. In the end, when we apply our methods, we use both approaches (the first-order Hamilton-Jacobi equations and the second-order geodesic equations) in our calculations to show that for physically relevant problems, the results obtained are equivalent up to good numerical precision. In a typical situation, we therefore apply the Hamilton-Jacobi method presented here, and when accuracy is the main factor, we fall back to solving the full geodesic equations.

We now discuss the motion of particles in curved space-time. Geodesic motion in a space-time characterized by a metric g_ij is governed by the Hamilton-Jacobi equation $$2\frac{\partial S}{\partial \tau }=g^{ij}\frac{\partial S}{\partial x^i}\frac{\partial S}{\partial x^j},$$(17)

where g^ij is the inverse metric and S denotes the Hamilton principal function. For the two Killing vectors t^α = (1, 0, 0, 0) (asymptotic time symmetry) and ϕ^α = (0, 0, 0, 1) (axisymmetry about the rotational axis) in rotating space-time, the Frobenius theorem implies the existence of a family of two surfaces orthogonal to these vectors (see, e.g., Friedman & Stergioulas 2013). This means that there are surfaces of constant t and ϕ in our space-time, yielding two constants of motion, namely energy E and the z-component of the angular momentum, L_z. We then seek a solution of Eq. (17) in the form $$S=\frac{1}{2}{\delta }_1\tau -Et+L_z\varphi +S_{\overline{r}}(\overline{r})+S_{\theta }(\theta ),$$(18)

where δ₁ is related to the rest-mass of the particle we study. With the metric function Eq. (4), this becomes $$\begin{array}{ll}{\delta }_1{\overline{r}}^2{\overline{B}}^2{e}^{-2\overline{\nu }}= \hfill & {\overline{r}}^2{\overline{B}}^2{e}^{-2\overline{\zeta }}{\partial }_{\overline{r}}{S}_{\overline{r}}{S}^2-e^{-4\overline{\nu }}{\overline{B}}^2{\overline{r}}^2{(E-L_z\overline{\omega })}^2\hfill \\ \hfill & +{\overline{B}}^2{e}^{-2\overline{\zeta }}{\partial }_{\theta }{S}_{\theta }^2+\frac{L_z^2}{{\sin }^2\theta }.\hfill \end{array}$$(19)

After reorganizing terms and introducing a simplifying notation $e^{\overline{\zeta }}/\overline{B}\equiv e^{\overline{\mu }}$, we obtain $$\begin{array}{ll}\hfill & e^{-2\overline{\mu }}{\partial }_{\theta }{S}_{\theta }^2+\frac{L_z^2}{{\sin }^2\theta }=\hfill \\ \hfill & {\overline{B}}^2{e}^{-2\overline{\nu }}{\overline{r}}^2(e^{2(\overline{\nu }-\overline{\zeta })}{\partial }_{\overline{r}}{S}_{\overline{r}}^2-{\delta }_1-e^{-2\overline{\nu }}{(E-L_z\overline{\omega })}^2).\hfill \end{array}$$(20)

The individual terms in Eq. (20) only depend on r or only on θ, that is, the dependence on r and θ is separable if $\overline{\nu }$, $\overline{B}$, and $\overline{\zeta }$ only depend on r and μ only depends on θ. This is the case to first order in the stellar rotation rate $\widehat{\Omega }$ because $e^{\overline{\mu }}=1+\mathcal{O}({\widehat{\Omega }}^2)$, in addition to $\overline{\nu }={\overline{\nu }}_0(\overline{r})+\mathcal{O}({\widehat{\Omega }}^2)$ and $\overline{B}={\overline{B}}_0(\overline{r})+\mathcal{O}({\widehat{\Omega }}^2)$ (see Table 1). However, to second-order in $\widehat{\Omega }$, the individual terms in Eq. (20) depend on both r and θ, that is, the dependence on r and θ is not separable. For geodesics, however, these higher-order deviations only contribute very close to the actual NS surface, and neglecting them enables us to obtain accurate approximations of the photon path.

When we now assume separability, we can introduce a separation variable $\mathcal{C}$ known as Carter’s constant (third constant of motion) in order to solve the differential Eq. (20). By noting that the conjugate momenta correspond to the first derivatives of S with respect to the generalized coordinates, we can write the components of four-momentum p as $$\begin{array}{ll}{p}_t\hfill & =-E,\hfill \\ p_{\overline{r}}\hfill & =\pm e^{\overline{\zeta }-2\overline{\nu }}{\left({\delta }_1{e}^{2\overline{\nu }}+{(E-L_z\overline{\omega })}^2-\frac{\mathcal{C}}{{\overline{B}}^2{e}^{-4\overline{\nu }}{\overline{r}}^2}\right)}^{1/2},\hfill \\ p_{\theta }\hfill & =\pm e^{\overline{\mu }}{\left(\mathcal{C}-\frac{L_z^2}{{\sin }^2\theta }\right)}^{1/2},\hfill \\ p_{\varphi }\hfill & =L_z.\hfill \end{array}$$

Similarly, the components of a local tetrad frame are $$\begin{array}{ll}{p}^{(t)}\hfill & =-p_{(t)}=-e_{(t)}^{\widehat{\mu }}{p}_{\widehat{\mu }}=-e^{-\overline{\nu }}{p}_t,\hfill \\ p^{(\overline{r})}\hfill & =p_{(\overline{r})}=e_{(\overline{r})}^{\widehat{\mu }}{p}_{\widehat{\mu }}=e^{-\overline{\zeta }+\overline{\nu }}{p}_{\overline{r}},\hfill \\ p^{(\theta )}\hfill & =p_{(\theta )}=e_{(\theta )}^{\widehat{\mu }}{p}_{\widehat{\mu }}=\frac{1}{\overline{r}}{e}^{-\overline{\zeta }+\overline{\nu }}{p}_{\theta },\hfill \\ p^{(\varphi )}\hfill & =p_{(\varphi )}=e_{(\varphi )}^{\widehat{\mu }}{p}_{\widehat{\mu }}=\frac{1}{e^{-\overline{\nu }}\overline{B}\overline{r}\sin \theta }{p}_{\varphi },\hfill \end{array}$$

where $e_{(a)}^{\widehat{\mu }}$ with index $a=t,\overline{r},\theta$ and ϕ are the tetrads of metric Eq. (4). Since we only consider null geodesics (i.e., photons), we now set δ₁ = 0.

2.4 Photon ray tracing

We now consider radiation that is emitted from the surface of the star at an emission point (r_e, θ_e, ϕ_e), as seen in the static frame. The radiation travels along a geodesic with a specific intensity I_E as measured by an observer comoving with the emission point. It is observed at an image plane situated at a radial distance r, with r →∞. We then wish to calculate the projected image of the star at this image plane.

First we set up the coordinate system so that the plane of observation is toward ϕ = 0 and θ = i, where i is the angle of inclination (see Fig. 1). The geodesic will be emitted with a four-momentum ${\underset{\_}{p}}_{\text{e}}$, and if it is eventually observed at the image plane at infinity, it will have a final four-momentum of $(E,{\widehat{p}}_r,0,0)$, purely in the radial direction. Likewise, the components of the position must satisfy $$\begin{array}{ll}\theta \hfill & \to i,\hfill \\ \varphi \hfill & \to 0,\hfill \end{array}$$

as r→∞. The change in the time and angular components along the geodesic can be written as $$\begin{array}{ll}dt\hfill & =\frac{p^t}{p^{\overline{r}}}\text{d}\overline{r},\hfill \\ d\theta \hfill & =\frac{p^{\theta }}{p^{\overline{r}}}\text{d}\overline{r},\hfill \\ d\varphi \hfill & =\frac{p^{\varphi }}{p^{\overline{r}}}\text{d}\overline{r},\hfill \end{array}$$

yielding a total change of angles Δθ and Δϕ when integrating from r_e to ∞. The condition for being observed is then $$\begin{array}{ll}{\theta }_{\text{e}}+\Delta \theta \hfill & =i,\hfill \\ {\varphi }_{\text{e}}+\Delta \varphi \hfill & =0.\hfill \end{array}$$

The projected image of the star on the image plane can then be described by two celestial coordinates: abscissa $\widehat{x}$ and ordinate ŷ. Making use of the tetrad components Eqs. (25)–(28), we obtain (Chandrasekhar 1998) $$\widehat{x}={\left(\frac{r{p}^{(\varphi )}}{p^{(t)}}\right)}_{r\to \infty }=\frac{1}{\sin i}\frac{L_z}{E},$$(36)

and $$\widehat{y}={\left(\frac{r{p}^{(\theta )}}{p^{(t)}}\right)}_{r\to \infty }=\frac{\sqrt{\mathcal{C}-\frac{L_z^2}{{\sin }^{2}i}}}{E}.$$(37)

Here it is useful to transform into a polar coordinate system on the image plane, as Eqs. (36) and (37) strongly suggest a more intuitive form if this is done. In this system we use as coordinates the radial distance from the center point, or the impact parameter b, and the polar rotation angle χ. We take χ to increase clockwise from the projected spin axis of the neutron star, with χ = 0 corresponding to the projected direction from the south to the north pole of the neutron star. We then express the impact parameter b and the angle χ via L_z and $\mathcal{C}$ as $$b=\frac{\sqrt{\mathcal{C}}}{E},$$(38)

and $$\sin \chi =\frac{1}{\sin i}\frac{L_z}{\sqrt{\mathcal{C}}}.$$(39)

Here, the nature of Carter’s constant as a generalized squared angular momentum is apparent. The constants of motion, combined with the geodesic null condition p^μp_μ = 0, allow us to solve p^θ and p^ϕ in terms of $\overline{r}$. As the final step, we can substitute the four-momentum components and image-plane coordinates into Eqs. (31)–(33) and solve the system of three first-order differential equations (in terms of t, θ, and ϕ) with $\overline{r}$ as a variable.

2.5 Redshift and emission angle

It is the most convenient to define all radiative processes in the co-rotating frame of the star. We denote variables defined in a frame that is momentarily comoving with the stellar surface with a prime. On the other hand, our distant observer is stationary and moving along the timelike Killing vector. Hence, we need to transform between stationary and rotating frames by using the four-velocity of the star’s fluid. To make a connection to the theory of special relativity, it is convenient to define two frames: a corotating rest frame of the fluid K′, and a non-rotating static frame K. Laws of physics for the radiative transfer take the usual form in the K′. A stationary observer, on the other hand, is in the non-rotating frame K from where the fluid is seen to move relativistically. We therefore need to transform between these two frames. In addition, we need to take into account that a particle released from infinity with zero angular momentum will acquire non-zero angular velocity in the direction of the star’s rotation as a result of the dragging of inertial frames.

Four-velocity of a stationary observer with zero angular momentum (so-called ZAMO) is $o^{\alpha }=N_o(t^{\alpha }+\overline{\omega }{\varphi }^{\alpha })$, where the normalization factor $N_o=e^{-\overline{\nu }}{(1-v_{\omega }{}^2)}^{-1/2}$ is obtained from o_αo^α = −1, and the velocity of the frame is $$v_{\omega }=\overline{\omega }\overline{B}{e}^{-2\overline{\nu }}\overline{r}\sin \theta ,$$(40)

indicating that the angular velocity of the ZAMO as measured by an inertial observer at infinity is $o^{\varphi }/o^t=d\varphi /dt=\overline{\omega }$.

The four-velocity s^α of a circular flow can be defined using the timelike and rotational Killing vectors as s^α = N_s(t^α + Ωϕ^α), where the normalization factor is defined as $N_s=e^{-\overline{\nu }}{(1-v_{\text{z}}{}^2)}^{-1/2}$ determined by s_αs^α = −1. Here the velocity $$v_{\text{z}}=(\Omega -\overline{\omega })\overline{B}{e}^{-2\overline{\nu }}\overline{r}\sin \theta ,$$(41)

can be identified as the three-velocity measured in the frame of the ZAMO, an observer rotating with a velocity of v_ω.

The total redshift is then given by an inner-product between a photon u^α and a four-velocity of the star’s fluid s^α. With these definitions, the redshift is $$1+z=-s_{\alpha }{u}^{\alpha }=e^{-\overline{\nu }} {\delta }^{-1},$$(42)

consisting of the gravitational part $e^{\overline{\nu }}$ and of the Doppler-like factor $$\delta =\frac{\sqrt{1-v_{\text{z}}{}^2}}{1-\Omega L_z}.$$(43)

To compute the emission angle, we again have to take the rotating frame into account. This can be done by introducing a projection operator $${\overline{h}}_{ab}=g_{ab}+s_a{s}_b,$$(44)

which projects four-vectors from the non-rotating frame to the rotating frame where the radiative processes are defined. Our metric tensor for the rotating observer is then $${\overline{h}}_{ab}d{x}^{a}d{x}^b=e^{2(\overline{\zeta }-\overline{\nu })}(d{\overline{r}}^2+{\overline{r}}^{2}d\theta )+\frac{{\overline{B}}^2{e}^{-2\overline{\nu }}{\overline{r}}^2{\sin }^2\theta }{1-v_Z^2}{(d\varphi -\Omega dt)}^2.$$(45)

As a definition, we can take the emission angle to be the angle between photon and a space-like surface normal vector $n_{\alpha }=N_n(0,\cos \gamma ,\overline{r}\sin \gamma ,0)$, with normalization $N_n=e^{\overline{\zeta }-\overline{\nu }}$ fulfilling n_α n^α = −1. The line element in spherical coordinates is $ds=dr\widehat{r}+rd\theta \widehat{\theta }+r\sin \theta d\varphi \widehat{\varphi }$, and by combining this with Eq. (16), we obtain the presented surface normal. When projected to the rotating frame, we can then obtain the angle from the generalized dot-product definition between two vectors as $$\cos {\alpha }^{\prime }=\frac{{\overline{h}}_{ab}{n}^a{u}^b}{{({\overline{h}}_{ab}{n}^a{n}^a)}^{1/2}{({\overline{h}}_{ab}{u}^a{u}^b)}^{1/2}}.$$(46)

Using the metric defined by Eq. (45) and a photon with components Eqs. (21)−(24), we obtain $$\cos {\alpha }^{\prime }=\delta e^{2\overline{\nu }-\overline{\zeta }}\left[p_{\overline{r}}\cos \gamma +\frac{p_{\theta }}{\overline{r}}\sin \gamma \right].$$(47)

For the non-rotating observer, we similarly obtain $$\cos \alpha =\cos {\alpha }^{\prime }{\delta }^{-1},$$(48)

by setting Ω → 0. Here it suffices to notice that we are only interested in the emission angle value at the surface of the star, that is, $\overline{r}=\overline{R}(\theta )$. The result here is identical to the emission angle obtained with a special-relativistic approach, using flat-space trigonometry and Lorentz-boosting with δ to the rotating frame (see, e.g., Poutanen & Beloborodov 2006).

2.6 Corotating coordinates

Next we define some quantities for a corotating observer located at the surface of the star. This helps us to connect the previously presented backward-in-time method to the methods where light rays are propagated from the star to the image plane (forward-in-time methods).

Transforming from the observer’s non-rotating frame K to the fluid rest frame K′ is easily done using the previously defined projection operator ${\overline{h}}_{ab}$ given by Eq. (45). We next express this projection operator in the normal coordinate system. Using the Eqs. (10) and (11), we obtain a longitudinally Lorentz-boosted metric tensor $$h_{ab}d{x}^{a}d{x}^b=e^{-2\overline{\nu }}d{r}^2+\frac{e^{2\overline{\zeta }}}{{\overline{B}}^2}{r}^{2}d{\theta }^2+{\gamma }_{\text{L}}^2{r}^2{\sin }^2\theta {(d\varphi -\Omega dt)}^2.$$(49)

The result agrees with the Schwarzschild coordinate system up to first order in rotation because $e^{\overline{\zeta }}/\overline{B}\approx 1+\mathcal{O}({\widehat{\Omega }}^2)$. This notation can be further simplified by defining a new azimuthal angular coordinate as ϕ′ ≡ ϕ −Ωt that is to be used by the rotating observer². It is important to note here that because of the rotation, the azimuthal angle ϕ and the time t are now coupled for the comoving observer. This new normalized longitudinal coordinate ensures that both the rotating and the stationary observer agree that a circle drawn around the star has 2π radians. The expression for the new longitudinal coordinate is also seen to be Lorentz-stretched by a factor of ${\gamma }_{\text{L}}={(1-v_{\text{z}}{}^2)}^{-1/2}$. The corotating observer can then use this projection operator to define space-like vectors orthogonal to their world line.

Additionally, it is useful to consider another projection operator m that will project from the 3D space to the 2D surface of the star. We can define this projection as $$m_{ab}=g_{ab}+s_a{s}_b-n_a{n}_b=h_{ab}-n_a{n}_b,$$(50)

where n_a is the unit normal to the surface. From here, it is easy to verify that it is perpendicular to the surface of the star as m_ab n^a = 0 and to the velocity of the surface as m_abs^a = 0. For simplicity, we now consider a spherical star so that the surface normal reduces to $n^a=e^{\overline{\nu }}{\delta }_r^a$, where ${\delta }_a^b$ is the Kronecker delta function. Then m can be expressed as $$m_{ab}d{x}^{a}d{x}^b=\frac{e^{2\overline{\zeta }}}{{\overline{B}}^2}{r}^{2}d{\theta }^2+{\gamma }_{\text{L}}^2{r}^2{\sin }^2\theta {(d\varphi -\Omega dt)}^2.$$(51)

This is effectively a 2D metric tensor of the star’s surface that is perpendicular to the world line of the corotating observer.

Definitions here are also relevant for the so-called Schwarzschild+Doppler (S+D) approximation (see, e.g., Poutanen & Beloborodov 2006). In the S+D approximation, the observer’s polar coordinate plane (b, χ) is connected to the corotating spherical coordinates (ϕ′, θ′) of the star. We note that this connection is done algebraically in the S+D method and so there is no need to determine the full path of the ray using partial differential equations as the problem reduces to calculating the so-called lensing integral alone. The S+D calculations are done in a special relativistic framework where quantities are defined in a corotating coordinate system that are then Lorentz-boosted into the static non-rotating frame. We now show the correct expression of this change of frame so that the results match the backward-in-time method (see also Cadeau et al. 2007, for an alternative derivation). In the usual Schwarzschild metric, the impact parameter b can be obtained as a function of the emission angle α given by Eq. (48) as $$b=\frac{1}{\sqrt{1-u}}R\sin \alpha ,$$(52)

by setting cosγ → 1, sin γ → 0 (spherical star), and ω → 0 along with $e^{-2\nu }\to {(1-u)}^{-1/2}$ (Schwarzschild metric, i.e., $\mathcal{O}({\widehat{\Omega }}^1)$). Here we use the emission angle as measured by the non-rotating observer in frame K. In this case, the solid angle is then simplified to $$d{\Omega }_o=bdbd\chi =\frac{1}{1-u}\frac{d\cos \alpha }{d\cos \psi }{R}^2\cos \alpha \frac{d\chi d\cos \psi }{D^2},$$(53)

where ψ is the lensing angle. Using the spherical symmetry of the Schwarzschild metric, we can directly connect the ψ and θ along with χ and ϕ to obtain $$d{\Omega }_o=\frac{1}{1-u}\frac{d\cos \alpha }{d\cos \psi } \cos \alpha \frac{dS}{D^2},$$(54)

where dS = R² sinθdθdϕ is the area element for the non-rotating static observer in K. Using dΩ_o to compute the flux, we would then obtain the observed (received) flux valid for the observer in frame K.

Integration over some finite-sized features, on the other hand, is done on the surface of the star in the corotating frame K′. This results in a mixing of the different frames because we would like the integration to happen simultaneously for the observer in K, that is to say, we need a connection between t, ϕ, and ϕ′. To do this, we need to define the differential area element in the corotating frame. This can be obtained by considering the 2D metric tensor m of the surface given by Eq. (50). Using this, we can derive the corotating differential area element as $$d{S}^{\prime }=\sqrt{\det m_{ab}}d{x}^{a}d{x}^b=\frac{e^{\overline{\zeta }}}{\overline{B}}{\gamma }_{\text{L}}{R}^2\sin \theta d\theta (d\varphi -\Omega dt).$$(55)

This means that the rotation will result in a stretching of the area element by a factor of γ_L, whereas the quadrupole moments will deform it by a factor of $e^{\overline{\zeta }}/B$. This result is obtained purely from differential geometry.

Next we work only in the Schwarzschild metric to draw a direct connection to the S+D approximation. Using Eq. (55), we obtain $$d{S}^{\prime }={\gamma }_{\text{L}}{R}^2\sin \theta d\theta (d\varphi -\Omega dt)={\gamma }_{\text{L}}{R}^2\sin {\theta }^{\prime }d{\theta }^{\prime }d{\varphi }^{\prime },$$(56)

as given in the corotating coordinates defined as θ′≡ θ and ϕ′≡ ϕ −Ωt. This differs by a factor of γ_L from an incorrect result that would be obtained by erroneously assuming that dS′ = R² sinθ′dθ′dϕ′. To transform from the non-rotating K frame to the corotating frame K′, we can use the Lorentz invariance of the photon beam cross-section given as (Terrell 1959; Lind & Blandford 1985) $$\cos \alpha dS=\cos {\alpha }^{\prime } d{S}^{\prime }.$$(57)

The connection between the emission angles cosα′ and cosα is also known from Eqs. (47) and (48) and is seen to be a simple Doppler boost factor δ. We then obtain $$dS=\delta d{S}^{\prime }.$$(58)

Finally, the total observed angular size is then seen to be $$d{\Omega }_o=\frac{{\gamma }_{\text{L}}\delta }{1-u}\frac{d\cos \alpha }{d\cos \psi } \cos \alpha \frac{R^2\sin {\theta }^{\prime }d{\theta }^{\prime }d{\varphi }^{\prime }}{D^2}.$$(59)

2.7 Emission

The observed (i.e., received) flux at photon energy E from a small area on an image plane is $$d{F}_E=I_{E}d{\Omega }_o,$$(60)

where I_E is the specific intensity of the radiation at infinity, and dΩ_o is the solid angle subtended by the element as measured by the observer. The total flux is then the integral of these elements over the image plane. As a final step, this observed flux has to be connected to the actual emerging radiation.

From Eq. (42), the relation between the emergent energy E′ to the observed energy E is $E/E^{\prime }={(1+z)}^{-1}$. The connection between the monochromatic observed and local intensity is then (see, e.g., Misner et al. 1973; Rybicki & Lightman 1979) $$I_E={\left(\frac{E}{E^{\prime }}\right)}^3 {I^{\prime }}_{E^{\prime }}({\alpha }^{\prime }),$$(61)

where ${I^{\prime }}_{E^{\prime }}({\alpha }^{\prime })$ is the intensity computed in the frame comoving with the emitting area. The radiation here is emitted in the direction of the angle α′ defined in the local rotating frame. Integrating over the energies, we obtain the bolometric intensity $$I={\left(\frac{E}{E^{\prime }}\right)}^4 I^{\prime }({\alpha }^{\prime }).$$(62)

The total (monochromatic) flux as a function of the observer’s time F_E (t) can then be obtained by integrating over the whole image, $$F_E(t)={\displaystyle \int I_E}(t) \text{d}{\Omega }_o={\displaystyle \int {\displaystyle \int \frac{{I^{\prime }}_{E^{\prime }}(t_{*},{\alpha }^{\prime })}{{(1+z)}^3}}} \frac{b\text{d}b\text{\hspace{0.17em}d}\chi }{D^2},$$(63)

where t_* = t −Δt is the time when the photon was emitted as measured in the non-rotating coordinate system. This can be computed when we know the total travel time Δt against some reference photon, for example, the one with the shortest path to the observer.

All of these quantities on the (non-rotating) spherical coordinate system (θ, ϕ) are then mapped to the observer’s polar image coordinates (b, χ) via ray tracing. The original longitudinal coordinate of the emission is easily obtained from ϕ_e = ϕ − t_*Ω because both t_* and Ω are defined for a distant observer, and change in the azimuthal coordinate is Lorentz invariant. This allows us to connect the observables that our distant observer will see to the local rest frame of the gas where most of the physical processes are naturally defined.

2.8 Angular distribution of radiation

2.8.1 Blackbody radiation

For pulse profile calculations, the simplest angular distribution of radiation is the isotropic radiation. Here we consider blackbody emission described by the specific intensity $$B_{\text{E}}(T)=5.04\times {10}^{22}\frac{E^3}{\text{exp}(E/T)-1}\text{\hspace{0.17em}erg\hspace{0.17em}}{\text{s}}^{-1}{\text{cm}}^{-2}{\text{keV}}^{-1}{\text{sr}}^{-1},$$(64)

where T and E are given in keV.

2.8.2 Atmosphere dominated by electron scattering

Next we consider beamed radiation. For simplicity, we continue to assume that the spectral distribution is given by the Planck function, but now we assume that the angular distribution corresponds to that given by coherent electron scattering in a plane-parallel, semi-infinite (optical depth τ →∞) atmosphere. This beaming pattern is described by the so-called Hopf function H(μ). Introducing a variable μ ≡ cosα′, the result is $$I_{\text{E}}(\mu )=B_{\text{E}}(T)\frac{H(\mu )}{2{\alpha }_1},$$(65)

where $${\alpha }_n={\displaystyle {\int }_0^{1}H}(\mu ){\mu }^{n}d\mu ,$$(66)

are the moments of the function H(μ), which is a solution of the Ambartsumian-Chandrasekhar integral equation (see, e.g. Chandrasekhar 1960; Sobolev 1963) $$H(\mu )=1+\mu H(\mu ){\displaystyle {\int }_0^1\frac{\Psi (\eta )}{\mu +\eta }}H(\eta )d\eta .$$(67)

Here Ψ(μ) is the characteristic function, which depends on the scattering law considered. For Rayleigh (Thomson) scattering, it is $$\Psi (\mu )=\frac{3}{16}(3-{\mu }^2).$$(68)

Given Ψ, the integral Eq. (67) can then be iteratively solved, for example, by computing $$\begin{array}{ll}{H}_{n+1}(\mu )=\frac{1}{2}{H}_n(\mu )+\frac{1}{2}(\hfill & \sqrt{1-2{\displaystyle {\int }_0^1\Psi }(\eta )\text{d}\eta }\hfill \\ \hfill & +{\displaystyle {\int }_0^1\frac{\eta \Psi (\eta )}{\mu +\eta }}{H}_n(\eta )\text{d}\eta {)}^{-1},\hfill \end{array}$$(69)

with a starting guess of $$H_0(\mu )=1+2.3\mu -0.3{\mu }^2.$$(70)

We note that Eq. (65) is physically inconsistent (blackbody radiation must be isotropic), but we adopt this spectrum and beaming pattern for the sake of illustration. Electron-scattering atmospheres can produce spectra that have spectral shapes similar to a Planck function, but they are much less efficient. Eq. (65) can describe such emission approximately, but only if it is preceded by an efficiency factor that depends on the color-correction factor f_c as $f_{\text{c}}^{-4}\approx 0.15$ (see, e.g., Suleimanov et al. 2011, 2012). We also note that even the simple polynomial expansion Eq. (70) has an accuracy better than <2%, and it can therefore be used in approximate solutions with a corresponding first moment of α₁ = 1.19167.

2.9 Method of solution

In practice, when calculating the observed time-dependent emission, we have to

• 1.

  set up the image plane;

• 2.

  propagate the geodesics using either the full equations of motion or the approximate Hamiltonian-Jacobi result,

• 3.

  compute the actual number of photons received now that we know the connection between (b, χ) and (θ, ϕ_e).

We trace the photons from the image plate at infinity all the way to the surface of the star by solving the first-order differential Eqs. (31)–(33) or the full geodesic equations backward in time. In our numerical computations, we place the image plane at ~10⁵R_e. For the Hamilton-Jacobi formalism, we make a simple variable substitution $\tilde{x}=1/r$ that helps us by stretching the step when far from the star and shortening it when approaching the star surface. Here an adaptive step size second-order Heun Runge-Kutta integrator is used with the forward Euler method as the predictor and trapezoidal method as the corrector. All photons that travel more than 1.05 R_e away from the star after a U-turn are terminated and considered to have missed the star. The full geodesic equations are solved using the ARCMANCER code (see Pihajoki et al. 2017, and the related equations therein).

Our image plate is defined using a polar coordinate system with a radial coordinate b (i.e., the impact parameter) and an angular coordinate χ. We also employ a non-equidistant grid in both coordinates to accommodate the extra resolution needed around the edges of the star. The radial coordinate b is defined using a Gauss-Laguerre abscissa (i.e., e^−b-weighted), and the angle coordinate is weighted with a simple sinusoidal function so that the resolution is increased around the top and bottom parts where χ = 0 or π, which is near the location of the poles (see Fig. 2). By ray tracing, we then obtain a mapping between the image plane and the surface of the star, defined on a grid. Arbitrary positions in the image plane are obtained by a quadratic interpolation in (b, χ) space.

For pulse profile calculations where only a small part of the star is emitting, we first search a crude location of the spot on the image plane and then impose a fine subgrid around it in order to accurately calculate the flux from this small patch. The subgrid itself is defined either in a polar grid (by constraining minimum and maximum χ and b) or in a Cartesian image grid (by constrainingminimum and maximum $\widehat{x}$ and ŷ) depending on the total area covered in the observer’s sky. To calculate the total flux F(t), we then integrate this small subgrid by using an adaptive multidimensional integration. The algorithm is based on a tensor product of nested Clenshaw-Curtis quadrature rules and is implemented using the CUBATURE package³.

Such a general way to treat the problem of course also has its disadvantages. Ray tracing photons in general is a computationally very expensive problem. For fast calculations, other more approximate ways exist to solve the problem, such as the oblate Schwarzschild method, where the symmetries of the Schwarzschild space-times are extensively used and the ray tracing reduces to lensing angle integrals (see, e.g., Poutanen & Beloborodov 2006; Morsink et al. 2007). We emphasize that our focus is not to compete with these methods in speed, but to verify their results using a more general description of the problem.

Fig. 2 Example of a non-equidistant polar grid used in our ray tracing with Nr = 20 and Nχ= 30 points. Thered dashed line corresponds to the outline of the actual oblate star that is covered with a chessboard pattern.

3 Applications and verification

Next we present some applications of the framework to some simple physical problems related to neutron stars to showcase possible applications of the code. The examples are also meant to act as a further verification, as we provide a comparison with existing literature results, when possible. Most notably, we perform an extensive comparison against AMP pulse profiles computed with a forward-in-time method as presented in Poutanen & Beloborodov (2006) and AlGendy & Morsink (2014). Our code serves as a great cross-verification tool for these types of special relativistic formulations because our framework is fully general relativistic and propagates photons backward in time from the observer to the surface.

Fig. 3 Formation of an image in curved space-time. Left panel: 3D visualization of the photon trajectories in the curved space-time starting from the star and ending at the observer’s image plane. Right panel: image that the observer sees using the Cartesian $\widehat{x}$ and ŷ coordinates. For illustrative purposes, the neutron star surface is covered with a chessboard pattern.

3.1 Images of neutron stars

As a first application of the code, we can determine the photon trajectories using the ray tracing algorithm and produce an image of the neutron star as as seen by the observer. This also shows how we can connect the Cartesian coordinates $\widehat{x}$ and ŷ of the observer to the coordinates ϕ_e and θ_e of the star. The left panel in Fig. 3 shows the trajectory of the photons in 3D, using a ẑ -coordinate in addition to the Cartesian image plane coordinates $\widehat{x}$ and ŷ. The star is chosen to have R = 12 km, M = 1.5 M_⊙, and has a spherical shape, whereas the observer is located at the equator with an inclination of i = 90°. Here the photons originate from the image plane located at ẑ = 20 and are then propagated backward in time until they intersect with the surface of the star (center of the star located at ẑ = 0), visualized with a spherical see-through wire-grid frame. The right panel of the figure shows the projected image as seen by a distant observer. Here the star is covered with a chessboard pattern to show how the ϕ_e and θ_e coordinates on the neutron star surface are seen by the distant observer. In case of no rotation (f = 0 Hz), the image outline is verified to be mirror symmetric with respect to reflection along the $\widehat{x}=0$ vertical axis and along the ŷ = 0 horizontal axis of the image.

3.2 Accuracy of the split Hamilton-Jacobi propagator

Next, we study the feasibility of the split Hamilton-Jacobi propagator by comparing to results from the general-purpose geodesic solver ARCMANCER (Pihajoki et al. 2017). The comparison solver directly solves the Lagrangian equations of motion of the geodesic in an arbitrary user-given metric. We compute the change in the Jacobi constant as computed with both sides of Eq. (20), and the change in the value of the Hamiltonian of the geodesic for two different neutron star configurations: R_e = 12 km, M = 1.6 M_⊙, and ν = 400 Hz, and R_e = 15 km, M = 1.4 M_⊙, and ν = 600 Hz. The observer inclinations are i = 15°, 45°, and 75°. The error in the Hamiltonian or in the Jacobi constant reflects the error in the photon path. To study the effect this has on the actual observables, we compare the values of the photon redshifts z obtained either with the full numerical propagator or with the split Hamilton-Jacobi propagator. Deviations in this value as a function of the location then reflect the error not only in z itself, but also in the ϕ_e and θ_e coordinates. The results are shown in Fig. 4.

In general, we see that the assumption of separability, as measured by the variation in the Jacobi constant, is good to a level of 10⁻³–10⁻², except for geodesics that hit the center of the neutron star from the observer’s point of view. However, the examples were deliberately chosen to be extreme, and the approximation of separability is much better for more slowly rotating neutron stars. Furthermore, the relative error in the observed redshift is always smaller than 7 × 10⁻³, even for these extreme cases.

This small error has two reasons. First, the splitting of the Hamilton-Jacobi equation is an excellent approximation because the quadrupole moment produces a deviation in the metric only very close to the star. When photons are propagated from a distant location to the stellar surface, the effect on the trajectory is negligible. Second, the splitting does not affect the redshift calculations as we use the exact form of 1 + z as given by Eq. (42). When these two aspects of the method are combined, we observe the excellent agreement against the results obtained from the full geodesic equations. Last, we note that in cases of extreme rotation or when high precision is required, the geodesic propagation can easily be made using the ARCMANCER instead of Eqs. (31)–(33), while using the rest of the results in this paper for the actual radiation computations.

3.3 Line profiles

Next, we study the energy-dependence of the stellar flux by computing the observed energy distribution F(E) of photons emitted from the stellar surface at a single energy E′ as measuredin the comoving frame. In order to minimize any source of error, we use ARCMANCER in this section to solve the geodesics in all of the subsequent calculations. Because of the variation of the redshift across the surface of the star that is caused by Doppler boosting and because of the oblate shape of the star, the observer sees a range of energies (Özel & Psaltis 2003; Bhattacharyya et al. 2006; Chang et al. 2006). Following Bauböck et al. (2015), we are interested in constraining the convolution (smearing) kernel $\mathcal{G}(E,E^{\prime })$ defined via $$F(E)={\displaystyle \int I^{\prime }}(E^{\prime })\mathcal{G}(E,E^{\prime })\text{d}{E}^{\prime },$$(71)

where we have dropped the time and angle dependency of the specific intensity I′ and have explicitly written all quantities to be functions of the energy. It follows from Eq. (63) that the actual flux of photons with an observed energy E and emitted energy E′ is then $$\begin{array}{ll}F(E)\hfill & ={\displaystyle \iint \frac{I^{\prime }(E^{\prime })}{{(1+z)}^3}}\frac{b\text{d}b\text{d}\chi }{D^2}\hfill \\ \hfill & ={\displaystyle \iiint \frac{I^{\prime }(E^{\prime })}{{(1+z)}^4}} \delta \left(E-\frac{E^{\prime }}{1+z}\right)\frac{b\text{d}b\text{d}\chi }{D^2}\text{d}{E}^{\prime },\hfill \end{array}$$(72)

where δ(x) is the Dirac delta function. The convolution kernel, or the so-called line profile, we are after is then $$\mathcal{G}(E,E^{\prime })={\displaystyle \iint \frac{1}{{(1+z)}^4}} \delta \left[E-\frac{E^{\prime }}{1+z}\right]\frac{b\text{d}b\text{d}\chi }{D^2}.$$(73)

Examples of line profiles are shown in Fig. 5 for different space-times and star configurations. The flux is normalized so that the emitted bolometric flux is unity, that is, the area encapsulated by the profile is one. In each case, the star is taken to have ν = 700 Hz, R_e = 10 km, and M = 1.4 M_⊙. The observer inclination is i = 20° and emission is taken to be isotropic, for simplicity. This figure shows line profiles for spherical and oblate stars, assuming for simplicity that the exterior space-time is the Schwarzschild space-time, as well as results for oblate stars in the appropriate second-order exterior space-time, that is, a space-time with the appropriate quadrupole moment given by relations Eq. (7)–(9).

The line profiles computed using the Schwarzschild metric with a spherical star appear to be smooth and asymmetrical with an enhancement toward higher energies caused by the relativistic Doppler boosting (see, e.g., Özel & Psaltis 2003). For an oblate star, the increased redshift of the regions near the pole shifts the peak toward lower energies. The resulting line profile is fairly symmetric (see, e.g., Bauböck et al. 2013b). However, when a physically more realistic metric with a non-zero quadrupole moment is used, the high-energy part of the line profile is further enhanced. This again results in an asymmetrical line profile. When the value of the quadrupole moment is increased to unphysically high levels, the line profile develops a narrow peak in the high-energy part. This effect highly depends on the observer inclination relative to the rotation axis of the star, however.

We now study the line profile shape in full detail using the BENDER code. In order to fully map the change in the line profile shape as a function of observer inclination, we calculated different profiles for three different cases: M = 1.1 M_⊙, 1.5 M_⊙, and 1.8 M_⊙. Here we consider only rapidly spinning stars and hence set the spin to 600 Hz, which is close to the maximum value observed for AMPs. For each mass, the equatorial radius R_e and observer inclination i were taken to span the full range from 10 to 16 km, and 0 to 90°, respectively. Examples of the observed line profiles are shown in Fig. 6 for i = 5°, 10°, 20°, 40°, 60°, and 90°. From here it is easy to see that the profile appears to be smooth at almost all observer viewing angles. Only at i ≲ 5°, a sharp spike starts to form. In this case, however, the actual observed width of the profile is already below 0.03 × E, whereas the spike is as narrow as 0.01 × E. For a spectral energy feature at around 10 keV, therefore, a resolution of 0.1 keV would be needed to resolve it.

We can also try and quantify the observed effect more thoroughly by introducing the full width at half-maximum (FWHM) of the profile (i.e., the width of the profile at F_max∕2). In addition, we consider the full width at tenth-maximum (FWTM) that reflects the total width of the profile (i.e., the width of the profile at F_max∕10). These values are shown for different radii and observer inclinations in Fig. 7. They are also a useful measure of how the rotation would smear the observed spectra: the FWTM gives a quantitative estimate of how widely smeared any narrow feature, such as a line or an absorption edge, would be observed. The FWHM, on the other hand, quantifies the energy resolution needed to resolve the exact effects from rotation itself. Finally, we can also use their ratio to describe the shape of the line profile: the narrower (and hence localized) the line profile feature, the smaller this fraction. For a narrow peak we expect an FWHM/FWTM of around ~ 10 %. This ratio is shown in the bottom panels of Fig. 7. For a star rotating at 600 Hz, a narrow line feature is visible only for observers with inclinations in a very narrow range, that is, 3° < i < 6°, regardless of the mass or radius of the star.

These results can be compared to the results reported in Bauböck et al. (2013b). Here the line profiles using the Schwarzschild exterior metric are seen to match our calculations well. For the profiles computed using a metric that includes corrections up to second order in $\widehat{\Omega }$, we see a clear deviation. Most notably, the line profiles we compute only contain narrow features in a very restricted range of observer inclinations 3° < i < 6°. However, Bauböck et al. (2013b) found narrow spectral features with observer inclinations i ≲ 30° with similar neutron star parameters. A possible reason for this discrepancy can be traced back to how the value of the quadrupole moment is computed. The value of our quadrupole moment is derived from the scaling relations, whereas Bauböck et al. (2013b) set their value of q by hand. They used the Hartle-Thorne metric (Hartle & Thorne 1968) in their calculations, where the quadrupole moment is given by q_inv = −j²(1 + η), where j is the dimensionless spin parameter (see Eq. (5); a in the notation of Bauböck et al. 2013b), and η is the strength of the deviation from a spherical potential (η = 0 reduces to Kerr-like space-time). Bauböck et al. (2013b) selected η = 3.3, which is a typical value given by the FPS EoS (Lorenz et al. 1993) for a star with M ≈ 1.4 M_⊙ (see Laarakkers & Poisson 1999). For the angular momentum they adopted j = 0.357, again as given by the FPS EoS at ν ≈ 700 Hz. With this value, their quadrupole moment is then q_inv ≈−0.548. The radius they imposed was R = 10 km.

However, we note that by selecting an individual EoS and setting the star’s mass and angular momentum, the radius of the star is already determined for physically realistic parameter combinations. In their case, the FPS EoS would yield a considerably different radius of R_e ≈ 11.8 km (Cook et al. 1994; Laarakkers & Poisson 1999). Moreover, Bauböck et al. (2013b) did not include a contribution from the pressure quadrupole moment β in their coordinate-invariant quadrupole expression (Pappas & Apostolatos 2012). On the other hand, we obtain for M = 1.4 M_⊙, R_e = 10 km, and ν = 700 Hz the following values: j ≈ 0.275, q ≈ −0.268, and β ≈ 0.010, which then give q_inv ≈−0.255. Hence, the value of the q_inv used in Bauböck et al. (2013b) is approximately twice that of a physically realistic neutron star with R_e = 10 km.

In general, the quadrupole moment is larger for a stiffer EoS, because a stiffer EoS produces a larger star and the quadrupole moment scales with the square of the radius. For us, this scaling is taken into account by relations Eqs. (7) and (8), which are obtained by fitting a large library of EoSs (see Bauböck et al. 2013a; AlGendy & Morsink 2014) to yield a consistent quadrupole moment at any given mass, radius, and angular velocity. This scaling also hides the difference between the non-rotating and rotating radii because it is formulated using the equatorial radius R_e. Alternatively, the corresponding R₀ of a non-rotating configuration might be considered, for which R₀ ≤ R_e for any given $\widehat{\Omega }$. This distinction between rotating and non-rotating radii is important as EoS modeling for the cold dense matter inside neutron stars is typically done assuming non-rotating radii. In this particular comparison, our j and q_inv are therefore smaller because we require that the radius be 10 km. The line profile emerging from a such a star is shown with a red solid line in Fig. 5and is not seen to develop a narrow core. For an oblate star, we need to artificially increase q by a factor of 4, so that q = −1.07, in order for the line profile to produce a narrow core, as seen in the red dashed line in Fig. 5. In conclusion, we are only able to reproduce the narrow peak with a large observer inclination of 20°, shown in Fig. 2 of Bauböck et al. (2013b), by artificially increasing the value of q.

Bauböck et al. (2013a) subsequently revised their calculations and recomputed their observed line profile in Hartle-Thorne metric with values of the quadrupole moment q originating from a similar physically consistent empirical parameterization. In this case, Bauböck et al. (2013a) still found a narrow spectral feature in the line profile for a neutron star with R_e = 10 km, M = 1.4 M_⊙, and ν = 700 Hz, similar to the parameters used in Bauböck et al. (2013b; see Fig. 5, Bauböck et al. 2013a). However, the observer inclination was not specified. Based on the results we present in Figs. 6 and 7, however we can say that the formation of a narrow peak for these neutron star parameters is only possible in a very limited range of observer inclinations, of around i ~ 5°.

Fig. 4 Errors in ray tracing two neutron stars with Re = 12 km, M = 1.6 M⊙, ν = 400 Hz or Re = 15 km, M = 1.4 M⊙, ν = 600 Hz, and with observer inclinations i = 15°, 45° and 75°, solving the full geodesic equation vs. the split Hamilton-Jacobi equation. The leftmost panels show the maximum variation in the Hamiltonian H of the geodesic, while the two center panels show the maximum variation in the Jacobi constant, computed either with the left side (C1) or the right side (C2) of Eq. (20). Rightmost panels: relative error in the redshift computed by solving the split Hamilton-Jacobi equation compared with the redshift computed by solving the full geodesic equation.

Fig. 5 Line profiles from a star with R = 10 km, M = 1.4 M⊙, and a rotation frequency of 700 Hz seen by an observer at an inclination i = 20° computed by solving the geodesic equations using ARCMANCER. The black line shows the profile of a spherical star with a Schwarzschild exterior space-time. The blue line represents the profile of an oblate star with Re = 10 km and a Schwarzschild exterior space-time. The red solid line denotes the profile of an oblate star with an exterior space-time that has a non-zero quadrupole moment (see text). The red dashed line shows the profile of an oblate star with a quadrupole moment that has been artificially increased by a factor of 4 (see text).

Fig. 6 Exact shapes of line profiles for different neutron stars spinning at 600 Hz. Observer inclinations span a range from i = 5° (red), 10° (blue), 20° (green), 40° (purple), and 60° (orange), to 90° (black). The energy in the horizontal axis is scaled with the compactness $\sqrt{1-u}=\sqrt{1-2GM/R_{\text{e}}}$ that would be expected from the gravitational redshift alone. Likewise, the flux in the vertical axis is always normalized with the peak flux to better show the evolution of the profile shape.

3.4 AMP pulse profiles

From here on, we move to time-dependent ray tracing problems by considering pulse profiles from AMPs. Here a hot spot on the stellar surface is emitting, and the star is seen to rotate with a frequency of Ω. The internal accuracy of the calculations is only set by the error tolerance of the numerical integration of the flux. Hereafter we use a relative tolerance of 5 × 10⁻³. The results here were obtained using both the split Hamilton-Jacobi propagator and ARCMANCER, and they match given the numerical error tolerance. For simplicity, we only show the results obtained by the split Hamilton-Jacobi method in the following discussion.

The definition of the differential surface element is given by Eq. (56) and hence correctly includes the γ_L factor. This comparison is crucial in order to verify that all of the physics is correctly incorporated in the formulations of the given methods. Results between the two methods are therefore expected to agree up to the numerical tolerance.

First, a general comparison of the ray tracing algorithm with the S+D approximation was made using the Schwarzschild metric. For simplicity, only spherical stars were considered here. The main parameters were the stellar mass M = 1.6 M_⊙, the stellar radius R = 12 km, the observer inclination i = 60°, and the colatitude of the spot θ_s = 50°. The effective radiation temperature was set to T_eff = 2 keV. The distance to the source was assumed to be D = 10 kpc. We defined a circular spotwith an angular radius ρ of either 1 or 30 degrees. Here the spot size is defined using its angular size in the corotating coordinate system. The angular distribution of the radiation corresponds either to an isotropic blackbody with constant intensity or to an atmosphere dominated by electron scattering, that is, the Hopf profile.

The light curves are computed in 128 time bins. Zero time t = 0 corresponds to the moment when the spot center crosses the plane defined by the spin axis and the direction to the observer. We computed curves for the five following quantities: monochromatic photon flux (ph cm⁻² s⁻¹ keV⁻¹) at the observer energies E = 2,  6, and 12 keV, bolometric photon flux (ph cm⁻² s⁻¹), and the bolometric energy flux (erg cm⁻² s⁻¹).

The comparison of these light curves is shown in Fig. 8 for a slowly rotating star (1 Hz) and in Fig. 9 for a fast-rotating star (400 Hz). In practice, comparing the profiles for slow rotation only tests our ray tracing routines because special relativistic effects (Doppler boosting, angle aberration, and so on) are negligible. The overall agreement of the two different methods is excellent, and from here, a baseline accuracy of about <0.2% relative error is obtained for the mapping of quantities between image plane and stellar surface. No large deviation between isotropic and Hopf profile is detected either, indicating a good agreement in our emission angle computations and formulation. Similarly, when rotation is increased and special relativistic effects start to play a role, we are usually able to reproduce the pulse profiles down to < 0.3% relative error, except near ϕ_e ~ 0. Here the tilt of the spot increases, and even though the absolute error remains the same, the relative error grows because the observed flux is increasingly lower for a more inclined spot. This situation is numerically expensive when integrating the observed flux from the NS image. In this case, we set set a bound on the number of flux integrand evaluations (typically ~ 10⁷ function calls) that in effect set an absolute error for the flux. It is then only in these rare cases that our integrator does not respect the relative error criteria set by us.

Next we compare emission from oblate stars. The surface here is defined using the radius function Eq. (12), but the star is still embedded in a symmetric Schwarzschild space-time. The parameters we used are an equatorial radius R_e = 12 km (in the usual Schwarzschild metric), a neutron star mass of M = 1.4 M_⊙, an extreme rotational frequency ν = 700 Hz, an observer inclination i = 45°, and a spot angular size of ρ = 10°. The effective temperature of the radiation was again taken to be T_eff = 2 keV, and the distance to be D = 10 kpc. Here the spot size is defined in a corotating spherical coordinate system on top of a unit sphere and is then projected onto the oblate inclined surface. To trace the effects of the changing surface, we considered the spot in three different locations at colatitudes of θ_s = 18°, 45°, and 90°. Additionally, we considered a spot with angle-dependent emission intensity. This was done using the electron-scattering atmosphere with θ_s = 45°.

A comparison of the oblate light curves is shown in Fig. 10. Again we obtain an excellent agreement with the forward-in-time method, with similar errors as in the spherical case (relative error < 0.2%). This agreement is of course expected since our method is general and does not depend on the shape of the emitting surface. The only large deviation is again seen when the spot is viewed from an extreme angle for θ_s = 90°, just before the occultation. We therefore conclude that the two methods, forward-in-time and backward-in-time, agree all the way up to the numerical tolerance.

After the comparisons, we calculate as a last example full skymaps of the emerging radiation from rapidly rotating oblate AMPs, as shown in Fig. 11. The emission from the AMP is shown for all possible observers and is mapped to the vertical axis using the observer’s inclination angle i. The horizontal axis of the map is the usual pulse phase. The brightness of the skymap is proportional to the received bolometric photon number flux. Taking a slice of the skymap at one particular value of i produces the light curve as seen by the observer at that inclination. The calculations here were made for an extreme case of an NS with R_e = 15 km, M = 1.6 M_⊙, and ν = 600 Hz. We considered a spot size of ρ = 10°, with varying colatitudes ranging from near the pole at θ_s = 10° to the equator at θ_s = 90°. We considered the cases of one spot and two antipodal spots with isotropic beaming and blackbody emission with T_eff = 2 keV. As a result, we can see that full occultations are only observed with one spot. From here it is also easy to see the variation in phase of the flux maxima and minima when the viewing angle of the observer changes. The effect becomes most prominent with two spots located at θ_s ~ 50−70° (second antipodal spot at θ_s,2 ~ 110−130°), and the minimum is seen to change from around phase of 0.2 all the way to 0.4.

We also show the corresponding strength of the observed pulsations for each inclination and spot colatitude combinations in Fig. 12. Here the color of the image corresponds to the pulse fraction computed by $$f_p=\frac{F_{\text{max}}-F_{\text{min}}}{F_{\text{max}}+F_{\text{min}}},$$(74)

where F_min and F_max are the minimum and maximum values in the bolometric light curves, respectively. From here the symmetry between θ_s and i becomes obvious as the lines of constant amplitude appear almost symmetric against switching between x and y axis.

Fig. 7 Line profile full width at tenth-maximum (top row) and full width at half-maximum (middle row) as a function of radius and observer inclination computed for different neutron star configurations of M = 1.1 M⊙, 1.5 M⊙, and 1.8 M⊙ spinning at 600 Hz. Additionally, the bottom row shows the ratio of these two, which can be used to quantify the width of the spiky part of the line profile. Note the different inclination scale on the bottom row that is used to highlight the region i < 10°, where the profile evolves from a smooth to a spiky shape.

Fig. 8 Light-curve comparisons for Schwarzschild space-time with a slowly rotating spherical star (R = 12 km, M = 1.6 M⊙, ν = 1 Hz, i = 60°, θs = 50°, ρ = 1°, and Teff = 2 keV) emitting according to a blackbody or Hopf profile with a spot size of either 1 or 30 degrees. The black solid line shows the pulse profiles computed using the S+D approximation (forward-in-time method; see text), and the red dashed line is a profile computed with the code presented here. Lower panel: residuals as Δ = (modelS+D∕model − 1) × 100%.

Fig. 9 Light-curve comparisons for Schwarzschild space-time with a rapidly rotating spherical star (ν = 400 Hz). The other parameters and symbols are the same as in Fig. 8.

Fig. 10 Light-curve comparisons for oblate Schwarzschild space-times with three different spot colatitudes: θs = 18°, 45°, and 90°. Additionally, the bottom row shows the comparison of the pulse profile for an electron-scattering atmosphere for θs = 45°. The parameters of the star, the spot, and the observer are Re = 12.0 km, M = 1.4 M⊙, ν = 700 Hz, i = 45°, and ρ = 10°. The other parameters and symbols are the same as in Fig. 8.

Fig. 11 Skymaps of the emitted radiation as produced by a rapidly rotating oblate AMP with one or two antipodal spots. The star is taken to have Re = 15 km, M = 1.6 M⊙, and ν = 600 Hz. The spot sizes are ρ = 10°, and the emission is coming from an isotropic blackbody with Teff = 2 keV. The red curves enclose the area where occultation is observed.

Fig. 12 Pulse fractions observed from a rapidly rotating oblate AMP. The parameters of the star and the spot(s) are the same as in Fig. 11.

4 Summary

We have presented a detailed study of radiation emerging from and near rotating compact objects. A framework of formulae for solving this problem was derived in a fully general relativistic manner. The formulae were then specialized to the context of rotating neutron stars.

First, we gave a detailed description of the second order in rotation space-time metric in Sect. 2.1. The space-time we used has a non-zero coordinate-invariant mass quadrupole moment q_inv. The components q and β of q_inv are defined via approximate relations for a wide span of neutron star masses, radii, and spins, following AlGendy & Morsink (2014). When the rotation increases, the star also starts to deviate from a sphere because the gravitational force weakens on the equator because the centrifugal force increases. An approximate relation for the resulting oblate spheroidal shape of the star was again obtained (Morsink et al. 2007; AlGendy & Morsink 2014) and was implemented in Sect. 2.2 for an easy but consistent description of the surface.

Second, we derived a new approximate ray tracing approach using the so-called split Hamilton-Jacobi method (also known as super-Hamiltonian method). This derivation was presented and discussed in detail in Sects. 2.3 and 2.4. Insteadof using the geodesic equation that is a second-order differential equation, we separated the Hamilton-Jacobi equation using a third constant of motion known as Carter’s constant. The method is exact up to first order in rotation (Kerr-like space-time with frame-dragging effects) but remains sufficiently accurate also for second order in rotation because deviations caused by the quadrupole moments are small. Formulating the components of the four-momentum vector like this has the useful feature that the polarization of the radiation can easily be taken into account (see, e.g., Chandrasekhar 1998; Dexter 2016)

Third, we gave a thorough description of the calculations related to the actual emission of the radiation. Effects such as redshift, Doppler boosting, and emission angle of the photon were discussed in a fully general relativistic manner in Sects. 2.5 and 2.7. In the special relativistic formulation (see, e.g., Poutanen & Beloborodov 2006), the calculations were made in a flat space-time and were then Lorentz-boosted to the rotating relativistic frame. In Sect. 2.6 we presented a derivation of the solid-angle element that we defined using a rotating coordinate system. The purpose of this was to clarify some common misunderstandings in the literature of how this transformation from the corotating to the static coordinate system can be achieved. We also briefly discussed in Sect. 2.8 the actual intensity of the emerging radiation and presented an iterative method to solve the Chandrasekhar-Ambartsumian integral, along a new approximate polynomial expansion that is related to the angle-dependent electron scattering atmosphere. We then described the actual intensity of the emerging radiation and used as a simple model the angle-dependent electron-scattering atmosphere presented in Sect. 2.8. Numerical methods for solving all of the presented equations were then laid out in Sect. 2.9.

Finally, in Sect. 3 we presented some applications of the framework. A connection to previous work in the literature was also made, when possible. As a first simple example, we showed how the image of an NS is formed in curved space-time. Next, we studied the energy-dependent emission by considering the emerging line profiles. Most notably, we concluded that when a consistent formulation is used to describe how the increasing eccentricity of the star is coupled to the related quadrupole moments, the resulting line profiles develop a narrow core only at an observer inclination of i ~ 5°. Otherwise, the smearing kernels are smooth functions. The effect of the rotational smearing on the observed energy spectra can be estimated using the FWHM and FWTM of the kernels, which we computed for all observer viewing angles. We then studied AMP pulse profiles extensively and thoroughly, and we compared our results to results obtained using existing special relativistic methods found in the literature. Here the agreement between the two methods was found to be excellent when the correct differential surface area element presented in Sect. 2.6 was used. Last, we computed full skymaps of the radiation that emerges from rapidly rotating AMPs, taking into account the oblate shape of the star and the quadrupole moments of the space-time metric.

Acknowledgments

We sincerely thank various people that helped us to improve the formulation of the theory and the actual paper by engaging the authors in helpful discussions: Most notably we would like to thank Juri Poutanen, Fred Lamb, Jean in ’t Zand, Cole Miller, and Sharon Morsink for their help and comments. We also offer special thanks to Tuomo Salmi for producing the pulse profiles for the code comparison. JN acknowledges support from the University of Turku Graduate School in Physical and Chemical Sciences. PP acknowledges support from the Academy of Finland, grant no. 1274931. This work benefited from discussions at the International Symposium on Neutron Stars in the Multi-Messenger Era in May 2016, supported by the National Science Foundation under Grant No. PHY-1430152 (JINA Center for the Evolution of the Elements).

References

All Tables

All Figures

Fig. 1 Geometry of the system. Here we show the underlying spherical coordinate system with ϕ and θ coordinates along with the Cartesian (x, y, z) coordinate system. In addition, the observer’s $\widehat{x}-\widehat{y}$ image plane is shown as viewed from an inclination angle i. The star is taken to rotate rapidly around the y -axis, which leads to an oblate (squeezed) shape for the emitting surface, and the radial vector r and the surface normal n therefore start to differ from each other.

In the text

	Fig. 2 Example of a non-equidistant polar grid used in our ray tracing with Nr = 20 and Nχ= 30 points. Thered dashed line corresponds to the outline of the actual oblate star that is covered with a chessboard pattern.
In the text

	Fig. 3 Formation of an image in curved space-time. Left panel: 3D visualization of the photon trajectories in the curved space-time starting from the star and ending at the observer’s image plane. Right panel: image that the observer sees using the Cartesian $\widehat{x}$ and ŷ coordinates. For illustrative purposes, the neutron star surface is covered with a chessboard pattern.
In the text

Fig. 4 Errors in ray tracing two neutron stars with Re = 12 km, M = 1.6 M⊙, ν = 400 Hz or Re = 15 km, M = 1.4 M⊙, ν = 600 Hz, and with observer inclinations i = 15°, 45° and 75°, solving the full geodesic equation vs. the split Hamilton-Jacobi equation. The leftmost panels show the maximum variation in the Hamiltonian H of the geodesic, while the two center panels show the maximum variation in the Jacobi constant, computed either with the left side (C1) or the right side (C2) of Eq. (20). Rightmost panels: relative error in the redshift computed by solving the split Hamilton-Jacobi equation compared with the redshift computed by solving the full geodesic equation.

In the text

Fig. 5 Line profiles from a star with R = 10 km, M = 1.4 M⊙, and a rotation frequency of 700 Hz seen by an observer at an inclination i = 20° computed by solving the geodesic equations using ARCMANCER. The black line shows the profile of a spherical star with a Schwarzschild exterior space-time. The blue line represents the profile of an oblate star with Re = 10 km and a Schwarzschild exterior space-time. The red solid line denotes the profile of an oblate star with an exterior space-time that has a non-zero quadrupole moment (see text). The red dashed line shows the profile of an oblate star with a quadrupole moment that has been artificially increased by a factor of 4 (see text).

In the text

Fig. 6 Exact shapes of line profiles for different neutron stars spinning at 600 Hz. Observer inclinations span a range from i = 5° (red), 10° (blue), 20° (green), 40° (purple), and 60° (orange), to 90° (black). The energy in the horizontal axis is scaled with the compactness $\sqrt{1-u}=\sqrt{1-2GM/R_{\text{e}}}$ that would be expected from the gravitational redshift alone. Likewise, the flux in the vertical axis is always normalized with the peak flux to better show the evolution of the profile shape.

In the text

Fig. 7 Line profile full width at tenth-maximum (top row) and full width at half-maximum (middle row) as a function of radius and observer inclination computed for different neutron star configurations of M = 1.1 M⊙, 1.5 M⊙, and 1.8 M⊙ spinning at 600 Hz. Additionally, the bottom row shows the ratio of these two, which can be used to quantify the width of the spiky part of the line profile. Note the different inclination scale on the bottom row that is used to highlight the region i < 10°, where the profile evolves from a smooth to a spiky shape.

In the text

Fig. 8 Light-curve comparisons for Schwarzschild space-time with a slowly rotating spherical star (R = 12 km, M = 1.6 M⊙, ν = 1 Hz, i = 60°, θs = 50°, ρ = 1°, and Teff = 2 keV) emitting according to a blackbody or Hopf profile with a spot size of either 1 or 30 degrees. The black solid line shows the pulse profiles computed using the S+D approximation (forward-in-time method; see text), and the red dashed line is a profile computed with the code presented here. Lower panel: residuals as Δ = (modelS+D∕model − 1) × 100%.

In the text

	Fig. 9 Light-curve comparisons for Schwarzschild space-time with a rapidly rotating spherical star (ν = 400 Hz). The other parameters and symbols are the same as in Fig. 8.
In the text

Fig. 10 Light-curve comparisons for oblate Schwarzschild space-times with three different spot colatitudes: θs = 18°, 45°, and 90°. Additionally, the bottom row shows the comparison of the pulse profile for an electron-scattering atmosphere for θs = 45°. The parameters of the star, the spot, and the observer are Re = 12.0 km, M = 1.4 M⊙, ν = 700 Hz, i = 45°, and ρ = 10°. The other parameters and symbols are the same as in Fig. 8.

In the text

	Fig. 11 Skymaps of the emitted radiation as produced by a rapidly rotating oblate AMP with one or two antipodal spots. The star is taken to have Re = 15 km, M = 1.6 M⊙, and ν = 600 Hz. The spot sizes are ρ = 10°, and the emission is coming from an isotropic blackbody with Teff = 2 keV. The red curves enclose the area where occultation is observed.
In the text

	Fig. 12 Pulse fractions observed from a rapidly rotating oblate AMP. The parameters of the star and the spot(s) are the same as in Fig. 11.
In the text