© The Authors 2026

1 Introduction

Rotation-induced modulations of emission from neutron-star surfaces known as X-ray pulse profiles encode the surface temperature distribution and the global parameters: mass, M, and radius, R. The emission is often dominated by one or more localized, hot regions (i.e., hotspots). Photons from these regions experience gravitational redshift and strong-field light bending, as well as special-relativistic Doppler boosting and time delays. These effects leave signatures in the energy-resolved waveform, encoding the star’s compactness, our primary interest, as well as other global and geometric parameters (Pechenick et al. 1983; Beloborodov 2002; Poutanen & Beloborodov 2006; Morsink et al. 2007; Cadeau et al. 2007; Psaltis & Özel 2014; Watts et al. 2016; Nättilä & Pihajoki 2018; Poutanen 2020). Although this picture is conceptually clean, to model observed pulse profiles accurately, and to extract the hotspot configuration and mass-radius (M-R) relation, we must incorporate the realistic atmospheric beaming pattern, including composition and magnetic-field effects (Zavlin et al. 1996; Ho & Lai 2001, 2003; Suleimanov et al. 2012; Salmi et al. 2020, 2023), rotation-induced oblateness (Miller & Lamb 1998; Poutanen & Gierliński 2003; Morsink et al. 2007; AlGendy & Morsink 2014; Nättilä & Pihajoki 2018; Jakab & Morsink 2025), hotspot geometry (Poutanen & Beloborodov 2006; Viironen & Poutanen 2004; Bogdanov et al. 2008), interstellar absorption (Verner et al. 1996; Wilms et al. 2000), and the detailed instrument response (Gendreau et al. 2016; Prigozhin et al. 2016; Remillard et al. 2022; Markwardt et al. 2024). With these ingredients and practical approximations established, reliable forward modeling procedures can be carried out across the full parameter space, enabling inference of M and R using the X-ray pulse profile modeling (PPM) technique (Bogdanov et al. 2019b, 2021; Miller et al. 2019; Riley et al. 2019; Miller et al. 2021; Riley et al. 2021; Salmi et al. 2024a; Choudhury et al. 2024a; Vinciguerra et al. 2024; Dittmann et al. 2024).

Thermal emission from neutron-star surfaces peaks in the soft X-ray band (∼0.1–2 keV) (Zavlin et al. 1996; Potekhin et al. 2015) and, thus, obtaining precise pulse profile measurements requires high throughput at these energies, sub-microsecond time tagging, and good energy resolution. NASA’s Neutron Star Interior Composition Explorer (NICER) provides this capability: a 0.2–12 keV bandpass, large effective area at soft X-rays, ∼100 ns absolute timing, and ∼10² eV energy resolution near 1 keV (Gendreau et al. 2016; Prigozhin et al. 2016). NICER’s primary pulse profile targets are rotation-powered millisecond pulsars (MSPs), whose highly coherent spin permits phase-folding of long observations to build energy-resolved waveforms with high signal-to-noise ratio. Statistical uncertainties then decrease roughly as the square root of the total photon counts (or exposure time). Accordingly, NICER has concentrated on a small set of bright, well-characterized MSPs that would be suitable for running the PPM (Bogdanov et al. 2019a; Guillot et al. 2019). Looking ahead, the upcoming enhanced X-ray Timing and Polarimetry mission (eXTP), currently planned for launch in the early 2030s, is designed to combine large-area timing, spectroscopy, and polarimetry, and is expected to make PPM an even more powerful probe of neutron-star dense matter (Zhang et al. 2025; Li et al. 2025; Watts et al. 2019).

The first landmark NICER PPM result came from PSR J0030+0451, an isolated millisecond pulsar with no prior radio mass measurement (Manchester et al. 2005; Lommen et al. 2006). Independent analyses of the same NICER dataset produced joint posteriors for the stellar mass and radius and favored nonantipodal emitting-region geometries, establishing the first NICER mass-radius constraints via PPM. Reanalyses have refined these constraints (Miller et al. 2019; Riley et al. 2019; Vinciguerra et al. 2024). NICER’s subsequent measurement for PSR J0740+6620, whose mass is tightly constrained by radio timing at ≃2 M_⊙ (Cromartie et al. 2020; Fonseca et al. 2021), yielded radius measurements at the high-mass end from independent analyses (Miller et al. 2021; Riley et al. 2021), with subsequent updates using additional NICER data (Salmi et al. 2024a; Dittmann et al. 2024) placing strong constraints on the dense-matter equation of state and disfavoring very soft models that predict markedly smaller radii (e.g., Raaijmakers et al. 2021b; Zhang & Li 2021). In both sources, the inferred hotspot locations are inconsistent with a simple centered dipole, suggesting more complex (e.g., offset or multipolar) surface-field configurations (Bilous et al. 2019; Pétri et al. 2025; Kalapotharakos et al. 2021; Chen et al. 2020; Huang & Chen 2025). Recent and ongoing NICER campaigns have targeted bright MSPs such as PSR J0437−4715, whose radio timing yields a precise mass of M = 1.418 ± 0.044 M_⊙, and PSR J0614−3329, thereby extending radius measurements toward the canonical ∼1.4 M_⊙ mass scale (Reardon et al. 2016; Reardon et al. 2024; Choudhury et al. 2024a; Mauviard et al. 2025).

Together with tidal-deformability constraints from gravitational-wave events (Abbott et al. 2017; Abbott et al. 2018; Abbott et al. 2019), these results have energized the nuclear-astrophysics community and motivated increasingly parametric and nonparametric equation-of-state inference frameworks that confront neutron-star structure with multi-messenger data (e.g., De et al. 2018; Landry & Essick 2019; Raaijmakers et al. 2019, 2020; Raaijmakers et al. 2021a; Essick et al. 2020; Greif et al. 2019; Legred et al. 2021; Huth et al. 2022; Annala et al. 2022; Huang et al. 2024a,c,b; Rutherford et al. 2024; Zhou & Huang 2025; Cartaxo et al. 2026; Huang & Sourav 2025; Huang 2025a,b). Since NICER’s initial PPM results, the applications and methodological advancements of PPM have expanded substantially. Open-source frameworks, such as X-PSI, have made fully Bayesian PPM accessible and reproducible, providing flexible hotspot parameterizations (Riley et al. 2023a,b).

Reanalyses of PSR J0030+0451 have shown that the choice of hotspot maps and the treatment of instrumental and background systematics can shift the inferred M–R posteriors, especially when joint NICER+XMM constraints are imposed (Vinciguerra et al. 2024). This approach has since been extended to additional NICER targets, notably PSR J0437−4715 (Choudhury et al. 2024a), and recent studies of PSR J1231−1411 (Salmi et al. 2024b). In parallel, key physical ingredients continue to be refined and several formalisms now propagate polarization predictions suitable for cross-mission tests (Heyl & Shaviv 2000; Viironen & Poutanen 2004; Loktev et al. 2020; Salmi et al. 2025). Methodological progress also includes faster and more accurate relativistic ray-tracing schemes (Nättilä & Pihajoki 2018; Poutanen 2020). Overall, PPM remains a rapidly advancing field.

Despite major advances, several bottlenecks still limit current PPM analyses. First, degeneracies among hotspot geometry and stellar compactness can broaden posterior inferences (Poutanen & Beloborodov 2006; Lo et al. 2013; Psaltis & Özel 2014; Zhao et al. 2024; Güver et al. 2026; Bootsma et al. 2025). One response is to adopt physics-informed hotspot parameterizations that reduce ad hoc geometric degrees of freedom; such approaches have recently been applied to PSR J0030+0451 and PSR J0437−4715 (Huang & Chen 2025). It has also been shown by Lamb et al. (2009a,b); Miller & Lamb (2015) that for NICER-like data, several plausible sources of systematic error do not yield apparently acceptable fits with radius biases exceeding 1 − σ.

Second, PPM is computationally intensive: accurate forward models must couple relativistic ray tracing, realistic atmospheric beaming, rotation-induced oblateness, and instrumental background and calibration systematics. Rigorous Bayesian inference (e.g., nested sampling or ensemble Markov chain Monte Carlo, i.e., MCMC) typically requires ∼10⁶–10⁸ likelihood evaluations in a high-dimensional parameter space (Speagle 2020; Buchner 2021) and incorporating radio constraints or joint NICER+XMM data further increases the cost (Miller et al. 2021; Salmi et al. 2024a; Choudhury et al. 2024a). Third, there is an accuracy-speed trade-off: coarse angular, phase, and energy grids can bias the likelihood and, thus, the posteriors for certain hotspot configurations; whereas the near-converged resolutions needed for robust inference (i.e., to render numerical errors subdominant to Poisson counting fluctuation) make forward modeling computationally demanding and can limit the scale of large-scale Bayesian analyses (Choudhury et al. 2024b; Bogdanov et al. 2019b, 2021; Salmi et al. 2024b; Vinciguerra et al. 2024). Recent cross-code tests and case studies identify extreme geometries where insufficient resolution measurably alters the waveform and likelihood computation, motivating higher-resolution validation runs (Choudhury et al. 2024b). At present, the relatively large observational uncertainties from NICER mean that resolution choices do not yet have a visible effect on shifting the inferred mass-radius. Nonetheless, these differences are expected to become important for future missions such as eXTP, for which PPM is a key dense-matter science driver (Watts et al. 2019; Li et al. 2025). In practice, a single high-fidelity likelihood evaluation can take ∼10, s or more (depending on resolution and model complexity), so many analyses adopt simplified hotspot parameterizations and/or multistage work-flows (coarse-resolution exploration followed by high-resolution verification) to keep compute costs tractable. The principal bottleneck remains computational speed.

In this work, we address these limitations by introducing a GPU-accelerated PPM framework written in C++/CUDA. The framework supports energy- and phase-resolved data. NICER-specific instrument response and background models are implemented as modular components that can be swapped for other observatories. From the outset, the forward model accommodates complex, physics-motivated surface-temperature distributions beyond circular, isothermal caps, while preserving high-fidelity physics. Conceptually, our forward model adopts the standard physical ingredients established in prior CPU frameworks (e.g., Riley et al. 2023b). Our contributions are: (i) a reengineered, GPU-optimized implementation of PPM for contemporary hardware; (ii) end-to-end reproducibility checks against existing benchmarks and codebases and a demonstration that the accuracy-speed trade-off can be broken by performing high-fidelity computation at unprecedented speed; and (iii) identification of a source of systematic error introduced by interpolating atmosphere lookup tables near grid boundaries.

Both the CPU and GPU implementations are released as open source to provide a reproducible foundation for future studies. In validation tests, the framework reproduces published waveform benchmarks to within ≲10⁻³ relative error while delivering 10³–10⁴× speedups over representative CPU baselines. This reduces a single likelihood evaluation from seconds-minutes to milliseconds on a single GPU, making full Bayesian inference, which was previously practical only on large clusters, feasible on personal hardware.

The increased computation speed enables high angular, phase, and energy resolution without a prohibitive cost, mitigating the resolution-induced likelihood biases noted in prior work (Choudhury et al. 2024b). Consequently, a much broader family of fine-detail hotspot maps becomes computationally tractable, including cases previously impractical due to compute budgets. We anticipate that this framework will facilitate more realistic geometry priors, wider parameter exploration, and multi-instrument or multiwavelength joint analyses in forthcoming PPM studies built on this codebase.

The paper is organized as follows. Section 2 summarizes the theoretical framework, highlights subtleties that affect numerical accuracy, presents a transparent computational recipe to facilitate full reproducibility, and details our GPU implementation. Section 3 provides the first set of validation tests, demonstrating that our CPU and GPU implementations reproduce theoretical benchmarks under standard OS assumptions (see Bogdanov et al. 2019b). Section 4 investigates interpolation within atmosphere tables, detailing our interpolation strategy and introducing two new test cases that underscore the need for careful treatment. Section 5 reports production-grade computations that incorporate instrument response, interstellar absorption, and realistic atmospheric beaming (i.e., all the essential ingredients of PPM forward computation) and shows that the framework remains fast and robust even for hotspot configurations previously regarded as impractical, thereby overcoming the long-standing accuracy-speed trade-off. Section 6 concludes with a brief summary and discussion. Appendices A–C supplement the main text with rederivations of some core formulas useful for PPM forward computation.

2 Model framework and computational recipe

In this section, we present the theoretical framework for X-ray PPM forward computation and provide a transparent, step-by-step recipe that covers the theoretical framework and implementation details to enable end-to-end reproducibility. The objective of PPM inference is to infer fundamental neutron-star properties from the energy- and phase-resolved X-ray waveforms. Given the statistical precision of modern data, we aim to limit forward-model inaccuracies to ≲1% in phase-resolved flux where feasible.

The relevant physics combines general relativity (GR) and special relativity (SR). In GR, strong gravity near the neutron star (NS) alters photon trajectories and energies: gravitational light bending increases the visible surface fraction, revealing regions beyond the local limb; gravitational redshift reduces photon energies; and time-of-flight delays in curved space-time introduce phase shifts. Rotation adds SR effects: Doppler boosting modulates the observed flux as a hotspot approaches or recedes from the observer, while rapid spin induces stellar oblateness. These effects are set by the global parameters (M, R, ν_∗) and by the temperature map on the neutron-star surface. Foundational works include Pechenick et al. (1983); Beloborodov (2002); Poutanen & Gierliński (2003); Miller & Lamb (1998); Poutanen & Beloborodov (2006); Morsink et al. (2007); Cadeau et al. (2007); Psaltis & Özel (2014); Watts et al. (2016); Nättilä & Pihajoki (2018); Poutanen (2020); Zavlin et al. (1996); Ho & Lai (2001, 2003); Suleimanov et al. (2012); Salmi et al. (2020, 2023); AlGendy & Morsink (2014); Viironen & Poutanen (2004); Bogdanov et al. (2008).

The theoretical waveform is then propagated through the astrophysical and instrumental layers that connect the source to the data, namely, interstellar absorption and the instrument response. We address each component explicitly in dedicated subsections.

Our approach builds on modeling practices summarized by Bogdanov et al. (2019b, 2021); Salmi et al. (2023) and the open-source X-PSI package (Riley et al. 2023b), which we use as standard benchmarks. We extend these frameworks in a self-contained manner by providing fully relativistic derivations and implementations, applying controlled approximations only when they improve efficiency without biasing inference.

For orientation, Section 2 is organized as follows. Section 2.1 defines the physical setup and notation (parameters and geometry). Section 2.2 specifies the spacetime metric and the oblate-surface model. Section 2.3 details photon propagation (light bending and time-of-flight delays). Section 2.4 incorporates the rotational SR Doppler effect. Section 2.5 models the surface emission and atmospheric treatment, integrating to obtain the observed flux. Section 2.6 discusses the default surface temperature map deployed in this work and the surface discretization strategy. Section 2.7 discusses interstellar absorption and instrument response. Section 2.8 provides the detailed road map of our computation recipe, and the GPU implementation is discussed in Section 2.9. The objective is a production-grade, ready-to-use codebase for future PPM studies.

