Nonlinear energy transfer and current sheet development in localized Alfven wavepacket collisions in the strong turbulence limit

TL;DR

This paper demonstrates through gyrokinetic simulations that localized Alfven wavepacket collisions in space plasmas effectively transfer energy to smaller scales and generate current sheets, even when kinetic effects are considered.

Contribution

It extends previous idealized models by showing that key properties of energy transfer and current sheet formation persist in more realistic wavepacket collisions including kinetic effects.

Findings

Wavepacket collisions facilitate perpendicular energy cascade.

Current sheets are formed during wavepacket overlap.

Nonlinear interactions cease after wavepackets separate.

Abstract

In space and astrophysical plasmas, turbulence is responsible for transferring energy from large scales driven by violent events or instabilities, to smaller scales where turbulent energy is ultimately converted into plasma heat by dissipative mechanisms. The nonlinear interaction between counterpropagating Alfven waves, denoted Alfven wave collisions, drives this turbulent energy cascade, as recognized by early work with incompressible magnetohydrodynamic (MHD) equations. Recent work employing analytical calculations and nonlinear gyrokinetic simulations of Alfven wave collisions in an idealized periodic initial state have demonstrated the key properties that strong Alfven wave collisions mediate effectively the transfer of energy to smaller perpendicular scales and self-consistently generate current sheets. For the more realistic case of the collision between two initially separated Alfven wavepackets, we use a nonlinear gyrokinetic simulation to show here that these key properties persist: strong Alfven wavepacket collisions indeed facilitate the perpendicular cascade of energy and give rise to current sheets. Furthermore, the evolution shows that nonlinear interactions occur only while the wavepackets overlap, followed by a clean separation of the wavepackets with straight uniform magnetic fields and the cessation of nonlinear evolution in between collisions, even in the gyrokinetic simulation presented here which resolves dispersive and kinetic effects beyond the reach of the MHD theory.

Figures

Equations

$$\frac{\partial {\mathbf{z}}^{\pm }}{\partial t}\mp {\mathbf{v}}_A\cdot \nabla {\mathbf{z}}^{\pm }=-{\mathbf{z}}^{\mp }\cdot \nabla {\mathbf{z}}^{\pm }-\nabla P/(n_{0i}{m}_i+n_{0e}{m}_e),$$

$$\delta {\hat{B}}_z(k_x,k_y,z)=\left[\frac{\delta {\hat{\overline{B}}}_z(k_x,k_y,{\beta }_i,T_i/T_e)}{\delta {\hat{\overline{B}}}_y(k_x,k_y,{\beta }_i,T_i/T_e)}\right]\delta {\hat{B}}_y(k_x,k_y,z)$$

$$\delta {\hat{B}}_j(k_x,k_y,z)=\sum _{k_z}\delta {\hat{B}}_j(k_x,k_y,k_z)e^{i{k}_{z}z},$$

$$\sum _{k_z}\delta {\hat{B}}_z(k_x,k_y,k_z)e^{i{k}_{z}z}=\sum _{k_z}\left[\frac{\delta {\hat{\overline{B}}}_z(k_x,k_y)}{\delta {\hat{\overline{B}}}_y(k_x,k_y)}\right]\delta {\hat{B}}_y(k_x,k_y,k_z)e^{i{k}_{z}z},$$

$$\delta {\hat{B}}_j(k_x,k_y,k_z)=\left[\frac{\delta {\hat{\overline{B}}}_z(k_x,k_y,{\beta }_i,T_i/T_e)}{\delta {\hat{\overline{B}}}_y(k_x,k_y,{\beta }_i,T_i/T_e)}\right]\delta {\hat{B}}_y(k_x,k_y,k_z).$$

$$w(z)=\exp \left[{\left(\frac{z-z_0}{{\Delta }_z}\right)}^p\right]+\exp \left[{\left(\frac{z-z_0+L_z}{{\Delta }_z}\right)}^p\right]+\exp \left[{\left(\frac{z-z_0-L_z}{{\Delta }_z}\right)}^p\right]$$

$$f(z)=\cos (k_{z}z+\delta ),$$

Full text

††volume: 650

{PACS}

1 Introduction

Turbulence plays an important role in facilitating particle energization and plasma heating in space and astrophysical plasmas, influencing the macroscopic evolution of many poorly understood systems, such as the solar corona and solar wind, planetary magnetospheres, and black hole accretion disks. The turbulent cascade mediates the transfer of energy from magnetic fields and plasma flows at large scales down to much smaller scales where dissipation mechanisms can effectively remove energy from the turbulent fluctuations, ultimately converting that energy to plasma heat. Understanding the details of this nonlinear turbulent cascade to small scales and of the mechanisms by which the turbulent energy is thermalized represents a grand challenge problem in heliophysics and astrophysics.

