Lane–Emden equation

Consider a self-gravitating, spherically symmetric fluid in hydrostatic equilibrium. Mass is conserved and thus described by the continuity equation

${\displaystyle {\frac {dm}{dr}}=4\pi r^{2}\rho }$

where ${\displaystyle \rho }$ is a function of ${\displaystyle r}$. The equation of hydrostatic equilibrium is

${\displaystyle {\frac {1}{\rho }}{\frac {dP}{dr}}=-{\frac {Gm}{r^{2}}}}$

where ${\displaystyle m}$ is also a function of ${\displaystyle r}$. Differentiating again gives

${\displaystyle {\begin{aligned}{\frac {d}{dr}}\left({\frac {1}{\rho }}{\frac {dP}{dr}}\right)&={\frac {2Gm}{r^{3}}}-{\frac {G}{r^{2}}}{\frac {dm}{dr}}\\&=-{\frac {2}{\rho r}}{\frac {dP}{dr}}-4\pi G\rho \end{aligned}}}$

where the continuity equation has been used to replace the mass gradient. Multiplying both sides by ${\displaystyle r^{2}}$ and collecting the derivatives of ${\displaystyle P}$ on the left, one can write

${\displaystyle r^{2}{\frac {d}{dr}}\left({\frac {1}{\rho }}{\frac {dP}{dr}}\right)+{\frac {2r}{\rho }}{\frac {dP}{dr}}={\frac {d}{dr}}\left({\frac {r^{2}}{\rho }}{\frac {dP}{dr}}\right)=-4\pi Gr^{2}\rho }$

Dividing both sides by ${\displaystyle r^{2}}$ yields, in some sense, a dimensional form of the desired equation. If, in addition, we substitute for the polytropic equation of state with ${\displaystyle P=K\rho _{c}^{1+{\frac {1}{n}}}\theta ^{n+1}}$ and ${\displaystyle \rho =\rho _{c}\theta ^{n}}$, we have

${\displaystyle {\frac {1}{r^{2}}}{\frac {d}{dr}}\left(r^{2}K\rho _{c}^{\frac {1}{n}}(n+1){\frac {d\theta }{dr}}\right)=-4\pi G\rho _{c}\theta ^{n}}$

Gathering the constants and substituting ${\displaystyle r=\alpha \xi }$, where

${\displaystyle \alpha ^{2}=(n+1)K\rho _{c}^{{\frac {1}{n}}-1}/4\pi G,}$

we have the Lane–Emden equation,

${\displaystyle {\frac {1}{\xi ^{2}}}{\frac {d}{d\xi }}\left({\xi ^{2}{\frac {d\theta }{d\xi }}}\right)+\theta ^{n}=0}$

Equivalently, one can start with Poisson's equation,

${\displaystyle \nabla ^{2}\Phi ={\frac {1}{r^{2}}}{\frac {d}{dr}}\left(r^{2}{\frac {d\Phi }{dr}}\right)=4\pi G\rho }$

One can replace the gradient of the potential using the hydrostatic equilibrium, via

${\displaystyle {\frac {d\Phi }{dr}}=-{\frac {1}{\rho }}{\frac {dP}{dr}}}$

which again yields the dimensional form of the Lane–Emden equation.

If ${\displaystyle n=0}$, the equation becomes

${\displaystyle {\frac {1}{\xi ^{2}}}{\frac {d}{d\xi }}\left(\xi ^{2}{\frac {d\theta }{d\xi }}\right)+1=0}$

Re-arranging and integrating once gives

${\displaystyle \xi ^{2}{\frac {d\theta }{d\xi }}=C_{1}-{\frac {1}{3}}\xi ^{3}}$

Dividing both sides by ${\displaystyle \xi ^{2}}$ and integrating again gives

${\displaystyle \theta (\xi )=C_{0}-{\frac {C_{1}}{\xi }}-{\frac {1}{6}}\xi ^{2}}$

