Part IV: Space Weather | Chapter 14

Solar Energetic Particles

Diffusive shock acceleration, first-order Fermi, Parker transport, and SEP event types

14.1 Diffusive Shock Acceleration (DSA)

Derivation 1: First-Order Fermi Acceleration at a Shock

DSA is the primary mechanism for accelerating particles to high energies at CME-driven shocks. Particles gain energy by bouncing back and forth across the shock.

Step 1. A particle crossing the shock from upstream to downstream sees the downstream plasma approaching at relative speed $\Delta v = v_1 - v_2 = v_1(1 - 1/s)$ where$s = \rho_2/\rho_1$ is the compression ratio.

Step 2. The average fractional energy gain per shock crossing cycle:

$$\frac{\Delta E}{E} = \frac{4}{3}\frac{v_1 - v_2}{c} = \frac{4}{3}\frac{v_1}{c}\frac{s-1}{s}$$

Step 3. The probability of escape per cycle is $P_{\text{esc}} = 4v_2/(cv_s)$ for relativistic particles. The resulting spectrum is a power law:

$$\boxed{f(p) \propto p^{-q}, \quad q = \frac{3s}{s-1}}$$

For a strong shock with $s = 4$ (adiabatic, $\gamma = 5/3$):$q = 4$, giving a differential energy spectrum $dN/dE \propto E^{-2}$. This is remarkably close to observed SEP spectral indices.

14.2 Parker Transport Equation

Derivation 2: Focused Transport in the Heliosphere

Step 1. The Parker transport equation describes the propagation of energetic particles in the heliosphere including diffusion, convection, and energy changes:

$$\boxed{\frac{\partial f}{\partial t} = \nabla \cdot (\kappa \cdot \nabla f) - \mathbf{v}_{\text{sw}} \cdot \nabla f + \frac{1}{3}(\nabla \cdot \mathbf{v}_{\text{sw}})\frac{\partial f}{\partial \ln p} + Q}$$

The terms represent: (1) spatial diffusion along and across $\mathbf{B}$, (2) convection with solar wind, (3) adiabatic energy changes, and (4) source term.

Step 2. The parallel mean free path depends on the turbulence spectrum:

$$\lambda_\parallel = \frac{3v}{8}\int_0^1 \frac{(1-\mu^2)^2}{D_{\mu\mu}} d\mu$$

where $D_{\mu\mu}$ is the pitch-angle diffusion coefficient from wave-particle interactions. Typical values: $\lambda_\parallel \sim 0.1\text{--}1$ AU for ~10 MeV protons.

14.3 SEP Event Types

Derivation 3: Impulsive vs Gradual Events

Step 1. Impulsive events: accelerated in flare reconnection regions. Characterized by electron-rich, ${}^3\text{He}$-rich (${}^3\text{He}/{}^4\text{He} \sim 1$, vs 0.0004 in solar wind), heavy-ion enhanced. Duration: hours.

Step 2. Gradual events: accelerated at CME-driven shocks. Proton-rich, composition similar to corona/solar wind. Duration: days. Can fill the inner heliosphere.

Impulsive Events

• Source: Flare reconnection
• Duration: hours
• ${}^3\text{He}$ enriched ($\times 10^3$)
• Fe/O enhanced
• Electron-rich
• Small spatial extent

Gradual Events

• Source: CME-driven shock
• Duration: days
• Normal composition
• Proton-rich
• Fills heliosphere
• Radiation hazard

14.4 Maximum Energy and Acceleration Time

Derivation 4: Hillas Criterion

Step 1. A particle can be accelerated only as long as its gyroradius is smaller than the accelerator size:

$$\boxed{E_{\max} = ZeBR \approx 300 Z B[\mu\text{T}] R[\text{AU}] \text{ MeV}}$$

For a CME shock with $B \sim 100$ nT, $R \sim 0.5$ AU:$E_{\max} \sim 15$ GeV for protons, consistent with ground-level enhancement (GLE) events.

14.5 Acceleration Timescale

Derivation 5: DSA Acceleration Time

Step 1. The characteristic time for DSA to accelerate a particle to momentum $p$:

$$\boxed{\tau_{\text{acc}} = \frac{3}{v_1 - v_2}\left(\frac{\kappa_1}{v_1} + \frac{\kappa_2}{v_2}\right) \ln\frac{p}{p_0}}$$

For typical parameters ($v_1 = 1000$ km/s, $\kappa \sim 10^{18}$ m$^2$/s):$\tau_{\text{acc}} \sim 10^4$ s for 10 MeV protons, consistent with the observed onset times of gradual SEP events.