Early research on incompressible magnetohydrodynamics (MHD) turbulence in the 1960s (Iroshnikov, 1963; Kraichnan, 1965) suggested that nonlinear interactions between counterpropagating Alfvén waves—or Alfvén wave collisions—support the turbulent cascade of energy from large to small scales. Following significant previous studies on weak incompressible MHD turbulence (Sridhar & Goldreich, 1994; Montgomery & Matthaeus, 1995; Ng & Bhattacharjee, 1996; Galtier et al., 2000), recent work has elucidated the mechanism of energy transfer in Alfvén wave collisions in the weakly nonlinear limit by computing an asymptotic analytical solution using incompressible MHD (Howes & Nielson, 2013), verifying that solution using nonlinear gyrokinetic simulations in the MHD limit of perpendicular scales larger than the Larmor radius, $k_{\perp}\rho_{i}\gg 1$ (Nielson et al., 2013), and confirming the results experimentally in the laboratory (Howes et al., 2012, 2013; Drake et al., 2013, 2014, 2016). The derivation of the analytical solution was possible using the idealized initial conditions of two overlapping, perpendicularly polarized Alfvén waves in a periodic geometry, and solving for the nonlinear evolution of the system. For initial plane Alfvén waves with wavevectors \mathbf{k}_{1}^{+}=k_{\perp}\mbox{\hat{\mathbf{x}}}-k_{\parallel}\mbox{\hat{\mathbf{z}}} and \mathbf{k}_{1}^{-}=k_{\perp}\mbox{\hat{\mathbf{y}}}+k_{\parallel}\mbox{\hat{\mathbf{z}}}, the nonlinear energy transfer is mediated by a nonlinearly generated, purely magnetic mode with wavevector \mathbf{k}_{2}^{(0)}=k_{\perp}\mbox{\hat{\mathbf{x}}}+k_{\perp}\mbox{\hat{\mathbf{y}}}, which can be interpreted as an oscillating shear in the magnetic field along which the Alfvén waves propagate (Maron & Goldreich, 2001; Howes & Bourouaine, 2017). The nonlinear interaction between $\mathbf{k}_{1}^{\pm}$ and $\mathbf{k}_{2}^{(0)}$ yields a secular transfer of energy from the $\mathbf{k}_{1}^{+}$ Alfvén wave to an Alfvén wave with \mathbf{k}_{3}^{+}=2k_{\perp}\mbox{\hat{\mathbf{x}}}+k_{\perp}\mbox{\hat{\mathbf{y}}}-k_{\parallel}\mbox{\hat{\mathbf{z}}}, and from the $\mathbf{k}_{1}^{-}$ Alfvén wave to an Alfvén wave with \mathbf{k}_{3}^{-}=k_{\perp}\mbox{\hat{\mathbf{x}}}+2k_{\perp}\mbox{\hat{\mathbf{y}}}+k_{\parallel}\mbox{\hat{\mathbf{z}}}. Since the energy is transferred to an Alfvén wave with a higher perpendicular wavenumber, $|\mathbf{k}_{3_{\perp}}^{\pm}|>|\mathbf{k}_{1\perp}^{\pm}|$, this interaction represents the fundamental mechanism by which turbulence transfers energy from larger to smaller scales.

Another important discovery about plasma turbulence followed from the finding that the nonlinear evolution of MHD turbulence simulations leads to the development of intermittent current sheets (Matthaeus & Montgomery, 1980; Meneguzzi et al., 1981), and that the dissipation of turbulent energy is found to be largely concentrated in the vicinity of these intermittent current sheets (Uritsky et al., 2010; Osman et al., 2011; Zhdankin et al., 2013). This finding has motivated significant recent efforts to seek evidence of the spatial localization of plasma heating by the dissipation of turbulence in current sheets through statistical analyses of solar wind observations (Osman et al., 2011; Borovsky & Denton, 2011; Osman et al., 2012; Perri et al., 2012; Wang et al., 2013; Wu et al., 2013; Osman et al., 2014) and numerical simulations (Wan et al., 2012; Karimabadi et al., 2013; TenBarge & Howes, 2013; Wu et al., 2013; Zhdankin et al., 2013). Although these works clearly demonstrate a connection between current sheets and plasma heating, the origin of these current sheets in plasma turbulence remains unknown: do they represent advected flux tube boundaries (Borovsky, 2008, 2010), or are they generated dynamically by the turbulence itself (Boldyrev et al., 2011; Zhdankin et al., 2012)? A significant breakthrough on this question was the discovery that Alfvén wave collisions in the strong turbulence limit naturally generate current sheets (Howes, 2016), making a connection for the first time between the nonlinear mechanism governing the transfer of energy to small scales and the self-consistent development of intermittent current sheets.

The analytical solution for energy transfer in weak Alfvén wave collisions (Howes & Nielson, 2013) and the simulations showing that strong Alfvén wave collisions naturally generate current sheets (Howes, 2016) were based on an idealized initial condition in which two finite-amplitude, plane Alfvén waves are initially overlapping in a periodic geometry before they begin to interact nonlinearly. Here we eliminate the unrealistic aspect of those studies by simulating the strong nonlinear interactions between two initially separated Alfvén wavepackets, with the aim to determine whether the general properties of Alfvén wave collisions found in the idealized case persist in this more realistic case of colliding wavepackets. Specifically, we focus here on answering two questions: (i) Do collisions between Alfvén wavepackets still mediate the transfer of energy to small perpendicular scales?; and (ii) Do Alfvén wavepacket collisions in the strongly nonlinear limit still lead to the development of intermittent current sheets? A companion paper (Verniero & Howes, 2017) will investigate whether the nonlinearly generated $(k_{x}/k_{\perp},k_{y}/k_{\perp},k_{z}/k_{\parallel})=(1,1,0)$ mode still mediates the energy transfer in the weakly collisional limit of Alfvén wavepacket collisions.

In §2, we describe the setup of this strong, localized Alfvén wavepacket collision simulation. The nonlinear evolution of this simulation is analyzed in §3, with particular emphasis on the perpendicular cascade of energy in §3.1, current sheet development in §3.2, evolution of the energy in §3.3, and general qualitative properties of localized Alfvén wavepacket collisions in §3.4. Conclusions are presented in §4.

2 Simulation

Here we employ the Astrophysical Gyrokinetics code AstroGK (Numata et al., 2010) to perform a gyrokinetic simulation of the nonlinear interaction between two initially separated, counterpropagating Alfvén wavepackets in the strongly nonlinear limit.

AstroGK evolves the perturbed gyroaveraged distribution function $h_{s}(x,y,z,\lambda,\varepsilon)$ for each species $s$, the scalar potential $\varphi$, the parallel vector potential $A_{\parallel}$, and the parallel magnetic field perturbation $\delta B_{\parallel}$ according to the gyrokinetic equation and the gyroaveraged Maxwell’s equations (Frieman & Chen, 1982; Howes et al., 2006). Velocity space coordinates are $\lambda=v_{\perp}^{2}/v^{2}$ and $\varepsilon=v^{2}/2$. The domain is a periodic box of size $L_{\perp}^{2}\times L_{z}$, elongated along the straight, uniform mean magnetic field \mathbf{B}_{0}=B_{0}\mbox{\hat{\mathbf{z}}}, where all quantities may be rescaled to any parallel dimension satisfying $L_{z}/L_{\perp}\gg 1$. Uniform Maxwellian equilibria for ions (protons) and electrons are chosen, with a realistic mass ratio $m_{i}/m_{e}=1836$. Spatial dimensions $(x,y)$ perpendicular to the mean field are treated pseudospectrally; an upwind finite-difference scheme is used in the parallel direction, $z$. Collisions employ a fully conservative, linearized collision operator with energy diffusion and pitch-angle scattering (Abel et al., 2008; Barnes et al., 2009).

