Department of Physics · UDESC

Seminar · 2026

Ray tracing in relativistic stars

From geometry
to observation

How spacetime curvature decides the image, the spectrum and the light curve that an X-ray telescope receives from a neutron star.

Outline

Four parts, one path

I

The star inside

What holds a neutron star up: the TOV equation, the equation of state and the mass–radius curve.

II

Relativistic ray tracing

Null geodesics in Schwarzschild: light bending, redshift and time delays — and how to follow the rays backwards, from the telescope to the surface.

III

Pulsars and neutron stars

The laboratories: pulsars and magnetars, hot spots, rotation, atmospheres — and what NICER has already measured.

IV

Recent results

My work: magnetar light curves, a massive low-field magnetar, accretion columns — and tests of alternative theories of gravity.

The premise

Everything we know arrives as photons

Compact stars are natural laboratories of strong-field gravity — but we know them only through the photons that left them and crossed a curved spacetime.

metric→null geodesics→image · spectrum · light curve→telescope→M, R, field, gravity

The direct problem

Given the star and the metric, predict what is observed. That is ray tracing.

The inverse problem

Given what is observed, infer mass, radius, spot geometry and magnetic field.

The test

If the metric is not Schwarzschild, the light curves change: observation tests gravity.

I

Part one

The star
inside

Before following the light, the source: what holds more than a solar mass within 12 km — the pressure of the densest matter we know, against the gravity of general relativity.

anatomyTOV equilibriumequation of statemass and radius

Anatomy · a TOV solution with the SLy equation of state

From the crust to the unknown core

to scale · the same r axis above and below

Atmosphere and ocean: centimetres to metres of H, He or C — where the light of Part II comes from.
Outer crust: a lattice of ever more neutron-rich nuclei in a sea of electrons, up to 4×10¹¹ g/cm³, when neutrons start to drip out of the nuclei.
Inner crust: nuclei and free, superfluid neutrons; at the bottom, “nuclear pasta”. Up to ~ρ₀/2.
Outer core: neutrons, with some protons, electrons and muons, up to ~2ρ₀.
Inner core: up to ~6ρ₀. Hyperons? Free quarks? The equation of state decides.

ρ₀ = 2.7×10¹⁴ g/cm³ is the density inside an atomic nucleus. TOV profile with SLy (Douchin & Haensel 2001, in the piecewise polytropes of Read et al. 2009).

The anatomy of a neutron star, computed here: the TOV equation solved with a realistic equation of state, SLy. Above, a slice of the star to scale; below, the density as a function of radius, on the same axis. The density falls from almost 10¹⁵ grams per cubic centimetre at the centre — three times that of an atomic nucleus — to the surface. The crust is a shell less than a kilometre thick: outside, a lattice of nuclei in a sea of electrons; below 4×10¹¹, neutrons start to drip out of the nuclei, and at the bottom of the crust the nuclei turn into nuclear pasta. Then the core: neutrons with a few protons, electrons and muons; and, at the centre, what nobody knows — hyperons, free quarks. Change the mass: a 2 solar-mass star is smaller, has its centre at six times nuclear density and a 400-metre crust.

Tolman · Oppenheimer & Volkoff · 1939

Hydrostatic equilibrium in general relativity

dPdr = −Gmρr²Newton · dm/dr = 4πr²ρ

dPdr = −Gmρr² (1 + Pρc²) (1 + 4πr³Pmc²) (1 − 2Gmrc²)⁻¹ the TOV equation · ρc² = ε, the energy density

pressure has inertia

Pressure is energy, and energy resists acceleration like mass. At the centre of this star: ×1.15

pressure gravitates

Pressure is also a source of gravity and adds to the mass. At the centre: ×1.46

space curves

Radial distances stretch, and the factor grows outwards. At the surface: ×1.55

More pressure, more gravity: squeezing harder is not enough, and a maximum mass appears — which Newton does not have.

central density, ρ_c · SLy

The equation that holds the star up. In Newton, hydrostatic equilibrium says that pressure falls inwards at exactly the rate needed to hold the weight of the layers above. Tolman, Oppenheimer and Volkoff, in 1939, wrote the same thing in general relativity, and three factors appear — they light up one at a time. First, pressure has inertia: pressure is energy. Second, pressure also gravitates: it adds to the mass as a source of gravity. Third, geometry: radial distances stretch. In a 1.4 solar-mass star the three are 15, 46 and 54 per cent — not small corrections. And the consequence is in the plot: with the same equation of state, in Newton the mass grows without end with the central density; in the TOV, more pressure means more gravity, and the curve has a maximum, 2.05 solar masses for SLy. Beyond it the star is unstable and collapses into a black hole.

The equation of state

How hard matter resists being squeezed

Read et al. 2009

The equation of state gives the pressure as a function of density, P(ρ): what the TOV needs to close.
Up to ~ρ₀ laboratory nuclear physics constrains the curve; above, each model extrapolates: two- and three-nucleon forces (WFF1, APR4), Skyrme (SLy), relativistic Brueckner–Hartree–Fock (ENG, MPA1), relativistic mean field (MS1; H4, with hyperons).
Stiff — more pressure at the same density — gives a bigger, heavier star; soft, a smaller, lighter one.
Limits: sound cannot outrun light, and 2 M☉ pulsars rule out the softest.

c_s²/c² = ΓP/(ε + P). The dot marks the centre of the maximum-mass star; beyond it, the curve is faded. The jumps come from the polytrope fit, which changes Γ at 10^14.7 and 10¹⁵ g/cm³; near the centre of the heaviest stars some fits exceed c.

Tolman–Oppenheimer–Volkoff · equation of state

From nuclear physics to the M–R curve

piecewise polytropes · observational constraintshover over a curve

dPdr = −G(ε + P)(m + 4πr³P/c²)c² r(r − 2Gm/c²)relativistic hydrostatic equilibrium (TOV)

Each equation of state P(ε) gives a mass–radius curve and a maximum mass. Measuring the mass and radius of the same star picks between them.

The lines of constant compactness look ahead to Part II: below R = 1.76 R_S the whole surface is in view; at 1.5 R_S lies the photon sphere.

Piecewise polytropes of Read, Lackey, Owen & Friedman (2009); constraints (68%): NICER — Riley et al. 2019, 2021; Miller et al. 2019; Choudhury et al. 2024; PSR J0348+0432 (Antoniadis et al. 2013); GW170817 (Abbott et al. 2018).

II

Part two

Relativistic
ray tracing

Near a mass, light does not travel in a straight line: it follows the null geodesics of the metric. Ray tracing means following those geodesics backwards — from the telescope to the source.

null geodesicslight bendingredshiftthe image of a compact starthe recipe

SMACS 0723 (JWST) · NASA, ESA, CSA, STScI · public domain

25 November 1915

Einstein and his field equations

Albert Einstein, 1921 — **Albert Einstein**1921

G_μν = 8πGc⁴ T_μνcurvature = energy and momentum content

d²x^μdλ² + Γ^μ_αβ dx^αdλ dx^βdλ = 0geodesic

ds² = g_μν dx^μdx^νthe metric is the field

“Spacetime tells matter how to move; matter tells spacetime how to curve.”
John A. Wheeler

The classical tests

Sobral, 1919: light bends

effect	GR predicts
perihelion of Mercury (1915)	43″ / century
light grazing the Sun (1919)	1.75″
redshift, Pound–Rebka (1959)	2.5×10⁻¹⁵ over 22.5 m
Shapiro delay (1964)	~200 μs
GPS clocks	+38 μs / day

Sobral expedition (Ceará, Brazil): 1.98 ± 0.12″. Príncipe: 1.61 ± 0.30″. The sky over north-eastern Brazil confirmed Einstein.

Sobral, 29 May 1919 · the measurement

Two plates, the same sky

simulated plate · true scale

width of the star image

prediction

Night: the Hyades photographed in July, before dawn, with the same telescope — the Sun far away.

Eclipse: the same field on 29 May, with the Moon covering the Sun. Every star moves away from the Sun, radially.

δ = 4GMc²b = 1.75″ R_☉btwice the Newtonian value

On the whole plate nothing moves: 1″ is a tenth of a pixel. In the loupes, at true scale, as under the micrometer, the star jumps by a fraction of its own width — and the centre of a blob can be measured far better than its width.

