Particle acceleration in an MHD-scale system of multiple current sheets

We investigate particle acceleration in an MHD-scale system of multiple current sheets by performing 2D and 3D MHD simulations combined with a test particle simulation. The system is unstable for the tearing-mode instability, and magnetic islands are produced by magnetic reconnection. Due to the interaction of magnetic islands, the system turns into a turbulent state. The 2D (3D) case yields both $-5/3$ ($-11/3$ and $-7/3$) power-law spacetra for magnetic and velocity fluctuations. Particles are efficiently energized by the generated turbulence, and it forms a power-law tail with an index of $-2.2$ and $-4.2$ in the energy distribution function for the 2D and 3D case, respectively. We find more energetic particles outside magnetic islands rather than inside. We observe super-diffusion in the 2D ($\sim t^{2.27}$) and 3D ($\sim t^{1.2}$) case in the energy space of energetic particles.


INTRODUCTION
Magnetic reconnection at a current sheet is a fundamental process in plasma physics [1,2,3]. Magnetic reconnection can be characterized as a topological change of anti-parallel magnetic fields where the frozen-in condition is broken. The reconnected magnetic field drags plasma away due to the tension force.
The outflow speed roughly corresponds to the Alfvén speed. As a result of magnetic reconnection, two separated plasmas are mixed together.
It is believed that magnetic reconnection is capable of generating energetic particles [4]. Several mechanisms have been proposed so far: (1) Speiser (meandering) motion across anti-parallel magnetic fields directly accelerates particles by the inductive electric field [5], (2) particles gain energy due to the conservation of the first adiabatic at the pileup region of magnetic field [6], (3) Fermi-type acceleration occurs due to the contraction of the magnetic islands [7,8,9]. Several kinetic simulations show the existence of non-thermal particles forming a power-law tail in the energy distribution function associated with the evolution of magnetic reconnection [10,11,12,13,14,15,16].
The system of multiple current sheets has been paid attention to in recent years. There are several places forming multiple current sheets in the heliosphere. For instance, heliospheric current sheets (HCSs) [17] are usually stable in the solar wind, but they are compressed at the heliospheric termination shock and can be unstable in the heliosheath. Spacecraft observations across the sector boundaries often find multiple thin current sheets inside a HCS, and these can be interpreted as the folding of individual magnetic flux tubes [18,19,20]. Not only in the heliosphere, but pulsar winds also have a similar structure, and it is believed that the interaction of the current sheets with the pulsar termination shock produces energetic particles and is responsible for conversion of Poynting dominated outflows to the observed radiation via energetic particles produced by the interaction [21,22,23,24,25].
Once those current sheets become unstable, it is thought that the system produces several magnetic islands due to magnetic reconnection and then turns into a turbulent state. Zhang and Ma [26], Akramov and Baty [27] performed MHD simulations of double current sheets and showed that growing magnetic islands interact with each other and then the system appears to be turbulence. Gingell et al. [28], Burgess et al. [29] performed 3D hybrid kinetic simulations of multiple current sheets. The system is unstable for the tearing-mode and drift-kink instability, these instabilities lead the system to a turbulent state with a power-law spectrum with an index of −7/3 for magnetic fluctuations.
Particle acceleration among multiple magnetic islands has been proposed as an efficient acceleration process. Zank et al. [30,31], le Roux et al. [32] model a focused Parker transport equation including the effects of electric field induced by magnetic island reconnection and magnetic island contraction to understand the flux of anomalous cosmic rays observed by Voyager spacecraft, which continuously increases in the downstream of the heliospheric termination shock. The model successfully reproduces the observed flux and shows that the energy spectrum becomes harder because of acceleration by magnetic islands in the downstream of a shock wave [31,33]. This has been also observed at interplanetary shock waves at 5 au [34,35].
However, recent kinetic simulations of multiple current sheets for a non-relativistic plasma did not show very efficient particle acceleration as expected by models. Drake et al. [36] performed 2D full PIC simulations of multiple current sheets and observed particle energization by a few decades in energy, but a power-law energy distribution did not form. 3D hybrid kinetic simulations was done by Burgess et al. [29], and apparent particle acceleration of ions and pickup ions was not found. Nakanotani et al. [37] investigated the interaction of current sheets with a shock wave and found ion flux increase associated with the evolution of the tearing-mode instability of current sheets in the downstream of the shock wave.
However, the power-law index of the energy spectrum unchanged associated with the generation of multiple islands due to the tearing-mode instability. Note that particle acceleration in multiple current sheets of a relativistic electron-positron plasma has been shown efficient [38].
One question arises: how efficient is particle acceleration on a larger scale, say, MHD scale? Recently, Arnold et al. [16] show that electrons are efficiently accelerated by Fermi acceleration due to the coalescence of magnetic islands using MHD simulations combined with guiding-center approximated electrons and including kinetic effects of energetic electrons. They pointed out that standard PIC simulations yield only a short power-law tail which extends a decade in energy because of the limitation of the simulation size. Therefore, this can be a reason why particle acceleration in the previous studies of multiple current sheets is not efficient as expected. We pursue to answer if particle acceleration on a larger scale is efficient or not.
In this study, we combine MHD simulations and test particle simulations to investigate particle acceleration. This method has been used for several investigations of particle acceleration in magnetic reconnection and turbulence for non-relativistic [39,40,41,42] and relativistic particles [43,44]. Although they ignore a feedback effect from energetic particles to MHD simulations, they provide valuable insights for particle acceleration in the MHD scale, which can not be easily obtained in kinetic simulations due to computational limitation. A similar idea has been applied for test-particle electrons in hybrid kinetic simulations [45,46]. We perform 2D and 3D MHD simulations of multiple current sheets combined with test particle simulations. This paper is organized as follows. In Section 2, we describe the scheme of an MHD simulation combined with a test particle simulation and initial conditions. Section 3 shows results of 2D and 3D simulations that present the evolution of multiple current sheets, particle acceleration, and particle diffusion in the energy space. The last section provides the discussion and conclusion to derive that particle acceleration in MHD-scale multiple current sheets is indeed efficient in both 2D and 3D systems.

