Temperature-dependent Saturation of Weibel-type Instabilities in Counter-streaming Plasmas

TL;DR

This study uses advanced simulations to show that the temperature of counter-streaming plasmas significantly affects the saturation and magnetic field generation of Weibel-type instabilities, impacting astrophysical and laboratory plasma magnetization.

Contribution

First 2X2V continuum Vlasov-Maxwell simulations reveal temperature-dependent saturation effects in Weibel instabilities in counter-streaming plasmas.

Findings

Colder plasmas exhibit significantly reduced temperature anisotropy.

Magnetization levels are orders of magnitude lower in colder plasmas.

Energy conversion efficiency depends strongly on plasma temperature.

Abstract

We present the first 2X2V continuum Vlasov-Maxwell simulations of interpenetrating, unmagnetized plasmas to study the competition between two-stream, Oblique, and filamentation modes in the weakly relativistic regime. We find that after nonlinear saturation of the fastest-growing two-stream and Oblique modes, the effective temperature anisotropy, which drives current filament formation via the secular Weibel instability, has a strong dependence on the internal temperature of the counter-streaming plasmas. The effective temperature anisotropy is significantly more reduced in colder than in hotter plasmas, leading to orders of magnitude lower magnetization for colder plasmas. A strong dependence of the energy conversion efficiency of Weibel-type instabilities on internal beam temperature has implications for determining their contribution to the observed magnetization of many astrophysical and laboratory plasmas.

Equations

$$\frac{\partial f_s}{\partial t}+\vec{v}\cdot \frac{\partial f_s}{\partial \vec{x}}+\frac{q_s}{m_s}(\vec{E}+\vec{v}\times \vec{B})\cdot \frac{\partial f_s}{\partial \vec{v}}=0,$$

$$f_{0,e}(v_x,v_y)=\frac{1}{2\pi v_{th}^2}{e}^{-\frac{v_x^2}{2v_{th}^2}}\left[e^{-\frac{(v_y-u_d{)}^2}{2v_{th}^2}}+e^{-\frac{(v_y+u_d{)}^2}{2v_{th}^2}}\right],$$

$$D_{ij}=\frac{{\omega }^2}{c^2}\left(k_i{k}_j-k^2{\delta }_{ij}\right)+{ϵ}_{ij},$$

$${ϵ}_{ij}=\left(1-\sum _s\frac{{\omega }_{p,s}^2}{{\omega }^2}\right){\delta }_{ij}+\sum _s\frac{{\omega }_{p,s}^2}{{\omega }^2}{\int }_{-\infty }^{\infty }{v}_i{v}_j\frac{\vec{k}\cdot {\nabla }_{\vec{v}}{f}_{0,s}}{\omega -\vec{k}\cdot \vec{v}}{d}^{2}v.$$

$$D=\left(\begin{array}[]{cc}D_{11}&D_{12}\\ D_{21}&D_{22}\end{array}\right)$$

$$D_{11}=1-\frac{{\omega }_{pe}^2}{4k^2{v}_{th}^2}\left(Z^{\prime }({\xi }_{+})+Z^{\prime }({\xi }_{-})\right)$$

$$D_{12}=D_{21}=\frac{{\omega }_{pe}^2{u}_d\cos \theta }{4\omega k{v}_{th}^2}\left(Z^{\prime }({\xi }_{+})-Z^{\prime }({\xi }_{-})\right)$$

$$D_{22}=1-\frac{{\omega }_{pe}^2}{{\omega }^2}-\frac{k^2{c}^2}{{\omega }^2}-\frac{{\omega }_{pe}^2\left(u_d^2{\cos }^2\theta +v_{th}^2\right)}{4{\omega }^2{v}_{th}^2}\left(Z^{\prime }({\xi }_{+})+Z^{\prime }({\xi }_{-})\right),$$

$$B_z(t=0)=\sum _{n_x,n_y=0}^{16,16}{\tilde{B}}_{n_x,n_y}\sin (\frac{2\pi n_{x}x}{L_x}+\frac{2\pi n_{y}y}{L_y}+{\tilde{\varphi }}_{n_x,n_y}).$$

