**AST4310, Autumn 2020, Python version**

This project was originally written by Robert J. Rutten, and converted to notebook by Tiago M. D. Pereira, using contributions from Luc Rouppe van der Voort, Lluís Mas Ribas, and Henrik Eklund.

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

In [1]:

```
%%javascript # Allow equation numbers
MathJax.Hub.Config({
TeX: { equationNumbers: { autoNumber: "AMS" } }
});
```

In [2]:

```
%matplotlib inline
import numpy
import matplotlib.pyplot as plt
from astropy import units
from astropy import constants
from astropy.visualization import quantity_support
from IPython.display import set_matplotlib_formats
quantity_support()
set_matplotlib_formats('svg')
plt.rc('legend', frameon=False)
```

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

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

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

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

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

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

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.

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.

* 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. You can use T=5000 K.
* 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?

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`

class with a method `read_atom()`

:

In [3]:

```
class Atom:
"""
Reads atomic data, calculates level populations according to Boltzmann's law,
and ionisation fractions according to Saha's law.
"""
def __init__(self, atomfile=None):
"""
Parameters
----------
atomfile : string, optional
Name of file with atomic data. If not present, atomic data needs
to be loaded with the .read_atom method.
"""
self.loaded = False
if atomfile:
self.read_atom(atomfile)
def read_atom(self, filename):
"""
Reads atom structure from text file.
Parameters
----------
filename: string
Name of file with atomic data.
"""
tmp = numpy.loadtxt(filename, unpack=True)
self.n_stages = int(tmp[2].max()) + 1
# Get maximum number of levels in any stage
self.max_levels = 0
for i in range(self.n_stages):
self.max_levels = max(self.max_levels, (tmp[2] == i).sum())
# Populate level energies and statistical weights
# Use a square array filled with NaNs for non-existing levels
chi = numpy.empty((self.n_stages, self.max_levels))
chi.fill(numpy.nan)
self.g = numpy.copy(chi)
for i in range(self.n_stages):
nlevels = (tmp[2] == i).sum()
chi[i, :nlevels] = tmp[0][tmp[2] == i]
self.g[i, :nlevels] = tmp[1][tmp[2] == i]
# Put units, convert from cm-1 to Joule
chi = (chi / units.cm).to('aJ', equivalencies=units.spectral())
# Save ionisation energies, saved as energy of first level in each stage
self.chi_ion = chi[:, 0].copy()
# Save level energies relative to ground level in each stage
self.chi = chi - self.chi_ion[:, numpy.newaxis]
self.loaded = True
```

To solve this exercise, it is recommended you add a few methods to class `Atom`

to 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 a temperature range of 100 - 175,000 K and $P_e$ = 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?

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.