solving the lane-emden equation

---

It is possible to solve the Lane-Emden equation analytically, i.e. to derive a simple mathematical expression for the dimensionless density as a function of dimensionless radius , for only three values of the polytropic index, n. These are

n = 0,     = 1 - (² / 6),

n = 1,     = sin / ,    and

n = 5,     = 1 / [1+(² / 3)]^0.5.

Solutions for all other values of n must be obtained numerically, i.e. using a computer to plot versus , as will now be described.

We start by expressing the Lane-Emden equation in the form

d² / d² = - [ (2 / ) d / d ] - ⁿ.

The numerical integration technique we will adopt consists of stepping outwards in radius from the centre of the star and evaluating the density at each radius. At each radius, the value of density, _i+1, is given by the value of density at the previous radius, _i , plus the amount the density changes by over the step. Denoting the step length as , we therefore have

_i+1 = _i + [ (d / d)_i+1].

The rate of change of density with radius, d / d , is an unknown in the above equation. To determine its value at each radius, we can use the same technique and write

(d / d)_i+1 = (d / d)_i + ( d² / d²).

The d² / d² term in the above equation is given by the rearranged form of the Lane-Emden equation given above, so we obtain

(d / d)_i+1 = (d / d)_i - ( [ (2 / _i) . (d / d)_i ] + _iⁿ ) . .

The numerical integration for a particular value of n can now proceed as follows. Starting at the centre, at which the values of , d / d and are known (because the boundary conditions tell us that d / d = 0 and = 1 at = 0), determine (d / d)_i+1. This value can then be used to determine _i+1. The radius is then incremented by adding to and the process repeated until the surface of the star is reached (when becomes negative). The results of such a numerical integration for five different values of n are shown in figure 22. Note that in this integration, was set to 0.001 and the starting value for was 0.00001 rather than 0 (in order to avoid an infinity in the calculations).

 

Figure 22:

Numerical solutions to the Lane-Emden equation for (left-to-right) n = 0, 1, 2, 3, 4 and 5.

Figure 22 shows that decreasing the polytropic index results in a stellar model in which the mass is more and more centrally condensed. In addition, it can be seen that the density of the n = 5 polytrope never falls below zero and, by inspection of the analytical solution given above, will do so at an infinite radius. The n = 5 polytrope (and n    5 polytropes in general) thus represents the non-physical case of an infinite star. For the n < 5 polytropes, the solution drops below zero at a finite value of , and hence the root (i.e. the radius of the polytrope, _R) can be determined using a linear interpolation between the points immediately before and after becomes negative. The roots for a range of polytropic indices are listed in table 2.

 

Table 2:

Numerical solutions to the Lane-Emden equation, showing the polytropic index n and the values of and d / d at the surface of the polytrope. It can be seen that the roots of the n = 0 and 1 solutions are in agreement with the roots of the analytical solutions listed above (which can be derived by setting =0 and rearranging for ).

n	R	- d / d = R
0	2.45	3.33 x 10-1 x R
1	3.14	1.01 x 10-1 x R
2	4.35	2.92 x 10-2 x R
3	6.90	6.14 x 10-3 x R
4	15.00	5.33 x 10-4 x R

How do these polytropic models compare to the results of a detailed solution of the equations of stellar structure? In order to make this comparison, which we shall do by comparing an n=3 polytropic model of the Sun (known as the Eddington Standard Model) with the so-called Standard Solar Model (SSM) [Bahcall 1998; Physics Letters B, 433, 1], we must convert the dimensionless radius, , and density, , to actual radius (in m) and density (in kg m^-3). We must also determine how the mass, pressure and temperature vary with radius. We do this as follows.

First, we must determine the value of the scale factor, . At the surface of the n=3 solar polytrope, where = 0, we can write

= R / _R ,

where R is the radius of the star (in this case, the radius of the Sun) and _R is the value of at the surface (i.e. the root of the Lane-Emden equation for n=3, as listed in table 2).

Next, we determine the mass as a function of radius. The rate of change of mass with radius is given by the equation of mass conservation

dM / dr = 4 r² .

By integrating and substituting r = and = _c  ⁿ into the above equation, we obtain

M = 4 r² dr = 4 ³ _c ^{2  n} d.



To evaluate this equation we need to rewrite the Lane-Emden equation in the form

² d / d = - ^{2 n} d.



Combining the last two equations, we obtain an equation for the mass of the star as a function of radius

M = - 4 ³ _c ² d / d.

The value of d / d as a function of is known from the numerical integration procedure outlined above. The value of _c, however, is not known. To determine it, we combine the equation for M above with a formula for the mean density of the Sun, i.e.

_av = 3M / 4r³ = - 3 _c (1 / ) d / d _{= R ,}

where the last term (between the vertical lines) means that it is evaluated at the surface of the star and can be calculated from the figures for n=3 given in table 2. We can therefore use this equation to determine _c, which in turn allows us to determine the variation of M with . This can be transformed to the variation of M with r using the formula r = . The result is plotted in the second panel from the top in figure 23.

 

Figure 23:

Comparison of numerical solution for an n = 3 polytrope of the Sun versus the Standard Solar Model.

It is now straightforward to determine the variation of density with radius using the equation

= _{c ⁿ.}

The result is plotted in the top panel of figure 23.

We can determine the variation of pressure with radius using the polytropic equation of state

P = K ^(n+1)/n = K _c ^{(n+1)/n n+1}.

To evaluate this equation we only require K, which we can obtain from the definition of scale length we adopted when deriving the Lane-Emden equation:

= [(n+1)K / 4 G _c^(n-1)/n] ^0.5.

The result is plotted in the third panel from the top in figure 23.

We can determine the variation of temperature with radius by equating the equation of state of an ideal gas with the polytropic equation of state

P = kT / m_H = K _c ^{(n+1)/n n+1},

which upon rearrangement and substitution of   = _{c ⁿ} gives

T = m_HK _c ^1/n / k.

The result is plotted in the bottom panel of figure 23. For comparison, we also plot the profiles of temperature, density, pressure and mass from the SSM, which is the most up-to-date solution to the equations of stellar structure currently available. It can be seen that the polytrope model does remarkably well, considering how simple the physics is - we have used only the mass and radius of the Sun and an assumption about the relationship between internal pressure and density as a function of radius. The agreement is particularly good at the core of the star, where the polytrope gives a central density of 7.65 x 10⁴ kg m^-3, a central pressure of 1.25 x 10¹⁶ N m^-2 and a central temperature of 1.18 x 10⁷ K, in comparison with the SSM values of 1.52 x 10⁵ kg m^-3, 2.34 x 10¹⁶ N m^-2 and 1.57 x 10⁷ K. Only in the outer regions of the Sun, where convection takes place, do the two solutions deviate significantly from one another.

---

©Vik Dhillon, 13th November 2012