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

Charles Hermite (1822 - 1901)

Sunday morning in Paris was opportunity to walk to the Montparnasse cemetery and pay respect to some of my personal heroes buried there. While the graves of Beckett, Cortázar and Poincaré attract quite many attention and visitors, it is less known that the great French mathematician Charles Hermite is buried there as well. Not only that he was a mentor to Henri Poincaré, Henri Padé, Thomas Stieltjes and Mihajlo Petrovic Alas, but his work on interpolation and function approximation is at the very core of the modern numerical methods used in computational fluid dynamics and radiative transfer (even when this is not so obvious or properly acknowledged). His work from 1878 ("Sur la formule d'interpolation de Lagrange") should be read by anyone interested in function approximation. The name on his grave stone are barely readable nowadays.

Robert Jelle Rutten (1942 - 2022)

So long, my teacher and friend. You will be hugely missed.

Effective Lande g-factor (in IDL)

The effective Lande g-factor (Shenstone and Blaire, 1929) is defined for a spectral line. It shows how much the $\sigma$ components of the normal Zeeman triplet are separated from the line center.  The higher the Lande factor is, the line is more sensitive to the magnetic field. It is computed as a combination of g-factors for each of the involved atomic levels. Derivation of the equation can be found, for example, in Landi Degl'Innocenti (1982SoPh...77..285L, see it for more details and for an alternative expression useful in the polarized radiative transfer).

For each of the levels (u for upper, l for lower) in LS coupling holds:
$$g_\mathrm{LS} =\frac{3}{2}+\frac{S (S+1) - L (L+1)}{2J(J+1)}, $$ 
where $L$, $S$ and $J$ are the orbitatl, spin and total angular momentum quantum numbers of the atomic level (note that it is undefined for $J=0$). The effective Lande g-factor is then defined as:
$$\bar{g} = \frac{1}{2}(g_\mathrm{l} +g_\mathrm{u}) + \frac{1}{4}(g_\mathrm{l} - g_\mathrm{u})(J_\mathrm{l}(J_\mathrm{l}+1) - J_\mathrm{u}(J_\mathrm{u}+1)).$$

Here is my IDL code for computing the effective Lande g-factor of a spectral line in LS coupling. At the input the code takes the quantum number of the lower and the upper level of a transition. The numbers may be specified either as an array [L, S, J] or in spectroscopic notation.

Example:

IDL> PRINT, LANDE_FACTOR(lower  = [3., 1, 2.], upper = '3p1.0')
  Upper term: 3p1.0 => l = 1.00000s = 1.00000j = 1.00000
          0.250000

IDL> PRINT, LANDE_FACTOR(lower  = '5F1.0', upper = '5d0.0')
  Lower term: 5f1.0 => l = 3.00000s = 2.00000j = 1.00000
  Upper term: 5d0.0 => l = 2.00000s = 2.00000j = 0.00000
          0.00000

Download: lande_factor.pro

Kurucz' Atlas in GNU Linux

I have just found about great effort made by Sbordone, Bonifacio & Castelli to make the Atlas code of Bob Kurucz available at GNU Linux.

There is a webpage: http://atmos.obspm.fr/ providing the code and the documentation.

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.

1D Solar Models in IDL: Spruit's Convection Model

The semi-empirical models of the solar atmosphere rely on the observed intensities and, therefore, they cannot say much about the invisible convection zone below the surface. Spruit (1974SoPh...34..277S) constructed a 1D model of the solar convection zone using the mixing-length theory with 4 free parameters. The model is constructed so that it matches the HSRA model (Gingerich et al, 1971SoPh...18..347G) of the solar atmosphere. The model of Spruit is sometimes used to initiate the convection simulations with 3D numerical codes.

Fig 1. The HSRA model atmosphere (dashed) and Spruit's model of convective zone (solid). This is reproduced after Fig.3 of Spruit's paper.

1D Solar Atmosphere Models in IDL: Gingerich et al (HSRA)