The boundary conditions ${\displaystyle \theta (0)=1}$ and ${\displaystyle \theta '(0)=0}$ imply that the constants of integration are ${\displaystyle C_{0}=1}$ and ${\displaystyle C_{1}=0}$. Therefore,

${\displaystyle \theta (\xi )=1-{\frac {1}{6}}\xi ^{2}}$

When ${\displaystyle n=1}$, the equation can be expanded in the form

${\displaystyle {\frac {d^{2}\theta }{d\xi ^{2}}}+{\frac {2}{\xi }}{\frac {d\theta }{d\xi }}+\theta =0}$

One assumes a power series solution:

${\displaystyle \theta (\xi )=\sum _{n=0}^{\infty }a_{n}\xi ^{n}}$

This leads to a recursive relationship for the expansion coefficients:

${\displaystyle a_{n+2}=-{\frac {a_{n}}{(n+3)(n+2)}}}$

This relation can be solved leading to the general solution:

${\displaystyle \theta (\xi )=a_{0}{\frac {\sin \xi }{\xi }}+a_{1}{\frac {\cos \xi }{\xi }}}$

The boundary condition for a physical polytrope demands that ${\displaystyle \theta (\xi )\rightarrow 1}$ as ${\displaystyle \xi \rightarrow 0}$. This requires that ${\displaystyle a_{0}=1,a_{1}=0}$, thus leading to the solution:

${\displaystyle \theta (\xi )={\frac {\sin \xi }{\xi }}}$

We start from with the Lane–Emden equation:

${\displaystyle {\frac {1}{\xi ^{2}}}{\frac {d}{d\xi }}\left(\xi ^{2}{\frac {d\theta }{d\xi }}\right)+\theta ^{5}=0}$

Rewriting for ${\displaystyle {\frac {d\theta }{d\xi }}}$ produces:

${\displaystyle {\frac {d\theta }{d\xi }}={\frac {1}{2}}\left(1+{\frac {\xi ^{2}}{3}}\right)^{3/2}{\frac {2\xi }{3}}={\frac {\xi ^{3}}{3\left[1+{\frac {\xi ^{2}}{3}}\right]^{3/2}}}}$

Differentiating with respect to ξ leads to:

${\displaystyle \theta ^{5}={\frac {\xi ^{2}}{\left[1+{\frac {\xi ^{2}}{3}}\right]^{3/2}}}+{\frac {3\xi ^{2}}{9\left[1+{\frac {\xi ^{2}}{3}}\right]^{5/2}}}={\frac {9}{9\left[1+{\frac {\xi ^{2}}{3}}\right]^{5/2}}}}$

Reduced, we come by:

${\displaystyle \theta ^{5}={\frac {1}{\left[1+{\frac {\xi ^{2}}{3}}\right]^{5/2}}}}$

Therefore, the Lane–Emden equation has the solution

${\displaystyle \theta (\xi )={\frac {1}{\sqrt {1+\xi ^{2}/3}}}}$

when ${\displaystyle n=5}$. This solution is finite in mass but infinite in radial extent, and therefore the complete polytrope does not represent a physical solution. Chandrasekhar believed for a long time that finding other solution for ${\displaystyle n=5}$ "is complicated and involves elliptic integrals".

In 1962, Sambhunath Srivastava found an explicit solution when ${\displaystyle n=5}$.[2] His solution is given by

${\displaystyle \theta ={\frac {\sin(\ln {\sqrt {\xi }})}{{\sqrt {\xi }}[3-2\sin ^{2}(\ln {\sqrt {\xi }})]}},}$

and from this solution, a family of solutions ${\displaystyle \theta (\xi )\rightarrow {\sqrt {A}}\,\theta (A\xi )}$ can be obtained using homology transformation. Since this solution does not satisfy the conditions at the origin (in fact, it is oscillatory with amplitudes growing indefinitely as the origin is approached), this solution can be used in composite stellar models.

It is known that if ${\displaystyle \theta (\xi )}$ is a solution of the Lane–Emden equation, then so is ${\displaystyle C^{2/n+1}\theta (C\xi )}$.[7] Solutions that are related in this way are called homologous; the process that transforms them is homology. If one chooses variables that are invariant to homology, then we can reduce the order of the Lane–Emden equation by one.

A variety of such variables exist. A suitable choice is

${\displaystyle U={\frac {d\log m}{d\log r}}={\frac {\xi ^{3}\theta ^{n}}{\varphi }}}$

and

${\displaystyle V={\frac {d\log P}{d\log r}}=(n+1){\frac {\varphi }{\xi \theta }}}$

We can differentiate the logarithms of these variables with respect to ${\displaystyle \xi }$, which gives

${\displaystyle {\frac {1}{U}}{\frac {dU}{d\xi }}={\frac {1}{\xi }}(3-n(n+1)^{-1}V-U)}$

and

${\displaystyle {\frac {1}{V}}{\frac {dV}{d\xi }}={\frac {1}{\xi }}(-1+U+(n+1)^{-1}V)}$.

Finally, we can divide these two equations to eliminate the dependence on ${\displaystyle \xi }$, which leaves

${\displaystyle {\frac {dV}{dU}}=-{\frac {V}{U}}\left({\frac {U+(n+1)^{-1}V-1}{U+n(n+1)^{-1}V-3}}\right)}$

This is now a single first-order equation.

The homology-invariant equation can be regarded as the autonomous pair of equations

${\displaystyle {\frac {dU}{d\log \xi }}=-U(U+n(n+1)^{-1}V-3)}$

and

${\displaystyle {\frac {dV}{d\log \xi }}=V(U+(n+1)^{-1}V-1).}$

The behaviour of solutions to these equations can be determined by linear stability analysis. The critical points of the equation (where ${\displaystyle dV/d\log \xi =dU/d\log \xi =0}$) and the eigenvalues and eigenvectors of the Jacobian matrix are tabulated below.[8]

${\displaystyle {\begin{array}{|l|l|l|}\hline {\text{Critical point}}&{\text{Eigenvalues}}&{\text{Eigenvectors}}\\\hline (0,0)&3,-1&(1,0),(0,1)\\(3,0)&-3,2&(1,0),(-3n,5+5n)\\(0,n+1)&1,3-n&(0,1),(2-n,1+n)\\\left({\dfrac {n-3}{n-1}},2{\dfrac {n+1}{n-1}}\right)&{\dfrac {n-5\pm \Delta _{n}}{2-2n}}&(1-n\mp \Delta _{n},4+4n)\\\hline \end{array}}}$