plasma_parameters.tex

\documentclass[11pt,a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amssymb}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{siunitx}
\usepackage{geometry}
\geometry{margin=1in}
\usepackage{pythontex}
\usepackage{hyperref}
\usepackage{float}

\title{Fundamental Plasma Parameters:\\From Debye Shielding to Wave Propagation}
\author{Computational Plasma Physics Laboratory}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
We present a comprehensive computational analysis of fundamental plasma parameters spanning laboratory, space, and fusion plasmas. Starting from the Debye length and plasma frequency, we derive characteristic length and time scales that govern collective behavior. We compute the plasma parameter $N_D$ to verify the plasma approximation, analyze Coulomb collision frequencies and Spitzer conductivity, and examine magnetized plasma dynamics through Larmor radii and cyclotron frequencies. Wave propagation characteristics including Langmuir, ion acoustic, and Alfvén waves are calculated across parameter regimes spanning electron densities from $10^{6}$ to $10^{26}$ m$^{-3}$ and temperatures from 0.1 eV to 10 keV. Our results demonstrate that the dimensionless ratios of these fundamental scales determine whether plasmas are collisional or collisionless, magnetized or unmagnetized, and which wave modes can propagate.
\end{abstract}

\section{Introduction}

Plasma behavior is governed by a hierarchy of characteristic length and time scales that emerge from the interplay between electromagnetic forces, thermal motion, and magnetic confinement \cite{chen1984introduction, goldston1995introduction}. The most fundamental of these is the Debye length $\lambda_D$, which characterizes the distance over which electrostatic perturbations are shielded by collective electron redistribution \cite{debye1923theory}. When the number of particles in a Debye sphere $N_D = n \lambda_D^3$ greatly exceeds unity, collective effects dominate individual particle interactions, and the system exhibits true plasma behavior \cite{spitzer1956physics}.

The plasma frequency $\omega_{pe}$ sets the fundamental time scale for electrostatic oscillations, arising when electrons are displaced from a uniform ion background \cite{tonks1929oscillations}. The ratio of collision frequency to plasma frequency determines whether the plasma is collisional (fluid-like) or collisionless (kinetic), fundamentally altering the physics of transport and wave propagation \cite{krall1973principles}.

In magnetized plasmas, charged particles gyrate around magnetic field lines with cyclotron frequency $\omega_c$ and Larmor radius $r_L$, introducing anisotropy and enabling new wave modes including electromagnetic Alfvén waves that transport energy along field lines at the Alfvén velocity \cite{alfven1942existence}. The magnetization parameter $\omega_c/\nu_{coll}$ determines whether particles complete many gyro-orbits between collisions, a critical distinction for confinement in fusion devices \cite{wesson2011tokamaks}.

This report computes these fundamental parameters across the full range of natural and laboratory plasmas, from the tenuous solar wind to dense fusion cores, demonstrating how dimensionless ratios constructed from these scales predict plasma regime transitions and govern collective phenomena.

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
from scipy import constants
from scipy.optimize import fsolve

# Configure matplotlib for publication-quality plots
plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.size'] = 11

# Physical constants
epsilon_0 = constants.epsilon_0  # Permittivity of free space [F/m]
e = constants.e  # Elementary charge [C]
m_e = constants.m_e  # Electron mass [kg]
m_p = constants.m_p  # Proton mass [kg]
k_B = constants.k  # Boltzmann constant [J/K]
c = constants.c  # Speed of light [m/s]
mu_0 = constants.mu_0  # Permeability of free space [H/m]

# Function definitions for plasma parameters
def debye_length(n_e, T_e):
    """
    Compute electron Debye length.

    Parameters:
    -----------
    n_e : float or array
        Electron density [m^-3]
    T_e : float or array
        Electron temperature [eV]

    Returns:
    --------
    lambda_D : float or array
        Debye length [m]
    """
    T_e_joules = T_e * e  # Convert eV to Joules
    return np.sqrt(epsilon_0 * T_e_joules / (n_e * e**2))

def plasma_frequency(n_e):
    """
    Compute electron plasma frequency.

    Parameters:
    -----------
    n_e : float or array
        Electron density [m^-3]

    Returns:
    --------
    omega_pe : float or array
        Plasma frequency [rad/s]
    """
    return np.sqrt(n_e * e**2 / (epsilon_0 * m_e))

def plasma_parameter(n_e, T_e):
    """
    Compute number of particles in Debye sphere.

    Parameters:
    -----------
    n_e : float or array
        Electron density [m^-3]
    T_e : float or array
        Electron temperature [eV]

    Returns:
    --------
    N_D : float or array
        Number of particles in Debye sphere (dimensionless)
    """
    lambda_D = debye_length(n_e, T_e)
    return (4.0/3.0) * np.pi * n_e * lambda_D**3

def coulomb_logarithm(n_e, T_e):
    """
    Compute Coulomb logarithm for electron-ion collisions.

    Parameters:
    -----------
    n_e : float or array
        Electron density [m^-3]
    T_e : float or array
        Electron temperature [eV]

    Returns:
    --------
    ln_Lambda : float or array
        Coulomb logarithm (dimensionless)
    """
    # Approximate formula valid for most plasmas
    return 23.0 - 0.5 * np.log(n_e / 1e6) + 1.25 * np.log(T_e)

def collision_frequency(n_e, T_e):
    """
    Compute electron-ion collision frequency (90-degree deflection).

    Parameters:
    -----------
    n_e : float or array
        Electron density [m^-3]
    T_e : float or array
        Electron temperature [eV]

    Returns:
    --------
    nu_ei : float or array
        Collision frequency [Hz]
    """
    T_e_joules = T_e * e
    ln_Lambda = coulomb_logarithm(n_e, T_e)
    # Spitzer collision frequency
    return (n_e * e**4 * ln_Lambda) / (12 * np.pi**1.5 * epsilon_0**2 * m_e**0.5 * T_e_joules**1.5)

def mean_free_path(n_e, T_e):
    """
    Compute electron-ion mean free path.

    Parameters:
    -----------
    n_e : float or array
        Electron density [m^-3]
    T_e : float or array
        Electron temperature [eV]

    Returns:
    --------
    lambda_mfp : float or array
        Mean free path [m]
    """
    v_thermal = np.sqrt(2 * T_e * e / m_e)
    nu_ei = collision_frequency(n_e, T_e)
    return v_thermal / nu_ei

def spitzer_conductivity(T_e):
    """
    Compute Spitzer electrical conductivity.

    Parameters:
    -----------
    T_e : float or array
        Electron temperature [eV]

    Returns:
    --------
    sigma : float or array
        Electrical conductivity [S/m]
    """
    T_e_joules = T_e * e
    # Spitzer conductivity (SI units)
    return 1.96e4 * (T_e)**1.5  # Simplified form in S/m

def cyclotron_frequency(B, m):
    """
    Compute cyclotron frequency.

    Parameters:
    -----------
    B : float or array
        Magnetic field strength [T]
    m : float
        Particle mass [kg]

    Returns:
    --------
    omega_c : float or array
        Cyclotron frequency [rad/s]
    """
    return e * B / m

