Lecture 5 : Polytropes concluded

Lecture 5 : Polytropes concluded

1  Properties of Polytropes

Recall from the previous lecture, that there are three known analytical solutions to the Lane-Emden equation. For n=5, the Lane-Emden equation has no roots, and thus represents the non-physical case of an infinite star. The two other cases are n=0 and n=1.





Problem 5.1 Show that for n=0 the root is x0=Ö6 and for n=1, x1=p.





We will now solve the Lane-Emden equation, for a star of mass M and radius R.

1.1  Scaling length, a

The radius, R is just the first root of x and is denoted xR, which can be obtained numerically or analytically, depending on n. The root gives us the scaling length a = R/xR.

1.2  Mass and density

The mass can be obtained as follows. First, multiply the L-E equation by x2. This yields


x2  df

dx
= - ó
õ
x

0 
x2 fn dx.
(1)

Now note that the internal mass as a function of radius is given by


M(r) = ó
õ
r

0 
4 pr2 rdr = 4 pa3 rc ó
õ
x

0 
x2 fn dx.
(2)





Problem 5.2 Substitute Eqn 5.2 into 5.1 to show that the mass is given by


M = - 4 pa3 rc ê
ê
x2  df

dx
ê
ê


x = xR 
(3)

where rc is the central density





The mean density [`(r)] is just the mass M, divided by the volume V = [ 4/3] pR3, or [`(r)] = [ 3M/(4 pR3)]





Problem 5.3 Show that the central density is related to the mean density by


-
r
 
= - 3 rc ê
ê
 1

x
 df

dx
ê
ê


x = xR 
(4)





This gives us the central density from the known mean density of the star. The density as a function of radius is then just


r(x) = rc fn(x).
(5)

1.3  Pressure and Temperature

From the definition of the pressure equation


P = K rg
(6)

where g = 1 + [ 1/n], we obtain


P = K rc1 + [ 1/n] fn+1.
(7)

We now only require K to get the density as a function of radius in the star. We can get this from our definition of the scaling length a, or we can integrate the equation of hydrostatic equilibrium to obtain the pressure directly. The central pressure Pc is found to be


Pc =  M2 G

R4
é
ë
4 p(n+1) ê
ê
 df

dx
ê
ê
2

x = xR 
ù
û
-1

 
.
(8)

From this one obtains K,


K =  Pc

rcg
.
(9)

The temperature, T comes from the pressure P and the assumed equation of state for the plasma


P = NkT
(10)

where N is the number density of particles. For a plasma of mean molecular weight m, this is


P =  rk T

mmH
.
(11)

The central temperature, Tc, is then


Tc =  Pc mmp

k rc
.
(12)





Problem 5.4 Find an expression for the temperature as a function of x in terms of known quantites (e.g. mass, M, radius R, central density Pc, central temperature, Tc and the solution to the L-E equation, f(x)).





1.4  Evaluation of quantities at the stellar surface

Several equations above depend on the quantity df/dx evaluated at the surface (i.e. x = r/a). We would like to evaluate these numerically using our computer program. Consider the sequence of evaluated points along the solution to the L-E equation illustrated in figure 5.1.

roots.gif
Figure 1: Numerical solution to the L-E equation for n=1 plotted close to the root f(x) = 0. Lower panel: The points before and after the solution crosses through f = 0 are shown as (xn,fn) and (xn+1,fn+1). The values of the derivative of the density function f¢ = df/dx, at the same points before and after the solution are shown in the upper panel

Denoting the solutions before and after the sequence crosses through f = 0, as (xn,fn) and (xn+1,fn+1), we can make a linear interpolation between these points to estimate the root of the equation xR


xR = xn+1 + (xn+1-xn)  fn+1

(fn-fn+1)
.
(13)

the value of the derivative at the root at x = xR is denoted f¢R (see figure 5.1) and is


f¢R = f¢n+1 -  (xn+1-xR)

(xn+1-xn)
(f¢n+1-f¢n)
(14)





Problem 5.5 A new version of the fortran program laneem.f is listed in the appendix. This version prints out the root of the L-E equation and the value of f¢ at the root. Modify this program to print out -x2f¢ as well. Compute this quantity at the stellar surface, and compare your results with the data in Table 5.1.





