Properties of the Firstorder Fermi acceleration in fast magnetic reconnection driven by turbulence in collisional MHD flows
Abstract
Fast magnetic reconnection may occur in different astrophysical sources, producing flarelike emission and particle acceleration. Currently, this process is being studied as an efficient mechanism to accelerate particles via a firstorder Fermi process. In this work we analyse the acceleration rate and the energy distribution of test particles injected in threedimensional magnetohydrodynamical (MHD) domains with largescale current sheets where reconnection is made fast by the presence of turbulence. We study the dependence of the particle acceleration time with the relevant parameters of the embedded turbulence, i.e., the Alfvén speed , the injection power and scale (). We find that the acceleration time follows a powerlaw dependence with the particle kinetic energy: , with for a vast range of values of . The acceleration time decreases with the Alfvén speed (and therefore with the reconnection velocity) as expected, having an approximate dependence , with for particles reaching kinetic energies between , respectively. Furthermore, we find that the acceleration time is only weakly dependent on the and parameters of the turbulence. The particle spectrum develops a highenergy tail which can be fitted by a hard powerlaw already in the early times of the acceleration, in consistency with the results of kinetic studies of particle acceleration by magnetic reconnection in collisionless plasmas.
keywords:
acceleration of particles — magnetic reconnection — magnetohydrodynamics — methods: numerical1 Introduction
Magnetic reconnection occurs when two magnetic fluxes of opposite polarity encounter each other. Under finite magnetic resistivity conditions a current sheet is formed at the discontinuity surface, where the field lines annihilate. Direct evidence of magnetic reconnection in astrophysical and space environments like the solar corona and the Earth magnetotail indicate that in some circumstances reconnection can be very fast, with rates which are a substantial fraction of the Alfvén speed .
Fast reconnection breaks the magnetic field topology releasing magnetic energy explosively thus explaining the bursty emission, for instance, in solar flares. Relativistic particles are always observed in connection with these flares suggesting that magnetic reconnection can lead to direct particle acceleration (see e.g., the reviews de Gouveia Dal Pino, Kowal, & Lazarian, 2014; de Gouveia Dal Pino & Kowal, 2015; Uzdensky, 2011, and references therein).
In analogy to diffusive shock acceleration (DSA), in which particles confined between the upstream and downstream flows undergo a firstorder Fermi acceleration, de Gouveia dal Pino & Lazarian (2005) (hereafter GL05) proposed a similar process occurring within the current sheet where trapped particles bounce back and forth between the converging magnetic fluxes of opposite polarity in the largescale reconnection region. The particles gyrorotate around a reconnected magnetic field (see Figure 2b in Kowal, de Gouveia Dal Pino, & Lazarian, 2011), gaining energy due to collisions with magnetic irregularities at a rate (where is the reconnection speed) implying a firstorder Fermi process with an exponential energy growth after several round trips (GL05, de Gouveia Dal Pino & Kowal 2015). A similar process was also invoked by Drake et al. (2006) who investigated particles accelerated inside twodimensional contracting magnetic islands (or loops). In Kowal, de Gouveia Dal Pino, & Lazarian (2011) it has been demonstrated the equivalence between the two mechanisms for driving firstorder Fermi acceleration. This process has been extensively tested numerically mainly through twodimensional (2D) particleincell (PIC) simulations of collisionless electronion or electronpositron plasmas (e.g., Drake et al., 2006; Zenitani & Hoshino, 2001, 2007, 2008; Lyubarsky & Liverts, 2008; Drake et al., 2010; ClausenBrown & Lyutikov, 2012; Cerutti et al., 2014; Li et al., 2015), and more recently also through threedimensional (3D) PIC simulations (Sironi & Spitkovsky, 2014; Guo et al., 2015, 2016). However, these simulations can probe acceleration only at the kinetic scales of the plasma, of a few hundreds of the inertial length (, where is the plasma frequency). To assess the firstorder Fermi process in the large scales of the collisional MHD flows commonly observed in astrophysical systems, Kowal, de Gouveia Dal Pino, & Lazarian (2011, 2012) have also successfully tested it in 2D and 3D MHD simulations injecting test particles in the reconnection domain.
Currently, fast magnetic reconnection is regarded as a potentially important mechanism to accelerate particles not only in the solar system context (e.g., Drake et al., 2006, 2009; Gordovskyy, Browning, & Vekstein, 2010; Gordovskyy & Browning, 2011; Zharkova et al., 2011; Lazarian & Opher, 2009; Drake et al., 2010; Lazarian & Desiati, 2010; Li et al., 2015), but also beyond it, in galactic and extragalactic environments such as jetaccretion disk systems (e.g., de Gouveia dal Pino & Lazarian, 2005; de Gouveia Dal Pino, Piovezan, & Kadowaki, 2010; de Gouveia Dal Pino et al., 2010; Giannios, 2010; del Valle et al., 2011; Kadowaki, de Gouveia Dal Pino, & Singh, 2015; Khiali, de Gouveia Dal Pino, & del Valle, 2015; Khiali, de Gouveia Dal Pino, & Sol, 2015), pulsar winds and GRBs (e.g., Lazarian et al., 2003; Zenitani & Hoshino, 2007; Zhang & Yan, 2011; Uzdensky, 2011; ClausenBrown & Lyutikov, 2012; Cerutti et al., 2014; Sironi & Spitkovsky, 2014; Guo et al., 2014, 2015; Singh, Mizuno, & de Gouveia Dal Pino, 2016). It has been also related to the production of ultrahighenergy cosmic rays (e.g., de Gouveia Dal Pino & Lazarian, 2000, 2001; Kotera & Olinto, 2011). Besides, the accelerated particles may produce detectable nonthermal emission in a wide range of energies, specially at gamma rays (e.g., del Valle et al., 2011; Vieyro & Romero, 2012; Cerutti et al., 2014; Khiali, de Gouveia Dal Pino, & del Valle, 2015; Khiali, de Gouveia Dal Pino, & Sol, 2015; Kadowaki, de Gouveia Dal Pino, & Singh, 2015; Singh, de Gouveia Dal Pino, & Kadowaki, 2015) or neutrinos (e.g., Khiali & de Gouveia Dal Pino, 2016), therefore, studies of the acceleration rate and the particle powerlaw index are fundamental for understanding and modelling this emission.
As remarked above, in order to obtain an efficient acceleration process, reconnection has
to be fast. In collisioneless plasmas, this is usually ensured by kinetic instabilities
or by the Hall effect (in the case of an electronion plasma), both relevant only at plasma kinetic scales. In largescale collisional MHD systems, fast reconnection can be driven either by anomalous resistivity (Parker, 1979; Biskamp, 1997; Zenitani, Hesse, & Klimas, 2009) or by turbulence (Lazarian & Vishniac, 1999; Kowal et al., 2009, 2012; Lazarian et al., 2012)
In a weak turbulent medium, the wandering of the magnetic field lines allows for many simultaneous events of reconnection to happen at the same time. Moreover, the reconnected flux is more efficiently removed due to turbulence which broadens the outflow channel (see Figure 1 in Kowal et al., 2009, for example). These two factors make such reconnection fast. According to Lazarian & Vishniac (1999), , where and are the injection velocity and scale of the turbulence, respectively. It is easy to see that for the upper limit, i.e. and , the maximum reconnection rate is . Both features, the simultaneous reconnection events and the broadened reconnection layer, are very important for accelerating particles, as demonstrated in Kowal, de Gouveia Dal Pino, & Lazarian (2011, 2012).
In this work we extend the earlier numerical studies of Kowal, de Gouveia Dal Pino, & Lazarian (2011, 2012) of the acceleration of test particles in collisional,
nonrelativistic
In the next section we summarize the main aspects of the theory of firstorder Fermi acceleration within current sheets with fast reconnection driven by turbulence. In Sec. 3, we describe the numerical methodology used in this work to perform the calculations. In Sec. 4.1 we show the computed acceleration time for different models. In Sec. 4.2 we analyse the distribution of the accelerated particles. In Sec. 5 we discuss the results and draw our conclusions.
2 Acceleration model
As remarked, in this work we explore the firstorder Fermi acceleration process within largescale current sheets in collisional MHD domains with fast reconnection driven by embedded turbulence. The overall aspects of this mechanism have been thoroughly discussed in several papers and recent reviews (GL05, Kowal, de Gouveia Dal Pino, & Lazarian, 2011, 2012; de Gouveia Dal Pino & Kowal, 2015) and here we just summarize the main assumptions.
The mechanism of firstorder Fermi acceleration operating within largescale current sheets first proposed by GL05, with trapped particles moving back and forth between the two converging reconnecting magnetic flux tubes, undergoing collisions with magnetic fluctuations and suffering a net energy gain
(1) 
after each round trip, is completely general and works either in collisional or collisionless fluids in two and threedimension domains (Kowal, de Gouveia Dal Pino, & Lazarian, 2011). Furthermore, it is equivalent to the firstorder particle acceleration process within twodimensional converging magnetic islands proposed by (Drake et al., 2006), as demonstrated in Kowal, de Gouveia Dal Pino, & Lazarian (2011).
The expression (1) indicates that after several round trips the particle energy must increase exponentially. In addition, it clearly shows that reconnection has to be fast () in order to make the overall process efficient. For
instance, in the surroundings of relativistic sources
(GL05, de Gouveia Dal Pino, Piovezan, & Kadowaki, 2010; Giannios, 2010; Lazarian, 2005)
Making a simple approximation that particles would escape from the acceleration zone with a similar rate as in shock acceleration, GL05 predicted an analytical powerlaw energy distribution for the accelerated particles which is actually compatible, e.g., with observed synchrotron powerlaw spectra slope in microquasars. Relaxing the escape assumption, Drury (2012) found that the particles may leave the acceleration zone with a similar rate as in shock acceleration, which results in if the compression ratio between the outflow and inflow densities in the reconnection site is large. It is worth noting that the Lazarian & Vishniac (1999) fast reconnection is built upon an incompressible turbulent theory and the numerical simulations employed here which confirm this (Kowal et al. 2009, 2011, 2012a) were done in the nearly incompressible case, where the inflow/outflow density ratio is close to unity. Indeed, the variations of density in our computational domain are smaller than 5% and are very sensitive to the sonic Mach number. Since our sound speed is 4.0, we have very small sonic Mach number, of the order of 0.25 or smaller. Therefore, Drury’s (2012) result does not apply in our context and some anisotropy between the parallel and perpendicular components of the particle velocities is necessary in order to ensure efficient particle acceleration. In Kowal, de Gouveia Dal Pino, & Lazarian (2011, 2012), it has been found that the anisotropy arises from particles preferentially being accelerated either in the parallel or the perpendicular direction in the beginning of the process (see also discussion in de Gouveia Dal Pino, Kowal, & Lazarian, 2014; de Gouveia Dal Pino & Kowal, 2015). As in Kowal, de Gouveia Dal Pino, & Lazarian (2012), in Section 3 we depict results of numerical simulations of particle acceleration by magnetic reconnection in large current sheets in the presence of turbulence where the evolution of both particle velocity components is tracked separately in order to further address this point.
2.1 Fast reconnection driven by turbulence
In this work we focus on fast reconnection driven by turbulence as described by Lazarian & Vishniac (1999) (see also Eyink, Lazarian, & Vishniac, 2011) and tested numerically through threedimensional simulations by Kowal et al. (2009, 2012). As mentioned in Introduction, in this model, the magnetic field wandering is the key process that induces fast, independent of Ohmic resistivity, magnetic reconnection. The presence of turbulence enables the wandering of the magnetic lines and therefore, reconnection events occurring simultaneously making reconnection very fast. For extremely weak turbulence the reconnection rate reduces to the well known slow SweetParker rate. For further reading on this fast reconnection model we refer to the recent reviews by Lazarian et al. (2012, 2014, 2016).
The predicted dependence of magnetic reconnection on the properties of turbulence, i.e. on the injection power and injection scale of turbulence , is given by (Lazarian & Vishniac, 1999; Kowal et al., 2009):
(2) 
but the numerical tests by Kowal et al. (2009, 2012) indicate a weaker dependence on , i.e. . In this work we consider the estimated relation between and the turbulence to explore the dependence of the firstorder Fermi particle acceleration on the turbulent fast reconnection parameters.
2.2 Particle acceleration within current sheets with fast reconnection driven by turbulence
In a SweetParker reconnection configuration with reconnection made artificially fast by enhanced numerical resistivity, particles are accelerated through a firstorder Fermi process, as predicted in GL05 and demonstrated numerically in Kowal, de Gouveia Dal Pino, & Lazarian (2011) and Kowal, de Gouveia Dal Pino, & Lazarian (2012). When turbulence is included within the current sheet (see Kowal, de Gouveia Dal Pino, & Lazarian, 2012) the acceleration process is improved. As mentioned earlier, this is because the presence of turbulence allows the formation of a thick volume filled with multiple simultaneously reconnecting magnetic lines, making the process intrinsically threedimensional. Charged particles trapped within this volume suffer several headon scatterings with the contracting magnetic fluctuations, which significantly improves the acceleration.
The local effective accelerating electric field is , where is the particle velocity, is the local velocity of plasma fluctuations, and is the local magnetic field (see below).
In the next section we briefly describe the numerical method employed to construct the scenario above and inject test particles, as in Kowal et al. (2012a).
3 Numerical Methodology
Following Kowal, de Gouveia Dal Pino, & Lazarian (2012), we inject test particles (5,000  10,000 protons) into a frozenintime 3D MHD domain with a largescale current sheet containing embedded weakly stochastic turbulence. This domain is built by integrating numerically the collisional MHD equations, which are appropriate for the description of most macroscopic astrophysical flows, until the turbulence in the current sheet reaches a steadystate. Considering that the macroscopic MHD dynamical times are much longer than those of the test particles, using a single snapshot to accelerate the particles is a reasonable assumption (see Sec. 4.1.2).
For each test particle we numerically solve the relativistic equation of motion:
(3) 
where , , and are the particle mass, charge, Lorentz factor, and particle velocity, respectively. and are the electric and magnetic fields, respectively. The electric field is obtained from Ohm’s law:
(4) 
where is the flow velocity, is the current density and is the Ohmic resistivity coefficient. The contribution to the electric field due to the resistivity is neglected, since its effects on the particle acceleration for realistic resistivities are smaller. Therefore, the equation of motion is
(5) 
This equation is integrated in time for the particles using a 4 order RungeKutta method, and the fields are interpolated up to second order (see Kowal, de Gouveia Dal Pino, & Lazarian, 2012).
Particles are initiated with random positions and velocity directions and a thermal velocity distribution. The thermal distribution speed is fixed to for all cases. No radiative loss is included, so particles lose or gain energy only by the interactions with the plasma fluctuations.
3.1 Initial Conditions
Distinctly from Kowal, de Gouveia Dal Pino, & Lazarian (2012) who considered only one MHD model of magnetic reconnection, here we explore different conditions for the embedded turbulence in the current sheet. The models studied are described in Table 1.
All models begin with a Harris current sheet , initialized using the magnetic vector potential , with a uniform guide field . The initial density profile is set from the condition of uniform total pressure, and the initial velocity is set to zero. The reconnection is initiated by a small perturbation to the magnetic vector potential = , with and (see also Kowal et al. (2009) for more details).
The velocity and magnetic field are expressed in Alfvén speed units, defined by the antiparallel component of the magnetic field, and the initial unperturbed density . The initial strength of the antiparallel component of the magnetic field is set to in all models, and the guide field to . The distance unit is defined by the length of the computational box in the direction. The box has size and . The time is measured in units of . The computational box has a grid of for all models. For more details see Kowal et al. (2009) and Kowal et al. (2012).
Turbulence is driven using a method described in Alvelius (1999). The driving is made by injecting discrete velocity fluctuations in the Fourier space concentrated around the wave vector , corresponding to the injection scale . The amplitude of driving is solely determined by its power , the number of driven Fourier components, and the time step of driving. In all models the turbulence is subAlfvénic. See Kowal et al. (2009) for further details. When the turbulence is developed and reach a steadysate, we injected test particles in this domain and followed their trajectories according to Eq. (5). The turbulence injection power and length scale for the studied models are given in Table 1.
A map illustrating the trajectory of an injected test particle in the 3D magnetic reconnection domain is depicted in Figure 1.
Model  I  II  III  IV  V 

