—

Boltzmann plot of the first electronic excited states of carbon.

The equilibrium is computed via Cantera at a given temperature and pressure.

From Boltzmann equation, we have:

\[n_\text{u} = n_0 \frac{g_\text{u} \exp\left(-\frac{E_\text{u}}{k_\text{B} T}\right)}{Z(T)}\]

Dividing by the total density, and rearranging terms, we obtain:

\[\log\left(\frac{x_\text{u}}{g_\text{u}}\right) = \log\left(\frac{x_\text{0}}{Z(T)}\right) -\frac{E_\text{u}}{k_\text{B} T}\]

Tags: kinetics equilibrium reverse reaction CH4 C2H2 Janev

# This is an option for the online documentation, so that each image is displayed separately.
# sphinx_gallery_multi_image = "single"

Import the required libraries.

import cantera as ct
import matplotlib.pyplot as plt
import numpy as np
from adjustText import adjust_text

import rizer.misc.units as u
from rizer.misc.ct_utils import load_goutier2025_mechanism
from rizer.misc.plt_utils import get_text, set_mpl_style

set_mpl_style(nb_columns=1)

Define mechanisms and parameters.

gas = load_goutier2025_mechanism("gas_with_electronic_excited_states")
assert "C(1D)" in gas.species_names
assert "H(n=2)" in gas.species_names

# Temperatures in Kelvin.
temperatures = np.array(
    [
        # 2000,
        # 5000,
        10000,
        15000,
        20000,
        30000,
        # 40000,
        # 50000,
    ]
)
# Pressure in Pascals.
P = ct.one_atm
# Initial mole fraction (here, we start with pure methane).
X_0 = "CH4:1"

# Energy and degenaracy of the first electronic excited states of carbon.
# cf. https://physics.nist.gov/cgi-bin/ASD/energy1.pl?de=0&spectrum=C+I&units=1&format=0&output=0&page_size=15&multiplet_ordered=0&conf_out=on&term_out=on&level_out=on&unc_out=1&j_out=on&g_out=on&lande_out=on&perc_out=on&biblio=on&temp=&submit=Retrieve+Data
energies_C_eV = np.array([0.0, 1.2637284, 2.6840136, 4.1826219])
gs_C = np.array([9, 5, 1, 5])
# Energy and degenaracy of the first electronic excited states of hydrogen.
# https://physics.nist.gov/cgi-bin/ASD/energy1.pl?de=0&spectrum=H+I&units=1&format=0&output=0&page_size=15&multiplet_ordered=0&conf_out=on&term_out=on&level_out=on&unc_out=1&j_out=on&g_out=on&lande_out=on&perc_out=on&biblio=on&temp=&submit=Retrieve+Data
energies_H_eV = np.array([0.0, 10.1988358, 12.0875052, 12.7485393])
gs_H = np.array([2, 8, 18, 32])

Plot the Boltzmann plot.

fig, ax = plt.subplots()

texts = []
for T in temperatures:
    # Equilibrate the mixture.
    gas.TPX = T, P, X_0
    gas.equilibrate("TP")
    # Get the mole fraction for carbon and hydrogen excited states.
    X_C = gas["C", "C(1D)", "C(1S)", "C(5So)"].X
    X_H = gas["H", "H(n=2)", "H(n=3)", "H(n=4)"].X

    # Plot log(x/g) vs energy for carbon excited states.
    points = ax.scatter(energies_C_eV, np.log(X_C / gs_C), marker="o", s=300)

    # Get the color of the last line (here it is a scatter plot with same color for all points).
    color = points.get_facecolor()[0]

    # Plot log(x/g) vs energy for hydrogen excited states.
    ax.scatter(energies_H_eV, np.log(X_H / gs_H), marker="s", s=300, color=color)

    # Fit data to `y=b-a*x` for carbon excited states.
    coefs_C = np.polynomial.polynomial.polyfit(
        x=energies_C_eV, y=np.log(X_C / gs_C), deg=1
    )
    ax.plot(
        energies_C_eV,
        np.polynomial.polynomial.polyval(energies_C_eV, coefs_C),
        ls="--",
        lw=4,
        color=color,
    )

    # Fit data to `y=b-a*x` for hydrogen excited states.
    coefs_H = np.polynomial.polynomial.polyfit(
        x=energies_H_eV, y=np.log(X_H / gs_H), deg=1
    )
    ax.plot(
        energies_H_eV,
        np.polynomial.polynomial.polyval(energies_H_eV, coefs_H),
        ls="-",
        lw=4,
        color=color,
    )

    # Get back the temperature from the fit for carbon excited states.
    # y = b - a*x
    a_per_eV = -coefs_C[1]  # eV^(-1)
    a_per_J = a_per_eV / (u.eV_to_J)  # J^(-1)
    # a = 1/(k_B * T)
    T_fit_C = 1 / (a_per_J * u.k_b)  # K
    print(f"Temperature from fit for carbon excited states: {T_fit_C:.2e} K")
    # Annotate at 2 eV.
    texts.append(
        get_text(
            2,
            float(np.polynomial.polynomial.polyval(2, coefs_C)),
            f"{round(int(T_fit_C), -1):,} K (C*)".replace(",", " "),
            ax=ax,
        )
    )

    # Get back the temperature from the fit for hydrogen excited states.
    # y = b - a*x
    a_per_eV = -coefs_H[1]  # eV^(-1)
    a_per_J = a_per_eV / (u.eV_to_J)  # J^(-1)
    # a = 1/(k_B * T)
    T_fit_H = 1 / (a_per_J * u.k_b)  # K
    print(f"Temperature from fit for hydrogen excited states: {T_fit_H:.2e} K")
    # Annotate at 7 eV.
    texts.append(
        get_text(
            7,
            float(np.polynomial.polynomial.polyval(7, coefs_H)),
            f"{round(int(T_fit_H), -1):,} K (H*)".replace(",", " "),
            ax=ax,
        )
    )

ax.set_xlabel(r"$\varepsilon$ [eV]")
ax.set_ylabel(r"$\log\left(\frac{x}{g}\right)$")
ax.set_title("Boltzmann plot of the first 4 electronic excited states of C and H.")
adjust_text(texts, avoid_self=False)
plt.show()

Temperature from fit for carbon excited states: 1.00e+04 K
Temperature from fit for hydrogen excited states: 1.00e+04 K
Temperature from fit for carbon excited states: 1.50e+04 K
Temperature from fit for hydrogen excited states: 1.50e+04 K
Temperature from fit for carbon excited states: 2.00e+04 K
Temperature from fit for hydrogen excited states: 2.00e+04 K
Temperature from fit for carbon excited states: 3.00e+04 K
Temperature from fit for hydrogen excited states: 3.00e+04 K