What is the fraction of ionized hydrogen atoms under certain approximation

Given the ionization reaction $2H^0 <=> H^+ + H^-$ where:

$H^+$ is the ionized state with electron occupancy n = 0 and energy 0

$H^0$ is the doubly degenerate neutral hydrogen atom with electron occupancy n = 1 and energy $-Delta$.

$H^-$ is the ionized state with electron occupancy n = 2 and energy $-(epsilon + delta)$

I want to calculate the fractions of $H^+$ and $H^-$ ions using the condition $ langle n rangle = 1$ (which “is a valid approximation in the case of the solar photosphere since the free electron concentration is less than the hydrogen concentration by a factor greater than $10^3$”). Further, I want to show that at low temperature, this fraction is $exp{-beta(Delta-epsilon)/2}$.

In summary, my approach has been to

1. Calculate the grand partition function $mathcal{Z}$

2. Use the condition $langle n rangle = frac{alpha}{mathcal{Z}} frac{partial mathcal{Z}}{partial alpha}$ ($alpha equiv e^{beta zeta}$) to determine an approximate value for the chemical potential $zeta$.

3. Evaluate the fraction of $H^+$ and $H^-$ using $P(n,i) = frac{exp{beta(nzeta – E_i)}}{sum exp{beta(nzeta – E_i)}}$

The problem I am having is the third part. I expect the fraction of $H^+$ and the fraction of $H^-$ to be the same, which is not the case according to my calculation. Also, my solution does not show that at low T the fraction is $exp{-beta(Delta-epsilon)/2}$.

In greater detail, I have: $$mathcal{Z} = sum exp{beta(nzeta – E_i)} = 2 + alpha e^{beta Delta} + alpha^2 e^{beta (Delta + epsilon)}$$

Then $langle n rangle = frac{alpha}{mathcal{Z}} frac{partial mathcal{Z}}{partial alpha}= frac{alpha e^{beta Delta} + 2 alpha^2 e^{beta epsilon}}{2 + alpha e^{beta Delta} + alpha^2 e^{beta epsilon}}=1 $ implies $alpha = pm (2e^{beta epsilon})^{1/2}$ where I discard the negative solution.

Then since I assume the gas is ideal $zeta = k_B T ln{alpha} = frac{1}{2} k_B T ln{2} – frac{epsilon}{2}$

Using this value for the chemical potential, I plug in to the grand canonical distribution:

For the fraction of $H^+$ I expect:

$$P(n=0) = frac{1}{2 + alpha e^{beta Delta} + alpha^2 e^{beta epsilon}}$$

and for the fraction of $H^-$ I expect:

$$P(n=2) = frac{e^{beta (2 zeta + epsilon)}}{2 + alpha e^{beta Delta} + alpha^2 e^{beta epsilon}}$$

These two probabilities are not the same. However, I don’t understand how that could be, since I imagine the dissociation equation necessarily implies equal fractions of $H^+$ and $H^-$. Further, since neither of the probabilities agrees with the low temperature limit I am trying to show, I imagine that I am misunderstanding how to calculate the fraction of each ionization state.

My questions:

• Does the dissociation equation imply equal fractions of $H^+$ and $H^-$?

• Is the fraction of $H^+$/$H^-$ proportional to the Gibbs factor divided by the grand partition function (i.e. $P(n=0)$/$P(n=2)$ as I’ve written it)?