Table 1: Numerical solutions to the L-E equation, showing the index n, first root xR, a quantity related to the mass (see text) and the ratio of the central density rc to the mean density [`(r)].
n xR-xR2 f¢R rc/[`(r)]
0 2.45 4.90 1.00
1 3.14 3.14 3.30
2 4.35 2.41 11.40
3 6.90 2.02 54.2
4 14.96 1.80 622.4

2  Polytrope for a perfect gas

When a small element of gas is compressed so that no heat is conducted into or out of the element, the change is described as adiabatic.

If the gas element is fully ionised gas and contains N particles, the total energy E is related to the temperature T by the Boltzmann constant k


E =  3

2
N k T.
(15)

If the gas obeys the perfect gas law, then we have


P =  N k T

V
(16)

and so for a perfect gas,


E =  3

2
PV.
(17)

In differential form this is just


dE =  3

2
(V dP + P dV).
(18)

The first law of thermodynamics states that the internal energy change in the gas is given by


dE = -P dV.
(19)

Combining these results we get that


 dP

P
= -  5

3
 dV

V
(20)

which may be integrated easily to yield


P µ V-5/3.
(21)

The mass of the gas element remains unchanged (the number of particles N is constant), and so the density is just inversely proportional to the volume. Thus


P µ r5/3.
(22)

This is a polytrope of index n=1.5.

More generally, the energy of the gas can take the form


E =  b

2
N k T
(23)

where b is the number of degrees-of-freedom in the internal state of the gas. This leads to


P µ r(1+2/b).
(24)

The polytropic index thus increases with the number of degrees-of-freedom. A limiting case is the isothermal gas, in which the temperature T is everywhere constant, so that P = rT µ r. The isothermal gas thus corresponds to (1+2/b) = 1, or b = ¥. For such a gas sphere, the polytropic index is n = ¥ >> 5. Recall that for all polytropes with n > 5, the star is unbounded (infinite mass and radius), so it is clear that stars cannot be isothermal spheres.

3  Polytrope for convective stars

In the context of a star, convection is the process by which heat is transferred by the motion of gas cells of lesser or greater density than the local density.

Convection is an efficient method of transporting heat out of the star, and is typically much more efficient than conduction. As a result, the heat energy within a rising or falling element does not change very much during the physical motion, and can be considered to take place adiabatically. As a result of the mixing of gas elements due to turbulence, the gas may be all on the same adiabat, which is to say that the equation of state satisfies


P µ r(1+2/b)
(25)

and all fluid elements have the same constant of proportionality.

Convection is particularly important in low mass stars, about 1/3 of the mass of the Sun. The Sun itself has a convective shell occupying the outer » 30% of its radius, while the inner part is radiative. We will look further into the balance of energy transport by radiation and convection in later lectures. As a result of the Sun not being fully convective, it is better described by a polytrope of n=3, rather than by n=3/2.

4  Eddington Standard Model

Consider that the internal pressure of the star comes from two sources, the gas pressure Pg and the photon (radiation) pressure, PR. As we shall see in the next lecture, the radiation pressure is proportional to the temperature T4, while the gas pressure is proportional to T.


PR =  aT4

3
(26)

and


Pg =  rk T

mmH
.
(27)

If we assume that the ratio of gas to radiation pressure is uniform throughout the star, we can express the total pressure P as


P = PR + Pg = (1-b)P + bP
(28)

and the ratio is


PR : Pg = 1-b: b.
(29)

Eliminating T, we get


P = (  k

mmH
)4/3 (  3(1-b)

ab4
)1/3r4/3.
(30)

This is a polytrope with g = 4/3, or alternatively n=3, and is called the Eddington Standard Model.

5  Solar polytrope

le.gif
Figure 2: Comparison of numerical solution for an n=3 polytrope of the Sun versus the Standard Solar Model (SSM). The SSM is the best available model for the Solar interior taking into account much more physics than we have looked at so far. The agreement is quite good, with the biggest disagreements being in the outer regions, where convection plays a significant role.

From the above arguments, we have good reason to model the Sun as an n=3 polytrope. Using the fortran program given below, one can solve the L-E equation for f = f(x). Various quantities at the surface of the star are listed in Table 5.1. For n=3, we have that xR = 6.90, -xRf¢x = xR = 2.02 and thus f¢x = xR = -0.042. We have, for the Sun, MO = 1.99 ×1030 kg and RO = 6.96 ×1010 cm.

