NGC 7319: GTC vs JWST

The first pictures from James Webb Space Telescope are simply mind blowing! Here is a comparison of NGC 7319 (from Stephan's quintet of galaxies) from Gran Telescopio Canarias (still the largest single aperture optical telescope) and from JWST. One should keep in mind that these are different instruments observing and different wavelengths. The alignment is also not perfect.

Quadratic splines: Part I

Some basic properties of parabola in the context of the second-order Hermite interpolation equation. 

There is a one-dimensional dataset with $n$ points: $x = (x_0, \dots, x_{n-1})$ and $y = (y_0, \dots, y_{n-1})$, such that $x_{i} > x_{i-1}$ for any $i$. We introduce the following labels: $$h_{i-1} = x_i - x_{i-1}$$ $$\delta_{i-1} = \frac{y_i - y_{i-1}}{x_i - x_{i-1}}$$ $h$ is always positive, while $\delta$ can have any value. On each segment, we approximate the unknown function with a parabola. On the $[i-1, i]$ segment: \begin{equation}f(x) = a_0 + a_1 (x - x_{i-1}) + a_2 (x - x_{i-1})^2\end{equation} The first and the second derivative of the parabola are: \begin{equation}f'(x) = a_1 + 2 a_2 (x - x_{i-1})\end{equation} \begin{equation}f''(x) = 2 a_2 \end{equation} The Hermite interpolation theorem provides us with a tool to determine the $a$ coefficients on the $[i-1, i]$ segment if we know $y_{i-1}$, $y_i$ and one of the first derivative in the end points, either $y'_{i-1}$ or $y'_i$. Let's assume that we know ("know" in the sense "we can compute") $y'_{i-1}$. Then there are three equations that the parabola must obey: \begin{equation}f(x_{i-1}) = y_{i-1} = a_0\end{equation} \begin{equation}f'(x_{i-1}) = y'_{i-1} = a_1\end{equation} \begin{equation}f(x_{i}) = y_i = a_0 + a_1 h_{i-1} + a_2 h_{i-1}^2\end{equation} The solution is then trivial: \begin{equation}a_0 = y_{i-1}\end{equation} \begin{equation}a_1 = y'_{i-1}\end{equation} \begin{equation}a_2 = \frac{y_{i} - y_{i-1} - y'_{i-1} h_{i-1}}{h_{i-1}^2}\end{equation} Once the $a$ coefficients are known, the parabola is fully determined. That also means that the derivative in the point $i$ is fixed: $$f'(x_{i}) = y'_i = a_1 + 2 a_2 (x_i - x_{i-1}) = y'_{i-1} + 2 \frac{y_{i} - y_{i-1} - y'_{i-1} h_{i-1}}{h_{i-1}^2} h_{i-1}$$ $$y'_i = \frac{2(y_{i} - y_{i-1}) - y'_{i-1} h_{i-1}}{h_{i-1}} $$ \begin{equation}y'_i = 2\delta_{i-1} - y'_{i-1}\end{equation} Local extremum

The extremum of this parabola is located at $x_m$: $$f'(x_m) = 0 = a_1 + 2 a_2 (x_m - x_{i-1})$$ \begin{equation}x_m = x_{i-1} - \frac{a_1}{2a_2}\end{equation} Monotonicity

The function $f(x)$ is monotone on the $[i-1, i]$ segment if it does not have an extremum on it, i.e. if $$x_m < x_{i-1}\;\;\;\;\;\;\;\;\;\;\mathrm{or}\;\;\;\;\;\;\;\;\;\; x_m > x_i$$ The first condition readily translates into: $$x_{i-1} - \frac{a_1}{2a_2} < x_{i-1}$$ $$ \frac{a_1}{2a_2} > 0$$ or, when we substitute $a_1$ and $a_2$: $$ \frac{y'_{i-1} h_{i-1}^2}{2(y_{i} - y_{i-1} - y'_{i-1} h_{i-1})} > 0$$ $$ \frac{y'_{i-1} h_{i-1}}{2(\delta_{i-1} - y'_{i-1})} > 0$$ As by definition $h_{i-1} > 0$ $$ \frac{y'_{i-1}}{2(\delta_{i-1} - y'_{i-1})} > 0$$ If $y'_{i-1} > 0$ then $\delta_{i-1} - y'_{i-1} > 0$ and $\delta_{i-1} > y'_{i-1}$. If $y'_{i-1} < 0$, then $\delta_{i-1} < y'_{i-1}$. Therefore, the first condition ($x_m < x_{i-1}$) reduces to: \begin{equation}\boxed{|\delta_{i-1}| > |y'_{i-1}|}\end{equation} The second condition leads to: $$x_{i-1} - \frac{a_1}{2a_2} > x_{i}$$ $$-x_{i} + x_{i-1} - \frac{a_1}{2a_2} > 0$$ $$h_{i-1} + \frac{a_1}{2a_2} < 0$$ $$h_{i-1} + \frac{y'_{i-1} h_{i-1}}{2(\delta_{i-1} - y'_{i-1})} < 0$$ $$\frac{2\delta_{i-1} - y'_{i-1}}{2(\delta_{i-1} - y'_{i-1})} < 0$$ There are two cases. First, if $2\delta_{i-1} - y'_{i-1} > 0$, i.e. $2\delta_{i-1} > y'_{i-1}$. then it must be $\delta_{i-1} - y'_{i-1}< 0$, i.e.$\delta_{i-1} < y'_{i-1}$. Therefore, the first case leads to $2\delta_{i-1} > y'_{i-1} > \delta_{i-1}$. Secondly, if $2\delta_{i-1} - y'_{i-1} < 0$, i.e. $2\delta_{i-1} < y'_{i-1}$, it must be $\delta_{i-1} - y'_{i-1}> 0$, i.e. $\delta_{i-1} > y'_{i-1}$. So the second case leades to $\delta_{i-1} > y'_{i-1}> 2\delta_{i-1}$. The first case is possible only when both $\delta_{i-1}$ and $y'_{i-1}$ are positive and the second case is possible only when both are negative. Therefore, together, the two cases reduce the second condition ($x_m > x_{i}$) to: \begin{equation}\boxed{|2\delta_{i-1}| > |y'_{i-1}| > |\delta_{i-1}|}\end{equation} Together, the two conditions combine into \begin{equation}\boxed{|2\delta_{i-1}| > |y'_{i-1}|}\end{equation} A much faster way to reach the same result is to realize that the the derivatives $y'_{i-1}$ and $y'_{i}$ must be of the same sign if the function is monotone on the $[i-1, i]$ segment, i.e. $$y'_{i-1} \cdot (2\delta_{i-1} - y'_{i-1}) > 0$$ which obviously leads to $|2\delta_{i-1}| > |y'_{i-1}|$.

Therefore, for a Hermite parabola defined by $y_{i-1}$, $y_i$ and $y'_{i-1}$ to be monotone on the interval $(i-1, i)$, the derivative $y'_{i-1}$ must be of the same sign as the linear slope $\delta_{i-1}$ and smaller in magnitude than $2\delta_{i-1}$.

Convexity

For the function $f$ to be convex on $(i-1, i)$ (if parabola is convex on a segment, it's convex everywhere) the condition is $f''(x) > 0$. In the case of our parabola: $$f''(x) = 2a_2 > 0$$. $$\frac{y_{i} - y_{i-1} - y'_{i-1} h_{i-1}}{h_{i-1}^2} > 0$$ or $$y_{i} - y_{i-1} - y'_{i-1} h_{i-1} > 0$$ \begin{equation}\delta_{i-1} > y'_{i-1}\end{equation} This result is obvious $\delta_{i-1} = y'_{i-1}$ is the straight line going through the end-points of the segment. For the parabola to be concave, it must hold $f''(x) < 0$, i.e. \begin{equation}\delta_{i-1} < y'_{i-1}\end{equation} 

 

XXXI Canary Islands Winter School: Computational Fluid Dynamics in Astrophysics (Photos)

Two weeks of interesting lectures, intensive hands-on exercises, discussions and meeting new people... A couple of photos by Claudio Dalla Vecchia that perfectly capture the spirit of the school. 

XXXI Canary Islands Winter School: Computational Fluid Dynamics in Astrophysics

Registration is now open for the XXXI Canary Islands Winter School of Astrophysics, which will be held this year in La Laguna, Tenerife, Spain from 19-28 November, on the topic of "Computational Fluid Dynamics in Astrophysics". Confirmed invited lecturers are

Fernando Moreno-Insertis (Fluid dynamics, conservation laws and waves)
Ake Nordlund (Stellar and planetary formation)
Maarit Käpylä (Dynamics of solar and stellar convection zones)
Matthias Rempel (Solar Atmosphere)
Sascha Husa (General relativity and compact objects)
Stefanie Walch-Gassner (Interstellar matter)
Tom Theus (Galaxy evolution)

The school will consist of lectures, seminars tutorials, and social events including visits to the observatories in Tenerife and La Palma. It is aimed at advanced graduate students in astrophysics, as well as postdoctoral researchers, with a strong interest in the computational fluid dynamics and it applications astrophysics.

For further information and to register, please visit

http://www.iac.es/winterschool/2019/

Claudio Dalla Vecchia
Nikola Vitas
Elena Khomenko

Three-dimensional simulations of solar magneto-convection including effects of partial ionization

Astronomy & Astrophysics Volume 618 (October 2018) is just out with a figure from our paper on the cover!

The paper (Khomenko, Vitas, Collados & de Vicente 2018: A&A, ADS) describes the effects of the partial ionization on the structure, dynamics and energy balance of the low chromosphere.

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.

How to display magnetic field lines in IDL?

It is a common problem in visualization of magnetic fields. If we assume that the potential $A$ of the magnetic field is knows, the field lines are, by definition, iso-$A$ curves. The quickest way to display them is by using the CONTOUR procedure. If the potential field is given as a 2D variable potential, then:

IDL> CONTOUR, potential, levels = levels, /xs, /ys

gives:




However, the information here is not complete without showing the actual direction of the field along the lines. It is easy to do it in IDL:

IDL> CONTOUR, potential, levels = levels, /xs, /ys
IDL> CONTOUR, potential, levels = levels, /xs, /ys, path_xy = c
IDL> FOR i = 1, N_ELEMENTS(c)/2-1, 50 DO $     ARROW, c[0, i-1], c[1, i-1], c[0, i], c[1, i], /norm, /solid, hsize = 5         


produces:

It is much better now, but the arrow heads from different field lines make it crowded. If we are interested in individual field lines, we can add color to this: