Crab Pulsar Magnetosphere

A pure WebGL/GLSL visualization that integrates the oblique rotator of PSR B0531+21, the pulsar wind termination shock, bow shock, jets, striped-wind current sheet, and supernova remnant into a single three-dimensional emission-extinction volume

1. Overview

The core of this project is not to overlay a 2D nebula image on top of a pulsar. Instead, within a single perspective ray cast from the camera, the outer supernova shell, pulsar wind nebula, termination shock torus, expanding bow shock, bipolar jets, magnetospheric current sheet, and central neutron star are all computed in the same world coordinates and the same depth order.

As a result, tilting the camera rotates all structures together, and foreground gas naturally attenuates the torus and jets in the background. The approach of making only the central pulsar 3D while compositing the surrounding nebula as a screen-fixed 2D background is not used here.

In this project, the Kali spherical inversion is used not as an independent cosmic background, but solely as a sparse density modulator for the supernova shell.

2. Scene Composition and Density Model

Ultimate is not a compositor that adds together several finished images; it is a single volume sampler that, at a given position $\mathbf x\in\mathbb R^3$, returns the emission coefficient $\mathbf j(\mathbf x)$ and the extinction coefficient $\sigma_t(\mathbf x)$.

The reason these model boundaries matter is that returning each structure to a separate screen layer would again decouple tilt, perspective, and occlusion.

3. Observed Values, Model Values, and Visual Scale Factors

A distinction is made between physically measured values and values chosen to make structures readable on screen.

The neutron star, the light cylinder (thousands of km in scale), the termination shock (roughly a light-year in scale), and the supernova remnant (several light-years in scale) cannot all be displayed simultaneously at true proportions on a single screen. Timing constants and spatial scale factors are not the same kind of quantity. The spatial values are world-space layout choices intended to show meaningful relationships between structures, not physical units.

4. Camera and Ray-Volume Intersection

4.1 Perspective Camera

Screen coordinates are divided by the vertical resolution to correct for the aspect ratio.

Starting from the camera position $\mathbf o$ looking toward the origin, the forward vector $\mathbf f$, right vector $\mathbf e_r$, and up vector $\mathbf e_u$ are constructed, and the ray direction is computed.

Implementation — **shaders/ult.frag**, **main()**

vec2 screen = (fragCoord - 0.5 * iResolution.xy) / iResolution.y;
vec3 forward = normalize(-rayOrigin);
vec3 right = normalize(cross(forward, vec3(0.0, 1.0, 0.0)));
vec3 up = normalize(cross(right, forward));
vec3 rayDirection = normalize(
    forward * 1.48
    + right * screen.x * uZoom
    + up * screen.y * uZoom
);

Here $z$ is the Field Zoom from the UI. The default camera distance is $3.18$ and the default pitch is $0.27\,\mathrm{rad}$; dragging changes yaw and pitch. This single camera is applied exactly once to all structures, from the outer shell down to the stellar surface.

4.2 Bounding Sphere Intersection

The ray is

and the intersection with an origin-centered sphere of radius $R$ is found from the two roots of the following quadratic equation.

Implementation — **shaders/ult.frag**, **raySphere()**

float projection = dot(rayOrigin, rayDirection);
float discriminant = projection * projection
                   - dot(rayOrigin, rayOrigin) + radius * radius;
if (discriminant < 0.0) return vec2(-1.0);
float root = sqrt(discriminant);
return vec2(-projection - root, -projection + root);

Pixels with $\Delta<0$ do not sample the volume at all. Only pixels that intersect the sphere advance through the $[s_{near},s_{far}]$ interval in 52 steps.

5. Emission-Extinction Integration

The position at each step is $\mathbf x_k=\mathbf r(s_k)$, and the volume sampler returns vec4(emission.rgb, extinction).

The current discrete approximation of the continuous radiative transfer equation is as follows.

The step alpha in GLSL is

Implementation — **shaders/ult.frag**, 52-step integration in **main()**

float alpha = 1.0 - exp(-sampleValue.a * stepSize * 2.15);
accumulated += transmittance * sampleValue.rgb * stepSize * 1.36;
transmittance *= 1.0 - alpha;
if (transmittance < 0.012) break;

1 - alpha is exactly $\exp(-2.15\sigma_t\Delta s)$. The values 1.36 and 2.15 are not real spectroscopic coefficients; they are rendering scale factors that separate the optically thin synchrotron emission from the weak extinction.

Optically thin synchrotron emission accumulates more strongly than extinction. This separation resolves the problem in earlier implementations where the outer shell acted like a thick fog that swallowed the central torus and jets entirely.