2.1 Geometrical setup

We modeled a rotating neutron star of mass M, equatorial radius R_e, and spin frequency ν_∗ (angular frequency of Ω ≡ 2πν_∗). At the equator of a typical NS, the dimensionless surface speed is $${\beta }_{\text{eq}}\equiv \frac{v_{\text{eq}}}{c}\approx 0.1–0.3,$$(1)

so both GR and SR effects are necessary for sub-percent-level waveform accuracy.

Photons emitted from the oblate stellar surface experience light bending, gravitational redshift, time-of-flight delays, and Doppler boosting. To place these ingredients in a coherent framework, we adopted a star-centered coordinate system with the observer at infinity along a fixed line of sight. The stellar surface was discretized into small elements (i.e., surface patches) and we traced photons from each patch to the observer, integrating the contributions over the surface. Rather than restricting emission to one or a few circular, isothermal spots, our algorithm is designed from the outset to accept full-surface temperature maps, making it suitable for nontrivial configurations (e.g., accretion-heated distributions or physics-motivated hotspot maps). This relaxes symmetry assumptions that are typical for spot-based models; consequently, we emphasize computational optimizations to ensure the tractability of full-surface maps.

Our geometric notation follows Loktev et al. (2020) for consistency with prior literature. A schematic of the setup is shown in Figure 1. For clarity, the key unit vectors and angles used in this work are summarized in Tables 1 and 2.

Physically, i denotes the observer’s inclination; (θ, ϕ) specify the location of an emitting surface element; α is the emission angle measured from the local radial direction in the near-star static frame; and ψ is the angle between the surface radial vector and the observer’s line of sight. In GR, regions with ψ > π/2 remain visible due to gravitational light bending; the maximum visible ψ is determined by the compactness (defined below).

For a spherical star the surface normal coincides with the radial unit vector, $\widehat{n}=\widehat{r}$; hence, η = 0 and σ = α. For an oblate star, rotational flattening tilts the surface normal away from ȓ; hence, η ≠ 0. We quantify this with the dimensionless surface-slope parameter, $$f\equiv \frac{1}{R\sqrt{1-u}}\frac{dR}{d\theta },u\equiv \frac{R_{\text{S}}}{R},$$(2)

with R_S ≡ 2GM/c² the Schwarzschild radius. We have $$\cos \eta =\widehat{r}\cdot \widehat{n}=\frac{1}{\sqrt{1+f^2}},\sin \eta =\frac{f}{\sqrt{1+f^2}}.$$(3)

The emission angle with respect to the local normal is $$\begin{array}{l}\cos \sigma =\cos \alpha \cos \eta +\frac{\sin \alpha \sin \eta }{\sin \psi }\hfill \\ \times (\cos i\sin \theta -\sin i\cos \theta \cos \varphi ).\hfill \end{array}$$(4)

These definitions lay the foundation for computing the observed flux. This subsection provides a robust geometric framework that connects the static surface description with the fully relativistic treatment used in our forward model.

Fig. 1 Demonstration of the geometric setup used in our X-ray PPM forward computation, shown in the lab-frame coordinate system $\widehat{x}-\widehat{y}=\widehat{z}$. A surface patch (orange square) on the stellar surface at colatitude θ and rotational phase ϕ, viewed along the observer’s line-of-sight unit vector, $\widehat{k}$. Definitions of all unit vectors are given in Table 1 and angle conventions are summarized in Table 2.

2.2 Spacetime and oblate-Schwarzschild (OS) approximation

We go on to specify the exterior spacetime and the rotation-induced stellar shape. Throughout this work, we adopted the oblate-Schwarzschild (OS) approximation: we traced the photon geodesics in the Schwarzschild spacetime of mass, M, while representing the stellar surface as rotationally flattened using the empirical shape prescription of Morsink et al. (2007), as refined by AlGendy & Morsink (2014) and Silva et al. (2021). This extends the classic “Schwarzschild+Doppler” (S+D) treatment (Poutanen & Beloborodov 2006) by including surface oblateness, while neglecting frame dragging and the external quadrupole moment of the metric. The Schwarzschild exterior metric is expressed as $$d{s}^2=-\left(1-\frac{R_{\text{S}}}{r}\right)c^{2}d{t}^2+{\left(1-\frac{R_{\text{S}}}{r}\right)}^{-1}d{r}^2$$(5) $$+r^{2}d{\theta }^2+r^2{\sin }^2\theta d{\varphi }^2.$$(6)

The g_tt component determines gravitational redshift and time dilation, while spatial curvature bends null geodesics. For a spherical star, the stellar radius is constant, R = R_e, while the comoving-frame area element is $$d{S}^{\prime }=R_0^{2}d\cos \theta d\varphi .$$(7)

Rotation deforms the star because centrifugal support reduces the effective gravity near the equator. We modeled the radius profile as $$\begin{array}{l}R=R_{\text{e}}(1+o_2{\cos }^2\theta ),\hfill \\ o_2={\text{Ω}}^2(-0.788+1.030x),\hfill \end{array}$$(8)

where o₂ is the shape coefficient and the dimensionless frequency, $\overline{\Omega }$, and compactness parameter, x, are defined as $$\overline{\text{Ω}}\equiv 2\pi {\nu }_{*}\sqrt{\frac{R_{\text{e}}^3}{GM}},x\equiv \frac{GM}{c^2{R}_{\text{e}}}.$$(9)

The quadratic term encodes polar flattening (o₂ < 0) and remains accurate for ν_∗ ≲ 700 Hz and R_S/R_e ≲ 1, with fractional errors ≲1% when higher-order terms are neglected. For oblate star, the corresponding surface area element is expressed as $$d{S}^{\prime }=R^2\sqrt{1+f^2}d\cos \theta d\varphi ,$$(10)

with f defined in Equation (2). This ensures that the element measures the proper area of the tilted surface. The effective surface gravity varies with latitude due to the combined effects of oblateness and rotation. Following AlGendy & Morsink (2014); Bogdanov et al. (2019b), we have $$\begin{array}{l}g(\theta )/g_0=1+{\sin }^2\theta [{\overline{\text{Ω}}}^2(-0.791+0.776x)\hfill \\ +{\overline{\text{Ω}}}^4(-1.315x+2.431x^2)+{\overline{\text{Ω}}}^6(-1.172x)]\hfill \\ +{\cos }^2\theta [{\overline{\text{Ω}}}^2(1.138-1.431x)\hfill \\ +{\overline{\text{Ω}}}^4(0.653x-2.864x^2)+{\overline{\text{Ω}}}^6(0.975x)\hfill \\ -{\overline{\text{Ω}}}^4(13.47x-27.13x^2)]\hfill \\ +\left|\cos \theta \right|{\overline{\text{Ω}}}^4(13.47x-27.13x^2).\hfill \end{array}$$(11)

Here, the reference gravity g₀ of the corresponding spherical star is $$g_0=\frac{GM}{R_{\text{e}}^2}\frac{1}{\sqrt{1-R_{\text{S}}/R_{\text{e}}}}.$$(12)

At the equator (θ = π/2), gravity is reduced, potentially affecting atmospheric structure and beaming, whereas at the poles it is enhanced.

The OS approximation is well suited to the millisecond pulsars targeted by NICER (typical spins ν_∗ ≲ 400 Hz) and delivers ≲1% waveform accuracy in the tests we consider. Errors remain modest up to ∼600–700 Hz, which motivates recent refinements of the OS approximation (Jakab & Morsink 2025). At higher spins (≳ 700–800 Hz), frame dragging and the external metric’s quadrupole terms become non-negligible, and Hartle–Thorne or fully relativistic rotating-star spacetime might be required. These choices fix the geometric setting for the rest of the calculation: null geodesics in Schwarzschild provide the mapping from the local emission direction and emission time at the surface to the observed angle and phase, naturally incorporating gravitational light bending and time-of-flight delays.

2.3 Photon propagation: light bending and time delays

The dominant GR effects relevant to PPM are gravitational light bending, gravitational redshift, and time-of-flight delays in the curved spacetime around the neutron star. Gravitational light bending deflects photon trajectories away from straight lines, increasing the visible fraction of the surface; regions with ψ > π/2 that would be hidden in Newtonian geometry can remain visible. Time-of-flight delays introduce phase shifts between photons emitted at different surface locations and directions, modifying the waveform shape and its harmonic content. Accurate treatment of these effects is essential because modeling errors can bias the inferred compactness GM/(Rc²) and hence the derived M and R, rather than merely inflating statistical uncertainties.

Gravitational light bending is quantified by the deflection angle ψ for a photon emitted at radius, R, with the emission angle, α, measured from the local radial direction. For outgoing rays (0 ≤ α < π/2), it is given by the Schwarzschild null-geodesic integral, $$\psi (u,\alpha )={{\displaystyle \int }}_R^{\infty }\frac{dr}{r^2}{\left[\frac{1}{b^2}-\frac{1}{r^2}\left(1-\frac{R_{\text{S}}}{r}\right)\right]}^{-1/2},$$(13)

where the impact parameter is $b=R\sin \alpha /\sqrt{1-u}\text{\hspace{0.17em}and\hspace{0.17em}}u=R_S/R$. This integral derives from the null geodesic equation in the Schwarzschild metric, reflecting how gravity pulls photons inward.

For inward emission (α > π/2), the photon path passes through periastron before escaping via $$\psi (u,\alpha )=2{\psi }_{\text{max}}-\psi (u,\pi -\alpha ),$$(14)

where ψ_max ≡ ψ(u_p, π/2) and the periastron compactness is u_p = u/ρ, with $$\rho =\frac{-2\sin \alpha }{\sqrt{3(1-u)}}\cos \left(\frac{1}{3}\left[\arccos \left(\frac{3u\sqrt{3(1-u)}}{2\sin \alpha }\right)+2\pi \right]\right).$$(15)

Using Equations (13)–(15), we precomputed two-dimensional (2D) lookup tables that map $\left(u,\cos \psi \right)\mapsto \left(\cos \alpha ,d\cos \alpha /d\cos \psi \right)$; these speed up ray tracing while retaining a high level of accuracy.

To calculate the time delay for a spherical star, we define zero delay for a photon emitted from the point directly beneath the observer (along $\widehat{k}$). The gravitational and geometric delay for a photon emitted at the same radius R but at angle α < π/2 is $$\frac{c\text{Δ}{t}_p}{R}(u,\alpha )=\frac{1}{R}{{\displaystyle \int }}_R^{\infty }\frac{dr}{1-R_{\text{S}}/r}\left\{{\left[1-\frac{b^2}{r^2}\left(1-\frac{R_{\text{S}}}{r}\right)\right]}^{-1/2}-1\right\}.$$(16)

For an oblate star, an additional delay arises from the difference between the actual emission radius, R, and the equatorial radius, R_e, $$c\delta t(R)=R_{\text{e}}-R+R_{\text{S}}\ln \left(\frac{R_{\text{e}}-R_{\text{S}}}{R-R_{\text{S}}}\right).$$(17)

For direct emission (α < π/2), the total delay is ∆t = ∆t_p + δt. For inward emission (α > π/2), the photon first moves inward and then outward; the total delay is $$\begin{array}{l}\frac{c\text{Δ}t}{R}(u,\alpha )=\frac{c}{R}[2\text{Δ}{t}_p(u_p,\pi /2)-\text{Δ}{t}_p(u,\pi -\alpha )\hfill \\ +2\delta t(R_p)-\delta t(R)],\hfill \end{array}$$(18)

where R_p = ρR is the periastron radius. Precomputing $\left(u,\cos \alpha \right)\mapsto \left(c\Delta t/R\right)$ yields fast, stable phase corrections for rotating stars. For a detailed explanation and discussion about the derivation above, we refer to Salmi et al. (2018, Section 2.1 and Appendix A). Interpolating from these three tables during likelihood evaluations avoids repeated geodesic integrations and is therefore essential for Bayesian inference, which typically requires ∼10⁶–10⁸ forward-model evaluations.

2.4 Rotation: Relativistic Doppler effect

We incorporate rotation by superimposing SR surface effects on the full GR propagation described above: we Lorentz-boosted the emitted photons in comoving frame into the local static frame on the oblate surface and then trace them along Schwarzschild null geodesics (e.g., AlGendy & Morsink 2014; Bogdanov et al. 2019b). For millisecond pulsars with spins ν_∗ ∼ 300 Hz (typical equatorial β_eq ≡ v/c ∼ 0.1), this “OS + Doppler” treatment provides an excellent approximation, whereas higher spins may necessitate Hartle–Thorne or fully relativistic rotating-star metrics. Special-relativistic Doppler boosting generates pulse asymmetry and higher harmonics, which are crucial for recovering the hotspot geometry on the neutron-star surface.

The surface velocity, accounting for the gravitational-redshift factor in the metric, is $$\beta =\frac{2\pi {\nu }_{*}R\sin \theta }{c\sqrt{1-u}},\gamma =\frac{1}{\sqrt{1-{\beta }^2}},$$(19)

where the $\sqrt{1-u}$ factor accounts for the gravitational-redshift correction to the rotating frequency.

The relativistic Doppler factor $$\delta =\frac{1}{\gamma (1-\beta \cos \xi )},$$(20)

modifies the observed photon energy, E_obs, in static frame relative to the emitted energy in the comoving frame, E_emit, as E_obs = δE_emit. This factor depends on the velocity, β, of the emitting region and the angle, ξ, of the velocity vector with respect to the emitting direction. It tends to amplify energy and intensity for photons from approaching patches (cos ξ > 0) and diminishes them for receding ones. The angles in the comoving frame (where the patch is at rest) are transformed via Doppler aberration, $$\mu =\cos {\sigma }^{\prime }=\delta \cos \sigma ,\cos {\alpha }^{\prime }=\delta \cos \alpha .$$(21)

Quantities with a prime (′) are measured in the comoving frame. A detailed derivation of the above result is presented in Appendix A.

The rapid rotation of the neutron star means that the emission from a surface hotspot will be continuously Doppler-shifted throughout the star’s rotational phase. As the hotspot moves across the stellar disk towards the observer, its emission is progressively blueshifted, reaching maximum blueshift as it approaches the line of sight. Conversely, as it recedes, the emission becomes increasingly redshifted. This modulation of photon energy directly impacts the observed flux and the shape of the pulse profile. Therefore, any realistic model of a pulse profile must incorporate the relativistic Doppler effect to correctly interpret the observed light curves and accurately constrain fundamental neutron-star parameters. These SR effects couple directly to the intrinsic atmospheric beaming, as discussed in the next section.

2.5 Surface emission and atmospheric beaming

