Reconnection and particle acceleration in interacting flux ropes – II. 3D effects on test particles in magnetically dominated plasmas

Abstract

1 INTRODUCTION

Magnetic reconnection causes a magnetic field to rapidly and violently rearrange its topology. This topological change affects plasma energetics and is one of the processes controlling energy exchange between different plasma system constituents (Kulsrud 1998; Priest & Forbes 2000). The main process of interest here is the conversion of energy available in the magnetic field into non-thermal particle distributions in magnetically dominated plasmas (low plasma-β). This phenomenon causes violent energy releases in a wide range of astrophysical events, including various kinds of flares and bursts of high energy (UV, X-ray and gamma-ray) (Uzdensky 2016). Solar flares are the most prominent and well-studied classical examples of reconnection (Masuda et al. 1994; Krucker et al. 2010). Observations reveal that 10–50 per cent of magnetic energy is converted into energetic charged particles (Lin & Hudson 1976) and that particles develop a power-law energy distribution containing energy of the same order as the converted magnetic energy (Krucker et al. 2010; Oka et al. 2015). In some observations of solar flares, the emission has no distinguishable thermal part and almost all electrons are accelerated to non-thermal energies (Krucker et al. 2010; Krucker & Battagli 2014). Electrons in solar flares can reach energies up to 5–50 MeV, while protons gain energies up to several GeV (Evenson et al. 1984; Aschwanden 2002). Reconnection has also been proposed as a mechanism for powerful flares and particle acceleration in more extreme settings like pulsar wind nebulae (Michel 1982; Lyubarsky & Kirk 2001; Porth et al. 2016), gamma-ray bursts (McKinney & Uzdensky 2012), magnetospheres of magnetars (Lyutikov 2003, 2006; Meng et al. 2014; Elenbaas et al. 2016) and in coronae and jets of accreting black holes and active galactic nuclei (Goodman & Uzdensky 2008; de Gouveia Dal Pino, Piovezan & Kadowaki 2010; Wilkins et al. 2015).

Flares from astrophysical objects require energy from macroscopic scales to be transferred to the microscopic scales on which particles are accelerated. The change of topology of the magnetic field configuration on large scales is well described by magnetohydrodynamics (MHD) and the dissipation at small scales is described by kinetic (particle) theory. The energetics of the plasma can be split into the part relevant at the fluid level plus non-thermal particle distributions. The MHD approach covers the overall scales and energetics of the system but does not give any information on particle dynamics. Kinetic approaches fully describe the microscopic scale, but are too costly to cover full astrophysical systems. In this work, we treat electrons and ions as test particles embedded in a thermal (MHD) plasma. Particle acceleration associated with reconnection and shocks in magnetically dominated Newtonian plasmas in the solar corona is studied extensively with test particle approaches (Rosdahl & Galsgaard 2009; Gordovskyy, Browning & Vekstein 2010; Gordovskyy & Browning 2011a,b; Gordovskyy et al. 2014; Zhou et al. 2015; Pinto et al. 2016; Zhou et al. 2016; Ripperda et al. 2017; Threlfall, Neukirch & Parnell 2017) and even in relativistic plasmas in the context of pulsar wind nebulae (Porth et al. 2016; Giacchè & Kirk 2017). The test particles are guided by MHD fields without giving feedback to these fields.

To initiate reconnection in the MHD plasma, we perturb an equilibrium of two adjacent, antiparallel and repelling current channels (Richard, Sydora & Ashour-Abdalla 1990; Keppens, Porth & Xia 2014; Ripperda et al. 2017). Translation and rotation of the currents cause a plasma disruption. Ripperda et al. (2017) showed that reconnection occurring due to this tilt instability is an efficient source of highly energetic particles in 2.5D settings. In 3D configurations, the kink instability interacts with the tilt instability, redistributing the poloidal magnetic field. Current channels undergoing a tilt or tilt–kink instability typically show a growth phase on Alfvénic time-scales in which kinetic energy grows exponentially. This energy is expected to be released and transferred to charged particles via magnetic reconnection. In the stellar corona context, the kink instability is one of the well-explored routes to initiate flares. Accessing a kink instability via antiparallel, repelling flux ropes is intimately connected to the coalescence instability of two attracting and merging flux ropes. The current filaments eventually attracting or repelling typically form as a result of turbulent plasma processes or instabilities as the tearing mode. Both the tilt and coalescence instability have been studied extensively in 2D configurations (see e.g. Longcope & Strauss 1993; Strauss & Longcope 1998; Marliani & Strauss 1999; Ng et al. 2008). Here we propose that such antiparallel currents may also form in stellar coronae where they play a role in magnetic island interactions and can cause particle acceleration. The evolution that is observed shows that repelling currents lead to sudden release of magnetic energy and current dissipation. The repelling islands typically cause localized reconnection, thin current sheet development and strongly curved magnetic fields. On a large scale, the current interaction can lead to rapid transfer of significant amounts of energy from the currents to the particles. Our model is representative for the top parts of adjacent flux ropes as seen in extreme ultraviolet observations of the solar corona. If two of these loops develop antiparallel currents, the tilt instability route to reconnection is accessible (Keppens et al. 2014). This model was studied in 2D by Richard et al. (1990) and extended to 3D by Keppens et al. (2014) and Ripperda et al. (2017). Here we investigate the effect of an additional kink instability on particle acceleration as well as the effects of boundary conditions, initial conditions and resistivity models on particle dynamics during the full 3D MHD evolution. We use high-resolution MHD results of Ripperda et al. (2017) for the evolution of the repelling current channels. For test particle simulations, we make use of the latest addition to the mpi-amrvac code (Porth et al. 2014) to dynamically evolve test particle populations during MHD evolution. We apply the guiding centre approximation in which particle gyration is neglected and only the particle velocity parallel to the magnetic field is evolved. This approach is valid in typical non-relativistic, magnetically dominated plasmas where the particle energy density is less than the energy density of the underlying MHD fluid. The numerical methods employed are described in detail in Section 2. In Section 3, 2.5D test particle simulations are discussed, and in Section 4 we discuss 3D simulations and the effect of the kink instability on test particle acceleration and energetics. We compare the results of the guiding centre approximation to solutions obtained from the full particle equation of motion in Appendix A3.

2 NUMERICAL SETUP

From solving the set of resistive, compressible MHD equations we obtain the magnetic and electric fields |$\boldsymbol {B}$| and |$\boldsymbol {E}$| and the total current density |$\boldsymbol {J} = \nabla \times \boldsymbol {B}$|⁠. The MHD data are scaled to CGS units before being used in the test particle calculations. To analyse the energetics and acceleration of electrons and protons in the MHD simulations, we follow the orbits of test particles in the MHD flow, similar to the approach of Porth et al. (2016) and Ripperda et al. (2017). In 3D simulations particles are injected uniformly in the domain. To analyse the behaviour of particles in the reconnection zones specifically, several runs are performed with a fraction of 0.99 of the ensemble uniformly distributed in space in a rectangular block, encapsulating the two (displaced) current channels and the areas with the largest current density, x ∈ [−1L, 1L], y ∈ [−2L, 2L], z ∈ [−3L, 3L]. The other fraction of 0.01 of the ensemble is uniformly distributed over the full domain x ∈ [−3L, 3L], y ∈ [−3L, 3L], z ∈ [−3L, 3L], including the surrounding background. In 2.5D simulations, with translational invariance in the z-direction for MHD evolution, the particles are distributed in the (x, y)-plane in accordance with the 3D simulations, at z = 0. In case I3De (see Table 1) the current channels are twice as long as in the reference cases, meaning that the particles are distributed over the domain x ∈ [−3L, 3L], y ∈ [−3L, 3L], z ∈ [−6L, 6L] and the MHD resolution is 300 × 300 × 600.

Equation (4) is advanced with a second-order symplectic Boris scheme (e.g. Birdsall & Langdon 1991). Each particle is advanced with an adaptive, individual time-step, ensuring that a single gyration is resolved by at least 60 steps, to ensure numerical stability. Equations (5)–(7) are advanced with a fourth-order Runge-Kutta scheme with adaptive time stepping. Here, the particle time-step δt is determined based on its parallel acceleration a = dv_‖/dt and velocity |$v = \sqrt{(v_{\Vert })^2 + (v_{\perp })^2}$| as the minimum of δr/v and v/a, where δr is the grid step. This grid step is restricted such that a particle cannot cross more than one cell of the MHD grid in one time-step. The fields |$\boldsymbol {E}$| and |$\boldsymbol {B}$|⁠, and for the GCA equations their spatial derivatives, are obtained at the particles position via linear interpolations in space and time between the fluid steps limited by the Courant-Friedrichs-Lewy, condition. The particles gyroradius is also calculated at every time-step and compared to the typical cell size to monitor the validity of the guiding centre approximation.

In the (x, y)-plane we employ open boundary conditions, in which the particles leaving the physical domain are destroyed. In 2.5D simulations we limit the length of the flux ropes in the z-direction, which is invariant for MHD fields. A particle crossing an artificially set boundary, at z = 3L or z = −3L (consistent with the z-boundaries in 3D simulations), is destroyed. For each destroyed particle, a new particle is injected at the opposite z-boundary with a thermal velocity from a Maxwellian distribution. This thermal bath boundary condition limits the length of the flux rope to 6L and counteracts particles accelerating indefinitely in the invariant z-direction. In 3D configurations, we have periodic boundary conditions for MHD fields, where particles leaving a z-boundary are periodically injected at the opposite z-boundary, consistent with MHD. In specific cases (see Table 1), we employ a similar boundary condition as in 2.5D configurations where particles leaving a z-boundary are destroyed. For each destroyed particle, a new particle is injected at the opposite boundary with a thermal velocity according to a Maxwellian. In this way, a large enough ensemble of particles is retained at all times, to achieve accurate statistics. Consequently, the length of the flux rope is limited to 6L or equivalently ∼0.1 solar radii. This boundary condition realistically mimics the injection of thermal particles in a (curved) flux rope in the corona of a star.

In Table 1, we list all particle simulations including the dimension, the type of particles (indicated by e⁻ for electrons and p⁺ for protons), the number of particles N_tot, the initial spatial distribution (which is always uniform in the domain given) and the equations of motion solved for the particles (guiding centre approximation, indicated with GCA or full equations of motion indicated with Lorentz). We indicate the typical length l the particles can travel in the current channels in the z-direction. Infinite length corresponds to periodic boundary conditions. A finite value corresponds to distance between the two opposite z-boundaries at which a particle is typically destroyed and a thermal particle is injected, respectively. We also mention the resistivity used in the particles equations of motion η_p (either 10⁻⁹ or equal to the resistivity applied for the MHD evolution η_MHD = 10⁻⁴) and the magnetic Reynolds number R_m typical for the simulation parameters. In the last two columns, we show the particles maximum kinetic energy |$\mathcal {E}_{{\rm kin,max}}/(m_0 c^2) = \gamma _{{\rm max}} -1$| and the particles maximum energy |$\mathcal {E}_{{\rm max}} = \gamma m_0 c^2$| in MeV for each run.

3 RESULTS IN 2.5D CONFIGURATIONS

Flares are strongly transient phenomena and particles accelerate during such an event. Recently, Ripperda et al. (2017) combined MHD and test particle methods to investigate proton and electron acceleration in static MHD snapshots of repelling flux tubes in 2.5D. Here particle dynamics are evolved simultaneously with MHD evolution. The 2.5D results presented by Ripperda et al. (2017) show very hard energy distributions and an inverted power-law spectrum due to indefinite acceleration in the infinitely long current channels. Here we suggest several solutions to obtain more realistic distributions with a power-law spectrum and energies in accordance with observations, both in 2.5D setups and in 3D setups. The perturbed equilibrium of adjacent and antiparallel currents develops a tilt instability in which the current channels start to displace and rotate. This instability is indicated by an exponential growth phase of the kinetic energy, that is reached after t ≈ 4t_S in our 2.5D simulation (see Fig. 1 ). After this phase, the non-linear regime is reached at t ≈ 6t_S, showing highly chaotic behaviour and magnetic energy is converted into kinetic energy. In 3D both phases are delayed until t ≈ 6t_S and t ≈ 8t_S, respectively, due to the magnetic tension caused by the kinking of the channels (see Fig. 1). We evolve particles during the whole evolution shown in Fig. 1; however, we are mainly interested in particles accelerating due to the tilt instability from t ≈ 6t_S onwards. Particle acceleration is quantified by means of energy distributions and pitch angle distributions. Electron energy distributions associated with solar flares typically have a high-energy tail that partially can be fitted with a power-law function |$f(\mathcal {E}) \propto (\mathcal {E})^{-p}$| with p ≥ 1. A longer time spent in the current channels corresponds to higher energies and harder energy distributions. In Ripperda et al. (2017), it is shown that interacting flux ropes in 2.5D configurations are an efficient mechanism to accelerate particles. However, the spectra found are very hard and the power-law slope is even inverted (p < 0) compared to what is expected based on observations, even on very short time-scales Δt ≪ 0.1t_S. The maximum energies found also exceed electron energies of 5–50 MeV associated with solar flares. The main cause mentioned is the 2.5D character of the setup, meaning that the current channels have an infinite length and hence, particles can accelerate indefinitely. Here the length of the flux ropes is limited to 6L by applying a thermal bath at z = ±3L, in accordance with the periodic boundaries in 3D simulations. In Fig. 2, we show the kinetic energy [|$\mathcal {E}_{{\rm kin}}/(m_0 c^2) = \gamma - 1$|] distribution (left-hand panel) and the pitch angle [|$\alpha = \arctan (v_{\Vert }/v_{\perp })$|] distribution (right-hand panel) counted by particle number, for electrons in case A2De. The spectra are coloured by the MHD time t_S, from magenta to red, with magenta corresponding to early times and red to late times. The initial distributions are depicted by a dashed black line. Initially the electrons are distributed in the regions of the (displaced) current channels −1L ≤ x ≤ 1L; −2L ≤ y ≤ 2L from a Maxwellian with thermal speed |$v_{{\rm th,e}} = v_{{\rm th,p}} = \sqrt{2 k_{\rm B} T \rho _0/(m_{{\rm p}} p_0)}$|⁠, to improve statistics of particles in reconnection regions. We have confirmed that a simulation with an initially uniform electron distribution gives similar results, with the same maximum energy γ_max but a smaller fraction of particles in the high-energy tail. The limited flux rope length binds the kinetic energy to γ − 1 ≲ 10², corresponding to |$\mathcal {E}_{{\rm max}} \approx 50$| MeV, which is three orders of magnitude smaller than in the case with infinitely long flux ropes (Ripperda et al. 2017). Acceleration in the current channels at early times causes the slope of the spectrum to be inverted, p < 0. Acceleration in the direction parallel to the magnetic field is dominant over particle drifts as can be seen from the pitch angle distribution in the right-hand panel of Fig. 2 that is strongly peaked around α = 0. For protons (case B2Dp in Table 1) the kinetic energy spectra look fairly similar, with two major differences. The maximum kinetic energy is limited by γ − 1 ≲ 10⁻², corresponding to E_max ≈ 957 MeV, due to the mass difference between electrons and protons, and less protons have left the domain through the open x- and y-boundaries. The power-law index of the high-energy tail is inverted, p < 0. In the pitch angle spectra an asymmetry with respect to α = 0 can be observed, due to electrons and protons accelerating in opposite directions. This observation is visible in all simulations carried out in this work. The tendency for particle to accelerate along the magnetic field lines, and hence obtain a very small pitch angle, is in agreement with the findings of Gordovskyy et al. (2014). However, curvature acceleration resulting in increasing parallel velocity is neglected in their setup and particle collisions are incorporated. This asymmetry develops directly after t = 0, when protons accelerate parallel to the magnetic field. Therefore, there are more particles with a pitch angle slightly larger than zero. For electrons this asymmetry is present as well, with more particles with a pitch angle slightly smaller than zero. However, because electrons develop a larger parallel velocity than protons, the α = 0 peak is sharper for electrons and the asymmetry around α = 0 is less pronounced than for protons. For protons the α = 0 peak is represented by a peak at cos (ζ) = 1 and for electrons at cos (ζ) = −1 if we define cos (ζ) = v_‖/v with ζ the angle between the velocity vector and the magnetic field vector.

