2.1. Kinetic theory

The Weibel instability arises from initially considering a Maxwellian plasma with anisotropic temperatures, wherein the initial equilibrium plasma VDF is assumed to be of the form

(2.1) \begin{equation} f_{s,0} = \left (\frac {m_s}{2\pi k_B}\right )^{3/2}\frac {1}{\sqrt {T_{s,\parallel }T_{s,\perp }^2}}\exp \left [-\frac {m_sv_x^2}{2k_BT_{s,\parallel }}-\frac {m_s\left (v_y^2+v_z^{2}\right )}{2k_BT_{s,\perp }}\right ], \end{equation}

where $m_s$ is the mass of species $s$ , $k_B$ is the Boltzmann constant, $v_i$ is the species velocity in the $i{\text {th}}$ direction, i.e. $i=\{1,2,3\}=\{x,y,z\}$ , and we have assumed, without loss of generality, that the plasma species has temperature $T_{s,\parallel }=T_{s,x}$ in the $x$ ( $\parallel$ ) direction, different from the temperature $T_{s,\perp }$ , where $\perp$ refers to the directions perpendicular to $x$ (i.e. $y$ and $z$ ), and that we are in the reference frame where the species is at rest (i.e. there is no bulk velocity).

The electromagnetic perturbation analysis of this distribution function is presented in the literature (Weibel Reference Weibel1959) and textbooks (Krall et al. Reference Krall, Trivelpiece and Gross1973), and is reproduced here for completeness. The dispersion relation is found by finding the zeros of the dielectric function (Krall et al. Reference Krall, Trivelpiece and Gross1973) for perturbations in the $x$ direction:

(2.2) \begin{align} D_{\text {K}}(k_x,\omega ) &= 1-\frac {k_x^2c^2}{\omega ^2}-\sum _s\frac {\omega _{s,p}^2}{\omega ^2}\left [1+\int \frac {v_{\perp}^2\left (\partial f_{s,0}/\partial v_x\right )}{v_x-\omega /k_x}\text{d}\textbf {v}\right ] \nonumber \\ &= 1-\frac {k_x^2c^2}{\omega ^2}-\sum _s\frac {\omega _{s,p}^2}{\omega ^2}\left [1+\frac {T_{s,\perp }}{2T_{s,\parallel }}Z'\left (\xi _{s}\right )\right ]=0, \end{align}

where $\omega _{s,p} = \sqrt {n_se^2/(m_s\epsilon _0)}$ is the species plasma frequency, where $n_s$ is the species number density, $e$ is the fundamental charge and $\epsilon _0$ is the permittivity of free space, $\xi _{s} = \omega /\left (k_x\sqrt {2T_{s,\parallel }/m_s}\right )$ , and $Z'(\xi )$ is the first derivative of the plasma dispersion function (Fried & Conte Reference Fried and Conte2015). Note that $\xi$ may be interpreted as the wave speed normalised by the thermal velocity; or, in the case for imaginary $\omega$ , a normalised growth rate. It is oftentimes useful to employ the following asymptotic expansions for $Z'(\xi )$ :

(2.3) \begin{equation} Z'(\xi )\simeq \begin{cases} {-2i\sqrt {\pi }\xi e^{-\xi ^2}+\frac {1}{\xi ^2}+\frac {3}{2\xi ^4}+\frac {15}{4\xi ^6}+\cdots } \quad \text {for } |\xi |\gg 1, \\ \\ {-2i\sqrt {\pi }\xi e^{-\xi ^2}-2+4\xi ^2-\frac {8\xi ^4}{3}+\cdots } \quad\;\; \text {for } |\xi |\ll 1. \end{cases} \end{equation}

2.2. Ten-moment fluid – Maxwell theory

To determine the fluid dispersion relation, we begin with the collisionless 10-moment fluid equations for a charged fluid using the primitive variables (Hakim Reference Hakim2008; Kuldinow et al. Reference Kuldinow, Yamashita, Mansour and Hara2024b ):

(2.4) \begin{equation} \frac {\text{d}\rho }{\text{d}t}+\rho \frac {\partial u_i}{\partial x_i}=0, \end{equation}

(2.5) \begin{equation} \frac {\text{d}u_i}{\text{d}t}+\frac {1}{\rho }\frac {\partial p_{ij}}{\partial x_j}=\frac {q}{m}\left (E_i+\epsilon _{ijk}u_jB_k\right ), \end{equation}

(2.6) \begin{equation} \frac {\text{d}p_{ij}}{\text{d}t} +p_{ij}\frac {\partial u_k}{\partial x_k} +p_{jk}\frac {\partial u_i}{\partial x_k} +p_{ki}\frac {\partial u_j}{\partial x_k} +\frac {\partial q_{ijk}}{\partial x_k} =\frac {q}{m}\left (\epsilon _{ikm}p_{jk}B_m+\epsilon _{jkm}p_{ik}B_m\right ), \end{equation}

where $\frac {\text{d}}{\text{d}t} = \frac {\partial }{\partial t}+u_i\frac {\partial }{\partial x_i}$ is the Lagrangian derivative, $\rho$ is the species mass density, $u_i$ is the species bulk velocity, $p_{ij}$ is the species pressure tensor, $q_{ijk}$ is the species heat flux tensor, $E_i$ and $B_i$ are the electric and magnetic fields, respectively, the species has charge $q$ and mass $m$ , and $\epsilon _{ijk}$ is the Levi–Civita symbol, employing the Einstein summation convention.

Note that conventional 5-moment models, which solve for $\rho, u_i$ and a total energy $\mathcal {E}$ , require closure for the isotropic pressure, $p$ , and deviatoric stress tensor, $\tau _{ij}$ . The pressure is usually closed by the equation of state, which is usually described using a polytropic index, and the stress tensor is typically described using a constitutive model (e.g. viscosity and shear). However, the 10-moment model explicitly formulates the pressure tensor evolution, which eliminates the use of an equation of state for the pressure and now only requires closure in the heat flux (third velocity moment).

For the derivation of a Weibel-type instability, we will assume a one-dimensional slab geometry in the $x$ direction, i.e. $\partial _y = \partial _z = 0$ , and we will only need equations for $u_y$ and $p_{xy}$ , which will lead to only to field components $B_z$ and $E_y$ . It is to be noted that perturbations to the density, $\rho$ , longitudinal velocity, $u_x$ , and diagonal pressure elements, $\{p_{xx},p_{yy},p_{zz}\}$ , only couple to the growth of the Weibel instability through second-order terms. For linear perturbations, only the transport equations for $u_y$ and $p_{xy}$ need to be considered, which can be written from (2.5) and (2.6) as

(2.7) \begin{equation} \frac {\text{d}u_y}{\text{d}t}+\frac {1}{\rho }\frac {\partial p_{xy}}{\partial x} = \frac {q}{m}\left (E_y+u_zB_x-u_xB_z\right ), \end{equation}

(2.8) \begin{equation} \frac {\text{d}p_{xy}}{\text{d}t} +2p_{xy}\frac {\partial u_x}{\partial x} +p_{xx}\frac {\partial u_y}{\partial x} +\frac {\partial q_{xxy}}{\partial x} =\frac {q}{m}\left [p_{xz}B_x+\left (p_{yy}-p_{xx}\right )B_z-p_{yz}B_y\right ]. \end{equation}

These equations must now be coupled with the Maxwell equations for the evolution of the electric and magnetic fields. In particular, we take linear perturbations of Faraday’s law and Ampère’s law for the magnetic field in only the $z$ direction:

(2.9) \begin{equation} \frac {\partial B_z}{\partial t} = -\frac {\partial E_y}{\partial x}, \end{equation}

(2.10) \begin{equation} -\frac {\partial B_z}{\partial x} = \mu _0\sum _s\frac {q_s\rho _s}{m_s}u_{s,y}+\mu _0\epsilon _0 \frac {\partial E_y}{\partial t}, \end{equation}

where $\mu _0$ is the vacuum permeability. Now, we may take perturbations around an initially uniform, stationary, field-less, but anisotropic initial condition:

(2.11) \begin{equation} \rho = \rho _0+\rho '\varpi, \quad u_i = u_i'\varpi, \end{equation}

(2.12) \begin{equation} p_{xx} = p_{xx,0}+p_{xx}'\varpi, \quad p_{yy} = p_{yy,0}+p_{yy}'\varpi, \quad p_{zz} = p_{zz,0}+p_{zz}'\varpi, \end{equation}

(2.13) \begin{equation} p_{xy} = p_{xy}'\varpi, \quad p_{xz} = p_{xz}'\varpi, \quad p_{yz} = p_{yz}'\varpi, \end{equation}

(2.14) \begin{equation} q_{xxy} = 0,\quad E_i = E_i'\varpi, \quad B_i = B_i'\varpi, \quad \end{equation}

where quantities subscripted with $0$ are invariant in space, primed quantities are assumed small and the linear perturbations are assumed as $\varpi = \exp \left [i\left (k_xx-\omega t\right )\right ]$ . The assumption that $q_{xxy} = 0$ is justified by comparing with the theory presented by Kuldinow et al. (Reference Kuldinow, Yamashita, Mansour and Hara2024b ), where $q_{xxy} \propto \partial _yT$ . Since we are assuming no variation in the $y$ direction, the estimate for the fluid heat flux vanishes. Note that this argument is independent of whether we are in the linear or nonlinear regime. For collisionless plasmas, spectral or non-local heat flux models, like those developed by Hammett & Perkins (Reference Hammett and Perkins1990) and Sharma et al. (Reference Sharma, Hammett, Quataert and Stone2006), may provide better agreement with kinetic results. The magnetic field suppresses cross-field transport (Mitchner & Kruger Jr. Reference Mitchner and Kruger1973), also limiting the importance of heat flux. However, studying the effect of heat flux closures on the behaviour of the fluid Weibel instability is reserved for future work.

Substituting these forms of the variables into (2.7)–(2.10) and keeping only terms to first order, we obtain

(2.15) \begin{equation} -i\omega u_y'+\frac {ik_xp_{xy}'}{\rho _0}=\frac {q}{m}E_y', \end{equation}

