Kinetic plasma turbulence during the nonlinear stage of the Kelvin-Helmholtz instability
Using a full kinetic, implicit particle-in-cell code, iPiC3D, we studied the properties of plasma kinetic turbulence, such as would be found at the interface between the solar wind and the Earth magnetosphere at low latitude during northwards periods. In this case, in the presence of a magnetic field B oriented mostly perpendicular to the velocity shear, turbulence is fed by the disruption of a Kelvin-Helmholtz vortex chain via secondary instabilities, vortex pairing and non-linear interactions.
We found that the magnetic energy spectral cascade between ion and electron inertial scales, and , is in agreement with satellite observations and other previous numerical simulations; however, in our case the spectrum ends with a peak beyond due to the occurrence of the lower hybrid drift instability. The electric energy spectrum is influenced by effects of secondary instabilities: anomalous resistivity, fed by the development of the lower hybrid drift instability, steepens the spectral decay and, depending on the alignment or anti-alignment of B and the shear vorticity, peaks due to ion-Bernstein waves may dominate the spectrum around . These waves are generated by counter-streaming flow structures, through flux freezing also responsible for reconnection of the in-plane component of the magnetic field, which then generates electron pressure anisotropy and flattening of the field-aligned component of the electron distribution function.
With a large separation between the scales of system dynamics and dissipation, astrophysical and space flows are often in a turbulent state, at the same time, they consist dominantly of fully ionised plasmas. A major implication of the previous sentence is the multiscale aspect of our descriptions and calculations.
Plasma species are subject to a variety of physical processes, to which one can associate a typical length scale (and a correlated time scale): collisional length scales, the scales at which waves and particle motions decouple (inertial length) and magnetisation scales (Larmor radius). Depending on temperature, density and field strength, these plasma scales can be ordered in various ways. In a turbulent flow, characterised by the nonlinear transport of energy between scales, an injection scale and dissipative scales are also to be taken into account. While in a fluid description the dissipation scale is a single length set by collisions, in weakly collisional plasmas several energy redistribution paths can coexist, and energy dissipation is the result of a complex interaction of micro-instabilities.
In Astronomy, one would like to think in terms of the larger (observable) scales and events, far above the necessary resolution set by the plasma-kinetic equations - and often above the scales associated with fluid turbulence as well. However, to correctly predict the transport coefficients used in such a macro desciption, one has to look at the actual interaction scales. In this paper, we focus on the energy and momentum transfer between solar wind and the magnetosphere near the equator when the earth ans solar wind magnetic fields are approximately in the same direction Southwood74 (); Chen93 (); Hasegawa04 (); Pope09 (); Masters09 (); Sundberg10 (). This plasma shear flow configuration is unstable to the Kelvin-Helmholtz instability and has been extensively studies by means of fluid simulations(Nagano79, ; Miura82, ; Nakamura10, ; Henri12, ) However, the scales on which the instability operates, are below the fluid limit and sheath width, growth rates or transport coefficients are governed by plasma kinetic laws.
Typically plasma kinetic simulations are limited to describing very small scale physical phenomena, as explicit calculation imposes severe constraints on the resolution, enforcing us to resolve the electron associated temporal and spatial scales. Working with an implicit formulation of the equations, those resolution restrictions are loosened brackbill1982implicit (), allowing us to solve larger systems at an affordable computational cost. The current work attempts to bridge the range from the largest electron governed scales, becoming available in the latest satellite data Alexandrova08 (); Safrankova13 (), to scales where one could consider the plasma as a single turbulent fluid.
There are different ways to self-consistently drive turbulent flows that have been used in numerical simulations: using an external driver to inject energy in the system gressel2008 (), starting from the decay of a large amplitude perturbation delsordo2010 () or letting an instability evolve nonlinearly dalziel1999 (). In this work, we use this last approach, namely the nonlinear stage of the magnetized Kelvin-Helmholtz instability. The initial shear flow is intended here as a source of free energy to drive a turbulent stage in order to study the properties of magnetized plasma turbulence.
Our objectives here are threefold: first we wish to observe the consequences of working with a reduced resolution, enabled by using an implicit scheme, by comparing with simulations fully resolving kinetic scales Karimabadi13 (). Second, we are aiming at a sufficient separation of scales to achieve a clear energy cascade scaling, comparable to observational data. And finally, we want to assess the observed kinetic effects in terms of more macroscopic behaviour.
The paper is divided as follows. After a brief review of plasma turbulence studies in section II, the model used in this study is described in section III and the outcome of this full kinetic simulation is presented in section IV. Results are analysed in section V. Our findings are concisely summarised in section VI.
Ii Review of plasma turbulence studies
Let us first have a brief look at the general nature of turbulence in magnetised plasmas to see what kind of behaviour we may expect in our own simulations, starting from the outer scale.
While our simulation does not have any collisions, given a sufficient domain size, it should be possible to retrieve the Alfvénic cascade observed in 2D magneto-hydrodynamic (MHD) turbulence, as it does not rely on collisions and exists deep into the kinetic regimeAlex2009 ().
In MHD context,’Turbulence’ typically refers to the strong, balanced turbulence cascade of magnetic and specific kinetic energy in an incompressible plasma. Strong turbulence implies dominance of the non-linear term over the dissipative term, as opposed to weak turbulence GaltierNazarenko2000 (); Galtier2006 (); Mininni07 (); PerezBoldyrev08 (); BoldyrevPerez09 (); Saur02 (). Incompressibility results in a straightforward coupling between the evolution of velocities and magnetic fields in the MHD and electron-MHD (EMHD) subspectra, allowing us to describe the problem in Elsässer variables , i.e. using linear combinations of the (electron) velocity and the magnetic field. Studies into the effects of compressibility ChoLaz02 (); ChoLaz03 (); Alex2009 (); Kowal10 (), concluded that the slow compressible modes passively exchange energy with the incompressible Alfvénic modes, while the fast modes are energetically less relevant, leaving the spectral scaling exponent unchanged. Balanced turbulence is defined by equation of the Elsässer energies , implying the absence of cross-helicity. According to general consensus Lithwick07 (); BeresLaz08 (); PodestaBhat10 (); PerezBoldyrev09 () unbalanced turbulence does not change the scaling, only relative amplitude of the Elsässer energies.
While 3D MHD Alfvénic turbulence is dividing the scientific world GS95 (); Haugen04 (); Mininni07 (); PodestaBhat10 (); BeresLaz10 (); PerezBoldyrev10 (); Boldyrev11 (); Alexandrova13 (), the scaling power of -3/2 for its 2D counterpart is well established from theory Iroshnikov63 (); Kraichnan65 (), in numericsBiskamp1989 (); Politano1989 (); Biskamp2001 () and in experiments Sommeria1986 ().
For magnetic and electron kinetic energy in the EMHD regime a -7/3 scaling (kinetic Alfvén turbulence) was predicted and found in EMHD simulations, for 2D Biskamp96 () as well as 3D Biskamp99 () and in observations Bale05 (). However, more recent solar wind observations suggest a steeper -8/3 spectrum Alexandrova08 (); Safrankova13 (), for this mismatch several possible explanations existHowes11 (); BoldyrevPerez12 (). Explicit PIC simulations Karimabadi13 () suggest that the -8/3 spectrum also exists in 2D. Corresponding to these different explanations, an electric field energy scaling with exponent -1/3Howes11 () or -2/3BoldyrevPerez12 () is expected. At scales below the electron inertial length, the magnetic field spectrum drops strongly due to Landau damping Howes11 ().
Iii Numerical setup
To drive a fully developed turbulent state starting from a shear layer requires an evolution on a time scale of many hundred of ion gyro-periods, much larger than the typical time scales that can be covered by standard Particle-in-Cell (PIC) methods. To overcome this problem and cover such a large period of time, we make use of the fully kinetic, fully electromagnetic Particle-in-Cell code iPIC3D iPIC3D (), which implements the moment implicit method to suppress numerical instabilities when using large simulation time steps brackbill1982implicit (). In the following, all quantities are normalized to ion quantities: the ion gyro-frequency, , the ion inertial length, or , and the Alfvén velocity . A reduced ion-to-electron mass ratio is used for computational reasons.
We consider a 2D physical space with 3D vector fields corresponding, in phase space, to a 2D-3V configuration. The size of the numerical box is using grid points, for a spatial resolution of (the electron inertial length). The initial sheared velocity field is characterized by a double shear layer where the velocity varies from to . Such a double shear layer is used in order to impose periodic boundary conditions, thus avoiding spurious effects driven by the boundary conditions. The shear layers are equispaced and located at and . The velocity profile reads:
where the maximum velocity field strength is corresponding to a velocity jump . We take a shear scale length of a few ion skin depths (or ion inertial lengths), namely . Outside from the shear layers where the system is homogeneous, the initial magnetic field is , where and ; the initial ion and electron thermal velocities, and respectively, are isotropic; the initial density is uniform and equal to one. Quasineutrality is imposed everywhere at the beginning of the simulation. The plasma beta is , so that the ion inertial length and gyroradius are roughly of same order.
One of the main difficulties in the kinetic modeling of shear flows is the choice of the initial conditions since only a few kinetic equilibria are known for shear flow configurations CaiStoreyNeubert1990PhFlB (). On the other hand, when using a force balance MHD-like equilibrium as initial condition, the tensor pressure reacts very fast at scales not far from the ion kinetic scale Henrietal2013POP () thus introducing strong fluctuations in the system. For these reasons, to mimic situations where the shear length is a few ion Larmor radii wide, such as the magnetosheath - magnetosphere boundary, we chose to implement an “extended two-fluid equilibrium”, that retains first order corrections in terms of Finite Larmor Radii (FLR) effects Cerrietal2013POP (). Although it is not a proper kinetic equilibrium, this setup has been proven to be particularly efficient in preventing the generation of unwanted artifacts driven by the kinetic reaction to an initial fluid-like setup Henrietal2013POP (); Cerrietal2013POP (). This consideration is particularly important in the range of parameters used in this study. The profiles of the thermal velocities are set according to the profiles of the diagonal elements of the pressure tensor in the “Cerri equilibrium”Cerrietal2013POP () and the particles are loaded with a weight corresponding to the equilibrium density profile.
In each cell, 100 particles are loaded for each species, for a total of about 1 billion particles. The simulation ran for a few million CPU hours on 8192 cores on FERMI at CINECA until the system reaches a time .
The numerical box dimension has been chosen so that the fastest growing mode is m=2, resulting in the formation of two Kelvin-Helmholtz vortices in each shear layer at the beginning of the nonlinear phase. The vortices are eventually disrupted during the nonlinear evolution of the system, forming two turbulent layers. We hereafter concentrate on the properties of plasma turbulence at the turbulent stage of the full kinetic simulation.
iv.1 Effect of vorticity orientation
One would expect a MHD setup with the same initial conditions for magnetic and flow fields to behave symmetric in the y-direction. In figure 1 we see that especially for the electric field there are strong differences between the two turbulent layers, while for the ion velocity there is almost no difference.
The differences become more clear when we study the y-dependence by averaging over x and time. In figure 2 we see a mean density and magnetic field decrease in the upper part of the domain. Figure 3 shows a mean electric field corresponding to , modified in the top turbulent layer. We find that the energy in the fluctuations (with a finite wavelength in x) in the top layer is much larger than the energy associated with the mean field, while in the bottom layer the x-dependent contribution to the electric energy is weak.
The vorticity of the shear flow in the top layer is in the same direction as the magnetic field. This means that, for the electrons, the Lorentz force is pointing in the same direction as the acceleration by the flow, enhancing their rotational motions. The resulting electron currents also create a magnetic field that reduces the initially imposed one. The magnetic field and flow are largely frozen-in, making the density follow the behaviour of B (). In the bottom layer, vorticity and B are anti-aligned and the forces on the electrons counteract, approximately cancelling out. This results in asymmetry in the electron currents between the top and bottom of the domain, seen in figure 3 on the right. The electric energy shows similar asymmetric behaviour.
To explain the similarity, this collisionless plasma would need some kind of resistivity. A comparison of the spatial behaviour of electric field and current, Figure 4, does show a strong correlation.
Resistivity can be measured through taking an average of the (collisionless) momentum balance of the electrons:
As the electron mass is small, we can drop the first term on the left, and the remaining terms are dominated by the electric field. Estimating the diffusivity (actually a tensor) as the ratio , we find, shown in figure 5, localised values around unity, surrounded by regions of weaker diffusivity (by approximately a factor 10). The most diffusive regions coincide with strong electric fields and small scale fluctuations.
iv.2 Spectral analysis
We studied the 2D spectra of two subdomains within the upper and bottom turbulent layer. The non-periodic boundaries in y were mitigated by the use of a Hamming window, though for high angles with the x-axis some noise exists.
In figure 6 we show the spectral scaling of the components of the magnetic field for different orientations. The behaviour is largely isotropic and rather similar for both layers. Scaling between the ion and electron scales is close to the expected value of -8/3, found also in Haynes et al. Haynes2014 () and Karimabadi et al.Karimabadi13 (). Even though the domain size is on the order of , the resulting forcing scales provided by the instability are insufficiently large to accommodate a fluid like energy cascade. One can see a slight bump in the spectrum of , stronger for high angles, between and .
Contrary to the magnetic field spectrum, and as could already be anticipated from the previous section, the k-space behaviour of the electric field is vastly different between the top and bottom layer In the lower layer spectrum in figure 7, we see two main features: a strong double peak, a bit beyond , and a bump, a bit beyond . The bump is clearest in where, like in the magnetic spectrum, it is strongest for high angles. The double peak is strongest in , where it favours small angles. The scaling between the and is obscured by these two disturbances - scaling indications on the figures are for reference, rather than fits. In the upper layer, traces of the double peak introduce noise on the measurement, though the, again largely isotropic, scaling is far steeper than expected for a collisionless plasma. The slope between -2 and -3 reminds of the magnetic spectrum and suggests the electric field follows from currents.
The final set of spectra included here, figure 8 shows the electron velocity and density behaviour. The ion density matches the electrons closely, while all ion velocity components drop rather steeply around down to noise levels, as in Karimabadi et al.Karimabadi13 (). The velocity spectrum is largely dominated by the features already observed in the electric field: a double peak in along the -axis a bit beyond in the bottom layer (parallel to its strongest appearance in E) and a bump in (orthogonal to the orientation of its strongest appearances in both E and B), along the -axis a bit beyond in both layers. The lowest four panels in figure 8 indicate the density scaling, we again notice an increase around the bump location in the other quantities. Only by looking at the scaling of the density difference between electrons and ions, do we recover the peaks at .
iv.3 Velocity distribution functions
We calculated the velocity distribution function on a 64 by 128 grid, corresponding to the processor distribution. In figure 9, showing the spatial variation of the kurtosis for the distribution function along respectively and , we see that for the electrons neither orientation is fully Maxwellian. The distribution function is shown in figure 10 for different locations in the plane. The velocity distribution in the plane is approximately isotropic and appears Maxwellian with a high energy tail of varying strength, while is a flat-top distribution. An approximation of the temperature ratio between the in and out of the plane dependances of the distribution function, using the corresponding variances (see figure 11) shows regions with large ratios of parallel to perpendicular temperature for the electrons. The ion velocity distribution is approximately Maxwellian and does not show strong anisotropies.
v.1 Pressure anisotropy and electrostatic waves
Visually we observed, and through the spectra we confirmed the presence of electrostatic waves with preferred propagation along . The waves are longitudinal and largely perpendicular to the magnetic field. They occur in regions with electron pressure anisotropy (the ion pressure is approximately isotropic), namely where the pressure component along B is larger, by a factor of about 5, than the pressure components perpendicular to B, see figure 12.
The double peak in the electric energy spectrum, figure 7, suggests the presence of the first and second proton cyclotron harmonic Bernstein waves. For propagation at high but less than perpendicular angle with the magnetic field these can result in a double peaked spectrum, becoming a single flat peak for orthogonal propagation LiHabbal2001 (); SahraouiBelmontGoldstein2012 (). The spectral peaks for propagation along - generally at a higher angle with B are less pronounced and suggest a combination of the two spectra mentioned.
A cut through Kelvin-Helmholtz-rolls will show alternating flow directions. It was suggested ReynoldsGanguli1998 () that such a geometry would become unstable to ion-Bernstein waves for layer separations on the order of the ion Larmor radius.
With approximately frozen-in magnetic fields, alternating flows correspond to alternating magnetic fields and hence reconnection. In earlier studies Egedal2002 (); EgedalLeDaughton2012 () magnetic reconnection is linked to flattening of the electron distribution and electron pressure anisotropy of the same order as observed here.
While in the current stage of the instability evolution, there are no clean rolls, we do find similar velocity and field structures, see figure 13, in regions with strong anisotropy, ion-Bernstein-waves and a flat topped electron velocity distribution. A contour plot for one of the regions of interest for the relevant quantities is given in figure 14.
v.2 Lower hybrid drift instability
In ,, and we found the generation of waves with k-vector preferrably along y and norm between and , shown in figure 15. This signature agrees with the lower hybrid drift instability (LHDi).
The LHDi is typically associated with anomalous resistivity, which we already proposed in the asymmetry section, figure 3, and anticipated from the steep E-spectrum, figure 7. The LHDi identified regions indeed correspond to the stronger diffusivity locations in our estimate through Ohm’s law, recall figure 5. The strongly localised diffusivity agrees with earlier simulations of the LHDi by Innocenti & Lapenta (2007 innocenti2007 (), and unpublished results).
While our domain size was of the order , the turbulent injection scale, set by the Kelvin-Helmholtz instability, was only a bit larger than and we were not able to resolve a fluid like spectral cascade.
For the magnetic field, we were able to reproduce the kinetic Alfvén wave spectrum as seen in existing simulations and observations.
The behaviour of the electric field is significantly different from what one would expect from a cascade in a collisionless plasma. The deviation is caused by two separate physical phenomena occuring on respectively ion and electron scales.
In the region where the shear layer vorticity and magnetic field are anti-aligned, around , an in-plane, alternating flow structure excites Ion-Bernstein waves, resulting in a peaked electric energy spectrum. This flow structure is also tied to the magnetic field structure, causing reconnection of the in-plane magnetic field, which, in turn, generates anisotropy in the electron velocity.
Troughout the whole domain, near , we find the lower hybrid drift instability (LHDI), showing up in both the magnetic and the electric spectrum. Beside the peak at the resonance wave length, the LHDI also affects the electric field spectrum through the generation of anomalous resistivity, resulting in a far steeper slope.
Acknowledgements.The research leading to these results has received funding from the European Commission’s Seventh Framework Programme (FP7/2007-2013) under the grant agreement SWIFF (project n 263340, www.swiff.eu). We acknowledge the access to the HPC resources of CINECA made available within the Distributed European Computing Initiative by the PRACE-2IP, receiving funding from the European Community’s Seventh Framework Programme (FP7/2007-2013) under grant agreement n RI-283493 and project number 2012071282. This work was supported by the Italian Super-computing Center CINECA under the ISCRA initiative.
-  D. J. Southwood, Plan. Space Sci. 22, 483 (1974).
-  S.-H. Chen, M. G. Kivelson, J. T. Gosling, R. J. Walker, and A. J. Lazarus, J. Geophys. Res. 98, 5727 (1993).
-  H. Hasegawa, M. Fujimoto, T.-D. Phan, H. Rème, A. Balogh, M. W. Dunlop, C. Hashimoto, and R. TanDokoro, Nature (London) 430, 755 (2004).
-  S. A. Pope, M. A. Balikhin, T. L. Zhang, A. O. Fedorov, M. Gedalin, and S. Barabash, Geophys. Res. Lett. 36, 7202 (2009).
-  A. Masters, N. Achilleos, C. Bertucci, M. K. Dougherty, S. J. Kanani, C. S. Arridge, H. J. McAndrews, and A. J. Coates, Plan. Space Sci. 57, 1769 (2009).
-  T. Sundberg, S. A. Boardsen, J. A. Slavin, L. G. Blomberg, and H. Korth, Plan. Space Sci. 58, 1434 (2010).
-  H. Nagano, Plan. Space Sci. 27, 881 (1979).
-  A. Miura and P. L. Pritchett, J. Geophys. Res. 87, 7431 (1982).
-  T. K. M. Nakamura, H. Hasegawa, and I. Shinohara, Phys. Plasmas 17, 042119 (2010).
-  P. Henri, F. Califano, M. Faganello, and F. Pegoraro, Phys. Plasmas 19, 072908 (2012).
-  J. U. Brackbill and D. W. Forslund, Journal of Computational Physics 46, 271 (1982).
-  O. Alexandrova, V. Carbone, P. Veltri, and L. Sorriso-Valvo, Astrophys. J. 674, 1153 (2008).
-  J. Šafránková, Z. Němeček, L. Přech, and G. N. Zastenker, Phys. Rev. Lett. 110, 025004 (2013).
-  O. Gressel, D. Elstner, U. Ziegler, and G. Rüdiger, Astron. Astrophys. 486, L35 (2008).
-  F. Del Sordo, S. Candelaresi, and A. Brandenburg, Phys. Rev. E 81, 036401 (2010).
-  S. B. Dalziel, P. F. Linden, and D. L. Youngs, J. Fluid Mech. 399, 1 (1999).
-  H. Karimabadi, V. Roytershteyn, M. Wan, W. H. Matthaeus, W. Daughton, P. Wu, M. Shay, B. Loring, J. Borovsky, E. Leonardis, S. C. Chapman, and T. K. M. Nakamura, Phys. Plasmas 20, 012303 (2013).
-  A. A. Schekochihin, S. C. Cowley, W. Dorland, G. W. Hammett, G. G. Howes, E. Quataert, and T. Tatsuno, Astrophys. J. Suppl. 182, 310 (2009).
-  S. Galtier, S. V. Nazarenko, A. C. Newell, and A. Pouquet, J. Plasma Phys. 63, 447 (2000).
-  S. Galtier, J. Plasma Phys. 72, 721 (2006).
-  P. D. Mininni and A. Pouquet, Phys. Rev. Lett. 99, 254502 (2007).
-  J. C. Perez and S. Boldyrev, Astrophys. J. Lett. 672, L61 (2008).
-  S. Boldyrev and J. C. Perez, Phys. Rev. Lett. 103, 225001 (2009).
-  J. Saur, H. Politano, A. Pouquet, and W. H. Matthaeus, Astron. Astrophys. 386, 699 (2002).
-  J. Cho and A. Lazarian, Phys. Rev. Lett. 88, 245001 (2002).
-  J. Cho and A. Lazarian, Mon. Not. Roy. Astron. Soc. 345, 325 (2003).
-  G. Kowal and A. Lazarian, Astrophys. J. 720, 742 (2010).
-  Y. Lithwick, P. Goldreich, and S. Sridhar, Astrophys. J. 655, 269 (2007).
-  A. Beresnyak and A. Lazarian, Astrophys. J. 682, 1070 (2008).
-  J. J. Podesta and A. Bhattacharjee, Astrophys. J. 718, 1151 (2010).
-  J. C. Perez and S. Boldyrev, Phys. Rev. Lett. 102, 025003 (2009).
-  P. Goldreich and S. Sridhar, Astrophys. J. 438, 763 (1995).
-  N. E. Haugen, A. Brandenburg, and W. Dobler, Phys. Rev. E 70, 016308 (2004).
-  A. Beresnyak and A. Lazarian, Astrophys. J. Lett. 722, L110 (2010).
-  J. C. Perez and S. Boldyrev, Astrophys. J. Lett. 710, L63 (2010).
-  S. Boldyrev, J. C. Perez, J. E. Borovsky, and J. J. Podesta, Astrophys. J. Lett. 741, L19 (2011).
-  O. Alexandrova, C. H. K. Chen, L. Sorriso-Valvo, T. S. Horbury, and S. D. Bale, Space Sci. Rev. 178, 101 (2013).
-  P. S. Iroshnikov, Astron. Zhurnal 40, 742 (1963).
-  R. H. Kraichnan, Phys. Fluids 8, 1385 (1965).
-  D. Biskamp and H. Welter, Physics of Fluids B 1, 1964 (1989).
-  H. Politano, A. Pouquet, and P. L. Sulem, Physics of Fluids B 1, 2330 (1989).
-  D. Biskamp and E. Schwarz, Phys. Plasmas 8, 3282 (2001).
-  J. Sommeria, J. Fluid Mech. 170, 139 (1986).
-  D. Biskamp, E. Schwarz, and J. F. Drake, Phys. Rev. Lett. 76, 1264 (1996).
-  D. Biskamp, E. Schwarz, A. Zeiler, A. Celani, and J. F. Drake, Phys. Plasmas 6, 751 (1999).
-  S. D. Bale, P. J. Kellogg, F. S. Mozer, T. S. Horbury, and H. Reme, Phys. Rev. Lett. 94, 215002 (2005).
-  G. G. Howes, J. M. Tenbarge, W. Dorland, E. Quataert, A. A. Schekochihin, R. Numata, and T. Tatsuno, Phys. Rev. Lett. 107, 035004 (2011).
-  S. Boldyrev and J. C. Perez, Astrophys. J. Lett. 758, L44 (2012).
-  S. Markidis, G. Lapenta, and R. Uddin, Mathematics and Computers in Simulation 80, 509 (2010).
-  D. Cai, L. R. O. Storey, and T. Neubert, Physics of Fluids B 2, 75 (1990).
-  P. Henri, S. S. Cerri, F. Califano, F. Pegoraro, C. Rossi, M. Faganello, O. Šebek, P. M. Trávníček, P. Hellinger, J. T. Frederiksen, A. Nordlund, S. Markidis, R. Keppens, and G. Lapenta, Phys. Plasmas 20, 102118 (2013).
-  S. S. Cerri, P. Henri, F. Califano, D. Del Sarto, M. Faganello, and F. Pegoraro, Phys. Plasmas 20, 112112 (2013).
-  C. T. Haynes, D. Burgess, and E. Camporeale, Astrophys. J. 783, 38 (2014).
-  X. Li and S. R. Habbal, J. Geophys. Res. 106, 10669 (2001).
-  F. Sahraoui, G. Belmont, and M. L. Goldstein, Astrophys. J. 748, 100 (2012).
-  M. A. Reynolds and G. Ganguli, Phys. Plasmas 5, 2504 (1998).
-  J. Egedal, Phys. Plasmas 9, 1095 (2002).
-  J. Egedal, A. Le, and W. Daughton, Phys. Plasmas 20, 061201 (2013).
-  M. E. Innocenti and G. Lapenta, Plasma Phys. Contr. Fusion 49, 521 (2007).