Figure 1.

Kinetic energy density evolution of the MHD fluid for 2.5D with effective resolution 2400² and 3D with effective resolution 300³, at all times integrated over a single current channel as identified by an advected tracer.

Open in new tabDownload slide

Figure 2.

2.5D kinetic energy distribution counted by particle number (left-hand panel) and pitch angle distribution (right-hand panel) for case A2De with 20 000 electrons from t = 0 to t = 10, solved for the guiding centre approximation with η_p = 10⁻⁴, with thermal bath applied and thermal velocity |$v_{{\rm th,p}}=\sqrt{2 k_{\rm B} T \rho _0/(m_{{\rm p}} p_0)}$|⁠. Time is measured in units of L/c_S, see the colour bar at the right, with L = 10 × 10⁶ m and c_S the speed of sound. The initial Maxwellian is depicted with a dashed, black line, showing the thermal part of the distribution.

Open in new tabDownload slide

4 RESULTS IN 3D CONFIGURATIONS

In 3D the perturbation in the velocity field consists of a z-component and dependence. This introduces variations in the z-direction that is invariant in 2.5D setups. The repelling and rotating current channels develop an additional kink instability that causes reconnection. Strong and thin current sheets develop at the boundaries and in between the two repelling islands (see Fig. 3 for the total current density magnitude in 2.5D and in 3D at t = 9t,). In Ripperda et al. (2017) it is shown that reconnection in this setup is indicated by a non-zero resistive electric field parallel to the magnetic field. This resistive electric field is plotted in Fig. 4 , with selected reconnecting magnetic field lines, initially at t = 0t_S and far into the non-linear regime at t = 9t_S. Particles mainly accelerate parallel to the resistive electric field and hence parallel to the magnetic field (Ripperda et al. 2017). In this section we investigate the effect of the kink instability on particle acceleration as well as the influence of the initial particle velocity distribution, the length of the flux ropes and resistivity for both electrons and protons. The effect of the initial spatial distribution and total number of particles on particle statistics are reported in Appendix A1 and A2. The effect of gyration is monitored in specific cases in Appendix A3.

Figure 3.

The total current density magnitude |$|\boldsymbol {J}|$| with a linear colour scale saturated to show values between [0, 50] in 2.5D (left-hand panel) and between [0, 10] in 3D (middle for a top view and right-hand panel for a side view). The secondary islands are visible in the 2.5D case, whereas thin current sheets are visible both in 2.5D and in 3D, indicated by a strong current magnitude.

Open in new tabDownload slide

Figure 4.

3D view of the parallel, resistive electric field E_|| with slices cut through the three axes, initially at t = 0t_S and in the non-linear regime at t = 9t_S, respectively. The colour scale is saturated at 0.001. Selected magnetic field lines in the current channel area are shown, coloured by their total current density value, saturated at j = 10. The box size is 6L × 6L × 6L with L = 10⁹ cm.

Open in new tabDownload slide

4.1 Effect of the kink instability on particle distributions

The effect of the kink instability is best visible for electrons. We show the kinetic energy spectra and the pitch angle spectra obtained from guiding centre simulations, with thermal velocity |$v_{{\rm th,e}} = v_{{\rm th,p}} = \sqrt{2 k_{\rm B} T \rho _0/(m_{{\rm p}} p_0)}$| with initially a fraction of 0.99 of the particles in the current channel area and flux ropes with length 6L in Fig. 5 for 20 000 electrons (case A3De). The main observable difference due to the kink instability is the development of a medium energy tail in the electron energy distribution with 10⁻⁵ ≤ γ − 1 ≤ 1 starting at t ≈ 8t_S in the non-linear regime. A second noticeable effect is the redistribution of electrons in the thermal distribution γ − 1 ≤ 10⁻⁵ after t ≈ 8t_S. This is attributed to the curvature of the magnetic field due to the kink, expelling particles from the current channels. These particles lose their energy in the ambient medium and either remain there and eventually leave the domain through the open boundaries, or they are caught again in either of the two current channels and re-accelerate. The electron energy distributions develop a high-energy peak at γ − 1 ≈ 60 (∼31 MeV) and the high-energy tail has an inverted power-law index p < 0. The difference in maximum energy compared to 2.5D results can be explained by the peak current reached in the MHD evolution. The current sheet is narrower in 2.5D due to a higher resolution (Ripperda et al. 2017), resulting in a larger peak current that accelerates particles to higher energy.

Figure 5.

3D distributions for case A3De with 20 000 electrons with η_p = 10⁻⁴ and 6L for the length of the current channels. Compared to 2.5D case A2De with equivalent settings in Fig. 2 a medium energy tail, visible as a bump in between the two peaks, after t ≈ 8t_S. The particles are also redistributed in the thermal distribution at late times.

Open in new tabDownload slide

In Fig. 6 we show kinetic energy (left-hand panel) and pitch angle distributions (right-hand panel) for case B3Dp, with 20 000 protons with thermal speed |$v_{{\rm th,p}} = \sqrt{2 k_{\rm B} T \rho _0/(m_{{\rm p}} p_0)}$| in the same setup as electrons in case A3De. The proton distributions are very similar to the 2.5D results. The maximum kinetic energy in the high-energy peak for protons is γ − 1 ≲ 2 × 10⁻² or |$\mathcal {E} \lesssim 957$| MeV and the high-energy tail has an inverted power-law index p < 0, similar to 2.5D case B2Dp. The pitch angle distributions are strongly peaked around α = 0. However, the asymmetry is less pronounced in 3D. This is attributed to randomization of pitch angles due to 3D effects in the non-linear regime after t ≈ 8t_S. The particles form two distinct populations. Protons and electrons moving inside the current channels develop a very high parallel velocity and a pitch angle very close to zero (positive for protons and negative for electrons). This is visible in the pitch angle distributions through the asymmetry around α = 0 at early times (before the linear growth phase of the tilt–kink instability at t ≈ 5t_S). This population of particles is also observed in the kinetic energy distributions at γ_max − 1 ≈ 2 × 10⁻² for protons and γ_max − 1 ≈ 6 × 10¹ for electrons. During the non-linear phase of the tilt–kink evolution (t ≳ 5t_S), particles are expelled from the channels due to the kink and reach less high energies (γ_max − 1 ≈ 2 × 10⁻³ for protons and γ_max − 1 ≈ 2 × 10¹ for electrons). In the pitch angle distributions, this is observed through the more symmetric and less high peak around α = 0. This is in accordance with Gordovskyy & Browning (2011b) concluding that the width of the peak in the pitch angle distribution depends on the (maximum) particle energy. For higher maximum energy, the particle species show a narrower peak. This is explained by taking into account that particles accelerate mostly along the magnetic field and drift velocities are negligible before reconnection occurs. Proton pitch angle distributions are wider than electron pitch angle distributions because electrons reach a larger parallel velocity and therefore the asymmetry at α = 0 is more pronounced for protons.

