MAGNETIC FIELD GENERATION IN CORE-SHEATH JETS VIA THE KINETIC KELVIN–HELMHOLTZ INSTABILITY

In this simulation study we use a core-sheath plasma jet structure instead of the counter-streaming plasma setups used in previous simulations by Alves et al. (2012, 2014), Grismayer et al. (2013a, 2013b), and Liang et al. (2013, 2014). The basic setup and illustrative results are shown in Figure 1. In our setup, a jet core with velocity v_core in the positive x direction resides in the middle of the computational box. The upper and lower quarters of the numerical grid contain a sheath plasma that can be stationary or moving with velocity v_sheath in the positive x direction (Nishikawa et al. 2013a, 2013b). This model is similar to the setup in our RMHD simulations (Mizuno et al. 2007) that used a cylindrical jet core. However, here we represent the jet core and sheath as plasma slabs. Initially, the system is charge and current neutral. The simulations have been performed using a numerical grid with (L_x, L_y, L_z) = (1005Δ, 205Δ, 205Δ) (simulation cell size: Δ = 1) and periodic boundary conditions in all dimensions. The jet and sheath (electron) plasma number density measured in the simulation frame is n_jt = n_am = 8. The electron skin depth, λ_s = c/ω_pe = 12.2Δ, where ω_pe = (e²n_am/₀m_e)^1/2 is the electron plasma frequency and the electron Debye length for the ambient electrons λ_D is 1.2Δ. The jet-electron thermal velocity is v_{jt, th, e} = 0.014c in the jet reference frame, where c is the speed of light. The electron thermal velocity in the ambient plasma is v_{am, th, e} = 0.03c, and ion thermal velocities are smaller by (m_i/m_e)^1/2. Simulations were performed using an electron–positron (e^±) plasma or an electron–proton (e⁻ − p⁺ with m_p/m_e = 1836) plasma for jet Lorentz factors of 1.5, 5.0, and 15.0 with the sheath plasma at rest (v_sheath = 0).

Zoom In Zoom Out Reset image size

An illustration of the development of the velocity shear surfaces is also shown in Figure 1 for e⁻ − p⁺ and e^± plasmas with v_core = 0.9978 c(γ_core = 15). For the e⁻ − p⁺ case, a nearly DC magnetic field is generated at the shear surfaces. The B_y magnetic field component is generated with negative values (blue) at z = 150Δ and positive values (red) at z = 50Δ. Additionally, a B_z (and B_x) magnetic field component, shown by the small arrows in Figures 1(b) and (c), is generated at the shear surfaces by current filaments (see Section 3). On the other hand, for the e^± case, a relatively long wavelength (∼100Δ) AC magnetic field is generated at the shear surfaces. Note the alternating B_y > 0 (red) and B_y < 0 (blue) in Figure 1(c) along the flow direction. These results are similar to those found by Liang et al. (2013a, 2013b). However, due to the 2D nature of their simulations and a counter-streaming setup, there are some differences in the structure.

We consider a sharp velocity shear surface at z = 0 with "jet" plasma at z > 0 and "ambient" plasma at z < 0 with flow in jet and/or ambient plasma in the x direction. Here the y direction is infinite. Following Gruzinov (2008), Alves et al. (2012, 2014), and Grismayer et al. (2013b), we assume uniform initial conditions on either side of the velocity shear surface, infinitely massive ions, and perturbations to the initial conditions of the following form:

$$\begin{eqnarray} \nonumber n(x,z,t) & = & n_0(z) + n_1(x,z,t) \\ \nonumber {\bf v}(x,z,t) & = & v_{\rm x0}(z) + {\bf v}_1(x,z,t)\\ {\bf E}(x,z,t) & = & E_{\rm x1}(x,z,t) + E_{\rm z1}(x,z,t)\\ \nonumber {\bf B}(x,z,t) & = & B_{\rm y1}(x,z,t)\\ \nonumber {\bf J}(x,z,t) & = & J_{\rm x1}(x,z,t) + J_{\rm z1}(x,z,t) . \end{eqnarray} \tag{ 1 }$$

We extend their results to a non-counter-streaming setup. Here we make the assumption that v_x0 > 0 is constant over the domain z > 0 but can take any constant positive or negative value v_x0≷0 over the domain z < 0. With these assumptions we look at density, velocity, current and electric field perturbations along the flow, x axis, that are also a function of the normal, z axis, to the velocity shear surface. The magnetic field perturbations are transverse to the flow, y axis, and parallel to the shear surface. It is assumed that perturbations are of the form

$$\begin{equation} f_1(x,z,t) = f_1(z) e^{i(kx-\omega t)}, \end{equation} \tag{ 2 }$$

and the wavevector k ≡ k_x is parallel to the flow direction. Thus we are considering a velocity shear surface that is infinite transverse to the flow direction and perturbations are independent of y, i.e., k_y = 0.

Derivation of the dispersion relation proceeds as in Alves et al. (2014) and the dispersion relation can be written in the following form¹³:

$$\begin{eqnarray} &&{-}[ k^2 + \omega ^2_{\rm p+}/c^2 - \omega ^2/c^2 ]^{1/2} \left[ {\omega _{\rm p+}^2/\gamma _{+}^2 \over (\omega - kv_{+})^2} - 1 \right]\nonumber\\ &&\quad\times \left[ \left({\omega _{\rm p+}^2 - \omega _{\rm p-}^2 \over c^2} \right) + \left({\omega ^2 \over c^2} - k^2 - {\omega _{\rm p+}^2 \over c^2}\right) \right] \nonumber\\ &&\quad+[ k^2 + \omega ^2_{\rm p-}/c^2 - \omega ^2/c^2 ]^{1/2} \left[ {\omega _{\rm p-}^2/\gamma _{-}^2 \over (\omega - kv_{-})^2} - 1 \right]\nonumber\\ &&\quad\times \left[ \left({\omega _{\rm p+}^2 - \omega _{\rm p-}^2 \over c^2} \right) - \left({\omega ^2 \over c^2} - k^2 - {\omega _{\rm p-}^2 \over c^2}\right) \right] = 0 , \end{eqnarray} \tag{ 3 }$$