When photons escape from surface hotspots, a simple starting point is to model the emergent spectrum as a blackbody in the comoving frame. For quantitative inference, however, physically motivated neutron-star atmosphere models are required. Thin (likely accreted) hydrogen or helium layers produce angle-and energy-dependent specific intensities (“beaming”) and spectra that deviate from a blackbody spectrum, typically appearing harder and exhibiting energy- and composition-dependent limb darkening or brightening. For rotation-powered millisecond pulsars with magnetic fields B ∼ 10⁸–10⁹ G, weakly magnetized hydrogen atmospheres are commonly adopted; whereas strongly magnetized models are required at higher fields. We followed this atmosphere-based approach throughout, drawing on established calculations and recent updates (Zavlin et al. 1996; Ho & Lai 2001, 2003; Suleimanov et al. 2012; Salmi et al. 2020, 2023; Choudhury et al. 2024b).

According to (Bogdanov et al. 2019b, Equation (31)), the number flux of photons can be calculated as $$\begin{array}{l}df(E)=\sqrt{1-u}{\delta }^3\mu \frac{d\cos \alpha }{d\cos \psi }\frac{I(E^{\prime },\mu )}{E}\hfill \\ \times \frac{\gamma R^2\sqrt{1+f^2}d\cos \theta d\varphi }{D^2},\hfill \end{array}$$(22)

where df denotes the photon count rate per unit detector area and per unit observed photon energy; I is the atmospheric specific intensity; E and E′ denote the observed and emitted photon energies, respectively; and D is the distance between the star and the observer. A detailed derivation of this expression is provided in Appendix C. Integrating over all visible surface elements, d cos θ dϕ, yields the full pulse profile.

For realistic beaming, we followed the procedure of Bogdanov et al. (2021). The characteristic field strength at which magnetic effects become non-negligible in the atmosphere is $B_0\equiv e^3{m}_e^{2}c/{\hslash }^3\simeq 2\times {10}^9\text{\hspace{0.17em}G}$. Magnetic effects are negligible in this context because surface fields are smaller than B₀. We also fixed the composition and ionization. After accretion ceases, gravitational settling stratifies the outer layers, leaving the lightest element at the photosphere; for recycled MSPs, this is typically hydrogen, with helium as a plausible variant. Over the NICER band and at T_eff ∼ (0.5–2) × 10⁶ K, differences in angular beaming between H and He are modest (a few percent below 1 keV and at most tens of percent near 2 keV). However, at grazing angles (µ ≲ 0.5, σ′ ≳ 60^◦) these differences can become substantial (see Salmi et al. 2023). We therefore adopted a pure, nonmagnetic H atmosphere as the baseline and treat a pure He atmosphere as a systematic variant. At these temperatures the neutral fraction is small and for µ ≳ 0.5 fully ionized opacities reproduce the beaming of partially ionized solutions to within a few percent over 0.5–2 keV. Accordingly, we use fully ionized models for I_ν(µ). We note that non-local thermodynamic equilibrium (NLTE) and Compton corrections are also negligible in this regime (sub-percent below 2 keV). These choices keep the beaming physics self-consistent with the magnetic-field assumption above and isolate composition and ionization effects to small, quantifiable systematics in PPM inference. Throughout, we adopt the standard assumption of a deep-heated atmosphere in which external energy is deposited below the photosphere. Alternative heating geometries can alter the emergent spectrum and beaming (e.g., Bauböck et al. 2019; Salmi et al. 2020) and full consistency with the magnetosphere would require modeling the return-current particle velocities and their stopping column depths, which we defer to future work for improvement.

To compute realistic beaming I_ν(µ), we did not solve the radiative transfer during PPM inference. Instead, we used precomputed, high-resolution lookup tables of $\log \left(I{T}_{\text{eff}}^{-3}\right)$ as a function of log T_eff, log g, µ, log(E/kT_eff) The hydrogen grid spans log T_eff ∈ [5.1, 6.8] (helium: [5.1, 6.5]), log g ∈ [13.7, 14.7], µ ∈ [10⁻⁶, 1] with extra nodes near µ ≃ 0, 1, and log(E/kT_eff) ∈ [−1.3, 2.0]. Production analyses use the fully ionized H grid. An identically sampled He grid is propagated as a composition systematic. For the full details of this treatment, we refer to Bogdanov et al. (2021) for further discussion.

With the atmosphere model established, the near-surface radiative physics is specified. Two additional components complete the forward model used for data comparison: (i) the surface temperature map and its numerical discretization on the oblate star; and (ii) the propagation to the detector, including interstellar absorption and the instrument response. The next two subsections develop these ingredients in detail. We now transition from the GR+SR forward model to its observational implementation, which converts model photon fluxes into measured counts.

2.6 Surface temperature maps and surface discretizations

Unlike most published PPM studies, which, for computational reasons, typically discretize only parametrized hotspot regions, our framework natively discretizes the entire stellar surface temperature map. Prior work has largely employed one- or two-spot models with simple geometries (e.g., circular caps) and single or dual temperature components. While effective for many applications, such prescriptions cannot capture spatially complex temperature structures within a spot or extended patterns across the surface. Full-surface discretization, in contrast, accommodates nontrivial temperature distributions and global morphologies, which are required in several physical scenarios (e.g., accreting neutron stars with extended heating patterns and nonuniform temperature distribution; see Das et al. 2022; Das et al. 2025).

In our codebase, the primary input is a user-specified full-surface map T_eff(θ, ϕ) defined on an oblate surface mesh in colatitude, θ, and azimuth, ϕ. The mesh resolution is user-controlled and can be increased until interpolation errors are sub-percent in the phase–energy domain.

For benchmarking with results from the literature, we also provide parameterizations that reproduce standard geometric hotspot models used in previous PPM works. For example, a single circular cap is defined by its center, (θ₀, ϕ₀), angular radius, ζ, and spot temperature ,T_spot; all mesh elements with angular separation ≤ ζ from (θ₀, ϕ₀) are assigned T_spot, while elements outside the cap are set to 0, indicating that no photons are emitted. The same strategy extends to antipodal two-cap configurations and other predefined shapes.

Our default physics-motivated generator for whole-surface temperature maps follows Huang & Chen (2025), which derives return-current heating patterns from a force-free magnetosphere with an off-centered dipole. In that study, the implementation used the X-PSI (Riley et al. 2023a,b) full-surface module and achieved per-map generation times of order one second on a CPU, sufficient for static tests and small-scale inference but not optimized for large-scale PPM inference and for testing different sampler setups for robustness. In the present framework, we adopt this physics-motivated hotspot module as the default for hotspot configurations in future PPM inference. For the full algorithmic details and validation of the hotspot configuration, see Huang & Chen (2025).

2.7 Interstellar absorption and instrument response

Photons emerging from the stellar surface are attenuated by interstellar absorption, which we model as a multiplicative factor applied to the intrinsic spectrum. Accounting for this energy-dependent attenuation is essential for linking theoretical pulse profiles to the observed data and enabling reliable parameter inference. In this work, we model line-of-sight attenuation with a multiplicative transmission factor applied to the emitted photon spectrum, $$f_{\text{ISM}}(E;N_{\text{H}})=\exp [-{\sigma }_{\text{ISM}}(E)N_{\text{H}}]f_0(E),$$(23)

where f₀(E) is the source-frame spectral flux, N_H is the hydrogen column, and σ_ISM(E) is the energy-dependent cross section from the tbnew prescription (modern gas-phase abundances and grain physics) (Wilms et al. 2000). For efficiency, we use a precomputed attenuation vector on the NICER energy grid at a fiducial $N_{\text{H}}^{\text{ref}}$ and exploit the linear scaling in the exponent to obtain transmission for arbitrary N_H. The MSPs targeted for PPM have very low columns (∼10²⁰ cm⁻²), so extinction is small over 0.3–2 keV and, being phase independent, has negligible impact on geometric inferences; we nevertheless carry N_H as a nuisance parameter and propagate it through the likelihood.

Another component that connects the theoretical pulse profile to the observed data is the instrument response. Here, our instrument response follows the same strategy implemented in previous NICER works as a standard pipeline. Comparison to NICER data is performed in detector space by convolving the attenuated model spectrum with the instrument response, $$C(N)={{\displaystyle \int }}_0^{\infty }R(N,E)f_{\text{ISM}}(E;N_{\text{H}})dE,$$(24)

where C(N) is the predicted count rate in channel N (counts s⁻¹) and R(N, E) is the full response (RMF×ARF), namely, the probability of detecting a photon of energy, E, in channel, N, multiplied by the effective area. Energy-resolved pulse profiles are generated by evaluating f₀(E) from the atmospheric beaming tables for each visible surface element and phase bin, applying the ISM transmission, and convolving with R(N, E). We verified our forward model by reproducing standard XSPEC results for simple spectra with the same response files and by crosschecking independent implementations. All production fits use the identical RMF+ARF set and channelization as the data to ensure consistency.

2.8 Computation recipe

In this section, we outline the step-by-step computational procedure for generating a phase-resolved spectrum from a rotating neutron star. The process (i) discretizes the stellar surface, (ii) precomputes gravitational light bending, (iii) computes the observed flux from each surface element, and (iv) convolves the result with models for interstellar absorption and the instrument response.

2.8.1 Surface discretization

To accurately model the emission from the star’s surface, we first divide it into a finite number of small patches or pixels. For this, we employ the HEALPix (Hierarchical Equal Area isoLatitude Pixelization) scheme (Górski et al. 2005). This method is particularly advantageous because it divides the sphere into patches of equal area, ensuring that each patch contributes proportionally to the total flux and simplifying the geometric calculations. The resolution of this grid is controlled by the parameter N_side. Let N_θ be the number of rings and i_θ = 1,…, N_θ the ring index. Each ring is further discretized into N_ϕ azimuthal bins with i_ϕ = 1, …, N_ϕ. The center coordinates (θ, ϕ) of each surface patch are computed from (i_θ, N_θ, i_ϕ, N_ϕ) according to the chosen discretization scheme.

For 1 ≤ i_θ < N_side and 1 ≤ i_ϕ ≤ 4i_θ, we have $$\cos \theta (i_{\theta },i_{\varphi },N_{\text{side}})=1-\frac{i_{\theta }^2}{3N_{\text{side}}^2},$$(25) $$\varphi (i_{\theta },i_{\varphi },N_{\text{side}})=\frac{\pi }{2i_{\theta }}\left(i_{\varphi }-\frac{1}{2}\right).$$(26)

For N_side ≤ i_θ ≤ 2N_side and 1 ≤ i_ϕ ≤ 4N_side, we have $$\cos \theta (i_{\theta },i_{\varphi ,}{N}_{\text{side}})=\frac{4}{3}-\frac{2i_{\theta }}{3N_{\text{side}}},$$(27) $$\varphi (i_{\theta },i_{\varphi },N_{\text{side}})=\frac{\pi }{2N_{\text{side}}}\left(i_{\varphi }-\frac{(i_{\theta }-N_{\text{side}}+1)\mathrm{mod}2}{2}\right).$$(28)

The southern hemisphere is a mirror image of the northern hemisphere, i.e., for 2N_side < i_θ ≤ 4N_side − 1 and 1 ≤ i_ϕ ≤ 4 min(4N_side − i_θ, N_side), we have $$\cos \theta (i_{\theta },i_{\varphi },N_{\text{side}})=-\cos \theta (4N_{\text{side}}-i_{\theta },i_{\varphi },N_{\text{side}}).$$(29) $$\varphi (i_{\theta },i_{\varphi },N_{\text{side}})=\varphi (4N_{\text{side}}-i_{\theta },i_{\varphi },N_{\text{side}}).$$(30)

In the HEALPix discretization scheme, the surface of the sphere is partitioned into $12N_{\text{side}}^2$ equal-area patches, so each patch covers a solid angle of $d\cos \theta d\phi =\pi /\left(3N_{\text{side}}^2\right)$.

In modeling complex emission geometries, we identify the surface patches that belong to the emitting region and designate them as active patches, whose total number is N_active. This discretization can introduce a small scaling bias at sharp boundaries because the pixelization may slightly over- or under-represent the true emitting area. A more careful edge treatment could mitigate this bias and potentially reduce computational load, but given that performance is satisfactory under the current setup, we defer such refinements to a future code release.

2.8.2 Lensing table generation

As discussed in Section 2.3, computing gravitational light bending and the associated time-of-flight delays for every surface patch is computationally expensive. To accelerate the calculation, we precompute three lookup tables: (i) the cosine of the emission angle, cos α; (ii) the lensing factor, d cos α/d cos ψ; and (iii) the time delay, ∆t.

These tables are indexed by grids defined over the key physical parameters. We use a helper function to define linearly spaced arrays: $$A=\text{Linspace}(x_{\text{min}},x_{\text{max}},x_{\text{num}}),$$(31)

to denote an array A of size x_num, whose elements are given by $$A[i]=x_{\text{min}}+(x_{\text{max}}-x_{\text{min}})\frac{i}{x_{\text{num}}-1},$$(32)

where i = 0, …, x_num − 1. The grids for the compactness (u), the cosine of the final bending angle (cos ψ), and the cosine of the emitting angle (cos α) are defined as $$\begin{array}{l}u=\text{Linspace}(u_{\text{min}},u_{\text{max}},N_u),\hfill \\ c_{\psi }=\text{Linspace}(c_{\psi ,\text{min}},c_{\psi ,\text{max}},N_{c_{\psi }}),\hfill \\ c_{\alpha }=\text{Linspace}(c_{\alpha ,\text{min}},c_{\alpha ,\text{max}},N_{c_{\alpha }}).\hfill \end{array}$$

For all calculations, we adopted the following empirically chosen values, selected to ensure all lookups remain within bounds for typical neutron-star parameters: $$\begin{array}{l}{u}_{\text{min}}=0.3,u_{\text{max}}=0.5,N_u=256,\hfill \\ c_{\psi ,\text{min}}=-0.9,c_{\psi ,\text{max}}=1.0,N_{c_{\psi }}=512,\hfill \\ c_{\alpha ,\text{min}}=-0.2,c_{\alpha ,\text{max}}=1.0,N_{c_{\alpha }}=512.\hfill \end{array}$$

For α < π/2, $$\psi (u,\alpha )=\frac{2\sin \alpha }{\sqrt{1-u}}{{\displaystyle \int }}_0^1\frac{xdx}{\sqrt{{\cos }^2\alpha +x^{2}q}},$$(33) $$q=[2-x^2-u{(1-x^2)}^2/(1-u)]{\sin }^2\alpha ,$$(34)

and for α > π/2, $$\psi (u,\alpha )=2\psi (u/\rho ,\pi /2)-\psi (u,\pi -\alpha ).$$(35)