How it is measured. Two photographs of the same patch of sky: the Hyades cluster, in Taurus. One at night, in July, when the Sun is already far away; the other during the eclipse, with the Sun right in front. The panel alternates between them, like a blink comparator, and everything is at true scale, with no exaggeration. On the whole plate nothing moves: one arcsecond is a tenth of a pixel on this screen. But nobody measured by eye: the plate went to a machine with a micrometer, which is what the loupes show — six arcseconds of sky around a star, with the wires at the night position. The star image is about two arcseconds wide, from turbulence and the optics, and during the eclipse it jumps by almost half of that, always away from the Sun, the symbol on the rim of the loupe. Near the Sun it jumps more, far away less: the deflection falls as 1/b. On the plate of the 4-inch lens, with a 5.8-metre focal length, one arcsecond is 28 micrometres. The centre of a blob can be measured much better than its width, and the fit for scale and rotation uses all the stars together. Click a star to put it in a loupe. Switch to Newton: the jump drops by half. That was the difference Sobral had to tell apart, and it measured 1.98″.

Einstein’s calculation · November 1915

Where the 1.75″ comes from

impact parameter, b/R_☉

δ = 4GM_☉c²b = 2r_sbb: from the centre of the Sun to the line of the incoming ray

[GMc²b] = (m³ kg⁻¹ s⁻²)(kg)(m² s⁻²)(m) = 1dimensionless: the angle comes out in radians

b = R_☉: δ = 5.91 km6.96×10⁵ km = 8.49×10⁻⁶ rad4GM_☉/c² = 5.91 km: twice the Sun’s r_s

1 rad = 180π × 3600″ ≈ 206,265″ → δ ≈ 1.75″Newton, with light as a particle: half of it, 0.87″

The calculation behind the 1.75. For a ray passing far from a mass M, general relativity gives the deflection δ = 4GM/c²b — twice what a Newtonian calculation gives treating light as a particle. b is the impact parameter: the distance between the centre of the Sun and the line the ray would follow without gravity, the perpendicular in the drawing. The deflection falls as 1/b, the plot below: at the limb, 1.75 arcseconds; at twice the solar radius, where κ² Tauri was, half of that. The result comes out in radians, because GM/c²b is dimensionless. 4GM☉/c² is 5.9 kilometres, twice the Sun’s Schwarzschild radius; divided by the 696 thousand kilometres of the solar radius, it gives 8.5 millionths of a radian; times 206,265 arcseconds per radian, 1.75. Sobral measured 1.98 ± 0.12 and Príncipe 1.61 ± 0.30; Newton would give 0.87. Drag b on the plot.

1859 – 1915 · the perihelion of Mercury

Mercury’s orbit does not close

perihelion advance	per century
observed (relative to the stars)	574″
other planets, by Newton	531″
unexplained remainder	43″

Δφ = 6πGM_☉c² a(1 − e²) = 0.104″ per orbitGR · Einstein, November 1915 · 415 orbits per century → 43.0″

Le Verrier found the remainder in 1859 and proposed a hidden planet, Vulcan. It was never found. Einstein got the 43″ with no adjustable parameter.

“For a few days, I was beside myself with joyous excitement.”
Einstein to Ehrenfest, January 1916

December 1915 · Russian front

Schwarzschild and the first exact solution

ds² = −(1 − r_sr) c²dt² + dr²1 − r_s/r + r²dΩ²r_s = 2GM/c² · valid outside any non-rotating spherical body

1 r_s

event horizon

1.5 r_s

photon sphere

3 r_s

ISCO: last stable orbit

2.95 km

r_s of the Sun (Earth: 8.9 mm)

He wrote to Einstein on 22 December 1915 while serving as an artillery officer. He died in May 1916.

What we need from the metric

Curvature does three things to a photon

1 · It bends the path

The photon follows a null geodesic. Near the surface of a neutron star the bending is tens of degrees: we see more than half of the star.

2 · It changes the energy

A photon emitted at r with local energy E reaches infinity with

E_∞ = E √(1 − R_S/r)gravitational redshift

3 · It delays

Different paths take different times (Shapiro delay): photons from each point of the star arrive tens of microseconds apart.

Ray tracing puts all three into a single calculation, pixel by pixel, for each instant of the rotation.

In a neutron star, R ≈ 2.5–3 R_S: none of the three effects is a small correction.

And all three depend only on the metric outside the star — which is why a light curve carries information about gravity.

Null geodesics in Schwarzschild

Two constants and a potential

The metric depends on neither t nor φ: energy and angular momentum are conserved.

E = (1 − R_Sr) ṫ

L = r² φ̇

For light only the ratio matters, the impact parameter b = L/E:

1b² = ṙ²L² + V(r), V = 1r²(1 − R_Sr)a “potential” from geometry alone

(dudφ)² = 1b² − u²(1 − R_Su)the orbit · u = 1/r

The maximum of V lies at r = 1.5 R_S: the photon sphere. Below b_c = (3√3/2) R_S the ray comes from the horizon — or from the star.

the potential and the orbit · same bdrag b

impact parameter, b/R_S

Schwarzschild geodesics

How light bends

u″ + u = 32 R_S u²u = 1/r · orbit equation for light

b > 3√3 GM/c² ≈ 2.6 R_S: the light is deflected and escapes.
b = 2.6 R_S: unstable orbit on the photon sphere, r = 1.5 R_S.
b < 2.6 R_S: captured.

The term on the right is the General Relativity correction. Without it, light would travel in a straight line.

Event Horizon Telescope · 2019 and 2022

The shadows

The ring is ~5 R_S across: it is the photon sphere, at 1.5 R_S, seen through the lens of the black hole itself — shadow diameter = 2√27 GM/c² ≈ 5.2 R_S.

Compactness

Where GR stops being a correction

the radius and the Schwarzschild radius of the same mass

Surface compared with infinity: z = (1 − r_s/R)^−1/2 − 1, g = (GM/R²)(1 − r_s/R)^−1/2 · deflection from the exact orbit integral (weak field: 2r_s/R) · ISCO, photon sphere and horizon: a 1.4 M☉ star of that size · the objects advance on their own until the first click.

The yardstick of the seminar is compactness: the radius of the object compared with the Schwarzschild radius of its mass, the size at which it would become a black hole. On the left, a zoom in powers of ten: the whole Earth would have to fit in a ball 9 millimetres in radius, a marble; the Sun, in 3 kilometres. That is why relativity is a correction here: on the Earth the surface clock loses 60 microseconds a day, and at the Sun grazing light is deflected by 1.75 arcseconds — the Sobral value. In a white dwarf, r_s/R is 5×10⁻⁴. In a neutron star the radius is only three times r_s: the escape velocity is more than half the speed of light, the redshift 23.5%, a grazing ray is deflected by tens of degrees and we see well over half of the surface. There, no series expansion will save us: the full metric is needed. And the scale carries on to the photon sphere and the horizon.

The angle that matters

From the surface to infinity: ψ(α)

A photon leaves the surface at an angle α to the normal. By the time it reaches the observer, it has turned through ψ around the centre:

b = R sin α√(1 − R_S/R)

ψ(b) = ∫₀^1/R du√(b⁻² − u²(1 − R_Su))exact

1 − cos α ≈ (1 − cos ψ)(1 − R_SR)Beloborodov (2002)

Without gravity, ψ = α and the limb is at 90°. With gravity, the limb moves further round: we see behind the star. The approximation is excellent down to R ≈ 3 R_S and fails precisely where the secondary image appears.

ψ(α): exact, Beloborodov and flat spacedrag the compactness

stellar radius, R/R_S

The naive idea: from the star to the camera

Many rays, hardly any pixels

Every point of the surface emits in every direction; each ray is followed along its geodesic, outwards.
Only those that pass through the camera aperture become pixels: here, 2 in every 10,000.
Completing this image would take ~10⁸ rays: hours of computing.
With a real telescope, 1 m² at 1 kpc, 1 photon in every 10⁴⁰ arrives.

3D count: uniform point on the surface, isotropic outward direction, exact geodesic to the camera, at 13 R_S. Only a sample is drawn, in the plane; the rays that enter the camera are all shown.

The method: from the camera to the star

One pixel, one ray

The rays leave the camera, one per pixel: the path of light travelled backwards. The direction of the pixel gives the impact parameter b.
The geodesic is integrated until the ray escapes to the sky or hits the surface at some point, at an angle α.
The pixel takes the colour of the point the light left from — corrected for redshift and for the emission angle.
Every ray computed becomes a pixel: 3,136 rays, 3,136 pixels, against the 10⁸ of the previous slide.

With spherical symmetry, ψ depends only on b: one table ψ(b), computed once, serves every pixel. Each ray is rotated into the plane of the drawing.