4 Results
4.1 Acceleration time
As mentioned above and from Eq. (2), the reconnection rate in the presence of turbulence, for a fixed value of the antiparallel magnetic field component ( ), depends only on and its injection scale (). The acceleration rate is naturally expected to depend on (de Gouveia Dal Pino & Kowal, 2015). In what follows we study the dependence of the acceleration time on , and .
For computing the mean acceleration time as a function of energy, we calculate where is the number of particles with energy between and , and is the acceleration time of particle ; the corresponding standard deviation is .
Dependence of the acceleration time on and the particle energy
For fixed different values of the turbulence injection parameters we study the particle acceleration as a function of .
In Figure
2 we show the kinetic energy evolution in logarithmic diagrams for , and
for model I. The initial particles distribution is the same in
all cases.
Particles are injected with a thermal distribution with a temperature corresponding to the sound speed of the isothermal MHD model (see Kowal, de Gouveia Dal Pino, & Lazarian, 2011).
The maps depict energies for both the parallel (red) and the perpendicular (blue) particle velocity components
In the absence of losses, the maximum energy that particles can reach
(accelerated by the first order mechanism) depends on the size of the
acceleration region (the effective thickness of the reconnection layer) and on the value of
the magnetic field. Indeed, in order to be confined the particle Larmor radius should
be smaller than the thickness of the acceleration region
(6) 
hence it changes with (or ).
If we take the thickness of the acceleration region as the width of the current sheet predicted from mass conservation, , where is the length of the reconnection layer which in our case is the size of the computational domain, then according to LazarianVishniac theory (e.g., Eyink, Lazarian, & Vishniac, 2011). However, in the 3D simulations employed here (see also, Kowal, de Gouveia Dal Pino, & Lazarian, 2012) the turbulence was injected in a region around the current sheet of thickness , therefore, the maximum energy should be constrained by rather than by and in this case . This is consistent with the values of the maximum energies inferred from Figure 1 which increase with .

