Part II: Solar Atmosphere | Chapter 6

The Chromosphere

Temperature minimum, chromospheric heating, spicules, and the Saha equation

6.1 The Temperature Minimum Region

Derivation 1: Radiative Equilibrium Temperature Minimum

Above the photosphere, the temperature initially decreases to a minimum of about 4400 K at a height of ~500 km before rising dramatically into the chromosphere and corona. This minimum marks the boundary between radiative cooling and non-radiative heating.

Step 1. In radiative equilibrium, the temperature at low optical depths follows:

$$T^4(\tau) = \frac{1}{2} T_{\text{eff}}^4 \left(1 + \frac{3}{2}\tau\right)$$

Step 2. As $\tau \to 0$, $T \to T_{\text{eff}} / 2^{1/4} \approx 4860$ K. The observed minimum of 4400 K is even lower because the atmosphere is not in strict radiative equilibrium at these heights—CO molecular cooling plays an important role.

Step 3. Above the temperature minimum, some non-radiative heating mechanism must deposit energy to raise the temperature. The chromospheric energy budget requires:

$$\boxed{F_{\text{heat}} \approx 4 \times 10^3 \text{ W/m}^2 \text{ (quiet Sun)} \sim 2 \times 10^4 \text{ W/m}^2 \text{ (active regions)}}$$

This is only about $10^{-4}$ of the total solar flux, but maintaining the chromospheric temperature against strong radiative losses (mainly in Ca II, Mg II, and H$\alpha$) requires a continuous supply of mechanical or magnetic energy.

6.2 The Saha Equation for Ionization

Derivation 2: Saha Ionization Equilibrium

The ionization state of the chromospheric plasma determines the opacity, emission, and energy balance. The Saha equation gives the ionization fraction in thermal equilibrium.

Step 1. Consider the ionization reaction $A \rightleftharpoons A^+ + e^-$. In chemical equilibrium, the chemical potentials balance:$\mu_A = \mu_{A^+} + \mu_e$.

Step 2. Using the quantum statistical mechanics result for ideal gases with$\mu = -k_BT \ln(Z/N)$ and the thermal de Broglie wavelength$\Lambda = h/\sqrt{2\pi m k_B T}$:

$$\frac{n_{i+1} n_e}{n_i} = \frac{2 g_{i+1}}{g_i} \frac{1}{\Lambda_e^3} \exp\left(-\frac{\chi_i}{k_B T}\right)$$

Step 3. Substituting the electron thermal wavelength:

$$\boxed{\frac{n_{i+1} n_e}{n_i} = \frac{2 g_{i+1}}{g_i} \left(\frac{2\pi m_e k_B T}{h^2}\right)^{3/2} \exp\left(-\frac{\chi_i}{k_B T}\right)}$$

Step 4. For hydrogen ($\chi = 13.6$ eV), at the temperature minimum (4400 K): the ionization fraction is $x \sim 10^{-4}$. By 8000 K (upper chromosphere):$x \sim 0.5$. The sharp transition creates a dramatic change in opacity.

The partial ionization of hydrogen in the chromosphere is responsible for the extreme non-LTE effects that make chromospheric spectral line formation notoriously difficult to model.

6.3 Chromospheric Heating Mechanisms

Derivation 3: Acoustic Wave Dissipation

Acoustic waves generated by turbulent convection steepen into shocks as they propagate upward into the decreasing-density chromosphere.

Step 1. An acoustic wave with velocity amplitude $v$ carries an energy flux:

$$F_{\text{ac}} = \frac{1}{2} \rho v^2 c_s$$

Step 2. As the wave propagates upward, conservation of energy flux requires the velocity amplitude to grow: $v \propto \rho^{-1/2}$. The wave steepens into a shock when $v \sim c_s$, which occurs at a height:

$$h_{\text{shock}} \sim H_P \ln\left(\frac{c_s}{v_0}\right) \approx 2H_P \approx 300 \text{ km}$$

Step 3. The dissipation rate per unit volume for a weak shock of Mach number$M = v/c_s$ is:

$$\boxed{Q = \frac{\rho c_s^3}{12 H_P} (M - 1)^3 \quad \text{(weak shock heating)}}$$