with velocities v_±, associated Lorentz factors γ_±, and plasma frequencies $\omega _{\rm p\pm }^2 \equiv {4 \pi n_{\pm }e^2 / \gamma _{\pm } m_{\rm e}}$ appropriate to z₊ > 0 and z₋ < 0, respectively.

2.2.1. Analytic Solutions

Our generalization of previously published work to allow motion of the z < 0 plasma in the ±x direction, i.e., v₋≷0, allows comparison with existing velocity shear surface counter-streaming solutions and also allows for velocity shear surface solutions representing a high speed "jet" plasma moving through an already relativistic "ambient" plasma. In particular, our generalization provides velocity shear surface solutions appropriate to spine-sheath AGN jet scenarios (Mizuno et al. 2007; Hardee et al. 2007; Walg et al. 2013; Clausen-Brown et al. 2013; Murphy et al. 2013 and references therein) or the "needles-in-a-jet" or "jet-in-a-jet" scenarios proposed in the blazar AGN context (e.g., Nalewajko et al. 2011 and references therein). To avoid confusion we change the notation used in Equation (3) to n_jt = n₊, n_am = n₋, v_jt = v₊, v_am = v₋, γ_jt = γ₊ and γ_am = γ₋. We also use the definition

$$\begin{equation*} \omega _{\rm p}^2 \equiv {4 \pi n_{\rm e} e^2 \over \gamma ^3 m_{\rm e}} \end{equation*}$$

keeping the Lorentz factor cubed in the denominator as this represents the frequency for plasma oscillations parallel to the direction of motion. We make these changes and rewrite Equation (3) more compactly as

$$\begin{eqnarray} &&(k^2c^2 + \gamma _{\rm am}^2\omega _{\rm p,am}^2 - \omega ^2)^{1/2}(\omega -kv_{\rm am})^2 [(\omega - kv_{\rm jt})^2 - \omega _{\rm p,jt}^2]\nonumber \\ &&\quad+ (k^2c^2 + \gamma _{\rm jt}^2\omega _{\rm p,jt}^2 - \omega ^2)^{1/2}(\omega - kv_{\rm jt})^2\nonumber\\ &&\quad\times [(\omega - kv_{\rm am})^2 - \omega _{\rm p,am}^2] = 0. \end{eqnarray} \tag{ 4 }$$

For counter-streaming velocities v_am = −v_jt = −v₀ and equal densities n_jt = n_am = n₀, Equation (4) becomes (e.g., Gruzinov 2008)

$$\begin{equation} (k^2c^2 + \gamma _0^2\omega _{\rm p0}^2 - \omega ^2)^{1/2} \lbrace 2(\omega ^2 - k^2v_0^2)^2 - 2\omega _{\rm p0}^2(\omega ^2 + k^2v_0^2) \rbrace = 0, \end{equation} \tag{ 5 }$$

with a solution $\omega ^2 = k^2c^2 + \gamma _0^2\omega _{\rm p0}^2$ that can be identified with transverse electromagnetic waves (the electric E_z and magnetic B_y field components are transverse to the wavevector k = k_x) and solutions to

$$\begin{equation} \omega ^4 - (2k^2v_0^2 + \omega _{\rm p0}^2)\omega ^2 + (k^2v_0^2 - \omega _{\rm p0}^2)k^2v_0^2 = 0, \end{equation} \tag{ 6 }$$

given by

$$\begin{equation} \omega ^2 = {{\omega _{\rm p0}^2}\over {2}} \left[ \left(1 + 2 {{k^2v_0^2}\over {\omega _{\rm p0}^2}}\right) \pm \left(1 + 8 {{k^2v_0^2}\over {\omega _{\rm p0}^2}}\right)^{1/2} \right], \end{equation} \tag{ 7 }$$

which can be identified with longitudinal electrostatic plasma oscillations (the electric E_x field component is parallel to the wavevector k_x). The purely real solution, "+" sign in Equation (7), in the limit $k^2v_0^2/\omega _{\rm p0}^2 \ll 1$ is $\omega ^2 \sim \omega ^2_{\rm p0}$ and in the limit $k^2v_0^2/\omega _{\rm p0}^2 \gg 1$ is $\omega ^2 \sim k^2v^2_0$. The second solution, "−" sign in Equation (7), is purely imaginary when $k^2v_0^2/\omega _{\rm p0}^2 < 1$, has a maximum growth rate $\omega ^2 = - \omega _{\rm p0}^2/8$ when $k^2v_0^2/\omega _{\rm p0}^2 =3/8$, is purely real when $k^2v_0^2/\omega _{\rm p0}^2 > 1$, and in the limit $k^2v_0^2/\omega _{\rm p0}^2 \gg 1$ becomes $\omega ^2 \sim k^2v^2_0$. This second solution is identical to the classic electrostatic two-stream instability associated with interpenetrating counter-streaming equal density relativistic plasmas. Note the difference in the "transverse" plasma frequency $\gamma _0^2\omega _{\rm p0}^2 = 4 \pi n_0 e^2/ \gamma _0 m_{\rm e}$ associated with transverse waves and the "longitudinal" plasma frequency $\omega _{\rm p0}^2 = 4 \pi n_0 e^2/ \gamma _0^3 m_{\rm e}$ associated with longitudinal waves. If densities in jet and ambient plasmas are unequal, n_jt ≠ n_am, and we normalize by the "longitudinal" plasma frequency $\omega _{\rm p,jt} = 4 \pi n_{\rm jt} e^2/ \gamma _0^3 m_{\rm e}$ and define ω' ≡ ω/ω_{p, jt}, k' ≡ kv₀/ω_{p, jt} and β₀ ≡ v₀/c Equation (4) reduces to Equation (29) in Alves et al. (2014)

