AST4310, Autumn 2020, Julia version

# Project 1: Basic Spectral Line Formation ("Cecilia Payne")¶

This project was originally written by Robert J. Rutten, and converted to notebook by Tiago M. D. Pereira.

#### Header and imports¶

The cells below contain some code to label equations in Markdown and some recommended julia imports to solve the exercises.

In :
macro javascript_str(s) display("text/javascript", s); end

javascript"""
MathJax.Hub.Config({
TeX: { equationNumbers: { autoNumber: "AMS" } }
});"""

In :
using Unitful
using UnitfulRecipes   # To make unitful play with plots, otherwise need to ustrip(temp)
using Plots
import PhysicalConstants.CODATA2018: c_0, k_B, m_e, h
import DelimitedFiles


## 1. Background¶

In this project you will explain the spectral-type sequence that is studied morphologically in the first exercise and that is summarised in the Figure below. A selection of stellar spectrograms illustrating the Harvard spectral sequence. These example spectra are printed positively, with the absorption lines dark on a bright background. Wavelengths in Ångström (1 Å = 0.1 nm = 10$^{-10}$ m). The peak brightness shifts from left to right from the "early-type" stars (O and B) to the "late-type" stars (G and lower). The sun has spectral type G2 V and is a late-type star. The early-type stars display the hydrogen Balmer lines prominently, but these become weak in solar-type spectra in which the Ca$^+$ H and K resonance lines are strongest. The M dwarfs on the bottom display strong molecular bands. From Novotny (1973).

You so re-act the work of Cecilia Payne at Harvard. Her 1925 PhD thesis was called "undoubtedly the most brilliant PhD thesis ever written in astronomy" by Otto Struve. Its opening sentences are:

The application of physics in the domain of astronomy constitutes a line of investigation that seems to possess almost unbounded possibilities. In the stars we examine matter in quantities and under conditions unattainable in the laboratory. The increase in scope is counterbalanced, however, by a serious limitation - the stars are not accessible to experiment, only to observation, and there is no very direct way to establish the validity of laws, deduced in the laboratory, when they are extrapolated to stellar conditions.

Extrapolation of terrestrial physics laws is precisely what Payne did in her thesis. She applied the newly derived Saha distribution for different ionisation stages of an element to stellar spectra, finding that the empirical Harvard classification represents primarily a temperature scale. Her work crowned efforts of Saha, Russell, Fowler, Milne, Pannekoek and others along the same lines. It illustrates that detailed physics, in this case atomic physics, is usually needed to explain cosmic phenomena. Cecilia Payne (1900 - 1979) was educated at Cambridge (England) by Milne and Eddington. She went to the US in 1923 and spent the rest of her career at Harvard (in the other Cambridge near Boston). Her 1925 thesis was the first one in astronomy at Harvard University and remains highly readable as a wide review of stellar spectroscopy at the time. The main conclusion was that stellar composition does not change much from star to star. Russell had already suggested so a decade earlier, but her thesis, the first Harvard Observatory Monograph, brought the point home. Copied from Hearnshaw (1986).

### 1.1 Payne’s line strength diagram¶

The key graph in Payne's thesis (page 131, earlier published in Payne (1924)) is reprinted below. Clearly, the observed behavior in the upper panel is qualitatively explained by the computed behavior in the lower panel. We will recompute the latter. The strengths of selected lines along the spectral sequence. Upper panel: variations of observed line strengths with spectral type in the Harvard sequence. The latter is plotted in reversed order on a non-linear scale that was obtained by making the peaks coincide with the corresponding peaks in the lower panel. The y-axis units are eye estimates on an arbitrary scale. Lower panel: Saha-Boltzmann predictions of the fractional concentration $N_{r,s}/N$ of the lower level of the lines indicated in the upper panel, each labeled with its ionisation stage, on logarithmic y-axis scales that are specified per species at the bottom, against temperature $T$ along the x axis given in units of 1000 K along the top.The pressure was taken constant at $P_e = N_e k T =$ 13.1 Pa. From Novotny (1973), who took it from Payne (1924).

### 1.2 The Boltzmann and Saha laws¶

