PlasmaBO

Installation

using Pkg
Pkg.add("PlasmaBO")

Problem Definition

The core problem is solving the linear dispersion relation for waves in a magnetized plasma:

\[\det(\mathbf{D}(ω, \mathbf{k})) = \left|\mathbf{K}(ω, \mathbf{k})+\left(k k-k^2 \mathbf{I}\right) \frac{c^2}{ω^2}\right|=0\]

where $\mathbf{D}$ is the dispersion tensor, $ω$ is the complex frequency, $\mathbf{k}$ is the wave vector, $\mathbf{K}=\mathbf{I}+\mathbf{Q}$, $\mathbf{Q}=-\frac{\mathbf{\sigma}}{i ω \epsilon_0}$, and

\[\mathbf{\sigma}=-i \sum_s \frac{q_s^2 n_s}{m_s} \sum_{n=-\infty}^{\infty} \int_{-\infty}^{\infty} \int_0^{\infty} \frac{2 \pi v_{\perp} d v_{\perp} d v_{\|}}{ω-n ω_{c s}-k_{\|} v_{\|}} \mathbf{\Pi}_s\]

with

\[\mathbf{\Pi}_s=\left[\begin{array}{ccc} A_s \frac{n^2 v_{\perp}}{\mu_s^2} J_n^2 & i A_s \frac{n v_{\perp}}{\mu_s} J_n J_n^{\prime} & B_s \frac{n v_{\perp}}{\mu_s} J_n^2 \\ -i A_s \frac{n v_{\perp}}{\mu_s} J_n J_n^{\prime} & A_s v_{\perp} J_n^{\prime 2} & -i B_s v_{\perp} J_n J_n^{\prime} \\ A_s \frac{n v_{\|}}{\mu_s} J_n^2 & i A_s v_{\|} J_n J_n^{\prime} & B_s v_{\|} J_n^2 \end{array}\right],\]

where $\mu_s=\frac{k_{\perp} v_{\perp}}{ω_{c s}}, A_s=\left(1-\frac{k_{\|} v_{\|}}{ω}\right) \frac{\partial f_{f_0}}{\partial v_{\perp}}+\frac{k_{\|} v_{\perp}}{ω} \frac{\partial f_{s 0}}{\partial v_{\|}}$, and $B_s=\frac{n ω_{c s} v_{\|}}{ω v_{\perp}} \frac{\partial f_{s 0}}{\partial v_{\perp}} +\left(1-\frac{n ω_{ω s}}{ω}\right) \frac{\partial f_{s_0}}{\partial v_{\|}}$.

This formulation is valid for non-relativistic, arbitrary gyrotropic distributions.

Input:

Plasma Parameters: For each species $s$, density $n_s$, temperature $T_s$ (anisotropic $T_{\parallel}, T_{\perp}$), drift velocity $v_{d,s}$, charge $q_s$, and mass $m_s$.
- Note: When arbitrary velocity distributions are provided, temperature are not needed as they are encoded in the distribution.
Coordinates & Fields:
- \[z\]
  : direction parallel to the background magnetic field ($B_0$).
- \[x\]
  : one perpendicular direction (any perpendicular direction is equivalent).
- Velocities: vz (parallel), vx (perpendicular), thermal speeds vtz, vtx.
Wave Parameters: Wave vector $\mathbf{k} = (k_\perp, 0, k_\parallel)$, where $θ$ is the propagation angle between $k$ and $B_0$, and $k_∥ = k \cos θ$, $k_⊥ = k \sinθ$.

Output:

Complex Frequency: $ω = ω_r + i\gamma$.
- \[ω_r\]
  : Real frequency (oscillation).
- \[\gamma\]
  : Growth rate (instability if $\gamma > 0$, damping if $\gamma < 0$).

Methodology

Solving the kinetic dispersion relation is traditionally challenging because:

Transcendental Equations: The dispersion tensor involves the Plasma Dispersion Function $Z(\zeta)$, making the equation $D(ω) = 0$ transcendental and highly nonlinear.
Root Finding Difficulty: Traditional iterative methods (like Newton-Raphson) require good initial guesses. Without them, it's easy to miss modes or converge to the wrong root.
Multimode Nature: Plasmas support multiple wave modes simultaneously (e.g., Alfvén, Whistlers, Langmuir), making it hard to track them all.

