EMDEN'S EQUATION

by Reinaldo Baretti Machín  

  Free counter and web stats

 

References :

1. Introduction to Nonlinear Differential and Integral Equations by Harold T. Davis ,  page  371 .

2.  http://mathworld.wolfram.com/Lane-EmdenDifferentialEquation.html

3.  http://www.jgiesen.de/astro/standardmodel/index.htm

4.  Polytropic Stars ,   http://www.astrophysicsspectator.com/topics/stars/Polytrope.html

 

 

ABSTRACT:  Emden's equation is solved numerically by finite differences. Its physical application to stars  is reviewed within the context of classical gravitation and thermodynamics.

 

 

 

Part I - numerical solution

 

Emden's equation has the form

   

                                      d2 y/dx2 +(2/x) dy/dx = - k yn                                         .  (1)

 

The physical applications of most interest, to be explained in part II,  impose initial conditions y(0)=1 ,  (dy/dx)0 = 0.

 

Assume for the moment that y ~ dimensionless , x ~ length ~ L . Then the factor  k ~ 1/L2 . Also (2/x) ~ 1 / L ,thus the factor 2 is dimensionless.  As given in (1) there is one scale of length    Lscale = 1/k1/2 .

Let k= 1. Lscale = 1/k1/2 = 1  and we should pick   ∆x << Lscale =1 .

 

The finite difference solution to (1) is  

 

                            yn =  2yn-1 - yn-2  + ( ∆x  )2 { -(2/xn-1)( yn-1 -yn-2 )/ ∆x   - k(yn-1 ) n  }  . (2)

 

Eq(2) is solved with the Fortran code given below for some values of the power n ;  n=0 , n=1 and n=2.

Physical considerations have to be brought to select the appropriate values of n or what is related to it ,the coefficient  γ.

                                                     

 

Fig. 1    Solution corresponding to n=0.

Fig 2.    Solution corresponding to n=1.

 

Fig. 3    Solution corresponding to n=2.

FORTRAN code

c Emden's equation HT Davis NL DE page 371
c d^2phi/dr^2+(2/r)d phi/dr= -a^2*phi^n Phi~gravitational potential
c Ic phi(0)=1 . phi(dr)=1.
data a,nstep ,n ,rf/1.,1000,1,8./
phi0=1.
phi1=1.
dr=rf/float(nstep)
kp=int(float(nstep)/80.)
kount=kp
print 100,0.,phi0
do 10 i=2,nstep
r=dr*float(i)
phi2=2.*phi1-phi0 +dr**2*(-(2./(r-dr))*(phi1-phi0)/dr
$-a**2*phi1**n )
if(i.eq.kount)then
print 100 ,r , phi2
kount=kount+kp
endif
phi0=phi1
phi1=phi2
10 continue
100 format('r,phi=',2(3x,e10.3))
stop
end
 

Part II - Derivation of Emden's equation

 

Fig. 4  Classical model of a star where the hydrostatic gravitational pressure is balanced by the kinetic pressure.

To begin there will be three variables pressure (P) ,  gravitational potential ( Φ)  and  density ρ . Three equations are needed to relate

them.

The gravitational  field  g(r) ~ m/s2 is defined by the negative of the gradient of the gravitational potential Φ(r) ,

                                            g(r) =  - dΦ /dr                                                             .      (3)

The units of Φ  ~ (m/s)2 .

At a distance (r) from the center of the star the hydrostatic  pressure decreases as we go from r to r +dr  by

        dP =  -weight of shell of width dr/ area = - g(r) ρ(r) 4 π r2 dr /(4 π r2 )

             = - g(r) ρ(r)  dr  = ρ(r)  d Φ     ~ Pascals ~ N/m2                                              .  (4)

Poisson equation is

                      ∂2Φ/∂r2 +     (2/r ) (∂Φ /∂r) = -4π G ρ(r)   ~ 1/s2  (G=6.67E-11 Nm2 /kg2 )  . (5)

The third equation is the adiabatic law from thermodynamics . The star as a whole on a "short" time span is not exchanging heat with the universe. No radiation pressure is taken into account. Of course this is only an approximation.

From thermodynamics in an adiabatic process pressure and volume satisfy

                                      P V γ  = constant    , where  the adiabatic index in general has a range    4/3  ≤  γ ≤ 5/3   ( see ref. 4) .

The lower value 5/3 is used for the case of a fully ionized plasma and 4/3 is used when the pressure is supported in part by radiation.

But       V= mass/ ρ           hence  ,

                                               P = K ργ                                        .   (6)

K is an empirical constant and so is γ . The choice of γ may not be unique. Notice that  if pressure and density are known at some radial distance  r1  then the constant K can be  determined  ;

                                            K = P1 / ρ(r1 )                              .      (7)

