2. Energy invariants and the approximation of a narrow spectrum

Consider a non-relativistic collisionless plasma with a uniaxial anisotropic distribution of electrons, $f(t, \boldsymbol {r}, \boldsymbol {v})$, where the effective temperature along the $z$-axis is higher than in the transverse directions, in which temperatures are equal: $T_z > T_x = T_y$ (i.e. axial symmetry, including the case of a uniform external magnetic field applied along $z$). Here the temperature $T_j = m_\mathrm {e} \langle v_j^2 \rangle - m_\mathrm {e} \langle v_j \rangle ^2$ for each Cartesian component $j = x, y, z$; angular brackets denote averaging in velocity space: $\langle \cdots \rangle = (n_0 / n) \int \cdots f \,{\rm d}^3 \boldsymbol {v}$, where $n = n_0 \int f \, {\rm d}^3 \boldsymbol {v}$ stands for plasma number density, $n_0 = n(t = 0)$ its initial value and $m_\mathrm {e}$ is the electron mass. According to the Vlasov–Maxwell equations with massive ions (considered below to be immobile), such an anisotropic distribution of electrons is subject to the Weibel instability (Weibel Reference Weibel1959; Davidson Reference Davidson1989) that leads to aperiodic growth of magnetic field perturbations with wavevectors $\boldsymbol {k}$ oriented, in general, at various angles to the anisotropy $z$ axis (see, e.g., Vagin & Uryupin Reference Vagin and Uryupin2014; Silva, Afeyan & Silva Reference Silva, Afeyan and Silva2021). If the plasma is initially homogeneous, $n_0 = \textrm {const.}$, and only the transverse modes of an extraordinary type are excited, in which the projection $k_z = 0$ and the electric field is directed along the $z$ axis, then the plasma remains homogeneous along this axis at all times, so that the projection of the generalized momentum onto the $z$ axis, $P_z$, is conserved on particle trajectories. It can be shown that its square, $P_z^2$, averaged over velocities and the orthogonal plane $xy$, is also independent on time. This condition, together with the conservation of the total energy of particles and fields, allows one to formulate two invariants first introduced by Davidson & Hammer (Reference Davidson and Hammer1971).

Disregarding the potential electric field and the change in plasma number density ($n = \textrm {const.}$), as well as the displacement current, which, for a wide range of plasma parameters, in the Weibel modes appears to be small compared with the current of charged particles, we can write these invariants in the following form:

(2.1)\begin{gather} \left[ \frac{n T_z}{2} + n \frac{ m_\mathrm{e} }{2} V_z^2 \right] + \sum\limits_{k_x,k_y} \frac{| B_{k} |^2}{8{\rm \pi}} \left( 2 + \frac{ \omega_\mathrm{p} ^2}{c^2} k^{-2} \right) = \textrm{const.} , \end{gather}

(2.2)\begin{gather}\left[ \frac{n T_x}{2} + \frac{n T_y}{2} \right] - \sum\limits_{k_x,k_y} \frac{| B_{k} |^2}{8{\rm \pi}} \left( 1 + \frac{ \omega_\mathrm{p} ^2}{c^2} k^{-2} \right) = \textrm{const.} \end{gather}

Here $V_z (t, x, y) = \langle v_z \rangle$ is the mean velocity of electrons along the $z$ axis, square brackets denote averaging in the plane $xy$, $[ \cdots ] = S^{-1} \iint \cdots \, {{\rm d} x} \, {{\rm d} y}$, $\omega _\mathrm {p} = (4 {\rm \pi}e^2 n / m_\mathrm {e} )^{1/2}$ is the plasma frequency, $e$ is the elementary charge, $c$ is the speed of light in vacuum, and we introduced the discrete modes (Fourier harmonics) $\boldsymbol {B}_{\boldsymbol {k}} = S^{-1} \iint _S \! \boldsymbol {B} \exp (- i k_x x - i k_y y) \, {{\rm d} x} \, {{\rm d} y}$ of the magnetic field $\boldsymbol {B}(t,x,y)$, defined in the simulation domain $S$ with periodic boundary conditions. The spatial spectrum of the magnetic field growing from the noise is averaged over the azimuthal angle in the plane $xy$ in accordance with the axial symmetry of the initial distribution function, so that the modes $\boldsymbol {B}_k$ are characterized only by the radial wavenumber $k$, although the summation in (2.1) and (2.2) takes into account all directions of the wavevectors. It is also presumed that, during the instability, the mean particle velocity in the plane $xy$ remains zero, $\langle v_{x,y} \rangle = 0$.

