# plasma_plots.theory.kinetic

*module*

Kinetic dispersion relations of unmagnetized plasmas: Langmuir waves and Landau damping,
ion-acoustic waves, beam-plasma, two-stream and bump-on-tail instabilities, and the Weibel
instability.

> **Units**
>
> Electrostatic results are normalized to the reference electrons: time in 1/ω_pe, length in the
> Debye length λ_De and velocity in the thermal speed v_the = √(T_e/m_e), so that
> ω_pe = v_the/λ_De = 1. Densities are relative to the reference electron density, charges in e and
> masses in m_e. A Maxwellian species of thermal speed v_th = √(T/m) drifting at u enters through
> ζ = (ω − k u)/(√2 |k| v_th). The electromagnetic Weibel relation ([`weibel()`][weibel]) uses c
> instead: k in ω_pe/c and speeds in c.
> 
> All frequencies are complex for ``exp(i(kx − ωt))``: a positive imaginary part is a growth rate,
> a negative one a damping rate. Functions take scalars or arrays that broadcast against each
> other, and return a Python complex for scalar input, else a complex numpy array.

**Examples**

```pycon
>>> omega = langmuir(0.5)
>>> round(omega.real, 6), round(omega.imag, 6)
(1.415662, -0.153359)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L1-L1)

## plasma_plots.theory.kinetic.PROTON_ELECTRON_MASS_RATIO

*attribute* · *module attribute*

```python
PROTON_ELECTRON_MASS_RATIO = 1836.15267343
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L34-L34)

## plasma_plots.theory.kinetic.Maxwellian

*class* · *dataclass*

```python
class Maxwellian
```

A drifting Maxwellian species, in the normalized units of this module.

f(v) = n/(√(2π) v_th) exp(−(v − u)²/(2 v_th²)). The fields may also be arrays that broadcast
against the wavenumbers.

**Attributes**

- `density` (`float`) — The density n_s, relative to the reference electron density. Default: ``1.0``.
- `charge` (`float`) — The charge q_s in units of e. Default: ``-1.0`` (electrons).
- `mass` (`float`) — The mass m_s in units of m_e. Default: ``1.0``.
- `thermal_speed` (`float`) — The thermal speed v_th,s = √(T_s/m_s) in units of v_the, positive. Default: ``1.0``.
- `drift` (`float`) — The drift speed u_s in units of v_the. Default: ``0.0``.

**Examples**

```pycon
>>> Maxwellian().plasma_frequency
1.0
>>> beam = Maxwellian(density=0.1, thermal_speed=0.5, drift=4.5)
>>> round(float(beam.plasma_frequency), 4)
0.3162
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L94-L177)

### plasma_plots.theory.kinetic.Maxwellian.charge

*attribute* · *class attribute* · *instance attribute*

```python
charge: float = -1.0
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L124-L124)

### plasma_plots.theory.kinetic.Maxwellian.density

*attribute* · *class attribute* · *instance attribute*

```python
density: float = 1.0
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L123-L123)

### plasma_plots.theory.kinetic.Maxwellian.drift

*attribute* · *class attribute* · *instance attribute*

```python
drift: float = 0.0
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L127-L127)

### plasma_plots.theory.kinetic.Maxwellian.mass

*attribute* · *class attribute* · *instance attribute*

```python
mass: float = 1.0
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L125-L125)

### plasma_plots.theory.kinetic.Maxwellian.thermal_speed

*attribute* · *class attribute* · *instance attribute*

```python
thermal_speed: float = 1.0
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L126-L126)

### plasma_plots.theory.kinetic.Maxwellian.plasma_frequency

*property*

```python
plasma_frequency
```

The plasma frequency ω_ps = √(n_s q_s²/m_s), in units of ω_pe.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L130-L134)

### plasma_plots.theory.kinetic.Maxwellian.ions

*method* · *classmethod*

```python
def ions(mass_ratio=PROTON_ELECTRON_MASS_RATIO, temperature_ratio=1.0, charge=1.0, drift=0.0)
```

Create quasi-neutral ions for the reference electrons.

**Parameters**

