Large-scale Compression Acceleration during Magnetic Reconnection in a Low-$\beta$ Plasma

In solar flares and other astrophysical systems, a major challenge for solving particle acceleration problem associated with magnetic reconnection is the enormous scale separation between kinetic scales and observed reconnection scale. Because of this, it has been difficult to draw any definite conclusions by just using kinetic simulations. Particle acceleration model that solves energetic particle transport equation can capture the main acceleration physics found in kinetic simulations, and thus provide a practical way to make observable predictions and directly compare model results with observations. Here we study compression particle acceleration in magnetic reconnection by solving Parker (diffusion-advection) transport equation using velocity and magnetic fields from two-dimensional high-Lundquist-number magnetohydrodynamics (MHD) simulations of a low-$\beta$ reconnection layer. We show that the compressible reconnection layer can give significant particle acceleration, leading to the formation of power-law particle energy distributions. We analyze the acceleration rate and find that the acceleration in the reconnection layer is a mixture of first order and second order Fermi processes. When including a guide field, we find the spectrum becomes steeper and both the power-law cutoff energy and maximum particle energy decrease as plasma becomes less compressible. This model produces 2D particle distribution that one can use to generate radiation map and directly compare with solar flare observations. This provides a framework to explain particle acceleration at large-scale astrophysical reconnection sites, such as solar flares.


