Build a star

Writing ρ = ρ_c θ^n and measuring the radius in units of α² = (n + 1)K ρ_c^(1/n − 1)/(4πG) turns hydrostatic equilibrium and mass conservation into the Lane-Emden equation,

(1/ξ²) d/dξ (ξ² dθ/dξ) + θ^n = 0, with θ(0) = 1 and θ′(0) = 0.

The solver integrates it by fourth-order Runge-Kutta from a series start, in steps of 0.01 that grow with ξ, and halves the last step until the first zero, ξ₁, is bracketed within 10⁻⁹. At n = 0 and n = 1 the equation has closed forms, θ = 1 − ξ²/6 and θ = sin ξ/ξ; at n = 5 the solution (1 + ξ²/3)^(−1/2) never reaches zero.

The caption under the solution reports the comparison, and the tabulated values of Chandrasekhar (1939) for n = 1.5, 3 and 4 are reproduced to the printed digits (ξ₁ = 3.65375, 6.89685 and 14.97155).

With M and R the constants follow: ρ_c = ρ̄ ξ₁/(3|θ′(ξ₁)|) with ρ̄ = 3M/(4πR³), P_c = GM²/(4πR⁴(n + 1)θ′(ξ₁)²), and the ideal gas with the mean molecular weight of a fully ionised mixture, μ = 4/(3 + 5X − Z), gives T_c = μ m_u P_c/(k ρ_c). Inside, ρ = ρ_c θ^n, P = P_c θ^(n+1) and T = T_c θ at r = R ξ/ξ₁. Radiation pressure, degeneracy and any change of composition with depth are left out: a polytrope has one μ throughout.

The energy generation is the pp chain plus the CNO cycle in the approximations of Kippenhahn, Weigert and Weiss (2012), eqs. (18.63) and (18.65), which fit the NACRE rates of Angulo et al. (1999): ε_pp = 2.57 × 10⁴ ψ f₁₁ g₁₁ ρ X² T₉^(−2/3) exp(−3.381 T₉^(−1/3)) with g₁₁ = 1 + 3.82 T₉ + 1.51 T₉² + 0.144 T₉³ − 0.0114 T₉⁴, and ε_CNO = 8.24 × 10²⁵ g₁₄ X_CNO X ρ T₉^(−2/3) exp(−15.231 T₉^(−1/3) − (T₉/0.8)²) with g₁₄ = 1 − 2.00 T₉ + 3.41 T₉² − 2.43 T₉³, both in erg g⁻¹ s⁻¹ with ρ in g cm⁻³ and T₉ = T/10⁹ K. The page takes ψ = 1 (the pp I branch alone), f₁₁ = 1 (no screening) and X_CNO = 0.66 Z, the share of C, N and O among the metals in the solar mixture of Asplund et al. (2009).

The luminosity is L = ∫ 4πr² ρ ε dr over the model, by Simpson's rule on 400 intervals.

The model luminosity rises steeply with n because the central temperature does, so one index matches the observed luminosity: the dot on the curve and the first readout.

Draw repeats that match with M, R and L sampled from the published uncertainties (a normal distribution; a split normal where the error bars are asymmetric; in log L where the paper quotes it that way) and collects the indices in the histogram with their median and 16th to 84th percentiles, until 2,000 draws or Pause; the small triangle under its axis marks the index at the central values. The nominal Sun and the MESA model carry no uncertainties, so there is nothing to draw for them.

For the Sun and the MESA preset the dots are the interior of a 1 M☉, [Fe/H] = 0 MESA model at 2.76 Gyr, the stored profile nearest the Sun's age (see the data notes), and the shaded band is its convective envelope, which begins at 0.74 R. The polytrope has a single index everywhere; the model's envelope follows a nearly adiabatic run of pressure with density that a core-fitted index cannot share, and its composition is not uniform: the hydrogen fraction, the dotted curve in the energy panel, falls from 0.70 at the surface to 0.50 at the centre, where the pp chain has been at work.

The two small panels show the solution θ(ξ) against the three exact cases, and the energy generation split into its two chains with the fraction of the luminosity produced inside each radius. Hover over the figures for values; a click on the luminosity curve sets n.