Part IV: Space Weather | Chapter 13

Coronal Mass Ejections

CME morphology, flux rope models, Lorentz force driving, and the drag-based propagation model

13.1 CME Morphology and Properties

Derivation 1: CME Kinetic Energy and Mass

CMEs are large-scale eruptions of magnetized plasma from the solar corona. A typical CME has mass$m \sim 10^{12}\text{--}10^{13}$ kg and speed 400-2000 km/s.

Step 1. The kinetic energy of a CME:

$$\boxed{E_k = \frac{1}{2}mv^2 \sim \frac{1}{2}(10^{13})(10^6)^2 = 5 \times 10^{24} \text{ J}}$$

Step 2. The magnetic energy of the erupting flux rope (with $B \sim 0.01$ T, $R \sim 10^8$ m):

$$E_B = \frac{B^2}{2\mu_0}\frac{4}{3}\pi R^3 \sim 10^{25} \text{ J}$$

The three-part structure of a CME as seen in coronagraph images: bright leading edge (compressed sheath), dark cavity (flux rope), and bright core (erupting prominence material).

13.2 Flux Rope Model

Derivation 2: Toroidal Flux Rope Equilibrium

Step 1. The upward "hoop force" on a toroidal current-carrying flux rope:

$$F_{\text{hoop}} = \frac{\mu_0 I^2}{4\pi R}\left(\ln\frac{8R}{a} - 1 + \frac{l_i}{2}\right)$$

where $R$ is the major radius, $a$ is the minor radius, and $l_i$ is the internal inductance.

Step 2. This is balanced by the downward magnetic tension of the overlying strapping field $B_{\text{ext}}$:

$$\boxed{F_{\text{tension}} = \frac{B_{\text{ext}}^2}{2\mu_0} \cdot 2\pi R}$$

Step 3. The torus instability occurs when the external field decreases sufficiently fast with height. The critical decay index:

$$n = -\frac{d\ln B_{\text{ext}}}{d\ln h} > n_{\text{crit}} \approx 1.5$$

When the flux rope rises to a height where $n > 1.5$, the overlying field can no longer restrain it, and a CME erupts. This provides a clear, testable criterion for CME initiation.

13.3 Lorentz Force Driving

Derivation 3: Lorentz Self-Force on an Expanding Flux Rope

Step 1. The $j \times B$ force in the radial direction for a flux rope with axial field $B_z$ and azimuthal field $B_\phi$:

$$\boxed{F_L = \frac{1}{\mu_0}\left(-\frac{B_\phi^2}{r} - \frac{d}{dr}\frac{B_\phi^2 + B_z^2}{2}\right)}$$

The first term is the pinch (inward tension of the azimuthal field) and the second is the magnetic pressure gradient. For a CME, the net Lorentz force is initially outward (dominated by the pressure gradient), accelerating the ejecta to speeds of 500-2000 km/s within $\sim 2\text{--}5 R_\odot$.

13.4 Drag-Based Propagation Model

Derivation 4: Aerodynamic Drag in the Solar Wind

Beyond $\sim 20 R_\odot$, the Lorentz force becomes negligible and the CME interacts aerodynamically with the ambient solar wind.

Step 1. The equation of motion:

$$\boxed{\frac{d^2R}{dt^2} = -\gamma(v - v_{\text{sw}})|v - v_{\text{sw}}|}$$

where $\gamma = c_d A \rho_{\text{sw}} / (M_{\text{CME}} + M_{\text{virtual}})$ is the drag parameter, typically $\gamma \sim 0.2\text{--}2 \times 10^{-7}$ km$^{-1}$.

Step 2. For constant $\gamma$ and $v_{\text{sw}}$, this has an analytical solution. Fast CMEs decelerate and slow CMEs accelerate toward the solar wind speed (~400 km/s).

Transit times to Earth range from ~15 hours (extreme events like the Carrington event, 1859) to ~5 days (slow CMEs). The drag-based model (DBM) is widely used for operational space weather forecasting.

13.5 ICME Structure at 1 AU

Derivation 5: Magnetic Cloud Fitting

Step 1. An interplanetary CME (ICME) at 1 AU often shows a magnetic cloud signature modeled as a force-free cylinder (Lundquist solution):

$$B_z(r) = B_0 J_0(\alpha r)$$$$\boxed{B_\phi(r) = B_0 H J_1(\alpha r)}$$