$$\bar{A}={\int }_0^{L_y}{\int }_0^{L_x}\frac{\int (v_y-{\bar{v}}_y{)}^{2}f(x,y,\vec{v})\mathrm{d}\vec{v}}{\int (v_x-{\bar{v}}_x{)}^{2}f(x,y,\vec{v})\mathrm{d}\vec{v}}\mathrm{d}x\mathrm{d}y,$$

$$f_{0,e}(v_x,v_y)=\frac{m_e}{2\pi \sqrt{k_B{T}_x{k}_B{T}_y}}{e}^{-\frac{-m_e{v}_x^2}{2k_{B}Tx}-\frac{m_e{v}_y^2}{2k_{B}Ty}},$$

Full text

††slugcomment: (Received December 16, 2018; Revised February 5, 2019; Accepted February 6, 2019 )

1. Introduction

Electromagnetic instabilities of the Weibel-type have been extensively studied as mechanisms for generating the sub-equipartion, long-lived magnetic fields whose effects are inferred or observed in laser-driven plasmas (Fox et al., 2013; Huntington et al., 2015), gamma ray bursts (Medvedev & Loeb, 1999), pulsar wind outflows (Kazimura at al., 1998), and cosmological scenarios (Schlickeiser & Shukla, 2003; Lazar et al., 2009). Interpenetrating, unmagnetized plasma flows common to the initial conditions of these systems support a diversity of competing instabilities spanned by electrostatic two-stream (TS) and electromagnetic Oblique and filamentation modes (Bret, 2009). Obtaining the efficiency $\epsilon_{B}=E_{M}/E_{K}$ of the conversion of initial kinetic energy $E_{K}=\langle\frac{1}{2}\rho v^{2}\rangle$ into magnetic energy $E_{M}=\langle B^{2}/2\mu_{0}\rangle$ by all of the simultaneously present instabilities is critical for estimating the role that Weibel-type generated magnetic fields play in both laboratory and astrophysical contexts.

Magnetic field generation in counter-streaming plasmas is associated with the filamentation instability (FI; Fried (1959)), while in plasmas with anisotropic bi-Maxwellian velocity distributions, it is attributed to the Weibel instability (WI; Weibel (1959)). However, the secular WI can also be present in counter-streaming plasmas after saturation of a faster growing TS mode relaxes the initial plasma velocity distribution to an anisotropic bi-Maxwellian (Schlickeiser & Shukla, 2003). In this Letter, we refer to these modes as Weibel-type instabilities, e.g. the FI and WI.

All of these modes grow on plasma frequency time scales and therefore require kinetic numerical simulations to study their evolution and competition in the nonlinear regime. Both subrelativistic and relativistic PIC simulations in 2X and 3X (Fonseca et al., 2003; Silva et al., 2003; Sakai et al., 2004; Nishikawa et al., 2003, 2005; Kato & Takabe, 2008; Takamoto et al., 2018) have shown promising efficiencies of $0.25\%\lesssim\epsilon_{B}\lesssim 2\%$ in the long-term quasi-steady state. Typical simulations proceed first forming electron Alfvén limited current filaments on small scales that coalesce into larger filaments, followed by a similar process for ions (Frederiksen et al., 2004; Hededal et al., 2004). We note that recent 2X PIC simulations have found the merging process to disrupt the coherent filament structures, reducing $\epsilon_{B}$ (Kumar et al., 2015; Takamoto et al., 2018). Filament merging and break-up events along with various electron trapping and acceleration mechanisms have been found to heat electrons, possibly explaining the origin of nonthermal radiation observed in several astrophysical scenarios.

Continuum Vlasov-Maxwell codes are an alternative approach that has historically been lucrative in the 1X2V setting (Califano et al., 1998), confirming magnetic trapping as the mechanism responsible for saturation (Davidson et al., 1972) and finding $\epsilon_{B}\sim 1\%$ for weakly relativistic drift velocities $0.1<u_{d}/c<0.6$. However, continuum Vlasov-Maxwell methods have not been employed as widely as the PIC method due to the expense in directly discretizing the distribution function on a high dimensional phase space grid. Nevertheless, advancements in algorithms and computing power have led to a resurgence in the application of continuum Vlasov-Maxwell codes to problems requiring kinetic modeling, motivated by the advantages of the direct-discretization noise-free representation of the distribution function. For example, Cagas et al. (2017) found that for colder beams, electrostatic potentials due to secondary streaming induced by the FI play an equally dominant role as magnetic trapping in saturation of the 1X FI.

