Derivation of Saha equation through "volume accessible to electrons"

Question

I'm self-learning some introductory plasma physics and I'm analyzing the derivation of the Saha equation for a simple first ionization plasma, with the aim of providing some insights on the involved orders of magnitude. The development I see from various sources, as far as I can understand, should go along the following steps:

• In the gas we have $N_b$ neutral atoms, $N_a$ ionized atoms and $N_e$ electrons
• We can write the Boltzmann ratio between the population $B$ (atoms) and $A+e^-$ (ions and free electrons) as $$\frac{N_a}{N_b}=\frac{g_a}{g_b}\exp\left(-\frac{E_a-E_b}{kT}\right)$$ and, neglecting the kinetic energy of the atoms and ions, $E_a=E_b+I+\frac{p_e^2}{2m}$
• We are dealing with a continuum of electron momentum states and we can more precisely write, thinking $\frac{dN_a}{N_b}=\frac{g_a(p)}{g_b}\exp{\left(-\frac{I}{kT}\right)}\exp{\left(-\frac{p_e^2}{2mkT}\right)}\, dp$ , $$\frac{N_a}{N_b}=q\frac{4\pi V}{h^3}\int{g_a(p)\exp{\left(-\frac{p_e^2}{2mkT}\right)}\, dp}$$ where we have incorporated the ionization contribution $I$ and every discrete degeneracy of $g_a$ and $g_b$ in the coefficient $q$
• Now the step where the question arises: all the sources say that $V=\frac{1}{n_e}$ (reciprocal of electron density) so that, after solving the usual integral and using the densities $n_a$ and $n_b$, we obtain $$\frac{n_a n_e}{n_i}=q\frac{\left(2\pi mkT\right)^\frac{3}{2}}{h^3}$$

So, what tells us that $V=\frac{1}{n_e}$? I assume that the measure on phase space is entirely included in the $h^3$ factor, so that I'm not at all convinced by this treatment and precisely in quantifying the space volume as described. Please note that such a way of working it out is presented here and here. On the other side I found more convincing the methods recalling chemical potential or other thermodynamic concepts as in this lecture. So, is this method rigorously correct? If yes, please explain with more detail, because it's completely unclear to me. Also, it may be good to know a book source where that method is used.

---

Addition on UTC March, 30 2023: I've found two books which detail somewhat the derivation

1. Hubeny & Mihalas, "Theory of Stellar Atmospheres", 2015, p. 93
2. Smirnov, "Principles of Statistical Physics", 2006, p. 75

Please note that the two books disagree in some way on the justification for using the relation $V=\frac{1}{n_e}$, one imposing a single electron in the volume, the other invoking the exchange symmetry. This enforces my idea that, albeit in widespread use, the continuum state counting procedure isn't commonly understood.

Answer 1

After thinking a lot, it came to me the idea that it's all about a statistical effect of electron distribution. Let's take the number of energy levels $g_s$ around an energy value and evaluate in how many ways $N_e$ electrons can be arranged among these levels. It is obviously the first step for deriving the Fermi-Dirac statistics, even that's not the point here: $$w_e(N_e)=\frac{g_s!}{N_e!(g_s-N_e)!}$$ Now, in order to evaluate the ratio of the statistical weights, one has to evaluate $$\frac{w_e(N_e+1)}{w_e(N_e)}=\frac{g_s-N_e}{N_e+1}\approx\frac{g_s}{N_e}$$ where the last approximation works well in the classical (non-degenerate) case and with $N_e\gg1$. Now, with all the caution about the meaning of infinitesimal quantities in statistical physics, let's consider the energy shell sufficiently small so that contains $dg_s$ levels and $dN_e$ electrons. Neglecting spin, we'll rewrite the expression in the original question as below (now $g_a$ represents only the ion discrete contribution with $q$ absorbing that discrete weight and the ionization part). $$\frac{N_a}{N_b}=\frac{w_e(N_e+1)}{w_e(N_e)}\frac{g_a}{g_b}\exp\left(-\frac{E_a-E_b}{kT}\right)=q\frac{dg_s}{dN_e}\exp{\left(-\frac{p^2}{2mkT}\right)}$$ Finally, we write explicitly $dg_s(p)=\frac{V 4\pi p^2\, dp}{h^3}$, work with the densities $dn_e=\frac{dN_e}{V}$, $\frac{n_a}{n_b}=\frac{N_a}{N_b}$ and separate the differentials: $$\frac{n_a}{n_b}\,dn_e=q\frac{4\pi}{h^3}p^2\exp{\left(-\frac{p^2}{2mkT}\right)}\,dp$$ Now it's only matter of integrating.