Skip to main content
Sandbox Physics
Planetary atmospheres from thermal tails to winds

AC13 · L3

Build a wind through its sonic point

Sandbox Physics ·

A smooth accelerating isothermal wind must cross the sound speed where gravity and the spherical pressure term balance. Continuity fixes its density shape; a lower-boundary density fixes its mass flow. Heating must then be checked separately.

Use continuity to eliminate density

Let rr increase outward from a fixed spherical gas reservoir of radius RbR_{\rm b}. Radial bulk velocity uu is positive outward, mass density ρ\rho is in kilograms per cubic meter and P=ρc2P=\rho c^2 is pressure. Temperature and particle mass stay fixed, so c2=kT/mc^2=kT/m is constant. Gravity points inward. Assume steady, spherical, inviscid flow, with no mass added along the path.

The mass inside each thin shell stays constant, so inflow and outflow through its faces match. This gives the same mass per second through every sphere. Differentiate its logarithm to replace the density gradient in the radial momentum equation.

M˙=4πr2ρu,dln⁡ρdr=−2r−1ududr\dot M=4\pi r^2\rho u,\qquad\frac{d\ln\rho}{dr}=-\frac2r-\frac1u\frac{du}{dr}
The two geometric powers come from the sphere area. Acceleration and geometric expansion both lower the density.
ududr=−c2dln⁡ρdr−GMr2⟹(u−c2u)dudr=2c2r−GMr2u\frac{du}{dr}=-c^2\frac{d\ln\rho}{dr}-\frac{GM}{r^2}\quad\Longrightarrow\quad\left(u-\frac{c^2}{u}\right)\frac{du}{dr}=\frac{2c^2}{r}-\frac{GM}{r^2}
The left side is advective acceleration. Substituting continuity moves the velocity-gradient pressure contribution to the left.

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

Pass smoothly through the critical point

At u=cu=c the coefficient of the velocity derivative vanishes. A finite slope therefore requires the right side to vanish at the same point, giving rs=GM/(2c2)r_{\rm s}=GM/(2c^2). This condition locates the critical point; it does not yet choose the direction of the crossing.

Write x=r/rsx=r/r_{\rm s} and M=u/c\mathcal M=u/c. Close to the critical point, put M=1+a(x−1)\mathcal M=1+a(x-1) in the dimensionless momentum equation. Its leading terms are 2a2(x−1)=2(x−1)2a^2(x-1)=2(x-1), so a=±1a=\pm1. The positive sign is the accelerating wind. Numerically, its regular expansion supplies a starting value just to either side of the singular point.

(M−1M)dMdx=2(1x−1x2),dudr∣s=crs\left(\mathcal M-\frac1{\mathcal M}\right)\frac{d\mathcal M}{dx}=2\left(\frac1x-\frac1{x^2}\right),\quad\left.\frac{du}{dr}\right|_{\rm s}=\frac c{r_{\rm s}}
The original zero divided by zero is replaced by the regular positive slope.
M=1+δ−δ22+δ34+O(δ4),δ=x−1\mathcal M=1+\delta-\frac{\delta^2}{2}+\frac{\delta^3}{4}+O(\delta^4),\qquad\delta=x-1
This local expansion picks the accelerating branch before numerical integration begins.

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

Choose the correct root and normalize the gas

Integrate the dimensionless equation once, choosing the constant so that x=1x=1 gives M=1\mathcal M=1. With w=M2w=\mathcal M^2, the result below usually has two positive roots. The accelerating solution takes the subsonic root inside and the supersonic root outside. Taking the smaller root everywhere creates a different solution.

The browser solves for ln⁡w\ln w on separate negative and positive intervals. An independent reference uses the principal Lambert-W branch inside and its other real branch outside. Once velocity is known, divide the continuity expression at radius rr by the same expression at the reservoir to recover density.

w−ln⁡w=4ln⁡x+4x−3w-\ln w=4\ln x+\frac4x-3
Both sides equal one at the sonic point. The sign of log w selects the branch.
n(r)=nbubu(r)(Rbr)2,M˙P=4πRb2mnbubn(r)=n_{\rm b}\frac{u_{\rm b}}{u(r)}\left(\frac{R_{\rm b}}r\right)^2,\qquad\dot M_{\rm P}=4\pi R_{\rm b}^2mn_{\rm b}u_{\rm b}
The base number density sets the mass normalization. It is independent of the stellar heating fraction.

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

Work through the default atmosphere

Take GM=4(3.986004418×1014) m3 s−2GM=4(3.986004418\times10^{14})\,\mathrm{m^3\,s^{-2}}, Rb=2(6.371×106) mR_{\rm b}=2(6.371\times10^6)\,\mathrm m, T=4000 KT=4000\,\mathrm K and nb=1015 m−3n_{\rm b}=10^{15}\,\mathrm{m^{-3}}. With the adopted particle mass 1.6735×10−27 kg1.6735\times10^{-27}\,\mathrm{kg}, the sound speed is 5.7446 km s−15.7446\,\mathrm{km\,s^{-1}}. The critical radius is 2.41574×107 m2.41574\times10^7\,\mathrm m or 1.89589Rb1.89589R_{\rm b}.