INTRODUCTION
Energy conversion and particle acceleration in strongly magnetized plasmas are important processes that hold the key for understanding many explosive solar and astrophysical high-energy phenomena (Zweibel & Yamada 2009;Lin 2011). Magnetic reconnection is a major mechanism that drives the release of magnetic energy and nonthermal particle acceleration by reorganizing the topology and connectivity of magnetic field lines. One of the best examples for magnetic reconnection and the associated particle acceleration is solar flares. Observations have suggested that magnetic reconnection converts 10% to 50% of the magnetic energy (up to ∼ 10 33 ergs) into plasma kinetic energy within 1 − 10 minutes. During the process, a large amount of electrons in the flare region (> 10 36 electrons) are accelerated into a power-law energy spectrum f (ε) ∝ ε −s with spectral index s ∼ 3 to more than s = 9 with a medium about 5 (Lin & Hudson 1976;Krucker et al. 2010;Oka et al. 2013Oka et al. , 2015Effenberger et al. 2017). The acceleration of ions in a flare region can be as efficient as that of electrons. This is suggested by RHESSI 's observation on the correlation between electrongenerated hard X-ray flux and ion-generated γ-ray flux (Shih et al. 2009). In-situ solar energetic particle (SEP) observation has also shown that electron and ion spectra often resemble power-law distributions (Mason et al. 2012). How such efficient particle acceleration occurs over a large-scale reconnection region remains an important unsolved problem in reconnection study.
During solar flares, large-scale magnetic reconnection is in the weakly collisional (high Lundquist number) regime and is likely to have magnetic structures with a range of spatial scales. One attractive scenario emerged in the past decade is the plasmoid-dominated reconnection, where a hierarchy of plasmoids develop in a macroscopic reconnection layer (Shibata & Tanuma 2001;Loureiro et al. 2007; Bhattacharjee et al. 2009;Comisso et al. 2016) and naturally bring the current sheet from the macroscopic scale to kinetic scale (Daughton et al. 2009;Ji & Daughton 2011). It is therefore important to study particle acceleration in magnetic reconnection using a multi-scale approach. For magnetic reconnection at kinetic scales, kinetic simulations provide first-principle description for particle acceleration, but the domain size is limited due to the demanding computational expense.
The standard approach to study particle acceleration on large scales is to solve energetic particle transport equation (e.g. Parker 1965;Zank et al. 2014), but this has not been applied in reconnection study until recently (see below for a more detailed discussion). Instead, test-particle simulations have been widely used to study particle acceleration during reconnection on large scales. Below we review the previous theories and numerical simulations on particle acceleration in magnetic reconnection.
Particle-in-cell (PIC) kinetic simulation has been popular in modeling particle acceleration during magnetic reconnection, as it includes the full range of plasma physics. Previous kinetic simulations have extensively studied several acceleration mechanisms, such as direct acceleration close to the reconnection X-point (Hoshino et al. 2001;Drake et al. 2005;Fu et al. 2006;Oka et al. 2010;Egedal et al. 2012Egedal et al. , 2015Wang et al. 2016), Fermi-type acceleration in contracting magnetic islands (Drake et al. 2006;Oka et al. 2010), and acceleration in island-merging regions (Oka et al. 2010;Liu et al. 2011;Drake et al. 2013;Nalewajko et al. 2015). By summing over particle guiding-center motions, several recent studies have identified curvature drift along the motional electric field as the major particle acceleration mechanism (Dahlin et al. 2014;Guo et al. 2014Guo et al. , 2015Li et al. 2015Li et al. , 2017 in the weak guide-field case. However, because of the enormous scale separation between kinetic scales (ion skin depth ∼ 10 − 100 m) and scale of the observed reconnection region (∼ 10 7 m), it has been difficult to draw any definite conclusion and compare solar flare observations to modeling results. To overcome this major difficulty and solve particle acceleration problem in solar flare reconnection, one has to come up with a description for the acceleration of particles in macroscopic fluid scale.
Test-particle simulations are widely used in studying particle acceleration during solar flares. Both full particle orbits and the particle guiding-center motions have been calculated in background electric and magnetic fields provided by MHD simulations. Under the guiding-center approximation, one can solve particle motions in realistic scales by removing the high-frequency gyromotions. Test-particle method usually generates hard power-law energy spectra (Onofri et al. 2006;Gordovskyy et al. 2010a,b;Zhou et al. 2015Zhou et al. , 2016 that can extend to tens of keV for electrons and tens of MeV for protons but may be too hard to explain the observations (power-law index for electrons 1 < s < 2).
Acceleration due to parallel electric field is usually the dominant particle acceleration mechanism found in the these simulations. This is likely due to the large anomalous resistivity and coarse grids used in these simulations, resulting in much broader current layers and much larger resistive electric field than that in real systems. Furthermore, the large anomalous resistivity is not supported by current 3D PIC simulations of reconnection layers (Roytershteyn et al. 2012;Liu et al. 2013;Le et al. 2018). One can avoid this problem by ignoring the parallel electric field completely (Zhou et al. 2015;Birn et al. 2017), leading to particle energy spectra that are close to solar flare observations. But this method still does not take into account the effect of wave-particle interaction that scatters particles and changes the acceleration processes.
The standard approach to solve large-scale particle acceleration and transport problem is to use energetic particle transport theory, which has been widely used in studying shock acceleration and cosmic ray transport. The primary acceleration mechanism is due to adiabatic compression and is included in the Parker transport equation (Parker 1965;Blandford & Eichler 1987). Various other acceleration mechanisms (e.g. fluid shear and fluid acceleration) could also be included in the transport theory (Earl et al. 1988;Zank 2014). It is worthwhile noting that the acceleration due to curvature drift and gradient drift found to be important in earlier kinetic simulations has been included in the transport theory (Jones 1990;le Roux & Webb 2009). Several studies have attempted to develop similar transport theories (or reduced kinetic equations) for studying particle acceleration during reconnection (Drake et al. 2006(Drake et al. , 2013Egedal et al. 2013;Zank 2014;le Roux et al. 2015;Montag et al. 2017). These studies include previously studied particle acceleration mechanisms, such as parallel reconnection electric field and contracting and merging magnetic islands. While some of the studies assume that the reconnection layer is incompressible and only consider incompressible effects (e.g. Drake et al. 2006Drake et al. , 2013, other recent studies emphasized both compressible and incompressible effects le Roux et al. 2015;Montag et al. 2017). Recently, for the first time, Li et al. (2018) used fully kinetic simulations to show that compression energization dominates the acceleration of high-energy particles in reconnection with a weak guide field (< 20% of the reconnecting component), and the compression and shear effects are comparable when the guide field is moderate (∼ 0.5 times of the reconnecting magnetic field component). Meanwhile, some recent MHD simulations also suggest that the reconnection layer is compressible especially when the plasma β is low and the guide field is weak (Birn et al. 2012;Provornikova et al. 2016). These simulation results suggest that one may study particle acceleration in a large-scale solar flare reconnection site using the transport theory. Drury (2012) considered reconnection acceleration by assuming the reconnection region as a black box with a certain compression ratio r between the upstream and downstream regions. He found that compression acceleration leads to a power-law spectrum f (p) ∝ p −χ and the spectral index depends on the compression ratio in a similar way as in diffusive shock acceleration χ = −3r/(r − 1). For nonrelativistic particles, the spectral index s for energy spectrum f (ε) is related to χ by s = (χ−1)/2.
As discussed above, reconnection layer in the weakly collisional regime may have magnetic structures in various scales. It is worthwhile studying whether power-law energy spectrum can still develop and how the spectral features depend on key plasma parameters of the reconnection layer. The goal of this paper is to study large-scale compression acceleration during magnetic reconnection in the plasmoid-dominated regime.
In this paper, we solve Parker (diffusion-advection) transport equation using the background velocity and magnetic fields from high-Lundquist-number MHD simulations of a low-β reconnection layer. We assume that electrons and protons are already energetic and can interact with the background magnetic fluctuation existed in the reconnection region. In Section 2, we describe the MHD simulations and stochastic integration method for solving the Parker transport equation. In Section 3, we present our simulation results. We show that particles are significantly accelerated by the compression reconnection layer in the plasmoid-dominated regime. The acceleration leads to formation of power-law energy distribution for both electrons and protons. The power-law index, cutoff energy and the maximum energy depend on the guide-field strength and the diffusion model. This model also produces 2D particle distribution that one can use to generate radiation map and directly compare with observations. This provides a framework to explain particle acceleration at large-scale reconnection sites, such as solar flares. In Section 4, we discuss the conclusions and the implications based on our simulation results.

MHD Simulations
We carry out simulations of magnetic reconnection using the Athena MHD code (Stone et al. 2008).
We use a third-order piecewise parabolic reconstruction, the Harten-Lax-van Leer Discontinuities (HLLD) Riemann solver, the MUSCL-Hancock (VL) integrator, and the constrained transport (CT) algorithm to ensure the divergence-free state of the magnetic field. The code solves the resistive MHD equations ∂ρ ∂t + ∇ · (ρv) = 0, where ρ is the mass density, v is the velocity, e is the total energy density, B is the magnetic field, p is the gas pressure, γ (=5/3) is the adiabatic index, j is the current density, and η is the resistivity.
Unless specified otherwise, we normalize the simulations by L 0 = 5000 km (the simulation box size is 10 4 km) and v A = 1000 km s −1 , which are the typical parameters of the reconnection site of a solar flare. We assume the normalized magnetic field B 0 = 50 G and particle density is 1.2 × 10 10 cm −3 .
We choose η = 10 −5 and the same box sizes L x = L y = 2 in all simulations, resulting a Lundquist and y ∈ [0, 2]. The simulations start from two current sheets with where B 0 = 1.0 is the strength of the reconnecting magnetic field, B g is the strength of the guide field, x 1 = 0.5 and x 2 = 1.5 are the x-positions of the current sheets, and λ = 0.005 is the half-thickness of the current sheet. The grid sizes are n x × n y = 8192 × 4096, so we can resolve the initial current sheet by at least 10 cells. We employ an initial magnetic flux perturbation to speed up the reconnection onset.
where ψ 0 is the amplitude of the perturbation. Initially the total pressure (gas pressure + magnetic pressure) is uniform in the simulation box. We choose ψ 0 = 10 −4 so that the initial density variation is under 2.6%. The initial plasma β = 2p/B 2 ≈ 0.1. We choose periodic boundary conditions along both x and y directions. We perform 4 simulations with B g = 0, 0.2, 0.5 and 1.0. The initial

Solving Parker Transport Equation
We then solve Parker's transport equation where f (x i , p, t) is the particle distribution function as a function of the particle position x i , momentum p (isotropic momentum assumed), and time t; κ is the spatial diffusion coefficient tensor, v is the bulk plasma velocity, and Q is the source. Note that the particle drift v d is out of the simulation plane and is not considered here. The diffusion coefficient tensor is given by where κ and κ ⊥ are the parallel and perpendicular diffusion coefficients. κ can be calculated from the quasilinear theory (Jokipii 1971). Assuming that magnetic turbulence is well-developed and has an isotropic power spectrum P ∼ k −5/3 , the resulting κ ∼ p 4/3 when the particle gyroradius is much smaller than the correlation length of turbulence. In particular, we use the following expression for κ (Giacalone & Jokipii 1999), where v is the particle speed, L c is the turbulence correlation length, Ω 0 is the particle gyrofrequency, We assume that the correlation length L c is equal to simulation box size/30 ≈ 333 km, which is the largest eddy size in a reconnection-driven turbulence as shown by 3D MHD simulations of magnetic reconnection (Huang & Bhattacharjee 2016). We assume the average magnetic field B 0 = 50 G and σ 2 = 1. Then, κ = 1.5 × 10 14 cm 2 s −1 for 10 keV protons and 4.0 × 10 14 cm 2 s −1 for 1 keV electrons, corresponding to 0.003κ 0 and 0.008κ 0 using simulation units. Test-particle simulations have suggested that κ ⊥ /κ is about 0.02-0.04 and is nearly independent of particle energy (Giacalone & Jokipii 1999). There are also observational evidence suggesting κ ⊥ /κ can be quite large (e.g., Zhang et al. 2003;Dwyer et al. 1997). Here we examine the effect of κ ⊥ /κ by adopting three different perpendicular diffusion κ ⊥ /κ = 0.01, 0.05 and 1.0.
The Parker transport equation can be solved by integrating the stochastic differential equation corresponding to the Fokker-Planck form of the transport equation (Zhang 1999;Florinski & Pogorelov 2009;Pei et al. 2010;Kong et al. 2017). Neglecting the source term Q in Equation (3) and assuming which is equivalent to a system of stochastic differential equations (SDEs) of Ito type where σ α µ σ α ν σ = 2κ µν , dW is the normalized distributed random number with mean 0 and variance √ ∆t, and ∆t is the time step for stochastic integration. This corresponds to a Wiener process.
Numerical approximation is often-used for the Wiener process to replace the normal distribution.
We use a uniform distribution in [− √ 3, √ 3] in the code. For a two-dimensional problem, The parameters used at particle locations are calculated from v x , v y , B x , B y , ∇ · v, ∂B x /∂x, ∂B x /∂y, ∂B y /∂x, ∂B y /∂y, which are all obtained from the MHD simulations. We interpolate these parameters to the particle positions and then calculate other required parameters: κ and κ ⊥ can be functions of B x , B y and B, so ∂κ /∂x, ∂κ /∂y, ∂κ ⊥ /∂x, and ∂κ ⊥ /∂y still depend on the derivatives ∂B x /∂x, ∂B x /∂y, ∂B y /∂x, ∂B y /∂y. The detailed expressions depend on the diffusion model to choose.
In this work, we use a derivative-free Milstein method (Burrage et al. 2004) to solve the stochastic differential equation. It is different from the usual method due to one more term, which makes it a higher-order method.
where X corresponds to spatial positions x, y and particle momentum p in our simulation. f (X t , t) is the deterministic term; g(X t , t) is the probabilistic term; h is the time step; N (0, 1) indicates a normal distribution, which substituted with uniform distribution [− √ 3, √ 3] in our simulations to speed up the computation. For a 1D problem, the particle moves a distance satisfying l 2 and l x should be much smaller than the spatial variation scale of the fields. In this work, we assume ∆x 2 < ∆x 2 and choose ∆t so that l x δ x , where δ x is the grid size. For our 2D problems, we choose the following criteria to determine the time step.

Compression in a Reconnection Layer
As reconnection evolves, the current sheet becomes thinner and eventually unstable to the plasmoid instability (Loureiro et al. 2007;Bhattacharjee et al. 2009;Comisso et al. 2016). Figure 1 shows the time evolution of the out-of-plane current density j z and plasma density ρ. At t = 2.5τ A , where τ A is the Alfvén crossing time L y /v A , the current sheet just starts to break into magnetic islands (Figure 1 (a) and (d)). These magnetic islands tend to contract due to magnetic tension force and merge with each other and form larger islands (t = 7.5τ A and 10τ A ). Figure 1 (b) and (c) show that new islands are continuously generated in the unstable current sheet. During these processes, the maximum plasma density increases from 1.0 to 3.0 or higher (Figure 1 (e) and (f)). The regions with enhanced density are concentrated in magnetic islands, reconnection exhausts, and inflow regions around the top and bottom sides of the magnetic islands. Due to the mass conservation in the simulation domain, density decreases in the inflow regions close to the reconnection layer and some regions in the islands. Particles can be accelerated and decelerated when crossing these regions. We expect that the net effect will be acceleration because on average, the density increases as particles move from the inflow to the outflow regions.
The enhanced plasma density suggests that the plasma in reconnection layer is compressed. To further examine the plasma compressibility, Figure 2 shows the time evolution of the density distributions f (ρ) for different runs. Plasma density evolves to have a broad distribution from a nearly uniform value ρ 0 initially. The distributions constantly change as the simulation evolves, suggesting that the reconnection layer is very dynamic. Take the run with B g = 0 for example, ρ/ρ 0 reaches about 6 and then decreases to about 4, suggesting that the compressed plasma in the reconnection layer can expand at late stage. Due to the mass conservation, Figure 2 shows significant distribution with ρ/ρ 0 < 1. The guide field plays an important role in controlling the plasma compressibility. As B g increases, the maximum density decreases from about 6 when B g = 0.0 to 2.7 when B g = 1.0.
This result is consistent with previous MHD simulations (Birn et al. 2012;Provornikova et al. 2016).
Note that f (ρ) for B g = 0.2 is close to the case with B g = 0, indicating that a weak guide field is not dynamically important here. This is because the magnetic pressure from the guide field component is only 0.04 times that of the reconnecting component. The broad f (ρ) and the nonuniform spatial distribution of ρ indicate that not all particles can "see" the entire density transition and that particle energy spectrum might not be a simple function of the compression ratio as that predicted by diffusion-advection analysis in a planar current sheet (Drury 2012). Figure 1. The out-of-plane current density j z and plasma density ρ at t = 2.5τ A , 7.5τ A , and 10τ A for half of the simulation box (x = 1.0 − 2.0), where τ A is the Alfvén crossing time L y /v A . The initial plasma density ≈ 1.0.

Particle Acceleration due to Compression: Constant Diffusion Coefficients
The onset time for fast reconnection varies with the guide field. Since we are mostly interested in the phase when the plasmoid instability is developed, we start solving the acceleration of energetic particles by injecting pseudo particles in the simulation when strong reconnection electric field emerges. Figure 3 shows the time evolution of the maximum value of reconnection electric field |ηj z |, where η is the resistivity. |ηj z | max starts growing at different time as the guide field varies. For runs with B g = 0 and 0.2, the rise time is almost the same. For runs with higher B g , it takes longer for |ηj z | max to grow. Based on this result, we inject particles at 2τ A when B g = 0 or 0.2, at 2.5τ A when B g = 0.5, and at 5τ A when B g = 1.0. For all the simulation cases, we continue to run the simulation for 10τ A and solve the transport equation.
We now discuss the results of energetic particle acceleration. Figure 4 shows the final particle energy spectra when using constant diffusion coefficients. We show two sets of simulations, one for protons with an initial energy 10 keV and κ = κ ⊥ = 0.003κ 0 (Figure 4 (a)) and the other for electrons with an initial energy 1 keV and κ = κ ⊥ = 0.008κ 0 (Figure 4 (b)). The eventual particle energy spectra resemble power-law distributions. When the guide field is weak, the power-law distributions extend several orders of magnitude in energy. As the guide field gets stronger, the power-law spectra become steeper and shorter, indicating particle acceleration is more efficient in reconnection with a weaker guide field. The spectra are close to each other for cases with B g = 0 and B g = 0.2. This is because the compressibility of the two cases are close to each other ( Figure 2). When B g increases to 1.0, particle spectrum becomes very steep with f (ε) ∼ ε −8.45 for protons and f (ε) ∼ ε −10.1 for electrons, and the maximum energy is less than 10 times of the initial particle energy. These results show that the guide-field strength is critical for particle acceleration during magnetic reconnection. When the guide field is weak, the plasma is strongly compressed in the reconnection layer, leading to an energy spectrum harder than that of the strong guide field case. This trend for the relation between the spectral index and the compressed plasma density is in agreement with Drury (2012), except that the spectral index also has a weak dependence on the diffusion coefficient.
To examine the nature of particle acceleration in a reconnection layer, we then study how the particle acceleration rate depends on the flow speed, which is about the Alfvén speed v A in a reconnection layer. We add another three simulations with fixed κ = κ ⊥ = 1.5 × 10 14 cm 2 /s and L 0 = 5000 km but different v A from 300 km/s to 10 4 km/s for the MHD run with B g = 0. The normalized κ and κ ⊥ then change from 0.01κ 0 to 3 × 10 −4 κ 0 . For each pseudo-particle, we calculate the acceleration rate dp/dt = ∆p/∆t for each short time interval ∆t = 0.0005τ A . Then, we statistically calculate f (ε) Figure 4. Particle energy distributions for cases with constant diffusion coefficients. p is particle momentum. ε indicates particle energy and is normalized by the initial particle energy ε 0 . The dashed lines indicate power-law fittings. For (a), we assume that the particles are protons with an initial energy 10 keV and κ = κ ⊥ = 0.003κ 0 . For (b), we assume that the particles are electrons with an initial energy 1 keV and the acceleration rate for all particles in the system. Figure 5 (a) shows the distributions of dp/dt , and v A0 is 300 km/s in our normalization. The measured acceleration rate is close to zero near the injected momentum since most of injected pseudo-particles are outside of the reconnection layer in the beginning. At higher energies, the acceleration rate becomes a power-law like distribution as a function of momentum dp/dt = C(p/p 0 ) α . The acceleration rate index α is 1.06 -1.10, which is expected as particles gain energy through the compression term −p∇ · v/3 in Parker transport equation. Figure 5 (a) also shows that the acceleration rate increases when the Alfvén speed gets larger. To further study the scaling of the acceleration rate with respect to v A , we fit C as a function of v A in Figure 5 (b). We find the acceleration rate normalization km/s in our normalization, suggesting that the acceleration mechanism is a mixture of 1st-order Fermi mechanism (∝ v A ) and 2nd-order Fermi mechanism (∝ v 2 A ). This is because particle can experience both compression and expansion in the reconnection layer. But, on average, particles experience a net compression as they move into the reconnection layer, where plasma is strongly compressed as shown in Figure 1. Since the reconnection layer is dynamically evolving, the acceleration rate is time-dependent as well. Figure 5 (c) and (d) show time evolution of the acceleration rate index α and the acceleration rate normalization C. The α index fluctuates throughout the simulation. For the three cases with stronger acceleration, the power-law index fluctuates around 1.1. For the case with v A = 300 km/s, the index is larger, which is likely due to statistical errors as only a small number of particles are accelerated to high energy. Figure 5 (d) shows the acceleration rate generally decreases as the simulation evolves, which is likely because reconnection becomes saturated in the late stage.

Particle Acceleration due to Compression: Energy Dependent Diffusion Coefficients
The constant and isotropic diffusion coefficient is a simplified assumption. In reality, κ usually depends on particle momentum. According to the quasi-linear theory (Equation 5), κ ∼ p 4/3 for nonrelativistic particles propagating in magnetic turbulence with a Kolmogorov power spectrum.
The diffusion coefficient in directions parallel and perpendicular to the magnetic field can be quite different and previous test-particle calculations give a perpendicular diffusion coefficient about a few percent of the parallel diffusion. Figure 6 shows the final energy spectra when we use energy dependent κ = κ 1 (p/p 0 ) 4/3 (κ 1 = 0.008κ 0 for electrons and 0.003κ 0 for protons) with 3 different κ ⊥ /κ : κ ⊥ = κ , κ ⊥ = 0.05κ , and κ ⊥ = 0.01κ . The figure shows several trends. First, particles still develop powerlaw energy spectra, but power-law energy range is shorter and the spectra roll over at certain energies depending on the diffusion model. The maximum particle energies are lower compared with the case with constant κ because high-energy particles can escape from the acceleration regions much easier due to their larger diffusion coefficients. Second, as the ratio κ ⊥ /κ decreases, the spectra become harder and the maximum energy is higher. The spectra change dramatically for cases with B g = 1.0.
The power-law index s changes from s ∼ 8.5 to s ∼ 4 for protons and from s ∼ 12 to s ∼ 4.5.
This is because when cross-field diffusion gets smaller, particles could stay in the acceleration regions for a longer time. Third, the maximum energies get close for cases with weak or moderate guide field (B g ≤ 0.5) even though the power-law part is steeper for cases with higher guide field. Finally, in all the cases, protons can be accelerated to hundreds of keV and electrons can be accelerated to (a) dp/dt as a function of particle momentum. Note that we have normalized the simulation time t in all runs with L 0 /v A0 , which is 300 km/s in our normalization, sot = tv A0 /L 0 . We accumulate dp/dt and particle number n p in each momentum bin every 0.0005τ A from 2τ A to 12τ A and calculate dp/dt = ( dp/dt)v A /(n p v A0 ). The solid lines are simulation data and the dashed lines are the power-law fittings Cp α , where C is the acceleration rate normalization and α is the acceleration rate index.
The accelerated particles are not uniformly distributed in simulations. Figure 7 shows the spatial distributions of high-energy electrons (9 − 36 keV) for the simulation using the MHD run with B g = 0, κ = 0.008κ 0 , and κ ⊥ = 0.01κ . At an earlier time (t = 7.5τ A ), high-energy electrons are mostly in the island at y ∼ 1.4, the top side of the large island at y ∼ 0.5, and the island merging region at y ∼ 1.65, suggesting that these regions are efficient at accelerating particles. As the simulation evolves, high-energy particles are advected with reconnection outflow and also diffuse to broader regions. Close to the end of the simulation (t = 10τ A ), high-energy particles become more uniform but their distribution still peaks at the two ends of the large magnetic island and in the reconnection exhausts. This geometry is similar to the above-the-loop-top hard X-ray sources observed in solar flares Oka et al. 2015). The confinement of high-energy electrons could potentially explain hard X-ray emission observed by RHESSI. Figure 7. Spatial distributions of high-energetic particles for the MHD run without a guide field at t = 7.5τ A , 8.8τ A and 10.0τ A . Here, we assume that the particles are electrons with an initial energy ε 0 =1 keV, κ = 0.008κ 0 , and κ ⊥ /κ = 0.01. We choose particles with energy 9.0 ≤ ε/ε 0 < 36.0.

Trajectories of Pseudo-Particles
To further illustrate how particles are accelerated, Figure 8 shows a representative pseudo-particle trajectory in the case with a constant κ = 0.003κ 0 and without a guide field. We mark three red dots to indicate the three major acceleration phases, including reconnection exhaust, contracting islands, and island-merging regions. Initially, the particle slowly gets advected into the reconnection layer. It gains energy in a short period of time (7τ A < t < 8τ A ) when the particle diffuses across the reconnection current sheet, where the background plasma is highly compressed. The particle is then trapped in a magnetic island and gains more energy but the rate of energy increase becomes lower. This is because particles can lose energy when they cross expanding regions of the magnetic island ( Figure 1). In the late phase, the small island merges with the large island and the particle gets accelerated and decelerated multiple times but still gains more energy on average.

DISCUSSION AND CONCLUSION
In this work, we have studied particle acceleration in a large-scale reconnection site through solving the Parker energetic particle transport equation using velocity and magnetic fields from high-Lundquist number MHD simulations of magnetic reconnection. We found that compression in the reconnection layer leads to significant particle acceleration and the formation of power-law energy distributions for both electrons and ions. As the guide field becomes stronger, the power-law distribution gets steeper, and the energy rollover of the power-law distribution and the maximum particle energy decrease. The power-law index for electrons is about 2.4 − 13.1, depending on the guide field strength, which is close to the range found in solar flare observations Oka et al. 2018). The strong dependence of particle acceleration on the guide field may be tested in observations (e.g., Qiu et al. (2010)). When the perpendicular spatial diffusion is much smaller than the parallel diffusion, we found the maximum electron energy reaches ∼ 100 keV and the maximum proton energy reaches a few MeV. Detailed analysis shows that the acceleration rate ∝ v 1.36 A , indicating a mixture of 1st-order Fermi and 2nd-order Fermi processes.
Our simulations also generate 2D spatial distributions of energetic particles. We found the energetic particles are concentrated in reconnection exhausts and magnetic islands. If combined with a radiation model, the 2D distributions could be used to make predicted radiation map that is comparable with hard X-ray observation by RHESSI and FOXSI and microwave imaging by radio observatories such as Very Large Array (VLA) and Expanded Owens Valley Solar Array (EOVSA) (Gary et al. 2018).
Our results are consistent with Drury (2012), which shows that the spectral index depends on compressibility of the reconnection layer. But we found that the spectral index is not just a simple expression of the compression ratio between the outflow and inflow regions. This is likely due to the complex structures (e.g. magnetic islands) and multiple compression and expansion regions formed in the reconnection layer. We found in our simulations that the particle energy spectra depend on the diffusion model, especially the ratio of perpendicular diffusion coefficient and parallel diffusion coefficient. Particle diffusion processes depend on the properties of turbulence in the reconnection region such as turbulence spectrum, the turbulence amplitude, the correlation length, and the turbulence anisotropy, which are still under active research (Huang & Bhattacharjee 2016;Beresnyak 2017;Kowal et al. 2017;Loureiro & Boldyrev 2017a,b;Mallet et al. 2017;Boldyrev & Loureiro 2017;Comisso et al. 2018;Walker et al. 2018;Dong et al. 2018). We expect a better understanding of these turbulence properties and hence the particle diffusion processes in a reconnection layer in the near future.
While fluid compression is the only acceleration mechanism considered in this study, incompressible effects (e.g. fluid shear) could also accelerate particles (Drake et al. 2006;Zank et al. 2014;le Roux et al. 2015;Li et al. 2018), potentially leading to stronger particle acceleration. Quantifying how other mechanisms change the particle spectral shape and maximum energies may be important for future studies.
Our 2D simulations have a few limitations. First, the periodic boundary conditions allows the large island to grows the system size, while in a solar flare, the largest island is likely to be ejected out of the reconnection layer and cannot grow to the system size, thus the current boundary conditions might lead to stronger particle acceleration. Second, the 2D configuration prevents the field variation along the out-of-plane direction, which might affect compression energization that depends on the divergence of fluid velocity. Third, we use a plasma β = 0.1 instead of a lower plasma beta which may be present for solar flares, due to technical difficulties when doing high-Lundquist-number simulations.
Lower plasma β might lead to stronger compression and hence stronger particle acceleration.
To conclude, we find that fluid compression in a reconnection layer leads to significant particle acceleration and the formation of power-law energy distributions for both electrons and ions. The compressibility of the reconnection region, which depends on the guide field, determines the spectral index and cutoff energy of the power-law distribution, and the maximum particle energy. The diffusion coefficient and its anisotropy also influence the key features of the nonthermal particle spectra.
Our analysis shows that the acceleration in the reconnection layer is a mixture of 1st-order Fermi and 2nd-order Fermi processes. Our model includes the acceleration mechanism derived from fully kinetic PIC simulations (Li et al. 2018), and also applies to a macroscopic reconnection layer like in a solar flare. The resulting time-dependent spatial and energy distributions of energetic particles can provide explanations to observed energetic particle emissions in solar flares and other astrophysical regimes.