Part IV: Space Weather | Chapter 16

Sun-Earth Connection

Parker spiral, sector structure, magnetopause standoff, Dungey cycle, and ring current injection

16.1 The Parker Spiral

Derivation 1: Spiral Angle of the Interplanetary Magnetic Field

As the solar wind carries the frozen-in magnetic field outward while the Sun rotates, the field lines form an Archimedean spiral in the equatorial plane.

Step 1. In the rotating frame of the Sun (angular velocity $\Omega_\odot$), a radially emitted plasma parcel maintains a fixed azimuthal position. In the inertial frame, the azimuthal displacement at radius $r$ is:

$$\phi(r) = \phi_0 - \frac{\Omega_\odot (r - R_\odot)}{v_{\text{sw}}}$$

Step 2. Since the magnetic field is frozen to the plasma, the field line shape follows the plasma path. The radial and azimuthal components of $\mathbf{B}$ satisfy:

$$\frac{B_\phi}{B_r} = \frac{r \, d\phi}{dr} = -\frac{\Omega_\odot r}{v_{\text{sw}}}$$

Step 3. The spiral angle $\psi$ (angle between $\mathbf{B}$and the radial direction):

$$\boxed{\tan\psi = \frac{|B_\phi|}{B_r} = \frac{\Omega_\odot r}{v_{\text{sw}}}}$$

Step 4. At Earth (1 AU = $1.496 \times 10^{11}$ m):

$$\tan\psi = \frac{2.87 \times 10^{-6} \times 1.496 \times 10^{11}}{4 \times 10^5} \approx 1.07 \quad \Longrightarrow \quad \psi \approx 45^\circ$$

The garden-hose angle of ~45 degrees at 1 AU is well confirmed by in-situ measurements. For fast wind (700 km/s): $\psi \approx 28^\circ$; for slow wind (350 km/s): $\psi \approx 51^\circ$.

16.2 Heliospheric Current Sheet and Sector Structure

Derivation 2: Magnetic Sector Boundaries

Step 1. The solar magnetic dipole axis is tilted with respect to the rotation axis by angle $\alpha$. The heliospheric current sheet (HCS) separates regions of opposite magnetic polarity and is warped according to:

$$\boxed{\theta_{\text{HCS}}(\phi) = \frac{\pi}{2} + \arcsin\left(\sin\alpha\sin\left(\phi - \frac{\Omega r}{v_{\text{sw}}}\right)\right)}$$

Step 2. As Earth orbits the Sun (or equivalently, the Sun rotates), it passes through the warped HCS, alternating between "toward" and "away" magnetic sectors. The number of sectors depends on the complexity of the coronal field.

During solar minimum, $\alpha$ is small and the HCS is nearly flat (two-sector structure). During solar maximum, $\alpha$ can reach 70+ degrees, creating a highly warped "ballerina skirt" with four or more sectors.

16.3 Magnetopause Standoff Distance

Derivation 3: Pressure Balance at the Magnetopause

Step 1. The magnetopause is the boundary where the solar wind dynamic pressure balances the Earth's magnetic pressure:

$$\rho_{\text{sw}} v_{\text{sw}}^2 = \frac{B_E^2}{2\mu_0}\left(\frac{R_E}{r}\right)^6$$

Step 2. Solving for the standoff distance (sub-solar point):

$$\boxed{r_{\text{mp}} = R_E\left(\frac{B_E^2}{2\mu_0 \rho_{\text{sw}} v_{\text{sw}}^2}\right)^{1/6}}$$

Step 3. For typical solar wind ($n = 5$ cm$^{-3}$,$v = 400$ km/s, $B_E = 3.1 \times 10^{-5}$ T):

$$r_{\text{mp}} \approx 10.3 \, R_E$$

During intense storms with $n = 20$ cm$^{-3}$ and $v = 800$ km/s, the magnetopause can be compressed to $\sim 6 R_E$, exposing geosynchronous satellites ($6.6 R_E$) to the solar wind.

16.4 The Dungey Cycle

Derivation 4: Magnetic Reconnection-Driven Convection

The Dungey cycle (1961) describes the circulation of magnetic flux through the magnetosphere driven by reconnection at the magnetopause and in the magnetotail.

Step 1. When the IMF has a southward component ($B_z < 0$), reconnection occurs at the dayside magnetopause:

$$\Phi_{\text{reconnect}} = E_{\text{sw}} \ell_{\text{merging}} \approx v_{\text{sw}} B_z \ell$$

Step 2. The cross-polar cap potential, which drives magnetospheric convection:

$$\boxed{\Phi_{pc} \approx v_{\text{sw}} B_z \ell_{\text{eff}} \sim 50\text{--}200 \text{ kV (typical storm)}}$$

Step 3. The cycle consists of:

Dayside reconnection opens closed field lines
Open flux is swept antisunward over the poles by solar wind
Open flux accumulates in the magnetotail lobes
Tail reconnection at the X-line closes open flux
Newly closed flux convects earthward, energizing plasma
Flux circulates back to the dayside, completing the cycle

The Dungey cycle time is $\sim 1\text{--}4$ hours. During storms with sustained southward IMF, continuous reconnection drives intense convection, ring current injection, and auroral activity.

16.5 Ring Current Injection

Derivation 5: Adiabatic Particle Energization

Step 1. As flux tubes convect earthward in the Dungey cycle, they compress. Particles trapped on these flux tubes are adiabatically energized through conservation of the first adiabatic invariant:

$$\mu = \frac{p_\perp^2}{2mB} = \text{const}$$

Step 2. As a particle moves from $L_1$ to $L_2 < L_1$, the magnetic field increases as $B \propto L^{-3}$. The perpendicular energy increases:

$$\boxed{\frac{E_\perp(L_2)}{E_\perp(L_1)} = \frac{B(L_2)}{B(L_1)} = \left(\frac{L_1}{L_2}\right)^3}$$

Step 3. A 1 keV ion at $L = 10$ that is convected to $L = 4$gains energy by a factor of $(10/4)^3 \approx 16$, reaching ~16 keV. This populates the ring current at $L = 3\text{--}7$ with 10-200 keV ions (mainly$\text{H}^+$ and $\text{O}^+$).

The ring current particles gradient-curvature drift azimuthally (ions westward, electrons eastward), producing a net westward current that depresses the surface magnetic field (measured as Dst). The ring current decays through charge exchange with exospheric neutrals on a timescale of$\sim 1\text{--}10$ days.

Numerical Simulation

