2.1. The 10-moment collisionless fluid system

Our simulations use the 10-moment collisionless fluid model of the Gkeyll continuum plasma physics solver framework. The 10-moment model employed here uses the following set of non-relativistic moment equations:

(2.1) \begin{align} \frac {\partial n}{\partial t} + \frac {\partial }{\partial x_j} (n u_j) &= 0,\\[-10pt]\nonumber \end{align}

(2.2) \begin{align} m n \frac {\partial u_j}{\partial t} + m n u_k \frac {\partial u_j}{\partial x_k} + \frac {\partial p_{jk}}{\partial x_k} &= n q [ E_j + \epsilon _{jkl} u_k B_l ],\\[-10pt]\nonumber \end{align}

(2.3) \begin{align} \frac {\partial p_{jk}}{\partial t} + \frac {\partial }{\partial x_l} (p_{jk} u_l) + p_{jl} \frac {\partial u_k}{\partial x_l} + p_{kl} \frac {\partial u_j}{\partial x_l} + \frac {\partial q_{jkl}}{\partial x_l} &= \frac {q}{m} [ \epsilon _{jlm} p_{lk} B_m + \epsilon _{klm} p_{lj} B_m ], \end{align}

where we have suppressed the species index $\alpha$ , corresponding to mass and charge $m_\alpha$ and $q_\alpha$ , but otherwise $n_\alpha =\int f_\alpha\, \text{d}^3v$ is the density, where $f_\alpha (\boldsymbol{x},\boldsymbol{v},t)$ is the particle distribution function, and the fluid flow velocity vector, pressure tensor and heat flux tensor components are defined by $u_{j,\alpha }=n_\alpha ^{-1} \int v_j f_\alpha\, \text{d}^3v$ , $p_{jk,\alpha }=m_\alpha \int (v_j-u_j)(v_k-u_k) f_\alpha\,\text{d}^3v$ and $q_{jkl,\alpha }=m_\alpha \int (v_j-u_j)(v_k-u_k)(v_l-u_l) f_\alpha\, \text{d}^3v$ , respectively. Furthermore, $E_j$ and $B_j$ denote the electric and magnetic field components, which are evolved according to Maxwell’s equations, $\partial _t \boldsymbol{B}=-\boldsymbol{\nabla }\times \boldsymbol{E}$ and $\partial _t \boldsymbol{E}=c^2\boldsymbol{\nabla }\times \boldsymbol{B}-\epsilon _0^{-1}\sum _\alpha q_\alpha n_\alpha \boldsymbol{u}_\alpha$ . The solver uses a locally implicit scheme for the source terms, while solving the homogeneous portion of the equations by the wave propagation method (LeVeque Reference LeVeque1997). This approach allows certain unimportant scales to be under-resolved without leading to numerical difficulties. In particular, our simulations strongly under-resolve the Debye scale, which would not be possible in most particle-in-cell kinetic codes. This provides us greater flexibility in the plasma parameter regimes considered.

Importantly, the full pressure tensor is evolved for each species, which allows some pressure anisotropy driven instabilities (Gary Reference Gary1993) to be retained, such as the firehose instability (Walters et al. Reference Walters, Klein, Lichko, Juno and TenBarge2024) and, with particular relevance to the current work, the Weibel instability. Previous work by Kuldinow & Hara (Reference Kuldinow and Hara2025) employing the 10-moment system with a zero heat flux closure ( $q_{jkl}=0$ ) demonstrated that the model can predict electron Weibel instability growth rates in the long wavelength and the strongly driven limits, and accurately reproduces the high wavenumber cutoff, as well as the nonlinear saturation level, of the instability.

Here, we use a closure that drives the system towards an isotropic pressure (Wang et al. Reference Wang, Hakim, Bhattacharjee and Germaschewski2015)

(2.4) \begin{equation} \frac {\partial q_{jkl,\alpha }}{\partial x_l} = k_{0,\alpha } \frac {v_{\text{th},\alpha }}{\sqrt {2}} (p_{jk,\alpha } - p_\alpha \delta _{jk}), \end{equation}