In the derivation of (2.1) and (2.2) the homogeneity of the number density and non-relativistic character of the plasma are used in the approximate equality $[ a_z^2 \langle n_0 \gamma ^{-1} \rangle ] \approx n [ a_z^2 ] = n \sum k^{-2} |B_k|^2$, where $a_z (t, x, y)$ is the vector potential, $\gamma (v)$ is the Lorentz factor of electrons. Relations similar to (2.1) and (2.2) for the case of a relativistic plasma with the temperatures $T_x = T_y > T_z$ and the modes with wavevectors along the $z$ axis, when the electric field is orthogonal to this axis, are given in Yang, Arons & Langdon (Reference Yang, Arons and Langdon1994). These relations for whistlers in a non-relativistic plasma in the same geometry are considered in Ossakow, Ott & Haber (Reference Ossakow, Ott and Haber1972).

Making use of the relation between the mode amplitudes of the magnetic field and the mean velocity of electrons, $k c | B_{k} | \approx 4 {\rm \pi}e n | V_{z k} |$, that follows from the Maxwell equation disregarding the displacement current and the contribution of ions to the Weibel currents, from (2.1) we have

(2.3)\begin{equation} \left[ \frac{n T_z}{2} \right] + \sum\limits_{k_x,k_y} \frac{| B_{k} |^2}{8{\rm \pi}} \left( \frac{c^2}{ \omega_\mathrm{p} ^2} k^2 + 2 + \frac{ \omega_\mathrm{p} ^2}{c^2} k^{-2} \right) = \textrm{const.} \end{equation}

Despite the fact that during the development of the thermal anisotropy-driven Weibel instability many different modes govern the dynamics of the total field (Stockem et al. Reference Stockem, Dieckmann and Schlickeiser2010; Borodachev et al. Reference Borodachev, Garasev, Kolomiets, Kocharovsky, Martyanov and Nechaev2017), let us assume that at every moment of time the energy is concentrated only in a narrow spectral range of wavevector lengths near $\tilde {k}$ and replace each summation by a single term. (This assumption is justified by simulations, in particular presented in § 3.) Then from (2.2) and (2.3) we finally obtain the following expression for the spatially averaged energy density of the magnetic field $w_B = \sum _{ { \tilde {k}_x,\tilde {k}_y }} | B_{ \tilde {k} } |^2 / (8 {\rm \pi})$:

(2.4)\begin{equation} \frac{ {2} w_B (t)}{[ n_0 T_{z 0} ]} = \frac{ {2} \tilde{k} ^2}{1 + \tilde{k} ^2} \frac{( A_{0} + 1 )^{-1} ( A_{0} - A )}{A + 2 \tilde{k} ^2 + 3} \leq \frac{{2}}{2 \sqrt{6} + 5} \frac{A_{0}}{A_{0} + 1} . \end{equation}

Here $\tilde {k} = k c / \omega _\mathrm {p}$ is the dimensionless wavenumber characterizing the spectrum. When the spectrum is narrow, $\tilde {k}$ can be chosen as corresponding to the spectral peak. For a broader spectrum it is reasonable to define $\tilde {k}$ as an averaged inverse square, so that $\tilde {k} ^{-2} = (8 {\rm \pi}w_B)^{-1} \sum _{k_x,k_y} k^{-2} | B_{k} |^2$. This definition naturally follows from (2.2) and (2.3) and performs surprisingly well (see § 3), probably because the term $\propto k^{-2}$ is principal in the sums at the late stage of the instability when the spectrum shifts towards smaller $k$. Next, $A = [ T_z ] / [ T_\perp ] - 1$ is the average plasma anisotropy parameter, and the quantities with index $0$ refer to the initial time moment. Although it was assumed that initially $T_x = T_y$, in the definition of the anisotropy parameter $A$ we use the effective transverse temperature $T_\perp = (T_x + T_y) / 2$, keeping in mind possible generalizations of the relation (2.4) to a case when the axial symmetry is absent.