The functional forms along $z$ of the initial Alfvén wavepackets used in this simulation are shown in Fig. 1. Initially at $z<0$ is a wavepacket with a magnetic field perturbation $\delta B_{y}$ polarized in the $y$ direction, a unipolar variation in $z$ (red), and a perpendicular structure with wavenumber $\mathbf{k}_{\perp}^{-}=(k_{x}L_{\perp},k_{y}L_{\perp})=(1,0)$. The eigenfunction dictating the different field components and perturbed distribution functions for this Alfvén wavepacket is determined by solving the linear, collisionless, gyrokinetic dispersion relation (Howes et al., 2006) for an Alfvén mode with the chosen perpendicular Fourier wavevector $(k_{x},k_{y})$ to obtain the complex eigenfrequency $\omega$, the complex Fourier coefficients for the eigenfunctions of the electromagnetic potentials $\hat{\phi}$, $\hat{A}_{\parallel}$ and $\delta\hat{B}_{\parallel}$, and the complex perturbed gyrokinetic distribution functions for the ions $\hat{h}_{i}(v_{\parallel},v_{\perp})$ and electrons $\hat{h}_{e}(v_{\parallel},v_{\perp})$, where the hat symbol denotes the $(k_{x},k_{y})$ Fourier coefficient (Nielson et al., 2013). The procedure for localizing this Alfvén wavepacket in the $z$ direction is described in Appendix A. Through this procedure, this unipolar Alfvén wavepacket in Fig. 1 propagates in the $+z$ direction. The other Alfvén wavepacket, initially at $z>0$, has a magnetic field perturbation $\delta B_{x}$ polarized in the $x$ direction, a dipolar structure in $z$ (blue), and a perpendicular structure with wavenumber $\mathbf{k}_{\perp}^{+}=(k_{x}L_{\perp},k_{y}L_{\perp})=(0,1)$; the eigenfunction specified by the same procedure dictates that this dipolar wavepacket propagates in the $-z$ direction. Note that these unsymmetric initial conditions were chosen to further test the complexities in Alfvén wave collisions. In particular by having a unipolar wavepacket, the collision does not shear and unshear in an oscillatory manner as it does in the interaction between two dipolar wavepackets.

The plasma parameters for this strong Alfvén wavepacket collision simulation are ion plasma beta $\beta_{i}=1$ and ion-to-electron temperature ratio $T_{i}/T_{e}=1$. To study the nonlinear evolution of this Alfvén wavepacket collision in the limit $k_{\perp}\rho_{i}\ll 1$, we choose a perpendicular simulation domain size $L_{\perp}=40\pi\rho_{i}$ with simulation resolution $(n_{x},n_{y},n_{z},n_{\lambda},n_{\varepsilon},n_{s})=(64,64,128,32,32,2)$. Therefore, the initial Alfvén wavepackets have perpendicular wavevectors $\mathbf{k}_{\perp}^{-}=(k_{x}\rho_{i},k_{y}\rho_{i})=(0.05,0)$ for the unipolar wavepacket and $\mathbf{k}_{\perp}^{+}=(k_{x}\rho_{i},k_{y}\rho_{i})=(0,0.05)$ for the dipolar wavepacket, so both waves have the same initial perpendicular wavenumber $k_{\perp}^{\pm}\rho_{i}=0.05$, but are polarized perpendicular to each other. The fully resolved perpendicular range in this dealiased pseudospectral method covers $0.05\leq k_{\perp}\rho_{i}\leq 1.05$. Here the ion thermal Larmor radius is $\rho_{i}=v_{ti}/\Omega_{i}$, the ion thermal velocity is $v_{ti}^{2}=2T_{i}/m_{i}$, the ion cyclotron frequency is $\Omega_{i}=q_{i}B_{0}/(m_{i}c)$, and the temperature is given in energy units. The parallel length of the simulation domain is $L_{z}$, extending over the range [$-L_{z}/2$,$L_{z}/2$]. Note that the simulation domain is triply periodic, so when a wavepacket exits the domain at $z=\pm L_{z}/2$, it re-enters at the opposite end at $z=\mp L_{z}/2$, enabling these two wavepackets to undergo successive collisions with each other. The linearized Landau collision operator (Abel et al., 2008; Barnes et al., 2009) is employed with collisional coefficients $\nu_{i}=\nu_{e}=10^{-3}$, yielding weakly collisional dynamics with $\nu_{s}/\omega\ll 1$.

The amplitude of the initial wavepackets is parameterized by the nonlinearity parameter (Goldreich & Sridhar, 1995), defined by taking the ratio of the magnitudes of the linear to the nonlinear terms in the incompressible MHD equations (Howes & Nielson, 2013; Nielson et al., 2013). In terms of Elsasser variables, defined by $\mathbf{z}^{\pm}=\mathbf{u}\pm\delta\mathbf{B}/\sqrt{4\pi(n_{0i}m_{i}+n_{0e}m_{e})}$, the nonlinearity parameter is defined by

$\chi^{\pm}\equiv|\mathbf{z}^{\mp}\cdot\nabla\mathbf{z}^{\pm}|/|\mathbf{v}_{A}\cdot\nabla\mathbf{z}^{\pm}|$, where $\chi^{\pm}$ characterizes the strength of the nonlinear distortion of the $\mathbf{z}^{\pm}$ Alfvén wave by the counterpropagating $\mathbf{z}^{\mp}$ Alfvén wave. For the particular initial Alfvén wavepackets shown in Fig. 1, the nonlinearity parameter simplifies to $\chi^{\pm}=2k_{\perp}\delta B_{\perp}^{\mp}/(k_{\parallel}B_{0})$. With the $\mathbf{z}^{\pm}$ wavepackets having parallel wavenumbers of approximately $k_{\parallel}a_{0}=\mp 3$, where $a_{0}=L_{z}/2\pi$, the amplitude of the unipolar wavepacket $(\delta B_{\perp}^{-}/B_{0})(a_{0}/\rho_{i})\simeq 60$ gives $\chi^{+}=2$ and the amplitude of the dipolar wavepacket $(\delta B_{\perp}^{+}/B_{0})(a_{0}/\rho_{i})\simeq 40$ gives $\chi^{-}=1.3$. Strong, critical balanced turbulence (Goldreich & Sridhar, 1995) corresponds to a nonlinearity parameter of $\chi\sim 1$, so this simulation falls into the desired limit of strong Alfvén wavepacket collisions.