Frontiers
We combine an MHD simulation with a test particle simulation to investigate particle acceleration in a system of multiple current sheets. We solve the following compressible ideal-MHD equations, (1) ∂ t e + ∇ · (hV + E × B) = 0; (3) where ρ is the plasma density, V plasma velocities, P plasma pressure, B magnetic fields, E electric fields, η artifitical magnetic resistivity, and J current densities. γ is an adiabatic index, and we set γ = 5/3.
We use an MHD scheme proposed by Kawai [47]. The first spatial derivative is calculated by the sixthorder compact scheme, and time integration is done by the third-order total variation diminishing (TVD) Runge-Kutta scheme [48]. The artifitial magnetic resistivity has the following form [47], where C η is an dimensionless and arbitrary parameter, c s the local sound speed, χ l referes to the Cartesian coordinates in the l-direction, and ∆χ l is the lobal grid spacing in the l-direction. Here, we set ∆ = ∆x 2 + ∆y 2 + ∆z 2 , which gives larger amount of magnetic resistivity compared to the form in [47].
The overbar denotes an approximate truncated Gaussian filter [49]. We use a fourth-order explicit scheme [50] for the fourth derivative. The magnetic resistivity with this form automatically localizes in regions where the current density has a strong gradient, such as current sheets. Therefore, the resistivity tends to damp turbulence less than a constant magnetic resistivity. We also introduce an artificial bulk viscosity and mass diffusivity to capture a shock wave and contact discontinuity correctly [47]. We note that the divergence-free condition (∇ · B = 0) is satisfied at around machine accuracy (∼ 10 −13 ) since we use a central-type finite difference scheme [51,47].
We set multiple current sheets in a periodic box. We assume the force-free condition for the current sheets [52,53,54], where B 0 is the in-plane magnetic field, d is the distant between two neighboring current sheets, L 0 the half thickness of a current sheet, and B g the background magnetic field. The plasma density and pressure are set to be uniform. We add small fluctuations δA z in the z−componet of the vector potential to initiate magnetic reconnection at current sheets for 2D and 3D simulations, where δA 0 is a constant value, and φ 2D and φ 3D are random phases for 2D and 3D simulations, respectively.
We use δA 0 = 0.05 and 0.02 for the 2D and 3D case, respectively. We confirmed that the overall evolution of the current sheets was similar as uniform random fluctuations were used and, therefore, it does not depend on the choice of initial fluctuations.
Simulation parameters used in the MHD simulations are as follows. We use L 0 as the unit length of the simulation and the alfvén speed v A0 defined by B = B 2 0 + B 2 g as the unit speed so that L 0 = 1 and v A0 = 1. We also set the uniform plasma density to ρ 0 = 1. The size of the simulation box is L with the grid number N x × N y × N z = 1024 × 256 × 256 for 2D and 3D simulations, respectively. The total plasma beta (β = β i + β e ) corresponds to 1. Here, β i and β e are the ion and electron plasma beta, respectively. We set the parameter C η = 2 for both 2D and 3D simulations. We put 4 current sheets in the box (d = 10L 0 ). The Courant-Friedrichs-Lewy (CFL) number is 0.5 and 0.25 for 2D and 3D simulations, respectively. In this study, we only consider cases without the background magnetic field (B g = 0).