$$\begin{eqnarray} &&\left(\gamma _0^2 {n_{\rm am} \over n_{\rm jt}} + k^{\prime 2}/\beta _0^2- \omega ^{\prime 2}\right)^{1/2} [ (\omega ^{\prime } + k^{\prime })^2 - (\omega ^{\prime 2} - k^{\prime 2})^2 ] \nonumber \\ &&+(\gamma _0^2 + k^{\prime 2}/\beta _0^2 - \omega ^{\prime 2})^{1/2} \left[ {n_{\rm am} \over n_{\rm jt}}(\omega ^{\prime } - k^{\prime })^2 -(\omega ^{\prime 2} - k^{\prime 2})^2 \right] = 0,\nonumber\\ \end{eqnarray} \tag{ 8 }$$

albeit with Lorentz factors in the leading terms resulting from our "longitudinal" as opposed to the "transverse" plasma frequency normalization used in Alves et al. (2012, 2014). Note that the transverse electromagnetic wave solution is not allowed for unequal densities on opposite sides of the velocity shear surface.

Analytic solutions of the dispersion relation, Equation (4), allowing for different densities and velocities on either side of the velocity shear surface, can be found in the low (kc ≪ ω_p) and high (kc ≫ ω_p) wavenumber limits. In the low wavenumber limit, Equation (4) can be written as

$$\begin{eqnarray} &&(\gamma _{\rm am}^2\omega _{\rm p,am}^2 \,{-}\, \omega ^2)^{1/2} (\omega ^2 - \omega _{\rm p,jt}^2)(\omega - kv_{\rm am})^2 + (\gamma _{\rm jt}^2\omega _{\rm p,jt}^2 - \omega ^2)^{1/2}\nonumber\\ &&\quad \times (\omega ^2 - \omega _{\rm p,am}^2)(\omega - kv_{\rm jt})^2 \sim 0. \end{eqnarray} \tag{ 9 }$$

A complex solution to Equation (9) can be written as

$$\begin{eqnarray} \omega &\sim& {(\gamma _{\rm am}\omega _{\rm p,jt}kv_{\rm am} + \gamma _{\rm jt}\omega _{\rm p,am}kv_{\rm jt}) \over (\gamma _{\rm am}\omega _{\rm p,jt} + \gamma _{\rm jt}\omega _{\rm p,am})}\nonumber\\ && \pm i {(\gamma _{\rm am}\omega _{\rm p,jt} \gamma _{\rm jt}\omega _{\rm p,am})^{1/2} \over (\gamma _{\rm am}\omega _{\rm p,jt} + \gamma _{\rm jt}\omega _{\rm p,am})}k(v_{\rm jt}-v_{\rm am}). \end{eqnarray} \tag{ 10 }$$

In Equation (10) the real part gives the phase (drift) velocity and the imaginary part gives the temporal growth rate and directly shows the dependence of the growth rate on the velocity difference across the shear surface. Note that for counter-streaming v_am = −v_jt the phase (drift) velocity is zero provided densities are equal on either side of the velocity shear.

In the low wavenumber limit where v_am = 0 and γ_am = 1 relevant to the numerical simulations, Equation (10) becomes

$$\begin{equation} \omega \sim {(\gamma _{\rm jt}\omega _{\rm p,am}/\omega _{\rm p,jt}) \over (1 + \gamma _{\rm jt}\omega _{\rm p,am}/\omega _{\rm p,jt})}kv_{\rm jt} \pm i {(\gamma _{\rm jt}\omega _{\rm p,am}/\omega _{\rm p,jt})^{1/2} \over (1 + \gamma _{\rm jt}\omega _{\rm p,am}/\omega _{\rm p,jt})}kv_{\rm jt}. \end{equation} \tag{ 11 }$$

Here we see that the phase velocity (drift speed) v_ph ≡ ω/k → v_jt as $\gamma _{\rm jt}\omega _{\rm p,am}/\omega _{\rm p,jt} = \gamma _{\rm jt}^{5/2} (n_{\rm am}/n_{\rm jt})^{1/2}$ increases and the temporal growth rate is maximized when $\gamma _{\rm jt}\omega _{\rm p,am}/\omega _{\rm p,jt} = \gamma _{\rm jt}^{5/2} (n_{\rm am}/n_{\rm jt})^{1/2} = 1$. In the limit where $\gamma _{\rm jt}\omega _{\rm p,am}/\omega _{\rm p,jt} =\gamma _{\rm jt}^{5/2} (n_{\rm am}/n_{\rm jt})^{1/2} \gg 1$ relevant to the numerical simulations the phase (drift) velocity is comparable to the jet speed and the low wavenumber growth rate scales with $\gamma _{\rm jt}^{-5/4}$. The low wavenumber limit complex solution is similar in form to the hydrodynamic KHI solution at low wavenumbers (Hardee et al. 2007). In addition to the complex solution, one purely real solution for v_am = 0 and γ_am = 1 and with $\gamma _{\rm jt}^2\omega _{\rm p,am}^2 > \omega _{\rm p,jt}^2$ relevant to the numerical simulations is given by

$$\begin{equation} \omega ^2 \sim \left[ { \gamma _{\rm jt}^2\omega _{\rm p,am}^2 - \omega _{\rm p,jt}^2 \over \omega _{\rm p,am}^2 - (2-\gamma _{\rm jt}^2)\omega _{\rm p,jt}^2}\right] \omega _{\rm p,jt}^2, \end{equation} \tag{ 12 }$$

and in the high jet Lorentz factor limit where $\omega _{\rm p,am}^2 > \gamma _{\rm jt}^2\omega _{\rm p,jt}^2$ (recall that $\omega _{\rm p,jt}^2 \equiv 4 \pi n_{\rm jt} e^2/ \gamma ^3_{\rm jt} m_{\rm e}$) becomes $\omega ^2 \sim \gamma _{\rm jt}^2 \omega _{\rm p,jt}^2 = 4 \pi n_{\rm jt} e^2/ \gamma _{\rm jt} m_{\rm e}$.

In the high wavenumber limit Equation (4) becomes

