Tuesday 11 July 2017

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}


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}.$$


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:
c_p = c_V + \frac{k}{m_\mathrm{H}}\,(1+x)\,\frac{(1 + D\,H)^2}{1-D}.
This expression can be written explicitly in terms of $x$ as:
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]

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

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}


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})}$$

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.

Tuesday 4 July 2017

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

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         

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:

Tuesday 11 April 2017

Deep-learning about horizontal velocities at the solar surface

The velocity fields are of great importance for understanding dynamics and structure of the solar atmosphere. The line of sight velocities are coded in the wavelength shifts of the spectral lines, thanks to the Doppler effect, and relatively easy to measure. On the other hand, the orthogonal ("horizontal") components of the velocity vector are impossible to measure directly.

The most popular method for estimating the horizontal velocities is so-called local correlation tracking (LCT, November & Simon, 1988). It is based on comparing successive images of the solar surface in the continuum light and transforming their differences into information about the horizontal fields. However, the LCT algorithm suffers from several limitations.

In a paper by Andres Asensio Ramos and Iker S. Requerey (with a small contribution from my side) accepted by A&A and published on Arxiv some weeks ago (2017arXiv170305128A) this problem is tackled by the deep-learning approach. A deep fully convolutional neural network is trained on synthetic observations from 3D MHD simulations of the solar photosphere and then applied to the real observation with the IMaX instrument on board the SUNRISE balloon (Martinez Pillet et al, 2011; Solanki, 2010). The method is validated using simulation snapshots of the quiet sun produced with the MANCHA code that I have been developing in the last couple of years.

Monday 20 March 2017

K-means clustering

The problem of clustering is a rather general one: If one has $m$ observations or measurements in $n$ dimensional space, how to identify $k$ clusters (classes, groups, types) of measurements and their centroids (representatives)?

The k-means method is extremely simple, rather robust and widely used in it numerous variants. It is essentially very similar (but not identical) to Lloyd's algorithm (aka Voronoi relaxation or interpolation used in computer sciences).


Let's use the following indices: $i$ counts measurements, $i \in [0, m-1]$; $j$ counts dimensions, $j \in [0, n-1]$; $l$ counts clusters, $l \in [0, k-1]$.

Each measurement in $n$-dimensional space is represented by a vector $x_i = \{x_{i, 0}, \dots x_{i, n-1}\}$, where index $i$ is counting different measurements ($i = 0, \dots, m-1$). The algorithm can be summarized as:

1. Choose randomly $k$ measurements as initial cluster centers: $c_0, ..., c_{k-1}$. Obviously, each of the clusters is also $n$-dimensional vector.

2. Compute Euclidean distance $D_{i, l}$ between every measurement $x_i$ and every cluster center $c_l$:
$$D_{i, l} = \sqrt{\sum_{j=0}^{n-1} (x_{i, j} - c_{l, j})^2}.$$
3. Assign every measurement $x_i$ to the cluster represented by the closest cluster center $c_l$.

4. Now compute new cluster centers by simply averaging all the measurements in each cluster.

5. Go back to 2. and keep iterating until none of the measurements changes its cluster in two successive iterations.

This procedure is initiated randomly and the result will be slightly different in every run. The result of clustering (and the actual number of necessary iteration) significantly depends on the initial choice of cluster centers. The easiest way to improve the algorithm is to improve the initial choice, i.e. to alter only the step 1. and then to iterate as before. There are to simple alternatives for the initialization.

Wednesday 1 February 2017

First observation of linear polarization in the forbidden [OI] 630.03 nm line

In a new paper (de Wijn, Socas-Navarro & Vitas, 2017, ApJ, 836, 29D) we present the first results of our observations of a sunspot and an active region using the SP/SOT instrument on board the Hinode satellite. The novelty in our observation is a trick that we used to double the standard wavelength range observed by the instrument. Thanks to that, we were able to see the sun not only in the two iron lines at 630.2 nm, but also in four other lines. One of those is particularly interesting: the forbidden ground-based line of neutral oxygen ([OI] 630.03 nm). It is one of only few oxygen lines in the solar spectrum and probably the best diagnostics of the solar oxygen abundance. For the first time ever we observed the linear polarization in this line! As an M2 (magnetic dipole) transition, it is predicted by the theory (Landi degl'Innocenti and Landi, 2004, Section 6.8) that this line produces the linear polarization signal with the opposite sign to the lines produced by E2 transitions. It is also the first time that linear polarization in M2 and nearby E2 lines is measured simultaneously, so that the flip in sign is obvious (see the left-most spectral line in the red circle in the Figure; in linear polarization it has "W" shape, while all other lines in the wavelength range have "M" shapes). This result may bring new light to the ongoing debate on the solar oxygen crisis.


More details of this unique observation will appear soon in a follow-up publication.