Using the same method as Cagas et al. (2017), this Letter presents the first 2X2V continuum Vlasov-Maxwell simulations of interpenetrating plasma flows aimed at incorporating the competition of the two-stream and Oblique modes with the FI. We find that the magnetization parameter $\epsilon_{B}$ is strongly dependent on the internal temperature of the beams in the weakly relativistic regime, wherein the growth rates of the various modes are comparable. Saturation of TS and Oblique modes in colder beams tends to scatter the anisotropy in the particle velocity distribution that drives the later formation of the current filaments via the secular WI. This behavior results in a complete lack of magnetization, $\epsilon_{B}\lesssim 10^{-5}$, in the nonlinear saturation phase of colder beams and typically reported values of $\epsilon_{B}\sim 10^{-2}$ for hotter beams.

2. Problem Setup

We work in the framework of the Vlasov-Maxwell system. The Vlasov equation for species $s$ is given by

[TABLE]

where $f_{s}$=$f_{s}(x,y,v_{x},v_{y},t)$ is the species distribution function and $\vec{E}=E_{x}(x,y)\hat{x}+E_{y}(x,y)\hat{y}$ and $\vec{B}=B_{z}(x,y)\hat{z}$ are the electric and magnetic fields that are evolved using Maxwell’s equations.

The initial plasma flow is modeled as two uniform density electron streams with opposite flow speeds $u_{d}$ in the y direction, allowing study of the purely reactive ($\omega=\mathrm{Im}[\omega]>0$) form of the FI. The beams are neutralized by a uniform background of ions. The distribution function at $t=0$ for electrons is given by

[TABLE]

where $v_{th}=\sqrt{k_{B}T_{e}/m_{e}}$.

The counter-streaming motion provides the free energy source that drives the TS, Oblique, and filamentation instabilities, which convert kinetic energy into electric and magnetic energy. TS and FI modes generate primarily electrostatic and magnetic energy, respectively, while Oblique modes qualitatively develop as a mixture of the two instabilities.

To study this problem in the nonlinear regime, we perform continuum kinetic simulations in two periodic configuration space dimensions and two velocity dimensions (2X2V) using the Gkeyll framework (Juno et al., 2018). Gkeyll uses the discontinuous Galerkin method for spatial discretization and a strong-stability preserving Runge-Kutta method for time stepping to solve the full Vlasov-Maxwell system in the non-relativistic regime.

3. Linear Theory

We define the angle $\theta$ of the wave vector $\vec{k}$ from the x-axis so that $\theta=0^{\circ}$ corresponds to pure FI modes, $\theta=90^{\circ}$ to pure TS modes, and any angle in between to Oblique modes. Linearizing the Vlasov-Maxwell system gives the general dispersion matrix

[TABLE]

with

[TABLE]

We find that using $\vec{k}$ aligned coordinates gives the simplest representation of $D$, i.e., rotating the original coordinates by angle $\theta$. Substitution of $f_{0,e}$ into the general dispersion matrix and neglecting the ion contribution then gives

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $Z(\xi_{\pm})$ is the plasma dispersion function with $\xi_{\pm}=\frac{\omega\pm ku_{d}\sin\theta}{\sqrt{2}kv_{th}}$. Eigenmodes of the system correspond to solving $\text{det}(D)=0$ for $\omega$ with corresponding eigenvectors satisfying $R^{T}DR\vec{E}=0$, where $R$ is the matrix for rotation by angle $-\theta$.

At any given $v_{th}$ and $u_{d}$, there will be a fastest-growing mode at some $k=|\vec{k}|$ and $\theta$. As seen in Figure 1, the fastest-growing mode transitions from $\theta=0^{\circ}$ to $\theta=90^{\circ}$ in the weakly relativistic region $0.1c\lesssim u_{d}\lesssim 0.5c$. The TS growth rate dominates that of the FI for $u_{d}\lesssim 0.1c$, because the FI is an electromagnetic instability, generally growing by a factor of $O(u_{d}/c)$ slower. While for $u_{d}>0.5c$, the FI growth rate dominates; however, at these flow velocities relativistic effects become significant and the non-relativistic dispersion relation from equation (5) provides only qualitative results.

In the transition region, most of the modes have growth rates roughly within a factor of 3-4 of each other. While the fastest-growing mode will reach saturation first, the other modes can continue to grow approximately at their linear rates, because the nonlinear saturation stage of the dominant mode may not have had enough time to significantly alter the property of the initial distribution function driving the instability. For example, TS saturation may scatter the distinct counter-streaming motions that drive the FI, but an effective temperature anisotropy capable of driving the secular WI may remain. This competition leads to interaction between the nonlinear stages of many of the modes, hence the final steady state can be different than simply the steady state of the fastest-growing mode alone.

4. Simulation Results

To study the simultaneous competition of TS, Oblique, and FI modes, we initialize a bath of electric and magnetic fluctuations with random amplitudes and phases, e.g.,

[TABLE]

$E_{x}(t=0)$ and $E_{y}(t=0)$ are initialized similarly, and they are given equal average energy densities, $\langle\epsilon_{0}E_{x}^{2}/2\rangle=\langle\epsilon_{0}E_{y}^{2}/2\rangle=\langle B_{z}^{2}/2\mu_{0}\rangle\approx 10^{-7}E_{K}$.

We present the results of four simulations, shown by the red crosses in Figure 1, where we fix the drift velocity at $u_{d}=0.1c$, but vary the temperature of the beams by choosing $v_{th}/u_{d}\in\{0.1,0.175,0.25,0.5\}$. The box sizes, respectively, are $L_{x}/d_{e}\in\{2.7,3.8,4.4,7.7\}$ and $L_{y}/d_{e}\in\{3.1,4.0,4.8,6.3\}$. Box sizes $L_{x}=2\pi/k_{FI}^{max}$ and $L_{y}=2\pi m/k_{TS}^{max}$ are chosen to be roughly equal $L_{x}\approx L_{y}$ while fitting a single fastest-growing wavelength of FI and an integer number, $m\approx k_{TS}^{max}/k_{FI}^{max}$, of TS wavelengths. All simulations used the Serendipity polynomial basis of order 2 for phase space discretization (Arnold & Awanou, 2011; Juno et al., 2018). The evolution of the magnetic

($\epsilon_{B}$) and electric ($\epsilon_{E}$) energy densities normalized to the initial kinetic energy is shown in Figure 2 as solid lines for 2X2V with random initial modes and as dashed lines for 1X2V with pure 1X FI initial modes at the same beam parameters.

In all 2X2V simulations, a combination of TS and multiple Oblique modes exponentially grow the fastest during the linear phase, as predicted by the linear growth rates in Figure 1, and then are the first to saturate. The presence of electromagnetic Oblique modes is indicated in Figure 2 by the exponential growth of both electric and magnetic energy, as pure TS modes would appear only as growing electric energy. The simultaneously present FI modes grow much more slowly, and during the linear phase are comparable to the pure FI modes in 1X2V shown as dashed lines in Figure 2. Following saturation, the nonlinear interaction of potential wells formed by the saturation of two-stream and Oblique modes, the tilted current filaments of Oblique modes, and the vertical current filament formation associated with the potentially still-growing FI become dominant.

