Recompute a cluster sightline from saved inputs

Download three files: the saved inputs, the script and the saved outputs. The script needs Python 3 and its standard library only. The outputs were produced with Python 3.9.6. The physics is the two-state mixing of Photon and axion as two mixing states and the random-cell model of Random domains, polarization and stellar cooling.

Step 1: read the inputs and their units

Open the inputs file first. The units block fixes every unit before any number appears: photon energy in keV, coupling in inverse GeV, field in microgauss, electron density per cubic centimetre, axion mass in eV, lengths in kpc and angles in radians. Seven energies from 0.6 to 9 keV are listed once and shared by all cases.

There are three cases. The uniform case holds one field direction split into 12 cells of 8 kpc. The tangled case is the Lab's default sightline, 24 cells with fixed but irregular angles. The resonance case sets the axion mass equal to the plasma frequency of the gas. Each case stores the cell angles explicitly. The Lab draws its angles from a seeded generator. Saving the angles lets anyone rerun the case with any random-number library.

Step 2: run the script

python3 -I recompute_sightline.py sightline-inputs.json my-outputs.json

The -I flag keeps Python from importing modules from the current folder. Compare your file with the saved outputs. They should agree to the last few digits. The Lab's own transfer matrices, computed in a different language and a different way, agree with these outputs to better than $10^{-10}$ in survival probability. That agreement is checked by a repository test.

Step 3: what one cell does

Each cell rotates the state into the field frame, evolves the parallel photon and the axion in the eigenbasis of their two-by-two matrix, gives the perpendicular photon its plasma phase, and rotates back. This is the central function:

def evolve_cell(state, plasma, axion, mixing, length, angle):
    ax, ay, a = state
    c, s = math.cos(angle), math.sin(angle)
    parallel = c * ax + s * ay
    perpendicular = (-s * ax + c * ay) * cmath.exp(1j * plasma * length)
    theta = 0.5 * math.atan2(2 * mixing, plasma - axion)
    mean = 0.5 * (plasma + axion)
    half_split = 0.5 * math.hypot(plasma - axion, 2 * mixing)
    ct, st = math.cos(theta), math.sin(theta)
    upper = (ct * parallel + st * a) * cmath.exp(1j * (mean + half_split) * length)
    lower = (-st * parallel + ct * a) * cmath.exp(1j * (mean - half_split) * length)
    parallel, a = ct * upper - st * lower, st * upper + ct * lower
    return [c * parallel - s * perpendicular, s * parallel + c * perpendicular, a]

The three arguments plasma, axion and mixing are $\Delta_{\rm pl}$, $\Delta_a$ and $\Delta_{a\gamma}$ in inverse kiloparsecs. The script derives them from CODATA constants. One kiloparsec is $1.5637\times10^{26}\,\mathrm{eV^{-1}}$, one tesla is $195.35\,\mathrm{eV^2}$ and the plasma frequency is $3.713\times10^{-11}$ eV at one electron per cubic centimetre. These three values are written into the outputs under derivedConstants, so a reader can see which constants produced a given file.

Unpolarized light is handled by running the two pure polarizations separately and averaging their density matrices with weight one half each. The linear polarization degree is

Step 4: the checks that must pass

Probability. Photon survival plus axion fraction must equal one. The largest deviation over all 21 saved rows is $3.6\times10^{-15}$, which is rounding error.

Uniform case against the formula. For one uniform field the survival of unpolarized light is $1-P_{\gamma\to a}/2$ with the two-state probability. Every uniform row agrees with this closed form to about $10^{-15}$. At 6.6 keV both give 0.5572. At 3 keV both give 0.9729, which sits near a node of the oscillation.

Resonance flatness. With the axion mass equal to the plasma frequency, survival is 0.7761 at all seven energies.

Monte Carlo mean. The tangled case's saved angles are one realization. To test the ensemble rule, the script draws 2000 new sets of 24 uniform angles with Python's own generator, seed 79, and averages the axion fraction at 3 keV:

A difference of 1.3 standard errors is consistent with chance. A pull above about three would point to a bug in either the propagation or the formula.

Step 5: four failure cases

The outputs also store four deliberate mistakes, each with a number that shows the damage.

The Euler case is the most common in practice. Any integrator that does not preserve the norm will slowly create or destroy photons, and the error grows with path length. The exact cell transfer avoids this. A fourth-order Runge–Kutta integrator with a 2 pc step, used as an independent check in the repository tests, agrees with it to $10^{-6}$.

Step 6: recompute a sightline exported from the Lab

The Lab's export button writes a JSON record with format sandboxphysics.axion-photon-conversion.v1. Its input block stores logarithms for coupling, density and mass, with keys such as logCouplingPerGev. Convert each with ten to that power and copy the domains list unchanged. Then add the result as a new case in your inputs file. The export also stores spectrum.simulatedSurvival and, separately, spectrum.idealReference. Compare your recomputed rows with the simulated curve. The ideal reference is a different object: the closed formula in uniform mode, or the ensemble mean in tangled mode.

Step 7: read the frozen limits with their assumptions

The Lab's stellar ruler shows four published limits. They are copied as printed and not recomputed, so each carries its authors' assumptions.

The X-ray value is eq. (5) of Reynolds et al., 2020. The gamma-ray value is from the abstract of Ajello et al., 2016, which we did not read in full. The CAST and stellar values are eqs. (89.20) and (89.16) of the PDG review of axions, 2025. The stellar mass edge in the Lab is set conservatively at 1 keV, because Ayala et al., 2014 discuss particles below a few keV but print no sharp cutoff.

The confidence levels differ, 95 percent against 99.7 percent, so the numbers are not directly interchangeable. Reynolds and colleagues give $6.3\times10^{-13}\,\mathrm{GeV^{-1}}$ for their field model A. Changing the field model alone moves the X-ray limit by tens of percent.

A checklist for the next sightline

1. Fix units for every input and record them in the file.
2. Save field angles explicitly, not only a seed.
3. Check normalization on every row.
4. Check a uniform case against the closed formula and a resonance case for flatness.
5. Compare a mean over many realizations with the ensemble rule, and report the pull.
6. Sample energies finely enough to resolve one realization's ripples.
7. Label every curve as a simulation, an ideal reference or an observation, and keep published limits with their confidence level, mass range and field model.

Try it in the Lab: Export a sightline from the Lab and compare it with the saved outputs.