The spectrum of an optically thin plasma#

utu.spectrum computes the emission lines of a hot, thin plasma from the CHIANTI atomic database, which it reads through fiasco. A line is worth an intensity only once the plasma emitting it is described, so this page builds one: a column of material spread over a range of temperatures at a fixed pressure. Out of it comes the spectrum that column would show a spectrograph, which is how an instrument is designed around the lines it means to observe.

The ions here are named rather than discovered. ions() will find every ion with a line in a range of wavelengths, but the answer depends on which ions the database at hand describes, and the one these pages are built against describes only a few. Naming them keeps this page saying the same thing wherever it is run.

import astropy.units as u
import fiasco
import matplotlib.pyplot as plt
import named_arrays as na
import numpy as np

import utu

# fiasco says which approximations it has fallen back on for each ion,
# which is worth reading once and not seven times
fiasco.log.setLevel("ERROR")

Naming ions#

CHIANTI names an ion by its element and its charge state, counting from one, so that the neutral atom is 1. Spectroscopists write the same thing with a Roman numeral, and spectroscopic() translates between the two. It takes a name or an array of them.

ion = ["He 1", "O 3", "O 4", "O 5", "Ne 7", "Mg 10", "Si 11"]

[utu.spectrum.spectroscopic(i) for i in ion]

['He I', 'O III', 'O IV', 'O V', 'Ne VII', 'Mg X', 'Si XI']

The plasma#

The plasma is described by three arrays which share an axis of temperature.

The density follows from the temperature and a pressure, since the corona is near enough to hydrostatic that the two are not independent. Giving the density the axis name of the temperature is what says so: utu.spectrum reads a density on the axis of the temperature as one density per temperature, and a density on an axis of its own as every density at every temperature.

The emission measure is how much material there is at each temperature, \(\int n_e n_H \, dh\) through the column. It is the product of the densities of electrons and of hydrogen rather than the square of the first, since fiasco leaves the ratio of the two out of its contribution functions. A real one is measured; this one is a Gaussian in \(\log T\) centred where the transition region emits, holding \(10^{27}\,\mathrm{cm}^{-5}\) of material in all.

axis = "temperature"

temperature = na.geomspace(1e4, 1e7, axis=axis, num=61) * u.K

pressure = 1e15 * u.K / u.cm ** 3
density = pressure / temperature

weight = np.exp(-np.square((np.log10(temperature / u.K) - 5.8) / 0.3) / 2)
emission_measure = 1e27 / u.cm ** 5 * weight / weight.sum(axis)

emission_measure.sum(axis)

ScalarArray(
    ndarray=1.e+27 1 / cm5,
    shape={},
)

Contribution functions#

The contribution function \(G(T)\) of a line is the power it radiates per unit emission measure, and it is sharply peaked: a line is emitted only by plasma near the temperature at which its ion exists at all. This is what makes a line a thermometer, and what makes the choice of a passband a choice about which part of the atmosphere to look at.

contribution_function() returns every line of one ion at once, as a function of wavelength, so that a line and its strength cannot come apart.

result = utu.spectrum.contribution_function(
    ion=fiasco.Ion("O 5", temperature.ndarray),
    density=density,
    axis_temperature=axis,
)

result.outputs.shape

{'temperature': 61, 'line': 2995}

The strongest of them is the line at 629.7 Å, the line ESIS was built to observe, and it is emitted at a quarter of a million kelvin.

index = np.argmax(result.outputs.max(axis), axis="line")

fig, ax = plt.subplots(figsize=(6, 3), constrained_layout=True)
na.plt.plot(temperature, result.outputs[index], ax=ax, axis=axis)
ax.set_xscale("log")
ax.set_xlabel(f"temperature ({temperature.unit:latex_inline})")
ax.set_ylabel(f"$G(T)$ ({result.outputs.unit:latex_inline})")
ax.set_title(
    f"{utu.spectrum.spectroscopic('O 5')} "
    f"{result.inputs[index].ndarray:latex_inline}"
);
../_images/spectrum_4_0.png

The lines#

lines() does this for every ion given to it, weights each contribution function by the emission measure, sums over temperature, and returns what is left brightest first. The wavelength and the ion of a line are components of the inputs, and its intensity is the output, so that sorting or selecting lines carries all three together.

result = utu.spectrum.lines(
    temperature=temperature,
    density=density,
    emission_measure=emission_measure,
    wavelength_min=550 * u.AA,
    wavelength_max=680 * u.AA,
    ions=ion,
)