3 Evolution of the Nonlinear Interaction

The basic evolution of this strong Alfvén wavepacket collision simulation is illustrated by three-dimensional contour plots of the parallel current density $j_{z}$ associated with each of the interacting wavepackets, shown in Fig. 2 and displayed as a movie in Fig. 3. Time is normalized in terms of the time for a single Alfvén wavepacket collision, $T_{c}$, during which the initially separated Alfvén wavepackets approach each other along $z$ (with the $+z$ direction from left to right in Fig. 2), overlap and interact nonlinearly, and then move away from each other after the collision. At $t=0$, the unipolar wavepacket is centered at $z=-L_{z}/4$, and the dipolar wavepacket at $z=+L_{z}/4$; during a collision time $T_{c}$, each wavepacket propagates at the Alfvén velocity $v_{A}$ over a parallel distance $L_{z}/2$. The midpoint of each collision occurs at $t/T_{c}=0.5,1.5,2.5$.

In Fig. 2(a) are plotted isocontours of the normalized parallel current density $j_{z}/j_{0}$ for each of the wavepackets at $t=0$, where the current density is normalized by $j_{0}=n_{0}q_{i}v_{ti}L_{\perp}/L_{z}$. It is clear that the unipolar wavepacket (at $z<0$, left side) initially has only perpendicular variation in the $x$ direction, while the dipolar wavepacket (at $z>0$, right side) initially has only perpendicular variation in the $y$ direction. Fig. 2(b) shows the midpoint of the first collision occurring at $z=0$ and $t/T_{c}=0.5$, showing a significantly more complicated perpendicular structure parallel current density $j_{z}$.

After the first collision at $t/T_{c}=1.0$ in Fig. 2(c), the unipolar wavepacket, now at $z=+L_{z}/4$, has gained some variation in the $y$ direction, and the dipolar wavepacket, now at $z=-L_{z}/4$, has developed variation in the $x$ direction. Each wavepacket has been distorted by passing through, and interacting nonlinearly with, the counterpropagating wavepacket. Mathematically, when expressed in terms of Fourier modes, the strong Alfvén wavepacket collision has mediated the nonlinear transfer of energy from the two initial perpendicular Fourier modes to other Fourier modes with larger values of $k_{\perp}$, as shown quantitatively in §3.1. Therefore this visualization clearly shows the nonlinear cascade of energy to smaller scales in strong Alfvén wavepacket collisions, the fundamental building block of astrophysical plasma turbulence, a key result of this study.

Because the simulation domain is periodic in the $z$ direction, the Alfvén wavepackets undergo a second collision at the boundary of the domain $z=\pm L_{z}/2$ at $t/T_{c}=1.5$, shown in Fig. 2(d), followed by a third collision at $z=0$ at $t/T_{c}=2.5$, shown in Fig. 2(f). Below we explore in more detail the cascade of energy to smaller perpendicular scales, the development of current sheets, the evolution of the energy in perpendicular Fourier modes, and key properties of localized Alfvén wavepacket collisions.

3.1 The Perpendicular Cascade of Energy

To explore the nonlinear cascade of smaller perpendicular scales in this strong Alfvén wavepacket collision simulation, we plot in Fig. 4 the perpendicular magnetic energy of fluctuations in perpendicular Fourier space integrated over $z$, $E_{B_{\perp}(k_{x},k_{y})}\equiv\int_{-L_{z}/2}^{L_{z}/2}dz|\delta B_{\perp}(k_{x},k_{y})|^{2}/8\pi$ in arbitrary units. In Fig. 4(a) is the perpendicular magnetic energy of the initial Alfvén wavepackets at $t=0$, showing all of the energy is contained within the three perpendicular Fourier modes $(k_{x}\rho_{i},k_{y}\rho_{i})=(1,0)$, $(0,1)$, and $(-1,0)$. Note that AstroGK uses a reality condition imposed on the complex Fourier coefficients such that the magnetic field can be described using only Fourier modes in the upper half-plane $k_{y}\geq 0$. The reality condition requires $\delta\hat{B}(k_{x},k_{y})=\delta\hat{B}^{*}(-k_{x},-k_{y})$, so if there is power in the $(1,0)$ mode, there must be equivalent power in the $(-1,0)$ mode, as seen in Fig. 4(a).

After the first collision at $t/T_{c}=1$, Fig. 4(b) shows definitively that the strong nonlinear interactions have transferred significant energy to many other modes in the perpendicular Fourier plane $(k_{x},k_{y})$. Note that the perpendicular wavenumber is given by $k_{\perp}=\sqrt{k_{x}+k_{y}}$, so that modes further away from the origin represent smaller scale fluctuations in the perpendicular plane. The plot after the second collision at $t/T_{c}=2$ in Fig. 4(c) shows that successive collisions continue to facilitate the perpendicular cascade of energy. These results demonstrate definitively that the finding from the idealized periodic case—that Alfvén wave collisions mediate the nonlinear transfer of energy to smaller perpendicular scales—indeed persists in the more realistic case of localized Alfvén wavepacket collisions, answering the first key question posed in the introduction.

3.2 Current Sheet Development