We created two 2D arrays, C_α and LF, both of a size (N_u, $N_{C_{\psi }}$), to store precomputed lookup tables for the cosine of the local emission angle cos α and the lensing factor d cos α/d cos ψ. For each u[i_u] in the array u, we set α_cur = 0, calculated ψ_cur = ψ(u[i_u], α_cur), and appended (α_cur, ψ_cur) to a list {(α, ψ), . . .}. We then increased α_cur by dα = π/2048 (an empirically chosen constant; it can be decreased for higher accuracy) until ψ_cur > arccos c_ψ,min. With a simple transformation, we converted the list {(α, ψ), . . .} to {(cos ψ, cos α), . . .}; using numerical differentiation, we obtain the list {(cos ψ, d cos α/d cos ψ), . . .}. Then, for each $c_{\psi }\left[i_{C_{\psi }}\right]$ in the array c_ψ, we interpolated $c_{\psi }\left[i_{C_{\psi }}\right]$ within these two lists to obtain the corresponding $C_{\alpha }\left[i_u,i_{C_{\psi }}\right]$ and $LF\left[i_u,i_{C_{\psi }}\right]$. Next, we moved on to the calculation of the time delay. For α < π/2, $$\frac{c\text{Δ}{t}_p}{R}(u,\alpha )=\frac{2{\sin }^2\alpha }{1-u}{{\displaystyle \int }}_0^1\frac{xdx}{\sqrt{{\cos }^2\alpha +x^{2}q}\left(1+\sqrt{{\cos }^2\alpha +x^{2}q}\right)}.$$(36)

For α > π/2, $$\frac{c\text{Δ}{t}_p}{R}(u,\alpha )=2\frac{c\text{Δ}{t}_p}{R}\left(u/\rho ,\pi /2\right)-\frac{c\text{Δ}{t}_p}{R}(u,\pi -\alpha ).$$(37)

We used a 2D array TD of size (N_u, $N_{C_{\alpha }}$) to store the lookup table of time delay as a function of (u, c_α). For each grid point $\left(u\left[i_u\right],c_{\alpha }\left[i_{C_{\alpha }}\right]\right)$, we directly calculated the time delay using the integral above and stored the result, expressed as $$TD[i_u,i_{c_{\alpha }}]=\frac{c\text{Δ}{t}_p}{R}(u[i_u],\arccos c_{\alpha }[i_{c_{\alpha }}]).$$(38)

The light-bending integral (Equation (33)) and the time-delay integral (Equation (36)) are given in Salmi et al. (2018, Appendix A).

2.8.3 Number flux calculation

This step computes the observed flux from each surface patch and assembles the pulse profile. We used a widely employed “shift-and-add” technique: rather than simulating every active patch across all rotation phases, we first computed a detailed light curve for a single pseudo-patch located at ϕ = 0 on a given colatitude ring θ. This template light curve was then phase-shifted and added to the total to account for the actual longitudes of the active patches on that ring.

For each ring indexed by i_θ at a colatitude, θ, we determined which patches are active and collect the azimuthal angles, ϕ, of the active patches into an array, ϕ_a, of a length, $N_{{\phi }_{\alpha }}$. If $N_{{\phi }_{\alpha }}=0$, the ring contains no active patches and is skipped. Otherwise, we prepared four arrays to store the pseudo-patch light-curve data at ϕ = 0. These were reused later when accumulating contributions from the active patches: $$\begin{array}{cccc}F[i_{fp}],& z[i_{fp}],& {\phi }_o[i_{fp}],& \mu [i_{fp}]\end{array}$$(39)

where $i_{\text{fp}}=0,\dots ,N_{\text{fp}}$ and $N_{\text{fp}}$ is the number of fine rotation-phase points, this parameter controls precision. When the pseudo-patch rotates to azimuth $\varphi \left(i_{\text{fp}}\right)=2\pi i_{\text{fp}}/N_{\text{fp}}$, we perform the following for each $i_{\text{fp}}\in \left[0,N_{\text{fp}}\right]$:

1. Initialize invisible-range indices $i_{\text{inv\_min}}=N_{\text{fp}}$ and i_{inv_max} = 0.

2. Compute cos ψ = sin i sin θ cos ϕ + cos i cos θ.

   • Use the lookup table to obtain cos α and the lensing factor d cos α/d cos ψ.

   • Evaluate cos σ from Equation (4): $$\cos \sigma =\cos \alpha \cos \eta +\frac{\sin \alpha \sin \eta }{\sin \psi }\times (\cos i\sin \theta -\sin i\cos \theta \cos \varphi )$$(40)

   • If cos σ < 0, the point is invisible (and if $i_{\text{fp}}=0$, the entire ring is invisible). Update $i_{\text{inv\_min}}=\min \left(i_{\text{inv\_min}},i_{\text{fp}}\right)$ and $i_{\text{inv\_max}}=\max \left(i_{\text{inv\_max}},i_{\text{fp}}\right)$. Otherwise, set $$\begin{array}{l}F[i_{fp}]=\sqrt{1-u}{\delta }^3\cos {\sigma }^{\prime }\frac{d\cos \alpha }{d\cos \psi }\frac{\gamma R^2\sqrt{1+f^2}d\cos \theta d\varphi }{D^2},\hfill \\ z[i_{fp}]=\frac{1}{\delta \sqrt{1-u}},\hfill \\ {\phi }_o[i_{fp}]=\frac{\varphi }{2\pi }+{\nu }_{\ast }\text{Δ}t,\hfill \\ \mu [i_{fp}]=\delta \cos \sigma .\hfill \end{array}$$(41)

If a ring is partially visible, we must handle the phase segment where it is hidden. For these invisible points $(i_{\text{inv\_min}}\le i_{\text{fp}}\le i_{\text{inv\_max}})$, we set the flux and apparent emission angle to zero. To ensure continuity and avoid numerical artifacts, the redshift and observed phase are smoothly filled using linear interpolation between the values at the points of disappearance (i_{inv_min} − 1) and reappearance (i_{inv_max} + 1), namely, $$F[i_{fp}]=0,\mu [i_{fp}]=0,$$(42) $$\begin{array}{ll}z\left[i_{{\text{f}}_{\text{p}}}\right]=\hfill & L(i_{{\text{f}}_{\text{p}}};i_{\text{inv\_min}}-1,i_{\text{inv\_max}}+1,\hfill \\ \hfill & z\left[i_{\text{inv\_min}}-1\right],z\left[i_{\text{inv\_max}}+1\right]),\hfill \end{array}$$(43) $$\begin{array}{l}{\phi }_o[i_{fp}]=L(i_{fp};i_{inv\text{\_}min}-1,i_{inv\text{\_}max}+1,\\ \begin{array}{cc}& {\phi }_o[i_{inv\text{\_}min}-1],{\phi }_o[i_{inv\text{\_}max}+1])\end{array}.\end{array}$$(44)

where L denotes linear interpolation: $$\text{L}(\text{x};\text{x0},\text{x1},\text{y0},\text{y1})\text{=y0+}(\text{y1-y0})\frac{\text{x-x0}}{\text{x1-x0}}.$$(45)

After preparing these calculations, we used the “shift-and-add” technique to accumulate all active patches’ contributions to the total waveform. We defined the phase points to be output as an array, φ_out, of a size, $N_{{\phi }_{\text{out}}}$, the observed photon energies as an array, E_o, of a size $N_{E_o}$, and the total output photon number flux as a matrix, F_t, of a specific size $\left(N_{{\phi }_{\text{out}}},N_{E_0}\right)$. For each $\varphi ={\varphi }_a\left[i_{{\varphi }_a}\right]$ with $i_{{\varphi }_a}=0,\dots ,N_{{\varphi }_a}-1$, each $\phi ={\phi }_{\text{out}}\left[i_{{\phi }_{\text{out}}}\right]$ with $i_{{\phi }_{\text{out}}}=0,\dots ,N_{{\phi }_{\text{out}}}-1$, and each $E=E_0\left[i_{E_0}\right]$ with $i_{E_0}=0,\dots ,N_{E_0}-1$. Next, we can compute $${\phi }_P=(\phi +\frac{\phi }{2\pi }+1-{\phi }_0[0])\mathrm{mod}1)+{\phi }_0[0].$$(46) $$\begin{array}{l}{F}_p=Interp({\phi }_p;{\phi }_o,F),\\ z_p=Interp({\phi }_p;{\phi }_o,z),\\ {\mu }_p=Interp({\phi }_p;{\phi }_o,\mu ),\end{array}$$(47)

where Interp(x₀; x, y) denotes building an interpolation function from arrays x and y and evaluating it at x₀. The contribution to the total flux is found by scaling, F_p, with the specific intensity of the atmospheric emission, I(z_pE, µ_p), which depends on the redshifted energy and the apparent emission angle. This contribution is then added to the final light-curve matrix, F_t, $$F_t[i_{{\phi }_{out}},i_{E_o}]+=F_p\frac{I(z_{p}E,{\mu }_p)}{E},$$(48)

2.8.4 Interstellar absorption and instrument response

The photon number flux calculated above is a theoretical quantity. To compare it with real data, we must model its interaction with the interstellar medium (ISM) and the response of the detector. We can define logarithmically spaced arrays using the helper function $$A=Logspace(x_{min},x_{max},x_{num})$$(49)

which represents an array A of length x_num whose elements are given by $$A[i]=\exp \left(\ln x_{min}+(\ln x_{max}-\ln x_{min})\frac{i}{x_{num}-1}\right),$$(50)

where i = 0, …, x_num − 1.

Throughout all calculations, we use the grid $$E=Logspace(0.1\text{\hspace{0.17em}keV},\text{}5.215\text{\hspace{0.17em}keV},N_{E_o}=128).$$(51)

We modeled the absorption by the ISM using a precomputed table from the tbnew model. The table provides the attenuation fraction as a function of energy for a reference hydrogen column density $\left(N_{\text{H}}^{\text{ref}}=4\times {10}^{18}{\text{\hspace{0.17em}cm}}^{-2}\right)$. We used col1 to represent the energy column in the table and col2 to represent the attenuation-fraction column. Our simulation uses N_H = 2 × 10¹⁹ cm⁻². Assuming the optical depth scales linearly with N_H, the new attenuation is the original value raised to the power of the ratio $N_{\text{H}}/N_{\text{H}}^{\text{ref}}\text{\hspace{0.17em}(i}.\text{e}.,=5)$ (i.e., = 5). We created an array ATT by interpolating the table to our energy grid and applying this scaling; ATT is of a size denoted as $N_{E_0}$ and elements expressed as $$ATT[i_{E_o}]=Interp{(E_o[i_{E_o}];col1,col2)}^5.$$(52)

Next, we convolve the signal with the instrumental response of the detector (e.g., NICER). This requires two standard data files: an ARF file and an RMF file. The ARF file provides the effective area of the detector across different energies and contains three columns that we denote E_l, E_h, and A_e all with size N_e. The RMF file describes how photons of a given true energy are recorded in different detector output channels and contains a matrix, R, of a specific size (N_e, N_c). We multiplied the effective area from the ARF with the RMF to obtain a fine-resolution response matrix, R_f, with a specific size (N_e, N_c): $$R_f[i_e,i_c]=A_e[i_e]\times R[i_e,i_c].$$(53)

