Solution (source code)

= Solution

Ignoring <helium>, <charge neutrality> gives $n_e=n_p$, while <baryon> conservation gives $n_B=n_p+n_H$. Defining the <hydrogen ionization fraction> $X_e=n_e/n_B$ yields $n_p=n_e=X_en_B$ and $n_H=(1-X_e)n_B$. Substituting in the abundance relation gives
$$
\frac{1-X_e}{X_e^2}=n_B\left(\frac{2\pi}{m_eT}\right)^{3/2}e^{I/T}.
$$
Use $n_B=\eta n_\gamma$ and $n_\gamma=2\zeta(3)T^3/\pi^2$ to obtain the <Saha equation>:
$$
\boxed{\frac{1-X_e}{X_e^2}=\frac{2\zeta(3)}{\pi^2}\eta\left(\frac{2\pi T}{m_e}\right)^{3/2}e^{I/T}.}
$$
If the positive right side is denoted by $A(T)$, the physical root is
$$
\boxed{X_e(T)=\frac2{1+\sqrt{1+4A(T)}}.}
$$
This form avoids subtracting nearly equal numbers when $A$ is small, and makes the limits transparent: $X_e\simeq1-A$ for $A\ll1$, whereas $X_e\simeq A^{-1/2}$ for $A\gg1$. It is an equilibrium prediction; a kinetic <cosmological recombination> treatment is needed after reactions cease to keep up with expansion.