def larmor_radius(v_perp, B, m):
    """
    Compute Larmor radius.

    Parameters:
    -----------
    v_perp : float or array
        Perpendicular velocity [m/s]
    B : float or array
        Magnetic field strength [T]
    m : float
        Particle mass [kg]

    Returns:
    --------
    r_L : float or array
        Larmor radius [m]
    """
    omega_c = cyclotron_frequency(B, m)
    return v_perp / omega_c

def thermal_larmor_radius(T, B, m):
    """
    Compute thermal Larmor radius.

    Parameters:
    -----------
    T : float or array
        Temperature [eV]
    B : float or array
        Magnetic field strength [T]
    m : float
        Particle mass [kg]

    Returns:
    --------
    r_L : float or array
        Thermal Larmor radius [m]
    """
    v_thermal = np.sqrt(2 * T * e / m)
    return larmor_radius(v_thermal, B, m)

def alfven_velocity(B, n):
    """
    Compute Alfvén velocity.

    Parameters:
    -----------
    B : float or array
        Magnetic field strength [T]
    n : float or array
        Ion density [m^-3]

    Returns:
    --------
    v_A : float or array
        Alfvén velocity [m/s]
    """
    return B / np.sqrt(mu_0 * n * m_p)

def ion_acoustic_velocity(T_e, T_i, gamma_e=1, gamma_i=3):
    """
    Compute ion acoustic velocity.

    Parameters:
    -----------
    T_e : float or array
        Electron temperature [eV]
    T_i : float or array
        Ion temperature [eV]
    gamma_e : float
        Electron adiabatic index (1 for isothermal, 3 for adiabatic)
    gamma_i : float
        Ion adiabatic index

    Returns:
    --------
    c_s : float or array
        Ion acoustic velocity [m/s]
    """
    return np.sqrt((gamma_e * T_e * e + gamma_i * T_i * e) / m_p)

\end{pycode}

\section{Debye Shielding and Plasma Frequency}

The Debye length arises from a balance between electrostatic potential energy and thermal kinetic energy. When a test charge is introduced into a plasma, electrons redistribute to screen the perturbation over a characteristic distance:
\begin{equation}
\lambda_D = \sqrt{\frac{\epsilon_0 k_B T_e}{n_e e^2}} = 743 \sqrt{\frac{T_e[\text{eV}]}{n_e[\text{m}^{-3}]}} \, \text{m}
\end{equation}

The electron plasma frequency represents collective oscillations when charge displacement creates a restoring electric field:
\begin{equation}
\omega_{pe} = \sqrt{\frac{n_e e^2}{\epsilon_0 m_e}} = 56.4 \sqrt{n_e[\text{m}^{-3}]} \, \text{rad/s}
\end{equation}

For a system to exhibit collective plasma behavior, the number of particles in a Debye sphere must greatly exceed unity:
\begin{equation}
N_D = n_e \lambda_D^3 = \frac{4\pi}{3} n_e \left(\frac{\epsilon_0 k_B T_e}{n_e e^2}\right)^{3/2} \gg 1
\end{equation}

\begin{pycode}
# Parameter space for different plasma regimes
n_e_range = np.logspace(6, 26, 200)  # 10^6 to 10^26 m^-3

# Representative temperatures
temperatures = {
    'Solar Corona': 100,  # eV
    'Tokamak Edge': 10,   # eV
    'Tokamak Core': 5000, # eV
    'Ionosphere': 0.1     # eV
}

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 12))

# Plot 1: Debye length vs density for different temperatures
for label, T_e in temperatures.items():
    lambda_D = debye_length(n_e_range, T_e)
    ax1.loglog(n_e_range, lambda_D * 1e6, linewidth=2.5, label=f'{label} ({T_e} eV)')

# Mark specific plasma regimes
regimes = {
    'Ionosphere': (1e11, 0.1),
    'Solar Wind': (1e7, 10),
    'Solar Corona': (1e14, 100),
    'Glow Discharge': (1e16, 2),
    'Tokamak Edge': (1e19, 10),
    'Tokamak Core': (1e20, 5000),
    'ICF Compression': (1e26, 1000)
}

for regime, (n, T) in regimes.items():
    ld = debye_length(n, T) * 1e6
    ax1.plot(n, ld, 'o', markersize=8, color='black')
    ax1.annotate(regime, xy=(n, ld), xytext=(10, 10),
                textcoords='offset points', fontsize=8,
                bbox=dict(boxstyle='round,pad=0.3', facecolor='yellow', alpha=0.5))

ax1.set_xlabel(r'Electron Density $n_e$ [m$^{-3}$]', fontsize=12)
ax1.set_ylabel(r'Debye Length $\lambda_D$ [$\mu$m]', fontsize=12)
ax1.set_title(r'Debye Shielding Length Across Plasma Regimes', fontsize=13, fontweight='bold')
ax1.legend(loc='upper right', fontsize=9)
ax1.grid(True, alpha=0.3, which='both')
ax1.set_xlim([1e6, 1e27])
ax1.set_ylim([1e-6, 1e6])

# Plot 2: Plasma frequency and period
omega_pe = plasma_frequency(n_e_range)
f_pe = omega_pe / (2 * np.pi)
tau_pe = 1 / f_pe

ax2.loglog(n_e_range, f_pe, 'b-', linewidth=2.5, label=r'$f_{pe}$ [Hz]')
ax2_twin = ax2.twinx()
ax2_twin.loglog(n_e_range, tau_pe * 1e12, 'r--', linewidth=2.5, label=r'$\tau_{pe}$ [ps]')

ax2.set_xlabel(r'Electron Density $n_e$ [m$^{-3}$]', fontsize=12)
ax2.set_ylabel(r'Plasma Frequency $f_{pe}$ [Hz]', fontsize=12, color='b')
ax2_twin.set_ylabel(r'Plasma Period $\tau_{pe}$ [ps]', fontsize=12, color='r')
ax2.set_title(r'Electron Plasma Oscillation Timescale', fontsize=13, fontweight='bold')
ax2.tick_params(axis='y', labelcolor='b')
ax2_twin.tick_params(axis='y', labelcolor='r')
ax2.grid(True, alpha=0.3, which='both')
ax2.legend(loc='upper left', fontsize=10)
ax2_twin.legend(loc='lower right', fontsize=10)

# Plot 3: Plasma parameter N_D
T_e_range = np.logspace(-1, 4, 200)  # 0.1 to 10,000 eV
N_D_grid = np.zeros((len(n_e_range), len(T_e_range)))

for i, n in enumerate(n_e_range):
    for j, T in enumerate(T_e_range):
        N_D_grid[i, j] = plasma_parameter(n, T)

# Create contour plot
levels = np.logspace(-2, 12, 15)
cs = ax3.contourf(T_e_range, n_e_range, N_D_grid, levels=levels,
                  cmap='RdYlGn', norm=plt.matplotlib.colors.LogNorm())
ax3.contour(T_e_range, n_e_range, N_D_grid, levels=[1], colors='black',
            linewidths=3, linestyles='--')

# Mark regime boundaries
ax3.annotate(r'$N_D < 1$: Not a Plasma', xy=(0.2, 1e8), fontsize=11,
            bbox=dict(boxstyle='round', facecolor='red', alpha=0.7))