For plane Alfvén wave collisions in the idealized periodic case, plots of the parallel current density $j_{z}$ in the $(x,y)$ plane perpendicular to the equilibrium magnetic field demonstrate that strong Alfvén wave collisions self-consistently generate current sheets (Howes, 2016). These current sheets extended the full parallel length of the original Alfvén waves, with widths in the perpendicular plane of approximately the perpendicular wavelength of the original interacting Alfvén waves, but with a much smaller thickness in the perpendicular plane. Here we use plots of the $(x,y)$ plane perpendicular to the equilibrium magnetic field to determine whether this current sheet formation persists in the more realistic case of strong collisions between localized Alfvén wavepackets.

First, we examine the development of the parallel current density $j_{z}$ in the $(x,y)$ plane during the course of each collision. In Fig. 5, we plot the normalized parallel current density $j_{z}/j_{0}$ at $z=0$ (a) early in the first collision at $t/T_{c}=0.3$ and (b) later in the same collision at $t/T_{c}=0.5$, where we remind the reader that the midpoint of each collision occurs at $t/T_{c}=0.5,1.5,2.5$. The nearly circular $j_{z}$ pattern in Fig. 5(a) is largely due to the linear superposition of the current of each of the initial Alfvén wavepackets (the top row of Fig. 6 shows the current patterns of each of the initial Alfvén wavepackets). The nonlinear interaction of the Alfvén wavepackets during the first collision leads to the thinning of these initial circular current patterns into a more sheet-like morphology in Fig. 5(b). The second row of Fig. 5 shows the evolution of the second collision in the plane $z=L_{z}/2$ (at the parallel boundary of the periodic domain) at (c) $t/T_{c}=1.3$ and (d) $t/T_{c}=1.5$. Here we see the current sheets becoming yet more elongated and intense through the course of the second collision. The third row shows plots of $j_{z}/j_{0}$ during the third collision at $z=0$, with further thinning and intensification of the current sheets. Note that the three-dimensional isocontours of $j_{z}$ in Fig. 2 shows that the parallel extent of the current sheets is the same as the parallel length of the original wavepackets.

These results show clearly that strong collisions of localized Alfvén wavepackets indeed self-consistently generate intermittent current sheets, definitively answering the second question in the introduction.

Another important question, however, is whether the current sheets that develop during the wavepacket collisions persist within each wavepacket after the collision. In Fig. 6, we show the initial pattern at $t=0$ of the normalized parallel current density $j_{z}/j_{0}$ for (a) the unipolar wave at $z=-L_{z}/4$ with only variation in the $x$ direction and (b) the dipolar wave at $z=+L_{z}/4$ with only variation in the $y$ direction. After the first collision, we show $j_{z}$ for (c) the unipolar wave at $z=+L_{z}/4$ and $t/T_{c}=0.98$ and (d) the dipolar wave at $z=-L_{z}/4$ and $t/T_{c}=0.96$. At this time when the wavepackets are no longer overlapping in $z$, one can clearly see the distortion of current density pattern, due to the previous collision, has lead to a thinning and intensification of the current density into a more sheet-like morphology that persists after the wavepackets have separated. After the second collision, we plot (e) the unipolar wave at $z=-L_{z}/4$ and $t/T_{c}=2.1$ and (f) the dipolar wave at $z=+L_{z}/4$ and $t/T_{c}=2.0$. Here we can see that the current sheets that develop as a consequence of the strong Alfvén wavepacket collisions indeed persist within the wavepackets after the collision is over, showing that these interactions in the more realistic case may explain the ubiquitous observations of intermittent current sheets in plasma turbulence.

3.3 Evolution of Energy

A simple physical interpretation of the evolution of the energy in this localized Alfvén wavepacket collision can be developed using the equations of incompressible MHD for guidance. These equations, written in terms of Elsasser variables, take the form

[TABLE]

and $\nabla\cdot\mathbf{z}^{\pm}=0$. Here $\mathbf{v}_{A}=\mathbf{B}_{0}/\sqrt{4\pi(n_{0i}m_{i}+n_{0e}m_{e})}$ is the Alfvén velocity due to the equilibrium field \mathbf{B}_{0}=B_{0}\mbox{\hat{\mathbf{z}}} where $\mathbf{B}=\mathbf{B}_{0}+\delta\mathbf{B}$, $P$ is total pressure (thermal plus magnetic), and $n_{0i}m_{i}+n_{0e}m_{e}$ is mass density. Recall that the Elsasser variables are defined by $\mathbf{z}^{\pm}=\mathbf{u}\pm\delta\mathbf{B}/\sqrt{4\pi(n_{0i}m_{i}+n_{0e}m_{e})}$, representing waves that propagate up or down the mean magnetic field. The nonlinear term, $\mathbf{z}^{\mp}\cdot\nabla\mathbf{z}^{\pm}$, governs the nonlinear interactions of the counterpropagating Alfvén wave collisions.

In Fig. 7, we plot the temporal evolution of the magnetic energy $E_{B_{\perp}(k_{x},k_{y})}$ in some of the lowest perpendicular Fourier modes $(k_{x},k_{y})$, illustrating some of the key properties of Alfvén wavepacket collisions. The first main point is that there is no nonlinear transfer of energy among Fourier modes until the two wavepackets begin to overlap in space along $z$ (first vertical dashed line). As expected from theoretical considerations of incompressible MHD turbulence (Iroshnikov, 1963; Kraichnan, 1965; Sridhar & Goldreich, 1994; Montgomery & Matthaeus, 1995; Ng & Bhattacharjee, 1996; Galtier et al., 2000; Howes & Nielson, 2013), the nonlinear interaction is zero unless both $z^{-}$ (the unipolar wavepacket) and $z^{+}$ (the dipolar wavepacket) are both non-zero at the same point in space, as can be seen by inspection of the nonlinear term $\mathbf{z}^{\mp}\cdot\nabla\mathbf{z}^{\pm}$. Between collisions, when the wavepackets have no overlap in $z$, there is no nonlinear transfer of energy among different $(k_{x},k_{y})$ modes, as clearly seen in Fig. 7. The second main point is that the perpendicular Fourier modes $(k_{x}L_{\perp},k_{y}L_{\perp})$ of the initial Alfvén wavepackets, $(1,0)$ and $(0,1)$ generally lose energy to nonlinearly generated perpendicular Fourier modes during each collision. This is the transfer of energy to smaller perpendicular scales, illustrated by the plots of perpendicular magnetic energy in $(k_{x},k_{y})$ Fourier space in Fig. 4. A more quantitative examination contrasting the nonlinear transfer of energy between the idealized plane Alfvén wave collisions in a periodic geometry and localized Alfvén wavepackets collisions, in both the weakly and strongly nonlinear limits, will be presented in a subsequent paper (Verniero & Howes, 2017).