Frontiers
At the same time, we solve the following equation of motion in a normalized form for non-relativistic particles using the standard Buneman-Boris method, Here, α = T 0 Ω c where T 0 is the characteristic time scale of the MHD simulation and Ω c is the cyclotron frequency of particles. The parameter α is an arbitrary and user-specified parameter since the system of the ideal MHD is scale-free, and we set α = 500. The same normalization used in the MHD simulation is applied in the equation of motion, and the particle energy is normalized by where m p is the particle mass. The total number of particles is N p = 5 * 1024 * 256 and 1024 * 256 * 256 for 2D and 3D simulations, respectively. We distribute particles uniformly in space, and they have the Maxwell distribution for velocities with a temperature of T p = 0.25. Here, we assume the equal temperature between ions and electrons. We introduce sub-cycles when calculating the equation of motion with a time step of ∆t p = ∆t M HD /250 where ∆t M HD is the time step calculated in the MHD simulation since the MHD time step can be larger than the cyclotron period. Although we do not have to specify if the test particle simulation in the MHD simulation is for electrons or ions, the parameter α = 500 can be appropriate for ions rather than electrons since α may become much larger for electrons in the scale of our interest [41].

2D case
Multiple current sheets evolve into a turbulent state. Fig. 1 shows the time evolution of the current density J z from t = 0 to 300. The black lines show the magnetic field lines There are 4 current sheets located equidistant at the initial time. The added initial fluctuations onset the tearing-mode instability, and we can see that magnetic reconnection occurs in the current sheets at t = 50. Since the phase of the fluctuations is random, the location of magnetic reconnection is also random. As the simulation proceeds, magnetic islands produced by the magnetic reconnection grow in size and merge with each other in the same current sheet. When the size of magnetic islands is roughly equal to or larger than the initial current sheet distance (10L 0 ), magnetic islands start interacting (t = 150). We observe that regions outside the magnetic islands become turbulent. At the later time (t = 300), the size of merging islands becomes around 20L 0 , and the system turns into turbulence.
The turbulence exhibits a −5/3 power-law in the magnetic and velocity fluctuations. confirm that the later stage of the system is in a turbulent state.
Non-thermal particles are produced during the evolution of the multiple current sheets into turbulence.
where N κ is the number of particles, k B the Boltzmann constant, T κ the kappa temperature, Γ the The location of energetic particles depend on the stage of the evolution of multiple current sheets. Fig. 4 shows the time evolution of the energy density defined by, where E min is the minimum energy and set to E min = 4, so that we count only energetic particles. These panels correspond to different times, t = 0, 50, 100, 150, 300 from top to bottom. The white lines are the magnetic fields lines. Note that the color scales are different at each time. It is obvious that there is no energetic particles at the initial time. After the onset of magnetic reconnection (t = 50), there are some Frontiers energetic particles produced along a current sheet. This acceleration is typical for magnetic reconnection [8,16]. As magnetic islands grow in size, we can see that energetic particles are trapped inside of magnetic islands. At t = 150, when magnetic islands interact with each other, it seems that energetic particles locate among magnetic islands rather than trapped inside. This is more evident at the end of the simulation (t = 300), the energy density outside of magnetic islands is much higher than the inside. Therefore, we can conclude that energetic particles are firstly accelerated inside current sheets and trapped inside of magnetic islands, then released and further accelerated after magnetic islands start interacting with each other. Note that Hoshino [38] also observed that energetic particles locate outside of magnetic islands in the simulation of multiple current sheets.
Particles are efficiently accelerated by turbulence. In Fig. 5, the left-top panel shows the time evolution of the energy of a typically accelerated particle. The shaded regions denoted by (a)-(c) correspond to the other panels in Fig. 5 The color scale in the panels (a)-(c) represents the particle energy. There are three major acceleration events, the first one is at t = 130 and the acceleration is a quick energization. As seen in the panel (a), the particle is accelerated by a reconnection outflow of a single current sheet. When the particle enters a current sheet, it is kicked and moves along the outflow. The second (t = 155) and third (t = 270) accelerations are formally similar and accelerated by turbulence. As mentioned early, the turbulence is produced by the interaction of magnetic islands, and it starts from T ∼ 150. The motion of the particle appears stochastic in the panels (b) and (c), and the acceleration time is gradual compared to the first acceleration. The slopes of the two acceleration times are consistent. The particle energy finally reaches E = 60. The particle trajectory indicates that, at first, a particle is energized in a single current sheet and then is further accelerated by turbulence produced due to the interaction of magnetic islands.
The diffusion of energetic particles in the energy space is super-diffusive. Fig. 6 shows the mean square displacement (MSD) of energy < ∆E 2 > of energetic particles [57,58]. We only consider particles whose energy is larger than E = 4 since the motion of lower-energy particles may significantly change the MSD [58]. The definition of < ∆E 2 > is as follows, where N p is the number of energetic particles (E > 4). Here, ∆E(t) is the displacement of a particle energy, ∆E(t) = E(t) − E(0) where E(0) is the initial particle energy. However, since the number of energetc particles are a few until t = 20 and ∼ 10 4 at T ∼ 150 (not shown here), we only consider the time after t = 150. The MSD of energy can be fitted by a power-law < ∆E 2 >∝ t a E with a power-law index of a E = 2.27. This indicates that the energy transport is super-diffusive. Note that the index a r < 1 is sub-diffusion and a r = 1 is normal diffusion.