Figure 3 shows the parallel and perpendicular mean acceleration times, as functions of the kinetic energy for , for the same MHD fast reconnection model (model I). The evolution of both particle velocity components, parallel and perpendicular to the magnetic field, is tracked separately, in order to detect anisotropies that may be important for particle confinement, as mentioned above. The energy is in units.
The acceleration time in Figure 2 can be fitted by a power law in energy, , where can be estimated for a given energy band as .
Figure 3 presents index as a function of the kinetic energy for the parallel and perpendicular acceleration for different values of of model I. In both cases, there is a peak in around at very low energies (). From Figure 2, we see that this peak occurs when particles are undergoing the initial slow drift acceleration (due to the spatially varying magnetic field) just before they enter the reconnection region (see also Figure 5 in Kowal, de Gouveia Dal Pino, & Lazarian 2011). When particles enter the fast (exponential) acceleration regime in the reconnection zone (Figure 2), drops quickly to values below (at energies ), therefore increasing the acceleration rate. In this regime, remains within the range from 0.2 to 0.6, until the energy reaches the saturation value and starts to increase again very smoothly as particles undergo further drift acceleration, outside the reconnection region, in agreement with Figure 2.
Figure 5 shows the averaged value of as a function of for the parallel and perpendicular components for model I. In these diagrams, in order to evaluate only in the regime of faster (firstorder Fermi) acceleration, we have selected the energy ranges for averaging starting at up to , , , and . The values of span from 0.2 to 0.6. For the indices are almost independent of the averaging energy band, as in Figure 4. For , the turbulence inside the reconnection region is in the relativistic regime and the MHD simulations considered in this work are expected to be no longer valid. This explains the increase of seen in Figure 5 for the smallest values of and thus a decrease of the acceleration rate where we expected the opposite. For around 60 and above, the index starts to strongly depend on the averaging energy band making the error bars much larger, and suggesting a trend of saturation of the value of in this low (and reconnection velocity) region. We further notice that the blue lines in these diagrams, which give the value of averaged in the low energy band (between 0.03 and 1) is dominated by the region in Figure 2 right before the entering into the zone of fast acceleration. Since the differences of the values estimated for different energy bands are still smaller than error bars, we can conclude that the characteristic values of in the firstorder Fermi regime lie between for a large range of values of , both for the parallel and perpendicular acceleration components.
Now, going back to Figure 3, we can estimate analytically the minimum acceleration time as (e.g., de Gouveia Dal Pino & Kowal, 2015):
(7) 
In Figure 3 we show the value of calculated for , and 0.1 . We clearly see that the acceleration times we obtained for the model with are longer than this lower limit which is attained only for the maximum energy in the upper right part of the diagram, as one expect since the expression above is actually valid for the maximum energy (see Eq. 6).
The above results are compatible with the notion that the acceleration time must decrease with (see also Eq. 1). To explore more quantitatively this issue, we show in Figure 6 the acceleration time as a function of for different proton energies , and . For computing in physical units we consider Alfvén time , with , the same time units as in the MHD simulations. The acceleration time decreases with for a fixed energy . Figure 6 also shows the best fitted function . For the three example energies considered here we get , with , and , respectively. The acceleration time in a DSA process scales with , where is the shock velocity; in the acceleration by reconnection it seems that there is a stronger dependence on , at least at high energies.
Dependence of the acceleration time on the environment evolution
We have also computed the acceleration time of the particles injected in three different dynamical time steps of the MHD turbulent reconnection site for model I, after the turbulence reached a statistical steady state. The time scale of the MHD environment is assumed to be much longer than the particles time scales, therefore we expect no significant changes in the particles evolution and acceleration when considering different dynamical times of the same MHD reconnection model (Kowal, de Gouveia Dal Pino, & Lazarian, 2012). From the computed acceleration times for the three simulated time steps of Model I, for both parallel and perpendicular components, we conclude that there is no significant change as expected. We have found the same behaviour for all values of (as seen in Figure 3). Also, a similar power law with an index 0.4 at relativistic energies is observed for these three dynamical times.
Dependence of the acceleration time on
The reconnection rate increases with the injection power (see Eq. 2) and therefore, the acceleration rate is also expected to increase for larger . We computed the mean acceleration time for models II and III and compared them with model I, in order to test the dependence of the acceleration efficiency on the turbulent power. Figure 7 shows mean acceleration times as functions of for , for both velocity components. We clearly see that for energies smaller than the saturation value, i.e., in the region of fast acceleration, the acceleration efficiency increases with (corresponding to shorter acceleration time , as expected from Eq. (1)), although the differences within an order of magnitude are, in general, encompassed by the error bars. For models II and III the powerlaw index in the relation is a little steeper, with . At energies larger than the saturation value (around for models II and III and for model I), the maximum acceleration times are similar in all models.
The weak dependence seen of the acceleration time with can be attributed to the fact that particles are accelerated while being scattered by the turbulent magnetic fluctuations in the plasma (). For larger injection power the velocity (and magnetic) fluctuations are stronger and therefore, the time the particles remain confined in turbulent patches gaining energy may be longer than in a case with smaller injection power .
Dependence of the acceleration time on
According to Eq. (2), the reconnection rate is also expected to increase with injection scale of the turbulence . We computed the mean acceleration time for models IV and V and compared with model I to test this dependence. Figure 8 shows the mean acceleration time as a function of for , for both, parallel and perpendicular, acceleration components with transparent areas corresponding to the estimation error of . For energies smaller than , there is no dependence on the injection scale, which could be attributed to the fixed size of the turbulent region in all models. However, for larger energies we observed the expected behaviour, where the efficiency of acceleration increases, or the acceleration time decreases with the decrease of . Although the differences are again within an order of magnitude at most (before the saturated energies are reached) and generally encompassed by the error bars. For models IV and V, the slopes of the relation in the region of fast growth are also slightly steeper than in model I, with .
Comparing Figures 7 and 8, we note that the acceleration time for a given has a slightly stronger dependence on than on , particularly at the highest energies near the saturation of the fast growth. According to Eq. (2) one might expect the opposite, i.e., a stronger dependence on . Nevertheless, we have also remarked already that the numerical simulations of Kowal et al. (2009) predict a weaker dependence of on than that predicted in Eq. (2) and thus we might expect that this would reflect in the results for the acceleration time as well.
4.2 Distribution of the accelerated particles
The protons initially have a thermal distribution, with the same thermal speed for all models. The particles begin to gain energy slowly until the Fermi mechanism starts operating. In Figure 9 the distribution of particles at 0.07 after the injection is shown for the parallel and perpendicular components of velocity; the color lines represent different values of (from to ); only accelerated particles are plotted. The initial distribution is the same for all values of . The number of particles accelerated in the direction perpendicular to the local magnetic field is larger than that of the particles accelerated in the parallel direction. For the largest higher energies have already been reached and the distributions are wider in energy. In the most efficient acceleration case (the highest ) the distribution already exhibits a powerlaw tail due to the Fermi acceleration (see Kowal, de Gouveia Dal Pino, & Lazarian, 2011, 2012).
In Figure 10 the evolution of the proton distribution during the acceleration is shown, for the parallel (red lines) and perpendicular velocity components (blue lines) to the local magnetic field. The seven time steps exhibited in the figure correspond to . We show the distributions for the 30, 90 and 200 cases of model I. The evolution of the particle distribution is similar for all values of . Initially, the particles have a normal distribution (not shown here) and, as particles gain energy, they start to populate higher values of energy. At some point (see the second curve from left to right in each diagram) the distribution becomes flatter than a normal distribution at higher energies. Due to our numerical setup, particles never stop accelerating as they are continuously reinjected into the system and therefore, the distribution shifts to higher and higher energies. As remarked, even after the particles attain the maximum energies allowed in the reconnection region by the firstorder Fermi process, they suffer further drift acceleration at a smaller rate. Thus, in the absence of radiative losses or a escape from the acceleration region, the particles may gain energy continuously. Most of the particles reach very high energies as time evolves. Their distribution may even get a positive powerlaw index in very evolved times (as we can see in the figure), but of course in real systems we should expect an interruption of this acceleration process with the escape of the particles from the finite volume of the acceleration zone and/or due to radiative losses.
It is interesting to note that for all cases shown in Figure 10 there are more particles being accelerated in the perpendicular than in the parallel direction to the local magnetic field. This anisotropy ensures the success of the acceleration process (see discussion in de Gouveia Dal Pino & Kowal, 2015, and references there in). The effective electric field accelerating a particle (with velocity u) is given by Eq. (5) and may lead to both parallel and normal acceleration directions to the local mean field depending on the original direction of the interacting magnetic fluctuation of the turbulent flow (v), or in other words on the resulting direction of (uv) with respect to the local B (given by the sum of the mean plus the fluctuating component). Nevertheless, the fact that we see in Figure 10 more particles in the perpendicular direction than in the parallel direction is an indication that the effective electric field is predominant in the normal directions on average.
During the total integrated time considered in Figure 10, , the particles in models with larger reach higher energies. Since the size of the system is independent of and is the same for all models, naturally models with larger ratios and thus larger acceleration rate will accelerate particles to higher energies (at the same time interval). We should point out that in the cases shown in this Figure, the maximum energy achievable by the firstorder mechanism (at the saturation of the fast growth) has already been reached at (see Figure 2).
As stressed above, since in our simulations the particles are continuously accelerated and there is no physical mechanism to allow them to escape, it is not possible to obtain the actual distribution of the accelerated particles. Nonetheless, we can make some estimate of the powerlaw index of the distribution soon after the particles start to populate the highenergy tail.
In Figure 11 we show the total number (not only the accelerated ones) of particles as a function of energy for two different time steps. Each figure corresponds to the cases 30, 90 and 200 of Model I. The initial normal distribution () is shown in gray dashed lines. The earliest time step plotted in each case corresponds to the approximate time when a highenergy powerlaw tail starts to form (i.e., when particles reach kinetic energies larger than , according to Figure 2); the second time corresponds to a little longer time step. In each time step we can distinguish two components in the distributions: a normal one (shown in solid gray line) + a powerlaw tail that fits the high energies starting at a certain energy we denote as . We see that the first powerlaw index at the earlier times can be fitted by to . However, these values would be a little smaller (corresponding to slightly steeper spectra) if we had taken the spectra at a little earlier time. Therefore, these values can be taken as approximate ones. The second powerlaw at later times in all the cases is flatter due to the effects discussed above and therefore, they must be taken only as illustrative of the limitations of the method. Of course, in realistic systems, the presence of physical particle escape from the acceleration zone, radiative losses and dynamical feedback of the accelerated particles into the plasma will result in steeper spectrum in the late times too (). We have also estimated the percentage of particles with energies (populating the highenergy region of the distribution). In all the cases this percentage increases from the first to the second time step.
Dependence of the particles spectrum on
For all models of Table 1, the particle distribution behaves in the same way as in model I above; as particles reach higher energies they start to develop a powerlaw. The anisotropic distributions with the production of more accelerated particles in the perpendicular than in the parallel direction are also present in models II and III (for the studied cases 50, 100, and 150).
A considerable dependence on is found in the evolution of the distribution of the accelerated particles. At each given time, model I which has the largest injection turbulence power, more particles accelerate and reach higher energies than in other models of Table 1. The differences increase with energy, in consistency with Figure 7. Model III is the less efficient as expected. Models I and II reach relativistic energies and exhibit a highenergy powerlaw tail at , while model III only shows this behaviour at a later time.
Dependence of the particles spectrum on
We also compared the particles energy distribution of Models IV and V with Model I, differing only in the turbulence injection wavenumber (according to Table 1), for the cases and 150. As in Model I, the other models also exhibit similar anisotropy with more particles being accelerated in the perpendicular than in the parallel direction to the local magnetic field. We also found a weaker dependence in the evolved spectrum with than with , as before (see Figures 7 and 8). The differences between the distributions become more relevant at the high energies near the saturation of the fast growth region, with the distribution in Model I slightly flatter than in Models IV and V, respectively, but within the error bars the distributions are very similar.
Effects of the initial particle energy distribution
In all the models studied so far, we considered the same initial thermal particle distribution with a thermal speed . This distribution corresponds to an initial energy . We have also computed the acceleration of particles for a model with the same turbulence properties of model I and with , but changing the initial particle distribution, considering two other cases: one with a larger thermal speed () and the other with a smaller one () than Model I. When comparing the acceleration rates of these three cases we find that after some time they all behave identically, i.e., the initially less energetic particle distribution takes longer time to enter the magnetic reconnection acceleration zone, but once the particles start to be accelerated, they evolve similarly to the other distributions, particularly after reaching values of the Lorentz factor .
Figure 12 depicts the particles distribution after a time, for the three cases. The best fitted total (normal + powerlaw) distribution is also shown for each model. We see that they all produce a similar powerlaw, independently of the initial thermal energy of particles.
5 Summary and Discussion
In this work we have investigated the firstorder Fermi acceleration of particles within largescale current sheets with fast magnetic reconnection driven by turbulence, using 3D collisional MHD simulations with the injection of test thermal particles, following the same approach as in (Kowal, de Gouveia Dal Pino, & Lazarian, 2012). We extended here this earlier study (see also Kowal, de Gouveia Dal Pino, & Lazarian, 2011; de Gouveia Dal Pino, Kowal, & Lazarian, 2014; de Gouveia Dal Pino & Kowal, 2015) by examining the effects of the parameters of the reconnection on the effective acceleration rate and the evolution of the spectrum of the particles. We considered models with different values of and different turbulence injection scale and power.
The main results can be summarized as follows.