Acoustic wave heating can account for the quiet-Sun chromospheric energy requirements, but is insufficient for active-region chromospheres. Magnetic mechanisms (Alfven wave dissipation, magnetic reconnection, current dissipation) are needed in magnetically dominated regions.

6.4 Chromospheric Diagnostics: H$\alpha$ and Ca II

Derivation 4: H$\alpha$ Line Formation

Step 1. H$\alpha$ (656.3 nm) is the $n=3 \to n=2$transition of hydrogen. The opacity depends on the population of the $n=2$ level, which in LTE is given by the Boltzmann equation:

$$\frac{n_2}{n_1} = \frac{g_2}{g_1} \exp\left(-\frac{E_{21}}{k_BT}\right) = 4 \exp\left(-\frac{10.2 \text{ eV}}{k_BT}\right)$$

Step 2. At T = 6000 K, $n_2/n_1 \sim 5 \times 10^{-9}$—only a tiny fraction of hydrogen atoms are in the $n=2$ state. But the H$\alpha$opacity is proportional to $n_2$, and the column density is enormous, making the line optically thick.

Step 3. In the chromosphere, H$\alpha$ is strongly affected by non-LTE effects. The source function departs from the Planck function because:

$$S_\nu = \frac{(1-\epsilon)J_\nu + \epsilon B_\nu}{1} \neq B_\nu$$

where $\epsilon$ is the photon destruction probability and $J_\nu$ is the mean intensity. In H$\alpha$, $\epsilon \ll 1$, meaning the source function is dominated by scattering rather than thermal emission. Ca II H&K lines (393.4, 396.8 nm) are the strongest chromospheric lines, with emission reversals in their cores indicating the temperature rise.

6.5 Spicules and Chromospheric Dynamics

Derivation 5: Spicule Dynamics as Rebound Shocks

Spicules are jet-like features that extend from the chromosphere into the corona, reaching heights of 5,000-10,000 km with velocities of 20-100 km/s. About $10^5$ spicules are present on the Sun at any time.

Step 1. A simple model treats spicules as gas driven upward by a pressure pulse at the base. The equation of motion for a column of chromospheric gas:

$$\rho \frac{dv}{dt} = -\frac{dP}{dz} - \rho g$$

Step 2. For a pressure pulse $\Delta P$ driving material to height$h$ against gravity, energy conservation gives:

$$h = \frac{v_0^2}{2g} = \frac{(\Delta P / \rho)}{2g}$$

Step 3. For $v_0 \approx 25$ km/s and $g = 274$ m/s$^2$:

$$\boxed{h = \frac{(25 \times 10^3)^2}{2 \times 274} \approx 1140 \text{ km}}$$

This underestimates the observed heights, suggesting that additional driving forces (magnetic tension, Alfven waves) play a role. Type II spicules, discovered by Hinode, are faster (50-100 km/s), shorter-lived (~10-150 s), and may contribute to coronal heating by injecting hot plasma ($\sim 10^5$ K) into the corona.

Numerical Simulation

Chromosphere: Temperature Profile, Saha Ionization, Wave Steepening, Spicule Dynamics

Python

