Part III: Solar Magnetic Activity | Chapter 12

Solar Flares

CSHKP model, magnetic reconnection, Sweet-Parker and Petschek rates, and particle acceleration

12.1 The Standard Flare Model (CSHKP)

Derivation 1: Flare Energy from Magnetic Field Annihilation

The CSHKP (Carmichael-Sturrock-Hirayama-Kopp-Pneuman) model describes the standard two-ribbon flare geometry: a vertical current sheet forms above a magnetic arcade, reconnection occurs, and energy is released into the corona and chromosphere.

Step 1. The available magnetic free energy is the excess over the potential field:

$$E_{\text{free}} = \frac{1}{2\mu_0}\int (B^2 - B_p^2) \, dV$$

Step 2. For a volume $V = L^3$ with field $B \sim 0.03$ T (300 G):

$$\boxed{E_{\text{free}} \sim \frac{B^2 L^3}{2\mu_0} \sim \frac{(0.03)^2 \times (3 \times 10^7)^3}{2 \times 4\pi \times 10^{-7}} \sim 10^{25} \text{ J (X-class flare)}}$$

This is $\sim 10^{32}$ erg, consistent with observations of large flares. The flare classification is based on peak soft X-ray flux: C ($10^{-6}$ W/m$^2$), M ($10^{-5}$), X ($10^{-4}$).

12.2 Sweet-Parker Reconnection

Derivation 2: Sweet-Parker Reconnection Rate

Sweet (1958) and Parker (1957) derived the reconnection rate for a steady-state resistive current sheet of length $L$ and thickness $\delta$.

Step 1. Mass conservation: plasma flows in at speed $v_{\text{in}}$through area $L$ and out at $v_{\text{out}}$ through area $\delta$:

$$\rho v_{\text{in}} L = \rho v_{\text{out}} \delta$$

Step 2. The outflow is driven by the magnetic tension, accelerating plasma to the Alfven speed: $v_{\text{out}} = v_A = B/\sqrt{\mu_0\rho}$.

Step 3. The current sheet thickness is set by resistive diffusion balancing the inflow:$\delta = \eta / v_{\text{in}}$.

Step 4. Combining these:

$$\boxed{\frac{v_{\text{in}}}{v_A} = \frac{1}{\sqrt{R_m}} = \frac{1}{\sqrt{Lv_A/\eta}}}$$$$\boxed{\delta = \frac{L}{\sqrt{R_m}}}$$

For the solar corona: $R_m \sim 10^{12}$, so $v_{\text{in}}/v_A \sim 10^{-6}$. This gives a reconnection time $\tau \sim L/(v_{\text{in}}) \sim 10^6 L/v_A \sim 10^8$ s — far too slow to explain flares (which release energy in minutes). This is the "reconnection rate problem."

12.3 Petschek Reconnection

Derivation 3: Fast Reconnection via Standing Shocks

Petschek (1964) proposed that the current sheet need not extend over the full system length. Instead, slow-mode standing shocks form, processing the inflowing plasma over a much larger area.

Step 1. The Petschek geometry has a compact diffusion region of length $l \ll L$with two pairs of standing slow-mode shocks extending from the ends.

Step 2. The maximum reconnection rate:

$$\boxed{\frac{v_{\text{in}}}{v_A} \sim \frac{\pi}{8\ln R_m} \sim 0.01\text{--}0.1}$$

This is fast enough to explain flare timescales. Modern numerical simulations with anomalous resistivity, plasmoid instability, and Hall MHD confirm reconnection rates of 0.01-0.1 of the Alfven speed, consistent with Petschek-like fast reconnection.

12.4 Particle Acceleration

Derivation 4: DC Electric Field Acceleration

Step 1. In the reconnection region, a DC electric field $E \sim v_{\text{in}}B$accelerates charged particles along the current sheet:

$$E \sim v_{\text{in}} B \sim 0.01 v_A B = 0.01 \times 10^6 \times 0.01 = 100 \text{ V/m}$$

Step 2. Over a length $\ell$, a particle gains energy:

$$\boxed{\varepsilon = qE\ell \sim 100 \times 10^7 \text{ eV} = 1 \text{ GeV (for } \ell = 10^7 \text{ m)}}$$

