Evolution of three-dimensional Relativistic Ion Weibel Instability: Competition with Kink Instability

TL;DR

This paper investigates the nonlinear behavior of the relativistic Weibel instability in collisionless plasmas through simulations, revealing its stability against kink instability and its potential to generate seed magnetic fields for MHD dynamos.

Contribution

It provides a detailed analysis of the nonlinear saturation of the relativistic Weibel instability and demonstrates its stability against kink instability in the MHD regime, supporting its role in magnetic field generation.

Findings

Weibel instability's nonlinear phase characterized by filament current densities.

Relativistic kink instability does not affect Weibel filaments up to 1000 ω_{p,i}^{-1}.

Stable magnetic fields from Weibel instability can seed MHD dynamos.

Abstract

In this paper, we report our recent findings on the relativistic Weibel instability and its nonlinear saturation by performing numerical simulations of collisionless plasmas. Analysis of the obtained numerical results revealed that the nonlinear phase of the Weibel instability can be described by characteristic phases based on the Weibel filaments' current density in terms of particle and Alfv\'en limit currents. We also analyzed the relativistic kink instability based on the energy principle in the magnetohydrodynamic (MHD) regime, and found that the Weibel filaments do not suffer from the kink-type instability in the MHD regime up to 1000 $\omega_{p,i}^{-1}$. This finding allowed a magnetic field to be sustained by relativistic Weibel instability that was stable enough to be a seed for MHD dynamos.

Tables

Table 1. Table 1: List of the parameters. γ 𝛾 \gamma is the initial Lorentz factor of the particles, and M , m 𝑀 𝑚 M,m are ion and electron masses, respectively; n 0 subscript 𝑛 0 n_{0} is the initial particle number per cell; L x , L ⟂ subscript 𝐿 x subscript 𝐿 perpendicular-to L_{\rm x},L_{\perp} are the numerical box size in parallel and perpendicular to the initial beam direction.

Name	$\gamma$	$M/m$	$n_0$	$L_x{\omega }_{p,i}/c$	$L_{⟂}{\omega }_{p,i}/c$
runA	2	100	20	18	15.36
runAa	2	100	20	18	6
runAb	2	100	20	18	1
runB	5	100	20	18	15.36
runB1	5	4	20	18	15.36
runB2	5	25	20	18	15.36
runB3 (2D,x-z)	5	1836	20	49	95.6
runBa (TMK18)	5	25	10	20	51.2
runC	20	100	20	18	15.36

Equations

$$I_{\mathrm{A}}=I_0⟨\gamma {\beta }_{||}⟩,$$

$$I_{\mathrm{P}}=\pi R^2{J}_{\mathrm{P}},J_{\mathrm{P}}\equiv e\mu n_s⟨{\beta }_{||}⟩c,$$

$${ϵ}_{\mathrm{B}}\equiv \frac{B^2}{8\pi Mn{c}^2(\gamma -1)},$$