Figure 6.

3D distributions for case B3Dp with 20 000 protons with η_p = 10⁻⁴ and 6L for the length of the current channels. Compared to equivalent case A3De for electrons in Fig. 5, protons develop a less high maximum Lorentz factor and no medium-energy tail forms. The asymmetry due to the acceleration of electrons and protons in opposite directions is more visible for the heavier protons in the pitch angle distribution on the right-hand panel.

Open in new tabDownload slide

From the peaked pitch angle distributions in the right-hand panels of Figs 5 and 6, we conclude that parallel acceleration is dominant. However, which effect causes this is not clear by just comparing the particle drifts from equation (5). To analyse which acceleration mechanism causes this peak, the electrons are split in high-energy particles, with γ − 1 ≥ 1 (⁠|$\mathcal {E} \gtrsim 0.5$| MeV), and low-energy particles, with γ − 1 < 1. There are four contributions to the parallel acceleration in the momentum equation (6): the first two |$m_0\gamma \boldsymbol {u_E}\cdot (v_{\Vert }(\boldsymbol {\hat{b}}\cdot \nabla )\boldsymbol {\hat{b}})$| and |$m_0\gamma \boldsymbol {u_E}\cdot ((\boldsymbol {u_E}\cdot \nabla )\boldsymbol {\hat{b}})$| due to the change of direction of the magnetic field (curvature and polarization effects respectively). The third qE_‖ due to resistive electric field and the last one |$-\mu _{\rm r} \boldsymbol {\hat{b}}\cdot \nabla [B(1-E^{2}_{\perp }/B^2)^{1/2}]/\gamma$| is the mirror deceleration effect. In Fig. 7 the spatial distribution of high-energy electrons (left-hand panel) and the low-energy electrons (right-hand panel) is shown, coloured by (γ − 1) representing the particles kinetic energy (top panels) and by the curvature term |$\gamma \boldsymbol {u_E}\cdot (v_{\Vert }(\boldsymbol {\hat{b}}\cdot \nabla )\boldsymbol {\hat{b}})$| in the guiding centre momentum equation (6) (bottom panel). This curvature term is found to be dominant at all times due to the initially curved magnetic field and the kink instability adding further curvature (see also further on in Section 4.3 in Fig. 13). In the left-hand channel (seen from the top), the particles move upwards, because of the direction of the magnetic field and hence the current density; in the right-hand channel, they move downwards. The particles are assigned a thermal speed every time they cross a z-boundary and therefore all fast particles are on the left of Fig. 7 and all the slow particles are on the right, with the thermal particles at the foot points. The magnetic curvature is equally divided between left and right channels. However the fast particles in the left-hand panel of Fig. 7 travelled for a longer time in the current channels and obtained a larger acceleration and energy, mainly from the curvature. The curvature is stronger on the outside of the channels than on the inside and the fastest particles (indicated in red in the left-hand top panel) are therefore located at the foot points of the channels, on the outside, where also the current density is largest (Ripperda et al. 2017).

Figure 7.

Spatial distribution of electrons at t = 9t_S in case A3De, with thermal bath applied, particle resistivity η_p = η_MHD = 10⁻⁴. In the panels on the left-hand-side high-energy electrons, with γ − 1 ≥ 1, are shown. In the top panel, electrons are coloured by their kinetic energy E_kin/(m₀c²) = γ − 1. In the bottom left-hand panel the same high-energy electrons are coloured by the dominant contribution of the acceleration mechanism, the first term in equation (6), the magnetic curvature magnitude |$|\gamma \mathbf {u_E}\cdot (v_{\Vert }(\mathbf {\hat{b}}\cdot \nabla )\mathbf {\hat{b}})|$|⁠. In the right-hand panels the low-energy electrons, with γ − 1 < 1 are shown, coloured by their kinetic energy in the top panel and by the magnitude of the magnetic curvature acceleration in the bottom panel. The particles are projected in the y, z-plane with the line of sight along the x-direction. The particles follow magnetic field lines shown in 4 and at late times they reside mostly in the current channels, thus indicating the structure of the two current filaments.

Open in new tabDownload slide

4.2 Individual particle dynamics

To analyse how individual particles energize in the magnetic field of two interacting flux ropes, we look at electron trajectories in run A3De. It is interesting to look at the energy evolution of a particle that is initially travelling in the current channel but is then expelled from the current channel. In Fig. 8, we show the trajectory of such a particle from a side-view, coloured by its Lorentz factor. The electron cycles several times through the current channels until t = 9.177t_S. From then onwards the trajectory is portrayed by thick squares coloured by the time t counted in t_S, until t = 9.192t_S. This time interval corresponds to the red rectangle in the left-hand and the middle panel in Fig. 9, where the evolution of the Lorentz factor is shown. On the right-hand axis of the middle panel of Fig. 9, the z-position of the particle is depicted with a magenta dashed line to indicate where the particle is on the trajectory in Fig. 8. The particle is first injected with a thermal speed at the bottom boundary at z = −3. The particle accelerates and travels upwards until it is mirrored by the magnetic field at t ≈ 9.179t_S. Then it starts travelling downwards towards z = −3 and decelerates. At t ≈ 9.181 the particle reaches z = −3, it is injected at z = 3 with a thermal speed. It is then mirrored again towards z = 3. At t ≈ 9.181 it leaves at z = 3 and is injected with a thermal speed at z = −3. It accelerates until it reaches a Lorentz factor of γ ≈ 2 along the same trajectory upwards through the current channel. At t = 9.1835 it is expelled from the kinking current channel. Outside the current channel the particle moves downwards, following the magnetic field, until it reaches z = −3 at t ≈ 9.187t_S. It is injected again with a thermal speed at the top boundary at z = 3 and accelerates by following a field line outside the current channel. The particle Lorentz factor increases with time spent in the flux tube and the z-position increases (or decreases depending on the orientation) linearly in the current channels. In the right-hand panel of Fig. 9, the evolution of the magnitude of the drift terms from the right-hand-side of equation (5) is shown. The parallel velocity is dominant and quickly approaches the speed of light (black dashed line) when the particle accelerates.

Figure 8.

Side view in the y, z-plane of the trajectory of an electron that is expelled from the kinking current channel in case A3De in the time interval [8.80t_S, 9.25t_S]. The trajectory is coloured by the particles Lorentz factor. Between 9.177t_S and 9.192t_S, when the particle is expelled from the flux rope the trajectory is marked with thicker squares and coloured by the time t counted in t_S indicated by the second legend. Thermal bath boundary conditions are applied at the top and bottom boundaries, where the fast particle is destroyed and a thermal particle is injected at the opposite boundary from a Maxwellian with thermal speed |$v_{th,e} = \sqrt{2 k_{\rm B} T \rho _0/(m_{{\rm p}} p_0)}$| and γ_th,e = 1. It is clearly visible that the particle leaves the current channel at t ≈ 9.183t_S.

Open in new tabDownload slide

Figure 9.

Temporal evolution of the Lorentz factor of the particle in Fig. 8 zoomed into the time interval t = 9.177t_S to t = 9.192t_S, (the period coloured by time t in Fig. 8). To make a link with the complicated particles trajectory in Fig. 8 the z-position of the particle is plotted with a dashed magenta line on the right-hand axis of the left-hand panel. At t ≈ 9.179t_S the particle is mirrored and starts decelerating downwards in the channel. The particle is expelled from the current channel at t ≈ 9.183t_S. At t ≈ 9.181t_S and t ≈ 9.187t_S, the particle reaches a boundary and is injected at the opposite boundary with a thermal velocity by the thermal bath. In the right-hand panel the temporal evolution of the magnitude of all drift terms on the right-hand-side of equation (5) are shown, normalized to the speed of light (indicated by the black dashed line), during the same period.