where $p_\alpha =\delta _{lm}p_{lm,\alpha }/3$ is the scalar pressure and $v_{\text{th},\alpha }=\sqrt {2 p_\alpha /(n_\alpha m_\alpha )}$ is the thermal speed. The parameters $k_{0,\alpha }$ adjust the strength of the isotropisation for each species. This isotropisation closure is not derived to capture a specific kinetic process (in contrast to, e.g. the Hammett–Perkins closure (Hammett & Perkins Reference Hammett and Perkins1990) for Landau damping) and the parameters $k_{0,\alpha }$ should therefore be regarded as free. However, by studying the parameters’ effect in a controlled manner, they can be used to approximate selected aspects of kinetic physics. In particular, as discussed in § 2.2, they regulate the damping of flows and magnetic fields in the unmagnetised limit, allowing one to emulate effective kinetic dissipation. In addition, as we show in Appendix A, the closure parameter can also be tuned to reproduce key characteristics of the electron Weibel instability. While other closures for the heat flux are available in Gkeyll, such as non-local (Ng et al. Reference Ng, Hakim, Bhattacharjee, Stanier and Daughton2017) or gradient-based (Ng et al. Reference Ng, Hakim, Wang and Bhattacharjee2020), we choose the local relaxation closure for its explicit connection to dissipation and our ability to determine effective Reynolds numbers from the choice of the closure parameter.

2.2. Effect of closure parameters and effective Reynolds numbers

Unlike in resistive-MHD, where the concepts of dynamo theory have mostly been developed, in the 10-moment system, or a collisionless kinetic system for that matter, there is no explicit kinematic viscosity $\nu$ and magnetic diffusivity $\eta$ . In addition, the corresponding damping rates of flow and magnetic perturbations can, in general, have different wavenumber ( $k$ ) dependence to the $\propto k^2$ diffusive behaviour. Thus, the usual definitions of the fluid and magnetic Reynolds numbers, $\textrm{Re}=u_0 L/\nu$ and $\textrm{Rm}=u_0 L/\eta$ , respectively, are less trivial to evaluate,Footnote ⁴ and they have different implications concerning the scale separation between the outer scale and the fluctuation cutoff scales.

To put our results in context, it is nevertheless useful to quantify the magnetic and fluid damping in terms of analogous dimensionless quantities, comparing inverse time scales of processes creating and dissipating magnetic and fluid kinetic energy on the outer scale $L$ . The fluid Reynolds number compares inertial processes ${\sim} u_0/L$ with (viscous) flow damping ${\sim} \nu /L^2$ , while the magnetic Reynolds number compares field line stretching by the flow ${\sim} u_0/L$ with resistive dissipation of magnetic fields ${\sim} \eta /L^2$ . Thus, a plausible generalisation of these quantities is

(2.5) \begin{equation} \textrm{Re} = \frac {u_0/L}{\left .\gamma _U^-\right |_{k_0}} \qquad \text{and} \qquad \textrm{Rm} = \frac {u_0/L}{\left .\gamma _B^-\right |_{k_0}}, \end{equation}

where $\left .\gamma _U^-\right |_{k_0}$ and $\left .\gamma _B^-\right |_{k_0}$ are the damping rates of the mass flow and the magnetic field, respectively, taken at the wavenumber corresponding to the box scale $k=k_0\equiv 2\pi /L$ .

In the following, we will show that in the 10-moment system employing the heat flux closure given in (2.4), the damping of the mass flow is governed by the ion closure parameter $k_{0,i}$ and the magnetic field damping by the electron closure parameter $k_{0,e}$ . Specifically, the damping rate of the mass flow scales approximately as $\gamma _U^-\propto k^2/k_{0,i}$ , that is, with the same $k$ -dependence as the diffusive flow damping in a viscous fluid. However, the damping rate of the magnetic field scales approximately as $\gamma _B^-\propto k^4/k_{0,e}$ , which is consistent with a hyperdiffusive damping. As we shall see, such separation of the electron and ion physics requires the mass ratio to be sufficiently large.