ax3.annotate(r'$N_D \gg 1$: Collective Behavior', xy=(50, 1e18), fontsize=11,
            bbox=dict(boxstyle='round', facecolor='green', alpha=0.7))

ax3.set_xscale('log')
ax3.set_yscale('log')
ax3.set_xlabel(r'Electron Temperature $T_e$ [eV]', fontsize=12)
ax3.set_ylabel(r'Electron Density $n_e$ [m$^{-3}$]', fontsize=12)
ax3.set_title(r'Plasma Parameter $N_D = n_e \lambda_D^3$', fontsize=13, fontweight='bold')
cbar = plt.colorbar(cs, ax=ax3)
cbar.set_label(r'$N_D$ (dimensionless)', fontsize=11)

# Plot regime points
for regime, (n, T) in regimes.items():
    N_D_val = plasma_parameter(n, T)
    if N_D_val > 1:
        ax3.plot(T, n, 'o', markersize=6, color='blue', markeredgecolor='black')

# Plot 4: Ratio of Debye length to system size
# For tokamak-like parameters
n_tokamak = np.logspace(18, 21, 100)
T_tokamak = np.logspace(0, 4, 100)
system_size = 1.0  # 1 meter characteristic size

lambda_D_normalized = np.zeros((len(n_tokamak), len(T_tokamak)))
for i, n in enumerate(n_tokamak):
    for j, T in enumerate(T_tokamak):
        lambda_D_normalized[i, j] = debye_length(n, T) / system_size

levels_norm = np.logspace(-8, 0, 9)
cs2 = ax4.contourf(T_tokamak, n_tokamak, lambda_D_normalized,
                   levels=levels_norm, cmap='viridis',
                   norm=plt.matplotlib.colors.LogNorm())

ax4.set_xscale('log')
ax4.set_yscale('log')
ax4.set_xlabel(r'Temperature $T_e$ [eV]', fontsize=12)
ax4.set_ylabel(r'Density $n_e$ [m$^{-3}$]', fontsize=12)
ax4.set_title(r'Debye Length Normalized to System Size (1 m)', fontsize=13, fontweight='bold')
cbar2 = plt.colorbar(cs2, ax=ax4)
cbar2.set_label(r'$\lambda_D / L$ (dimensionless)', fontsize=11)

# Mark fusion-relevant regime
ax4.axvspan(1000, 10000, alpha=0.2, color='red', label='Fusion Core')
ax4.axhspan(1e19, 1e21, alpha=0.2, color='blue', label='Fusion Density')
ax4.legend(loc='lower left', fontsize=9)

plt.tight_layout()
plt.savefig('plasma_parameters_debye_plasma_freq.pdf', dpi=300, bbox_inches='tight')
plt.close()

# Store key values for table
debye_values = {}
for regime, (n, T) in regimes.items():
    debye_values[regime] = {
        'n_e': n,
        'T_e': T,
        'lambda_D': debye_length(n, T),
        'omega_pe': plasma_frequency(n),
        'N_D': plasma_parameter(n, T)
    }
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.98\textwidth]{plasma_parameters_debye_plasma_freq.pdf}
\caption{Fundamental plasma parameters across density-temperature space. \textbf{Top Left:} Debye shielding length decreases with increasing density as more particles are available to screen perturbations, ranging from millimeters in the ionosphere to nanometers in compressed fusion fuel. Representative regimes span eight orders of magnitude in density. \textbf{Top Right:} Plasma frequency and corresponding oscillation period, showing the fundamental timescale for electrostatic oscillations increases from femtoseconds at solid-density to microseconds in space plasmas. \textbf{Bottom Left:} Plasma parameter $N_D$ in density-temperature space, with the critical boundary $N_D = 1$ (dashed black line) separating strongly-coupled systems from weakly-coupled plasmas where collective effects dominate. Most natural and laboratory plasmas occupy the $N_D \gg 1$ regime (green) where mean-field theory applies. \textbf{Bottom Right:} Ratio of Debye length to system size for fusion-relevant parameters, demonstrating that $\lambda_D/L \ll 1$ ensures quasi-neutrality and justifies the plasma approximation for macroscopic devices.}
\end{figure}

\section{Collisional Processes and Transport}

Coulomb collisions between charged particles determine transport coefficients and distinguish collisional from collisionless regimes. The 90-degree deflection collision frequency is:
\begin{equation}
\nu_{ei} = \frac{n_e e^4 \ln\Lambda}{12\pi^{3/2} \epsilon_0^2 m_e^{1/2} (k_B T_e)^{3/2}}
\end{equation}

where $\ln\Lambda$ is the Coulomb logarithm accounting for cumulative small-angle scattering. The electron-ion mean free path is $\lambda_{mfp} = v_{th,e} / \nu_{ei}$, where $v_{th,e} = \sqrt{2k_B T_e / m_e}$ is the thermal velocity.

The Spitzer electrical conductivity for a fully ionized plasma is:
\begin{equation}
\sigma_{Spitzer} = 1.96 \times 10^4 \, T_e^{3/2} \, \text{S/m}
\end{equation}

with $T_e$ in eV. This strong temperature dependence ($\sigma \propto T^{3/2}$) arises because faster electrons scatter less frequently and have higher momentum.

\begin{pycode}
# Calculate collision properties
T_e_coll = np.logspace(-1, 4, 200)  # 0.1 to 10,000 eV

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 12))

# Plot 1: Collision frequency vs temperature for different densities
densities = [1e18, 1e19, 1e20, 1e21]  # m^-3
colors = ['blue', 'green', 'orange', 'red']

for n, color in zip(densities, colors):
    nu_ei = collision_frequency(n, T_e_coll)
    omega_pe_val = plasma_frequency(n)

    ax1.loglog(T_e_coll, nu_ei, linewidth=2.5, color=color,
              label=f'$n_e = {n:.0e}$ m$^{{-3}}$')
    ax1.axhline(omega_pe_val, linestyle='--', linewidth=1.5, color=color, alpha=0.5)

ax1.set_xlabel(r'Electron Temperature $T_e$ [eV]', fontsize=12)
ax1.set_ylabel(r'Collision Frequency $\nu_{ei}$ [Hz]', fontsize=12)
ax1.set_title(r'Electron-Ion Collision Frequency (dashed: $\omega_{pe}$)', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3, which='both')
ax1.set_ylim([1e6, 1e16])

# Plot 2: Collisionality parameter nu/omega
for n, color in zip(densities, colors):
    nu_ei = collision_frequency(n, T_e_coll)
    omega_pe_val = plasma_frequency(n)
    collisionality = nu_ei / omega_pe_val

    ax2.loglog(T_e_coll, collisionality, linewidth=2.5, color=color,
              label=f'$n_e = {n:.0e}$ m$^{{-3}}$')

ax2.axhline(1, linestyle='--', linewidth=2, color='black', label=r'$\nu/\omega = 1$')
ax2.fill_between(T_e_coll, 1e-6, 1, alpha=0.2, color='blue', label='Collisionless')
ax2.fill_between(T_e_coll, 1, 1e6, alpha=0.2, color='red', label='Collisional')