Flares accelerate electrons to 10-100 keV (producing hard X-rays via bremsstrahlung) and protons to MeV-GeV energies. The nonthermal emission carries 10-50% of the total flare energy.

12.5 Thermal and Nonthermal Emission

Derivation 5: The Neupert Effect

The Neupert effect (1968) is the empirical observation that the soft X-ray (SXR) light curve is approximately the time integral of the hard X-ray (HXR) light curve.

Step 1. Accelerated electrons stream down magnetic field lines and produce HXR bremsstrahlung when they hit the dense chromosphere (thick-target model).

Step 2. The deposited energy heats chromospheric plasma, which expands upward (chromospheric evaporation), filling coronal loops with hot ($\sim 10^7$ K) plasma that emits in soft X-rays.

Step 3. The SXR emission measure grows as the integral of the energy deposition rate:

$$\boxed{F_{\text{SXR}}(t) \propto \int_0^t F_{\text{HXR}}(t') \, dt'}$$

This is a key diagnostic linking the acceleration (impulsive phase, HXR) to the thermal response (gradual phase, SXR), supporting the thick-target model of flare energy deposition.

Numerical Simulation

Solar Flares: Reconnection Rates, Neupert Effect, Current Sheet, Energy Spectrum

Python

script.py150 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle('Solar Flare Physics', fontsize=16, fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

for ax in axes.flat:
    ax.set_facecolor('#0a0a0a')
    ax.tick_params(colors='white')
    ax.xaxis.label.set_color('white')
    ax.yaxis.label.set_color('white')
    ax.title.set_color('white')
    for spine in ax.spines.values():
        spine.set_color('#333')

# --- Panel 1: Reconnection rates ---
ax1 = axes[0, 0]
Rm = np.logspace(4, 14, 500)

SP_rate = 1.0 / np.sqrt(Rm)
Petschek_rate = np.pi / (8 * np.log(Rm))
fast_rate = 0.01 * np.ones_like(Rm)

ax1.loglog(Rm, SP_rate, color='#60a5fa', linewidth=2.5, label='Sweet-Parker')
ax1.loglog(Rm, Petschek_rate, color='#f43f5e', linewidth=2.5, label='Petschek')
ax1.loglog(Rm, fast_rate, color='#fbbf24', linewidth=2, linestyle='--', label='Observed rate ~0.01')
ax1.axvline(x=1e12, color='white', linestyle=':', alpha=0.3, label='Solar corona R_m')
ax1.set_xlabel('Magnetic Reynolds number R_m')
ax1.set_ylabel('Reconnection rate v_in / v_A')
ax1.set_title('Reconnection Rate vs R_m')
ax1.set_ylim(1e-8, 1)
ax1.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: Neupert effect ---
ax2 = axes[0, 1]
t = np.linspace(0, 30, 500)  # minutes

# HXR burst (impulsive)
HXR = 0.8 * np.exp(-((t - 5)/1.5)**2) + 0.5 * np.exp(-((t - 8)/2)**2) + 0.3 * np.exp(-((t - 11)/1)**2)
HXR = np.maximum(HXR, 0)

# SXR = integral of HXR (Neupert effect) + cooling
from scipy.integrate import cumulative_trapezoid
SXR_raw = np.zeros_like(t)
SXR_raw[1:] = cumulative_trapezoid(HXR, t)
# Add radiative cooling
tau_cool = 15  # minutes
SXR = SXR_raw * np.exp(-t / tau_cool)
SXR = SXR / np.max(SXR)
HXR_norm = HXR / np.max(HXR)

ax2.plot(t, HXR_norm, color='#f43f5e', linewidth=2.5, label='HXR (nonthermal)')
ax2.plot(t, SXR, color='#60a5fa', linewidth=2.5, label='SXR (thermal)')
ax2.fill_between(t, 0, HXR_norm, alpha=0.15, color='#f43f5e')
ax2.set_xlabel('Time [minutes]')
ax2.set_ylabel('Normalized flux')
ax2.set_title('Neupert Effect: HXR and SXR Light Curves')
ax2.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: Current sheet geometry ---
ax3 = axes[1, 0]

# Draw Sweet-Parker current sheet
sheet_L = 3.0
sheet_delta = 0.2
ax3.fill_between([-sheet_L, sheet_L], [-sheet_delta, -sheet_delta],
                 [sheet_delta, sheet_delta], color='#fbbf24', alpha=0.3)
ax3.plot([-sheet_L, sheet_L], [sheet_delta, sheet_delta], color='#fbbf24', linewidth=2)
ax3.plot([-sheet_L, sheet_L], [-sheet_delta, -sheet_delta], color='#fbbf24', linewidth=2)

# Inflow arrows
for y_off in [-1.5, -1.0, 1.0, 1.5]:
    ax3.annotate('', xy=(0, sheet_delta * np.sign(y_off)),
                xytext=(0, y_off),
                arrowprops=dict(arrowstyle='->', color='#60a5fa', linewidth=2))

# Outflow arrows
for x_off in [-sheet_L, sheet_L]:
    ax3.annotate('', xy=(x_off * 1.3, 0),
                xytext=(x_off, 0),
                arrowprops=dict(arrowstyle='->', color='#f43f5e', linewidth=2))

# Magnetic field lines
for y in [0.4, 0.8, 1.2]:
    x_line = np.linspace(-4, -0.5, 50)
    ax3.plot(x_line, np.full_like(x_line, y), color='#4ade80', linewidth=1)
    ax3.plot(x_line, np.full_like(x_line, -y), color='#4ade80', linewidth=1)
    ax3.plot(-x_line, np.full_like(x_line, y), color='#a78bfa', linewidth=1)
    ax3.plot(-x_line, np.full_like(x_line, -y), color='#a78bfa', linewidth=1)

ax3.text(0, 1.8, 'v_in', fontsize=10, color='#60a5fa', ha='center')
ax3.text(4.0, 0.2, 'v_out ~ v_A', fontsize=10, color='#f43f5e', ha='center')
ax3.text(0, 0, 'delta', fontsize=9, color='#fbbf24', ha='center', va='center')
ax3.annotate('L', xy=(sheet_L, -0.4), fontsize=10, color='#fbbf24')
ax3.set_xlim(-5, 5)
ax3.set_ylim(-2, 2)
ax3.set_aspect('equal')
ax3.set_title('Sweet-Parker Current Sheet')
ax3.set_xlabel('x')
ax3.set_ylabel('y')

# --- Panel 4: Flare energy spectrum ---
ax4 = axes[1, 1]
E_flare = np.logspace(24, 33, 200)  # Joules

# Power-law frequency distribution
alpha_freq = 1.8
dN_dE = E_flare**(-alpha_freq)
dN_dE = dN_dE / dN_dE[0] * 1e4

ax4.loglog(E_flare, dN_dE, color='#f43f5e', linewidth=2.5)
ax4.axvline(x=1e24, color='#4ade80', linestyle='--', alpha=0.5)
ax4.axvline(x=1e28, color='#fbbf24', linestyle='--', alpha=0.5)
ax4.axvline(x=1e32, color='#fb923c', linestyle='--', alpha=0.5)

ax4.text(1e24, 1e-2, 'Nanoflares', fontsize=8, color='#4ade80', rotation=90, va='bottom')
ax4.text(1e28, 1e-2, 'Microflares', fontsize=8, color='#fbbf24', rotation=90, va='bottom')
ax4.text(1e32, 1e-2, 'X-class', fontsize=8, color='#fb923c', rotation=90, va='bottom')

ax4.set_xlabel('Flare energy [J]')
ax4.set_ylabel('dN/dE [relative]')
ax4.set_title('Flare Frequency-Energy Distribution')
ax4.annotate('dN/dE ~ E^-1.8', xy=(1e29, 1e-5), fontsize=10, color='white')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("Plot saved to output.png")

print("\n" + "=" * 60)
print("SOLAR FLARES")
print("=" * 60)
print("\n--- Reconnection Rates ---")
print("  Sweet-Parker: v_in/v_A ~ 1/sqrt(R_m) ~ 10^-6")
print("  Petschek: v_in/v_A ~ pi/(8*ln(R_m)) ~ 0.01-0.03")
print("  Observed: ~0.01-0.1 v_A")
print("\n--- Flare Energetics ---")
print("  Nanoflare: ~10^17 J (10^24 erg)")
print("  C-class: ~10^22 J (10^29 erg)")
print("  M-class: ~10^23-10^24 J")
print("  X-class: ~10^25 J (10^32 erg)")
print("\n--- Particle Acceleration ---")
print("  Electrons: 10-100 keV (HXR bremsstrahlung)")
print("  Protons: 1 MeV - 1 GeV")
print("  Total nonthermal energy: 10-50% of flare energy")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

12.6 Full Sweet-Parker Reconnection Derivation

From Mass and Energy Conservation

Step 1: Geometry. Consider a steady-state current sheet of length $L$ and half-thickness $\delta$. Plasma flows in at speed $v_{\text{in}}$ from top and bottom, and exits at speed $v_{\text{out}}$ from the two ends.

Step 2: Mass conservation. Incompressible flow in the 2D sheet:

$$\rho\,v_{\text{in}}\,(2L) = \rho\,v_{\text{out}}\,(2\delta) \quad\Rightarrow\quad v_{\text{in}} L = v_{\text{out}}\delta \quad\ldots(1)$$

Step 3: Energy conservation. The outflow is accelerated by the magnetic tension (the $j\times B$ force). The upstream magnetic energy density is converted to kinetic energy:

$$\frac{B^2}{2\mu_0} = \frac{1}{2}\rho\,v_{\text{out}}^2 \quad\Rightarrow\quad v_{\text{out}} = v_A = \frac{B}{\sqrt{\mu_0\rho}} \quad\ldots(2)$$

Step 4: Ohmic diffusion balance. Inside the sheet, the inflow speed must balance the rate of magnetic diffusion across the sheet thickness:

$$v_{\text{in}} = \frac{\eta}{\delta} \quad\Rightarrow\quad \delta = \frac{\eta}{v_{\text{in}}} \quad\ldots(3)$$

Step 5: Combining. From (1): $\delta = v_{\text{in}}L/v_A$. Equating with (3):

$$\frac{\eta}{v_{\text{in}}} = \frac{v_{\text{in}}L}{v_A} \quad\Rightarrow\quad v_{\text{in}}^2 = \frac{\eta\,v_A}{L}$$

Step 6. Defining $R_m = Lv_A/\eta$:

$$\boxed{\frac{v_{\text{in}}}{v_A} = \frac{1}{\sqrt{R_m}} \approx 10^{-6} \text{ for coronal } R_m\sim 10^{12}}$$$$\boxed{\delta = \frac{L}{\sqrt{R_m}} \approx 1 \text{ m for } L=10^7 \text{ m}}$$

The reconnection timescale $\tau_{SP} = L/v_{\text{in}} = \sqrt{\tau_A\tau_d}$ is the geometric mean of the Alfven time and diffusion time. For the corona: $\tau_{SP} \sim 10^8$ s (~3 years), far too slow for flares observed to release energy in minutes.

12.7 Why Petschek Reconnection Is Fast

Standing Slow-Mode Shocks and Logarithmic Scaling

Step 1. In the Petschek model, the resistive diffusion region has length $l^* \ll L$. Four standing slow-mode shocks extend from the corners of this compact region, subtending an angle$\theta$ to the inflow direction.

Step 2. Most of the energy conversion occurs at the shocks, not in the diffusion region. The opening angle $\theta$ determines the inflow speed: $v_{\text{in}} \approx v_A\sin\theta$.

Step 3. The maximum rate occurs when the diffusion region length $l^*$ shrinks to the minimum allowed by the Sweet-Parker scaling within that sub-region. Self-consistently:

$$l^* \sim L\exp\left(-\frac{\pi v_A}{4 v_{\text{in}}}\right)$$

Step 4. Maximizing the inflow rate subject to this constraint:

$$\boxed{\frac{v_{\text{in}}}{v_A}\bigg|_{\max} \approx \frac{\pi}{8\ln R_m} \approx 0.01\text{--}0.03 \text{ for coronal } R_m}$$

The logarithmic dependence on $R_m$ makes the rate only weakly sensitive to the resistivity, explaining why flare timescales are similar despite large uncertainties in the effective resistivity. Modern simulations with plasmoid instability achieve rates of ~0.01, consistent with observations.

12.8 Power-Law Spectrum from First-Order Fermi Acceleration

Deriving $f(E) \propto E^{-\delta}$ from Statistical Arguments

Step 1. At each interaction with converging magnetic mirrors (or equivalently, each complete Fermi cycle), a particle gains fractional energy $\epsilon = \Delta E/E$:

$$E_n = E_0(1+\epsilon)^n \quad\text{after } n \text{ interactions}$$

Step 2. The probability of remaining in the acceleration region after $n$interactions is $P_n = (1-P_{\text{esc}})^n$ where $P_{\text{esc}}$ is the escape probability per cycle.

Step 3. Eliminating $n$: from $n = \ln(E/E_0)/\ln(1+\epsilon)$:

$$N(>E) \propto P_n \propto \left(\frac{E}{E_0}\right)^{\ln(1-P_{\text{esc}})/\ln(1+\epsilon)}$$

Step 4. For small $\epsilon$ and $P_{\text{esc}}$:

$$\boxed{f(E) = \frac{dN}{dE} \propto E^{-\delta}, \quad \delta = 1 + \frac{P_{\text{esc}}}{\epsilon}}$$

For the first-order Fermi process at a shock with compression ratio $r$:$\epsilon \propto (r-1)/r$ and $P_{\text{esc}} \propto 1/r$, giving$\delta = (r+2)/(2(r-1))$ for the energy spectrum. For a strong shock ($r=4$):$\delta = 2$, producing the classic $E^{-2}$ spectrum. In flares, the observed hard X-ray spectral index of 3-5 corresponds to an electron spectrum of$\delta \approx 2\text{--}4$, broadly consistent.

12.9 CSHKP Flare Model Diagram

The standard CSHKP model shows the geometry of a two-ribbon flare with reconnection in a vertical current sheet above the magnetic arcade.

CSHKP standard flare model: magnetic reconnection at the X-point in the current sheet releases energy. Accelerated electrons stream down field lines producing HXR at chromospheric footpoints (ribbons). Heated plasma fills post-flare loops visible in SXR and EUV.

Extended Simulation: Reconnection Rates & GOES Light Curves

Extended: Reconnection Models, GOES Flare Simulation, Current Sheets, Fermi Spectra

Python

script.py141 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle('Extended: Reconnection Rates and Flare Diagnostics', fontsize=16, fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

# --- Panel 1: Detailed reconnection rate comparison ---
ax1 = axes[0, 0]
Rm = np.logspace(3, 14, 500)

SP = 1.0 / np.sqrt(Rm)
Petschek = np.pi / (8 * np.log(Rm))
# Plasmoid-mediated (for Rm > Rm_crit ~ 10^4)
plasmoid = np.where(Rm > 1e4, 0.01 * np.ones_like(Rm), SP)
# Collisionless (Hall MHD)
hall = 0.1 * np.ones_like(Rm)

ax1.loglog(Rm, SP, color='#60a5fa', linewidth=2.5, label='Sweet-Parker: 1/sqrt(Rm)')
ax1.loglog(Rm, Petschek, color='#f43f5e', linewidth=2.5, label='Petschek: pi/(8*ln(Rm))')
ax1.loglog(Rm, plasmoid, color='#fbbf24', linewidth=2, linestyle='--', label='Plasmoid instability: ~0.01')
ax1.loglog(Rm, hall, color='#4ade80', linewidth=2, linestyle=':', label='Hall MHD: ~0.1')

# Solar corona
ax1.axvspan(1e10, 1e14, alpha=0.05, color='white')
ax1.text(3e11, 2e-1, 'Solar corona', fontsize=9, color='white', ha='center')

ax1.set_xlabel('Magnetic Reynolds Number Rm')
ax1.set_ylabel('Reconnection rate v_in / v_A')
ax1.set_title('Reconnection Rate Models')
ax1.set_ylim(1e-8, 1)
ax1.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: GOES X-ray light curve simulation ---
ax2 = axes[0, 1]
t = np.linspace(-10, 120, 1000)  # minutes relative to flare onset

# Simulate a complex flare with multiple bursts
# SXR (1-8 angstrom)
sxr_base = 1e-7  # C1 background
bursts = [
    (0, 5, 1e-4),    # main burst
    (8, 3, 5e-5),    # secondary peak
    (15, 2, 2e-5),   # tertiary
]
sxr = np.full_like(t, sxr_base)
for t0, sigma, amp in bursts:
    rise = np.where(t > t0, amp * (1 - np.exp(-(t - t0) / sigma)), 0)
    decay_t = t0 + 3 * sigma
    decay = np.where(t > decay_t, amp * np.exp(-(t - decay_t) / 20), rise)
    combined = np.where(t > decay_t, decay, rise)
    sxr = sxr + combined

# HXR (25-50 keV)
hxr_base = 0.1
hxr = hxr_base + 100 * np.exp(-((t - 2)/2)**2) + 60 * np.exp(-((t - 9)/1.5)**2) + 30 * np.exp(-((t - 16)/1)**2)

ax2.semilogy(t, sxr, color='#f43f5e', linewidth=2.5, label='SXR 1-8 A (GOES)')
ax2_twin = ax2.twinx()
ax2_twin.plot(t, hxr, color='#fbbf24', linewidth=2, alpha=0.7, label='HXR 25-50 keV (RHESSI)')

# GOES classification lines
for level, name, col in [(1e-6, 'C', '#4ade80'), (1e-5, 'M', '#fbbf24'), (1e-4, 'X', '#f43f5e')]:
    ax2.axhline(y=level, color=col, linestyle=':', alpha=0.3)
    ax2.text(115, level*1.2, name, fontsize=9, color=col)

ax2.set_xlabel('Time [minutes from onset]')
ax2.set_ylabel('SXR Flux [W/m^2]', color='#f43f5e')
ax2_twin.set_ylabel('HXR counts [relative]', color='#fbbf24')
ax2.set_title('Simulated GOES X-ray Light Curve')
ax2.set_xlim(-10, 120)
ax2.tick_params(axis='y', labelcolor='#f43f5e')
ax2_twin.tick_params(axis='y', labelcolor='#fbbf24')

# --- Panel 3: Current sheet aspect ratio ---
ax3 = axes[1, 0]
Rm_range = np.logspace(4, 14, 200)
delta_over_L_SP = 1.0 / np.sqrt(Rm_range)

# Petschek: much thicker effective region
delta_over_L_Pet = np.pi / (8 * np.log(Rm_range))

ax3.loglog(Rm_range, delta_over_L_SP, color='#60a5fa', linewidth=2.5, label='Sweet-Parker: delta/L = 1/sqrt(Rm)')
ax3.loglog(Rm_range, delta_over_L_Pet, color='#f43f5e', linewidth=2.5, label='Petschek: l*/L (effective)')

ax3.set_xlabel('Magnetic Reynolds Number Rm')
ax3.set_ylabel('Current sheet aspect ratio delta/L')
ax3.set_title('Current Sheet Thinning')
ax3.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# Annotate physical thicknesses
ax3.annotate('delta ~ 1 m\n(SP, coronal Rm)', xy=(1e12, 1e-6), fontsize=9, color='#60a5fa',
            ha='center')

# --- Panel 4: Fermi power-law spectrum ---
ax4 = axes[1, 1]
E = np.logspace(0, 3, 500)  # keV

for delta_spec, col, lbl in [(2.0, '#f43f5e', 'delta=2 (strong shock, r=4)'),
                               (3.0, '#fb923c', 'delta=3 (r=2)'),
                               (4.0, '#fbbf24', 'delta=4 (r=1.5)'),
                               (5.0, '#4ade80', 'delta=5 (r=1.33)')]:
    spectrum = E**(-delta_spec) * np.exp(-E / 500)
    ax4.loglog(E, spectrum / spectrum[0], color=col, linewidth=2.5, label=lbl)

ax4.set_xlabel('Energy [keV]')
ax4.set_ylabel('f(E) [normalized]')
ax4.set_title('Fermi Acceleration: Power-Law Spectra')
ax4.set_ylim(1e-10, 10)
ax4.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')
ax4.text(50, 1e-3, 'dN/dE ~ E^{-delta}', fontsize=11, color='white')