Magnetic Field Saturation of the Ion Weibel Instability in Interpenetrating Relativistic Plasmas

The magnetic field is a fundamental physical quantity found everywhere in the universe. It plays essential roles in galaxy spiral arms (Fletcher et al. 2011), high-energy astrophysical phenomena (Shakura & Sunyaev 1976; Piran 1999; Bühler & Blandford 2014), star formation (Basu & Mouschovias 1994; Machida et al. 2008), and cyclic dynamos at the interior of stars and planets (Glatzmaier et al. 1999; Hotta et al. 2016). Because of its ubiquity, the origin of the magnetic field is of great interest in various research fields.

The Weibel instability arises from anisotropy in the plasma velocity distribution function (Fried 1959; Weibel 1959). Because the instability feeds the free energy of the anisotropy and converts it into magnetic energy, it is one of the promising mechanisms responsible for the origin of the magnetic field. Following early theoretical works (Medvedev & Loeb 1999; Frederiksen et al. 2004; Hededal et al. 2004; Kato 2005; Medvedev et al. 2006; Silva 2006; Achterberg et al. 2007; Achterberg & Wiersma 2007; Bret 2009), numerical simulations have revealed its wide applicability. The generated magnetic field plays crucial roles in collisionless shock formation (Kato 2007; Kato & Takabe 2008; Spitkovsky 2008), associated particle accelerations (Sironi et al. 2013; Matsumoto et al. 2015), and afterglows of gamma-ray bursts (Nousek et al. 2006; Tomita & Ohira 2016). Recent experiments using high-power laser facilities have succeeded in showing detailed images of the nonlinear structure characterized by filaments of the current (Fox et al. 2013; Huntington et al. 2015; Park et al. 2015; Huntington et al. 2017).

Kato (2005) proposed a theoretical model describing magnetic field evolution in cylindrical currents of the electron Weibel instability, and found that the maximum magnetic field could be obtained when the current reached the so-called Alfvén current (Alfvén 1939; Lawson 1958; Honda 2000),

$$\begin{eqnarray}&&{I}_{{\rm{A}},{\rm{s}}}={I}_{0}\langle \gamma {\beta }_{| | }\rangle ,\end{eqnarray} \tag{ 1 }$$

where ${I}_{0}\equiv {m}_{s}{c}^{3}/{q}_{s}$, m_s, and q_s are the mass and charge of the particle of each species s, respectively, and c is speed of light. ${\beta }_{| | }$ is the magnitude of the three-velocity in the unit of light speed along the beam direction, and γ is the Lorentz factor. The angle brackets $\langle \rangle$ denote an average over the beam volume. This gives an upper limit of the cylindrical current in which the Larmor radius of particles in the self-generated magnetic field is comparable with the filament size, and thus determines the magnetic field saturation.

The saturation mechanism of the ion Weibel instability is not well understood despite its importance in understanding attainable maximum magnetic energy. Using particle-in-cell (PIC) simulations, Frederiksen et al. (2004) reported on the early evolution of the ion Weibel filaments and found that the generated magnetic field keeps growing even after saturation of the electron Weibel instability. In contrast, Kumar et al. (2015) reported that the Weibel filaments finally disappear in the two-dimensional (2D) case because of strong electric fields generated by either the Buneman instability or by an oblique unstable mode (Bret 2009). More recently, an analytical model of coalescence of the ion Weibel filaments was proposed for non-relativistic plasmas based on the kinetic quasi-linear theory of the instability, and was validated by 2D and three-dimensional (3D) PIC simulations in the weak nonlinear regime (Ruyer et al. 2015, 2016).

In this Letter, we present the evolution and saturation of the Weibel instability in interpenetrating relativistic ion-electron plasmas by means of 3D PIC simulations. Large-scale, long-term simulations enabled us to elucidate for the first time the coalescence of ion-scale current filaments and the resulting magnetic field saturation after reaching the ion Alfvén current limit.

To explore the nonlinear evolution of the Weibel instability, we used a fully kinetic electromagnetic PIC code (Matsumoto et al. 2017), which implements the quadratic function for the particle weighting and the implicit finite-difference time-domain (FDTD) scheme for the Maxwell equation. Simulation runs were initialized by two cold unmagnetized counter flows with a bulk Lorentz factor of γ = 5 in the laboratory frame. The beams were set in the x-direction with the periodic boundary condition in all directions. The cold counter-streaming condition induced a very large anisotropy of the velocity distribution in the system. Consequently, the relativistic Weibel instability grew rapidly in the whole region of the simulation domain. In the following, time and spatial length were normalized to the inverse of the non-relativistic ion plasma frequency ${\omega }_{p,i}^{-1}$ and the ion inertial length c/ω_p,i using the initial values, respectively. Ten particles per cell were used for each species of ion and electron with the mass ratio⁴ of M/m = 25. The cell size was set as Δ = 0.1 c/ω_p,e, equivalent to $0.045\ c/{\tilde{\omega }}_{p,e}$ in terms of the proper plasma frequency ${\tilde{\omega }}_{p,e}$, to resolve the initial electron Weibel instability accurately. The numerical Cherenkov instability in the relativistic plasma flows was suppressed by carefully choosing a magical time step size of ${\rm{\Delta }}t={\rm{\Delta }}/c(=0.1{\omega }_{{pe}}^{-1}$) specific to our semi-implicit PIC code (Ikeya & Matsumoto 2015). The simulation domain size was taken as L_x × L_y × L_z = 20.0 c/ω_p,i × 51.2 c/ω_p,i × 51.2 c/ω_p,i. Such very large-scale PIC simulations were performed on the Japanese K computer using 16,384 nodes (131,072 processor cores) and 200 TB of physical memory.