The Harvard-Smithsonian Reference Atmosphere (HSRA) by Gingerich et al (1971, 1971SoPh...18..347G) is another widely used semi-empirical model of the idealized plane-parallel solar atmosphere in hydrostatic equilibrium. This model also includes the chromosphere where the hydrogen ionization is solved using the statistical equilibrium equations. The model is still often used as the reference solar atmosphere in 1D and its number of citations steadily grows (> 820 so far). The paper is very clearly written and it's a very recommendable reading especially if you use this model. Beside the model itself, there is several tables with computed intensities, brightness temperature, optical depths at different wavelengths, etc.

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.

Kourganoff graphs

Scharzschild equation provides a tool to easily find and analyze the solution of the radiative transfer cases for some simple source functions. First of all, the lambda operator is clearly linear, so\begin{equation}
\Lambda[a\,f(t) + b\,g(t)] = a\,\Lambda[f(t)] + b\,\Lambda[g(t)].
\end{equation}

Formal solution and Schwarzschild equation

The radiation transport equation (RTE) for time-independent plane-parallel atmosphere has well known form (e.g. Rutten, 2003):
\begin{equation}
\mu\,\frac{\mathrm{d}I_{\nu}}{\mathrm{d} \tau_{\nu}} = I_{\nu} - S_{\nu},
\end{equation}
 where $I_{\nu}$ is specific intensity, $S_{\nu}$ source function, $\tau_{\nu}$ is the radial optical depth and $\mu$ is the cosine of the angle between the line of sight and vertical axis. Mathematically, this is a linear first order differential equation with constant coefficients. Physically, it describes propagation of radiation through an interacting medium. Properties of the medium are hidden in the source function and the optical depth. However, to make the first step toward the solution of RTE, we will assume that the optical depth is independent variable and that the source function is a know function of it.

Exponential integrals

The formal solution of the radiative transfer equation can be expressed using exponential integrals. These are mathematical functions defined as:
\begin{equation}
E_n(x) = \int_1^{\infty}\, \exp (-yx)\frac{\mathrm{d}x}{y^n}, \;\;\;\;n = 1, 2, ...
\end{equation}
There is several important properties of these functions:

Carlson's Quadrature Set B

Beside the quadrature Set A that is described and computed in the previous blog, Carlson (1963) described an alternative set of quadratures, Set B. This set corresponds to odd approximations. It "relates to Set A in the way Double-P relates to Single-P." (ibid). The principles on which this set is constructed are same as in the case of Set A. However, the equations slightly differ. First we determine the $W_l$ coefficients from the system:
$$ \sum_{l=1}^{n/2-1} W_l = \frac{n-3}{3}, $$ $$ W_l^2 = W_1^2 + (l-1)\Delta^\prime, $$ where $\Delta^\prime = 2/n$. This system is again solved by Newton-Raphson formula. The $W_l$ values are then used to get $w_l$ (as with Set A).

Carlson's Quadrature Set A

When the radiative transfer equation is solved in three dimensions, it is necessary to discretize the angular dependence of the specific intensity. In principle, the continuous integral over solid angle is replaced by a discrete sum over selected rays. This is the often used technique known as the discrete ordinates method (DOM). However, there is no unique solution for the discretization of the solid angle.

The state-of-art radiative magnetohydrodynamical codes (eg. MURaM and BIFROST) use a quadrature known as Set A by Carlson (1963). The advantage of this particular quadrature is that it is invariant with respect to 90^o axis rotation.

Terminology

The quadrature is defined by directional cosines $\mu_l$ and their weights $w_l$. They are combined into sets that define rays $\Omega_m = (\mu_x, \mu_y, \mu_z)$ and their weights $w_m$. The order of quadrature $n$ specifies the number of $\mu$ levels in the [-1, 1] interval. Since the quadrature is symmetric in respect to axis, it is sufficient to find values of the angles in half of that range. Thus the number of $\mu$ levels is $l_{max} = n/2$. The number of rays per octant in this scheme is $n(n+2)/8$. Each ray is labeled by three numbers indexing $\mu$ level along one of the axes. The sum of indices of a ray is the same of all rays, $l_{max}+2$. For example, when $n=8$ and $l_{max}=4$, possible indices are 1, 2, 3 and 4. Their possible variations are those that have sum of 6. The number of these variations (and the rays) is $n(n+2)/8 = 10$. The number of different classes (variation without repetition with the specified sum) is $l_{max}-1 = 3$. In our example: (4, 1, 1), (3, 2, 1) and (2, 2, 2).