Using (6)  we get       dP =  γ K ργ-1 dρ  and equating to (4)  gives

                                                      γ K ργ-1 dρ  = ρ dΦ                    . (8)

Integration gives     

                                                     Φ(r) = [γ/(γ-1)] K ρ (r) γ-1  + C   .  (9)

The zero of potential is chosen at the surface r=R ,where ρ(R) =0  , so Φ(R) ≡ 0  + C   , hence C=0.  Φ(r) is a positive quantity.  It is not the usual choice, recall that for a point particle Φ(r)= -Gm/r , the potential is zero at infinity. Instead of  γ the parameter n is introduced  by  , γ = 1 +1/n  .

Solving for the density one gets

                                                          ρ(r) = [ (n+1)K ]-n Φ(r) n     . (10)

Substitution in (5) leads to an equation containing only the potential Φ

                         ∂2Φ/ ∂r2  +  (2/r ) (∂Φ /∂r) = -4π G [ (n+1)K ]-n Φ(r) n  = - a2 Φ(r) n  ,   (11)

where  a2 = 4π G [ (n+1)K ]-n . The dimensions of the constant a  will be dependent on the power n (called the polytropic index)  which in turn is a function of the specific heats ratio γ.

To retrieve Emden's equation from (11), introduce Φ(r) = Φ0 y   and   r = x/( a Φ0n/2 ) ≡k x ,  where small  k=1/(a Φ0n/2 ) . The

small k has units of length.  Then Φ0 =  Φ(r=0) yielding     

                                                               d2 y/dx2 +(2/x) dy/dx = -  yn       .               (12)   

The values of y of physical  interest  are those between 0 and 1.  For at x=0 , r=0 ; y=1 and Φ= Φ0 . On the other hand at     r=R     

Φ(R) =0     and y=0   .   

Let x0 the first root of the solution of (12) -for a particular index n-  then the radius of the sphere is

                                                         R = x0 /  (a Φ0n/2 )                                   .  (13) 

Ata radius r            

                                                            dΦ/dr = -GM(r) /r2      .                              (14)   

Or in terms of x and y and  Φ0             

                                                       (Φ0 /k) (dy/dx) x  = - GM(x)/(kx)2 .                   (15)

Eq(15) serves to calculate  Φ0 , the potential at the center of the star. Let x= x0 ( the first root) , then M(x0 ) = Mtotal and

the derivative (dy/dx)x0  can be read from the numerical solution. notice that k also depends on a power of Φ0 . It turns out that

                                Φ0 = {  a [ (-GM/x02 ) / (dy/dx)x0 ] }2/(2-n)               .              (16)

Three applications of Emden's equation will be discussed.  It is used to determine  the density ρ(r) and central pressure P0 and temperature T0 .The basic data are the coefficient γ the radius R and total mass M.

At a distance r from the center the density is

                                    ρ(r) = 3 M(r)/(4π r3 ) =  (3/4π ) (1/G) {-(1/r)dΦ/dr }  ,        (17-a)

                                          = {3Φ0 /(4πk2 G)} ( -(1/x) dy/dx )                   .           (17-b)  

 

At the center                                       ρ(0)= {3Φ0 /(4πk2 G)} { -(1/x) dy/dx )x=0 }       ,  (18)

but  -(1/x) (dy/dx )0 = 1/3 for all solutions as seen from our numerical results given below. Thus the central density is

                                                                ρ(0)=  Φ0 /(4πk2 G)}                               .    (19)

Numerical results , see below, show that (1/x)(dy/dx)x=0 = -1/3.  

Another way of getting ρ(0) is the following . From (17) we get the mean density as a function of the derivative at the surface r=R or equivalently  x=x0 ,

                               ρmean = 3 M(R)/(4π R3 ) =    {3Φ0 /(4πk2 G)} ( -(1/x) dy/dx ) x=x0    .  (20)

Dividing  eq. (19) by eq. (20) gives ,

                                             ρ(0)/ρmean  = {  (-(3/x)dy/dx ) x=x0 } -1                            . (21)

                                          

  y(0)=1.,(dy/dx)i=0.n=0
r,phi,D= 0.120E-01 0.100E+01 -0.333E+00
r,phi,D= 0.160E-01 0.100E+01 -0.333E+00
r,phi,D= 0.200E-01 0.100E+01 -0.333E+00
r,phi,D= 0.240E-01 0.100E+01 -0.333E+00
r,phi,D= 0.280E-01 0.100E+01 -0.333E+00
r,phi,D= 0.320E-01 0.100E+01 -0.333E+00
r,phi,D= 0.360E-01 0.100E+01 -0.333E+00
r,phi,D= 0.400E-01 0.100E+01 -0.333E+00
    

 y(0)=1.,(dy/dx)i=0.n=1
