Numeric Model Comparison: Chang-Cooper vs Analytic¶

This example compares the analytic synchrotron model (model="sync") with the full numeric Chang-Cooper electron distribution solver (model="numeric") to validate consistency and highlight where the numeric model differs.

Background¶

The default "sync" model uses the standard analytic broken-power-law synchrotron spectrum (Sari, Piran & Narayan 1998). The "numeric" model instead constructs the steady-state electron energy distribution on a log-spaced grid and computes synchrotron emissivity and self-absorption by integrating the exact spectral kernel \(F(\nu/\nu_c)\) over the distribution. This naturally captures:

• Exact spectral shapes from the Dermer (2009) synchrotron kernel, rather than broken power-law approximations
• Synchrotron self-absorption computed directly from the electron distribution derivative \(d/d\gamma[N/\gamma^2]\)
• Smooth spectral transitions at the characteristic breaks \(\nu_m\) and \(\nu_c\)
• Pair production at very high energies (optional, via include_pp=1)

The analytic and numeric models agree to \(\sim 10\%\) at optically thin frequencies, with differences arising from the broken power-law approximation in the analytic model. At radio frequencies where SSA matters, the two models use different absorption treatments and can differ by a factor of \(\sim 2\).

Physical parameters¶

We use a standard on-axis top-hat jet:

Computing the models¶

import numpy as np
from blastwave import FluxDensity_tophat

DAY = 86400.0

P = {
    "Eiso":    1e53,
    "lf":      300.0,
    "theta_c": 0.1,
    "A":       0.0,
    "n0":      1.0,
    "eps_e":   0.1,
    "eps_b":   0.01,
    "p":       2.2,
    "theta_v": 0.0,
    "d":       1000.0,
    "z":       0.2,
}

t = np.geomspace(0.01 * DAY, 1000.0 * DAY, 200)

# Analytic synchrotron (Sari+98 broken power law)
F_sync = {
    "100 GHz": FluxDensity_tophat(t, 1e11 * np.ones_like(t), P,
                                   tmin=1.0, tmax=1e8, model="sync"),
    "1 keV":   FluxDensity_tophat(t, 2.418e17 * np.ones_like(t), P,
                                   tmin=1.0, tmax=1e8, model="sync"),
}

# Chang-Cooper numeric
F_num = {
    "100 GHz": FluxDensity_tophat(t, 1e11 * np.ones_like(t), P,
                                   tmin=1.0, tmax=1e8, model="numeric"),
    "1 keV":   FluxDensity_tophat(t, 2.418e17 * np.ones_like(t), P,
                                   tmin=1.0, tmax=1e8, model="numeric"),
}

Light curves and residuals¶

The left panel shows light curves at three representative frequencies. Solid lines are the analytic sync model; dashed lines are the numeric model. At optically thin frequencies (100 GHz near the spectral peak, optical, X-ray), the two models track each other closely across three decades in time.

The right panel shows fractional residuals. The numeric model typically agrees with the analytic Sari+98 formula to within \(\sim 10\text{--}30\%\), with the sign and magnitude of the offset depending on the spectral regime:

• Near the spectral peak (\(\nu \sim \nu_m\)): the numeric model can be \(\sim 20\text{--}30\%\) higher because the exact \(F(x)\) kernel is broader than the broken power-law approximation
• Above the cooling break (\(\nu \gg \nu_c\)): the numeric model is \(\sim 10\%\) lower, reflecting the difference between the exact spectral integral and the Sari+98 coefficient

Broadband SED comparison¶

The spectral energy distribution (SED) highlights differences near the characteristic breaks \(\nu_m\), \(\nu_c\), and the SSA turnover.

Key features visible in the SEDs:

• At low frequencies, the numeric model shows SSA self-absorption (the \(\nu^{5/2}\) rise) while the bare sync model does not include SSA. The SSA turnover is computed directly from the electron distribution.
• Near the spectral breaks (\(\nu_m\), \(\nu_c\)), the numeric model produces smooth transitions rather than the sharp power-law junctions of the analytic model.
• At high frequencies, both models follow the same \(\nu^{-p/2}\) power law with \(\sim 10\%\) normalization difference.

Discussion¶

The "numeric" model constructs the steady-state cooled electron distribution analytically on a log-spaced \(\gamma\) grid:

• Slow cooling (\(\gamma_c > \gamma_m\)): \(N(\gamma) \propto \gamma^{-p}\) for \(\gamma_m \le \gamma \le \gamma_c\), steepening to \(\gamma^{-(p+1)}\) above \(\gamma_c\)
• Fast cooling (\(\gamma_c < \gamma_m\)): \(N(\gamma) \propto \gamma^{-2}\) for \(\gamma_c \le \gamma \le \gamma_m\), steepening to \(\gamma^{-(p+1)}\) above \(\gamma_m\)

The emissivity is then computed by integrating the exact Dermer (2009) synchrotron kernel \(F(\nu/\nu_c(\gamma))\) over this distribution, and SSA is computed from the distribution derivative. This avoids the broken power-law approximations of the analytic model while remaining fast (a single grid evaluation per frequency, no PDE time-stepping).

The \(\sim 10\text{--}30\%\) differences between the models are expected and well-understood:

1. The Sari+98 broken power-law approximation uses asymptotic spectral slopes matched at sharp break frequencies. The exact integral smoothly interpolates between regimes.
2. The pitch-angle averaging factor (PITCH_ANGLE_AVG \(\approx 0.72\)) used in the analytic model is an approximate correction; the numeric model computes the perpendicular emissivity, leading to small normalization offsets.
3. SSA in the numeric model is computed self-consistently from the electron distribution, while the analytic SSA models use a Rayleigh-Jeans blackbody approximation.

Full script¶

The complete analysis script is at examples/numeric_comparison.py.