In thermodynamical equilibrium (TE) macroscopic equipartition laws hold with the gas temperature as the major parameter. These are the Kirchhoff, Planck, Wien and Stefan-Boltzmann laws for radiation, and the Maxwell, Saha and Boltzmann laws for matter. In this exercise we are concerned with the latter two. They describe the division of the particles of a specific element over its different ionisation stages and over the discrete energy levels within each stage. For example, the Saha law specifies the distribution of iron particles between neutral iron (Fe), once-ionised iron (Fe$^+$), twice-ionised iron (Fe$^{2+}$), etc., whereas the Boltzmann law specifies the sub-distribution of the iron particles per ionisation stage over the discrete energy levels that each of the Fe, Fe$^+$, Fe$^{2+}$ etc. stages may occupy. (In astronomy one doesn't write ions as Fe$^{3+}$ but rather as Fe IV. More precisely: Fe I is the spectrum of neutral iron Fe, Fe II the spectrum of once-ionised iron Fe$^+$, etc.)

The figure below illustrates the energy level structure of neutral hydrogen. Energy level diagram for hydrogen. They approach the ionisation threshold at $\chi_H$ = 2.18 aJ for $n \rightarrow \infty$. The principal quantum number $n$ equals the level counter $s$ in this simple structure. The fine structure of each level (splitting in $2\,n^2$ sublevels) is not shown. For each of the first four hydrogen series the principal bound-bound transitions between bound levels are marked by vertical lines with the name and the wavelength of the corresponding spectral line. The series limits $(n=\infty$) are also marked. A bound-free ionisation/recombination transition is added to the Balmer series. The amount of energy above the ionisation threshold represents the kinetic energy that is gained or lost. A free-free transition (radiative encounter between a bare proton and a free electron) is also marked. The bound-free and free-free transitions contribute to stellar continua, while the bound-bound transitions produce the hydrogen lines. The Lyman lines are in the ultraviolet, the Balmer lines are in the visible and the Paschen and Brackett lines are in the infrared. Some Balmer lines are present in the stellar spectrograms shown above. The solar Balmer $\alpha$ line is usually called H$\alpha$.

#### 1.2.1 Boltzmann law¶

In TE the partitioning of a specific atom or ion stage over its discrete energy levels ("excitation equilibrium")

\begin{equation} \frac{n_{r,s}}{N_r} = \frac{g_{r,s}}{U_r} \mathrm{e}^{-\chi_{r,s}/kT}, \label{eq:5.16a} \end{equation}

with $T$ the temperature, $k$ the Boltzmann constant, $n_{r,s}$ the number of particles per m$^3$ in level $s$ of ionisation stage $r$, $g_{r,s}$ the statistical weight of that level, and $\chi_{r,s}$ the excitation energy of that level measured from the ground state $(r,1)$, $N_r \equiv \sum_s n_{r,s}$ the total particle density in all levels of ionisation stage $r$, and $U_r$ its partition function defined by

\begin{equation} U_r \equiv \sum_s g_{r,s} \mathrm{e}^{-\chi_{r,s}/kT}. \label{eq:5.16} \end{equation}

Thus, the neutral stage has $r=1$, each ground state is at $s=1$, and each ground state has excitation energy $\chi_{r,1}=0$ and ionisation energy to the next stage $\chi_r$. A radiative deexcitation between levels $(r,s)$ and $(r,t)$, with level $s$ 'higher' than level $t$, releases a photon with energy $\chi_{r,s} - \chi_{r,t} = h\nu = h c/\lambda$, with $h$ the Planck constant, $\nu$ the photon frequency, $c$ the velocity of light and $\lambda$ the wavelength. The excitation energy $\chi_{r,s}$ is the energy difference between the excited level $(r,s)$ and the ground state $(r,1)$. Astronomers usually call it "excitation potential" and measure it from the ground state up (Physicists prefer "binding energy" from the ground state of the next ion down in wavenumbers (m$^{-1}$)) in electron volt, with 1 eV corresponding to $1.6022 \times 10^{-19}$ J. For example, the H I Balmer $\alpha$ line results from photonic transitions between levels $n=2$ and $n=3$ of neutral hydrogen, with $\chi_{1,3} = 12.09$ eV, $\chi_{1,2} = 10.20$ eV and wavelength $\lambda = hc/(\chi_{1,3} - \chi_{1,2})= 656.5$ nm.

The number densities $n_{r,s}$ and $n_{r,t}$ are called "level populations" and are usually measured per m$^3$.

The statistical weights $g_{r,s}$ measure the degeneracy of levels due to magnetic fine splitting. The latter occurs only in the presence of an external magnetic field; in its absence, magnetic fine-structure levels coincide and may accommodate more particles than allocated per single level by the Pauli exclusion principle.
The weights measure such excess.
For example, neutral hydrogen atoms have $g_{1,1}=2$ for their ground state because the electron and proton spins can be parallel or anti-parallel (the fine-structure transition between the two states produces the 21 cm radio line from interstellar gas).

#### 1.2.2 Saha law¶

In TE the particle partitioning over the various ionisation stages of an element ("ionisation equilibrium") is given by the Saha distribution:

\begin{equation} \frac{N_{r+1}}{N_r} = \frac{1}{N_e} \frac{2U_{r+1}}{U_r} \left(\frac{2 \pi m_e kT}{h^2}\right)^{3/2} \mathrm{e}^{-\chi_r/kT}, \label{eq:5.17} \end{equation}

with $N_e$ the electron density, $m_e$ the electron mass, $\chi_r$ the threshold energy needed to ionise stage $r$ to stage $r+1$, and $U_{r+1}$ and $U_r$ the partition functions of ionisation stages $r+1$ and $r$ defined by (\ref{eq:5.16}). The ionisation energy $\chi_r$ is the minimum photon energy that is absorbed at ionisation or emitted at recombination in a bound-free interaction. The factor two represents the statistical weight of the freed electron, which has $g_e = 2$ due to the two orientations that its spin may take. The scaling with $1/N_e$ says that ionisation is easier if there is room for the resulting free electron or, reversedly, that recombination from stage $r+1$ to stage $r$ requires catching a free electron. The kinetic energy of the free electron contributes the $(\ldots)^{3/2}$ term through the Maxwell velocity distribution.

#### 1.2.3 Saha-Boltzmann Populations of Hydrogen¶

The level energies and statistical weights of hydrogen can be calculated exactly using the relations:

\begin{align} \tag{4} g_{1,s} &= 2s^2 \\ \tag{5} \chi_{1, s} &= 13.598 (1- 1/s^2)\; \mathrm{eV} \end{align}

By far the largest jump in energy is from the ground level to the first excited stage, as illustrated in the level diagram earlier. The first excited level $s=2$ is at such high excitation energy that its small Boltzmann factor makes its population negligible in comparison with the ground state population. However, the increase with $s$ at large values of $s$ shows that the hydrogen partition function would get infinite if too many levels are included in the summation, because $g_{1,s} \sim s^2$ while $\chi_{1,s} \rightarrow 13.598$ eV. Actually, all atoms and ions share in this behavior at very high excitation energy since they all get to be "hydrogenic" in nature when the valence electron sits in a nearly detached orbit. This singularity has been a cause of much debate, but real atoms are not worried by it. They are never alone and loose their identity through interactions with neighbours long before they grow as large as this. (They are never alone where TE holds. In intergalactic space they may be lonely but they won't sit in their high levels.) A reasonable cut-off value to the orbit size is set by the mean atomic interdistance $N^{-1/3}$ (page 260 of Rybicky & Lightman);

\begin{equation} \tag{6} s_\mathrm{max} \approx \sqrt{\frac{r}{a_0}} \,\,\, N^{-1/6}, \end{equation}

giving $s_\mathrm{max} \approx 100$ for hydrogen at $N_\mathrm{H} = 10^{18}$ m$^{-3}$. With such a cut-off the partition functions are generally not much larger than the ground state weights at all temperatures of interest - those at which the pertinent stage of ionisation is not devoid of population anyhow.

### Exercise 1: The Boltzmann and Saha laws¶

* Inspect the hydrogen energy level diagram in section 1.2. Which transitions correspond to the hydrogen lines in the image with stellar spectrograms (section 1)? Which transitions share lower levels and which share upper levels? * Payne's basic assumption was that the strength of the absorption lines observed in stellar spectra scales with the population density of the lower level of the corresponding transition. Why would she think so? (It is not correct, but generally stellar absorption lines do get stronger at larger lower-level population. In this exercise we follow her example and assume that the scaling is linear.) * Use this expectation to give initial rough estimates of the strength ratios of the $\alpha$ lines in the the H I Lyman, Balmer, Paschen and Brackett series. * Explain from equations (1) and (3) why the Saha and Boltzmann distributions behave differently for increasing temperature. * Speculate how ionisation can fully deplete a stage (e.g. all atoms can transition from neutral to ionised) while excitation puts only a few atoms in levels just below the ionisation level. Hint: what is the limit of the Saha and Boltzman ratios for infinite temperature?

### Exercise 2: Saha-Boltzmann populations of a simplified Ca atom¶

Using a model atom, we can calculate numerically the quantities given above in the Boltzmann and Saha laws. A model atom is simply a collection of the atomic data necessary: level energies $\chi_{r,s}$, ionisation energies $\chi_r$, and the statistical weights $g_{r,s}$. For temperature, density, and pressure we will use values similar to those on a stellar atmosphere. We start with a simplified model of calcium, with three levels and four states per level.

Why calcium? It is an important element, and later we will use it to make a comparison with hydrogen, by far the most abundant element in the universe. In this directory you will find model atoms for several elements (*_atom.txt). The atom file format is: one level per line, with the energy in wavenumbers (cm$^{-1}$), followed by the statistical weight and ionisation stage. In this exercise you will use Ca_atom.txt. To read the file format into a data structure you can manipulate, you can use the following Atom structure and read_atom() function:

In :
struct Atom{T <: Real}
g::Array{T}
chi::Array{<:Unitful.Energy{T}}
chi_ion::Array{<:Unitful.Energy{T}}
nstages::Int
max_levels::Int
end

"""
Reads atom structure from text file.

# Parameters

filename: string
Name of file with atomic data.
"""
function read_atom(filename)
tmp = DelimitedFiles.readdlm(filename, comments=true)
nstages = Int(maximum(tmp[:, 3])) + 1
max_levels = 0
for i=1:nstages
# Count all ocurrences when x = level number
max_levels = max(max_levels, count(x -> x == i - 1, tmp[:, 3]))
end
# Populate level energies and statistical weights
chi = zeros(Float64, nstages, max_levels)
g = deepcopy(chi)
for i=1:nstages
nlevels = count(x -> x == i - 1, tmp[:, 3])
chi[i, 1:nlevels] = tmp[findall(x -> x == i - 1, tmp[:, 3]), 1]
g[i, 1:nlevels] = tmp[findall(x -> x == i - 1, tmp[:, 3]), 2]
end
chi = (c_0 * h * chi * u"cm^-1") .|> u"aJ"  # unit conversion from cm^-1 to aJ
chi_ion = deepcopy(chi[:, 1])
chi .-= chi_ion
chi[chi .< 0u"aJ"] .= 0u"aJ"  # Reinstate original missing values as zero
return Atom(g, chi, chi_ion, nstages, max_levels)
end

Out:
read_atom

To solve this exercise, it is recommended you write functions to manipulate the structure Atom and calculate the partition function and populations according to Saha-Boltzmann.

* Using the simplified Ca atom (Ca_atom.txt), compute the partition functions $U_r$ for T=5000, 10000, and 20000 K. What can you say about the temperature dependence of $U_r$? * Plot a "Payne curve" for the simplified Ca atom using the same temperature range (100 - 175.000 K) and electron pressure (100 Pa) * Make a separate figure with a Payne curve for an element of your choice. You can choose one of the existing model atoms, or you can use the [NIST atomic spectra database](https://physics.nist.gov/PhysRefData/ASD/levels_form.html) to build a model for any atom you'd like. How does it compare with Ca?

### Exercise 3: Solar Ca$^+$K versus H$\alpha$¶

The figure below compares the solar spectrum of H$\alpha$, the principal line in the Balmer sequence (transition fron $n=$3 to $n=$2, $\lambda =$ 656.5 nm, see diagram in section 1.2) with the spectrum of Ca$^+$K (transition from $n=$2 to $n=$1, $\lambda =$ 393.5 nm), known in astronomical notation by Ca II K. (The K is an extension from Draper to the original alphabetic solar spectrum feature list by Fraunhofer, who only named Ca II H in his alphabetic feature naming from red to blue.)

The main line in the lefthand panel (Ca$^+$K) is much stronger than the one at right (H$\alpha$). Since stars as the sun are mostly made up of hydrogen, it may come as a surprise that a calcium line is much stronger than a hydrogen line. However, by now (in your role of being Cecilia Payne) it is clear to you that line strength ratios between different elements do not only depend on their abundance ratio but also on the temperature. We will quantify this dependence for this solar line pair. Two sections of the solar spectrum displaying the strongest lines in the visible region (except that Ca$^+$ H at $\lambda =$ 396.7 nm is very similar to its doublet twin Ca$^+$ K). The lefthand segment contains the central part of the Ca$^+$ K line in the violet region of the spectrum, the righthand segment the H$\alpha$ line in the red. Each segment is 3 nm wide. The vertical axis measures disk-averaged solar intensity so that these plots portray the solar spectrum as if it came from a non-resolved distant star. The units are kW m$^{-2}$ nm$^{-1}$, obtained by integrating the solar intensity over the solar disk (solid angle). The line crowding is much larger in the violet than in the red. Most of the numerous "blends" (overlapping lines) that are superimposed on the wings of Ca$^+$ K are due to iron, an element with extraordinary rich energy level structure. The strong blend at 394.5 nm is due to aluminium atoms. Very close to line center the Ca$^+$ K line displays two tiny emission peaks which betray the presence of magnetic fields in a way that is not understood. Stars that are more active than the sun have much higher peaks in their Ca$^+$ K line cores. The damping wings start just outside these peaks and extend much further than the plotted range. The H$\alpha$ line in the righthand panel is much weaker. The cores of both lines are formed in the solar chromosphere, the regime where magnetic fields take over from the gas pressure in dominating solar fine structure. There are hundreds of studies of cool-star chromospheres using the Ca$^+$ K line core because its emission is an indicator of stellar magnetism. The Ca$^+$ K line is also much used in recent studies of solar chromospheric dynamics. The temperature sensitivity of H$\alpha$ makes it complementary to Ca$^+$ K as an atmospheric diagnostic, especially of the magnetic field structures and their not-understood heating processes in the upper chromosphere. These show up as thread-like brightenings and darkenings arranged as flower petals in images taken in H$\alpha$ line-center radiation. These plots are made from the "Kitt Peak Solar Flux Atlas" compiled by Robert Kurucz, made with the Fourier Transform Spectrometer at Kitt Peak, a 1-metre Michelson interferometer with the unsurpassed spectral resolution of $\lambda / \Delta \lambda =$ 10$^6$.

To compute the relative strength between the Ca$^+$ K and H$\alpha$ lines you will use the code developed so far. Your Atom class can be used to compute the populations of any model atom, and in the current directory we have model atoms for H and Ca. Another necessary ingredient is the Ca abundance relative to hydrogen. In the Sun, this is about 2x10$^{-6}$. You will need to compare the number of absorbers on the lower level of each line: for Ca$^+$ K this is the ground state of the singly-ionised atom, while for H$\alpha$ it is the first excited state, neutral atom.

* Explain qualitatively why the solar Ca$^+$ K line is much stronger than the solar H$\alpha$ line, even though hydrogen is not ionised in the solar photosphere and low chromosphere ($T \approx$ 4000 - 6000 K) where these lines are formed, and even though the solar Ca/H abundance ratio is only $N_\mathrm{Ca}/N_\mathrm{H} = 2 \times 10^{−6}$. Assume again that the observed line strength scales with the lower-level population density (which it does, although nonlinearly through a "curve of growth" as you will see in Project 2). * Prove your explanation by computing and plotting the expected strength ratio of these two lines as function of temperature for $P_e = 10^2$ dyne cm$^{-2}$. Make use of H_atom.txt and Ca_atom.txt. * Plot the relative population changes $(\Delta n_\mathrm{Ca} / \Delta T) / n_\mathrm{Ca}$ and $(\Delta n_\mathrm{H} / \Delta T) / n_\mathrm{H}$ for the lower levels of Ca$^+$K and H$\alpha$, using $\Delta T=$ 1 K. Around $T=$ 5600 K the Ca$^+$K curve dips down to very small values; the H$\alpha$ curve does that around $T=$ 9500 K. Thus, for $T \approx$ 5600 K the temperature sensitivity of Ca$^+$K is much smaller than the temperature sensitivity of H$\alpha$. Explain each flank of the two population curves and the dips in the two temperature sensitivity curves. (The dips can be diagnosed by overplotting the variation with temperature of each population in relative units.) * Find at which temperature the hydrogen in stellar photospheres with $P_e =$ 10 Pa is about 50% ionised. Plot the neutral and ionised fractions of hydrogen as a function of temperature.