The nonlinear dynamics forming or disrupting current filaments are strongly temperature dependent. Figure 3 shows spatial snapshots of $E_{y}$ and $B_{z}$ for the (upper panels) hot, $v_{\mathrm{th}}/u_{d}=0.5$, and (lower panels) cold, $v_{\mathrm{th}}/u_{d}=0.1$, cases during and long after the time of saturation. In the post-saturation phase of the hot case, the horizontally varying $B_{z}$ associated with vertical current filaments is dominantly visible and steady in time (Figure 3d), while in the cold case $B_{z}$ is spatially disorganized and strongly fluctuating in time (Figure 3h). We find that a qualitatively similar pattern occurs for $u_{d}/c\in\{1/30,0.25,0.4\}$. Note that $u_{d}/c=1/30$ is the upper limit of intergalactic plasma drift velocities relevant to cosmological scenarios and cosmological magnetic field generation (Miniati, 2002; Schlickeiser & Shukla, 2003; Lazar et al., 2009).

In all cases, the electric energy quickly grows, saturates, and then rapidly decreases. In Figure 4 we plot different parts of the phase space before and after the disruption of the electric energy for the same hot and cold cases shown in Figure 3. In the hot case, near the time of saturation, clear 2D electron tubes form (Figure 3a) due to the fastest-growing wavelength of the TS $y$-varying $E_{y}$ extending in the $x$ direction, which appear as 1X electron holes, or BGK modes, in the $Y$ versus $V_{y}$ cut at a fixed $X$ value shown in Figure 4a. The 2D electron tubes are slightly skewed due to the simultaneously growing, albeit slower, Oblique modes. TS saturation creates a plateau in the $V_{y}$ direction of the velocity distribution function (Figure 4d), leaving an effective temperature anisotropy from which the now-secular Weibel instability can grow. In Figure 5, we find that the spatially averaged temperature anisotropy, $\bar{A}$, drops from $\bar{A}=5$ at $t\omega_{pe}=40$ just before saturation to $\bar{A}\approx 3$ at $t\omega_{pe}=100$ soon after, where

[TABLE]

and $\bar{v}_{i}=\int v_{i}f(x,y,\vec{v})\mathrm{d}\vec{v}/\int f(x,y,\vec{v})\mathrm{d}\vec{v}$. On the order of an electron bounce period, $T_{b}$, in the BGK mode after saturation, filaments from the previous FI and now growing secular WI disrupt all of the electron tubes, leading to the observed sudden decrease in stored, what is primarily, electrostatic energy, where $T_{b}\sim 2\times 2\pi/(k_{TS}^{max}u_{d})\sim 20\omega_{pe}^{-1}$. The remaining temperature anisotropy of $\bar{A}\approx 3$ at $t\omega_{pe}=100$ drives further Weibel growth, and visible filaments appear around $t\omega_{pe}=150$ (not shown), where $\epsilon_{B}$ grows again. The magnetic energy then saturates around $\epsilon_{B}\sim 10^{-2}$ near the Alfvén-limited regime, $\rho_{g}\sim m_{e}u_{d}/(eB_{z})\sim 7d_{e}\sim L_{x}$, and enters a steady-state oscillation at the magnetic bounce frequency, agreeing closely with PIC and 1X2V simulations.

In the cold case, many Oblique modes with $\theta\gtrsim 30^{\circ}$ have the same or faster growth rates than TS and therefore saturate at similar times, creating 2D electron tubes with far more spatial structure and with tilted current filaments passing through (Figure 3e and 4e). The phase space has sharp structures with the main two beams rotating in the $V_{x}-V_{y}$ plane (Figure 4g), followed by the beams fragmenting and developing fine structure, as shown in the distribution functions in Figure 4h and 4f at much later time. The spatially averaged temperature anisotropy sharply drops (Figure 5) from $\bar{A}=101$ at $t\omega_{pe}=20$ immediately prior to saturation, to $\bar{A}\approx 1.2$ at $t\omega_{pe}=50$ soon after. However, unlike the hot case, $\bar{A}$ remains constant through the rest of the simulation. Filaments do not form after disruption of the BGK modes due to the lack of a significant residual effective temperature anisotropy to drive the secular WI, and the end result is a bath of primarily electrostatic fluctuations. The lack of current filaments does not provide a steady structure to sustain magnetic fields, thus $\epsilon_{B}$ decays to values orders of magnitude smaller than that observed in PIC and 1X2V simulations.