PlasmaBO.jl overcomes these challenges using a Matrix Eigenvalue Method (Xie [1], Xie and Xiao [2]):

Rational Approximation: We approximate the Plasma Dispersion Function $Z(\zeta)$ using a multi-pole (J-pole) expansion:
math Z(\zeta) \approx \sum_{j=1}^J \frac{b_j}{\zeta - c_j}
Linearization: This substitution allows us to transform the nonlinear dispersion relation into a linear matrix eigenvalue problem of the form:
math \mathcal{M} \mathbf{V} = ω \mathbf{V}
where $\mathcal{M}$ is a sparse matrix constructed from plasma parameters and auxiliary variables.
Global Solution: By computing the eigenvalues of $\mathcal{M}$, we obtain all roots $ω$ of the dispersion relation simultaneously.

Key Advantages:

No Initial Guesses Required: The solver finds all unstable and damped modes automatically.
Robostness: Effectively handles complex kinetic effects and arbitrary velocity distributions (via the BO-HH expansion).
Efficiency: Optimized for finding multiple modes in complex plasmas without manual intervention.

Features

Hermite-Hermite (BO-HH) expansion solver for arbitrary/analytic distributions
- Maxwellian / BiMaxwellian
Analytic Product-Bi-Kappa (BO-PBK) solver
Multi-fluid solver
Integration with ChargedParticles.jl for specifying different species, e.g., Maxwellian("O-18 3+", n, Tpara)

Usage Examples

Check out the ring beam instability example for detailed usage instructions, also see firehose instability example for using with arbitrary velocity distributions, BO-PBK example for using with kappa distributions (BO-PBK), cold plasma example for comparing kinetic and fluid solvers, dispersion surface tracking example for 2D scanning and mode tracking, and wave polarization and handedness example for analyzing wave modes in cold plasma.

Solvers

PlasmaBO.solve — Function

solve(species, B0, kx, kz, alg = BOHH; kw...)
solve(species, B0, ks, θs, alg = BOHH; kw...)

Convenience interface for solving a single wavevector (kx, kz) or dispersion scan over multiple wavevectors (ks, θs).

Keyword arguments are passed to the algorithm constructor. Defaults to BOHH solver.

When should you use each?

Use BO-PBK if the plasma is reasonably described by PBK and you care about κ-dependent instability thresholds, loss-cone physics, and you want speed + robustness.

Use BO-Arbitrary if the distribution is numerical (e.g., from spacecraft data) and cannot be parameterized analytically (e.g., has plateaus, shoulders, or holes). You are willing to trade efficiency for universality.

References

Xie et al. [3], Xie [4]

[1]: H. Xie. Efficient Framework for Solving Plasma Waves with Arbitrary Distributions. Physics of Plasmas 32, 60702 (2025).
[2]: H. Xie and Y. Xiao. PDRK: A General Kinetic Dispersion Relation Solver for Magnetized Plasma*. Plasma Science and Technology 18, 97 (2016).
[3]: H.-s. Xie, R. Denton, J.-s. Zhao and W. Liu. BO 2.0: Plasma Wave and Instability Analysis with Enhanced Polarization Calculations (Mar 2021), arXiv:2103.16014 [physics].
[4]: H.-s. Xie. PDRF: A General Dispersion Relation Solver for Magnetized Multi-Fluid Plasma. Computer Physics Communications 185, 670–675 (2014).
[5]: T. Umeda, S. Matsukiyo, T. Amano and Y. Miyoshi. A Numerical Electromagnetic Linear Dispersion Relation for Maxwellian Ring-Beam Velocity Distributions. Physics of Plasmas 19, 72107 (2012).
[6]: S. P. Gary, C. W. Smith, M. A. Lee, M. L. Goldstein and D. W. Forslund. Electromagnetic ion beam instabilities. The Physics of fluids 27, 1852–1862 (1984).
[7]: T. Cattaert, M. A. Hellberg and R. L. Mace. Oblique Propagation of Electromagnetic Waves in a Kappa-Maxwellian Plasma. Physics of Plasmas 14, 82111 (2007).
[8]: W. Bai and H. Xie. BO-PBK: A Comprehensive Solver for Dispersion Relations of Obliquely Propagating Waves in Magnetized Multi-Species Plasma with Anisotropic Loss-Cone Drift Product-Bi-Kappa Distribution (Dec 2025), arXiv:2512.06901 [physics].

API Reference