To reduce the banding artifacts of fixed-step sampling, a hash jitter is added to the first sample position of each pixel.

float stepSize = (farDistance - nearDistance) / float(VOLUME_STEPS);
float jitter = hash(fragCoord + fract(iTime)) * stepSize;
float sampleDistance = nearDistance + jitter;

Integration terminates early when transmittance drops below $T<0.012$. Background stars are added only after the volume integration, so the foreground nebula density correctly attenuates the background.

6. 3D Density Field of the Supernova Shell

6.1 Anisotropic Ellipsoid

The outer remnant is defined not as a sphere but as an ellipsoidal distance with axis ratios $(1.22,0.72,0.96)$.

$\mathbf p_a$ is a slowly advected position. Different phases and frequencies are assigned to the three axes so that the outer material does not rotate rigidly like a flat disk.

Implementation — **shaders/ult.frag**, **sampleUltimateVolume()**

float delayedFlow = time * 0.055 - radius * 3.7;
vec3 advectedPoint = point + 0.045 * vec3(
    sin(2.8 * point.y + delayedFlow),
    sin(2.6 * point.z + delayedFlow + 2.1),
    sin(3.0 * point.x + delayedFlow + 4.2)
);
float ellipsoidRadius = length(advectedPoint / vec3(1.22, 0.72, 0.96));
float shellBand = exp(-pow((ellipsoidRadius - 0.91) / 0.12, 2.0));

$\omega_{vis}=0.055$ and $k=3.7$ are advection coefficients in visual world-space time and are not physical SI units.

6.2 FBM and Ridged Noise

Large gas clumps are handled by a 3-octave FBM, while sharp filaments are handled by ridged noise.

Implementation — **shaders/ult.frag**, **fbm3()**·**ridged3()**

for (int octave = 0; octave < 3; octave++) {
    value += amplitude * noise3(point);
    point = point.yzx * 2.03 + vec3(5.2, 1.3, 8.1);
    amplitude *= 0.48;
}

float ridge = 1.0 - abs(2.0 * noise3(point) - 1.0);
value += amplitude * ridge * ridge;

A fixed 3D basis rotateNoiseSpace() is applied first so that the value-noise lattice of the first octave does not appear as a square grid or voxel band. Domain warping, which feeds the low-frequency cloud value back into the coordinates, transforms the filaments from a cell grid into curved gas flows.

6.3 Limited Use of Kali Spherical Inversion

The Kali-type spherical inversion is iterated three times.

Implementation — **shaders/ult.frag**, **kaliDensity()**

for (int fold = 0; fold < 3; fold++) {
    folded = abs(folded)
           / max(dot(folded, folded), 0.075)
           - 0.659;
}

However, this result is not rendered as an independent fractal universe. The orbit-length energy is used only as a sparse shell density modulator to create irregular voids and condensations.

7. Pulsar Wind Nebula: Torus, Bow Shock, and Jets

All points are projected into the pulsar frame using the rotation axis $\mathbf s$ and its two perpendicular basis vectors $\mathbf e_1,\mathbf e_2$.

Implementation — **shaders/ult.frag**, **sampleUltimateVolume()**

vec3 frame = vec3(
    dot(point, equatorA),
    dot(point, equatorB),
    dot(point, spinAxis)
);
float equatorialRadius = length(frame.xy);
float azimuth = atan(frame.y, frame.x);

The equatorial radius is defined as $\rho=\sqrt{p_{f,x}^2+p_{f,y}^2}$ and the axial distance as $z_s=p_{f,z}$.

7.1 Termination Shock Torus

float torusDistance = length(vec2(
    equatorialRadius - WORLD_TERMINATION_SHOCK,
    frame.z
));
float torus = exp(-pow(torusDistance / 0.025, 2.0));

Azimuthal knot noise and approach-side weighting are added so the result does not become a perfectly uniform neon ring.

7.2 Expanding Bow Shock

The three bow shocks are staggered in age by $1/3$ of a cycle from one another.

float age = mod(time / 4.6 + float(wispIndex) / 3.0, 1.0);
float wispRadius = WORLD_TERMINATION_SHOCK + age * 0.33;
float wispDistance = length(vec2(equatorialRadius - wispRadius, frame.z));
float arc = exp(-pow(wispDistance / (0.012 + age * 0.020), 2.0));
float wisp = arc * fragment * pow(1.0 - age, 1.7);