How it is done in practice. You start from the camera, not from the star: from each pixel a ray goes out, travelling backwards along the path light would take; the direction of the pixel gives the impact parameter. That ray is integrated: either it escapes, and the pixel is sky, or it hits the surface at some point, with some emission angle, and the pixel takes the colour of that point. On the left, the rays leaving the camera, first one by one, then by the hundreds; on the right, the image building up. Every ray computed becomes a pixel: 3,136 rays for 3,136 pixels, against the 10⁸ of the previous slide. In a spherical metric the angle ψ depends only on b, so a table ψ(b) handles every pixel — it is what the labs on my website use, and what the panel on the next slide does in real time.

What a compact star looks like

We see behind the star

the same star, more and more compact · exact integraldrag the compactness and the inclination

radius, R/R_S

axis inclination

On the left, the star as it would look without gravity; on the right, ray-traced. The grid is the same; gravity bends the light around the star and reveals the far hemisphere.

The image is larger: apparent radius R/√(1 − R_S/R).
For R < 1.76 R_S the whole surface is visible — the far pole becomes a ring.
Below 1.5 R_S the star would lie inside the photon sphere.

1959 · Pound & Rebka · Jefferson tower, Harvard

Light loses energy as it climbs

14.4 keV γ rays from ⁵⁷Fe (Mössbauer effect). Pound & Rebka, PRL 4, 337 (1960): ~10%; Pound & Snider, Phys. Rev. 140, B788 (1965): ~1%. Every 12 s source and absorber swap places, as in the experiment.

Pound and Rebka, 1959, in the 22.5-metre tower of the Jefferson laboratory at Harvard. Gravity shifts the frequency of light by gh/c², 2.5 parts in 10¹⁵. To measure it they used the Mössbauer effect of iron-57: a 14.4 keV line narrow enough to see that shift. Light that falls arrives with slightly more energy and is no longer absorbed; to compensate, they moved the absorber away at 0.74 micrometre per second — about 2.7 millimetres per hour. In the middle diagram, on top is the line of the arriving photon, shifted by gravity to +0.74; below, the line the absorber accepts, which slides by Doppler as it moves. When the two match, absorption is maximal — at v = gh/c = 0.74 μm/s. Every 12 seconds source and absorber swap places and the match moves to −0.74, just as they did to cancel systematics. On the right, the same law at a 1.4 solar-mass neutron star: with no exaggeration, 500 nanometres leave green and arrive orange, at 618 nm. It is the same factor neutrinos suffer on leaving the proto-neutron star.

Liouville · what reaches the telescope

The spectrum that climbs out of the well

Along a ray, I_ν/ν³ is conserved. With g = √(1 − R_S/R):

I_obs(E) = g³ I_em(E/g)specific intensity

T_∞ = g T

R_∞ = R/g

A black body stays a black body, only cooler. The luminosity seen from afar is L_∞ = g² L.

A spectral line of known energy E₀ arrives at g E₀: measuring the line means measuring the compactness.

thermal spectrum with an absorption linedrag the compactness

radius, R/R_S

local temperature, kT

Part II in summary

The ray tracing recipe

Choose the exterior metric: Schwarzschild; Hartle–Thorne or Kerr with rotation; a different one for a different theory of gravity.
Build the table ψ(b) and the delay Δt(b) by integrating the null geodesics once.
For each pixel (or each surface element): b → emitting point, angle α, redshift g, delay.
Intensity: I_obs = g³ I(E/g, α′), with the atmosphere’s beaming and the Doppler shift from rotation.
Add up: image, flux — and, phase by phase, the light curve and the phase-resolved spectrum.

Exact or approximate?

Exact integral (or elliptic integrals): always valid, secondary images included.
Beloborodov (2002) and Poutanen & Beloborodov (2006): closed formulae for ψ and for the delay, excellent for R ≳ 3 R_S.
Oblate Schwarzschild (Morsink et al. 2007): a star flattened by rotation, for millisecond pulsars.

This is the engine behind NICER’s mass and radius measurements — and behind the results of Part IV.

III

Part three

Pulsars and
neutron stars

The laboratories: more compact than anything short of a black hole, spinning, magnetised — and with their surface visible in X-rays.

pulsars and magnetarshot spotsrotation and atmosphereNICER

Pulsar (illustration) · NASA · public domain

Neutron stars in numbers

A solar mass in a city

Neutron star over Manhattan, to scale — A neutron star ~20 km across over Manhattan, to scale.

1.1–2.1 M☉

measured masses

11–13 km

radius (NICER, gravitational waves)

2–3 R_S

compactness: the regime of Part II

10⁸–10¹⁵ G

surface magnetic field

Central density a few times that of the atomic nucleus, ~10¹⁵ g/cm³. Spin periods from 1.4 ms to tens of seconds.

Surface gravity is ~10¹¹ times the Earth’s; the escape velocity, 0.5 to 0.7 c.

Cambridge · 28 November 1967

The first pulsar: “LGM-1”

Jocelyn Bell Burnell — **Jocelyn Bell**in 1967, PhD student

**Antony Hewish**supervisor · Nobel Prize 1974

Paper chart of the first pulsar, 1967 — The original chart, August 1967: the “scruff” pencilled in as CP 1919.

Bell noticed a bit of “scruff” that kept coming back on the radio telescope’s paper chart: pulses every 1.337 s, as regular as clockwork. Nickname: LGM-1, for Little Green Men. It was a spinning neutron star.

Pulsar theory

What is measured, and what is inferred

ATNF catalogue · 2,855 pulsars with measured Ṗhover over the points

observable	expression
period	P
characteristic age	t_c = P/2Ṗ
magnetic field	B = (3Ic³PṖ/8π²R⁶)^1/2
rotational energy lost	Ė_rot = −4π²IṖ/P³

Canonical neutron star

M = 1.4 M☉ · R = 10 km · I = 10⁴⁵ g cm²

B = 3.2×10¹⁹ (PṖ)^1/2 G

|Ė_rot| = 3.9×10⁴⁶ Ṗ/P³ erg/s

Only P and Ṗ are measured. The rest comes from assuming a spinning magnetic dipole — the straight lines in the diagram.

Transient activity

A magnetar, spinning

event	energy
short bursts (~0.1 s)	10³⁹ – 10⁴¹ erg
outbursts (months)	10⁴¹ – 10⁴³ erg
giant flares	10⁴⁴ – 10⁴⁶ erg

The field lines, anchored in the crust, get twisted over time. When the crust can no longer take the stress, it cracks — a starquake — and the magnetic energy comes out all at once.

Giant flares in the Galaxy and the Magellanic Clouds: 1979, 1998 and 2004.

The surface in X-rays

Why hot spots appear

Millisecond pulsars

Particles accelerated in the magnetosphere flow back and heat the polar caps to ~10⁶ K. They are NICER’s targets.

Magnetars

Field decay and the currents along twisted field lines heat regions of the crust — which change after each outburst.

Isolated cooling stars

The magnetic field makes heat conduction anisotropic: poles hotter than the equator (the “Magnificent Seven”).

NICER on the International Space Station — **NICER** · on the Space Station since 2017 · ~100 ns time resolution

XMM-Newton observatory — **XMM-Newton** · ESA, since 1999 · the largest area in soft X-rays

Chandra observatory — **Chandra** · NASA, since 1999 · the sharpest angular resolution

Ray tracing a hot spot

One spot, one light curve

static star · exact integral, Beloborodov or flat space

radius, R/R_S

spot ↔ axis, i

observer ↔ axis, j

angular radius of the spot

Two spots · Beloborodov’s classes

Antipodal, or not

two spots: one at each magnetic pole

magnetic axis ↔ rotation, i

observer ↔ rotation, j

offset of the second spot

Rotation and atmosphere

Doppler, aberration, delays and beaming

15° spot · 1.4 M☉, 12 km, spinning

spin frequency

Doppler: the surface moves at up to 0.1–0.2 c; the approaching side is brighter (bolometric flux ∝ δ⁵).
Aberration: the emission angle changes in the frame of the surface.
Delays: tens of μs between points on the star — fractions of the period at hundreds of Hz.
Beaming: an atmosphere does not emit equally in all directions; the profile depends on it.

Together they make the pulse asymmetric and create harmonics — and the asymmetry measures R separately from M.

NICER · pulse-profile modelling

Mass and radius from the pulse shape