script.py135 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle('Chromospheric Physics', fontsize=16, fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

for ax in axes.flat:
    ax.set_facecolor('#0a0a0a')
    ax.tick_params(colors='white')
    ax.xaxis.label.set_color('white')
    ax.yaxis.label.set_color('white')
    ax.title.set_color('white')
    for spine in ax.spines.values():
        spine.set_color('#333')

# --- Panel 1: VAL-C Temperature profile ---
ax1 = axes[0, 0]
h = np.linspace(-100, 2500, 500)  # height in km

# VAL-C model (approximate)
T = np.where(h < 0, 6400 + 8 * h,
    np.where(h < 500, 6400 - 4.0 * h,
    np.where(h < 900, 4400,
    np.where(h < 2000, 4400 + 4.5 * (h - 900),
    4400 + 4.5 * 1100 + 100 * (h - 2000)))))

ax1.plot(h, T, color='#f43f5e', linewidth=2.5)
ax1.axhline(y=4400, color='#fbbf24', linestyle='--', alpha=0.5, label='T_min = 4400 K')
ax1.axvspan(-100, 0, alpha=0.1, color='orange', label='Photosphere')
ax1.axvspan(0, 500, alpha=0.1, color='red')
ax1.axvspan(500, 2200, alpha=0.1, color='magenta', label='Chromosphere')

ax1.set_xlabel('Height [km]')
ax1.set_ylabel('Temperature [K]')
ax1.set_title('VAL-C Temperature Profile')
ax1.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: Saha ionization fraction ---
ax2 = axes[0, 1]
T_saha = np.linspace(3000, 15000, 500)
k_B = 8.617e-5  # eV/K
chi = 13.6  # eV
m_e = 9.109e-31
h_planck = 6.626e-34
k_B_SI = 1.381e-23

# Saha equation for hydrogen (assume n_e from self-consistent solution)
# n_e * n_p / n_H = Saha_coeff
# For fixed total n_H + n_p, solve quadratic
n_total = 1e17  # cm^-3 (typical chromospheric)

Saha_coeff = 2.0 * (2 * np.pi * m_e * k_B_SI * T_saha / h_planck**2)**1.5 * np.exp(-chi / (k_B * T_saha))
Saha_coeff = Saha_coeff * 1e-6  # convert to cm^-3

# Solve n_e^2 / (n_total - n_e) = Saha_coeff
# n_e^2 + Saha * n_e - Saha * n_total = 0
x_ion = (-Saha_coeff + np.sqrt(Saha_coeff**2 + 4 * Saha_coeff * n_total)) / (2 * n_total)
x_ion = np.clip(x_ion, 1e-10, 1.0)

ax2.semilogy(T_saha, x_ion, color='#fb923c', linewidth=2.5)
ax2.axvline(x=4400, color='#fbbf24', linestyle='--', alpha=0.5, label='T_min')
ax2.axvline(x=6000, color='#f43f5e', linestyle='--', alpha=0.5, label='Photosphere')
ax2.axvline(x=10000, color='#a78bfa', linestyle='--', alpha=0.5, label='Upper chrom.')
ax2.set_xlabel('Temperature [K]')
ax2.set_ylabel('Ionization fraction x')
ax2.set_title('Hydrogen Ionization (Saha Eq.)')
ax2.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: Acoustic shock steepening ---
ax3 = axes[1, 0]
z = np.linspace(0, 2000, 300)  # km
H_P = 150  # km

# Wave amplitude growth
v0 = 1.0  # km/s at z=0
v_amp = v0 * np.exp(z / (2 * H_P))
c_s = 7.0  # km/s

# Clamp at shock formation
shock_height = 2 * H_P * np.log(c_s / v0)
v_amp_clipped = np.where(z < shock_height, v_amp, c_s + 0.5 * (z - shock_height) / H_P)

ax3.plot(z, v_amp, color='#f43f5e', linewidth=2.5, label='Wave amplitude v(z)')
ax3.axhline(y=c_s, color='#60a5fa', linestyle='--', linewidth=1.5, label='Sound speed c_s')
ax3.axvline(x=shock_height, color='#fbbf24', linestyle=':', linewidth=1.5, label='Shock formation')
ax3.fill_between(z, 0, np.minimum(v_amp, c_s), alpha=0.15, color='#f43f5e')
ax3.set_xlabel('Height [km]')
ax3.set_ylabel('Velocity amplitude [km/s]')
ax3.set_title('Acoustic Wave Steepening')
ax3.set_ylim(0, 15)
ax3.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 4: Spicule trajectory ---
ax4 = axes[1, 1]
g = 0.274  # km/s^2
t = np.linspace(0, 600, 300)  # seconds

for v0_sp, col, lbl in [(25, '#f43f5e', 'Type I (25 km/s)'),
                          (50, '#fb923c', 'Type II (50 km/s)'),
                          (100, '#fbbf24', 'Type II fast (100 km/s)')]:
    h_sp = v0_sp * t - 0.5 * g * t**2
    h_sp = np.maximum(h_sp, 0)
    ax4.plot(t, h_sp / 1000, color=col, linewidth=2.5, label=lbl)

ax4.axhline(y=10, color='white', linestyle=':', alpha=0.3, label='Typical spicule height')
ax4.set_xlabel('Time [s]')
ax4.set_ylabel('Height [Mm]')
ax4.set_title('Spicule Ballistic Trajectories')
ax4.set_ylim(0, 25)
ax4.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("Plot saved to output.png")

print("\n" + "=" * 60)
print("CHROMOSPHERE")
print("=" * 60)
print("\n--- Temperature Structure ---")
print("  T_min ~ 4400 K at h ~ 500 km")
print("  Upper chromosphere T ~ 10,000-20,000 K")
print("  Transition region: T rises from 20,000 to 10^6 K over ~100 km")
print("\n--- Heating Requirements ---")
print("  Quiet Sun: ~4 x 10^3 W/m^2")
print("  Active region: ~2 x 10^4 W/m^2")
print("\n--- Spicules ---")
print("  Type I: v ~ 25 km/s, h ~ 5 Mm, lifetime ~ 5 min")
print("  Type II: v ~ 50-100 km/s, h ~ 5-10 Mm, lifetime ~ 10-150 s")
print("  Number on Sun: ~10^5 at any time")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6.6 The Saha Equation — Full Statistical-Mechanical Derivation

We derive the Saha equation from the grand canonical ensemble, connecting microscopic quantum mechanics to macroscopic ionization balance.

Derivation from Partition Functions

Step 1. Consider the ionization reaction$X_i \rightleftharpoons X_{i+1} + e^-$. In thermal equilibrium the chemical potentials satisfy $\mu_i = \mu_{i+1} + \mu_e$.

Step 2. For an ideal gas at temperature $T$, the chemical potential of species $s$ is:

$$\mu_s = -k_B T \ln\!\left(\frac{g_s}{\lambda_s^3 \, n_s}\right)$$

where $g_s$ is the internal partition function (statistical weight),$\lambda_s = h / \sqrt{2\pi m_s k_B T}$ is the thermal de Broglie wavelength, and $n_s$ is the number density.

Step 3. Substituting the chemical potential balance and noting that the ionization energy $\chi_i$ shifts the zero of the internal energy of the$(i+1)$-th ion:

$$\ln\!\left(\frac{g_i}{\lambda_i^3 n_i}\right) = \ln\!\left(\frac{g_{i+1}}{\lambda_{i+1}^3 n_{i+1}}\right) + \ln\!\left(\frac{2}{\lambda_e^3 n_e}\right) - \frac{\chi_i}{k_B T}$$

The factor of 2 for the electron comes from its two spin states ($g_e = 2$).

Step 4. Since $m_i \approx m_{i+1}$ (the electron mass is negligible compared to the ion mass), $\lambda_i \approx \lambda_{i+1}$, and these cancel. Exponentiating:

$$\boxed{\frac{n_{i+1}\,n_e}{n_i} = \frac{2\,g_{i+1}}{g_i}\,\frac{1}{\lambda_e^3}\,\exp\!\left(-\frac{\chi_i}{k_BT}\right) = \frac{2\,g_{i+1}}{g_i}\left(\frac{2\pi m_e k_BT}{h^2}\right)^{\!3/2}\exp\!\left(-\frac{\chi_i}{k_BT}\right)}$$

The Saha ionization equation

Application to Chromospheric Hydrogen

Step 5. For hydrogen ($\chi = 13.6$ eV, $g_1 = 2$for the ground state, $g_2 = 1$ for the proton), in a pure hydrogen plasma with$n_e = n_p$ and total number density $n_{\text{tot}} = n_H + n_p$, the ionization fraction $x = n_p / n_{\text{tot}}$ satisfies:

$$\frac{x^2}{1 - x}\,n_{\text{tot}} = \Phi(T), \qquad \Phi(T) = \left(\frac{2\pi m_e k_BT}{h^2}\right)^{3/2}\exp\!\left(-\frac{13.6\text{ eV}}{k_BT}\right)$$

Step 6. This is a quadratic in $x$:$x^2 + (\Phi / n_{\text{tot}})x - \Phi / n_{\text{tot}} = 0$, with the physical root:

$$x = \frac{-\Phi + \sqrt{\Phi^2 + 4\Phi\,n_{\text{tot}}}}{2\,n_{\text{tot}}}$$

At the temperature minimum ($T = 4400$ K, $n_{\text{tot}} \sim 10^{17}$ cm$^{-3}$):$x \sim 10^{-4}$ (almost neutral). At $T = 10{,}000$ K:$x \sim 0.5$. This dramatic transition drives the extreme non-LTE effects of chromospheric spectral line formation.

6.7 Radiative Equilibrium Temperature Minimum — Detailed Derivation

Step 1. In a grey atmosphere, the Eddington solution for the mean intensity is:

$$J(\tau) = \frac{3F}{4\pi}\left(\tau + \frac{2}{3}\right)$$

where $F = \sigma T_{\text{eff}}^4$ is the bolometric flux. In radiative equilibrium$S = J$, so $B(\tau) = J(\tau)$ and$\sigma T^4 / \pi = J$.

Step 2. As we move outward to the surface ($\tau \to 0$):

$$T^4(\tau \to 0) = \frac{3}{4}T_{\text{eff}}^4 \cdot \frac{2}{3} = \frac{1}{2}T_{\text{eff}}^4$$

Step 3. Therefore the minimum temperature in a grey radiative equilibrium atmosphere is:

$$\boxed{T_{\min} = \frac{T_{\text{eff}}}{2^{1/4}} = \frac{5778}{1.189} \approx 4860 \text{ K}}$$

Grey atmosphere temperature minimum (Milne-Eddington)

Step 4. The observed minimum is even lower ($\approx 4400$ K) because:

• Non-grey effects: UV opacity is much larger than IR opacity, so UV photons escape from shallower (cooler) layers
• CO molecular cooling: CO ro-vibrational bands near 4.7 $\mu$m are efficient radiators, cooling the gas below the grey prediction
• Line blanketing: thousands of spectral lines redistribute energy from UV to longer wavelengths

Above the temperature minimum, radiative equilibrium breaks down entirely. Non-radiative heating (acoustic waves, magnetic dissipation) must supply energy to produce the observed chromospheric temperature rise.

Diagram: Chromospheric Temperature Profile

Advanced Simulation: Saha Ionization and Chromospheric Emission

Saha Ionization Fraction vs Temperature and Chromospheric Emission Spectrum

Python

script.py125 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
fig.suptitle('Saha Ionization & Chromospheric Emission Spectrum',
             fontsize=15, fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

# === Panel 1: Hydrogen ionization fraction vs T for different densities ===
ax1 = axes[0]
T_arr = np.linspace(3000, 25000, 800)
k_B_eV = 8.617e-5    # eV / K
chi_H = 13.6          # eV
m_e = 9.109e-31       # kg
h_p = 6.626e-34       # J s
k_B_SI = 1.381e-23    # J / K

densities = [1e15, 1e16, 1e17, 1e18]
colors_d = ['#60a5fa', '#4ade80', '#f43f5e', '#fbbf24']
labels_d = ['n_tot = 10^15 cm^-3', 'n_tot = 10^16', 'n_tot = 10^17', 'n_tot = 10^18']

for n_tot_cgs, col, lbl in zip(densities, colors_d, labels_d):
    n_tot = n_tot_cgs * 1e6   # convert to m^-3
    Phi = (2.0 * np.pi * m_e * k_B_SI * T_arr / h_p**2)**1.5 * np.exp(-chi_H / (k_B_eV * T_arr))
    # Phi is in m^-3
    ratio = Phi / n_tot
    x_ion = (-ratio + np.sqrt(ratio**2 + 4.0 * ratio)) / 2.0
    x_ion = np.clip(x_ion, 1e-12, 1.0)
    ax1.semilogy(T_arr, x_ion, color=col, linewidth=2, label=lbl)

ax1.axvline(x=4400, color='#fbbf24', linestyle='--', alpha=0.4, linewidth=1)
ax1.text(4500, 1e-1, 'T_min', color='#fbbf24', fontsize=9)
ax1.axvline(x=10000, color='#f43f5e', linestyle='--', alpha=0.4, linewidth=1)
ax1.text(10100, 1e-1, 'Upper chrom.', color='#f43f5e', fontsize=9)
ax1.axhline(y=0.5, color='white', linestyle=':', alpha=0.2)
ax1.set_xlabel('Temperature [K]')
ax1.set_ylabel('Ionization fraction x = n_p / n_tot')
ax1.set_title('Hydrogen Ionization (Saha Equation)')
ax1.set_ylim(1e-10, 2)
ax1.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white', loc='lower right')

# === Panel 2: Schematic chromospheric emission spectrum ===
ax2 = axes[1]
wavelength = np.linspace(380, 700, 2000)  # nm

# Photospheric continuum (approximate blackbody shape at 5778 K)
h_c = 6.626e-34
c_light = 3e8
k_Bol = 1.381e-23
lam_m = wavelength * 1e-9
B_phot = 2 * h_c * c_light**2 / lam_m**5 / (np.exp(h_c * c_light / (lam_m * k_Bol * 5778)) - 1)
B_phot = B_phot / B_phot.max()

# Chromospheric emission lines (Gaussian profiles on top)
def emission_line(lam, lam0, width, strength):
    return strength * np.exp(-0.5 * ((lam - lam0) / width)**2)

chrom = np.zeros_like(wavelength)
# H-alpha 656.3 nm
chrom += emission_line(wavelength, 656.3, 0.5, 0.8)
# H-beta 486.1 nm
chrom += emission_line(wavelength, 486.1, 0.4, 0.35)
# H-gamma 434.0 nm
chrom += emission_line(wavelength, 434.0, 0.35, 0.2)
# Ca II K 393.4 nm
chrom += emission_line(wavelength, 393.4, 0.6, 0.95)
# Ca II H 396.8 nm
chrom += emission_line(wavelength, 396.8, 0.5, 0.7)

ax2.fill_between(wavelength, 0, B_phot * 0.3, color='#fb923c', alpha=0.15, label='Photospheric continuum')
ax2.plot(wavelength, chrom, color='#f43f5e', linewidth=2, label='Chromospheric emission')

# Label the lines
ax2.annotate('H-alpha\n656.3 nm', xy=(656.3, 0.82), xytext=(620, 0.95),
            fontsize=8, color='#fbbf24',
            arrowprops=dict(arrowstyle='->', color='#fbbf24', lw=0.8))
ax2.annotate('Ca II K\n393.4 nm', xy=(393.4, 0.97), xytext=(410, 1.05),
            fontsize=8, color='#4ade80',
            arrowprops=dict(arrowstyle='->', color='#4ade80', lw=0.8))
ax2.annotate('Ca II H\n396.8 nm', xy=(396.8, 0.72), xytext=(430, 0.78),
            fontsize=8, color='#4ade80',
            arrowprops=dict(arrowstyle='->', color='#4ade80', lw=0.8))
ax2.annotate('H-beta', xy=(486.1, 0.37), xytext=(500, 0.5),
            fontsize=8, color='#a78bfa',
            arrowprops=dict(arrowstyle='->', color='#a78bfa', lw=0.8))

ax2.set_xlabel('Wavelength [nm]')
ax2.set_ylabel('Relative intensity')
ax2.set_title('Chromospheric Emission Spectrum (schematic)')
ax2.set_ylim(0, 1.15)
ax2.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("Plot saved to output.png")
print("")
print("=== Saha Equation Results (n_tot = 10^17 cm^-3) ===")

T_check = [4400, 6000, 8000, 10000, 15000, 20000]
n_tot_check = 1e17 * 1e6
for Tc in T_check:
    Phi_c = (2.0 * np.pi * m_e * k_B_SI * Tc / h_p**2)**1.5 * np.exp(-chi_H / (k_B_eV * Tc))
    ratio_c = Phi_c / n_tot_check
    xc = (-ratio_c + np.sqrt(ratio_c**2 + 4.0 * ratio_c)) / 2.0
    print("  T = " + str(Tc) + " K:  x = " + str("%.2e" % xc))

print("")
print("=== Key Chromospheric Emission Lines ===")
print("  Ca II K  393.4 nm  (strongest chromospheric line)")
print("  Ca II H  396.8 nm")
print("  H-alpha  656.3 nm  (Balmer series)")
print("  H-beta   486.1 nm")