ax2.set_xlabel(r'Electron Temperature $T_e$ [eV]', fontsize=12)
ax2.set_ylabel(r'Collisionality $\nu_{ei}/\omega_{pe}$', fontsize=12)
ax2.set_title(r'Collisionality Parameter: Fluid vs Kinetic Regimes', fontsize=13, fontweight='bold')
ax2.legend(fontsize=9, loc='upper right')
ax2.grid(True, alpha=0.3, which='both')
ax2.set_ylim([1e-3, 1e3])

# Plot 3: Mean free path vs Debye length
for n, color in zip(densities, colors):
    mfp = mean_free_path(n, T_e_coll)
    ld = debye_length(n, T_e_coll)
    ratio = mfp / ld

    ax3.loglog(T_e_coll, ratio, linewidth=2.5, color=color,
              label=f'$n_e = {n:.0e}$ m$^{{-3}}$')

ax3.axhline(1, linestyle='--', linewidth=2, color='black')
ax3.set_xlabel(r'Electron Temperature $T_e$ [eV]', fontsize=12)
ax3.set_ylabel(r'$\lambda_{mfp} / \lambda_D$', fontsize=12)
ax3.set_title(r'Mean Free Path Normalized to Debye Length', fontsize=13, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3, which='both')

# Plot 4: Spitzer conductivity
sigma_spitzer = spitzer_conductivity(T_e_coll)

ax4.loglog(T_e_coll, sigma_spitzer, linewidth=3, color='purple')
ax4.set_xlabel(r'Electron Temperature $T_e$ [eV]', fontsize=12)
ax4.set_ylabel(r'Spitzer Conductivity $\sigma$ [S/m]', fontsize=12)
ax4.set_title(r'Electrical Conductivity: $\sigma \propto T_e^{3/2}$', fontsize=13, fontweight='bold')
ax4.grid(True, alpha=0.3, which='both')

# Mark reference materials
copper_conductivity = 5.96e7  # S/m at room temperature
ax4.axhline(copper_conductivity, linestyle='--', linewidth=2, color='orange',
           label='Copper (300 K)')

# Find where plasma exceeds copper conductivity
T_exceed_copper = T_e_coll[sigma_spitzer > copper_conductivity][0]
ax4.plot(T_exceed_copper, copper_conductivity, 'ro', markersize=10)
ax4.annotate(f'$T_e = {T_exceed_copper:.1f}$ eV\n(plasma exceeds Cu)',
            xy=(T_exceed_copper, copper_conductivity),
            xytext=(20, -30), textcoords='offset points',
            fontsize=10, bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.7),
            arrowprops=dict(arrowstyle='->', lw=2))

ax4.legend(fontsize=10, loc='lower right')

plt.tight_layout()
plt.savefig('plasma_parameters_collisions.pdf', dpi=300, bbox_inches='tight')
plt.close()

# Calculate collision parameters for regimes
collision_data = {}
for regime, (n, T) in regimes.items():
    nu_ei = collision_frequency(n, T)
    omega_pe_val = plasma_frequency(n)
    mfp = mean_free_path(n, T)
    ld = debye_length(n, T)
    sigma = spitzer_conductivity(T)

    collision_data[regime] = {
        'nu_ei': nu_ei,
        'collisionality': nu_ei / omega_pe_val,
        'mfp': mfp,
        'mfp_over_ld': mfp / ld,
        'sigma': sigma
    }
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.98\textwidth]{plasma_parameters_collisions.pdf}
\caption{Collisional processes and transport properties. \textbf{Top Left:} Electron-ion collision frequency decreases as $T^{-3/2}$ due to reduced Coulomb scattering cross-section at higher velocities. Dashed horizontal lines indicate plasma frequencies $\omega_{pe}$ for each density; when $\nu_{ei} < \omega_{pe}$ the plasma oscillates many times between collisions. \textbf{Top Right:} Collisionality parameter $\nu_{ei}/\omega_{pe}$ separates collisional regimes ($> 1$, red shading) where fluid theory applies from collisionless regimes ($< 1$, blue shading) requiring kinetic treatment. Most fusion plasmas are deeply collisionless. \textbf{Bottom Left:} Ratio of electron mean free path to Debye length, demonstrating that collisions typically occur over distances much larger than the screening length, validating the Debye-Hückel approximation. \textbf{Bottom Right:} Spitzer conductivity increases as $T_e^{3/2}$, exceeding copper's conductivity above $\sim 8$ eV. Hot fusion plasmas are excellent electrical conductors despite their low density.}
\end{figure}

\section{Magnetized Plasma Dynamics}

In the presence of a magnetic field $\mathbf{B}$, charged particles execute helical trajectories with cyclotron frequency:
\begin{equation}
\omega_c = \frac{eB}{m}
\end{equation}

The Larmor radius (gyroradius) for a particle with perpendicular velocity $v_\perp$ is:
\begin{equation}
r_L = \frac{v_\perp}{\omega_c} = \frac{m v_\perp}{eB}
\end{equation}

For a thermal distribution, the characteristic thermal Larmor radius is $r_{L,th} = v_{th} / \omega_c$, where $v_{th} = \sqrt{2k_B T/m}$.

The magnetization parameter $\omega_c \tau_{coll} = \omega_c / \nu$ determines whether particles complete many gyro-orbits between collisions. When $\omega_c \tau_{coll} \gg 1$, transport is highly anisotropic with diffusion predominantly along field lines.

\begin{pycode}
# Magnetic field strengths
B_range = np.logspace(-6, 2, 200)  # 1 microTesla to 100 Tesla

# Temperature range
T_range = np.logspace(-1, 4, 200)  # 0.1 to 10,000 eV

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 12))

# Plot 1: Larmor radius vs B field for different temperatures
for T in [1, 10, 100, 1000, 5000]:
    r_L_electron = thermal_larmor_radius(T, B_range, m_e) * 1e3  # mm
    r_L_proton = thermal_larmor_radius(T, B_range, m_p) * 1e3    # mm

    ax1.loglog(B_range, r_L_electron, linewidth=2, label=f'e$^-$, {T} eV')
    ax1.loglog(B_range, r_L_proton, linewidth=2, linestyle='--', label=f'p$^+$, {T} eV')

# Mark typical magnetic fields
B_markers = {
    'Earth': 5e-5,
    'Tokamak': 5,
    'Stellarator': 2,
    'Z-pinch': 10,
    'Magnetar': 1e8  # Off scale but conceptual
}

for name, B_val in B_markers.items():
    if B_val < 100:
        ax1.axvline(B_val, linestyle=':', alpha=0.5, color='gray')
        ax1.text(B_val, 1e2, name, rotation=90, va='bottom', fontsize=8)

ax1.set_xlabel(r'Magnetic Field $B$ [T]', fontsize=12)
ax1.set_ylabel(r'Thermal Larmor Radius $r_L$ [mm]', fontsize=12)
ax1.set_title(r'Larmor Radius: Electron (solid) vs Proton (dashed)', fontsize=13, fontweight='bold')
ax1.legend(fontsize=8, ncol=2, loc='upper right')
ax1.grid(True, alpha=0.3, which='both')
ax1.set_xlim([1e-6, 1e2])
ax1.set_ylim([1e-4, 1e4])

# Plot 2: Cyclotron frequency
omega_ce = cyclotron_frequency(B_range, m_e) / (2 * np.pi)  # Hz
omega_ci = cyclotron_frequency(B_range, m_p) / (2 * np.pi)  # Hz