Note that the obtained relation is valid for any initial values of the anisotropy parameter and arbitrary distribution functions that could change in a complex way at the nonlinear stage of the Weibel instability – as long as such a change does not violate the system homogeneity along the $z$ axis and is slow enough so that the displacement current remains small.

The maximum of the left-hand side of the inequality (2.4) is reached at $\tilde {k} ^4 = (A + 3)/2$, which corresponds to the implicit optimization in terms of the distribution function that provides the greatest transfer of the thermal energy of electron longitudinal motion into the turbulent magnetic field energy at a given decrease in the anisotropy parameter from $A_0$ to $A$. The absolute maximum is reached at a complete isotropization, i.e. $A = 0$, and is written on the right-hand side of the inequality (2.4). Therefore, the inequality shows that regardless of the initial particle distribution, even with an infinitely high anisotropy parameter, $A_{0} \rightarrow \infty$, the Weibel instability of transverse modes cannot provide a magnetization level, $2 w_B / [ n_0 T_{z 0} ]$, higher than approximately $0.2$.

In the equations above we neglected the contribution of ions to the Weibel currents and the change in the ion kinetic energy, which for a multicomponent plasma is actually valid only on a limited time interval (see, e.g., Borodachev et al. Reference Borodachev, Garasev, Kolomiets, Kocharovsky, Martyanov and Nechaev2017). However, the account for the presence of ions (or other electron fractions) is possible within the framework of the approximations made and leads to the relation

(2.5)\begin{equation} \frac{w_B (t)}{[ K_{z0} ]} = \frac{2 \tilde{k} ^{2} ( \tilde{A} _{0} + 1 )^{-1} ( \tilde{A} _{0} - \tilde{A} )}{ \tilde{A} (1+ \tilde{k} ^{2}) + 5 \tilde{k} ^{2} + 3} ,\end{equation}

where $K_{x,y,z} = 0.5 \sum _\alpha n_\alpha m_\alpha \langle v_{x,y,z}^2 \rangle _\alpha$ is the energy density of the kinetic motion of all particles along the $x$, $y$ and $z$ axes, respectively, $n_\alpha$ and $m_\alpha$ are the number density and masses of particles of a sort $\alpha$, angular brackets denote averaging in velocity space over the distribution function of particles of a sort $\alpha$, and the anisotropy parameter is redefined as $\tilde {A} = 2 [K_z] / [K_x + K_y] - 1$.

Derivation and analysis of the relation (2.5) or its analogues for the problems with different geometry of anisotropy under the assumption of a narrow spatial spectrum of turbulence, created by currents of all particles, will be given elsewhere. As an example, let us point out the case of a disk-type anisotropy of the distribution function, when $T_x = T_y > T_z$ and the Weibel instability is determined by the modes with wavevectors parallel to the $z$ axis; then the equalities in (2.4) and (2.5) remain valid, if we replace $A$ (or $\tilde {A}$) with $4 A + 3$ (and similarly $A_0$ or $\tilde {A} _{0}$).

3. Numerical simulation results and verification of the analytical relation

We verify the relation (2.4) by particle-in-cell simulations with the kinetic relativistic code EPOCH (Arber et al. Reference Arber, Bennett, Brady, Lawrence-Douglas, Ramsay, Sircombe, Gillies, Evans, Schmitz, Bell and Ridgers2015) in two-dimensional geometry for different sets of plasma parameters. Let us demonstrate the validity of the relation by three examples. For all three simulations we choose as the initial condition a spatially homogeneous bi-Maxwellian distribution function of electrons, $f (t = 0, \boldsymbol {v}) \propto \prod _{j = x,y,z} \exp [ - m_\mathrm {e} v_j^2 / (2 T_{j 0}) ]$. The initial electron number density is equal to $n_0 = 10^6$ cm$^{-3}$, ions are immobile. For the simulation (a) we use the temperatures $T_{z0} = 11$ keV and $T_{\perp 0} = 1$ keV, so that the initial anisotropy parameter is $A_0 = 10$, and external fields are absent. For the simulation (b), $T_{z0} = 3$ keV and $T_{\perp 0} = 2.7$ keV, $A_0 = 10$, and the external magnetic field $B_{0z} = 0.53$ G along the $z$ axis is applied. For the simulation (c), $T_{z0} = 3.4$ keV and $T_{\perp 0} = 2.7$ keV, $A_0 = 0.25$, external fields are absent. Two-dimensional simulation domain $S$ is represented by a square with periodic boundary conditions and a side length of (a) $800 r_ {\mathrm {e}}$ divided into $800 \times 800$ cells with $1000$ particles per cell, (b) $1600 r_ {\mathrm {e}}$ divided into $800 \times 800$ cells with $400$ particles per cell, (c) $1600 r_ {\mathrm {e}}$ divided into $400 \times 400$ cells with $1200$ particles per cell. Here $r_ {\mathrm {e}} = T_{z0}^{1/2} ( 8 {\rm \pi}e^2 n_0 )^{-1/2}$ is the Debye length, and time is normalized to $2 {\rm \pi}/ \omega _\mathrm {p}$ below.