where $J_0, J_1$ are Bessel functions, $H = \pm 1$ is the chirality, and $\alpha = 2.405/R_0$. Typical parameters at 1 AU: $B_0 \sim 20$ nT,$R_0 \sim 0.1$ AU. The smooth rotation of the magnetic field vector through the cloud is the defining signature used for identification.

Numerical Simulation

CMEs: Drag-Based Propagation, Velocity Profiles, Magnetic Cloud, Torus Instability

Python

script.py121 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.special import j0, j1

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle('Coronal Mass Ejections', 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: Drag-based model ---
ax1 = axes[0, 0]
v_sw = 400  # km/s
gamma = 1e-7  # km^-1

def drag_model(t, y):
    R, v = y
    dv = -gamma * (v - v_sw) * abs(v - v_sw)
    return [v, dv]

for v0, col, lbl in [(2000, '#f43f5e', 'v0=2000 km/s'),
                      (1000, '#fb923c', 'v0=1000 km/s'),
                      (500, '#fbbf24', 'v0=500 km/s'),
                      (250, '#4ade80', 'v0=250 km/s')]:
    R_sun_km = 6.957e5
    sol = solve_ivp(drag_model, [0, 400], [20*R_sun_km, v0],
                   max_step=0.5, dense_output=True, t_eval=np.linspace(0, 400, 1000))
    AU_km = 1.496e8
    ax1.plot(sol.t, sol.y[0] / AU_km, color=col, linewidth=2, label=lbl)

ax1.axhline(y=1.0, color='white', linestyle=':', alpha=0.3, label='1 AU')
ax1.set_xlabel('Time [hours]')
ax1.set_ylabel('Distance [AU]')
ax1.set_title('Drag-Based CME Propagation')
ax1.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')
ax1.set_xlim(0, 200)

# --- Panel 2: CME velocity profiles ---
ax2 = axes[0, 1]
for v0, col, lbl in [(2000, '#f43f5e', 'Fast CME'),
                      (1000, '#fb923c', 'Average CME'),
                      (250, '#4ade80', 'Slow CME')]:
    sol = solve_ivp(drag_model, [0, 400], [20*6.957e5, v0],
                   max_step=0.5, t_eval=np.linspace(0, 300, 500))
    ax2.plot(sol.y[0] / (6.957e5), sol.y[1], color=col, linewidth=2.5, label=lbl)

ax2.axhline(y=v_sw, color='white', linestyle='--', alpha=0.5, label='v_sw = 400 km/s')
ax2.set_xlabel('Distance [R_sun]')
ax2.set_ylabel('Speed [km/s]')
ax2.set_title('CME Speed vs Distance')
ax2.set_xlim(20, 215)
ax2.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: Magnetic cloud (Lundquist solution) ---
ax3 = axes[1, 0]
r = np.linspace(0, 1, 200)  # normalized to cloud radius
alpha = 2.405

Bz = j0(alpha * r)
Bphi = j1(alpha * r)
B_total = np.sqrt(Bz**2 + Bphi**2)

ax3.plot(r, Bz, color='#f43f5e', linewidth=2.5, label='B_z (axial)')
ax3.plot(r, Bphi, color='#60a5fa', linewidth=2.5, label='B_phi (azimuthal)')
ax3.plot(r, B_total, color='#fbbf24', linewidth=2, linestyle='--', label='|B| total')
ax3.axhline(y=0, color='white', linestyle=':', alpha=0.3)
ax3.set_xlabel('r / R_0')
ax3.set_ylabel('B / B_0')
ax3.set_title('Magnetic Cloud: Lundquist Solution')
ax3.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 4: Torus instability ---
ax4 = axes[1, 1]
h = np.linspace(0.1, 3, 200)  # height in units of source depth

# External field decay profiles
for n_val, col, lbl in [(0.5, '#4ade80', 'n=0.5 (stable)'),
                         (1.0, '#fbbf24', 'n=1.0 (marginal)'),
                         (1.5, '#fb923c', 'n=1.5 (critical)'),
                         (2.0, '#f43f5e', 'n=2.0 (unstable)')]:
    B_ext = h**(-n_val)
    ax4.plot(h, B_ext / B_ext[0], color=col, linewidth=2.5, label=lbl)

ax4.axhline(y=0, color='white', linestyle=':', alpha=0.3)
ax4.set_xlabel('Height [normalized]')
ax4.set_ylabel('B_ext / B_ext(0)')
ax4.set_title('Torus Instability: Field Decay Index')
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("CORONAL MASS EJECTIONS")
print("=" * 60)
print("\n--- CME Properties ---")
print("  Mass: 10^12 - 10^13 kg")
print("  Speed: 200 - 3000 km/s")
print("  Kinetic energy: 10^23 - 10^25 J")
print("  Rate: ~3/day (solar max), ~0.5/day (solar min)")
print("\n--- Transit Times to Earth ---")
print("  Fast (2000 km/s): ~20 hours")
print("  Average (1000 km/s): ~40 hours")
print("  Slow (400 km/s): ~100 hours")
print("  Carrington event (1859): ~17 hours")
print("\n--- Torus Instability ---")
print("  Critical decay index n_crit ~ 1.5")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

13.6 Flux Rope Equilibrium: Grad-Shafranov Equation

Force-Free Flux Rope in Cylindrical Geometry

Step 1. For an axisymmetric flux rope with translational symmetry along $z$, the magnetic field can be written in terms of a flux function $\psi(r,z)$:

$$\mathbf{B} = \frac{1}{r}\nabla\psi\times\hat{z} + B_z(r)\hat{z}$$

Step 2. For a force-free cylindrical configuration ($\partial/\partial z = 0$,$\partial/\partial\phi = 0$), the force-free condition $\nabla\times\mathbf{B}=\alpha\mathbf{B}$in cylindrical coordinates gives:

$$-\frac{dB_z}{dr} = \alpha B_\phi, \qquad \frac{1}{r}\frac{d(rB_\phi)}{dr} = \alpha B_z$$

Step 3. For constant $\alpha$, these combine into a Bessel equation:

$$\frac{d^2B_z}{dr^2} + \frac{1}{r}\frac{dB_z}{dr} + \alpha^2 B_z = 0$$

Step 4. The solution is the Lundquist flux rope:

$$\boxed{B_z(r) = B_0 J_0(\alpha r), \quad B_\phi(r) = B_0 H J_1(\alpha r)}$$

where $H = \pm 1$ is the handedness (chirality) of the flux rope, $J_0, J_1$are Bessel functions, and $\alpha = 2.405/R_0$ ensures $B_z$ vanishes at the boundary. This is the most widely used model for fitting magnetic clouds observed at 1 AU.

13.7 Torus Instability: Critical Decay Index

Deriving $n_{\text{crit}} = 3/2$

Step 1. A current ring of current $I$ at height $R$ above the photosphere experiences an upward hoop force:

$$F_{\text{hoop}} = \frac{\mu_0 I^2}{4\pi R}\left(\ln\frac{8R}{a} - 1 + \frac{l_i}{2}\right) \approx \frac{\mu_0 I^2}{4\pi R}\Lambda$$

Step 2. The external (strapping) field exerts a downward Lorentz force per unit length:

$$F_{\text{strap}} = I B_{\text{ext}}(R) \cdot 2\pi R$$

Step 3. Equilibrium: $F_{\text{hoop}} = F_{\text{strap}}$. Now perturb$R \to R + \delta R$ and check stability:

$$\frac{dF_{\text{hoop}}}{dR} \sim -\frac{F_{\text{hoop}}}{R}, \qquad \frac{dF_{\text{strap}}}{dR} = I\frac{d(B_{\text{ext}}\cdot 2\pi R)}{dR}$$

Step 4. Instability occurs when the restoring force decreases faster than the driving force, i.e., when $|dF_{\text{strap}}/dR| < |dF_{\text{hoop}}/dR|$. Defining the decay index:

$$n = -\frac{d\ln B_{\text{ext}}}{d\ln R}$$

Step 5. The stability condition becomes:

$$\boxed{n > n_{\text{crit}} = \frac{3}{2} \quad\text{(unstable, CME eruption)}}$$

Observations confirm that CMEs erupt when the flux rope apex reaches the height where$n \approx 1.5$. Typical eruption heights are 50,000-100,000 km above the photosphere. The torus instability provides the most successful quantitative criterion for CME initiation.

13.8 Drag-Based Propagation: Analytical Solution

Full Equation of Motion and Arrival Time

Step 1. The complete equation including gravity and drag:

$$\frac{d^2R}{dt^2} = -\gamma(v - v_{\text{sw}})|v - v_{\text{sw}}| - \frac{GM_\odot}{R^2}$$

For $R \gg R_\odot$, gravity is negligible and with constant $\gamma$ and $v_{\text{sw}}$:

Step 2. Define $w = v - v_{\text{sw}}$. Then $dw/dt = -\gamma w|w|$. For a fast CME ($w > 0$):

$$\frac{dw}{w^2} = -\gamma\,dt \quad\Rightarrow\quad w(t) = \frac{w_0}{1 + \gamma w_0 t}$$

Step 3. Integrating for position:

$$\boxed{R(t) = R_0 + v_{\text{sw}}\,t + \frac{1}{\gamma}\ln(1 + \gamma\,w_0\,t)}$$

Step 4. Setting $R = 1$ AU and solving for $t$ gives the arrival time. For a CME with $v_0 = 2000$ km/s, $\gamma = 10^{-7}$ km$^{-1}$: arrival in ~20 hours. For $v_0 = 500$ km/s: arrival in ~80 hours.

The drag parameter $\gamma = c_d A\rho_{\text{sw}}/(M + M_v)$ depends on the drag coefficient ($c_d \sim 1$), CME cross-section $A$, solar wind density, and the "virtual mass" (entrained ambient medium). Typical values: $\gamma \sim 0.2\text{--}2\times10^{-7}$km$^{-1}$.

13.9 CME Three-Part Structure

Coronagraph images reveal the classic three-part structure of a CME: bright front, dark cavity, and bright core.

Classic three-part CME structure as seen in coronagraph images: bright compressed sheath (leading edge), dark cavity (low-density flux rope with helical magnetic field), and bright core (cool, dense prominence material).

Extended Simulation: CME Propagation & Earth Arrival

Extended: CME Arrival Prediction, Decay Index, Speed Evolution at 1 AU

Python

script.py159 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle('CME Propagation and Arrival Time Prediction', fontsize=16, fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

AU_km = 1.496e8
R_sun_km = 6.957e5

# --- Panel 1: Analytical vs numerical drag model ---
ax1 = axes[0, 0]
v_sw = 400  # km/s
gamma = 1e-7

def drag_ode(t, y):
    R, v = y
    w = v - v_sw
    dv = -gamma * w * abs(w)
    return [v, dv]

for v0, col, lbl in [(2500, '#f43f5e', 'v0=2500 (Carrington-class)'),
                      (1500, '#fb923c', 'v0=1500 km/s'),
                      (800, '#fbbf24', 'v0=800 km/s'),
                      (300, '#4ade80', 'v0=300 km/s (slow)')]:
    R0 = 20 * R_sun_km
    sol = solve_ivp(drag_ode, [0, 300], [R0, v0],
                   max_step=0.5, t_eval=np.linspace(0, 300, 500))

# Analytical solution
    w0 = v0 - v_sw
    t_arr = sol.t
    if abs(w0) > 1:
        sign_w = np.sign(w0)
        R_analytical = R0 + v_sw * t_arr * 3600 + sign_w / gamma * np.log(1 + gamma * abs(w0) * t_arr * 3600)
    else:
        R_analytical = R0 + v_sw * t_arr * 3600

ax1.plot(sol.t, sol.y[0] / AU_km, color=col, linewidth=2, label=lbl)
    ax1.plot(t_arr, R_analytical / AU_km, '--', color=col, linewidth=1, alpha=0.5)

ax1.axhline(y=1.0, color='white', linestyle=':', alpha=0.3)
ax1.text(5, 1.03, '1 AU (Earth)', fontsize=9, color='white')
ax1.set_xlabel('Time [hours]')
ax1.set_ylabel('Distance [AU]')
ax1.set_title('CME Propagation (solid=numerical, dashed=analytical)')
ax1.set_xlim(0, 150)
ax1.set_ylim(0, 2)
ax1.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: Arrival time vs initial speed ---
ax2 = axes[0, 1]
v0_range = np.linspace(200, 3000, 200)
arrival_times = []

for v0 in v0_range:
    R0 = 20 * R_sun_km
    try:
        sol = solve_ivp(drag_ode, [0, 500], [R0, v0], max_step=1.0,
                       events=lambda t, y: y[0] - AU_km, dense_output=True)
        if sol.t_events[0].size > 0:
            arrival_times.append(sol.t_events[0][0])
        else:
            arrival_times.append(np.nan)
    except Exception:
        arrival_times.append(np.nan)

ax2.plot(v0_range, arrival_times, color='#f43f5e', linewidth=2.5)
ax2.axhline(y=17, color='#fbbf24', linestyle='--', alpha=0.5, label='Carrington event (~17 hr)')
ax2.axhline(y=96, color='#4ade80', linestyle='--', alpha=0.5, label='Typical slow CME (~4 days)')
ax2.set_xlabel('Initial CME speed [km/s]')
ax2.set_ylabel('Arrival time at 1 AU [hours]')
ax2.set_title('CME Arrival Time vs Initial Speed')
ax2.set_ylim(0, 200)
ax2.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: Decay index profile ---
ax3 = axes[1, 0]
h = np.linspace(10, 200, 300)  # height in Mm

# Different active region configurations
configs = [
    ('Simple bipole', 1.5, '#f43f5e'),
    ('Quadrupolar', 2.5, '#fb923c'),
    ('Dipole + nearby AR', 1.0, '#fbbf24'),
]

for name, n_base, col in configs:
    # n increases with height
    n_profile = n_base * (h / 100)**0.5
    ax3.plot(h, n_profile, color=col, linewidth=2.5, label=name)

ax3.axhline(y=1.5, color='white', linestyle='--', linewidth=2, alpha=0.5, label='n_crit = 1.5')
ax3.fill_between(h, 1.5, 4, alpha=0.05, color='#f43f5e')
ax3.text(150, 2.5, 'Unstable\n(CME eruption)', fontsize=10, color='#f43f5e', ha='center')
ax3.text(30, 0.5, 'Stable\n(confined)', fontsize=10, color='#4ade80', ha='center')
ax3.set_xlabel('Height above photosphere [Mm]')
ax3.set_ylabel('Decay index n')
ax3.set_title('Torus Instability: Decay Index Profiles')
ax3.set_ylim(0, 4)
ax3.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 4: CME speed at 1 AU vs initial speed ---
ax4 = axes[1, 1]
v_final = []
for v0 in v0_range:
    R0 = 20 * R_sun_km
    try:
        sol = solve_ivp(drag_ode, [0, 500], [R0, v0], max_step=1.0,
                       events=lambda t, y: y[0] - AU_km)
        if sol.t_events[0].size > 0:
            idx = -1
            v_final.append(sol.y_events[0][0][1])
        else:
            v_final.append(np.nan)
    except Exception:
        v_final.append(np.nan)

ax4.plot(v0_range, v_final, color='#f43f5e', linewidth=2.5, label='v at 1 AU')
ax4.plot([200, 3000], [200, 3000], 'w--', linewidth=1, alpha=0.3, label='v_1AU = v_0 (no drag)')
ax4.axhline(y=v_sw, color='#4ade80', linestyle=':', alpha=0.5, label='v_sw = 400 km/s')
ax4.set_xlabel('Initial CME speed [km/s]')
ax4.set_ylabel('Speed at 1 AU [km/s]')
ax4.set_title('CME Speed Evolution: Initial vs Arrival')
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("EXTENDED: CME PROPAGATION AND ARRIVAL")
print("=" * 60)
print("\n--- Drag-Based Model ---")
print("  R(t) = R0 + v_sw*t + (1/gamma)*ln(1 + gamma*w0*t)")
print("  gamma ~ 0.2-2 x 10^-7 km^-1")
print("\n--- Torus Instability ---")
print("  n_crit = 3/2 (for thin current ring)")
print("  Typical eruption height: 50-100 Mm")
print("\n--- Arrival Time Examples ---")
print("  Carrington (v0~2500): ~15-17 hours")
print("  Fast CME (v0~1500): ~25 hours")
print("  Average (v0~800): ~50 hours")
print("  Slow (v0~300): ~120+ hours")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Solar Flares Next: Solar Energetic Particles →

Coronal Mass Ejections

13.1 CME Morphology and Properties

Derivation 1: CME Kinetic Energy and Mass

13.2 Flux Rope Model

Derivation 2: Toroidal Flux Rope Equilibrium

13.3 Lorentz Force Driving

Derivation 3: Lorentz Self-Force on an Expanding Flux Rope

13.4 Drag-Based Propagation Model

Derivation 4: Aerodynamic Drag in the Solar Wind

13.5 ICME Structure at 1 AU

Derivation 5: Magnetic Cloud Fitting

Numerical Simulation

CMEs: Drag-Based Propagation, Velocity Profiles, Magnetic Cloud, Torus Instability

13.6 Flux Rope Equilibrium: Grad-Shafranov Equation

Force-Free Flux Rope in Cylindrical Geometry

13.7 Torus Instability: Critical Decay Index

Deriving \(n_{\text{crit}} = 3/2\)

13.8 Drag-Based Propagation: Analytical Solution

Full Equation of Motion and Arrival Time

13.9 CME Three-Part Structure

Extended Simulation: CME Propagation & Earth Arrival

Extended: CME Arrival Prediction, Decay Index, Speed Evolution at 1 AU