- `mass_ratio` (`float`) (default: `PROTON_ELECTRON_MASS_RATIO`) — m_i/m_e. Default: the proton mass ratio, 1836.15.
- `temperature_ratio` (`float`) (default: `1.0`) — T_e/T_i. Default: ``1.0``.
- `charge` (`float`) (default: `1.0`) — The charge number Z; the density is 1/Z. Default: ``1.0``.
- `drift` (`float`) (default: `0.0`) — The drift speed in units of v_the. Default: ``0.0``.

**Returns**

- (`Maxwellian`) — Ions with thermal speed √(T_i/m_i) = 1/√(mass_ratio · temperature_ratio) in units of v_the.

**Examples**

```pycon
>>> ions = Maxwellian.ions(mass_ratio=100.0, temperature_ratio=4.0)
>>> ions.thermal_speed, round(float(ions.plasma_frequency), 4)
(0.05, 0.1)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L136-L177)

## plasma_plots.theory.kinetic.beam_plasma_cold

*function*

```python
def beam_plasma_cold(k, beam_speed, beam_density=0.1, plasma_density=1.0, all_roots=False)
```

Compute the cold beam-plasma frequencies: the roots of a quartic.

A cold electron beam of density n_b at speed v_b through a cold plasma of density n_p
(immobile ions): 1 = n_p/ω² + n_b/(ω − k v_b)², i.e.
ω⁴ − 2kv_b ω³ + (k²v_b² − n_p − n_b) ω² + 2kv_b n_p ω − n_p k²v_b² = 0. For n_b ≪ n_p the
maximum growth is γ ≈ (√3/2)(n_b/(2 n_p))^(1/3) √n_p at k v_b ≈ √n_p.

**Parameters**

- `k` (`float or array_like`) — The wavenumber in units of 1/λ_De (any length unit L works, with speeds in ω_pe L).
- `beam_speed` (`float or array_like`) — v_b.
- `beam_density` (`float or array_like`) (default: `0.1`) — n_b. Default: ``0.1``.
- `plasma_density` (`float or array_like`) (default: `1.0`) — n_p. Default: ``1.0``.
- `all_roots` (`bool`) (default: `False`) — Return all four roots, sorted by real part. Default: ``False``.

**Returns**

- (`complex or numpy.ndarray`) — The root with the largest imaginary part where it is unstable, else the real root closest to the slow beam mode k v_b − √n_b. With ``all_roots``, an extra last axis of length 4.

> **References**
>
> C. K. Birdsall and A. B. Langdon, Plasma Physics via Computer Simulation (IOP, 1991),
> sec. 5.10.
> R. J. Briggs, Electron-Stream Interaction with Plasmas (MIT Press, 1964).

**Examples**

```pycon
>>> omega = beam_plasma_cold(1.0, beam_speed=1.0, beam_density=0.001)
>>> round(omega.real, 4), round(omega.imag, 4)
(0.9588, 0.066)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L770-L829)

## plasma_plots.theory.kinetic.bohm_gross

*function*

```python
def bohm_gross(k)
```

Compute the Bohm–Gross frequency ω = √(1 + 3k²) of Langmuir waves.

The fluid (adiabatic, γ = 3) Langmuir frequency, the real part of [`langmuir()`][plasma_plots.theory.kinetic.langmuir] for
k λ_De ≪ 1.

**Parameters**

- `k` (`float or array_like`) — The wavenumber in units of 1/λ_De.

**Returns**

- (`complex or numpy.ndarray`) — ω in units of ω_pe (real-valued).

> **References**
>
> D. Bohm and E. P. Gross, "Theory of plasma oscillations. A. Origin of medium-like
> behavior", Phys. Rev. 75, 1851 (1949).

**Examples**

