(1) the equation of state can be approximated by the ideal gas law:
\begin{equation}p = n\,kT,\end{equation}
where $p$ is the gas pressure, $n$ is the total number density, $T$ is the temperature and $k$ is the Boltzmann constant.
For partially ionized gas:
(2) the ionization fraction can be approximated by its equilibrium value that is set by Saha's ionization formula:
\begin{equation}\frac{n_{i+1}n_{\mathrm{e}}}{n_i} = \frac{1}{\Phi(T)},\end{equation}
$$\frac{1}{\Phi(T)} = \frac{2 U_{i+1}}{U_i} \left(\frac{2 \pi m_e k T}{h^2}\right)^{3/2} \mathrm{e}^{-\chi/kT},$$
where $U_i$ and $U_{i+1}$ are the partition functions of the ionization stages $i$ and $i+1$, $m_\mathrm{e}$ is the electron mass, $h$ is Planck's constant and $\chi$ is the ionization energy.
For a mixture of gases, it means that each component can be approximated by ideal gas, so that the partial pressure of the $i$ component is:
\begin{equation}p_j = n_j\,kT, \end{equation}
and the total pressure is equal to the sum of the partial pressures (Dalton's law):
\begin{equation}p = \sum p_\mathrm{j},\end{equation}
where the summation goes over all components. The components in this case are the free electrons and all the species of the neutral and charged atoms and molecules that are present in the plasma. In the general case the system of the equations is non-linear and contains large number of equations (one for each atomic or molecular specie in every relevant ionization stage). In the theory of stellar atmosphere it is often necessary to solve this system when two parameters are know.
Nevertheless, there are some trivial solutions. The most notorious is the example of the atmosphere composed of pure hydrogen when the H2 molecules and the negative hydrogen ions are neglected. In that case there are only three types of particles present: protons, neutral H-atoms and free electrons. The rest of this blog refers to that case only.
Solution for given temperature and pressure
This simple exercise is described in Mihalas (1970, p.73, 1970stat.book.....M). The total number of particles in this case is:
\begin{equation}n = n_\mathrm{e} + n_\mathrm{H} + n_\mathrm{H^+} = n_\mathrm{e} + n_\mathrm{H}^\mathrm{tot},\label{eq:pureh1}\end{equation}
where $n_\mathrm{H}^\mathrm{tot} = n_\mathrm{H} + n_\mathrm{H^+}$ is the total number of neutral and ionized H atoms (the total number of H nuclei). The ratio of $n_\mathrm{H^+}$ and $n_\mathrm{H}$ is given by Saha's equation:
\begin{equation}\frac{n_\mathrm{H^+} n_{\mathrm{e}}}{n_\mathrm{H}} = \frac{1}{\Phi(T)}. \label{eq:pureh2}\end{equation}
In addition we know that the number of electrons must be equal to the number of protons: \begin{equation}n_\mathrm{H^+} = n_\mathrm{e}.\label{eq:pureh3}\end{equation}
If we eliminate $n_\mathrm{H}$ and $n_\mathrm{H^+}$ from Eq.$\ref{eq:pureh2}$ and Eq.$\ref{eq:pureh3}$, then Eq.$\ref{eq:pureh1}$ becomes quadratic equation for the electron pressure:
\begin{equation}n = 2 n_\mathrm{e} + n_\mathrm{e}^2 \Phi(T),\label{eq:pureh4}\end{equation}
with the solution:
\begin{equation}n_{\mathrm{e}} = \left(\sqrt{1 + n\,\Phi(T)}-1 \right)\,\frac{1}{\Phi(T)}.\label{eq:pureh5}\end{equation}
To express this solution in terms of the total pressure $p$ and the electron pressure $p_\mathrm{e}$, we assume that the ideal gas law is valid not only for the mixture, but for all of its components as well. Thus: $p_\mathrm{e} = n_\mathrm{e}\,kT$ and the solution for the electron pressure becomes:
\begin{equation}p_{\mathrm{e}} = \left(\sqrt{1 + \frac{p}{kT}\,\Phi(T)}-1 \right)\,\frac{kT}{\Phi(T)}.\label{eq:pureh6}\end{equation}
(In Mihalas' book there is a typo in this equation.) The partial pressure of the hydrogen atoms (neutral and ionized) is obviously:
$$p_\mathrm{H}^\mathrm{tot} = p - p_\mathrm{e}.$$Once we know the electron pressure, it is easy to find other variables.
The ionization fraction is defined as:
\begin{equation}x \equiv \frac{n_\mathrm{H^+}}{n_\mathrm{H}^\mathrm{tot}}= \frac{n_\mathrm{e}}{n_\mathrm{H}^\mathrm{tot}}.\label{eq:x}\end{equation}
Density
The mass density of a mixture specified by the number densities of its components is:
\begin{equation}\rho \equiv \sum_i n_i\,m_i = n_\mathrm{H}\,m_\mathrm{H}+n_\mathrm{H^+}\,m_\mathrm{H^+}+n_\mathrm{e}\,m_\mathrm{e},\label{eq:rho}\end{equation}
where $m_\mathrm{H}$, $m_\mathrm{H^+}$ and $m_\mathrm{e}$ are masses of one neutral hydrogen atom, one proton and one electron. Since $n_\mathrm{H^+} = n_\mathrm{e}$ and $m_\mathrm{H} \approx m_\mathrm{H^+} \gg m_\mathrm{e}$, the density is
$$\rho \approx m_\mathrm{H}\,n_\mathrm{H}^\mathrm{tot}.$$
Alternatives
If the density is given instead of the total pressure, then the system becomes:
\begin{eqnarray}\rho &=& n_\mathrm{H}\,m_\mathrm{H} + n_\mathrm{e}(m_\mathrm{H} + m_\mathrm{e}),\nonumber\\ \Phi\,n_\mathrm{e}^2 &=& n_\mathrm{H},\nonumber\end{eqnarray}
yielding again to a quadratic equation for the electron density:
$$\Phi\,m_\mathrm{H}\,n_\mathrm{e}^2 +(m_\mathrm{H} + m_\mathrm{e}) n_\mathrm{e} - \rho = 0.$$
If the electron pressure is given instead of the density of the pressure, then the solution for the total pressure is $$p = p_\mathrm{e}(\Phi\,p_\mathrm{e} + 2kT),$$ what is equivalent to Eq.$\ref{eq:pureh4}$.
Mean molecular weight
The mean molecular weight is defined as
$$\mu = \frac{\rho}{n\,m_\mathscr{A}},$$
where $m_\mathscr{A}$ is the atomic mass unit in grams. If all gas in neutral,
$$\mu = \frac{m_\mathrm{H}\,n_\mathrm{H}}{n_\mathrm{H}\,m_\mathscr{A}} = A_\mathrm{H},$$
with $A_\mathrm{H}$ being the atomic mass of hydrogen ($\approx 1.008$). If all atoms are ionized, the density remains unchanged, but the number of particles is doubled and, therefore, $\mu \approx 0.5$. (The mean molecular weight is always roughly equal to the total number of nuclei per particle.)
Specific internal energy
The internal energy is equal to the sum of the kinetic part and the part due to the ionization (in principle, there should be another term due to the excitation here, but we neglect it in this example):
$$E_\mathrm{int} = \frac{3}{2}N\,kT + N_\mathrm{e}\,\chi.$$
The specific internal energy pre volume is then $$\varepsilon_V = \frac{E_\mathrm{int}}{V} = \frac{3}{2}n\,kT + n_\mathrm{e}\,\chi,$$ and the specific internal energy per mass is simply
$$\varepsilon_m = \frac{\varepsilon_V}{\rho}.$$
Specific heat capacity at constant volume (density)
The heat capacity at constant volume $C_V$ is defined as: $$C_V = \left(\frac{\partial E_\mathrm{int}}{\partial T}\right)_V.$$ The heat capacity per mass at constant volume $c_V$ is defined as: \begin{equation}c_V = \left(\frac{\partial \varepsilon_m}{\partial T}\right)_\rho = \left[\frac{\partial }{\partial T}\left(\frac{3}{2}\frac{p}{\rho} + \frac{n_\mathrm{e}}{\rho}\,\chi\right) \right]_\rho = \frac{3}{2}\frac{1}{\rho}\left(\frac{\partial p }{\partial T}\right)_\rho + \frac{\chi}{\rho} \left(\frac{\partial n_\mathrm{e}}{\partial T}\right)_\rho. \label{eq:cv}\end{equation} From the defition of the ionization fraction (Eq.$\ref{eq:x}$), it directly follows: $$\left(\frac{\partial n_\mathrm{e}}{\partial T}\right)_\rho = \left(\frac{\partial x n_\mathrm{H}^\mathrm{tot}}{\partial T}\right)_\rho = n_\mathrm{H}^\mathrm{tot}\left(\frac{\partial x }{\partial T}\right)_\rho.$$
The gas pressure is now \begin{equation}p = n\,kT = (n_\mathrm{H}^\mathrm{tot} + n_\mathrm{e})\,kT = n_\mathrm{H}^\mathrm{tot}(1 + x)\,kT = \frac{\rho}{m_\mathrm{H}}(1 + x)\,kT,\label{eq:p}\end{equation} and thus \begin{equation}\left(\frac{\partial p }{\partial T}\right)_\rho = n_\mathrm{H}^\mathrm{tot} k \left[(1 + x) + T \left(\frac{\partial x}{\partial T}\right)_\rho \right].\label{eq:dpdt}\end{equation}The heat capacity for the constant volume becomes:\begin{eqnarray}c_V &=& \frac{n_\mathrm{H}^\mathrm{tot}}{\rho}\left[\frac{3}{2} k \left((1 + x) + T \left(\frac{\partial x}{\partial T}\right)_\rho \right) + \chi \left(\frac{\partial x }{\partial T}\right)_\rho\right].\nonumber\\ &=& \frac{n_\mathrm{H}^\mathrm{tot}}{\rho}\left[\frac{3}{2} k (1 + x) + \left(\frac{3}{2} kT + \chi\right) \left( \frac{\partial x }{\partial T}\right)_\rho\right].\label{eq:cv0}\end{eqnarray}
To find $\partial x/\partial T$, let's rewrite Saha's equation in terms of $x$:
\begin{equation}\frac{n_\mathrm{H^+} n_{\mathrm{e}}}{n_\mathrm{H}} = \frac{x^2}{1-x}\,n_\mathrm{H}^\mathrm{tot}\end{equation} \begin{equation}\frac{ x^2}{1-x} = \frac{m_\mathrm{H}}{\rho} \left( \frac{2\pi m_\mathrm{e} \,kT}{h^2}\right)^{3/2} \mathrm{e}^{-\frac{\chi}{kT}} = q_\rho\,T^{3/2} \mathrm{e}^{-\frac{\chi}{kT}},\label{eq:sahax}\end{equation} where $q_\rho$ contains all the constants. Now the derivative of Eq.$\ref{eq:sahax}$ with respect to $T$ and constant $\rho$ is\begin{eqnarray}\frac{2x(1-x)+x^2}{(1-x)^2}\,\left(\frac{\partial x}{\partial T}\right)_\rho &=& q_\rho \left[\frac{3}{2}T^{1/2} + T^{3/2}\left(\frac{\chi}{kT^2}\right)\right]\,\mathrm{e}^{-\frac{\chi}{kT}},\nonumber\\\frac{x(2-x)}{(1-x)^2}\,\left(\frac{\partial x}{\partial T}\right)_\rho &=& \frac{1}{T} \left(q_\rho T^{3/2} \,\mathrm{e}^{-\frac{\chi}{kT}}\right) \left(\frac{3}{2} + \frac{\chi}{kT}\right),\nonumber\\
\left(\frac{\partial x}{\partial T}\right)_\rho &=& \frac{(1-x)^2}{x(2-x)}\,\frac{x^2}{1-x}\,\frac{1}{T} \left(\frac{3}{2} + \frac{\chi}{kT}\right),\nonumber\\ \left(\frac{\partial x}{\partial T}\right)_\rho &=& \frac{x(1-x)}{(2-x)}\, \frac{1}{T}\, \left(\frac{3}{2} + \frac{\chi}{kT}\right).\label{eq:dxdt}\end{eqnarray} At this point it is worth to introduce to two new variables to simplify the notation:\begin{equation}D = D(x) \equiv \frac{x\,(1-x)}{(1+x)\,(2-x)},\end{equation}and\begin{equation}H = H(T) \equiv \frac{3}{2} + \frac{\chi}{kT}.\end{equation}
Note that $D\rightarrow 0$ when $x\rightarrow 0$ or $x\rightarrow 1$. Note as well that $1 - D$ reduces to:
$$1-D = \frac{2}{(1+x)\,(2-x)}.$$
Equation $\ref{eq:dxdt}$ we now write simply as: \begin{equation}\left(\frac{\partial x}{\partial T}\right)_\rho = (1 + x)\,\frac{1}{T}\,D\,H.\label{eq:dxdt1}\end{equation}
Substituting Eq.$\ref{eq:dxdt1}$ into Eq.$\ref{eq:cv0}$ yields to:
\begin{equation}c_V = \frac{n_\mathrm{H}^\mathrm{tot}}{\rho}\left[\frac{3}{2} k (1 + x) + \left(\frac{3}{2} kT + \chi\right)\, (1 + x)\,\frac{1}{T}\,D\,H\right]\end{equation}
\begin{equation} c_V= \frac{k}{m_\mathrm{H}}\,(1 + x)\, \left[\frac{3}{2} + D\,H^2\right],\label{eq:cv7}\end{equation}where we used $n_\mathrm{H}^\mathrm{tot}/\rho = 1/m_\mathrm{H}$. Explicitly in terms of $x$ the specific heat at the constant density is:
\begin{equation} c_V= \frac{k}{m_\mathrm{H}}\,(1 + x)\, \left[\frac{3}{2} + \frac{x\,(1-x)}{(1+x)\,(2-x)}\,\left(\frac{3}{2} + \frac{\chi}{kT}\right)^2\right].\label{eq:cv8}\end{equation}
Specific heat capacity at constant pressure
The next thing to derive is an equivalent expression for the heat capacity per mass at constant pressure. We start from a relation that can easily be derived from the first principle of thermodynamics:
\begin{equation}c_p = c_V + \frac{T}{\rho^2}\,\left(\frac{\partial p}{\partial T}\right)_\rho^2\left(\frac{\partial p}{\partial \rho}\right)_T^{-1}.\end{equation}
First we rewrite it as:
\begin{eqnarray}c_p &=& c_V + \frac{T}{\rho^2}\,\left[\frac{p}{T} \left(\frac{\partial \ln p}{\partial \ln T}\right)_\rho\right]^2\,\left[\frac{p}{\rho}\,\left(\frac{\partial \ln p}{\partial \ln \rho}\right)_T\right]^{-1}.\nonumber\\ &=& c_V + \frac{k}{m_\mathrm{H}}\,(1+x)\,\left(\frac{\partial \ln p}{\partial \ln T}\right)_\rho^2\,\left(\frac{\partial \ln p}{\partial \ln \rho}\right)_T^{-1}.\end{eqnarray}
The two derivatives we find by differentiating the logarithm of Eq.$\ref{eq:p}$:
$$\ln p = \ln (k/m_\mathrm{H}) + \ln\rho + \ln T + \ln (1 + x),$$
$$\mathrm{d} \ln p = \mathrm{d} \ln\rho + \mathrm{d} \ln T + \mathrm{d} \ln (1 + x),$$
yielding to:
\begin{equation}\left(\frac{\partial \ln p}{\partial \ln T}\right)_\rho = 1 + \frac{1}{1+x}\left(\frac{\partial x}{\partial \ln T}\right)_\rho = 1 + D\,H,\label{eq:dlnp_rho}\end{equation}
\begin{equation}\left(\frac{\partial \ln p}{\partial \ln \rho}\right)_T = 1 + \frac{1}{1+x}\left(\frac{\partial x}{\partial \ln \rho}\right)_T = 1 - D.\label{eq:dlnp_t}\end{equation}
Therefore, for $c_p$ we can write:
\begin{equation}
c_p = c_V + \frac{k}{m_\mathrm{H}}\,(1+x)\,\frac{(1 + D\,H)^2}{1-D}.
\end{equation}
This expression can be written explicitly in terms of $x$ as:
\begin{eqnarray}
c_p &=& \frac{k}{m_\mathrm{H}}\,(1+x) \left[\frac{3}{2} + \frac{1 + 2 D\,H + D\,H^2}{1 - D} \right]\nonumber\\
&=& \frac{k}{m_\mathrm{H}}\,(1+x) \left[\frac{3}{2} + \frac{(2-x)(1+x) + 2(1-x)x \,H + (1-x)x\,H^2}{2} \right]\nonumber\\
&=& \frac{k}{m_\mathrm{H}}\,(1+x) \left[\frac{5}{2} + \frac{(1-x)x}{2}( 1 + 2 \,H + x\,H^2) \right]\nonumber\\ &=& \frac{k}{m_\mathrm{H}}\,(1+x) \left[\frac{5}{2} + \frac{(1-x)x}{2}(1 +H)^2 \right]\nonumber\\
&=& \frac{k}{m_\mathrm{H}}\,(1+x) \left[\frac{5}{2} + \frac{(1-x)x}{2}\left(\frac{5}{2} +\frac{\chi}{kT}\right)^2 \right]
\end{eqnarray}
The ratio of the heat capacities (sometimes labeled with $\gamma$) for the partially ionized pure hydrogen is:
\begin{equation}\gamma \equiv \frac{c_p}{c_V} = 1 + \frac{(1+ D\,H)^2}{(1 - D)\,(3/2 + D\,H^2)}
.\label{eq:cpcv}\end{equation}
Mayer's relation
Mayer's relation for ideal gas states that $c_p = c_v + nk$. In a gas mixture, this should apply to each component of free particles that contribute to the internal energy. Before we had: \begin{equation} c_p = c_V + \frac{k}{m_\mathrm{H}}\,(1+x)\,\frac{(1 + D\,H)^2}{1-D}.
\end{equation} Therefore, for Mayer's equation to be true, we have to show that: $$n = \frac{1}{m_\mathrm{H}}\,(1+x)\,\frac{(1 + D\,H)^2}{1-D} $$
Adiabatic exponents
In addition we derive the explicit expression of the adiabatic exponents $\Gamma$ of the partially ionized pure hydrogen gas. The exponents are defined as:
\begin{equation}\gamma \equiv \Gamma_1 \equiv \left(\frac{\partial \ln p}{\partial \ln \rho}\right)_S,\end{equation}
\begin{equation}\nabla_\mathrm{ad}\equiv\frac{\Gamma_2 -1}{\Gamma_2} \equiv \left(\frac{\partial \ln T}{\partial \ln p}\right)_S,\label{eq:G1}\end{equation}
\begin{equation}\gamma\,\nabla_\mathrm{ad} \equiv \Gamma_3 -1 \equiv \left(\frac{\partial \ln T}{\partial \ln \rho}\right)_S.\end{equation}
Comment on the meaining of $\Gamma$'s. Note that both $\gamma$ in Eq.$\ref{eq:cpcv}$ and in Eq.$\ref{eq:G1}$ is not the same. It is easy to express the three adiabatic coefficients as function of the gradients in Eqs.$\ref{eq:dlnp_rho}$ and $\ref{eq:dlnp_t}$:
\begin{equation}\Gamma_1 =\left( \frac{\partial \ln p}{\partial \ln \rho}\right)_T,\label{eq:G1_1}\end{equation}
\begin{equation}\frac{\Gamma_2 -1}{\Gamma_2} = \frac{k}{m_\mathrm{H}}\,\frac{1}{c_p} \left( \frac{\partial \ln p}{\partial \ln T}\right)_\rho \left( \frac{\partial \ln p}{\partial \ln \rho}\right)_T^{-1},\label{eq:G2_1}\end{equation}
\begin{equation}\Gamma_3 - 1 = \frac{k}{m_\mathrm{H}}\,\frac{1}{c_v} \left( \frac{\partial \ln p}{\partial \ln T}\right)_\rho.\label{eq:G3_1}\end{equation}
After substituting the results for the gradients Eqs.$\ref{eq:dlnp_rho}$ and $\ref{eq:dlnp_t}$, the final expression of the adiabatic exponents for the partially ionized pure hydrogen plasma are:
\begin{equation}\Gamma_1= \frac{c_p}{c_v}\,(1 - D),\end{equation}
\begin{equation}\frac{\Gamma_2 -1}{\Gamma_2} = \frac{c_p - c_V}{c_p}\,\frac{1}{1 + D\,H},\end{equation}
\begin{equation}\Gamma_3 - 1= \frac{c_p - c_V}{c_V}\,\frac{1-D}{1 + D\,H}.\end{equation}
Entropy
Let's now compute the specific entropy per mass $s$. It can be shown that (see my notes on computing entropy for $T$ and $\rho$ as independent variables: $$ds = \frac{1}{T} \left[\left(\frac{\partial \varepsilon_m}{\partial \rho}\right)_T - \frac{p}{\rho^2}\right] d\rho + \frac{1}{T}\left(\frac{\partial \varepsilon_m}{\partial T}\right)_\rho dT $$ To solve for these derivatives analytically, we need to find $\varepsilon_m = \varepsilon_m(T, \rho)$: $$\varepsilon_m = \frac{1}{\rho} \left[\frac{3}{2}n\,kT + n_\mathrm{e}\,\chi\right]$$ and its derivatives: $$\left(\frac{\partial \varepsilon_m}{\partial \rho}\right)_T = -\frac{1}{\rho^2} \left[\frac{3}{2}n\,kT + n_\mathrm{e}\,\chi\right] + \frac{3}{2} \frac{kT}{\rho} \left(\frac{\partial n}{\partial \rho}\right)_T + \frac{\chi}{\rho}\left(\frac{\partial n_\mathrm{e}}{\partial \rho}\right)_T$$ and $$\left(\frac{\partial \varepsilon_m}{\partial T}\right)_\rho =$$ Where we still need to express $n$ and $n_e$ as functions of $T$, $\rho$. By definition: $$\rho = n_\mathrm{H}\,m_\mathrm{H} + n_\mathrm{e}(m_\mathrm{H} + m_\mathrm{e})$$ and $$n = n_\mathrm{H} + 2 n_\mathrm{e}$$ while Saha gives us: $$\Phi\,n_\mathrm{e}^2 = n_\mathrm{H}$$
As we saw before, the last two equations combine into a quadratic equation with a solution: $$n_{\mathrm{e}} = \left(\sqrt{1 + n\,\Phi(T)}-1 \right)\,\frac{1}{\Phi(T)}.$$ And for $\rho(n)$ we get explicitly: $$\rho = (n - n_\mathrm{e}) \,m_\mathrm{H} + n_\mathrm{e} m_\mathrm{e}$$ $$\rho = \left(n - \left(\sqrt{1 + n\,\Phi(T)}-1 \right)\,\frac{1}{\Phi(T)}\right) \,m_\mathrm{H} + \left(\sqrt{1 + n\,\Phi(T)}-1 \right)\,\frac{1}{\Phi(T)} m_\mathrm{e}$$ $$\rho = n m_\mathrm{H} - \frac{\sqrt{1 + n\,\Phi(T)}-1 }{\Phi(T)} (m_\mathrm{H} -m_\mathrm{e})$$ Therefore: $$\left(\frac{\partial \rho }{\partial n}\right)_T = m_\mathrm{H} - (m_\mathrm{H} -m_\mathrm{e}) \frac{1}{2 \sqrt{1 + n\,\Phi(T)}}$$ $$\left(\frac{\partial n }{\partial \rho}\right)_T = \frac{2 \sqrt{1 + n\,\Phi(T)}}{2 m_\mathrm{H}\, \sqrt{1 + n\,\Phi(T)} - (m_\mathrm{H} -m_\mathrm{e})}$$
References
The general thermodynamical expressions and their derivation can be found in any textbook on the thermodynamics or on the stellar structure, e.g. in Chandrasekhar, S. (1967, 1967aits.book.....C). There is a brief mentioning of the equation of state for the pure hydrogen in Mihalas, D. (1970, p.73, 1970stat.book.....M). Biermann, L. (1942, 1942ZA.....21..320B) was the first to derive the explicit expressions for the $c_p$ and $c_V$ (see Lobel, A. (2001, 2001ApJ...558..780L)). Knapp, G (lecture notes) gives most of the derivations above neatly and tidily. Hansen, C.J. et al (2004, 2004sipp.book.....H) derive several cases including the pure hydrogen and discuss the solution.
IDL implementation
IDL function that evaluates the equation of state for pure hydrogen for given temperature and gas pressure can be downloaded from here: eos_pureh.pro.