ax2.loglog(B_range, omega_ce, linewidth=3, color='blue', label='Electron')
ax2.loglog(B_range, omega_ci, linewidth=3, color='red', label='Proton')

ax2.set_xlabel(r'Magnetic Field $B$ [T]', fontsize=12)
ax2.set_ylabel(r'Cyclotron Frequency $f_c$ [Hz]', fontsize=12)
ax2.set_title(r'Cyclotron Frequency: $f_c = eB/(2\pi m)$', fontsize=13, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3, which='both')

# Mark frequency bands
ax2.axhspan(1e6, 1e9, alpha=0.1, color='green', label='RF range')
ax2.axhspan(1e9, 3e11, alpha=0.1, color='yellow', label='Microwave')

# Plot 3: Magnetization parameter (omega_c / nu)
n_mag = 1e20  # Fixed density for tokamak
T_mag = np.logspace(0, 4, 200)

magnetization_electron = np.zeros((len(B_range), len(T_mag)))
magnetization_ion = np.zeros((len(B_range), len(T_mag)))

for i, B in enumerate(B_range):
    for j, T in enumerate(T_mag):
        omega_ce_val = cyclotron_frequency(B, m_e)
        omega_ci_val = cyclotron_frequency(B, m_p)
        nu_ei = collision_frequency(n_mag, T)

        magnetization_electron[i, j] = omega_ce_val / nu_ei
        magnetization_ion[i, j] = omega_ci_val / nu_ei

# Plot electron magnetization
levels_mag = np.logspace(-2, 6, 9)
cs = ax3.contourf(T_mag, B_range, magnetization_electron,
                  levels=levels_mag, cmap='plasma',
                  norm=plt.matplotlib.colors.LogNorm())
ax3.contour(T_mag, B_range, magnetization_electron, levels=[1],
           colors='white', linewidths=3, linestyles='--')

ax3.set_xscale('log')
ax3.set_yscale('log')
ax3.set_xlabel(r'Temperature $T_e$ [eV]', fontsize=12)
ax3.set_ylabel(r'Magnetic Field $B$ [T]', fontsize=12)
ax3.set_title(r'Electron Magnetization $\omega_{ce}/\nu_{ei}$ ($n_e = 10^{20}$ m$^{-3}$)',
             fontsize=13, fontweight='bold')
cbar = plt.colorbar(cs, ax=ax3)
cbar.set_label(r'$\omega_{ce}/\nu_{ei}$', fontsize=11)

# Mark fusion regime
ax3.axvspan(1000, 10000, alpha=0.2, color='cyan')
ax3.axhspan(1, 10, alpha=0.2, color='cyan')

# Plot 4: Ratio of Larmor radius to system size
system_size_fusion = 1.0  # 1 meter
r_L_normalized = thermal_larmor_radius(1000, B_range, m_e) / system_size_fusion

ax4.loglog(B_range, r_L_normalized, linewidth=3, color='purple', label='Electron (1 keV)')

r_L_normalized_ion = thermal_larmor_radius(1000, B_range, m_p) / system_size_fusion
ax4.loglog(B_range, r_L_normalized_ion, linewidth=3, color='orange',
          linestyle='--', label='Proton (1 keV)')

ax4.axhline(0.01, linestyle=':', linewidth=2, color='green',
           label=r'$r_L/L = 0.01$ (good confinement)')
ax4.axhline(1, linestyle=':', linewidth=2, color='red',
           label=r'$r_L/L = 1$ (poor confinement)')

ax4.set_xlabel(r'Magnetic Field $B$ [T]', fontsize=12)
ax4.set_ylabel(r'$r_L / L$ (dimensionless)', fontsize=12)
ax4.set_title(r'Larmor Radius Normalized to System Size (1 m)', fontsize=13, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3, which='both')
ax4.set_xlim([1e-6, 1e2])
ax4.set_ylim([1e-6, 1e4])

# Mark typical fusion fields
for name, B_val in [('Earth', 5e-5), ('Tokamak', 5), ('Stellarator', 2)]:
    if B_val < 100:
        ax4.axvline(B_val, linestyle=':', alpha=0.5, color='gray', linewidth=1)

plt.tight_layout()
plt.savefig('plasma_parameters_magnetized.pdf', dpi=300, bbox_inches='tight')
plt.close()

# Calculate magnetization parameters
magnetic_data = {}
magnetic_fields = {
    'Earth Magnetosphere': 5e-5,
    'Tokamak (ITER)': 5.3,
    'Stellarator (W7-X)': 2.5,
    'RFP': 0.5,
    'Laser Plasma': 100
}

for field_name, B in magnetic_fields.items():
    for regime, (n, T) in regimes.items():
        if 'Tokamak' in regime or 'Solar' in regime:
            omega_ce_val = cyclotron_frequency(B, m_e)
            omega_ci_val = cyclotron_frequency(B, m_p)
            nu_ei = collision_frequency(n, T)
            r_L_e = thermal_larmor_radius(T, B, m_e)
            r_L_i = thermal_larmor_radius(T, B, m_p)

            key = f"{regime}_{field_name}"
            magnetic_data[key] = {
                'B': B,
                'omega_ce': omega_ce_val,
                'omega_ci': omega_ci_val,
                'magnetization_e': omega_ce_val / nu_ei if nu_ei > 0 else np.inf,
                'r_L_e': r_L_e,
                'r_L_i': r_L_i
            }
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.98\textwidth]{plasma_parameters_magnetized.pdf}
\caption{Magnetized plasma parameters and confinement metrics. \textbf{Top Left:} Thermal Larmor radius decreases with magnetic field strength, with electrons (solid) gyrating $\sqrt{m_p/m_e} \approx 43$ times tighter than protons (dashed). Vertical gray lines mark typical field strengths from Earth's magnetosphere ($50 \mu$T) to tokamaks ($\sim 5$ T). At 5 T and 1 keV, electrons have sub-millimeter gyroradii while ions have centimeter-scale orbits. \textbf{Top Right:} Cyclotron frequencies span radio (MHz) to microwave (GHz) bands, determining resonant heating schemes used in fusion devices. \textbf{Bottom Left:} Electron magnetization parameter $\omega_{ce}/\nu_{ei}$ in temperature-field space for fusion density $n_e = 10^{20}$ m$^{-3}$. White dashed contour at unity separates weakly magnetized (fluid-like) from strongly magnetized (guiding-center) regimes. Cyan shading indicates typical tokamak operating space where $\omega_{ce}/\nu_{ei} \sim 10^3$. \textbf{Bottom Right:} Normalized Larmor radius determines confinement quality; $r_L/L < 0.01$ (green line) indicates good magnetic confinement, while $r_L/L \sim 1$ (red line) implies particles reach walls before completing orbits. Multi-Tesla fields are required to confine keV plasmas in meter-scale devices.}
\end{figure}

\section{Wave Propagation in Plasmas}

Plasmas support a rich spectrum of wave modes arising from coupling between electromagnetic fields, pressure gradients, and magnetic tension. The three fundamental modes are:

\textbf{1. Langmuir Waves (Electron Plasma Waves):} Electrostatic oscillations at the plasma frequency, representing collective electron motion against a stationary ion background. Dispersion relation:
\begin{equation}
\omega^2 = \omega_{pe}^2 + \frac{3}{2} k^2 v_{th,e}^2
\end{equation}

\textbf{2. Ion Acoustic Waves:} Electrostatic sound waves in which ion inertia provides the restoring force and electron pressure provides stiffness. Propagation velocity:
\begin{equation}
c_s = \sqrt{\frac{\gamma_e k_B T_e + \gamma_i k_B T_i}{m_i}}
\end{equation}

For $T_e \gg T_i$ and isothermal electrons ($\gamma_e = 1$), this simplifies to $c_s \approx \sqrt{k_B T_e / m_i}$.

\textbf{3. Alfvén Waves:} Transverse electromagnetic waves propagating along magnetic field lines with velocity determined by magnetic tension:
\begin{equation}
v_A = \frac{B}{\sqrt{\mu_0 \rho}} = \frac{B}{\sqrt{\mu_0 n_i m_i}}
\end{equation}

These waves are crucial for energy transport in magnetized plasmas and play a central role in solar corona heating and magnetospheric dynamics.

\begin{pycode}
# Wave propagation calculations
n_wave = np.logspace(16, 22, 200)  # Density range for wave analysis

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 12))