As each bow shock ages, its width broadens and its brightness decreases as $(1-a_i)^{1.7}$. An azimuthal sinusoid clips part of each ring so that the bow shocks appear as incomplete arcs, as seen in observations, rather than as closed circles.

7.3 Fixed-Axis Jet

The jet width increases gradually along the axis.

float alongJet = frame.z;
float jetBend = 0.018 * sin(12.0 * alongJet - time * 0.55)
              * smoothstep(0.10, 0.72, abs(alongJet));
vec2 bentJet = frame.xy
             - vec2(jetBend, -jetBend * 0.45) * sign(alongJet);
float jetRadius = length(bentJet);
float jetWidth = 0.020 + 0.040 * abs(alongJet);

A slow kink is applied to the centerline, but the owner of the axis is always $\mathbf s$. An implementation in which the jet precesses like a propeller around the rotating magnetic axis $\mathbf m(t)$ is not used.

8. Oblique Rotator and Magnetosphere

8.1 Magnetic Axis

Given the rotation axis $\mathbf s$, its perpendicular basis $\mathbf e_1,\mathbf e_2$, and the magnetic axis inclination $\alpha$, the magnetic moment is

Implementation — **shaders/pulsar-model.glsl**, **magneticAxis()**

return normalize(
    cos(MAGNETIC_OBLIQUITY) * spinAxis
    + sin(MAGNETIC_OBLIQUITY)
    * (cos(spinPhase) * basisA + sin(spinPhase) * basisB)
);

The two magnetic poles are always exactly opposite, at $+\mathbf m$ and $-\mathbf m$. The error of using the observed $0.40016$ cycle separation between two brightness peaks as the spatial angle between the two poles is avoided.

8.2 Light Cylinder and Spin-Down Luminosity

Implementation boundary — **shaders/pulsar-model.glsl**

const float LIGHT_CYLINDER = 0.105;
const float MAGNETIC_OBLIQUITY = 0.87266463;

These two physical expressions serve as the structural reference for the model. Because the current shader does not integrate $c$, $\mu$, and $\Omega$ in SI units to compute luminosity, no placeholder code corresponding to $L$ is presented. What is implemented is only the normalized boundary based on $R_{LC}$ and the geometry of $\alpha=50^\circ$.

The WORLD_LIGHT_CYLINDER=0.105 on screen is not a direct conversion of the physical length of $R_{LC}$; it is a visual radius used to read the near-zone and wind-zone boundaries.

8.3 Closed Dipole Field

The interior of the light cylinder uses the dipole flux invariant.

Implementation — **shaders/ult.frag**, **sampleUltimateVolume()**

vec3 moment = magneticAxis(
    spinPhase - radius / WORLD_LIGHT_CYLINDER * 0.12,
    FIELD_PA
);
float magneticCosine = dot(radialDirection, moment);
float sinThetaSquared = max(0.016, 1.0 - magneticCosine * magneticCosine);
float radiusUnits = radius / WORLD_LIGHT_CYLINDER;
float fluxInvariant = sinThetaSquared / max(radiusUnits, 0.08);
float closed = 1.0 - smoothstep(0.88, 1.12, radiusUnits);
float fieldLines = pow(
    max(0.0, 0.5 + 0.5 * cos(10.0 * fluxInvariant)),
    16.0
) * closed;

The shader converts the contours of $\Psi$ into thin emission bands and attenuates the closed-field region as it approaches the light cylinder.

8.4 Retarded Striped-Wind Current Sheet

Outside the light cylinder, the shader uses the null surface of the retarded split-monopole current sheet rather than a vacuum dipole lighthouse beam.

In the code,

is computed, and

Implementation — **shaders/ult.frag**, **sampleUltimateVolume()**

float lightCylinderUnits = radius / WORLD_LIGHT_CYLINDER;
vec3 retardedAxis = magneticAxis(
    spinPhase - lightCylinderUnits,
    FIELD_PA
);
float sheetCoordinate = dot(retardedAxis, radialDirection);
float sheetWidth = 0.028 + 0.005 * min(lightCylinderUnits, 5.0);
float sheet = exp(-pow(sheetCoordinate / sheetWidth, 2.0))
            * smoothstep(0.86, 1.14, lightCylinderUnits)
            * exp(-radius * 2.0);

a finite-width 3D emitting surface is constructed. This structure is not a 2D spiral line drawn on screen; it is an actual volume surface whose front and back faces intersect along the ray.

9. Pulse Light Curve — Not a Pulsating Sphere

The wrapped distance around the phase boundary is