Figure 1 shows the evolution of the normalized magnetic energy density, $w_B / ( n_0 T_{z 0} )$, given by the relation (2.4), in which the evolving anisotropy parameter $A$ and the indicated wavenumber $\tilde {k}$ from the simulation are substituted. Despite the complex character of the magnetic turbulence dynamics and the difference of the shape of the electron distribution from a bi-Maxwellian one (see below), the analytically found relation agrees well with the simulation data. This agreement is much better than that given by the previous quasiadiabatic theory (Montes & Peyraud Reference Montes and Peyraud1972; Lemons, Winske & Gary Reference Lemons, Winske and Gary1979) assuming the self-similar evolution of the initial bi-Maxwellian particle distribution function and resulting in an unreasonable substitution $\tilde {k} ^2 = A$ in the relation (2.4). Small discrepancies of approximately ${\sim }5$ % between the expression (2.4) and the simulation with the initial anisotropy $A_0 = 0.25$ are apparently related to the insufficient narrowness of the spatial spectrum of the magnetic turbulence. Note, that the width of the spectrum divided by the wavenumber of its maximum is greater than $1$ for the entire simulation, but even in that case (2.4) gives a decent approximation and outperforms significantly the quasiadiabatic theory.

Figure 1. Normalized magnetic field energy density, ${2} w_B [ n_0 T_{z 0} ]^{-1} 10^2$, obtained from the simulation (black solid curve) and found using the relation (2.4) analytically (red dashes) in dependence on time. The blue dash–dotted curve corresponds to the quasiadiabatic approximation (Montes & Peyraud Reference Montes and Peyraud1972; Lemons et al. Reference Lemons, Winske and Gary1979), which is equivalent to the substitution $\tilde {k} ^2 = A$ in the relation (2.4). Panels correspond to simulations with (a) $A_0 = 10$ and $B_{0z} = 0$, (b) $A_0 = 10$ and $B_{0z} = 0.53$ G, (c) $A_0 = 0.25$ and $B_{0z} = 0$.

Let us discuss in more detail simulation (a) with $A_0 = 10$, $B_{0z} = 0$. Calculations show that in the course of the instability, the electron currents, emerging in the plasma, flow mainly along the $z$ axis and generate magnetic field with lines lying in the $xy$-plane (see figure 2). At the linear stage of the instability, the electron number density slightly increases (no more than on 5 %) in the regions where these currents are localized and equally decreases in the regions of the magnetic field localization. However, after the magnetic field growth saturates, a noticeable modulation of electron number density vanishes. Meanwhile the electron velocity distribution considerably differs from a bi-Maxwellian and, on average, gradually becomes more isotropic (figure 3), so that the anisotropy parameter decreases monotonically to the value $A \approx 0.3$ at the moment of dimensionless time $t = 160$, when the simulation ends.

Figure 2. Distribution of the electron number density relative modulation, $n / n_0 - 1$, in the simulation domain at the moment of time $t \approx 15$ obtained in the simulation with $A_0 = 10$, $B_{0z} = 0$. Several magnetic field lines are shown in blue. The axes represent the coordinates $x$ and $y$ in units of the Debye length $r_ {\mathrm {e}}$.

Figure 3. Evolution of the electron anisotropy parameter $A = [T_z] / [ T_\perp ] - 1$ (solid curve) and the normalized magnetic field energy density ${2} w_B [ n_0 T_{z 0} ]^{-1} 10^2$ (dashes) obtained in the simulation with $A_0 = 10$, $B_{0z} = 0$.