Open in new tabDownload slide

4.3 Effect of resistivity on particle distributions

A resistivity of η_p = 10⁻⁴ results in a magnetic Reynolds number of R_m = 2 · 10⁵, which is three order of magnitude lower than in the solar corona. To moderate the acceleration parallel to the magnetic field, and therewith the peaked pitch angle spectra and the hard, inverted energy spectra, a resistivity model is proposed that lowers the resistivity to realistic solar corona values. To restrict the contribution of a resistive electric field, the resistivity in the GCA equations is set to η_p = 10⁻⁹ resulting in a Magnetic Reynolds number R_m = 2 · 10¹⁰ in the solar corona regime. The energy spectra for case C3De with η_p = 10⁻⁹ and periodic (infinite) flux ropes are shown in the left-hand panel of Fig. 10 and the pitch angle spectra in the right-hand panel. A high-energy peak develops at t ≈ 5t_S with maximum Lorentz factor γ − 1 ≲ 10³, resulting in |$\mathcal {E}_{{\rm max}} \approx 500$| MeV for electrons and a medium energy part with an inverted slope. The pitch angle is still strongly peaked at α = 0 and the energy spectra still show an inverted slope.

Figure 10.

Distributions for case C3De with 20 000 electrons in 3D MHD with lowered resistivity η_p = 10⁻⁹ and infinite length of the current channels (i.e. no thermal bath boundary condition applied). Compare to case A3De with higher resistivity and thermal bath in Fig. 5.

Open in new tabDownload slide

In case D3Dp 20 000 protons are evolved in similar settings as in case C3De for electrons. We find a final kinetic energy distribution (see left-hand panel of Fig. 11) with maximum Lorentz factor γ − 1 ≤ 2 · 10⁻² (⁠|$\mathcal {E} \lesssim 957$| MeV) and an inverted power-law slope. The pitch angle spectrum (see right-hand panel of Fig. 11) is peaked around α = 0 but the peak is less dominant than in case B3Dp, where η_p = η_MHD = 10⁻⁴ but the length of the flux ropes is limited to 6L.

Figure 11.

Distributions for case D3Dp with 20 000 protons uniformly distributed in 3D MHD with lowered resistivity η_p = 10⁻⁹ and infinite length of the current channels (i.e. no thermal bath applied). Compare to proton case B3Dp with higher resistivity and thermal bath in Fig. 6.

Open in new tabDownload slide

The high-energy particles with γ − 1 ≥ 1 (⁠|$\mathcal {E} \gtrsim 0.5$| MeV) in case C3De are depicted in Fig. 12, coloured by the kinetic energy (γ − 1) in the left-hand panel and by the magnitude of the curvature acceleration |$m_0\gamma \mathbf {u_E}\cdot (v_{\Vert }(\mathbf {\hat{b}}\cdot \nabla )\mathbf {\hat{b}})$| in the right-hand panel. Unlike for case A3De (compare to Fig. 7), the fast particles are not just located at the foot points, but distributed over the whole area of the current channels, with the fastest particles in the middle (see the left-hand top panel). This area corresponds to the region with the strongest current density (and hence resistive electric field) and magnetic field curvature [see Fig. 4 for the resistive electric field with magnetic field lines in this setup and Ripperda et al. (2017) for more detail on MHD results]. The slow particles, with γ − 1 < 1, are residing outside the current channels. Particles mostly accelerate in the regions with strongest curvature, at the outside of the kinked channels. This can also be seen from the distribution of the acceleration terms in the bottom panel of Fig. 13. For case C3De the distribution of curvature acceleration at t = 9t_S shows a tail consisting of the fast particles (in magenta with crosses as indicators). Compared to case A3De (Fig. 7) there is a clearer separation between slow and fast particles and the contribution of curvature acceleration to the fast particles. In case A3De (Fig. 7), all particles are accelerated by the magnetic curvature and the thermal bath prevents a tail to arise in the acceleration distribution.

Figure 12.

Spatial distribution of electrons at t = 9t_S in case C3De, without thermal bath, particle resistivity η_p = 10⁻⁵ × η_MHD = 10⁻⁹. In the left-hand panel high-energy electrons, with γ − 1 ≥ 1, are coloured by their kinetic energy E_kin/(m₀c²) = γ − 1. In the right-hand panel the same high-energy electrons are coloured by the dominant contribution of the acceleration mechanism, the first term in equation (6), the magnetic curvature magnitude |$|\gamma \mathbf {u_E}\cdot (v_{\Vert }(\mathbf {\hat{b}}\cdot \nabla )\mathbf {\hat{b}})|$|⁠. The particles are projected in the y, z-plane with the line of sight along the x-direction. The particles follow magnetic field lines shown in 4 and at late times they reside mostly in the current channels, thus indicating the structure of the two current filaments.

Open in new tabDownload slide

Figure 13.

Comparison of the number distribution of the four acceleration terms in the right-hand-side of equation (6) for 20 000 electrons in case C3De (dotted lines with crosses as indicators) and case A3De (solid lines with dots as indicators) before non-linear phase at t = 6t_S (top panel) and during non-linear phase at t = 9t_S (bottom panel). In case C3De periodic boundary conditions are applied (i.e. infinitely long flux ropes) and η_p = 10⁻⁵ × η_MHD = 10⁻⁹ and in case A3De the flux rope length is 6L and resistivity is set η_p = η_MHD = 10⁻⁴. The magnetic mirror term [fourth term in equation (6)] is indicated by a green line, the resistive electric field acceleration [third term in equation (6)] by a black line and the two terms attributed to the change of direction of the magnetic field by red [polarization term, the second term in equation 6] and magenta [magnetic curvature term, the first term in equation (6)].

Open in new tabDownload slide

The limited length of the flux ropes and the realistic particle resistivity separately counteract the strong parallel acceleration and the hard spectra found by Ripperda et al. (2017). These two solutions are compared by quantifying the contributions of the four separate terms in the momentum equation (6). We can determine whether the peaked distributions have a physical cause or are due to a too high MHD resistivity set for computational purpose. In Fig. 13 the distribution of the four terms contributing to the particle momentum in the momentum equation (6) are shown for case C3De and A3De at t = 6, before the current channels start kinking and at t = 9, in the non-linear regime, respectively. We see that at t = 6 the acceleration mechanisms show a similar distribution, compared between runs with thermal bath and η_p = 10⁻⁴ (case A3De) and the runs with periodic boundary conditions and η_p = 10⁻⁹ (case C3De), except the term due to the resistive, parallel electric field. For this term, |$q E_{\Vert }/m_0 = q \eta _{\rm p} \boldsymbol {J} \cdot \boldsymbol {\hat{b}}/m_0$| we can see the direct effect of decreased resistivity. At t = 9 however, in the bottom panel of Fig. 13, the distributions of the curvature acceleration contribution |$\gamma \boldsymbol {u_E}\cdot (v_{\Vert }(\boldsymbol {\hat{b}}\cdot \nabla )\boldsymbol {\hat{b}})$| and the acceleration due to polarization |$\gamma \boldsymbol {u_E}\cdot ((\boldsymbol {u_E}\cdot \nabla )\boldsymbol {\hat{b}})$| are shifted by two orders of magnitude for case C3De with η_p = 10⁻⁹. For the case with thermal bath (A3De), they are not. The mirror effect |$-\mu _{\rm r} \boldsymbol {\hat{b}}\cdot \nabla [B(1-E^{2}_{\perp }/B^2)^{1/2}]/\gamma$| is negligible in both cases C3De and A3De. The resistive acceleration qE_‖/m₀ has not changed with respect to the spectrum before the channels started kinking at t = 6t_S. Conclusively, the resistive electric field is not the main cause for the hard spectra obtained. The periodic boundary conditions and therewith the indefinite acceleration in the z-direction, in the current channels seem to be dominant and this effect is counteracted by thermal bath boundary conditions. However, even with a thermal bath applied, an inverted slope is observed in the energy spectra.

