HW #5

1. Boltzmann equation

<< Graphics`Graphics`

We start with the Boltzmann equation
\(TraditionalForm\`\(d\ n\)\/\(d\ t\) + 
          3  H\ n = \(-⟨σ\ v⟩\) \((n\^2 - n\_eq\^2)\)\).
We defined the yield Y = n/s and rewrite the equation in terms of x = m/T instead of time.  Because the entropy density scales as (R^(-3)), it satisfies the equation
((d s)/(d t)) + 3  H s = 0
and hence the 3  H term cancels for the yield,
\(TraditionalForm\`\(d\ Y\)\/\(d\ t\) = \(-⟨σ\ v⟩\) \(s(Y\^2 - Y\_eq\^2)\)\).
To rewrite the time variable by x, we need the relationship between time and temperature in radiation dominated universe.  The expansion rate is 
\(TraditionalForm\`H\^2 = \(\(\(8  π\)\/3\) \(G\_N\) ρ\_R = \(\(8  π\)\/3\) \(G\_N\) \(g\_*\) \(π\^2\/30\) T\^4\)\).
Using \(TraditionalForm\`H = \(\(R\&·\)\/R = \(d\ R\)\/\(R\ d\ t\)\)\), and noting T ∝ (R^(-1)) and hence ((d R)/R) = -((d T)/T), we find d t = -(1/H) ((d T)/T).  Because H(T) = H(m) ((T^2)/(m^2)),
d t = -((m^2)/(H(m))) ((d T)/(T^3)).  Then the Boltzmann equation becomes 
\(TraditionalForm\`\(d\ Y\)\/\(d\ T\) = \(m\^2\/\(H(
                m)\)\) \(1\/T\^3\) ⟨σ\ v⟩ \(s(Y\^2 - Y\_eq\^2)\)\)
Using s(T) = s(m) ((T^3)/(m^3)), it becomes
\(TraditionalForm\`\(d\ Y\)\/\(d\ T\) = \(\(m\^2\/\(H(m)\)\) \(\(s(
                  m)\)\/m\^3\) ⟨σ\ v⟩ \((Y\^2 - 
                Y\_eq\^2)\) = \(\(s(m)\)\/\(H(
                  m)\)\) \(1\/m\) ⟨σ\ v⟩ \((Y\^2 - Y\_eq\^2)\)\)\).
Finally, d T = d (m/x) = -m ((d x)/(x^2)), and
\(TraditionalForm\`\(d\ Y\)\/\(d\ x\) = \(-\(1\/x\^2\)\) \(\(s(
                m)\)\/\(H(
                m)\)\) ⟨σ\ v⟩ \((Y\^2 - Y\_eq\^2)\)\)
\(TraditionalForm\`\(\(=\)\(\(-0.368\) \(\(g\_*\)\^\(1/
                2\)\) \(1\/x\^2\) 
          m\ \(M\_Pl\) ⟨σ\ v⟩ \((Y\^2 - Y\_eq\^2)\)\)\)\)
Here, (Subscript[M, Pl]) = (((8  π (Subscript[G, N])))^(-1/2)) = 
          2.44×(10^18) GeV.

N[((((1/3) ((π^2)/30)))^(-1/2)) 
        Zeta[3]/(π^2)]

To integrate the Boltzmann equation, we need to have an expression for the equilibrium yield.  Once the particle is non-relativistic, the difference in statistics is not important.  The number density is 
\(TraditionalForm\`n\_eq = \(∫\(\(\(d\^3\) 
                    p\)\/\((2  π)\)\^3\) 
              e\^\(\(-β\)\ E\) = \(∫\(\(\(d\^3\) 
                      p\)\/\((2  π)\)\^3\) 
                e\^\(\(-β\)\ \((m + p\^2/2  m)\)\) = \(\(e\^\(\(-β\)\ m\)\)(\(m\ T\)\/\(2  π\))\)\^\(3/2\)\)\)\).
The yield is
\(TraditionalForm\`Y\_eq = \(n\_eq\/s = \(\(1\/\(g\_*\)\) \(π\^2\/\(ζ(
                    3)\)\) \(1\/T\^3\) \(\(e\^\(\(-β\)\ m\)\)(\(m\ T\)\/\(2  π\))\)\^\(3/2\) = \(\(π\^2\/\(\(g\_*\) \(ζ(
                        3)\)\)\) \(\(e\^\(-x\)\)(x\/\(2  π\))\)\^\(3/2\) = 0.521\ \(\(g\_*\)\^\(-1\)\) \(x\^\(3/2\)\) e\^\(-x\)\)\)\)\).

N[((π^2)/(Zeta[3] (((2  π))^(3/2))))]

The boundary condition is such that Y = (Subscript[Y, eq]) at x = 0.  Unfortunately for x = 0 (or T ≫ m), the non-relativistic approximation we made to work out (Subscript[Y, eq])above is no longer good.  Because the result is not too sensitive to the initial condition, we set Y = (Subscript[Y, eq])at x = 1.  We will verify later that indeed the result is insensitive to this choice.  We use GeV unit for everything.

S-wave

Mathematica unfortunately seems to have trouble dealing with big numbers such as (Subscript[M, Pl]).  We can help it by solving for y = ((s(m))/(H(
                m))) ⟨σ v⟩ 
          Y
((d y)/(d x)) = -(1/(x^2)) ((y^2) - 
              ((Subscript[y, eq])^2)),
\(TraditionalForm\`y\_eq = \(\(\(s(m)\)\/\(H(
                  m)\)\) ⟨σ\ v⟩ \(1\/s\) \(\(e\^\(\(-β\)\ m\)\)(\(m\ T\)\/\(2  π\))\)\^\(3/2\) = \(\((\(1\/3\) π\^2\/30)\)\^\(\(-1\)/
                  2\)\) \(M\_Pl\) \(m\/T\^3\) ⟨σ\ v⟩ \(\(e\^\(\(-β\)\ m\)\)(\(m\ T\)\/\(2  π\))\)\^\(3/2\)\)\)
=0.192 (Subscript[M, Pl]) 
          m ⟨σ v⟩ (x^(3/2)) 
          (e^(-x))

N[((((1/3) ((π^2)/30)))^(-1/
              2)) ((((1/(2  π))))^(3/2))]

Mathematica still balks when I let it integrate all the way from x = 1 to 10000.  I break it up into 1 < x < 50 and the rest.

solution1 = 
      NDSolve[{y'[
                x] == -(1/(x^2)) (y[x]^2 - (((0.192 (Subscript[M, Pl]) m σ (x^(3/2)) (E^(-x))))^2)), 
            y[1] == 0.192 (Subscript[M, Pl]) m σ (1^(3/2)) 
                (E^(-1))} /. {m -> 1000, 
            g -> 100, σ -> (10^(-10)), 
            (Subscript[M, Pl]) -> 2.44 (10^18)}, y, {x, 1, 50}]

Evaluate[y[50] /. solution1]

solution2 = 
      NDSolve[{y'[
                x] == -(1/(x^2)) (y[x]^2 - (((0.192 (Subscript[M, Pl]) m σ (x^(3/2)) (E^(-x))))^2)), 
            y[50] == 42.144129023721824`} /. {m -> 1000, 
            g -> 100, σ -> (10^(-10)), 
            (Subscript[M, Pl]) -> 2.44 (10^18)}, y, {x, 50, 10000}]