$$\begin{eqnarray} &&(k^2c^2 - \omega ^2)^{1/2}[2(\omega -kv_{\rm am})^2(\omega - kv_{\rm jt})^2 - \omega _{\rm p,jt}^2 (\omega - kv_{\rm am})^2\nonumber\\ &&\quad - \omega _{\rm p,am}^2(\omega - kv_{\rm jt})^2] \sim 0. \end{eqnarray} \tag{ 13 }$$

Here the solution with ω² ∼ k²c² for electromagnetic waves formally exists only when the plasma frequencies on either side of the velocity shear are equal, i.e., γ_amω_{p, am} = γ_jtω_{p, jt}. Two additional solutions are found from

$$\begin{equation} 2(\omega -kv_{\rm am})^2(\omega - kv_{\rm jt})^2 - \omega _{\rm p,jt}^2 (\omega - kv_{\rm am})^2 - \omega _{\rm p,am}^2(\omega - kv_{\rm jt})^2 \sim 0, \end{equation} \tag{ 14 }$$

where for v_am = 0, as in the simulations, the solutions become

$$\begin{equation} \omega \sim kv_{\rm jt} \pm \omega _{\rm p,jt}/ \sqrt{2}\; {\rm and}\; \omega \sim \pm \omega _{\rm p,am}/ \sqrt{2} \end{equation} \tag{ 15 }$$

and correspond to electrostatic plasma oscillations on either side of the velocity shear surface.

Zoom In Zoom Out Reset image size

2.2.2. Numerical Solution of the Dispersion Relation for v_am = 0

A numerical solution to the dispersion relation for the Lorentz factors γ_jt = (a) 1.5, (b) 5.0, and (c) 15.0 used in the simulations is shown in Figure 2. In all cases the jet and ambient medium are assigned equal number densities determined in the ambient (simulation) frame. Thus, the plasma frequency ratios for the three cases are $\omega _{\rm p,am}/\omega _{\rm p,jt} = \gamma ^{3/2}_{\rm jt} =$ (a) 1.84, (b) 11.18, and (c) 58.09. At small and large wavenumbers the numerical solutions agree with the analytic low (Equations (11) and (12)) and high (Equation (15)) wavenumber solutions almost exactly. The low wavenumber complex solution, Equation (11), provides an excellent estimate up to wavenumbers within a factor of two of the maximally unstable wavenumber, k ≡ k*. The large wavenumber solutions, Equation (15), provide excellent estimates at wavenumbers more than a factor of two above the maximum marginally unstable wavenumber, k ≡ k^max.

The numerical solutions show that the maximum growth rates are $\omega _{\rm I}^*/\omega _{\rm p,jt} =$ (a) 0.472, (b) 0.934, and (c) 1.464 at wavenumbers k*c/ω_{p, jt} = (a) 2.344, (b) 7.079, and (c) 27.542, and wavelengths λ*(ω_{p, jt}/c) = (a) 2.68, (b) 0.888, and (c) 0.228. The maximum marginally unstable wavenumber is k^maxc/ω_{p, jt} = (a) 3.715, (b) 9.550, and (c) 33.113. If we scale the growth rate and wavelength at maximum growth to the ambient plasma frequency we obtain maximum growth rates $\omega _{\rm I}^*/\omega _{\rm p,am} =$ (a) 0.256, (b) 0.079, and (c) 0.025 at wavelengths λ* = (a) 4.93(c/ω_{p, am}), (b) 9.92(c/ω_{p, am}), and (c) 13.25(c/ω_{p, am}) and the minimum marginally unstable wavelength is λ^min = (a) 3.11(c/ω_{p, am}), (b) 7.35(c/ω_{p, am}), and (c) 11.0(c/ω_{p, am}).

Numerical solution of the dispersion relation suggests that

$$\begin{eqnarray} \omega _{\rm I}^* &\sim& 0.4 \gamma _{\rm jt}^{1/2}\omega _{\rm p,jt},\quad {\rm and} \nonumber \\ k^*v_{\rm jt} &\sim& {(1 + \gamma _{\rm jt}\omega _{\rm p,am}/\omega _{\rm p,jt})\over (\gamma _{\rm jt}\omega _{\rm p,am}/\omega _{\rm p,jt})^{1/2}} \omega _{\rm p,jt} \end{eqnarray} \tag{ 16 }$$

provide an excellent zeroth order estimate for the maximum temporal growth rate and a reasonable zeroth order estimate for the wavenumber at maximum growth. The maximum temporal growth rate estimate using Equation (16) lies within 6% of the numerically determined values. The maximally growing wavenumber estimate using Equation (16) ranges from (a) 20% above to (c) 5% below the numerically determined values where the Equation (16) estimate has been obtained by using $\omega _{\rm I}^* \sim \omega _{\rm p,jt}$ in Equation (11).

It is important to note that the maximum temporal growth rate $\omega ^*_{\rm I} \propto \gamma _{\rm jt}^{-1}$ and does not decrease as rapidly with Lorentz factor as the equal density counter-streaming maximum growth rate for which $\omega ^{\max }_{\rm I} \propto \gamma _{\rm jt}^{-3/2}$. It should be noted that this analysis assumes that the ion mass is infinite and thus may not be applicable to electron–positron cases.

Zoom In Zoom Out Reset image size

Current densities in the flow direction, J_x, associated with the two velocity shear surfaces are shown in Figure 3. In the e⁻ − p⁺ cases, J_x fluctuates in strength but not direction on either side of the velocity shear surfaces and currents run parallel to the velocity shear surface. Currents are negative on the ambient side and positive on the jet side of the velocity shear surfaces. The fluctuations are most easily seen in the blue-black on the ambient side of the velocity shear surfaces and the maximum and minimum current amplitude is J_x ∼ (a) ±6.26, (c) ±8.53, and (e) ±11.77. In the e^± cases, oblique current filaments grow at the velocity shear boundaries and the maximum and minimum current amplitude is J_x ∼ (b) ±30.1, (d) ±94.7, and (f) ±12.8. The e^± current fluctuations lead to much larger variation in the magnetic field component, B_y, associated with the velocity shear surfaces. These panels make it clear that fluctuations have the shortest spacing for the cases with γ_jt = 1.5, fluctuation spacing is about two times larger for the cases with γ_jt = 5, and also about two times larger for the cases with γ_jt = 15. In the e⁻ − p⁺ cases, current filaments are found along the velocity shear in the very early stages and outside the jet at later times (see Figure 10(d)).

