Part II: Solar Atmosphere | Chapter 8

Solar Wind Physics

Parker's isothermal wind, the critical point, fast and slow wind, and MHD wind models

8.1 Parker's Isothermal Solar Wind

Derivation 1: The Parker Wind Equation (Full Derivation)

Eugene Parker (1958) showed that a static corona is impossible—it must expand as a supersonic wind. This was one of the most consequential theoretical predictions in space physics, confirmed just four years later by Mariner 2. The derivation is a classic of mathematical physics.

Historical Context

Before Parker, Chapman (1957) had shown that the corona has extremely high thermal conductivity, producing temperatures of $\sim 10^5$ K even at Earth's orbit. Biermann (1951) had inferred the existence of a continuous solar corpuscular radiation from observations of comet tail deflections. Parker unified these ideas by showing that the hydrodynamic equations demand a supersonic outflow—the solar wind.

Step 1: Governing Equations. We start with the steady-state, spherically symmetric, hydrodynamic equations for a fully ionized hydrogen plasma. The continuity equation ensures mass conservation:

$$\frac{1}{r^2}\frac{d}{dr}(\rho v r^2) = 0 \quad \Longrightarrow \quad \rho v r^2 = \text{const} = \frac{\dot{M}}{4\pi}$$

where $\dot{M}$ is the mass loss rate and $v(r)$ is the radial velocity. This tells us that $\rho = \dot{M} / (4\pi r^2 v)$—if the velocity increases outward, the density must decrease faster than $r^{-2}$.

Step 2: Momentum Equation. The radial momentum equation (Euler equation with gravity):

$$\rho v \frac{dv}{dr} = -\frac{dP}{dr} - \frac{GM_\odot \rho}{r^2}$$

Step 3: Isothermal Assumption. For an isothermal corona at temperature $T$, the equation of state is $P = \rho c_s^2$ where the isothermal sound speed is$c_s = \sqrt{k_BT / (\mu m_H)}$. Taking the derivative:

$$\frac{dP}{dr} = c_s^2 \frac{d\rho}{dr}$$

Step 4: Eliminating Density. From the continuity equation, we can express$d\rho/dr$ in terms of $dv/dr$. Differentiating $\rho v r^2 = \text{const}$:

$$\frac{1}{\rho}\frac{d\rho}{dr} = -\frac{1}{v}\frac{dv}{dr} - \frac{2}{r}$$

Step 5: Combining. Substituting $dP/dr = c_s^2 d\rho/dr$ into the momentum equation and using the density derivative from Step 4:

$$\rho v \frac{dv}{dr} = -c_s^2 \rho\left(-\frac{1}{v}\frac{dv}{dr} - \frac{2}{r}\right) - \frac{GM_\odot\rho}{r^2}$$

Expanding and rearranging:

$$v\frac{dv}{dr} + \frac{c_s^2}{v}\frac{dv}{dr} = \frac{2c_s^2}{r} - \frac{GM_\odot}{r^2}$$

Step 6: The Parker Equation. Factoring out $dv/dr$:

$$\boxed{\left(v - \frac{c_s^2}{v}\right)\frac{dv}{dr} = \frac{2c_s^2}{r} - \frac{GM_\odot}{r^2}}$$

The Parker solar wind equation

This is a first-order ODE with a critical point where both sides vanish simultaneously.

Step 7: Critical Point Analysis. The right-hand side vanishes at the critical radius:

$$r_c = \frac{GM_\odot}{2c_s^2}$$

At this radius, the left-hand side must also vanish. The factor $(v - c_s^2/v)$equals zero when $v = c_s$. Therefore, the wind passes through the sound speed at the critical radius. Numerically:

$$r_c = \frac{6.674 \times 10^{-11} \times 1.989 \times 10^{30}}{2 \times (1.66 \times 10^5)^2} \approx 5.8 R_\odot \quad \text{(for } T = 10^6 \text{ K)}$$

Step 8: Topology of Solutions. The Parker equation can be integrated analytically. Defining $u = v/c_s$ and $x = r/r_c$:

$$u^2 - \ln u^2 = 4\ln x + \frac{4}{x} + C$$

where $C$ is an integration constant. There are four classes of solutions:

• Class I (Parker wind): Subsonic at $r < r_c$, passes through sonic point, supersonic for $r > r_c$. $C = -3$.
• Class II (accretion): Supersonic infall that decelerates through the sonic point.
• Class III (breeze): Always subsonic, velocity increases then decreases. Not physical (finite pressure at infinity).
• Class IV: Always supersonic. Not physical (requires supersonic base).

Only the Class I (Parker wind) solution satisfies both boundary conditions: subsonic near the Sun and vanishing pressure at infinity. This transonic solution was Parker's key insight. The predicted velocity at 1 AU is $v \approx 300\text{--}400$ km/s for$T \approx 1\text{--}2$ MK, in excellent agreement with observations.

8.2 The Critical Point and Mass Loss Rate

Derivation 2: Mass Loss Rate from the Parker Solution

Step 1. At the critical point, $v(r_c) = c_s$ and$\rho(r_c)$ is determined by the coronal base conditions. The density at $r_c$ is related to the base density by hydrostatic equilibrium below the critical point:

$$\rho(r_c) \approx \rho_0 \exp\left(-\frac{GM_\odot}{c_s^2}\left(\frac{1}{R_\odot} - \frac{1}{r_c}\right)\right)$$

Step 2. Since $r_c = GM_\odot / (2c_s^2)$, the exponent simplifies to:

$$\frac{GM_\odot}{c_s^2 R_\odot}\left(1 - \frac{R_\odot}{r_c}\right) = \frac{2r_c}{R_\odot}\left(1 - \frac{R_\odot}{r_c}\right) = 2\left(\frac{r_c}{R_\odot} - 1\right)$$

Step 3. The mass loss rate is:

$$\boxed{\dot{M} = 4\pi r_c^2 \rho(r_c) c_s \approx 4\pi r_c^2 \rho_0 c_s \exp\left(-2\frac{r_c}{R_\odot} + 2\right)}$$

The observed solar mass loss rate is $\dot{M}_\odot \approx 2 \times 10^{-14} M_\odot$/yr = $1.3 \times 10^9$ kg/s. This is negligible for stellar evolution but drives the heliosphere and space weather. The exponential sensitivity to $r_c/R_\odot$(hence to $T$) explains why the wind properties are so sensitive to coronal conditions.

8.3 Fast and Slow Solar Wind

Derivation 3: Energy Equation and Terminal Velocity

The isothermal Parker wind predicts $v \sim 400$ km/s, but observations show two distinct wind components.

Step 1. The Bernoulli integral for the Parker wind:

$$\frac{v^2}{2} + c_s^2 \ln\rho - \frac{GM_\odot}{r} = \text{const}$$

Step 2. At large distances ($r \to \infty$), the gravitational and logarithmic terms become negligible, and the terminal velocity is:

$$\boxed{v_\infty^2 \approx 4c_s^2 \ln\left(\frac{r_c}{R_\odot}\right) - 2\frac{GM_\odot}{R_\odot} + v_0^2}$$

Step 3. Observational characteristics:

Fast Wind

• Speed: 600-800 km/s
• Source: Coronal holes
• Density: ~3 cm$^{-3}$ at 1 AU
• $T_p > T_e$ (ion-heated)
• Low variability

Slow Wind

• Speed: 300-400 km/s
• Source: Streamer belt / boundaries
• Density: ~10 cm$^{-3}$ at 1 AU
• $T_p \approx T_e$
• Highly variable, compositionally distinct

8.4 The Alfven Radius

Derivation 4: Alfven Radius and Angular Momentum Loss

Step 1. The Alfven radius $r_A$ is where the wind speed equals the Alfven speed: $v(r_A) = v_A(r_A)$. Beyond this point, the wind can no longer communicate back to the Sun via Alfven waves.

Step 2. For a radial magnetic field $B_r \propto r^{-2}$ and$\rho v r^2 = \text{const}$:

$$v_A(r) = \frac{B_r(r)}{\sqrt{\mu_0 \rho(r)}} = \frac{B_0 R_\odot^2}{r^2 \sqrt{\mu_0 \rho_0 v_0 R_\odot^2 / (v r^2)}} = B_0 R_\odot^2 \sqrt{\frac{v}{\mu_0 \dot{M}/(4\pi)}} \frac{1}{r^2}$$

Step 3. The angular momentum loss rate per unit mass from the Weber-Davis model:

$$\boxed{\dot{J} = \dot{M} \Omega r_A^2}$$

The wind removes angular momentum as if it were co-rotating with the Sun out to $r_A$. Parker Solar Probe crossed the Alfven surface at $r_A \approx 13\text{--}20 R_\odot$ in April 2021, directly measuring the sub-Alfvenic solar wind for the first time.

8.5 Weber-Davis MHD Wind Model

Derivation 5: MHD Wind with Rotation and Magnetic Field

Weber and Davis (1967) extended Parker's model to include solar rotation and the azimuthal magnetic field, producing a more complete description of the solar wind.

Step 1. The steady-state MHD equations in the equatorial plane include the induction equation for the azimuthal field:

$$B_\phi = B_r \frac{r(v_\phi - \Omega r)}{v_r - v_{Ar}} = -B_r \frac{r\Omega(1 - r^2/r_A^2)}{v_r(1 - M_A^{-2})}$$

where $M_A = v_r / v_{Ar}$ is the radial Alfven Mach number.

Step 2. Two conserved quantities along streamlines are the specific angular momentum and specific energy:

$$\mathcal{L} = r v_\phi - \frac{r B_r B_\phi}{\mu_0 \rho v_r} = \Omega r_A^2 \quad \text{(angular momentum)}$$

$$\mathcal{E} = \frac{v^2}{2} + h - \frac{GM_\odot}{r} - \frac{\Omega r B_r B_\phi}{\mu_0 \rho v_r} \quad \text{(energy)}$$

Step 3. The MHD wind has three critical points (sonic, Alfven, fast magnetosonic):

• Slow critical point: $v_r = c_{\text{slow}}$ (near $r \sim 2 R_\odot$)
• Alfven critical point: $v_r = v_A$ (at $r = r_A \sim 15 R_\odot$)
• Fast critical point: $v_r = c_{\text{fast}}$ (at $r \sim 20 R_\odot$)

