---
title: "Recompute a cluster sightline from saved inputs"
description: "Rerun uniform and random-field sightlines with a standard-library Python script, check normalization and a Monte Carlo mean, and inspect four failure cases."
canonical_url: "https://sandboxphysics.com/stories/axion-photon-mixing/recompute-a-cluster-sightline/"
language: "en"
page_type: "story"
alternate_languages:
  zh-CN: "https://sandboxphysics.com/zh/stories/axion-photon-mixing/recompute-a-cluster-sightline/"
last_modified: "2026-11-12"
---

[Searching for axions with cosmic magnetic fields](https://sandboxphysics.com/stories/axion-photon-mixing/)

# Recompute a cluster sightline from saved inputs

Sandbox Physics · 2026-11-12 · L3 · Method

By the end of this guide you will have rerun three saved photon–axion sightlines with a short Python script that shares no code with the Lab. You will check that probability is conserved, that the uniform case matches the closed formula, and that a Monte Carlo mean over random field directions matches the ensemble rule. You will also see four ways the calculation goes wrong. Everything here is an offline simulation of model inputs. None of it is an observation.

Download three files: the [saved inputs](https://sandboxphysics.com/stories/axion-photon-mixing/materials/sightline-inputs.json), the [script](https://sandboxphysics.com/stories/axion-photon-mixing/materials/recompute_sightline.py) and the [saved outputs](https://sandboxphysics.com/stories/axion-photon-mixing/materials/sightline-outputs.json). 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](https://sandboxphysics.com/stories/axion-photon-mixing/photon-and-axion-as-two-states/) and the random-cell model of [Random domains, polarization and stellar cooling](https://sandboxphysics.com/stories/axion-photon-mixing/random-domains-polarization-and-stars/).

## 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

$$
\Pi_L=\frac{\sqrt{(\rho_{xx}-\rho_{yy})^2+4(\mathrm{Re}\,\rho_{xy})^2}}{\rho_{xx}+\rho_{yy}}.
$$

## 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:

**Table 5.1** Monte Carlo mean against the random-angle rule, at 3 keV.

| Quantity | Value |
| --- | --- |
| Mean axion fraction from 2000 realizations | $0.2434$ |
| Standard error of that mean | $0.0027$ |
| Product rule over 24 cells | $0.2399$ |
| Difference in standard errors | $1.3$ |

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.

**Table 5.2** Saved failure cases and what they show.

| Mistake | Measured effect | Lesson |
| --- | --- | --- |
| Plasma term with the wrong sign at resonance | Survival spread over energy rises from $2\times10^{-15}$ to $0.013$ | The resonance needs $\Delta_{\rm pl}$ and $\Delta_a$ to have the same sign |
| Twelve energies instead of 400 | Interpolation misses by up to $0.28$ | One realization has narrow ripples |
| Sixteen forward-Euler steps in one strong cell | Total probability $1.055$ | Euler steps are not unitary |
| Treating one sightline as the mean | Largest gap $0.22$ over 0.5 to 10 keV | Compare realizations with realizations |

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.

**Table 5.3** Published limits shown in the Lab and the assumptions they rest on.

| Limit | Value | Holds for | Rests on |
| --- | --- | --- | --- |
| Chandra spectrum of NGC 1275 | $7.9\times10^{-13}\,\mathrm{GeV^{-1}}$, 99.7% | $m_a<10^{-12}$ eV | Field model B of the Perseus cluster |
| Fermi-LAT spectrum of NGC 1275 | $5\times10^{-12}\,\mathrm{GeV^{-1}}$, 95% | 0.5 to 5 neV | Six years of LAT data and an assumed Perseus field model |
| CAST helioscope | $5.8\times10^{-11}\,\mathrm{GeV^{-1}}$, 95% | $m_a<0.02$ eV | The solar axion flux from the standard solar model |
| Horizontal-branch star counts | $6.6\times10^{-11}\,\mathrm{GeV^{-1}}$, 95% | up to a few keV | Stellar models and the helium abundance |

The X-ray value is eq. (5) of [Reynolds et al., 2020](https://arxiv.org/abs/1907.05475v4). The gamma-ray value is from the abstract of [Ajello et al., 2016](https://arxiv.org/abs/1603.06978v1), 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](https://pdg.lbl.gov/2025/reviews/rpp2025-rev-axions.pdf). The stellar mass edge in the Lab is set conservatively at 1 keV, because [Ayala et al., 2014](https://arxiv.org/abs/1406.6053v2) 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](https://sandboxphysics.com/labs/axion-photon-conversion/): Export a sightline from the Lab and compare it with the saved outputs.

## Sources

1. [Reynolds et al., 2020](https://arxiv.org/abs/1907.05475v4)
2. [Ajello et al., 2016](https://arxiv.org/abs/1603.06978v1)
3. [PDG review of axions, 2025](https://pdg.lbl.gov/2025/reviews/rpp2025-rev-axions.pdf)
4. [Ayala et al., 2014](https://arxiv.org/abs/1406.6053v2)

## Keep reading

Searching for axions with cosmic magnetic fields · Reading path

### From two mixing states to a published limit

1. [Article **How could an invisible particle change a distant spectrum?**](https://sandboxphysics.com/stories/axion-photon-mixing/invisible-particles-and-distant-spectra/)
2. [Article **Photon and axion as two mixing states**](https://sandboxphysics.com/stories/axion-photon-mixing/photon-and-axion-as-two-states/)
3. [Article **Mixing in a cluster magnetic field and plasma**](https://sandboxphysics.com/stories/axion-photon-mixing/mixing-in-cluster-fields-and-plasma/)
4. [Article **Random domains, polarization and stellar cooling**](https://sandboxphysics.com/stories/axion-photon-mixing/random-domains-polarization-and-stars/)
5. [You are here **Recompute a cluster sightline from saved inputs**](https://sandboxphysics.com/stories/axion-photon-mixing/recompute-a-cluster-sightline/)
6. [Article **Reading the Chandra NGC 1275 axion limit**](https://sandboxphysics.com/stories/axion-photon-mixing/reading-the-chandra-ngc-1275-limit/)