```pycon
>>> bohm_gross(0.0)
(1+0j)
>>> round(bohm_gross(0.5).real, 4)
1.3229
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L485-L513)

## plasma_plots.theory.kinetic.bump_on_tail

*function*

```python
def bump_on_tail(k, beam_density=0.1, beam_speed=4.5, beam_thermal_speed=0.5, bulk_density=None)
```

Compute the most unstable kinetic root of a bump-on-tail distribution.

A Maxwellian bulk (density n_0, thermal speed 1, at rest) plus a Maxwellian beam (n_b, v_b,
v_th,b) over immobile ions. Newton is started from the Bohm–Gross frequency and from phase
speeds across the beam (k·v for v from v_b − 3 v_th,b to v_b + v_th,b), and the root with the
largest imaginary part is returned: the growing Langmuir/beam mode where the bump makes the
distribution unstable (phase speed on its rising flank), the least damped root found where
not. For k < 0 the roots are those of |k| mirrored, ω → −conj ω.

**Parameters**

- `k` (`float or array_like`) — The wavenumber in units of 1/λ_De, nonzero.
- `beam_density` (`float or array_like`) (default: `0.1`) — n_b. Default: ``0.1``.
- `beam_speed` (`float or array_like`) (default: `4.5`) — v_b in units of v_the. Default: ``4.5``.
- `beam_thermal_speed` (`float or array_like`) (default: `0.5`) — v_th,b in units of v_the, positive. Default: ``0.5``.
- `bulk_density` (`float or array_like`) (default: `None`) — n_0. Default: ``1 − n_b``, so the total electron density is 1.

**Returns**

- (`complex or numpy.ndarray`) — ω in units of ω_pe; ``nan`` if no guess converged.

> **References**
>
> T. M. O'Neil and J. H. Malmberg, "Transition of the dispersion roots from beam-type to
> Landau-type solutions", Phys. Fluids 11, 1754 (1968).

**Examples**

```pycon
>>> omega = bump_on_tail(0.3)
>>> round(omega.real, 4), round(omega.imag, 4)
(1.0012, 0.1981)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L942-L1010)

## plasma_plots.theory.kinetic.electrostatic_dielectric

*function*

```python
def electrostatic_dielectric(omega, k, species=None, derivative=0)
```

Compute the electrostatic dielectric function ε(ω, k) = 1 + Σ_s χ_s of Maxwellian species.

Its zeros are the electrostatic (Langmuir, ion-acoustic, beam) modes.

**Parameters**

- `omega` (`complex or array_like`) — The complex frequency in units of ω_pe.
- `k` (`float or array_like`) — The wavenumber in units of 1/λ_De, nonzero.
- `species` (`Maxwellian or sequence of Maxwellian`) (default: `None`) — The mobile species. Default: the reference electrons, ``Maxwellian()``, with immobile ions.
- `derivative` (`(0, 1)`) (default: `0`) — 0 for ε, 1 for ∂ε/∂ω. Default: ``0``.

**Returns**

- (`complex or numpy.ndarray`) — ε or ∂ε/∂ω, broadcast over ``omega``, ``k`` and the species fields.

> **See Also**
>
> [`susceptibility()`][plasma_plots.theory.kinetic.susceptibility] : One species' contribution χ_s.

**Examples**

```pycon
>>> abs(electrostatic_dielectric(langmuir(0.3), 0.3)) < 1e-10
True
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L236-L274)

## plasma_plots.theory.kinetic.ion_acoustic

*function*

```python
def ion_acoustic(k, temperature_ratio=10.0, mass_ratio=PROTON_ELECTRON_MASS_RATIO)
```

Compute the kinetic ion-acoustic root of Maxwellian electrons and ions.

The zero of ε = 1 + χ_e + χ_i with the reference electrons and ions of
[`Maxwellian.ions()`][plasma_plots.theory.kinetic.Maxwellian.ions], damped by electron and ion Landau damping. It is continued from
k λ_De = 0.01 and T_e/T_i = max(τ, 30), seeded with the fluid frequency and the weak damping
γ/ω = −√(π/8) (√(m_e/m_i) + τ^(3/2) exp(−τ/2 − 3/2)) (Chen, eq. 7.144).

**Parameters**

- `k` (`float or array_like`) — The wavenumber in units of 1/λ_De, positive (ω(−k) = −conj ω(k)).
- `temperature_ratio` (`float or array_like`) (default: `10.0`) — τ = T_e/T_i. Default: ``10.0``.
- `mass_ratio` (`float or array_like`) (default: `PROTON_ELECTRON_MASS_RATIO`) — μ = m_i/m_e. Default: the proton mass ratio, 1836.15.

**Returns**

- (`complex or numpy.ndarray`) — ω in units of ω_pe; ``nan`` where the root was lost.

> **See Also**
>
> [`ion_acoustic_fluid()`][plasma_plots.theory.kinetic.ion_acoustic_fluid] : The fluid estimate.

> **References**
>
> F. F. Chen, Introduction to Plasma Physics and Controlled Fusion, 3rd ed. (Springer, 2016),
> sec. 7.6.
> B. D. Fried and R. W. Gould, "Longitudinal ion oscillations in a hot plasma", Phys. Fluids 4,
> 139 (1961).

**Examples**

```pycon
>>> omega = ion_acoustic(0.1, temperature_ratio=10.0, mass_ratio=100.0)
>>> round(omega.real, 5), round(omega.imag, 5)
(0.01178, -0.00076)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L655-L716)