LogLogPlot[Evaluate[y[x] /. solution1], {x, 1, 50}, 
      PlotRange -> {{1, 10000}, {1, (10^11)}}]

LogLogPlot[Evaluate[y[x] /. solution2], {x, 50, 10000}, 
      PlotRange -> {{1, 10000}, {1, (10^11)}}]

Show[%, %%]

Compared to the equlibrium values,

LogLogPlot[
      0.192 (Subscript[M, Pl]) m σ (x^(3/2)) 
          (E^(-x)) /. {m -> 1000, 
          g -> 100, σ -> (10^(-10)), 
          (Subscript[M, Pl]) -> 2.44 (10^18)}, {x, 1, 1000}, 
      PlotRange -> {{1, 10000}, {1, (10^11)}}, 
      PlotStyle -> RGBColor[1, 0, 0]]

Show[%, %%]

the solution is basically right on the equilibrium until x~20, and then becomes approximately constant afterwards.  This is why we don't expect the result to be sensitive to the initial condition as long as it starts on the equilibrium for x < 10 or so.  Verify this point by taking the ratio

LogLinearPlot[
      0.192 (Subscript[M, Pl]) m σ (x^(3/2)) 
          (E^(-x))/Evaluate[y[x] /. solution1[[1]]] /. {m -> 1000, 
          g -> 100, σ -> (10^(-10)), 
          (Subscript[M, Pl]) -> 2.44 (10^18)}, {x, 1, 50}, 
      PlotStyle -> RGBColor[1, 0, 0]]