4.4 Realistic conditions for solar flares

In case E3De we combine a limited flux rope length and resistivity for realistic solar corona conditions. In a flux rope of length 6L and a resistivity of η_p = 10⁻⁹ we evolve 20 000 electrons, with the aim to restrict a high-energy peak and the accompanying inverted power-law index for the high-energy tail. The pitch angle distributions (right-hand panel of Fig. 14) are still dominated by α = 0, but there are more particles with a non-zero pitch angle, compared to cases A3De and C3De. The high-energy peak observed in cases A3De and C3De has now disappeared in the kinetic energy distributions (left-hand panel of Fig. 14). The slope of the high-energy tail at times t ≳ 7t_S (coloured green to red), when the channels start kinking, is still inverted. However, at this time there are not many particles left in the thermal part of the distribution. The maximum kinetic energy is bounded by γ − 1 ≲ 5 (⁠|$\mathcal {E} \lesssim 3$| MeV), due to the thermal bath and the length of the channels. How many particles remain thermal is affected by the initial velocity distribution. The particles maximum energy is bounded by the time the particles spend in the flux rope, and hence by the length of the flux rope.

Figure 14.

Distributions for case E3De with 20 000 electrons in 3D MHD with lowered resistivity η_p = 10⁻⁹ and a limited length of 6L for the length of the current channels (i.e. thermal bath applied). Compare to cases A3De with higher resistivity (Fig. 5) and C3De with infinitely long current channels and lowered resistivity η_p = 10⁻⁹ (Fig. 10).

Open in new tabDownload slide

4.5 Effect of initial velocity distribution on particle distributions

Up to now it is assumed that all particles, electrons and protons, have the typical MHD velocity |$v_{{\rm th,p}}=\sqrt{2 k_{\rm B} T \rho _0/(m_{{\rm p}} p_0)}$| as thermal speed. Assuming energy equipartition between electrons and protons results in a larger thermal speed for electrons |$v_{{\rm th,e}}= \sqrt{m_{\rm p}/m_{\rm e}} \times v_{{\rm th,p}}$|⁠. In case G3De we explore the effect of a larger initial electron velocity with flux rope length 6L and particle resistivity set as η_p = 10⁻⁹. Electrons are initialized from a Maxwellian with thermal speed |$v_{{\rm th,e}} = \sqrt{m_{\rm p}/m_{\rm e}} \times v_{{\rm th,p}}$| in accordance with energy equipartition and they are uniformly distributed over the spatial domain −3L ≤ x ≤ 3L; −3L ≤ y ≤ 3L; −3L ≤ z ≤ 3L. Few particles accelerate inside the current channels at early times t < 6t_S in the linear phase (see Fig. 15). The maximum electron energy is limited by the length of the current channels and the thermal bath to γ − 1 ≲ 7 or |$\mathcal {E} \lesssim 4.1$| MeV, which is in the range of observed electron energies coming from solar flares. New thermal particles are injected from the thermal bath for every destroyed high-energy particle leaving a periodic boundary, maintaining a thermal distribution dominant in number of particles. The high-energy tail that develops during the exponential growth phase of the tilt–kink instability after t ≥ 6t_S has a power-law distribution with spectral index p > 1 (indicated by the dotted black line in the left-hand panel of Fig. 15). The high-energy particles are dominant neither in number of particles nor in energy content due to decreased particle resistivity and the thermal bath. The pitch angle distributions in the right-hand panel of Fig. 15 are initially (nearly) uniform and the peak due to particles accelerated parallel to the magnetic field at α = 0 develops only after t ≥ 6t_S. The peak at α = 0 is less dominant than in all other electron cases.

Figure 15.

Distributions for 20 000 electrons in case G3De with an initially uniform spatial distribution with thermal velocity |$v_{{\rm th,e}}=\sqrt{2 k_{\rm B} T \rho _0/(m_{{\rm p}} p_0)}\sqrt{m_{{\rm p}}/m_{{\rm e}}}$| and lowered resistivity η_p = 10⁻⁹ and 6L for the length of the current channels (i.e. thermal bath applied). Compare to case E3De (Fig. 14) with an initially smaller thermal velocity v_th,e = v_th,p for electrons. As a guide for the eye a typical power-law distribution of the form dN/d(γ − 1) ∝ (γ − 1)⁻¹ is plotted (dotted black line), corresponding to equal energy content in each decade of γ − 1.

Open in new tabDownload slide

In case F3Dp we explore the same configuration for protons, initiated uniformly from a Maxwellian with v_th,p. For protons assuming energy equipartition or a generic fluid velocity as thermal speed results in the same initial energy distribution. The results are similar to electron case G3De with the difference that the maximum energy is limited to γ − 1 ≲ 3 · 10⁻² or |$\mathcal {E} \lesssim 941$| MeV (see the left-hand panel of Fig. 16) due to the mass difference and consequently the lower thermal speed. Few particles accelerate inside the current channels at early times t < 6t_S in the linear phase. In the non-linear phase from t > 6t_S onwards, a high-energy tail develops with a power-law distribution with index p > 1 (indicated by the dotted black line in the left-hand panel of Fig. 16). Because of the lower maximum energy reached, less protons leave the domain through the open x, y-boundaries, compared to electrons in case G3De and the thermal distribution remains dominant at all times. The pitch angle distribution (see the right-hand panel of Fig. 16) remains (nearly) uniform till t = 6t_S and afterward it is peaked around α = 0 but at least an order of magnitude smaller than in all other proton cases. In Fig. 17 contributions of the four separate terms in the momentum equation (6) are quantified. Comparing to Fig. 13 for cases A3De and C3De the effect of a finite flux rope length and low resistivity combined is shown for the four acceleration mechanisms. We see that at t = 9 the magnetic curvature is the dominant acceleration mechanism. There is no peak in the distribution. The magnetic curvature acceleration is three orders of magnitude lower than in case C3De (with η_p = 10⁻⁹ and periodic flux ropes) and one order of magnitude lower than in case A3De (with finite flux rope length and η_p = 10⁻⁴). The other three acceleration mechanisms are at least three orders of magnitude smaller than the curvature acceleration. The curvature acceleration term in the GCA momentum equation (6) is proportional to the parallel velocity of the particle. A particle following a (curved) field line accelerates parallel to the field line. Even in the initially straight current channel the field lines are curved due to the initial magnetic field distribution in equation (2). Conclusively, a realistic magnetic Reynolds number and finite flux rope length result in a realistic maximum energy reached, power-law index of the high-energy tail of the kinetic energy distribution and the shape of the pitch angle distributions. Case G3De for electrons and case F3Dp for protons give the most realistic results for particle acceleration due to interacting flux ropes in the solar corona. In the next section, we explore the effect of the finite length of a flux rope for the settings of case G3De.

Figure 16.

Distributions for 20 000 protons in case F3Dp with an initially uniform spatial distribution and lowered resistivity η_p = 10⁻⁹ and 6L for the length of the current channels. Compare to cases B3Dp (Fig. 6) with infinitely long current channels and lowered resistivity η_p = 10⁻⁹ and D3Dp (Fig. 10) with thermal bath applied and higher resistivity η_p = 10⁻⁴. As a guide for the eye, a typical power-law distribution of the form dN/d(γ − 1) ∝ (γ − 1)⁻¹ is plotted (dotted black line), corresponding to equal energy content in each decade of γ − 1.