Figure 1 shows 3D profiles of the magnetization parameter _B ≡ B²/8π Mn₀c² (γ − 1), where B is the magnetic field strength, representing the ratio of the generated magnetic energy to the initial beam kinetic energy. At a characteristic time stage of the electron Alfvén current ($t=20\ {\omega }_{p,i}^{-1}$), the simulation domain is filled with filamentary structures elongated in the background beam direction with a typical size of 10 times the electron inertia length (Figure 1(a)). These filaments continue to coalesce into larger scales (inverse cascade) beyond a time stage of $t=280\ {\omega }_{p,i}^{-1}$ at which the ion Alfvén current is reached (Figure 1(b)).

Temporal evolution of the mean current in a filament is shown in Figure 2(a), and is compared with the Alfvén current (Equation (1)) for each species of the ion (s = i) and electron (s = e) using the initial beam speed. Here, the current is $I=\langle | J| {\rangle }_{{yz}}\pi {R}^{2}$, where J is the current density and R was measured as a half wavelength at the peak of the Fourier-transformed current density in the y–z plane, and $\langle {\rangle }_{{yz}}$ denotes the area average in the y–z plane. The current reached the electron Alfvén current at $t=20\ {\omega }_{p,i}^{-1}$, after which it continuously grew and finally reached the ion Alfvén current at $t=280\ {\omega }_{p,i}^{-1}$. The current did not exceed the Alfvén current limit up to $560\ {\omega }_{p,i}^{-1}$ in the present simulation run.

Zoom In Zoom Out Reset image size

The magnetic field energy (_B averaged in the simulation domain) also saturated around the time when the current reached the ion Alfvén current and sustained at _B ∼ 0.015 (Figure 2(b)). This strong magnetic field level was maintained until the simulation run ended at $t=560\ {\omega }_{p,i}^{-1}$. Comparisons with 2D simulations exhibited distinct differences in the saturation phase even using the same parameters. In the case of the 2D in-plane current (simulations in the x–z plane), _B evolved similarly to the 3D case until $t=100\ {\omega }_{p,i}^{-1}$. After that, however, _B decreased rapidly. In the case of the 2D out-of-plane current (simulations in the y–z plane), _B saturated at $t=10\ {\omega }_{p,i}^{-1}$, and then _B decreased monotonically.

Regarding particle acceleration, we did not observe non-thermal particles clearly during the saturation process. The particle velocity distribution functions for protons and electrons are both essentially thermal.

Figures 3(a) and (b) present 2D profiles of the current density for 3D and 2D in-plane cases. It clearly shows that the current filaments in the 2D in-plane case are much more fragmented than in the 3D case, whereas in the 3D case, fluctuations between filaments can be found. The difference stems from the geometrical restriction in the 2D in-plane case, in which a current filament first collides with a filament with the different polarity in the course of coalescence with the nearby filament with the same polarity. This results in enhanced electrostatic field excitation by the beam instability (Kumar et al. 2015). Indeed, the breakup of the filaments in the 2D in-plane case occurred right after the peak of the electrostatic field ($t=10\ {\omega }_{p,i}^{-1}$) as shown in Figure 3(c). In contrast, oppositely signed filaments in the 3D case can pass each other without collisions during the coalescence of the filaments with the same polarity, resulting in a weak excitation of the electrostatic mode (blue curve in Figure 3(c)). These results suggest that the breakup of current filaments, and thus the rapid decay in magnetic energy, are characteristic only in the 2D in-plane case (Haugbølle 2011).

Zoom In Zoom Out Reset image size

Geometrically, it could be expected that the current filaments in the 2D out-of-plane case would evolve similarly to the 3D case because the current filaments can avoid collision with the neighboring oppositely flowing ones. Figure 2(b), however, shows that the magnetic energy monotonically decreased after saturation at $t=10\ {\omega }_{p,i}^{-1}$, and finally became much smaller than in the 3D case. Figure 4 shows snapshots of _B for the 3D (Figure 4(a)) and the 2D out-of-plane (Figure 4(b)) cases in the final evolution stage at $t=560\ {\omega }_{p,i}^{-1}$. In the late nonlinear stage, the currents were no longer isolated cylindrical beams (Kato 2005), but were rather connected to each other to form a network structure. It was found that strong magnetic field areas in the 3D case were distributed broadly at the interface of the currents, whereas in the 2D out-of-plane case they were confined to a very narrow layer. The local maximum values, however, were similar to each other (_B ∼ 0.1).

