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.
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.
The mass can be obtained as follows. First, multiply the L-E equation by x2. This yields
| (1) |
Now note that the internal mass as a function of radius is given by
| (2) |
Problem 5.2 Substitute Eqn 5.2 into 5.1 to show that the mass is given by
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
|
This gives us the central density from the known mean density of the star. The density as a function of radius is then just
| (5) |
From the definition of the pressure equation
| (6) |
where g = 1 + [ 1/n], we obtain
| (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
| (8) |
From this one obtains K,
| (9) |
The temperature, T comes from the pressure P and the assumed equation of state for the plasma
| (10) |
where N is the number density of particles. For a plasma of mean molecular weight m, this is
| (11) |
The central temperature, Tc, is then
| (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)).
|
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.
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
| (13) |
the value of the derivative at the root at x = xR is denoted f¢R (see figure 5.1) and is
| (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.
|
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 |
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
| (15) |
If the gas obeys the perfect gas law, then we have
| (16) |
and so for a perfect gas,
| (17) |
In differential form this is just
| (18) |
The first law of thermodynamics states that the internal energy change in the gas is given by
| (19) |
Combining these results we get that
| (20) |
which may be integrated easily to yield
| (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
| (22) |
This is a polytrope of index n=1.5.
More generally, the energy of the gas can take the form
| (23) |
where b is the number of degrees-of-freedom in the internal state of the gas. This leads to
| (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.
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
| (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.
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.
| (26) |
and
| (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
| (28) |
and the ratio is
| (29) |
Eliminating T, we get
| (30) |
This is a polytrope with g = 4/3, or alternatively n=3, and is called the Eddington Standard Model.
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
| (31) |
and thus the scale length of the polytrope, a is
| (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.
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?
|
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