Numerical solution of Emden–Chandrasekhar equation

In astrophysics, the Emden–Chandrasekhar equation is a dimensionless form of the Poisson equation for the density distribution of a spherically symmetric isothermal gas sphere subjected to its own gravitational force, named after Robert Emden and Subrahmanyan Chandrasekhar.^[1][2] The equation was first introduced by Robert Emden in 1907.^[3] The equation^[4] reads

where ${\displaystyle \xi }$ is the dimensionless radius and ${\displaystyle \psi }$ is the related to the density of the gas sphere as ${\displaystyle \rho =\rho _{c}e^{-\psi }}$, where ${\displaystyle \rho _{c}}$ is the density of the gas at the centre. The equation has no known explicit solution. If a polytropic fluid is used instead of an isothermal fluid, one obtains the Lane–Emden equation. The isothermal assumption is usually modeled to describe the core of a star. The equation is solved with the initial conditions,

${\displaystyle \psi =0,\quad {\frac {d\psi }{d\xi }}=0\quad {\text{at}}\quad \xi =0.}$

The equation appears in other branches of physics as well, for example the same equation appears in the Frank-Kamenetskii explosion theory for a spherical vessel. The relativistic version of this spherically symmetric isothermal model was studied by Subrahmanyan Chandrasekhar in 1972.^[5]

Derivation

For an isothermal gaseous star, the pressure ${\displaystyle p}$ is due to the kinetic pressure and radiation pressure

${\displaystyle p=\rho {\frac {k_{B}}{WH}}T+{\frac {4\sigma }{3c}}T^{4}}$

whe

• ${\displaystyle \rho }$ is the density
• ${\displaystyle k_{B}}$ is the Boltzmann constant
• ${\displaystyle W}$ is the mean molecular weight
• ${\displaystyle H}$ is the mass of the proton
• ${\displaystyle T}$ is the temperature of the star
• ${\displaystyle \sigma }$ is the Stefan–Boltzmann constant
• ${\displaystyle c}$ is the speed of light

The equation for equilibrium of the star requires a balance between the pressure force and gravitational force

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

where ${\displaystyle r}$ is the radius measured from the center and ${\displaystyle G}$ is the gravitational constant. The equation is re-written as

${\displaystyle {\frac {k_{B}T}{WH}}{\frac {1}{r^{2}}}{\frac {d}{dr}}\left(r^{2}{\frac {d\ln \rho }{dr}}\right)=-4\pi G\rho }$

Actual solution and asymptotic solution

Introducing the transformation

${\displaystyle \psi =\ln {\frac {\rho _{c}}{\rho }},\quad \xi =r\left({\frac {4\pi G\rho _{c}WH}{k_{B}T}}\right)^{1/2}}$

where ${\displaystyle \rho _{c}}$ is the central density of the star, leads to

${\displaystyle {\frac {1}{\xi ^{2}}}{\frac {d}{d\xi }}\left(\xi ^{2}{\frac {d\psi }{d\xi }}\right)=e^{-\psi }}$

The boundary conditions are

${\displaystyle \psi =0,\quad {\frac {d\psi }{d\xi }}=0\quad {\text{at}}\quad \xi =0}$

For ${\displaystyle \xi \ll 1}$, the solution goes like

${\displaystyle \psi ={\frac {\xi ^{2}}{6}}-{\frac {\xi ^{4}}{120}}+{\frac {\xi ^{6}}{1890}}+\cdots }$

Limitations of the model

Assuming isothermal sphere has some disadvantages. Though the density obtained as solution of this isothermal gas sphere decreases from the centre, it decreases too slowly to give a well-defined surface and finite mass for the sphere. It can be shown that, as ${\displaystyle \xi \gg 1}$,

${\displaystyle {\frac {\rho }{\rho _{c}}}=e^{-\psi }={\frac {2}{\xi ^{2}}}\left[1+{\frac {A}{\xi ^{1/2}}}\cos \left({\frac {\sqrt {7}}{2}}\ln \xi +\delta \right)+O(\xi ^{-1})\right]}$

where ${\displaystyle A}$ and ${\displaystyle \delta }$ are constants which will be obtained with numerical solution. This behavior of density gives rise to increase in mass with increase in radius. Thus, the model is usually valid to describe the core of the star, where the temperature is approximately constant.^[6]

Singular solution

Introducing the transformation ${\displaystyle x=1/\xi }$ transforms the equation to

${\displaystyle x^{4}{\frac {d^{2}\psi }{dx^{2}}}=e^{-\psi }}$

The equation has a singular solution given by

${\displaystyle e^{-\psi _{s}}=2x^{2},\quad {\text{or}}\quad -\psi _{s}=2\ln x+\ln 2}$

Therefore, a new variable can be introduced as ${\displaystyle -\psi =2\ln x+z}$, where the equation for ${\displaystyle z}$ can be derived,

${\displaystyle {\frac {d^{2}z}{dt^{2}}}-{\frac {dz}{dt}}+e^{z}-2=0,\quad {\text{where}}\quad t=\ln x}$

This equation can be reduced to first order by introducing

${\displaystyle y={\frac {dz}{dt}}=\xi {\frac {d\psi }{d\xi }}-2}$

then we have

${\displaystyle y{\frac {dy}{dz}}-y+e^{z}-2=0}$

Reduction

There is another reduction due to Edward Arthur Milne. Let us define

${\displaystyle u={\frac {\xi e^{-\psi }}{d\psi /d\xi }},\quad v=\xi {\frac {d\psi }{d\xi }}}$

then

${\displaystyle {\frac {u}{v}}{\frac {dv}{du}}=-{\frac {u-1}{u+v-3}}}$

Properties

• If ${\displaystyle \psi (\xi )}$ is a solution to Emden–Chandrasekhar equation, then ${\displaystyle \psi (A\xi )-2\ln A}$ is also a solution of the equation, where ${\displaystyle A}$ is an arbitrary constant.
• The solutions of the Emden–Chandrasekhar equation which are finite at the origin have necessarily ${\displaystyle d\psi /d\xi =0}$ at ${\displaystyle \xi =0}$

See also

• Lane–Emden equation
• Frank-Kamenetskii theory
• Chandrasekhar's white dwarf equation

References