Hot spots of PSR J0030+0451 inferred by NICER — PSR J0030+0451 · the solutions of the two teams: Amsterdam (Riley+ 2019, left) and Maryland (Miller+ 2019). In both, the spots are far from antipodal.

pulsar	mass (M☉)	radius (km)
PSR J0030+0451 · Riley+ 2019	1.34^+0.15_−0.16	12.7^+1.1_−1.2
PSR J0030+0451 · Miller+ 2019	1.44^+0.15_−0.14	13.0^+1.2_−1.1
PSR J0740+6620 · Riley+ 2021	2.07^+0.07_−0.07	12.4^+1.3_−1.0
PSR J0437−4715 · Choudhury+ 2024	1.42 ± 0.04	11.4^+1.0_−0.6

NICER time-stamps each photon to ~100 ns. Folding millions of photons at the period gives the pulse profile; ray tracing, with M, R, the geometry and the spots left free, fits that profile by Bayesian inference.

J0030’s spots came out non-antipodal and elongated: the field is not a centred dipole.

IV

Part four

Recent
results

Ray tracing applied to real data — magnetars, a massive magnetar, accretion columns — and the same tool used to test gravity beyond Einstein.

SGR J1745−29003XMM J1852+0033accretion columnsscalar-tensor gravitybraneworld

Magnetar (illustration) · ESO/L. Calçada · CC BY 4.0

2017 – 2026

Light curves as an instrument

2017

A&A 599, A87: rotation-powered SGRs and AXPs, with realistic neutron stars

2020

ApJ 889, 165: SGR J1745−2900 and a multipolar field, from the light curve

2023

Galaxies 11, 18: pulse profiles of accretion columns

2024

JHEAp 42, 52: 3XMM J1852+0033, a massive low-field magnetar

2024

EPJC 84, 1337: light curves as a test of scalar-tensor gravity

2024

PRD 110, 103014: stratified rotation and the moment of inertia

2025

PRD 112, 083010: quark stars and the QCD equation of state

2026

arXiv:2607.09836: compact stars in a high-tension braneworld

2026

JHEAp (submitted): the column of 1A 0535+262 through its NuSTAR spectrum

The common thread: reading in the light of a compact star its mass, its radius, its field — and the gravity that curves it.

de Lima, Coelho, Pereira, Rodrigues & Rueda · ApJ 889, 165 (2020)

SGR J1745−2900: more than two spots

Chandra light curves of SGR J1745-2900 with three-spot fits — Chandra light curves (2013 and 2016) and the three-spot fits.

Positions of the three spots in 2013 and 2016 — The three spots: two at the poles, one in the southern hemisphere.

The magnetar 0.1 pc from Sgr A*; 3.5 years of Chandra data.
Model: M, R, up to three spots of any size and position, Schwarzschild light bending; fitted with a genetic algorithm.
Two or three spots describe the data; with three, they are large (16–67°).
A field with non-trivial poloidal and toroidal structure — the same message NICER gave for J0030.

1.4 ± 0.1

M☉

10.9 ± 1.5

km

de Lima et al. · JHEAp 42, 52 (2024) · 1/2

3XMM J1852+0033: ray tracing and inference

One of the five pulse profiles (XMM-Newton, 2008) and the best fit: dashed, the contribution of each spot.

A transient magnetar near the remnant Kes 79, P = 11.56 s; until then there were only upper limits on the field.
Five XMM-Newton epochs over three weeks in 2008.
Ray tracing in Schwarzschild with circular spots of any size (Turolla & Nobili 2013) and magnetised atmospheres (NSMAX: H, C, O), with anisotropic emission.
Bayesian inference (MCMC): mass, radius, field and spot positions shared by all epochs; size and temperature free to vary.
Tested by injection and recovery of simulated data.

Best model

Three small spots, two in the southern hemisphere, and a carbon atmosphere: χ²/dof = 1.22, well ahead of two spots or of hydrogen.

de Lima et al. · JHEAp 42, 52 (2024) · 2/2

A massive, low-field magnetar

2.09^+0.16_−0.09 M☉

mass (1σ): among the highest ever estimated

12.0 ± 1.4 km

radius (2σ)

~8×10¹¹ G

surface field: log B = 11.89^+0.19_−0.93, below the critical value

Posterior distributions of mass, radius and field — Posteriors of M, R and log B.

Best-fit geometry: three spots — The best-fit geometry.

Interpretation: accretion of supernova material buried the field, which is now re-emerging; the carbon atmosphere would come from diffusive nuclear burning of H/He.

Compactness GM/Rc² ≈ 0.26: favours equations of state of intermediate stiffness.

Silva, de Lima, Coelho, Stecchini, de Liz & de Araujo · Galaxies 11, 18 (2023)

When the light comes from a column

3D mesh of the star with accretion columns — The star and the columns (truncated cone and cylinder) as a 3D mesh.

Sketch of the ray tracing with and without a periastron — Photons emitted above the surface may pass through a **periastron** before escaping.

1A 0535+262 (NuSTAR, 2018): observed profile and the two-column model.

In X-ray binaries, matter falls along the field lines and forms columns above the poles: the emission comes out sideways (fan beam).
Exact ray tracing in Schwarzschild, tabulated, including for points above the surface; Doppler included.
Genetic algorithm for the priors, then MCMC.

1.53^+0.52_−0.40

M☉

12.1^+1.1_−0.9

km

Silva, de Lima, Coelho & de Araújo · submitted to JHEAp (2026) · the model

Inside the column, through the spectrum

The same star as on the previous slide, now through its spectrum: NuSTAR, 3–44 keV, L ≈ 3.7×10³⁶ erg/s — 16 times below the critical luminosity.
Becker & Wolff (2022) model: a conical column; photons gain energy from the hot electrons and from the infalling gas (thermal and bulk Comptonization) and escape through the walls.
12 parameters, 4 independent MCMC chains of 256 walkers. The posterior has two modes; the dominant one wins by Δln 𝒫 ≈ 55.

v(y) = −v_ff [k_∞²(y³ − 1) + k₀²y⁴]^1/2y = r/R★ · free fall high up, braked near the surface

A submitted paper, led by Rodrigo Silva, with Jaziel Coelho and José Carlos de Araújo: the same star as on the previous slide, 1A 0535+262, now through its spectrum rather than its pulse profile. A NuSTAR observation at a luminosity sixteen times below the critical one: the subcritical regime. The model is Becker and Wolff’s: a conical column over the pole, gas falling along the field lines and braking near the surface, and photons gaining energy by scattering off the hot electrons and off the infalling gas itself, until they escape through the walls — the fan beam. The animation uses the geometry and velocity profile of the fit: an 11-degree cone, the top at four and a half stellar radii, and gas reaching the surface at a tenth of the speed of light. There are twelve parameters, sampled with four independent MCMC chains; the posterior has two modes, and the dominant one is far more probable.

Silva, de Lima, Coelho & de Araújo · submitted to JHEAp (2026) · the companion

A Be star: a disc at the equator

Illustration of a Be/X binary: the Be star with its disc, a stream of gas and the neutron star with an accretion disc — A Be/X binary: the Be star’s disc, the gas torn from it and the accretion disc around the neutron star.

Be stars: hot, of type B and late O, spinning close to break-up and losing gas from the equator, which forms a disc.
The “e” in O9.7–B0 IIIe stands for the hydrogen emission lines, such as Hα, which come from the disc; a double-peaked profile shows it rotates.
In 1A 0535+262 the companion is the neutron star, on a 111-day orbit that crosses the edge of the disc — the next slide.

Silva, de Lima, Coelho & de Araújo · submitted to JHEAp (2026) · the system

Where the column’s gas comes from

1A 0535+262: a neutron star pulsing every 103 s, with B ≈ 5×10¹² G, on an eccentric orbit — e ≈ 0.45, 111 days — around HD 245770, a Be star (O9.7–B0 IIIe), at 1.8 kpc.
The Be star spins fast and sheds a disc of gas from its equator.
At periastron the neutron star crosses the edge of the disc and steals gas: an X-ray outburst every orbit.
Close to it, the field cuts the disc at the magnetospheric radius and leads the gas to the poles: the columns.

In 2020 a giant outburst reached 1.2×10³⁸ erg/s; in quiescence, ~7×10³⁴ erg/s. The orbit has the system’s eccentricity and period; the stars, the disc and the close-up are schematic.

