Colliding Alfv\'enic wave packets in MHD, Hall and kinetic simulations

TL;DR

This paper investigates the collision of Alfvénic wave packets using MHD, Hall, and kinetic simulations, revealing dispersive and kinetic effects, non-Maxwellian features, and solar wind-like proton distributions.

Contribution

It extends the Parker-Moffatt problem analysis by incorporating Hall and kinetic models, providing detailed insights into dispersive and kinetic phenomena during wave packet collisions.

Findings

Regions with strong temperature anisotropy form during collisions

Proton distribution functions develop beams along magnetic field lines

Dispersive and kinetic effects are significant in wave interactions

Abstract

The analysis of the Parker-Moffatt problem, recently revisited in Pezzi et al. (2016), is here extended by including the Hall magnetohydrodynamics and two hybrid kinetic Vlasov-Maxwell numerical models. The presence of dispersive and kinetic features is studied in detail and a comparison between the two kinetic codes is also reported. Focus on the presence of non-Maxwellian signatures shows that - during the collision - regions characterized by strong temperature anisotropy are recovered and the proton distribution function displays a beam along the direction of the magnetic field, similar to some recent observations of the solar wind.

Equations

$${\partial }_t\rho +\nabla \cdot (\rho \mathbf{u})=0$$

$${\partial }_t\mathbf{u}+(\mathbf{u}\cdot \nabla )\mathbf{u}=-\frac{\tilde{\beta }}{2\rho }\nabla (\rho T)+\frac{1}{\rho }\left[(\nabla \times \mathbf{B})\times \mathbf{B}\right]$$

$${\partial }_t\mathbf{B}=\nabla \times \left[\mathbf{u}\times \mathbf{B}-\frac{\widetilde{ϵ}}{\rho }(\nabla \times \mathbf{B})\times \mathbf{B}\right]$$

$${\partial }_{t}T+(\mathbf{u}\cdot \nabla )T+(\gamma -1)T(\nabla \cdot \mathbf{u})=0$$

$${\partial }_{t}f+\mathbf{v}\cdot \nabla f+\frac{1}{\widetilde{ϵ}}$$

$$\frac{\partial \mathbf{B}}{\partial t}=-\nabla \times \mathbf{E};$$

$$ϵ(x,y,t)=\frac{1}{n}\sqrt{\int {\left[f(\mathbf{x},\mathbf{v},t)-f_M(\mathbf{x},\mathbf{v},t)\right]}^{2}d\mathbf{v}}$$

Full text

{PACS}