Numerical Simulation

SEPs: DSA Spectra, Intensity Profiles, Energy Gain, Composition

Python

script.py138 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 Energetic Particles', 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: DSA power law spectrum ---
ax1 = axes[0, 0]
E = np.logspace(-1, 4, 500)  # MeV

for s, col, lbl in [(4.0, '#f43f5e', 's=4 (strong shock, q=4)'),
                     (3.0, '#fb923c', 's=3 (q=4.5)'),
                     (2.5, '#fbbf24', 's=2.5 (q=5)')]:
    q = 3 * s / (s - 1)
    gamma_spec = (q + 2) / 2  # dN/dE spectral index
    dNdE = E**(-gamma_spec)
    # Exponential cutoff at E_max
    E_max = 1000 * (s / 4)**2  # MeV
    dNdE *= np.exp(-E / E_max)
    dNdE = dNdE / dNdE[0]
    ax1.loglog(E, dNdE, color=col, linewidth=2.5, label=lbl)

ax1.set_xlabel('Energy [MeV]')
ax1.set_ylabel('dN/dE [relative]')
ax1.set_title('DSA Energy Spectra')
ax1.set_ylim(1e-12, 10)
ax1.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: SEP intensity-time profiles ---
ax2 = axes[0, 1]
t = np.linspace(0, 72, 500)  # hours

# Impulsive event
imp_rise = np.exp(-((t - 2)/0.5)**2) * np.exp(-t/3)
imp_rise = np.maximum(imp_rise, 0)

# Gradual event
grad_rise = np.where(t > 5, (t - 5)**0.5 * np.exp(-(t - 5)/20), 0)
grad_rise = grad_rise / np.max(grad_rise)

# Shock passage spike
shock_spike = 0.3 * np.exp(-((t - 40)/2)**2)

ax2.semilogy(t, imp_rise * 100 + 0.1, color='#f43f5e', linewidth=2.5, label='Impulsive (electrons)')
ax2.semilogy(t, (grad_rise + shock_spike) * 100 + 0.1, color='#60a5fa', linewidth=2.5,
            label='Gradual (protons)')
ax2.axvline(x=40, color='#fbbf24', linestyle=':', alpha=0.5, label='Shock arrival')
ax2.set_xlabel('Time [hours after flare]')
ax2.set_ylabel('Intensity [relative]')
ax2.set_title('SEP Intensity-Time Profiles')
ax2.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: Fermi acceleration cartoon ---
ax3 = axes[1, 0]
# Particle trajectory bouncing between upstream and downstream
np.random.seed(42)
x = [0]
y = [0]
E_particle = [1]
for i in range(50):
    # Upstream diffusion
    dx = np.random.normal(-2, 1)
    dy = np.random.normal(0, 0.5)
    x.append(x[-1] + dx)
    y.append(y[-1] + dy)
    # Cross shock (at x=0)
    x.append(0)
    y.append(y[-1] + np.random.normal(0, 0.3))
    E_particle.append(E_particle[-1] * 1.05)
    E_particle.append(E_particle[-1])
    # Downstream diffusion
    dx = np.random.normal(1, 0.5)
    dy = np.random.normal(0, 0.5)
    x.append(x[-1] + dx)
    y.append(y[-1] + dy)
    E_particle.append(E_particle[-1])
    # Cross shock back
    x.append(0)
    y.append(y[-1] + np.random.normal(0, 0.3))
    E_particle.append(E_particle[-1] * 1.05)

crossings = np.arange(len(E_particle))
ax3.plot(crossings[:100], E_particle[:100], color='#f43f5e', linewidth=1.5)
ax3.set_xlabel('Number of shock interactions')
ax3.set_ylabel('Energy [relative]')
ax3.set_title('Energy Gain via DSA')
ax3.axhline(y=1, color='white', linestyle=':', alpha=0.3)

# --- Panel 4: Composition comparison ---
ax4 = axes[1, 1]
elements = ['H', 'He', '3He', 'C', 'O', 'Ne', 'Fe']
sw_abund = [1e6, 5e4, 0.02, 0.3, 0.5, 0.1, 0.05]
imp_abund = [1e5, 5e4, 20, 1, 0.5, 0.3, 0.5]
grad_abund = [1e6, 5e4, 0.03, 0.4, 0.6, 0.12, 0.06]

x_pos = np.arange(len(elements))
width = 0.25