r,phi,D= 0.120E-01 0.100E+01 -0.333E+00
r,phi,D= 0.160E-01 0.100E+01 -0.333E+00
r,phi,D= 0.200E-01 0.100E+01 -0.333E+00
r,phi,D= 0.240E-01 0.100E+01 -0.333E+00
r,phi,D= 0.280E-01 0.100E+01 -0.333E+00
r,phi,D= 0.320E-01 0.100E+01 -0.333E+00
r,phi,D= 0.360E-01 0.100E+01 -0.333E+00
r,phi,D= 0.400E-01 0.100E+01 -0.333E+00

 

 y(0)=1.,(dy/dx)i=0.n= 3
r,phi,D= 0.120E-01 0.100E+01 -0.333E+00
r,phi,D= 0.160E-01 0.100E+01 -0.333E+00
r,phi,D= 0.200E-01 0.100E+01 -0.333E+00
r,phi,D= 0.240E-01 0.100E+01 -0.333E+00
r,phi,D= 0.280E-01 0.100E+01 -0.333E+00
r,phi,D= 0.320E-01 0.100E+01 -0.333E+00
r,phi,D= 0.360E-01 0.100E+01 -0.333E+00
r,phi,D= 0.400E-01 0.100E+01 -0.333E+00

 

We proceed to get formulas for the central pressure P0  and the central temperature T0 . From (10) write

                                                             Φ0 = (n+1) K ρ0(1/n)                         (22 )

and from (6)                                           P0 = K ρ0 (1/n + 1)                        .    (23)

Dividing eq (21) by (2) gives the central pressure formula

                                                     P0  = ρ0 Φ0 /(n+1)                             .     (24)

The central temperature is essentially related to P0 by the equation of state,

                                      P0 =  (kBoltz. / mhydrogen ) ρ0 T0

                                          =  ρ0 Φ0 /(n+1)               ,

thus                           

                                           T0 = ( mhydrogen / kBoltz.)  Φ0 /(n+1)             .        (25)

In SI units ( mhydrogen / kBoltz.)  = 1.21 E-4  kelvin -(s/m)2   . If we assume that radiation pressure contributes by the same amount as the hydrostatic pressure , the constant is multiplied by a factor of  2 and we rewrite T0 as

                                      T0 = 2 ( mhydrogen / kBoltz.)  Φ0 /(n+1)

                                           = 2.42 E-4  Φ0 /(n+1)  ~ kelvin   when Φ0 ~ (m/s)2 .   (26)

Example I :   Example taken from page 376 reference 1. 

Th bright component in Capella has M=8.60 E 30 kg  , R =9.55E 9 m and it assumed that γ =4/3 .

Find the central values of  density ρ0  , gravitational potential  Φ0 , pressure P0  and temperature T0 .

The plot of y vs  x for n=3 is given  next in fig 5.

                               Fig 5.  Solution corresponding to n=3 , γ = 4/3, x0 = 6.92 .


Fig 6. Plot of   { -(3/x) y' }^-1 needed to calculate ρ(0) from (21)   . At x0 its value is about 55; x0 = 6.92 .

 

To calculate ρ0 , Φ0 , P0 and T0 we used FORTRAN code  #2 .

There is quite an  entanglement  of formulas , so there may be more than one way to go about it.

a) calculation of ρmean  :

                        ρmean = 3M/(4π R3 ) = 2.36 kg/m3

b) calculation of   ρ0 

From (21)  and using the values of the derivatives in Fig. 6

                                 ρ(0) = ρmean  {  (-(3/x)dy/dx ) x=x0 } -1

                                        = (2.36) (55) = 130 kg/m3

c)  calculation of   k  ( don't confuse with capital K)

From the definition above, see the lines before eq(12) ,      r = x/( a Φ0n/2 ) ≡k x . So

                                      k =  R/x0  ~ meter

                                         = 9.55E 9/6.92 = 1.38E9 m

d) calculation of Φ0 , from eq (19) 

                                                     Φ0 = 4π k2 G ρ0

                                                          = 4π ( 1.38 E9 )2 (6.67E-11) (130)

                                                          = 2.07E 11 (m/s)2 

e) calculation of P0 , from (24)

                                                      P0  = ρ0 Φ0 /(n+1) 

                                                          = 6.73E+12 pascal ( N/m2 )    . 

f) calculation of T0

                                                       T0  = 2.42 E-4  Φ0 /(n+1)

                                                            = 1.25E+7 kelvin    .

 

Results using FORTRAN code # 2.
 rhozero ~(kg/m^3)= 129.646973
