Opacity for realistic 3D MHD simulations of cool stellar atmospheres

The first paper of Andrea Perdomo Garcia is just submitted for publication in Astronomy & Astrophysics, and out on arxiv.org/abs/2306.03744. The paper is all about computing the opacities for realistic modelling of cool stellar atmospheres. It is divided in three unities. First (Section 3) it describes the computation of detailed monochromatic opacity including millions of atomic and molecular spectral lines and millions of wavelength points. For this the code SYNSPEC (Hubeny and Lanz, 2011, 2017a, b) is used. Then (Section 4) the monochromatic opacities are used to construct opacity distribution function which reduces the number of wavelength points from millions to thousands. The results are compared in detail with ones produced by Kurucz. Some striking similarities and some warning differences are found. Finally (Section 5), the opacity distribution function to construct opacity bins. This method, originally proposed by Nordlund (1982) is the key ingredient for realistically simulating stellar atmospheres in 3D as it reduced the problem further, from thousands of wavelength points to only a few. However, the method depends on a choice of some free parameters. In our paper the possible choices are carefully analyzed and some interesting conclusions are offered. 

In Sect.3 there are two figures (Figs.2 and 3) that I find very useful and illustrative. The monochromatic opacity (Fig.2) and the radiative heating rate (Fig.3) are shown as 2D functions of wavelength (X-axis) and height in the atmosphere (Y-axis) for four different cool stars (all with solar metalicity). Optical depths in the continuum and continuum+lines are overplotted.

(Andrea is the final year PhD student at Instituto de Astrofisica de Canarias and Univeridad de La Laguna, supervised by Manolo Collados Vera and myself. Stay tuned, more cool stuff is coming out from her research this year.)

Ideal EOS for partially ionized hydrogen plasma

Plasma in the solar photosphere can be safely approximated by the ideal gas (non-interacting, randomly moving point-like particles) in the local thermodynamical equilibrium (TE). That means that

(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.

 

Ionization fraction

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.

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);

Notes on some basic quantities, units and constants: Part II

Mass fractions

The mass fraction of H, He and the metals is a quick way to specify the chemical composition of the stellar plasma. The mass fraction of hydrogen is defined as ratio of the mass of hydrogen particles and the total mass of all particles (in a given volume): $$ X = \frac{M_\mathrm{H}}{M}= \frac{\rho_\mathrm{H}}{\rho}, $$ where $\rho_\mathrm{H}$ and $\rho$ are the respective (mass) densities. The hydrogen mass density is equal to the hydrogen number density ($n_\mathrm{H}$) times the mass of one hydrogen particle $m_{\mathrm{H}} = A_{\mathrm{H}}\,m_{\mathscr{A}}$, thus: $$ X = \frac{M_\mathrm{H}}{M}= \frac{n_\mathrm{H}\,m_{\mathrm{H}}}{\rho}= \frac{n_\mathrm{H}\,A_{\mathrm{H}}\,m_{\mathscr{A}}}{\rho}. $$ Similar to that, for helium and the metals we define: $$ Y = \frac{M_\mathrm{He}}{M}= \frac{n_\mathrm{He}\,m_{\mathrm{He}}}{\rho}= \frac{n_\mathrm{He}\,A_{\mathrm{He}}\,m_{\mathscr{A}}}{\rho}, $$ $$ Z = \frac{M_\mathrm{metals}}{M}= \frac{n_\mathrm{metals}\,m_{\mathrm{metals}}}{\rho} = \frac{\sum_{i=3} n_\mathrm{i}\,A_{\mathrm{i}}\,m_{\mathscr{A}}}{\rho},. $$

Notes on some basic quantities, units and constants: Part I

There is a definition of mole as a unit that every student learns at very elementary level. Although mole is one of the seven base units of the International System, it is a bit specific and sometimes creates confusion. Here is my attempt to clarify the concept of mole and to derive in one place some useful relations used in the radiative transfer and atmospheric modeling. The first version of these notes I wrote for a group of students at the University of Belgrade many years ago. I still use them as a personal reminder. 