Authors should not enter PACS codes directly on the manuscript, as these must be chosen during the online submission process and will then be added during the typesetting process (see http://www.aip.org/pacs/ for the full list of PACS codes)

1 Introduction

The interaction of two oppositely propagating Alfvénic wave packets has been studied for more than half a century. This interaction has been proposed as an elementary step in the analysis of magnetohydrodynamics (MHD) turbulence (Elsässer, 1950; Iroshnikov, 1964; Kraichnan, 1965; Dobrowolny et al., 1980a, b; Velli et al., 1989; Sridhar & Goldreich, 1994; Goldreich & Sridhar, 1995; Ng & Bhattacharjee, 1996; Matthaeus et al., 1999; Galtier et al., 2000; Verdini et al., 2009; Howes & Nielson, 2013; Nielson et al., 2013). Indeed, in the framework of ideal incompressible magnetohydrodynamics (MHD), large amplitude perturbations in which the magnetic ${\bf b}$ and bulk velocity ${\bf u}$ fluctuations are either perfectly correlated, or perfectly anti-correlated, are solutions of the governing equations. To induce nonlinear couplings among the fluctuations, and therefore to excite turbulence, it is necessary to simultaneously consider magnetic fluctuations ${\bf b}$ and velocity fluctuations ${\bf u}$ that have an arbitrary sense of correlation. This may be accomplished by superposing the two senses of correlation, in Alfvén units, ${\bf u}=+{\bf b}$ and ${\bf u}=-{\bf b}$.

Based on these considerations, Moffatt (1978) and Parker (1979) analyzed the collision of large-amplitude Alfvén wave packets in the framework of the incompressible MHD, observing that their interaction is limited to the time interval in which they overlap. During this temporal window wave packets can transfer energy and modify their spatial structure; however, after the collision, packets return to undisturbed propagation without further interactions. The Parker–Moffatt problem has been recently revisited in (Pezzi et al., 2016) (hereafter, Paper I) with the motivation to extend its description to the realm of the kinetic plasmas. In fact, the scenario described by Parker & Moffatt is potentially applicable to astrophysical plasmas such as solar wind (Belcher & Davis, 1971; Bruno et al., 1985; Verdini et al., 2009) or solar corona (Matthaeus et al., 1999; Tomczyk et al., 2007), where Alfvénic perturbations represent one of the main components of fluctuations. However, since such systems often exhibit compressive activity as well as dispersion and kinetic signatures (Alexandrova et al., 2008; Bruno & Carbone, 2013; Sahraoui et al., 2007; Gary et al., 2010; Marsch, 2006; Valentini et al., 2011; Servidio et al., 2012; Valentini et al., 2014; Servidio et al., 2015; He et al., 2015; Roberts et al., 2016; Lion et al., 2016; Perrone et al., 2016), it is of considerable interest to include these features in the analysis of the Parker & Moffatt problem.

In particular, in Paper I, it has been found that during the wave packets interaction, as prescribed by Parker & Moffatt, nonlinear coupling processes cause the magnetic energy spectra to evolve towards isotropy, while energy transfers towards smaller spatial scales. Moreover, the new ingredients introduced with the kinetic simulation (Hall and kinetic effects) play a significant role and several features of the evolution in the Vlasov case differ with respect to the MHD evolution. Here we extend that study to discern the role of dispersive and genuinely kinetic effects, supplementing the previously considered MHD and Vlasov simulations, by introducing also an Hall MHD simulation. Moreover, we also examine this basic problem by means of a hybrid Particle-in-Cell simulation (HPIC), which allows comparison of two different numerical approaches (HVM and HPIC), which refer to the same physical model. We may anticipate that, in the HPIC case, the system dynamics at small scales is affected by the presence of particles thermal noise and only the features related to large spatial scales are properly recovered during the evolution of the two wave packets. Based on this consideration, we employ mainly the HVM simulation to highlight the presence of kinetic effects during the wave packets interaction. In particular, during the collision of the wave packets, the proton velocity distribution function (VDF) exhibits a beam along the background magnetic field direction, similar to some solar wind observations (He et al., 2015). We note that the present paper compares results from four different models in the context of a single physical problem, and is therefore also a contribution in the spirit of the “Turbulence Dissipation Challenge” that has been recently discussed in the space plasma community Parashar et al. (2015b).

The paper is organized as follows: in Sec. 2 the theoretical model and the numerical codes are presented. In Sec. 3, we compare the several simulations by focusing on the description of some fluid-like diagnostics. Sec. 4, examines kinetic signatures in the HVM simulation. Finally we conclude in Sec. 5 by summarizing our results.

2 Theoretical models and numerical approaches

As discussed above, here we approach the problem concerning the interaction of two Alfvénic wave packets by means of fluid and hybrid numerical simulations. For problems such as this, the system dimensionality is fundamental: in fact, a proper description should consider a three-dimensional physical space (i.e. three-dimensional wave vectors), where both parallel and perpendicular cascades are taken into account (Howes, 2015; Parashar et al., 2015a, 2016). However, dynamical range of the spatial scales (wave numbers) represented in the model is equally important to capture nonlinear couplings during the wave packet interaction. Furthermore, performing a kinetic HVM simulation which contemporaneously includes a full $3D$–$3V$ phase space and while also retaining a good spatial resolution is too demanding for the present High Performance Computing capability. Given that numerous runs are required to complete a study such as the present one, a fully 3D approach would be prohibitive. Therefore we adopt a $2.5D$ physical space, where vectorial fields are three-dimensional but their variations depend only on two spatial coordinates ($x$ and $y$). It is worth noting that $2.5D$ captures the qualitative nature of many processes very well even though there might be some quantitative differences for some processes (Karimabadi et al., 2013; Wan et al., 2015; Li et al., 2016).

The fluid models here considered are MHD and Hall MHD, whose dimensionless equations are:

[TABLE]

In Eqs. (1)–(4) spatial coordinates ${\bf x}=(x,y)$ and time $t$ are respectively normalized to $\tilde{L}$ and $\tilde{t}_{A}=\tilde{L}/\tilde{c}_{A}$. The magnetic field ${\bf B}={\bf B}_{0}+{\bf b}$ is scaled to the typical magnetic field ${\tilde{B}}$, while mass density $\rho$, fluid velocity ${\bf u}$, temperature $T$ and pressure $p=\rho T$ are scaled to typical values ${\tilde{\rho}}$, ${\tilde{c}}_{A}={\tilde{B}}/(4\pi{\tilde{\rho}})^{1/2}$, ${\tilde{T}}$ and ${\tilde{p}}=2\kappa_{B}{\tilde{\rho}}{\tilde{T}}/m_{p}$ (being $\kappa_{B}$ the Boltzmann constant and $m_{p}$ the proton mass), respectively. Moreover, ${\tilde{\beta}}={\tilde{p}}/({\tilde{B}}^{2}/8\pi)$ is a typical value for the kinetic to magnetic pressure ratio; $\gamma=5/3$ is the adiabatic index and ${\tilde{\epsilon}}={\tilde{d}}_{p}/{\tilde{L}}$ (being ${\tilde{d}}_{p}={\tilde{c}}_{A}/{\tilde{\Omega}}_{cp}$ the proton skin depth) is the Hall parameter, which is set to zero in the pure MHD case. Details about the numerical algorithm can be found in Vásconez et al. (2015); Pucci et al. (2016).

On the other hand, hybrid Vlasov-Maxwell simulations have been performed by using two different numerical codes: an Eulerian hybrid Vlasov-Maxwell code (HVM) (Valentini et al., 2007) and a hybrid Particle-in-cell (HPIC) code (Parashar et al., 2009). For both cases protons are described by a kinetic equation,while electrons are a Maxwellian, isothermal fluid. In the Vlasov model, an Eulerian representation of the Vlasov equation for protons is numerically integrated. In PIC method, the distribution function is Monte-Carlo discretized and the Newton-Lorentz equations are updated for the “macro-particles”. Electromagnetic fields, charge density and current density are computed on a grid (Birdsall & Langdon, 2004; Dawson, 1983).

Dimensionless HVM equations are:

[TABLE]

where $f=f({\bf x},{\bf v},t)$ is the proton distribution function. In Eqs. (6)–(7), velocities ${\bf v}$ are scaled to the Alfvén speed ${\tilde{c}_{A}}$, while the proton number density $n=\int f\,d^{3}v$, the proton bulk velocity ${\bf u}=n^{-1}\int{\bf v}f\,d^{3}v$ and the proton pressure tensor ${\Pi}_{ij}=n^{-1}\int(v-u)_{i}\,(v-u)_{j}f\,d^{3}v$, obtained as moments of the distribution function, are normalized to $\tilde{n}=\tilde{\rho}/m_{p}$, $\tilde{c}_{A}$ and ${\tilde{p}}$, respectively. The electric field $\bf{E}$, the current density ${\bf j}={\bf\nabla}\times{\bf B}$ and the electron pressure $P_{e}$ are scaled to $\tilde{E}=(\tilde{c}_{A}\tilde{B})/c$, $\tilde{j}=c\tilde{B}/(4\pi\tilde{L})$ and $\tilde{p}$, respectively. Electron inertia effects have been considered in Ohm’s law to prevent numerical instabilities (being $m_{e}/m_{p}=0.01$, where $m_{e}$ is the electron mass), while no external resistivity $\eta$ is introduced. A detailed description of the HVM algorithm can be found in Valentini et al. (2007). On the other hand, the hybrid PIC run has been performed using the P3D hybrid code (Zeiler et al., 2002) and all the numerical and physical parameters are the same as the HVM run. The code has been extensively used for reconnection and turbulence (e.g. Parashar et al., 2009; Malakit et al., 2009).

In both classes of performed simulation (fluid and kinetic), the spatial domain $D(x,y)=[0,8\pi]\times[0,2\pi]$ is discretized with $(N_{x},N_{y})=(1024,256)$ in such a way that $\Delta x=\Delta y$ and spatial boundary conditions are periodic. For the HVM run, the velocity space is discretized with an uniform grid with $51$ points in each direction, in the region $v_{i}=[-v_{max},v_{max}]$ (being $v_{max}=2.5{\tilde{c}_{A}}$) and velocity domain boundary conditions assume $f=0$ for $|v_{i}|>v_{max}$ ($i=x,y,z$); while, in the HPIC case, the number of particles per cell is $400$. Moreover $\beta_{p}=2v_{th,p}^{2}/{\tilde{c}_{A}}^{2}={\tilde{\beta}}/2=0.5$ (i.e. $v_{max}=5v_{th,p}$), ${\bf u_{e}}={\bf u}-{\tilde{\epsilon}}{\bf j}/n$, ${\tilde{\epsilon}}=9.8\times 10^{-2}$, $k_{d_{p}}={\tilde{\epsilon}}^{-1}\simeq 10$ and $k_{d_{e}}=\sqrt{m_{p}/m_{e}}\times{\tilde{\epsilon}}^{-1}\simeq 100$. The background magnetic field is mainly perpendicular to the $x-y$ plane: ${\mathbf{B}_{0}}=B_{0}(\sin\theta,0,\cos\theta)$, where $\theta=\cos^{-1}\left[\left(\mathbf{B_{0}}\cdot{\hat{\bf z}}\right)/B_{0}\right]=6^{\circ}$ and $B_{0}=|\mathbf{B_{0}}|$.

In the initial conditions, ions are isotropic and homogeneous (Maxwellian velocity distribution function in each spatial point) for both kinetic simulations.

Large amplitude magnetic ${\bf b}$ and bulk velocity ${\bf u}$ perturbations are introduced, while no density perturbations are taken into account (which implies nonzero total pressure fluctuations). Initial perturbations consist of two Alfvénic wave packets with opposite velocity-magnetic field correlation. The packets are separated along $x$ and, since $B_{0,x}\neq 0$, they counter-propagate. The nominal time for the collision, evaluated with respect to the center of each wave packet, is $\tau\simeq 58.9$.

The magnetic field perturbation $\mathbf{b}$ has been built in such a way that $\mathbf{B_{0}}\cdot\mathbf{b}=0$ is satisfied in each spatial point. Then the velocity field perturbation $\mathbf{u}$ is built by imposing that $\mathbf{u}$ and $\mathbf{b}$ are correlated (anti-correlated) for the wave packet which moves against (along) $B_{0x}$. A detailed discussion about the properties of the initial perturbations can be found in Paper I. The condition $B=|{\bf B}|=const$ is not satisfied by our initial perturbations, while this condition would be a requirement in defining a large amplitude Alfvén mode in the context of a compressible MHD model. This suggests that pressure and density fluctuations are generated during the wave packet evolution.

The perturbations intensity is $\langle b\rangle_{rms}/B_{0}=0.2$, therefore the Mach number is $M_{s}=\langle u\rangle_{rms}/v_{th,p}=0.4$. The intensity of fluctuations with respect to the in-plane field $B_{0x}$ is quite strong, with a value of about $2$. This last parameter can be associated with $\tau_{nl}/\tau_{coll}$ (characteristic nonlinear time $\tau_{nl}$; characteristic collision time $\tau_{coll}$), whose value gives insight about the type of turbulence which could be generated. Here $\tau_{nl}/\tau_{coll}\simeq 0.5$, hence nonlinear effects are important to approach a strong turbulence scenario.

3 Numerical results: a comparison between several codes

In this Section we focus on the description of the results of the four different simulations (MHD, HMHD, HVM and HPIC) by focusing on some “fluid”-like diagnostics which help to understand the system dynamics and, also, to compare the numerical codes.

Figure 1 reports a direct comparison between the simulations, showing the contour plots of the out-of-plane component of the current density $j_{z}=({\bf\nabla}\times{\bf B})\cdot{\hat{\bf z}}$. Vertical columns from left to right in Fig. 1 refer to MHD, HMHD, HVM and HPIC simulations, respectively; while each horizontal row refers to a different time instant: $t=29.5$ (a), $t=\tau=58.9$ (b), $t=70.7$ (c) and $t=98.2$ (d).

Significant differences are recovered in the MHD case with respect to the HMHD, HVM and HPIC runs. While the MHD evolution is symmetric with respect to the center of the $x$ direction, in the other cases this symmetry is broken also before the wave packets interaction due to the presence of dispersive effects which differentiate the propagation along and across the background magnetic field. Moreover, during the wave packets overlap [Fig. 1(b)], smaller scales structures are formed in the HMHD and the HVM cases with respect to the pure MHD evolution, while the HPIC run - despite it recovers several significant features of the wave packets interaction - suffers the presence of particles thermal noise, which has been artificially smoothed out in Fig. 1.

After the collision [Fig. 1 (c) and (d)], the difference between the MHD and the other simulations becomes stronger. In particular, some vortical structures at the center of the spatial domain are recovered in the HMHD and HVM cases, in contrast to the pure MHD case. Moreover, the Vlasov simulation exhibits some secondary ripples in front of each wave packet whose nature could be related to some wave-like fluctuations. These secondary, low-amplitude ripples are not recovered in the other simulations: in fact, they cannot be appreciated in the HPIC run where the noise prevents the formation of such structures while, in the Hall simulation, they are only roughly visible. The nature of these low-amplitude ripples is compatible with a KAW-like activity and will be reported in detail in a separate paper.

In order to compare models and codes, we display, in Fig. 2, the temporal evolution of the energy variations $\Delta E$. Black, red and blue lines indicate respectively the kinetic $\Delta E_{kin}$, thermal $\Delta E_{th}$ and magnetic $\Delta E_{B}$ energy variations, while each panel from (a) to (d) refers to the MHD, HMHD, HVM and HPIC runs, respectively. The evolution of $\Delta E_{kin}$ and $\Delta E_{B}$ is quite comparable in all the performed simulations and, in the temporal range where the wave packets collide, magnetic and kinetic energy is exchanged. On the other hand, the evolution of the thermal energy $\Delta E_{th}$ differs in the HPIC case compared to the other simulations. Indeed, $\Delta E_{th}$ remains quite close to zero for all the simulations except for the HPIC run, where it grows almost linearly for the presence of numerical noise. It is worth to note that, as the number of particles increases, the evolution of $\Delta E_{th}$ would get closer to the one obtained in the MHD, HMHD and HVM simulations.

The scenario described by Moffatt and Parker is also based on the property, in ideal incompressible MHD, that two wave packets separately conserve energy, which is equivalent to conservation of both total energy and cross helicity $\sigma_{c}$. It is natural therefore to examine evolution of cross helicity, as well as the evolution of the residual energy $\sigma_{r}$, which gives information about the relative strength of magnetic fluctuations and the fluid velocity fluctuations. Figures 3(a–b) show the temporal evolution of normalized residual energy $\sigma_{r}$ (a) and the normalized cross-helicity $\sigma_{c}$ (b). These quantities are defined as follows: $\sigma_{r}=(e^{u}-e^{b})/(e^{u}+e^{b})=2e^{u}/(e^{+}+e^{-})$, where $e^{r}=e^{u}-e^{b}$, $e^{\pm}=\langle(\mathbf{z^{\pm}})^{2}\rangle/2$ ($\mathbf{z^{\pm}}=\mathbf{u}\pm\mathbf{b}$), $e^{u}=\langle\mathbf{u}^{2}\rangle/2$ and $e^{b}=\langle\mathbf{b}^{2}\rangle/2$; $\sigma_{c}=(e^{+}-e^{-})/(e^{+}+e^{-})=2e^{c}/(e^{u}+e^{b})$, being $e^{c}=\langle\mathbf{u}\cdot\mathbf{b}\rangle/2$. In each panel of Fig. 3, black, dashed blue, dashed green and red lines refer to MHD, HMHD, HVM and HPIC cases, respectively.

Figure 3 (a) shows the evolution of the normalized residual energy $\sigma_{r}$, which is similar in all the simulations. In particular $\sigma_{r}\simeq 0$ in the initial stage,it oscillates during the wave packets overlapping, and finally it returns to $\sigma_{r}\simeq 0$ after the collisions. The $\sigma_{r}$ oscillations are well correlated with the oscillations of $\Delta E_{B}$ and $\Delta E_{kin}$ seen in Fig. 2.

Deeper insights are revealed by the evolution of the cross-helicity $\sigma_{c}$. Indeed, for ideal incompressible MHD, the cross helicity remains constant, and for this initial condition, $\sigma_{c}=0$. Here, $\sigma_{c}$ is well-preserved in the MHD run despite this simulation being compressible. This means that the compressible effects, introduced here by the fact that initial perturbations are not pressured balanced, are not strong enough to break the $\sigma_{c}$ invariance. On the other hand, for the remaining simulations (HMHD, HVM and HPIC), $\sigma_{c}$ is not preserved: i) it shows a jump around $t=\tau=58.9$, due to the presence of kinetic and dispersive effects, and ii) there is an initial growth of $\sigma_{c}$ followed by a relaxation phase. It seems also significant to point out that, the initial growth of $\sigma_{c}$ occurs faster in the kinetic cases compared to the HMHD one. This may reflect the fact that, the initial condition evolves differently in the Hall MHD simulation compared to the kinetic runs.