LogLinearPlot[Evaluate[y[x] /. solution2], {x, 50, 10000}, 
      PlotRange -> {20, 30}]

Therefore, using the notation in the problem, Y(∞) = ((H(
                m))/(s(m) (Subscript[σ, 0]))) y(∞), and y(∞) is about when the abundance starts to deviate significantly from the equilibrium value.  Because this behavior is mainly due to the exponential dropoff of (Subscript[Y, eq]), it is expected to be rather insensitive to the chioce of m and (Subscript[σ, 0]).  Indeed, by varying m, we find

Table[{(10^t), Clear[yint1, yint2, solution1, solution2, solution3]; 
        solution1 = 
          NDSolve[{y'[
                    x] == -(1/(x^2)) (y[x]^2 - (((0.192 (Subscript[M, Pl]) m σ (x^(3/2)) (E^(-x))))^2)), 
                y[1] == 0.192 (Subscript[M, Pl]) m σ (1^(3/2)) 
                    (E^(-1))} /. {m -> (10^t), 
                g -> 100, σ -> (10^(-10)), 
                (Subscript[M, Pl]) -> 2.44 (10^18)}, y, {x, 1, 20}]; 
        yint1 = Evaluate[y[20] /. solution1[[1]]]; 
        solution2 = 
          NDSolve[{y'[
                    x] == -(1/(x^2)) (y[x]^2 - (((0.192 (Subscript[M, Pl]) m σ (x^(3/2)) (E^(-x))))^2)), 
                y[20] == yint1} /. {m -> (10^t), 
                g -> 100, σ -> (10^(-10)), 
                (Subscript[M, Pl]) -> 2.44 (10^18)}, y, {x, 20, 1000}]; 
        yint2 = Evaluate[y[1000] /. solution2[[1]]]; 
        solution3 = 
          NDSolve[{y'[
                    x] == -(1/(x^2)) (y[x]^2 - (((0.192 (Subscript[M, Pl]) m σ (x^(3/2)) (E^(-x))))^2)), 
                y[1000] == yint2} /. {m -> (10^t), 
                g -> 100, σ -> (10^(-10)), 
                (Subscript[M, Pl]) -> 2.44 (10^18)}, y, {x, 1000, 10000}]; 
        Evaluate[y[10000] /. solution3[[1]]]}, {t, 0.5, 3.5, 
        0.1}]

Only a mild variation for a very wide range of m.

(caution: the version of Mathematica I have seems to have a strange memory effect that always gives the first entry wrong.  I get the reasonable results by quitting the kernel once and re-run the command.)

P-wave