Zoom In Zoom Out Reset image size

In order to make a more exact comparison with dispersion relation solutions, Figure 4 shows fluctuation in B_y relative to the average, 〈B_y〉, along one-dimensional (1D) cuts made parallel to the x axis at y = 100Δ for z/Δ = 52, 54, and 56. It is important to realize that the computational box is periodic in the x direction and only an integer number of wavelengths can fit in the computational box. In the e⁻ − p⁺ and e^± γ_jt = 1.5 cases, variation in fluctuation spacing along the x axis allows λ ∼ 50Δ ± 5Δ, i.e., λ = 1000Δ/(20 ± 1). While fluctuation amplitudes are over 10 times larger for the e^± case, to our measurement accuracy, the e⁻ − p⁺ and e^± fluctuation wavelengths are equal. In the γ_jt = 5 cases we find λ ∼ 100Δ ± 10Δ, i.e., λ = 1000Δ/(10 ± 1). Again fluctuation amplitudes are over 10 times larger for the e^± case, but to our accuracy, e⁻ − p⁺ and e^± wavelengths are equal. In the γ_jt = 15 cases, we again find that λ ∼ 100Δ ± 10Δ and fluctuation amplitudes are about 10 times larger for the e^± case.

Comparison of the observed oscillations with the theoretically predicted fastest growing wavelength, and minimum marginally unstable wavelength associated with each Lorentz factor, suggests the following interpretation. The observed oscillation wavelength for the γ_jt = 1.5 cases becomes λ^obs ∼ (4.1 ± 0.4)(c/ω_{p, am}), recall that c/ω_{p, am} = 12.2Δ, and the observed wavelength lies between the wavelengths λ* = 4.93(c/ω_{p, am}) and λ^min = 3.11(c/ω_{p, am}), predicted theoretically. The observed oscillation wavelength for the γ_jt = 5 cases becomes λ^obs ∼ (8.2 ± 0.8)(c/ω_{p, am}) and also lies between the wavelengths λ* = 9.92(c/ω_{p, am}) and λ^min = 7.35(c/ω_{p, am}), predicted theoretically. Thus both e⁻ − p⁺ and e^± γ_jt = 1.5 and γ_jt = 5 cases show compelling evidence for fluctuation wavelengths near to the theoretically predicted fastest growing wavelength. Note that the predicted minimum e-folding time of $\tau _{\rm e}^* \sim 3.9 \omega _{\rm p,am}^{-1}$ allows for ∼51 e-foldings at $t = 200 \omega _{\rm p,am}^{-1}$ for the γ_jt = 1.5 cases. For the γ_jt = 5 cases the predicted minimum e-folding time of $\tau _{\rm e}^* \sim 12.7 \omega _{\rm p,am}^{-1}$ allows for ∼20 e-foldings at $t=250 \omega _{\rm p,am}^{-1}$.

For the γ_jt = 15 cases, the theoretically predicted fastest growing wavelength and minimum marginally unstable wavelength are λ* = 13.25(c/ω_{p, am}) and λ^min = 11.0(c/ω_{p, am}), respectively. The observed oscillation wavelength of λ^obs ∼ (8.2 ± 0.8)(c/ω_{p, am}) is somewhat shorter than the predicted minimum marginally unstable wavelength but is consistent with wave growth within the predicted unstable wavelength range. We note that the minimum e-folding time of $\tau _{\rm e}^* \sim 40 \omega _{\rm p,am}^{-1}$ has only allowed for ∼7 e-foldings at $t = 300 \omega _{\rm p,am}^{-1}$ for the γ_jt = 15 cases. This is likely insufficient time for the electrostatic mode to be fully developed and fluctuation wavelengths in these cases may be influenced by the transverse current filament structure that is discussed in Section 3.

Figure 5 shows the time evolution of magnetic and electric field energy for the six simulation cases. In general, total field energy growth appears saturated for the γ_jt = 1.5 cases is still slowly growing for the γ_jt = 5 cases and remains more rapidly growing for the γ_jt = 15 cases at the simulation end time. The growth rate clearly decreases as the Lorentz factor increases but the growth time does not appear to have increased by the factor of ∼3 (γ_jt = 5) and ∼10 (γ_jt = 15) relative to the γ_jt = 1.5 cases as suggested by the maximum growth rate found from the dispersion relation solutions. In the e⁻ − p⁺ cases, the total magnetic field energy exceeds the total electric field energy by factors of ∼4 (γ_jt = 1.5) to ∼10 (γ_jt = 15). In the e^± cases, the total electric field energy is more comparable to the magnetic field energy with the total magnetic field energy exceeding the total electric field energy by only factors of ≳ 1 (γ_jt = 1.5) to ∼4 (γ_jt = 15).

Zoom In Zoom Out Reset image size

In all cases the total magnetic field energy is primarily from the B_y magnetic field component. In the e⁻ − p⁺ cases, the total electric field energy is primarily from the E_z component and secondarily from the E_x component. Field energy associated with the E_y and B_z > B_x field components is from one to three orders of magnitude smaller. In the e⁻ − p⁺ cases field energy associated with the E_z component at first grows rapidly but then is overtaken by the growth of the field energy associated with the B_y component with an accompanying slower growth and lesser field energy associated with the E_x field component. In the e^± cases, the total electric field energy is now primarily from the E_y component, secondarily from the E_z component, and only thirdly from the E_x component. The electric field energy shows rapid initial growth, much more rapid than for the e⁻ − p⁺ cases, that is eventually overtaken by the growth in the magnetic energy associated primarily with the B_y magnetic field component. Again there is not much energy associated with the B_x component, but there is now a significant amount of energy in the B_z field component relative to the B_y field component, unlike in the e⁻ − p⁺ cases.