In order to understand the role of the Hall physics, we computed the normalized generalized cross helicity $\sigma_{g}=2e^{g}/(e^{u}+e^{b})$, where $e^{g}=0.5\ \langle\mathbf{u}\cdot\mathbf{b}+{\tilde{\epsilon}}\mathbf{\omega}\cdot\mathbf{u}/2\rangle$, and $\omega={\bf\nabla}\times\mathbf{u}$, which is an invariant of incompressible HMHD (Turner, 1986; Servidio et al., 2008). Figure 3(c) displays the temporal evolution of $\sigma_{g}(t)$ for the MHD (black), the HMHD (dashed blue), HVM (dashed green) and HPIC (red) simulations. Note that the evolution of $\sigma_{g}$ is trivial for the MHD simulation where, since $\tilde{\epsilon}=0$, $\sigma_{g}=\sigma_{c}$. Moreover, it can be easily appreciated that, for the HMHD case, $\sigma_{g}$ is almost preserved and does not exhibit any significant variation due to the collision itself, even though it shows a slight increase in the initial stages of the simulation followed by a decay towards $\sigma_{g}=0$ [similar to the growth of $\sigma_{c}$ recovered in Fig. 3(b)]. On the other hand the two kinetic cases, which exhibit a similar behavior, show a fast growth of $\sigma_{g}$ in the initial stage of the simulations followed by a decay phase [similar to the growth of $\sigma_{c}$ recovered in Fig. 3(b)]; then, during the collision, $\sigma_{g}$ significantly increases. We may explain the evolution of $\sigma_{c}$ and $\sigma_{g}$ as follows. In the MHD run, compressive effects contained in the initial condition as well as compressible activity generated during the evolution are not strong enough to break the invariance of $\sigma_{c}$ (i.e. of $\sigma_{g}$). Instead, in the Hall MHD simulation, the first break of the $\sigma_{c}$ invariance observed in the initial stage of the simulation cannot be associated with the Hall effect since also $\sigma_{g}$ is not preserved in this temporal region and $\sigma_{c}$ and $\sigma_{g}$ have a similar evolution. On the other hand, the jump recovered in $\sigma_{c}$ around $t\simeq\tau=58.9$ is significantly related to the Hall physics. In fact, since $\sigma_{g}$ does not exhibit a similar jump at $t\simeq\tau$, we argue that the physics which produces the growth of $\sigma_{c}$ is the Hall physics (which is taken into account in the invariance of $\sigma_{g}$). Finally, the production of both $\sigma_{c}$ and $\sigma_{g}$ recovered in the kinetic simulations cannot be completely associated with the Hall effect (which, of course, is still present) but kinetic and compressive effects may have an important role.