Where the gas falling into the column comes from. 1A 0535+262 is a Be/X-ray binary: the neutron star, pulsing every 103 seconds, moves on a very eccentric 111-day orbit around a Be star — a hot star spinning so fast that it sheds a disc of gas from its equator. The orbit in the animation has the right eccentricity and period, with Keplerian timing: slow far away, fast at periastron. It is at periastron that the neutron star crosses the edge of the disc and tears gas off it; the gas spirals into a small disc around it and falls in. Below, the light curve: one X-ray outburst per passage. In the close-up, side-on, what happens near the star: the disc is cut at the magnetospheric radius, where the 10¹² gauss field takes over, and the gas slides down the field lines to the poles — the columns of the previous slide. Sometimes a giant outburst comes, like the one in 2020, which reached 10³⁸ erg per second; the paper’s observation is from a quiet phase.

Silva, de Lima, Coelho & de Araújo · submitted to JHEAp (2026) · the Becker & Wolff model

One equation for the photons in the column

v ∂f∂R = ε3R² ∂(R²v)∂R ∂f∂ε + 1R² ∂∂R(R²c3n_eσ_∥ ∂f∂R) + n_eσ̄cm_ec² 1ε² ∂∂ε[ε⁴(f + kT_e ∂f∂ε)] − ft_esc + Q(R, ε) steady-state transport equation for the photon distribution, f(R, ε) · Becker & Wolff (2022), conical geometry

advection

the infalling gas carries the photons along, with the v(y) of the previous slide

bulk Comptonization

the flow converges: scattering off the infalling gas, photons gain energy

diffusion

photons move along the column, with the strong-field cross section σ_∥

thermal Comptonization

the Kompaneets term: electrons at T_e drive the spectrum towards ~kT_e

escape

photons leave through the walls in a time t_esc: the fan beam

sources

seed photons: bremsstrahlung, cyclotron and blackbody from the base

How it is solved

f_G = Σ_n C_n g_n(τ) (εkT_e)^κ−4 e^−ε/2kT_e M_κ,μ(ε_minkT_e) W_κ,μ(ε_maxkT_e)

Separation of variables: eigenfunctions g_n(τ) in optical depth and Whittaker functions in energy. The spectrum is this Green’s function convolved with the sources, adding the wall and the top of the column. Magnetic cross sections with vacuum polarization; cyclotron and iron lines, empirical. The implementation was checked equation by equation.

What the MCMC samples: 12 parameters

αvelocityξdiffusionψbulk/thermalθ₀openingk₀impacty_topheightT_eelectronsṀaccretionBfieldd_Feironf_syssystematicsDdistance

Fixed: M = 1.4 M☉, R★ = 10 km, N_H and the cyclotron line, outside the band. B has a prior from the line (46.3 keV). The likelihood is Gaussian, with the systematic term f_sys, and charges the energy budget L_acc = GMṀ/R★ within a factor of 2.

The details of the model. Becker and Wolff write a transport equation for the photon distribution in the column, f of R and ε, and solve it analytically. Term by term — the terms light up one at a time, and a click pins one. On the left, advection: the infalling gas carries the photons along. On the right, first bulk Comptonization: the flow converges, and photons scattering off the infalling gas gain energy, as in first-order Fermi acceleration. Then diffusion along the column, with the cross section reduced by the field. Then the Kompaneets term, thermal Comptonization, which drives the spectrum towards the electron temperature. Escape through the walls, which gives the fan beam. And the sources: bremsstrahlung, cyclotron and the blackbody from the base of the column. The solution separates variables — spatial eigenfunctions and Whittaker functions in energy — and the spectrum is the Green’s function convolved with the sources. In the fit this becomes twelve parameters: the temperature, the accretion rate, the cone geometry, the field, with a prior from the cyclotron line, and a free systematic uncertainty.

Silva, de Lima, Coelho & de Araújo · submitted to JHEAp (2026) · results

A hot, tall column — and a missing piece

NuSTAR spectrum of 1A 0535+262 with the Becker and Wolff model and residuals — The NuSTAR spectrum and the posterior (median and 1σ and 2σ bands), with the residuals. Inset: the iron line.

8.56×10⁷ K

electron temperature, to 1% · kT ≈ 7.4 keV

4.5 R★

top of the column, ~35 km above the surface

0.070

ψ, bulk/thermal: thermal Comptonization dominates

log θ₀ − 1.36 log Ṁ
= −22.32 ± 0.10

the spectrum fixes the combination: more Ṁ, a wider column

A missing piece

The field comes out 2.3σ above that of the cyclotron line — which would then sit at 55 keV, not the measured 46 — and the residuals move together (runs test: z = −3.7). Candidate: recombination between excited Landau levels, seen in this source and absent from the model.

Dominant posterior mode. χ²_ν = 1.00 with a free systematic uncertainty is not a goodness-of-fit test: the residuals are the diagnostic.

The results. What the spectrum measures well is the electron temperature of the column: 8.6 times 10 to the 7 kelvin, to one per cent — kT of 7.4 keV — and the four chains agree with each other. Thermal Comptonization dominates: the parameter measuring bulk against thermal comes out at 0.07. Geometry and accretion rate are degenerate: more infalling matter can be compensated by a wider column, and the spectrum fixes the combination of the two. The top of the emitting region sits at four and a half stellar radii, some 35 kilometres above the surface — tall for a subcritical source, where a radiation-dominated shock is not expected. And the fit shows a missing piece: the magnetic field comes out stronger than the cyclotron line implies, and the residuals are correlated. One candidate is recombination emission between excited Landau levels, seen in this source and absent from the model. A reduced χ² of one proves nothing here, because the systematic uncertainty is free.

Brans & Dicke 1961 · Damour & Esposito-Farèse 1993

What if gravity is not Einstein’s?

What “coupling to matter” means

g̃_μν = A²(φ) g_μνthe metric matter feels

Matter measures distances and times with g̃: its rulers and clocks change with φ. Light does not notice — a conformal change leaves null geodesics alone — but the radius matter measures is R = A_sR_E.

α(φ) = d ln Adφthe strength of the coupling

Matter sources the field, □φ = −4πGαT, and feels an extra force, ∝ α∇φ: between two bodies, G_eff = G(1 + α²).

In scalar-tensor theories, besides the metric there is a scalar field φ that couples to matter. In the Solar System the coupling is tiny (α² < 10⁻⁵, from the Cassini probe); in a neutron star spontaneous scalarisation can occur: with A = e^βφ²/2 and β ≲ −4.35, φ = 0 becomes unstable and the field grows on its own, like spontaneous magnetisation.

Outside the star, the metric is no longer Schwarzschild: it is the Just metric, with a scalar charge Q:

ds² = −(1 − ar)^b/adt² + (1 − ar)^−b/adr² + (1 − ar)^{1 − b/a}r²dΩ²b = 2GM/c² · a/b = √(1 + Q²) · Einstein frame

Two stars with the same mass and the same radius — “doppelgängers” — bend light in different ways. The light curve can tell the theories apart.

Now gravity. In scalar-tensor theories, gravity has a scalar field in addition to the metric. Brans and Dicke proposed the first one in 1961; Damour and Esposito-Farèse showed in 1993 that, in neutron stars, the field can grow spontaneously — scalarisation — even if in the Solar System everything looks like general relativity. Outside a scalarised star the metric is the Just metric, with a scalar charge Q. The question of our work: do two stars with the same mass and the same radius, one in general relativity and the other scalarised, have different light curves? What “coupling to matter” means: matter does not measure distances and times with the metric g, but with g̃ = A²(φ) g — rulers and clocks made of matter change with the field. The strength of that coupling is α, the derivative of ln A: matter sources the field, and the field pulls on matter with an extra force, so the effective gravity between two bodies becomes G(1 + α²). Light does not notice the conformal change, but the radius matter measures is A_s times the Einstein-frame one — and it is this A_s, 0.95 or 1.05, that changes on the next slide. In the Solar System the Cassini probe limits α² to below 10⁻⁵; with A = exp(βφ²/2), however, if β is more negative than −4.35, inside a neutron star φ = 0 becomes unstable and the field grows, like the spontaneous magnetisation of a ferromagnet.

Ottoni, Coelho, de Lima, Pereira & Rueda · EPJC 84, 1337 (2024)

The same star, different gravity

M = 2.1 M☉ · one 10° spot · ι = θ = 90° · static star · exact integral

compactness, GM/Rc²

scalar charge, Q