# Plot 1: Phase velocities of different waves
T_wave = 100  # eV
B_wave = 1.0  # Tesla

omega_pe_wave = plasma_frequency(n_wave)
f_pe_wave = omega_pe_wave / (2 * np.pi)

# Ion acoustic velocity (assuming T_i = 0.1 * T_e)
c_s = ion_acoustic_velocity(T_wave, 0.1 * T_wave, gamma_e=1, gamma_i=3)

# Alfven velocity
v_A_array = alfven_velocity(B_wave, n_wave)

# Electron thermal velocity
v_th_e = np.sqrt(2 * T_wave * e / m_e)

# Phase velocity of Langmuir wave at k*lambda_D = 0.3
k_lambda_D = 0.3
lambda_D_array = debye_length(n_wave, T_wave)
k_array = k_lambda_D / lambda_D_array
v_ph_langmuir = omega_pe_wave / k_array

ax1.loglog(n_wave, v_A_array * 1e-3, linewidth=3, color='blue', label=r'Alfv\'en $v_A$')
ax1.axhline(c_s * 1e-3, linewidth=3, color='green', label=r'Ion Acoustic $c_s$')
ax1.axhline(v_th_e * 1e-3, linewidth=3, color='red', linestyle='--', label=r'Electron Thermal $v_{th,e}$')
ax1.loglog(n_wave, v_ph_langmuir * 1e-3, linewidth=2.5, color='purple',
          label=r'Langmuir $v_{ph}$ ($k\lambda_D = 0.3$)')

ax1.set_xlabel(r'Density $n$ [m$^{-3}$]', fontsize=12)
ax1.set_ylabel(r'Velocity [km/s]', fontsize=12)
ax1.set_title(f'Wave Phase Velocities ($T_e = {T_wave}$ eV, $B = {B_wave}$ T)',
             fontsize=13, fontweight='bold')
ax1.legend(fontsize=11, loc='best')
ax1.grid(True, alpha=0.3, which='both')

# Speed of light reference
ax1.axhline(c * 1e-3, linewidth=2, color='black', linestyle=':', alpha=0.5, label='Speed of light')

# Plot 2: Langmuir wave dispersion
k_range = np.logspace(-3, 2, 300)  # k * lambda_D
omega_langmuir = np.zeros((len(k_range), 3))  # Three temperatures

T_dispersion = [10, 100, 1000]  # eV
n_dispersion = 1e20  # m^-3

for idx, T_disp in enumerate(T_dispersion):
    omega_pe_disp = plasma_frequency(n_dispersion)
    lambda_D_disp = debye_length(n_dispersion, T_disp)
    v_th_e_disp = np.sqrt(2 * T_disp * e / m_e)

    omega_langmuir[:, idx] = np.sqrt(omega_pe_disp**2 + 1.5 * (k_range / lambda_D_disp * v_th_e_disp)**2)

    ax2.loglog(k_range, omega_langmuir[:, idx] / omega_pe_disp, linewidth=2.5,
              label=f'$T_e = {T_disp}$ eV')

ax2.axhline(1, linestyle='--', linewidth=2, color='black', alpha=0.5)
ax2.set_xlabel(r'$k\lambda_D$ (dimensionless)', fontsize=12)
ax2.set_ylabel(r'$\omega / \omega_{pe}$ (dimensionless)', fontsize=12)
ax2.set_title(r'Langmuir Wave Dispersion: $\omega^2 = \omega_{pe}^2 + \frac{3}{2}k^2 v_{th}^2$',
             fontsize=13, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3, which='both')
ax2.set_xlim([1e-3, 1e2])

# Plot 3: Alfvén velocity across plasma regimes
B_alfven = np.logspace(-6, 1, 200)  # 1 microT to 10 T
n_alfven_values = [1e18, 1e19, 1e20, 1e21]

for n_alf in n_alfven_values:
    v_A_sweep = alfven_velocity(B_alfven, n_alf)
    ax3.loglog(B_alfven, v_A_sweep * 1e-3, linewidth=2.5, label=f'$n = {n_alf:.0e}$ m$^{{-3}}$')

# Mark regimes where v_A is comparable to other velocities
ax3.axhline(c_s * 1e-3, linestyle=':', linewidth=2, color='green', alpha=0.5,
           label=f'$c_s$ ({T_wave} eV)')
ax3.axhline(1e4, linestyle=':', linewidth=2, color='red', alpha=0.5,
           label='$10^4$ km/s')

ax3.set_xlabel(r'Magnetic Field $B$ [T]', fontsize=12)
ax3.set_ylabel(r'Alfv\'en Velocity $v_A$ [km/s]', fontsize=12)
ax3.set_title(r'Alfv\'en Velocity: $v_A = B/\sqrt{\mu_0 n m_p}$', fontsize=13, fontweight='bold')
ax3.legend(fontsize=10, loc='lower right')
ax3.grid(True, alpha=0.3, which='both')

# Plot 4: Plasma beta parameter
T_beta = np.logspace(-1, 4, 200)
B_beta = 5.0  # Tesla (tokamak field)
n_beta = 1e20  # m^-3

# Thermal pressure
p_thermal = n_beta * k_B * T_beta * e  # Pascals

# Magnetic pressure
p_magnetic = B_beta**2 / (2 * mu_0)

# Beta parameter
beta = p_thermal / p_magnetic

ax4.loglog(T_beta, beta, linewidth=3, color='darkblue')
ax4.axhline(1, linestyle='--', linewidth=2, color='black', label=r'$\beta = 1$')
ax4.fill_between(T_beta, 1e-6, 1, alpha=0.2, color='blue', label=r'$\beta < 1$ (Low-$\beta$)')
ax4.fill_between(T_beta, 1, 1e2, alpha=0.2, color='red', label=r'$\beta > 1$ (High-$\beta$)')

ax4.set_xlabel(r'Temperature $T$ [eV]', fontsize=12)
ax4.set_ylabel(r'Plasma Beta $\beta = p_{thermal}/p_{magnetic}$', fontsize=12)
ax4.set_title(f'Plasma Beta Parameter ($n = {n_beta:.0e}$ m$^{{-3}}$, $B = {B_beta}$ T)',
             fontsize=13, fontweight='bold')