In order to explore the role of small scales into the dynamics of colliding wave packets, we computed the averaged mean squared current density $\langle j_{z}^{2}\rangle$ as a function of time. This quantity indicates the presence of small scale activity (such as production of small scale current sheets), and is reported in Figure 4 for all the simulations. As in the previous figures, black, blue dashed, green dashed and red lines refer to the MHD, HMHD, HVM and HPIC cases, respectively. All models show a peak of $\langle j_{z}^{2}\rangle(t)$ around the collision time $t\simeq\tau$ due to the collision of wave packets. After the collision, some high-intensity current activity persists in all the simulations, which present a qualitatively similar evolution of $\langle j_{z}^{2}\rangle(t)$.

Other quantities that provide physical details about our simulations are $\epsilon_{\rho}=\langle\delta\rho^{2}\rangle$ (compressibility) and the enstrophy $\epsilon_{\omega}=\langle\mathbf{\omega}^{2}\rangle/2$ (fluid vorticity $\omega$). Note that $\delta\rho=\rho-\langle\rho\rangle$. Figure 5 reports the temporal evolution of $\epsilon_{\rho}$ (a) and $\epsilon_{\omega}$ (b) for all the runs. Black, blue dashed, green dashed and red lines indicate respectively the MHD, HMHD, HVM and HPIC cases. The $\epsilon_{\rho}$ evolution shows that density fluctuations peak around $t\simeq 63.8$ and $t\simeq 83.4$. The two peaks are respectively associated with the collision between the packets and with the propagation of magneto-sonic fluctuations generated by the initial strong collision with provide a sort of “echo” of the original interaction. Moreover, from the initial stage of the simulations, $\epsilon_{\rho}$ exhibits some modulations, which are produced by the absence of a pressure balance in the initial condition. In fact, as packets start to evolve, low-amplitude fast perturbations (clearly visible in the density contour plots, not shown here) propagate across the box and collide faster compared to the “main” wave packets themselves. Moreover, by comparing the different simulations, one notices that, for $t\lesssim 20$, kinetic and Hall runs tend to produce a similar evolution of $\epsilon_{\rho}$, slightly bigger compared to the MHD case. Then, around $t\simeq 20$, the HMHD run displays a stronger compressibility with respect to the kinetic cases. This difference is probably due to the presence of kinetic damping phenomena which occur in the kinetic cases.