and the main pulse and interpulse are approximated by two Gaussians.

Implementation — **shaders/pulsar-model.glsl**, **wrapPhase()**·**pulseProfile()**

float wrapPhase(float d) {
    d = abs(fract(d));
    return min(d, 1.0 - d);
}

float mainPulse = exp(-0.5 * pow(mainDistance / MAIN_SIGMA, 2.0));
float interPulse = exp(-0.5 * pow(interDistance / INTER_SIGMA, 2.0));
return mainPulse + INTERPULSE_AMP * interPulse;

This function does not change the star's size or position. The caustic brightness variation that occurs when the rotating magnetospheric emission pattern sweeps past the line of sight is applied only to the small polar cap and the central magnetospheric emission.

Rendering the actual $29.58936\,\mathrm{Hz}$ rotation in real time would yield only about two frames per rotation on a 60 Hz display. Therefore,

const float CRAB_FREQUENCY = 29.58936;
const float PULSAR_SLOWDOWN = 60.0;
float phase = fract(time * CRAB_FREQUENCY / PULSAR_SLOWDOWN);

a 60× slowdown is applied to make the structure readable. PULSAR_SLOWDOWN=60 is an observational time scale factor, not a physical quantity.

10. Central Neutron Star and Noisy Corona

The neutron star is not blurred into the volume density; instead it is intersected precisely as a separate small sphere. For a star of radius $R_*=0.012$, the front-face depth of the sphere is found, and the moment a raymarch sample reaches that depth, the surface emission is inserted and the transmittance behind it is reduced. This allows gas in front of the star to occlude it, while the star itself occludes the magnetosphere behind it.

Keeping a large white sphere makes an unresolved pulsar look like a planet. The current radius has been reduced from the previous $0.035$ to $0.012$, retaining only the minimum visual exaggeration needed to keep it visible on screen.

The noise is not used to roughen the neutron star surface like rock. Instead, three terms are combined in the corona within radius $0.14$.

1. A compact envelope of the form $\exp[-(r/0.058)^2]$

2. Magnetic pole concentration of the form $|\hat{\mathbf r}\cdot\mathbf m|^9$

3. Localized brightness holes formed by rotated ridged noise

Implementation — **shaders/ult.frag**, **sampleUltimateVolume()** and **main()**

float polarFlux = pow(abs(dot(radialDirection, coreMoment)), 9.0);
float returnCurrent = exp(-pow(dot(radialDirection, spinAxis) / 0.16, 2.0));
float coreEnvelope = exp(-pow(radius / 0.058, 2.0));
float coreTexture = 0.20 + 0.80 * ridged3(
    rotateNoiseSpace(point * 56.0)
    + vec3(time * 0.09, -time * 0.05, time * 0.07)
);
float corePlasma = coreEnvelope * coreTexture
                 * (0.16 + 0.62 * polarFlux + 0.22 * returnCurrent);

const float STAR_RADIUS = 0.012;
vec2 starHit = raySphere(rayOrigin, rayDirection, STAR_RADIUS);

As a result, the center reads not as a perfect white ball but as a small opaque core surrounded by fragmented, anisotropic plasma.

11. Color, Exposure, and Tone Mapping

All emission accumulates in linear light, and ACES approximate tone mapping is applied exactly once at the end.

Implementation — **shaders/ult.frag**, **acesToneMap()** and **main()**

vec3 numerator = color * (2.51 * color + 0.03);
vec3 denominator = color * (2.43 * color + 0.59) + 0.14;
return clamp(numerator / denominator, 0.0, 1.0);

// main()
vec3 color = acesToneMap(accumulated * uIntensity * 1.34);
color = pow(max(color, 0.0), vec3(0.91));
color += (hash(fragCoord + fract(iTime)) - 0.5) / 255.0;

A mild gamma shaping and dithering of less than $1/255$ are then added. Because each structure's image is not individually tone-mapped before compositing, color and energy ordering are preserved even where the torus and shell overlap.

The palette references the physical structures listed below, but it is not the direct result of integrating spectral emission lines.

• Cyan-white: synchrotron interior, approach-side torus, jets

• Blue to magenta: current sheets and magnetosphere of opposite polarities

• Red to orange: condensed filaments in the supernova remnant

12. Performance Design

• Only pixels that intersect the bounding sphere undergo the 52-step raymarch.

• Shell noise is evaluated only in regions where shellBand > 0.002.

• The noisy corona is evaluated only where $r<0.14$.

• The torus and jets also execute their expensive functions only within their own narrow distance bounds.

