---
title: "Solve and check a critical atmospheric wind"
description: "Run a standalone Python critical-wind benchmark, compare a regular ODE integration with Lambert-W branches and diagnose failed physical assumptions."
canonical_url: "https://sandboxphysics.com/stories/atmospheric-escape/solve-a-wind/"
language: "en"
page_type: "story"
alternate_languages:
  zh-CN: "https://sandboxphysics.com/zh/stories/atmospheric-escape/solve-a-wind/"
last_modified: "2026-10-28"
---

[Planetary atmospheres from thermal tails to winds](https://sandboxphysics.com/stories/atmospheric-escape/)

AC13 · L3

# Solve and check a critical atmospheric wind

Sandbox Physics · 2026-10-28

Find the critical radius numerically, integrate the accelerating solution from both sides of it, and test its conserved quantities. Then compare the power requirement and collision condition before interpreting a computed rate.

## Before you begin

Read Build a wind through its sonic point. The basic replay needs Python 3.10+ and its standard library; the optional analytic comparison uses the pinned NumPy/SciPy versions below. No repository checkout is required.

## Save the input before running anything

Download the scripts, input, provenance, instructions and frozen outputs listed below into one ordinary folder. The input defines five synthetic single-component atmospheres: default wind, a strongly bound tail, a sparse atmosphere, an XUV power shortfall and a critical point below the reservoir. No observed planet profile or MESA evolution output enters this experiment.

The JSON fixes SI constants, boundary radius, density, prescribed temperature, XUV irradiance and usable power fraction. The script hashes the exact input bytes into the output. Keep the downloaded input unchanged for the first run; use a separately named copy for experiments. The provenance file lists the byte counts and SHA-256 values and distinguishes authored teaching inputs from the papers that motivate the equations.

```
python3 -W error replay-wind.py --input escape-input.json --out wind-replay-rerun.json
# Optional analytic reference, in a virtual environment:
python3 -m venv .venv
source .venv/bin/activate
python3 -m pip install -r requirements-reference.txt
python3 reference-lambert.py --out lambert-reference-rerun.json
```

Reading: Cranmer 2004, New views of the solar wind with the Lambert W function, v1 · Volkov et al. 2011, Thermally driven atmospheric escape, arXiv v1

## Bracket the place where both terms vanish

The critical radius is found by bisection on $2c^2/r-GM/r^2$ between $0.01R_{\rm b}$ and $100R_{\rm b}$. At every step retain the half-interval with a sign change. A missing sign change raises an error. All five fixed examples are bracketed; arbitrary new inputs must be checked.

Compare the result with $GM/(2c^2)$. For the default it is $2.41574013\times10^7\,\mathrm m$. This root only locates the zero of the numerator. The velocity still has to pass through $u=c$ on the accelerating branch. Directly dividing numerator by denominator at that point produces zero divided by zero.

Reading: Cranmer 2004, New views of the solar wind with the Lambert W function, v1

## Start with a regular slope, not a singular division

Integrate $q=\ln\mathcal M$ instead of velocity. This keeps the state finite even for extremely small subsonic speeds. Start at $x=1\pm10^{-5}$ with the local expansion below, then integrate inward or outward using fourth-order Runge–Kutta. Targets closer than $10^{-4}$ to the sonic point use the expansion directly.

The maximum radial step is $0.0015\min(x,1)$. Each requested radius is integrated independently from the regular sonic neighborhood. The output also halves this maximum step at seven radii from $0.03r_{\rm s}$ to $12r_{\rm s}$, and saves the change in log Mach number. This checks numerical sensitivity of the ODE path; its agreement with a separate analytic solution tests a different implementation.

$$
\frac{dq}{dx}=\frac{2(1/x-1/x^2)}{e^{2q}-1},\quad q=\delta-\delta^2+\frac{13}{12}\delta^3+O(\delta^4)
$$

The exponential-minus-one function avoids loss of precision in the small denominator away from the exact critical point. The series handles the critical neighborhood.

Reading: Cranmer 2004, New views of the solar wind with the Lambert W function, v1

## Compare with two analytic branches

Set $D=x^{-4}e^{3-4/x}$. The implicit Parker relation can be inverted using the Lambert-W function, defined by $W(z)e^{W(z)}=z$. In the interval of interest it has two real branches. The branch names in SciPy are zero and minus one.

The optional script evaluates the principal branch below the sonic radius and the other real branch above it. At exactly $x=1$ it assigns the regular limit explicitly. It compares eleven radii with the independently integrated ODE. At twice the sonic radius, choosing the principal branch gives a subsonic root instead of the intended accelerating wind. Both can satisfy an implicit equation, which is why a small algebraic residual alone cannot identify the physical branch.

$$
\mathcal M^2=-W_0(-D)\quad(x<1),\qquad\mathcal M^2=-W_{-1}(-D)\quad(x>1)
$$

The full accelerating solution switches real branches at their common sonic endpoint.

Reading: Cranmer 2004, New views of the solar wind with the Lambert W function, v1

## Test the flow and the approximations separately

Recover density from the base boundary and calculate $4\pi r^2\rho u$ at all saved radii. Also compute $u^2/2+c^2\ln(n/n_{\rm b})-GM/r$, the isothermal momentum integral. The first check tests the density normalization; the second responds to errors in the independently integrated velocity shape.

The frozen default has relative mass-flux spread of order ten to the minus fifteen and a momentum-integral spread below ten to the minus ten in sound-speed-squared units. Those are numerical residuals. They do not include an uncertainty in temperature, the collision cross section or the omission of cooling.

The table and curve below use the public replay output. The strongly bound case gives an illustrative Jeans rate near ten milligrams per second, but its formal Parker sonic Knudsen number is enormous. The shortfall case has a well-behaved numerical wind and a power requirement about 752 times the available XUV power. The buried-critical case is already supersonic at the reservoir. The output retains each failed interpretation beside its numerical results.

Independent regularized ODE result, default synthetic inputs. Horizontal axis is radius divided by the reservoir radius, logarithmic; the dashed line marks the numerically located critical radius.

Actual synthetic results from replay-wind.py: fixed neutral-H-like particle mass and effective collision cross section, prescribed temperature, two-Earth-radius reservoir, outer radius twelve times the reservoir. Inputs include the full settings for each case. Formal wind rates are listed even when their physical assumptions fail.

| Case | Sonic radius / base | Formal Parker rate · kg/s | Sonic Knudsen | Required / available power |
| --- | --- | --- | --- | --- |
| Default wind | $1.9$ | $7.70\times10^{6}$ | $4.02\times10^{-3}$ | $2.28\times10^{-1}$ |
| Rare tail | $1.52\times10^{1}$ | $6.77\times10^{-4}$ | $1.83\times10^{8}$ | $2.05\times10^{-9}$ |
| Sparse gas | $5.06$ | $7.90\times10^{2}$ | $9.05\times10^{1}$ | $2.80\times10^{-4}$ |
| Heat shortfall | $1.26$ | $1.84\times10^{7}$ | $1.37\times10^{-3}$ | $7.52\times10^{2}$ |
| Buried critical point | $2.11\times10^{-1}$ | $7.11\times10^{7}$ | $7.26\times10^{-5}$ | $3.78$ |

Reading: Volkov et al. 2011, Thermally driven atmospheric escape, arXiv v1 · Murray-Clay, Chiang & Murray2009, Atmospheric escape from hot Jupiters, v2

## Make one controlled modification

Copy the input and double only the base number density in the default case. Predict the wind rate, the sonic Knudsen number and the required-to-available power ratio before running. The respective changes should be factors of two, one half and two. The velocity shape stays fixed because its equation contains gravity and sound speed, not the density normalization.

For a second experiment, change only the declared outer radius from twelve to twenty-four base radii. Re-run the mechanical-energy comparison. More work and acceleration must be supplied as gas travels farther. Keep track of this outer boundary whenever you quote an isothermal heat requirement.

The public files may be reused as the scientific teaching material described in their README. They contain standalone numerical analysis, not a copy of the protected interactive Lab. To study a particular planet, replace the prescribed temperature model with an energy/chemistry calculation and establish appropriate lower and outer boundaries.

Reading: Murray-Clay, Chiang & Murray2009, Atmospheric escape from hot Jupiters, v2 · Cranmer 2004, New views of the solar wind with the Lambert W function, v1

## Exercises and explanations

1. Your implicit equation residual is tiny, but the speed falls below the sound speed beyond the critical radius. Which test catches the mistake?

   Read the explanation

   Check the branch and the sign of acceleration. At x>1 the accelerating transonic solution must use the supersonic Lambert-W branch. The principal branch can satisfy the same implicit relation, so equation residual alone is insufficient.

2. The output input hash differs from the frozen reference. Can you compare all last digits directly?

   Read the explanation

   First identify the input change. A different boundary, density or temperature is a different calculation. Restore the fixed input for a numerical reproduction, then use a separate modified run for a physical experiment.

## Try the prediction

### [Compare the saved wind with the interactive one](https://sandboxphysics.com/labs/planetary-atmospheric-escape/)

1. Reset the Lab and export its calculation. Compare the default critical radius, mass flow, power ratio and sonic Knudsen number with the Python output.
2. Choose Collisions run out and identify which mathematical curves remain computable after the fluid interpretation fails.
3. Change only the base density and test the predicted factors before trying a different temperature.

[Previous: Build a wind through its sonic point](https://sandboxphysics.com/stories/atmospheric-escape/critical-wind/) [Next: Read the competing explanations of the radius valley](https://sandboxphysics.com/stories/atmospheric-escape/radius-valley/)

Save these files together. The numerical replay uses the Python3.10+ standard library. Only the optional Lambert-W comparison needs the pinned NumPy/SciPy dependencies.

- [Independent critical ODE replay](https://sandboxphysics.com/stories/atmospheric-escape/downloads/replay-wind.py)
- [Fixed synthetic atmosphere inputs](https://sandboxphysics.com/stories/atmospheric-escape/downloads/escape-input.json)
- [Frozen numerical results and failed regimes](https://sandboxphysics.com/stories/atmospheric-escape/downloads/wind-replay.json)
- [Independent Lambert-W comparison](https://sandboxphysics.com/stories/atmospheric-escape/downloads/reference-lambert.py)
- [Frozen analytic comparisons](https://sandboxphysics.com/stories/atmospheric-escape/downloads/lambert-reference.json)
- [Optional pinned dependencies](https://sandboxphysics.com/stories/atmospheric-escape/downloads/requirements-reference.txt)
- [Versions, units, permissions and hashes](https://sandboxphysics.com/stories/atmospheric-escape/downloads/provenance.json)
- [Standalone instructions and scope](https://sandboxphysics.com/stories/atmospheric-escape/downloads/README.md)

## Sources and calculation scope

1. [Cranmer 2004, New views of the solar wind with the Lambert W function, v1](https://arxiv.org/pdf/astro-ph/0406176v1)

   7June2004; AJP 72,1397–1403. Full text, SectionsIII–IV, equations12–16: isothermal momentum equation, critical point and the two real Lambert-W branches. The solar ionized-hydrogen sound-speed convention is replaced by our explicitly fixed single-particle mass.

2. [Volkov et al. 2011, Thermally driven atmospheric escape, arXiv v1](https://arxiv.org/pdf/1009.5110v1)

   26September2010 preprint, published ApJL729,L24(2011). Full text, Section 2 and Figures1–5: hard-sphere mean free path, radius-based versus local scale-height Knudsen numbers, kinetic escape and conditional Jeans corrections. Their heating below the lower boundary differs from a globally maintained isothermal wind. The adopted Lab collision cross section is a teaching input, not their molecular-nitrogen value.

3. [Murray-Clay, Chiang & Murray2009, Atmospheric escape from hot Jupiters, v2](https://arxiv.org/pdf/0811.0006v2)

   18November2008 revision, ApJ 693,23(2009). Full text, Section 3.2, equations19/22 and Figure 7: absorbing area, energy-limited scaling and the cooling/recombination limitations of a fixed efficiency. This source develops hot-Jupiter radiation hydrodynamics; the small-planet illustration uses only the explicitly stated global power comparison.