The spatial spectrum of magnetic field, as expected, has an axial symmetry, therefore the spectrum can be completely described by the distribution of the power spectral density of the field, $|B_k|^2$, over the radial wavenumber $k$, obtained by averaging the original two-dimensional spectrum over wavevector directions in the $xy$-plane. In figure 4(a) this distribution is shown at different moments of dimensionless time: at the beginning ($t = 10$) and in the middle ($t = 15$) of the linear stage of instability; at the stage of its saturation ($t = 40$); and also for several later moments at the nonlinear stage of magnetic field decay. The characteristic time of instability development (inverse linear growth rate) equals approximately $2$.

Figure 4. (a) Power spectrum of magnetic field $|B_k|^2$ (in arbitrary units) in dependence on the normalized radial wavenumber at successive moments of dimensionless time (indicated in the legend) in the simulation with $A_0 = 10$ and no external magnetic field. (b) Comparison of the spectra obtained in the simulations with $A_0 = 0.25$ at the time moment $t = 800$ (blue dash–dotted line), $A_0 = 10$ and $B_{0z} = 0.53$ G at $t = 500$ (red dashes), $A_0 = 10$ and $B_{0z} = 0$ at $t = 80$ (black solid line). The time moments chosen correspond to the nonlinear stage after the saturation.

In all three simulations the Weibel instability spectrum has a prominent peak (the ratio of its width to the position of the maximum is of the order of $1$) justifying the assumption made before (2.4). This is shown in figure 4(b) in which the spectra obtained in different simulations are compared at some late time moment after the saturation of the instability.

According to figure 4(a), the dynamics of magnetic turbulence, even in the considered axially symmetric two-dimensional case, is rather complex. Starting from the moment $t \approx 40$, when the average magnetic field saturates and begins to decay, the evolution of the spectrum becomes self-similar: the power laws for its long- and short-wavelength slopes establish $|B_k|^2 \propto k^{6}$ and $|B_k|^2 \propto k^{-7}$, respectively, and do not change further with time (see figure 4a). Moreover, the growth of long-wavelength modes and the decay of short-wavelength ones at the far nonlinear stage also follow power laws in time, e.g. $|B_k|^2 \propto t^{-5/2}$ for the decaying modes (not shown in this paper). The evolution of the full spectrum is shown in figure 5. It could be seen that the wavenumber $\tilde {k}$, corresponding to the maximum power of magnetic field, decreases with time approximately following a square root law, $\tilde {k} \propto t^{-1/2}$, which was demonstrated for collisionless plasma earlier (see Silva et al. (Reference Silva, Schoeffler, Vieira, Hoshino, Fonseca and Silva2020), Zhou et al. (Reference Zhou, Wu, Loureiro and Uzdensky2021) and also figures 3 and 5 in Borodachev et al. Reference Borodachev, Garasev, Kolomiets, Kocharovsky, Martyanov and Nechaev2017). However, in different problems of Weibel turbulence, various power laws of self-similar spectrum evolution (for the slopes of instantaneous spectra, time-dependence of individual modes or the dominant wavenumber) were reported; see, e.g., Romanov et al. (Reference Romanov, Bychenkov, Rozmus, Capjack and Fedosejevs2004), Stockem et al. (Reference Stockem, Dieckmann and Schlickeiser2010), Nechaev et al. (Reference Nechaev, Garasev, Kocharovsky and Kocharovsky2020) and Kuznetsov et al. (Reference Kuznetsov, Nechaev, Garasev and Kocharovsky2023) for a bi-Maxwellian plasma and Gruzinov (Reference Gruzinov2001), Frederiksen et al. (Reference Frederiksen, Hededal, Haugbølle and Nordlund2004), Medvedev et al. (Reference Medvedev, Fiore, Fonseca, Silva and Mori2005) and Dieckmann et al. (Reference Dieckmann, Lerche, Shukla and Drury2007) for various two-stream systems.

Figure 5. Evolution of the spatial spectrum of magnetic field (shades of grey represent power spectral density in arbitrary units) obtained in the simulation with $A_0 = 10$, $B_{0z} = 0$. Blue line correspond to the power-law function $\propto t^{-1/2}$.