## plasma_plots.theory.kinetic.ion_acoustic_fluid

*function*

```python
def ion_acoustic_fluid(k, temperature_ratio=10.0, mass_ratio=PROTON_ELECTRON_MASS_RATIO, adiabatic_index=3.0)
```

Compute the fluid ion-acoustic frequency with Boltzmann electrons and adiabatic ions.

ω² = k² (T_e/(1 + k² λ_De²) + γ_i T_i)/m_i, i.e. in normalized units
ω = k √((1/(1 + k²) + γ_i/τ)/μ) with τ = T_e/T_i and μ = m_i/m_e; for k λ_De ≪ 1,
ω = k c_s with c_s = √((T_e + γ_i T_i)/m_i).

**Parameters**

- `k` (`float or array_like`) — The wavenumber in units of 1/λ_De.
- `temperature_ratio` (`float or array_like`) (default: `10.0`) — τ = T_e/T_i. Default: ``10.0``.
- `mass_ratio` (`float or array_like`) (default: `PROTON_ELECTRON_MASS_RATIO`) — μ = m_i/m_e. Default: the proton mass ratio, 1836.15.
- `adiabatic_index` (`float`) (default: `3.0`) — γ_i of the ions. Default: ``3.0`` (one-dimensional adiabatic compression).

**Returns**

- (`complex or numpy.ndarray`) — ω in units of ω_pe (real-valued).

> **References**
>
> F. F. Chen, Introduction to Plasma Physics and Controlled Fusion, 3rd ed. (Springer, 2016),
> sec. 4.6.

**Examples**

```pycon
>>> round(
...     ion_acoustic_fluid(0.1, temperature_ratio=10.0, mass_ratio=100.0).real,
...     5,
... )
0.01136
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L604-L652)

## plasma_plots.theory.kinetic.landau_damping_weak

*function*

```python
def landau_damping_weak(k)
```

Compute the weak-damping approximation of Langmuir waves: Bohm–Gross plus Landau damping.

ω = √(1 + 3k²) − i √(π/8) k⁻³ exp(−1/(2k²) − 3/2), the textbook expansion for k λ_De ≪ 1
(Chen, eq. 7.133). It overestimates the damping at moderate k (by ~60 % at k = 0.3); use
[`langmuir()`][plasma_plots.theory.kinetic.langmuir] for the exact root.

**Parameters**

- `k` (`float or array_like`) — The wavenumber in units of 1/λ_De, nonzero.

**Returns**

- (`complex or numpy.ndarray`) — ω in units of ω_pe.

> **References**
>
> F. F. Chen, Introduction to Plasma Physics and Controlled Fusion, 3rd ed. (Springer, 2016),
> sec. 7.5.

**Examples**

```pycon
>>> omega = landau_damping_weak(0.2)
>>> round(omega.real, 4), round(omega.imag, 7)
(1.0583, -6.51e-05)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L516-L548)

## plasma_plots.theory.kinetic.langmuir

*function*

```python
def langmuir(k)
```

Compute the exact kinetic Langmuir root: the least-damped zero of ε(ω, k) near Bohm–Gross.

Electrons (Maxwellian, v_th = 1) with immobile ions. The root is continued in k from
k λ_De = 0.1, where the weak-damping approximation is accurate.

**Parameters**

- `k` (`float or array_like`) — The wavenumber in units of 1/λ_De.

**Returns**

- (`complex or numpy.ndarray`) — ω in units of ω_pe, with ω(−k) = ω(k) and ω(0) = 1; ``nan`` where the root was lost.

