Equation of state: Vardya - Mihalas - Wittmann

These are my notes on equation of state derived initially by Vardya (1965), described by Mihalas (1967) and popularized by Wittmann (1974). It is widely used (or at least present as an option) in spectral synthesis codes like SIR or NICOLE. It is also prepared for the MANCHA code with several modifications and additionally computed quantities. However, the equation is derived in this particular form to be solvable on computing resources half a century ago. From today's perspective, this formulation is rather obsolete. While the results of the VMW EOS are still largely reliable, various tricks introduced to control numerical stability limit its usability to rather restricted range of the pressure and temperature.

Here I derive a simple equation of state for the solar atmosphere following the classical work of Vardya (1965), Mihalas (1967) and Wittmann (1974). The EOS is based on the Saha ionisation equilibrium and the instanteneous chemical equilibrium for the molecules. The main ingredient is hydrogen. It's included as atomic hydrogen (H), negative hydrogen ion (H-), positive hydrogen ion (H+), and as H2 and H2+ molecules. For all other atoms the neutral and the first two ionisation stages are included.

The equations were first published by Vardya (1967). Mihalas (1967) gave a simple numerical algorithm for an efficient solution of the system. Wittmann (1974) copied the equations and add a corrective factor that provides numerical stability at high temperature.

The derivation here follows Mihalas. However, in the original derivation there is a couple of inconsistencies that obscure the procedure. Here I write the equations in a correct and consistent way.

Definitions

Let's first define the pressures:
$p_{\mathrm{H}}$ - partial pressure of the neutral H atoms;
$p_{\mathrm{H^+}}$ - partial pressure of the positive H ions (protons);

$p_{\mathrm{H^-}}$ - partial pressure of the negative H ions;

$p_{\mathrm{H_2}}$ - partial pressure of the H2 molecules;

$p_{\mathrm{H_2^+}}$ - partial pressure of the H2+ molecules;

$p_{\mathrm{e}}$ - electron pressure.

These pressure are proportional to the corresponding number densities:

$$p_{\mathrm{H}} = n_{\mathrm{H}} k T;$$

$$p_{\mathrm{H^+}} = n_{\mathrm{H^+}} k T;$$

$$p_{\mathrm{H^-}} = n_{\mathrm{H^-}} k T;$$

$$p_{\mathrm{H_2}} = n_{\mathrm{H_2}} k T;$$

$$p_{\mathrm{H_2^+}} = n_{\mathrm{H_2^+}} k T;$$

$$p_{\mathrm{e}} = n_{\mathrm{e}} k T.$$

Particle conservation

The total pressure due to the presence of H in the plasma is then:

\begin{equation} p_{\mathrm{H^{tot}}} = p_{\mathrm{H}} + p_{\mathrm{H^+}} + p_{\mathrm{H^-}}+p_{\mathrm{H_2}}+p_{\mathrm{H_2^+}}. \end{equation}

However, if we write this equation in terms of the number density, it becomes: 

\begin{equation} [n_{\mathrm{H^{tot}}}]_\mathrm{particles} = n_{\mathrm{H}} + n_{\mathrm{H^+}} + n_{\mathrm{H^-}}+n_{\mathrm{H_2}}+n_{\mathrm{H_2^+}}. \end{equation}

On the other hand, the particle (nuclei) conservation requires that 

\begin{equation}[n_{\mathrm{H^{tot}}}]_\mathrm{nuclei} = n_{\mathrm{H}} + n_{\mathrm{H^+}} + n_{\mathrm{H^-}}+2 n_{\mathrm{H_2}}+2 n_{\mathrm{H_2^+}}, \end{equation}

and

\begin{equation} \pi_{\mathrm{H^{tot}}} = p_{\mathrm{H}} + p_{\mathrm{H^+}} + p_{\mathrm{H^-}}+2 p_{\mathrm{H_2}}+ 2 p_{\mathrm{H_2^+}}, \end{equation} where $\pi_{\mathrm{H^{tot}}}$ is the total pressure due to the hydrogen particles as if they were all free (atoms or ions), i.e. as if there were no molecules formed. This pressure is clearly fictitious when that conditions is not fulfilled.