ax4.bar(x_pos - width, sw_abund, width, color='#4ade80', alpha=0.8, label='Solar wind')
ax4.bar(x_pos, imp_abund, width, color='#f43f5e', alpha=0.8, label='Impulsive SEP')
ax4.bar(x_pos + width, grad_abund, width, color='#60a5fa', alpha=0.8, label='Gradual SEP')
ax4.set_yscale('log')
ax4.set_xticks(x_pos)
ax4.set_xticklabels(elements, color='white')
ax4.set_ylabel('Relative abundance')
ax4.set_title('SEP Composition')
ax4.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='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 ENERGETIC PARTICLES")
print("=" * 60)
print("\n--- DSA Theory ---")
print("  Strong shock (s=4): spectral index q = 4, dN/dE ~ E^-2")
print("  Acceleration time (10 MeV): ~hours")
print("  Maximum energy: ~GeV (protons)")
print("\n--- Event Types ---")
print("  Impulsive: flare-accelerated, electron/3He-rich, hours")
print("  Gradual: shock-accelerated, proton-rich, days")
print("  Ground-level events (GLEs): ~70 since 1942")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

14.6 Detailed DSA Power-Law Derivation

Deriving $f(p) \propto p^{-3r/(r-1)}$ from Compression Ratio

Step 1. Consider a plane shock with upstream velocity $u_1$ and downstream velocity $u_2 = u_1/r$ where $r$ is the compression ratio. In the shock frame, particles scatter isotropically in each frame.

Step 2. A particle crossing from upstream to downstream sees the downstream gas approaching at speed $\Delta u = u_1 - u_2 = u_1(1-1/r)$. The fractional momentum gain per crossing (averaging over pitch angles):

$$\left\langle\frac{\Delta p}{p}\right\rangle_{\text{up}\to\text{down}} = \frac{2}{3}\frac{\Delta u}{v}$$

Step 3. Similarly crossing back from downstream to upstream gives the same gain. Per complete cycle (up$\to$down$\to$up):

$$\left\langle\frac{\Delta p}{p}\right\rangle_{\text{cycle}} = \frac{4}{3}\frac{u_1 - u_2}{v} = \frac{4}{3}\frac{u_1}{v}\frac{r-1}{r}$$

Step 4. The escape probability per cycle. Only particles moving away from the shock downstream can escape. The flux of particles crossing the shock from downstream to upstream vs escaping downstream:

$$P_{\text{esc}} = \frac{4u_2}{v} = \frac{4u_1}{rv}$$

Step 5. After $n$ cycles: $p = p_0(1+\epsilon)^n$ and$N \propto (1-P_{\text{esc}})^n$. Eliminating $n$:

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

Step 6. For small $\epsilon$ and $P_{\text{esc}}$:$\ln(1-P_{\text{esc}})/\ln(1+\epsilon) \approx -P_{\text{esc}}/\epsilon$. Computing the ratio:

$$\frac{P_{\text{esc}}}{\epsilon} = \frac{4u_1/(rv)}{(4/3)(u_1/v)(r-1)/r} = \frac{3}{r-1}$$

Step 7. Therefore the differential spectrum:

$$\boxed{f(p) = \frac{dN}{dp} \propto p^{-q}, \quad q = \frac{3r}{r-1}}$$

For a strong adiabatic shock ($\gamma=5/3$): $r=4$, $q=4$, giving a differential energy spectrum $dN/dE \propto E^{-2}$. This is remarkably universal: it depends only on the compression ratio, not on the details of the scattering or the shock speed.

14.7 Parker Transport Equation: Full Derivation

Physical Origin of Each Term

Step 1. Start from the focused transport equation for the phase-space density$f(\mathbf{x}, p, t)$ of energetic particles in the solar wind:

$$\boxed{\frac{\partial f}{\partial t} = \underbrace{\nabla\cdot(\kappa\cdot\nabla f)}_{\text{diffusion}} - \underbrace{\mathbf{v}_{\text{sw}}\cdot\nabla f}_{\text{convection}} + \underbrace{\frac{1}{3}(\nabla\cdot\mathbf{v}_{\text{sw}})\,p\frac{\partial f}{\partial p}}_{\text{adiabatic cooling}} + \underbrace{Q}_{\text{source}}}$$

Term 1: Spatial diffusion. The diffusion tensor $\kappa_{ij}$ has parallel and perpendicular components relative to $\mathbf{B}$:

$$\kappa_\parallel = \frac{v\lambda_\parallel}{3}, \qquad \kappa_\perp \approx \frac{\kappa_\parallel}{1 + (\lambda_\parallel/r_g)^2}$$

Term 2: Convection. The solar wind carries particles outward at $v_{\text{sw}} \sim 400$ km/s.

