GitHub

PlasmaBO

DOI version

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:

  1. Transcendental Equations: The dispersion tensor involves the Plasma Dispersion Function $Z(\zeta)$, making the equation $D(ω) = 0$ transcendental and highly nonlinear.
  2. 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.
  3. 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]):

  1. 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}

  2. 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.

  3. 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, and dispersion surface tracking example for 2D scanning and mode tracking.

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.

See also: BOHH, BOPBK, BOFluid

source

BO-PBK (BOPBK) is an analytic, distribution-aware eigen-solver optimized for kappa plasmas, whereas BO-Arbitrary (BOHH) is a universal but numerically heavier framework that approximates any distribution at the cost of efficiency and low-κ accuracy.

PlasmaBO.BOPBK — Type
BOPBK(; N = 2)

Dispersion solver using the PBK matrix formulation.

N controls the truncation order of the cyclotron harmonic index used to build the dispersion matrix.

source
PlasmaBO.BOHH — Type
BOHH(; N = 2, J = 8)

Dispersion solver using the Hermite-Hankel (HH) matrix formulation.

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.BOFluid — Type
BOFluid()

Dispersion solver using the multi-fluid electromagnetic matrix formulation.

source

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]
[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]

API Reference

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-Hankel (HH) matrix formulation.

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.BOPBK — Type
BOPBK(; N = 2)

Dispersion solver using the PBK matrix formulation.

N controls the truncation order of the cyclotron harmonic index used to build the dispersion matrix.

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{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.funAn — Method
funAn(n, a, d, m, p)

Perpendicular integral for Jn²: ∫ Jn²(ay) exp(-(y-d)²) (y-d)^m y^p dy.

source
PlasmaBO.funBn — Method
funBn(n, a, d, m, p)

Perpendicular integral for Jn·Jn': ∫ Jn(ay) Jn'(ay) exp(-(y-d)²) (y-d)^m y^p dy.

source
PlasmaBO.funCn — Method
funCn(n, a, d, m, p)

Perpendicular integral for (Jn')²: ∫ [Jn'(ay)]² exp(-(y-d)²) (y-d)^m y^p dy.

source
PlasmaBO.funIn — Method
funIn(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 — Function
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.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)
) 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)

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

Algorithm

  1. Compute a0_{l,m} via numerical integration:

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

  2. Normalize by distribution integral

  3. Convert to power-law basis: a{l,m} = hermitea0toa(a0_{l,m})

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.

See also: BOHH, BOPBK, BOFluid

source
PlasmaBO.track — Method
track(solution, point)

Track a single dispersion branch across parameter space from initial point.

The algorithm:

  1. Start at the k-point nearest to the initial point
  2. Track bidirectionally using interpolation to predict ω at each k
  3. Select nearest eigenvalue to prediction

Example

initial = BranchPoint(0.03, 0.1721im)
k_branch, ω_branch = track(solution, initial)

See also: BranchPoint

source