Polytrope

Two of the basic equations of stellar structure, hydrostatic equilibrium (recall P=ρgr and g=GM/r²):

dP = ρgdr = ρGMdr/r² ...(1)

and conservation of mass (recall M=ρV and dV=4πr²dr)

dM = 4π r² ρdr ...(2)

Rewrite (1) as M = (-r²/Gρ)dP/dR then (2) can be written

d/dr [-(r²/Gρ)dP/dr] = 4πr²ρ

d/dr [(r²/ρ)dP/dr] = -4πr²ρG

(1/r²)d/dr[(r²/ρ)dP/dr] = -4πGρ ...(3)

Which can be closed if we have, eg,
(remember the ratio of specific heats γ=c_p/c_v>1?)

P = Kρ^γ = Kρ^1+1/n ...(4a) and dP/dr=K(1+1/n)ρ^1/n ...(4b)

where n = 1/(γ-1) is called the polytropic index.
(Remember how γ=5/3 for an ideal monatomic gas?)

Eliminating P and dP/dr in (3) we have

(1/r²)d/dr[K (n+1)/n ρ^{1/n -1} r² dρ/dr] = -4πGρ ...(5)

or, since d/dr ρ^1/n = (1/n)ρ^{1/n -1} dρ/dr,

(1/r²)d/dr[K(n+1)r²dρ^1/n/dr] = -4πGρ ...(6)

We now have a simple OrdinaryDifferentialEquation with boundary conditions for a physical thing like a star

dP/dr=0 at the center and ρ=0 at the surface.

Now we introduce a dimensionless radius and density

r/R=ξ and ρ/ρ_c=θⁿ ie, put r=Rξ, dr=Rdξ, ρ= ρ_cθⁿ then (6) becomes

(1/ξ²) d/dξ [ξ² dθ/dξ] = -[(4πGR²)/K(n+1)]ρ_c^1-1/n θⁿ

Choose the scaling factor R=[K(n+1)/(4πG)]^½ρ_c^(1-n)/2n

and the equation becomes ξ^-2 d/dξ[ξ² dθ/dξ] = -θⁿ

which can be written as the Lane-Emden equation:

d²θ/dξ² + 2/ξ dθ/dξ + θⁿ = 0 ...(7)

for a polytropic model of index n. Given n the ODE can be integrated outwards with boundary conditions at the center ξ=0, θ=1, dθ/dξ=0.

For n<5 θ will fall to zero for some ξ=ξ_r which can be considered the surface of the polyreope. For n=1, 3 and 5 the solution is analytic:

n=0, θ = 1-ξ²/6,

n=1, θ = sinξ / ξ,

n=5, θ = (1+ξ²/3)^-1/2.

For other values numerical integration is used, rearrange the Land-Emden equation as

d²θ/dξ² = -2/ξ dθ/dξ - θⁿ
Step outwards in ξ_i+1=ξ+Δξ and evaluate θ(ξ_i+1) = θ(ξ)+dθ/dξ Δξ
where

dθ/dξ_i+1=dθ/dξ_i+d²θ/dξ²_i Δξ
and from te rearranged L-E equation we get

dθ/dξ_i+1=dθ/dξ_i-(2/ξ dθ/dξ_i+θⁿ)Δξ

Try Δξ=0.001 and to avoid blowup at ξ=0 for the first few steps use an expansion like

θ(ξ)=1-ξ²/6+nξ⁴/120-n(8n-5)ξ⁶/15120+....
or even more terms from HERE.

And HERE is an XLS starter kit... drag the columns down as far as neccesary to get θ=0 -- almost 7000 for an n=3 polytrope!

For a relativistic degenerate electron gas Pe goes as ρ^4/3 like an n=3 polytrope.
The n=3 case also applies to a star in which the gas pressure is a constant fraction of the radiative pressure (ie the star is in radiative equilibrium, Eddington's "standard model").

For a nonreletavistic degenerate electron gas Pe goes as ρ^5/3 like an n=1.5 polytrope.

The n=1 case is a self-gravitating isothermal gas which can be used to approximate the distribution of stars in a globular cluster.
²