Going back to
\(TraditionalForm\`\(d\ Y\)\/\(d\ x\) = \(-\(1\/x\^2\)\) \(\(s(
                m)\)\/\(H(
                m)\)\) ⟨σ\ v⟩ \((Y\^2 - Y\_eq\^2)\)\)
\(TraditionalForm\`\(\(=\)\(\(-0.368\) \(\(g\_*\)\^\(1/
                2\)\) \(1\/x\^2\) 
          m\ \(M\_Pl\) ⟨σ\ v⟩ \((Y\^2 - Y\_eq\^2)\)\)\)\)
we now substitute \(TraditionalForm\`⟨σ\ v⟩ = \(σ\_0\) x\^\(-1\)\).
((d Y)/(d x)) = -(1/(x^3)) ((s(
                m))/(H(m))) (Subscript[σ, 0])((Y^2) - ((Subscript[Y, eq])^2))
Mathematica again seems to have trouble dealing with big numbers such as (Subscript[M, Pl]).  We can help it by solving for y = ((s(m))/(H(m))) (Subscript[σ, 0]) Y
((d y)/(d x)) = -(1/(x^3)) ((y^2) - 
              ((Subscript[y, eq])^2)),
\(TraditionalForm\`y\_eq = \(\(\(s(m)\)\/\(H(
                  m)\)\) ⟨σ\ v⟩ \(1\/s\) \(\(e\^\(\(-β\)\ m\)\)(\(m\ T\)\/\(2  π\))\)\^\(3/2\) = \(\((\(1\/3\) π\^2\/30)\)\^\(\(-1\)/
                  2\)\) \(M\_Pl\) \(m\/T\^3\) ⟨σ\ v⟩ \(\(e\^\(\(-β\)\ m\)\)(\(m\ T\)\/\(2  π\))\)\^\(3/2\)\)\)
=0.192 (Subscript[M, Pl]) 
          m ⟨σ v⟩ (x^(3/2)) 
          (e^(-x))

N[((((1/3) ((π^2)/30)))^(-1/
              2)) ((((1/(2  π))))^(3/2))]

Mathematica still balks when I let it integrate all the way from x = 1 to 10000.  I break it up into 1 < x < 50 and the rest.

solution1 = 
      NDSolve[{y'[
                x] == -(1/(x^3)) (y[x]^2 - (((0.192 (Subscript[M, Pl]) m σ (x^(3/2)) (E^(-x))))^2)), 
            y[1] == 0.192 (Subscript[M, Pl]) m σ (1^(3/2)) 
                (E^(-1))} /. {m -> 1000, 
            g -> 100, σ -> (10^(-10)), 
            (Subscript[M, Pl]) -> 2.44 (10^18)}, y, {x, 1, 50}]

Evaluate[y[50] /. solution1]

solution2 = 
      NDSolve[{y'[
                x] == -(1/(x^3)) (y[x]^2 - (((0.192 (Subscript[M, Pl]) m σ (x^(3/2)) (E^(-x))))^2)), 
            y[50] == 941.0060810113391`} /. {m -> 1000, 
            g -> 100, σ -> (10^(-10)), 
            (Subscript[M, Pl]) -> 2.44 (10^18)}, y, {x, 50, 10000}]

LogLogPlot[Evaluate[y[x] /. solution1], {x, 1, 50}, 
      PlotRange -> {{1, 10000}, {1, (10^11)}}]

LogLogPlot[Evaluate[y[x] /. solution2], {x, 50, 10000}, 
      PlotRange -> {{1, 10000}, {1, (10^11)}}]

Show[%, %%]

Compared to the equlibrium values,

LogLogPlot[
      0.192 (Subscript[M, Pl]) m σ (x^(3/2)) 
          (E^(-x)) /. {m -> 1000, 
          g -> 100, σ -> (10^(-10)), 
          (Subscript[M, Pl]) -> 2.44 (10^18)}, {x, 1, 1000}, 
      PlotRange -> {{1, 10000}, {1, (10^11)}}, 
      PlotStyle -> RGBColor[1, 0, 0]]

Show[%, %%]

the solution is basically right on the equilibrium until x~20, and then becomes approximately constant afterwards.  This is why we don't expect the result to be sensitive to the initial condition as long as it starts on the equilibrium for x < 10 or so.  Verify this point by taking the ratio

LogLinearPlot[
      0.192 (Subscript[M, Pl]) m σ (x^(3/2)) 
          (E^(-x))/Evaluate[y[x] /. solution1[[1]]] /. {m -> 1000, 
          g -> 100, σ -> (10^(-10)), 
          (Subscript[M, Pl]) -> 2.44 (10^18)}, {x, 1, 50}, 
      PlotStyle -> RGBColor[1, 0, 0]]