The longitudinal kKHI mode discussed in detail in Section 2.2 would lead to growth in the E_x, E_z, and B_y field components. At least this is approximately in agreement with what we find for the e⁻ − p⁺ cases, although the growth time does not increase as rapidly with the Lorentz factor as predicted. On the other hand, in the e^± cases, we find the electric field energy dominated by the E_y component and a significant amount of magnetic field energy in the B_z component. The fact that the growth of the total magnetic and electric field energies is not decreasing with the Lorentz factor as rapidly as predicted by the longitudinal dispersion relation and that, particularly in the e^± cases, significant magnetic and electric field components develop that are not described by the analysis in Section 2.2, which provides evidence for additional processes operating in the velocity shear region that have not been captured by a longitudinal dispersion relation, and, in particular, the magnetic and electric fields imply the presence of growing transverse modes.

Zoom In Zoom Out Reset image size

Figure 6 shows the structure of the B_y component of the magnetic field in the y–z plane (jet flows out of the page) at the midpoint of the simulation box, x = 500Δ, and 1D cuts along the z axis showing the magnitude and direction of the magnetic field components at the midpoint of the simulation box, x = 500Δ and y = 100Δ for the e⁻ − p⁺ case and the e^± case at simulation time $t = 300\,\omega _{\rm pe}^{-1}$, both with γ_jt = 15. Comparison of the transverse structures in the y direction at the velocity shear surfaces shown in panels (a) and (d) with the parallel structures in the x direction shown in Figure 1 in panels (b) and (c) shows that the fluctuations transverse to the jet in the y direction are much more rapid than fluctuations along the jet in the x direction. In the e⁻ − p⁺ case, magnetic fields appear relatively uniform at the velocity shear surfaces along the transverse y direction just as was seen at the velocity shear surfaces along the parallel x direction, with almost no transverse fluctuations visible in the magnetic field (small fluctuations in the y direction over distances on the order of ∼10Δ are visible in the currents in Figure 7(b)), whereas small longitudinal mode fluctuations in the x direction occur over distances ∼100Δ. For the electron–positron case, the magnetic field alternates in both the y and z directions and these transverse fluctuations occur over distances on the order of ∼10Δ, whereas longitudinal mode fluctuations in the x direction occur over distances ∼100Δ.

Zoom In Zoom Out Reset image size

The 1D cuts show that the B_y field component dominates in the e⁻ − p⁺ case, that the B_y field component is about an order of magnitude smaller for the e^± case, and that the B_z component is significant for the e^± case, as already indicated in Figure 5. The 1D cuts also show that there is magnetic field sign reversal on either side of the maximum that is relatively small for the e⁻ − p⁺ case but is much more significant for the e^± case, which can be seen also in Figure 6(d). More details are revealed by the enlargement of the region contained in the squares. For the e⁻ − p⁺ case, the generated relatively uniform DC magnetic field is symmetric about the velocity shear surface, e.g., note that B_y > 0 immediately around the shear surface and B_y < 0 in the jet and ambient plasmas at somewhat larger distances from the shear surface. On the other hand, for the e^± case the generated AC magnetic field resides largely on the jet side of the velocity shear surface.

Figure 7 shows how the J_x current structure in a small y–z plane, responsible for the magnetic field structure shown in Figure 6. Motion of electrons and/or positrons across the shear surface produces the electric currents shown also in Figure 7 by the arrows. Relativistic jet flow is out of the page and in the e⁻ − p⁺ case positive (red/orange) and negative (blue/black) current flows along the jet and the sheath side of the velocity shear surfaces, respectively. Positive currents are stronger than the negative currents, leading to the generation of the B_y magnetic field component, shown in Figures 6(a)–(c). In the e^± case, a complex current structure appears on the jet side of the velocity shear surface. The associated magnetic fields are then folded and twisted by vortical plasma motions. The vortices appear like "islands" in the magnetic field. In the currents, it is possible to see that the transverse fluctuation scale is similar in the e⁻ − p⁺ and e^± cases, but the structures are considerably different.

It seems likely that the development of transverse filamentary structure has influenced the longitudinal structure studied in Section 2. In general, we find that the kKHI grows on timescales t∝γ_jt, albeit growth also depends on the density ratio across the velocity shear. Once particles have scattered across the velocity shear via kKHI or thermal motions, structure associated with interpenetrating relativistic plasmas can develop. For k ≡ k_x, there is the beam space charge instability (Bludeman et al. 1960 and see also Equations (8) and (9) in Hardee & Rose 1978) with

$$\begin{eqnarray} \omega &=& kv_{\rm jt} \pm i {\omega _{\rm p,jt}\over \omega _{\rm p,am}} kv_{\rm jt},\quad {\rm and} \nonumber \\ \omega _{\rm I}^{\max } &=& 0.69 (\omega _{\rm p,jt}^2 \omega _{\rm p,am})^{1/3}\, {\rm at}\; kv_{\rm jt} = \omega _{\rm p,am} - 0.4(\omega _{\rm p,jt}^2\omega _{\rm p,am})^{1/3}.\nonumber\\ \end{eqnarray} \tag{ 17 }$$

For k ≡ k_y there is the ordinary mode (filamentation) instability with

$$\begin{equation} \omega _{\rm R} = 0,\; {\rm and}\; \omega _{\rm I} = {{\omega _{\rm p,am}}\over {(\omega _{\rm p,am}^2 + \omega _{\rm p,jt}^2)^{1/2}}} (\gamma _{\rm jt} \omega _{\rm p,jt}) {v_{\rm jt} \over c}. \end{equation} \tag{ 18 }$$

