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 3, 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 ) \]

Simulated isotropic galaxy map (+ dipole)

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

Linear fit to density vs. declination (Secrest+21)

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!
Take data-generating process \(f_M: \mathbf{\theta} \to \mathbf{x}\).
Use a neural network to learn \(\alertmath{\mathcal{L}}\).
This is SBI.

Normalising Flows

Goal: describe the \(d\)-dimensional base distribution \({\mathcal{L}(\mathbf{x}\mid\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).

Results

Bayes Factors

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