To establish how the damping of the mass flows and magnetic fields depend on the closure parameters, we have performed one-dimensional simulations initialised with sinusoidally varying, divergence-free perturbations of the respective quantities. The self-consistent exponential decay of these perturbations was used to fit flow and magnetic decay rates, $\gamma _U^-$ and $\gamma _B^-$ . In particular, for the mass flow decay studies, the only non-vanishing component of the (equal) ion and electron flows are initially $u_{\alpha z} = u_{\alpha z 0} \sin (k x)$ . For the magnetic decay studies, only the electrons have a non-zero initial velocity of the same form with a flow amplitude $u_{ez0}=B_{y0}k/(e\mu _0 n_e)$ , corresponding to a magnetic field perturbation $B_y = B_{y0} \cos (k x)$ . These simulations use the same linear box size as our dynamo simulations (namely, $L$ that defines $k_0=2\pi /L$ used for normalisation), while employing ten times higher resolution for improved accuracy.

Figure 1.

Dependence of damping rates on closure parameters (circle markers, main figure) and wavenumber (diamond markers, inset figure). In both of these scans, star markers correspond to the baseline case. Dashed lines indicate fitted relevant power-law behaviour. (a) Mass flow damping rate as a function of $k_{0,i}/k_0$ and $k/k_0$ . (b) Magnetic field damping rate as a function of $k_{0,e}/k_0$ and $k/k_0$ .

Figure 1 illustrates the numerically obtained dependencies of the damping rates in the vicinity of our baseline parameters (to be detailed in § 2.3). Over certain ranges, both $\gamma _U^-$ and $\gamma _B^-$ are approximately inversely proportional to their respective closure parameter: $\gamma _U^-\propto k_{0,i}^{-1}$ and $\gamma _B^-\propto k_{0,e}^{-1}$ . For the flow damping, this dependence can be understood by a linearised analysis, with the assumptions that in the momentum equation, the $m_in_i\partial _t u_{iz}$ and $\partial _x p_{ixz}$ terms nearly balance, while in the relevant $xz$ component of the pressure tensor equation, the generation of this off-diagonal pressure component by the term $p_{ixx}\partial _x u_{iz}$ is nearly compensated by its destruction by the isotropising closure term $\propto v_{\textrm{th},i}k_{0,i}p_{ixz}$ . This behaviour yields a dependence of the damping rate $\gamma _U^-\propto v_{\textrm {th},i} k^2/k_{0,i}$ , consistent with the scalings found in simulations around the baseline point. A somewhat more involved calculation of the magnetic damping rate, where Maxwell’s equations are also invoked in addition to the fluid ones, yields a similar regime $\gamma _B^-\propto v_{\textrm {th},e} k^4 \delta _e^2/k_{0,e}$ , where $\delta _\alpha =c/\omega _{p\alpha }$ denotes the inertial length and $\omega _{p\alpha }=\sqrt {n_\alpha q_\alpha ^2/(\epsilon _0 m_\alpha )}$ the plasma frequency of species $\alpha$ . We note that this differs from the exact kinetic damping of magnetic fields (Mikhailovskii Reference Mikhailovskii1980; Pusztai et al. Reference Pusztai, Juno, Brandenburg, TenBarge, Hakim, Francisquez and Sundström2020), that is, of the form $\gamma _B^-\propto v_{\textrm {th},e} k^3 \delta _e^2$ ; thus, our closure introduces an extra $k/k_{0,e}$ factor.

Clearly, the $1/k_{0,\alpha }$ scaling cannot hold for arbitrarily low values of the closure coefficients $k_{0,\alpha }$ , because at low values, the damping rates saturate at a finite value. This behaviour can be seen for $\gamma _B^-$ in figure 1(b). In fact, box-scale magnetic perturbations at our baseline parameters are in the transition region to this saturation (while higher $k$ values are typically in the $\gamma _B^-\propto k_{0,e}^{-1}$ regime). However, increasing $k_{0,e}$ brings our simulations away from this saturation, providing a reduced magnetic damping, which is essential for the interpretation of the dynamo results. Concerning flow damping, we observe the damping rate to level off also at high values of $k_{0,i}$ , and our baseline case is in this vicinity. It is of no concern though, as we do not perform scans in the ion closure parameter.

Figure 2.

Mass-ratio dependence of the damping rate scaling exponents, defined through the power law relations $\gamma _U^-\propto k_{0,i}^{\alpha _{U,i}}k_{0,e}^{\alpha _{U,e}}$ , $\gamma _B^-\propto k_{0,i}^{\alpha _{B,i}}k_{0,e}^{\alpha _{B,e}}$ .

The effect of the mass ratio on the scaling of the damping rates is illustrated in figure 2 in terms of the exponents defined by