PlasmaBO.BOFluid
PlasmaBO.BOHH
PlasmaBO.BranchPoint
PlasmaBO.FluidSpecies
PlasmaBO.FluidSpecies
PlasmaBO.JPoleCoefficients
PlasmaBO.PBK.BOPBK
PlasmaBO.funIn
PlasmaBO.get_jpole_coefficients
PlasmaBO.handedness
PlasmaBO.hermite_H
PlasmaBO.hermite_a0_to_a
PlasmaBO.hermite_coefficients_matrix
PlasmaBO.hermite_expansion
PlasmaBO.polarization_ratio
PlasmaBO.solve
PlasmaBO.track

PlasmaBO.BOFluid — Type

BOFluid()

Dispersion solver using the multi-fluid electromagnetic matrix formulation.

source

PlasmaBO.BOHH — Type

BOHH(; N = 2, J = 8)

Dispersion solver using the Hermite-Hermite (HH) basis expansion.

N controls the truncation order of the cyclotron harmonic index. J controls the truncation order of the number of poles for Z-function approximation

This method transforms the dispersion relation into a matrix eigenvalue problem using J-pole approximation for the plasma dispersion function, allowing simultaneous computation of all wave modes.

source

PlasmaBO.BranchPoint — Type

BranchPoint{T}
BranchPoint(k, ω)

Initial point specification for dispersion branch tracking.

Fields

k: Wave vector at initial point
ω: Frequency at initial point
locator: Callable used to locate the seed eigenvalue at the initial point

Example

# Track unstable branch starting at k=0.03 with γ=0.1721
point1 = BranchPoint(0.03, 0.1721im)     # Auto: track by Im

# Track real frequency mode at k=0.5 with ωᵣ=2.5
point2 = BranchPoint(0.5, 2.5)  # Auto: track by Re

source

PlasmaBO.FluidSpecies — Type

FluidSpecies(n, Tz, Tp = Tz; vdz=0.0, gamma_z=1.0, gamma_p=1.0, particle=:p, ...)

source

PlasmaBO.FluidSpecies — Type

FluidSpecies{T}

Fluid species parameters for the multi-fluid dispersion relation solver.

Fields

q: Charge
m: Mass
n: Number density (m⁻³)
Tz: Parallel temperature (eV)
Tp: Perpendicular temperature (eV)
vdz: Parallel drift velocity (in units of c)
gamma_z: Parallel polytrope exponent (default: 1.0 for isothermal)
gamma_p: Perpendicular polytrope exponent (default: 1.0 for isothermal)

source

PlasmaBO.JPoleCoefficients — Type

JPoleCoefficients{T}

J-pole approximation coefficients for the plasma dispersion function. Z(ζ) ≈ Σⱼ bⱼ/(ζ - cⱼ)

source

PlasmaBO.funIn — Method

funIn(n)
funIn(T, n)

Normalization integral I_n = ∫ v^n exp(-v²) dv / √π.

Returns Γ((n+1)/2)/√π for even n, 0 for odd n.

source

PlasmaBO.get_jpole_coefficients — Method

get_jpole_coefficients(J=8)

Get J-pole approximation coefficients for the plasma dispersion function.

Supported values: J ∈ {4, 6, 8, 10, 12, 16, 20, 24, 28, 32}

source

PlasmaBO.handedness — Method

handedness(v; kx=0.0, kz=0.0, threshold=1e-5)
handedness(v, ω; kx=0.0, kz=0.0, threshold=1e-5)
handedness(v, ω, kx, kz; threshold=1e-5)

Determine the physical polarization handedness of the wave mode represented by the eigenvector v with complex frequency ω and wavevector (kx, kz) relative to the background magnetic field direction (+z direction).

We assume the wave phase convention is e^{i(k_z z - ω t)} with wave propagation along +z. This physical handedness correctly projects the electric field perpendicular to the wavevector k to support oblique modes, and factors in the sign of real(ω) to determine the real-time rotation sense of the physical field.

Returns:

:R (Right-handed/whistler-like/electron-sense)
:L (Left-handed/Alfvén-like/ion-sense)
:linear (Linearly polarized or near-linear within the threshold)

source

PlasmaBO.hermite_H — Method

hermite_H(n::Int, x)

Compute physicists' Hermite polynomial H_n(x) using recurrence relation.

Recurrence Relations

H_0(x) = 1
H_1(x) = 2x
H{n+1}(x) = 2x*Hn(x) - 2n*H_{n-1}(x)

