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(n\zeta - E_i)}}{\sum \exp{\beta(n\zeta - 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(n\zeta - 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)?

The $H^0$ (n=1) level is doubly degenerate, rather than the $H^+$ (n=0) level. Thus $\mathcal{Z}=1 + 2\alpha e^{\beta \Delta} + \alpha^2 e^{\beta(\Delta + \epsilon)}$. It follows from the technique outlined above that $\langle n \rangle = \frac{2 \alpha e^{\beta \Delta} + 2 \alpha^2 e^{\beta(\Delta + \epsilon)}}{1 + 2 \alpha e^{\beta \Delta} + \alpha^2 e^{\beta(\Delta + \epsilon)}}$ and thus $\alpha = e^{-\beta(\Delta + \epsilon)/2}$.
Since $\alpha \equiv e^{\beta \zeta}$ then $\zeta = -(\Delta + \epsilon)/2$ and we can solve for the fraction of ions under the approximation $\langle n \rangle =1$. Indeed the fractions of $H^+$ and $H^-$ are equal:
$P=\frac{1}{2 + 2 e^{\beta(\Delta - \epsilon)/2}} \approx e^{-\beta (\Delta - \epsilon)/2}$ at low temperatures.