Open in new tabDownload slide

Figure 17.

The number distribution of the four acceleration terms in the right-hand-side of equation (6) for 20 000 electrons in case G3De at t = 9t_S. Compare to Fig. 13 for cases A3De and C3De. The length of the flux rope is 6L and resistivity is set as η_p = 10⁻⁵ × η_MHD = 10⁻⁹. The magnetic mirror term [fourth term in equation (6)] is indicated by a green line, the resistive electric field acceleration [third term in equation (6)] by a black line and the two terms attributed to the change of direction of the magnetic field by red [polarization term in equation (6)] and magenta [magnetic curvature term in equation (6)].

Open in new tabDownload slide

4.6 Effect of the length of the flux ropes on particle distributions

In case of periodic boundary conditions, particles are either expelled from the flux rope or they travel through the flux rope until the simulation ends. When they are expelled from the flux rope they thermalize in the ambient medium until they leave the domain or they are caught into the flux ropes again. The particles causing the high-energy peak are in the current channels for a long time and typically cycle through the current channel many times. These fast electrons travelling in a flux rope typically accelerate up to γ_max ≈ 1 · 10³ if the flux rope is infinitely long (i.e. if no thermal bath is applied as in case C3De for electrons). During the simulation an electron with a parallel velocity v_‖ ≈ c travels a maximum distance of Δl ≈ c × 10t_S ≈ 2.6 · 10¹³ cm, or over 4000 cycles through the flux ropes, if it is not expelled at some point. This is equal to Δl ≈ 3.6 · 10² solar radii in 10 soundcrossing times t_S. The medium-energy tail in case C3De (see Fig. 10) consists of particles with an average Lorentz factor of γ − 1 ≈ 10⁻². This is equal to a parallel velocity of v_‖ ≈ 0.14c. Even these particles can travel a maximum distance of Δl ≈ 0.14c × 10t_S ≈ 3.4 · 10¹² cm, or 52 solar radii through the flux ropes. Thermal bath boundary conditions limit this to one cycle, or Δl = 6L ≈ 8.6 · 10⁻² solar radii. To analyse how the length of the flux rope affects the energy distribution, case H3De has settings similar to case G3De, with the total length of the current channels 12L (or 1.7 · 10⁻¹ solar radii) instead of 6L. All particles are initially in the current channel area. To obtain equal accuracy the resolution is also doubled in the z-direction. The shape of the distributions for case H3De in Fig. 18 is similar to case G3De (Fig. 15). The maximum energy reached by the particles in the flux ropes with length 12L is approximately doubled compared to case G3De (γ_max − 1 ≲ 2 · 10¹ or |$\mathcal {E} \lesssim 11$| MeV). Initially few particles accelerate in the flux tubes and the energy distribution in the left-hand panel remains largely Maxwellian. In the non-linear phase, after t = 6t_S a non-thermal power-law tail forms due to the tilt–kink instability and subsequent reconnection. The pitch angle distribution in the right-hand panel remains nearly flat until t = 6t_S, after which the peak around α = 0 becomes dominant due to particles accelerating along the magnetic field.

Figure 18.

Distributions for 20 000 electrons in case H3De with lowered resistivity η_p = 10⁻⁹, initial thermal speed |$v_{{\rm th,e}}=\sqrt{2 k_{\rm B} T \rho _0/(m_{{\rm p}} p_0)}\sqrt{m_{{\rm p}}/m_{{\rm e}}}$| and now double the length 12L of the current channels (i.e. thermal bath applied). Compare to case G3De (Fig. 15) with half the length (6L) of the current channels. As a guide for the eye, a typical power-law distribution of the form dN/d(γ − 1) ∝ (γ − 1)⁻¹ is plotted (dotted black line), corresponding to equal energy content in each decade of γ − 1.

Open in new tabDownload slide

5 CONCLUSIONS

We find that reconnection in the low plasma-β regime drives efficient energy conversion from magnetic energy to kinetic energy and that electrons and protons are efficiently accelerated to non-thermal energy distributions. We observe two populations of high-energy particles. A high-energy peak of particles trapped inside the current channels at early times, where they accelerate efficiently along the magnetic field. Electrons reach maximum energies between ∼10 and ∼500 MeV. The maximum energy depends strongly on initial velocity distribution, plasma resistivity and the length of the flux tubes. A second population consists of electrons accelerating in the reconnection zones at late times in the non-linear phase. These electrons generate a high-energy tail in between the peak and the Maxwellian part of the distribution due to particles accelerating, in quite a narrow range between ∼0.5 and ∼10 MeV. Protons reach maximum energy of approximately 1000 MeV in all cases.

We proposed two solutions to limit indefinite particle acceleration due to infinitely long flux ropes as found in Ripperda et al. (2017) for 2.5D simulations of low plasma-β environments and extended this setup to full 3D simulations for electrons and protons. One solution is to apply a thermal bath along the length-direction of the flux ropes at the periodic boundaries (case A3De for electrons and B3Dp for protons). This solution effectively limits the length of the flux ropes to realistic scales and destroys fast particles leaving the domain. For every particle ending up in the thermal bath, a new thermal particle is injected at the opposite periodic boundary. This assures to keep a steady total number of particles. The particles maximum kinetic energy grows linearly with the length of the flux ropes and the time spent in the flux rope. In a flux rope of length 6L the maximum particle energy is limited to γ − 1 ≲ 100, or |$\mathcal {E} \lesssim 50$| MeV for electrons and γ − 1 ≲ 2 · 10⁻², or |$\mathcal {E} \lesssim 960$| MeV for protons. The majority of particles reach a non-thermal energy before they reach the end of the flux rope. A high-energy tail forms in all cases with limited flux rope length. The shape of the high-energy tail and maximum energy depend on the initial conditions and resistivity set for the simulation.

The second solution counteracts acceleration in the direction of the resistive electric field, parallel to the magnetic field, as proposed by Zhou et al. (2016). A particle resistivity is applied that is typically lower than the MHD resistivity and results in realistic magnetic Reynolds number for solar corona conditions (case C3De for electrons and D3Dp for protons). This particle resistivity effectively lowers the electric field that is felt by the particles, but does not affect the MHD evolution and the formation of strong and thin current sheets. It limits the number of particles that accelerate to high energy. Despite having less particles in the high-energy peak, the maximum particle energy is not limited due to the infinite length of the (periodic) flux ropes. Electrons accelerate to γ − 1 ≲ 1000 or |$\mathcal {E} \lesssim 500$| MeV and protons to γ − 1 ≲ 2 · 10⁻² or |$\mathcal {E} \lesssim 960$| MeV. A high-energy distribution forms in all cases and its shape and maximum energy depends on the initial conditions and the resistivity set.

Despite limited maximum particle energy an inverted power-law index (p < 0) is found for the non-thermal distribution in the cases with low resistivity or limited flux rope length. This is due to the high-energy peak developing from particles accelerating in the current channels at early times. These high-energy particles dominate the energy census at late times. To limit the number of particles in the high-energy tail we combined limiting the flux rope length and the resistivity to realistic values (cases E3De, G3De and H3De for electrons and F3Dp for protons). This results in setups with realistic magnetic Reynolds number and a flux rope length limited to a fraction of a solar radius. The maximum energy grows with the length of the flux ropes and for electrons |$\mathcal {E} \lesssim 4$| MeV for length 60 Mm and |$\mathcal {E} \lesssim 11$| MeV for length 120 Mm and for protons |$\mathcal {E} \lesssim 1$| GeV. The shape of the slope of the high-energy tail formed depends on the initial velocity distribution. Assuming energy equipartition between electrons and protons, we find high-energy power-law distributions |$f(\mathcal {E}) \sim (\mathcal {E})^{-p}$| with p ≥ 1, as opposed to the inverted spectra found if the flux rope length is not limited and the resistivity is not lowered to realistic solar corona values. The high-energy particles are dominant neither in number nor in energy, such that the test particle approximation is valid. In all these cases a part of the particle ensemble and its energy content is in the non-thermal distribution. These findings are in good agreement with the 2D kinetic PIC results of Li et al. (2015) for particle acceleration in force-free current sheets in low-β electron–proton plasmas. The non-thermal particle distributions found can explain the efficient electron acceleration in low plasma-β environments such as solar flares.