3D case
Multiple current sheets in the 3D simulation box become turbulent via the tearing-mode instability.
The 3D simulation has the same condition as the 2D simulation but it extends the simulation box to the z−direction by 40L 0 and we use a smaller value of the CFL number (c CF L = 0.25). Fig. 7 shows the snapshots of the current density J z at the different times t = 0, 100, 150, 200. As the same in the 2D simulation, four current sheets are located inside of the simulation parallel to the x − z plane at t = 0. Small fluctuations seen at t = 0 are because of the initial fluctuations defined by Eq. 14. The initial fluctuations initiate the tearing-mode instability, and magnetic reconnection proceeds at the current sheets. The location of the magnetic reconnection is also random as well as the 2D simulation. Although magnetic islands grow in size after the onset of magnetic reconnection, the shape of magnetic islands is not as clear as magnetic islands in the 2D simulation. This is because the magnetic reconnection occurs at random on the current sheets and magnetic islands merge with each other in the 3D simulation. Therefore, the evolution of current sheets is much more complicated in the 3D simulation. Due to the interaction of the destabilized current sheets, the system appears turbulent at t = 150. At the end of the simulation (t = 200), small scale fluctuations are more visible than at t = 150, and the system transits to a highly-turbulent system. The turbulence appears isotropic since there is no background magnetic field.
The spectra of the magnetic and velocity fluctuations form a −11/3 and −7/3 power-law, respectively. Non-thermal particles are produced during the evolution of the multiple current sheets and form a powerlaw tail. Fig. 9 shows the energy distribution of test particles at different times corresponding to the color scale. After the onset of the magnetic reconnection, the existence of non-thermal particles is not as obvious as in the 2D simulation. The particles seem to be heated rather than accelerated. However, a power-law tail starts forming after the turbulence begins to be created (t ∼ 125). Super-diffusion of energetic particles is observed in the energy space. Fig. 10 shows the time evolution of the energy MSD of energetic particles. We consider only particles whose energy is larger than 4. The number of particles is a few until t = 50. Therefore, we focus on the later time. After magnetic islands start interacting with each other (t = 125), the MSD is fitted by ∝ t 1.2 . This indicates that the particle acceleration in the evolution of the turbulence is super-diffusive.