3.4 Wavepacket Collisions

The three-dimensional visualization of this localized Alfvén wavepacket collision in Fig. 2 presents a concise overview of many of the properties examined above. The series of panels (a) through (c) show clearly the distortion of the original wavepackets that arises during the collision and persists after the wavepackets have separated again. This is a physical visualization of the nonlinear cascade of energy to smaller scales, the key characteristic of turbulence in space and astrophysical plasmas. Although a little difficult to see clearly in the 3D projection in Fig. 2, the nonlinear evolution of this strong Alfvén wavepacket collision also leads to the self-consistent generation of a sheet-like morphology for the regions of intense parallel current density $j_{z}$. This thinning and intensifying of the current into sheets is more clearly seen in the perpendicular cross section at $z=L_{z}/2$ (the right-hand $z$ boundary of the simulation domain) in Fig. 2(d); the same cross section at $z=L_{z}/2$ is also shown in Fig. 5(d).

In addition to the isocontours of the normalized parallel current density $j_{z}/j_{0}$, we also plot the paths of a sample of magnetic field lines as they traverse the domain. Beginning at points on a $5\times 5$ grid at $z=-L_{z}/2$, we trace the field lines in the $+z$ direction. Although the perpendicular magnetic field fluctuation associated with each Alfvén wavepacket is small compared to the equilibrium magnetic field, $\delta B_{\perp}\ll B_{0}$, one can still see the distortion of the field lines as they pass through each wavepacket. In the space between wavepackets, both before and after each collision, the magnetic field lines are straight and uniform, further illustrating the point that nonlinear evolution ceases when the wavepackets are separated, even in a gyrokinetic simulation where dispersive and kinetic effects are resolved.

It is worthwhile also noting the important point that, even after each strong Alfvén wavepacket collision, the individual wavepackets continue to propagate along the equilibrium magnetic field and they remain localized in $z$. There is a very small spreading of the wavepacket, hardly noticeable in Fig. 2, due to the fact that Alfvén waves become dispersive at $k_{\perp}\rho_{i}\gtrsim 1$, with an increasing parallel group velocity. Therefore, some of the higher $k_{\perp}$ modes that are nonlinearly generated by the strong Alfvén wavepacket collision in this gyrokinetic simulation—a numerical approach that resolves these dispersive kinetic effects at $k_{\perp}\rho_{i}\gtrsim 1$—will propagate a little faster than Alfvén waves in the MHD limit $k_{\perp}\rho_{i}\ll 1$ that have a parallel group velocity at the Alfvén speed $v_{g\parallel}=v_{A}$, causing a slight spreading out of the localized wavepacket.

3.5 Relation to Similar Wavepacket Collision Simulations

A recent study by Pezzi et al. (2017) has used numerical simulations to tackle what they term the “Parker-Moffatt problem.” The basic question is inspired by the properties of incompressible MHD equations given in (1): do two initially separated, localized Alfvén wavepackets interact nonlinearly only during the time that they overlap, and then cease evolving after they have passed through each other and become separated again? The simulation results shown in Fig. 2, and examined in further detail in the subsections above, definitively answer this question—yes.

The results of the gyrokinetic simulation of strong Alfvén wavepacket collisions presented here, however, stand in stark contrast to the results of the study by Pezzi et al. (2017), which finds that after the collision, in both their compressible MHD and Hybrid Vlasov-Maxwell simulations, the wavepackets do not separate cleanly and continue to evolve, as illustrated in Figure 1 of their paper. What is the cause for the dramatic contradiction in these results?

We argue here that the limitation of the Pezzi et al. (2017) simulations to two spatial dimensions is the root of these striking differences. Due to the large computational cost of running the Hybrid Vlasov-Maxwell (HVM) code (Valentini et al., 2007) in three spatial dimensions, the authors chose to approximate the 3D problem by using an oblique 2D approach (Howes, 2014), where spatial variation is allowed in two dimensions in the $(x,y)$ plane, but the equilibrium magnetic field is slightly tilted away from the normal to the plane, \mathbf{B}_{0}=B_{0}(\sin\theta\mbox{\hat{\mathbf{x}}}+\cos\theta\mbox{\hat{\mathbf{z}}}) with $\theta=6^{\circ}$, yielding a small component of the equilibrium magnetic field in the $x$ direction. Since plane Alfvén waves only propagate if there is a nonzero component of the wavevector parallel to the equilibrium magnetic field, $k_{\parallel}\neq 0$, this trick enables the $x$ component of a wavevector to have a small component parallel to $\mathbf{B}_{0}$, enabling propagating Alfvén waves to be simulated.

The limitation to two spatial dimensions has two important effects: (i) it precludes the possibility of initializing two perpendicularly polarized, counterpropagating wavepackets (as initialized here), which is the essential set up of the Parker-Moffatt problem to be studied; and (ii) it significantly limits the nonlinear couplings that can arise, constraining the results in an unphysical way (Howes, 2014).