(2.16) \begin{equation} -i\omega p_{xy}'+ik_xp_{xx,0}u_y'=\frac {qB_z'}{m}\left (p_{yy,0}-p_{xx,0}\right ), \end{equation}

(2.17) \begin{equation} -i\omega B_z'=-ik_xE_y', \end{equation}

(2.18) \begin{equation} -ik_x B_z'=\mu _0\epsilon _0\sum _s\frac {q_s\rho _{s,0}}{m_s\epsilon _0}u_{s,y}'-i\mu _0\epsilon _0\omega E_y'. \end{equation}

Algebraic manipulations yield the following dispersion relation for the 10-moment Maxwell Weibel instability:

(2.19) \begin{equation} D_{\text {10M}}\left (k_x,\omega \right ) = 1-\frac {k_x^2c^2}{\omega ^2}-\sum _s\frac {w_{s,p}^2}{\omega ^2}\left (1+\frac {T_{s,\perp }}{2T_{s,\parallel }}\frac {1}{\xi _s^2-1/2}\right )=0. \end{equation}

Equation (2.19) is consistent with previous work (Sarrat et al. Reference Sarrat, Del Sarto and Ghizzo2016, Reference Sarrat, Del Sarto and Ghizzo2017) and care has been taken to express the dielectric function in a form reminiscent of the kinetic dispersion relation. By comparing with the kinetic dispersion relation, i.e. (2.2), we observe that for large $\xi$ , $D_{\text {K}}(k_x,\omega ) = D_{\text {10M}} (k_x,\omega )+O\left (\xi ^{-4}\right )$ . Recalling that $\xi$ is a normalised wavespeed, this indicates good agreement between the two models at large values of $\omega /k_x\propto \xi$ . It is discussed later that $\text {Re}(\omega )\approx 0\ll \text {Im}(\omega )$ for the Weibel instability, so large values of $\omega$ correspond to large Weibel growth rates.

2.3. Ten-moment fluid – Darwin approximation

In the Darwin approximation, the plasma is assumed to be non-relativistic, that is, the bulk velocity in (2.10) obeys $u\ll c$ , and also that the thermal speed $\sqrt {k_BT/m}\ll c$ . Thus, the displacement current is negligible in the frequencies of interest ( $\omega \ll ck_x$ ). The Darwin equations are written in terms of the electric scalar potential, $\phi$ , and the magnetic vector potential, $A_i$ , as

(2.20) \begin{equation} \frac {\partial ^2 \phi }{\partial x_k\partial x_k}=-\sum _s\frac {q_s\rho _s}{m_s\epsilon _0}, \end{equation}

(2.21) \begin{equation} \frac {\partial ^2 A_i}{\partial x_k\partial x_k}=\mu _0\epsilon _0\frac {\partial ^2\phi }{\partial t\partial x_i}-\sum _s\frac {\mu _0q_s\rho _su_{s,i}}{m_s}, \end{equation}

from which the electric and magnetic fields are found as

(2.22) \begin{equation} E_i = -\frac {\partial \phi }{\partial x_i}-\frac {\partial A_i}{\partial t},\hspace {30pt} B_i = \epsilon _{ijk}\frac {\partial A_k}{\partial x_j}. \end{equation}

In one spatial dimension, (2.20) and (2.21) reduce to

(2.23) \begin{equation} \frac {\partial ^2 \phi }{\partial x^2}=-\sum _s\frac {q_s\rho _s}{m_s\epsilon _0}, \end{equation}

(2.24) \begin{equation} \frac {\partial ^2 A_{\perp}}{\partial x^2}=-\sum _s\frac {\mu _0q_s\rho _su_{s,\perp}}{m_s}. \end{equation}

As the perturbations of $\rho$ are second-order effects, the perturbations of $\phi$ are also second order from (2.23). Similar to Maxwell’s equations, we only need the perturbations of $E_y$ and $B_z$ for the linear perturbation analysis. The linear perturbations of (2.22) and (2.24) yield

(2.25) \begin{equation} E_y' = i\omega A_y' = \frac {\omega }{k_x}\left (ik_xA_y'\right ) = \frac {\omega }{k_x}B_z' \end{equation}

and

(2.26) \begin{equation} B_z' = ik_xA_y' = \frac {i\mu _0\epsilon _0}{k_x}\sum _s\frac {q_s\rho _{s,0}}{m_s\epsilon _0}u_{s,y}', \end{equation}

respectively. By comparison to the Maxwell system in the previous section, (2.25) and (2.17) are identical, while (2.26) is equivalent to (2.18) but sans the displacement current term. Substituting the Darwin model ((2.25) and (2.26)) into the 10-moment fluid equations ((2.15) and (2.16)), we obtain the dispersion relation for the 10-moment Darwin Weibel instability:

(2.27) \begin{equation} D_{\text{10D}}\left (k_x,\omega \right ) = -\frac {k_x^2c^2}{\omega ^2}-\sum _s\frac {w_{s,p}^2}{\omega ^2}\left (1+\frac {T_{s,\perp }}{2T_{s,\parallel }}\frac {1}{\xi _s^2-1/2}\right )=0. \end{equation}

It is to be noted that the first term on the right-hand side of (2.19) is eliminated using the Darwin approximation in (2.27) under the assumption that $k_x^2c^2/\omega ^2\gg 1$ .

2.4. Magnetohydrodynamic approximation

Another common assumption in the plasmadynamics community is the magnetohydrodynamic (MHD) approximation. The plasma is considered as a single fluid, rather than a sum of individual fluids. The equations are then written in terms of the total mass density, $\bar{\rho}$ , mass velocity, $U_i$ , and net current density, $j_i$ :

(2.28) \begin{equation} \bar {\rho } = \sum _s\rho _s = \bar {n}\sum _sm_s, \end{equation}

(2.29) \begin{equation} U_i = \frac {1}{\bar {\rho }}\sum _s\rho _su_{s,i}, \end{equation}

(2.30) \begin{equation} j_i = \sum _s\frac {q_s\rho _su_{s,i}}{m_s}. \end{equation}

The MHD equations of motion are derived by summing or differencing the individual fluid species equations of motion to obtain equations for the bulk plasma quantities $\{\bar {\rho },U_i,j_i\}$ . The field dynamics are modified to use Ampère’s equation (neglecting the displacement current) and the generalised Ohm’s law, which is equivalent to the fluid momentum equations, assuming that $m_{\text{e}}/m_{\text{i}}\ll 1$ and $T_{\text{i}}\lesssim T_{\text{e}}$ , where $\text {e}$ and $\text {i}$ refer to electron and ion species, respectively. The evolution and field equations in MHD are written (Krall et al. Reference Krall, Trivelpiece and Gross1973)

(2.31) \begin{equation} \frac {\partial \bar {\rho }}{\partial t}+U_i\frac {\partial \bar {\rho }}{\partial x_i}+\bar {\rho }\frac {\partial U_i}{\partial x_i}=0, \end{equation}

(2.32) \begin{equation} \bar {\rho }\frac {\partial U_i}{\partial t}+U_j\frac {\partial U_i}{\partial x_j}=\epsilon _{ijk}j_jB_k-\frac {\partial p_{\text{e},ij}}{\partial x_j}, \end{equation}

(2.33) \begin{equation} E_i+\epsilon _{ijk}U_{j}B_k = -\frac {1}{e\bar {n}}\frac {\partial p_{\text{e},ij}}{\partial x_j}+\frac {m_{\text{e}}}{\bar {n}e^2}\frac {\text{d} j_{i}}{\text{d} t}, \end{equation}

(2.34) \begin{equation} \epsilon _{ijk}\frac {\partial B_k}{\partial x_j}=\mu _0j_i, \end{equation}

where the final term in (2.33), i.e. the unsteady component in Ohm’s law, provides feedback from the magnetic field to the currents and is important in recovering the correct Weibel instability dynamics. The electron inertia terms ( $\propto u\partial _xu$ ) are negligible for the linear perturbation theory assuming an initially stationary electron fluid, i.e. $\vec {u}_{\text{e},0} = 0$ . The generalised Ohm’s law in (2.33) can be inserted into Faraday’s equation to obtain

(2.35) \begin{equation} \frac {\partial B_i}{\partial t} = \epsilon _{ijk}\frac {\partial }{\partial x_j}\left (\epsilon _{klm}U_{l}B_m+\frac {1}{e\bar {n}}\frac {\partial p_{\text{e},kl}}{\partial x_l}-\frac {m_{\text{e}}}{\bar {n}e^2}\frac {\partial j_{k}}{\partial t}\right ). \end{equation}

This conclusion is similar in form to previous work discussing the role of the so-called ‘canonical battery effect’ in magnetic reconnection (Yoon & Bellan Reference Yoon and Bellan2019) and magnetogenesis (Laishram, Yun & Yoon Reference Laishram, Yun and Yoon2024).

Taking perturbations of (2.31)–(2.35) and coupling them to the pressure evolution equation, (2.16), yields the linear dispersion relation for the 10-moment MHD Weibel instability:

(2.36) \begin{equation} D_{10\text{MHD}}(k_x,\omega ) = -\frac {k_x^2c^2}{\omega ^2}-\frac {\omega _{\text{e},p}^2}{\omega ^2}\left ( 1 + \frac {T_{\text{e},\perp }}{2T_{\text{e},\parallel }} \frac {1}{\xi _{\text{e}}^2 - 1/2}\right )=0. \end{equation}

Here, it can be seen that (2.36) is consistent with (2.27) when only accounting for the electron susceptibility, $s=\text {e}$ . Without the unsteady component to Ohm’s law, the parenthetical term in (2.36) becomes $1-T_{\perp }/T_{\parallel }-2\xi ^2$ , which is inconsistent with the kinetic solution for all $\xi$ . Including the unsteady component, however, achieves an equivalent dispersion relation to the 10-moment model with Darwin equations. Thus, to accurately capture the Weibel instability using the MHD approximation, the unsteady component must be retained in the generalised Ohm’s law. Furthermore, our theory suggests that for extended MHD-type simulations, including a full pressure tensor but not retaining the unsteady component to Ohm’s law may lead to unphysical instability growth from the incorrect dispersion relation. A full derivation and discussion of (2.36) is presented in Appendix A.

3.1. Electron Weibel growth rate

Figure 1 shows the electron Weibel growth rate, i.e. $\tilde {\gamma } = \text {Im}\left (\tilde {\omega }\right )$ , as a function of $\tilde {k}_x$ and $T_{\text{e},\perp }/T_{\text{e},\parallel }$ for $\tilde {T}_{\text{e},\parallel }=6.25\times 10^{-4}$ or $\sqrt {\tilde {T}_{\text{e},\parallel }/\tilde {m}_{\text{e}}}=0.025$ for $\tilde {k}_x\in \left [0.01,10\right ]$ and $T_{\text{e},\perp }/T_{\text{e},\parallel }\in \left [1,100\right ]$ . Here, the ion susceptibility is neglected. There is negligible difference ( ${\lt}1\,\%$ for the entire plot) between the Maxwell and Darwin fluid growth rates for the parameters chosen, as shown in figures 1(b) and 1(c), respectively, since $\sqrt {\tilde {T}_{\text{e},\parallel }/\tilde {m}_{\text{e}}}\ll \tilde {c}=1$ . While good qualitative agreement can be observed between all the models, the 10-moment fluid dispersion relations (figures 1(b,c)) predict slightly higher growth rates than the kinetic dispersion relation (figure 1(a)). It is to be noted that the real frequency $\omega _r=\text {Re}\left (\tilde {\omega }\right )\simeq 0$ in each of these plots, because (2.2), (2.19) and (2.27) are functions of $\omega ^2$ in the high- $\xi$ limit and thus admit purely imaginary solutions.

Figure 1. Theoretical growth rates, $\tilde {\gamma }=\gamma /\omega _{\text{e},p}$ for the electron Weibel instability with $\tilde {T}_{\text{e},\parallel }=6.25\times 10^{-4}$ , showing the cutoff wavenumber (dashed line) and using the (a) kinetic dispersion relation, (2.2), (b) 10-moment Maxwell dispersion relation, (2.19), and (c) 10-moment Darwin dispersion relation, (2.27).

For every wavenumber $\tilde {k}_x$ , we note that there is a cutoff temperature ratio, $T_\perp /T_{\parallel }$ , under which the instability does not grow. Or conversely, for a given temperature ratio, there is a maximum wavenumber above which there is no instability. The cutoff condition occurs when $\xi =0$ solves the dispersion relation. Multiplying (2.2) by $\omega ^2$ and taking the limit of $\{\omega, \xi \}\rightarrow 0$ , noting that $Z'(0)=-2$ , since $Z'(\xi ) = -2\left [1+\xi Z(\xi )\right ]$ (see also (2.3)), we find that the cutoff condition can be written as

(3.2) \begin{equation} (k_x)_{\text{cutoff}} = \sqrt {\sum _s\frac {\omega _{s,p}^2}{c^2}\left (\frac {T_{s,\perp }}{T_{s,\parallel }} - 1\right )} = \frac {\omega _{\text{e},p}}{c}\sqrt {\sum _s\frac {m_{\text{e}}}{m_s}\left (\frac {T_{s,\perp }}{T_{s,\parallel }} - 1\right )}, \end{equation}

which describes the cutoff $k_x$ for a given $T_{\perp }/T_{\parallel }$ for each species. For the electron Weibel instability, we can also write the cutoff temperature ratio for a given $k_x$ as

(3.3) \begin{equation} \left (\frac {T_\perp }{T_{\parallel }}\right )_{\text{cutoff}}=1+\left (\frac {k_xc}{\omega _p}\right )^2. \end{equation}

Furthermore, by performing a similar treatment to the above with the 10-moment dispersion relations, (2.19) and (2.27), one obtains the exact same cutoff condition as the kinetic result. The cutoff condition is plotted as a dotted line in figure 1.

In figure 2, we present plots of the kinetic and fluid growth rates for various choices of the temperature ratio to provide a quantitative comparison. For a particular $k_x$ (e.g. $\tilde {k}_x\lt 0.3$ ), as the temperature ratio increases, the agreement between the kinetic and fluid predictions converge. This is due to the fact that the growth rate $\gamma$ (and $\xi \propto \gamma$ ) increases with $T_\perp /T_{\parallel }$ and thus the fluid dispersion relation more closely approximates the kinetic one. It can be seen from figure 2 that the cutoff condition occurs as predicted in (3.2).

Figure 2. Plot of $\tilde {\gamma }_{\text {K}}$ (solid lines; kinetic) and $\tilde {\gamma }_{\text {F}}$ (dashed lines; 10-moment Darwin, which is almost identical to 10-moment Maxwell) for the electron Weibel instability as a function of the wavenumber for various values of the temperature ratio $T_\perp /T_{\parallel }$ .

We also note that the fluid models provide good order-of-magnitude estimates for the maximum growth rate and the wavenumber at which the growth rate is largest. Using the 10-moment Darwin (10D) theory (2.27) and solving for $d\omega /dk=0$ , we can analytically determine the wavenumber of the maximum growth rate, i.e. $\left (k_x\right )_{\text{max}}$ , and the maximum growth rate, $\gamma _{\text{max}}$ :

(3.4) \begin{equation} \left (k_x\right )_{\text{max,10D}}=\frac {\omega _p}{c}\left (\sqrt {\frac {T_\perp }{T_{\parallel }}}-1\right )^{1/2}, \quad \gamma _{\text{max,10D}}=\frac {\omega _p}{c}\sqrt {\frac {T_\parallel }{m}}\left (\sqrt {\frac {T_\perp }{T_{\parallel }}}-1\right ). \end{equation}

To determine these values for the full kinetic dispersion relation, one would need to solve the plasma dispersion function, for which there is no closed-form solution.

To quantify the discrepancy between the kinetic and fluid dispersion relations, figure 3 presents the relative error of the linear growth rates between the fluid and kinetic theories. The relative error is defined as

\begin{equation*} \Delta = \left |\frac {\gamma _{\text {F}}-\gamma _{\text {K}}}{\gamma _{\text {K}}}\right |=\left |\frac {\gamma _{\text {F}}}{\gamma _{\text {K}}}-1\right |, \end{equation*}

where $\gamma _{\text {K}}$ is is the growth rate calculated from the kinetic theory, (2.2), and $\gamma _{\text {F}}$ is the growth rate calculated from the fluid theory, (2.19). The relative error is greatest near to the cutoff wavenumber, where both growth rates go to zero. However, the error decreases linearly with $T_\perp /T_{\parallel }$ at large temperature ratios (e.g. $T_\perp /T_\parallel \gt 3$ and $\tilde {k}_x \lt 0.3)$ . By combining (2.27) and (3.3), we can estimate for the fluid dispersion relations for $T_\perp /T_\parallel \gg 1$ :

(3.5) \begin{equation} \xi ^2 = \frac {1}{2}\left [1-\frac {T_\perp /T_{\parallel }}{\left (T_\perp /T_{\parallel }\right )_{\text{cutoff}}}\right ] \simeq -\frac {1}{2}\frac {T_\perp /T_{\parallel }}{\left (T_\perp /T_{\parallel }\right )_{\text{cutoff}}}, \end{equation}

where $\left (T_\perp /T_{\parallel }\right )_{\text{cutoff}}$ is given in (3.3). Furthermore, we can relate the relative error in the growth rate between the kinetic and fluid dispersion relations to the relative error in the estimate of the dispersion function:

(3.6) \begin{equation} \left |\frac {\delta \gamma }{\gamma }\right |=\left |\frac {\delta \xi }{\xi }\right |\simeq \left |\frac {\delta Z'}{2Z'}\right |, \end{equation}

where the final equality is obtained by estimating $Z'\simeq 1/\xi ^2$ . By comparing the first two terms of the kinetic $Z'$ (2.3) with the 10-moment Darwin dispersion relation (2.27), we find that

(3.7) \begin{equation} \left |\frac {\delta Z'}{Z'}\right | = \left |\frac {1}{Z'(\xi )} \left ( Z'(\xi )-\frac {1}{\xi ^2-1/2}\right )\right |\simeq \frac {1}{\xi ^2}. \end{equation}

Substituting (3.5) into (3.6) and (3.7), we find that the relative error in the growth rate between the kinetic and fluid dispersion relations can be written as

(3.8) \begin{equation} \Delta \approx \frac {\left (T_\perp /T_{\parallel }\right )_{\text{cutoff}}}{T_\perp /T_{\parallel }} \end{equation}

for large temperature ratios, $T_\perp /T_\parallel \gg 1$ . It is for this reason that in figure 3, we see, for instance, an error of approximately $10\,\%$ for $T_\perp /T_{\parallel }=10$ at small $k_x$ away from the cutoff wavenumber, e.g. $k_x \leqslant \left (k_x\right )_{\text{cutoff}}/10$ .

Figure 3. Relative discrepancy between kinetic and 10-moment Maxwell fluid growth rates $\Delta = |\gamma _F/\gamma _K-1| (\%).$ The dashed line denotes the cutoff wavenumber.

It is to be noted that, by combining (3.3), (3.4) and (3.8), the relative error in the growth rate estimate at the wavenumber of maximum growth can be obtained as

(3.9) \begin{equation} \Delta _{\text{max}} = \frac {1+\left (k_{\text{max}}c/\omega _p\right )^2}{T_\perp /T_\parallel } = \left (\frac {T_\perp }{T_\parallel }\right )^{-1/2}. \end{equation}

3.2. Ion Weibel growth rate

The previous analysis considered hot electrons with anisotropic temperature and cold (immobile) ions, that is, ions do not participate in the dynamics or the growth rate. The ion Weibel instability arises when there is a strong temperature anisotropy in the ions, rather than electrons. Due to their heavy mass, ions often require longer times to suffer sufficient collisions to come to equilibrium and thus can often retain more anisotropy than electrons (Ng et al. Reference Ng, Hakim, Wang and Bhattacharjee2020). Ion populations with different drift velocities can act as a large effective temperature anisotropy, leading to ion Weibel-like instabilities (Chang, Wong & Wu Reference Chang, Wong and Wu1990).

Figure 4 shows the ion Weibel growth rate, i.e. $\tilde{\gamma} = \text {Im}(\tilde{\omega} )$ , as a function of $\tilde {k}_x$ and $T_{{\text{i}},\perp }/T_{{\text{i}},\parallel }$ ; the cutoff condition as defined in (3.3) is plotted as a dotted line. We assume a hydrogen plasma (i.e. $m_{\text{i}}/m_{\text{e}} = 1836$ ) and set $T_{\text{e},\parallel } = T_{\text{e},\perp }$ , so that the electron Weibel instability does not occur. Here, $\sqrt {\tilde {T}_{\text{e},\parallel }/\tilde {m}_{\text{e}}}=\sqrt {\tilde {T}_{\text{e},\perp }/\tilde {m}_{\text{e}}}=0.025$ and $\sqrt {\tilde {T}_{{\text{i}},\parallel }/\tilde {m}_{\text{i}}}=0.001$ for $\tilde {k}_x\in \left [0.1,300\right ]$ and $T_{{\text{i}},\perp }/T_{{\text{i}},\parallel }\in \left [1,10^5\right ]$ . This study includes both warm electrons and ions where $T_{{\text{i}},\parallel }/T_{\text{e},\parallel } \approx 2.94$ . We note that the ion Weibel growth rate ( $\gamma _{\text {iW}}$ ) is much smaller than the electron Weibel growth rate ( $\gamma _{\text {eW}}$ ); namely, $\gamma _{\text {eW}}/\omega _{\text{e},p}=0.1$ is achieved in the electron Weibel instability for electron temperature ratios of approximately $T_{\text{e},\perp }/T_{\text{e},\parallel } = 25$ ; however, $\gamma _{\text {iW}}/\omega _{{\text{i}},p}=0.1$ is achieved for ion temperature ratios $T_{{\text{i}},\perp }/T_{{\text{i}},\parallel } \gtrsim 10^4$ . Similar to the electron Weibel instability, the discrepancy between the Maxwell and Darwin fluid growth rates is $\lt 1\,\%$ for the parameters chosen here, since the thermal speeds of each species are much less than the speed of light. The small value of $\gamma$ also ensures that $\left (k_xc/\omega \right )^2\gg 1$ and that the Darwin approximation is valid.

Figure 4. Theoretical growth rates, $\tilde {\gamma }=\gamma /\omega _{\text{e},p}$ , for the ion Weibel instability as a function of the ion temperature ratio and wavenumber. Here, $\tilde {T}_{\text{e},\parallel }=\tilde {T}_{\text{e},\perp }=6.25\times 10^{-4}$ and $\tilde {T}_{{\text{i}},\parallel }=1\times 10^{-6}$ . (a) Kinetic dispersion relation, (2.2). (b) Ten-moment Maxwell dispersion relation, (2.19). (c) Ten-moment Darwin dispersion relation, (2.27).

In figure 5, we present the kinetic and fluid growth rates for the various choices of the ion temperature ratio, $T_{{\text{i}},\perp }/T_{{\text{i}},\parallel }$ , to provide a quantitative comparison. Note that the behaviour here is different from the electron Weibel instability in figure 2. Most notably, the fluid growth rates (dashed lines) do not converge to the kinetic growth rate (solid lines) as $T_{{\text{i}},\perp }/T_{{\text{i}},\parallel }$ increases for all $k_x$ , while the cutoff condition obtained from the kinetic and 10-moment dispersion relations is identical, as explained in § 3.1. The discrepancy even at small $k_x$ can be explained by noting that the kinetic dispersion relation reduces to the fluid theory only when both $\xi _{\text{e}}\gg 1$ and $\xi _{\text{i}} \gg 1$ . For the electron Weibel instability, indeed $\xi _{\text{e}} \gg 1$ , resulting in the convergence of kinetic and fluid theories at small $\tilde {k}_x$ , as shown in figure 2. However, for the ion Weibel instability, the values of $\xi$ are reduced compared with the electron Weibel instability because (1) $|\tilde {\gamma }|=|\tilde {\omega }|\lt 10^{-2}$ is small (note: $\omega _r=0$ ) and (2) the growth rates are positive for higher $\tilde {k}_x$ . These factors lead to $\xi _{{e}}\ll 1$ , contributing to the discrepancies between the fluid and kinetic theories. Nonetheless, the fluid theories predict the correct cutoff condition, and we note qualitatively good agreement in the magnitude of the ion Weibel instability growth rate obtained from the 10-moment dispersion relation with the kinetic theory due to $\xi _{\text{i}}\gtrsim 1$ .

Figure 5. Plot of $\tilde {\gamma }_{\text {K}}$ (solid lines; kinetic) and $\tilde {\gamma }_{\text {F}}$ (dashed lines; 10-moment Darwin, which is almost identical to 10-moment Maxwell) for the ion Weibel instability as a function of the wavenumber for various values of the temperature ratio $T_{{\text{i}},\perp }/T_{{\text{i}}, \parallel }$ .

5.1. Kinetic model

We use an implicit charge-, canonical momentum- and energy-conserving Darwin electromagnetic particle-in-cell (EM-PIC) model, described by Chen & Chacón (Reference Chen and Chacón2014). The model has been tested against Weibel instability theory and other standard electromagnetic test cases.

The EM-PIC model employs a Crank–Nicolson discretisation for the particle update:

(5.1) \begin{equation} \frac {x^{\nu +1}_p-x_p^\nu }{\Delta \tau ^\nu } = v_{p,x}^{\nu +\frac {1}{2}}, \end{equation}

(5.2) \begin{equation} \frac {v^{\nu +1}_{p,i}-v_{p,i}^\nu }{\Delta \tau ^\nu } = \frac {q_p}{m_p}\left (E_{p,i}^{\nu +\frac {1}{2}}+\epsilon _{ijk}v_{p,j}^{\nu +\frac {1}{2}}B_{p,k}^{\nu +\frac {1}{2}}\right ), \end{equation}

where $\{x_p,v_{p,i}\}$ is the particle position and velocity, $\Delta \tau ^{\nu }$ is the $\nu {\text {th}}$ substep within $\Delta t$ , $E_{i,p}^{\nu +\frac {1}{2}} = E_{i}\left (x_{p}^{\nu +\frac {1}{2}},t^{\nu +\frac {1}{2}}\right )$ and $B_{i,p}^{\nu+\frac{1}{2}} = B_i\left(x_p^{\nu +\frac{1}{2}},t^{\nu +\frac{1}{2}}\right)$ are the electric and magnetic fields interpolated to the particle position at the half-substep time, and following a Crank–Nicolson-type discretisation, $X^{\nu +\frac {1}{2}}=\left (X^{\nu +1}+X^{\nu }\right )/2$ for all quantities $X$ . The particle velocity is updated using the Boris method. Note that both the particle positions and velocities are defined on the full-substep times, and the velocity update, (5.2), requires the fields at $x^{\nu +\frac {1}{2}}$ , which itself is dependent on $v^{\nu +\frac {1}{2}}$ , and thus $v^{\nu +1}$ . Thus, the particle evolution equations, (5.1) and (5.2), are solved implicitly using Picard iteration. The particle updates are coupled with the electromagnetic field updates using the Darwin model:

(5.3) \begin{equation} \epsilon _0\frac {E_{x,\ell +\frac {1}{2}}^{n+1}-E_{x,\ell +\frac {1}{2}}^{n}}{\Delta t}+\sum _s\bar {j}_{s,x,\ell +\frac {1}{2}}^{n+\frac {1}{2}}-\left \langle \bar {j}^{n+\frac {1}{2}}_x\right \rangle =0, \end{equation}

(5.4) \begin{equation} \frac {1}{\mu _0}\frac {A_{\perp, \ell -1}^{n+\frac {1}{2}}-2A_{\perp, \ell }^{n+\frac {1}{2}}+A_{\perp, \ell +1}^{n+\frac {1}{2}}}{\Delta x^2}+\sum _s\bar {j}_{s,\perp, \ell }^{n+\frac {1}{2}}-\left \langle \bar {j}_\perp ^{n+\frac {1}{2}}\right \rangle =0, \end{equation}

where the current components $\bar {j}_{s,i,l}$ are orbit averaged over the substeps for each species, $s$ . Here,

(5.5) \begin{equation} \bar {j}^{n+\frac {1}{2}}_{x,\ell +\frac {1}{2}} = \frac {1}{\Delta t\Delta x}\sum _p\sum _\nu q_pv^{\nu +\frac {1}{2}}_{p,x}S_m\left (x_p^{\nu +\frac {1}{2}}-x_{\ell +\frac {1}{2}}\right )\Delta \tau ^{\nu }, \end{equation}

(5.6) \begin{equation} \bar {j}^{n+\frac {1}{2}}_{\perp, \ell } = \frac {1}{\Delta t\Delta x}\sum _p\sum _\nu q_pv^{\nu +\frac {1}{2}}_{p,\perp }S_l\left (x_p^{\nu +\frac {1}{2}}-x_{\ell }\right )\Delta \tau ^{\nu }, \end{equation}

where $S_m$ and $S_l$ are B-spline shape functions used for consistency with the field-to-particle interpolation to achieve energy conservation, $\langle \bar {j}_i\rangle =\sum _{s,\ell } \bar {j}_{s,i,\ell }/N_x$ is the spatial average of the net current density (i.e., adding both ion and electron current densities) and is necessary in periodic geometries to enforce periodicity of the fields (Chen, Chacón & Barnes Reference Chen, Chacón and Barnes2011; Chen & Chacón, Reference Chen and Chacón2014), and $N_x$ is the number of cells in the $x$ direction.

5.2. Fluid model

The 10-moment plasma model solves the coupled set of fluid equations, which are written in one dimension using conservative variables:

(5.7) \begin{equation} \frac {\partial \textbf {U}}{\partial {t}} + \frac {\partial {\textbf {F}}}{\partial {x}} = \textbf {S}, \end{equation}

where the vector of conservative variables $\textbf {U} = \{\rho, \rho u_i, \rho u_iu_j+p_{ij}\}=\{\rho, \Gamma _i,S_{ij}\}$ , the vector of conservative fluxes $\textbf {F} = \{\rho u_x, \rho u_iu_x+p_{ix}, \rho u_iu_ju_x+u_ip_{jx}+u_jp_{xi}+u_xp_{ij}\}$ and the vector of conservative sources $\textbf {S} = \{0,\rho E_i+\epsilon _{ijk}\rho u_jB_k,\rho u_i E_j+\rho u_j E_i +\epsilon _{ikl}\left (p_{jk}+\rho u_ju_k\right )B_l +\epsilon _{jkl}\left (p_{ik}+\rho u_i u_k\right )B_l\}$ . Note that this conservative form of the equations is consistent with the primitive form which was used in the derivation of the dispersion relations in § 2. The fluid equation is coupled with the Darwin equations, given in (2.23) and (2.24).

5.2.1. Numerical scheme

For the numerical solver, we employ a finite volume approach, in which the state variables are defined at the cell centres (Kuldinow et al. Reference Kuldinow, Yamashita and Hara2024a ,Reference Kuldinow, Yamashita, Mansour and Hara b ; Sahu, Mansour & Hara Reference Sahu, Mansour and Hara2020). The fluid equations in one dimension are discretised using a Crank–Nicolson scheme:

(5.8) \begin{equation} \frac {\textbf {U}_\ell ^{n+1}-\textbf {U}_\ell ^n}{\Delta t} + \frac {\textbf {F}_{\ell +\frac {1}{2}}^{n+\frac {1}{2}}-\textbf {F}_{\ell -\frac {1}{2}}^{n+\frac {1}{2}}}{\Delta x} = \textbf {S}_\ell ^{n+\frac {1}{2}}, \end{equation}

while the Darwin equations are similarly discretised:

(5.9) \begin{equation} \epsilon _0\frac {E_{x,\ell }^{n+1}-E_{x,\ell }^{n}}{\Delta t}+\sum _s j_{s,x,\ell }^{n+\frac {1}{2}}-\left \langle j_x^{n+\frac {1}{2}}\right \rangle =0, \end{equation}

(5.10) \begin{equation} \frac {1}{\mu _0}\frac {A_{\perp, \ell -1}^{n+\frac {1}{2}}-2A_{\perp, \ell }^{n+\frac {1}{2}}+A_{\perp, \ell +1}^{n+\frac {1}{2}}}{\Delta x^2}+\sum _s j_{s,\perp, \ell }^{n+\frac {1}{2}}-\left \langle j_\perp ^{n+\frac {1}{2}}\right \rangle =0, \end{equation}

where $\ell$ is the cell index and $\ell +\frac {1}{2}$ are cell interface positions, the superscript $n$ denotes the time step, $j_{s,i} = q_s\Gamma _{s,i}/m_s$ and $\langle j_i\rangle = \sum _{s,\ell } j_{s,i,\ell }/N_x$ . All quantities are defined at cell centres; interfacial values obtained using the monotonic upwinding scheme for conservation laws (MUSCL) with the van Leer harmonic limiter (van Leer Reference van Leer1979) are used to calculate the flux quantities.

Equations (5.8)–(5.10) can be written in the form

(5.11) \begin{equation} \mathcal {F}\left (\mathcal {M}^{n+1},E_x^{n+1},A_\perp ^{n+1}\right )=\mathcal {F}\left (\textbf {V}^{n+1}\right )=\textbf {R}\Rightarrow 0, \end{equation}

where $\mathcal {M}=\{\rho, \Gamma _i,S_{ij}\}$ is the vector of fluid moments at all grid points and all quantities at time $t=t^n$ are assumed to be known for the nonlinear function $\mathcal {F}$ , and $\textbf {R}$ is the residual which we want to vanish. The root of (5.11) is solved by a preconditioned Jacobian-free Newton Krylov iteration (Kelley Reference Kelley1995; Knoll & Keyes Reference Knoll and Keyes2004). Given an initial guess for the solution at iteration $\kappa$ , $\textbf {V}^{n+1,(\kappa )}$ , we can calculate the residual and update the solution vector as

(5.12) \begin{equation} \textbf {V}^{n+1,(\kappa +1)}=\textbf {V}^{n+1,(\kappa )}+\delta \textbf {V}^{(\kappa )}, \end{equation}

where the update to the solution vector satisfies

(5.13) \begin{equation} \left [\frac {\partial \mathcal {F}}{\partial \textbf {V}}\right ]^{(\kappa )}\delta \textbf {V}^{(\kappa )}=-\textbf {R}^{(\kappa )}. \end{equation}

The inversion of the Jacobian matrix is handled through the preconditioned generalised minimum residual method (GMRES), and the Jacobian-vector product is estimated by a finite difference:

(5.14) \begin{equation} \left [\frac {\partial \mathcal {F}}{\partial \textbf {V}}\right ]^{(\kappa )}\delta \textbf {V}^{(\kappa )} \simeq \frac {1}{\epsilon }\left [\mathcal {F}\left (\textbf {V}^{(\kappa )}+\epsilon \delta \textbf {V}^{(\kappa )}\right )-\mathcal {F}\left (\textbf {V}^{(\kappa )}\right )\right ] \end{equation}

for a small value of $\epsilon \simeq 10^{-8}$ . The iteration ( $\kappa$ ) is performed until convergence, determined by the magnitude of the residual $|\textbf {R}|\lt \text {tol}=10^{-8}$ .

5.2.2. GMRES preconditioner

The GMRES preconditioner modifies (5.13) as

(5.15) \begin{equation} M^{-1}\left (\left [\frac {\partial \mathcal {F}}{\partial \textbf {V}}\right ]^{(\kappa )}\delta \textbf {V}^{(\kappa )}\right )=M^{-1}\left (-\textbf {R}^{(\kappa )}\right ), \end{equation}

where $M^{-1}$ is a function that is an approximate inverse of $\frac {\partial \mathcal {F}}{d\textbf {V}}$ , i.e. $M^{-1}(\textbf {X})\simeq \left (\frac {\partial \mathcal {F}}{d\textbf {V}}\right )^{-1}\textbf {X}$ , allowing for faster convergence of the GMRES method. The function $M^{-1}$ is determined by solving a simplified form of (5.13) and (5.14). To simplify and make the equations analytically invertible, we neglect the fluxes, $\textbf {F}$ , in (5.8) and neglect the effect of magnetic field coupling on the transverse momenta. Specifically, writing out the components of the simplified Newton step, we find

(5.16) \begin{equation} -R_{\rho, \ell } = \frac {\delta \rho _\ell }{\Delta t}, \end{equation}

\begin{align*} \hspace{-45pt}-R_{\Gamma _x,\ell } = \frac {\delta \Gamma _{x,\ell }}{\Delta t} -\frac {q}{2m}\Bigg ( \delta \rho _\ell E_{x,\ell }^{n+\frac {1}{2}}+\rho _\ell ^{n+\frac {1}{2}}\delta E_{x,\ell }\hspace{45pt} \end{align*}

(5.17) \begin{align} +{\kern2pt}\delta \Gamma _{y,\ell }B_{z,\ell }^{n+\frac {1}{2}}-\delta \Gamma _{z,\ell }B_{y,\ell }^{n+\frac {1}{2}} +\Gamma _{y,\ell }^{n+\frac {1}{2}}\delta B_{z,\ell }-\Gamma _{z,\ell }^{n+\frac {1}{2}}\delta B_{y,\ell } \Bigg ), \end{align}

(5.18) \begin{equation} -R_{\Gamma _\perp, \ell } = \frac {\delta \Gamma _{\perp, \ell }}{\Delta t}-\frac {q}{2m}\left (\delta \rho _\ell E_{\perp, \ell }^{n+\frac {1}{2}}+\rho _\ell ^{n+\frac {1}{2}}\delta E_{\perp, \ell }\right ), \end{equation}

\begin{align*} \hspace{-20pt}-R_{S_{ij},\ell } = \frac {\delta S_{ij,\ell }}{\Delta t}-\frac {q}{2m}\Bigg [ \delta \Gamma _{i,\ell } E_{j,\ell }^{n+\frac {1}{2}}+\Gamma _{i,\ell }^{n+\frac {1}{2}}\delta E_{j,\ell } +\delta \Gamma _{j,\ell } E_{i,\ell }^{n+\frac {1}{2}}+\Gamma _{j,\ell }^{n+\frac {1}{2}}\delta E_{i,\ell } \hspace{20pt}\end{align*}

(5.19) \begin{align} \hspace{5pt}+{\kern2pt}\epsilon _{ikl}\left (\delta S_{jk,\ell }B_{l,\ell }^{n+\frac {1}{2}}-\delta S_{jk,\ell }B_{l,\ell }^{n+\frac {1}{2}}\right ) +\epsilon _{jkl}\left (S_{ik,\ell }^{n+\frac {1}{2}}\delta B_{l,\ell }-S_{jk,\ell }^{n+\frac {1}{2}}\delta B_{l,\ell }\right ) \Bigg ],\hspace{-5pt} \end{align}

(5.20) \begin{equation} -R_{E_x,\ell } = \epsilon _0\frac {\delta E_{x,\ell }}{\Delta t}+\sum _s\frac {q_s\Gamma _{x,\ell }}{m_s}-\langle j_x\rangle, \end{equation}

(5.21) \begin{equation} -R_{A_\perp, \ell } = \frac {1}{\mu _0}\frac {\delta A_{\perp, \ell -1}-2\delta A_{\perp, \ell }+\delta A_{\perp, \ell +1}}{\Delta x^2}+\sum _s\frac {q_s\Gamma _{\perp, \ell }}{m_s}-\langle j_\perp \rangle, \end{equation}

where $R_{Q,\ell }$ is the residual of the equation for quantity $Q$ , and $\delta Q$ are the components of $\delta \textbf {V}$ . In addition, the field quantities are updated as

(5.22) \begin{equation} \delta E_{\perp, \ell } = -\frac {\delta A_{\perp, \ell }}{\Delta t}, \end{equation}

(5.23) \begin{equation} \delta B_{i,\ell } = \epsilon _{ixk}\frac {\delta A_{k,\ell +1}-\delta A_{k,\ell -1}}{2\Delta x}. \end{equation}

Equations (5.18), (5.21) and (5.22) form a closed set, equivalent to a tridiagonal matrix-vector product, through which the preconditioner updates $\delta A_\perp$ can be obtained and thus, $\delta E_\perp$ , $\delta \Gamma _\perp$ and $\delta B_i$ . Once those updates are known, (5.16), (5.17) and (5.20) are closed, and can be solved to obtain $\delta \rho$ , $\delta \Gamma _x$ and $\delta E_x$ . Finally, the previously determined quantities can be be substituted into (5.19) to solve for $\delta S_{ij}$ . This approximate inversion provides an approximate solution to (5.13), which is used as a preconditioner to the GMRES solution of the nonlinear Newton step, (5.15).

6.1. Initialisation

The set-up is an electron Weibel instability; all presented quantities are normalised as described in (3.1); we drop the tilde on all variables for the remainder of this section. We assume a periodic domain of length $L_x = 32$ (i.e. the domain wavenumber is $k_{L} = 2\pi / L_x \approx 0.196$ ) with $N_x=128$ uniformly spaced cells, using the normalised units defined above. The kinetic simulations employ 2000 particles per cell and both models employ a time step of $\Delta t = 1$ (recall that time is normalised with $\omega _{\text{e}, p}^{-1}$ ). The initial electron distribution function is assumed to be a shifted bi-Maxwellian distribution:

(6.1) \begin{equation} f_{\text{e}}(x,v_i;\,t=0) = \frac {1}{\sqrt {2\pi } v_{\text{e,th},\parallel }v_{\text{e,th},\perp }^2}\exp \left [\frac {v_x^2}{2v_{\text{e},\text{th},\parallel }^2}+\frac {\left [v_y-u_{\text{e},y}(x)\right ]^2+v_z^2}{2v_{\text{e},\text{th},\perp }^2}\right ], \end{equation}

where $v_{\text{e,th},i}=\sqrt {T_{\text{e},i}}$ is the thermal velocity in the $i$ direction and $u_{\text{e},y}(x)=\delta \cos (k_L x)$ is a small initial transverse velocity perturbation, where $\delta = 1\times 10^{-5}$ and $k_{L}=2\pi /L_x$ . We consider two test cases where $T_\perp /T_\parallel = 2.56$ and $T_\perp /T_\parallel = 25.6$ . For both test cases, $T_\parallel = T_x = 6.25\times 10^{-4}$ . As shown in figure 1, for $T_\perp /T_\parallel = 2.56$ , $\gamma _{\text{K}} = 4\times 10^{-3}$ and $\gamma _{\text{10M}} = 6\times 10^{-3}$ ( $|\xi |\sim 0.6$ ), approximately a $50\,\%$ discrepancy between the kinetic and 10 M dispersion relations. However, for $T_\perp /T_\parallel = 25.6$ , $\gamma _{\text{K}} = 2.30\times 10^{-2}$ and $\gamma _{\text{10M}} = 2.38\times 10^{-2}$ ( $|\xi |\sim 3.3$ ), which results in only an approximately $3\,\%$ discrepancy, despite the moderate size of $\xi$ .

In both cases, the initial conditions seed the fundamental wavenumber ( $m=1$ ), determined by the domain size, which is not the wavenumber of maximum growth, $k_{\text{max}}$ , predicted by (3.4). For $T_\perp /T_\parallel = 2.56$ , the growth rate is maximum at $k_{\text{max}} = 0.775$ , which is 3.9 times larger than $k_{L}$ and the growth rate evaluated for $k_{L}$ is 3.7 times smaller than the maximum growth rate. Additionally, for $T_\perp /T_\parallel = 25.6$ , the growth rate is maximum at $k_{\text{max}} = 2.01$ , which is $10.3$ times larger than $k_{L}$ , resulting in a growth rate 4.3 times smaller than $\gamma _{\text{max}}$ . The initial conditions are kept ‘quiet’, i.e. noise-less, so that the fundamental $k_L$ remains dominant through saturation. It is to be noted that in simulations with multimode or noisy initial conditions, the wavenumber of maximum growth, as predicted by the linear theory, is dominant.

The kinetic simulations each took approximately 4.5 hours, while the fluid simulations each took approximately 0.5 hours, using a single processor. It is however noted that the comparison of computational cost is highly dependent on the number of particles per cell in the kinetic simulation, and likewise, how important finite particle noise is to the system dynamics.

6.2. Comparison of results to linear and saturation theory

Figure 6 shows the total magnetic field energy, $\int _{\Omega _x}\frac {1}{2\mu _0}|B|^2 \text{d}x=\frac {\Delta x}{2\mu _0}\sum _\ell |B_\ell |^2$ , as a function of time for the two $T_\perp /T_\parallel$ cases.

Figure 6. Temporal variation of magnetic field energy obtained from kinetic and fluid models. Comparisons are made to linear growth rate predictions, (2.2) and (2.27), and estimates of saturation field strength, (4.2), for (a) $T_\perp / T_\parallel = 2.56$ and (b) $T_\perp / T_\parallel = 25.6$ .

In figure 6(a), we find good agreement between the simulation results and the linear growth rates obtained from both kinetic and fluid theories for $T_\perp /T_\parallel = 2.56$ . It is to be noted that both computational models use the Darwin equations, thus neglecting the displacement current. For the non-relativistic set-ups ( $c\gg \omega /k$ ) presented in this work, the discrepancy between the theoretical Maxwell and Darwin growth rates is less than $1\,\%$ , and does not drastically affect the saturation properties. We also find that, despite the different growth rates, the results obtained from the kinetic (PIC) and 10-moment fluid models saturate at approximately the same magnetic field energy, consistent with the field strength predicted by the saturation theory, as shown in (4.2).

Figure 6(b) shows good agreement between the simulation results and theory for the higher anisotropy case, i.e. $T_\perp /T_\parallel = 25.6$ . The difference between the two predicted growth rates obtained from PIC and 10-moment dispersion relations agree to within $2.5\,\%$ , and again the results from both models saturate at nearly the same field strength, similar to that predicted by theory.

6.3. Comparison of magnetic field profiles

To provide a qualitative comparison of the fluid and kinetic results outside of the global metric of total magnetic field energy, we present the contour plots of the magnetic field in space and time. However, the magnetic field spans many scales due to the linear exponential growth phase and admits positive and negative values. Thus, for visualisation purposes, we will plot the symmetric logarithm of the magnetic field:

(6.2) \begin{equation} \text {symlog}\left [B(x),C\right ] = \text {sgn}\left [B(x)\right ]\log _{10}\left [1+C\left |B(x)\right |\right ]. \end{equation}

Unlike the common logarithm, the symmetric logarithm is continuous at zero but has the same asymptotic behaviour. For the remainder of this study, we will employ $C=10^{3}$ . This comes at the cost of losing information at scales smaller than $10^{-3}$ , for this choice of the constant. However, since we will be directing our attention towards the saturation properties of the magnetic field, this will not affect our discussion. Note that for $|B|\gg 10^{-3}$ , $\text {symlog}\left [B(x),10^3\right ]\simeq \text {sgn}\left [B(x)\right ]\left [3+\log _{10}\left |B(x)\right |\right ]$ .

6.3.1. Low temperature ratio: $T_\perp /T_{\parallel }=2.56$

Figure 7 shows $\text {symlog}\left [B(x),10^3\right ]$ for the initial temperature ratio $T_\perp /T_{\parallel }=2.56$ . Recall that for this lower temperature ratio, the growth rates of the kinetic and fluid descriptions of the Weibel instability differ by approximately 50 %, as shown in figure 6(a). The kinetic result, shown in figure 7(a), does not reach the peak magnetic field until $t\simeq 1250$ , while the 10-moment result in figure 7(b) saturates early at approximately $t\simeq 800$ . This is because the growth rates are estimated as ${\gamma }_{\text{K}} = 4\times 10^{-3}$ and $\gamma _{\text{10M}} = 6\times 10^{-3}$ for kinetic and 10-moment dispersion relations, respectively. We note that the kinetic simulation, due to wave mixing and non-thermal particles, generates high- $k_x$ waves that look like ripples in the contour plots. These high- $k_x$ structures born from individual particle dynamics eventually come to dominate the dynamics at saturation. In contrast, the 10-moment simulations capture neither the non-thermal particles nor significant deviations from Maxwellian; near saturation, the fluid model harmonics of the initial $m=1$ mode also saturate and we observe less structured behaviour. Despite such a difference, the linear growth and nonlinear saturation is captured using the 10-moment fluid model.

Figure 7. Symmetric logarithm of the magnetic field profile (6.2) for $T_\perp /T_{\parallel }=2.56$ : (a) kinetic simulation and (b) 10-moment fluid simulation.

6.3.2. High temperature ratio: $T_{\perp }/T_{\parallel }=25.6$

Figure 8 shows $\text {symlog}\left [B(x),10^3\right ]$ for the initial temperature ratio $T_\perp /T_{\parallel }=25.6$ . Recall that for this higher temperature ratio, the growth rates of the kinetic and fluid descriptions of the Weibel instability agree well, as shown in figures 1 and 3. Furthermore, we can see that the profiles of the kinetic and fluid models at saturation $\left (t \simeq 275\right )$ in figure 8 show approximately $m=3$ modes ( $k_x=3k_{x,0}$ , where $k_{x,0}=2\pi /L_x$ ) before both models predict that higher $k_x$ modes also saturate, leading to fine structure in the magnetic field profiles.

Figure 8. Symmetric logarithm of the magnetic field profile (6.2) for $T_\perp /T_{\parallel }=25.6$ : (a) kinetic simulation and (b) 10-moment simulation.

As with the lower temperature case, in figure 7, we note many small-amplitude, high- $k_x$ ripples in the kinetic simulation, which grow to share nearly the same magnetic field energy as the fundamental mode at saturation. However, the fluid simulation is more ordered in the sense that after saturation, we observe a large-amplitude $m=3$ mode near $t = 310$ , then an $m=5$ mode near $t=350$ , rather than the spectrum visible in the kinetic simulation. Again, this is most likely due to there being no damping mechanism in the fluid simulation (e.g. phase mixing), while the kinetic simulation necessarily exhibits particle–wave resonance and non-thermal diffusion. However, despite the fact that we are not capturing these inherently kinetic effects, it is remarkable that the 10-moment fluid model can accurately capture the linear growth and saturation level.