(2.6) \begin{equation} \gamma _U^-\propto k_{0,i}^{\alpha _{U,i}}k_{0,e}^{\alpha _{U,e}}k^{\alpha _{U,k}}, \qquad {} \gamma _B^-\propto k_{0,i}^{\alpha _{B,i}}k_{0,e}^{\alpha _{B,e}}k^{\alpha _{B,k}}. \end{equation}

We consider a starting point with parameters $k = k_0$ and $k_{0,i} = 13.6 k_0$ , as in our baseline case, and set $k_{0,e} = 13.6 k_0$ to ensure consistency at unit mass ratio. The quantities $k_{0,i}$ , $k_{0,e}$ and $k$ are then perturbed about their baseline values in magnetic-field and mass-flow damping tests to determine the above-mentioned exponents. In all simulations, the electron mass is held fixed while the ion mass is varied.

The exponents characterising the $k$ dependence remain nearly constant over the entire mass-ratio scan, with $\alpha _{U,k}\approx 2$ and $\alpha _{B,k}\approx 4$ , exhibiting maximum deviations of only $1.2\,\,\%$ and $0.05\,\%$ , respectively. As these variations are negligible, the corresponding exponents are omitted from the figure. At unit mass ratio, the effects of ion (or, in this case, positron) and electron closures on the damping of a given quantity are identical, as expected. As the mass ratio increases, $\alpha _{U,i}$ and $\alpha _{B,e}$ approach $-1$ , while $\alpha _{U,e}$ and $\alpha _{B,i}$ approach $0$ . In other words, at sufficiently large mass ratio, magnetic-field damping depends solely on $k_{0,e}$ , whereas mass-flow damping depends solely on $k_{0,i}$ , with cross-dependencies becoming negligible. This decoupling between electron and ion closures is already well converged at the mass ratio used in this study, $m_i/m_e = 100$ . Finally, we note that over the mass-ratio range considered here, the sums $\alpha _{B,i} + \alpha _{U,i}$ and $\alpha _{B,e} + \alpha _{U,e}$ remain close to $-1$ , with maximum deviations of $1.9\,\%$ and $0.8\,\%$ , respectively.

Since the $k_{0,\alpha }\rightarrow 0$ limit of the 10-moment system is not singular, at sufficiently low values of the closure parameters $k_{0,\alpha }$ , these scalings must be broken and the damping rates saturate towards constant values. As seen in figure 1(b), we are in the vicinity of such a saturation concerning $\gamma _B^-(k_{0,e})$ at our baseline parameters; however, the scan in $k_{0,e}$ we will perform moves deeper into the $\gamma _B^-\propto k_{0,e}^{-1}$ regime. Concerning the mass flows, we work in the $\gamma _U^-\propto k^2/k_{0,i}$ regime and we keep $k_{0,i}$ fixed.

2.3. Simulation settings

To be relevant to the fast turbulent dynamo in collisionless astrophysical media, our simulations must have sustained three-dimensional flows. In reality, the flows are the result of gravitational collapse in the presence of dark matter, while our simulations include forced flows of non-gravitating baryonic matter in the plasma state. To produce a fast dynamo at a relatively low value of the critical magnetic Reynolds number, we drive our flows towards a simulation box-scale ( $L$ ) Galloway–Proctor (GP) model flow (Galloway & Proctor Reference Galloway and Proctor1992). More specifically, both the electrons and ions are accelerated by the same spatio-temporally varying acceleration $\boldsymbol{a}_\alpha (\boldsymbol{x},t)=C_f \boldsymbol{U}_{\textrm{GP}}(\boldsymbol{x},t)/t_{\text{str},i}$ , which is proportional to the GP flow,

(2.7) \begin{equation} \boldsymbol{U}_{\textrm {GP}}(\boldsymbol{x},t) = u_0 \left [ \begin{matrix} A \sin (k_0 z+\sin \omega t)+C \cos (k_0 y+\cos \omega t) \\[3pt] A \cos (k_0 z+\sin \omega t) \\[3pt] C \sin (k_0 y+\cos \omega t) \end{matrix} \right ], \end{equation}