Although Pezzi et al. (2017) state that they set up two Alfvénic wavepackets, several pieces of evidence suggest that their initial wavepackets are indeed not Alfvénic . First, they initialize the fluctuations by imposing the constraint that fluctuations are strictly transverse to the equilibrium magnetic field $\mathbf{B}_{0}\cdot\delta\mathbf{B}=0$. Although this is a necessary condition for an Alfvén wave, it is not sufficient: the polarization of the magnetic field fluctuation for the Alfvén wave with wavevector $\mathbf{k}$ must also be oriented in the direction given by $\mathbf{B}_{0}\times\mathbf{k}$ (Maron & Goldreich, 2001; Howes & Nielson, 2013; Howes, 2014). Any component of the magnetic field fluctuation in the other direction of the plane given by $\mathbf{B}_{0}\cdot\delta\mathbf{B}=0$ represents not an Alfvén wave, but a different mode: in incompressible MHD, it would belong to the pseudo-Alfvén wave (Maron & Goldreich, 2001; Howes & Nielson, 2013; Howes, 2014), and in compressible MHD, it would represent some combination of the fast and slow magnetosonic modes (Cho & Lazarian, 2003). That their initialization in fact excites some magnetosonic modes is supported by two results in their paper. First, their initialization excites density fluctuations even before the wavepackets collide, a feature they attribute to magnetic pressure fluctuations due to finite-amplitude magnetic fluctuations, but that could simply be because their initialization directly included these linear magnetosonic modes. Second, even before the waves collide in the MHD simulation, shown in the upper left panel of their Figure 1, one can see small disturbances well ahead of and behind each wavepacket, suggesting that magnetosonic modes, which have a different group velocity from the Alfvén mode, are dispersively spreading out the wavepacket, a linear effect.

The second impact of the limitation to two spatial dimensions is that the nonlinear couplings possible in the oblique 2D configuration are unphysically constrained relative to the three dimensional case, as explained in detail in Howes (2014). This constraint can artificially enhance the nonlinear coupling to magnetosonic modes relative to the coupling to Alfvénic modes, dramatically altering the nonlinear evolution.

Our perspective is that any investigation of the so-called Parker-Moffatt problem is most cleanly performed in three spatial dimensions. Our results show definitively that when two localized Alfvén wavepackets collide, the interaction occurs only while they overlap, ceasing when the wavepackets separate cleanly, leaving distorted wavepackets that have essentially the same finite extent along the equilibrium magnetic field. Before concluding, it is also worthwhile emphasizing that our results where obtained with a gyrokinetic simulation that resolves all of the low-frequency—not including cyclotron and fast magnetosonic physics (Howes et al., 2006)—compressible and kinetic effects of a weakly collisional plasma.

4 Conclusion

Our results here show that many of the fundamental properties of Alfvén wave collisions, originally characterized analytically and numerically in the idealized case of a periodic Alfvén wave collision, persist under the more realistic conditions of a localized Alfvén wavepacket collision. Specifically, we have demonstrated that strong localized Alfvén wavepacket collisions are effective in mediating the nonlinear cascade of energy to small perpendicular scales, as demonstrated in Fig. 4. Furthermore, we have shown that strong localized Alfvén wavepacket collisions also self-consistently generate current sheets that persist even between collisions, as shown in Figures 5 and 6, confirming a robust mechanism to explain the ubiquitous current sheets observed in numerical simulations of plasma turbulence (Wan et al., 2012; Karimabadi et al., 2013; TenBarge & Howes, 2013; Wu et al., 2013; Zhdankin et al., 2013) and inferred from measurements of solar wind turbulence (Osman et al., 2011; Borovsky & Denton, 2011; Osman et al., 2012; Perri et al., 2012; Wang et al., 2013; Wu et al., 2013; Osman et al., 2014).

The overall evolution of localized Alfvén wavepacket collisions involves strong nonlinear interactions occurring while the wavepackets overlap, followed by a clean separation of the wavepackets with straight uniform magnetic fields in between and the cessation of nonlinear evolution, as visualized in Fig. 2. The wavepackets remain localized along the equilibrium magnetic field before and after each wavepacket collision. It is important to emphasize that these characteristics are predicted based on the properties of the incompressible MHD equations, but are found even in the gyrokinetic simulation presented here which resolves dispersive and kinetic effects beyond the reach of the MHD theory. That these important properties of Alfvén wave collisions persist in the realistic localized wavepacket case further supports the contention that Alfvén wave collisions represent the fundamental building block of astrophysical plasma turbulence.

Previous analytical and numerical work in the limit of weakly nonlinear, periodic Alfvén wave collisions, (Howes & Nielson, 2013; Nielson et al., 2013) has shown that a nonlinearly generated mode, which is not a solution of the linear dispersion relation, serves to mediate the energy transfer to small perpendicular scales. Future work will address the question of whether this nonlinearly generated mode still plays a key role in the more realistic case of localized Alfvén wavepacket collisions (Verniero & Howes, 2017).

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. 1048957, NSF PHY-10033446, NSF CAREER AGS-1054061, and DOE DE-SC0014599. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. We especially thank one of their staff members, Dave Pugmire, for invaluable assistance via email. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575, through NSF XSEDE Award PHY090084.

Appendix A Windowing Function for Localized Alfvén Wavepacket Collisions

A.1 Parallel Dependence of Alfvén Eigenfunctions

Initializing Alfvén wavepackets localized along the equilibrium magnetic field using the AstroGK code (Numata et al., 2010) takes advantage of the remarkable feature of gyrokinetics that the linear collisionless gyrokinetic dispersion relation, when suitably normalized, is independent of the parallel wavenumber of fluctuations, as explained below.

For a fully ionized proton-electron plasma with Maxwellian equilibrium velocity distributions and zero net drift velocities, the complex frequency that is the solution of the Vlasov-Maxwell linear dispersion relation (Stix, 1992; Quataert, 1998) can be expressed as a function of five dimensionless parameters, $\omega/\Omega_{i}=\tilde{\omega}_{VM}(k_{\parallel}\rho_{i},k_{\perp}\rho_{i},\beta_{i},T_{i}/T_{e},v_{ti}/c)$ (Howes et al., 2006). In the gyrokinetic limit of non-relativistic ($v_{ti}/c\ll 1$), low-frequency ($\omega/\Omega_{i}\ll 1$), anisotropic ($k_{\parallel}/k_{\perp}\ll 1$) fluctuations, the complex eigenfrequency of the linear collisionless gyrokinetic dispersion relation has just three dimensionless parameters when normalized by $k_{\parallel}$, $\omega/(k_{\parallel}v_{A})=\overline{\omega}_{GK}(k_{\perp}\rho_{i},\beta_{i},T_{i}/T_{e})$ (Howes et al., 2006). The fact that the suitably normalized complex eigenfrequency and complex Fourier coefficients of the eigenfunction in gyrokinetics is independent of the parallel wavenumber $k_{\parallel}$ means that the phase and amplitude relationships among the different components of the eigenfunction are independent of the parallel wavenumber or coordinate.

