2
$\begingroup$

I have a differential equation which is, $$ \frac{d X_e}{dz} = \frac{C_r(z)}{H(z)}\times \frac{1}{(1+z)}\Big[\frac{38\eta T^3}{25\pi^2}\alpha (T)X_e^2 - b(T)(1 - X_e) \Big] $$ Where $$ \alpha(T) = \langle \sigma v \rangle \approx 9.8 \Big( \frac{\alpha_s}{m_e} \Big)^2 \Big( \frac{E_I}{T} \Big)^{1/2} \ln\Big( \frac{E_I}{T} \Big)$$

$$ b(T) = \alpha(T) \Big( \frac{m_e T}{2\pi} \Big)^{3/2} e^{-E_I/T}$$ $$C_r = \frac{\Gamma_{2s} + 3\mathcal{P}\Gamma_{2p}}{\Gamma_{2s} +3\mathcal{P}\Gamma_{2p} + \beta }$$ $$\beta = b(T)e^{3E_I/4T}$$

with $\Gamma_{2s} = 8.227 s^{-1}$ and $3\mathcal{P}\Gamma_{2p} = \frac{675}{2432}\Big(\frac{E_I}{T}\Big)^3 \frac{H(T)}{(1 - X_e)\zeta(3)}$. Also, $$ H(z) = \sqrt{\Omega_m} H_0 (1 + z)^{3/2} \Big( 1 + \frac{1 + z}{1 + z_{eq}} \Big) $$ $\Omega_m = 0.309$, $H_0 = 1.5\times 10^{-33}$ and $z_{eq} = 3395$. Finally $T = 0.2329meV\times (1 + z)$

This is peeble's equation in Recombination. I want to solve this in between $z= 10$ to $z= 1800$ with the initial condition $X_e[1800] = 1$. But I am unable to do show. The is an error (see at the end). The code is:

Subscript[E, I] = 13.6;
Subscript[m, e] = 5.109989461*10^5;
\[Eta] = 6*10^-10;
Subscript[z, eq] = 3395;
Subscript[H, 0 ] = 1.5*10^-33;
prefac\[Alpha] = (9.8/(137*Subscript[m, e])^2);
prefac\[Beta] = (Subscript[m, e]/(2*Pi))^(3/2);
Subscript[\[CapitalOmega], m] = 0.3089;

All constants are here.