LogLinearPlot[Evaluate[y[x] /. solution2], {x, 50, 10000}, 
      PlotRange -> {700, 900}]

Therefore, using the notation in the problem, Y(∞) = ((H(
                m))/(s(m) (Subscript[σ, 0]))) y(∞), and y(∞) is about 2 ((Subscript[x, f])^2) where (Subscript[x, f]) is when the abundance starts to deviate significantly from the equilibrium value.  Because this behavior is mainly due to the exponential dropoff of (Subscript[Y, eq]), it is expected to be rather insensitive to the chioce of m and (Subscript[σ, 0]).  Indeed, by varying m, we find

Table[{(10^t), Clear[yint1, yint2, solution1, solution2, solution3]; 
        solution1 = 
          NDSolve[{y'[
                    x] == -(1/(x^3)) (y[x]^2 - (((0.192 (Subscript[M, Pl]) m σ (x^(3/2)) (E^(-x))))^2)), 
                y[1] == 0.192 (Subscript[M, Pl]) m σ (1^(3/2)) 
                    (E^(-1))} /. {m -> (10^t), 
                g -> 100, σ -> (10^(-10)), 
                (Subscript[M, Pl]) -> 2.44 (10^18)}, y, {x, 1, 20}]; 
        yint1 = Evaluate[y[20] /. solution1[[1]]]; 
        solution2 = 
          NDSolve[{y'[
                    x] == -(1/(x^3)) (y[x]^2 - (((0.192 (Subscript[M, Pl]) m σ (x^(3/2)) (E^(-x))))^2)), 
                y[20] == yint1} /. {m -> (10^t), 
                g -> 100, σ -> (10^(-10)), 
                (Subscript[M, Pl]) -> 2.44 (10^18)}, y, {x, 20, 1000}]; 
        yint2 = Evaluate[y[1000] /. solution2[[1]]]; 
        solution3 = 
          NDSolve[{y'[
                    x] == -(1/(x^3)) (y[x]^2 - (((0.192 (Subscript[M, Pl]) m σ (x^(3/2)) (E^(-x))))^2)), 
                y[1000] == yint2} /. {m -> (10^t), 
                g -> 100, σ -> (10^(-10)), 
                (Subscript[M, Pl]) -> 2.44 (10^18)}, y, {x, 1000, 10000}]; 
        Evaluate[y[10000] /. solution3[[1]]]}, {t, 0.5, 3.5, 
        0.1}]

Only a mild variation for a very wide range of m.

(caution: the version of Mathematica I have seems to have a strange memory effect that always gives the first entry wrong.  I get the reasonable results by quitting the kernel once and re-run the command.)

2. (Subscript[Ω, DM])

S-wave

Once we have the yield, it is easy to convert it to the current energy density.  
(Subscript[ρ, DM]) = m Y(∞) (Subscript[s, 0]) = 
          m ((H(m))/(s(m) (Subscript[σ, 0]))) y(∞) 
            (Subscript[s, 0])
while (Subscript[ρ, c]) = ((3 ((Subscript[H, 0])^2))/(8  π (Subscript[G, N]))) = 
          3  ((Subscript[H, 0])^2) ((Subscript[M, Pl])^2), (Subscript[x, f]) = y(∞), and 