ax4.legend(fontsize=11)
ax4.grid(True, alpha=0.3, which='both')

# Mark fusion-relevant temperature
T_fusion = 10000  # eV (10 keV)
beta_fusion = (n_beta * k_B * T_fusion * e) / (B_beta**2 / (2 * mu_0))
ax4.plot(T_fusion, beta_fusion, 'ro', markersize=10)
ax4.annotate(f'Fusion core\n$\\beta = {beta_fusion:.3f}$', xy=(T_fusion, beta_fusion),
            xytext=(30, -30), textcoords='offset points',
            fontsize=10, bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.7),
            arrowprops=dict(arrowstyle='->', lw=2))

plt.tight_layout()
plt.savefig('plasma_parameters_waves.pdf', dpi=300, bbox_inches='tight')
plt.close()

# Calculate wave properties for table
wave_data = {}
for regime, (n, T) in regimes.items():
    if 'Tokamak' in regime or 'Solar' in regime:
        omega_pe_val = plasma_frequency(n)
        c_s_val = ion_acoustic_velocity(T, 0.1 * T)

        # Assume B field for each regime
        if 'Tokamak' in regime:
            B_val = 5.0
        elif 'Solar' in regime:
            B_val = 0.01
        else:
            B_val = 1e-5

        v_A_val = alfven_velocity(B_val, n)

        wave_data[regime] = {
            'f_pe': omega_pe_val / (2 * np.pi),
            'c_s': c_s_val,
            'v_A': v_A_val,
            'B': B_val
        }
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.98\textwidth]{plasma_parameters_waves.pdf}
\caption{Wave propagation in plasmas across multiple modes. \textbf{Top Left:} Characteristic velocities for different wave types at 100 eV and 1 T. Alfvén velocity (blue) decreases with increasing density as more mass must be accelerated by magnetic tension. Ion acoustic velocity (green) is density-independent and scales as $\sqrt{T_e/m_i}$. Langmuir wave phase velocity (purple) exceeds thermal velocity due to thermal pressure corrections to the cold plasma dispersion. All remain well below the speed of light (black dotted). \textbf{Top Right:} Langmuir wave dispersion showing transition from cold plasma limit ($\omega \approx \omega_{pe}$ for $k\lambda_D \ll 1$) to kinetic regime where thermal effects dominate ($\omega \propto k$ for $k\lambda_D \gg 1$). Higher temperatures extend the transition region. \textbf{Bottom Left:} Alfvén velocity spans four orders of magnitude from $\sim 10$ km/s in dense fusion plasmas to $> 10^4$ km/s in tenuous magnetospheres. Intersection with ion acoustic velocity (green dotted at 440 km/s for 100 eV) marks regime boundaries for MHD wave coupling. \textbf{Bottom Right:} Plasma beta parameter $\beta = p_{thermal}/p_{magnetic}$ determines whether magnetic field (low-$\beta$, blue) or thermal pressure (high-$\beta$, red) dominates. Tokamaks operate at $\beta \sim 0.05$ (marked point at 10 keV), requiring strong magnetic confinement. The $\beta = 1$ boundary (black dashed) separates magnetic-pressure-dominated from kinetic-pressure-dominated regimes.}
\end{figure}

\section{Comprehensive Parameter Tables}

\begin{pycode}
# Generate comprehensive table of plasma parameters

print(r'\begin{table}[H]')
print(r'\centering')
print(r'\caption{Fundamental Plasma Parameters Across Regimes}')
print(r'\begin{tabular}{@{}lcccccc@{}}')
print(r'\toprule')
print(r'Regime & $n_e$ [m$^{-3}$] & $T_e$ [eV] & $\lambda_D$ [m] & $\omega_{pe}/(2\pi)$ [Hz] & $N_D$ & $\nu_{ei}/\omega_{pe}$ \\')
print(r'\midrule')

for regime, (n, T) in sorted(regimes.items()):
    ld = debye_length(n, T)
    f_pe = plasma_frequency(n) / (2 * np.pi)
    N_D = plasma_parameter(n, T)
    nu_ei = collision_frequency(n, T)
    omega_pe_val = plasma_frequency(n)
    coll_ratio = nu_ei / omega_pe_val

    print(f"{regime} & {n:.1e} & {T:.1e} & {ld:.2e} & {f_pe:.2e} & {N_D:.2e} & {coll_ratio:.2e} \\\\")

print(r'\bottomrule')
print(r'\end{tabular}')
print(r'\end{table}')

print()
print()

# Table 2: Magnetized plasma parameters
print(r'\begin{table}[H]')
print(r'\centering')
print(r'\caption{Magnetized Plasma Parameters}')
print(r'\begin{tabular}{@{}lccccc@{}}')
print(r'\toprule')
print(r'Configuration & $B$ [T] & $T_e$ [eV] & $r_{L,e}$ [mm] & $r_{L,i}$ [mm] & $\omega_{ce}/(2\pi)$ [GHz] \\')
print(r'\midrule')

mag_configs = [
    ('Tokamak (ITER)', 5.3, 10000),
    ('Stellarator (W7-X)', 2.5, 2000),
    ('RFP', 0.5, 500),
    ('Solar Corona', 0.01, 100),
    ('Earth Magnetosphere', 5e-5, 1),
]

for name, B, T in mag_configs:
    r_L_e = thermal_larmor_radius(T, B, m_e) * 1e3  # mm
    r_L_i = thermal_larmor_radius(T, B, m_p) * 1e3  # mm
    f_ce = cyclotron_frequency(B, m_e) / (2 * np.pi * 1e9)  # GHz

    print(f"{name} & {B:.2e} & {T:.0f} & {r_L_e:.3f} & {r_L_i:.2f} & {f_ce:.2f} \\\\")

print(r'\bottomrule')
print(r'\end{tabular')
print(r'\end{table}')

print()
print()

# Table 3: Wave velocities
print(r'\begin{table}[H]')
print(r'\centering')
print(r'\caption{Wave Propagation Velocities}')
print(r'\begin{tabular}{@{}lcccc@{}}')
print(r'\toprule')
print(r'Regime & $c_s$ [km/s] & $v_A$ [km/s] & $v_{th,e}$ [km/s] & $v_A/c_s$ \\')
print(r'\midrule')

wave_configs = [
    ('Tokamak Core', 1e20, 5000, 5.0),
    ('Solar Corona', 1e14, 100, 0.01),
    ('Solar Wind', 1e7, 10, 1e-5),
]

for name, n, T, B in wave_configs:
    c_s_val = ion_acoustic_velocity(T, 0.1 * T) * 1e-3  # km/s
    v_A_val = alfven_velocity(B, n) * 1e-3  # km/s
    v_th_e_val = np.sqrt(2 * T * e / m_e) * 1e-3  # km/s
    ratio = v_A_val / c_s_val

    print(f"{name} & {c_s_val:.1f} & {v_A_val:.1f} & {v_th_e_val:.0f} & {ratio:.2f} \\\\")

print(r'\bottomrule')
print(r'\end{tabular}')
print(r'\end{table}')