Consider, for example, the procedure for initializing an Alfvén wavepacket localized along the direction parallel to the equilibrium magnetic field \mathbf{B}_{0}=B_{0}\mbox{\hat{\mathbf{z}}}. For a perpendicular Fourier component of that wavepacket with a wavevector $(k_{x},k_{y})$, we specify the complex perpendicular Fourier coefficient of the $y$ component of the perturbed magnetic field by $\delta\hat{B}_{y}(k_{x},k_{y},z)$. To obtain the complex perpendicular Fourier coefficient of any other component of the eigenfunction, for example the parallel magnetic field perturbation, $\delta\hat{B}_{z}(k_{x},k_{y},z)$, we simply multiply the ratio of the suitably normalized linear eigenfunction components,

[TABLE]

The key property of gyrokinetics that we exploit is the fact that the ratio in brackets, determined from the linear eigenfunction, is independent of the parallel wavenumber $k_{z}$ and parallel coordinate $z$. For simplicity, we suppress the dependence in the plasma parameters $(\beta_{i},T_{i}/T_{e})$ below, as they are constants for a given calculation.

To demonstrate that the linear eigenfunction derived for a given perpendicular wavevector $(k_{x},k_{y})$ can be used to construct a wavepacket of arbitrary functional form in $z$, we note that the Fourier transform in $z$ of a component of the fluctuating magnetic field can be expressed as

[TABLE]

where the index $j$ denotes any component $x$, $y$, or $z$. By substituting (3) into the left side of (2) for $\delta\hat{B}_{z}(k_{x},k_{y},z)$ and right side of (2) for $\delta\hat{B}_{y}(k_{x},k_{y},z)$, we obtain

[TABLE]

where we have exploited the fact that the factor in brackets is independent of $k_{z}$ to bring it inside of the summation on the right-hand side. The orthogonality of basis functions $e^{ik_{z}z}$ means that the corresponding terms in the sum must be equal, yielding the result

[TABLE]

Therefore, the phase and amplitude relationships between different components of the eigenfunction are the same whether the parallel dependence is expressed in coordinate space $z$ or Fourier space $k_{z}$.

The procedure for initializing a localized wavepacket in AstroGK with an arbitrary functional form in the parallel direction $z$ therefore follows: (i) the form of the wavepacket in the parallel direction $z$ for the perpendicular Fourier component of the parallel vector potential $\delta\hat{A}_{\parallel}(k_{x},k_{y},z)$ is specified; (ii) the linear collisionless gyrokinetic dispersion relation $\overline{\omega}_{GK}(k_{\perp}\rho_{i},\beta_{i},T_{i}/T_{e})$ is solved for the other suitably normalized electromagnetic field potentials $\hat{\overline{\phi}}(k_{x},k_{y})$ and $\delta\hat{\overline{B}}_{\parallel}(k_{x},k_{y})$ as well as the perturbed ion and electron gyrokinetic distribution functions, $\hat{\overline{h}}_{i}(v_{\parallel},v_{\perp},k_{x},k_{y})$ and $\hat{\overline{h}}_{e}(v_{\parallel},v_{\perp},k_{x},k_{y})$; and (iii) the ratio of the suitably normalized eigenfunctions are used to compute the form of the wavepacket in the parallel direction $z$ for the perpendicular Fourier component of all the other components of the eigenfunction, for example $\delta\hat{B}_{\parallel}(k_{x},k_{y},z)$. This completes the specification of a wavepacket of arbitrary parallel functional dependence for the chosen linear gyrokinetic wave mode.

A.2 Specification of Windowing Function

To localize the Alfvén wavepacket, we use a windowing function $w(z)$ for the complex Fourier wave amplitude of each perpendicular Fourier mode. The windowing function in $z$ for a given perpendicular Fourier mode is

[TABLE]

where the three window function parameters are the center of the window $z_{0}$, the width of the window $\Delta_{z}$ and the power $p$ of the exponential. Note that the second and third terms are needed to ensure periodicity in $z$ when the window does not fall to zero at the $z$ limits of the domain (it is important, since only two terms are included, that the window does fall to zero by at least an additional distance $L_{z}$ beyond each $z$ boundary).

The sinusoidal variation in $z$ of each perpendicular Fourier mode is given by

[TABLE]

where $k_{z}$ is the wavenumber along the equilibrium magnetic field and $\delta$ is a phase-shift for that mode. The total initial waveform along $z$ is then given by the product of the sinusoidal mode and the windowing function, $w(z)f(z)$.

For the localized Alfvén wavepacket collision simulation presented in this paper, the parameters for the dipolar wave (which propagates in the $-z$ direction) are $k_{z}a_{0}=-3$, $\delta=0$ $z_{0}=\pi/2a_{0}=L_{\parallel}/4$, $\Delta_{z}=1/2a_{0}$, and the default power $p=2$. For the unipolar wave (which propagates in the $+z$ direction), the parameters are $k_{z}a_{0}=+1$, $\delta=\pi/2$ $z_{0}=-\pi/2a_{0}=-L_{\parallel}/4$, $\Delta_{z}=1/2a_{0}$, and the default power $p=2$. Note that the non-zero phase $\delta$ for this unipolar wavepacket means it has an initial form $-\sin(k_{z}z)$.

Note that the aspect ratio of the characteristic perpendicular length scale to the characteristic parallel length scale $a_{0}/\rho_{i}\equiv\epsilon\ll 1$ is the small gyrokinetic expansion parameter. The parallel domain scale is given by $k_{\parallel}a_{0}=1$, where the parallel domain length $L_{\parallel}=2\pi/k_{\parallel}$, so the domain length can also be expressed as $L_{\parallel}=2\pi a_{0}$. All parallel length scales in AstroGK are normalized by $a_{0}$, whereas all perpendicular length scales are normalized by $\rho_{i}$, or equivalently dimensionless perpendicular wavenumber $k_{\perp}\rho_{i}$.