\(TraditionalForm\`\(H(m)\)\^2 = \(1\/\(3 
                 M\_Pl\^2\)\) \(g\_*\) \(π\^2\/30\) m\^4\), \(TraditionalForm\`s(m) = \(g\_*\) \(\(2  π\^2\)\/45\) m\^3\),
and hence
\(TraditionalForm\`Ω\_DM = 
        m\ \ \(\@\(\(\(g\_*\)\/\(3  M\_Pl\^2\)\) \(π\^2\/30\) 
                m\^4\)\) \(45\/\(2 \( π\^2\) \(g\_*\) 
                m\^3\)\) \(x\_f\/σ\_0\) \(s\_0\) 
          1\/\(3 \( H\_0\^2\) M\_Pl\^2\)\)
\(TraditionalForm\`\(\(=\)\(0.252\ \(\(x\_f\) s\_0\)\/\(\(\(g\_*\)\^\(1/2\)\) \(M\_Pl\^3\) \(H\_0\^2\) σ\_0\)\)\)\)

N[(1/3) ((((1/3) ((π^2)/30)))^(1/
              2)) (((((2  (π^2))/45)))^(-1))]

We found (Subscript[x, f]) ≈ 23 for the choice of parameters.  

From HW #3, we found 
s = (Subscript[s, γ]) + 
            (Subscript[s, ν]) = (1 + (21/22)) 
              (Subscript[s, γ]) = (43/22) 2 (( 2  (π^2))/45) ((Subscript[T, 0])^3) = 
              2890 (cm^(-3))

(1 + (21/22)) 2 ((2  (π^2))/45) (hbarc^(-3)) 
        ((Subscript[T, 0])^3) /. {hbarc -> 0.1973 (10^(-4)), 
        (Subscript[T, 0]) -> 2.725*8.617 (10^(-5))}

It is useful to re-express the Hubble constant in the GeV unit (sounds crazy): 
\(TraditionalForm\`H\_0 = \(100\ h\ km\ \(sec\^\(-1\)\) 
            Mpc\^\(-1\) = \(\(\(100\ h\ km\ sec\^\(-1\)\)\/\(3.00×\(10\^5\) km\ sec\^\(-1\)\)\) \(0.1973\ GeV\ fm\)\/\(\(10\^6\) 3.086\ \(10\^16\) m\) = 2.131\ \(10\^\(-42\)\) GeV\ h\)\)\)

(100/(3 (10^5))) ((0.1973 (10^(-15)))/((10^6) 3.086 (10^16)))

Therefore (Subscript[Ω, DM]) (h^2) is given by

0.251646 (((Subscript[x, f]) (Subscript[s, 0]) (hbarc^3))/((g^(1/2)) ((Subscript[M, Pl])^3) ((Subscript[H, 0])^2) (Subscript[σ, 0]))) /. {(Subscript[x, f]) -> 23, (Subscript[s, 0]) -> 2890, 
        hbarc -> 0.1973 (10^(-13)), (Subscript[M, Pl]) -> 2.44 (10^18), 
        (Subscript[H, 0]) -> 2.13 (10^(-42)), 
        g -> 100, (Subscript[σ, 0]) -> (10^(-10))}

This is too big.  

To get the realistic value (Subscript[Ω, DM]) (h^2) = 0.12, we want (Subscript[σ, 0]) = 
        1.6×(10^(-9)) (GeV^(-2)).  This is the range of cross sections of our interest.

1.9492533091097788`/0.12

P-wave

This time, y(∞) is enhanced to 775.  Therefore (Subscript[Ω, DM]) (h^2) is given by

0.251646 ((y (Subscript[s, 0]) (hbarc^3))/((g^(1/2)) ((Subscript[M, Pl])^3) ((Subscript[H, 0])^2) (Subscript[σ, 0]))) /. {y -> 775, (Subscript[s, 0]) -> 2890, 
        hbarc -> 0.1973 (10^(-13)), (Subscript[M, Pl]) -> 2.44 (10^18), 
        (Subscript[H, 0]) -> 2.13 (10^(-42)), 
        g -> 100, (Subscript[σ, 0]) -> (10^(-10))}

This is too big. 

To get the realistic value (Subscript[Ω, DM]) (h^2) = 0.12, we want (Subscript[σ, 0]) = 
        5.5×(10^(-8)) (GeV^(-2)) would do much better.  This is the range of cross sections of our interest.

65.6813615026121`/0.12

In either case, the cross section is of the order of electroweak scale, (Subscript[σ, 0])~((π (α^2))/(m^2))with m~TeV.

Optional

To be followed.