This panel redoes the paper’s calculation. A star of 2.1 solar masses with compactness 0.284 — the value at which, in general relativity, the critical deflection angle is exactly 180 degrees: the point opposite the observer sits right at the edge of visibility. A 10-degree spot passes behind the star. In general relativity, the curve has a smooth minimum. In scalar-tensor theory 1, with A_s = 0.95, the critical angle drops to 0.90π: an invisible zone of 18 degrees opens up behind the star, and the spot is eclipsed — the flux goes to zero. In theory 2, with A_s = 1.05, the angle exceeds π and a zone of multiple images appears: a lensing peak in the middle of the minimum. Same mass, same radius, qualitatively different curves. Play with the scalar charge: at Q = 0 the three coincide.

Ottoni et al. · EPJC 84, 1337 (2024) · results

Up to 80% difference in flux

Critical deflection angle versus compactness — The critical angle versus M/R: the GR line and the scalar-tensor band.

Flux at 700 Hz and relative difference — Spinning at 700 Hz, with Doppler and delays: flux and difference relative to GR.

The compactness at which ψ_c = π ranges from 0.275 to 0.305 in scalar-tensor theory, against 0.284 in GR.
The flux from a spot differs from GR by up to ~80%.
The radius changes by up to ~10% for the same stellar parameters: today, this limits Q ≲ 0.5.
With NICER, J0740+6620 has C = 0.21–0.28: not yet enough to separate the theories.

Outlook

Radii with smaller uncertainties — eXTP, NewAthena, and gravitational waves from Einstein Telescope and LISA — may constrain scalar-tensor theories through the pulse shape.

Randall & Sundrum 1999 · Shiromizu, Maeda & Sasaki 2000

The Universe on a membrane

In a braneworld, our Universe is a four-dimensional hypersurface — the brane — inside a space with one extra dimension, the bulk. Matter and light are stuck on the brane; only gravity propagates in the bulk.

The idea was put forward to attack the hierarchy problem (why gravity is so weak) and the cosmological-constant problem and accelerated expansion.

G_μν = 8πT_μν + 48πλ S_μν − E_μνShiromizu, Maeda & Sasaki 2000 · κ₅⁴ = 48π/λ (G = c = 1)

S_μν: local corrections, quadratic in the matter (∝ T²). In a fluid they add ρ²/2λ to the density: they matter when ρ approaches the brane tension λ.
E_μν: the bulk Weyl tensor projected onto the brane — nonlocal: it cannot be computed from brane quantities alone.
For λ → ∞ and E_μν → 0, General Relativity is recovered.

In the paper: an effective 4D closure, inspired by the Sahni–Shtanov scenario (2003) — not a complete 5D bulk solution.

What a braneworld is. The idea, from Randall and Sundrum in 1999, is that our four-dimensional Universe is a membrane — the brane — embedded in a space with one more dimension, the bulk. All Standard Model fields, matter and light, are stuck on the brane; only gravity spreads into the bulk. It was proposed to explain why gravity is so much weaker than the other forces and, later, for dark energy. Shiromizu, Maeda and Sasaki showed what Einstein's equations look like from inside the brane: two new terms appear. A local one, quadratic in the matter, with coefficient 48π/λ, where λ is the brane tension: in a fluid it adds ρ²/2λ to the density, and it only matters when ρ is comparable to λ. And a nonlocal one, the projection of the bulk Weyl tensor: the gravitational field of the extra dimension seen from here, which cannot be computed from what lives on the brane alone. In our paper we use an effective, four-dimensional version inspired by the Sahni–Shtanov scenario.

The modified TOV equations · SLy equation of state

Same nuclear physics, different gravitation

The equation of state of the matter (SLy) does not change. What changes is how much it gravitates: the star now feels an effective density and pressure.

ρ_eff = ρ + ρ²2λ + ρ_𝒰
p_eff = p + pρλ + ρ²2λ + p_𝒰local term (λ) · Weyl sector (𝒰)

dmdr = 4πr² ρ_eff
dpdr = −(ρ + p) m + 4πr³p_effr(r − 2m)matter is conserved: the (ρ + p) factor is the usual one

The local term strengthens attraction and shrinks the star — but only when ρ approaches λ. At the centre of a neutron star ρ ≈ 2×10⁻³ km⁻²; with λ ~ 10⁴ km⁻², ρ/λ ~ 10⁻⁷.

modified TOV · SLy

log₁₀(λ / km⁻²)

α_𝒰

w_𝒰

How this enters a star. The nuclear physics does not change: the equation of state is still SLy, the same as in general relativity. What changes is how much the matter gravitates. The mass equation gets an effective density, and the numerator of the TOV equation an effective pressure — each with a local term, controlled by the tension λ, and a Weyl term. The (ρ + p) factor stays as usual because ordinary matter is conserved. The panel solves these equations here, in the browser, with SLy. The "local term only" button sets λ = 10⁻³ km⁻²: the effective density goes up, gravity squeezes harder, and the maximum mass collapses. But for that λ has to be comparable to the central density of the star, about 2×10⁻³ km⁻². With the λ the data prefer, 10⁴, the local term is ten million times smaller than ρ — it is the Weyl sector that moves the curve. "Best fit" is the paper's best sample: a maximum mass of 2.28 and a radius of 13.2 km, against 2.05 and 11.7 in general relativity.

What comes from the bulk

The Weyl sector: screening or extra attraction

E_μν is the imprint of the bulk's gravitational field on the brane. There it acts as an effective fluid — in cosmology, the “dark radiation”. But the brane does not determine this fluid: a closure is needed.

ρ_𝒰 = α_𝒰 ρ

p_𝒰 = w_𝒰 ρ_𝒰

Isotropic closure; the anisotropic Weyl stress (Δ_𝒰) is left for future work.

if…	effect on the star
α_𝒰 < 0	ρ_eff = ρ(1 + α_𝒰) < ρ: screening — less mass gravitates, larger and heavier stars
α_𝒰 > 0	more attraction: smaller and lighter stars
w_𝒰	only appears as w_𝒰α_𝒰ρ in the effective pressure: it matters little when α_𝒰 is small, and it can compensate α_𝒰

M–R curves with λ = 10⁴ km⁻² · SLydashed: GR

The Weyl sector is the interesting part. The projected Weyl tensor is the gravitational field of the bulk seen from the brane, and on the brane it behaves as an extra fluid — in cosmology it shows up as “dark radiation”. The problem is that it is not determined by the brane equations: a closure has to be chosen. Ours is the simplest one: a Weyl density proportional to the matter density, with coefficient α, and a pressure with equation of state w. On the left of the panel, the fan in α with w = 0: negative α reduces the density that gravitates — a screening — and the star grows; positive α squeezes it. On the right, the fan in w with α = −0.15: w only enters multiplied by α, so different pairs of α and w give similar curves. That degeneracy is what leaves w poorly determined by the data.

dos Santos Júnior, de Lima & Moraes · arXiv:2607.09836 (2026, submitted) · results

High tension, mildly negative Weyl coupling

Data and inference

GW170817 (radii from the tidal deformability) + NICER: PSR J0740+6620 and PSR J1231−1411 — both of them results of pulse-profile ray tracing. Ensemble MCMC (4 × 32 walkers × 3,000 steps, 288,000 samples), priors centred at large λ and at the GR values.

	median	68%	GR
log₁₀(λ/km⁻²)	3.98	2.54 – 5.50	∞
α_𝒰	−0.15	−0.23 – 0.15	0
w_𝒰	0.24	−0.78 – 0.61	0
M_max (M☉)	2.30	2.22 – 2.44	2.05
R_1.4 (km)	13.31	12.74 – 13.85	11.76

Braneworld mass–radius band against GR — The posterior M–R band (68% and 95%), GR+SLy (dashed) and the constraints used.

λ ~ 10⁴ km⁻²: the local term vanishes (ρ/λ ~ 10⁻⁷). The posterior nearly repeats the prior: the data are compatible with high tension, they do not measure λ.
α_𝒰 < 0 in 83% of the samples: a mild screening. It is the Weyl sector, not the local term, that pushes the curve to the right.
A secondary branch (17%): α_𝒰 > 0 and w_𝒰 < 0 compensate — smaller R_1.4, larger M_max, compactness near the imposed bound.
M_max ≈ 2.3 M☉: consistent with PSR J0952−0607 (2.35 ± 0.11 M☉). GW190814 (2.50–2.67 M☉) only fits in the 95% tail.
The message: heavier stars are not a generic prediction of braneworlds — they depend on the sign of the Weyl sector.