> **See Also**
>
> [`landau_damping_weak()`][plasma_plots.theory.kinetic.landau_damping_weak] : The small-k approximation.

> **References**
>
> L. D. Landau, "On the vibrations of the electronic plasma", J. Phys. USSR 10, 25 (1946).
> J. Canosa, "Numerical solution of Landau's dispersion equation", J. Comput. Phys. 13, 158
> (1973).

**Examples**

```pycon
>>> omega = langmuir([0.3, 0.4, 0.5, 1.0])
>>> np.round(omega, 4)
array([1.1598-0.0126j, 1.2851-0.0661j, 1.4157-0.1534j, 2.0459-0.8513j])
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L551-L598)

## plasma_plots.theory.kinetic.maximum_growth

*function*

```python
def maximum_growth(function, k_range, samples=64, tol=1e-08)
```

Find the wavenumber of maximum growth rate of a dispersion relation.

``function`` is sampled at ``samples`` points across ``k_range``, and the largest Im ω is
refined by golden-section search between the neighbours of the best sample.

**Parameters**

- `function` (`callable`) — ``function(k)`` returning the complex ω for an array of wavenumbers, e.g. [`bump_on_tail()`][plasma_plots.theory.kinetic.bump_on_tail] or ``lambda k: two_stream(k, 3.0, 0.3)``.
- `k_range` (`(float, float)`) — The interval of wavenumbers searched.
- `samples` (`int`) (default: `64`) — The number of samples of the initial scan. Default: ``64``.
- `tol` (`float`) (default: `1e-08`) — The absolute tolerance in k of the refinement. Default: ``1e-8``.

**Returns**

- (`tuple of (float, complex)`) — The wavenumber of maximum growth and ω there.

**Raises**

- `ValueError` — If ``function`` returns no finite value on the scan.

**Examples**

```pycon
>>> k, omega = maximum_growth(lambda k: two_stream_cold(k, 1.0), (0.01, 1.4))
>>> round(k, 5), round(omega.imag, 6)  # √(3/8), √0.5/2
(0.61237, 0.353553)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L1013-L1071)

## plasma_plots.theory.kinetic.solve_dispersion

*function*

```python
def solve_dispersion(function, k, guess, derivative=None, continuation=True, tol=1e-11, maxiter=60)
```

Find a complex root ω(k) of a dispersion relation for every wavenumber.

Damped Newton iteration in ω, with the analytic derivative when one is given and a
central difference otherwise. With ``continuation``, the wavenumbers are solved in the given
order and each root seeds the next (linearly extrapolated from the last two), which follows one
branch along a k-grid; otherwise every k starts from ``guess`` and all are solved at once.

**Parameters**

- `function` (`callable`) — ``function(omega, k)``, analytic in ω, whose zero is sought, e.g. ``lambda w, k: electrostatic_dielectric(w, k, species)``. It must broadcast over arrays.
- `k` (`float or array_like`) — The wavenumbers; one-dimensional with ``continuation``, else any shape.
- `guess` (`complex or array_like`) — The starting frequency: a scalar for the first k with ``continuation``, else anything that broadcasts against ``k``.
- `derivative` (`callable`) (default: `None`) — ``derivative(omega, k)``, the ω-derivative of ``function``. Default: a central difference.
- `continuation` (`bool`) (default: `True`) — Continue the root along ``k``. Default: ``True``.
- `tol` (`float`) (default: `1e-11`) — Converged once the Newton step is below ``tol · |ω|``. Default: ``1e-11``.
- `maxiter` (`int`) (default: `60`) — The maximum number of Newton steps per wavenumber. Default: ``60``.

**Returns**

- (`complex or numpy.ndarray`) — ω for every k, ``nan`` where Newton did not converge (no exception is raised); with ``continuation`` the next wavenumber is then seeded by the last converged root.

**Raises**

- `ValueError` — With ``continuation``, if ``k`` has more than one dimension or ``guess`` is not a scalar.

**Examples**

The Langmuir branch followed from k = 0.2 to 0.5:

```pycon
>>> k = np.linspace(0.2, 0.5, 7)
>>> omega = solve_dispersion(electrostatic_dielectric, k, guess=1.06)
>>> round(complex(omega[-1]).real, 4), round(complex(omega[-1]).imag, 4)
(1.4157, -0.1534)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L361-L453)

## plasma_plots.theory.kinetic.susceptibility

*function*

```python
def susceptibility(omega, k, species, derivative=0)
```

Compute the electrostatic susceptibility χ_s(ω, k) of one Maxwellian species.

χ_s = (ω_ps²/(k² v_th,s²)) [1 + ζ_s Z(ζ_s)], with ζ_s = (ω − k u_s)/(√2 |k| v_th,s) and Z the
plasma dispersion function (Landau's continuation, valid for any complex ω). For a negative k
the drift enters as k u_s, so χ_s(ω, −k) is χ_s(ω, k) with the drift reversed.

**Parameters**

- `omega` (`complex or array_like`) — The complex frequency in units of ω_pe.
- `k` (`float or array_like`) — The wavenumber in units of 1/λ_De, nonzero.
- `species` (`Maxwellian`) — The species.
- `derivative` (`(0, 1)`) (default: `0`) — 0 for χ_s, 1 for ∂χ_s/∂ω. Default: ``0``.

**Returns**

- (`complex or numpy.ndarray`) — χ_s or ∂χ_s/∂ω, broadcast over ``omega``, ``k`` and the species fields.

**Raises**

- `ValueError` — If ``derivative`` is not 0 or 1.

> **References**
>
> B. D. Fried and S. D. Conte, The Plasma Dispersion Function (Academic Press, 1961).
> T. H. Stix, Waves in Plasmas (AIP, 1992), ch. 8.

**Examples**

The cold limit χ → −ω_p²/ω²:

```pycon
>>> round(susceptibility(10.0, 0.01, Maxwellian()).real, 6)
-0.01
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L183-L233)

## plasma_plots.theory.kinetic.two_stream

*function*

```python
def two_stream(k, beam_speed, thermal_speed, beam_density=0.5)
```

Compute the kinetic purely growing root of two warm counter-streaming electron beams.

Two Maxwellian electron beams of density n_b and thermal speed v_t at ±v_b over immobile
ions. By symmetry ε(iγ, k) is real, and the two-stream instability is purely growing,
ω = iγ: it exists where ε(0, k) < 0, and γ is the largest zero of ε(iγ, k), found by a scan
and bisection. Where the beams are stable, the least damped root reached by Newton from
phase speeds between 0 and v_b is returned instead, with Re ω ≥ 0 (−conj ω is a root too):
past the marginal wavenumber the growing root iγ meets its mirror −iγ near γ = 0 and turns
into a weakly damped oscillating pair, so γ(k) stays continuous. For v_t → 0 the growth
approaches [`two_stream_cold()`][plasma_plots.theory.kinetic.two_stream_cold].

**Parameters**

- `k` (`float or array_like`) — The wavenumber in units of 1/λ_De (ω(−k) = ω(k)).
- `beam_speed` (`float or array_like`) — v_b in units of v_the.
- `thermal_speed` (`float or array_like`) — v_t of each beam in units of v_the, positive.
- `beam_density` (`float or array_like`) (default: `0.5`) — n_b, each beam's density. Default: ``0.5``, so the total density is 1.

**Returns**

- (`complex or numpy.ndarray`) — ω in units of ω_pe: iγ where unstable; ``nan`` where stable and Newton found no root.

> **See Also**
>
> [`two_stream_cold()`][plasma_plots.theory.kinetic.two_stream_cold] : The cold limit.

> **References**
>
> T. H. Stix, Waves in Plasmas (AIP, 1992), sec. 8.11.
> C. K. Birdsall and A. B. Langdon, Plasma Physics via Computer Simulation (IOP, 1991),
> sec. 5.10.

**Examples**