$$R\propto \{\begin{array}{ll}{t}^2& (\text{if}M/m\lesssim 25),\\ t^{0.76}& (\text{if}M/m=100).\end{array}$$

$$\mathcal{E}=\int d{x}^3{T}_{00}=\int d{x}^3\left[\rho h{\gamma }^2-p+\frac{B^2+E^2}{8\pi }\right].$$

$$\mathcal{E}=\int d{x}^3\left[\rho \gamma +\rho \gamma (\gamma -1)+\frac{({\gamma }^2-1)\Gamma +1}{\Gamma -1}p+\frac{B^2+E^2}{8\pi }\right].$$

$$\mathcal{E}=\int d{x}^3\left[\rho h{\dot{\xi }}^2+\rho +\frac{p}{\Gamma -1}+\frac{B^2+E^2}{8\pi }\right]+O({\dot{\xi }}^3).$$

$$4\pi J=\frac{1}{r}\frac{d}{dr}(rB)c.$$

$$B(r<R)=\frac{\sqrt{2}{I}_0}{c}\frac{r}{R^2},$$

$$B(r>R)=\frac{\sqrt{2}{I}_0}{cr},$$

$$\nabla p_g=\mathbf{J}\times \mathbf{B}.$$

$$\delta W^{(2)}=\frac{1}{2}\int d{x}^3\left[\frac{Q^2}{4\pi }+\mathbf{J}\cdot (\xi \times \mathbf{Q})+\Gamma p_{\mathrm{g}}(\nabla \cdot \xi {)}^2+(\xi \cdot \nabla p_{\mathrm{g}})(\nabla \cdot \xi )\right],$$

$$\xi =\mathrm{Re}\left[({\xi }_r,i{\xi }_{\theta },i{\xi }_z)e^{i(m\theta +kz)}\right],$$

$${\xi }_r$$

$${\xi }_{\theta }$$

$${\xi }_z$$

$$\delta W^{(2)}=\frac{\pi }{2}{\int }_0^{\infty }g{\xi }^{2}dr,$$

$$g=\frac{1}{r}\frac{m^2{B}^2}{k^2{r}^2+m^2}+\frac{m^2{B}^2}{r}-2\frac{B}{r}\frac{d}{dr}(rB)+m^2\frac{d}{dr}\left[\frac{B^2}{k^2{r}^2+m^2}\right].$$

$$g=-\frac{k^{2}r}{k^2{r}^2+1}{B}_{\theta }^2.$$

$$\delta W^{(2)}$$

$$F(k,R,R_{\mathrm{c}})$$

$$F(k,R,R_{\mathrm{c}})\to \{\begin{array}{ll}\left[1/2+2\ln (R_{\mathrm{c}}/R)\right]k^2& (\text{if}kR\ll 1)\\ 2-R^2/R_{\mathrm{c}}^2& (\text{if}kR\gg 1)\end{array},$$

$${\Gamma }_{\mathrm{grow}}\lesssim \sqrt{\frac{\delta W^{(2)}}{K}},$$

$$K\equiv \int \rho h({\xi }_r^2+{\xi }_{\theta }^2+{\xi }_{\varphi }^2).$$

$$\rho h\frac{{\partial }^2}{\partial t^2}\xi =-\nabla p+\mathbf{j}\times \mathbf{B},$$

$$\frac{\pi {\xi }^2}{2k^2}\left[{\rho }_0{h}_{\mathrm{out}}\{5k^2{R}_{\mathrm{c}}^2-3\ln (k^2{R}_{\mathrm{c}}^2+1)\}-4\Delta p\{5k^2{R}^2-3\ln (k^2{R}^2+1)\}\right],$$

$${\Gamma }_{\mathrm{grow}}={\left[\frac{I_0^2{k}^2}{R^2{c}^2}\frac{1-\ln [1+k^2{R}^2]/k^2{R}^2]+k^2{R}^2\left\{\ln [1+1/k^2{R}^2]-\ln [1+1/k^2{R}_{\mathrm{c}}^2]\right\}}{{\rho }_0{h}_{\mathrm{out}}\{5k^2{R}_{\mathrm{c}}^2-3\ln (k^2{R}_{\mathrm{c}}^2+1)\}-4\Delta p\{5k^2{R}^2-3\ln (k^2{R}^2+1)\}}\right]}^{1/2}$$

$${\Gamma }_{\mathrm{grow}}\to \{\begin{array}{ll}\frac{c_{\mathrm{A}}k}{2\sqrt{2}}\left(\frac{R}{R_{\mathrm{c}}}\right)\sqrt{\frac{1+4\ln (R/R_{\mathrm{c}})}{1+\frac{4\Delta p}{{\rho }_0{h}_0}\frac{R_{\mathrm{c}}^2-R^2}{R_{\mathrm{c}}^2}}}& (\text{if}kR\ll 1),\\ \frac{c_{\mathrm{A}}}{\sqrt{10}{R}_{\mathrm{c}}}\sqrt{\frac{2-(R/R_{\mathrm{c}}{)}^2}{1+\frac{4\Delta p}{{\rho }_0{h}_0}\frac{R_{\mathrm{c}}^2-R^2}{R_{\mathrm{c}}^2}}}& (\text{if}kR\gg 1).\end{array}$$

$${\Gamma }_{\mathrm{grow}}\simeq 0.15\frac{c_{\mathrm{A}}}{R},$$

$${\Gamma }_{\mathrm{grow},\mathrm{lab}}\simeq 0.15\frac{c_{\mathrm{A}}}{{\gamma }_{\mathrm{fluid}}R},$$

$$N_{\mathrm{TPL}}\equiv \frac{L_{\mathrm{BD}}^2}{4\pi S_{+/-}},$$

Full text

editorials, notices — miscellaneous — catalogs — surveys

††journal: ApJ

1 Introduction

Weibel instability converts the anisotropy of plasma temperature into a magnetic field (Weibel, 1959; Fried, 1959), and is considered as an origin of magnetic field in various phenomena. In particular, it is expected that magnetic field turbulence resulting from relativistic Weibel instability is responsible for the relativistic collisionless shock formation (Kato & Takabe, 2008; Spitkovsky, 2008), which is one of the most probable candidates for accelerating the ultra-high-energy cosmic-rays (Hillas, 1984) through diffusive shock acceleration (Kirk & Duffy, 1999; Achterberg et al., 2001) and wake field acceleration (Tajima & Dawson, 1979; Hoshino, 2008; Iwamoto et al., 2017, 2018). And Weibel instability is also regarded as providing the necessary magnetic field to explain the observed afterglow of the gamma-ray burst (Nousek et al., 2006).

For the reasons mentioned above, the relativistic Weibel instability has been studied extensively under various situations and parameters (Kato, 2005; Califano et al., 2006; Achterberg & Wiersma, 2007; Achterberg et al., 2007; Pétri & Kirk, 2007). Importantly, recent developments in high-power laser facilities allow direct investigations of the Weibel instability in laboratory. In particular, the Omega and NIF laser experiments have recently performed several projects which are considered to have observed the early phase development of nonrelativistic Weibel filaments (Huntington et al., 2015; Park et al., 2015; Huntington et al., 2017). Although these recent laser experiments have not yet been successful in realizing collisionless shock formation, we expect that the next-generation laser facilities would succeed in achieving fully developed relativistic collisionless shock via the relativistic Weibel instability.

In this paper, we report the further analysis of our previous work on the relativistic Weibel instability (Takamoto et al., 2018) (TMK18 in the following). The Weibel instability is a very popular mechanism in plasma instability and has been studied extensively, including analytical and numerical approaches. However, these work were limited to either only the early phase or to the two-dimensional (2D) system. In our work, large-scale and long-time simulations of the relativistic Weibel instability were performed. In Section 2, the numerical setup is explained. In Section 3, numerical results and our analyses are presented. In Section 4, the kink instability in the Weibel filaments is discussed. Section 5 discusses the saturation mechanism of the Weibel instability, and Section 6 summarizes our results.

2 Numerical Setups

In this paper, the temporal evolution of the relativistic ion-electron Weibel instability is modeled by a collisionless plasma; it is integrated by a three-dimensional (3D) relativistic electromagnetic particle-in-cell (PIC) code developed by Matsumoto et al. (2017), which is a vectorized, hybrid-parallel 3D PIC simulation code with a quadratic particle weighting function. As considered in TMK18, two cold, $T=0$, unmagnetized beam counter-flows are modeled in our 2D and 3D simulation boxes. The beam is set in the x-direction; in each cell, half of the particles are set to move in the positive x-direction, and the other particles move in the negative x-direction. The periodic boundary is considered in all directions. Consequently, a relativistic ion-electron Weibel instability grows in the whole numerical box region because of the very large anisotropy. For a unit of length and time, the nonrelativistic ion inertial length $c/\omega_{p,i}$ and the inverse of the ion plasma frequency $\omega_{p,i}$ are considered. To resolve the initial electron Weibel instability accurately, the cell size is set as $\Delta=0.1c/\omega_{p,e}$, equivalent to $0.07c/\tilde{\omega}_{p,e}$ in the case of $\gamma=2$, $0.045c/\tilde{\omega}_{p,e}$ in the case of $\gamma=5$, and $0.022c/\tilde{\omega}_{p,e}$ in the case of $\gamma=20$, where $\tilde{\omega}_{p,e}$ is the relativistic electron plasma frequency. To reduce the numerical Cherenkov instability, the CFL number is set to unity, $\Delta t=\Delta/c$, which is a magical time-step size specific to our semi-implicit PIC code (Ikeya & Matsumoto, 2015). Simulation parameters of the initial beam bulk Lorentz factor $\gamma$, the ion-to-electron mass ratio $M/m$, number of particles per cell per species $n_{0}$, and the system size ($L_{\rm x}$, $L_{\rm y}=L_{\rm z}=L_{\perp}$) are summarized in table 1 for different simulation runs.

3 Results

3.1 Theoretical Consideration of Current Evolution

In this section, we briefly review the current evolution of Weibel filaments reported by Kato (2005) in the case of a 2D (out-of-plane current) study.

In the early nonlinear phase, the Weibel instability is known to generate current filaments because of the generated magnetic field, that is, two lines of current flowing in the same direction attract each other, but those in the opposite direction repel each other, resulting in an accumulation of currents.

In the later nonlinear phase, Kato (2005) pointed out that there are two possible upper limits of current; One is the Alfvén limit current (Alfvén, 1939), which results in the magnetic field for which the gyro-radius of the particles is equivalent to the current filaments’ radius. In Gaussian units, this is given as:

[TABLE]

where $I_{0}\equiv m_{s}c^{3}/q_{s}$, $m_{s}$ and $q_{s}$ are the mass and charge of the particles, respectively, and $s$ means ion or electrons. $\beta_{||}$ is the magnitude of the three-velocity in units of light velocity along the filaments’ direction, and $\gamma$ is the Lorentz factor. The angle brackets $\langle\rangle$ denote the average over the beam volume. This gives a theoretically maximal magnetic field (Alfvén, 1939; Medvedev & Loeb, 1999). The other limit is the particle limit current, which is the maximum current $I_{\rm P}$ carried by all the particles in the beam, defined as:

[TABLE]

where $R$ is the radius of filaments, $n_{s}$ is the density of particles, and $\mu=1$ for electron-proton plasmas, and $\mu=2$ for electron-positron plasma.

Following these upper current limits, Kato (2005) pointed out that there are two paths of evolution depending on the initial anisotropy of particle momentum: i.e., the Alfvén limit ($I_{\rm A}<I_{\rm P}$), and the particle limit ($I_{\rm A}>I_{\rm P}$). The Alfvén limit occurs when the initial anisotropy is small or the beams have a nonrelativistic velocity. In this case, the currents by Weibel instability saturates at the Alfvén current, and increase its radius without changing its current, resulting in a gradually decreasing magnetic field. In contrast, the particle limit occurs when the initial anisotropy is large or for beams with relativistic velocity. In this case, currents continue to evolve by coalescing until reaching the Alfvén current.

Kato (2005) studied these saturation mechanisms using a 2D relativistic PIC simulation, and reported a detailed analysis. However, it is still unclear whether these results can be applied in the 3D case because there are several important instabilities that occur only when a 3D case is considered, such as the kink instability. In the following, we investigate the saturation mechanism of the Weibel instability in a 3D space.

3.2 Temporal Evolution of Magnetic Field and Current

Figure 1 shows the temporal evolution of the magnetic field generated by Weibel instability, which is measured by

[TABLE]

where $M$ is the ion mass, $n$ is the initial ion number density in laboratory frame, $c$ is the light velocity, and $\gamma$ is the Lorentz factor of the initial ion. The term $\epsilon_{\rm B}$ measures the energy density of the magnetic field with respect to the initial ion kinetic energy density, and the plotted values are averaged over the entire volume. The left panel shows the Lorentz factor dependence of the generated magnetic field in terms of time scale in the particle comoving frame, plotting the runs A, B, C ($\gamma=2,5,20$). It shows that the generated magnetic is insensitive to the plasma Lorentz factor. This indicates that the generated magnetic field can be described by Ampère’s law: $B\sim 4\pi nq(c/\tilde{\omega}_{\rm p,i})$ which makes the magnetization parameter $\epsilon_{\rm B}$ a constant value 111 We consider that the differences in the linear phase reflect those in the spatial resolution because of the Lorentz contraction. .

The right panel shows the mass dependence of the magnetic field evolution. Although they show a very different evolution before $10\omega_{\rm p,i}^{-1}$, the saturation values after $10\omega_{\rm p,i}^{-1}$ take a similar value around 0.01 to 0.02. This can be understood from the current evolution discussed in the following.

Figure 2 is the temporal evolution of current density in the x-direction in the case of runBa, where the plotted values are measured at $x=L_{\rm x}/2$ averaged over this plane 222 In Figure 2, both the particle and the Alfvén limit currents were calculated using Equations (1) and (2). The physical variables were measured at a plane on $x=L_{\rm x}/2$ and were averaged over that plane. Concerning the $\gamma$ factor in the Alfvén limit current, the initial Lorentz factor was used.

.

This indicates that the temporal evolution can be divided into six phases, i.e.: (1) a linear electron Weibel growth phase, (2) nonlinear phase of electron Weibel instability following the electron particle limit current, (3) a linear ion Weibel growth phase, (4) a nonlinear ion Weibel growth phase following the ion particle limit current, (5) a coalescing phase of Weibel filaments, (6) a saturation phase due to either the Alfvén limit current (runBa) or limited numerical box size. Phase (3) can be observed more clearly in run B3, where $\epsilon_{\rm B}$ initially increases exponentially because of the growth in the electron Weibel instability (1) and (2) , and changes its growth rate around $t=0.3\omega_{p,i}^{-1}$. In the other runs, although the phases (1) and (2) exist, the smaller value of mass ratio causes these regions to degenerate into the linear growth phase of ion Weibel instability. Around $t=4\omega_{p,i}^{-1}$, the current density becomes larger than the electron particle limit, and instead starts to follow the ion particle limit around $t=10\omega_{p,i}^{-1}$, which coincides with the first peak of $\epsilon_{\rm B}$ other than runB3 whose large mass ratio allows the electron particle limit to show a clear peak. After that, the coalescence of the filaments starts, and the filaments increase their radius as discussed later. At this phase, the current density does not follow the ion particle limit current density and its value starts to decrease because the current reaches the electron Alfvén limit and electrons can no longer participate in increasing current. Finally, at $t\gtrsim 300\omega_{p,i}^{-1}$ of runBa, the current density also reaches the ion Alfvén current, which also limits the ion contribution to the current, resulting in a reduced current density. The $\epsilon_{\rm B}$ reaches the saturation value at around $0.02$. In our 3D runs, we did not observe the saturation because of the kink instability, which will be discussed in detail in Section 4.

3.3 Temporal Evolution of Radius

The left panel of Figure 3 is the temporal evolution of the ion filament radii with their initial Lorentz factor varying in the case of $M/m=100$ (runs A,B,C). In this paper, the filament radius is determined as half of the peak wavelength of the 2D energy spectrum of the ion velocity in a plane perpendicular to the initial beam direction. This indicates that the evolution of the radius is nearly independent of the initial Lorentz factor though the timestep of the initial rapid increase of their radius becomes slower as the Lorentz factor increase because of the longer growth time. The temporal evolution of the radius $R$ in the coalescence phase can be fitted by: $R\propto t^{0.76}$.

The right panel of Figure 3 is the temporal evolution of the ion filament radii with varying their mass ratio $M/m$ in the case of $\gamma=5$ (runs B,B1,B2,Ba). It shows that the temporal evolution of the filament radius can be in general described by a power law, $R\propto t^{n}$. The temporal evolution of the radius becomes shallower in the later phase around $t>100\omega_{p,i}$. This is because each filament covers more than 10 % of the numerical box in this phase and the effects of the periodic boundary cannot be neglected. In particular, it indicates that the radius curve becomes steeper when the mass ratio is small, $M/m\lesssim 25$, and shallower when the mass ratio is large, $M/m=100$, which can be written as:

[TABLE]

3.4 Temporal Evolution of Anisotropy

The temporal evolution of temperature of ions and electrons of runBa (TMK18) is presented in Figure 4. Initially, they only have the parallel temperature from counter-beam distribution, and a very small perpendicular temperature from numerical heating 333 The parallel and perpendicular temperatures are measured in terms of the initial beam direction (x-direction), but the not magnetic field.

. The figure shows that the parallel temperature of the ions and electrons starts to reduce when they reach the particle-limit current. This is because the particle-limit phase is equivalent to the filament coalescence phase if the initial anisotropy is sufficiently large. Interestingly, the temperature evolution of ions and electrons after the particle limit shows a different behavior. Electrons become nearly isotropic immediately after reaching the particle-limit current, and continue to be isotropic.

the ion perpendicular temperature also starts to grow during electron’s particle limit because of turbulent magnetic field generated from electron current filaments. However, the ions do not become isotropic even after reaching ion particle-limit current. We consider that this is because the magnetic at this stage field was too weak to bend the ions’ trajectory; In contrast, the electrons were magnetized because of its small mass. This allows their temperature to grow nearly isotropically even after ion filament coalescing phase, which can also be observed in the right panel of Figure 4 where electrons’ perpendicular temperature increasing with ion filament radius in the power law relation 444 Note that the electron temperature is not perfectly isotropic but perpendicular temperature is a little larger. We consider that this is due to the nature of the current filament coalescing processes that occur only with current filaments flowing in the same direction. This allows strong heating in the perpendicular direction to the current flow.

. The fluctuations in beam direction (x-direction) play an important role in the electrons’ perpendicular temperature evolution. In TMK18, it was reported that the electron’s perpendicular temperature stayed cold in the case of 2D simulation with an out-of-plane beam direction. The difference comes from the existence of a fluctuation in the beam direction, which scatters electrons through the electrostatic field resulting from these fluctuations and causes effective heating. This fact claims that we need greater caution when using the 2D approximation of plasma simulations, such as PIC, and applying the obtained result to the 3D case.

4 Relativistic Magnetohydrodynamic (MHD) Kink Instability

In this section, the kink instability in this simulation is considered. We derive the MHD kink instability in Weibel filaments using the energy principle (Newcomb, 1960; Milosavljević & Nakar, 2006). First, we discuss the relativistic energy principle following nonrelativistic work by Kulsrud (2005). Note that our treatment is not covariant but considers only the fluid comoving frame, similar to a work by Takamoto (2018). In this section, we consider an MHD plasma for simplicity. Possible kinetic effects will be discussed in Section 5.

First, the total energy in a relativistic MHD plasma can be written as:

[TABLE]

Assuming the ideal equation of state: $h=1+(p/\rho)[\Gamma/(\Gamma-1)]$, this can be rewritten as:

[TABLE]

We consider an energy change from a static state, so that the Lorentz factor $\gamma$ can be replaced by $1+\dot{\xi}^{2}/2$ where $\xi$ is a small displacement. The above equation reduces to:

[TABLE]

This shows that we can use the same discussion by Kulsrud (2005) formally by just modifying the kinetic energy from $\rho_{0}\dot{\xi}^{2}/2$ to $\rho h\dot{\xi}^{2}$ 555 This expression seems different from the nonrelativistic kinetic energy by a factor of 2. This is because the nonrelativistic density is measured in the laboratory frame, and is given as $\rho_{\rm NR}\simeq\rho(1+\dot{\xi}^{2}/2)$, as discussed by Landau & Lifshitz (1959). In addition, we should subtract the rest mass energy ($\rho\gamma$) from the energy density. This results in the nonrelativistic expression.

. Note that the additional enthalpy $h$ formally includes the effect of the Lorentz factor of the electrons and positively charged particles moving with a relativistic velocity. Hence, it is formally able to consider the direction of particle counterflows if we consider the anisotropic temperatures $T_{||},\ T_{\perp}$ that are responsible for the Weibel instability. In this section, however, we only consider an isotropic temperature for simplicity.

Following Milosavljević & Nakar (2006), we consider a cylindrical and symmetric current filament. Using the toroidal magnetic field ${\bf B}=B(r){\bf e}_{\theta}$ and the axial current density ${\bf J}=J(r){\bf e}_{z}$, the Ampère’s law can be written as:

[TABLE]

In the following, we assume a homogeneous current $I_{0}$ in the filament, $r<R$, where $R$ is the filament radius. From Equation (8), the root-mean-square value of the magnetic field is given as:

[TABLE]

and the magnetic field outside the current is

[TABLE]

where the reduction of $\sqrt{2}$ is due to the assumption of the presence of a strongly fluctuating magnetic field.

The static equilibrium is kept by the pressure balance:

[TABLE]

In the following, we consider a small perturbation on this static equilibrium state, and the growth of an MHD kink instability.

In general, the second-order potential energy change by a displacement $\xi$ is given as:

[TABLE]

where ${\bf Q}\equiv\nabla\times(\xi\times{\bf B})$. Assuming

[TABLE]

Newcomb (1960) found that the minimum $\delta W$ can be obtained if

[TABLE]

where $\xi$ is a function of $r$. The condition $\nabla\cdot\xi=0$ is satisfied. In this case, the $\delta W$ can be written as:

[TABLE]

where

[TABLE]

For the kink instability, we consider the case $m=1$. We also assume that the trial displacement vector $\xi$ is a constant, which allows us to reduce the complexity of the analysis. Although the obtained growth rate is not an exact value of the growth rate, obtaining an analytic expression of the growth rate provides a fruitful physical insight into this instability in the relativistic regime. The above procedure reduces Equation (18) to

[TABLE]

To proceed to the calculation of Equation (18), Equations (9) and (10) are substituted into the above equations. In this paper, we consider multiple filaments resulting from the Weibel instability and there are finite distances between them. Therefore, we set a new parameter $R_{\rm c}$ as half the distance between these filaments, and we use this as an upper value of integration in Equation (17). We obtain

[TABLE]

In the limit of a small and large value of $k$, the function $F(k,R,R_{\rm c})$ becomes

[TABLE]

and the detailed dependence on $k$ is given in Figure 5. Equations (20) and (21) show that $\delta W^{(2)}$ is always negative, which means that these filaments are unstable for the MHD kink instability for any wavelength.

Next, we calculate the kinetic energy. According to Kulsrud (2005), we can estimate the upper limit of the growth rate $\Gamma_{\rm grow}$ of perturbations as the following expression:

[TABLE]

where $K$ is the kinetic energy given as:

[TABLE]

Note that the kinetic energy includes the density $\rho$ and the pressure $p_{\rm gas}$ in its enthalpy for the following reason. The equation of motion of RMHD can be written as:

[TABLE]

in a second-order accuracy in $\xi$. This means that the discussion deriving the growth rate of the nonrelativistic case in Kulsrud (2005) can also be formally applied in the relativistic case if we replace the kinetic energy from $\rho\dot{\xi}^{2}/2$ to $\rho h\dot{\xi}^{2}$.

The top panel of Figure 6 shows the profiles of the ion density and pressure along a magenta line listed in the bottom panel of Figure 6. They show that the density is almost uniform in all regions, and pressure is constant with different values according to whether it is inside or outside of the filaments 666 The pressure in the bottom panel of Figure 6 is measured in the ion rest frame. Inside the filaments, the ions maintain the initial beam velocity due to the collisionless property and the temperature in the x-direction becomes cold. In contrast, in the contact region, there are two ion components with opposite flow directions, which results in a very high temperature in the x-direction. , which we set as $\rho_{0},p_{\rm in},p_{\rm out}$, respectively. Substituting Equations (14), (15), and (16), the kinetic energy becomes

[TABLE]

where $h_{\rm out}$ is the enthalpy in the contact region of the filaments, and $\Delta p\equiv p_{\rm out}-p_{\rm in}$. From Equations (20), (21), and (26), the growth rate is given as:

[TABLE]

In the long and short wavelength limit, this reduces to

[TABLE]

Figure 7 is a plot of Equation (27) when $R_{\rm c}=2R$ and $p_{\rm out}=2p_{\rm in}=2Mnc^{2}$. It shows that the growth rate reaches its maximum value around $k\simeq c/\omega_{p,i}$, the value of which can be well reproduced by the short wavelength limit in Equation (28).

From our simulation results, we set $R_{\rm c}\simeq 2R$, $p_{\rm out}\simeq 2p_{\rm in}\sim 2Mnc^{2}$, the growth rate in the range of wavelength shorter than $R$ becomes

[TABLE]

and this shows that the growth of the MHD kink instability in current filaments with Alfvén current is around 10 times the Alfvén crossing time of filaments 777 Our assumption of a comoving boundary with the current filament is in fact invalid and it has a steep velocity shear, which would cause another instability, such as the Kelvin-Helmholtz instability. However, this was neglected here to avoid a too-complex analysis. In addition, the radius $R$ is the shortest limit of validity of MHD approximation. Hence, Equation (29) also corresponds to the maximum growth rate which is valid in the case of the MHD approximation.

. Note that our expression differs from Equation (14) of Milosavljević & Nakar (2006) because of the modification in Equation (25).

Finally, this analysis was performed in the fluid rest frame; therefore, the growth rate in the laboratory frame could be obtained as:

[TABLE]

where $\gamma_{\rm fluid}$ is the Lorentz factor of the background fluid. This is also similar to the expression obtained in the case of force-free approximation (Lyubarskii, 1999) although its Lorentz factor is not exactly equal to the fluid one but to the phase velocity of the perturbation.

5 Saturation Mechanisms – Suppression of MHD Kink Instability

Our simulations indicate that the Alfvén current limit and isotropization of particles are responsible for the saturation of relativistic Weibel instability. However, it has long been believed that the Weibel filaments would suffer from kink instability, resulting in saturation of the filament evolution. In Section 4, we calculated the growth rate of kink instability in the case of Weibel filaments using the energy principle given by Equation (30). Left panel of Figure 8 plots the growing time of MHD kink instability, the inverse of Equation (30), in the cases of runBa. Note that the Alfvén velocity was calculated considering the effect of the relativistic pressure, which reduces the Alfvén velocity more than expected from $\epsilon_{\rm B}$. This indicated that the MHD kink instability cannot grow in both cases because the growth time is much longer than the simulation timescale.

In addition, we consider that the formation of a nontrivial 2D structure of current filaments in the plane perpendicular to the beam direction is also important. The bottom panel in Figure 6 shows that the Weibel filaments do not exist as isolated cylindrical filaments which was assumed for the analysis of kink instability, but as a network of multifilaments. This nontrivial network structure would cancel the nonlinear evolution of kink instability, and keep the filaments stable. To characterize the topology of the current filament network quantitatively, we define the following topological variable:

[TABLE]

where $L_{\rm BD}$ is the length of the boundary between the current filaments, and $S_{+/-}$ is the area of current filaments with a positive/negative current. The above new variable reduces to the number of filaments if the filaments are isolated cylindrical currents with a similar radius : $L_{\rm BD}\rightarrow\Sigma 2\pi R=2\pi NR$ and $S\rightarrow\Sigma\pi R^{2}=\pi NR^{2}$. This means that $N_{\rm TPL}$ decreases as filaments merge and a network of the filaments is formed. The right panel of Figure 8 is the temporal evolution of $N_{\rm TPL}$ and the number of pseudocircles measured through a typical length scale the value of which is shown in the right panel of Figure 3, which counts the number of filaments if the filaments are isolated 888Here, $S_{+/-}$ were measured by counting the area where the current density is larger than 2.5% of the averaged current density. $L_{\rm BD}^{2}$ were measured by counting the area with a current density of less than 2.5% of the averaged current density, and divided by $\pi R^{2}$, where R is the value plotted in Figure 3. The threshold 2.5% is comparable to the level of numerical fluctuations. . The figure shows that the $N_{\rm TPL}$ is nearly always less than the pseudocircles, indicating that the filaments form a network and the nonlinear evolution of kink instability cannot grow due to the cancellation in the network structure 999 Although the effective $N_{\rm TPL}$ becomes partly less than the actual value, we consider that this is due to the numerical error for calculating $N_{\rm TPL}$ from our numerical simulation results.

.

Although the above discussion assumes an MHD approximation for the kink instability, it is still possible that some kinetic effects play important roles because the radius of cylindrical Alfvén currents is comparable to the Larmor radius of the generated magnetic field. Ruyer & Fiuza (2018) indicated that the resonance between the particles’ gyromotion and the current filament motion can induce a kinetic kink-type instability. The unstable wavelength is predicted as comparable with one gyrorotation length: $2\pi v_{\perp,i}/\omega_{B}=2\pi\sqrt{\gamma/2\epsilon_{B}(\gamma-1)}(v_{\perp,i}/c)(c/\tilde{\omega}_{p,i})=2\pi\sqrt{\gamma^{2}/2\epsilon_{B}(\gamma-1)}(v_{\perp,i}/c)(c/\omega_{p,i})$ where $\omega_{B}=eB/\gamma Mc$ is the Larmor frequency and a cold plasma is assumed. The growing time is predicted as: $\sqrt{\gamma}c/(\kappa v_{\perp,i}\omega_{p,i})=c/(\kappa v_{\perp,i}\tilde{\omega}_{p,i})$, where $\kappa$ is the screening factor of current by hot electrons. The relativistic modification from the Lorentz dilation was applied in the above expressions. In our simulation (runBa), the growing time reduced to around $11\omega_{p,i}^{-1}$ and the unstable wavelength around $90c/\omega_{p,i}$, where we used $\gamma=5$, $\epsilon_{B}=0.015$, $\kappa=0.2$, and $v_{\perp,i}=c$. This means that it might be able to grow, but our numerical box is not large enough to check it, although it might be suppressed by a similar effect from the network structure discussed above. Unfortunately, the present numerical resource does not allow us to check the above mechanism with a kinetically sufficient resolution ($\Delta\lesssim 0.1c/\omega_{p,e}$), and we leave this as our future work 101010 This topic has now been actively studied, and the recent developments can be found in, e.g., Alves et al. (2018) and Nishikawa et al. (2019) which discuss the importance of the kinetic Kelvin-Helmholtz instability, the application to AGN jets, and so on.

.

6 Summary

In this paper, a detailed explanation of our previous work TMK18 was provided. It was found that the magnetization parameter $\epsilon_{\rm B}$ was nearly insensitive to the initial beam Lorentz factor and mass ratio, saturating around 0.01 to 0.03. It was also found that the evolution of the relativistic Weibel instability can be divided into six phases which are closely related with the beam plasma current limit: particle and Alfvén limit currents. The important finding was that the relativistic MHD kink instability was analyzed using the energy principle method. The obtained growth time is much longer than our simulation time. This indicates that the relativistic MHD kink instability may not suppress the relativistic Weibel instability, and the resulting magnetic field can be a seed of a large scale MHD instability.

We would like to thank Takanobu Amano, Masahiro Hoshino 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). This work is 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.).