Without applying the two solutions proposed the whole ensemble of particles contains so much energy that kinetic feedback of the particles to the electromagnetic fields cannot be neglected. The assumption that the energy content of the particles is much lower than that of the fluid is invalid. A kinetic description allows particles to lose their energy to the fields through kinetic instabilities and particle–field interaction that is ignored in the test particle approximation. However at solar length scales, considering that reconnection occurs globally in the simulation box, PIC simulations are extremely demanding numerically. In cases with realistic magnetic Reynolds numbers and solar length scales for the flux ropes the number of accelerated particles is small compared to the total number of particles. The fields induced by these particles should not substantially affect the MHD evolution and reconnection.

We applied a guiding centre approximation, ignoring the gyration of particles and reducing computing time. Via the guiding centre approximation we demonstrated that magnetic curvature is the leading acceleration mechanism in all cases. The magnetic field inside the current channels is curved initially and the kink instability introduces more curvature, making the curvature acceleration dominant at all times. The curvature acceleration is proportional to the velocity of a particle parallel to the magnetic field, enhancing the effect of the curvature. Particles are mainly accelerated parallel to the magnetic field and hence the resistive electric field, resulting in pitch angle distributions dominated by a peak at α = 0 in all cases. The width and height of the pitch angle peak depends on the maximum parallel particle velocity and hence on particle energy. The α = 0 peak shows an electron/proton asymmetry, previously observed by Gordovskyy et al. (2014), caused by protons predominantly moving with positive parallel velocity and electrons mostly with negative parallel velocity along the magnetic field.

The guiding centre approximation is valid since the gyration radius remains much smaller than typical cell size in all runs. For protons, the gyration might not be negligible due to a larger mass and hence a larger gyroradius. To monitor the validity of our approach, results are compared to a run where proton gyration is not neglected and the full equations of motion (4) are solved in Appendix A3 (case K3Dp). Applying the GCA has little to no effect on the energy distribution, compared to a case where particle gyration is fully resolved. The GCA is accurate in magnetically dominated, Newtonian plasmas under solar corona conditions. Results are also confirmed for both an ensemble of 200 000 particles (Appendix A1, case I3De) and for an initially uniform spatial distribution (Appendix A2, case J3De) to show that obtained statistics are accurate. Since kinetic feedback of the particles on the fields is neglected, the accumulated current from energetic particles is unimportant and it is not necessary to increase the total number of particles.

Acknowledgements

This research was supported by projects GOA/2015-014 (2014–2018 KU Leuven) and the Interuniversity Attraction Poles Programme by the Belgian Science Policy Office (IAP P7/08 CHARM). OP is supported by the ERC synergy grant ‘BlackHoleCam: Imaging the Event Horizon of Black Holes’ (Grant No. 610058). The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation Flanders (FWO) and the Flemish Government – department EWI. BR likes to thank Lorenzo Sironi, Fabio Bacchini, Norbert Magyar, Jannis Teunissen, Kirit Makwana, Matthieu Leroy and Dimitris Millas for fruitful discussions and comments.

REFERENCES

APPENDIX A: VALIDATION TESTS

A1 Number of particles

To assure that the results are statistically accurate, a setup similar to case A3De is evolved with an initial distribution of 200 000 electrons for case I3De. Minor differences are visible in Fig. A1, mainly due to smoother distributions. The energy range and global features of both the kinetic energy distributions and the pitch angle distribution are in agreement with results for 20 000 electrons. Evolving more particles has no effect on the physical results, other than improving statistics, since there is no feedback of the particles on the fields, nor(( any particle–particle interaction.

Figure A1.

Kinetic energy distribution as in Fig. 5 for case I3De with 200 000 electrons in 3D MHD with η_p = 10⁻⁴ and 6L for the length of the current channels.

Open in new tabDownload slide

A2 Spatial distribution

Based on the findings of Ripperda et al. (2017) for 2.5D configurations, a fraction of 0.99 of the particles is initially distributed in the area of the current channels −1L ≤ x ≤ 1L; −2L ≤ y ≤ 2L; −3L ≤ z ≤ 3L. In this way using 20 000 particles results in accurate statistics, similar to a run with 200.000 particles, and computational resources are mainly used on particles that are likely to accelerate. In 3D configurations, a larger area of the box is filled with reconnecting magnetic field due to the kink instability. To monitor whether it is accurate to distribute most particles in the region of the current channels in 3D configurations, a setup similar to case E3De with thermal bath and η_p = 10⁻⁹ is run for case J3De, now with an initially uniform spatial distribution. Electrons are randomly and uniformly initialized in the whole box −3L ≤ x ≤ 3L; −3L ≤ y ≤ 3L; −3L ≤ z ≤ 3L. The resulting energy distributions in Fig. A2 show the same shape and trend and the maximum energy γ − 1 ≲ 5 is equal to the maximum energy obtained in case E3De in Fig. 14. The main difference is that there are less particles in the high-energy tail. However, there are not more particles remaining in the thermal distribution. The particles in the ambient medium, outside the current channels, leave the domain guided by the field lines on the time-scales considered. Also the pitch angle is still peaked around α = 0, but the peak is an order of magnitude lower than in case E3De. Distributing particles in the area of the current channels rather than in the ambient medium improves statistics on high-energy particles, but does not alter the physical features of the obtained spectra. This is in accordance with the 2.5D results of Ripperda et al. (2017).

Figure A2.

Kinetic energy distribution as in Fig. 14 for case J3De with 20 000 electrons in 3D MHD, uniformly distributed in the domain with η_p = 10⁻⁹ and 6L for the length of the current channels.

Open in new tabDownload slide

A3 Gyration

Protons have a larger gyroradius than electrons due to their mass and therefore the GCA is most likely to break down in proton simulations. To monitor the validity of the GCA, the full particle equations of motion (4) including gyromotion are evolved with a Boris scheme for protons in the same MHD background. In Fig. A3 we show the kinetic energy distribution for 20 000 protons in a setup where particle gyration is resolved. Case D3Dp is chosen for comparison to confirm that the high-energy peak found is not a numerical artefact. The energy distribution shows a very similar shape. A peak develops due to particles accelerating in the current channels and the maximum energy reached is limited by the length of the current channels and the thermal bath boundary conditions, such that the maximum energy is similar to case D3Dp, with γ − 1 ≤ 4 · 10⁻² (⁠|$\mathcal {E} \lesssim 976$| MeV). The high-energy tail shows an inverted power-law index, developing at similar time as in case D3Dp. Besides minor differences in the thermal part of the spectrum, due to the random assignment of a thermal velocity component v_x, v_y and v_z rather than v_‖ and v_⊥ in the GCA setups, the global aspects are similar to the GCA results in Fig. 11. For electrons the differences are expected'; to be even smaller, due to a smaller mass and hence completely negligible gyroradius.

Figure A3.

Kinetic energy distribution as in Fig. 11 for case K3Dp with 20 000 protons in 3D MHD, solved with the full equations of motion, with η_p = 10⁻⁹ and infinite length of the current channels.

Open in new tabDownload slide