The results. We used three datasets: the masses and radii of GW170817, inferred from the tidal deformability, and two NICER pulsars — J0740 and J1231 —, whose radii come, themselves, from pulse-profile ray tracing, the subject of this seminar. The sampling is an ensemble MCMC with 288 thousand samples. The brane tension comes out around 10⁴ km⁻², but honesty is needed: the posterior nearly repeats the prior; what the data say is that high tension is compatible, not that we measured λ. The Weyl coupling comes out mildly negative, in 83% of the samples, which gives the screening: the M–R curve moves to the right, with a maximum mass of 2.3 and a radius of 13.3 km. The parameter w is barely constrained. There is a secondary branch, with 17% of the samples, where positive α and negative w compensate. The maximum mass is consistent with the pulsar J0952, at 2.35; the secondary of GW190814 only enters in the tail. And the main message: in braneworlds, heavier stars do not come for free — they depend on the sign of the Weyl sector. Numbers from arXiv version 1; the paper is under review.

Itoh 1970 · Bodmer 1971 · Witten 1984 · the strange matter hypothesis

What if iron is not the bottom of the well?

At extreme densities neutrons may dissolve into quarks. The Bodmer–Witten hypothesis goes further: matter made of u, d and s quarks would be the ground state of matter — with an energy per baryon below that of ⁵⁶Fe, 930 MeV.

Pauli: each flavor has its own Fermi sea. A third flavor opens new levels and the top (μ) drops — even after paying for the s mass.

Nuclei don't decay because that would take ~A simultaneous weak d → s conversions.
And u–d matter alone must be unstable — otherwise iron itself would have turned into quarks.
If it holds, “neutron stars” may be quark stars, bound by the strong force.

the stability window · energy per baryon at p = 0

a₄ (interaction)

a₂ (10³ MeV²)

The strange matter hypothesis of Itoh, Bodmer and Witten. In a nucleus, quarks are locked inside protons and neutrons. But if matter were squeezed until the nucleons dissolved, free u and d quarks could turn some of themselves into s quarks through the weak interaction. Why? Pauli: each flavor has its own Fermi sea. With three flavors, the same number of quarks fits into lower levels and the energy drops — even after paying for the s mass. The hypothesis is that this three-flavor matter has a lower energy per baryon than iron: it would be the ground state of matter. Nuclei do not decay because that would take dozens of simultaneous weak conversions. And u–d matter alone must be unstable, or iron would already have turned into quarks. The panel shows this with the model of our paper: the amber curve is the energy per baryon of uds matter, the lilac one that of ud matter. The turquoise band is the window where uds lies below iron and ud above it: there, quark stars exist. To the left, iron would decay; to the right, strange matter only appears under pressure, inside hybrid stars.

Alford, Braby, Paris & Reddy 2005 · a QCD-informed equation of state

Dense QCD in three numbers

ω = −3a₄4π² μ⁴ + 3a₂4π² μ² + Bthermodynamic potential · p = −ω · μ: quark chemical potential

a₄ = 1 − c	the strong interaction. a₄ = 1 means free quarks; perturbative QCD suggests c ≳ 0.3
a₂ = m_s² − 4Δ²	the s-quark mass against the gap Δ of color superconductivity (CFL phase) — it can be negative
B	the bag constant: the cost of opening the QCD vacuum

At the surface p = 0, but the density is not: ρ_s ≈ 4B ~ 10¹⁴–10¹⁵ g/cm³. The star is self-bound: at small masses M ∝ R³, unlike a neutron star.

TOV + the tidal equation with the surface jump (Δy = −4πR³ρ_s/M), solved here. The paper's example gives 2.19 M☉, 11.97 km and Λ = 852 — the panel reproduces it.

TOV + tides · u, d, s quarks

a₄

a₂ (10³ MeV²)

B^1/4 (MeV)

To build stars we need an equation of state. The one by Alford and collaborators starts from QCD at very high density and adds corrections: a potential with three parameters. a4 measures the strong interaction: it is 1 for free quarks and about 0.7 if perturbative QCD is right. a2 combines the strange quark mass and the color-superconducting gap — in the CFL phase, where quarks of different colors and flavors pair up, a2 is m_s² minus four times the gap squared, and it can go negative if the gap is large. B is the bag constant, the cost of opening a bubble of perturbative vacuum. The panel solves the TOV and tidal equations for this matter. Look at the shape of the curve: it starts at the origin. A small quark star is a droplet of strange matter held by the strong force, and its radius grows as the cube root of the mass — the opposite of a neutron star, dashed. “Free quarks” does not reach 2.13 solar masses; the paper's example passes every test; and the “large gap” case, with negative a2, passes too — and goes through the compact object in HESS J1731, the pink point.

Zastrow, Pereira, de Lima & Horvath · PRD 112, 083010 (2025) · the method

Thousands of stars, drawn at random

Each point (a₂, a₄, B) inside the stability window is a whole family of stars: TOV for M and R, the tidal equation for Λ. The paper drew 263,287 points by Latin hypercube sampling; here the browser draws its own, live.

	accepted range	where it comes from
M_max	2.13 – 2.33 M☉	pulsars in binaries (Romani et al. 2022) · GW170817 + GRB 170817A (Rezzolla et al. 2018)
R_1.4	10.00 – 13.25 km	NICER radii — ray tracing — and gravitational waves (compilation by Pang et al. 2023)
Λ_1.4	< 970	tidal deformability from GW170817

Amber: the families that pass all three tests. The dashed lines show how far the accepted ones reach — and they converge to the paper's limits.

quasi-random sampling · live

The method of the paper by João Zastrow, at UDESC. The parameter space has three dimensions; each point inside the stability window defines an equation of state and, with it, a whole family of stars — maximum mass, radius at 1.4 solar masses and tidal deformability. We drew 263 thousand points by Latin hypercube sampling, which covers the space evenly without the cost of an MCMC; the price is that there are no credible intervals, only limits. The constraints are conservative: maximum mass between 2.13 and 2.33, radius between 10 and 13.25 km — and those radii come from NICER, that is, from ray tracing — and Λ below 970, from GW170817. The panel redoes the method here: it draws points, solves the stars and paints in amber the ones that pass. Switch the axes: maximum mass against a4 shows the upper cut near 0.8; against B, the bag-model shape, M_max ∝ B^−1/2. With a few thousand points the limits already approach those of the paper.

Zastrow, Pereira, de Lima & Horvath · PRD 112, 083010 (2025) · results

Interacting quarks, a dense surface

combined constraints	range
a₄	0.18 – 0.81
a₂ (10³ MeV²)	−202 – 42
B^1/4 (MeV)	114 – 185
R_1.4 (km)	10.00 – 12.26
Λ_1.4	270 – 970
ρ_surface (10¹⁵ g/cm³)	0.29 – 0.60
c_s/c	0.51 – 0.85

The sound speed exceeds c/√3 ≈ 0.58 — the conformal limit of QCD at infinite density — when a₂ < 0.

Allowed mass-radius curves for quark stars, with XMMU J173203.3−344518 — All M–R curves that pass the constraints, with XMMU J173203.3−344518 (orange).

a₄ ≲ 0.81: the matter departs by at least 20% from free quarks. Strong interactions cannot be neglected, at any mass.
Δ ≲ 230 MeV (with m_s = 100 MeV), from the lower limit on a₂ — in line with the independent bound of ≲ 216 MeV (Kurkela et al. 2024).
The surface jumps from vacuum to 1.1–2.2 ρ_sat: a number that enters the tides, the oscillation modes and the star's mountains.
R_1.4 = 10.0–12.3 km, like a neutron star: the radius alone cannot tell them apart; the tides can.
XMMU J1732−344 (0.77 M☉, 10.4 km, in HESS J1731−347) falls on the curves effortlessly — a quark-star candidate, if the atmosphere model behind its radius is right.

The results. Combining the three constraints, a4 stays below 0.81: quark matter in a star departs by at least 20% from a free gas, at any mass. The lower limit on a2, minus 202 thousand MeV², with a 100 MeV strange quark, gives a color gap of at most 230 MeV — and an independent group, with another method, found 216. The bag constant lands in a narrow range, which fixes the surface density between 1.1 and 2.2 times nuclear saturation: that jump is a boundary condition for tides, quasi-normal modes and ellipticity. The radius of a 1.4 solar-mass quark star lies between 10 and 12.3 km — the same as a neutron star; what can tell them apart is the deformability. And the compact object in HESS J1731, with 0.77 solar masses and 10.4 km — too light to come out of an ordinary supernova — falls naturally on these curves. But careful: that radius comes from spectral fitting, and depends on the atmosphere model and the temperature map — the same issues we discussed for hot spots.