The enstrophy $\epsilon_{\omega}$ is displayed in Fig. 5(b). All the runs exhibit a similar evolution of $\epsilon_{\omega}$ up to the wave packet collisions. Then MHD and HMHD cases exhibit a quite similar level of $\epsilon_{\omega}$, slightly bigger compared to the one recovered in the HVM and HPIC cases, where probably kinetic damping does not allow the formation of strong vortical structures at small scales.

It is interesting to compare different simulations also by looking at power spectral densities (PSDs). Figures 6 show the magnetic energy PSD integrated along $k_{y}$ $E_{b,y}(k_{x})=\sum_{k_{y}}E_{b}(k_{x},k_{y})$ (left column) and along $k_{x}$ $E_{b,x}(k_{y})=\sum_{k_{x}}E_{b}(k_{x},k_{y})$ (right column); while each row respectively refers to $t=29.5$ (top row), $t=\tau=58.9$ (center row) and $t=98.2$ (bottom row). The cyan dashed line shows the $k^{-5/3}$ slope for reference while, in each panel, black, blue, dashed green and red lines indicate respectively MHD, HMHD, HVM and HPIC simulations. Moreover, to compare the two wave-number directions, gray lines in each panel report the corresponding PSD obtained from the MHD run, reduced in the other direction [for example, in the top row left panel, the gray line refers to $E_{b,x}(k_{y})$ for the MHD simulation while other curves in th same panel report $E_{b,y}(k_{x})$]. It is interesting to note that, at $t=29.5$, all the simulations exhibit a steep spectrum in $E_{b,y}(k_{x})$, related to the initial condition, while the difference in power between $E_{b,y}(k_{x})$ and $E_{b,x}(k_{y})$ - the latter being significantly smaller than the former - tends to reduce in the final stages of the simulations. This suggest, as described in Paper I, the presence of nonlinear couplings which cause spectra to become more isotropic. Moreover the comparison between the different simulations indicates that the dynamics at large scales is described in a quite similar way by all runs, while at small scales some differences are revealed. In particular, the HPIC simulation is affected by particles noise while the HVM run contains more energy at small scales compared to MHD and HMHD. Note that the presence, in MHD and HMHD runs, of an explicit resistivity prevents the population of small scales.