Examples

hermite_H(0, 1.5)  # Returns 1.0
hermite_H(1, 1.5)  # Returns 3.0 (= 2*1.5)
hermite_H(2, 1.5)  # Returns 7.0 (= 2*1.5*3.0 - 2*1*1.0)

source

PlasmaBO.hermite_a0_to_a — Method

hermite_a0_to_a(a0lm)

Convert normalized Hermite basis coefficients: a0lm -> alm

Transforms expansion:

f(z,x) = Σ_{l,m} a0_{l,m} * ρ_l(z) * u_m(x)

to:

f(z,x) = Σ_{l,m} a_{l,m} * g_l(z) * h_m(x)

where:

ρl, um are normalized Hermite functions
g_l(z) ∝ z^l * exp(-z²/2)
h_m(x) ∝ x^m * exp(-x²/2)

source

PlasmaBO.hermite_coefficients_matrix — Method

hermite_coefficients_matrix(nmax)

Compute coefficient matrix for expanding Hermite polynomials in power basis given the maximum order nmax.

Returns matrix cHn[n+1, k+1] such that:

H_n(x) / √(2^n n! √π) = Σ_k cHn[n+1, k+1] * x^k / √(2^{n-k} n! √π)

This is used to convert between normalized Hermite basis and power-law basis.

Returns

cHn: (nmax+1) × (nmax+1) coefficient matrix

source

PlasmaBO.hermite_expansion — Method

hermite_expansion(
    fv::AbstractMatrix{T},
    vz, vx, vtz, vtp;
    Nz = 16, Nx = 16,
    dz = zero(T), dx = zero(T),
    atol = zero(T)
) where {T}

Compute Hermite expansion coefficients from gridded distribution function.

Expands arbitrary 2D distribution f(vparallel, vperp) as:

f(vz, vx) = Σ_{l=0}^{Nz} Σ_{m=0}^{Nx} a_{l,m} * ρ_l(vz) * u_m(vx)

where ρl and um are normalized Hermite basis functions. Nz, Nx are the maximum parallel and perpendicular Hermite indices, respectively.

Arguments

fv: Distribution function values on (vz, vx) grid
vz: Parallel velocity grid (1D or matches fv size)
vx: Perpendicular velocity grid (1D or matches fv size)
vtz: Parallel thermal velocity (default: 1.0)
vtp: Perpendicular thermal velocity (default: 1.0)
dz: Parallel drift velocity (default: 0.0)
dx: Perpendicular drift velocity (default: 0.0)
atol: Absolute tolerance for zeroing small coefficients

Returns

Named tuple with:

alm: (Nz+1) × (Nx+1) coefficient matrix in power-law basis
a0lm: (Nz+1) × (Nx+1) coefficient matrix in Hermite basis

Notes

Grid-quadrature noise in high-order coefficients is amplified enormously by the power-basis moment weights (absolute noise as small as 1e-21 at l ≈ 40 produces 1e-3 relative eigenvalue errors). Therefore, atol = 1e-10 is recommended whenever Nz or Nx exceeds ~16; genuine coefficients of smooth VDFs are far above this level.

Algorithm

Compute a0_{l,m} via numerical integration:

a0_{l,m} = ∫∫ f(z,x) ρ_l(z) u_m(x) dz dx * (2/(Lz*Lx))

Normalize by distribution integral
Convert to power-law basis: a{l,m} = hermitea0toa(a0_{l,m})

source

PlasmaBO.polarization_ratio — Method

polarization_ratio(v, kx, kz)
polarization_ratio(v; kx=0.0, kz=0.0)

Return the transverse polarization ratio P = Ey / E_perp of the eigenvector v, where E_perp is the transverse electric field component perpendicular to the wavevector k in the x-z plane: E_perp = Ex * cos(θ) - Ez * sin(θ).

If both kx and kz are zero (or omitted), it defaults to parallel propagation (kx = 0.0, kz = 1.0), where the ratio simplifies to Ey / Ex.

source

PlasmaBO.solve — Function

solve(species, B0, kx, kz, alg = BOHH; kw...)
solve(species, B0, ks, θs, alg = BOHH; kw...)

Convenience interface for solving a single wavevector (kx, kz) or dispersion scan over multiple wavevectors (ks, θs).

Keyword arguments are passed to the algorithm constructor. Defaults to BOHH solver.