Equation (18) is formally found in the limit $k^2c^2 \gg 2\omega _{\rm p,am}^2 + \gamma _{\rm jt}^2\omega _{\rm p,jt}^2 + \omega _{\rm p,jt}^2$ and is like Equation (11) in Hardee & Rose (1978), which assumed ω_{p, am} ≫ ω_{p, jt}, but now with Ω = eB/mc = 0 and allowing for ω_{p, jt} ∼ ω_{p, am} to reveal the density dependence. We have adopted the present notation in Equations (17) and (18) and note that they were derived originally in the context of electron–positron jet and ambient plasmas. They should also apply to electron–proton plasmas where the ions are assumed infinitely massive.

We see from the above that the longitudinal beam space charge instability grows on timescales t∝γ_jt that are comparable to the kKHI. On the other hand, filamentary structure associated with the transverse ordinary mode instability grows on timescales $t \propto \gamma _{\rm jt}^{1/2}$ and thus can grow faster than kKHI longitudinal structure for large Lorentz factors. While there is excellent agreement between observed longitudinal structure scales and theoretical prediction for the two lower Lorentz factors, such is not the case for the high Lorentz factor simulation. Here we believe that more rapid growth of transverse structure in the high Lorentz factor case has overwhelmed slower growth of the longitudinal kKHI and led to a longitudinal length scale that is less than the minimum unstable wavelength associated with the kKHI.

Zoom In Zoom Out Reset image size

Figure 8 shows how the J_x current structure in a small y–z plane square around the velocity shear surface (see Figure 7 for location) and the magnetic field strength and position relative to the shear surfaces along 1D cuts in the z direction change as a function of the Lorentz factor for the e^± cases at time $t = 300 \omega _{\rm pe}^{-1}$. Not presented here are the e⁻ − p⁺ cases, as they show very little change in the amplitude or the width of the amplified magnetic field region as a function of the Lorentz factor. This result may indicate a difference from the theory and counter-streaming simulation results (Grismayer et al. 2013a, 2013b; Alves et al. 2014) in which amplitude and width scaled with $\sqrt{\gamma _0}$ or just that our higher Lorentz factor e⁻ − p⁺ cases have not reached saturation. In the e^± cases, we see that the currents are located on the jet side of the shear surface for Lorentz factors γ_jt = (a) 15 and (b) 5, but are located on both sides of the shear surface for (c) γ_jt = 1.5. The maximum and minimum current density amplitudes are (a) ±6.26, (b) ±77.0, and (c) ±19.5, and the maximum magnetic field strength is smaller by about an order of magnitude for the (a) γ_jt = 15 case and smaller by about a factor of three for the (c) γ_jt = 1.5 case compared to the (b) γ_jt = 5 case. Here we do find an increase in the maximum field strength from the γ_jt = 1.5 case to the γ_jt = 5 case as suggested by the theory and counter-streaming results but with a decrease in the total shear layer width instead of the expected increase in shear layer width. The γ_jt = 15 simulation is not near saturation so cannot be directly compared to the lower Lorentz factor cases. However, it is clear that the increased inertia of the more relativistically moving jet plasma inhibits motion of jet electrons across the shear surface and affects the shear structure significantly compared to counter-streaming simulations in which both plasmas have the same inertia.

Temporal development of the total magnetic field energy shown in Figure 5 shows that the still slowly growing magnetic field energy for the (b) γ_jt = 5 case is comparable to the (c) γ_jt = 1.5 case at time $t = 300 \omega _{\rm pe}^{-1}$ and should become greater at later times. Figure 5 also shows a more rapid growth of the total magnetic field energy for the (a) γ = 15 case at this time. These results suggest that the differences in the current structure and the magnetic field strength and location may indicate a temporal development attributable to growth rate differences in addition to inertial differences. In fact, for the case with γ_jt = 1.5 at $t = 200 \omega _{\rm pe}^{-1}$, the current filaments with maximum and minimum values ±42.6 are located nearer to the velocity shear than at the later time of $t = 300 \omega _{\rm pe}^{-1}$ shown in Figure 8(c), where the maximum and minimum values are ±19.5. Thus, the differences seen in Figure 8 from high to low Lorentz factors may provide an indication of the temporal development of the current structure from fewer (γ_jt = 15) to more (γ_jt = 1.5) e-folding times.

Figure 9 provides a 3D display of the currents and magnetic fields at the velocity shear surface for the e^± case with γ_jt = 5 at $t = 250 \omega _{\rm pe}^{-1}$. The 3D display reveals current filaments with length in the x direction (see Figure 3(d)) much longer than the spacing in the y–z plane (see Figure 8(b)). Strong positive (red) and negative (blue) current filaments wrapped by magnetic field lines seen in 3D are seen in a 2D slice shown previously in Figure 8(b) (albeit here at an earlier time) in the x component of the current density (positive (red) and negative (blue)) and magnetic field (arrows) in the y–z plane. The positive and negative current filaments seen in 2D are now revealed to twist around each other with the longitudinal wavelength λ^obs ∼ 100Δ seen in Figure 3(d) and Figure 4(d). The total magnetic energy isosurface shows a concentration of the magnetic field around the current filaments.

Zoom In Zoom Out Reset image size

Figure 10 provides a 3D display of the currents and magnetic fields at the velocity shear surface for the e^± cases at the same time, $t = 300 \omega _{\rm pe}^{-1}$, as the 2D slices shown in Figure 8, and also shows the current and magnetic field structure at the velocity shear surface for the e⁻ − p⁺ case with γ_jt = 5. For the e^± cases, as indicated by the 2D slices, the 3D structure shows a current carrying region that thickens as the Lorentz factor decreases and at low Lorentz factor appears on both sides of the velocity shear surface. The 3D structure suggests a single layer of current filaments at γ_jt = 15 that broadens to a dual layer of current filaments at γ_jt = 5. At γ_jt = 1.5 current filaments on both sides of the velocity shear layer are much less well defined. A comparison between the γ_jt = 5 e^± (Figure 10(b)) and e⁻ − p⁺ (Figure 10(d)) cases shows the very different current and magnetic field structures at the velocity shear surface. For the e⁻ − p⁺case, the magnetic field is very uniform and largely confined to the velocity shear surface just below a strong positive (red) current sheet on the jet side. Outside the velocity shear surface we see a weaker negative (blue/green) current sheet, and further outside a filamented weak negative (blue/green) current region.