We thus obtain


a3 rc = -  MO

4 px2R f¢R
= 8 ×1031 g
(31)

and thus the scale length of the polytrope, a is


a = 1.009 ×1010 cm.
(32)

Table 2 shows further quantities : the central density, rc, the mean density [`(r)], the central pressure, Pc and the central temperature Tc (using the ideal gas law) and finally the constant of the pressure law, K (i.e. P = K r4/3). The profiles of temperature, density and pressure are shown in figure 5.2 compared to the ``Standard Solar Model'', due to Bahcall (Physics Letters B, 433, 1, 1998). 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 beween internal pressure and density as a function of radius.

Table 2: Properties of an n=3 polytropic model of the Sun
rc 80 g cm-3
[`(r)] 0.02 rc » 1.4 g cm-3
Pc 1.2 ×1017  dyne cm-2
Tc 12 ×106  K
K 3.8 ×1014





Problem 5.6 Extra credit: Get the standard solar model from
http://www.sns.ias.edu/~ jnb/SNdata/solarmodels.html,
and use the output from the fortran program along with GNUPLOT to make plots similar to figure 5.2. Try fitting polytropes other than n=3. Can one make a better polytrope model that the one shown in figure 5.2?





References

[]
J. Bahcall, Standard Solar Model. http://www.sns.ias.edu/~ jnb/SNdata/solarmodels.html

[]
R. Bowers and T. Deeming. 1984. Astrophysics I, Stars, Jones and Bartlett Publishers.

[]
Press, W. Introduction to Astrophysics, Chapter 5, 1997, http://cfatab.harvard.edu/whp/ay45toc.html

Appendix

The following program is a modification of last lecture's laneem.f. It computes various values at the surface of the star (i.e. the root of the L-E equation) by using linear interpolation between the final two points (the first of which is before the root, and the second after it). See figure 5.1.

        program lane_emden_roots

c**************************************************************
c* program to solve the Lane-Emden equation numerically       *
c*                                                            *
c* the equation to solve is                                   *
c* d^2x/dt^2 = - ((2/t)*(dx/dt) + x**n)                       *
c*                                                            *
c* where x is the dimensionless density function (psi)        *
c* and t is the dimensionless length (xi)                     *
c* initial conditions are      x = 1                          *
c*                         dx/dt = 0                          *
c*                                                            *
c* and the power index is n.                                  *
c* analytical solutions are known for n=0, 1, 5.              *
c*                                                            * 
c* Physical quantities are computed at the root, where x=0    *
c*                                                            *
c* CF. Oct 1999.                                              *
c**************************************************************

        real x, dxdt, n, dt, xlast, dxdtlast
        character*50 name

* index n of the polytrope
        
        write(6,*) 'Enter n, polytrope index '
        read(5,*) n
        
        write(6,*) 'Enter name of results file'
        read(5,*) name

* open the file for output of results

        open(20,file=name,status='new')

* variables are 
* dimensionless length : t 
* dimensionless density function : x

* specify initial values 

        x    = 1.0   
        dxdt = 0.0
        t    = 0.0001

* "radius step" is in length variable, t
* the smaller the step the more accurate the calculation
* and the longer it takes!

        dt = 0.001

* now the MAIN program loop

        do 100 i = 1, 10000

* calculate the next dxdt                  
          dxdt = dxdt - ( 2.0*dxdt/t + x**n )*dt
        
* calculate the next x    
          x = x + dxdt*dt 

* increment the dimensionless radius (move outward in the star)
          t = t + dt
        
* write out the results

          write(20,*) t,x,dxdt

* check if the solution has become negative yet  
          if (x.lt.0.0) goto 999

* save the current values of x and t, for determining the root later
          xlast = x
          dxdtlast = dxdt

100     continue

999     continue

* calculate the position of the root, using the current x (which 
* is the first value < 0, and the previous x (stored as xlast) 
* which is the last calculated value which was still > 0

        t1 = t + (x/(xlast-x))*dt
        dphidt1 = dxdt - (t-t1)*(dxdt-dxdtlast)/dt

        write(6,*) 'first root of phi at xi=',t1 
        write(6,*) 'at this point dphi/dt = ',dphidt1

        end





File translated from TEX by TTHgold, version 3.11.
On 13 Apr 2003, 11:03.