Term 3: Adiabatic cooling. In the expanding solar wind ($\nabla\cdot\mathbf{v}_{\text{sw}} = 2v_{\text{sw}}/r > 0$), particles lose energy as they do work against the expanding flow. This term causes spectral steepening at low energies.

Step 2. For a radial solar wind in 1D (neglecting perpendicular diffusion):

$$\frac{\partial f}{\partial t} = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\kappa_\parallel\frac{\partial f}{\partial r}\right) - v_{\text{sw}}\frac{\partial f}{\partial r} + \frac{2v_{\text{sw}}}{3r}\,p\frac{\partial f}{\partial p} + Q$$

This is the standard equation solved in SEP transport codes. The parallel mean free path$\lambda_\parallel$ ranges from 0.01 AU (scatter-dominated) to 1 AU (nearly scatter-free) for 10-100 MeV protons, controlling the shape of the SEP intensity-time profile.

14.8 SEP Event Timeline: Impulsive vs Gradual

The diagram shows the characteristic intensity-time profiles of impulsive (flare) and gradual (CME-driven shock) SEP events.

Comparison of impulsive and gradual SEP events. Impulsive events: rapid onset (minutes), short duration (hours), electron-rich, ${}^3\text{He}$-enriched. Gradual events: slower onset, days-long duration, proton-rich, with ESP (energetic storm particle) enhancement at shock arrival.

Extended Simulation: DSA Spectra & SEP Intensity Profiles

Extended: DSA Spectra, Spectral Index, SEP Transport, Energy-Dependent Profiles

Python

