How the Pulsar Simulator Works

Outside a non-rotating, spherical star, spacetime is described by the Schwarzschild metric, as around the black hole on our companion page. Every light ray stays in a plane through the centre, and its path depends on one number, the impact parameter b. A photon leaving radius r at angle α from the outward radial direction, as measured by an observer at rest there, has

On the black hole page each pixel's ray is integrated step by step. Here the star blocks every ray that would come close to the photon sphere (we keep R ≥ 3.6 M), so a faster route is possible. We tabulate two functions once, at start-up, and build every ray from them.

Swept angle. The angle a photon sweeps around the centre as it travels outward from r to infinity is

Travel time. The coordinate time for the same trip diverges, so we tabulate a regularized version, from which every finite time difference follows:

Both integrals have an integrable singularity where the ray is tangential (α → π/2). The substitution w = u(1 − t²), t = (cos α / √k) sinh z, with k = 2 sin²α (1 − 3u)/(1 − 2u), makes the integrand smooth, and 60-point Gauss–Legendre quadrature then gives F to 3 × 10⁻⁶ rad. The table holds F and S = T + r cos α + 2 ln r, which stays finite as r → ∞, on a 512 × 512 grid in α and u ∈ [0, 1/3.2]. A Web Worker builds it in about a second and a half, and the graphics card reads it as a floating-point texture.

Any ray from the table. For a ray between radius r and a camera at r_c, with the camera seeing it at angle θ (sin θ = b √(1 − 2u_c)/r_c):

with u_p the root of u²(1 − 2u) = 1/b² below 1/3, and the same combinations of T for arrival times. In flat space F(α, u) = α, and both formulas reduce to straight lines.

Where a point appears. A point at radius r, at angle ψ from the direction to the camera, has two images. The primary image is made by a ray that sweeps Φ = ψ and appears on the point's side of the star. The secondary image is made by a ray that sweeps 2π − ψ around the other side. Φ increases monotonically with the emission angle (we check this numerically). So each image is found by bisection on α between 0 and the grazing angle α_lim = π − arcsin(b_R √(1 − 2u) u), where b_R = R/√(1 − 2/R). If Φ(α_lim) is less than the target, the star hides that image. Rays that would loop round the star more than once are blocked by its surface when R > 3 M, so these two images are the only ones.

The graphics card runs this bisection (26 steps) for every vertex of every field line, every frame, in a separate render pass that writes screen positions into a texture. The line pass then draws each segment as a screen-space quad. Segments are dropped where an image jumps across the axis behind the star, which is where the primary and secondary images swap sides.

The camera is a static observer at r_c looking at the star's centre. A pixel at angle θ from the centre has b = r_c sin θ / √(1 − 2/r_c). If b < b_R, its ray reaches the surface, leaving it at emission angle sin α = b √(1 − 2/R)/R and at angle

from the sub-camera point, in the direction of the pixel. θ is computed as arctan(tan(¹/₂FOV) |ndc|) rather than from a dot product, because a star 3,000 M away subtends milliradians and arccos loses precision near 1 in single precision. The surface point is looked up at its retarded time, rotation phase φ_obs − Ω Δt, with Δt from the T table relative to the ray from the sub-camera point.

Frequency shift. The surface element moves at speed β = ΩR sin θ_spin / √(1 − 2/R), as measured by a static observer there. With k the photon's direction in that observer's frame and v the surface velocity,

The last factor is the static camera's own gravitational blueshift. Because I_ν/ν³ is conserved along a ray, a blackbody at temperature T is seen as a blackbody at gT, scaled by the beaming function evaluated at the comoving emission angle α′: either 1 (isotropic) or (1 + 2.06 cos α′)/1.687, the electron-scattering limb darkening of Chandrasekhar (1960), normalized to the same total flux. Colors come from the CIE 1931 color-matching functions, as on the black hole page. The hot spots are drawn at 8,000 K and the rest of the surface at a chosen fraction of that, as a visible stand-in for the X-ray temperatures. The “brightness” color mode instead shows g⁴I′ on a false-color scale.

Hot spots. Each spot is a circle of fixed angular radius in the co-rotating frame. By default the radius is that of the open-field polar cap of a flat-space dipole, sin²θ_pc = R/R_LC with R_LC = c/Ω (Goldreich & Julian 1969). Plasma-filled solutions have more open magnetic flux, so their polar caps are somewhat wider (Contopoulos et al. 1999).

Rotating dipole in vacuum. A point magnetic dipole m(t) turning at Ω about +y produces the field (Jackson 1999, §9.3, in its magnetic form)