To summarize this section, we compared our numerical codes by analyzing some global diagnostics and we conclude that the Moffatt-Parker scenario is quite well satisfied by MHD. However, other intriguing characteristics are observed when one moves beyond the MHD treatment. Moreover, the comparison between kinetic codes suggests that HVM and HPIC simulations display qualitatively similar features at large scales. However, when one aims to analyze the dynamics at small scales, HPIC simulations suffers from thermal particle noise. Magnetic energy spectra differ in the HPIC case compared with the HVM case. Moreover, by comparing the $j_{z}$ contour plots, one can easily appreciate how HPIC simulation is affected by particle noise. Based on these considerations, we continue the analysis of the kinetic features produced in Alfvén wave packet collisions by focusing only on the HVM case.

4 Kinetic features recovered during the wave packets interaction

We begin a description of the kinetic signatures present in the Vlasov simulation by looking at the temperature anisotropy. Fig. 7 reports the contour plots of the temperature anisotropy $T_{\perp}/T_{\parallel}$, where the parallel and perpendicular directions are evaluated in the local magnetic field frame (LBF), at four time instants: $t=29.5$ (a), $t=\tau=58.9$ (b), $t=70.7$ (c) and $t=98.2$ (d). Clearly temperature anisotropy is present even before the main wave packets collision, due to the fact that the initial wave packets are not linear eigenmodes of the Vlasov equation and, hence, their dynamical evolution leads to an anisotropy production. Moreover, a more careful analysis suggests that the left wave packet tends to produce regions where $T_{\perp}/T_{\parallel}<1$ close to the packet itself (which, as it can be appreciated in Fig. 1, is localized around $x\simeq 9.5$), while the right wave packet (localized around $x=15.7$) is characterized by $T_{\perp}/T_{\parallel}>1$. The presence of different temperature anisotropies is related to the broken symmetry with respect to the center of the $x$ direction. Indeed, the dynamics of the wave packets is different if they move parallel or anti-parallel to $B_{0,x}$. This produces the different temperature anisotropy recovered in the top panel of Fig. 7.