We defined the array, e_i, with a length, N_Eo, whose elements are given by $$e_I\text{[J]=}\{\begin{array}{cc}\text{1,}& ifi=j,\\ \text{0,}& \text{otherwise}.\end{array}$$(54)

Define the downsampling matrix DS with size ($N_{E_0}$, N_e) by $$DS[i_{E_o},i_e]={{\displaystyle \int }}_{E_l}^{E_h}\text{Interp}(x;E_o,e_{i_e})dx.$$(55)

Combining ISM attenuation, ATT, the downsampling matrix, DS, and the fine instrument response, R_f, yields the final, comprehensive response matrix, RSP. This matrix directly converts the theoretical, energy-resolved photon flux into the expected counts per detector channel. The resulting RSP has a specific size ($N_{E_0}$, N_c) and is given by $$RSP[i_{E_o},i_c]={\displaystyle \sum _{i_e=0}^{N_e-1}A}TT[i_{E_o}]\times DS[i_{E_o},i_e]\times R_f[i_e,i_c],$$(56)

2.9 GPU implementation

The NVIDIA GPU architecture is designed for massively parallel computation, executing thousands of threads simultaneously. The core processing unit is the Streaming Multiprocessor (SM), each of which contains hundreds of simpler arithmetic units known as CUDA cores. When a computational task is launched, it is organized into a grid of thread blocks; one or more blocks can be assigned to an SM for execution. The efficiency of this architecture depends on its memory hierarchy. All threads can access a large, high-capacity global memory, which is analogous to CPU random-access memory (RAM) but has high latency, making frequent access a primary performance bottleneck. To mitigate this, each SM is equipped with a small, extremely fast on-chip shared memory. This memory is private to the threads within a single block and acts as a programmable cache. In a typical GPU algorithm for a problem such as the one described below, a block of threads first cooperatively loads the necessary subset of data from global memory into this low-latency shared memory. All subsequent computations can then be performed by accessing the shared memory, reducing memory overhead and unlocking the full parallel processing power of the CUDA cores before writing the final results back to global memory.

In our waveform generation process (as described in Section 2.8), we treated each ring independently. In the GPU implementation, we assigned each ring’s computation to a dedicated SM. For each SM, the four helper arrays F, z, φ_o, and µ were staged in shared memory, whereas the number-flux matrix, F_t, the lensing lookup table, the (optional) numerical atmosphere lookup table, and the RSP matrix reside in global memory. When a kernel is launched, the host (CPU) provides all configuration data describing the task, including the star’s mass, radius, and spin frequency (M, R, ν_∗); the observer’s inclination (i); the hotspot’s shape, position, and size ({θ_i, ϕ_i, ζ_i, . . .}); and the ring indices it operates on (i_θ, N_θ). For each active surface patch, output phase, and output energy, the kernel computes the corresponding number flux and immediately writes the result to the matrix F_t stored in global memory using Equation (48). After all tasks on the SMs were complete, we multiplied the matrix, F_t, by the matrix, RSP, which combines ISM absorption and the detector’s instrument response. We then copy the final result from GPU global memory to CPU host memory for analysis and plotting.

3 OS approximation benchmarks

In this section, we describe how we tested the OS approximation and reproduced the theoretical benchmarks established in Bogdanov et al. (2019b) that were generated by the Illinois– Maryland (IM) group code. At this stage, no atmosphere model or instrument-specific modeling was included; all tests assumed blackbody emission. We therefore evaluated the code accuracy using blackbody spectra and prescribed beaming laws.

Following Bogdanov et al. (2019b), we verified both S+D and OS. The S+D scheme treats special-relativistic effects exactly at the surface, assuming a spherical star and using the exterior Schwarzschild metric; it is very fast and accurate for slowly rotating stars (see also Poutanen & Beloborodov 2006). The OS scheme incorporates the star’s spin-induced oblateness while retaining an exterior Schwarzschild spacetime; it is now the standard in X-ray PPM. Direct comparisons against ray tracing in numerically computed rotating-star metrics show that OS waveforms agree with the exact results to within ≲0.1% for spin frequencies ν_∗ ≲ 300 Hz, whereas S+D can incur percent-level errors at higher spins. Hence, we adopted OS as the default for our production runs (Bogdanov et al. 2019b; Morsink et al. 2007) and implemented the OS approximation in our GPU codebase.

We followed the hotspot geometries and spacetime choices in Bogdanov et al. (2019b). Specifically, we tested both point-like (ζ_spot = 0.01 rad) and extended (ζ_spot = 1 rad) circular hotspots. Beaming patterns included (i) isotropic emission and (ii) simple parametric laws ∝ cos² σ and ∝ sin² σ used for verification in Bogdanov et al. (2019b, Tables 3–4), where σ is the angle between the surface normal and the emission direction. In this section, we assume a canonical neutron-star mass M = 1.4 M_⊙, equatorial radius R_e = 12 km, distance D = 200 pc, and a surface temperature kT = 0.35 keV. We considered monochromatic pulse profiles evaluated at E_obs = 1 keV with a photon number flux of 1 cm⁻² s⁻¹ keV⁻¹, matching the reference tests in Bogdanov et al. (2019b).

We implemented both CPU and GPU versions of our code-base. For validation, we compare our waveforms phase-by-phase against the IM reference waveforms used in Bogdanov et al. (2019b). Consistent with the NICER validation protocol, our target is a fractional accuracy better than 0.1% across phases under likelihood-relevant weighting, where phase bins contribute in proportion to their expected counts. In this sense, we achieve the target accuracy, with larger relative deviations confined to ingress/egress when the hotspot lies behind the stellar limb and the flux is orders of magnitude below the pulse maximum. Because Poisson likelihood contributions scale with the expected counts, these low-flux bins carry little statistical weight in typical datasets. Moreover, the IM treatment about the stellar limb, as implemented in Bogdanov et al. (2019b), is a useful reference but not ground truth. Accordingly, we defer any certification of limb (small-µ) influence on the waveform computation to a dedicated stress test in Sect. 4, Since these systematic errors become significantly more consequential under realistic atmosphere model. These localized discrepancies have a negligible impact on inference because likelihood evaluations are generally insensitive to tiny, low-count portions of the profile.

Across all OS test cases (OS1a–OS1l; see Figure 2 for a subset), our CPU and GPU implementations reproduce the IM reference waveforms within the 0.1% target at essentially all phases outside eclipse ingress/egress. The few outliers occur precisely where the spot rotates out of and back into view; as noted above, these phases have minimal photon counts and do not affect practical parameter estimation. This level of agreement meets the precision standard required for theoretical PPM forward computation and provides a reliable foundation for extensions such as energy-resolved profiles with atmospheric beaming and instrument responses. In the next section, we incorporate a realistic atmosphere model and identify an important interpolation artifact in the widely adopted atmospheric lookup table. We then introduce two diagnostic test cases to quantify its impact and present a more detailed implementation that mitigates the artifact.

Fig. 2 Monochromatic (Eobs = 1 keV) waveforms under the oblate-Schwarzschild (OS) approximation for the OS1 test suite (subset shown), computed with our CPU (blue) and GPU (orange) implementations. The theoretical benchmarks from Bogdanov et al. (2019b) are shown in black. Residuals are taken with respect to the benchmark; red dashed lines indicate ±0.1%. Test-case definitions follow Bogdanov et al. (2019b).

4 Test cases on atmosphere interpolation

In the previous section, we explain how we verified that our implementation reproduces published theoretical benchmarks under the OS approximation and that the resulting light curves agree with canonical test problems. Before turning to production-grade studies, we highlight another common source of systematic error that is missing from much of the existing literature: interpolation within precomputed neutron-star atmosphere tables.

To avoid repeatedly solving the radiative-transfer problem on the fly, modern forward models use predefined lookup tables of the emergent specific intensity. This step is essential for realistic PPM: replacing a blackbody with a physically motivated atmosphere is required for meaningful comparison with NICER observations. In practice, many groups use nonmagnetized atmosphere grids for hydrogen or helium compositions, often referred to as NSX-H and NSX-He, which tabulate I_ν(µ; T_eff, g, E) as a function of photon energy E, emission-angle cosine µ = cos σ′, effective temperature T_eff, and surface gravity g (e.g., Ho & Lai 2001, 2003). In what follows, we examine how interpolation choices on these grids can affect the modeled waveforms, and we propose additional test cases designed to demonstrate the numerical behavior of interpolation on the NSX atmosphere lookup tables. In our experiments, evaluating the NSX atmosphere tables and performing the associated interpolation constitute the dominant computational cost in pulse profile calculations. More importantly, different interpolation schemes can introduce distinct numerical artifacts (e.g., overshooting near table boundaries), which may bias predicted fluxes and spectra quite significantly under our proposed test cases.

Several widely used codes adopt different interpolation strategies for these tables. As reported in Bogdanov et al. (2019b, 2021); Choudhury et al. (2024b), the IM code employs quadratic and quartic polynomial interpolation, X-PSI primarily uses cubic Lagrange interpolation, and the Alberta implementation relies on linear schemes (see Bogdanov et al. 2019b, 2021; Choudhury et al. 2024b for details). Because IM is not open source, its precise interpolation details are not publicly documented. In contrast, the publicly available X-PSI code uses cubic Lagrange interpolation for atmosphere lookups (Riley et al. 2023b; Bogdanov et al. 2021). Given the structure of the tabulated NSX H/He atmosphere grids (Ho & Lai 2001, 2003), these choices can introduce method-dependent artifacts, underscoring the need for careful treatment.

A commonly underappreciated numerical issue in atmosphere-table usage is polynomial overshoot near grid boundaries. Because the widely used nonmagnetized NSX H/He grids have finite spacing and regions of rapid curvature, selecting an interpolation scheme that avoids spurious oscillations while preserving physical constraints (e.g., non-negativity of the specific intensity) is nontrivial. In principle, as the grid becomes sufficiently dense, low-order (e.g., linear) interpolation is typically the most robust and economical. A comprehensive, controlled comparison of interpolation orders and monotone schemes, ideally alongside finer atmosphere grids or direct radiative-transfer recalculations, is deferred to future work.

In our experiments with NSX-H grids, cubic Lagrange interpolation can exhibit boundary instabilities (overshoot/undershoot) that depend on local spacing and curvature, particularly as the emission-angle cosine approaches an endpoint (µ → 0; grazing emission). To enforce non-negativity, the open-source X-PSI implementation clips any negative interpolants to zero.¹ While this safeguard prevents unphysical negative intensities, it still cannot compute the predicted flux downward properly in phases where many rays sample boundary regions.

Motivated by these considerations, we adopt a hybrid interpolation scheme in this work. For the temperature and gravity axes of the NSX-H table, log T_eff and log g (the NSX-H tables tabulate these quantities in logarithmic units), we employ linear interpolation. The grids are sufficiently dense, making linear interpolation both robust and accurate in practice. For the emission-energy axis, we use cubic Lagrange interpolation in log(E/kT_eff); for the emission-angle cosine µ, we use cubic Lagrange interpolation in the interior but switch to linear interpolation within boundary zones at µ → 0. Because a cubic Lagrange interpolant uses four points, we designate the outermost four grid nodes on the small-µ end as “protected” regions in which linear interpolation is enforced.

In a preliminary numerical test, we directly interpolated the NSX-H atmosphere tables and examined the resulting specific intensities. We found that the boundary behavior as µ → 0 (i.e., near the stellar limb) is particularly susceptible to interpolation overshoot, occasionally producing unphysical negative intensities. This sensitivity arises from the steep angular dependence of the tabulated intensities close to the limb and the discrete structure of the NSX-H grids. To illustrate how such overshoot can affect otherwise typical observing configurations, we present two deliberately chosen test cases.

Both tests use a single circular hotspot with uniform temperature T = 10⁶ K, representing the simplest hotspot configuration. The hotspot has an angular radius of 0.01 rad (i.e., a compact, highly localized patch), ensuring that viewing-geometry effects dominate, minimizing complex-shape-induced systematics and avoiding dilution from extended or complex emission regions. In both cases, the neutron-star mass and radius are fixed at M = 1.4 M_⊙ and R = 12 km. To probe rapid rotation while remaining within the OS regime, we set the spin frequency to 600 Hz. The same version of the NSX-H atmosphere tables is used throughout, and the source distance is fixed at 150 pc. The colatitude of the hotspot center is θ_c = 1.983 rad for test 1 and θ_c = π − 0.01 rad for test 2; the observer colatitude is i = 0.01 rad for test 1 and i = 0.964 rad for test 2. A complete specification of both configurations is summarized in Table 3, and their geometries are illustrated in Figure 3.

The two tests differ only in viewing geometry and hotspot location. They are constructed to create more phases when the spot approaches or recedes beyond the visible limb, driving rays with µ → 0, where unconstrained cubic interpolation tends to overshoot. Although these examples intentionally maximize the visibility of the effect, the underlying issue is generic: for many reasonable geometries, a rotating hotspot will periodically sample small-µ rays (as it grazes the limb or re-enters the line of sight), making careful treatment near µ → 0 essential.

Importantly, these setups are not exotic. Test 1 is consistent with a southern magnetic pole in a canonical inclined dipole or an off-centered dipole configuration. Test 2 is likewise plausible for an off-centered dipole whose south pole lies near the rotation pole. Variants such as a narrow ring hotspot at the same colatitude (e.g., arising from higher multipolar fields) can be analyzed analogously and would exhibit similar sensitivity to the limb region.

In each test case, we evaluated four strategies for interpolating the NSX table. Strategy (a) applies cubic interpolation uniformly across all four dimensions. Strategy (b) also uses cubic Lagrange interpolation but enforces a non-negativity safeguard by clipping negative values to zero. Strategy (c) employs cubic interpolation uniformly in all dimensions with linear interpolation enforced at the µ → 0 boundary. Strategy (d), our baseline in this work, uses linear interpolation in two dimensions and cubic interpolation in the remaining two, while reverting the nominally cubic dimensions to linear near the µ → 0 boundary to suppress overshoot. All results are produced with the same CPU codebase and a high-resolution configuration.

Figure 4 presents phase–energy-resolved X-ray pulse profiles for test 1. For a small viewing angle (i = 0.01), the profile exhibits mild asymmetry over spin phase. Under uniform cubic interpolation (panel a), we observe pronounced negative intensities (large undershoots) centered near phases ∼0.25, 0.50, and 0.75, predominantly in the soft X-ray band. These unphysical negatives would severely corrupt the inferred pulse profile and any derived likelihood. For the X-PSI code (which we have not run ourselves), a preliminary replication (T. Salmi, priv. comm.) reports a smaller negative flux (occurring at a different location) than we find here. We have not independently verified the phase–energy locations of these occurrences and therefore do not draw conclusions about where in the profile they arise.

Introducing non-negativity clipping (panel b) removes the undershoots but creates extended zero-intensity “holes” that still end up biasing the likelihood computation by artificially suppressing flux. In contrast, uniform cubic interpolation with linear treatment at the edges (panel c) yields a qualitatively reasonable profile: by construction, linear schemes at the boundaries strongly limit overshoot and thus avoid negative intensities. Our baseline mixed scheme (panel d) likewise produces a physically consistent profile; the boundary reversion to linear in the nominally cubic dimensions effectively eliminates the negative regions, directly implicating cubic overshoot at the table edges of the µ axis as the root cause of these artifacts.

Test 2 (Figure 3) is intentionally close in observables to test 1, but it is physically distinct: the hotspot is fixed at the rotational south pole while the line of sight changes which longitudes are viewed as the star rotates. This geometry is well motivated by centered or off-centered dipolar fields heating neutron-star hotspots. The same qualitative failure modes reappear under uniform cubic interpolation, as shown in Figure 5: three phase-localized negative regions in the soft band, which become zeroed “holes” after clipping, leading to substantial underestimation of photon counts. As before, both the uniform linear scheme and our mixed baseline recover stable, physically plausible profiles for this configuration.

Crucially, these cases are not isolated to a single demonstrated geometry. Small perturbations of the hotspot location, changes in mass–radius, or alternative hotspot shapes readily reproduce these negative regions when unconstrained cubic interpolation is used near table boundaries. Because Bayesian inference explores the full parameter space of a given model, it is effectively inevitable that such configurations will be encountered during inference. Unless interpolation is handled with care, mass–radius constraints risk hidden systematics and hence reduced credibility. We therefore recommend that existing pipelines incorporate these two test cases into their validation suites and adopt atmosphere-table interpolation strategies that control boundary overshoot (e.g., uniform linear or mixed schemes with boundary linearization, as implemented here).

After this detailed discussion of the implementation of the atmosphere model, we proceeded to production-grade tests. Specifically, we reproduced end-to-end computation cases that incorporate all necessary ingredients, such as our baseline atmosphere interpolation scheme, interstellar absorption, and convolution with the instrument response. These tests demonstrate that our treatment of the forward computation of the pulse profile is implemented correctly. With these components in place, the codebase is prepared for inference while achieving between three and four orders of magnitude speedup in pulse profile synthesis relative to existing CPU- and Python-based codebases.

Fig. 3 Demonstration of two geometric test cases for atmosphere interpolation. Test 1: point-like hotspot on the southern hemisphere; the line of sight is aligned with the spin axis and lies above the rotational north pole. Test 2: point-like hotspot at the south pole; the observer is located in the northern hemisphere.

Fig. 4 Simulated phase–energy-resolved X-ray pulse profiles for interpolation test 1, computed using four interpolation schemes applied to the NSX tables: (a) cubic interpolation in all dimensions; (b) as in (a), with negative intensities clipped to zero; (c) cubic interpolation in all dimensions with linear interpolation enforced at table boundaries; (d) the baseline configuration used in this work, cubic along two axes and linear along the remaining two, with linear interpolation enforced at table boundaries.

Fig. 5 Simulated phase–energy–resolved X-ray pulse profiles for interpolation test 2, computed using four interpolation schemes applied to the NSX tables: (a) cubic interpolation in all dimensions; (b) as in (a), with negative intensities clipped to zero; (c) cubic interpolation in all dimensions with linear interpolation enforced at table boundaries; (d) the baseline configuration used in this work, cubic along two axes and linear along the remaining two, with linear interpolation enforced at table boundaries.

5 Production-level computation result

In this section, we proceed by starting with theoretical benchmarks, which omit realistic radiative physics and other effects. We move on to production-grade tests that connect directly to observational data, which is done by incorporating a realistic neutron-star atmosphere, ISM absorption, and NICER’s instrument response. In this section, we compare the results of the CPU and GPU versions of the codebase. At previous section, we are not discussing the computation performance specifically, because the most time-consuming part of this PPM is the realistic atmosphere model table interpolation, so vast computation improvement is based on improving this atmosphere treatment strategy. Here in this section, we will compare both the modeling accuracy and the computation performance, using the benchmarks provided in Choudhury et al. (2024b). As we discuss in the previous section, we are using different strategy on surface discretizations, instead of using the discretization of only the hotspot region, we consider the whole star surface discretization as default option, so when comparing the modeling accuracy of the surface elements in the hotspot, we using the active grid number in the hotspot as the comparison. The goal of this section is to demonstrate that the extreme geometry tests of Choudhury et al. (2024b), including cases that are challenging to model accurately at standard resolutions in existing codes such as X-PSI and IM, can be computed at our standard setup with discrepancies well below the Poisson fluctuation level, while achieving 10⁴–10⁵× speedups relative to the Ultra-resolution X-PSI runs. In this section, CPU tests were performed on an Intel i7-13700KF with single core and GPU tests on an NVIDIA GeForce RTX 4080 (16 GB).

Following Choudhury et al. (2024b), we adopt the same five “extreme” waveform tests to compare PPM treatments. Four tests use a single-temperature hot spot. All shapes are built only from circular components, either emitting or masking, because some code in that test limited to circular primitives. Masking components simply hide the flux from the region they cover. We excluded the double-temperature test case from Choudhury et al. (2024b) because the compactness parameter was set to an extreme value, rendering the emission spot always visible through gravitational light bending and permitting multiple imaging–an effect not treated by our present code version. Complete parameter setup and the exact geometries (rings, crescents, and the bithermal configuration) are tabulated and illustrated in Choudhury et al. (2024b, Table 1 and Figure 1).

All tests share the same non-hotspot settings: distance D = 150 pc; observer inclination i = 2.391 rad measured from the rotational north pole (a viewing geometry under which the near pole remains visible across rotational phase), spin frequency at infinity ν = 300 Hz, and neutral hydrogen column density N_H = 2 × 10¹⁹ cm⁻², and exposure time of 10⁶ s. The inclination is chosen to mirror the viewing of a well-studied nearby millisecond pulsar PSR J0437−4715, and the spin frequency is set at 300 Hz because the OS approximation is accurate to ≲0.1% at and below this rate (Bogdanov et al. 2019b), while discrepancies grow to a few percent by ∼600 Hz (Cadeau et al. 2007; Jakab & Morsink 2025).

The four single-temperature tests use M = 1.4 M_⊙ and equatorial radius R_e = 12.0 km. Within the single-temperature set, two tests form nonconcentric rings by nesting a mask inside an emitter (with sizes varied to stress overlap logic and discretization), and two form elongated crescents by letting a mask protrude into an emitter to accentuate edge effects.

Atmosphere setup is uniform across tests. Each emitting component is modeled as a geometrically thin, fully ionized hydrogen atmosphere; specific intensities are drawn from NSX-H lookup tables (Ho & Lai 2001). The resulting spectra are then attenuated using the same N_H above. For any instrument-response handling and the full per-component quantities, we refer the reader to Choudhury et al. (2024b), whose setup we follow here.

5.1 Resolution tiers setup

We evaluate two resolution tiers designed for distinct purposes: (i) a standard tier that is sufficiently accurate for routine inference at practical cost and (ii) a high-resolution tier intended to match the fidelity of the “Ultra-style” reference settings in (Choudhury et al. 2024b). To isolate the effect of surface discretization, we vary only the active number of surface patches within each hotspot. All other precision controlled parameters are held fixed at high values comparable to those used in rigorous code-verification runs. This ensures that numerical errors associated with phase and energy discretization, along with per-patch ray statistics, are negligible. The remaining lever on accuracy and speed is the surface discretization.

Table 4 summarizes the two tiers detailed resolution setup. The “high-resolution” tier uses an active 256 × 256 surface tiling per hotspot, which is comparable to the surface discretization commonly used in “Ultra-level” resolution tier (e.g., in X-PSI). The “standard” tier uses 128 × 128 patches per hotspot and is intended to keep numerical systematics below typical Poisson counting uncertainties for representative NICER-like exposures, while offering a substantial speedup. Because our codebase differs in architecture and algorithms from other frameworks, a one-to-one mapping of resolution settings is not straightforward; by “comparable” we mean that our tiers achieve similar accuracy to the production-grade benchmarks reported by Choudhury et al. (2024b).

In our test runs, both CPU and GPU implementations produce phase-resolved, energy-resolved waveforms under identical physics and numerical settings. For validation, we compare our waveforms to a trusted reference (eg. X-PSI Ultra-resolution setup in Choudhury et al. (2024b)) and require fractional residuals ≲10⁻³ (i.e., ≲0.1%) over phases where the flux is non-negligible. Small, localized deviations may occur near the photon number extremes. These phases carry very low counts and have negligible impact on likelihood evaluations used in pulse profile inference.

The reason we need to compare our CPU and GPU version code internally is because porting a ray-tracing pipeline from CPU to GPU is nontrivial. Memory layouts, parallel reduction order, and instruction-level differences can lead to small floating-point discrepancies even for the same numerical algorithm. While CPU implementations can reach either tier, end-to-end, high-resolution scans of the full parameter space are typically computationally prohibitive within reasonable computation time, whereas the GPU implementation enables production-grade inference at the high-resolution setting. The standard tier of GPU setup is designed to be extremely fast on GPU but already become untractable on CPU even with taking advantage of modern C++ coding, while keeping numerical systematics subdominant to Poisson fluctuation for all these extreme hot-spot setup.

To demonstrate the robustness of our GPU code on the whole possible hotspot geometry parameter space, in the following subsections, we evaluate four geometries that stress the ray-tracing and integration pipeline: Ring-Eq, Ring-Polar, Crescent-Eq, and Crescent-Polar. Their precise definitions (spot shapes, locations, and masks) follow Choudhury et al. (2024b, Table 1 and Figure 1). Each case is discussed in turn using identical numerical algorithm and precision-controlled parameters so that differences reflect only hotspot geometry or resolution. To ensure full transparency and reproducibility, we have released the source code at GitHub repository². A versioned snapshot, including all scripts required to regenerate the figures for both CPU and GPU implementations at both resolution tiers, has been archived on Zenodo³.

Fig. 6 Energy-summed waveforms (left) and phase-summed spectra (right) for the Ring-Eq case. The bottom subpanels show the fractional difference (in %) of each configuration relative to the X-PSI Ultra-resolution benchmark (black dotted). CPU results are plotted with solid lines: standard resolution (blue solid) and high resolution (red solid). GPU results are plotted with dashed lines: standard resolution (blue dashed) and high resolution (red dashed). In the left panel, the legend reports χ2 as defined in Equation (57), computed jointly over energy and phase channels as a quantitative measure of the discrepancy between full energy-resolved pulse profile. In the right panel, the legend lists the runtime for each configuration. Gray dashed lines in the bottom subpanels indicate the expected Poisson fluctuations of the Ultra-resolution X-PSI benchmark.

Fig. 7 Energy-summed waveforms (left) and phase-summed spectra (right) for the Ring-Polar case. The bottom subpanels show the fractional difference (in %) of each configuration relative to the X-PSI Ultra-resolution benchmark (black dotted). CPU results are plotted with solid lines: standard resolution (blue solid) and high resolution (red solid). GPU results are plotted with dashed lines: standard resolution (blue dashed) and high resolution (red dashed). In the left panel, the legend reports χ2 as defined in Equation (57) computed jointly over energy and phase channels as a quantitative measure of the discrepancy between full energy-resolved pulse profile. In the right panel, the legend lists the runtime for each configuration. Gray dashed lines in the bottom subpanels indicate the expected Poisson fluctuations of the X-PSI Ultra-resolution benchmark.

Fig. 8 Energy-summed waveforms (left) and phase-summed spectra (right) for the Ring-Polar case. The bottom subpanels show the fractional difference (in %) of each configuration relative to the X-PSI Ultra-resolution benchmark (black dotted). CPU results are plotted with solid lines: standard resolution (blue solid) and high resolution (red solid). GPU results are plotted with dashed lines: standard resolution (blue dashed) and high resolution (red dashed). In the left panel, the legend reports χ2 as defined in Equation (57) computed jointly over energy and phase channels as a quantitative measure of the discrepancy between full energy-resolved pulse profile. In the right panel, the legend lists the runtime for each configuration. Gray dashed lines in the bottom subpanels indicate the expected Poisson fluctuations of the X-PSI Ultra-resolution benchmark.

Fig. 9 Energy-summed waveforms (left) and phase-summed spectra (right) for the Ring-Polar case. The bottom subpanels show the fractional difference (in %) of each configuration relative to the X-PSI Ultra-resolution benchmark (black dotted). CPU results are plotted with solid lines: standard resolution (blue solid) and high resolution (red solid). GPU results are plotted with dashed lines: standard resolution (blue dashed) and high resolution (red dashed). In the left panel, the legend reports χ2 as defined in Equation (57) computed jointly over energy and phase channels as a quantitative measure of the discrepancy between full energy-resolved pulse profile. In the right panel, the legend lists the runtime for each configuration. Gray dashed lines in the bottom subpanels indicate the expected Poisson fluctuations of the X-PSI Ultra-resolution benchmark.

5.2 Ring-Eq case

The Ring-Eq configuration is a ring-shaped hotspot centered on the stellar equator. It typically produces a single predominant pulse per rotation. Figure 6 compares energy-summed waveforms and phase-summed spectra across our two resolution tiers (standard and high) and across CPU and GPU backends, bench-marked against the X-PSI Ultra-resolution reference used in Choudhury et al. (2024b). The X-PSI Ultra-resolution configuration serves here (and throughout) as a high-fidelity standard for extreme geometries.

To quantify the agreement, we can compute a joint, 2D χ² statistic, $${\chi }^2=K{\displaystyle \sum _{i,j}\frac{(m_{ij}-d_{ij})}{m_{ij}},}$$(57)

where m_ij and d_ij denote model and reference counts in energy bin i and phase bin j, respectively. The normalization factor K is defined as $$K=\frac{{10}^6}{{\displaystyle \sum ij{m}_{ij}}}.$$(58)

The normalization constant is included to enable fair comparisons between different waveforms. Without it, Equation (57) implies that the overall χ² would scale with the total photon counts. Accordingly, we normalize all computations to a fixed total of 10⁶ counts, as suggested in Choudhury et al. (2024b). Under Poisson statistics and in the large-count limit, −2 ∆ ln L ≈ χ², so maintaining small χ² controls the corresponding change in log-likelihood ($\Delta \ln L\approx -\frac{1}{2}{\chi }^2$). Our goal is to demonstrate that both CPU and GPU implementations meet production-grade accuracy for pulse profile inference.

The energy-summed waveforms from all four runs (standard+high resolution × CPU+GPU) agree with the benchmark to within the expected Poisson fluctuations. Differences cluster by resolution rather than by backend: increasing the surface-tiling density (from 128 × 128 to 256 × 256 patches per hotspot) uniformly reduces fractional residuals, while switching from CPU to GPU at fixed resolution has no visible effect on accuracy. Small deviations concentrate near flux minima (when the ring rotates partially out of view), where counts are lowest. These phases carry negligible statistical weight in likelihood evaluations.

The phase-summed spectra exhibit the same pattern. Curves at a given resolution overlap closely regardless of backend, and the high-resolution tier yields uniformly smaller absolute fractional differences. This behavior is consistent with surface discretization being the dominant residual error once the per-patch rays, phase sampling, and energy sampling are held fixed at high settings.

Runtime measurements are shown in the legend of Figure 6. Within our codebase, increasing the surface-patch count by a factor of four (standard → high resolution) scales wall time by approximately the same factor (e.g., 5.63 s → 22.4 s on CPU), as expected for a largely linear compute-to-patch workload. Even with careful modern C++ implementation, exhaustive parameter-space exploration at high resolution is computationally prohibitive on CPU, whereas the GPU implementation achieves ∼2000× speedups over our CPU runs at identical settings, enabling production-grade inference even under the high-resolution setup. For context, the quoted X-PSI Ultra-resolution runtime for comparable resolution is ∼1–5 minutes; our GPU timings are 4–5 orders of magnitude faster while maintaining agreement at the χ² ∼ 10⁻² level.

5.3 Ring-Polar case

The Ring-Polar configuration is a small ring-shaped hotspot centered near the rotational pole. Because the emitting area is small, a spot at the same surface temperature produces fewer total counts than larger hotspot. Moreover, for typical neutron-star compactness (e.g., M ≈ 1.4 M_⊙, R_e ∼ 12 km), gravitational light bending keeps most of a polar ring at least partially in view over the rotation cycle for a wide range of observer inclinations, so the modulation amplitude is modest rather than strongly pulsed.

Figure 7 compares energy-summed waveforms and phase-summed spectra across our two resolution tiers (standard and high) and across CPU and GPU backends, using the X-PSI Ultra-resolution configuration as the high-fidelity reference. Agreement is quantified by the joint, 2D statistic χ². In our runs, energy-summed waveforms from all configurations are consistent with the benchmark at the Poisson level: the standard tier yields χ² ≈ 1.5 for both CPU and GPU, while the high-resolution tier reduces this to χ² ≈ 0.05. As in the other cases, differences group by resolution rather than by backend, indicating that surface discretization, not CPU vs. GPU, is the dominant factor governing residuals once other precision settings are held fixed.

The phase-summed spectra show the same pattern. At energies where the counts are highest, curves overlap closely; in the highest-energy channels (e.g., channel indices ≳ 250, where counts are lowest), fractional residuals can reach the ∼1–1.5% level without changing the total χ². This behavior is expected: low-count bins carry little statistical weight, so localized percent-level differences there have negligible impact on the likelihood computation.

Within a given backend, as in the Ring-Eq case, increasing the surface-patch count from 128 × 128 to 256 × 256 (a ×4 increase) scales the Ring-Polar runtime by approximately the same factor, as expected for a patch-parallel workload. Absolute runtime for the Ring-Polar case is comparable to, or slightly longer than, that for the equatorial ring. While the CPU implementation meets the accuracy target, exhaustive parameter scans at either resolution tier are computationally prohibitive on CPU within practical budgets. The GPU implementation achieves ∼10³× speedups at identical settings, enabling production-grade Bayesian inference at very high fidelity.

5.4 Crescent-Eq case

In the Crescent-Eq configuration, the hotspot is an elongated crescent located near the stellar equator. Because the crescent spans a substantial azimuthal extent, it remains visible over a large fraction of the rotational phase and yields comparatively high counts. As noted by Choudhury et al. (2024b), this geometry can be slightly more sensitive to code implementations and resolution choices, with percent-level differences sometimes appearing at phases where the flux approaches zero; such localized discrepancies are expected when counts are minimal and carry little statistical weight.

In our experiments (Figure 8), the energy-summed waveforms agree with the X-PSI Ultra-resolution benchmark at the ≲0.1% level across phase, and increasing the surface-patch density from the standard tier (128 × 128 patches per hotspot) to the high-resolution tier (256 × 256) produces no visually discernible change. The joint 2D statistic χ² defined earlier likewise remains very small: at the standard tier we obtain χ² ≈ 0.005, while at the high-resolution tier we find χ² ≈ 0.006. Differences of this magnitude are consistent with floating-point and discretization effects rather than systematic modeling errors; at such small values, monotonic improvements in χ² with resolution are not guaranteed.

The phase-summed spectra show a similar pattern. Curves overlap closely in energy ranges with substantial counts, while the highest-energy bins can exhibit slightly larger fractional residuals without materially affecting the total χ². Small backend-dependent trends (CPU vs. GPU) in these highest-energy bins are consistent with floating-point accumulation effects. At fixed resolution, results group by resolution rather than by hardware.

We attribute the robustness in this geometry, in part, to our full-surface discretization strategy: the entire stellar surface is finely meshed, and the crescent is selected by masking within that mesh. This approach represents sharp crescent edges cleanly across phase, reducing edge-alignment artifacts that can arise when only the hotspot boundary is discretized. While alternative pipelines employ different but valid tiling schemes (e.g., Bogdanov et al. 2019b, 2021), our results indicate that a full-surface mesh provides stable accuracy for elongated, edge-dominated spots.

Finally, the runtime remains highly favorable. On the GPU, the standard tier computes a single energy-resolved profile in ∼3.4 ms, enabling production-grade scans. On the CPU, the same configuration requires ∼7.16 s at the standard tier and ∼27 s at the high-resolution tier. Thus, although the CPU implementation attains the desired accuracy, only the GPU implementation provides the required speedup to support large-scale Bayesian inference.

The crescent geometry is particularly relevant for future studies incorporating higher-order magnetic multipoles, which can naturally yield noncircular, elongated hotspots. For example, Huang & Chen (2025) show that, in an off-centered dipole magnetosphere, incorporating a volumetric return-current region modifies the polar-cap temperature map, adding a crescent-like, lower-temperature component. These effects are modest for a centered (aligned) dipole but become more pronounced as the dipole offset increases. Thus, testing this crescent configuration under different placements is extremely valuable for future physics-motivated, hotspot-based PPM.

5.5 Crescent-Polar case

In the Crescent-Polar configuration, we place the same elongated hotspot near the rotational north pole and adopt a near-limb viewing geometry, such that only a small portion of the emitting region is directly visible at any given phase. In this edge-dominated regime, most detected photons originate from a narrow strip along the spot boundary where projected area, light bending, and atmospheric beaming combine to produce a steep intensity gradient. Consequently, the energy-summed waveform becomes highly sensitive to numerical resolution specifically at the spot edge. Across different codes, and even within the same code at different resolution settings, this geometry can produce visibly different waveforms. Although the largest fractional differences tend to occur in low-count bins, the associated mismatch metric, χ² (as defined in Equation (1) of Choudhury et al. (2024b) and here computed jointly over energy and phase channels), can still become large enough to change the likelihood by a factor exp(−χ²/2), potentially biasing inference if not controlled.

Figure 9 shows that our implementation closely matches the benchmark in this challenging configuration. Using both the CPU and GPU pipelines at the standard resolution tier yields χ² ≈ 1.96. In typical PPM analyses, a forward-model mismatch of order unity is negligible relative to the χ² variations that drive parameter updates and therefore has no visible impact on the inference. At the high-resolution tier, the discrepancy remains at χ² ≃ 1.82, indicating only marginal improvement in precision. The right panel shows that the spectral fractional differences decrease at high resolution, this is the primary source of the modest improvement from the standard to the high-resolution tier. By contrast, the left panel exhibits a persistent pattern in the fractional differences that generally stays below the Poisson fluctuation level; however, near half phase (the trough between the two peaks), the residuals reach up to ∼5%. This phase interval corresponds to when the crescent-shaped polar hotspot becomes faint and approaches the limb (µ → 0). As discussed in Section 4, cubic-interpolation overshoot in this regime, followed by the non-negativity clipping safeguard in X-PSI, leads to an underestimation of the local photon counts. While the clipping prevents the systematic from exceeding the Poisson fluctuation level, it nevertheless produces a real, phase-localized bias that is evident here and results in an overall underestimation of the counts. Comparing with Figure 7 of Choudhury et al. (2024b), the IM group’s results oscillate between under- and overestimating the photon counts relative to the X-PSI benchmark, suggesting that their atmosphere-model interpolation is sensitive to interpolation choices. In contrast, the Alberta codebase, which uses linear interpolation for the atmosphere model, is more robust at boundaries and does not suffer from overshoot-induced underestimation. Accordingly, it yields consistently higher counts at the limb. These comparisons reinforce that careful treatment of atmosphere-table interpolation is essential for reliable pulse profile computations. Consistent with our other tests, the phase-summed spectra indicate that the residual fractional differences are dominated by high-energy channels with very low counts.

From a performance standpoint, our code remains highly efficient in this case. On the CPU, the standard setup runs in ∼1.63 s and the high-resolution setup in ∼6.63 s. On the GPU, the corresponding run times are ∼1.77 ms and ∼6.9 ms, respectively, yielding speedups of order 10³ relative to our own CPU implementation for both resolution tiers. These gains are achieved without emulators or machine-learning surrogates; we use high-fidelity forward calculations with algorithmic refinements executed on the GPU.

This case underscores an important methodological point. Because edge-on, elongated spots can exist in physically motivated scenarios (e.g., modestly off-centered dipoles can place volumetric hotspot regions near a rotational pole; see Figures 3a–b of Huang & Chen 2025), there is no compelling physical argument to a priori forbid such geometries. Reliable inference therefore requires numerical treatments that remain accurate in edge-dominated regimes. By breaking the traditional speed-versus-accuracy trade-off through a GPU implementation, our code provides one practical pathway to control χ² to the ∼1 level even in this extreme configuration, while simultaneously achieving ∼10³–10⁴ speedups.

Our CPU implementation uses IEEE 754 double-precision floating-point numbers (64-bit, FP64), whereas the current GPU implementation on an RTX 4080 employs single precision (32-bit, FP32) for all computations. Because our forward model involves no long time integration, cumulative round-off error is limited; moreover, we adopt numerically stable formulations to control FP32 error. In the current benchmarks at matched resolution, the FP32 GPU results are indistinguishable from the FP64 CPU results, indicating that FP32 is adequate for our present configuration. From an engineering perspective, enabling FP64 in our GPU codebase is straightforward. Users who run this code on scientific GPU clusters may optionally switch to FP64 to improve numerical precision.

6 Discussion and conclusion

This work presents a GPU-accelerated X-ray pulse profile modeling (PPM) codebase that directly addresses the typical accuracy-speed trade-off in current analyses. We first outline the theoretical framework of pulse profile forward modeling, rigorously verifying and deriving all formulas essential to PPM forward computation (with more details in Appendices A–C). We then provide a detailed computational recipe for both CPU and GPU implementations and clearly document the porting of the CPU code to the GPU to maximize reproducibility.

We verified our implementation by first reproducing the theoretical light-curve benchmarks of Bogdanov et al. (2019b), demonstrating accurate computation of light curves for specified spacetime and hotspot geometries under the oblate-Schwarzschild (OS) approximation. The Schwarzschild + Doppler (S+D) approximation is also supported, although OS remains our default production setting.

We further examined a commonly ignored source of systematic error that arises when moving from theoretical benchmarks to models with realistic atmospheres. Because of the current structure of the widely used NSX model, several codebases interpolated within these precomputed neutron-star atmosphere lookup tables using high-order polynomial interpolation. Using two representative viewing-geometry test cases, we demonstrate that cubic Lagrange interpolation can produce significant overshoot, yielding unphysical negative specific intensities. One approach to mitigation is to clip these negative values to zero. However, in the test cases presented here, both strategies (cubic Lagrange interpolation with or without clipping) can lead to underestimation of the total photon flux at many rotational phases, particularly in the soft X-ray band, depending on implementation. Preliminary X-PSI tests indicate smaller negative flux values (T. Salmi, priv. comm.) and the phase–energy locations of any deficits may differ across implementations.

Moreover, these new test cases are neither extreme nor exotic. They plausibly represent the emission geometry of a southern polar cap arising from either a canonical centered dipole or an off-centered dipole configuration. The region of parameter space that exhibits this overshoot behavior is potentially broad; modest variations of the pulse profile model parameters readily produce similar behavior across different settings. In fact, for many common viewing and hotspot geometries, rotation will bring portions of the hotspot into and out of the line of sight near the limb, where these interpolation artifacts are most likely to arise. The interpolation artifact identified here suggests that previously reported PPM inferences that relied on high-order polynomial interpolation without boundary-aware treatment could contain an uncontrollable systematic error. Reexamining those analyses with a boundary-aware interpolation scheme, like the one proposed in this work, would be extremely valuable for assessing the impact on inferred parameters, especially M and R.

To mitigate interpolation artifacts, we adopted a mixed-order interpolation scheme. Within the interior of the atmospheric lookup tables, we used cubic Lagrange interpolation; but at the grazing-emission boundary, we reverted to linear interpolation to suppress overshoot. We also allow the interpolation order to differ across table dimensions, mixing linear and cubic along different axes to reduce computational cost while maintaining accuracy. The optimal configuration is table- and physics-dependent; constructing denser atmospheric grids would enable a more definitive validation of the interpolation scheme. A comprehensive evaluation is deferred to future work.

Building on these theoretical checks and the careful treatment of the atmosphere model, we turned to production-grade benchmarks, which include all of the necessary ingredients for PPM forward computation. We reproduced the suite of benchmarks in Choudhury et al. (2024b), incorporating a realistic atmospheric treatment, interstellar absorption, and the NICER instrument response. We implemented both CPU and GPU codes and examined two resolution configurations that differ primarily in the number of surface grid elements within the hotspot. Holding all other setup choices fixed allowed us to isolate how discretization influences runtime scaling while maintaining comparable physical fidelity.

The observed GPU performance indicates clear implications for inference, enabling higher-resolution forward models, broader parameter exploration, and more robust convergence diagnostics within modest computational budgets. With a modern C++ implementation, CPU runtimes at these extremely high resolutions, ranging from several seconds to tens of seconds per profile, remain prohibitive for full Bayesian analyses and robustness testing. In contrast, the GPU implementation presented here reduces per-profile computation time to the millisecond regime for the standard resolution and to ∼10 ms for the extremely high-resolution setting on a single RTX 4080 GPU. Because fully converged posterior exploration can require 10⁷−10⁸ likelihood evaluations (depending on the problem’s dimensionality), these timings make end-to-end inference feasible on a single GPU or modest GPU clusters.

Crucially, the speed improvement does not compromise accuracy. At the fine resolutions adopted here, our results closely match published benchmarks while retaining millisecond-scale runtimes. Even for geometries previously viewed as resolution-sensitive, our high-resolution GPU runs reproduce reference benchmarks, effectively removing the hard accuracy-speed bottleneck.

Beyond the common single- or two-temperature prescriptions used in many current inferences, our framework accepts arbitrary, user-specified global temperature maps on the stellar surface. This capability enables inference based on physics-motivated hotspots, such as those predicted by off-centered dipole or more general multipolar magnetic configurations (e.g., Huang & Chen 2025; Gralla et al. 2017; Lockhart et al. 2019), and naturally extends to accretion-powered millisecond pulsars, where global temperature distributions are more appropriate than bounded caps (Das et al. 2022; Das et al. 2025). Because purely geometric hotspot choices can introduce degeneracies and ad hoc modeling systematics, we advocate reducing arbitrary geometry in favor of models tied to the magnetospheric structure of the neutron star (e.g., controlled multipolar expansions) and potentially enabling multiwavelength constraints. Related efforts (e.g., Kalapotharakos et al. 2021; Olmschenk et al. 2025; Pétri et al. 2025) typically assume uniform-temperature caps; observational and theoretical indications suggest, however, that polar-cap regions are structured, with nonuniform temperatures and subregions dominating the X-ray emission (Spitkovsky 2006; Bai & Spitkovsky 2010; Gralla et al. 2017; Lockhart et al. 2019). Our code is designed to accommodate such more complex hotspot models.

While the present work emphasizes fast and accurate forward modeling, the inference layer is a natural next step. Our GPU optimization extends beyond waveform synthesis to the entire likelihood evaluation: phase- and energy-resolved model counts, background treatment, and Poisson log-likelihood accumulation run on GPU and parallelize naturally across energies, phase bins on different parameters. We support two complementary paths for sampling: (i) a thin Python/C++ interface that interoperates with established CPU-side samplers; for instance, MULTINEST (Feroz et al. 2009), ULTRANEST (Buchner 2021), and DYNESTY (Speagle 2020). We also kept the likelihood kernel resident on the GPU to minimize host–device overhead; and (ii) a fully GPU-resident parallel-tempered MCMC (PT-MCMC); and, resources permitting, a GPU-friendly nested-sampling engine) to handle multimodal posteriors and strong parameter correlations. Both modes can scale to multiple GPUs by sharding chains or live points across GPUs and batching likelihood calls within each device. Using this architecture together with physics-motivated temperature maps, we will reanalyze NICER sources, systematically varying samplers and settings to assess robustness and to deliver stable, reproducible constraints on neutron-star masses and radii.

