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

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[1.2094722121124053e-5 + 0.0623764416971194im]

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, θ)...  13%|███                    |  ETA: 0:00:08
Solving dispersion (k, θ)...  26%|██████                 |  ETA: 0:00:07
Solving dispersion (k, θ)...  39%|█████████              |  ETA: 0:00:06
Solving dispersion (k, θ)...  52%|████████████           |  ETA: 0:00:04
Solving dispersion (k, θ)...  65%|███████████████        |  ETA: 0:00:03
Solving dispersion (k, θ)...  78%|██████████████████     |  ETA: 0:00:02
Solving dispersion (k, θ)...  91%|█████████████████████  |  ETA: 0:00:01
Solving dispersion (k, θ)... 100%|███████████████████████| Time: 0:00:09

using CairoMakie

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