where $t_{\text{str},i}=L/v_{\text{th},i}$ is the ion streaming time, $u_0$ is the nominal flow speed, $k_0=2\pi /L$ is the lowest wavenumber fitting the simulation box, the constants $A$ and $C$ determining the stirring geometry are chosen to be unity here, while the constant $C_f$ sets the drive strength. In our baseline case, the representative flow speed is sub-sonic $u_0=0.35 c_s$ , with the sound speed $c_s=\sqrt {T_e/m_i}$ , which defines the outer turnover time $t_{\text{turn}}=L/u_0$ . As this is expected to be the characteristic time scale of dynamo growth, we use $t_{\text{turn}}$ for time normalisation with the nominal flow speed $u_0$ , while noting that the root mean square flow speed varies over time and is in general different from $u_0$ . We also set the stirring frequency of the flow to $\omega =2\pi /t_{\text{turn}}$ .

The simulations use a mass ratio of $m_i/m_e=100$ , and our baseline driving strength of $C_f=2$ and $k_{0,i}=13.6 k_0$ are chosen such that the flows remain subsonic throughout the simulation. As a reference, this value of $k_{0,i}$ corresponds to $\textrm{Re}= 0.535$ , while $k_{0,e}=1 k_0$ corresponds to $\textrm{Rm}= 307$ , according to numerical damping rate tests like those presented in § 2.2, yielding a large magnetic Prandtl number $\textrm{Pm}= \textrm{Rm}/\textrm{Re}=574$ , which is relevant for the ICM. Due to numerical limitations on achievable scale separations, our approach to the $\textrm{Pm}\gg 1$ limit sacrifices the fluid cascade and modelling the system at the viscous cutoff scale. This limitation results in smooth, large-scale flows with a rapidly decaying velocity spectrum, i.e. a form of turbulence that lacks an extended inertial range due to strong damping at small scales. Despite the absence of an extended inertial range, such flows can still sustain dynamo amplification through coherent stretching (Schekochihin et al. Reference Schekochihin, Cowley, Taylor, Maron and McWilliams2004). Moreover, it is worth emphasising that the dissipation in a collisionless system is not expected to be fixed: it evolves as the plasma becomes increasingly magnetised and is further modified by kinetic instabilities.

Unlike particle-in-cell simulations, the continuum solvers of Gkeyll are essentially noise-free, thus we need to prescribe magnetic perturbations from which the Weibel instability can grow. Note that due to the rapid growth of this instability, the fact that we need such a seed magnetic field is not an issue, as in astrophysical settings, the instability would produce macroscopic field strengths from arbitrarily small fluctuations on time scales shorter than those of the dynamo growth. Thus, as an initial condition, we prescribe a magnetic field of the form

(2.8) \begin{equation} B_{i}(t=0)=B_0 \sum _{j\ne i,n} b_{ij,n}\cos [n k_0 (x_j+L\phi _{ij,n})], \end{equation}

with random relative amplitudes $b_{ij,n}$ and phases $\phi _{ij,n}$ , both uniform on the unit interval, and a characteristic magnetic field strength $B_0$ . Here, we only seed the four lowest mode numbers $n$ . For our initial nominal field $B=B_0$ , the thermal electron Larmor radius $\rho _e=v_{\text{th},e} m_e/(e B)$ is $\rho _e=4.4 L$ , thus the system is non-magnetised. The system is non-relativistic, since $v_{\text{th},e}/c=0.0625$ .

A current consistent with this field, $\boldsymbol{j}=\mu _0^{-1}\boldsymbol{\nabla }\times \boldsymbol{B}$ , fully carried by the electrons, is also prescribed. In addition, there is an initial flow velocity, equal for the electrons and the ions, given by (2.7) for $t=0$ . The initial ratio of the box-integrated magnetic energy and kinetic energy in the flows is $3.7\times 10^{-6}$ . Our baseline box size is $L=70 \delta _i$ in each direction. With a spatial resolution of $N=140$ in each direction, such that the linear size of the finite volume is $\varDelta x=L/140$ , the ion inertial length is barely resolved $\delta _i=2\varDelta x$ , while the electron inertial length is under-resolved $\delta _e=0.2\varDelta x$ , along with the Debye length $\lambda _D\equiv \sqrt {\epsilon _0 T_e/(n_e e^2)}\approx 0.009\varDelta x$ . These scales are not a concern since the Debye scale is not playing an important role in the physics of the system, and the characteristic scales of the electron Weibel instability are significantly larger than $\delta _e$ due to the weak pressure anisotropy that develops in the system.