Sun-Earth Connection: Parker Spiral, Spiral Angle, Magnetopause, Adiabatic Energization

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('Sun-Earth Connection', 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: Parker spiral ---
ax1 = axes[0, 0]
Omega = 2.87e-6  # rad/s
R_sun_AU = 1.0 / 215.0

for v_sw, col, lbl in [(350, '#f43f5e', 'Slow wind (350 km/s)'),
                         (500, '#fb923c', 'Average (500 km/s)'),
                         (700, '#fbbf24', 'Fast wind (700 km/s)')]:
    v_AU = v_sw * 1e3 / 1.496e11  # AU/s
    r = np.linspace(R_sun_AU, 2.0, 500)  # AU
    phi = -Omega * (r - R_sun_AU) / v_AU

x = r * np.cos(phi)
    y = r * np.sin(phi)
    ax1.plot(x, y, color=col, linewidth=2, label=lbl)

# Additional spiral lines at different starting angles
    for dphi in [np.pi/2, np.pi, 3*np.pi/2]:
        x2 = r * np.cos(phi + dphi)
        y2 = r * np.sin(phi + dphi)
        ax1.plot(x2, y2, color=col, linewidth=0.8, alpha=0.4)

# Sun and Earth
ax1.plot(0, 0, 'o', color='#fbbf24', markersize=12)
theta_earth = 0
ax1.plot(np.cos(theta_earth), np.sin(theta_earth), 'o', color='#60a5fa', markersize=8)
ax1.annotate('Earth', xy=(1.05, 0.05), fontsize=8, color='#60a5fa')

ax1.set_xlim(-2, 2)
ax1.set_ylim(-2, 2)
ax1.set_aspect('equal')
ax1.set_xlabel('x [AU]')
ax1.set_ylabel('y [AU]')
ax1.set_title('Parker Spiral (Equatorial Plane)')
ax1.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white', loc='upper left')

# --- Panel 2: Spiral angle vs distance ---
ax2 = axes[0, 1]
r_AU = np.linspace(0.1, 5, 500)

for v_sw, col in [(350, '#f43f5e'), (500, '#fb923c'), (700, '#fbbf24')]:
    psi = np.degrees(np.arctan(Omega * r_AU * 1.496e11 / (v_sw * 1e3)))
    ax2.plot(r_AU, psi, color=col, linewidth=2.5, label='v_sw = ' + str(v_sw) + ' km/s')

ax2.axvline(x=1.0, color='white', linestyle=':', alpha=0.3, label='Earth')
ax2.axhline(y=45, color='white', linestyle='--', alpha=0.3)
ax2.set_xlabel('Distance [AU]')
ax2.set_ylabel('Spiral angle [deg]')
ax2.set_title('Parker Spiral Angle vs Distance')
ax2.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: Magnetopause standoff ---
ax3 = axes[1, 0]
v_range = np.linspace(300, 1000, 200)  # km/s
n_vals = [3, 5, 10, 20]  # cm^-3
colors_n = ['#4ade80', '#fbbf24', '#fb923c', '#f43f5e']
m_p = 1.673e-27
B_E = 3.1e-5
R_E = 6.371e6
mu_0 = 4 * np.pi * 1e-7

for n, col in zip(n_vals, colors_n):
    rho = n * 1e6 * m_p  # kg/m^3
    v = v_range * 1e3  # m/s
    P_dyn = 0.5 * rho * v**2  # Pascal (use 0.5 for efficiency, actual factor from empirical fit)
    r_mp = R_E * (B_E**2 / (2 * mu_0 * P_dyn))**(1.0/6.0)
    ax3.plot(v_range, r_mp / R_E, color=col, linewidth=2.5, label='n = ' + str(n) + ' cm^-3')

ax3.axhline(y=6.6, color='white', linestyle='--', alpha=0.3, label='GEO orbit (6.6 R_E)')
ax3.set_xlabel('Solar wind speed [km/s]')
ax3.set_ylabel('Magnetopause distance [R_E]')
ax3.set_title('Magnetopause Standoff Distance')
ax3.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 4: Adiabatic energization ---
ax4 = axes[1, 1]
L_start = np.linspace(6, 12, 200)
L_end_vals = [3, 4, 5]
colors_L = ['#f43f5e', '#fb923c', '#fbbf24']

E_init = 1.0  # keV

for L_end, col in zip(L_end_vals, colors_L):
    E_final = E_init * (L_start / L_end)**3
    ax4.semilogy(L_start, E_final, color=col, linewidth=2.5,
                label='L_final = ' + str(L_end))

ax4.axhline(y=50, color='white', linestyle=':', alpha=0.3)
ax4.annotate('Ring current energy range', xy=(6, 50), fontsize=8, color='white')
ax4.set_xlabel('Initial L-shell')
ax4.set_ylabel('Final energy [keV] (from 1 keV)')
ax4.set_title('Adiabatic Energization (mu conservation)')
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("SUN-EARTH CONNECTION")
print("=" * 60)
print("\n--- Parker Spiral ---")
psi_earth = np.degrees(np.arctan(Omega * 1.496e11 / (400e3)))
print("  Spiral angle at 1 AU (400 km/s): " + str(round(psi_earth, 1)) + " deg")
psi_fast = np.degrees(np.arctan(Omega * 1.496e11 / (700e3)))
print("  Spiral angle at 1 AU (700 km/s): " + str(round(psi_fast, 1)) + " deg")
print("\n--- Magnetopause ---")
n_typ = 5e6  # m^-3
rho_typ = n_typ * m_p
v_typ = 400e3
P_dyn = 0.5 * rho_typ * v_typ**2
r_mp = R_E * (B_E**2 / (2 * mu_0 * P_dyn))**(1.0/6.0)
print("  Standoff (n=5, v=400): " + str(round(r_mp/R_E, 1)) + " R_E")
print("\n--- Dungey Cycle ---")
print("  Cross-polar cap potential: 50-200 kV (storm)")
print("  Cycle time: 1-4 hours")
print("\n--- Ring Current ---")
print("  Energization factor (L=10 to L=4): " + str(round((10.0/4)**3, 1)) + "x")
print("  Main ions: H+ (80%), O+ (20% in storms)")
print("  Decay time: 1-10 days (charge exchange)")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

16.6 Full Parker Spiral Derivation

Step-by-Step from Frozen-in Field in the Rotating Frame

Step 1. In the frame co-rotating with the Sun at angular velocity $\Omega_\odot = 2.87\times10^{-6}$ rad/s, the solar wind is purely radial: $\mathbf{v}' = v_r\hat{r}$. The frozen-in condition means the magnetic field is parallel to the velocity in the co-rotating frame.

Step 2. Transform to the inertial frame. A point at $(r, \phi)$ in the rotating frame has velocity:

$$\mathbf{v} = v_r\hat{r} + \Omega_\odot r\sin\theta\,\hat{\phi}$$

Step 3. Since $\mathbf{B} \| \mathbf{v}'$ in the rotating frame, and$\mathbf{v}' = v_r\hat{r}$, we have $B_\phi' = 0$ in that frame. Transforming the electric field: $\mathbf{E}' = \mathbf{E} + \mathbf{v}_{\text{rot}}\times\mathbf{B}$. The frozen-in condition $\mathbf{E}' + \mathbf{v}'\times\mathbf{B} = 0$ gives:

$$\frac{B_\phi}{B_r} = -\frac{\Omega_\odot r\sin\theta}{v_r}$$

Step 4. With $B_r = B_0(R_\odot/r)^2$ from flux conservation and constant $v_r$:

$$B_\phi = -B_0\frac{R_\odot^2}{r^2}\cdot\frac{\Omega_\odot r\sin\theta}{v_r} = -B_0\frac{R_\odot^2\Omega_\odot\sin\theta}{rv_r}$$

Step 5. The spiral angle:

$$\boxed{\tan\psi = \frac{|B_\phi|}{B_r} = \frac{\Omega_\odot r\sin\theta}{v_{\text{sw}}}}$$

Step 6. The field line equation in the equatorial plane ($\theta = \pi/2$):

$$\frac{d\phi}{dr} = \frac{B_\phi}{rB_r} = -\frac{\Omega_\odot}{v_{\text{sw}}}$$

$$\boxed{r - r_0 = -\frac{v_{\text{sw}}}{\Omega_\odot}(\phi - \phi_0) \qquad\text{(Archimedean spiral)}}$$

At Earth ($r = 1$ AU, $v_{\text{sw}} = 400$ km/s):$\tan\psi = 2.87\times10^{-6}\times1.496\times10^{11}/4\times10^5 \approx 1.07$, giving $\psi \approx 47^\circ$. The total field strength at 1 AU:$B \approx 5\text{--}10$ nT, confirmed by in-situ measurements since Mariner 2 (1962).

16.7 Magnetopause Standoff Distance: Full Derivation

$R_{mp} = R_E(B_E^2/(2\mu_0\rho v^2))^{1/6}$

Step 1. At the sub-solar magnetopause, the solar wind is brought to rest (or diverted). The total external pressure:

$$P_{\text{ext}} = \frac{1}{2}\rho_{\text{sw}}v_{\text{sw}}^2 + P_{\text{th,sw}} + \frac{B_{\text{sw}}^2}{2\mu_0}$$

The thermal and magnetic pressures are typically small compared to the dynamic pressure, so:

$$P_{\text{ext}} \approx k_p\,\rho_{\text{sw}}\,v_{\text{sw}}^2$$

where $k_p \approx 0.88$ accounts for the bow shock.

Step 2. The internal magnetic pressure from the dipole field at distance $r$:

$$P_{\text{int}} = \frac{B^2(r)}{2\mu_0} = \frac{1}{2\mu_0}\left(\frac{B_E R_E^3}{r^3}\right)^2 = \frac{B_E^2}{2\mu_0}\left(\frac{R_E}{r}\right)^6$$

Step 3. Setting $P_{\text{ext}} = P_{\text{int}}$ and solving for $r$:

$$\boxed{R_{mp} = R_E\left(\frac{B_E^2}{2\mu_0\,k_p\,\rho_{\text{sw}}\,v_{\text{sw}}^2}\right)^{1/6}}$$

Step 4. Plugging in typical values ($n = 5$ cm$^{-3}$, $v = 400$ km/s):

$$\rho v^2 = 5\times10^6\times1.67\times10^{-27}\times(4\times10^5)^2 = 1.3\times10^{-9}\text{ Pa}$$

$$\frac{B_E^2}{2\mu_0} = \frac{(3.1\times10^{-5})^2}{2\times4\pi\times10^{-7}} = 0.38\text{ Pa}$$

$$R_{mp} = R_E\left(\frac{0.38}{1.3\times10^{-9}}\right)^{1/6} = R_E\times(2.9\times10^8)^{1/6} \approx 10.3\,R_E$$

The 1/6 power dependence means the magnetopause is rather insensitive to solar wind conditions. Even a factor of 10 increase in dynamic pressure only moves the boundary inward by a factor of$10^{1/6} \approx 1.47$. However, during extreme events ($n=50$ cm$^{-3}$,$v=1000$ km/s), $R_{mp}$ can shrink to $\sim 5 R_E$, exposing geosynchronous orbit ($6.6 R_E$) to the magnetosheath.

16.8 Dungey Cycle: Reconnection-Driven Convection

From Reconnection Rate to Substorm and Ring Current Injection

Step 1. The dayside reconnection rate determines the cross-polar cap potential:

$$\Phi_{pc} = E_{sw}\cdot\ell_{\text{eff}} = v_{\text{sw}}\,|B_z^{\text{south}}|\,\ell_{\text{eff}}$$

where $\ell_{\text{eff}} \sim 5\text{--}10\,R_E$ is the effective merging length. For $v_{\text{sw}} = 400$ km/s, $B_z = -10$ nT, $\ell = 7\,R_E$:

$$\Phi_{pc} = 4\times10^5\times10^{-8}\times7\times6.4\times10^6 \approx 180\text{ kV}$$

Step 2. This potential drives magnetospheric convection with an electric field$E = \Phi_{pc}/(2L_{pc})$ where $L_{pc}$ is the polar cap diameter. Flux tubes convect earthward on the nightside.

Step 3. The convection electric field maps to the equatorial plane, driving plasma inward and energizing it adiabatically. The energy injection rate into the ring current:

$$\boxed{\frac{dE_{RC}}{dt} \approx \epsilon\,v_{\text{sw}}\,B_z^2\,\ell^2/\mu_0 \sim 10^{12}\text{--}10^{13}\text{ W (during storms)}}$$

where $\epsilon \sim 0.1\text{--}0.2$ is the reconnection efficiency.

The Dungey cycle connects three key phenomena: (1) dayside reconnection opens magnetic flux, (2) the open flux is convected to the tail, loading energy (growth phase), and (3) tail reconnection releases the stored energy as substorms (expansion phase) and injects particles into the ring current. Sustained southward IMF keeps the cycle running, building up the ring current (storm main phase).

16.9 Parker Spiral: View from Above the Ecliptic

The Archimedean spiral pattern of the interplanetary magnetic field as seen from above the ecliptic plane.

Parker spiral magnetic field lines as seen from above the ecliptic. Red lines: field pointing away from the Sun. Blue lines: field pointing toward the Sun. The heliospheric current sheet separates the two sectors. At Earth (1 AU), the garden-hose angle is approximately 45 degrees for 400 km/s solar wind.

Extended Simulation: Parker Spiral & Magnetopause Standoff

Extended: Parker Spiral, IMF Components, Magnetopause Map, Dungey Potential

Python

script.py141 lines

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

fig = plt.figure(figsize=(14, 10))
fig.suptitle('Parker Spiral and Magnetopause Dynamics', fontsize=16, fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

# --- Panel 1: Parker spiral (polar plot) ---
ax1 = fig.add_subplot(221, polar=True)
ax1.set_facecolor('#0a0a0a')
ax1.tick_params(colors='white')
ax1.set_title('Parker Spiral (Equatorial Plane)', color='white', pad=15)

Omega = 2.87e-6  # rad/s
R_sun_AU = 1.0 / 215.0

for v_sw, col, lbl in [(300, '#f43f5e', '300 km/s (slow)'),
                         (400, '#fb923c', '400 km/s'),
                         (600, '#fbbf24', '600 km/s'),
                         (800, '#4ade80', '800 km/s (fast)')]:
    v_AU_per_s = v_sw * 1e3 / 1.496e11
    r = np.linspace(R_sun_AU, 3.0, 1000)
    phi = -Omega * (r - R_sun_AU) / v_AU_per_s

ax1.plot(phi, r, color=col, linewidth=1.5, label=lbl)
    # Additional spirals at different starting longitudes
    for dphi in [np.pi/2, np.pi, 3*np.pi/2]:
        ax1.plot(phi + dphi, r, color=col, linewidth=0.5, alpha=0.3)

ax1.set_rmax(3.0)
ax1.plot(0, 1.0, 'o', color='#60a5fa', markersize=6)  # Earth
ax1.plot(0, 0, 'o', color='#fbbf24', markersize=8)  # Sun
ax1.legend(fontsize=7, loc='upper right', facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: B components vs distance ---
ax2 = fig.add_subplot(222)
ax2.set_facecolor('#0a0a0a')
ax2.tick_params(colors='white')

r_AU = np.linspace(0.1, 5, 500)
B0_r = 3.5  # nT at 1 AU for Br
v_sw = 400  # km/s

Br = B0_r / r_AU**2  # radial component
Bphi = B0_r * Omega * r_AU * 1.496e11 / (v_sw * 1e3 * r_AU**2)
Btotal = np.sqrt(Br**2 + Bphi**2)
psi = np.degrees(np.arctan2(np.abs(Bphi), np.abs(Br)))

ax2.semilogy(r_AU, Br, color='#f43f5e', linewidth=2, label='|B_r| ~ r^{-2}')
ax2.semilogy(r_AU, Bphi, color='#60a5fa', linewidth=2, label='|B_phi| ~ r^{-1}')
ax2.semilogy(r_AU, Btotal, color='#fbbf24', linewidth=2.5, linestyle='--', label='|B| total')

ax2.axvline(x=1, color='white', linestyle=':', alpha=0.3)
ax2.text(1.05, 20, 'Earth', fontsize=9, color='white')
ax2.set_xlabel('Distance [AU]', color='white')
ax2.set_ylabel('|B| [nT]', color='white')
ax2.set_title('IMF Components vs Distance', color='white')
ax2.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')
for spine in ax2.spines.values():
    spine.set_color('#333')

# --- Panel 3: Magnetopause vs conditions (2D heatmap) ---
ax3 = fig.add_subplot(223)
ax3.set_facecolor('#0a0a0a')
ax3.tick_params(colors='white')

n_range = np.linspace(1, 50, 100)  # cm^-3
v_range_mp = np.linspace(300, 1200, 100)  # km/s
N, V = np.meshgrid(n_range, v_range_mp)

m_p = 1.673e-27
B_E = 3.1e-5
R_E = 6.371e6
mu_0 = 4 * np.pi * 1e-7

rho = N * 1e6 * m_p
P_dyn = 0.88 * rho * (V * 1e3)**2
R_mp = R_E * (B_E**2 / (2 * mu_0 * P_dyn))**(1.0/6.0) / R_E

im3 = ax3.pcolormesh(N, V, R_mp, cmap='RdYlGn', shading='auto', vmin=4, vmax=14)
plt.colorbar(im3, ax=ax3, label='R_mp [R_E]')
ax3.contour(N, V, R_mp, levels=[6.6], colors=['white'], linewidths=[2])
ax3.text(35, 900, 'R_mp = 6.6 R_E\n(GEO exposed)', fontsize=9, color='white', ha='center')

ax3.set_xlabel('Solar wind density [cm^-3]', color='white')
ax3.set_ylabel('Solar wind speed [km/s]', color='white')
ax3.set_title('Magnetopause Standoff Distance', color='white')
for spine in ax3.spines.values():
    spine.set_color('#333')

# --- Panel 4: Dungey cycle potential and energy injection ---
ax4 = fig.add_subplot(224)
ax4.set_facecolor('#0a0a0a')
ax4.tick_params(colors='white')

Bz_range = np.linspace(-30, 5, 200)  # nT
v_vals = [350, 500, 700]
colors_v = ['#fbbf24', '#fb923c', '#f43f5e']

for v, col in zip(v_vals, colors_v):
    # Reconnection potential (Milan et al.)
    Phi = np.where(Bz_range < 0,
                   0.1 * v * 1e3 * np.abs(Bz_range) * 1e-9 * 7 * 6.371e6 / 1e3,  # kV
                   10)  # quiet baseline kV
    ax4.plot(Bz_range, Phi, color=col, linewidth=2.5,
            label='v_sw = ' + str(v) + ' km/s')

ax4.axvline(x=0, color='white', linestyle=':', alpha=0.3)
ax4.set_xlabel('IMF Bz [nT]', color='white')
ax4.set_ylabel('Cross-polar cap potential [kV]', color='white')
ax4.set_title('Reconnection Rate vs Southward IMF', color='white')
ax4.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')
for spine in ax4.spines.values():
    spine.set_color('#333')

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: PARKER SPIRAL AND MAGNETOPAUSE")
print("=" * 60)
print("\n--- Parker Spiral ---")
print("  tan(psi) = Omega * r * sin(theta) / v_sw")
print("  At 1 AU (400 km/s): psi ~ 47 deg")
print("  B_r ~ r^{-2}, B_phi ~ r^{-1}")
print("  B transitions from radial to azimuthal at r ~ v_sw/Omega ~ 1 AU")
print("\n--- Magnetopause ---")
print("  R_mp = R_E * (B_E^2 / (2*mu_0*rho*v^2))^{1/6}")
print("  Typical: 10.3 R_E")
print("  Storm: 5-7 R_E")
print("  GEO orbit at 6.6 R_E can be exposed")
print("\n--- Dungey Cycle ---")
print("  Phi_pc ~ v_sw * Bz * l_eff ~ 50-200 kV (storms)")
print("  Cycle time: 1-4 hours")
print("  Energy injection: 10^12 - 10^13 W")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

16.10 Enhanced Parker Spiral: Numerical Examples and Distance Dependence

Spiral Angle at Multiple Heliospheric Distances

The Parker spiral angle $\psi$ depends on both heliocentric distance $r$ and solar wind speed $v_{\text{sw}}$. Recall:

$$\tan\psi(r) = \frac{\Omega_\odot\,r}{v_{\text{sw}}}$$

Numerical Example 1: Earth (1 AU). With $\Omega_\odot = 2.87 \times 10^{-6}$ rad/s,$r = 1.496 \times 10^{11}$ m, and $v_{\text{sw}} = 400$ km/s:

$$\tan\psi = \frac{2.87 \times 10^{-6} \times 1.496 \times 10^{11}}{4.0 \times 10^5} = 1.073 \quad\Longrightarrow\quad \psi \approx 47^\circ$$

Numerical Example 2: Mars (1.524 AU). At Mars orbit:

$$\tan\psi = \frac{2.87 \times 10^{-6} \times 1.524 \times 1.496 \times 10^{11}}{4.0 \times 10^5} = 1.636 \quad\Longrightarrow\quad \psi \approx 59^\circ$$

Numerical Example 3: Jupiter (5.2 AU). At Jupiter:

$$\tan\psi = \frac{2.87 \times 10^{-6} \times 5.2 \times 1.496 \times 10^{11}}{4.0 \times 10^5} = 5.58 \quad\Longrightarrow\quad \psi \approx 80^\circ$$

The spiral tightens dramatically with distance. At large $r$, $\psi \to 90^\circ$ and the field becomes nearly azimuthal. This is because $B_r \propto r^{-2}$ while $B_\phi \propto r^{-1}$:

$$\frac{|B_\phi|}{B_r} = \frac{\Omega_\odot r}{v_{\text{sw}}} \propto r \quad\Longrightarrow\quad B_\phi \text{ dominates at large } r$$

Speed Dependence at Fixed Distance

At 1 AU, the spiral angle varies strongly with solar wind speed:

Solar wind speed

Spiral angle at 1 AU

300 km/s (very slow)

$\psi \approx 55^\circ$

400 km/s (slow wind)

$\psi \approx 47^\circ$

500 km/s (average)

$\psi \approx 40^\circ$

600 km/s (fast wind)

$\psi \approx 35^\circ$

800 km/s (very fast)

$\psi \approx 28^\circ$

Radial distance where $\psi = 45^\circ$: Setting $\tan\psi = 1$:

$$\boxed{r_{45} = \frac{v_{\text{sw}}}{\Omega_\odot} = \frac{v_{\text{sw}}}{2.87 \times 10^{-6}} \approx 0.93 \text{ AU for } v_{\text{sw}} = 400 \text{ km/s}}$$

This critical distance $r_{45}$ marks where the field transitions from predominantly radial ($r \ll r_{45}$) to predominantly azimuthal ($r \gg r_{45}$). For fast wind (800 km/s),$r_{45} \approx 1.86$ AU, meaning the field remains more radial out to Mars. For slow wind (300 km/s),$r_{45} \approx 0.70$ AU, inside Venus's orbit.

16.11 Magnetopause Dynamics: Chapman-Ferraro Model and Reconnection Erosion

Full Chapman-Ferraro Derivation

Step 1. The Chapman-Ferraro model treats the magnetopause as a current sheet that confines the geomagnetic dipole field. The magnetopause current density is determined by the requirement that $\mathbf{B}_{\text{internal}} \neq 0$ inside and$\mathbf{B} = 0$ outside (in the idealized case). For a tangential discontinuity:

$$\mathbf{K} = \frac{1}{\mu_0}(\hat{n} \times \mathbf{B}_{\text{inside}})$$

where $\mathbf{K}$ is the surface current density (A/m) and $\hat{n}$ is the outward normal.

Step 2. The Chapman-Ferraro image dipole method: to ensure $B_n = 0$ at the boundary, place an image dipole of equal magnitude on the sunward side. The total field at the sub-solar magnetopause is then doubled:

$$B_{\text{mp}} = 2B_{\text{dipole}}(R_{\text{mp}}) = \frac{2B_E R_E^3}{R_{\text{mp}}^3}$$

Step 3. The magnetic pressure of this confined field must balance the total external pressure. Including the factor of 2 from the image dipole:

$$\frac{(2B_{\text{dipole}})^2}{2\mu_0} = k_p\,\rho_{\text{sw}}\,v_{\text{sw}}^2$$

Step 4. Solving for $R_{\text{mp}}$:

$$\frac{4B_E^2 R_E^6}{2\mu_0 R_{\text{mp}}^6} = k_p\,\rho_{\text{sw}}\,v_{\text{sw}}^2$$

$$\boxed{R_{\text{mp}} = R_E\left(\frac{2B_E^2}{\mu_0\,k_p\,\rho_{\text{sw}}\,v_{\text{sw}}^2}\right)^{1/6}}$$

The factor of 2 in the numerator (vs. the simpler form) comes from the image dipole enhancement. The parameter $k_p \approx 0.88$ accounts for the gasdynamic bow shock.

Reconnection Erosion of the Magnetopause

Step 5. When the IMF has a southward component ($B_z < 0$), dayside reconnection erodes magnetic flux from the dayside magnetosphere. This effectively weakens the internal magnetic pressure and moves the magnetopause earthward. The Shue et al. (1998) empirical model:

$$\boxed{R_{\text{mp}} = \left(11.4 + 0.013\,B_z\right)\left(\frac{P_{\text{dyn}}}{2}\right)^{-1/6.6} \quad [R_E]}$$

where $B_z$ is in nT and $P_{\text{dyn}}$ is in nPa. The linear $B_z$ term captures the erosion effect.

Step 6. Deriving the erosion magnitude $\Delta R_{\text{mp}}$ from reconnection. The rate of open flux production:

$$\frac{d\Phi_{\text{open}}}{dt} = E_{\text{rec}}\,\ell_{\text{X}} = \alpha_R\,v_{\text{sw}}\,|B_z|\,\ell_{\text{X}}$$

where $\alpha_R \sim 0.1$ is the reconnection efficiency and $\ell_X \sim 10\,R_E$is the X-line length. Each unit of open flux removed from the dayside weakens the confinement pressure. Over a time $\Delta t$:

$$\Delta\Phi_{\text{open}} = \alpha_R\,v_{\text{sw}}\,|B_z|\,\ell_X\,\Delta t$$

Step 7. The effective dipole moment is reduced: $M_{\text{eff}} = M_E - \Delta M$where $\Delta M / M_E \sim \Delta\Phi / \Phi_{\text{dayside}}$. Since$R_{\text{mp}} \propto M^{1/3}$:

$$\boxed{\Delta R_{\text{mp}} \approx -\frac{R_{\text{mp}}}{3}\frac{\Delta\Phi_{\text{open}}}{\Phi_{\text{dayside}}} \approx -\frac{R_{\text{mp}}}{3}\frac{\alpha_R\,v_{\text{sw}}\,|B_z|\,\ell_X\,\Delta t}{B_E R_E^2 / R_{\text{mp}}}}$$

For a strong storm with $|B_z| = 20$ nT, $v_{\text{sw}} = 500$ km/s, after 1 hour of sustained reconnection:

$$\Delta R_{\text{mp}} \approx -1.5\,R_E \quad \text{(from erosion alone, on top of pressure compression)}$$

Reconnection erosion and dynamic pressure compression act together. During extreme events like the March 1991 storm, the combined effect pushed the magnetopause inside geosynchronous orbit at$6.6\,R_E$, exposing multiple spacecraft to the magnetosheath.

16.12 Dungey Cycle: Detailed Derivations

Cross-Polar Cap Potential

Step 1. The solar wind imposes a motional electric field on the magnetosphere:

$$\mathbf{E}_{\text{sw}} = -\mathbf{v}_{\text{sw}} \times \mathbf{B}_{\text{IMF}}$$

For $\mathbf{v}_{\text{sw}} = v_{\text{sw}}\hat{x}$ and $\mathbf{B} = B_z\hat{z}$ (southward,$B_z < 0$):

$$E_y = -v_{\text{sw}}\,B_z = v_{\text{sw}}\,|B_z| > 0 \quad\text{(dawn-to-dusk)}$$

Step 2. This electric field is applied across an effective reconnection width$L_{\text{eff}}$. The cross-polar cap potential:

$$\boxed{\Phi_{\text{PC}} = E_y \cdot L_{\text{eff}} = v_{\text{sw}}\,|B_z|\,L_{\text{eff}}}$$

Step 3: Typical values. With $v_{\text{sw}} = 400$ km/s,$|B_z| = 5$ nT, $L_{\text{eff}} = 5\,R_E$:

$$\Phi_{\text{PC}} = 4 \times 10^5 \times 5 \times 10^{-9} \times 5 \times 6.371 \times 10^6 = 63.7 \text{ kV}$$

For a strong storm ($v_{\text{sw}} = 600$ km/s, $|B_z| = 20$ nT,$L_{\text{eff}} = 8\,R_E$):

$$\Phi_{\text{PC}} = 6 \times 10^5 \times 20 \times 10^{-9} \times 8 \times 6.371 \times 10^6 \approx 611 \text{ kV}$$

In practice, the polar cap potential saturates at ~150-250 kV due to Region-1 field-aligned current feedback (Siscoe et al., 2002). The empirical relationship:

$$\boxed{\Phi_{\text{PC,sat}} \approx \frac{\Phi_{\text{PC,linear}}}{1 + \Phi_{\text{PC,linear}}/\Phi_{\text{max}}} \quad\text{with } \Phi_{\text{max}} \approx 200 \text{ kV}}$$

Convection Electric Field and Two-Cell Pattern

Step 4. The convection electric field in the magnetosphere is derived from$\mathbf{E} = -\mathbf{v} \times \mathbf{B}$. In the equatorial plane, the dawn-to-dusk electric field drives $\mathbf{E} \times \mathbf{B}$ drift sunward:

$$v_{\text{conv}} = \frac{E}{B} = \frac{\Phi_{\text{PC}}}{L_{\text{tail}}\,B_{\text{eq}}}$$

Step 5. Mapping to the ionosphere using the dipole mapping factor$L_{\text{iono}} = L_{\text{eq}}\sqrt{R_E/L_{\text{eq}}R_E}$:

$$E_{\text{iono}} = \frac{\Phi_{\text{PC}}}{2R_{\text{PC}}} \approx \frac{100 \text{ kV}}{2 \times 1.5 \times 10^6 \text{ m}} \approx 33 \text{ mV/m}$$

Step 6. The resulting $\mathbf{E} \times \mathbf{B}$ drift in the ionosphere creates the characteristic two-cell convection pattern:

Dawn cell: Antisunward flow over the polar cap, return flow on the dawn flank. Counter-clockwise as viewed from above in the Northern Hemisphere.
Dusk cell: Antisunward flow over the polar cap, return flow on the dusk flank. Clockwise as viewed from above in the Northern Hemisphere.
Flow speeds: ~100-1000 m/s in the ionosphere, mapping to ~10-100 km/s in the equatorial magnetosphere.
Potential pattern: Dawn side is at negative potential (electron precipitation), dusk side at positive potential (ion precipitation).

The convection velocity in the ionosphere:

$$\boxed{v_{\text{iono}} = \frac{E_{\text{iono}}}{B_{\text{iono}}} = \frac{\Phi_{\text{PC}}/(2R_{\text{PC}})}{5 \times 10^{-5}} \approx 670 \text{ m/s for } \Phi_{\text{PC}} = 100 \text{ kV}}$$

SuperDARN radar observations directly measure this convection pattern. During quiet times ($\Phi_{\text{PC}} \sim 30$ kV), the pattern is weak and contracted. During storms ($\Phi_{\text{PC}} \sim 150$ kV), it expands to lower latitudes and drives rapid auroral zone convection at over 1 km/s.

16.13 Ring Current Injection: The Burton Equation

Derivation of the Dst Dynamics

Step 1. The Dst (Disturbance Storm-Time) index measures the symmetric depression of the horizontal magnetic field at Earth's surface caused by the ring current. The Dessler-Parker-Sckopke relation connects Dst to the total ring current energy:

$$\text{Dst} = -\frac{\mu_0}{4\pi}\frac{2E_{\text{RC}}}{B_E R_E^3}$$

Step 2. The Burton equation (Burton et al., 1975) models the time evolution of the pressure-corrected Dst index ($\text{Dst}^*$):

$$\boxed{\frac{d\,\text{Dst}^*}{dt} = Q(t) - \frac{\text{Dst}^*}{\tau}}$$

where $Q(t)$ is the injection function (nT/hr) and $\tau$ is the decay time constant.

Step 3. The pressure-corrected Dst removes the contribution of magnetopause currents:

$$\text{Dst}^* = \text{Dst} - b\sqrt{P_{\text{dyn}}} + c$$

where $b \approx 7.26$ nT/$\sqrt{\text{nPa}}$ and $c \approx 11$ nT (O'Brien and McPherron, 2000).

Step 4. The injection function $Q(t)$ depends on the solar wind electric field:

$$\boxed{Q(t) = \begin{cases} -4.4(v_{\text{sw}} B_s - E_c) & \text{if } v_{\text{sw}} B_s > E_c \\ 0 & \text{otherwise} \end{cases}}$$

where $B_s = \max(0, -B_z)$ is the southward IMF component, $E_c \approx 0.49$ mV/m is the threshold for injection, and $Q$ is in nT/hr.

Step 5. The decay time $\tau$ represents charge-exchange losses of ring current ions with exospheric neutrals. O'Brien and McPherron (2000) found an energy-dependent decay:

$$\boxed{\tau = 2.40\,e^{9.74/(4.69 + v_{\text{sw}} B_s)} \quad\text{hours}}$$

For quiet conditions, $\tau \approx 8$ hours. During intense driving, $\tau$ decreases to ~3-5 hours as higher-energy ions charge-exchange faster.

Step 6: Analytical solution. For constant injection $Q_0$ starting at $t = 0$with initial $\text{Dst}^*(0) = 0$:

$$\boxed{\text{Dst}^*(t) = Q_0\,\tau\left(1 - e^{-t/\tau}\right) \xrightarrow{t \gg \tau} Q_0\,\tau}$$

For $Q_0 = -30$ nT/hr and $\tau = 8$ hr:$\text{Dst}^* \to -240$ nT (intense storm).

After injection ceases, the recovery phase follows an exponential decay:$\text{Dst}^*(t) = \text{Dst}^*_{\min}\,e^{-(t - t_{\min})/\tau}$. A two-component decay with $\tau_{\text{fast}} \approx 3$ hr and $\tau_{\text{slow}} \approx 20$ hr better fits observations, reflecting the distinction between O$^+$ (fast decay) and H$^+$ (slow decay) contributions.

16.14 The Substorm Cycle: Loading-Unloading

Growth Phase, Expansion Phase, and Recovery

Step 1: Growth Phase (Loading). During southward IMF, dayside reconnection opens closed dipolar flux and transports it to the tail lobes. The lobe magnetic flux accumulates:

$$\frac{d\Phi_{\text{lobe}}}{dt} = \Phi_{\text{day}} - \Phi_{\text{night}}$$

where $\Phi_{\text{day}}$ is the dayside reconnection rate and $\Phi_{\text{night}}$is the (initially weak) tail reconnection rate.

Step 2. The lobe magnetic pressure increases as flux accumulates:

$$P_{\text{lobe}} = \frac{B_{\text{lobe}}^2}{2\mu_0} \propto \Phi_{\text{lobe}}^2$$

This increased pressure compresses the plasma sheet, thinning it from ~5 $R_E$ to ~0.5 $R_E$over the 30-90 minute growth phase.

Step 3. The loading timescale. The total open flux in the polar cap:

$$\Phi_{\text{PC}} = B_{\text{iono}}\,\pi\,R_{\text{PC}}^2 \approx 5 \times 10^{-5} \times \pi \times (15^\circ \times 111 \text{ km})^2 \approx 0.4 \text{ GWb}$$

The dayside reconnection voltage adds open flux at rate:

$$\frac{d\Phi_{\text{PC}}}{dt} = \Phi_{\text{PC,voltage}} \approx 60 \text{ kV} = 60 \text{ kWb/s}$$

Step 4. The time to significantly load the tail (e.g., double the polar cap flux):

$$\boxed{t_{\text{growth}} \approx \frac{\Delta\Phi_{\text{PC}}}{\Phi_{\text{PC,voltage}}} \approx \frac{0.4 \times 10^9}{60 \times 10^3} \approx 6700 \text{ s} \approx 110 \text{ min}}$$

In practice, substorm onset occurs after 30-90 minutes, before the tail is fully loaded, because instabilities in the thinned current sheet trigger explosive reconnection.

Expansion Phase (Unloading)

Step 5. The expansion phase begins when a near-Earth neutral line (NENL) forms at$\sim 20\text{--}30\,R_E$. The tail reconnection rate rapidly increases:

$$\Phi_{\text{night}} \gg \Phi_{\text{day}} \quad\text{(unloading)}$$

Step 6. The energy released during a substorm can be estimated from the lobe magnetic energy:

$$E_{\text{lobe}} = \frac{B_{\text{lobe}}^2}{2\mu_0}\,V_{\text{lobe}} \approx \frac{(20 \times 10^{-9})^2}{2 \times 4\pi \times 10^{-7}} \times (20\,R_E)^2 \times (100\,R_E) \approx 3 \times 10^{15} \text{ J}$$

This energy is partitioned among:

Plasma sheet heating and earthward convection (~40%)
Plasmoid ejection tailward (~30%)
Auroral precipitation and ionospheric Joule heating (~20%)
Ring current injection (~10%)

Step 7. The expansion phase duration is set by the Alfven transit time across the tail:

$$\boxed{t_{\text{expansion}} \approx \frac{L_{\text{tail}}}{v_A} \approx \frac{20\,R_E}{500 \text{ km/s}} \approx 4.3 \text{ min}}$$

The full expansion phase lasts 15-30 minutes as reconnection proceeds through the stored flux.

Recovery Phase

Step 8. The recovery phase begins when the tail reconnection rate drops below the dayside rate. Open flux is no longer being released faster than it is created. The polar cap contracts, the current sheet thickens, and magnetospheric convection returns to a quiet state over 1-2 hours.

During a geomagnetic storm, multiple substorm cycles occur in sequence, each injecting fresh particles into the ring current. The storm main phase consists of ~5-20 substorms over 6-24 hours, progressively building the ring current until the driving IMF turns northward.

16.15 Dungey Cycle: Magnetic Topology Diagram

The complete Dungey cycle showing dayside reconnection, antisunward convection over the poles, tail reconnection, and return flow through the flanks.

Dungey Cycle Steps:

Dayside reconnection — southward IMF merges with northward geomagnetic field, opening closed flux
Antisunward convection — solar wind drags open field lines over the polar caps
Lobe accumulation — open flux piles up in the north and south tail lobes
Tail reconnection — oppositely directed lobe fields reconnect at the neutral line, closing flux
Earthward return — newly closed flux tubes convect earthward, energizing trapped plasma
Flank return — closed flux circulates back to the dayside, completing the cycle (~1-4 hours)

Enhanced Simulation: Complete Sun-Earth Connection

Complete Sun-Earth Connection: Parker Spiral, IMF, Magnetopause, Dungey Cycle, Burton Dst, Substorms

Python

script.py312 lines

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

fig, axes = plt.subplots(2, 3, figsize=(20, 12))
fig.suptitle('Sun-Earth Connection: Complete Analysis', fontsize=18, fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

for ax in axes.flat:
    ax.set_facecolor('#0a0a0a')
    ax.tick_params(colors='white', labelsize=9)
    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('#444')

# ====== Panel 1: Parker spiral for 4 wind speeds (polar) ======
ax1 = fig.add_subplot(231, polar=True)
ax1.set_facecolor('#0a0a0a')
ax1.tick_params(colors='white', labelsize=8)
ax1.set_title('Parker Spiral\n(4 Wind Speeds)', color='white', pad=15, fontsize=11)

Omega = 2.87e-6
R_sun_AU = 1.0 / 215.0

wind_speeds = [
    (300, '#f43f5e', '300 km/s'),
    (400, '#fb923c', '400 km/s'),
    (600, '#fbbf24', '600 km/s'),
    (800, '#4ade80', '800 km/s'),
]

for v_sw, col, lbl in wind_speeds:
    v_AU_per_s = v_sw * 1e3 / 1.496e11
    r = np.linspace(R_sun_AU, 3.5, 1000)
    phi = -Omega * (r - R_sun_AU) / v_AU_per_s
    ax1.plot(phi, r, color=col, linewidth=2, label=lbl)
    for dphi in [np.pi/2, np.pi, 3*np.pi/2]:
        ax1.plot(phi + dphi, r, color=col, linewidth=0.4, alpha=0.25)

ax1.plot(0, 1.0, 'o', color='#60a5fa', markersize=7)
ax1.plot(0, 1.524, 's', color='#fb923c', markersize=5)
ax1.plot(0, 0, 'o', color='#fbbf24', markersize=10)
ax1.set_rmax(3.5)
ax1.legend(fontsize=7, loc='upper right', facecolor='#0a0a0a', edgecolor='#444', labelcolor='white')

# ====== Panel 2: IMF components vs distance ======
ax2 = axes[0, 1]
r_AU = np.linspace(0.1, 5.5, 500)
B0 = 3.5  # nT at 1 AU normalization

Br = B0 / r_AU**2
Bphi = B0 * Omega * r_AU * 1.496e11 / (400e3 * r_AU**2)
Btotal = np.sqrt(Br**2 + Bphi**2)

planets = [(0.72, 'Venus', '#fbbf24'), (1.0, 'Earth', '#60a5fa'),
           (1.52, 'Mars', '#fb923c'), (5.2, 'Jupiter', '#4ade80')]
for dist, name, col in planets:
    ax2.axvline(x=dist, color=col, linestyle=':', alpha=0.3, linewidth=1)
    ypos = B0 / dist**2 * 1.5
    ax2.text(dist + 0.05, ypos, name, fontsize=8, color=col, rotation=0)

# Mark crossover point
r_cross = 400e3 / (Omega * 1.496e11)
ax2.axvline(x=r_cross, color='white', linestyle='--', alpha=0.4, linewidth=1)
ax2.text(r_cross + 0.1, 0.3, 'B_r = B_phi', fontsize=8, color='white', alpha=0.6)

ax2.set_xlabel('Heliocentric distance [AU]')
ax2.set_ylabel('Magnetic field [nT]')
ax2.set_title('IMF Components vs Distance', fontsize=11)
ax2.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#444', labelcolor='white')
ax2.set_xlim(0, 5.5)
ax2.set_ylim(0.01, 100)

# ====== Panel 3: Magnetopause standoff 2D map ======
ax3 = axes[0, 2]
n_range = np.linspace(1, 60, 120)
v_range_mp = np.linspace(250, 1200, 120)
N, V = np.meshgrid(n_range, v_range_mp)

m_p = 1.673e-27
B_E = 3.1e-5
R_E = 6.371e6
mu_0 = 4 * np.pi * 1e-7

rho = N * 1e6 * m_p
P_dyn = 0.88 * rho * (V * 1e3)**2
R_mp = (B_E**2 / (2 * mu_0 * P_dyn))**(1.0/6.0)

im3 = ax3.pcolormesh(N, V, R_mp, cmap='RdYlGn', shading='auto', vmin=4, vmax=14)
cbar = plt.colorbar(im3, ax=ax3)
cbar.set_label('R_mp [R_E]', color='white', fontsize=9)
cbar.ax.tick_params(colors='white', labelsize=8)

contours = ax3.contour(N, V, R_mp, levels=[6.6, 8.0, 10.0, 12.0],
                        colors=['white', '#fbbf24', '#fb923c', '#f43f5e'], linewidths=[2, 1.5, 1.5, 1.5])
ax3.clabel(contours, fmt='%.1f', fontsize=8, colors='white')

ax3.plot(5, 400, '*', color='white', markersize=12)
ax3.text(7, 420, 'Typical', fontsize=8, color='white')
ax3.plot(20, 800, '*', color='#f43f5e', markersize=12)
ax3.text(22, 820, 'Storm', fontsize=8, color='#f43f5e')

ax3.set_xlabel('Density [cm^-3]')
ax3.set_ylabel('Speed [km/s]')
ax3.set_title('Magnetopause Standoff [R_E]', fontsize=11)

# ====== Panel 4: Cross-polar cap potential (Dungey cycle) ======
ax4 = axes[1, 0]
Bz_range = np.linspace(-30, 5, 300)

for v_sw, col, lbl in [(350, '#fbbf24', '350 km/s'), (500, '#fb923c', '500 km/s'), (700, '#f43f5e', '700 km/s')]:
    # Linear potential
    Phi_linear = np.where(Bz_range < 0,
                          v_sw * 1e3 * np.abs(Bz_range) * 1e-9 * 7 * R_E / 1e3,
                          10.0)
    # Saturated potential
    Phi_max = 200.0  # kV
    Phi_sat = Phi_linear / (1.0 + Phi_linear / Phi_max)
    ax4.plot(Bz_range, Phi_sat, color=col, linewidth=2.5, label='v_sw = ' + str(v_sw) + ' km/s')
    ax4.plot(Bz_range, Phi_linear, color=col, linewidth=1, linestyle=':', alpha=0.4)

ax4.axhline(y=200, color='white', linestyle='--', alpha=0.3, linewidth=1)
ax4.text(-28, 205, 'Saturation limit', fontsize=8, color='white', alpha=0.5)
ax4.axhline(y=50, color='white', linestyle=':', alpha=0.2)
ax4.text(-28, 55, 'Quiet-moderate', fontsize=8, color='white', alpha=0.4)
ax4.axvline(x=0, color='white', linestyle=':', alpha=0.3)

ax4.set_xlabel('IMF Bz [nT]')
ax4.set_ylabel('Cross-polar cap potential [kV]')
ax4.set_title('Dungey Cycle: Polar Cap Potential', fontsize=11)
ax4.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#444', labelcolor='white')
ax4.set_ylim(0, 300)

# ====== Panel 5: Burton equation Dst evolution ======
ax5 = axes[1, 1]

# Simulate a geomagnetic storm
hours = np.linspace(0, 72, 2000)
dt_hr = hours[1] - hours[0]
Dst_star = np.zeros_like(hours)

# Storm driver: southward IMF from t=6 to t=18 hours
vBs = np.zeros_like(hours)
# SSC at t=6, main phase 6-18h, recovery after
for i, t in enumerate(hours):
    if 6 <= t < 8:
        vBs[i] = 2.0  # mV/m, moderate onset
    elif 8 <= t < 12:
        vBs[i] = 5.0  # mV/m, intense main phase
    elif 12 <= t < 15:
        vBs[i] = 8.0  # mV/m, peak driving
    elif 15 <= t < 18:
        vBs[i] = 3.0  # mV/m, waning
    else:
        vBs[i] = 0.0

Ec = 0.49  # mV/m threshold
for i in range(1, len(hours)):
    if vBs[i] > Ec:
        Q = -4.4 * (vBs[i] - Ec)  # nT/hr
        tau = 2.40 * np.exp(9.74 / (4.69 + vBs[i]))
    else:
        Q = 0.0
        tau = 8.0
    Dst_star[i] = Dst_star[i-1] + (Q - Dst_star[i-1] / tau) * dt_hr

ax5.plot(hours, Dst_star, color='#f43f5e', linewidth=2.5, label='Dst* (Burton eq.)')

# Show driving on twin axis
ax5b = ax5.twinx()
ax5b.fill_between(hours, 0, vBs, color='#fb923c', alpha=0.2)
ax5b.plot(hours, vBs, color='#fb923c', linewidth=1.5, linestyle='--', label='vBs driver')
ax5b.set_ylabel('vBs [mV/m]', color='#fb923c', fontsize=9)
ax5b.tick_params(axis='y', colors='#fb923c', labelsize=8)
ax5b.set_ylim(0, 12)
ax5b.spines['right'].set_color('#fb923c')

ax5.axhline(y=0, color='white', linestyle=':', alpha=0.2)
ax5.axhline(y=-50, color='white', linestyle=':', alpha=0.15)
ax5.text(50, -45, 'Moderate storm', fontsize=8, color='white', alpha=0.4)
ax5.axhline(y=-100, color='white', linestyle=':', alpha=0.15)
ax5.text(50, -95, 'Intense storm', fontsize=8, color='white', alpha=0.4)

# Mark phases
ax5.axvspan(6, 8, alpha=0.05, color='#fbbf24')
ax5.axvspan(8, 18, alpha=0.08, color='#f43f5e')
ax5.axvspan(18, 72, alpha=0.03, color='#4ade80')
ax5.text(7, 10, 'SSC', fontsize=8, color='#fbbf24', ha='center')
ax5.text(13, 10, 'Main phase', fontsize=8, color='#f43f5e', ha='center')
ax5.text(40, 10, 'Recovery', fontsize=8, color='#4ade80', ha='center')

ax5.set_xlabel('Time [hours]')
ax5.set_ylabel('Dst* [nT]')
ax5.set_title('Burton Equation: Storm Dst Evolution', fontsize=11)
ax5.legend(fontsize=8, loc='lower left', facecolor='#0a0a0a', edgecolor='#444', labelcolor='white')
ax5.set_xlim(0, 72)

# ====== Panel 6: Substorm loading-unloading cycle ======
ax6 = axes[1, 2]

t_sub = np.linspace(0, 6, 1000)  # hours
# Simulate 3 substorm cycles
lobe_flux = np.zeros_like(t_sub)
al_index = np.zeros_like(t_sub)

flux_val = 0.0
for i, t in enumerate(t_sub):
    # Dayside reconnection adds flux continuously
    Phi_day = 60.0  # kWb/s = kV

# Substorm onsets at t=1.2, 2.8, 4.5 hours
    onsets = [1.2, 2.8, 4.5]
    Phi_night = 10.0  # quiet tail reconnection
    for onset in onsets:
        if onset <= t < onset + 0.3:
            # Expansion phase: rapid unloading
            Phi_night = 300.0 * np.exp(-(t - onset) / 0.15)
        elif onset + 0.3 <= t < onset + 1.0:
            # Recovery phase
            Phi_night = 30.0 * np.exp(-(t - onset - 0.3) / 0.5)

dt_s = (t_sub[1] - t_sub[0]) * 3600.0
    flux_val = flux_val + (Phi_day - Phi_night) * dt_s * 1e-3  # convert kWb to MWb-ish units
    if flux_val < 0:
        flux_val = 0.0
    lobe_flux[i] = flux_val

# AL index proxy
    for onset in onsets:
        if onset <= t < onset + 0.1:
            al_index[i] = -800 * (t - onset) / 0.1
        elif onset + 0.1 <= t < onset + 0.4:
            al_index[i] = -800 * np.exp(-(t - onset - 0.1) / 0.15)
        elif onset + 0.4 <= t < onset + 1.0:
            al_index[i] = min(al_index[i], -100 * np.exp(-(t - onset - 0.4) / 0.3))

ax6.plot(t_sub, lobe_flux / max(lobe_flux) * 100, color='#fb923c', linewidth=2.5, label='Lobe flux (normalized)')

ax6b = ax6.twinx()
ax6b.plot(t_sub, al_index, color='#f43f5e', linewidth=2, label='AL index proxy')
ax6b.set_ylabel('AL index [nT]', color='#f43f5e', fontsize=9)
ax6b.tick_params(axis='y', colors='#f43f5e', labelsize=8)
ax6b.spines['right'].set_color('#f43f5e')

# Mark substorm phases
for onset in [1.2, 2.8, 4.5]:
    ax6.axvline(x=onset, color='#f43f5e', linestyle='--', alpha=0.4, linewidth=1)
    ax6.text(onset - 0.5, 95, 'Growth', fontsize=7, color='#fb923c', ha='center')
    ax6.text(onset + 0.15, 95, 'Exp', fontsize=7, color='#f43f5e', ha='center')
    ax6.text(onset + 0.6, 95, 'Rec', fontsize=7, color='#4ade80', ha='center')

ax6.set_xlabel('Time [hours]')
ax6.set_ylabel('Lobe flux [% of peak]', color='#fb923c')
ax6.tick_params(axis='y', colors='#fb923c', labelsize=8)
ax6.set_title('Substorm Loading-Unloading Cycle', fontsize=11)
ax6.legend(fontsize=8, loc='upper left', facecolor='#0a0a0a', edgecolor='#444', labelcolor='white')
ax6b.legend(fontsize=8, loc='lower left', facecolor='#0a0a0a', edgecolor='#444', 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" + "=" * 70)
print("COMPLETE SUN-EARTH CONNECTION ANALYSIS")
print("=" * 70)

print("\n--- Enhanced Parker Spiral ---")
for v in [300, 400, 600, 800]:
    psi_earth = np.degrees(np.arctan(Omega * 1.496e11 / (v * 1e3)))
    psi_mars = np.degrees(np.arctan(Omega * 1.524 * 1.496e11 / (v * 1e3)))
    print("  v=" + str(v) + " km/s: psi(Earth)=" + str(round(psi_earth, 1)) + " deg, psi(Mars)=" + str(round(psi_mars, 1)) + " deg")

print("\n--- IMF Crossover ---")
r_cross_AU = 400e3 / (Omega * 1.496e11)
print("  B_r = B_phi at r = " + str(round(r_cross_AU, 2)) + " AU (for 400 km/s)")

print("\n--- Magnetopause (Chapman-Ferraro) ---")
n_vals = [5, 10, 20, 50]
v_vals = [400, 600, 800]
for n in n_vals:
    for v in v_vals:
        rho = n * 1e6 * m_p
        P = 0.88 * rho * (v * 1e3)**2
        Rmp = (B_E**2 / (2 * mu_0 * P))**(1.0/6.0) / R_E
        print("  n=" + str(n) + ", v=" + str(v) + ": R_mp = " + str(round(Rmp, 1)) + " R_E")

print("\n--- Burton Equation Storm ---")
print("  Peak Dst* = " + str(round(min(Dst_star), 0)) + " nT")
peak_idx = np.argmin(Dst_star)
print("  Time of minimum: " + str(round(hours[peak_idx], 1)) + " hours")
# Recovery: find when Dst recovers to -50
recovery_mask = (hours > hours[peak_idx]) & (Dst_star > -50)
if np.any(recovery_mask):
    t_rec = hours[recovery_mask][0]
    print("  Recovery to -50 nT: " + str(round(t_rec, 1)) + " hours")

print("\n--- Substorm Timescales ---")
print("  Growth phase: 30-90 min (typical 60 min)")
print("  Expansion phase: 15-30 min")
print("  Recovery phase: 30-60 min")
print("  Full cycle: ~2-3 hours")
print("  Energy per substorm: ~3e15 J")
print("=" * 70)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Space Weather Back to Course Overview →

Sun-Earth Connection

16.1 The Parker Spiral

Derivation 1: Spiral Angle of the Interplanetary Magnetic Field

16.2 Heliospheric Current Sheet and Sector Structure

Derivation 2: Magnetic Sector Boundaries

16.3 Magnetopause Standoff Distance

Derivation 3: Pressure Balance at the Magnetopause

16.4 The Dungey Cycle

Derivation 4: Magnetic Reconnection-Driven Convection

16.5 Ring Current Injection

Derivation 5: Adiabatic Particle Energization

Numerical Simulation

Sun-Earth Connection: Parker Spiral, Spiral Angle, Magnetopause, Adiabatic Energization

16.6 Full Parker Spiral Derivation

Step-by-Step from Frozen-in Field in the Rotating Frame

16.7 Magnetopause Standoff Distance: Full Derivation

\(R_{mp} = R_E(B_E^2/(2\mu_0\rho v^2))^{1/6}\)

16.8 Dungey Cycle: Reconnection-Driven Convection

From Reconnection Rate to Substorm and Ring Current Injection

16.9 Parker Spiral: View from Above the Ecliptic

Extended Simulation: Parker Spiral & Magnetopause Standoff

Extended: Parker Spiral, IMF Components, Magnetopause Map, Dungey Potential

16.10 Enhanced Parker Spiral: Numerical Examples and Distance Dependence

Spiral Angle at Multiple Heliospheric Distances

Speed Dependence at Fixed Distance

16.11 Magnetopause Dynamics: Chapman-Ferraro Model and Reconnection Erosion

Full Chapman-Ferraro Derivation

Reconnection Erosion of the Magnetopause

16.12 Dungey Cycle: Detailed Derivations

Cross-Polar Cap Potential

Convection Electric Field and Two-Cell Pattern

16.13 Ring Current Injection: The Burton Equation

Derivation of the Dst Dynamics

16.14 The Substorm Cycle: Loading-Unloading

Growth Phase, Expansion Phase, and Recovery

Expansion Phase (Unloading)

Recovery Phase

16.15 Dungey Cycle: Magnetic Topology Diagram

Enhanced Simulation: Complete Sun-Earth Connection

Complete Sun-Earth Connection: Parker Spiral, IMF, Magnetopause, Dungey Cycle, Burton Dst, Substorms