\end{pycode}

\section{Dimensionless Parameters and Regime Classification}

The physics of plasmas is governed by dimensionless ratios constructed from fundamental scales:

\begin{enumerate}
\item \textbf{Plasma Parameter:} $N_D = n_e \lambda_D^3$ determines coupling strength. For $N_D \gg 1$, collective effects dominate and mean-field theory applies. When $N_D \sim 1$, strong correlation effects become important.

\item \textbf{Collisionality:} $\nu_{ei} / \omega_{pe}$ distinguishes collisional ($> 1$) from collisionless ($< 1$) regimes. Collisionless plasmas require kinetic theory; collisional plasmas admit fluid descriptions.

\item \textbf{Magnetization:} $\omega_{ce} / \nu_{ei}$ determines whether particles complete gyro-orbits between collisions. For $\omega_{ce}/\nu_{ei} \gg 1$, transport is anisotropic and confined along field lines.

\item \textbf{Plasma Beta:} $\beta = p_{thermal} / p_{magnetic} = 2\mu_0 n k_B T / B^2$ measures thermal to magnetic pressure. Low-$\beta$ plasmas ($\beta \ll 1$) are magnetically dominated; high-$\beta$ plasmas are kinetically dominated.

\item \textbf{Normalized Larmor Radius:} $r_L / L$ determines gyrokinetic ordering. For $r_L/L \ll 1$, guiding center approximations are valid and magnetic confinement is effective.
\end{enumerate}

These ratios determine which physical processes are important and which theoretical approximations are valid.

\section{Conclusions}

This comprehensive analysis of fundamental plasma parameters demonstrates the multi-scale nature of plasma physics across twelve orders of magnitude in density (from $10^{6}$ m$^{-3}$ in the solar wind to $10^{26}$ m$^{-3}$ in compressed inertial fusion targets) and five orders of magnitude in temperature (from 0.1 eV in the ionosphere to 10 keV in tokamak cores).

The Debye length, ranging from nanometers to meters, sets the electrostatic shielding scale and determines when quasi-neutrality applies ($\lambda_D/L \ll 1$). The plasma frequency, spanning microhertz to petahertz, establishes the fundamental timescale for charge separation oscillations. All plasma regimes analyzed satisfy $N_D \gg 1$, confirming that collective effects dominate individual particle interactions.

Collision frequencies decrease sharply with temperature ($\nu_{ei} \propto T_e^{-3/2}$), rendering hot fusion plasmas deeply collisionless with $\nu_{ei}/\omega_{pe} \sim 10^{-3}$. This necessitates kinetic treatment of transport and turbulence. The Spitzer conductivity increases as $T_e^{3/2}$, exceeding copper's conductivity above 8 eV and reaching $10^7$ S/m in fusion cores.

In magnetized configurations, the electron Larmor radius scales as $r_{L,e} \propto \sqrt{T_e}/B$, requiring multi-Tesla fields to achieve $r_{L,e}/a < 0.01$ for effective confinement of keV plasmas. The magnetization parameter $\omega_{ce}/\nu_{ei} \sim 10^{3}$ in tokamaks ensures particles complete many gyro-orbits between collisions, enabling classical transport theory.

Wave analysis reveals that Alfvén velocities span $10^{1}$ to $10^{4}$ km/s depending on density and field strength, while ion acoustic velocities remain near $\sqrt{T_e/m_i} \sim 100$ km/s. The plasma beta in tokamaks is typically $\beta \sim 0.05$, indicating magnetic pressure dominance, while solar corona and magnetosphere plasmas can achieve $\beta \sim 1$.

These fundamental parameters provide the foundation for understanding collective phenomena including turbulence, waves, instabilities, and transport across all plasma regimes from astrophysical to laboratory systems.

\begin{thebibliography}{99}

\bibitem{chen1984introduction}
F.F. Chen, \textit{Introduction to Plasma Physics and Controlled Fusion}, 2nd ed., Plenum Press, New York (1984).

\bibitem{goldston1995introduction}
R.J. Goldston and P.H. Rutherford, \textit{Introduction to Plasma Physics}, Institute of Physics Publishing, Bristol (1995).

\bibitem{debye1923theory}
P. Debye and E. Hückel, ``The theory of electrolytes. I. Lowering of freezing point and related phenomena,'' \textit{Phys. Z.} \textbf{24}, 185 (1923).

\bibitem{spitzer1956physics}
L. Spitzer, Jr., \textit{Physics of Fully Ionized Gases}, Interscience Publishers, New York (1956).

\bibitem{tonks1929oscillations}
L. Tonks and I. Langmuir, ``Oscillations in ionized gases,'' \textit{Phys. Rev.} \textbf{33}, 195 (1929).

\bibitem{krall1973principles}
N.A. Krall and A.W. Trivelpiece, \textit{Principles of Plasma Physics}, McGraw-Hill, New York (1973).

\bibitem{alfven1942existence}
H. Alfvén, ``Existence of electromagnetic-hydrodynamic waves,'' \textit{Nature} \textbf{150}, 405 (1942).

\bibitem{wesson2011tokamaks}
J. Wesson and D.J. Campbell, \textit{Tokamaks}, 4th ed., Oxford University Press (2011).

\bibitem{landau1946oscillations}
L.D. Landau, ``On the vibrations of the electronic plasma,'' \textit{J. Phys. USSR} \textbf{10}, 25 (1946).

\bibitem{bohm1949minimum}
D. Bohm, E.H.S. Burhop, and H.S.W. Massey, ``The use of probes for plasma exploration in strong magnetic fields,'' in \textit{The Characteristics of Electrical Discharges in Magnetic Fields}, A. Guthrie and R.K. Wakerling, eds., McGraw-Hill (1949).

\bibitem{fried1961longitudinal}
B.D. Fried and S.D. Conte, \textit{The Plasma Dispersion Function}, Academic Press, New York (1961).

\bibitem{stix1992waves}
T.H. Stix, \textit{Waves in Plasmas}, American Institute of Physics, New York (1992).

\bibitem{bellan2006fundamentals}
P.M. Bellan, \textit{Fundamentals of Plasma Physics}, Cambridge University Press (2006).

\bibitem{freidberg2014plasma}
J.P. Freidberg, \textit{Ideal MHD}, Cambridge University Press (2014).

\bibitem{hazeltine2003plasma}
R.D. Hazeltine and J.D. Meiss, \textit{Plasma Confinement}, Dover Publications (2003).

\bibitem{nrl2019formulary}
NRL Plasma Formulary, Naval Research Laboratory, Washington DC (2019).

\bibitem{miyamoto2005plasma}
K. Miyamoto, \textit{Plasma Physics for Nuclear Fusion}, MIT Press (2005).

\bibitem{helander2002collisional}
P. Helander and D.J. Sigmar, \textit{Collisional Transport in Magnetized Plasmas}, Cambridge University Press (2002).

\bibitem{fitzpatrick2014plasma}
R. Fitzpatrick, \textit{Plasma Physics: An Introduction}, CRC Press (2014).

\bibitem{bittencourt2004fundamentals}
J.A. Bittencourt, \textit{Fundamentals of Plasma Physics}, 3rd ed., Springer (2004).

\end{thebibliography}

\end{document}