The acceleration time follows a powerlaw dependence with the particle energy, , with 0.2 which is weakly sensitive to the magnetic reconnection parameters of the injected turbulence, tested for a large range of values of .

The acceleration time dependence with the Alfvén velocity is , with for particle kinetic energies between , respectively and keeping the same trend approximately for larger energies (tested for model I).

For a given value of the ratio, the acceleration time is shorter for larger values of the turbulence injection parameters, i.e. and , as expected from theory. Nonetheless, the maximum differences between the models are generally less than an order of magnitude and are within the error bars due to the uncertainties in the evaluation of the acceleration times from the numerical simulations, so that we can conclude that these dependences are not relevant.

In all the cases studied here, the number of particles being accelerated in the perpendicular direction to the local magnetic field is larger than the ones being accelerated in the parallel direction. This unbalancing is important to ensure the effectiveness of the acceleration process (see below).

The particle spectrum of the accelerated particles develops a highenergy tail, which can be fitted by a hard powerlaw index , with to (or even a little smaller) in the early times of the acceleration and is independent of the initial thermal energy of the injected particle distribution.
These results have important implications for studies of particle acceleration specially in magnetically dominated regions of astrophysical environments like the surrounds of GRBs, black holes in AGNs and microquasars, and the relativistic jets associated to these sources. As remarked, most studies of firstorder Fermi particle acceleration by magnetic reconnection have been performed considering PIC simulations (e.g., Zenitani & Hoshino, 2007; Lyubarsky & Liverts, 2008; Drake et al., 2010; Cerutti et al., 2013; Sironi & Spitkovsky, 2014; Guo et al., 2014, 2015), which apply only to the kinetic scales of the flow. In order to probe the largescale properties of the acceleration by magnetic reconnection beyond the kinetic scales in collisional astrophysical systems like those mentioned above, an MHD description is required. The 2D and 3D studies undertaken by Kowal, de Gouveia Dal Pino, & Lazarian (2011, 2012) and in this work have explored exactly these macroscopic scales of the acceleration by magnetic reconnection and thus are complimentary to the former kinetic studies.
It should be noticed also that, contrary to what has been found in PIC simulations (e.g., Guo et al., 2014), Kowal, de Gouveia Dal Pino, & Lazarian (2011) demonstrated that the acceleration of energetic particles in 2D and 3D reconnection domains shows substantial differences, being more efficient in the second case. This justifies why our analysis here has focussed on 3D geometries of reconnection only.
The earlier collisional numerical studies (Kowal, de Gouveia Dal Pino, & Lazarian 2011, 2012; see also Dmitruk & Matthaeus 2003) and most of the ones presented here have neglected the time evolution of the MHD environment. This is in general expected to be valid since this time is much longer than the particle time scales, particularly when considering a firstorder Fermi process in a statistically steady state turbulent domain. In Sec. 4.1.2, we explored the particle acceleration considering different snapshots of the reconnection domain and found no significant changes, as predicted. Nonetheless, this evolution may be important when considering more realistic nonsteady flows and when calculating the spectra and loss effects in real astrophysical systems. Preliminary steps in this direction have been performed in studies like, e.g. Lehe, Parrish, & Quataert (2009); Khiali, de Gouveia Dal Pino, & del Valle (2015); Khiali & de Gouveia Dal Pino (2016); Khiali, de Gouveia Dal Pino, & Sol (2015).
Earlier analytical studies of the firstorder Fermi process in largescale current sheets (e.g., GL05; Drury, 2012) predicted that the acceleration time would be similar to that of shock acceleration, and the energy powerlaw spectrum of the accelerated particles could be either steeper or harder than the one predicted for shock acceleration and nearly independent on the reconnection velocity. These predictions, although based on rather simplified assumptions have been at least qualitatively confirmed by the results of this work. For a broad range of reconnection velocities represented by a fiducial parametric space encompassing , the acceleration time dependence with the kinetic particle energy is found to be , with .
Furthermore, the minimum analytically estimated acceleration time according to Eq. (7) is comparable to the values found in the simulations when the particles reach the maximum energy during the firstorder Fermi acceleration in the reconnection zone (the saturation energy). As we have seen, this maximum energy is attained when the particle Larmor radius becomes comparable to the size of the acceleration zone.
It is also remarkable that the powerlaw indices obtained for the particles distribution in the highenergy tail from our collisional MHD simulations in the large scales are comparable to the values obtained from the PIC simulations in the kinetic scales of the plasma (e.g., Zenitani & Hoshino, 2001; Drake, Swisdak, & Fermo, 2013; Sironi & Spitkovsky, 2014; Guo et al., 2014, 2015; Li et al., 2015).
We should stress that the acceleration process in magnetic reconnection sites with turbulence theory depends on , and ) that determines the firstorder energy gain; (ii) the thickness of the turbulent region which improves the particle scattering probability; and (iii) the strength and maximum scale of the velocity and magnetic field fluctuations within the turbulent region, which control the scattering mean free path (or time) which in turn depend on both and . Therefore, the overall acceleration process is very complex. In this work we analysed only the dependence of with and the turbulence injection parameters. Both, and the acceleration efficiency are clearly dominated by the dependence, as one should expect for any process driving the fast reconnection, though we also obtained some weak dependencies of the acceleration time with the turbulent parameters. For instance, the scattering should happen at scales equal or smaller than , this might be the reason why only the dependence on the injection power and not on injection scale is manifested at lower energies in our results (, compare Figures 7 and 8). Moreover, at these scales particles can be scattered many times on the same side of the current sheet, with the energy gain temporarily independent of the value of until they are scattered back across the magnetic discontinuity again to complete the first order Fermi cycle.
Having the points above in mind, we should remember that the turbulence is essentially the physical mechanism that drives fast reconnection in the largescale current sheets studied here. This is a potentially very important driving mechanism because turbulence is very common in astrophysical sources and environments. Nevertheless, the first order Fermi could in principle operate in current sheets with fast reconnection driven by other possible processes and the results should not differ substantially from the present ones. This is compatible with the results found above that show only a weak dependence of the acceleration rate with the parameters of the turbulence. This may also explain why our results are similar to those of kinetic PIC simulations, where the driving mechanisms of fast reconnection are generally very distinct.
It should be also stressed that the collisional MHD simulations shown here focussed on proton acceleration. Although applicable to electrons too, the numerical integration of the electron trajectories is much longer in MHD domains with test particles. Nevertheless, such tests are also needed. Hybrid simulations combining both the PIC and the MHD approach may be a good approach to this problem (e.g., Bai et al., 2015).
We further remark that we have tried to establish a link with the results of the PIC studies which probe only the kinetic scales up to 1000 skin depth scales. But in our collisional study only the injected particles with Larmor radii near the MHD scales are effectively accelerated. This limitation can be also solved using hybrid codes able to resolve both the kinetic and the MHD scales and make a smooth transition between them (de Gouveia Dal Pino & Kowal, 2015).
Another note is in order. This work should be distinguished from studies that examined particle acceleration in pure turbulent environments (which are not embedded in largescale current sheets, see e.g., Dmitruk & Matthaeus 2003; Zharkova et al. 2011; Dalena et al. 2014; Kowal, de Gouveia Dal Pino, & Lazarian 2012; de Gouveia Dal Pino & Kowal 2015; Brunetti & Lazarian 2016). For instance, Kowal, de Gouveia Dal Pino, & Lazarian (2012) have compared the two cases and concluded that in the cases with pure turbulence particle acceleration is probably dominated by a second order Fermi process, but further studies must be carried out in order to disentangle the processes.
Finally, cosmicray acceleration investigation in magnetic reconnection sites has still many challenges to overcome, particularly in collisional MHD and relativistic regimes. The present study has tried to advance a little in the first of these topics. With regard to the second one, i.e., the study of acceleration in relativistic domains of reconnection, there has been some recent advances both in collisionless descriptions (e.g., Cerutti et al., 2013; Sironi & Spitkovsky, 2014; Guo et al., 2014, 2015, and references therein), and in collisional relativistic MHD fast reconnection involving turbulence (e.g., de Gouveia Dal Pino & Kowal 2015; Lazarian et al. 2016; Singh, Mizuno, & de Gouveia Dal Pino 2016; Takamoto, Inoue, & Lazarian 2015; see also de Gouveia Dal Pino, Kowal, & Lazarian 2014 for a short review of both approaches). These are important issues to be explored further, specially for building more realistic models of flares and variability in the spectrum of compact sources to help in the interpretation of current high energy observations and in making predictions for upcoming new generation of instruments, like the Cherenkov Telescope Array (Actis et al., 2011; Acharya et al., 2013; Sol et al., 2013) and the ASTRI CTA MiniArray (Vercellone, 2016).
Acknowledgements
M. V. d.V. acknowledges CNPq/Twas for financial support. E.M.G.D.P. acknowledges partial support from the Brazilian agencies FAPESP (grant no. 2013/105595 and CNPq (grant no. 300083/947). G.K. acknowledges support from FAPESP (grants no. 2013/040732 and 2013/188150) and PNPD/CAPES (grant no. 1475088). M. V. d.V. thanks Reinaldo SantosLima for fruitful discussions on different topics addressed in this paper. The authors also acknowledge the anonymous referee for a careful review. This work has made use of the computing facilities of the Laboratory of Astroinformatics (IAG/USP, NAT/Unicsul), whose purchase was made possible by the Brazilian agency FAPESP (grant 2009/540064) and the INCTA. M. V. d.V. thanks the great hospitality of the Instituto de Astronomia, Geofísica e Ciências Atmosféricas (São Paulo University), where part of this work was developed.
Footnotes
 pubyear: 2016
 pagerange: Properties of the Firstorder Fermi acceleration in fast magnetic reconnection driven by turbulence in collisional MHD flows–Properties of the Firstorder Fermi acceleration in fast magnetic reconnection driven by turbulence in collisional MHD flows
 Alternative descriptions of fast reconnection in a collisional MHD scenario have been proposed also by Loureiro, Schekochihin, & Cowley (2007); Shibata & Tanuma (2001).
 Recent studies (Takamoto, Inoue, & Lazarian, 2015) indicate that in relativistic domains turbulent driven magnetic reconnection behaves similarly to the nonrelativistic case (see Lazarian et al., 2016; de Gouveia Dal Pino, Kowal, & Lazarian, 2014, for reviews).
 We note that in an earlier pioneering work, Kobak & Ostrowski (2000) also studied the role of MHD turbulence in the particle acceleration process in a volume with a reconnecting magnetic field. However, they did not consider a real turbulent cascade developed selfconsistently to affect the reconnection at the current sheet, as performed in Kowal, de Gouveia Dal Pino, & Lazarian (2011, 2012) and in the present work. They instead, employed a Monte Carlo method and mimicked the effects of the turbulence with smallamplitude pitch angle scatterings. Their approach did not allow them to detect any Fermi process. Besides, the limitations of their method did not allow them to explore the dependence of the acceleration rate with the parameters of the turbulence, or the Alfvén (and reconnection) velocity, as we do in the present work.
 Giannios (2010), in particular, repeated GL05 calculations for a relativistic MHD flow and obtained a net particle energy increase , where , which in the limit recovers the expression obtained by GL05.
 We note that the ”collisional” term in our study refers only to the macroscopic MHD flow where magnetic reconnection takes place. The test particles are injected in this flow and then are accelerated due to interactions with local fluctuations through the effective electric field as described by Eq. (5).
 The number of particles being accelerated preferentially in the perpendicular or parallel direction is counted considering: and in the parallel case, and and in the perpendicular one.
 We note that the regime of faster increase in the loglog diagrams of Figure 2, actually reflects the exponential growth of the energy with time for each of the thousands of accelerated particles (Figs. 4 and 5 of Kowal, de Gouveia Dal Pino, & Lazarian 2011 show this exponential growth of energy for each test particle more clearly). This increase, which in the loglog diagrams is provided by the collective behaviour of all particle sample  can be fitted by power laws with indices > 2 (see Kowal, de Gouveia Dal Pino, & Lazarian, 2012).
 We are not considering radiative loses, but for protons they may be important under extreme conditions.
 It should be stressed that, as in our model, the Fermi acceleration and resulting particle powerlaw spectrum obtained by Guo et al. (2014, 2015); Li et al. (2015) is due to the electric field produced by the magnetic fluctuations (), while in the case of Sironi & Spitkovsky (2014), it is argued that the acceleration is dominated by the resistive electric field component, which in our case is absent (see also Kowal, de Gouveia Dal Pino, & Lazarian 2012).
References
 Acharya B. S., et al., 2013, APh, 43, 3
 Actis M., et al., 2011, ExA, 32, 193
 Alvelius K., 1999, PhFl, 11, 1880
 Bai X.N., Caprioli D., Sironi L., Spitkovsky A., 2015, ApJ, 809, 55
 Barres de Almeida U., de Gouveia Dal Pino E. M., 2014, RMxAC, 44, 122
 Biskamp D., 1997, noma.book, 392
 Brunetti G., Lazarian A., 2016, MNRAS, 458, 2584
 Cerutti B., Werner G. R., Uzdensky D. A., Begelman M. C., 2014, ApJ, 782, 104
 Cerutti B., Werner G. R., Uzdensky D. A., Begelman M. C., 2013, ApJ, 770, 147
 ClausenBrown E., Lyutikov M., 2012, MNRAS, 426, 1374
 Dalena S., Rappazzo A. F., Dmitruk P., Greco A., Matthaeus W. H., 2014, ApJ, 783, 143
 de Gouveia Dal Pino E. M., Kowal G., Lazarian A., 2014, ASPC, 488, 8
 de Gouveia dal Pino E. M., Lazarian A., 2005, A&A, 441, 845 (GL05)
 de Gouveia Dal Pino E. M., Piovezan P. P., Kadowaki L. H. S., 2010, A&A, 518, A5
 de Gouveia Dal Pino E. M., Kowal G., 2015, ASSL, 407, 373
 de Gouveia Dal Pino E. M., Lazarian A., 2001, ApJ, 560, 358
 de Gouveia Dal Pino E. M., Lazarian A., 2000, ApJ, 536, L31
 de Gouveia Dal Pino E. M., Piovezan P., Kadowaki L., Kowal G., Lazarian A., 2010, HiA, 15, 247
 del Valle M. V., Romero G. E., LuqueEscamilla P. L., Martí J., Ramón SánchezSutil J., 2011, ApJ, 738, 115
 Dmitruk P., Matthaeus W. H., 2003, ApJ, 597, 1097
 Drake J. F., Cassak P. A., Shay M. A., Swisdak M., Quataert E., 2009, ApJ, 700, L16
 Drake J. F., Opher M., Swisdak M., Chamoun J. N., 2010, ApJ, 709, 963
 Drake J. F., Swisdak M., Che H., Shay M. A., 2006, Natur, 443, 553
 Drake J. F., Swisdak M., Fermo R., 2013, ApJ, 763, L5
 Drury L. O., 2012, MNRAS, 422, 2474
 Eyink G. L., Lazarian A., Vishniac E. T., 2011, ApJ, 743, 51
 Giannios D., 2010, MNRAS, 408, L46
 Gordovskyy M., Browning P. K., Vekstein G. E., 2010, ApJ, 720, 1603
 Gordovskyy M., Browning P. K., 2011, ApJ, 729, 101
 Guo F., Li H., Daughton W., Li X., Liu Y.H., 2016, PhPl, 23, 055708
 Guo F., Li H., Daughton W., Liu Y.H., 2014, PhRvL, 113, 155005
 Guo F., Liu Y.H., Daughton W., Li H., 2015, ApJ, 806, 167
 Kadowaki L. H. S., de Gouveia Dal Pino E. M., Singh C. B., 2015, ApJ, 802, 113
 Khiali B., de Gouveia Dal Pino E. M., 2016, MNRAS, 455, 838
 Khiali B., de Gouveia Dal Pino E. M., del Valle M. V., 2015, MNRAS, 449, 34
 Khiali B., de Gouveia Dal Pino E. M., Sol H., 2015, arXiv, arXiv:1504.07592
 Kobak T., Ostrowski M., 2000, MNRAS, 317, 973
 Kotera K., Olinto A. V., 2011, ARA&A, 49, 119
 Kowal G., Lazarian A., Vishniac E. T., OtmianowskaMazur K., 2012b, NPGeo, 19, 297
 Kowal G., de Gouveia Dal Pino E. M., Lazarian A., 2011, ApJ, 735, 102
 Kowal G., de Gouveia Dal Pino E. M., Lazarian A., 2012a, PhRvL, 108, 241102
 Kowal G., Lazarian A., Vishniac E. T., OtmianowskaMazur K., 2009, ApJ, 700, 63
 Lazarian A., 2005, AIPC, 784, 42
 Lazarian A., Desiati P., 2010, ApJ, 722, 188
 Lazarian A., Eyink G., Vishniac E., Kowal G., 2014, ASPC, 488, 23
 Lazarian A., Kowal G., Takamoto M., de Gouveia Dal Pino E. M., Cho J., 2016, ASSL, 427, 409
 Lazarian A., Opher M., 2009, ApJ, 703, 8
 Lazarian A., Petrosian V., Yan H., Cho J., 2003, astro, arXiv:astroph/0301181
 Lazarian A., Vishniac E. T., 1999, ApJ, 517, 700
 Lazarian A., Vlahos L., Kowal G., Yan H., Beresnyak A., de Gouveia Dal Pino E. M., 2012, SSRv, 173, 557
 Lehe R., Parrish I. J., Quataert E., 2009, ApJ, 707, 404
 Li X., Guo F., Li H., Li G., 2015, ApJ, 811, L24
 Loureiro N. F., Schekochihin A. A., Cowley S. C., 2007, PhPl, 14, 100703
 Lyubarsky Y., Liverts M., 2008, ApJ, 682, 14361442
 Parker E. N., 1979, cmft.book,
 Shibata K., Tanuma S., 2001, EP&S, 53, 473
 Singh C. B., de Gouveia Dal Pino E. M., Kadowaki L. H. S., 2015, ApJ, 799, L20
 Singh C. B., Mizuno Y., de Gouveia Dal Pino E. M., 2016, ApJ, 824, 48
 Sironi L., Spitkovsky A., 2014, ApJ, 783, L21
 Sol H., et al., 2013, APh, 43, 215
 Takamoto M., Inoue T., Lazarian A., 2015, ApJ, 815, 16
 Uzdensky D. A., 2011, SSRv, 160, 45
 Vercellone S., 2016, EPJWC, 121, 04006
 Vieyro F. L., Romero G. E., 2012, A&A, 542, A7
 Zenitani S., Hoshino M., 2008, ApJ, 677, 530544
 Zenitani S., Hoshino M., 2007, ApJ, 670, 702
 Zenitani S., Hoshino M., 2001, ApJ, 562, L63
 Zenitani S., Hesse M., Klimas A., 2009, ApJ, 696, 1385
 Zhang B., Yan H., 2011, ApJ, 726, 90
 Zharkova V. V., et al., 2011, SSRv, 159, 357