• Integration terminates early when transmittance satisfies $T<0.012$.

• All GLSL files are fetched in parallel and the program is compiled strictly at startup.

• Shader load or compilation failures are not silently replaced by a fallback screen.

js/shaders.js declares GLSL files as an explicit URL manifest so that both the project bundler and the browser can identify dependent files.

const SHADER_URLS = {
  vertex: new URL('../shaders/vertex.vert', import.meta.url),
  ultimate: new URL('../shaders/ult.frag', import.meta.url),
  pulsarModel: new URL('../shaders/pulsar-model.glsl', import.meta.url)
};

Combining a folder URL with dynamic filenames is not used.

13. Validation Criteria

The fact that the screen looks impressive alone does not establish that the model is correct. The following conditions are verified.

1. Common camera — when dragging or tilting, the shell, torus, bow shock, jets, and current sheet all move under the same perspective.

2. Depth occlusion — the foreground shell and diffuse gas attenuate background structures, and the star occludes the background magnetosphere.

3. Axis separation — the jets follow the fixed rotation axis, while only the magnetic axis and current sheet rotate.

4. Opposite magnetic poles — the polar caps are always located at $+\mathbf m$ and $-\mathbf m$.

5. Pulse semantics — the pulse toggle changes only emission intensity, not the star's radius or position.

6. Bow shock direction — the bow shocks propagate outward from the termination shock and are not locked as closed rings.

7. Error transparency — missing GLSL files or compilation failures are exposed via an error overlay.

8. Inactive model blocked — the legacy pulsar.frag is inaccessible from the official loader and UI.

14. Directory and Runtime Flow

spinning-pulsar-of-the-crab-nebula/
├── index.html                   # Ultimate 3D single entry point and control UI
├── style.css                    # Fullscreen canvas and off-canvas studio
├── README.md                    # Physics, mathematics, and rendering documentation
├── js/
│   ├── shaders.js               # Explicit GLSL URLs, model injection, and preset registry
│   ├── renderer.js              # WebGL compilation, uniforms, and 52-step render loop driver
│   ├── controls.js              # Drag/zoom/toggle and accessibility state
│   └── main.js                  # Initialization and explicit error boundary
└── shaders/
    ├── vertex.vert              # fullscreen triangle pair
    ├── pulsar-model.glsl        # Single owner of observational constants, pulse profile, and rotator frame
    ├── ult.frag                 # Canonical 3D volume renderer
    └── pulsar.frag              # Inactive legacy projection; not currently read by the loader

The runtime flow is as follows.

index.html
  → main.js
    → shaders.js ──fetch──> vertex.vert / pulsar-model.glsl / ult.frag
    → renderer.js ─compile/link──> Ultimate WebGL program
    → controls.js ─uniform state──> renderer.js

pulsar-model.glsl is injected exactly once at the /*__PULSAR_MODEL__*/ marker. If the marker is absent or appears more than once, initialization is halted.

15. UI and Controls

The GLSL Source tab in the panel shows the currently compiled ult.frag and the injected physics model, not a JavaScript copy.

16. Honestly Stated Limitations

• This is not a full MHD/PIC simulation. It combines mathematical conditions drawn from force-free and observational structures into a real-time visual model.

• Spatial scales are not proportional. Multiple visual scale factors are used so that the neutron star and the termination shock can be read simultaneously.

• Colors are a chosen palette. Actual spectral response, polarization, and Doppler spectrum are not computed.

• $\alpha=50^\circ$ is a representative value. It should not be interpreted as the sole confirmed direct measurement.

• The light curve is an empirical approximation. The two Gaussians summarize the outcome of caustic emission and do not directly integrate particle trajectories.

• Kali modulation is not a physical law. It is a procedural density tool for creating multi-scale voids in the remnant.

• The outer advection runs on visual time. It is not a physical time integration that compresses the actual centuries-long evolution of the supernova remnant into real time.

17. References and Lineage

Pulsar Magnetosphere and Wind

• Spitkovsky (2006), Time-dependent Force-free Pulsar Magnetospheres

• Bogovalov (1999), On the physics of cold MHD winds from oblique rotators

• Pétri (2022), radiative oblique magnetospheres

Crab Nebula Observations and Timing

• Chandra/Hubble Crab pulsar-wind movie

• Nakamori et al. (2025), JBO timing reference

• Jodrell Bank/HEASARC Crab ephemeris

Procedural Density Fields

• Kali — ShaderToy, spherical-inversion fractal family

• Inigo Quilez, procedural noise, SDF, and raymarching reference