Before proceeding to inference, we note an additional physical effect: multiple imaging due to strong-field gravitational light bending. For sufficiently compact stars, approximately when GM/(Rc²) ≳ 0.28, photon trajectories from a given surface element can reach the observer along multiple paths, yielding secondary images. Our current framework is not yet equipped to model these higher-order images and, thus, implementing this capability for production-grade inference is deferred to future work. For context, current NICER constraints for PSR J0740+6620 imply a central-value compactness likely below this threshold, although portions of the allowed parameter space may reach it; in those more compact configurations, this multiple-imaging effect could become relevant.

In summary, the codebase presented here provides a modern computational framework for PPM: it achieves production-grade accuracy at millisecond-scale runtimes. Our results are aligned with the majority of previous community tests and we propose two new atmosphere-interpolation tests. Our work substantially breaks the accuracy-speed bottleneck of Bayesian PPM studies. In parallel with advancing reanalyses of existing observational sources with physics-motivated hotspot models, we invite the community to use our GPU PPM codebase to reexamine the robustness of previous studies under diverse inference setups. For transparency and pedagogical reasons, we also provide a CPU reference implementation with a clear structure and high readability. Taken together, these efforts advance pulse profile modeling toward a more efficient, reliable, and reproducible regime.

Data availability

GPU-PPM Github repository: https://github.com/zhoutz/gpu_ppm. Zenodo repository of reproduction file of this work: https://doi.org/10.5281/zenodo.17315493.