T[z_] := 23.29*(1 + z)*10^-5;
\[Alpha][z_] := prefac\[Alpha]*Sqrt[Subscript[E, I]/T[z]]*Log[Subscript[E, I]/T[z]];
b[z_] := \[Alpha][z]* prefac\[Beta]*T[z]^(3/2)*E^(-Subscript[E, I]/T[z]);
H[z_] := Subscript[H, 0 ]*Sqrt[Subscript[\[CapitalOmega], m]]*(1 + z)^(3/2)*(1 + ((1 + z)/(1 + Subscript[z, eq])))
Subscript[\[CapitalGamma], s] = 8.227;
Subscript[\[CapitalGamma], p][z_] := (675/2432)*(Subscript[E, I]/T[z])^3*H[z]/((1 -Subscript[X, e][z])*Zeta[3]);
\[Beta][z_] := b[z]*E^(3 Subscript[E, I]/(4 T[z]));
Subscript[C, r][z_] := (Subscript[\[CapitalGamma], s] + Subscript[\[CapitalGamma], p][z])/(Subscript[\[CapitalGamma],s] + Subscript[\[CapitalGamma], p][z] + \[Beta][z]);
s = NDSolve[{Subscript[X, e]'[x] == (1/(H[x]*(1 + x)))*((38*\[Eta]*(T[x])^3*\[Alpha][x]*(Subscript[X, e][x])^2)/(23*Pi^2) - b[x]*(1 - Subscript[X, e][x])), Subscript[X, e][1800] == 1}, Subscript[X, e], {x, 10, 1800} ]

Can someone please help?. If the code is hard to read, here is an image. enter image description here

$\endgroup$
6
  • $\begingroup$ FYI: your Friedman equation is (most likely) wrong. To quickly see this, evaluate H[0] which doesn't give H0 as expected (ie the present day value of the Hubble parameter). $\endgroup$
    – Hans Olo
    Commented Dec 22, 2023 at 18:18
  • $\begingroup$ Actually I have not derived that .. It is from Baumann's book (eqn-3.169) $\endgroup$ Commented Dec 22, 2023 at 18:35
  • 1
    $\begingroup$ Ok, it's probably an expansion around the matter-radiation equality. Still it might be prudent to use a correct expression and check how the results changes $\endgroup$
    – Hans Olo
    Commented Dec 22, 2023 at 18:54
  • $\begingroup$ Sure thanks.. can you provide some source for the expression?.. $\endgroup$ Commented Dec 23, 2023 at 5:36
  • 1
    $\begingroup$ See the last equation in the section "Density Parameter" here: en.wikipedia.org/wiki/Friedmann_equations You can ignore curvature and use the equality redshift to relate the radiation to matter. I hope this helps. $\endgroup$
    – Hans Olo
    Commented Dec 23, 2023 at 10:55

1 Answer 1

5
$\begingroup$

The error message means MMA is trying to find the derivative at the initial condition, $x=1800$. Take a look at H[1800] or H[x] /. x->1800: $$\text{6.371926864638378$\grave{ }$*${}^{\wedge}$-29} \left(\frac{1801}{1800_{\text{eq}}+1}+1\right)$$

The problem is that MMA doesn't put in $z_{eq}$. The expert MMA programmers advise us to avoid the use of subscripts. For example see Item 3 of these pitfalls. Note that you have defined the $E_I = 13.6$. That looks like Euler's number with the imaginary unit as a subscript to MMA. That definition won't go away with Clear. See this answer for how to unset subscripts.

Your code runs fine without subscripts:

ClearAll["Global`*"]

EI = 13.6;
me = 5.109989461*10^5;
η = 6*10^-10;
zeq = 3395;
H0 = 1.5*10^-33;
prefacα = (9.8/(137*me)^2);
prefacβ = (me/(2*Pi))^(3/2);
Ωm = 0.3089;

T[z_] := 23.29*(1 + z)*10^-5;
α[z_] := prefacα*Sqrt[EI/T[z]]*Log[EI/T[z]];
b[z_] := α[z]*prefacβ*T[z]^(3/2)*E^(-EI/T[z]);
H[z_] := H0*Sqrt[Ωm]*(1 + z)^(3/2)*(1 + ((1 + z)/(1 + zeq)))
Γs = 8.227;
Γp[z_] := (675/2432)*(EI/T[z])^3*H[z]/((1 - Xe[z])*Zeta[3]);
β[z_] := b[z]*E^(3 EI/(4 T[z]));
Cr[z_] := (Γs + Γp[z])/(Γs + Γp[z] + β[z]);

ode = Xe'[x] == (1/(H[x]*(1 + x)))*
       ((38*η*(T[x])^3*α[x]*(Xe[x])^2)/(23*Pi^2) -
       b[x]*(1 - Xe[x]));

ode /. x -> 1800  /. Xe[1800] -> 1

(* Xe'[1800] == 1.67186 *)

And the solution is

s = NDSolveValue[{ode, Xe[1800] == 1}, Xe, {x, 10, 1800}];
Plot[s[x], {x, 10, 1800}, GridLines -> Automatic]

enter image description here

$\endgroup$
3
  • $\begingroup$ Thanks, Can you tell me one more thing? If I want to find the value of z where Xe is 0.5.. How to do that? $\endgroup$ Commented Dec 22, 2023 at 18:30
  • 2
    $\begingroup$ To find $z$ that solves $X_e(z) =.5$, evaluate FindInstance[{s[z] == .5, 10 < z < 1800}, z]. to get (* {{z -> 1376.57}} *). $\endgroup$
    – LouisB
    Commented Dec 22, 2023 at 19:18
  • $\begingroup$ Thank you once again for your help. $\endgroup$ Commented Dec 29, 2023 at 1:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.