The particle conservation in the pressure form can be rewritten as (Note that this is different from the first of Eq.(34) in Mihalas' article.):
\begin{equation}f_1 + f_2 + f_3 + 2 f_4 + 2 f_5 = 1, \end{equation}
or as
\begin{equation}f_1\left(1 + \frac{f_2}{f_1} + \frac{f_3}{f_1} + 2\frac{f_4}{f_1} + 2\frac{f_5}{f_1}\right) = 1, \end{equation}

where

\begin{eqnarray}
f_1&=&\frac{p_{\mathrm{H}}}{\pi_{\mathrm{H^{tot}}}},\nonumber\\
f_2&=&\frac{p_{\mathrm{H^+}}}{\pi_{\mathrm{H^{tot}}}},\nonumber\\
f_3&=&\frac{p_{\mathrm{H^-}}}{\pi_{\mathrm{H^{tot}}}},\nonumber\\
f_4&=&\frac{p_{\mathrm{H_2^+}}}{\pi_{\mathrm{H^{tot}}}},\nonumber\\
f_5&=&\frac{p_{\mathrm{H_2}}}{\pi_{\mathrm{H^{tot}}}}.
\end{eqnarray}
In addition we define:
\begin{equation}
f_{\mathrm{e}}=\frac{p_{\mathrm{e}}}{\pi_{\mathrm{H^{tot}}}}.
\end{equation}

Ionization and dissociation constants

The relations between the partial pressures are set by the Saha ionization formula and by the instantaneous chemical equilibrium. These relations are specified by the corresponding coefficients:
\begin{eqnarray}
K(H)&=&\frac{p_{\mathrm{H^+}}p_{\mathrm{e}}}{p_{\mathrm{H}}},
\nonumber\\
K(H^-)&=&\frac{p_{\mathrm{H}}p_{\mathrm{e}}}{p_{\mathrm{H^-}}},\nonumber\\
K(H_2)&=&\frac{p_{\mathrm{H}^2}}{p_{\mathrm{H_2}}},\nonumber
\nonumber\\
K(H_2^+)&=&\frac{p_{\mathrm{H}} p_{\mathrm{H^+}}}{p_{\mathrm{H_2^+}}}.
\end{eqnarray}
For the sake of simplicity we introduce new variables:
\begin{eqnarray}
g_2 &=& \frac{K(\mathrm{H})}{p_{\mathrm{e}}} ,\nonumber\\
g_3 &=& \frac{p_{\mathrm{e}}}{K(\mathrm{H})} ,\nonumber\\
g_4 &=& \frac{p_{\mathrm{e}}}{K(\mathrm{H_2^+})} ,\nonumber\\
g_5 &=& \frac{p_{\mathrm{e}}}{K(\mathrm{H_2}).}
\end{eqnarray}
The new variables can be easily related to the previously introduced $f_i$ ratios:
\begin{eqnarray}
g_2&=& \frac{K(H)}{p_{\mathrm{e}}}=\frac{p_{\mathrm{H^+}}}{p_{\mathrm{H}}} =\frac{f_2}{f_1}\nonumber\\
g_3&=&\frac{p_{\mathrm{e}}}{K(H^-)} = \frac{p_{\mathrm{H^-}}}{p_{\mathrm{H}}}=\frac{f_3}{f_1},\nonumber\\
g_4 &=& \frac{p_{\mathrm{e}}}{K(\mathrm{H_2^+})} = \frac{p_\mathrm{e}}{p_\mathrm{H}}\;\frac{p_\mathrm{H_2^+}}{p_\mathrm{H^+}} = \frac{f_\mathrm{e}}{f_1}\;\frac{f_4}{f_2}= \frac{f_\mathrm{e}}{f_1}\;\frac{f_4}{f_1}\;\frac{f_1}{f_2} = \frac{f_\mathrm{e}}{f_1}\;\frac{f_4}{f_1}\;\frac{1}{g_2},\nonumber\\
&\Rightarrow&  \frac{f_4}{f_1} = \frac{f_1}{f_\mathrm{e}}\,g_2\,g_4,\nonumber\\
g_5&=& \frac{p_{\mathrm{e}}}{K(\mathrm{H_2})} = \frac{p_\mathrm{e}}{p_\mathrm{H}}\;\frac{p_\mathrm{H_2^+}}{p_\mathrm{H}} = \frac{f_\mathrm{e}}{f_1}\;\frac{f_5}{f_1},\nonumber\\
 &\Rightarrow&  \frac{f_5}{f_1} = \frac{f_1}{f_\mathrm{e}}\,g_5.\nonumber\\
\end{eqnarray}


Now we can write the particle conservation as a function of only two variables $f_1$ and $f_{\mathrm{e}}$:
\begin{equation}f_1(1 + g_2 + g_3) + 2 \frac{f_1^2}{f_{\mathrm{e}}} (g_2 g_4 + g_5) = 1, \end{equation}

Charge conservation

Another equation required to close the system we find in the charge conservation law:
\begin{equation}
n_{\mathrm{e}} = n_{\mathrm{H^+}}-n_{\mathrm{H^-}}+n_{\mathrm{H_2^+}} + Q,
\end{equation}
where $Q$ stands for all the electrons freed by all single- and double-ionised elements other than hydrogen (index $i$ runs over the atomic elements other than hydrogen; index $j$ runs over the ionization stages):
\begin{equation}
Q = \sum_{z = 2}^{n_i} \alpha_i \sum_{i = j}^{2}n_{ji}.
\end{equation}

If we rewrite that in terms of pressures:

\begin{equation}
p_{\mathrm{e}} = p_{\mathrm{H^+}}-p_{\mathrm{H^-}}+p_{\mathrm{H_2^+}} + Q,
\end{equation}
or after dividing both sides with $\pi_{\mathrm{H^{tot}}}$
\begin{equation}
f_{\mathrm{e}} = f_2 - f_3 + f_4 + Q.
\end{equation}

Similar to the particle conservation law, the charge conservation  can be rewritten as a function solely of $f_1$ and $f_{\mathrm{e}}$:

\begin{equation} f_{\mathrm{e}} = f_1 (g_2 - g_3) +  \frac{f_1^2}{f_{\mathrm{e}}} g_2 g_4 + Q.\end{equation}

Solution

To solve the system we introduced the following substitutes:
\begin{eqnarray} a&=&(1 + g_2 + g_3),\nonumber\\ b&=&2\left(1 + \frac{g_2 g_4}{g_5}\right),\nonumber\\ c&=&g_5,\nonumber\\ d&=&g_2 - g_3,\nonumber\\ e&=&\frac{g_2 g_4}{g_5}. \end{eqnarray}

To simplify further derivation, the particle conservation we rewrite as 

\begin{equation}a\, f_1 + c\,b\, \frac{f_1^2}{f_{\mathrm{e}}} = 1, \end{equation}
and the charge conservation as
\begin{equation}f_{\mathrm{e}} = d\,f_1 + e\,c\, \frac{f_1^2}{f_{\mathrm{e}}} + Q.\end{equation}

Solve Eqs.(19...) and (20) for $f_{\mathrm{e}}$:
\begin{equation} f_{\mathrm{e}} = \frac{c\,b\,f_1^2}{1-a\, f_1}, \end{equation} \begin{equation} f_{\mathrm{e}}^2 = d\,f_1 f_{\mathrm{e}} + e\,c\,f_1^2 + Q f_{\mathrm{e}}, \end{equation} and substitute $F_{\mathrm{e}}$ from Eq.(22..) into Eq.(23):
\begin{equation} \frac{c^2\,b^2\,f_1^4}{(1-a\, f_1)^2} = \frac{cbd\,f_1^3}{(1-a\, f_1)} + e\,c\,f_1^2 + \frac{cb\,Qf_1^2}{(1-a\, f_1^2)} \end{equation} We multiply the last expression with $(1-a\, f_1)^2/(c\,f_1^2)$: \begin{eqnarray} &\,& cb^2\,f_1^2 = bd\,f_1(1-a\, f_1) + e(1-a\, f_1)^2 + b\,Q(1-a\, f_1),\nonumber\\ &\Rightarrow& cb^2\,f_1^2 = bd\,f_1 -abd\, f_1^2 + e - 2aef_1 + ea^2 f_1^2 + bQ -abQ\,f_1,\nonumber\\ &\Rightarrow& (cb^2+abd-ea^2)\,f_1^2 + (2ae-bd+abQ)\,f_1 - (e + bQ) = 0,\nonumber\\ &\Rightarrow& c_1\,f_1^2 + c_2\,f_1 + c_3 = 0,\end{eqnarray}

where
\begin{eqnarray}
c_1 &=& cb^2+abd-ea^2,\nonumber\\
c_2 &=& 2ae-bd+abQ,\nonumber\\
c_3 &=&- (e + bQ).
\end{eqnarray}
The solution of the quadratic equation gives the value of $f_1$:
\begin{equation}
f_1 = \frac{-c_2 \pm \sqrt{c_2^2 - 4c_1c_3}}{2c_1}.
\end{equation}

Which solution do we take?

Backward substitution

Once we know $F_1$, we can evaluate the other $F$ ratios:
\begin{equation}
f_2 = g_2 f_1,
\end{equation}
\begin{equation}
f_3 = g_3 f_1,
\end{equation}
Equation (6) can be written as:
\begin{equation}
a f_1 + (1 + e) f_5 = 1,
\end{equation}

thus for $f_5$ we have:

\begin{equation}
f_5 = \frac{1 - a f_1}{1 + e},
\end{equation}

and for $f_4$:

\begin{equation}
f_4 = \frac{g_2 g_4}{g_5} f_5 = e\,f_5.
\end{equation}
Finally, it is possible to compute $f_{\mathrm{e}}$:
\begin{equation}
f_{\mathrm{e}} = f_2 - f_3 + \frac{1}{2}f_4 + Q.
\end{equation}

The total pressure of hydrogen-nuclei is then (Eq.(9)):

\begin{equation} \pi_{\mathrm{H^{tot}}} = \frac{p_{\mathrm{e}}}{f_{\mathrm{e}}}, \end{equation}

and the partial pressures are

\begin{eqnarray} p_{\mathrm{H}} &=& f_1\pi_{\mathrm{H^{tot}}},\nonumber\\ p_{\mathrm{H^+}} &=& f_2\pi_{\mathrm{H^{tot}}},\nonumber\\ p_{\mathrm{H^-}} &=& f_3\pi_{\mathrm{H^{tot}}},\nonumber\\ p_{\mathrm{H_2^+}} &=& \frac{1}{2}\,\pi_{\mathrm{H_2^+}}= \frac{1}{2}f_4\pi_{\mathrm{H^{tot}}},\nonumber\\ p_{\mathrm{H_2}} &=& \frac{1}{2}\,\pi_{\mathrm{H_2}}= \frac{1}{2}f_5\pi_{\mathrm{H^{tot}}}. \end{eqnarray}

Total pressure

The total pressure of hydrogen particles is therefore

\begin{eqnarray} p_{\mathrm{H^{tot}}} &=& p_{\mathrm{H}} + p_{\mathrm{H^+}} + p_{\mathrm{H^-}}+p_{\mathrm{H_2}}+p_{\mathrm{H_2^+}} =
\nonumber\\ &=&\left[f_1 + f_2 + f_3 + \frac{1}{2}(f_4 + f_5)\right]\pi_{\mathrm{H^{tot}}}. \end{eqnarray}

Once we know $p_{\mathrm{H^{tot}}}$, the total gas pressure is
\begin{eqnarray}
p_{\mathrm{g}} &=&  p_{\mathrm{e}} +  p_{\mathrm{a}} =  p_{\mathrm{e}} +  p_{\mathrm{H^{tot}}}+  p_{\mathrm{nonH}}=p_{\mathrm{e}}\left[1 + \frac{p_{\mathrm{H^{tot}}} +  p_{\mathrm{nonH}}}{p_{\mathrm{e}}}\right] .
\end{eqnarray}