Acknowledgements

H.C. acknowledges support from NSF Grant AST-2308111. The authors extend gratitude to Sharon Morsink, Tuomo Salmi for providing valuable Benchmark, help on clarifying the theoretical details about PPM, Alexander Chen about the GPU implementation suggestions and H.C. would like to thank Anna Watts, Alexander Dittmann, Cole Miller, Roger Romani, Wenbin Lu, Roger Blanford for insightful discussion and comments.

References

Appendix A Relativistic Doppler effect

To account for the relativistic Doppler effect, we consider two inertial reference frames: the lab frame, in which the stellar surface element is in motion, and the comoving frame, in which the element is at rest. Throughout this section, we set light speed c = 1. A prime (′) denotes quantities measured in the comoving frame. The goal is to relate (i) the angle between the photon direction and the surface normal, and (ii) the solid-angle element, as seen in the two frames. These relations are what enter the specific-intensity transformations.

Suppose that, in the lab frame, a small emitting surface patch translates along the +ŷ direction with speed v, while its outward surface normal is aligned with $\widehat{z}$. A photon leaves the patch in direction $\widehat{k}$. We project $\widehat{k}$ onto the x-z plane. Define k_⊥ in x-z plane, $$k_{\perp }=\widehat{k}-\widehat{k}\cdot \widehat{y}{\widehat{k}}_{\perp }=\frac{k_{\perp }}{\left|k_{\perp }\right|}$$(A.1)