When the packets collide [Fig. 7 (b)], sheets characterized by a strong temperature anisotropy ($T_{\perp}/T_{\parallel}>1$) are recovered, correlated spatially with the current density structures. Then, at $t=70.7$ [Fig. 7 (c)], wave packets split again and a region, localized at $(x,y)\simeq(14.3,1.0)$, where the temperature anisotropy suddenly moves from values $T_{\perp}/T_{\parallel}<1$ towards $T_{\perp}/T_{\parallel}>1$ ones is present. We will show that this region also shows the presence of strong departures from the equilibrium Maxwellian shape. At the final stage of the simulations [Fig. 7 (d)], each wave packet continues traveling, accompanied by a persistent level of temperature anisotropy, which is, indeed, well correlated with the current structures.

It is interesting to point out that, beyond the presence of temperature anisotropies, regions characterized by a temperature nongyrotropy are also recovered. Many methods have been proposed by evaluating the temperature nongyrotropy (Aunai et al., 2013; Swisdak, 2016). Here we adopt the measure $D_{ng}$ (Aunai et al., 2013), which is proportional to the root-mean-square of the off-diagonal elements of the pressure tensor. Fig. 8 shows the contour plots of $D_{ng}$ at four time instants: $t=29.5$ (a), $t=\tau=58.9$ (b), $t=70.7$ (c) and $t=98.2$ (d). Moreover, as for the temperature anisotropy, the evolution of the two wave packets tends to produce temperature nongyrotropies also before the wave packets collision [Fig. 8(a)]. Then, during the collision [Fig. 8(b)–(c)], the temperature nongyrotropy becomes more intense and it is also quite well correlated with the current structures. At the final stage of the simulation [Fig. 8(d)], each wave packet is connoted by a level of nongyrotropy which is quite bigger compared to the value before the collision. The presence of nongyrotropic regions suggests that it is fundamental to retain a full velocity space where the VDF is let free to evolve and, eventually, distort.

To further support the idea that kinetic effects are produced during the interaction of the wave packets, we computed the $L^{2}$ norm difference (Servidio et al., 2012; Greco et al., 2012; Servidio et al., 2015):

[TABLE]

which measures the displacements of the proton VDF $f({\bf x},{\bf v},t)$ with respect to the associated Maxwellian distribution function $f_{M}({\bf x},{\bf v},t)$, built such that density, bulk speed and total temperature of the two VDFs are the same. Figure 9 reports the evolution of the $\epsilon_{max}(t)=max_{(x,y)}\epsilon(x,y,t)$ as a function of time. Clearly, as for the previous proxies of kinetic effects, also $\epsilon_{max}$ moves away from zero in the early phases of the simulation due to the fact that the initial-condition is not a Vlasov eigenmode. Moreover, after the first initial variation, $\epsilon_{max}$ is almost constant until wave packets interact. Then, during the collision, $\epsilon_{max}$ grows and reaches its maximum at $t=70.7$, quite later than the wave packet collisions. Then it decreases and saturates at a value about two times bigger than the value before the collision, thus suggesting, again, that there is “net” production of non-Maxwellian features in the VDF during the wave packets interaction.