phizero~(m/s)^2= 2.06962983E+011
pzero (pascals)= 6.70798132E+012
tzero~ kelvin= 12521167.
 

Example II :   Example taken from page 377 reference 1. 

Data for the sun is M= 1.99 E 30 kg , R=6.96 E 8 m  ,assume γ = 7/5 ,thus  n=1/(γ-1) = 5/2 .

We follow the same steps as in example I.  But we have to solve eq.(2) with n= 5/2 to find the expression

 -(3/x)dy/dx ) x=x0 .

We have near x0 , that  -(3/x)dy/dx ) x=x0 = 23.3   .

r,phi,D= 0.502E+01 0.282E-01 0.193E+02
r,phi,D= 0.508E+01 0.225E-01 0.200E+02
r,phi,D= 0.515E+01 0.170E-01 0.208E+02
r,phi,D= 0.521E+01 0.116E-01 0.216E+02
r,phi,D= 0.528E+01 0.641E-02 0.225E+02
r,phi,D= 0.535E+01 0.131E-02 0.233E+02
r,phi,D= 0.541E+01 0.1#RE+01 0.1#RE+01

Using Fortran code #2 with the new data,we obtain

rhomean,rhozero(kg/m^3)= 0.141E+04 0.328E+05
phizero~(m/s)^2 ,pzero (pascals),tzero(kelvin)= 0.466E+12 0.437E+16 0.3
22E+08

That is ,

   ρmean = 1.41 E3 kg/m3  ( higher or comparable to all planets beyond Mars)   ,  ρ0 = 3.28E4 kg/m3 .

Φ0 = 4.66 E11 (m/s)2 ,   P0 = 4.37E15 Pa    , T0 = 3.22 E 7 kelvin    .

 

 

Fortran code #2

c example of the calculation of central density , pressure and
c temperature of a star from Emden's equation
c see Nonlinear DE H. T. Davis Dover publications page 371
real K
data amass,radius,gama/8.60E30,9.55E9 ,1.33333 /
data xzero/6.92 /
data G/6.67E-11/
pi=2.*asin(1.)
ainvDx= 55.
rhom=amass/(4.*pi*radius**3/3)
an= 1./(gama-1.)
rhozero=rhom*ainvDx
K=radius/xzero
phizero=4.*pi*k**2*G*rhozero
Pzero=rhozero*phizero/(an+1.)
Tzero=2.42E-4*phizero/(an+1.)
print*,'rhozero (kg/m^3)=',rhozero
print*,'phizero~(m/s)^2=',phizero
print*,'pzero (pascals)=',pzero
print*,'tzero~ kelvin=',tzero
stop
end




                     

Fortran code (two examples with different n values)

c Emden's equation HT Davis NL DE page 371
c d^2phi/dr^2+(2/r)d phi/dr= -a^2*phi^n Phi~gravitaional potential
c Ic phi(0)=1 . phi(dr)=1.
real n
data a,nstep ,n ,rf/1.,1000,2.5,5.5/
phi0=1.
phi1=1.
dr=rf/float(nstep)
kp=int(float(nstep)/80.)
kount=kp
c print *,'y(0)=1.,(dy/dx)i=0.n=', n
print 100,0.,phi0,1.
do 10 i=2,nstep
r=dr*float(i)
phi2=2.*phi1-phi0 +dr**2*(-(2./(r-dr))*(phi1-phi0)/dr
$-a**2*phi1**n )
c if(i.gt.2.and.i.le.10)then
if(i.eq.kount)then
D=1./((-3./r)*(phi2-phi1)/dr)
print 100 ,r , phi2 ,D
kount=kount+kp
endif
phi0=phi1
phi1=phi2
10 continue
100 format('r,phi,D=',3(3x,e10.3))
stop
end
 

 

 

 

*****

c Emden's equation HT Davis NL DE page 371
c d^2phi/dr^2+(2/r)d phi/dr= -a^2*phi^n Phi~gravitational potential
c Ic phi(0)=1 . phi(dr)=1.
data a,nstep ,n ,rf/1.,1000,3,4./
phi0=1.
phi1=1.
dr=rf/float(nstep)
kp=int(float(nstep)/80.)
kount=kp
print *,'y(0)=1.,(dy/dx)i=0.n=', n
c print 100,0.,phi0,1.
do 10 i=2,nstep
r=dr*float(i)
phi2=2.*phi1-phi0 +dr**2*(-(2./(r-dr))*(phi1-phi0)/dr
$-a**2*phi1**n )
if(i.gt.2.and.i.le.10)then
D=(1./r)*(phi2-phi1)/dr
print 100 ,r , phi2 ,D
kount=kount+kp
endif
phi0=phi1
phi1=phi2
10 continue
100 format('r,phi,D=',3(3x,e10.3))
stop
end