HW #4

1. isothermal halo

The isothermal model of halo assumes the distribution 
\(TraditionalForm\`ρ(\(x\& -> \), \(v\& -> \)) = 
        N\ \(exp(\(1\/σ\^2\) \((Ψ - \(1\/2\) \(v\& -> \)\^2)\))\)\).
This ansatz is called isothermal because of its resemblance to the Boltzmann distribution
\(TraditionalForm\`ρ(\(x\& -> \), \(v\& -> \)) = \(N\ \(exp(\(1\/\(m\ σ\^2\)\) \((m\ Ψ - \(m\/2\) \(v\& -> \)\^2)\))\) = N\ e\^\(\(-E\)/k\ T\)\)\)
with the identifications k T = m (σ^2) and the single-particle energy E = (m/2) (((Overscript[v,  ->]) )^2) - 
          m Ψ.  (The griavitational potential energy is V = m Ψ.)  But this is not just a resemblance.  Any phase distribution function given in terms of the Hamiltonian is automatically a static solution to the Boltzmann equation, and furthermore an exponential form (Boltzmann distribution in thermal equilibrium) is known to be a particularly robust solution under small perturbations.  Therefore, it is a good guess for a stable configuration of the halo.

The spatial density is given by the integration upon the velocities, 
\(TraditionalForm\`ρ(
          r) = \(∫\(d\^3\) 
              v\ N\ \(exp(\(1\/σ\^2\) \((Ψ - \(1\/2\) \(v\& -> \)\^2)\))\) = \(\((2  π\ σ\^2)\)\^\(3/2\)\) 
            N\ e\^\(Ψ/σ\^2\)\)\).
With the boundary condition ρ(0) = (Subscript[ρ, 0]), (((2  π (σ^2)))^(3/2)) 
          N (e^(Ψ(0)/(σ^2))) = (Subscript[ρ, 0]), and hence 
ρ(r) = (Subscript[ρ, 0]) 
          (e^((Ψ(r) - Ψ(0))/(σ^2))).  Namely,
Ψ(
          r) = Ψ(0) + (σ^2) 
            log ((ρ(r))/(Subscript[ρ, 0])).
The Poisson equation is then
\(TraditionalForm\`\(1\/r\^2\) \(d\/\(d\ r\)\) \((\(r\^2\) \(d\/\(d\ r\)\) \(σ\^2\) 
              log\ \(ρ(r)\)\/ρ\_0)\) = \(-4\) π\ G\ \(ρ(
            r)\)\).
Now using the variables \(TraditionalForm\`ρ\&~ = ρ/ρ\_0\), \(TraditionalForm\`r = \(r\_0\) r\&~\), \(TraditionalForm\`r\_0 = \@\(9  σ\^2/4  π\ G\ ρ\_0\)\),
\(TraditionalForm\`\(σ\^2\/r\_0\^2\) \(1\/r\&~\^2\) \(d\/\(d\ r\&~\)\) \((\(r\&~\^2\) \(d\/\(d\ r\&~\)\) 
              log\ ρ\&~)\) = \(-4\) π\ G\ \(ρ\_0\) ρ\&~\)
or
\(TraditionalForm\`\(1\/r\&~\^2\) \(d\/\(d\ r\&~\)\) \((\(r\&~\^2\) \(d\/\(d\ r\&~\)\) 
              log\ ρ\&~)\) = \(\(-r\_0\^2\)\ \(\(4  π\ G\ ρ\_0\)\/σ\^2\) ρ\&~ = \(-9\)\ ρ\&~\)\)
Writing out the derivatives,
\(TraditionalForm\`2\ \(1\/r\&~\) ρ\&~'\/ρ\&~ + ρ\&~''\/ρ\&~ - \((ρ\&~'\/ρ\&~)\)\^2 = \(-9\)\ ρ\&~\)
subject to the boundary condition \(TraditionalForm\`\(ρ\&~\)(0) = 1\).  We also impose \(TraditionalForm\`ρ\&~' \((0)\) = 0\)to avoid a \(TraditionalForm\`1\/r\&~\)singularity in the above equation at the origin.  Then we find \(TraditionalForm\`ρ\&~'' \((0)\) = \(-3\)\).

solution = 
      NDSolve[{2 If[x == 0, -3, y'[x]/x]/y[x] + (y''[x]/y[x] - (((y'[x]/y[x]))^2)) == -9 y[x], 
          y[0] == 1, y'[0] == 0}, y, {x, 0, 100}]

<< Graphics`Graphics`

LogLogPlot[Evaluate[y[x] /. solution], {x, 0, 100}, 
      PlotRange -> {0.00001, 1}]

LogLogPlot[(2/9) (1/(1 + (x^2))), {x, 0.001, 100}, 
      PlotRange -> {0.00001, 1}, 
      PlotStyle -> RGBColor[1, 0, 0]]

Show[%, %%]

Indeed, the asymptotic behavior is given approximately by \(TraditionalForm\`ρ\&~ = \(2\/9\) 1\/\(1 + r\&~\^2\)\).  

This behavior can also be understood analytically.  Going back to the differential equation 
\(TraditionalForm\`2 
             r\&~\ ρ\&~'\/ρ\&~ + \(r\&~\^2\) ρ\&~''\/ρ\&~ - \(\(r\&~\^2\)(ρ\&~'\/ρ\&~)\)\^2 = \(-9\)\ \(r\&~\^2\) ρ\&~\)
and assume the asymptotic behavior \(TraditionalForm\`\(ρ\& ~~ \) A\ r\&~\^\(-n\)\).  Keeping only the leading terms in \(TraditionalForm\`1/r\&~\),
\(TraditionalForm\`2 \((\(-n\))\) + \((\(-n\))\) \((\(-n\) - 1)\) - 
          n\^2 = \(-9\) A\ r\&~\^\(2 - n\)\)
If n < 2, there is only a trivial solution A = 0.  If n > 2, -2 n + n(n + 1) - (n^2) = 0 and hence n = 0, a contradiction.  Therefore, n = 2, and then A = (2/9).

The rotation speed of stars embedded inside the halo is determined by the usual balance between the gravitational force and the centrifugal force, 
\(TraditionalForm\`m\ v\^2\/r = \ \(\(d\ V\)\/\(d\ r\) = \(\(-m\)\ \(d\ Ψ\)\/\(d\ r\) = \(-m\)\ \(σ\^2\) \(d\ ρ/d\ r\)\/ρ\)\)\).  
Therefore,
\(TraditionalForm\`v\^2 = \(\(-σ\^2\) \(r\/ρ\) \(d\ ρ\)\/\(d\ r\) = \(-σ\^2\) \(r\&~\) ρ\&~'\/ρ\&~\)\).
Note that the asymptotic rotation speed is found using the asymptotic solution obtained above,
\(TraditionalForm\`v\_∞\^2 = \(\(-σ\^2\) \(\(r\&~\)(\(-2\))\) r\&~\^\(-1\) = 2  σ\^2\)\), and hence (Subscript[v, ∞]) = (Sqrt[2]) σ.

The fact that the rotation speed approaches a constant instead of falling as (r^(-1/2)) is the most surprising feature of the data ("flat rotation curve"), which is reproduced in this model.

Plot[
      Sqrt[Evaluate[-x y'[x]/y[x] /. solution]], {x, 0.001, 100}, 
      PlotRange -> {0, 2}]

2. Gravitational Microlensing

It is presented in a separate PDF because I couldn't insert figures into Mathematica notebook.