result.shape

{'line': 1136}

The emission is isotropic, so the radiance of a line is its intensity spread over the whole sphere.

num = 8
brightest = result[{"line": slice(num)}]

radiance = brightest.outputs / (4 * np.pi * u.sr)
radiance = radiance.to(u.erg / u.s / u.cm ** 2 / u.sr)

name = utu.spectrum.spectroscopic(brightest.inputs.ion)
wavelength = brightest.inputs.wavelength

print(f"{'ion':>6s}{'wavelength':>14s}{'radiance':>14s}")
for i in range(num):
    j = {"line": i}
    print(f"{name[j].ndarray:>6s}"
          f"{wavelength[j].ndarray.to_value(u.AA):>14.2f}"
          f"{radiance[j].ndarray.value:>14.1f}")

   ion    wavelength      radiance
   O V        629.73        2652.8
  O IV        554.51         387.0
  O IV        554.08         150.4
  O IV        609.83         117.3
  Mg X        609.79         109.7
  O IV        555.26          77.9
  O IV        553.33          77.4
  O IV        608.40          63.3

The spectrum#

Drawn as a stem from zero for each line, this is the picture an instrument is designed against: where the light is, and how much of it there is beside everything else the passband will admit. stem() draws it and names the brightest lines, moving the labels apart so that none covers another or crosses a line it does not belong to.

It takes a spectrum in the shape lines() returns, so the radiance computed above has to be put back into one before it can be drawn.

fig, ax = plt.subplots(figsize=(7, 3), constrained_layout=True)
utu.spectrum.stem(
    na.FunctionArray(inputs=brightest.inputs, outputs=radiance),
    ax=ax,
)
ax.set_xlabel(f"wavelength ({u.AA:latex_inline})")
ax.set_ylabel(f"radiance ({radiance.unit:latex_inline})");
../_images/spectrum_7_0.png

Line profiles#

A line is not infinitely narrow. The ions emitting it move, so each photon is Doppler shifted by the velocity of its ion along the line of sight, and the line is spread over a range of wavelengths about as wide as the motions of the plasma: the thermal motion of its ions, any unresolved nonthermal motion, and a shift of the whole profile if the plasma is flowing. synthesize() spreads each line into a Gaussian with these properties, in bins of velocity around its rest wavelength. The plasma at each temperature gets a profile of its own, as wide as the thermal motion at that temperature, so the hot and the cool parts of a column can move differently.

It takes the contribution functions of an ion as contribution_function() returns them, along with the mass of the ion, which sets its thermal width. Here is the \(\mathrm{O\,V}\) line from above, from plasma with no nonthermal motion and from plasma with a nonthermal width of \(20\,\mathrm{km\,s^{-1}}\).

ion = fiasco.Ion("O 5", temperature.ndarray)

contribution = utu.spectrum.contribution_function(
    ion=ion,
    density=density,
    axis_temperature=axis,
)

index = np.argmax(contribution.outputs.max(axis), axis="line")

width_nonthermal = na.ScalarArray([0, 20] * u.km / u.s, axes="width")

result = utu.spectrum.synthesize(
    contribution=contribution[index],
    temperature=temperature,
    mass=ion.mass,
    emission_measure=emission_measure,
    velocity=na.linspace(-150, 150, axis="velocity", num=121) * u.km / u.s,
    width_nonthermal=width_nonthermal,
    axis_temperature=axis,
)

fig, ax = plt.subplots(figsize=(7, 3.5), constrained_layout=True)
na.plt.stairs(
    result.inputs.wavelength,
    result.outputs,
    ax=ax,
    axis="velocity",
    label=width_nonthermal,
)
ax.set_xlabel(f"wavelength ({u.AA:latex_inline})")
ax.set_ylabel(f"spectral radiance\n({result.outputs.unit:latex_inline})")
ax.legend(title="nonthermal width");
../_images/spectrum_8_0.png

The radiance of each bin is the integral of the profile across it, so broadening a line moves its light around without adding to it, and loses none of it as long as the bins span the line. Summed over the bins, both profiles give the radiance of the line in the table above.

width_bin = np.diff(result.inputs.wavelength, axis="velocity")

(result.outputs * width_bin).sum("velocity")

ScalarArray(
    ndarray=[2652.82184418, 2652.82184157] erg / (s sr cm2),
    shape={'width': 2},
)