Solve and check a critical atmospheric wind

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.

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.

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.

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