Beyond $r_A$, the magnetic field becomes predominantly azimuthal ($B_\phi \gg B_r$), winding into the Parker spiral. The Weber-Davis model also predicts the torque on the Sun, giving a spin-down timescale of $\sim 10^{10}$ years for solar-type stars.

Numerical Simulation

This simulation integrates the Parker wind equation numerically, showing all four solution classes, the transonic solution, velocity and density profiles, and the critical point topology.

Parker Solar Wind: Solution Topology, Velocity Profiles, and Critical Point Analysis

Python

script.py222 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.optimize import brentq

# ============================================================
# Parker Solar Wind: Complete Solution
# ============================================================

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle("Parker Solar Wind Model", 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")

# Constants
G = 6.674e-11
M_sun = 1.989e30
R_sun = 6.957e8
k_B = 1.381e-23
m_p = 1.673e-27
AU = 1.496e11

# --- Panel 1: Parker equation solution topology ---
ax1 = axes[0, 0]

# Dimensionless form: u^2 - ln(u^2) = 4*ln(x) + 4/x + C
# where u = v/c_s, x = r/r_c
x = np.linspace(0.2, 10, 1000)

# Parker wind solution (C = -3): passes through u=1 at x=1
C_parker = -3
rhs = 4 * np.log(x) + 4 / x + C_parker

# Solve for u at each x using Lambert W or direct inversion
# For the transonic solution, we need both branches
u_supersonic = np.zeros_like(x)
u_subsonic = np.zeros_like(x)

for i, xi in enumerate(x):
    rhs_val = 4 * np.log(xi) + 4 / xi + C_parker
    # u^2 - ln(u^2) = rhs_val
    # Need to solve this numerically
    def parker_eq(u, rhs_val=rhs_val):
        if u <= 0:
            return 1e10
        return u**2 - np.log(u**2) - rhs_val

try:
        u_supersonic[i] = brentq(parker_eq, 1.0, 20.0)
    except:
        u_supersonic[i] = np.nan
    try:
        u_subsonic[i] = brentq(parker_eq, 0.001, 1.0)
    except:
        u_subsonic[i] = np.nan

# Plot transonic solution
ax1.plot(x, u_supersonic, color="#f43f5e", linewidth=2.5, label="Parker wind (transonic)")
ax1.plot(x, u_subsonic, color="#f43f5e", linewidth=2.5)

# Breeze solutions (different C values)
for C_val, col, lbl in [(-2, "#60a5fa", "Breeze (C=-2)"),
                         (-1, "#4ade80", "Breeze (C=-1)"),
                         (-4, "#fbbf24", "Double-valued (C=-4)")]:
    u_up = np.zeros_like(x)
    u_dn = np.zeros_like(x)
    for i, xi in enumerate(x):
        rhs_val = 4 * np.log(xi) + 4 / xi + C_val
        try:
            u_up[i] = brentq(lambda u: u**2 - np.log(u**2) - rhs_val, 1.001, 30.0)
        except:
            u_up[i] = np.nan
        try:
            u_dn[i] = brentq(lambda u: u**2 - np.log(u**2) - rhs_val, 0.001, 0.999)
        except:
            u_dn[i] = np.nan

ax1.plot(x, u_up, color=col, linewidth=1.5, linestyle="--", alpha=0.7)
    ax1.plot(x, u_dn, color=col, linewidth=1.5, linestyle="--", alpha=0.7, label=lbl)

ax1.axhline(y=1, color="white", linestyle=":", alpha=0.3)
ax1.axvline(x=1, color="white", linestyle=":", alpha=0.3)
ax1.plot(1, 1, "wo", markersize=10, label="Critical point")
ax1.set_xlabel("r / r_c")
ax1.set_ylabel("v / c_s")
ax1.set_title("Parker Equation: Solution Topology")
ax1.set_ylim(0, 5)
ax1.set_xlim(0.2, 8)
ax1.legend(fontsize=7, facecolor="#0a0a0a", edgecolor="#333", labelcolor="white")

# --- Panel 2: Velocity profiles for different temperatures ---
ax2 = axes[0, 1]

for T_MK, col in [(0.5, "#60a5fa"), (1.0, "#f43f5e"), (1.5, "#fb923c"), (2.0, "#fbbf24")]:
    T = T_MK * 1e6
    c_s = np.sqrt(2 * k_B * T / m_p)  # isothermal for fully ionized H
    r_c = G * M_sun / (2 * c_s**2)

r_arr = np.linspace(R_sun, 250 * R_sun, 2000)
    x_arr = r_arr / r_c

v_arr = np.zeros_like(r_arr)
    for i, xi in enumerate(x_arr):
        rhs_val = 4 * np.log(xi) + 4 / xi - 3
        try:
            if xi >= 1:
                v_arr[i] = brentq(lambda u: u**2 - np.log(u**2) - rhs_val, 1.0, 30.0) * c_s
            else:
                v_arr[i] = brentq(lambda u: u**2 - np.log(u**2) - rhs_val, 0.001, 1.0) * c_s
        except:
            v_arr[i] = np.nan

ax2.plot(r_arr / R_sun, v_arr / 1e3, color=col, linewidth=2,
            label="T = " + str(T_MK) + " MK")

ax2.axhline(y=400, color="white", linestyle=":", alpha=0.3)
ax2.axvline(x=215, color="gray", linestyle="--", alpha=0.3)
ax2.text(216, 50, "1 AU", color="gray", fontsize=8)
ax2.set_xlabel("r / R_sun")
ax2.set_ylabel("Wind speed [km/s]")
ax2.set_title("Parker Wind: v(r) for Different T")
ax2.set_xlim(1, 250)
ax2.set_ylim(0, 800)
ax2.legend(fontsize=8, facecolor="#0a0a0a", edgecolor="#333", labelcolor="white")

# --- Panel 3: Density and velocity at 1 AU ---
ax3 = axes[1, 0]

T_range = np.linspace(0.5, 3.0, 50)  # MK
v_1AU = np.zeros_like(T_range)
n_1AU = np.zeros_like(T_range)

for i, T_MK in enumerate(T_range):
    T = T_MK * 1e6
    c_s = np.sqrt(2 * k_B * T / m_p)
    r_c = G * M_sun / (2 * c_s**2)
    x_AU = AU / r_c

rhs_val = 4 * np.log(x_AU) + 4 / x_AU - 3
    try:
        u_AU = brentq(lambda u: u**2 - np.log(u**2) - rhs_val, 1.0, 50.0)
        v_1AU[i] = u_AU * c_s / 1e3  # km/s

# Density from mass conservation: rho * v * r^2 = const
        # rho_0 = n_0 * m_p ~ 10^14 * 1.67e-27 kg/m^3
        n_0 = 1e14  # m^-3 at base
        rho_0 = n_0 * m_p
        # At critical point: rho_c * c_s * r_c^2 = rho_0 * v_0 * R_sun^2
        # This is approximate
        n_1AU[i] = n_0 * R_sun**2 * 1e3 / (AU**2 * v_1AU[i] * 1e3)  # cm^-3
        n_1AU[i] = n_1AU[i] * 1e-6  # rough scaling
    except:
        v_1AU[i] = np.nan
        n_1AU[i] = np.nan

ax3_twin = ax3.twinx()
ax3.plot(T_range, v_1AU, color="#f43f5e", linewidth=2.5, label="Wind speed")
ax3_twin.plot(T_range, n_1AU, color="#60a5fa", linewidth=2.5, label="Density")
ax3.set_xlabel("Coronal temperature [MK]")
ax3.set_ylabel("v at 1 AU [km/s]", color="#f43f5e")
ax3_twin.set_ylabel("n at 1 AU [relative]", color="#60a5fa")
ax3.set_title("Wind Properties vs Coronal Temperature")
ax3.tick_params(axis="y", labelcolor="#f43f5e")
ax3_twin.tick_params(axis="y", labelcolor="#60a5fa")

# --- Panel 4: Critical radius vs temperature ---
ax4 = axes[1, 1]
T_cr = np.linspace(0.3, 5.0, 200)
c_s_cr = np.sqrt(2 * k_B * T_cr * 1e6 / m_p)
r_c_cr = G * M_sun / (2 * c_s_cr**2) / R_sun

ax4.plot(T_cr, r_c_cr, color="#f43f5e", linewidth=2.5)
ax4.axhline(y=15, color="#fbbf24", linestyle="--", alpha=0.5, label="r_A ~ 15 R_sun (PSP)")
ax4.axhline(y=215, color="#4ade80", linestyle="--", alpha=0.5, label="1 AU")
ax4.set_xlabel("Coronal temperature [MK]")
ax4.set_ylabel("Critical radius [R_sun]")
ax4.set_title("Critical Radius vs Temperature")
ax4.set_yscale("log")
ax4.set_ylim(1, 300)
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("PARKER SOLAR WIND MODEL")
print("=" * 60)

for T_MK in [0.5, 1.0, 1.5, 2.0]:
    T = T_MK * 1e6
    c_s = np.sqrt(2 * k_B * T / m_p)
    r_c = G * M_sun / (2 * c_s**2)
    print("\n--- T = " + str(T_MK) + " MK ---")
    print("  Sound speed c_s = " + str(round(c_s/1e3, 1)) + " km/s")
    print("  Critical radius r_c = " + str(round(r_c/R_sun, 1)) + " R_sun")

x_AU = AU / r_c
    rhs_val = 4 * np.log(x_AU) + 4 / x_AU - 3
    try:
        u_AU = brentq(lambda u: u**2 - np.log(u**2) - rhs_val, 1.0, 50.0)
        print("  v(1 AU) = " + str(round(u_AU * c_s / 1e3, 0)) + " km/s")
    except:
        print("  v(1 AU): no solution found")

print("\n--- Observed Solar Wind ---")
print("  Fast wind: 600-800 km/s (coronal holes)")
print("  Slow wind: 300-400 km/s (streamer belt)")
print("  Mass loss rate: ~2 x 10^-14 M_sun/yr")
print("  Alfven radius: ~15-20 R_sun (Parker Solar Probe)")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

8.6 Parker Solar Wind — The Complete Derivation

This section provides the most detailed treatment of the Parker wind in the course, covering every algebraic step from the basic equations to the full solution topology.

Starting Point: The Fluid Equations

Step 1. We assume a steady, spherically symmetric, isothermal, non-magnetic, non-rotating outflow of a fully ionised hydrogen plasma. The three governing equations are:

$$\text{Mass:} \quad \frac{d}{dr}\!\left(\rho v r^2\right) = 0 \;\Longrightarrow\; \rho v r^2 = \frac{\dot{M}}{4\pi} = \text{const}$$

$$\text{Momentum:} \quad \rho v\frac{dv}{dr} = -\frac{dP}{dr} - \frac{GM_\odot\rho}{r^2}$$

$$\text{Equation of state:} \quad P = n k_B T = \frac{\rho k_B T}{\mu m_H} = \rho c_s^2, \quad c_s = \sqrt{\frac{k_B T}{\mu m_H}}$$

For fully ionised hydrogen with equal electron and proton temperatures,$\mu = 1/2$ (mean molecular weight), so$c_s = \sqrt{2k_BT/m_p}$.

Eliminating Density

Step 2. Logarithmically differentiate the continuity equation$\rho v r^2 = \text{const}$:

$$\frac{d\ln\rho}{dr} + \frac{d\ln v}{dr} + \frac{2}{r} = 0 \;\Longrightarrow\; \frac{1}{\rho}\frac{d\rho}{dr} = -\frac{1}{v}\frac{dv}{dr} - \frac{2}{r}$$

Step 3. Since $P = \rho c_s^2$ with $c_s$ constant (isothermal):

$$\frac{dP}{dr} = c_s^2\frac{d\rho}{dr} = c_s^2\rho\!\left(-\frac{1}{v}\frac{dv}{dr} - \frac{2}{r}\right)$$

Deriving the Parker Equation

Step 4. Substitute this into the momentum equation:

$$\rho v\frac{dv}{dr} = -c_s^2\rho\!\left(-\frac{1}{v}\frac{dv}{dr} - \frac{2}{r}\right) - \frac{GM_\odot\rho}{r^2}$$

Step 5. Dividing through by $\rho$ and expanding:

$$v\frac{dv}{dr} = \frac{c_s^2}{v}\frac{dv}{dr} + \frac{2c_s^2}{r} - \frac{GM_\odot}{r^2}$$

Step 6. Collecting the $dv/dr$ terms on the left:

$$\boxed{\left(v - \frac{c_s^2}{v}\right)\frac{dv}{dr} = \frac{2c_s^2}{r} - \frac{GM_\odot}{r^2} = \frac{2c_s^2}{r}\left(1 - \frac{r_c}{r}\right)}$$

The Parker solar wind equation, with $r_c = GM_\odot/(2c_s^2)$

Critical Point Analysis

Step 7. This ODE has the form $f(v)\,dv/dr = g(r)$. The right side vanishes at $r = r_c$. The left side factor$(v - c_s^2/v) = (v^2 - c_s^2)/v$ vanishes when $v = c_s$.

Step 8. At the critical point, $dv/dr$ is of the form 0/0. We use l'Hopital's rule or direct Taylor expansion. Write$v = c_s + \delta v$, $r = r_c + \delta r$ and expand to first order:

$$2c_s\,\delta v \cdot \frac{dv}{dr}\bigg|_c = \frac{2c_s^2}{r_c^2}\,\delta r$$

Since $\delta v = (dv/dr)|_c\,\delta r$:

$$\left.\frac{dv}{dr}\right|_c = \pm\frac{c_s}{r_c}$$

The $+$ sign gives the accelerating (wind/accretion) solution; the $-$ sign gives a decelerating solution. Only the wind solution ($+$, subsonic for $r < r_c$, supersonic for $r > r_c$) satisfies $P \to 0$ as $r \to \infty$.

Analytical Integration

Step 9. Define dimensionless variables $u = v/c_s$ and$x = r/r_c$. The Parker equation becomes:

$$\left(u - \frac{1}{u}\right)\frac{du}{dx} = \frac{2}{x}\left(1 - \frac{1}{x}\right)$$

Step 10. This separates:$\left(u - 1/u\right)du = 2(1/x - 1/x^2)\,dx$. Integrating both sides:

$$\frac{u^2}{2} - \ln u = 2\ln x + \frac{2}{x} + \frac{C}{2}$$

Multiplying through by 2 gives the standard form:

$$\boxed{u^2 - \ln u^2 = 4\ln x + \frac{4}{x} + C}$$

Integrated Parker equation in dimensionless form

The Four Classes of Solutions

Step 11. The transonic solution passes through $u = 1$ at$x = 1$. Setting $u = 1$, $x = 1$:$1 - 0 = 0 + 4 + C$, so $C = -3$.

• Class I — Parker wind ($C = -3$): Subsonic below $r_c$, passes through the sonic point, supersonic above. This is the physically correct solar wind solution. Velocity at 1 AU: $\sim 400$ km/s for $T = 10^6$ K.
• Class II — Accretion ($C = -3$): The time-reverse of the Parker wind. Supersonic infall decelerating through the sonic point. Relevant for Bondi accretion.
• Class III — Subsonic breeze ($C > -3$): Always subsonic. Velocity increases, reaches a maximum below $c_s$, then decreases. Unphysical because$P \to \text{finite}$ as $r \to \infty$, which cannot be matched to the ISM.
• Class IV — Supersonic everywhere ($C > -3$): Always supersonic. Unphysical because it requires $v > c_s$ at the coronal base.

Critical Radius: Numerical Values

$$r_c = \frac{GM_\odot}{2c_s^2} = \frac{GM_\odot \mu m_H}{2k_BT}$$

• $T = 0.5$ MK: $r_c \approx 23\,R_\odot$, $c_s \approx 91$ km/s
• $T = 1.0$ MK: $r_c \approx 11.5\,R_\odot$, $c_s \approx 128$ km/s
• $T = 1.5$ MK: $r_c \approx 7.7\,R_\odot$, $c_s \approx 157$ km/s
• $T = 2.0$ MK: $r_c \approx 5.8\,R_\odot$, $c_s \approx 181$ km/s

8.7 Weber-Davis MHD Wind — Extended Derivation

The Weber-Davis (1967) model extends Parker's solution by including solar rotation ($\Omega_\odot$) and the azimuthal magnetic field ($B_\phi$).

Conservation Laws

Step 1. In the equatorial plane, the induction equation for a steady axisymmetric flow with perfect conductivity ($\mathbf{E} + \mathbf{v} \times \mathbf{B} = 0$) gives:

$$v_r B_\phi - v_\phi B_r = -\Omega r B_r$$

which means field lines co-rotate with the Sun below the Alfven point. Solving for $B_\phi$:

$$B_\phi = B_r\frac{v_\phi - \Omega r}{v_r}$$

Step 2. The conserved specific angular momentum (including magnetic torque):

$$\boxed{\mathcal{L} = r v_\phi - \frac{r B_r B_\phi}{\mu_0 \rho v_r} = \Omega r_A^2}$$

Total specific angular momentum = that of rigid rotation out to the Alfven radius

The Alfven Radius

Step 3. Define the radial Alfven Mach number$M_A = v_r / v_{A,r}$ where$v_{A,r} = B_r / \sqrt{\mu_0 \rho}$. At $r = r_A$, $M_A = 1$ (the wind speed equals the Alfven speed).

Step 4. The $B_\phi$ expression has a singularity at $v_r = v_{A,r}$(the denominator vanishes). Regularity requires the numerator to also vanish at$r_A$, giving $v_\phi(r_A) = \Omega r_A$ — co-rotation at the Alfven point.

Angular Momentum Loss Rate

Step 5. The total angular momentum loss rate of the Sun is:

$$\boxed{\dot{J} = \dot{M}\,\Omega\,r_A^2}$$

The wind acts as a magnetic brake, removing angular momentum as if the wind were co-rotating rigidly out to $r_A$. With $r_A \approx 15\,R_\odot$:

$$\frac{\dot{J}}{\dot{M}} = \Omega r_A^2 \approx 2.9 \times 10^{-6} \times (15 \times 6.96 \times 10^8)^2 \approx 3.2 \times 10^{16} \text{ m}^2\text{/s}$$

Step 6. This angular momentum loss spins down the Sun on a timescale:

$$\tau_{\text{spin-down}} = \frac{J_\odot}{\dot{J}} \sim \frac{I_\odot \Omega}{\dot{M}\Omega r_A^2} = \frac{I_\odot}{\dot{M} r_A^2} \sim 10^{10} \text{ years}$$

This is consistent with the observed spin-down of solar-type stars (Skumanich 1972:$\Omega \propto t^{-1/2}$). Parker Solar Probe crossed the Alfven surface at$r_A \approx 13\text{--}20\,R_\odot$ in 2021, directly confirming this key prediction.

Diagram: Parker Wind Solution Topology

The four classes of Parker wind solutions meeting at the critical point where$v = c_s$ and $r = r_c$. Only the transonic Parker wind (Class I, solid red) satisfies the physical boundary conditions of a subsonic base and vanishing pressure at infinity.

Advanced Simulation: Parker Wind Velocity, Parameter Space, and Weber-Davis

Parker Wind Solution Topology, Parameter Space Explorer, and Weber-Davis MHD Wind

Python

script.py230 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.optimize import brentq

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle('Parker Wind & Weber-Davis MHD Wind',
             fontsize=15, 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')

G = 6.674e-11
M_sun = 1.989e30
R_sun = 6.957e8
k_B = 1.381e-23
m_p = 1.673e-27
AU = 1.496e11

def parker_eq(u, rhs_val):
    if u <= 0:
        return 1e10
    return u**2 - np.log(u**2) - rhs_val

# === Panel 1: Full solution topology (all 4 classes) ===
ax1 = axes[0, 0]
x = np.linspace(0.15, 6, 600)

# Transonic solution (C = -3)
u_sup = np.zeros_like(x)
u_sub = np.zeros_like(x)
for i, xi in enumerate(x):
    rhs = 4 * np.log(xi) + 4.0 / xi - 3.0
    try:
        u_sup[i] = brentq(parker_eq, 1.001, 30.0, args=(rhs,))
    except:
        u_sup[i] = np.nan
    try:
        u_sub[i] = brentq(parker_eq, 0.001, 0.999, args=(rhs,))
    except:
        u_sub[i] = np.nan

ax1.plot(x, u_sup, color='#f43f5e', linewidth=3, label='I: Parker wind')
ax1.plot(x, u_sub, color='#f43f5e', linewidth=3)
# The accretion solution is the same curve (reverse arrow sense)
ax1.plot(x, u_sup, color='#fb923c', linewidth=2, linestyle=':', alpha=0.0)
ax1.plot(x, u_sub, color='#fb923c', linewidth=2, linestyle=':', alpha=0.0)

# Breeze and supersonic solutions
for C_val, col, lbl, alpha_v in [(-2.0, '#4ade80', 'III: Breeze', 0.7),
                                   (-1.0, '#4ade80', '', 0.5),
                                   (-4.5, '#a78bfa', 'IV: Supersonic', 0.7),
                                   (-6.0, '#a78bfa', '', 0.5)]:
    u_a = np.zeros_like(x)
    u_b = np.zeros_like(x)
    for i, xi in enumerate(x):
        rhs = 4 * np.log(xi) + 4.0 / xi + C_val
        try:
            u_a[i] = brentq(parker_eq, 1.001, 30.0, args=(rhs,))
        except:
            u_a[i] = np.nan
        try:
            u_b[i] = brentq(parker_eq, 0.001, 0.999, args=(rhs,))
        except:
            u_b[i] = np.nan
    ax1.plot(x, u_a, color=col, linewidth=1.5, linestyle='--', alpha=alpha_v)
    lab = lbl if lbl else None
    ax1.plot(x, u_b, color=col, linewidth=1.5, linestyle='--', alpha=alpha_v, label=lab)

ax1.axhline(y=1, color='white', linestyle=':', alpha=0.2)
ax1.axvline(x=1, color='white', linestyle=':', alpha=0.2)
ax1.plot(1, 1, 'o', color='#fbbf24', markersize=10, zorder=5, label='Critical point')
ax1.set_xlabel('r / r_c')
ax1.set_ylabel('v / c_s')
ax1.set_title('Solution Topology (all 4 classes)')
ax1.set_ylim(0, 4)
ax1.set_xlim(0.15, 6)
ax1.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# === Panel 2: Parker wind v(r) for different T ===
ax2 = axes[0, 1]
for T_MK, col in [(0.5, '#60a5fa'), (1.0, '#f43f5e'), (1.5, '#fb923c'), (2.0, '#fbbf24'), (3.0, '#a78bfa')]:
    T = T_MK * 1e6
    c_s = np.sqrt(2 * k_B * T / m_p)
    r_c = G * M_sun / (2 * c_s**2)
    r_arr = np.linspace(R_sun, 250 * R_sun, 1500)
    x_arr = r_arr / r_c
    v_arr = np.zeros_like(r_arr)
    for i, xi in enumerate(x_arr):
        rhs = 4 * np.log(xi) + 4.0 / xi - 3.0
        try:
            if xi >= 1:
                v_arr[i] = brentq(parker_eq, 1.0, 30.0, args=(rhs,)) * c_s
            else:
                v_arr[i] = brentq(parker_eq, 0.001, 1.0, args=(rhs,)) * c_s
        except:
            v_arr[i] = np.nan
    ax2.plot(r_arr / R_sun, v_arr / 1e3, color=col, linewidth=2,
            label='T = ' + str(T_MK) + ' MK')

ax2.axvline(x=215, color='gray', linestyle='--', alpha=0.3)
ax2.text(220, 30, '1 AU', color='gray', fontsize=9)
ax2.axhline(y=400, color='white', linestyle=':', alpha=0.2)
ax2.axhline(y=700, color='white', linestyle=':', alpha=0.2)
ax2.set_xlabel('r / R_sun')
ax2.set_ylabel('Wind speed [km/s]')
ax2.set_title('Parker Wind: v(r) for Different Temperatures')
ax2.set_xlim(1, 250)
ax2.set_ylim(0, 900)
ax2.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# === Panel 3: v(1 AU) and r_c vs T — parameter space ===
ax3 = axes[1, 0]
T_range = np.linspace(0.3, 4.0, 100)
v_1AU = np.zeros_like(T_range)
r_c_arr = np.zeros_like(T_range)

for i, T_MK in enumerate(T_range):
    T = T_MK * 1e6
    c_s = np.sqrt(2 * k_B * T / m_p)
    r_c = G * M_sun / (2 * c_s**2)
    r_c_arr[i] = r_c / R_sun
    x_AU = AU / r_c
    rhs = 4 * np.log(x_AU) + 4.0 / x_AU - 3.0
    try:
        u_AU = brentq(parker_eq, 1.0, 50.0, args=(rhs,))
        v_1AU[i] = u_AU * c_s / 1e3
    except:
        v_1AU[i] = np.nan

ax3.plot(T_range, v_1AU, color='#f43f5e', linewidth=2.5, label='v at 1 AU')
ax3_twin = ax3.twinx()
ax3_twin.plot(T_range, r_c_arr, color='#60a5fa', linewidth=2, linestyle='--', label='r_c / R_sun')
ax3.set_xlabel('Coronal temperature [MK]')
ax3.set_ylabel('v(1 AU) [km/s]', color='#f43f5e')
ax3_twin.set_ylabel('r_c [R_sun]', color='#60a5fa')
ax3.tick_params(axis='y', labelcolor='#f43f5e')
ax3_twin.tick_params(axis='y', labelcolor='#60a5fa')
ax3.set_title('Parameter Space: v(1 AU) and r_c vs T')
ax3_twin.set_yscale('log')

# Observed range
ax3.axhspan(300, 400, alpha=0.1, color='#fb923c')
ax3.text(0.5, 350, 'Slow wind', color='#fb923c', fontsize=8)
ax3.axhspan(600, 800, alpha=0.1, color='#fbbf24')
ax3.text(0.5, 700, 'Fast wind', color='#fbbf24', fontsize=8)

# === Panel 4: Weber-Davis — azimuthal velocity and field ===
ax4 = axes[1, 1]
# Simplified Weber-Davis model
Omega = 2.9e-6  # rad/s
r_A = 15 * R_sun  # Alfven radius
r_range = np.linspace(R_sun, 250 * R_sun, 500)

# Parker wind velocity for T=1 MK
T_wd = 1.0e6
c_s_wd = np.sqrt(2 * k_B * T_wd / m_p)
r_c_wd = G * M_sun / (2 * c_s_wd**2)
v_r = np.zeros_like(r_range)
for i in range(len(r_range)):
    xi = r_range[i] / r_c_wd
    rhs = 4 * np.log(xi) + 4.0 / xi - 3.0
    try:
        if xi >= 1:
            v_r[i] = brentq(parker_eq, 1.0, 30.0, args=(rhs,)) * c_s_wd
        else:
            v_r[i] = brentq(parker_eq, 0.001, 1.0, args=(rhs,)) * c_s_wd
    except:
        v_r[i] = np.nan

# Alfven Mach number
M_A = np.where(r_range < r_A, r_range / r_A * np.sqrt(v_r / v_r[0]), v_r / v_r[0])
# Simplified: M_A ~ v_r * r^2 / (v_A0 * R_sun^2)

# Azimuthal velocity (approximate)
v_phi = Omega * r_range * np.where(r_range < r_A, 1.0, (r_A / r_range)**2)

# Parker spiral angle
psi = np.arctan2(Omega * r_range, v_r) * 180 / np.pi

ax4.plot(r_range / R_sun, v_phi / 1e3, color='#fb923c', linewidth=2.5, label='v_phi [km/s]')
ax4.plot(r_range / R_sun, v_r / 1e3, color='#f43f5e', linewidth=2, linestyle='--', label='v_r [km/s]')
ax4_twin = ax4.twinx()
ax4_twin.plot(r_range / R_sun, psi, color='#4ade80', linewidth=1.5, linestyle=':', label='Spiral angle')
ax4_twin.set_ylabel('Parker spiral angle [deg]', color='#4ade80')
ax4_twin.tick_params(axis='y', labelcolor='#4ade80')

ax4.axvline(x=15, color='#fbbf24', linestyle='--', alpha=0.4)
ax4.text(16, 300, 'r_A', color='#fbbf24', fontsize=10)
ax4.set_xlabel('r / R_sun')
ax4.set_ylabel('Velocity [km/s]')
ax4.set_title('Weber-Davis: v_r, v_phi, and Spiral Angle')
ax4.set_xlim(1, 250)
ax4.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white', loc='center right')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("Plot saved to output.png")
print("")
print("=== Parker Wind: Velocity at 1 AU ===")
for T_MK in [0.5, 1.0, 1.5, 2.0, 3.0]:
    T = T_MK * 1e6
    c_s = np.sqrt(2 * k_B * T / m_p)
    r_c = G * M_sun / (2 * c_s**2)
    x_AU = AU / r_c
    rhs = 4 * np.log(x_AU) + 4.0 / x_AU - 3.0
    try:
        u = brentq(parker_eq, 1.0, 50.0, args=(rhs,))
        v = u * c_s / 1e3
        print("  T = " + str(T_MK) + " MK: c_s = " + str(round(c_s/1e3)) + " km/s, r_c = " + str(round(r_c/R_sun, 1)) + " R_sun, v(1AU) = " + str(round(v)) + " km/s")
    except:
        print("  T = " + str(T_MK) + " MK: no transonic solution")
print("")
print("=== Weber-Davis Model ===")
print("  Alfven radius r_A ~ 15 R_sun")
print("  Angular momentum loss: J_dot = M_dot * Omega * r_A^2")
rA_m = 15 * R_sun
Jdot_per_Mdot = Omega * rA_m**2
print("  J_dot / M_dot = " + str("%.2e" % Jdot_per_Mdot) + " m^2/s")
print("  Parker spiral angle at 1 AU ~ " + str(round(np.degrees(np.arctan(Omega * AU / 400e3)))) + " degrees")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Corona Next: Magnetic Fields of the Sun →