Zoom In Zoom Out Reset image size

The broad strong magnetic field areas in the 3D case can be attributed to an electron gyro motion around the current interface in the nonlinear stage. Figure 4(c) shows the temporal evolution of the electron temperature perpendicular to the initial beam direction (x-direction) for the 3D (blue curve) and the 2D out-of-plane (green curve) cases. We found in the 3D result that electrons were heated preferentially in the perpendicular direction during coalescence of ion-scale filaments, and were finally energized to the relativistic temperature (k_B T_⊥ ∼ 10 mc²) at $t=560\ {\omega }_{p,i}^{-1}$, whereas the 2D result showed a heating process only in the very early stage, followed by a constant non-relativistic temperature (k_BT_⊥ ∼ 0.5 mc²).

In the 3D case, the non-uniform structures in the x-direction (Figure 3(a)) can cause the mixing of current filaments in the x-direction, resulting in significant heating of the electrons. This allows the broadening of the strong magnetic field area, and hence the total _B is 10 times larger than the 2D out-of-plane case. Figure 4(d) is the temporal evolution of the area where _B > 5 × 10⁻³, that is, the strongly magnetized region. In the 2D out-of-plane case, the strongly magnetized area quickly decreased after its peak in the early stage at $t=10\ {\omega }_{p,i}^{-1}$, whereas in the 3D case, it was sustained even in the late nonlinear stage. Note that the evolution of the strongly magnetized region in Figure 4(d) coincides with Figure 2(b), indicating that the different evolution of _B between the 3D and the 2D out-of-plane cases was the result of different filling factors of strongly magnetized regions in the simulation domain rather than the magnetic field strength itself.

We have presented the saturation mechanisms of the ion Weibel instability by means of long-term, large-scale PIC simulations. It was found that the Weibel-generated magnetic fields sustained for long time periods after reaching the ion Alfvén current limit in the 3D space; the Weibel filaments are stable during the filament merging process against secondary instabilities such as the Buneman and other electrostatic modes found in 2D in-plane simulation studies. We found that electron heating processes continue during the coalescence of ion-scale filaments in the late nonlinear stage and are crucial for the sustained magnetic field. Concerning the afterglows of gamma-ray bursts, the obtained value, _B ∼ 0.015, supports a model with moderately strong magnetic field (Medvedev & Loeb 1999). More importantly, our study showed that the magnetic field generated by the Weibel instability is stably sustained much longer than the plasma timescale, that is, at least for 560 ion plasma oscillations. Note that some recent studies of the timescale of the afterglow emission in the downstream reported the estimation (Sironi & Goodman 2007; Tomita & Ohira 2016) $7.4\times {10}^{7}\ {\omega }_{{\rm{p}},{\rm{i}}}^{-1}$ ${E}_{\mathrm{iso},53}^{1/3}{n}_{0}^{1/6}{{\rm{\Gamma }}}_{100}^{-5/3}$, where E_iso,53 is the isotropic equivalent energy normalized to 10⁵³ erg, n₀ is the number density of interstellar medium normalized to $1\,{\mathrm{cm}}^{-3}$, and Γ₁₀₀ is the Lorentz factor of the blastwave normalized to 100. Although that simulation timescale is still much shorter than this one, our results imply a possibility of a much more stable Weibel magnetic field than previously thought.

The present results indicate that collisionless shock simulations in 2D in-plane configurations underestimated the shock forming time, and overestimated the particle acceleration efficiency by the electrostatic field in the shock transition region. This is because the Buneman-type instability observed in the 2D in-plane case was suppressed in our 3D results, and the Weibel instability becomes dominant in relativistic and weakly magnetized counter flows (Bret 2009). Stochastic particle accelerations with Weibel-generated magnetic turbulence (Sironi et al. 2015) are more important than the direct acceleration by the electric field. As shown in the present work, the magnetic turbulence generated by the Weibel instability can be stronger in 3D than in 2D. Thus, this acceleration process can be more efficient in 3D than those observed in 2D simulations so far (Matsumoto et al. 2017).

Finally, we consider that the present results can provide a guideline for designing laboratory experiments using high-power laser facilities aiming to better understand magnetic field generation by the Weibel instability and collisionless shock formation. The experiment timescale of $\sim 200\ {\omega }_{p,i}^{-1}$ and the spatial scale size of $20\ c/{\omega }_{p,i}\times 50\ c/{\omega }_{p,i}\times 50\ c/{\omega }_{p,i}$ were required to allow current filaments to evolve into the ion Alfvén current and true magnetic field saturation.

We would like to thank T. Amano, M. Hoshino, and Y. Ohira for many fruitful comments and discussions. This research used computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (Project ID:hp160121, hp170125, hp170231) and Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan. This work was supported in part by the Postdoctoral Fellowships by the Japan Society for the Promotion of Science No. 201506571 (M.T.), and JSPS KAKENHI grant No, 17H02877 (Y.M.).