Zoom In Zoom Out Reset image size

The structures shown in Figures 9 and 10 are similar to those produced by the filamentation (Weibel-like) instability, associated with interpenetrating plasmas (see Equation (18)). We note that the change in structure from closely spaced current filaments of smaller diameter in a narrower region in the γ_jt = 15 case to the merged larger-diameter and less closely spaced filaments in a broader region in the γ_jt = 5 case shown in Figures 8 and 10 is like the expected temporal development for the filamentation instability as the number of e-folding times increases. Since in our simulations the $\gamma _{\rm jt} = 5\, \&\, 15$ cases have not reached saturation, one cannot say with certainty that the instability will not ultimately develop the structures seen in the saturated γ_jt = 1.5 case, in which the current filaments are probably fully developed by $t = 300 \omega _{\rm pe}^{-1}$. However, it seems more likely that the lack of significant current structure on the outside of the velocity shear surface in the two higher Lorentz factor e^± cases is a direct result of the increased inertia of the relativistically moving plasma.

The observed 2D and 3D structures indicate the development of longitudinal (electrostatic two-stream) and transverse (Weibel-like current filamentation, e.g., Nishikawa et al. 2005, 2006, 2009a) plasma instabilities usually associated with interpenetrating relativistic plasmas. Figure 11 shows the magnetic field components and the associated phase-space plots of electron motions in the z direction perpendicular to the velocity shear surface. These motions produce the mixing required to trigger interpenetrating plasma instabilities. Note that in e⁻ − p⁺ cases, mixing is almost completely associated with the electrons and in e^± cases both electrons and positrons participate in the mixing.

Zoom In Zoom Out Reset image size

The e⁻ − p⁺ γ_jt = 5 case at two different times illustrates the development of both the dominant B_y component of the magnetic field and the plasma mixing process. The magnetic field is initially strongest at the shear surface. The magnetic field strengthens and more deeply penetrates the ambient plasma with time. Slightly deeper penetration into the jet plasma also occurs with time. Ambient electrons (red dots) with v_z > 0 are moving toward and into the jet, and become more heated and penetrate deeper into the jet plasma with time. Relatively cold (note a very small thermal spread) jet electrons (blue dots) with v_z < 0 are moving outward and into the ambient plasma. These jet electrons, while remaining cold, penetrate deeper into the ambient plasma with time. At the simulation times presented, the electrons are mixed in space but not yet in velocity. Due to the relatively uniform DC magnetic field generated in the x and y directions, the phase-space plot shows a regular structure.

In the e^± γ_jt = 5 case, the dominant B_y component of the magnetic field is strongest at the shear surface and more deeply penetrates the jet plasma than the ambient plasma. Note the very different location of the magnetic field in this case versus the e⁻ − p⁺ case at the same simulation time. The electrons are less mixed spatially on the ambient side of the shear surface but with much more heating of the jet electrons than in the e⁻ − p⁺ case. Here most of the action resides on the jet side of the shear surface where the filamentation instability dominates the dynamics. Both ambient and jet electrons are accelerated in the strong AC magnetic and electric fields associated with the filamentation instability. Ambient electrons are more strongly heated than the jet electrons, but now there is a significant velocity mixing of the jet and ambient electrons and the ambient electrons penetrate into the jet more deeply than in the e⁻ − p⁺ case.

Zoom In Zoom Out Reset image size

Just as Figure 11 shows that the particle behavior near the velocity shear for the e^± and e⁻ − p⁺ cases is significantly different, Figure 12 shows that electron acceleration at the velocity shear also is significantly different in e^± and e⁻ − p⁺ cases. The v_x–v_y velocity component plots show that the electrons are accelerated in parallel and perpendicular directions (Figure 12(b)) in the e^± case and the electrons are accelerated only in the parallel direction (Figure 12(c)) in the e⁻ − p⁺ case.

Similar analysis for the two other jet Lorentz factors (not shown) shows that e⁻ − p⁺ cases with different Lorentz factor show the same parallel electron acceleration. On the other hand, there are some modest differences in the e^± cases for the two other jet Lorentz factors studied. These differences occur because in the γ_jt = 1.5 case the ambient and jet electrons are strongly mixed across the velocity shear surface, but in the γ_jt = 15 case, the ambient electrons are strongly mixed into the jet region, but jet electrons are only weakly mixed into the ambient plasma. This velocity and phase-space result is also revealed in Figures 8 and 10 that show a development of current filament structure on both sides of the velocity shear surface for the γ_jt = 1.5 case but only on one side of the velocity shear surface for the γ_jt = 15 case. In the γ_jt = 1.5 case, the mixing of jet and ambient electrons on both sides of the velocity shear is accompanied by mixing in v_x − v_y and acceleration in both parallel and perpendicular directions like the γ_jt = 5 case. In the γ_jt = 15 case, ambient and jet electrons are accelerated mainly in the perpendicular direction.

Our simulations do not follow the kKHI significantly past the saturation phase (see Section 3.1) and are thus too short to allow for significant electron acceleration in the self-generated fields of the shear flow instabilities. In fact, the strongest acceleration should occur after the growth of the kKHI fully saturates and electrons can probe the long-lived turbulent electric fields in the shear region. This has been demonstrated by Alves et al. (2014) for the electron–proton plasma, who discuss the acceleration process and estimate the maximum energy gain of an electron interacting in the electric field to be $\Delta \mathcal {E} _{\max } \propto mc^{2}\gamma _{0}^{4}$. Acceleration to super-thermal energies is thus possible for shear flows with relativistic Lorentz factors as shown by the PIC simulations in Alves et al. (2014). In our simulations, we trace the initial stages of electron acceleration. We note that the acceleration in the parallel direction observed in our e⁻ − p⁺ simulations resemble a process of electron surfing on the electric field structures in the shear region that provides acceleration mostly in the direction of the bulk plasma flow, as indicated in Alves et al. (2014).