Choose the subsonic root at the base, multiply its speed by density and area, and obtain M˙P=7.69715×106 kg s−1\dot M_{\rm P}=7.69715\times10^6\,\mathrm{kg\,s^{-1}}. At the sonic point the logarithmic density-gradient length is Lρ=rs/3L_\rho=r_{\rm s}/3: two terms come from area expansion and one from the regular acceleration. With the adopted collision cross section, Kns=0.00402\mathrm{Kn}_{\rm s}=0.00402. A threshold below 0.10.1 is a necessary continuum check used here.

The table below is generated from the public independent ODE replay, with the same synthetic input settings. Extra stored digits describe numerical reproducibility, not atmospheric accuracy. The two cold examples show why a mathematically defined Parker rate can be a poor physical description.

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.
CaseSonic radius / baseFormal Parker rate · kg/sSonic KnudsenRequired / available power
Default wind1.91.97.70×1067.70\times10^{6}4.02×10−34.02\times10^{-3}2.28×10−12.28\times10^{-1}
Rare tail1.52×1011.52\times10^{1}6.77×10−46.77\times10^{-4}1.83×1081.83\times10^{8}2.05×10−92.05\times10^{-9}
Sparse gas5.065.067.90×1027.90\times10^{2}9.05×1019.05\times10^{1}2.80×10−42.80\times10^{-4}
Heat shortfall1.261.261.84×1071.84\times10^{7}1.37×10−31.37\times10^{-3}7.52×1027.52\times10^{2}
Buried critical point2.11×10−12.11\times10^{-1}7.11×1077.11\times10^{7}7.26×10−57.26\times10^{-5}3.783.78

Supply the mechanical energy

Let incident high-energy irradiance FXUVF_{\rm XUV} be power per area perpendicular to the rays. A planet intercepts a disk, so absorbed usable power is Pav=ηπRXUV2FXUVP_{\rm av}=\eta\pi R_{\rm XUV}^2F_{\rm XUV}. The illustrative choice RXUV=RbR_{\rm XUV}=R_{\rm b} must be stated; a real absorption radius can differ from a transit radius. Dividing by the base binding energy per kilogram gives the conventional energy-limited scale.

For the independently prescribed-temperature wind, compare the mechanical-energy increase only to ro=12Rbr_{\rm o}=12R_{\rm b}. Constant temperature makes the specific enthalpy equal at the endpoints, so the net required heat is the kinetic plus gravitational increase. A globally isothermal wind needs continuing heat and an unbounded energy supply as the outer radius tends to infinity.

M˙EL=ηπRXUV2FXUVGM/Rb\dot M_{\rm EL}=\frac{\eta\pi R_{\rm XUV}^2F_{\rm XUV}}{GM/R_{\rm b}}
This is a power scale, separate from the density-normalized transport calculation.
Preq=M˙P[uo2−ub22+GM ⁣(1Rb−1ro)]P_{\rm req}=\dot M_{\rm P}\left[\frac{u_{\rm o}^2-u_{\rm b}^2}{2}+GM\!\left(\frac1{R_{\rm b}}-\frac1{r_{\rm o}}\right)\right]
The bracket has units of joules per kilogram. This is the heat that must be supplied along the computed flow to the finite outer radius.

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

Interpret a result before using it

At FXUV=100 W m−2F_{\rm XUV}=100\,\mathrm{W\,m^{-2}} and η=0.15\eta=0.15, the default gives M˙EL=6.11443×107 kg s−1\dot M_{\rm EL}=6.11443\times10^7\,\mathrm{kg\,s^{-1}} and Preq/Pav=0.228P_{\rm req}/P_{\rm av}=0.228. The available power could supply this mechanical increase. Temperature, ionization and cooling still need an energy equation.

Changing only irradiation changes this comparison, because the transport temperature has been independently prescribed. In The missing heating, the ratio rises to about 752: that mathematical wind needs another heat supply or a different temperature structure. If the critical point lies below the reservoir, the chosen branch is already supersonic at the lower boundary. If the sonic region is collisionless, fluid transport must be replaced by a kinetic treatment. These tests tell us which next calculation is needed.

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

Exercises and explanations

  1. Double the base density at fixed mass, radius and temperature. What happens to the Parker mass flow and sonic Knudsen number?

    Read the explanation

    The velocity profile stays fixed. Mass flow doubles and the mean free path halves, so the sonic Knudsen number halves. Required power doubles; available stellar power stays fixed.

  2. At fixed particle mass, double temperature. What happens to rsr_{\rm s} and λb\lambda_{\rm b}?

    Read the explanation

    Both halve, because sound speed squared is proportional to temperature. The entire velocity profile must be recomputed; the mass-loss rate does not follow from this simple factor alone.

  3. A calculation has constant mass flux to twelve digits but sonic Knudsen number90. Is it a verified physical wind?

    Read the explanation

    It is a numerically consistent fluid solution in a region where the fluid approximation fails. The conservation residual tests arithmetic within the equations; the collision criterion tests whether those equations describe the gas.

Try the prediction

Separate transport from heating

  1. Predict what changes when only base density doubles. Open the advanced controls and test it.
  2. Select The missing heating. Increase XUV power and compare the power ratio with the unchanged velocity curve.
  3. Select A buried critical point and identify the supersonic state at the reservoir.

Sources and calculation scope

  1. Cranmer 2004, New views of the solar wind with the Lambert W function, v1

    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

    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

    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.

FeedbackSomething broken, confusing, or missing on this page?Tell us
What kind of feedback?

Sent with this page’s address, its Lab ID, the site language and your device type. Please leave out personal details. Feedback is kept for up to 12 months.Privacy policy