We then parametrize the direction by a polar angle ξ measured from +ŷ direction and an azimuth ϕ measured in the x-z plane. With this choice, $$\widehat{k}\cdot \widehat{y}=\cos \xi ,{\widehat{k}}_{\perp }\cdot \widehat{x}=\cos \varphi ,{\widehat{k}}_{\perp }\cdot \widehat{z}=\sin \varphi ,$$(A.2)

so that the unit propagation vector has the component form $$\widehat{k}=(\sin \xi \cos \varphi ,\cos \xi ,\sin \xi \sin \varphi ).$$(A.3)

Along the null ray, the spatial coordinates advance at light speed in the direction $\widehat{k}$, hence $$x=t\sin \xi \cos \varphi ,y=t\cos \xi ,z=t\sin \xi \sin \varphi ,$$(A.4)

We now boost to the comoving frame with a standard Lorentz transformation along +ŷ: $$t^{\prime }=\gamma (t-vy),y^{\prime }=\gamma (y-vt),x^{\prime }=x,z^{\prime }=z,$$(A.5)

where the Lorentz factor is $$\gamma =\frac{1}{\sqrt{1-v^2}}.$$(A.6)

Next, we introduce the angle σ between the photon direction and the surface normal $\widehat{z}$. This is the local emission angle relevant for atmosphere beaming, in our geometry, $$\cos \sigma =\frac{z}{\sqrt{x^2+y^2+z^2}}=\sin \xi \sin \varphi ,$$(A.7)

and the same definition applied in the comoving frame gives σ′. Using the boosted coordinates and the null condition, we obtain $$\begin{array}{l}\cos {\sigma }^{\prime }=\frac{z^{\prime }}{\sqrt{{x^{\prime }}^2+{y^{\prime }}^2+{z^{\prime }}^2}}\\ =\frac{\sin \xi \sin \varphi }{\sqrt{{(\sin \xi \cos \varphi )}^2+{[\gamma (y/t-v)]}^2+{(\sin \xi \sin \varphi )}^2}}\\ =\frac{\cos \sigma }{\sqrt{{\sin }^2\xi +{\gamma }^2{(\cos \xi -v)}^2}}\\ =\frac{\cos \sigma }{\gamma \sqrt{(1-v^2){\sin }^2\xi +{\cos }^2\xi -2v\cos \xi +v^2}}\\ =\frac{\cos \sigma }{\gamma \sqrt{1-2v\cos \xi +v^2{\cos }^2\xi }}\\ =\frac{\cos \sigma }{\gamma (1-v\cos \xi )}\\ =\delta \cos \sigma ,\end{array}$$(A.8)

where we have identified the usual Doppler factor $$\delta =\frac{1}{\gamma (1-v\cos \xi )}.$$(A.9)

This expression shows explicitly how the emission angle with respect to the surface normal is modulated by the boost.

We now turn to the transformation of the solid angle. With our spherical coordinates we define $$d\text{Ω}=d\cos \xi d\varphi ,d{\text{Ω}}^{\prime }=d\cos {\xi }^{\prime }d{\varphi }^{\prime }$$(A.10)

Because the Lorentz boost is along y and leaves x and z unchanged, the azimuth in the x-z plane is unaffected: $${\varphi }^{\prime }=\text{arctan}\left(\frac{z^{\prime }}{x^{\prime }}\right)=\text{arctan}\left(\frac{z}{x}\right)=\varphi ,$$(A.11)

so that $$d\varphi =d{\varphi }^{\prime }.$$(A.12)

The polar cosine with respect to the boost axis is, by definition, $$\cos \xi =\frac{y}{\sqrt{x^2+y^2+z^2}},\text{a}$$(A.13)

and its correspondence in comoving frame is $$\cos {\xi }^{\prime }=\frac{y^{\prime }}{\sqrt{x^{\prime 2}+y^{\prime 2}+z^{\prime 2}}}=\frac{\gamma (y/t-v)}{\gamma (1-v\cos \xi )}=\frac{\cos \xi -v}{1-v\cos \xi }.$$(A.14)

This is the standard formula for angles measured relative to the boost direction. Differentiating with respect to cos ξ^′ gives the Jacobian of the transformation: $$\begin{array}{l}d\cos {\xi }^{\prime }=\frac{d\cos \xi (1-v\cos \xi )-(\cos \xi -v)(-vd\cos \xi )}{{(1-v\cos \xi )}^2}\\ =d\cos \xi \frac{1-v^2}{{(1-v\cos \xi )}^2}\\ ={\delta }^{2}d\cos \xi .\end{array}$$(A.15)

x^′2 + y^′2 + z^′2

which immediately implies the solid-angle relation $$d{\text{Ω}}^{\prime }={\delta }^{2}d\text{Ω}.$$(A.16)

Appendix B Time dilation approximation

In this section, we derive the mapping between emission time at a rotating stellar surface element and the reception time at a distant detector. The goal is to obtain a compact differential relation dt_emit ↔ dt_receive

Suppose a nongravitating spherical star which center locate at (x, y, z) = (0, 0, 0) with mass M and equatorial radius R_e, a distant observer is located in the +x direction. The location of a small beaming patch on the star surface while rotating we denote as r(t_emit)(all t are measured in observer’s frame). A photon emit from r, that finally hit the plane x = D, where D is a large number represent the star-observer distance, will undergo 2 stages.

The 1st stage is, from the photon been emitted, to it arrive the plane x = R, we use ∆t₁ to denote the the time elapsed in stage 1, $$t_{\text{hitR}}=t_{\text{emit}}+\text{Δ}{t}_1.$$(B.1)

Here, x = R is a fiducial plane placed beyond the stellar surface along $+\widehat{x}$.

The 2nd stage begins when photon hit x = R and ends when photon hit plane x = D, $$t_{\text{receive}}=t_{\text{hitR}}+\text{Δ}{t}_2,\text{Δ}{t}_2=\frac{D-R}{c}.$$(B.2)

Because the ray is taken to propagate parallel to $+\widehat{x}$ once it leaves the near zone, ∆t₂ is a constant independent of emission phase. This cleanly separates phase-dependent (Stage 1) from phase-independent (Stage 2) contributions.

In the 1st stage, we consider 2 consecutive emitted photon from the same patch, photon A and photon B, with emitting time relation $$t_{\text{emit}}^B=t_{\text{emit}}^A+d{t}_{\text{emit}},$$(B.3)

since the star is rotating, these 2 photon are not emit at the same location, so we have $\Delta t_1^A\ne \Delta t_1^B$. Since we are treating the non-gravitating case, we ignore the gravitational light bending effect, we assume the photon go straight line towards +x direction. From geometry we can derive $$\text{Δ}{t}_1=\frac{R-\widehat{x}\cdot r(t_{\text{emit}})}{c},$$(B.4)

so $$\frac{d(\text{Δ}{t}_1)}{d{t}_{\text{emit}}}=\frac{-\widehat{x}\cdot \dot{r}(t_{\text{emit}})}{c}=-\beta \cos \xi .$$(B.5)

where $\cos \xi =\widehat{x}\cdot \widehat{v}$ is the cosine of the angle between patch velocity and the photon emit direction, and $$\beta =\frac{v}{c}=\frac{\left|\dot{r}\right|}{c}.$$(B.6)

In stage 1, we get $$\frac{d{t}_{\text{hitR}}}{d{t}_{\text{emit}}}=1+\frac{d(\text{Δ}{t}_1)}{d{t}_{\text{emit}}}=1-\beta \cos \xi ,$$(B.7)

which expresses the local time-compression/expansion between emission and crossing of the x = R plane.

In stage 2, ∆t₂ is constant, so $$\frac{d{t}_{\text{receive}}}{d{t}_{\text{emit}}}=\frac{d{t}_{\text{hitR}}}{d{t}_{\text{emit}}}=1-\beta \cos \xi ,$$(B.8)

and the total mapping from emission to reception, in the non-gravitating straight-line limit, is purely kinematic. When the surface moves toward the observer (cos ξ > 0) the received cadence is compressed, when it recedes, it is dilated.

When considering a Schwarzschild gravitating star, the surface of star and the observer have different gravitational potential, so we have to include the gravitational redshift effect, which is $$\frac{d{t}_{\text{receive}}}{d{t}_{\text{emit}}}=\frac{1}{\sqrt{1-R_{\text{S}}/R}},$$(B.9)

where R_S is the Schwarzschild radius associated with M and R.

Combine both, we have $$d{t}_{\text{emit}}=\frac{\sqrt{1-R_{\text{S}}/R}}{1-\beta \cos \xi }d{t}_{\text{receive}},$$(B.10)

which is the working relation used in the main text. The numerator represents gravitational time dilation, and the denominator represents kinematic compression/expansion from motion along the line of sight. In the limiting cases β → 0 or R_S/R → 0, Equation (B.10) correctly reduces to purely gravitational or purely kinematic effects, respectively. If both vanish, dt_emit = dt_receive as expected.

Appendix C Blackbody radiation and flux formula

We begin with the Planck’s law of blackbody radiation, $$B(\nu ,T)=\frac{2h{\nu }^3}{c^2}\frac{1}{e^{h\nu /kT}-1},$$(C.1)

where B(ν, T) is the spectral radiance of a blackbody in thermal equilibrium at temperature T, with unit power per (unit solid angle × unit area normal to the propagation × unit frequency). Express it in more precise form, $$\frac{EdN}{dtd\text{Ω}\cos \sigma dSd\nu }=\frac{2h{\nu }^3}{c^2}\frac{1}{e^{h\nu /kT}-1},$$(C.2)

where E is the photon energy, dN the photon count, dΩ the differential solid angle around the propagation direction, cos σ dS the projected surface area element, and dν the frequency bin. Writing Planck’s law in this fully differential form will let us transform variables and frames cleanly.

The advantage of using this form is that we can easily change variables, suppose we want to change ν to E with $$E=h\nu ,dE=hd\nu ,$$(C.3)

we directly plug them in and get $$\frac{EdN}{dtd\text{Ω}\cos \sigma dSdE}=\frac{2E^3}{h^3{c}^2}\frac{1}{e^{E/kT}-1}.$$(C.4)

This recasts the blackbody intensity in per-energy units, which is convenient for applications on detector responses that are tabulated versus E.

We use prime(′) to denote quantities measured in the comoving frame of the emitting patch. For a small patch on the star surface emitting blackbody radiation, we rewrite the above equation in comoving frame, $$\frac{E^{\prime }d{N}^{\prime }}{d{t}^{\prime }d{\text{Ω}}^{\prime }\cos {\sigma }^{\prime }d{S}^{\prime }d{E}^{\prime }}=I_{bb}\equiv \frac{2{E^{\prime }}^3}{h^3{c}^2}\frac{1}{e^{E^{\prime }/k{T}^{\prime }}-1}.$$(C.5)

Here I_bb is the comoving specific intensity of a blackbody. The primes emphasize that both thermodynamic variables (T′) and emission angles (σ′) are defined in the instantaneous rest frame of the surface element. Our task below is to relate these comoving frame differentials to the ones measured by a distant observer.

We now calculate each term in the observer’s frame. According to Doppler effect and gravatitional redshift, $$E=\delta \sqrt{1-u}{E}^{\prime }.$$(C.6)

The observed photon energy is boosted by the special-relativistic Doppler factor δ and reduced by the gravitational redshift factor $\sqrt{1-u}$, with u ≡ R_S/R.

According to number conservation, each emit photon should arrive the observer, $$dN=d{N}^{\prime }.$$(C.7)

According to Equation (B.10), $$d{t}^{\prime }=\frac{\sqrt{1-u}}{1-\beta \cos \xi }dt=\sqrt{1-u}\delta \gamma dt.$$(C.8)

The comoving time element dt′ relates to the observer frame time interval dt by the same factors derived in Appendix B: gravitational time dilation (numerator) and kinematic arrival-time compression/expansion (denominator).

According to Equation (A.16), $$\begin{array}{l}d{\text{Ω}}^{\prime }={\delta }^{2}d\text{Ω}={\delta }^{2}d\cos \alpha d\lambda \hfill \\ ={\delta }^2\frac{d\cos \alpha }{d\cos \psi }d\cos \psi d\lambda ={\delta }^2\frac{d\cos \alpha }{d\cos \psi }\frac{dA}{D^2}.\hfill \end{array}$$(C.9)

where λ is the azimuthal angle around the outward radial vector, dA is the differential area of the detector, and D is the distance from star to the observer.

The first equality is the standard formula dΩ′ = δ²dΩ from Appendix A. We then (i) resolve dΩ into polar and azimuthal differentials around radial vector, (ii) introduce the gravitational light-bending Jacobian d cos α/d cos ψ that maps the apparent angle ψ at infinity to local emission angle α of spherical surface, and (iii) relate d cos ψdλ to detector’s collecting area dA via d cos ψdλ = dA/D².

Combine all these together, we have $$\frac{EdN}{dtdAdE}=\sqrt{1-u}{I}_{bb}{\delta }^3\gamma \cos {\sigma }^{\prime }\frac{d\cos \alpha }{d\cos \psi }\frac{d{S}^{\prime }}{D^2}.$$(C.10)

This is the energy flux per unit detector area and per unit energy. The prefactors have transparent origins: δ (from E), δ² (from dΩ′), and γ (from the dt′ mapping) combine into ${\delta }^3\gamma ,\sqrt{1-u}$ accounts for gravitational redshift, the cosine factor uses the comoving local emission angle σ′, and the Jacobian encodes Schwarzschild gravitational lensing. The dS′/D² factor represent the differential area of the emitting patch and the stellar distance information.