Issue 
A&A
Volume 615, July 2018



Article Number  A50  
Number of page(s)  19  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201630261  
Published online  11 July 2018 
Radiation from rapidly rotating oblate neutron stars
^{1}
Tuorla Observatory, Department of Physics and Astronomy, University of Turku,
Väisäläntie 20,
21500
Piikkiö,
Finland
^{2}
Nordita, KTH Royal Institute of Technology and Stockholm University,
Roslagstullsbacken 23,
10691
Stockholm,
Sweden
email: joonas.nattila@su.se
^{3}
University of Helsinki, Department of Physics,
Gustaf Hällströmin katu 2a,
00560
Helsinki,
Finland
Received:
15
December
2016
Accepted:
5
February
2018
A theoretical framework for emission originating from rapidly rotating oblate compact objects is described in detail. Using a HamiltonJacobi formalism, we show that special relativistic rotational effects such as aberration of angles, Doppler boosting, and time dilatation naturally emerge from the general relativistic treatment of rotating compact objects. We use the Butterworth–Ipser metric expanded up to the second order in rotation and hence include effects of light bending, framedragging, and quadrupole deviations on our geodesic calculations. We also give detailed descriptions of the numerical algorithms used and provide an opensource implementation of the numerical framework called BENDER. As an application, we study spectral line profiles (i.e., smearing kernels) from rapidly rotating oblate neutron stars. We find that in this metric description, the secondorder quadrupole effects are not strong enough to produce narrow observable features in the spectral energy distribution for almost any physically realistic parameter combination, and hence, actually detecting them is unlikely. The full width at tenthmaximum and full width at halfmaximum of the rotation smearing kernels are also reported for all viewing angles. These can then be used to quantitatively estimate the effects of rotational smearing on the observed spectra. We also calculate accurate pulse profiles and observer skymaps of emission from hot spots on rapidly rotating accreting millisecond pulsars. These allow us to quantify the strength of the pulse fractions one expects to observe from typical fastspinning millisecond pulsars.
Key words: methods: numerical / radiative transfer / stars: neutron / gravitation
© 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 spacetime 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^{1}, 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 spacetime 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 spaceborne Xray 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 nonrotating curved spacetime 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). Spacetimes 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 nextorder correction to a nonspherical 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 spacetime 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 HartleThorne metricvariant 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 Xray emission from accreting millisecond pulsars (AMPs; Wijnands & van der Klis 1998; Patruno & Watts 2012) and nuclearpowered 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 nuclearpowered millisecond pulsars, quasicoherent oscillations are observed during a thermonuclear type I Xray 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 secondorder corrections to the spacetime 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 Spacetime 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^{2} = 1, hence describing lengths in units of gravitational radius, where M is the mass of the compact object.
The exterior spacetime of a static, nonrotating, spherically symmetric mass is described by the wellknown Schwarzschild metric (1)
where r is the radial coordinate defined so that the area of a sphere at coordinate time t is 4πr^{2}, 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) (2)
where is the socalled isotropic radial coordinate, and we set . 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 (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 (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) (4)
where is the angular velocity of the local inertial frame, and the functions , and in the metric coefficients can be expanded in the powers of and ū (Butterworth & Ipser 1976). Here is the timedilation factor relating the proper time of the local observer to the time at infinity. Physical interpretation of follows from the fact that the proper circumference of a circle around the axis of symmetry is . Similarly, the interpretation of follows from the fact that acts as a conformal (angle preserving) factor of the spacetime. We also note that the time and space coordinates are connected in the isotropic metric via the term that also enters both the radial and angular terms. The zerothorder terms () of the series expansions are the familiar Schwarzschild metric coefficients expressed in isotropic coordinates (see Table 1).
The firstorder expansion in rotation () is qualitatively related to Kerr metric. In this case, we introduce an angular velocity term of the local inertial frame that accounts for the framedragging effects. It can be defined as (5)
where the dimensionless quantity and is the star’s angular momentum with moment of inertia I(M, Ω).
The secondorder expansion () corresponds to a similar approximation as the HartleThorne slowrotation spacetime (Hartle & Thorne 1968), which introduces two quadrupole moments into the metric. These secondorder 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 (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 coordinateinvariant quadrupole moment) are connected by what is called (approximative) universal relations. Bauböck et al. (2013a) defined empirical relations for these parameters in the HartleThorne 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 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)
We also note that both q and β are , whereas j is .
It is possible to transform between the Schwarzschild coordinate radius r and the isotropic coordinates using the relation (Friedman et al. 1986) (10)
The relation between the differentials of the two radial coordinates is (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 , 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.
Series expansion terms of the metric coefficients up to .
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 (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) (13)
and . The angle γ, defined as the angle between the radial unit vector r and the surface normal n is given by (15)
Then the normal to the surface can be defined using the radial vector r and the tangential vector θ as (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 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 HamiltonJacobi equation
In general, the motion of light rays in curved spacetime is governed by the secondorder geodesic equation. In this section we present an equivalent theoretical formalism based on HamiltonJacobi (sometimes also known as superHamiltonian) 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 firstorder HamiltonJacobi equations and the secondorder 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 HamiltonJacobi 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 spacetime. Geodesic motion in a spacetime characterized by a metric g_{ij} is governed by the HamiltonJacobi equation (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 spacetime, 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 spacetime, yielding two constants of motion, namely energy E and the zcomponent of the angular momentum, L_{z}. We then seek a solution of Eq. (17) in the form (18)
where δ_{1} is related to the restmass of the particle we study. With the metric function Eq. (4), this becomes (19)
After reorganizing terms and introducing a simplifying notation , we obtain (20)
The individual terms in Eq. (20) only depend on r or only on θ, that is, the dependence on r and θ is separable if , , and only depend on r and μ only depends on θ. This is the case to first order in the stellar rotation rate because , in addition to and (see Table 1). However, to secondorder in , 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 higherorder 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 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 fourmomentum p as
Similarly, the components of a local tetrad frame are
where with index and ϕ are the tetrads of metric Eq. (4). Since we only consider null geodesics (i.e., photons), we now set δ_{1} = 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 fourmomentum , and if it is eventually observed at the image plane at infinity, it will have a final fourmomentum of , purely in the radial direction. Likewise, the components of the position must satisfy
as r→∞. The change in the time and angular components along the geodesic can be written as
yielding a total change of angles Δθ and Δϕ when integrating from r_{e} to ∞. The condition for being observed is then
The projected image of the star on the image plane can then be described by two celestial coordinates: abscissa and ordinate ŷ. Making use of the tetrad components Eqs. (25)–(28), we obtain (Chandrasekhar 1998) (36)
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 as (38)
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 . As the final step, we can substitute the fourmomentum components and imageplane coordinates into Eqs. (31)–(33) and solve the system of three firstorder differential equations (in terms of t, θ, and ϕ) with as a variable.
2.5 Redshift and emission angle
It is the most convenient to define all radiative processes in the corotating 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 fourvelocity 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 nonrotating 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 nonrotating 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 nonzero angular velocity in the direction of the star’s rotation as a result of the dragging of inertial frames.
Fourvelocity of a stationary observer with zero angular momentum (socalled ZAMO) is , where the normalization factor is obtained from o_{α}o^{α} = −1, and the velocity of the frame is (40)
indicating that the angular velocity of the ZAMO as measured by an inertial observer at infinity is .
The fourvelocity 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 determined by s_{α}s^{α} = −1. Here the velocity (41)
can be identified as the threevelocity measured in the frame of the ZAMO, an observer rotating with a velocity of v_{ω}.
The total redshift is then given by an innerproduct between a photon u^{α} and a fourvelocity of the star’s fluid s^{α}. With these definitions, the redshift is (42)
consisting of the gravitational part and of the Dopplerlike factor (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 (44)
which projects fourvectors from the nonrotating frame to the rotating frame where the radiative processes are defined. Our metric tensor for the rotating observer is then (45)
As a definition, we can take the emission angle to be the angle between photon and a spacelike surface normal vector , with normalization fulfilling n_{α} n^{α} = −1. The line element in spherical coordinates is , 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 dotproduct definition between two vectors as (46)
Using the metric defined by Eq. (45) and a photon with components Eqs. (21)−(24), we obtain (47)
For the nonrotating observer, we similarly obtain (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, . The result here is identical to the emission angle obtained with a specialrelativistic approach, using flatspace trigonometry and Lorentzboosting 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 backwardintime method to the methods where light rays are propagated from the star to the image plane (forwardintime methods).
Transforming from the observer’s nonrotating frame K to the fluid rest frame K′ is easily done using the previously defined projection operator 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 Lorentzboosted metric tensor (49)
The result agrees with the Schwarzschild coordinate system up to first order in rotation because . This notation can be further simplified by defining a new azimuthal angular coordinate as ϕ′ ≡ ϕ −Ωt that is to be used by the rotating observer^{2}. 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 Lorentzstretched by a factor of . The corotating observer can then use this projection operator to define spacelike 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 (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_{ab}s^{a} = 0. For simplicity, we now consider a spherical star so that the surface normal reduces to , where is the Kronecker delta function. Then m can be expressed as (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 socalled 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 socalled 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 Lorentzboosted into the static nonrotating frame. We now show the correct expression of this change of frame so that the results match the backwardintime 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 (52)
by setting cosγ → 1, sin γ → 0 (spherical star), and ω → 0 along with (Schwarzschild metric, i.e., ). Here we use the emission angle as measured by the nonrotating observer in frame K. In this case, the solid angle is then simplified to (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 (54)
where dS = R^{2} sinθdθdϕ is the area element for the nonrotating 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 finitesized 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 (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 . 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 (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^{2} sinθ′dθ′dϕ′. To transform from the nonrotating K frame to the corotating frame K′, we can use the Lorentz invariance of the photon beam crosssection given as (Terrell 1959; Lind & Blandford 1985) (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 (58)
Finally, the total observed angular size is then seen to be (59)
2.7 Emission
The observed (i.e., received) flux at photon energy E from a small area on an image plane is (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 . The connection between the monochromatic observed and local intensity is then (see, e.g., Misner et al. 1973; Rybicki & Lightman 1979) (61)
where 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 (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, (63)
where t_{*} = t −Δt is the time when the photon was emitted as measured in the nonrotating 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 (nonrotating) 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 (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 planeparallel, semiinfinite (optical depth τ →∞) atmosphere. This beaming pattern is described by the socalled Hopf function H(μ). Introducing a variable μ ≡ cosα′, the result is (65)
are the moments of the function H(μ), which is a solution of the AmbartsumianChandrasekhar integral equation (see, e.g. Chandrasekhar 1960; Sobolev 1963) (67)
Here Ψ(μ) is the characteristic function, which depends on the scattering law considered. For Rayleigh (Thomson) scattering, it is (68)
Given Ψ, the integral Eq. (67) can then be iteratively solved, for example, by computing (69)
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. Electronscattering 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 colorcorrection factor f_{c} as (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} = 1.19167.
2.9 Method of solution
In practice, when calculating the observed timedependent emission, we have to
 1.
set up the image plane;
 2.
propagate the geodesics using either the full equations of motion or the approximate HamiltonianJacobi 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 firstorder differential Eqs. (31)–(33) or the full geodesic equations backward in time. In our numerical computations, we place the image plane at ~10^{5}R_{e}. For the HamiltonJacobi formalism, we make a simple variable substitution 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 secondorder Heun RungeKutta 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 Uturn 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 nonequidistant grid in both coordinates to accommodate the extra resolution needed around the edges of the star. The radial coordinate b is defined using a GaussLaguerre 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 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 ClenshawCurtis quadrature rules and is implemented using the CUBATURE package^{3}.
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 spacetimes 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 nonequidistant polar grid used in our ray tracing with N_{r} = 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 forwardintime method as presented in Poutanen & Beloborodov (2006) and AlGendy & Morsink (2014). Our code serves as a great crossverification 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 spacetime. Left panel: 3D visualization of the photon trajectories in the curved spacetime starting from the star and ending at the observer’s image plane. Right panel: image that the observer sees using the Cartesian 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 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 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 seethrough wiregrid 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 vertical axis and along the ŷ = 0 horizontal axis of the image.
3.2 Accuracy of the split HamiltonJacobi propagator
Next, we study the feasibility of the split HamiltonJacobi propagator by comparing to results from the generalpurpose geodesic solver ARCMANCER (Pihajoki et al. 2017). The comparison solver directly solves the Lagrangian equations of motion of the geodesic in an arbitrary usergiven 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 HamiltonJacobi 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^{−3}–10^{−2}, 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^{−3}, even for these extreme cases.
This small error has two reasons. First, the splitting of the HamiltonJacobi 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 energydependence 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 defined via (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 (72)
where δ(x) is the Dirac delta function. The convolution kernel, or the socalled line profile, we are after is then (73)
Examples of line profiles are shown in Fig. 5 for different spacetimes 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 spacetime is the Schwarzschild spacetime, as well as results for oblate stars in the appropriate secondorder exterior spacetime, that is, a spacetime 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 nonzero quadrupole moment is used, the highenergy 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 highenergy 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 halfmaximum (FWHM) of the profile (i.e., the width of the profile at F_{max}∕2). In addition, we consider the full width at tenthmaximum (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 , 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 HartleThorne metric (Hartle & Thorne 1968) in their calculations, where the quadrupole moment is given by q_{inv} = −j^{2}(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 Kerrlike spacetime). 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 coordinateinvariant 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 nonrotating and rotating radii because it is formulated using the equatorial radius R_{e}. Alternatively, the corresponding R_{0} of a nonrotating configuration might be considered, for which R_{0} ≤ R_{e} for any given . This distinction between rotating and nonrotating radii is important as EoS modeling for the cold dense matter inside neutron stars is typically done assuming nonrotating 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 HartleThorne 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 R_{e} = 12 km, M = 1.6 M_{⊙}, ν = 400 Hz or R_{e} = 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 HamiltonJacobi 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 (C_{1}) or the right side (C_{2}) of Eq. (20). Rightmost panels: relative error in the redshift computed by solving the split HamiltonJacobi 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 spacetime. The blue line represents the profile of an oblate star with R_{e} = 10 km and a Schwarzschild exterior spacetime. The red solid line denotes the profile of an oblate star with an exterior spacetime that has a nonzero 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 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 timedependent 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^{−3}. The results here were obtained using both the split HamiltonJacobi propagator and ARCMANCER, and they match given the numerical error tolerance. For simplicity, we only show the results obtained by the split HamiltonJacobi 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^{−2} s^{−1} keV^{−1}) at the observer energies E = 2, 6, and 12 keV, bolometric photon flux (ph cm^{−2} s^{−1}), and the bolometric energy flux (erg cm^{−2} s^{−1}).
The comparison of these light curves is shown in Fig. 8 for a slowly rotating star (1 Hz) and in Fig. 9 for a fastrotating 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^{7} 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 spacetime. 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 angledependent emission intensity. This was done using the electronscattering atmosphere with θ_{s} = 45°.
A comparison of the oblate light curves is shown in Fig. 10. Again we obtain an excellent agreement with the forwardintime 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, forwardintime and backwardintime, 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 (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 tenthmaximum (top row) and full width at halfmaximum (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 Lightcurve comparisons for Schwarzschild spacetime with a slowly rotating spherical star (R = 12 km, M = 1.6 M_{⊙}, ν = 1 Hz, i = 60°, θ_{s} = 50°, ρ = 1°, and T_{eff} = 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 (forwardintime method; see text), and the red dashed line is a profile computed with the code presented here. Lower panel: residuals as Δ = (model_{S+D}∕model − 1) × 100%. 
Fig. 9 Lightcurve comparisons for Schwarzschild spacetime with a rapidly rotating spherical star (ν = 400 Hz). The other parameters and symbols are the same as in Fig. 8. 
Fig. 10 Lightcurve comparisons for oblate Schwarzschild spacetimes with three different spot colatitudes: θ_{s} = 18°, 45°, and 90°. Additionally, the bottom row shows the comparison of the pulse profile for an electronscattering atmosphere for θ_{s} = 45°. The parameters of the star, the spot, and the observer are R_{e} = 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 R_{e} = 15 km, M = 1.6 M_{⊙}, and ν = 600 Hz. The spot sizes are ρ = 10°, and the emission is coming from an isotropic blackbody with T_{eff} = 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 spacetime metric in Sect. 2.1. The spacetime we used has a nonzero coordinateinvariant 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 socalled split HamiltonJacobi method (also known as superHamiltonian method). This derivation was presented and discussed in detail in Sects. 2.3 and 2.4. Insteadof using the geodesic equation that is a secondorder differential equation, we separated the HamiltonJacobi equation using a third constant of motion known as Carter’s constant. The method is exact up to first order in rotation (Kerrlike spacetime with framedragging 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 fourmomentum 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 spacetime and were then Lorentzboosted to the rotating relativistic frame. In Sect. 2.6 we presented a derivation of the solidangle 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 ChandrasekharAmbartsumian integral, along a new approximate polynomial expansion that is related to the angledependent electron scattering atmosphere. We then described the actual intensity of the emerging radiation and used as a simple model the angledependent electronscattering 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 spacetime. Next, we studied the energydependent 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 spacetime 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 MultiMessenger Era in May 2016, supported by the National Science Foundation under Grant No. PHY1430152 (JINA Center for the Evolution of the Elements).
References
 Agrawal, P. C. 2006, Adv. Space Res., 38, 2989 [NASA ADS] [CrossRef] [Google Scholar]
 AlGendy, M., & Morsink, S. M. 2014, ApJ, 791, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Bardeen, J. M., & Wagoner, R. V. 1971, ApJ, 167, 359 [NASA ADS] [CrossRef] [Google Scholar]
 Bauböck, M., Psaltis, D., Özel, F., & Johannsen, T. 2012, ApJ, 753, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Bauböck, M., Berti, E., Psaltis, D., & Özel, F. 2013a, ApJ, 777, 68 [NASA ADS] [CrossRef] [Google Scholar]
 Bauböck, M., Psaltis, D., & Özel, F. 2013b, ApJ, 766, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Bauböck, M., Özel, F., Psaltis, D., & Morsink, S. M. 2015, ApJ, 799, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Bhattacharyya, S., Miller, M. C., & Lamb, F. K. 2006, ApJ, 644, 1085 [NASA ADS] [CrossRef] [Google Scholar]
 Braje, T. M., Romani, R. W., & Rauch, K. P. 2000, ApJ, 531, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Butterworth, E. M., & Ipser, J. R. 1976, ApJ, 204, 200 [NASA ADS] [CrossRef] [Google Scholar]
 Cadeau, C., Leahy, D. A., & Morsink, S. M. 2005, ApJ, 618, 451 [NASA ADS] [CrossRef] [Google Scholar]
 Cadeau, C., Morsink, S. M., Leahy, D., & Campbell, S. S. 2007, ApJ, 654, 458 [NASA ADS] [CrossRef] [Google Scholar]
 Chan, C.k., Psaltis, D., & Özel, F. 2013, ApJ, 777, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1960, Radiative transfer (New York: Dover) [Google Scholar]
 Chandrasekhar, S. 1998, The Mathematical Theory of Black Holes, The International series of monographs on physics (Oxford: Oxford Univ. Press) [Google Scholar]
 Chang, P., Morsink, S., Bildsten, L., & Wasserman, I. 2006, ApJ, 636, L117 [NASA ADS] [CrossRef] [Google Scholar]
 Cook, G. B., Shapiro, S. L., & Teukolsky, S. A. 1994, ApJ, 424, 823 [NASA ADS] [CrossRef] [Google Scholar]
 Dexter, J. 2016, MNRAS, 462, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Dexter, J., & Agol, E. 2009, ApJ, 696, 1616 [NASA ADS] [CrossRef] [Google Scholar]
 Friedman, J., & Stergioulas, N. 2013, Rotating Relativistic Stars, Cambridge Monographs on Mathematical Physics (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Friedman, J. L., Ipser, J. R., & Parker, L. 1986, ApJ, 304, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Gendreau, K. C., Arzoumanian, Z., & Okajima, T. 2012, Space Telescopes and Instrumentation 2012: Ultraviolet to Gamma Ray, Proc. SPIE, 8443, 844313 [CrossRef] [Google Scholar]
 Hartle, J. B., & Thorne, K. S. 1968, ApJ, 153, 807 [NASA ADS] [CrossRef] [Google Scholar]
 Laarakkers, W. G., & Poisson, E. 1999, ApJ, 512, 282 [NASA ADS] [CrossRef] [Google Scholar]
 Lamb, F. K., Boutloukos, S., Van Wassenhove, S., et al. 2009a, ApJ, 706, 417 [NASA ADS] [CrossRef] [Google Scholar]
 Lamb, F. K., Boutloukos, S., Van Wassenhove, S., et al. 2009b, ApJ, 705, L36 [NASA ADS] [CrossRef] [Google Scholar]
 Lind, K. R., & Blandford, R. D. 1985, ApJ, 295, 358 [NASA ADS] [CrossRef] [Google Scholar]
 Lo, K. H., Miller, M. C., Bhattacharyya, S., & Lamb, F. K. 2013, ApJ, 776, 19 [Google Scholar]
 Lorenz, C. P., Ravenhall, D. G., & Pethick, C. J. 1993, Phys. Rev. Lett., 70, 379 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Miller, M. C., & Lamb, F. K. 1998, ApJ, 499, L37 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, M. C., & Lamb, F. K. 2015, ApJ, 808, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation (San Francisco: W.H. Freeman and Co.) [Google Scholar]
 Morsink, S. M., & Stella, L. 1999, ApJ, 513, 827 [NASA ADS] [CrossRef] [Google Scholar]
 Morsink, S. M., Leahy, D. A., Cadeau, C., & Braga, J. 2007, ApJ, 663, 1244 [NASA ADS] [CrossRef] [Google Scholar]
 Muno, M. P., Chakrabarty, D., Galloway, D. K., & Psaltis, D. 2002, ApJ, 580, 1048 [NASA ADS] [CrossRef] [Google Scholar]
 Narayan, R., Zhu, Y., Psaltis, D., & Saḑowski A. 2016, MNRAS, 457, 608 [NASA ADS] [CrossRef] [Google Scholar]
 Özel, F., & Psaltis, D. 2003, ApJ, 582, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Page, D. 1995, ApJ, 442, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Papitto, A., Torres, D. F., Rea, N., & Tauris, T. M. 2014, A&A, 566, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pappas, G., & Apostolatos, T. A. 2012, Phys. Rev. Lett., 108, 231104 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Patruno, A., & Watts, A. L. 2012, ArXiv eprints [arXiv: 1206.2727] [Google Scholar]
 Pechenick, K. R., Ftaclas, C., & Cohen, J. M. 1983, ApJ, 274, 846 [NASA ADS] [CrossRef] [Google Scholar]
 Pihajoki, P., Rantala, A., & Johansson, P. H. 2017, IAU Symp., 324, 347 [NASA ADS] [Google Scholar]
 Poutanen, J., & Beloborodov, A. M. 2006, MNRAS, 373, 836 [NASA ADS] [CrossRef] [Google Scholar]
 Poutanen, J., & Gierliński, M. 2003, MNRAS, 343, 1301 [NASA ADS] [CrossRef] [Google Scholar]
 Psaltis, D., & Johannsen, T. 2012, ApJ, 745, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Psaltis, D., & Özel, F. 2014, ApJ, 792, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Pu, H.Y., Yun, K., Younsi, Z., & Yoon, S.J. 2016, ApJ, 820, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Rybicki, G. B., & Lightman, A. P. 1979, Radiative processes in astrophysics (New York: WileyInterscience) [Google Scholar]
 Schnittman, J. D., & Krolik, J. H. 2013, ApJ, 777, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Shcherbakov, R. V., & McKinney, J. C. 2013, ApJ, 774, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Sobolev, V. V. 1963, A treatise on radiative transfer (Princeton: Van Nostrand) [Google Scholar]
 Soffitta, P., Barcons, X., Bellazzini, R., et al. 2013, Exp. Astron., 36, 523 [NASA ADS] [CrossRef] [Google Scholar]
 Suleimanov, V., Poutanen, J., & Werner, K. 2011, A&A, 527, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suleimanov, V., Poutanen, J., & Werner, K. 2012, A&A, 545, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Terrell, J. 1959, Phys. Rev., 116, 1041 [NASA ADS] [CrossRef] [Google Scholar]
 Vincent, F. H., Paumard, T., Gourgoulhon, E., & Perrin, G. 2011, Classical and Quantum Gravity, 28, 225011 [NASA ADS] [CrossRef] [Google Scholar]
 Watts, A. L. 2012, ARA&A, 50, 609 [NASA ADS] [CrossRef] [Google Scholar]
 Weinberg, N., Miller, M. C., & Lamb, D. Q. 2001, ApJ, 546, 1098 [NASA ADS] [CrossRef] [Google Scholar]
 Wijnands, R., & van der Klis, M. 1998, Nature, 394, 344 [NASA ADS] [CrossRef] [Google Scholar]
 Yagi, K., & Yunes, N. 2013, Science, 341, 365 [Google Scholar]
 Zhang, S. N., Feroci, M., Santangelo, A., et al. 2016, Proc. SPIE, 9905, 99051Q [Google Scholar]
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 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 nonequidistant polar grid used in our ray tracing with N_{r} = 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 spacetime. Left panel: 3D visualization of the photon trajectories in the curved spacetime starting from the star and ending at the observer’s image plane. Right panel: image that the observer sees using the Cartesian 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 R_{e} = 12 km, M = 1.6 M_{⊙}, ν = 400 Hz or R_{e} = 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 HamiltonJacobi 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 (C_{1}) or the right side (C_{2}) of Eq. (20). Rightmost panels: relative error in the redshift computed by solving the split HamiltonJacobi 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 spacetime. The blue line represents the profile of an oblate star with R_{e} = 10 km and a Schwarzschild exterior spacetime. The red solid line denotes the profile of an oblate star with an exterior spacetime that has a nonzero 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 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 tenthmaximum (top row) and full width at halfmaximum (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 Lightcurve comparisons for Schwarzschild spacetime with a slowly rotating spherical star (R = 12 km, M = 1.6 M_{⊙}, ν = 1 Hz, i = 60°, θ_{s} = 50°, ρ = 1°, and T_{eff} = 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 (forwardintime method; see text), and the red dashed line is a profile computed with the code presented here. Lower panel: residuals as Δ = (model_{S+D}∕model − 1) × 100%. 

In the text 
Fig. 9 Lightcurve comparisons for Schwarzschild spacetime 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 Lightcurve comparisons for oblate Schwarzschild spacetimes with three different spot colatitudes: θ_{s} = 18°, 45°, and 90°. Additionally, the bottom row shows the comparison of the pulse profile for an electronscattering atmosphere for θ_{s} = 45°. The parameters of the star, the spot, and the observer are R_{e} = 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 R_{e} = 15 km, M = 1.6 M_{⊙}, and ν = 600 Hz. The spot sizes are ρ = 10°, and the emission is coming from an isotropic blackbody with T_{eff} = 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.