Case: Firehose instability (kappa distribution)

This page demonstrates test firehose instability using Hermite-Hermite (HH) expansion for bi-kappa distribution. Here we first generate numerical samples of the distribution function using gen_fv2d. Then we use hermite_expansion to compute the Hermite expansion coefficients (this would also work for arbitrary distributions).

Simulation Parameters

• Background Magnetic Field: $B_0 = 0.1$ T
• Species:
  1. Electrons (Maxwellian)
     • Density: $n_e = 5 \times 10^{19}$ m$^{-3}$
     • Temperature: $T_e = 496.683$ eV
  2. Protons (Bi-Kappa)
     • Density: $n_p = 5 \times 10^{19}$ m$^{-3}$
     • Parallel Temperature: $T_{\parallel p} = 1986.734$ eV
     • Perpendicular Temperature: $T_{\perp p} = 993.367$ eV
     • Kappa Indices: $\kappa_{\parallel} = \kappa_{\perp} = 5.5$

using PlasmaBO

B0 = 0.1  # [Tesla]
θ = deg2rad(45)
n = 5.0e19

electron = Maxwellian(:e, n, 496.683)

κz = 5.5
κx = 5.5
proton = BiKappa2(5.e19, κz, κx, 1986.734, 993.367)
data = gen_fv2d(proton)

(vz = [-5.882268236282655e6 -5.882268236282655e6 … -5.882268236282655e6 -5.882268236282655e6; -5.852856895101242e6 -5.852856895101242e6 … -5.852856895101242e6 -5.852856895101242e6; … ; 5.82344555391983e6 5.82344555391983e6 … 5.82344555391983e6 5.82344555391983e6; 5.852856895101243e6 5.852856895101243e6 … 5.852856895101243e6 5.852856895101243e6], vx = [0.0 19729.727457323723 … 3.9262157640074207e6 3.9459454914647443e6; 0.0 19729.727457323723 … 3.9262157640074207e6 3.9459454914647443e6; … ; 0.0 19729.727457323723 … 3.9262157640074207e6 3.9459454914647443e6; 0.0 19729.727457323723 … 3.9262157640074207e6 3.9459454914647443e6], fv = [9.193937817053725e-27 9.166820152414298e-27 … 4.482658173520089e-35 4.214228270554484e-35; 9.779555794453612e-27 9.75071084035085e-27 … 4.7681860142764826e-35 4.482658173520089e-35; … ; 1.040537038788064e-26 1.0374679583761412e-26 … 5.073312387562397e-35 4.769513012461223e-35; 9.779555794453586e-27 9.750710840350824e-27 … 4.7681860142764703e-35 4.4826581735200773e-35], dvz = 29411.341181413278, dvx = 19729.727457323723, vtz = 588226.8236282655, vtx = 394594.54914647446, vdz = 0.0, vdx = 0.0)

alm = hermite_expansion(data.fv, data.vz, data.vx, data.vtz, data.vtx).alm
proton_param = HHSolverParam(proton, B0; alm = alm)

kn = 31.0613
k = kn / 4
wci = proton_param.wc
species = (proton_param, electron)

ωs = solve(species, B0, k .* sincos(θ)...; N = 2, J = 24)
ω_unstable = filter(ω -> isfinite(ω) && imag(ω) > 0.001 * wci, ωs)
println("Unstable modes (ω/ωci): ", ω_unstable ./ wci)

Unstable modes (ω/ωci): ComplexF64[3.627845807376187e-6 + 0.06238378093912506im]

Dispersion Curve Scan

julia> k_ranges = (0.05:0.02:0.5) .* kn;
julia> sol = solve(species, B0, k_ranges, θ; N = 2, J = 24);
Solving dispersion (k, θ)...  22%|█████                  |  ETA: 0:00:06
Solving dispersion (k, θ)...  65%|███████████████        |  ETA: 0:00:02
Solving dispersion (k, θ)...  91%|█████████████████████  |  ETA: 0:00:00
Solving dispersion (k, θ)... 100%|███████████████████████| Time: 0:00:04

using CairoMakie

let
    f, (ax, ax2) = plot(sol, kn, wci)
    ylims!(ax, [0, 1])
    ylims!(ax2, [-0.1, 0.1])
    f
end