script.py144 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('DSA Spectra and SEP Transport', fontsize=16, fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

# --- Panel 1: DSA spectra for different compression ratios ---
ax1 = axes[0, 0]
p = np.logspace(-1, 4, 500)  # momentum in units of p0

for r, col, lbl in [(4.0, '#f43f5e', 'r=4.0 (strong, q=4.0)'),
                     (3.5, '#fb923c', 'r=3.5 (q=4.2)'),
                     (3.0, '#fbbf24', 'r=3.0 (q=4.5)'),
                     (2.5, '#4ade80', 'r=2.5 (q=5.0)'),
                     (2.0, '#60a5fa', 'r=2.0 (q=6.0)')]:
    q = 3 * r / (r - 1)
    fp = p**(-q) * np.exp(-p / 1000)
    ax1.loglog(p, fp / fp[10], color=col, linewidth=2.5, label=lbl)

ax1.set_xlabel('Momentum p/p_0')
ax1.set_ylabel('f(p) [normalized]')
ax1.set_title('DSA: f(p) ~ p^{-3r/(r-1)}')
ax1.set_ylim(1e-15, 10)
ax1.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: Spectral index vs compression ratio ---
ax2 = axes[0, 1]
r_range = np.linspace(1.1, 7, 300)
q_range = 3 * r_range / (r_range - 1)
# Energy spectral index: gamma = (q+2)/2 for non-relativistic
gamma_nr = (q_range + 2) / 2

ax2.plot(r_range, q_range, color='#f43f5e', linewidth=2.5, label='Momentum index q')
ax2.plot(r_range, gamma_nr, color='#60a5fa', linewidth=2.5, label='Energy index gamma (NR)')

# Mark strong shock
ax2.axvline(x=4, color='#fbbf24', linestyle='--', alpha=0.5, label='r=4 (strong shock)')
ax2.axhline(y=4, color='#fbbf24', linestyle=':', alpha=0.3)

# Typical observed range
ax2.axhspan(3.5, 5.5, alpha=0.05, color='white')
ax2.text(5.5, 4.5, 'Typical\nobserved range', fontsize=9, color='white', ha='center')

ax2.set_xlabel('Compression ratio r')
ax2.set_ylabel('Spectral index')
ax2.set_title('DSA Spectral Index vs Compression Ratio')
ax2.set_xlim(1.5, 7)
ax2.set_ylim(1, 12)
ax2.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: SEP diffusive transport profiles ---
ax3 = axes[1, 0]
t = np.linspace(0.1, 72, 500)  # hours

# Solve diffusion-convection for different mean free paths
for lam_AU, col, lbl in [(0.05, '#f43f5e', 'lambda=0.05 AU (scatter-dom.)'),
                           (0.2, '#fb923c', 'lambda=0.2 AU'),
                           (0.5, '#fbbf24', 'lambda=0.5 AU'),
                           (1.0, '#4ade80', 'lambda=1.0 AU (scatter-free)')]:
    # Approximate solution: Reid-Axford profile
    # J(t) ~ (1/t) * exp(-r^2/(4*kappa*t))
    r_AU = 1.2  # path length along field line
    v_particle = 0.1 * 3e5  # 0.1c for ~5 MeV proton (km/s)
    kappa = v_particle * lam_AU * 1.496e8 / 3  # km^2/s
    t_sec = t * 3600

# Free-streaming arrival
    t_fs = r_AU * 1.496e8 / v_particle / 3600  # hours

J = np.where(t > t_fs * 0.5,
                 (1.0 / np.maximum(t_sec, 1)**1.5) * np.exp(-np.minimum(((r_AU * 1.496e8)**2) / (4 * kappa * np.maximum(t_sec, 1)), 500)),
                 0)
    J = J / np.nanmax(J)

ax3.semilogy(t, J + 1e-6, color=col, linewidth=2, label=lbl)

ax3.set_xlabel('Time [hours after injection]')
ax3.set_ylabel('Intensity [normalized]')
ax3.set_title('SEP Transport: Effect of Mean Free Path')
ax3.set_ylim(1e-4, 2)
ax3.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 4: Energy-dependent onset and decay ---
ax4 = axes[1, 1]
energies = [100, 30, 10, 3]  # MeV
colors_E = ['#f43f5e', '#fb923c', '#fbbf24', '#4ade80']

t = np.linspace(0, 48, 500)

for E_MeV, col in zip(energies, colors_E):
    # Higher energy -> faster arrival, faster decay
    v = np.sqrt(2 * E_MeV * 1.6e-13 / (1.67e-27))  # m/s
    t_onset = min(1.2 * 1.496e11 / v / 3600, 40)  # hours, capped
    tau_decay = max(10 / np.sqrt(E_MeV / 10), 1)  # hours

dt = t - t_onset
    rise = np.exp(-np.minimum(np.abs(dt - 2)**2 / 8, 100))
    decay = np.exp(-np.minimum(np.abs(dt) / tau_decay, 100))
    J = np.where(t > t_onset, rise * decay, 0)
    J = J / np.max(J) if np.max(J) > 0 else J
    ax4.semilogy(t, J + 1e-5, color=col, linewidth=2.5,
                label=str(E_MeV) + ' MeV')

ax4.set_xlabel('Time [hours]')
ax4.set_ylabel('Intensity [normalized]')
ax4.set_title('Energy-Dependent SEP Profiles')
ax4.set_ylim(1e-4, 2)
ax4.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')
ax4.text(3, 0.5, 'High energy:\nfaster onset', fontsize=9, color='#f43f5e')

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("EXTENDED: DSA SPECTRA AND SEP TRANSPORT")
print("=" * 60)
print("\n--- DSA Power Law ---")
print("  f(p) ~ p^{-3r/(r-1)}")
print("  r=4 (strong shock): q=4, dN/dE ~ E^{-2}")
print("  r=2 (weak shock): q=6, dN/dE ~ E^{-4}")
print("\n--- Parker Transport ---")
print("  Diffusion + convection + adiabatic cooling + source")
print("  lambda_parallel ~ 0.05-1 AU for 10 MeV protons")
print("  kappa = v*lambda/3")
print("\n--- Velocity Dispersion ---")
print("  100 MeV: v ~ 0.4c, onset ~ 12 min")
print("  10 MeV: v ~ 0.15c, onset ~ 30 min")
print("  1 MeV: v ~ 0.05c, onset ~ 100 min")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← CMEs Next: Space Weather →

Solar Energetic Particles

14.1 Diffusive Shock Acceleration (DSA)

Derivation 1: First-Order Fermi Acceleration at a Shock

14.2 Parker Transport Equation

Derivation 2: Focused Transport in the Heliosphere

14.3 SEP Event Types

Derivation 3: Impulsive vs Gradual Events

Impulsive Events

Gradual Events

14.4 Maximum Energy and Acceleration Time

Derivation 4: Hillas Criterion

14.5 Acceleration Timescale

Derivation 5: DSA Acceleration Time

Numerical Simulation

SEPs: DSA Spectra, Intensity Profiles, Energy Gain, Composition

14.6 Detailed DSA Power-Law Derivation

Deriving \(f(p) \propto p^{-3r/(r-1)}\) from Compression Ratio

14.7 Parker Transport Equation: Full Derivation

Physical Origin of Each Term

14.8 SEP Event Timeline: Impulsive vs Gradual

Extended Simulation: DSA Spectra & SEP Intensity Profiles

Extended: DSA Spectra, Spectral Index, SEP Transport, Energy-Dependent Profiles