Wising up to CatWISE

Using simulation-based inference to measure the cosmic dipole

Oliver Oayda

[email protected]

PhD Candidate

Sydney Institute for Astronomy

The University of Sydney

September 15, 2026

Simulation-based inference resolves systematics.

The Kinematic Dipole

Kinematic dipole as a sphere (Earth inside)

Counting Galaxies

\[ \require{color} \definecolor{secondarycolor}{RGB}{235,129,27} \newcommand{\alertmath}[1]{{\color{secondarycolor}{#1}}} N_{i} = \overline{N} (1 + \alertmath{\mathcal{D}} \cos \theta_i ) \]

Special Relativity

Aberration

Doppler boosting

// Fixed simulation inputs. Keeping these together makes the assumed physics
// and the plotted population straightforward to audit.
doppler_config = ({
  nSources: 60000,
  randomSeed: 73,
  differentialCountIndex: 2.5, // dN/dS_nu proportional to S_nu^(-2.5)
  minFlux: 1,
  maxFlux: 100,
  plotMaxFlux: 200,
  nBins: 52,
  spectralIndex: 1,
  cosTheta: 1,                 // small patch in the forward hemisphere
  fluxCut: 5
})

// Inverse-CDF samples from the truncated differential power law
// p(S_nu) proportional to S_nu^(-p), for S_min <= S_nu <= S_max.
doppler_rest_fluxes = {
  const p = doppler_config.differentialCountIndex;
  const exponent = 1 - p;
  const low = Math.pow(doppler_config.minFlux, exponent);
  const high = Math.pow(doppler_config.maxFlux, exponent);
  const rng = d3.randomLcg(doppler_config.randomSeed);

  return Array.from({ length: doppler_config.nSources }, () =>
    Math.pow(low + rng() * (high - low), 1 / exponent)
  );
}

// Observer-frame transformation requested for the forward patch:
//   gamma = 1 / sqrt(1 - beta^2)
//   delta = gamma (1 + beta cos(theta))
//   S'_nu = S_nu delta^(1 + alpha)
doppler_transform = {
  const beta = doppler_beta;
  const gamma = 1 / Math.sqrt(1 - beta * beta);
  const delta = gamma * (1 + beta * doppler_config.cosTheta);
  const fluxMultiplier = Math.pow(
    delta,
    1 + doppler_config.spectralIndex
  );

  return { beta, gamma, delta, fluxMultiplier };
}

doppler_bins = {
  const logMin = Math.log10(doppler_config.minFlux);
  const logMax = Math.log10(doppler_config.plotMaxFlux);
  const step = (logMax - logMin) / doppler_config.nBins;

  return Array.from({ length: doppler_config.nBins + 1 }, (_, i) =>
    Math.pow(10, logMin + i * step)
  );
}

doppler_histograms = {
  function histogram(values, series) {
    const counts = Array(doppler_config.nBins).fill(0);
    for (const value of values) {
      if (
        value < doppler_config.minFlux ||
        value >= doppler_config.plotMaxFlux
      ) continue;

      const index = Math.min(
        doppler_config.nBins - 1,
        d3.bisectRight(doppler_bins, value) - 1
      );
      counts[index] += 1;
    }

    const bars = counts
      .map((count, i) => ({
        x0: doppler_bins[i],
        x1: doppler_bins[i + 1],
        count,
        series
      }))
      .filter((d) => d.count > 0);

    // Close the stepped outline at both limits. On the logarithmic y-axis,
    // count = 1 is the visible baseline used for an empty boundary bin.
    const first = counts.findIndex((count) => count > 0);
    const last = counts.findLastIndex((count) => count > 0);
    const outline = first < 0
      ? []
      : [
          { x: doppler_bins[first], count: 1, series },
          ...counts.slice(first, last + 1).map((count, offset) => ({
            x: doppler_bins[first + offset],
            count: Math.max(1, count),
            series
          })),
          { x: doppler_bins[last + 1], count: 1, series }
        ];

    return { bars, outline };
  }

  const observedFluxes = doppler_rest_fluxes.map(
    (flux) => flux * doppler_transform.fluxMultiplier
  );

  const rest = histogram(doppler_rest_fluxes, "Rest frame");
  const observed = histogram(observedFluxes, "Observer frame");

  return {
    bars: [...rest.bars, ...observed.bars],
    outlines: [...rest.outline, ...observed.outline]
  };
}

Plot.plot({
  width: 500,
  height: 390,
  marginLeft: 62,
  marginRight: 18,
  marginTop: 34,
  marginBottom: 54,
  style: {
    fontFamily: "Fira Sans, Helvetica, sans-serif",
    fontSize: "13px"
  },
  color: {
    domain: ["Rest frame", "Observer frame"],
    range: ["#50a96d", "#ef7c3b"],
    legend: true
  },
  x: {
    type: "log",
    domain: [doppler_config.minFlux, doppler_config.plotMaxFlux],
    ticks: [1, 10, 100],
    label: "Flux density"
  },
  y: {
    type: "log",
    domain: [1, 12000],
    ticks: [1, 10, 100, 1000, 10000],
    label: "Count in each flux bin"
  },
  marks: [
    Plot.rectY(doppler_histograms.bars, {
      x1: "x0",
      x2: "x1",
      y1: 1,
      y2: "count",
      fill: "series",
      fillOpacity: 0.22,
      inset: 0
    }),
    Plot.line(doppler_histograms.outlines, {
      x: "x",
      y: "count",
      z: "series",
      stroke: "series",
      curve: "step-after",
      strokeWidth: 2.8
    }),
    Plot.ruleX([doppler_config.fluxCut], {
      stroke: "#222",
      strokeDasharray: "5,4",
      strokeOpacity: 0.65
    }),
    Plot.ruleY([1], {
      stroke: "#000",
      strokeWidth: 2.8
    }),
    Plot.ruleX([doppler_config.minFlux], {
      stroke: "#000",
      strokeWidth: 2.8
    }),
    Plot.text(
      [{ x: doppler_config.fluxCut, label: "Flux cut" }],
      {
        x: "x",
        text: "label",
        frameAnchor: "top",
        dy: 8,
        dx: 5,
        textAnchor: "start",
        fill: "#222"
      }
    )
  ]
})

viewof doppler_beta = {
  const root = html`<label class="doppler-beta-control">
    <span>Speed (c)</span>
    <input class="doppler-beta-slider" type="range" min="0" max="0.3" step="0.005" value="0">
    <output></output>
  </label>`;
  const input = root.querySelector("input");
  const output = root.querySelector("output");

  root.value = input.valueAsNumber;
  output.textContent = root.value.toFixed(3);
  input.addEventListener("input", () => {
    root.value = input.valueAsNumber;
    output.textContent = root.value.toFixed(3);
    root.dispatchEvent(new Event("input", { bubbles: true }));
  });

  return root;
}

The Amplitude Excess

The cosmic dipole should be consistent with the CMB dipole.

The CatWISE Sample

Correcting for the Bias

Null hypothesis test for the dipole amplitude

Null hypothesis test for the dipole direction

Is it good enough to make an ad hoc correction after the fact?

WISE’s Scanning Law

Obeys a scanning law over the survey’s lifetime

Photometric Errors

WISE coverage in W1 band for the Secrest+21 sample

Source photometric error (%) in W1 band for the Secrest+21 sample

Eddington Bias

\[ \require{color} \definecolor{tomato}{RGB}{255,88,62} \definecolor{cornflowerblue}{RGB}{89,138,234} {\color{tomato}{\textsf{Sampled flux}}} = {\color{cornflowerblue}{\textsf{True flux}}} + \textsf{Noise} \]

eddington_bins = {
  const xMin = 1;
  const xMax = 30;
  const nBins = 80;
  const logMin = Math.log10(xMin);
  const logMax = Math.log10(xMax);
  const step = (logMax - logMin) / nBins;

  return Array.from({ length: nBins + 1 }, (_, i) =>
    Math.pow(10, logMin + i * step)
  );
}

viewof eddington_cut_index = {
  const minIndex = d3.bisectLeft(eddington_bins, 1.2);
  const maxIndex = d3.bisectRight(eddington_bins, 4) - 1;
  const value = d3.bisectLeft(eddington_bins, 2);
  const root = html`<label style="display: flex; align-items: center; gap: 0.55em;">
    <span>Flux cut</span>
    <input type="range" min=${minIndex} max=${maxIndex} step="1" value=${value}>
    <span style="font-variant-numeric: tabular-nums; min-width: 2.8em;"></span>
  </label>`;
  const input = root.querySelector("input");
  const output = root.querySelector("span:last-child");

  root.value = input.valueAsNumber;
  output.textContent = eddington_bins[root.value].toFixed(2);
  input.addEventListener("input", () => {
    root.value = input.valueAsNumber;
    output.textContent = eddington_bins[root.value].toFixed(2);
    root.dispatchEvent(new Event("input", { bubbles: true }));
  });

  return root;
}

eddington_cut = eddington_bins[eddington_cut_index]

viewof eddington_noise = {
  const root = html`<label style="display: flex; align-items: center; gap: 0.55em;">
    <span>Noise</span>
    <input class="eddington-noise-slider" type="range" min="0" max="1.5" step="0.025" value="0.0">
    <span style="font-variant-numeric: tabular-nums; min-width: 2.8em;"></span>
  </label>`;
  const input = root.querySelector("input");
  const output = root.querySelector("span:last-child");

  root.value = input.valueAsNumber;
  output.textContent = root.value.toFixed(2);
  input.addEventListener("input", () => {
    root.value = input.valueAsNumber;
    output.textContent = root.value.toFixed(2);
    root.dispatchEvent(new Event("input", { bubbles: true }));
  });

  return root;
}

eddington_samples = {
  const n = 60000;
  const x = 3;
  const sLow = 1;
  const rng = d3.randomLcg(42);

  return Array.from({ length: n }, () => {
    const trueFlux = sLow * Math.pow(1 - rng(), -1 / x);
    return { trueFlux };
  });
}

eddington_histogram = {
  const xMin = 1;
  const xMax = 30;
  const bins = eddington_bins;
  const nBins = bins.length - 1;

  function histogram(values, series, color, opacity) {
    const counts = Array(nBins).fill(0);
    for (const value of values) {
      if (value < xMin || value >= xMax) continue;
      const i = Math.min(
        nBins - 1,
        d3.bisectRight(bins, value) - 1
      );
      counts[i] += 1;
    }
    return counts.map((count, i) => ({
      x0: bins[i],
      x1: bins[i + 1],
      count,
      series,
      color,
      opacity
    }));
  }

  const fluxError =
    eddington_noise === 0
      ? () => 0
      : d3.randomNormal.source(d3.randomLcg(314))(0, eddington_noise);
  const measuredSources = eddington_samples.map((d) => {
    const sampledFlux = d.trueFlux + fluxError();
    return {
      trueFlux: d.trueFlux,
      sampledFlux,
      included: sampledFlux > eddington_cut
    };
  });
  const includedSources = measuredSources.filter((d) => d.included);
  const trueFlux = eddington_samples.map((d) => d.trueFlux);
  const sampledFlux = includedSources.map((d) => d.sampledFlux);

  return [
    ...histogram(trueFlux, "True distribution", "#598aea", 0.3),
    ...histogram(sampledFlux, "Sampled above cut", "#ff583e", 0.4)
  ];
}

Plot.plot({
  width: 720,
  height: 410,
  marginLeft: 72,
  marginRight: 24,
  marginTop: 34,
  marginBottom: 58,
  style: {
    fontFamily: "Fira Sans, Helvetica, sans-serif",
    fontSize: "15px"
  },
  x: {
    type: "log",
    domain: [1, 30],
    ticks: [1, eddington_cut, 10, 20],
    tickFormat: (d) =>
      d === eddington_cut ? d.toFixed(2) : d3.format("~g")(d),
    label: "Flux density (brightness)"
  },
  y: {
    grid: false,
    label: "Count in each flux bin"
  },
  marks: [
    Plot.rectY(eddington_histogram, {
      x1: "x0",
      x2: "x1",
      y1: 0,
      y2: "count",
      fill: "color",
      fillOpacity: "opacity",
      inset: 0
    }),
    Plot.line(eddington_histogram, {
      x: "x0",
      y: "count",
      z: "series",
      curve: "step-after",
      stroke: "color",
      strokeWidth: 3,
      strokeOpacity: 0.95
    }),
    Plot.ruleX([eddington_cut], {
      stroke: "#222",
      strokeDasharray: "5,4",
      strokeOpacity: 0.55
    }),
    Plot.ruleY([0], {
      stroke: "#000",
      strokeWidth: 2.8
    }),
    Plot.ruleX([1], {
      stroke: "#000",
      strokeWidth: 2.8
    })
  ]
})

What if the noise varies over the sky? 🤔

This is just something we can simulate!

Simulations

\[m'_{\nu } \to m'_{\nu} + \alertmath{\Delta m'_{\nu}} \]

One Thing Left To Do…

\(\text{Posterior} = \dfrac{\text{Prior} \times \alertmath{\text{Likelihood}} }{\text{Evidence}}\)

Don't know \(\alertmath{\mathcal{L}}\)? No problem!
\(f_M\): parameters \(\to\) data.
Use a neural network to learn \(\alertmath{\mathcal{L}}\).
This is SBI.

Normalising Flows

Goal: describe the distribution \({\mathcal{L}(\boldsymbol{x}\mid\boldsymbol{\theta})}\) for data \(\boldsymbol{x}\) and parameters \(\boldsymbol{\theta}\)
We use a series of Masked Autoregressive Flows (MAFs).
A flow transforms a simple distribution into a complex one.

Composing individual bijections, adapted from Papamakarios et al. (2019).

The Nuts & Bolts

Algorithm: Sequential Neural Likelihood Estimator (SNLE)

Input: observation \(\boldsymbol{x}_0\), prior \( \pi\), simulator \(f_M\)

Output: posterior \( \mathcal{P}\), evidence \( \mathcal{Z}\)

\[ \definecolor{posteriorblue}{RGB}{0,181,235} \definecolor{evidencepurple}{RGB}{128,0,255} \definecolor{priorgreen}{RGB}{64,128,90} \definecolor{alertcolor}{RGB}{235,129,27} {\color{posteriorblue}\mathcal{P}(\boldsymbol{\theta}\mid\boldsymbol{x}_0)} = \frac{ {\color{priorgreen}\pi(\boldsymbol{\theta})} \times {\color{alertcolor}\mathcal{L}(\boldsymbol{x}_0\mid\boldsymbol{\theta})} }{ {\color{evidencepurple}\mathcal{Z}(\boldsymbol{x}_0)} } = \frac{ \mathcal{J}(\boldsymbol{x}_0, \boldsymbol{\theta}) }{ {\color{evidencepurple}\mathcal{Z}} } \]

Let proposal \(=\) prior

for each training round:

for each simulation this round:

sample \(\boldsymbol{\theta} \hookleftarrow\) proposal

simulate \(\boldsymbol{x} = f_M (\boldsymbol{\theta})\)

append the \( (\boldsymbol{\theta}, \boldsymbol{x}) \) pair

learn approx. likelihood \( \mathcal{L}_{\phi} \) using all pairs

compute posterior \( \mathcal{P} \) and evidence \(\mathcal{Z} \)

set proposal \(=\) posterior

return final posterior & evidence estimate

\[ \definecolor{priorgreen}{RGB}{64,128,90} \definecolor{alertcolor}{RGB}{235,129,27} \begin{aligned} \mathrm{Loss}(\phi) &= -\frac{1}{N} \sum_{i=1}^N \ln {\color{alertcolor}\mathcal{L}_{\phi}(x_i \mid \theta_i)} \\ &\approx -\mathbb{E}_{\mathcal{J}} \left[\ln {\color{alertcolor}\mathcal{L}_{\phi}}\right] \\ &= \mathbb{E}_{{\color{priorgreen}\pi}} \left[ D_{\mathrm{KL}} \left( {\color{alertcolor}\mathcal{L}} \,\|\, {\color{alertcolor}\mathcal{L}_{\phi}} \right) \right] + \mathrm{const.} \end{aligned} \]

Results

CatSIM Results

Confirms cosmic dipole tension.

Bayes Factors

Model	\(\ln \mathcal{B}\)
Dipole from Secrest+21	\(2.8\)	\(\pm\)	\(0.9\)
Dipole from Dam+23	\(1.9\)	\(\pm\)	\(0.8\)
Dipole with free parameters	\(0.0\)	\(\pm\)	\(0.0\)
CMB direction, free velocity	\(-2.4\)	\(\pm\)	\(0.9\)
CMB velocity & direction	\(-8.3\)	\(\pm\)	\(0.9\)
Free dipole but no extra error	\(-113.6\)	\(\pm\)	\(6.4\)