In order to appreciate the structure of $\epsilon$ in the spatial domain, the left panel of Figure 10 shows the contour plot of $\epsilon(x,y,t)$ at the time instant $t=70.7$ (when $\epsilon$ reaches its maximum value). The $\epsilon$ contours are correlated with the current structures and with the anisotropic/nongyrotropic regions. Moreover, a blob-like region where $\epsilon$ reaches its maximum is present. This area is associated with the region where the temperature anisotropy moves from $T_{\perp}/T_{\parallel}<1$ values towards $T_{\perp}/T_{\parallel}>1$ ones [See Fig. 7(c)]. In this area the VDF strongly departs from the Maxwellian. The right panel of Fig. 10 shows the three dimensional isosurface plot of the VDF at $t=70.7$ and in the spatial point $(x^{*},y^{*})$ where $\epsilon(x^{*},y^{*},t=70.7)=\max_{(x,y)}\epsilon(x,y,t=70.7)$. A strong beam, parallel to the local magnetic field direction, is observed in the VDF in Fig. 10. The drift speed of the beam is about ${\tilde{c}}_{A}$. The production of a beam due to the interaction of two wave packets has been also pointed out by He et al. (2015).

5 Conclusion

In this paper we have described the interaction of two Alfvénic wave packets by means of MHD, Hall MHD and hybrid kinetic simulations of the same physical configuration. Kinetic simulations have been performed with two different codes: an Eulerian Vlasov-Maxwell code (Valentini et al., 2007) and hybrid Particle-in-cell code (Parashar et al., 2009). By approaching the Parker & Moffatt problem within several physical frameworks, we have comparatively analyzed different effects (compressive, Hall and kinetic) which contribute to the general, complex puzzle.

The analysis performed in Paper I where MHD and Vlasov simulations were compared has been here extended by including the Hall MHD framework. In particular, we showed how, moving beyond the pure MHD treatment, dispersion as well as kinetic effects play an important role. Furthermore, the analysis of HMHD and kinetic runs allows to separate the presence of dispersive and purely kinetic features. It is also interesting to note that the compressive activity is different in the Hall case compared to the kinetic runs, indicating that some kinetic damping processes may be active in the Vlasov simulation. A separate type of comparison is afforded by comparative analysis of the HPIC and HVM runs, While these methods should describe, approximately, the same physics, i.e., the Vlasov treatment of collsionless plasma dynamics, the comparison between the different codes is interesting from a methodological perspective, and therefore represents a contribution to the Turbulent Dissipation Challenge (Parashar et al., 2015b). In particular, the two kinetic simulations performed are able to take into account the dynamics which occurs at large spatial scales and their comparison is successful in this range of scales. However the PIC runs lacks accuracy when smaller spatial scales are produced by the collision of the two wave packets, thus indicating that the Eulerian approach better describes the dynamics of the system at small spatial scales. Of course the comparison is expected to become better if the number of particle per cell in the PIC simulation gets bigger (Camporeale & Burgess, 2011; Franci et al., 2015).

Based on the last consideration, we have analyzed the production of kinetic signatures by focusing only on the HVM simulation. Several proxies which are routinely adopted for highlighting the presence of kinetic features indicate that wave packets tend to produce kinetic effects such as temperature anisotropies and nongyrotropies also before the main wave packets interaction. This is related to the fact that the initial condition, consisting of quasi-Alfvénic wave packets, is neither a Vlasov equilibrium nor a Vlasov eigenmode. Therefore it dynamically leads to the production of kinetic features.

The analysis of kinetic effects before and after the main wave packets collision indicates that some kinetic features are enhanced by the collision itself and each wave packet is significantly characterized by a strong degree of non-thermal signatures. In particular the presence of temperature nongyrotropies suggests that descriptions based on reduced velocity space assumptions may partially fail the description of such features. Finally, similarly to the observations of He et al. (2015), during the wave packet collision, a beam in the velocity distribution function is observed to form along the direction of the local magnetic field. This characteristic may connect our results with the general scenario of wave packets observed in natural plasmas such as the solar wind.

Research is supported by NSF AGS-1063439, AGS-1156094 (SHINE), AGS-1460130 (SHINE), and NASA grants NNX14AI63G (Heliophysics Grandchallenge Theory), and the Solar Probe Plus science team (ISIS/SWRI subcontract No. D99031L), and Agenzia Spaziale Italiana under the contract ASI-INAF 2015-039-R.O “Missione M4 di ESA: Partecipazione Italiana alla fase di assessment della missione THOR”. Kinetic numerical simulations here discussed have been run on the Fermi supercomputer at Cineca (Bologna, Italy), within the ISCRA-C project IsC37-COLALFWP (HP10CWCE72) and on the Newton parallel machine at University of Calabria (Rende, Italy).