We finally note that our simulations are performed in SI units with physical parameters comparable to those of Pusztai et al. (Reference Pusztai, Juno, Brandenburg, TenBarge, Hakim, Francisquez and Sundström2020). Namely, all physical constants have their actual value, except the ion mass that is $m_i=100 m_e$ ; both species have initially the same scalar temperature $T_e=T_i=1\,\textrm {keV}$ , the electron density is $n_e=2.3 \times 10^{28} \,\textrm {m}^{-3}$ , along with $k_0=2.58 \times 10^5\,\textrm {m}^{-1}$ . Furthermore, $B_0=1\,\textrm {T}$ . As the system has the scale invariance of a collisionless kinetic system, the fact that the plasma parameters (particularly the density) are not representative of the ICM is of no concern; we will present our results with practical normalisations.

2.4. Electron Weibel instability in non-magnetised collisionless plasma turbulence

The Weibel instability is a non-magnetised kinetic plasma instability that arises in the presence of an anisotropic velocity distribution. Transverse perturbations in the magnetic field deflect particles in a way that produces currents reinforcing the perturbation. As a result, it tends to generate filamentary current and magnetic field structures, and acts to restore pressure isotropy. It plays important roles in astrophysical and laboratory plasmas (e.g. collisionless shocks, gamma-ray bursts, laser-plasma experiments) by generating sizable magnetic fields in initially non-magnetised or weakly magnetised environments, thereby qualitatively affecting particle dynamics. The instability has both electron and ion variants, and additionally arises in pair-plasmas.

A large body of the Weibel literature is concerned with interpenetrating (often supersonic and relativistic) flow scenarios, due to its relevance for collisionless astrophysical shocks (Medvedev & Loeb Reference Medvedev and Loeb1999; Medvedev, Silva & Kamionkowski Reference Medvedev, Silva and Kamionkowski2006; Spitkovsky Reference Spitkovsky2007; Kato & Takabe Reference Kato and Takabe2008), such as supernova remnants, and laser-plasma experiments (Schoeffler et al. Reference Schoeffler, Loureiro, Fonseca and Silva2014; Huntington et al. Reference Huntington2015). Pressure anisotropies appearing in such scenarios are typically substantial. However, as demonstrated by Zhou et al. (Reference Zhou, Zhdankin, Kunz, Loureiro and Uzdensky2022), the Weibel instability also naturally occurs in non-magnetised collisionless plasma turbulence, since phase mixing in sheared flows creates pressure anisotropy $\varDelta _\alpha$ . For subsonic turbulence, as relevant for our topic, this anisotropy is small (see Appendix B for representative values from our simulations).

The intuitive and conventional definition $\varDelta _\alpha =p_{\|\alpha }/p_{\perp \alpha }-1$ implicitly assumes a pressure tensor with diagonal components $\{p_{\|\alpha }, p_{\perp \alpha },p_{\perp \alpha }\}$ ; this definition will need to be revised, since the pressure tensor in a non-magnetised system is non-gyrotropic and has, in general, three distinct eigenvalues. Thus, following Zhou et al. (Reference Zhou, Zhdankin, Kunz, Loureiro and Uzdensky2022), we use the definition $\varDelta _\alpha =\sqrt {\langle (p_{\textrm{max},\alpha }/p_{\textrm {min}, \alpha} )^2\rangle }-1$ , where $p_{\textrm {max},\alpha}$ is the highest eigenvalue of the pressure tensor, $p_{\textrm {min},\alpha}$ is the average of the two lower eigenvalues and the angle bracket denotes the spatial (box) average.

Since the pressure anisotropies of ions and electrons that develop in the sheared flow tend to be comparable, the electron Weibel instability will dominate the magnetic field growth as it has higher growth rates than the ion mode. Specifically, in terms of the above-mentioned anisotropy measure that is positive by construction, the fastest growing electron Weibel mode has a characteristic wavenumber $k\sim \varDelta _e^{1/2}/\delta _e$ and a growth rate of $\gamma _B \sim \varDelta _e^{3/2}\omega _{pe}v_{\text{th},e}/c$ . For $\varDelta _e \ll 1$ , the characteristic wavenumber is shifted to scales larger than the electron inertial length.