with every m evaluated at the retarded time t − r. This is an exact solution of Maxwell's equations in empty space. We check numerically that ∇·B = 0 and that B satisfies the wave equation. Deutsch (1955) solved the problem for a conducting sphere of radius R. Outside the star his field differs from the point dipole's by terms of relative size (R/R_LC)², under 4% even for the 716 Hz pulsar, and much less for slower ones. The pattern is stationary in the rotating frame, so field lines are traced once, from footpoints on the surface, by fourth-order Runge–Kutta on dx/ds = B̂ with steps of 2.5% of r (at most 0.06 R_LC), and then rotated rigidly each frame. Footpoints are chosen so that the static-dipole loop through each would cross the magnetic equator at radii spaced geometrically from 1.7 R to 2.4 R_LC.

Static dipole. The first term alone, turned rigidly. Lines are colored by whether their equatorial radius lies inside the light cylinder.

Plasma-filled sketch. A real magnetosphere is filled with plasma that carries currents (Goldreich & Julian 1969). The resulting “force-free” field has no closed-form solution for a tilted star; the standard results are numerical (Contopoulos et al. 1999 for the aligned case, Spitkovsky 2006 for tilted ones). The sketch is built from three pieces, only two of which are calculations:

• Inside the light cylinder, lines whose dipole loop closes within it (equatorial radius r_eq < R_LC) are drawn as static dipole loops.
• Far outside, the field is Bogovalov's (1999) exact oblique split monopole, the tilted generalization of Michel (1973). Field lines are Archimedean spirals on cones, φ = φ₀ − Ω(r − R), and lines from the two magnetic hemispheres are separated by a current sheet on the surface m̂(t − r + R)·n = 0. Each open line is assigned an asymptotic magnetic colatitude Θ by conserving its share of the open flux, 1 − cos Θ = R_LC/r_eq, which spreads the open flux evenly over each hemisphere as in the split monopole.
• Between 0.5 and 1.5 R_LC, each open line's colatitude and twist are blended smoothly (a cubic smoothstep) from the dipole to the monopole. This part is an interpolation, not a solution of any equation. The real field there contains a separatrix, a Y-shaped point where the closed zone meets the sheet, and currents that this sketch does not have.

The sheet is drawn out to 6 R_LC, about one wavelength of the striped wind (2πR_LC), as a translucent surface between 48 spiral spokes.

Gravity and the field. All three fields are flat-space solutions. In Schwarzschild spacetime a static dipole's vector potential becomes

(Wasserman & Shapiro 1983). For a fixed dipole moment, this makes the field at the surface of a 1.4-sun, 12 km star 1.36 times stronger in its radial component and 1.50 times stronger in its tangential component than the flat-space field. The field lines near the surface would therefore lean slightly more toward the equator. We have not included this. Frame dragging, which also alters the near field, is omitted too.

The observed flux at a given moment is the sum over the visible surface of the observed intensity times the solid angle each element subtends. For an observer at distance D → ∞, the light-bending map gives the solid angle directly (Pechenick et al. 1983; Beloborodov 2002):

where dS is the coordinate area of a fixed patch of the static sphere and I′ the bolometric intensity in the frame moving with the surface. ψ(α) and its derivative come from a 2,048-point evaluation of F(α, 1/R), resampled onto a uniform grid in ψ.

We integrate over a fixed grid on each spot in the co-rotating frame (14 rings by 40 sectors of equal area) at 384 rotation phases. Each contribution is placed at its arrival phase φ + Ω Δt by linear interpolation into 128 bins. This handles light-travel time correctly. Spreading contributions by arrival time automatically supplies the factor dt_emit/dt_obs that turns δ⁴ in the intensity into the δ⁵ of Poutanen & Gierliński (2003), so it is not applied twice. The rest of the surface, at uniform temperature, adds a constant (a uniform sphere looks the same at every phase). This is the “Schwarzschild + Doppler” approximation of Miller & Lamb (1998), without the oblateness correction used in NICER analyses (Bogdanov et al. 2019).

As a check, we also compute the profile by brute force. We trace 360,000 rays from a distant observer, each assigned an equal share of the image disk, and sum the observed intensity of the surface points they hit at their retarded times. The two methods agree to 0.15% of the peak (below).

How pulsars generate coherent radio emission is an open problem, but the geometry of the beam is well tested. We follow Blaskiewicz, Cordes & Wasserman (1991) and compute it numerically rather than with their first-order formulas:

1. Emission points. 28 × 120 points per pole on the sphere of radius r_em, inside the open-field boundary sin²θ_m = r_em/R_LC. Each point is weighted by a brightness pattern in q = θ_m/θ_open: a conal ring (Gaussian at q = 0.72, width 0.11), a core (Gaussian at 0, width 0.22), or uniform.
2. Beam direction. Charges slide along the static-dipole field line while co-rotating. Their velocity is v = κB̂ + Ω × r / √(1 − 2/r), with κ fixed by |v| = c. Highly relativistic charges radiate along v.
3. Polarization. Curvature radiation is polarized along the component of the charge's acceleration perpendicular to v. The acceleration is (dv/ds) κ + Ω × v: the change along the field line, by central differences, plus the rotation of the frame.
4. Light bending. With bending on, the ray's direction at infinity is the emission direction rotated within its plane through the centre by F(α, u_em) − α. The polarization vector is parallel-transported: its components perpendicular to and within that plane are preserved. The arrival time is T(α, u_em). With bending off, it is −r·k̂.
5. Observation. At each of 1,440 rotation phases, every emission direction within 5° of the line of sight contributes intensity and the Stokes parameters Q and U, weighted by a Gaussian of width 1°, at its arrival phase. The real beaming angle, about 1/γ, is far smaller. The width is a smoothing kernel that keeps the sampled curves continuous. The plotted angle is ½ arctan(U/Q), shown where the intensity exceeds 3% of the peak.