Pereira, Ottoni, Coelho, Rueda & de Lima · PRD 110, 103014 (2024)

When core and crust do not spin together

Almost every calculation treats a pulsar as a rigid body. But the interior may be superfluid, have a phase transition — a quark core — or gain accreted layers that spin differently, like the Sun's tachocline.

1r⁴ddr(r⁴j ω̄′) + 4j′r ω̄ = 0frame dragging in slow rotation · Hartle 1967 · ω̄ = Ω − ω

[ω̄]_R★ = [Ω], [ω̄′]_R★ = 0 · I = JΩ⁺the jump at R★, derived in the paper · Ω⁺: the surface spin, the one timing measures

J comes from the derivative of ω̄ at the surface, so it feels the jump: the same star, with the same Ω⁺, has a different moment of inertia.

Hartle + SLy · two-layer rotationlong dash: 1.0 M☉ · dotted: 2.0 M☉

mass

jump 𝒞 = [Ω]/ω̄⁻

R★ / R

The other paper is about rotation. Almost every pulsar model assumes the star spins as a rigid body. But the interior may be superfluid and weakly coupled to the crust; there may be a phase transition, with a quark core that, by angular momentum conservation, spins differently; or accreted layers from outside. The Sun has this: the tachocline separates an interior spinning rigidly from an envelope spinning differently. In general relativity, slow rotation is described by Hartle's equation for frame dragging. What the paper does is derive the jump conditions when Ω jumps at a radius R-star: omega-bar jumps with it, its derivative does not. Since the angular momentum J comes from the derivative at the surface, it feels the jump, and the moment of inertia, defined with the observed surface spin, changes. The panel solves this here with SLy. Drag R-star: if the jump is deep, below about 40% of the radius, almost nothing changes. Near the surface, the change in I reaches the change in Ω itself. The star on the left spins with both real rates, no exaggeration.

Pereira, Ottoni, Coelho, Rueda & de Lima · PRD 110, 103014 (2024) · results

Rigid rotation is an assumption — and it has a cost

Near the surface, ΔI/I → ΔΩ/Ω⁺; below ~40% of the radius the effect vanishes. It grows fast for R★ ≳ 0.6 R, and more in light stars.
Hybrid stars: the jump sits at the transition to the quark core, and only exceeds 5% with 𝒞 ~ 20% and R★ ≳ 0.8 R (figure).
In the radius: with a known mass, ΔR/R = ½ ΔI/I. Ten percent in I is 5% in R — the size of the NICER error bars, which come from ray tracing.
Direct measurement of I: the double pulsar J0737−3039A should give I to 11% by 2030 (today, I₄₅ = 1.15–1.48). Jumps of 10–20% in Ω cannot be told from rigid rotation.
Energy: with I 10–20% larger, five Fermi gamma-ray pulsars with L_γ > Ė fit back within the rotational energy.
I-Love-Q: the quadrupole depends on another function besides J — the universal relations weaken. And in CCOs, accreted layers may spin differently for up to ~10 thousand years.

Change in the moment of inertia of hybrid stars against mass — Hybrid stars (η = 0.39 and 0.77): ΔI/I with the jump at the phase transition, against mass.

pulsar (Fermi)	L_γ/Ė, I₀	+10%	+20%
J1057−5851	1.19	> 1	0.99
J1522−5735	1.09	0.99	0.91
J1650−4601	1.11	> 1	0.92

What this changes. Near the surface, the relative change in the moment of inertia tends to the relative change in rotation, and below about 40% of the radius it vanishes. In the hybrid stars of the figure, the jump sits at the transition to the quark core, which is deep: it only exceeds 5% with 20% jumps in rotation. Now the link with the rest of this talk: if the mass is known and the radius is inferred from the moment of inertia, a 10% error in I becomes 5% in the radius — the size of the NICER error bars, which come from ray tracing. The double pulsar should give the moment of inertia to 11% by 2030, and that will not separate a rigid star from one with 10 to 20% jumps. In the energy budget, an I larger by 10 to 20% lets some gamma-ray pulsars, which seemed to emit more than rotation allows, fit back within it — the table shows three, for 1.45 solar masses; for 2.04 there are two more. And the universal I-Love-Q relations, which assume rigid rotation, weaken.

Conclusions

From geometry to observation — and back

In a neutron star, light bending, redshift and delays are not corrections: they define the image and the light curve.
Backward ray tracing — with the exact integral when the compactness calls for it — links the metric to what a telescope measures.
Applied to Chandra, XMM-Newton and NuSTAR data, it gave the mass, radius, field and spot geometry of magnetars and of an accreting star — and, through its spectrum, the temperature and height of its column.
Since the light curve depends only on the exterior metric, it also tests gravity: scalar-tensor theories and braneworld.
And the interior matters: quark stars and stratified rotation change the M, R and I that ray tracing measures.

Next steps

Magnetised atmospheres and exact ray tracing together, for isolated cooling stars.
Fast rotation and metrics beyond Schwarzschild in the same fits.
Data from the next generation: eXTP, NewAthena, STROBE-X.

The panels in this seminar run in the browser: rafael-lima.pages.dev/talks.

Image credits

Freely licensed photos and figures

Albert Einstein — Ferdinand Schmutzer / Adam Cuerden · public domain

Karl Schwarzschild — unknown author · public domain

M87* (EHT) — Event Horizon Telescope Collaboration · CC BY 4.0

Sgr A* (EHT) — EHT Collaboration · CC BY 4.0

Jocelyn Bell, 1967 — Roger W Haworth · CC BY-SA 2.0

Antony Hewish — Woodruff T. Sullivan III / NRAO/AUI Archives (NRAO/AUI/NSF) · CC BY 4.0

Chart record of CP 1919, 1967 — Bill Thompson (Billthom) · CC BY-SA 4.0

Neutron star over Manhattan — NASA's Goddard Space Flight Center · public domain

1919 eclipse — ESO/Landessternwarte Heidelberg-Königstuhl/F. W. Dyson, A. S. Eddington, & C. Davidson · CC BY 4.0

SMACS 0723 (JWST) — NASA, ESA, CSA, and STScI · public domain

NICER — NASA · public domain

XMM-Newton — ESA/XMM-Newton · CC BY-SA 3.0 IGO

Chandra — NASA/CXC/NGST · public domain

Pulsar (illustration) — NASA · public domain

Magnetar (illustration) — ESO/L. Calçada · CC BY 4.0

Carl Brans — Cfbolner · CC BY-SA 4.0

Thibault Damour — Melirius · CC BY 3.0

PSR J0030+0451: hot-spot models — NASA's Goddard Space Flight Center · public domain

Paper figures (Part IV) — de Lima et al., ApJ 889, 165 (2020); de Lima et al., JHEAp 42, 52 (2024, CC BY 4.0 on arXiv); Silva et al., Galaxies 11, 18 (2023, CC BY 4.0); Ottoni et al., EPJC 84, 1337 (2024, CC BY 4.0); dos Santos Júnior, de Lima & Moraes, arXiv:2607.09836 (CC0); Pereira et al., PRD 110, 103014 (2024); Zastrow et al., PRD 112, 083010 (2025); Silva, de Lima, Coelho & de Araújo (submitted to JHEAp, 2026)

Simulated Sobral plate — positions from the Tycho-2 catalogue (Høg et al., A&A 355, L27, 2000), via VizieR/CDS

Be/X binary (illustration) — AI-generated

P–Ṗ diagram — ATNF pulsar catalogue 2.8.1 (Manchester et al., AJ 129, 1993, 2005)

Thank you

Questions?

Slides, animations and dashboards available at rafael-lima.pages.dev/talks

From geometryto observation

The star inside

Relativistic ray tracing

Pulsars and neutron stars

Recent results

The direct problem

The inverse problem

The test

The starinside

pressure has inertia

pressure gravitates

space curves

Relativisticray tracing

1 · It bends the path

2 · It changes the energy

3 · It delays

Exact or approximate?

Pulsars andneutron stars

Canonical neutron star

Millisecond pulsars

Magnetars

Isolated cooling stars

Recentresults

Best model

advection

bulk Comptonization

diffusion

thermal Comptonization

escape

sources

How it is solved

What the MCMC samples: 12 parameters

A missing piece

What “coupling to matter” means

Outlook

Data and inference

Next steps

Questions?

Notes

From geometry
to observation

The star
inside

Relativistic
ray tracing

Pulsars and
neutron stars

Recent
results