```pycon
>>> omega = two_stream(0.2, beam_speed=3.0, thermal_speed=0.3)
>>> round(omega.real, 6), round(omega.imag, 4)
(0.0, 0.3491)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L832-L939)

## plasma_plots.theory.kinetic.two_stream_cold

*function*

```python
def two_stream_cold(k, beam_speed, beam_density=0.5, all_roots=False)
```

Compute the cold symmetric two-stream frequencies in closed form.

Two cold electron beams of density n_b each at ±v_b over immobile ions:
1 = n_b/(ω − k v_b)² + n_b/(ω + k v_b)², a quadratic in ω² with the roots
ω² = k²v_b² + n_b ± √(4 n_b k² v_b² + n_b²). The "slow" root ω² < 0, i.e. the purely growing
ω = iγ, exists for k² v_b² < 2 n_b; the maximum growth is γ = √n_b/2 at k² v_b² = 3 n_b/4.

**Parameters**

- `k` (`float or array_like`) — The wavenumber in units of 1/λ_De (any length unit L works, with speeds in ω_pe L).
- `beam_speed` (`float or array_like`) — v_b, each beam's speed.
- `beam_density` (`float or array_like`) (default: `0.5`) — n_b, each beam's density (its ω_pb²). Default: ``0.5``, so the total density is 1.
- `all_roots` (`bool`) (default: `False`) — Return all four roots instead of the slow one. Default: ``False``.

**Returns**

- (`complex or numpy.ndarray`) — The slow root √(ω₋²) with Im ω ≥ 0: iγ where unstable, the real slow beam mode where stable. With ``all_roots``, an extra last axis of length 4: (√(ω₊²), −√(ω₊²), √(ω₋²), −√(ω₋²)).

> **References**
>
> C. K. Birdsall and A. B. Langdon, Plasma Physics via Computer Simulation (IOP, 1991),
> sec. 5.10.

**Examples**

```pycon
>>> omega = two_stream_cold(np.sqrt(3 / 8), beam_speed=1.0)
>>> round(omega.real, 6), round(omega.imag, 6)
(0.0, 0.353553)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L722-L767)

## plasma_plots.theory.kinetic.weibel

*function*

```python
def weibel(k, anisotropy, parallel_thermal_speed)
```

Compute the purely growing (or damped) root of the electron Weibel instability.

Transverse electromagnetic waves along k ∥ z in bi-Maxwellian electrons (T⊥ across k, T∥
along k) with immobile ions:

ω² − k²c² − ω_pe² + ω_pe² A [1 + ζ Z(ζ)] = 0,  A = T⊥/T∥,  ζ = ω/(√2 k v_∥),

with v_∥ = √(T∥/m_e). It reduces to light waves ω² = ω_pe² + k²c² for T → 0. The root is
ω = iγ, with γ > 0 exactly for k²c² < (A − 1) ω_pe² and γ < 0 (a damped, non-oscillating
mode) beyond; γ is found by bisection of the real function D(iγ), which decreases in γ. Near
the cutoff γ ≈ √(2/π) k v_∥ (A − 1 − k²c²/ω_pe²)/A.

**Units:** unlike the rest of this module, k in ω_pe/c, ω in ω_pe and speeds in c.

**Parameters**

- `k` (`float or array_like`) — The wavenumber in units of ω_pe/c (ω(−k) = ω(k), ω(0) = 0).
- `anisotropy` (`float or array_like`) — A = T⊥/T∥, positive; unstable for A > 1.
- `parallel_thermal_speed` (`float or array_like`) — v_∥ = √(T∥/m_e) in units of c, positive.

**Returns**

- (`complex or numpy.ndarray`) — ω = iγ in units of ω_pe.

> **References**
>
> E. S. Weibel, "Spontaneously growing transverse waves in a plasma due to an anisotropic
> velocity distribution", Phys. Rev. Lett. 2, 83 (1959).
> R. C. Davidson, D. A. Hammer, I. Haber and C. E. Wagner, "Nonlinear development of
> electromagnetic instabilities in anisotropic plasmas", Phys. Fluids 15, 317 (1972).
> N. A. Krall and A. W. Trivelpiece, Principles of Plasma Physics (McGraw-Hill, 1973), sec. 9.10.

**Examples**

```pycon
>>> round(weibel(1.0, anisotropy=4.0, parallel_thermal_speed=0.1).imag, 4)
0.061
>>> weibel(
...     np.sqrt(3.0), anisotropy=4.0, parallel_thermal_speed=0.1
... ).imag < 1e-12
True
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/theory/kinetic.py#L1077-L1148)