Beam display. The drawn beam uses the brightest ring of emission directions (q = 0.72, or the cap edge for a uniform beam), 40 per pole, from the same calculation. The cone draws each direction as a straight line out to R_LC. “Light in flight” draws where that light is at one instant. A photon emitted a time a ago from x along k̂ is now at Rot_y(x + ak̂, φ − Ωa). In the co-rotating frame this is the fixed curve Rot_y(x + ak̂, −Ωa), an Archimedean spiral with one turn per c P, drawn for ages up to 1.5 P. “Textbook beams” is not computed: it draws two rigid 4° cones along the magnetic axis, of the same length, for comparison only. The photon paths are straight lines from the emission point; bending of the paths themselves is left out, although their images on screen are lensed like everything else.

The dashed comparison curve is the rotating vector model (Radhakrishnan & Cooke 1969):

The sign and zero point of position angles vary between papers (Everett & Weisberg 2001). We measure the angle from the projected spin axis, with the sign matched to our numerical model in the low-altitude limit. For a beam emitted at height h, Blaskiewicz et al. predict to first order in h/R_LC that the pulse centre moves earlier by 2h/R_LC radians and the steepest point of the polarization swing later by the same amount, a total separation of 4h/R_LC.

Periods P and spin-down rates Ṗ are from the ATNF Pulsar Catalogue, version 2.6 (Manchester et al. 2005). The fastest, PSR J1748−2446ad at 716 Hz, was discovered by Hessels et al. (2006); the catalogue lists no spin-down rate for it, so no field or age is quoted. The “This star” panel quotes three standard derived quantities, which assume a moment of inertia I = 10⁴⁵ g cm²:

The field formula assumes a vacuum dipole at right angles to the spin axis, radiating L = (2/3) μ²Ω⁴ sin²χ / c³ (Ostriker & Gunn 1969). Plasma-filled magnetospheres lose energy at L ≈ μ²Ω⁴(1 + sin²χ)/c³ (Spitkovsky 2006). They spin down even when aligned, and the same Ṗ then implies a dipole moment smaller by a factor of 0.6 to 0.8, depending on χ. The catalogue value is a convention, not a measurement. Masses and radii of the three NICER pulsars are from Riley et al. (2019, 2021) and Choudhury et al. (2024). A later analysis of J0740+6620 with more data gives 12.49 km (Salmi et al. 2024).

Everyday numbers. The “In Everyday Terms” figures on the main page are for M = 1.4 suns and R = 12 km (R = 5.8 M). Mean density: M / (4πR³/3) = 3.9 × 10¹⁷ kg/m³, using the areal radius; a teaspoon is 5 mL, giving 1.9 × 10¹² kg, against about 5 × 10¹¹ kg for 8.1 billion people at a 62 kg average adult mass (an overestimate, since it counts children as adults). Surface gravity measured by an observer at rest on the surface: GM / (R² √(1 − 2GM/Rc²)) = 1.6 × 10¹² m/s²; impact speed from 1 m: √(2gh) ≈ 1,800 km/s. Clock rate at the surface: √(1 − 2GM/Rc²) = 0.81. Spin: 716 turns per second is 43,000 rpm, against roughly 20,000 rpm for a kitchen blender. Field: 10¹² G = 10⁸ T, against 3 T for a clinical MRI scanner.

The physics code runs unchanged in a test page, in headless Chrome. It is compared against an independent Python reference that uses adaptive quadrature (SciPy's quad), a direct ODE integration of the photon orbit equation u″ + u = 3u², and a separate ray-shooting solver for image positions. The two share no code.

The radio results are a first-order formula tested against a calculation that includes all orders. The computed separation runs 5% to 13% above 4h/R_LC, and splits unevenly: the polarization lag is 2.2h to 2.7h and the pulse lead 1.8h to 2.0h, against 2h each in the formula.

Two numbers we checked rather than assumed. Beloborodov's (2002) approximation, 1 − cos α = (1 − cos ψ)(1 − 2/R), is often quoted as accurate to about 1%. Against the exact sweep, its largest error in ψ is 0.49° at R = 7 M, 1.1° at 5.8 M and 2.4° at 5 M. At 4 M, the compactness of PSR J0740+6620, it fails near the limb by 27°. The simulator uses the exact table throughout. Second, the fraction of a 1.4-sun, 12 km star visible at once is 75.5% (ψ_max = 120.6°).