Counting "Elementary particles"

For the solar/stellar plasma, the "elementary" particles are atoms, ions (positive or negative), free electrons and molecules. In the very cool atmospheres there are dust particles as well. In the solar atmosphere dust can be completely neglected. Regarding the chemical composition, the atmospheric plasma is made out of hydrogen, helium and the metals (all other elements). There is no nuclear reactions and thus the total number density of nuclei per atomic specie is constant with time.

It is important to distinguish between the number of free atoms and the total number of atoms including those bound in the molecules. The former I denote as $N_{\mathrm{a}}^{\mathrm{free}}$, the latter as $N_{\mathrm{a}}^{\mathrm{tot}}$. The total number of atoms is identical to the number of atomic nuclei. The total number of molecules is $N_{\mathrm{m}}$.

The total number of particles $N$ is therefore:
$$N = N_\mathrm{e} + N_{\mathrm{a}}^{\mathrm{free}} + N_{\mathrm{m}},$$ or when there is no molecules $$N = N_\mathrm{e} + N_{\mathrm{a}}^{\mathrm{free}} = N_\mathrm{e} + N_\mathrm{H} + N_\mathrm{He} + N_\mathrm{\mathrm{metals}},$$ where $\mathrm{e}$, $\mathrm{m}$, $\mathrm{H}$ and $\mathrm{He}$ stand for the electrons, the molecules, hydrogen, helium and $\mathrm{metals}$ refer to all other elements together. The contribution of the metals can be further divided into the contributions of the individual elements.

Diatomic partition functions in IDL

Recently I posted a set of IDL routines for the atomic partition functions (PF) based on various polynomial fits that are available in the literature. Here I do the same for the diatomic molecules.

The routines for the diatomic PF use an internal catalog of 325 diatomic molecules. The ionized molecules are separate entries. The catalog is stored as an IDL structure in catalogue_of_diatomics.sav. The structure has the following tags for every molecule:

Name	Chemical formula, e.g. 'H2', 'CO', 'H2+'
Nc	Number of constituents, always 2 for diatomics
Nuc	Number of uniq constituents, e.g. 2 for CO, 1 for H2
Constituents	Constituent atoms, e.g. ['H', 'H'], ['C', 'O']
Charge	Charge, e.g. -1, 0, +1
Code	Kurucz's code for molecules, e.g. CO '0608.00'
D0	Dissociation energy in eV
Data	List of internally available data sets for this molecule
Type	Type of data: pf for partition function, kp for chemical equilibrium constant, eint for internal energy (one entry for each dataset)
References	List of ADS references (for each dataset)

Atomic partition functions in IDL

Many books on radiative transfer and spectroscopy introduce the partition function ($U$, PF) as a statistical weight or probability that an atom is in the ionization stage $r$ at the temperature $T$. It is defined as the sum over all excitation states ($i$) of their statistical weights ($g_i$) multiplied by an
exponential factor (cf. Rutten, 2003, p.30, Eq.289):
$$U_r \equiv \sum_i g_i \mathrm{e}^{-\frac{\chi_{r, i}}{kT}},$$ where $\chi_{r, i}$ is the excitation potential of the state $i$ in the ionization stage $r$, $k$ is the Boltzmann constant) and $T$ is the temperature.

The trouble with PF is that there is an infinite number of the energy levels. No matter how small their contribution becomes (the exponential factor is decreasing with the excitation potential), PF computed after this definition goes to infinity. Farther, that means that the probability of an atom to occupy any finite state would be zero. Since this is clearly unphysical, thus there is a need to define a reliable procedure to cut off the partition function. There are many approached described in literature. However, they are usually rather non-trivial and, therefore, computing partition function even for a limited range of temperature and for a specific atomic element can be a very tedious job. The results of these computations are often published as the coefficients of the best polynomial fits for a range of temperature values. The IDL procedure that you can find below is written to load these tables from several different sources (see below) and to interpolate them for given a temperature, an atomic specie and an ionization stage.