We have doubled and tripled the box sizes and found no qualitative difference in the results. In the hot case, filaments form at the same wavelengths $\lambda_{FI}^{max}=2\pi/k_{FI}^{max}$ as before, and then merge to larger scales. However, the merging process does not break up the filaments as was reported in the 2X case of recent PIC simulations (Kumar et al., 2015; Takamoto et al., 2018), although their beam drift velocities were relativistic, with relativistic $\gamma=5$. In the cold case, we still see no filament formation for larger box sizes. This result contradicts previous non-relativistic PIC simulations of Kato & Takabe (2008), in which they found filament formation and $\epsilon_{B}\sim 10^{-2}$ post-saturation for the same cold beam parameters. Note that in the cold case, our result is physically different than in Kumar et al. (2015) and Takamoto et al. (2018); they see the disruption of filaments after the merging of initial filaments forming at $\lambda_{FI}^{max}$ scales, while we do not observe formation of the $\lambda_{FI}^{max}$ scale filaments in the first place. Inherent particle noise in PIC simulations may potentially make beams effectively hotter than they should be, which weakens the effects in the cold beam regime that diminish the effective temperature anisotropy.

To disentangle the role of TS and Oblique modes versus temperature anisotropy, we also ran 2X2V simulations of the classic Weibel instability with a bi-Maxwellian initial distribution function:

[TABLE]

with $k_{B}T_{y}=m_{e}\left(u_{d}^{2}+v_{th}^{2}\right)$ and $k_{B}T_{x}=m_{e}v_{th}^{2}$ for the same anisotropy parameters $\bar{A}=101$ and $\bar{A}=5$. Current filament formation and merging were unaffected by temperature anisotropy, achieving $\epsilon_{B}\sim 10^{-2}$ in the nonlinear steady state for both cases. This result emphasizes the role that TS and Oblique mode saturation play in affecting the effective temperature anisotropy of the distribution function and in controlling the growth of current filaments via the secular WI.

5. Conclusions

We present the first 2X2V continuum-Vlasov Maxwell simulations of weakly relativistic interpenetrating plasma flows to show a strong dependence of current filament formation and resulting magnetization on the internal temperature of the beams. Hotter beams maintain an effective temperature anisotropy after saturation of the faster growing TS and Oblique modes, allowing for continuing current filament formation via the secular Weibel instability. The result is a final magnetization of $\epsilon_{B}\sim 10^{-2}$ due to current filaments reaching the Alfvén current limited regime and merging to larger scales. On the other hand, the effective temperature anisotropy for colder beams falls rapidly and remains constant after the saturation of TS and Oblique modes. Current filament formation does not follow, and instead, a bath of electrostatic fluctuations dominates with effectively no magnetization, $\epsilon_{B}\lesssim 10^{-5}$.

These results put testable constraints on the efficiency of magnetic field generation via Weibel-type instabilities in astrophysical and laboratory scenarios. Differences in observed magnetization for colliding plasmas of similar flow velocities could potentially be attributed to a difference in the beams’ internal temperatures. We have found this temperature dependence for the range $1/30\leq u_{d}/c\leq 0.4$, where $u_{d}/c=1/30$ corresponds to the upper bound of intergalactic plasma counter-streaming velocities relevant to cosmological scenarios. The source of disagreement between previous PIC simulations, which find filament formation, and our study, which does not, remains to be determined in the case of colder beams. The continuum Vlasov-Maxwell calculations presented in this Letter could be extended to include the evolution of an ion species, as well as the third dimension, while PIC calculations could be repeated with larger numbers of particles-per-cell to further mitigate the noise endemic to the PIC algorithm, even beyond the standard convergence analyses that have been performed in previous studies.

We are grateful for insights from conversations with Petr Cagas and William Fox. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant No. ACI-1548562, resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357, as well as the Eddy cluster at Princeton University. A. H. was supported by the U.S. Department of Energy under Contract No. DE-AC02-09CH11466 and Air Force Office of Scientific Research under grant No. FA9550-15-1-0193; J. J. was supported by a NASA Earth and Space Science Fellowship (Grant No. 80NSSC17K0428); and J.M. TenBarge was support by NSF SHINE award (AGS-1622306).