Ph361-Spring 2011: Stars

Introduction

This is an introductory course about stars, their structure, evolution and collective motion in galaxies and globular clusters. I will try to give you some exposure to the fundamentals, some familiarity with astronomical jargon and conventions and practice with doing calculations so that you will be in a good position to appreciate and follow modern research. However, I will not be able to give much of a discussion of that research. 

The two texts that are most useful for the two areas that are covered are:
Kippenhahn & Weigert - Stellar Structure and Evolution
Binney & Tremaine - Galactic Dynamics
These go far beyond what we will cover but you will not regret owning them!

Other texts that may be useful include:
Hansen & Kawaler - Stellar Interiors
Binney - Galactic Astronomy
Gray - Observation and Analysis of Stellar Photospheres
Bahcall - Neutrino Astrophysics
Chaisson & McMillan- Astronomy Today (A non-science major text but up to date on qualitative astronomy and terminology.)

1. Stellar Structure

Basics

Static Equilibrium

The equation of static equilibrium is ∇P=ρg=-ρ∇Φ. Note that isobars, equipotentials, and isobars coincide. The law of gravitation is -∇·g=(∇^2)Φ=4πGρ. These can be combined for spherical stars as (1/(r^2))(d/dr)((r^2)/ρ)(dP/dr)+4πGρ=0. To solve these, we need to know P(ρ) explicitly. Most stars determine P(ρ) implicitly using other physics and equations so this is not very useful in practice

Simple Stellar Models (KW2)

For a ``barytropic’’ star we can write y=∫(dP/ρ)- equal to the enthalpy when isentropic so that  (((d^2) y)/(dr^2))+(2/r)(dy/dr)+4πGρ(y)=0. 
For a uniform density star, P=(Subscript[P, 0])-((2((πGρ)^2)(r^2))/3), the radius R=(((((3(Subscript[P, 0]))/(2((πGρ)^2)))))^(1/2)), and the mass M=((((6/π)))^(1/2))(((((Subscript[P, 0])/(G (ρ^(4/3))))))^(3/2)). 
So,  (Subscript[P, 0]) = ((((π/6)))^(1/3)) (GM^(2/3))(ρ^(4/3))=((3(GM^2))/(8((πR)^4))).

Off[Solve::ifun];iP=ρ y[r]/.DSolve[{y''[r]+(2/r)y'[r]+4π iG ρ==0,y[0]==(iP0/ρ),y'[0]==0},y[r],r][[1]];
iR=Simplify[r/.Solve[iP==0,r][[2]],Assumptions->{iP0>0,iG>0}];
iM=((4π ρ (iR^3))/3);
m=iM((((r/R)))^3);Plot[{(x^3),1-(x^2)},{x,0,1},Frame->True,FrameLabel->{"(r/R)","(P/(Subscript[P, 0])),(m/M)","Uniform Density",""},RotateLabel->False]

Q1.1: How well does this model reproduce the central pressure of the Earth?

Helen Craig 
(Subscript[P, o])=(3/8) ((GM^2)/((πR)^4))and G=6.67e-11 (m^3)/(s^2) kg, M=5.98e24 kg, R=6.37e6m so
(Subscript[P, o]) = 172 GPa
From Handbook of Chemistry and physics, Pressure=361.7 GPa.
So for an order of magnitude, the model seems to give the right answer but it is also about half what it should be

Q1.2: Solve for the structure of a gaseous planet where P∝(ρ^2)approximately. What do you notice about the radius? Check with Jupiter, Saturn.

Rachel Reddick
Gas giant with P ∝ (ρ^2)
	Need to solve (1/(r^2))(d/dr)((r^2)/ρ)(dP/dr)+4πGρ=0
	In this setup, dP/dr = 2 a ρ dρ/dr
	So, it follows that:

G := 6.67*^-11 (*m^3/kg/s^2*)
Clear[G]
2 a/r^2 D[r^2D[ρ,r],r] + 4 Pi G == 0

So, set up to solve the differential equation here:
dens[r_]:=FullSimplify[ρ[r]/.DSolve[{2 a ρ''[r] + 4 a/r ρ'[r] + 4 Pi G ρ[r] == 0, ρ[0] == ρ0, ρ'[0]==0},ρ[r],r][[1]]]
dens[r]
dens[r_] := (((Sqrt[a]) ρ0 Sin[(((Sqrt[G]) (Sqrt[2 π]) r)/(Sqrt[a]))])/((Sqrt[G]) (Sqrt[2 π]) r))
FullSimplify[dens[Sqrt[Pi a /G/2]]]
myR = Sqrt[Pi a /G/2];
myM = Integrate[dens[r] 4 Pi r^2,{r,0,myR}] 
m[r_] := Integrate[dens[rp] 4 Pi rp^2 ,{rp,0,r}]/myM
m[r*myR]
Plot[((-π x Cos[π x]+Sin[π x])/π),{x,0,1},Frame->True,FrameLabel->{"r/R","m/M"}]

FullSimplify[dens[x*myR]/ρ0]
((Sin[π x])/(π x))
Plot[(Sin[Pi x]/(Pi x))^2,{x,0,1},Frame->True,FrameLabel->{"r/R","P(r)/P0"}]

FullSimplify[dens[x*myR]/ρ0]
((Sin[π x])/(π x))
Plot[(Sin[Pi x]/(Pi x))^2,{x,0,1},Frame->True,FrameLabel->{"r/R","P(r)/P0"}]

The strange this about this setup is that the radius of the planet is actually independent of its mass (!).

Jupiter’s mass is 1.9e27 kg; radius 71,000 km
Saturn’s mass is 5.68e26 kg; radius 60,000 km  (Note:  These are the equatorial radii.)

So -- despite a factor of four difference in mass, the two planets actually do have very similar radii.  Some of the difference is probably atributable to the fact that the two planets rotate at different rates.  Or there could be a small difference in the proportionality constant in the P∝ρ^2 relation.

Rotation

Rotating Stars (KW41,42, Tassoul 1978)

For non-rotating uniform density star, Φ=(Subscript[Φ, 0])+((2((πGρr)^2))/3). In order to continue to satisfy Poisson's equation for a slowly rotating starwe must add a term (Ar^2)(Ω^2)(Subscript[P, 2])(μ) . This is the leading correction. To match at the surface, the external potential must have the form Φ = -(GM/r) + ((((BΩ)^2) (R^5))/(r^3))(Subscript[P, 2])(μ). We must solve the problem in the rotating frame where there is no motion by introducing the centrifugal potential (Subscript[Φ, cen])=-(1/2) | Ω×r ((|)^2)=(1/3)(Ω^2)(r^2)((Subscript[P, 2])(μ)-1). The angular variation of the potential satisfies ((4πGρR)/3)δr+(A+(1/3))(Ω^2)(R^2)(Subscript[P, 2])=((4πGρR)/3)δr+(B+(1/3))(Ω^2)(R^2)(Subscript[P, 2])=0 implying A=B. The radial gradient must also match.
((4 πGρr)/3)+2((AΩ)^2)(Subscript[RP, 2])= ((4 ((πGρR)^3))/(3(r^2)))-3((BΩ)^2)(Subscript[RP, 2]), dropping the centrifugal rms which automatically cancel. Hence A=B=(1/2)  and δr = -((5 (Ω^2) R)/(8 πGρ))(Subscript[P, 2 ]). (Note this keeps the volume unchanged as ((∫_(Subscript[, -1]))^1)(Subscript[dμP, 2])(μ)=0.) So, (Subscript[R, eq])-(Subscript[R, pole])=((5 (Ω^2) (R^2))/(4g)).
Fast rotators are MacLaurin spheroids. 
To solve the more general problem, we use ellipsoidal coordinates.

Q1.3: What would be the answer if all the mass were concentrated at the center? Which is the better approximation for the Earth?

Q1.4: Show that ∇P=-ρ(∇Φ+Ω×(Ω×r)) for a differentially rotating star. Assume that P=P(ρ) and expand in cylindrical (s,z) coordinates. Show that Ω=Ω(s). Does the sun satisfy this rule?

Binary Stars (See KW2,41)

Restricted 3 - body problen

Much has learned from analyzing stars in binary systems. For the simple case, the orbit is circular, the stars co-rotate and their mass is assumed to be centrally concentrated. The potential in the rotating frame is Φ = -G[(Subscript[m, 1]) (({(((s + (Subscript[a, 1])))^2)+(z^2)+(c^2)})^(-1/2)) +(Subscript[m, 2]) (({(((s-(Subscript[a, 2])))^2)+(z^2)+(c^2)})^(-1/2))+(1/2)(ma^(-3))(s^2)], where c is a softening parameter. Extrema are “Lagrange” points, L1-L5.
Some of these potential considerations will re-surface when we discuss stellar dynamics.

Φ=-(m1/((((((s+((m2 a)/(m1+m2))))^2)+(z^2)+(c^2)))^(1/2)))-(m2/((((((s-((m1 a)/(m1+m2))))^2)+(z^2)+(c^2)))^(1/2)))-(((m1+m2)(s^2))/(2(a^3)))/.{m1->1,m2->.5,a->1,c->.1};
Φplt=ContourPlot[Φ,{s,-2,2},{z,-1,1},Contours->{-2.48,-2.65,-2.9},AspectRatio->Automatic,FrameLabel->{"s","z","Vertical equipotentials",""}]
Φ=-(m1/((((((s+((m2 a)/(m1+m2))))^2)+(y^2)+(c^2)))^(1/2)))-(m2/((((((s-((m1 a)/(m1+m2))))^2)+(y^2)+(c^2)))^(1/2)))-(((m1+m2)((s^2)+(y^2)))/(2(a^3)))/.{m1->1,m2->.5,a->1,c->.1};
Φplt=ContourPlot[Φ,{s,-2,2},{y,-2,2},Contours->{-2.09,-2.48,-2.65,-2.9},AspectRatio->Automatic,FrameLabel->{"s","y","Orbital equipotentials",""}]

Roche Limit

Mass transfer

As a result of orbital or stellar evolution, one star, mass m1, may expands to fill the equipotential surface up to the L1 saddle point. The volume enclosed by this critical equipotential surface, (Subscript[V, 1]), has an effective radius (Subscript[R, 1])=(((3V/4π))^(1/3))  where it can be shown that
(Subscript[R, 1])/a=0.38+0.2log((Subscript[m, 1])/(Subscript[m, 2])) , 0.3<(Subscript[m, 1])/(Subscript[m, 2])<20;  =0.462((((Subscript[m, 1])/m))^(1/3)) , (Subscript[m, 1])/(Subscript[m, 2])<0.8 .
The second formula, which is in practice the more useful, is equivalent to P≈10(((((Subscript[ρ, 1])/(1g (cm^(-3))))))^(-1/2))hr. This is the critical period for mass transfer to start.

plt1=Plot[.38+.2Log[10,q],{q,.3,20},PlotRange->{0,.7},Frame->True];
plt2=Plot[.46224((((q/(1+q))))^(1/3)),{q,0,.8},Frame->True];
Show[plt1,plt2,Frame->True,FrameLabel->{"(Subscript[m, 1])/(Subscript[m, 2])", "(Subscript[R, 1])/a","Roche Lobe radius",""},RotateLabel->False]

Q1.5 Verify the formula for the critical period from the quoted relation for (Subscript[R, 1])/a.

(Subscript[R, 1])/a=.462((((Subscript[m, 1])/m))^(1/3)) and (P^2)=4(π^2)(a^3)/((Subscript[Gm, 2]))

So (P^2)=4(π^2)/((Subscript[Gm, 2])) * (((Subscript[R, 1])^3)m)/(.462^3(Subscript[m, 1]))
= ((4(π^2)m3V)/((Subscript[Gm, 2])4(Subscript[πm, 1])((.462)^3)))
G=.865(cm^2)(g^(-1))(hr^(-2))
P^2 = 110.492m/((Subscript[m, 2])ρ)   and (Subscript[m, 1])/(Subscript[m, 2])<.8 so 1 < (((m/(Subscript[m, 2])))^(1/2))<1.34
and P≈10(ρ^(-1/2))

Q1.6* Verify the quoted relations for (Subscript[R, 1])/a  by integrating the Roche volumes numerically

Helen Craig
clear all;
close all;
numB = 500;
for indexB = 1 : 1 : numB
b = indexB/numB*20; % b is the ratio between masses
m (1) = 1/(1 + 1/b); % find m1 and m2
m (2) = 1/(1 + b);
[s, fval] = fminbnd (@(s) potential (s, m), -m (2), m (1) ); % find the point of maximum
barray (indexB) = b;				% potential between the masses
sarray (indexB) = 1 - (m (1) - s);
end
plot (barray, sarray)
hold on;
x = 0 : .01 : .8
y = .462*(1./(1 + 1./x));
plot (x, y, ' r') % plot the analytical solution
x = .3 : 1 : 20
y = .38 + .2*log (x);
plot (x, y, ' g')

Q1.7 We have assumed that (Subscript[m, 1]) is corotating. How do you think the criterion for mass transfer will change if it is not?

Question 1.7

Binary star system -- what happens to the mass transfer criterion if m1 is not corotating?

If m1 were the more massive object, dominating the system, then the center of rotation of the object and the binary is essentially the same.  In that case, if m1 were rotating more rapidly that its secondary orbited, then the radius at which material would cease to be bound would be smaller, so the radius of mass transfer should decrease in that case.  The opposite is true if it rotates more slowly, or somewhat retrograde -- the radius for mass transfer would increase.  Very retrugrade, and it would be smaller again.

In short, if m1 has a larger magnitude of angular momentum than the corotating case, the radius for mass transfer will decrease; smaller, and it will increase, to some maximum at zero angular momentum.  This should also apply to the more general case of possibly comparable masses, or when m1 is less massive than m2.