DISCUSSION AND CONCLUSION
Although the evolution of multiple current sheets is different in the 2D and 3D simulations, both cases yield turbulence at the end of the simulation. In the 2D case, current sheets are unstable for the tearingmode instabiliy, and magnetic islands are produced by magnetic reconnection. In the 3D case, magnetic reconnection occurs at random on current sheets (the x − z plane) and the evolution of magnetic islands differs along the z−direction. This makes the evolution of current sheets more complicated in the 3D case than in the 2D case. However, the system of both cases develops into a highly turbulent state at the end of the simulation.
The efficiency of particle acceleration in the 2D simulation is higher than that in the 3D simulation. While the power-law index of the energy distribution in the 2D case exhibits −2.2, it is −4.2 in the 3D case. This simply implies that the acceleration in the 2D case is more efficient than in the 3D case. The maximum energy of the 2D case (E max ∼ 300E 0 ) is also higher than that of the 3D case (E max ∼ 100E 0 ). We note that the power-law tails extend to the maximum energies. The index of the observed super-diffusion in the 2D case (2.27) is higher than that in the 3D case (1.2). This also indicates that the 2D acceleration is more efficient than the 3D acceleration. We interpret this because in the 2D case particles can be easily trapped in the turbulence more than in the 3D case.
Compared to previous studies, an MHD scale system of multiple current sheets is an efficient acceleration site. Former kinetic simulations [36,29,37] did not show a significant particle acceleration, for example, (1) no power-law tail and (2) acceleration by a factor of a few decades. However, as we have shown, the particle acceleration in the both 2D and 3D cases forms a power-law tail and is accelerated by a factor of more than 100. Therefore, we conclude that particle acceleration in an MHD-scale system of multiple current sheets is efficient. Although it is possible to directly verify this by extending kinetic simulations to MHD-scale, it may not be realistic due to the current computational power.
The particle acceleration observed in the 2D and 3D cases can be modeled by a fractional Fokker-Planck model, which is a generalization of a classical Fokker-Planck model. It is thought a super-diffusion in the energy space is an indication of efficient particle acceleration and can be related to the formation of a power-law tail [59,60,61,58]. There are several models for anomalous diffusions in the energy space as well as the real space using a fractional Fokker-Planck model to understand particle acceleration with the nature of anomalous diffusion often observed in space plasmas [62,59,63,60,64]. In a future study, we use a fractional Fokker-Planck model and compare it with several simulations by varying the background magnetic field.
We do not expect a plasma beta dependence on particle acceleration. Following the fact that plasma beta does not strongly affect the tearing-mode instability [65], We assume that multiple current sheets develop into a turbulent state in several values of plasma beta. Since the structure of magnetic reconnection appears to be turbulent in a low-beta plasma [66,67], we anticipate that particles can be still efficiently accelerated by the turbulence as well as we have seen in our simulations.
On the other hand, we expect some dependencies of particle acceleration on the strength of the background magnetic field. Several studies of magnetic reconnection show that particle acceleration becomes less efficient as the background magnetic field becomes strong [68,69,70,16]. This can be because particle motion for Fermi acceleration is limited by the background magnetic field. In the system of multiple current sheets with a strong magnetic field, an initial acceleration by a single magnetic reconnection becomes less efficient, and therefore the latter acceleration phase due to turbulence can be less efficient as well.
In conclusion, We have performed 2D and 3D MHD simulations of multiple current sheets combined with test particle simulations to investigate particle acceleration. In the both cases, multiple current sheets are unstable for the tearing-mode instability and transit to a turbulence state with a power-law spectra for magnetic and velocity fluctuations. We also observe the formulation of magnetic islands because of magnetic reconnection during the transition. Non-thermal particles are efficiently produced due to turbulence generated by the interaction of magnetic islands. Their energy distribution can be fitted by a Kappa distribution with a Kappa index (or power-law index) of −2.2 and −4.2 for the 2D and 3D case, respectively. The efficient acceleration is consistent with the observed super-diffusion in the energy space for the both cases, which can be modeled by a fractional Fokker-Planck model.

CONFLICT OF INTEREST STATEMENT
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

AUTHOR CONTRIBUTIONS
The Author Contributions section is mandatory for all articles, including articles by sole authors. If an appropriate statement is not provided on submission, a standard one will be inserted during the production process. The Author Contributions statement must describe the contributions of individual authors referred to by their initials and, in doing so, all authors agree to be accountable for the content of the work. Please see here for full authorship criteria.

FUNDING
Details of all funding sources should be provided, including grant numbers if applicable. Please ensure to add all necessary funding information, as after publication this is no longer possible.

ACKNOWLEDGMENTS
This is a short text to acknowledge the contributions of specific colleagues, institutions, or agencies that aided the efforts of the authors.

DATA AVAILABILITY STATEMENT
The simulation data shown in this paper are made available upon request by the corresponding author. Figure 1. Snapshots of the current density J z in the 2D case at different times, t = 0, 50, 100, 150, 300 from top to bottom. Black lines represent the magnetic field lines (contour lines for the vector potential A z ).