(I should take a look at the math for this, but I haven't had a chance yet.)

Compact Objects

White Dwarfs (See KW15,35, Shapiro & Teukolsky 1983)

Structure

Equation of State

A problem that is not artificial is the structure of a cold, fully degenerate white dwarf. First we need the density adopting the Pauli principle and assuming a degenerate gas of free electrons 
ρ=((8π (Subscript[μ, e])(Overscript[(Subscript[m, p]), ]) )/(3((Subscript[λ, c])^3)))(u^3)=1.9× (10^6)(u^3)g (cm^(-3))  where uc is the Fermi proper speed, (Subscript[λ, c])=2.4 (10^(-10))is the Compton wavelength and (Subscript[μ, e])=2 is the mean molecular weight per electron. 
The pressure is given by P=((8 (Subscript[πm, e]) (c^2) )/(3((Subscript[λ, c])^3)))((∫_0)^u)du' ((u')^4)(((1+((u')^2)))^(-1/2)) = 6.0×(10^22)(((u(1+(u^2)))^(1/2))(2(u^2)-3)+3(sinh^(-1))u) dyne (cm^(-2)). 
The specific heat ratio, Γ=(dlnP/(dlnρ))falls from (5/3)to (4/3)as it should as the Fermi speed becomes ultrarelativistic.

ρ=((8π 2 1.66 (10^(-24)))/(3 (((2.43 (10^(-10))))^3)))(u^3);
iP=((8π  9.11 (10^(-28))(((3 (10^10)))^2))/(3(((2.43 (10^(-10))))^3)))Integrate[(u^4)(((1+(u^2)))^(-1/2)),u];
ParametricPlot[{Log[10,ρ],Log[10,iP]},{u,.01,10},Frame->True,AspectRatio->Full,FrameLabel->{"log(ρ)","log(P)","Equation of State",""},RotateLabel->False];
Γ=((D[Log[iP],u])/(D[Log[ρ],u]));
Plot[Γ,{u,.01,10},Frame->True,FrameLabel->{"u","Γ","Specific Heat Ratio",""}]

Static Equilibrium

Note that dy=(((Subscript[m, e])(c^2))/((Subscript[μ, e])(Subscript[m, p])))dγ, where γ=(((1+(u^2)))^(1/2)) is the Fermi Lorentz factor.
Our differential equation then becomes (((d^2)γ)/((dξ)^2))+(2/ξ)((dγ)/(dξ))+((16π)/3)((((γ^2)-1))^(3/2))=0, where ξ =(r/(Subscript[r, c]))= ((((Subscript[μ, e])(Subscript[m, p]))/(Subscript[m, Pl]))) ((r/(Subscript[λ, c]))) = (r/(1.6 (10^9)cm)), (Subscript[m, Pl])=(((((ℏc)/G)))^(1/2))=2.2×(10^(-5))g is the Planck mass. 
The outer radius is R=(Subscript[ξ, m])(Subscript[r, c])where γ((Subscript[ξ, m]))=1. The mass is M=4((Subscript[πr, c])^3)(Overscript[(∫_0), (Subscript[ξ, m])])((dξξ)^2)ρ. The mass increases very slowly with increasing central density (Fermi energy) The limiting mass is ~1.45(Subscript[M, ⊙]). The basic reason why there is a limiting mass is that when the specific heat ratio becomes Γ~4/3, then the internal energy U scales as (R^(-1)) just like the gravitational energy and an equilibrium is not possible as it is when Γ~5/3 and U∝(R^(-2)) .  The limiting mass is the Chandrasekhar mass. More massive stars must collapse to make a neutron star or black hole. Note that for non-relativistic stars, P~(M^2)(R^(-4 ))~ ((((MR^(-3))))^(5/3))=>R~(M^(-1/3)).

rc=((2.43 (10^(-10))2.2 (10^(-5)))/(2 1.66 (10^(-24))));
ξ0=.001;γ0=2;ξm=.6;
wd=NDSolve[{ γ''[ξ]+(2/ξ)γ'[ξ]+((16π)/3)(((Max[0,((γ[ξ])^2)-1]))^(3/2))==0,γ[ξ0]==γ0,γ'[ξ0]==0},γ,{ξ,ξ0,ξm}][[1,1]];
ξm=.999ξ/.FindRoot[Evaluate[γ[ξ]/.wd]==1,{ξ,.2}];
iR=rc ξm
iM=((4π (rc^3))/(1.99 (10^33)))NIntegrate[(ξ^2)ρ/.u->(((((γ[ξ])^2)-1))^(1/2))/.wd,{ξ,ξ0,ξm}]
Plot[γ[ξ]/.wd,{ξ,ξ0,ξm}, Frame->True]

Q1.8 Compute the total (internal +gravitational) energy of a cold white dwarf.

Q1.9 Estimate the thermal energy for a hot, 0.6(Subscript[M, ⊙]) white dwarf with isothermal interior temperature T~(10^7)K. (Hint this will mostly be associated with the ions.)

Limitations

This description is quite accurate but allowance must be made for Coulomb effects, degenerate envelopes in young, hot stars, crystalline cores in old, cold stars, pycnonuclear reactions and electron capture in high mass stars.

Observations

White dwarfs are very common and observed with masses from 0.2 to 1.3(Subscript[M, ⊙])(median 0.6(Subscript[M, ⊙])) and surface temperatures from 4000-150000K depending upon age. They are very common as they are the fate of most low mass stars that have evolved off the main sequence. The most common composition is C-O  though He and heavier element stars can form.  White dwarfs are observed as Cataclysmic Variables in binary systems and can have periods as short as 5 min in ultracompact binaries. Accreting white dwarfs are believed to become Type 1a supernovae.

Q1.10* Sirius A - a 2 (Subscript[M, ⊙])A star with apparent magnitude -1.5 and absolute magnitude 1.4, the brightest star in the sky - has a parallax of 0.4” and a distance of 2.6 pc. It is a visual binary with a 50 yr period and semi major axis of 7.5”. Its1(Subscript[M, ⊙])  white dwarf companion has a radius of 0.0008(Subscript[R, ⊙]), a surface temperature of 25000K and a luminosity of 0.026(Subscript[L, ⊙]). Use your knowledge of astronomy and the mass-radius relation to demonstrate that this description is self-consistent.

Black Holes(KW37, Shapiro & Teukolsky 1983)

Schwarzschild Metric

A non-rotating hole is described by the Schwarzschild metric (ds^2)=-(1-((2m)/r))(dt^2)+((dr^2)/(1-2m/r))+(r^2)(((dθ)^2)+(sin^2)θ ((dϕ)^2)) in geometrical units where r is a circumferential radius and m is the gravitational mass which defines the scale. There is an event horizon at r=2m= 3(M/(Subscript[M, ⊙]))km.

Kerr Metric

A spinning hole is described by the Kerr metric which includes an additional parameter, the angular momentum in unit mass squared in geometrical units, a<m.  The event horizon is at r=m+((((m^2)-(a^2)))^(1/2)). The area theorem can be used to show that there is an irreducible mass (Subscript[m, i])= (([m(m+((((m^2) - (a^2)))^(1/2)))/2])^(1/2)). Note that m>(Subscript[m, i])>(2^(-1/2))m so that up to 0.29 of the rest mass is extractable in principle.

Observed Stellar Holes

Stellar black holes are most easily observed in binary systems We do not know the number of stellar holes in the Galaxy.

Neutron Stars(KW36, Shapiro & Teukolsky 1983)

TOV Equation

Consider a non-rotating, spherically symmetric neutron star. The exterior metric must be reducible to the Schwarzschild metric through Birkhoff's theorem. Let us continue to define r in the same way as for black holes. The metric can then be written in a standard form: (ds^2)=-(e^(2Φ))(dt^2)+((dr^2)/(1-2m/r))+(r^2)(((dθ)^2)+(sin^2)θ ((dϕ)^2)) where Φ=Φ(r) and m=m(r). 
For an observer at rest (u^α) =((e^(-Φ)),0,0,0), (Subscript[u, α]) =(-(e^Φ),0,0,0). Use the vanishing divergence of the stress energy tensor to obtain: (Subscript[P, ,r])=(wu^t)(Subscript[u, t,r]) = -(Subscript[wΦ, ,r]) or (dP/(dΦ))=-w=-(P+ρ)  
Next, evaluate the mixed Einstein tensor.

Clear[coord,metric,inversemetric,affine,ρ,m];
coord={t,r,θ,ϕ};
metric={{-Exp[2Φ[r]],0,0,0},{0,(((1-2m[r]/r))^(-1)), 0,0},{0,0,(r^2),0},{0,0,0,(r^2)((Sin[θ])^2)}};
inversemetric=FullSimplify[Inverse[metric]];
affine=Simplify[Table[(1/2)Sum[(inversemetric[[α,δ]])
(D[metric[[δ,β]],coord[[γ]] ]+
D[metric[[δ,γ]],coord[[β]] ]-D[metric[[β,γ]],coord[[δ]] ]),{δ,1,4}],
{α,1,4},{β,1,4},{γ,1,4}] ];
riemann=
Simplify[Table[D[affine[[α,β,δ]],coord[[γ]]]-D[affine[[α,β,γ]],coord[[δ]]]+Sum[affine[[α,ϵ,γ]]affine[[ϵ,β,δ]]-affine[[α,ϵ,δ]]affine[[ϵ,β,γ]],{ϵ,1,4}],{α,1,4},{β,1,4},{γ,1,4},{δ,1,4}]];
ricci=FullSimplify[Table[Sum[inversemetric[[α,δ]]riemann[[γ,δ,γ,β]],{γ,1,4},{δ,1,4}],{α,1,4},{β,1,4}]];
einstein=FullSimplify[ricci-(1/2)Tr[ricci]{{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}}]

so that 
(Subscript[(G^t), t]) = -2((Subscript[m, ,r])/(r^2)) =8π (Subscript[(T^t), t])=-8πρ => m = 4 π ((∫_0)^r)dr'((r')^2)ρ(r')
and 
(Subscript[(G^r), r]) =2(1-((2m)/r))(Subscript[Φ, 'r])-((2 m)/(r^3)) =8(Subscript[((πT)^r), r])=8πP  
so that
(Subscript[P, ,r])=-(Subscript[wΦ, ,r])=-(((m+4((πr)^3)P)(ρ+P))/((r^2)(1-2m/r))).
This is the Tolman-Oppenheimer-Volkoff equation. m is the interior mass energy and equals the m determined by eg orbits at large r at the surface of the star. m+4((πr)^3)P is the active gravitational mass and w the passive. The (1-2m/r) is kinematic. If the rest mass per baryon is μ and the proper density n, then 
(Subscript[m, rest])=4π∫dr (((r^2)(1-2m/r))^(-1/2))μ n. We can also define quasi-Newtonian internal energy density ρ-μn and gravitational  energy density -ρ[1-(((1-2m/r))^(1/2))] which agree with standard non-relativistic relations in the limit. At the surface, where P=0, r=R, m=M, g=(1/ρ)(dP/dr)=(M/(((R^2)(1=2M/R))^(1/2))).

Limiting Mass

Neutron stars are roughly10 km in radius - (10^(57/3))fm. They have a maximum mass for much the same reason as white dwarfs, although the physics is much less well constrained and speculative possibilities like strange quark stars cannot be ruled out. The effect of general relativity is to increase the critical exponent from Γ=4/3 to a value much closer to 2 as the star becomes denser. It is called the Oppenheimer-Volkoff mass ~2(Subscript[M, ⊙]). In 2010, Demorest et al measured a Shapiro delay in the millisecond pulsar PSR J1614-2230 and found its mass to be 2 (Subscript[M, ⊙]). This rules out soft candidate equations of state. 
Neutron stars have solid crusts, liquid oceans, superconducting and super fluid interiors. Magnetized neutron stars are radio and X-ray pulsars. They are formed in core collapse supernovae and there are thought to be about 3×(10^8)  in the Galaxy.

Q1.12 Adopt a simple EOS for high density neutron matter, (Subscript[P, 35])=0.02((Subscript[ρ, 35])^2), where the units are (10^35)erg (cm^(-3)). Obtain a rough mass-radius relation using the TOV equation. 
Question 1.12

Take (Subscript[P, 35])=0.02((Subscript[ρ, 35])^2), and get the mass-radius relation

In general (TOV):
D[mP[r],r] == - (m[r] + 4 Pi r^3 mP[r]) (ρ[r] + mP[r])/(r^2 (1-2m[r]/r))
Plugging in the equation of state as mP = k ρ^2, and noting that m[r] = 4 Pi Integral[r^2 ρ[r’],{r’,0,r}]
soln = DSolve[ {ρ'[r] == -1/(2k)(m[r]+4Pi r^3 k ρ[r]^2)(1+k ρ[r])/(r (r-2m[r])),
m'[r]==4 Pi r^2 ρ[r],
m[0]==0,ρ[0]==ρ0},{ρ[r],m[r]},r]
DSolve[{ρ'[r]==-(((1+k ρ[r]) (m[r]+4 k π (r^3) ((ρ[r])^2)))/(2 k r (r-2 m[r]))),m'[r]==4 π (r^2) ρ[r],m[0]==0,ρ[0]==ρ0},{ρ[r],m[r]},r]
Let’s try this again, working numerically, and include all the standard cgs units:
G = 6.67*^-8;
 
(*cm^3/kg/s^2*)
c = 3*^10; (*cm/s*)

msolar = 2*^33; (*g -- using to sorto out units*)
rscale = 100000; (*cm -- 1 km scale*)

nk = 2*^-37c^4;  (*erg/(g/cm^3)^2*)
nρ0 = 10^17;  (**central density -- free parameter -- g/cm^3*)

nrmin=.01;  (*km; numerical cutoff to get an actual answer*)
nrmax = 30; (*km; reasonable upper bound to cover max radius*)

mp0 = N[4 Pi (nrmin rscale)^2 nρ0]/msolar
m0 = N[4/3 Pi (nrmin rscale)^3 nρ0]/msolar

nsoln = NDSolve[{(r rscale)^2 m''[r ] msolar/rscale^2-2r rscale  m'[r ] msolar/rscale == -G/(2nk) (r rscale m[r ] msolar+nk/(4Pi)m'[r ]^2 msolar^2/rscale^2/c^2)(4Pi (r rscale)^2 + nk m'[r ] msolar/rscale/c^2)/(r rscale-2G m[r]msolar/c^2),
m'[nrmin]==4 Pi (nrmin rscale)^3 nρ0/msolar,m[nrmin]==4/3 Pi (nrmin rscale)^3 nρ0/msolar},m[r],{r,nrmin,nrmax}][[1]]
6.283185307179586`*^-10
2.0943951023931955`*^-7
{m[r]->InterpolatingFunction[{{0.01`,30.`}},"<>"][r]}
Now, plot the mass enclosed at each given radius, and density at each radius:
Plot[m[r]/.nsoln,{r,nrmin,nrmax},AxesLabel->{"r (km)","m(r) Msolar"}]
Plot[D[m[r]/.nsoln,r]/(4 Pi (rscale r)^2)*msolar/.r->myr,{myr,nrmin,nrmax},AxesLabel->{"r (km)","ρ(r) g/cm^3"}]



Evaluate[D[m[r]/.nsoln,r]]/.r->2*nrmin
myr/.FindRoot[ D[m[r]/.nsoln,r]/.r->myr,{myr,10*nrmin,nrmax}][[1]]
Evaluate[m[r]/.nsoln]/.r->%
2.540020299021696`*^-6
14.457135276876617`
1.8682812342174728`
Since this is all parametrized (as set up) by setting the central density, loop over the central density to get enough points to be worth plotting:
ρlist = {.5,.6,.7,.8,.9,1,1.25,1.5,2,2.5}*10^17;
mlist = {};
rlist = {};
pairlist = {};

For[i=0,i<Length[ρlist],i++,
nρ0= ρlist[[i+1]];
nsoln = NDSolve[{(r rscale)^2 m''[r ] msolar/rscale^2-2r rscale  m'[r ] msolar/rscale == -G/(2nk) (r rscale m[r ] msolar+nk/(4Pi)m'[r ]^2 msolar^2/rscale^2/c^2)(4Pi (r rscale)^2 + nk m'[r ] msolar/rscale/c^2)/(r rscale-2G m[r]msolar/c^2),
m'[nrmin]==4 Pi (nrmin rscale)^3 nρ0/msolar,m[nrmin]==4/3 Pi (nrmin rscale)^3 nρ0/msolar},m[r],{r,nrmin,nrmax}][[1]];
rlist = Append[rlist,myr/.FindRoot[ D[m[r]/.nsoln,r]/.r->myr,{myr,10*nrmin,nrmax}][[1]]];
mlist = Append[mlist,Evaluate[m[r]/.nsoln]/.r->rlist[[i+1]]];
pairlist=Append[pairlist,{mlist[[i+1]],rlist[[i+1]]}]]
rlist
mlist
{16.392365042749464`,15.9350417990007`,15.517413602751011`,15.134546105301675`,14.457135276876617`,13.146521483250018`,12.201567646676127`,11.489337314929216`}
{1.397921816345817`,1.5341458571215547`,1.64470218359452`,1.7347340776946303`,1.8682812342174728`,2.0377107426219605`,2.0913787225272866`,2.096050557271169`}
ListLinePlot[pairlist,AxesLabel->{"M (msolar)","R (km)"}]

Oscillations

Basic Ideas (KW38,39,40)

Radial Modes

The philosophy is quite simple. We write down the equations of motion and perform linear perturbations. We assume that the perturbations are adiabatic as in a sound wave and arrive at a second order ordinary differential equation for the radial displacement. This can be cast in Sturm-Liouville form. There is an infinite sequence of eigenmodes with increasing eigenfrequencies which depend upon the equilibrium structure of the star. It is also possible that the perturbations are non-adiabatic as they allow changes in opacity and energy generation. These can drive unstable modes as Cepheid variable stars.

Non - radial Modes

It is also possible to consider non-radial modes. Here there are two coupled partial differential equations if we ignore the change in the gravity - Cowling approximation - and four if we do not. The equations are separable and the angular parts involve spherical harmonics. There are two types of modes, p- modes driven by pressure and g-modes driven by buoyancy. There is also a single f mode for each l>1 which combines both features. Rotating stars break the degeneracy between different m values. This theory finds spectacular application in the sun.

White Dwarfs

The simplest application is to white dwarfs.

Q1.11* Solve for the fundamental p mode frequency of a 0.6(Subscript[M, ⊙])white dwarf.

2. The Sun

Fundamentals

Gas

Perfect Gas

Using standard kinetic theory/stat mech, <pv> =3kT=> P=(1/3)n<pv> =nkT=((ρkT)/(Subscript[μm, p])) where μ≃(((2X+3Y/4+Z/2))^(-1))=(((3X/2+Y/4+1/2))^(-1))≃0.6 for X=.73,Y=.25, Z=.02. So P=1.4×(10^(14 ))ρ (Subscript[T, 6])dyne (cm^(-2)) . nb middle-aged sun has more helium and larger μ.

μ=(((1.5 .73+.25 .25+.5))^(-1));
pg=((ρ 1.38 (10^(-16))(10^6)t6)/(μ 1.67 (10^(-24))));

Q2 .1* Show from first principles that the entropy per unit mass of a perfect gas of atomic hydrogen is (S/k)=(1/(Subscript[m, p]))ln(((e^(5/2))/((nλ)^3)))  where λ=(h/(((2(Subscript[πm, p])kT))^(1/2))). Hence show that the entropy of a fully ionized gas of fixed composition is (S/k)=(1/(Subscript[μm, p]))ln(((T^(3/2))/ρ))-(1/(Subscript[m, p]))[XlnX+(Y/4)ln(Y/32)+((1+X)/2)ln(((1+X)(1836^(3/2)))/2)]+const. (Note that the full expression must be used if the ionization changes.)

λ=6.63 (10^(-27))(((2π 1.38 (10^(-16))1.67 (10^(-24)) (10^6)t6))^(-1/2));
(1/(1.67 (10^(-24))))Log[ ((1.67 (10^(-24)))/( ρ (λ^3)))(E^(5/2))];

Helen

We can use dS=((-dF)/dT) and F=-kTlog(Z)
To find Z, use Z=(∏_j^(3N))(1/h)(∫_(-∞)^∞)(Subscript[dq, j])exp(-(Subscript[βp, j])/2(Subscript[m, j]))/N=(((L/h(Sqrt[2((πm)/β)])))^(3N))/N=(((L/λ))^(3N))/N!
calculate F:

F=-kTNlog(V/(λ^3))+kTlog(N!) and log(N!)≈Nlog(N)-N

So F = -NkT (log(V/Nλ^3)+1)=NkT(log(((ρλ)^3))-1)

Calculate S:

S=((-dF)/dT) and (d/dT)log(((ρλ)^3))=(d/dT)(logT^(-3/2))=-3/2T

S = -Nk (log((((ρλ)^3)))-1)+(3/(2T))NkT=k log(((e^(5/2))/((ρλ)^3)))

If μ=(((3X/2+Y/4+1/2))^(-1))

S=(1/(μm))log(((e^(5/2))/ρ)((((2πmμkT))^(3/2))/(h^3)))=(1/(μm))log(((T^(3/2))/ρ))-(1/(μm))log((μ^(3/2)))=(1/(μm))log(((T^(3/2))/ρ))-(3/2)(3X/2+Y/4+1/2)log(3X/2+Y/4+1/2)

Ionization Equilibrium (KW14)

Collisions maintain Local Thermodynamic Equilibrium. Let there be energy levels (Subscript[E, i]) (ignoring degeneracy which we can always break by perturbing the Hamiltonian). The occupancy of state i satisfies the Boltzmann distribution  (Subscript[n, i])∝(e^(-(Subscript[E, i])/kT)). Now consider an unexcited ion in ionization state r and compare its occupancy to that of an unexcited ion in state r+1, separated by ionization potential (Subscript[χ, r])(measured in eV) accompanied by a free electron with energy (p^2)/2(Subscript[m, e]) in a volume (d^3)p of  momentum space. The Boltzmann law gives ((Subscript[dn, r+1,0])/(Subscript[n, r,0]))= ((2(d^3)p)/((Subscript[n, e])(h^3)))(e^(-((Subscript[χ, r])+(p^2)/2m)/kT)). Integrating over momentum, we get the Saha equation
 (((Subscript[n, r+1,0])(Subscript[n, e]))/((Subscript[n, r]),0))= 2(((((2(Subscript[πm, e])kT)/(h^2))))^(3/2)) (e^(-(Subscript[χ, r])/kT)). 
If we need to include ions in excited states, then we simply apply the Boltzmann law again to obtain (Subscript[n, r])=(Subscript[n, r,0])(Subscript[Z, r]), where (Subscript[Z, r])=(∑_i)(e^(((Subscript[E, r,0])-(Subscript[E, r,i]))/kT)) is the partition function. 
So ((Subscript[n, r+1])/(Subscript[n, r]))=((2/(Subscript[n, e]))) (((Subscript[Z, r+1])/(Subscript[Z, r]))) (((((2 (Subscript[πm, e])kT)/(h^2))))^(3/2))(e^(-(Subscript[χ, r/kT])))=4.8((Subscript[n, e24])^(-1))(((Subscript[Z, r+1])/(Subscript[Z, r])))((Subscript[T, 6])^1.5)(e^(-((.0116χ)/(Subscript[T, 6])))). 
Note that when the gas is dilute, there are many free electron states available and the ionization occurs when χ~3-20 kT.

(2/((10^24)n24))(((((2π 9.1 (10^(-28))1.38 (10^(-16))(10^6)T6)/(((6.63 (10^(-27))))^2))))^1.5)(E^(-χ 1.6 (10^(-12))/(1.38 (10^(-16))(10^6)T6)));

Apply to solar photosphere where (Subscript[n, 24])∼(10^(-7)), (Subscript[T, 6])∼ 6×(10^(-3))  (Subscript[n, p])/(Subscript[n, H])∼

Q2 .2 Estimate the temperature at which hydrogen ionizes in a stellar atmosphere and explain in words why it is much less than χ/k.

Black Body Radiation

Planck Law

dE= (Subscript[I, νΩ])dνdΩdAdt =2(((ν^2)dνdΩ)/(c^3))nhνdAdt where (Subscript[I, νΩ]) is the intensity and n is the occupation number =((((e^(hν/kt))-1))^(-1)) in Bose equilibrium. (Subscript[I, νΩ])=(Subscript[B, νΩ])=((2((hν)^3))/(c^2))((((e^(hν/kt))-1))^(-1))the Planck distribution
Hence U = 8 π ((((kT/hc)))^3) kT ∫(((dxx^3)((e^x) - 1))^(-1)) =((8(π^5)(k^4))/(15(h^3)(c^3)))(T^4)
P = < (cos^2)θ>U=(1/3)U= ((4 σ)/(3c))(T^4) = (a/3)(T^4)=2.52×(10^9)((Subscript[T, 6])^4) dyne (cm^(-2))

rP=2.52 (10^9)(T6^4);σSB=5.7 (10^(-5));

Q2 .3 Show from first principles that the entropy of black body radiation is 4U/3T =74(T^3)k per unit volume

Helen
dS=(dU/T)=8(π^5)k/15(h^3)(c^3)   4((T^3)/T)dT=((8(π^5)(k^4))/(15(h^3)(c^3)))4(T^2)dT

so ((8(π^5)(k^4))/(15(h^3)(c^3))) 4((T^3)/3)=4U/3T

Further: k=1.38e-23 kg (m^2) (s^(-2)) (k^(-1)), h=6.626e-34(m^2)  kg /s, 3e8 m/s

so S=(((272)(1.38e-23))/((((6.62e-34))^3)(((3e8))^3)))(T^3)k=90(T^3)k

Sun

We use (Subscript[M, ⊙])=2.0×(10^33)g, (Subscript[R, ⊙])=7.0×(10^10)cm, (Subscript[T, ⊙])=5800 K, (Subscript[L, ⊙])=3.9×(10^33)erg (s^(-1)). ρ~100, (Subscript[T, 6])~10=>(Subscript[P, g])/(Subscript[P, r])~100.

sM=2.0 (10^33);sR=7.0 (10^10);sT=5800;sL=3.9 (10^33);4π (sR^2)σSB (sT^4)

Q2 .4 Estimate the mean density of the sun, the surface gravity, the escape velocity, and the flux of radiation at Earth. Ignore solar rotation.

Helen
ρ=(M/((4/3)((πR)^3)))=((2e30kg)/((4/3)((π700000m)^3)))=1.4e12kg/(m^3)

v=(Sqrt[2GM/r]) and G=6.7e-11(m^3)(kg^(-1))(s^(-2))
so v=2e7m/s

The flux for a frequency is (Subscript[f, ν])=(∫_0^(π/2))(Subscript[I, ν])cos(θ)dΩ=(Subscript[I, ν])π

The energy denity for a black body is (Subscript[u, ν])=((8((πν)^2))/(c^3))   ((hν)/(exp(hν/kT)-1))and (Subscript[I, ν])=((Subscript[cdu, ν])/(dΩ))=((2((hν)^3))/(c^2))(1/(exp(hν/kT)-1))
so total flux is
f=(∫_0^∞)(Subscript[f, ν])dν=((σT)^4) where σ=((2(π^5)(k^4))/(15(c^2)(h^3)))=5.7e-5(erg/(s (cm^2) (k^4)))

f=(6.7e-5)(5800)(erg/(s (cm^2)))=.3306(erg/(s (cm^2)))

Radiative Transfer(4/15,4/18)

Foundations

Equations of Radiative Transfer

Flux (Subscript[F, ν])=∫(Subscript[dΩcosθI, νΩ])  = (Subscript[πI, νΩ]) when isotropic in hemisphere. Pressure is (Subscript[P, ν])=(1/c)∫((dΩcos)^2)(Subscript[θI, νΩ])  = ((4π)/(3c))(Subscript[I, νΩ]) =((Subscript[U, ν])/3)when isotropic.
(Subscript[dI,  νΩ])= -(Subscript[κρdsI, νΩ]) +(Subscript[j, νΩ])ρds  where κ is opacity and (Subscript[j, νΩ]) is the emissivity. Conventional to treat scattering as absorption and re-emission.
((Subscript[dI,  νΩ])/(dτ))= -(Subscript[I, νΩ]) +((Subscript[j, νΩ])/κ) where τ=∫dsκρ is the optical depth. When the atoms, ions and electrons are maintained in LTE at temperature T, Kirchhoff’s law tells us that ((Subscript[j, νΩ])/κ)=(Subscript[B, ν])(T).  
If we make the diffusion approximation, with mfp <<R and isotropic re-emission,
(Subscript[F, ν])=- (c/(3κρ))((∂(Subscript[U, ν]))/(∂r))=-(c/(3κρ))((∂(Subscript[B, ν]))/(∂T))((∂T)/(∂r)). The flux is then F=-((4(acT^3))/(3(Subscript[κ, R])ρ))((∂T)/(∂r)), where ((Subscript[κ, R])^(-1)) =((∫((dν)/κ)((∂(Subscript[B, v]))/(∂T)))/(∫dν((∂(Subscript[B, v]))/(∂T))))is the Rosseland mean opacity. Hence L=-((16((πacr)^2)(T^3))/(3(Subscript[κ, R])ρ))(dT/dr). This allows for the photons that have longer mfp to carry more flux.

Q2.5 Suppose that κ=((Kν)^n).Evaluate (Subscript[κ, R]).

Helen

(1/(Subscript[k, R]))=(∫_0^∞)(k^(-1))udν/(∫_0^∞)udν and u=((8((πhν)^3))/(c^3))(1/(exp(hν/kT)-1))

(∫_0^∞)udν=((8(π^5)(k^4)(T^4))/(15(c^3)(h^3)))

(∫_0^∞)(k^(-1))udν = (1/k) ((8πh)/(c^3))((ν^(3-n))/(exp(hν/kT)-1))=(1/k)((8πh)/(c^3))((((kT/h)))^(4-n))Γ(5-n)ζ(5-n)

Here Γ(5-n)=(4-n)! is the gamma function and ζ(5-n) is the Riemann zeta function

ζ(1)=∞
ζ(2)=(π^2)/6
ζ(3)=1.20205...
ζ(4)=(π^4)/90
ζ(5)=1.036928...
ζ(6)=(π^6)/945
and so on...

Opacity

Einstein Relations

For two - level atoms in equilibrium, (Subscript[R, u->l])=(Subscript[N, u])A(1+n) =(Subscript[R, l->u])=(Subscript[N, l])An. The absorption coefficent must contain a correction for stimulated emission κ=(Subscript[κ, a])(1-(e^(-hν/kT))) ~(Subscript[κ, a])((hν)/kT).

Thomson Scattering

(Subscript[κ, T]) = ((Subscript[σ, T])/((Subscript[μ, e])(Subscript[m, p])))= 0.35(cm^2)(g^(-1))for solar composition.

Q2.6 Show that the luminosity of an electron scattering atmosphere is limited by ~(10^38)(M/(Subscript[M, ⊙])) erg (s^(-1)).

Helen 
(Subscript[F, grav])=(GMm/r)

We need the outward radiation pressure:

(Subscript[F, rad])=dp/dt=(Subscript[R, scat])hν/c where hν/c is momentum and (Subscript[R, scat])=nσc is the Thomson scattering on an electron.  n=(L/(4((πr)^2)chν)) is the density of photons at r with energy hν

(Subscript[F, rad])=((Lσ)/(4((πr)^2)c))

(Subscript[F, rad])=(Subscript[F, grav]) when 

L=((4πcGMm)/σ)=(((4π)(3e10)(6.7e-8)(2e33)(1.7e-24))/(6.7e-25)) (M/(Subscript[M, sun]))=(10^38)(erg/sec)(M/(Subscript[M, sun]))

Free - free Opacity

At low frequency, where quantum effects unimportant, electron oscillation energy (((e^2)(E^2))/((mν)^2)) gained every collision. κρc~-((dln((E^2)))/dt)~ (Subscript[n, e])(((Subscript[n, i])(((Ze^2)/kT)))^2)((((kT/(Subscript[m, e]))))^(1/2))((e^2)/((Subscript[m, e])(ν^2))) => (Subscript[κ, a])~(((Z^2)/A))((((Subscript[n, e])(e^6))/((((Subscript[m, p])((Subscript[m, e])kT))^(3/2))(ν^2)c)))((kT/(hν))). 
So, (Subscript[κ, Rff])~ <((Z^2)/A)>((((ρe)^6)(h^2))/(((Subscript[m, p])^2)((((Subscript[m, e])^(3/2))(kT))^(7/2))c))~ 70 ((Subscript[ρT, 6])^(-3.5)) (cm^2)(g^(-1))for solar composition. (Subscript[j, νΩ])~((2 ((hν)^3))/(c^2))(Subscript[κ, a]) (e^(-hν/kT))~(T^(-1/2))(e^(-hν/kT)) and j~(ρ^2)(T^(1/2)).

Bound - free opacity

Ionize not excite; dominated by metals as H ionized. σ~(((α^3)((χ/(hν))))^3)((((hc/χ)))^(2 )), ν>χ/h, typically ν~χ/h ~kT/h   (Subscript[κ, Rbf])~1.4×(10^(3 ))((Z/0.02))((Subscript[ρT, -6])^(-3.5)) (cm^2)(g^(-1))~20((Z/0.02))(Subscript[κ, Rff]) very approximately. Should use tables.

Bound - bound opacity

Opacity due to excitation of bound states. Important in envelopes of cool stars, especially with molecules and radicals.

Conduction

Electrons mfp~(([(ρ/(Subscript[m, p]))(((((e^2)/kT)))^2)])^(-1))<<photon mfp. Ignore on main sequence. However, very important in WD, NS as degeneracy increases mfp. Keeps cores isothermal.

Opacity Tables

Tables (from Deep) at http://opalopacity.llnl.gov/existing.html
 We use a uniform composition (X=.70,Y=0.28,Z=0.02) from Table 73. The data is organized by {log T, log ρ, log κ}. We use this table for linear interpolation, κ(ρ,(Subscript[T, 6])).

SetDirectory["/Users/rdb/Desktop/teaching/ph361/2011"];
κtxt=ReadList["kappa.txt",Number,RecordLists->False,RecordSeparators->","];
κlst1=Table[{κtxt[[3n-2]],κtxt[[3n-1]],κtxt[[3n]]},{n,1,Length[κtxt]/3}];
κlst=Sort[Join[κlst1,Flatten[Table[{3.7+.1i,1+.5j,Interpolation[κlst1,{3.7+.1i,1+.5j},InterpolationOrder->1]},{i,1,36},{j,1,2}],1],{{7.2,1,Interpolation[κlst1,{7.2,1},InterpolationOrder->1]}}]];

Opacity Function

κ[ρ_,t_]:=(10^(Interpolation[κlst,{Log[10,t]+6,Log[10,ρ]},InterpolationOrder->1]))
ListContourPlot[κlst,Contours->{-1,-.5,0,.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6},FrameLabel->{"Log T","log ρ","Log κ",""},RotateLabel->False]

The Sun

Estimates

U∼(M/(Subscript[m, p]))kT~(10^48)erg  ∼((GM^2)/R) (Virial theorem). (Subscript[τ, Kelvin - Helmholtz]) ∼ (U/L)∼10Myr <<(Subscript[τ, solar system])∼5Byr. =>nuclear reactions.
(Subscript[l, mfp])∼  (1/(κρ))∼ 3×(10^(-3))((Subscript[T, 6])^3.5)(ρ^(-2))cm  ∼0.1 cm; (Subscript[τ,  escape]) ∼((R^2)/((Subscript[l, mfp])c))~(10^(5 ))yr; ∼ (Subscript[τ, KH])((Subscript[U, rad])/U)as it should. This is very rough!

Q2.7 Prove the Virial Theorem for a γ=5/3 gas to show that the total energy is minus the internal energy.

Simple Model

As an example, let us guess that the density has a Gaussian profile, terminated at the photosphere where ρ=(10^(-17)). ρ=100(e^(-(((r/1.53))^2))). Measure lengths henceforth in units of (10^10)cm, masses in (10^30)g pressure in (10^15)dyne (cm^(-2)), and temperature in (10^6)K. Other quantities in cgs.

rs=7;ms=2000;ρ=100 (E^(-(((r/1.53))^2)));
Plot[ρ,{r,0,7 },Frame->True,PlotRange->{0,120},FrameLabel->{"(Subscript[r, 10])","ρ"},RotateLabel->False]

Compute the mass and gravity.

nds1:=NDSolve[{m'[r]==4π (r^2)ρ,m[0]==0.},m,{r,0,rs}];mass=nds1;g=((6.67 (10^(-8))(10^10)  Evaluate[m[r]/.mass])/(r^2))[[1]];
Plot[{Evaluate[m[r]/.mass],.01 g},{r,0,rs},Frame->True,PlotRange->{0,2500},FrameLabel->{"(Subscript[r, 10])","(Subscript[m, 30]),(Subscript[g, 2])"},RotateLabel->False]

Compute the pressure and temperature. Observe that the temperature is unrealistic in the outer parts of the star.

nds2:=NDSolve[{p'[r]==-(10^(-5))ρ g,p[rs]==.001},p,{r,.001,rs}];press=nds2;t=((Evaluate[p[r]/.press])/(.14ρ));
Plot[{Evaluate[p[r]/.press],10t},{r,.001,rs},Frame->True,PlotRange->{0,300},FrameLabel->{"(Subscript[r, 10])","(Subscript[P, 15]),(Subscript[T, 5])"},RotateLabel->False]

Q2.8 If we were going to correct our original density profile to ensure that the temperature would fall monotonicallhy with radius in the outer parts of tyhe star, should the density gradient fall off more or less rapidly with radius?

Convection

Schwarzschild Criterion

Adiabatic interchange of fluid elements unstable if upper one lighter/buoyant  than surroundings ie if entropy decreases. If μ constant, then ((Pρ)^(-5/3)), or ((Tρ)^(-2/3))or ((T^(5/2))/P)cannot decrease. Convection very efficient except at surface. So entropy constant in convection zone. (nb beware ionization zones).

Q2.9 Explain what happens to the temperture gradient in an ionization zone.

Mixing length theory

Heat flux F~(Subscript[C, P])ρvΔT, where v∼(((ℓgδρ/ρ))^(1/2)), where ℓ is the mixing length. Numerical experiments suggest ℓ∼ -(1/4)(dr/(dlnP )). So, F∼(Subscript[C, p])((ρ(g/T))^(1/2)) (((Δ|∇T|))^(3/2))(ℓ^2). Hence
((Δ | ∇T |)/(| ∇T |))∼(M^2)∼((((F/(Subscript[Pc, s]))))^(2/3))

Solar Atmosphere (Gray)

Grey Atmosphere

Constant κ. (Subscript[I, νΩ])->I etc. μ(dI/(dτ))=I-S, where the source function S is the Rayleigh-Jeans spectrum if LTE and optical depth is measured downwards and μ≡cosθ. Make moment expansion
(dF/(dτ))=cU-4πS, c (dP/(dτ))=F, If F const, P=((F(τ+k))/c)  Truncate with Eddington approximation U =3P=3(F/c)(τ+k)=((4π)/c)S
Integrating, I(0,μ)=(Overscript[(∫_0)((dτ)/μ)(e^(-τ/μ)), ∞])S and F=2π(Overscript[(∫_0), 1])dμμI(0,μ) = (3/2)(Overscript[(∫_0), 1])dμ((μ^2)+μk)F =>k=(2/3), U=3(F/c)(τ+(2/3)). 
S=(c/(4π))(aT^4)=(σ/π)(T^4) =>T=(((Subscript[T, e])((3/4)τ+(3/(2 ))))^(1/4)). Effective temperature is that at τ=2/3. Can use as surface boundary condition. Asymptotic temperature is 0.84(Subscript[T, e]). Model only approximate because κ declines with ν and Eddington closure inaccurate. Modern computational methods solve the PDE exactly for arbitrary opacity. Monte Carlo methods often used. Non-LTE, non-planar, non-static, non-stationary, non-photon, non-non-relativistic extensions can be important!

Limb-darkening

I(0,μ)=((3F)/(4π))(μ+2/3).  Observe lower T at limb. Qualitatively correct. nb have chromosphere and corona and ultimately get limb-brightening!

Continuum Opacity

Ly edge -91nm Ba edge 364nm, Pa edge 820nm... Solar photosphere too cool to have much hydrogen ionized. Indeed too little excited to n=3 to contribute much opacity. Most of the electrons come from metals. Most visual continuum opacity from (SuperMinus[H]) ion with χ=0.75 eV, σ ∼ (10^(-17))(cm^(-2)). (Subscript[n, (SuperMinus[H])])∼(10^10)(cm^(-3)).

Photospheric Conditions

(Subscript[n, H])∼(10^(-17)), ρ∼2×(10^(-7,)) T∼6000, P~2×(10^5), g∼3×(10^4). Combining  equations of radiative transport and pressure balance, P∼g/κ => κ∼0.15, which agrees with table and scaleheight H∼1/κρ~300km.

Absorption Lines

Mostly see dark (Fraunhofer) absorption lines from cooler regions above photosphere.
QM and semi-classical calculation gives Lorentzian/Breit-Wigner/Cauchy natural line profile σ(ν)=((ΓKf)/(2π(((δν)^2)+(Γ^2)/4)))where Γis FWHM f is the oscillator strength and K = ∫dνσ = (((πe)^2)/((Subscript[m, e])c)) =0.026(cm^2)(s^(-1)). Typical natural width of permitted line is ∼(10^(-12 ))cm.

Line Profiles

Conventionally express using equivalent width (Subscript[W, λ])=∫dλ(1-((Subscript[F, λ])/(Subscript[F, c]))). 
Doppler broadening due to thermal motions.   Width is 4×(10^(-5))((((Subscript[T, 4])/A))^(1/2))λ. Convolution with Lorentzian gives Voigt profile. 
Pressure broadening due to individual collisions interupting emission over decay time. Also produces Lorentzian profile.
Linear stark broadening due to electric field of many particles breaking degeneracy of l-levels in H - ΔE∝(r^(-2)). 
Quadratic Stark broadening due to  many particles creating fluctuating polarizing field ΔE∝(r^(-4)). 
For molecules rely on van der Waals forces, ΔE∝(r^(-6)). 
Also get line broadening due to “microturbulence” motion in photosphere.
Curve of (Subscript[W, λ])vs N, the column density of absorbers gives the curve of growth.

Q2.10 Justify the expressions for the natural and thermal line widths

Q2.11 Make a sketch of a generic curve of growth usiong the information above.

Solar Spectrum

Q 2.12 Why does this spectrum seem to emphasize the infrared? Sketch what it would look like if you plotted dF/dlogν.

Q 2.13 Try to identify as many lines as you can. Compare with a hot and cool star spectrum.

Energy Generation(KW18, HK6, Bahcall3)

Nuclear Physics

The binding energy per nucleon curve peaks at Fe at 8.5MeV, but is already 7MeV ∼ (10^(-5))  erg at He. It is reasonable to suspect that the most abundant element H burns into He inside main sequence stars. For the sun n∼(10^26)(cm^(-3)) , v∼300 km (s^(-1)), If the sun burns 0.1 (Subscript[M, ⊙]),  t ∼ ((10^(56-5))/(4 (10^33)))∼(10^(10 ))yr. Sun is half way through its main sequence life. σ~(1/nvt) ∼ (10^(-51)) (cm^2)∼(10^(-25)) (fm^(-2)) .  Suppression due to Coulomb repulsion. (Subscript[r, min])∼(e^2)/kT ∼1000fm.  (Subscript[λ, deB])∼300fm.

Tunneling

S-wave Schrodinger equation: (ψr)''=(k^2)ψr, k=((([2μ(V-E)])^(1/2))/ℏ). So, | ψr(|^2)∼(v^(-1))(e^(-2πη)); where 2πη=2((((2μ))^(1/2))/ℏ)∫((dr(((e^2)/r)-E))^(1/2))=4α (c/v) (Overscript[(∫_0), 1])((dr((r^(-1)) - 1))^(1/2)) = 2πα (c/v), where v is the relative velocity.

Integrate[((((r^(-1))-1))^(1/2)),{r,0,1}];

σ ∝ ((Subscript[λ, deB])^2) (e^(-2 πη))=((S(E))/E)(e^(-2 πη)) ; <σv> ∼((((2/π)))^(1/2))((((μ/kT)))^(3/2)) ((2S)/μ) ∫(dvve^(-(2πη+((μv)^2)/2kT)))   Exponent minimized when v=(Subscript[v, 0])=(((((2 παckT)/μ)))^(1/3))=600((Subscript[T, 6])^(1/3)) km (s^(-1)) and has value 34((Subscript[T, 6])^(1/3)).  Only suprathermal protons penetrate the barrier. Hence <σv>∼ (10^(-15))S(E) ((Subscript[T, 6])^(-2/3))(e^(-34((Subscript[T, 6])^(-1/3)))). When (Subscript[T, 6])~15, the suppression is ~(10^(-6 )). This treatment ignores many subtleties including effects of screening. The overall cross section measured by S(E) is very small as the process is weak.

PP Chains

Only stable intermediary is d. Made by p+p->d+(SuperPlus[e])+ν+0.42+1.02-0.26=1.2MeV. (An alternative is p+e+p->d+ν but this is unimportant in the sun.) 
Then d+p->(He^3)+γ  +5.5MeV, (He^3)+(He^3)->α+2p+12.9MeV. 
Adding,  4 p + 2 e -> α + 2 γ+2ν+(1/4)(2×1.2+2×5.5+12.9)=6.6 MeV/p.
Also PPII: (He^3)+α->(Be^7)+γ+1.6MeV, (Be^7)+e->(Li^7)+ν+0.86-0.80MeV, (Li^7)+p->2α+17.3MeV. PPIII: (Be^7)+p->(B^8)+γ,(B^8)->(Be^(8* ))+(SuperPlus[e])+ν, (B^(8* )) -> 2 α. ((Subscript[E, ν])<14MeV, <(Subscript[E, ν])≥7MeV.)  
PPI dominates (Subscript[T, 6])<15  - most of sun. PPIII dominates at high T.  (10^(-4 )) of power for sun.

Q2.14 What happens to primordial deuterium in the solar core?

CNO Cycle

At higher temperatures, Coulomb barriers of CNO can be overcome and they provide catalyzed routes to burning hydrogen. 
(C^12)+p->(N^13)+γ, (N^13) -> (C^13)+ (SuperPlus[e]) + ν, (C^13) + p -> (N^14) + γ, (N^14) + p ->(O^15)+γ, (O^15)->(N^15)+(SuperPlus[e])+ν, (N^15)+p->(C^12)+α. Also p-capture on (N^15), (O^16), (O^17).  
Cycle using existing CNO, sensitive to Z, T. Dominates (Subscript[T, 6])>18-20.

Helium Burning

At higher temperatures. α + α -> (Be^8), (Be^8)+α->(C^(12* )) enabled by resonance. Followed by more α capture to build up heavy elements.

Energy Generation

(Subscript[ϵ, ppI]) = 2.5 (10^6)((Subscript[T, 6])^(-2/3))(e^(-33.8((Subscript[T, 6])^(-1/3)))) ρ (X^2)erg (s^(-1))(g^(-1)),  (Subscript[ϵ, CNO]) = 8.7 (10^27)((Subscript[T, 6])^(-2/3))(e^(-152((Subscript[T, 6])^(-2/3)))) ρ XZ ((Subscript[T, 6])^16)erg (s^(-1))(g^(-1)). CNO dominates for (Subscript[T, 6])>17.

Plot[{2.5 (10^6)(t^(-2/3))(E^(-34(t^(-1/3)))),8.7 (10^27).02(t^(-2/3))(E^(-152(t^(-1/3))))},{t,10,20},Frame->True,FrameLabel->{"(Subscript[T, 6])","ϵ","PPI, CNO rates",""},RotateLabel->False]

Solar model(4/20,4/25)

Basic Equations

(dP/dm)=-(Gm/(4((πr)^4))); (dr/dm) = -(1/(4 ((πr)^2) ρ(P,T))); (dL/dm)=ϵ; (dT/dm)=-(GmT/(4((πr)^4)P))∇, where ∇=min[((3κLP)/(64((πσGmT)^4))),(dlnT/dlnP)] ignoring ionization, surface and composition effects as well as time dependence which can be important for low and high mass stars but is not an issue for the sun. Also important for vibrational instability eg Cepheids, Miras.

Boundary Conditions

r(0)=0, L(0)=0; P(M)=0;T(M)=0, for radiative envelopes. However, solar envelope is convective amnd use T(M)=(Subscript[T, e]). Treat atmosphere as isothermal/exponential P=g∫dzρ,  where ∫dzκρ=(2/3), g=(GM/(R^2))  so P = (2/3)(1 + ((dlogκ)/(dlogρ)))(g/(Subscript[κ, e])). If non-grey atmosphere, must us  (OverBar[ κ])=((∫(Subscript[dνL, ν])κ)/L)

Convection Zone

Should use mixing length theory in outer envelope where entropy decreases with radius. Express in terms ∇= (dlogT/dlogP)=(2/(5 )). Conserve sum of radiative and convective heat flux up to photosphere.

Integration

Shooting method:  Guess (Subscript[P, 0]), (Subscript[T, 0]) and integrate out. Guess R, L and integrate in. Match at interior point. Alternatively use relaxation method.

ZAMS Solar Model (Still under construction. Sorry.)

Sun

Mass in (Subscript[M, ⊙]); L in (Subscript[L, ⊙]);ρ  in g (cm^(-3)); r in (10^10)cm; P in (10^15)dyne (cm^(-2)); T in (10^6)K, κ in (cm^2)(g^(-1))interpolated from κlst. 
P=.14ρT; (Subscript[ϵ, ppI]) = 2.5 (10^6)ρ (X^2)((Subscript[T, 6])^(-2/3))(e^(-33.8((Subscript[T, 6])^(-1/3))))erg (s^(-1))(g^(-1))=6.8 (10^5)ρ ((Subscript[T, 6])^(-2/3)) (e^(-33.8 ((Subscript[T, 6])^(-1/3)))) (Subscript[L, ⊙])((Subscript[M, ⊙])^(-1));

mc=.001;mi=.8;tc=14;ρc=90;pc=.14 ρc tc;rc=(10^(-10))(((((3 2 (10^33)mc )/(4π ρc))))^(.333));lc=6.8 (10^5)ρc (tc^(-2/3))(E^(-34(tc^(-1/3))))mc;
nds=NDSolve[{p'[m]==-(((((2 (10^33)))^2)6.67 (10^(-8))m)/(4π (10^55)((r[m])^4))),r'[m]==((2 (10^3).14 t[m])/(4π ((r[m])^2)p[m])),t'[m]==-(((((2 (10^33)))^2)6.67 (10^(-8))m t[m])/(4π (10^55)((r[m])^4)p[m]))Min[.4,((10^(Interpolation[κlst,{Log[10,t[m]]+6,Log[10,p[m]]-Log[10,t[m]]-.854},InterpolationOrder->1])) ((3 1.95l[m](10^15)p[m])/(64π 5.67 (10^(-5))6.67 (10^(-8))m (10^24)((t[m])^4))))],l'[m]==((6.8 (10^5))/(.14 ))p[m]((t[m])^(-5/3))(E^(-34((t[m])^(-1/3)))),p[mc]==pc,r[mc]==rc,t[mc]==tc,l[mc]==lc},{p,r,t,l},{m,mc,mi}][[1]];

Simplifying and integrating out from the core,

mc=.001;mi=.8;tc=13.7;ρc=83;pc=.14 ρc tc;rc=7.80(((((mc )/(ρc))))^(.333));lc=6.8 (10^5)ρc (tc^(-2/3))(E^(-34(tc^(-1/3))))mc;

ndsout=NDSolve[{p'[m]==-((2120m)/((r[m])^4)),r'[m]==((22.3t[m])/( ((r[m])^2)p[m])),t'[m]==-((2120m t[m])/(((r[m])^4)p[m]))Min[.4,((7.69l[m]p[m])/(m ((t[m])^4)))(10^(Interpolation[κlst,{Log[10,t[m]]+6,Log[10,p[m]]-Log[10,t[m]]+.854},InterpolationOrder->1]))] ,l'[m]==4.9 (10^6)p[m]((t[m])^(-5/3))(E^(-34((t[m])^(-1/3)))),p[mc]==pc,r[mc]==rc,t[mc]==tc,l[mc]==lc},{p,r,t,l},{m,mc,mi}][[1]];
pplt=Plot[Evaluate[((p[m])/pc)/.ndsout],{m,mc,mi},Frame->True];rplt=Plot[Evaluate[.143r[m]/.ndsout],{m,mc,mi},Frame->True];
tplt=Plot[Evaluate[((t[m])/tc)/.ndsout],{m,mc,mi},Frame->True];lplt=Plot[Evaluate[l[m]/.ndsout],{m,mc,mi},Frame->True];
delplt=Plot[Evaluate[(( 7.69l[m]p[m])/(m ((t[m])^4)))(10^(Interpolation[κlst,{Log[10,t[m]]+6,Log[10,p[m]]-Log[10,t[m]]+.854},InterpolationOrder->1]))/.ndsout],{m,mc,mi},Frame->True]
Show[pplt,rplt,tplt,lplt,PlotRange->{0,1},FrameLabel->{"m/(Subscript[M, ⊙])","P/(Subscript[P, c]),r/(Subscript[R, ⊙]),T/(Subscript[T, c]),L/(Subscript[L, ⊙])"}]

Integrating in from the surface,

ms=1.;mi=.8;ps=4. (10^(-10));rs=6.;ts=6. (10^(-3));ls=.7;
ndsin=NDSolve[{p'[m]==-((2120m)/((r[m])^4)),r'[m]==((22.3t[m])/( ((r[m])^2)p[m])),t'[m]==-((2120m t[m])/(((r[m])^4)p[m]))Min[.4,((7.69l[m]p[m])/(m ((t[m])^4))) (10^(Interpolation[κlst,{Log[10,t[m]]+6,Log[10,p[m]]-Log[10,t[m]]+.854},InterpolationOrder->1])) ],l'[m]==0.,p[ms]==ps,r[ms]==rs,t[ms]==ts,l[ms]==ls},{p,r,t,l},{m,ms,mi}][[1]];
pplt=Plot[Evaluate[((p[m])/10)/.ndsin],{m,mi,ms},Frame->True];rplt=Plot[Evaluate[.143r[m]/.ndsin],{m,mi,ms},Frame->True];
tplt=Plot[Evaluate[((t[m])/tc)/.ndsin],{m,mi,ms},Frame->True];lplt=Plot[Evaluate[l[m]/.ndsin],{m,mi,ms},Frame->True];
Show[pplt,rplt,tplt,lplt,PlotRange->{0,1},FrameLabel->{"m/(Subscript[M, ⊙])","P/(Subscript[P, c]),r/(Subscript[R, ⊙]),T/(Subscript[T, c]),L/(Subscript[L, ⊙])"}]
ρplt=Plot[Evaluate[((p[m])/(.14t[m]))/.ndsin],{m,mi,ms},Frame->True];

Standard ZAMS Model

Following Solanki, Here is a standard solar model. The lower curves are for the zero age model, the upper curves are for the current sun with lower central hydrogen density and consequently higher temperature.

Logarithmic Temperature Variation

Approximate variation

sol={{0,0,13.7,160},{.2,1.2,10.,60.},{.4,1.7,8.,25},{.6,2.2,6.4,8.},{.8,2.9,4.6,2.5},{.9,3.6,3.5,.63},{1,6.2,.006,2 (10^(-10))}};
solp=Table[{sol[[i,1]],sol[[i,4]]},{i,1,Length[sol]}];pm=Interpolation[solp];
solpplt=Plot[((pm[m])/(pm[0])),{m,0,1},Frame->True,FrameLabel->{"m/(Subscript[M, ⊙])","P","Solanki",""},RotateLabel->False,PlotRange->{0,1}];
solr=Table[{sol[[i,1]],sol[[i,2]]},{i,1,Length[sol]}];rm=Interpolation[solr];
solrplt=Plot[((rm[m])/7),{m,0,1},Frame->True,FrameLabel->{"m/(Subscript[M, ⊙])","r","Solanki",""},RotateLabel->False,PlotRange->{0,1}];
solt=Table[{sol[[i,1]],sol[[i,3]]},{i,1,Length[sol]}];tm=Interpolation[solt];
soltplt=Plot[((tm[m])/(tm[0])),{m,0,1},Frame->True,FrameLabel->{"m/(Subscript[M, ⊙])","T","Solanki",""},RotateLabel->False,PlotRange->{0,1}];
solρplt=Plot[((pm[m]tm[0])/(tm[m]pm[0])),{m,0,1},Frame->True,FrameLabel->{"m/(Subscript[M, ⊙])","T","Solanki",""},RotateLabel->False,PlotRange->{0,1}];
Show[solpplt,solrplt,soltplt,solρplt]

Tests(4/27)

Solar Neutrinos

Test proposed in 1960’s by Bahcall et al. Davis experiment performed in Homestake mine S. Dakota to detect neutrinos that that can convert Cl in 600t (Subscript[C, 2])(Subscript[Cl, 4 ])into A which decays with half life ~30d. (Gallium and Cerenkov water detectors are also used. The most important flux is from  8B from <0.1(Subscript[R, ⊙]) in the energetically unimportant PPIII chain. (CNO neutrinos can also be important.) The measured rate was found to be 2 SNU (reactions per (10^36) target atoms) instead of 7.SNO used (Subscript[D, 2])O which also measures (Subscript[ν, μ]),(Subscript[ν, τ]) in NC reactions. Effect complicated by MSW effect for E>5MeV where electron induced dispersion mediates the rate of neutrino conversion. SNO and Kamionkande demonstrated neutrino oscillation - conversion to ν, τ neutroinos like a coupled oscillator.  Reactor (KamLAND) and atmospheric neutrino experiments are also important and lead to a choice of two mass hierarchies.

Q2.15 How many neutrinos, photons  will the sun radiate over its lifetime?

Helioseismology

Normal Modes

p modes (P ~ 5min) propagate in surface at sound speed, evanescent in atmosphere, g modes (P~160min) evanescent in convection zone, occasionally reported but never confirmed! Also f mode. Global p modes measure ((Subscript[c, s])^2)(r) and errors <(10^(-3)). Largest at tachocline.

Q2.16 Explain from first principles why the Period of solar p modes ~ 5 m.

Tomography

Also measure propagation times of disturbances due to solar flares (sunquakes) ion WKB limit. See back of sun. Explore vertical structure of sunspots. Sensitive to motion, magnetic field...

Rotation

Surface rotation varies with latitude 25-30d period from equator to pole. Interior variation on cones (not cylinders) down to tachocline and solid body rotation at intermediate value in core. Mixes composition. Source of solar dynamo. Reported variations?
http://www.youtube.com/watch?v=huDgiYHCGMs

Q2.17 Explain why the solar angular frequency swas expected to be constant on cylinders within the convection zone.

Asteroseismology

Just beginning. eg Rigil Kent - older sun. Illustrates that stars slow down with time.

Quadrupole moment

Measure oblateness - 10ppm. 
Perihelion advance of mercury -43” per century - contributed by solar quadrupole.
More generally apsidal motion in elliptical binary stars, proportional to mass quadrupole moment

Q2.18 Confirm the estimate of the solar oblateness.

3. Main Sequence

Hertzsprung - Russell Diagram (KW22)

Astronomy

Luminosity

m=-2.5log(F/(Subscript[F, 0])) in a given band or bolometrically. M is absolute magnitude which is apparent magnitude at 10pc. Vega(A0V, 2(Subscript[M, ⊙])2.5(Subscript[R, ⊙]), 40(Subscript[L, ⊙])10,000K, 8pc,  is roughly zero point for magnitdue scales 
Reddening a complication
Distances from parallax, membership of clusters galaxies, masses from binaries. 
Sun is 4.7 (bolometric) 4.8 in V (Yellow) and 5.5 in B(Blue) B-V measures temperature.
Luminosity Class ({{Ia, Supergiants}, {Ib, Supergiants}, {II, Bright Giants}, {III, Giants}, {IV, Subgiants}, {V, Dwarfs}, {VI, Subdwarfs}, {VII, White Dwarfs}}) 
Subdwarfs aredwarfs with  low metallicity and opacity.

Q3.1 What is the apparent magnitude of a (10^10)(Subscript[L, ⊙]) galaxy at z=0-.3

Q3.2 Show that the bolometric flux from a m=0 star is 2.5×(10^(-5))erg (cm^(-2))(s^(-1)) ,  and the luminosity of a M=0 star is 3×(10^35)erg (s^(-1)).

Helen
M=-2.5 Log(L/(Subscript[L, sun]))+(Subscript[M, sun])
(Subscript[L, sun])=3.8(x10^33)erg/s
(Subscript[M, sun])=4.83

If M=0 then 0=-2.5 Log(L/3.8(x10^33)erg/s)+4.83
L = 85.5*(3.8 x10^33)=3(x10^35)erg/s

m=0 is the flux from Vega so F=(L/(4((πd)^2)))=((40(Subscript[L, sun]))/(4((π(8pc))^2)))~2(x10^(-5))erg/(cm^(2 ))s

Reddening

In V-band, total extinction A(V)=(Subscript[m, V])-(Subscript[m, V])(0) ~((Subscript[N, H])/(2 (10^21)(cm^(-2))))=R(V)E(B-V) where R~3, and E(B-V) is selective extinction. Different types of dust in different environments.

Temperature

Effective temperature - Color in B-V - mass - luminosity - ({{Type, (Subscript[T, 3]), Lines, M, L, f}, {O, 25-40, HeII, 15-100, >30,000, small}, {B, 11-25, HeI, 2-15, 25-30,000, .0013}, {A, 7.5-11, H,Fe, 1.5-2, 5-25, .006}, {F, 6-7.5, CaII, Fe, 1.2-1.5, 1.5-5, .03}, {G, 5-6, CH, CN, 0.8-1.2, .6-1.5, .07}, {K, 3.5-5, Mols, 0.5-0.8, .1-.6, .12}, {M, 2-3.5, TiO, 0.08-0.5, <.1, .75}, {L, 1.3-2, Na, <0.08, □, □}, {T, .7-1.3, (Subscript[CH, 4]), <0.08, □, □}})Decimal subdivisions. Luminosity Class Sun is G2V.

Q3.3 Explain the relationship between B-V and temperature.

Helen
If we treat the star as a black body then (Subscript[λ, max])=3(x10^(-3))mK/T
A star with a smaller (Subscript[λ, max]) or one that is bluer will be hotter.  Roughly speaking, a hot star will have zero or negative B-V and a cold star will have positive B-V

Q3.4 Estimate the fraction of stars by spectral class in a magntiude limited sample

Q3.5 Estimate the Kelvin-Helmholtz lifetime for a 1000K brown dwarf.

Helen
The Kelvein-Helmholtz lifetime is given by (Subscript[τ, KH])=(U/L)=((3(M^2)G)/(5R)) (1/L)
For a brown dwarf we have M~.05(Subscript[M, sun])or 80(Subscript[M, J])=80(1.9(x10^27)kg)
L~(10^(-4))(Subscript[L, sun])=(10^(-4))(3.839(x10^26)J/s)
R~(Subscript[R, J])~71000000m
G=6.67(x10^(-11))(m^3)/kg (s^2)
(Subscript[τ, KH])=((3[80(1.9(x10^27))])^2)6.67(x10^(-11))/(5(71(x10^6))) (1/(3.839(x10^22)))
=3.3(x10^17)s=(10^10)yrs

ZAMS Composition

X=0.73, Y=0.25, Z=0.02 for sun

Stellar Populations

I Young, slow, disk high metals
II Old, fast, halo, low metals
III Primordial, very low metals

Main Sequence Models

Rough scaling L=(M^3.5), t=10(M^(-2.5))Gyr,  R=(M^(.8)) in solar units. Stars become more luminous and compact as they age.

ZAMS stars

Q3.6 How does the temperature scale with mass?

Helen
The Stefan-Boltzmann Law gives us L∝(T^4) and L∝(M^α)
 where α=3 for M>(Subscript[M, sun])
 and α=5 for M<(Subscript[M, sun])
 
 so (T^4)∝(M^3)or (T^4)∝(M^5)

4. Evolution

Late Evolution (5/4)

LowMass Stars(KW32)

M<0.8(Subscript[M, ⊙])

(Subscript[t, MS])   > 10Gyr. Most stars; Salpeter mass function N(>M)∝(M^(-1.35))

Q4.1 Use the scaling between mass and luminosity, to describe how the luminosity of a galaxy is contributed by stars of different mass.

~1(Subscript[M, ⊙])

Main Sequence  Contracting, isothermal, near degenerate core surrounded by H-burning shell with increasing luminosity on CNO cycle
Subgiant  Increasing (Subscript[P, shell]) causes the convective envelope to expand at constant L so that T ∝(R^(-1/2)). 
Giant  L ∝(R^2) at constant T until "tip" of giant branch where L>10000(Subscript[L, ⊙]), R > 100(Subscript[R, ⊙]). 
He Flash (Subscript[T, core]) ~ (10^8)K He burns quickly to C in degenerate shell without expansion. 
Horizontal Branch Core and envelope shrink ~100Myr. (Evolution complicated by slow, intense mass loss which is not well understood. Also mass transfer in close binaries.) 
Cepheids (young higher mass Pop I P>1d P-L relation) and RRLyrae (old lower mass Pop II P<1d)  pulsating stars reside here.
Supergiant C core contracts (Subscript[T, core])~3×(10^8)K, some C+α->O burning occurs envelope expands along "Asymptotic" Giant Branch (AGB)   
Hot C-O White Dwarf  Unstable He flashes and pulsational instability of the envelope eject mass to form planetary nebulae, illuminated by C+O+Ne WD. (Subscript[T, WD]) ~(10^5)K
Cool White Dwarf  cools over Hubble time
Nova WD: in binary system undergoes explosive H burning with mass loss and rapid increase in L to <(10^4)(Subscript[L, ⊙])followed by (~60d) decline; repeating
SN Ia - small minority build up mass to (Subscript[M, Ch])  and undergo deflagration -> detonation explosion, leaving no remnant, we think. Peak luminosity (5 ×(10^9)(Subscript[L, ⊙]) or (Subscript[M, V])∼-19.3) correlates with width of decay (~30d) - a measure of the 56Ni (which decays to 56Co   in 6d and then to 56Fe in 77d).  SN1a have been used as standard candle for cosmography.

Q4.2 Estimagte how much energy the sun will radiate on the main sequence, as a red giant and as a white dwarf

Globular Clusters

Test theory by observing isochrones. Density of stars ∝ time spent in phase. Turn off mass gives age as a function of metallicity. Note blue stragglers which may have been mass transfer binaries.

Q4.3 The luminosity of a star at the tip of the red giant branch is sometimes used as a distance indicator. Explain how this can be and outline some conerns you might have about the accuracy of this method.

~4(Subscript[M, ⊙])

Similar evolution to above except happens faster and core He burning is non-degenerate and steady; there is no asymptotic branch

Q4.4 Find out the main properties of 4 (Subscript[M, ⊙])nmain sequence; name some local examples. How long will they remin on hte MS?

High Mass Stars(KW31)

>8(Subscript[M, ⊙])

Core gets hot enough to burn C and then up α sequence to Si and Fe. degeneracy pressure  of (Subscript[M, Ch])core inadequate to support envelope (∼(10^14)cm) and get core collapse to nuclear density ->(10^8)cm, (10^(-6))ρ. γ((+)^56)Fe->13α+4n,  endothermic. 
Uncertainties due to mass loss, semiconvection (overstable pulsations driven by μ gradients that can mix shells). 
Energy release in heat and neutrinos during core collapse SN can expel out layers to core collapse SN. Radiant, kinetic, neutrino/binding energy  (10^(49,51,53)) erg
Type II hydrogen; Type Ib no H; Type Ic no He, progenitors of long GRBs
Movie showing 3D simulation of explosion of 15 (Subscript[M, ⊙])star in cm2.mpeg (Nordhaus et al 2010).
Current expectation is that 8-20(Subscript[M, ⊙])stars produce neutron stars; >20(Subscript[M, ⊙]) stars produce black holes. 
There are many underluminous and overluminous core collapse supernovae. Some may be attributable to pair production when the radiation becomes hot enough which may drive supernovae through violent pulsations. PPSN may leave no remnant.
Slides from Woosley.

Q4.5 Adoptiong the Salpeter mass function, make an estimate of the numbers of WD, NS and BH in our Galaxy. What worries do you have about these estimates?

Magnetars

A small minority of neutron stars are born with (10^(14-15))G field. Most neutron stars are believed to be born with modest spin periods. However, some are recycled as millisecond pulsars. A magnetar with a few ms spin period would have enough spin energy to drive a supernova explosion.

Q4.6 What fraction of the pressure inside a magnetar is magnetic?

Gamma Ray Bursts

GRBs comprise "Long" (>2s) and "Short" (<2s) bursts plus Soft Gamma Repeaters which are Magnetar "Novae".  Long bursts appear to be a rare type of core collapse supernovae, most plausibly involving formation of a black hole in a spinning core. Short bursts may be coalescing neutron star binaries. They create collimated relativistic outflows with Lorentz factors ~1000.

Q4.7 A Long GRB at redshift 5 emits (10^54)(((ΔΩ)/(4π)))  erg of energy in thte form of ∼100 MeV γ-rays in a beam with solid angle ΔΩ  in our direction. Estimate the number of photons collected by a gamma ray telescope with effective area 1000 (cm^2).

~100(Subscript[M, ⊙])

η Car and R136 are two prominent examples of stars that are argued to exceed 100 (Subscript[M, ⊙]). There is much mass loss before explosion.

Nucleosynthesis

Overview

H, He are mostly cosmological though additional He is created by H burning. Li, Be, B can be created by cosmic rays from C, N, O by spallation.
N is a by-product of the CNO cycle.  
α (C, O, Ne, Si) elements below the iron peak are created by further  burning just prior to core collapse in massive young stars. 
SN Ia create Fe group elements from older, low mass stars.

s, r -processes

Elements beyond Fe are mostly produced through two basic mechanisms.  
s-process During advanced evolution of massive stars when the rate of neutron capture is slow compared with time for β decays to occur. Accounts for about half the isotopes beyond Fe. n sources are mostly 13C+α(->^16)O+n, 22Ne+α(->^25)Mg+n.

r-process The opposite inequality obtains during a supernova explosion. The neutrons are produced rapidly especially by photodisintegration of 56Fe.

Q4.7 Look up the abundances of some of the isotopes shown in these diagrams and explain which of them are made by the s- and r-processes.

Early Evolution (5/16)

Interstellar Medium

Over ten percent of the disk is in the form of gas which is replensihed by mass loss, supernovae anmd novae. 
The interstellar medium is heated by supernovae and cooled by radiation
Most of the volume is in a hot (n~ .003(cm^(-3)),T~5×(10^5)K) phase.
Most of the mass in the disk is shared between atomic (HI) hydrogen  (n~1 (cm^(-3)), T ~(10^4)K) and molecular ((Subscript[H, 2]))(n~100(cm^(-3)), T~(10^2)K)phases. 
The characteristic pressure of these phases is ~(10^(-12)) dyne (cm^(-2)) and  is also exhibited by the cosmic rays ( ~10(()^(-9))  (cm^(-3))GeV protons), magnetic field (B ~3μG), interstellar starlight and the microwave background.

Molecular Clouds

The molecular gas is mostly found in ~ (10^(5-6))(Subscript[M, ⊙]) self-gravitating molecular clouds with  average density n∼(10^2)(cm^(-3))and peak density n∼(10^6)(cm^(-3)) . The closest is Orion. The characteristic collapse time should be (((Gρ))^(-1/2)) ∼10 Myr while they appear to survive for longer due to stellar feedback.  Molecular clouds are quite dusty and best observed in IR, mm and especially CO. Over 150 molecules have been identified.

Star Formation

Jeans’ Mass

A cloud collapses when its mass exceeds the Jeans' mass. This can be derived intuitively by observing that the dispersion relation for a sound wave is (ω^2)=(c^2)(k^2), where (c^2) is the square of the relevant sound speed, (dp/(dρ)),  isothermal or adiabatic. Just as with a plasma, except with opposite sign the effect of graviation will change this to (ω^2)=(c^2)(k^2)-4πGρ. When the gravitational time is shorter than the sound crossing time, the cloud will collapse. Conventionally we define a Jeans' length (Subscript[λ, J])=((((((πc)^2)/(Gρ))))^(1/2)) ∼0.4((c/(200(ms^(-1))))) ((((n/((10^3)(cm^(-3))))))^(-1/2))pc and a Jeans' mass (Subscript[M, J]) = (π/6)ρ ((Subscript[λ, J])^3)∼2((((c/(200(ms^(-1))))))^(3/2))((((n/((10^3)(cm^(-3))))))^(-1/2))(Subscript[M, ⊙]). (See BT for a more rigorous derivation.)

Four Stages of Star Formation

(a) There will be hierarchical fragmentation and roughly isothermal collapse until the fragments cannot cool and a photosphere forms. 
(b) Mass continues to accrete onto the star dependent in a poorly understood way on the ability of magnetic field to transport angular momentum. 
(c) The core temperature continues to rise until (Subscript[T, c])∼(10^6)K and D-burning occurs. This delays the onset of PPI and allows mass to continue to accrete as well as to flow out in a wind. 
(d) After the D is exhausted, H-burning begins and the star settles onto the MS. A circumstellar disk remains and planets will form.

Shu, Adams, Lizano(1987)

Hayashi Track

After burning D, a newly formed protostar will be quite bloated (R~100(Subscript[R, ⊙])), cool (T~3000K) and luminous (L∼1000(Subscript[L, ⊙])). It will be convective and collapse along the Hayashi track on the K-H timescale, ∼ a few Myr, radiating away its binding energy. The stars collapse initially at constant and then after H-burning commences, at constant L. The boundary of these tracks is the Hayashi limit. Stars should not be found to the right of it.

First Stars

The first (Pop III) stars formed in the early universe (z~10) will be devoid of metals (low opacity, no CNO). More massive stars up to ~ 400 (Subscript[M, ⊙]) can form and survive. They can evolve to undergo additional types of supernova explosion.

Woosley

5. Principles of Stellar Dynamics

Overview of Collisionless Stellar Systems

Stars and Dark Matter

The universe, clusters, groups and normal galaxies are dominated by dark matter; dwarf galaxies, galaxy cores and globular clusters less so. The Galaxy has a luminosity ∼2×(10^10)(Subscript[L, ⊙]),  mostly from an exponential disk with mass ∼5×(10^10)(Subscript[M, ⊙]) in a dark halo of mass ~(10^12)(Subscript[M, ⊙]) and size ∼100kpc (solar radius ~8kpc, solar speed ∼220 km (s^(-1))). The central BH mass is ∼4×(10^6)(Subscript[M, ⊙]). Our neighbor M31 is about twice our size, 700 kpc away and approaching us at 120 km (s^(-1)).

Dynamical Friction

The stars and dark matter rarely interact except through long range gravitational forces. To quantify this, note that a small deflection in a distant encounter with impact parameter b is (Subscript[δv, ⊥]) ~((2Gm)/bv). These deflections add stochastically and so the scattering rate is n∫2((πbdbv(((δv)/v)))^2)∼((8((πnG)^2)(m^2))/(v^3))ln(((Subscript[b, max])/(Subscript[b, min]))). For a self-gravitating cluster/galaxy containing N stars, NGm∼(v^2)R, where (Subscript[b, max])∼R; (Subscript[b, min])∼(Gm/(v^2)) and so the "relaxation" time for a cumulative π/2 collision is (N/(8 πlnN))(R/v). For the Galaxy, interior to the solar radius and ignoring the dark matter, N~(10^11), R∼8kpc, v∼220 km (s^(-1)) the collision time is (10^(17 ))yr; for a globular cluster with N~(10^5), v∼10 km (s^(-1)), R~ 3 pc, the collision time is several Gyr.  We will refine these estimates later but it demonstrates why stellar systems are usually regarded as collisionless in the first approximation. (A bigger problem is that many stars can evolve on relatively short timescales.)

Violent Relaxation

A much faster opportunity for changing the particle energies and orbits may arise during initial galaxy collapse and major mergers when the dynamical time and orbital period ∼(((Gρ))^(-1/2)) can be comparable and broad resonances can occur. This is thought to have created the stellar distributions in settled galaxies.

Potential Theory

Spheres (BT2.1, 2.2)

Poisson’s Equation

The first task is to solve Poisson’s equation, (∇^2)Φ=-∇·g=4πGρ for given ρ(x), whiuch may evolve. The simplest case is spherical symmetry.
g=-G((m(r))/(r^2))where m(r)=((∫_0)^r)4((πr)^('2))dr' ρ(r'). Note that only the interior mass contributes and the exterior field is that of the total mass concentrated at the origin - Newton’s theorems. 
So Φ(r) = -4 πG((∫_r)^∞)dr'((∫_0)^(r'))dr''(((((r'')/(r'))))^2)ρ(r’’) assigning the zero point at r=∞. Hence Φ(r)=-G[((m(r))/r)+((∫_r)^∞)((dm(r'))/(r'))].

Q 5.1 Prove and interpret this expression.

Gravitational Energy

W=(1/2)∫dmΦ, where the half avoids double counting.

Surface Density

If ρ=ρ(r), Σ(b) = (Overscript[(∫_(-∞)), (Underscript[, ∞])])dzρ =(Overscript[(Underscript[2∫, b]), ∞])drrρ(r)((((r^2)-(b^2) ))^(-1/2)). This is an Abel integral and the solution is given by multiplying by b((((b^2)-((r')^2)))^(-1/2)) and integrating over b to obtain ρ (r) = -(1/π) (Underscript[(∫^∞), r])((db((b^2)-(r^2)))^(-1/2))((dΣ)/db).

Q 5.2 Derive this expression and verify it for a NFW distribution, ρ∝(((r^(-1))(a+r))^(-2)). Abel integrals are common in stellar dynamics.

Q 5.3 Evaluate, the mass, potential energy, and circular velocity for the NFW distribution. (It is conventional to normalize this to the mass interior to the radius where the density is 200 times the contemporary critical density.)

Spheroids (BT2.5)

Oblate Spheroidal Coordinates

A useful solution to Laplace’s equation involves oblate spheroidal coordinates (r,z)->a e (cosh u sin v, sinh u cos v). The 2D metric is (ds^2)=(h^2)((du^2)+(dv^2)) where h=a e ((((cosh^2)u - (sin^2)v))^(1/2)). 
We can invert to obtain h=(((((((a^2)(e^2)-(r^2)+(z^2)))^2)+4(r^2)(z^2)))^(1/4)), u=(((cosh^(-1))((((a^2)(e^2)+(r^2)+(z^2)+(h^2))/(2(a^2)(e^2)))))^(1/2)), v=(((sin^(-1))((((a^2)(e^2)+(r^2)+(z^2)-(h^2))/(2(a^2)(e^2)))))^(1/2)),

e=.98;h[a_,r_,z_]:=(((((((a^2)(e^2)-(r^2)+(z^2)))^2)+4(r^2)(z^2)))^(1/4));
u[a_,r_,z_]:=ArcCosh[((((((a^2)(e^2)+(r^2)+(z^2)+((h[a,r,z])^2))/(2(a^2)(e^2)))))^(1/2))];v[a_,r_,z_]:=ArcSin[((((((a^2)(e^2)+(r^2)+(z^2)-((h[a,r,z])^2))/(2(a^2)(e^2)))))^(1/2))];

Laplace' s Equation

Laplace's equation separates and for Φ=Φ[u], (Subscript[d, u])(cosh u (Subscript[d, u])Φ)=0=> ((dΦ)/dM)=-(G/(a e)) × {(sin^(-1))e,(sin^(-1))(sech u)} , {interior, exterior} to a thin spheroidal shell of mass dM semi-major axis a and eccentricity e. The equipotentials are confocal ellipsoids with u=const.; the associated field lines are orthogonal hyperbolae with v=const. Note that on the shell, u=(sech^(-1))e, h=((a(1-(e^2)(sin^2)v))^(1/2)).

massplt=Plot[((((1-(r^2))(1-(e^2))))^(1/2)),{r,0,1},PlotStyle->{Green,Thickness[.01]},PlotRange->{{0,4},{0,4}},Frame->True,AspectRatio->Automatic,FrameLabel->{"r","z","Mass Shell",""},RotateLabel->False];
Φplt=Table[Plot[((((1-(r^2)/(a^2))((a^2)-(e^2))))^(1/2)),{r,0,a},PlotStyle->{Red,Thickness[.005]},PlotRange->{{0,4},{0,4}},Frame->True,AspectRatio->Automatic],{a,1.5,6,.5}];
gplt=Table[Plot[((((r^2)((Cot[v])^2)-(e^2)((Cos[v])^2)))^(1/2)),{r,((((1-(e^2)((Sin[v])^2))/(((Csc[v])^2)-(e^2))))^(1/2)),6},PlotStyle->{Thickness[.005]},PlotRange->{{0,4},{0,4}},Frame->True,AspectRatio->Automatic],{v,0,.5π,.1π}];
Show[massplt,Φplt,gplt];

Q 5.4 Verify that this solution does satisfy Laplace’s equation

Q 5.5 Find another coordinate transformation that admits separaable solutions to Laplace’s equation.

Potential

Introduce dM/da(a), so that integrating over similar mass spheroids, Φ =  (∫_0^((((r^2)+((z^2)/(1-(e^2)))))^(1/2)))(Subscript[da((( dΦ)/dM)), e])(dM/da)+ (∫_((((r^2)+((z^2)/(1-(e^2)))))^(1/2))^∞)(Subscript[da(((dΦ)/dM)), i])(( dM)/da).

dΦdm[a_,r_,z_]:=-((ArcSin[Sech[u[a,r,z]]])/(a e));dmda[a_]=((2a)/(((1+(a^2)))^2));
Φ[r_,z_]:=NIntegrate[dΦdm[a,r,z]dmda[a],{a,.01,((((r^2)+((z^2)/(1-(e^2)))))^(1/2))}]-((ArcSin[e])/e)NIntegrate[((dmda[a])/(a )),{a,((((r^2)+((z^2)/(1-(e^2)))))^(1/2)),∞}]
Φplt=ContourPlot[Φ[r,z],{r,0,4},{z,0,4},AspectRatio->Automatic,Contours->{-.8,-.6,-.4,-.2},ContourShading->False, ContourLabels->All,FrameLabel->{"r","z","Gravitational Potential and Field",""} ,RotateLabel->False];

Field

The field external to a shell is straightforwardly computed to be (dg/dM)= ((G sech u)/(a e h)). The field direction is given by (dr/dz)=tanhu tanv so that ((dg)/dM)=-(G/(h^2)) (tanhu sinv,cos v). The interior field vanishes.

dgdm[a_,r_,z_]:=-(({Tanh[u[a,r,z]]Sin[v[a,r,z]],Cos[v[a,r,z]]})/((h[a,r,z])^2));
g[r_,z_]:=NIntegrate[dgdm[a,r,z]dmda[a],{a,0,((((r^2)+((z^2)/(1-(e^2)))))^(1/2))}];
gplt=ListVectorPlot[Table[{{r,z},g[r,z]},{r,.001,4,.5},{z,.001,4,.5}],AspectRatio->Automatic];
Show[Φplt,gplt]

Field

The field external to a shell is straightforwardly computed to be (dg/dM)= ((G sech u)/(a e h)). The field direction is given by (dr/dz)=tanhu tanv so that ((dg)/dM)=-(G/(h^2)) (tanhu sinv,cos v). The interior field vanishes.

dgdm[a_,r_,z_]:=-(({Tanh[u[a,r,z]]Sin[v[a,r,z]],Cos[v[a,r,z]]})/((h[a,r,z])^2));
g[r_,z_]:=NIntegrate[dgdm[a,r,z]dmda[a],{a,0,((((r^2)+((z^2)/(1-(e^2)))))^(1/2))}];
gplt=ListVectorPlot[Table[{{r,z},g[r,z]},{r,.001,4,.5},{z,.001,4,.5}],AspectRatio->Automatic];
Show[Φplt,gplt]

Q 5.6 Confirm the transformation of g(u) to g(r,z).

Density

An element of distance along the field at the shell is ((a(1-(e^2)(sin^2)v))^(1/2))du. 
The thickness of the shell (between two similar ellipsoids is ds =(((1-(e^2)))^(1/2))(((1-(e^2)(sin^2)v))^(-1/2))da
The surface density in the shell is then, from Gauss’ theorem,  dΣ=(dM/(4((((πa)^2)(1-(e^2)(sin^2)v))^(1/2)))). 
Hence the density in the shell is ρ = ((dΣ)/ds) = (1/(4((((πa)^2)(1-(e^2)))^(1/2))))(dM/da)which is constant over the shell.

ρ[r_,z_]:=((dmda[a])/(4π (a^2)(((1-(e^2)))^(1/2))))/.a->((((r^2)+((z^2)/(1-(e^2)))))^(1/2));
ContourPlot[ρ[r,z],{r,0,4},{z,0,4},AspectRatio->Automatic,Contours->{.001,.0003,.0001,.00003,.00001},ContourShading->False,ContourLabels->All,FrameLabel->{"r","z","Density",""} ,RotateLabel->False];

Q 5.7 Confirm the four statements made above.

Orbits (BT3)

We can use an inefficient integrator to compute 3D orbits in arbitrary potentials.
If the density is uniform, the potential is harmonic and the angular frequency will be (((4πGρ))^(1/2)).  An elliptic orbit comprises two degenerate oscillators. 
We can generalize to axisymmetric and non-axisymmetric potentials. Nearly circular orbits can be handled perturbatively in the epicyclic approximation. They are generically non circular unlike classical epicycles! We will return to this when we consider rotating systems.

Φ=(1/2)((ξ^2)+(ψ^2)+(ζ^2));(*Φ=Log[(ξ^2)+(ψ^2)+(ζ^2)];*)gx=-(∂_ξ)Φ/.{ξ->x[t],ψ->y[t],ζ->z[t]};gy=-(∂_ψ)Φ/.{ξ->x[t],ψ->y[t],ζ->z[t]};gz=-(∂_ζ)Φ/.{ξ->x[t],ψ->y[t],ζ->z[t]};tmax=30.;sc=4.;
nds=NDSolve[{x''[t]==gx,y''[t]==gy,z''[t]==gz,x[0]==1,x'[0]==0,y[0]==2,y'[0]==1,z[0]==-1,z'[0]==-1},{x,y,z},{t,0,tmax}];
xyztplt=Plot[Evaluate[{x[t],y[t],z[t]}/.nds],{t,0,tmax}]
rzplt=ParametricPlot[Evaluate[{(((((x[t])^2)+((y[t])^2)))^(1/2)),z[t]}/.nds],{t,0,tmax},AspectRatio->Automatic,PlotRange->{{-sc,sc},{-sc,sc}},Frame->True,FrameLabel->{"R","z","Surface of Section",""},RotateLabel->False]
xyzplt=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.nds],{t,0,tmax},Boxed->False,AspectRatio->Automatic,PlotRange->{{-sc,sc},{-sc,sc},{-sc,sc}}]

Fluid Approximation (BT4.2)

Stars can be treated as a fluid and will satisfy fluid equation, nm ((d<v>)/dt)=-∇·P-nm∇Φ where the pressure tensor P=nm(<v⊗v>-<v>⊗<v>)which need not be isotropic.
If the distribution is spherically symmetric with no mean motions (P^(θθ))=(P^(ϕϕ)) =(P^tt)and ∇·P= (∂_r)(P^rr)+(2/r)((P^rr)-(P^tt)) . The equation of hydrostatic equilibrium then becomes
(∂_r)(((Subscript[nσ, r])^2)) + n((∂_r)Φ+2((((Subscript[σ, r])^2)-((Subscript[σ, θ])^2))/r))=0. 
This is easy to generalize to different coordinate systems and more general flows and pressure tensors, for example in the solar neighborhood.

x={r,θ,ϕ};g={{1,0,0},{0,(r^2),0},{0,0,(r^2)  ((Sin[θ])^2)}};ig=Inverse[g];aff=Simplify[Table[(1/2)Sum[ig[[i,l]]
(D[g[[l,j]],x[[k]] ]+D[g[[l,k]],x[[j]] ]-D[g[[j,k]],x[[l]] ]),{l,1,3}],{i,1,3},{j,1,3},{k,1,3}] ];press={{prr[r],0,0},{0,((ptt[r])/(r^2)),0},{0,0,((ptt[r])/((r^2)((Sin[θ])^2)))}};
divpress=Table[Sum[D[press[[i,j]], x[[j]]]+Sum[(aff[[i,k,j]]press[[k,j]]+aff[[j,k,j]]press[[i,k]]),{k,1,3}],{j,1,3}],{i,1,3}]

Kinetic Description (BT4)

A kinetic description is based on the Collisionless Boltzmann/Vlasov equation. Let f(p,q,t) be a distribution function in 6D phase space for a Hamiltonian system. Continuity becomes (∂_t)f+(∂_q)f(OverDot[q])+(∂_p)f(OverDot[p])= (∂_t)f + (OverDot[q])·(∂_q)f + (OverDot[p])·(∂_p)f=(Df/Dt)=0, using Hamilton’s equations. f is constant along a trajectory. Usually we set m=1 and write it as (∂_t)f + v·∇f -∇Φ·(∂_v)f=0.  We can add birth and death and collision terms to the rhs.

Q5.4 Derive the Jeans’equations from the Boltzmann equations by performing integrals like n=∫dvf etc

A constant of the motion in the {x,v} phase space satisfies C(x,v,t) =const. An integral satisfies I(x,v)=const. and can restrict the motion to a subspace or limited volume. Any function of the integrals solves the plain Vlasov equation - Jeans’ theorem. Functions of three integral usually suffice. It is often helpful to explore surfaces of section plotting points, say, where the star crosses the equator northward. This can identify non-analytic integrals

6. Spherical Stellar Systems

Globular Clusters (Paul Simeon 5/23)

Open clusters (OC)

Stellar type: Pop I stars
Location: galactic disk
Star count: 100 - 1000
Size: 1 - 10 pc
Median age: ~100 Myr (ongoing formation)
Number in Milky Way: ~100,000 (lower limit)

Globular clusters (GC)

Stellar type: Pop II stars
Location: galactic halo, distributed with spherical symmetry about our galaxy
Star count: (10^4) - (10^6), total mass ~ (10^4) - (10^6) (Subscript[M, ⊙]).
Median radius: ~10 pc
Median age: ~10 Gyr (formed long ago)
Number in Milky Way: 150 - 158 currently known with perhaps 10 - 20 more undiscovered
	Andromeda has ~500. M87 has ~13,000. 
Brightest GC in MW has (Subscript[M, V])≃-3, with ~(10^6) stars. The faintest have ~1000 stars, with (Subscript[M, V])≃-3. Most GCs have (Subscript[M, V])≃-7.
White dwarfs make up about 10% of the stellar population. Lots of pulsars too. (find estimate)
Core radius: radius at which the projected surface brightness of the GC is half its central value, ~1.5 pc, but can range from 1-10 pc.
Median radius: radius of the sphere that contains half of all the light emitted by the cluster, ~10 pc.
Tidal radius: radius within which the gravity of the GC dominates over the external gravitational field, ~35-50 pc.

Introduction for collisionless gravitating systems in steady state (imagine a globular cluster)

If we have enough stars, N≫1, the stars can be treated as a smooth fluid with a distribution function f(x,v,t) that gives the number of stars in a phase volume (d^3)x (d^3)v at a time t. However, don’t expect things to be physical when the volume in phase space or real space is too small to contain many stars. The distribution function of a fluid obeys the Boltzmann equation, or the Boltzmann transport equation,
((df(x,v,t))/(dt))=((∂f)/(∂t))+((∂f)/(∂x))·v+((∂f)/(∂v))·a
                    =((∂f)/(∂t))+((∂f)/(∂x))·v-((∂f)/(∂v))·∇ϕ=(Subscript[(((∂f)/(∂t))), collisions])=0.
The collision term on the RHS is zero for large systems of stars. Physical collisions are negligible, and gravitational collisions are significant only for t≥(Subscript[t, c])≈((N/(ln N)))(Subscript[t, dyn]), where (Subscript[t, dyn]) is the typical dynamical time scale of orbits in a smooth potential. If t≪(Subscript[t, c]), we can call this a collisionless system. The stars follow the gravitational potential, which is determined by 
(∇^2)ϕ=4π G ρ, and ρ(x,t)=m∫f(x,v,t)(d^3)v.

This equation for ρ assumes that all stars have the same mass m, which is completely untrue. One could expect different stellar populations to have different distributions, e.g., older, smaller Pop II stars form a halo while Pop I stars (which contain some young, large stars) are in the disk.

There are many possible solutions to these coupled equations, but let’s start with steady-state solutions, where f(x,v,t)=f(x,v).
The goal is to have f(x,v) and ρ(x,t) that are both consistent with each other and with a potential function. Modern computing allows one to solve these coupled equations for solutions. However, in the past people assumed a distribution function or a potential, and then calculated ρ from there. Here we will assume a distribution function that seems to work pretty well for our purposes.

Any orbit in a spherically symmetric potential necessarily allows four isolating integrals of motion: the energy per unit mass E=(1/2)(v^2)+ϕ (or equivalently the Hamiltonian) and the three components of angular momentum L. By Jeans theorem, any non-negative function of these integrals can serve as a distribution function of a stellar system. If our system is spherically symmetric in all of its properties, then f cannot depend on the direction of L but only its magnitude. So, any steady-state spherical system can be described by a distribution function f(E,L).

Binney & Tremaine and Padmanabhan use the convention of changing the sign of the potential, ψ=-ϕ+(Subscript[ϕ, 0]), which requires a shifted energy, ϵ=-E+(Subscript[ϕ, 0])=-(1/2)(v^2)+ψ. If we define our potential and energy this way, our solutions often take a form such that f or ρ are zero for ϵ≤0 or ψ≤0. Note that the new potential obeys (∇^2)ψ=-4π G ρ with the boundary condition that ψ->(Subscript[ϕ, 0]) as |x|->∞. I prefer to use the standard convention where particles are trapped in negative potential wells.

Systems with f(E)

The simplest models of stellar systems have f(x,v)=f(E)=f((1/2)(v^2)+ϕ).

King model

The King model is that of a truncated isothermal sphere. Close to the center of the cluster, where particles are tightly bound to the cluster, stars have a Maxwell-Boltzmann distribution of velocities, dn∝Exp[-((v^2)/(2(σ^2)))](d^3)v. This would correspond to a distribution function f(E)=((Subscript[ρ, 0])/(((2π (σ^2)))^(3/2)))Exp[(((E))/(σ^2))]. However, at large radii the weakly bound stars may be tidally stripped by a neighboring galaxy or non-uniform gravitational field (clusters are never in completely isolated environments). This tidal radius corresponds to a fixed energy, (Subscript[ϕ, 0])<0, at which a particle escapes the cluster, effectively reducing the escape speed. To thin out the stars at large radii, the modified distribution function takes the form  
f(E)=((Subscript[ρ, 0])/(((2π (σ^2)))^(3/2)))(Exp[(((E-(Subscript[ϕ, 0])))/(σ^2))]-1), E≤(Subscript[ϕ, 0])
f(E)=0, E>(Subscript[ϕ, 0])

(∇^2)ϕ=(1/(r^2))(d/(dr))((r^2)((dϕ)/(dr)))=4π G ρ 

Except for scale factors for the number of stars and the radius of the cluster, globular clusters appear to only differ in their central concentration, c=(Subscript[log, 10])(((Subscript[r, t])/(Subscript[r, 0]))). The core radius is defined as (Subscript[r, 0])=(Sqrt[((9(σ^2))/(4π G (Subscript[ρ, 0])))]), where (Subscript[ρ, 0]) is the central density. The tidal radius, (Subscript[r, t]), is defined by ϕ((Subscript[r, t]))=(Subscript[ϕ, 0]).

Clear[ρ,ρ0,σ,v,f,r,ϕ,ϕ0,gnewt];
σ=1;ρ0=1;ϕ0=-1;gnewt=0.5;
f=((ρ0)/(((2π (σ^2)))^(3/2)))(Exp[(((ϕ[r]+(1/2)(v^2)-ϕ0))/(σ^2))]-1);
ρ=Integrate[4π (v^2)f,{v,0,(Sqrt[2(ϕ0-ϕ[r])])},Assumptions->σ≥0]
(*Integrate[4π*(r^2)*gnewt*ρ ,{r,0,∞}]*)
(*(∫_0^∞)(1/(r^2))((∫_0^∞)4π*(r^2)*gnewt*ρ \[DifferentialD]r)\[DifferentialD]r*)
(*int1=Integrate[4π*(r^2)*gnewt*ρ,ϕ'[0],ϕ[0]==-10,ϕ[r],r,{r,0,∞}]*)
(*DSolve[{2r ϕ'[r]+(r^2)ϕ''[r]==4π*(r^2)*gnewt*ρ,ϕ'[0]==0,ϕ[0]==-10},ϕ[r],r]*)
(*NDSolve[{2r ϕ'[r]+r*r ϕ''[r]==4π*(r^2)*gnewt*(1/(3 (Sqrt[π])))(-6 DawsonF[(Sqrt[-1-ϕ[r]])]+2 (Sqrt[-1-ϕ[r]]) (5+2 ϕ[r])),ϕ'[0.01]==0,ϕ[0.01]==-10},ϕ[r],{r,0.01,10}]*)

(*Isothermal sphere*)
Clear[ρ,ρ0,σ,v,f,r,ϕ,ϕ0,gnewt];
σ=1;ρ0=1;ϕ0=-1;gnewt=0.5;
f=((ρ0)/(((2π (σ^2)))^(3/2)))Exp[(((ϕ[r]+(1/2)(v^2)-ϕ0))/(σ^2))];
ρ=Integrate[4π (v^2)f,{v,0,(Sqrt[2(ϕ0-ϕ[r])])},Assumptions->σ≥0]
phi:=NDSolve[{2r ϕ'[r]+r*r ϕ''[r]==4π*(r^2)*gnewt*ρ0 (Sqrt[ϕ0-ϕ[r]]) ((2/((Sqrt[π]) σ))-(((E^((-ϕ0+ϕ[r])/(σ^2))) Erf[((Sqrt[-ϕ0+ϕ[r]])/σ)])/(Sqrt[-ϕ0+ϕ[r]]))),ϕ'[0.1]==0,ϕ[0.1]==-10},ϕ[r],{r,0.1,100}];
Plot[phi,{r,1,10}]

Systems with f(E,L)

Relaxation time and the lifetime of a cluster

From the previous chapter, the “relaxation” time for a cumulative π/2 collision is (N/(8 πlnN))(R/v). For the Galaxy, interior to the solar radius and ignoring the dark matter, N~(10^11), R∼8kpc, v∼220 km (s^(-1)) the collision time is (10^(17 ))yr; for a globular cluster with N~(10^5), v∼10 km (s^(-1)), R~ 3 pc, the collision time is several Gyr (3 Gyr for these numbers). For small clusteres, the relaxation time could be just a few hundred million years. Over time scales comparable to the relaxation time, the high-velocity tail of the velocity distribution is repopulated, and those stars quickly escape the tidal radius. The cluster cools itself off this way and would eventually collapse. 

Some clusters appear to last longer than they should after many relaxation times because binary stars contribute a lot of energy to the cluster. A close binary star can have as much binding energy as the entire cluster, and up to half of the stars in a cluster are binaries. All of this latent energy can be added to the velocity of stars by stellar three-body collisions. This extra “heat” supports the cluster for much longer times.
For an estimate of the relative binding energies of binaries and clusters, we can use the virial theorem to get ((Subscript[BE, bin])/(Subscript[BE, gc]))∼(m/(N m))(((((Subscript[v, bin])/(Subscript[v, gc]))))^2). Velocities of binary stars can be up to 1000 km/s, so a single binary can provide significant kinetic energy to the stellar population.

Binaries in clusters

Binaries in clusters are classified as being hard or soft binaries if their binding energy |E|=((G (Subscript[m, 1])(Subscript[m, 2]))/(2a)) (where a is the semimajor axis of their obit) is large or small compared with the random kinetc energy of field stars, (1/2) (Subscript[m, f])(σ^2). Hard binaries are more likely to get harder and transfer energy to the field star, while soft binaries are more likely to get disrupted by a random encounter. The dynamics of hard and soft binaries should be treated differently.

Soft binaries

We treat soft binaries as a reduced mass orbiting a fixed mass. If a field star passes one star more closely than the other, the effect of the farther star is negligible. If the field star passes at a greater distance, then the center of mass velocity of the binary will change, not affecting the binding energy of the system. By the virial theorem, the magnitude of the kinetic energy of the stars of a soft binary is roughly the same as the binding energy, which is less than the kinetic energy of the average field star. As the system moves towards energy equipartition, the interactions with field stars will on average adding internal energy to the system. A soft binary may be disrupted by either ejection from a single interaction (called ionization) or by a series of more subtle interactions (called evaporation). Evaporation dominates for very soft binaries, and the two disruption processes are comparable when E≃m (σ^2).

It also turns out that in relaxed clusters the rate of evaporation of soft clusters is balanced by the formation rate from three-body interactions, resulting in an equilibrium number on the order of one soft binary per cluster, regardless of the number of stars in the cluster. Thus, soft binaries play no major role in the evolution of globular clusters.

Hard binaries

Hard binaries have binding energies |E|>(Subscript[m, f])(σ^2), and thus by virial theorem, the orbital velocities of the stars in a binary are much greater than the average passing star. When a field star interacts with a hard binary, either by exchange or flyby of the third star, the exiting star usually leaves with a speed comparable to the orbital speed; thus, the single star gained energy and the binary became harder. This is sometimes stated as Heggie’s law: on average, hard binaries become get harder and soft binaries get softer.

The rate of energy transfer from hard binaries becoming harder has been estimated from numerical studies to be <(OverDot[E])>≃-3.8 ((ρ (G^2)(m^3))/σ), equivalent to <ΔE>(Subscript[t, relax])≃-0.13 m (σ^2).

Compact objects can also capture stars if they pass close enough and at the right speed to resonate with stellar tidal oscillations. The internal energy dissipated in the stellar oscillation can reduce the kinetic energy of the star enough for it to be captured by the compact object. According to Padmanabhan, tidal capture is much more important than three-body encounters for forming binaries. Only when the number of stars is about 10 or fewer is the three-body process more dominant. B&T don’t seem to reach the same result, concluding the rate of tidal captures is uncertain.

X-ray binaries

1000 times more abundant in the densest GCs than in the Galaxy. These are compact objects accreting matter from a companion star. With stellar densities of GCs roughly 1000 times that of the galactic disk, it’s not surprising to find an abundance of these.

Millisecond pulsars