© The Authors 2025

1. Introduction

Solar flares cannot be explained without considering the effects of non-thermal electrons. The presence of these electrons is inferred directly through X-ray and radio emission, and indirectly through the optical, ultraviolet, and extreme ultraviolet radiation produced by the collisional heating they cause (Benz 2016). Magnetic reconnection is believed to be responsible for providing the energy that is required for the particles to reach non-thermal velocities, but the processes responsible for the particle acceleration are debated (Gordovskyy et al. 2019; Browning 2024). The observed signatures point towards a non-thermal power law energy distribution that varies significantly from flare to flare (Christe et al. 2008) and also as a function of space and time within a flare (Grigis & Benz 2004, 2006), indicating that the particle acceleration is sensitive to the large-scale magnetic configuration.

Magnetohydrodynamic (MHD) simulations are able to resolve large scales realistically (Hansteen et al. 2015), and they are able to approximate regions with magnetic reconnection using sophisticated resistivity models (see Færder et al. 2023 and references therein). Nevertheless, kinetic models are needed to describe the acceleration and transport of non-thermal particles, since MHD assumes thermal energy distributions.

Non-thermal particle beam models allow us to transport already accelerated particles in MHD simulations and deposit their energy along the way. In 3D, this was first done by Frogner et al. (2020) and was followed by more realistic beam propagation in Frogner & Gudiksen (2024). In Ruan et al. (2020) and Druett et al. (2024), the energy depositions are fed back into the MHD simulation, and the beams’ energy gains are subtracted from the MHD energy density. This enables self-consistent simulation of the energetics and atmospheric response of large flares where the accelerated particle beams contribute to chromospheric heating and upflows. However, particle beam models rely on the free parameters that define the initial beam velocity distribution. The study of particle acceleration helps constrain these parameters.

Non-thermal particle acceleration can occur on both large and small (kinetic) scales. Particle-in-cell (PIC) simulations by Drake et al. (2006) have shown that repeated first-order Fermi acceleration in contracting magnetic islands (or plasmoids) can be an efficient mechanism. In addition, PIC simulations of the magnetosphere by Egedal et al. (2012) indicate that electrons in reconnection exhausts are significantly accelerated by large-scale parallel electric fields. However, PIC simulations are too computationally expensive to include both the kinetic and large scales of flares on the Sun. A compromise is to surround PIC simulations with MHD boundaries, as done by Baumann et al. (2013). They found that parallel electric fields from the magnetic reconnection regions contributed the most to electron acceleration, and that the shape of the non-thermal power law depended mostly on the geometry and evolution of the magnetic field. Another approach is to simultaneously represent the electrons in the plasma as both a fluid and as particles. Doing this in 2D, while at the same time eliminating kinetic-scale parallel electric fields and approximating the particle motion with the guiding centre approximation (GCA), Arnold et al. (2021) simulated strong non-thermal acceleration via the first-order Fermi process, similar to Drake et al. (2006), but this time with a larger simulation domain, of the order of 10⁴ km. Finally, a third approach is to combine large-scale MHD simulations with test particles, which have been shown to accelerate electrons up to MeV energies (Gordovskyy et al. 2014) and can be used to synthesise observational signatures such as hard X-ray emission from MHD flare simulations (Pinto et al. 2016).

Test particles do not alter the external electromagnetic field in which they are embedded, nor do they interact with each other. This loss of self-consistency comes with the benefit of a reduction in computational cost, as the particle simulation reduces to a Monte Carlo problem and the external fields can be given by MHD or even analytically prescribed. The reduced cost makes it possible to simulate a very large number of particles, tracing them over very large distances and up to very high energies. From particle trajectories, one can study the characteristics of acceleration and infer statistical properties of the plasma (Marchand 2010). Stochastic approaches that simulate the interaction of test and background particles have frequently been used, starting with MacKinnon & Craig (1991) and explored by Fletcher (1995), Burge et al. (2014), Pinto et al. (2016), and Borissov et al. (2020). Without background collisions, the test particle simulations rely on two strong assumptions: collisionless motion and negligible particle feedback.

Collisionless motion is a valid approximation as long as the paths of the particles are shorter than the mean free path of the plasma. The non-thermal particles are believed to be accelerated in the thin and hot solar corona (Krucker et al. 2008), making collisionless motion a common assumption in test particle studies (Rosdahl & Galsgaard 2010; Gordovskyy et al. 2020; Zhao et al. 2021; Pallister et al. 2021; Bacchini et al. 2024; Mowbray et al. 2025). Excluding particle feedback has been suggested to artificially increase particle acceleration as it ignores the effect of return currents (Rosdahl & Galsgaard 2010). The exclusion is valid only as long as the electric current produced by the accelerated particles is negligible compared to the background currents. This depends on the fraction of non-thermal particles relative to the fraction of thermal particles. In some large flares, X-ray, microwave, and extreme ultraviolet measurements indicate that a large fraction, possibly all, of the electrons in a coronal volume can be accelerated (Krucker et al. 2010; Krucker & Battaglia 2014), but statistical studies show that in smaller flares the total energy in non-thermal particles is a smaller fraction of the peak thermal energy (Warmuth & Mann 2020) than in larger flares. From this we can infer that smaller flares accelerate a smaller portion of the particle population. This may restrict the use of test particles to smaller flares.

A source of smaller flares is magnetic reconnection in coronal, fan-spine topologies. Here, the current sheet at the magnetic null point can be a source of particle acceleration (Pallister et al. 2021). If the current sheet is stretched thin with sufficient resolution and appropriate plasma resistivity, the plasmoid instability will occur (Furth et al. 1963). Plasmoids of various sizes grow and merge with the exhaust of reconnection jets. This null point region resembles the flare model in Nishizuka & Shibata (2013), where plasmoid ejection and collision from the reconnection region drive particle acceleration.

We used a static snapshot from an MHD simulation of plasmoid-mediated reconnection in a 2D fan-spine magnetic topology as a background for test particle simulations, and studied the characteristics of emerging non-thermal electrons. The MHD simulation is of high enough resolution and large enough extent to resolve both plasmoid formation at the reconnecting null point and the large-scale fan-spine magnetic topology. We used the GCA for the particles’ equations of motion and simulated particle ensembles that are statistically representative of the MHD plasma. To ensure the ensembles contain non-thermal electrons, we used importance weighing to bias the particle sampling with knowledge gained from an initial test particle simulation.

2. Methods

The test particle simulation consisted of four parts: the equations of motion of the particles; the electric and magnetic fields that the forces acting on the particles depend on; a method for choosing the initial conditions of the particles; and a method for solving the equations of motion.

2.1. The particles’ equations of motion

A non-relativistic test particle embedded in an external magnetic (B) and electric field (E) is subject to the Lorentz force,

$$\begin{array}{c}\hfill \mathit{F}=q(\mathit{E}+\mathit{v}\times \mathit{B}),\end{array}$$(1)

where q and v are the charge and velocity of the particle, respectively. In typical solar coronal conditions, electrons have a very small Larmor radius and fast gyration due to their low mass-to-charge ratio and due to strong magnetic fields. Resolving their helical motion quickly becomes a computationally demanding task. To overcome this problem, we used the GCA, which averages the particle’s motion over its gyration. The set of equations that describe the motion of the guiding centre is valid as long as the electromagnetic field is sufficiently smooth over a gyration period. We considered the approximation valid whenever

$$\begin{array}{c}\hfill \frac{r_L}{L_B}\ll 1,\end{array}$$(2)

where r_L is the Larmor radius of the particle and L_B = B/|∇B| is the characteristic length of the external magnetic field.

Northrop (1961) derived the equations describing the time evolution of the position of the guiding centre R (in three dimensions) and the time evolution of its velocity component parallel to the magnetic field v_∥. Essentially, the motion of the guiding centre is split into a component parallel to the magnetic field and a component perpendicular to the magnetic field. The perpendicular component is represented by a set of drifts. Here we used the non-relativistic version of the equations (see Ripperda et al. 2018 for an extensive discussion of the equations). They are

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\mathit{R}}{\mathrm{d}t}& =v_{\Vert }\mathit{b}+{\mathit{v}}_E+{\mathit{v}}_{\mathrm{\nabla }B}+{\mathit{v}}_C+{\mathit{v}}_P\hfill \end{array}$$(3a)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{v}_{\Vert }}{\mathrm{d}t}& =\frac{1}{m}(q{E}_{\Vert }-\mu \mathit{b}·\mathrm{\nabla }B)+{\mathit{v}}_E·\frac{\mathrm{d}\mathit{b}}{\mathrm{d}t},\hfill \end{array}$$(3b)

where b = B/B is the magnetic field direction, E_∥ = E ⋅ b is the electric field component parallel to b, μ is the magnetic moment of the particle (which is assumed to be constant in the GCA), and m is the particle mass. The perpendicular drifts are

$$\begin{array}{ccc}\hfill {\mathit{v}}_E& =\frac{\mathit{E}\times \mathit{b}}{B}\hfill & \hfill (E\mathrm{cross}B\mathrm{drift}),\end{array}$$(4a)

$$\begin{array}{ccc}\hfill {\mathit{v}}_{\mathrm{\nabla }B}& =\frac{\mu }{\mathit{qB}}(\mathit{b}\times \mathrm{\nabla }B)\hfill & \hfill (\mathrm{gradient}B\mathrm{drift}),\end{array}$$(4b)

$$\begin{array}{ccc}\hfill {\mathit{v}}_C& =\frac{m{v}_{\Vert }}{\mathit{qB}}(\mathit{b}\times \frac{\mathrm{d}\mathit{b}}{\mathrm{d}t})\hfill & \hfill (\mathrm{curvature}\mathrm{drift}),\end{array}$$(4c)

$$\begin{array}{ccc}\hfill {\mathit{v}}_P& =\frac{m}{\mathit{qB}}(\mathit{b}\times \frac{\mathrm{d}{\mathit{v}}_E}{\mathrm{d}t})\hfill & \hfill (\mathrm{polarisation}\mathrm{drift}).\end{array}$$(4d)

The first term in Eq. (4) is the acceleration by parallel electric fields. The second term is responsible for magnetic mirroring, and the third term is first-order Fermi acceleration. The latter can be associated with the energy change due to a curvature drift along the electric field, since

$$\begin{array}{c}\hfill {\mathit{v}}_E·\frac{\mathrm{d}\mathit{b}}{\mathrm{d}t}=\frac{{\mathit{v}}_C·q\mathit{E}}{m{v}_{\Vert }}·\end{array}$$(5)

We defined the kinetic energy of the particle as the sum of three components, the parallel velocity, the speed of the gyration (v_⊥), and the magnitude of the vector sum of the drifts (v_Drifts):

$$\begin{array}{c}\hfill K=\frac{m}{2}(v_{\Vert }^2+v_{\mathrm{Drifts}}^2)+B\mu .\end{array}$$(6)

When we evaluated the material derivatives in Eq. (3), we neglected drifts other than v_E, which is usually the dominant one. The material derivative becomes

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\equiv \frac{\partial }{\partial t}+v_{\Vert }\mathit{b}·\mathrm{\nabla }+{\mathit{v}}_E·\mathrm{\nabla },\end{array}$$(7)

and we can then write the change of energy in the gyration as

$$\begin{array}{c}\hfill \frac{\mathrm{d}(B\mu )}{\mathrm{d}t}=\mu (\frac{\partial B}{\partial t}+v_{\Vert }\mathit{b}·\mathrm{\nabla }B+{\mathit{v}}_E·\mathrm{\nabla }B).\end{array}$$(8)

Here, the second term is again associated with magnetic mirroring and does not alter the total energy. The first and third terms represent betatron acceleration. This is a total energy increase or decrease as a result of an increasing or decreasing magnetic field strength, either through explicit time variance of the magnetic field or due to a gradient perpendicular to the magnetic field direction. Betatron acceleration may be associated with the gradient B drift since

$$\begin{array}{c}\hfill \mu {\mathit{v}}_E·\mathrm{\nabla }B={\mathit{v}}_{\mathrm{\nabla }B}·q\mathit{E}.\end{array}$$(9)

2.2. The electromagnetic field and its gradients

In the experiments described in this section, we assumed the timescale of the test particle motion to be much shorter than the characteristic timescales of the electric and magnetic fields. That is, E and B were assumed static, and the explicit time derivative in Eq. (7) becomes zero. We obtained the static electromagnetic field from a snapshot of an MHD simulation produced by the Bifrost code (Gudiksen et al. 2011). Bifrost solves the radiative MHD equations using a sixth-order, finite difference differential operator. The operator stencil is asymmetric so that the derivatives become staggered. To minimise the number of interpolations (which are of fifth order), Bifrost places its variables on a staggered grid. Bifrost is used to simulate the solar atmosphere realistically and with high performance and precision, all the way from the corona down to the photosphere, sometimes also including parts of the convection zone.

We selected a snapshot from the Bifrost simulation by Færder et al. (2024) that uses hyper-diffusion to model resistivity. This time series simulates a 2D slice of an idealised coronal fan-spine magnetic configuration. There is no gravity nor radiative cooling or heating. The initial magnetic field has a 70 μT negative parasitic polarity (at the lower boundary) embedded in a 30 μT positive background field (in the z-direction). The initial magnetic field and plasma velocity have no y components and thus always lie in the xz plane. Consequently, there are no electric field components parallel to the magnetic field since the electric current and field will always lie in the y-direction. The plasma is treated as an ideal, ionised solar abundance gas, with a mean molecular weight of 0.616 amu. The initial plasma density and temperature are homogeneous with values 3 ⋅ 10⁻¹³ kg m⁻³ and 0.61 MK, respectively. At the lower boundary, a constant horizontal velocity is imposed that has a peak value at the location of the inner spine of the magnetic field, advecting the inner spine parallel to the lower boundary. This stretches the current sheet at the null point, which eventually becomes so thin and long that the plasmoid instability is reached. From halfway in the time-evolving simulation, at a point where two plasmoids are present, we increased the resolution using linear interpolation, and ran the simulation over a small time span. We took the electromagnetic field and electron density from the final snapshot of this rerun to use as static fields in the test particle simulations.

Having the electromagnetic field on a grid, we needed interpolation to get the field values at an arbitrary position in space. For this purpose, we used cubic spline interpolation. We also needed to calculate the spatial derivatives ∇B, ∇b, and ∇v_E. To do this, we used forward mode automatic differentiation locally, at each particle position, using the Julia package ForwardDiff (Revels et al. 2016). Having methods for evaluating the electromagnetic field and its gradients at arbitrary positions, we only needed the particle initial conditions to start simulating their trajectories.

2.3. Particle sampling

The method for choosing the initial conditions is crucial for finding non-thermal particles and inferring properties of the plasma. We wanted to study the statistical nature of the non-thermal particles (e.g. the probability of a particle becoming non-thermal from a specific set of initial conditions), so we needed the initial conditions of the test particles to be a statistically accurate representation of the particles in the solar atmospheric plasma. In our case, the particles needed to reflect the electron number density and the plasma temperature.

2.3.1. Sampling initial positions and velocities

We sampled the initial positions from the electron density distribution using rejection sampling. The rejection sampling algorithm is described in Appendix A. The initial velocity was sampled from a Maxwell-Boltzmann distribution, where the distribution’s temperature corresponds to the MHD temperature at the initial position of the particle. To obtain the Maxwell-Boltzmann distribution, we sampled each of the three velocity components independently from a Gaussian distribution. Finally, from the initial velocity components (v) and the initial position (r), we evaluated the initial guiding centre position, the initial parallel velocity, the initial gyration velocity, and the magnetic moment according to the formulas:

$$\begin{array}{cc}\hfill \mathit{R}& =\mathit{r}+\frac{m}{qB(\mathit{r})}\{[\mathit{v}-{\mathit{v}}_E(\mathit{r})]\times \mathit{b}(\mathit{r})\},\hfill \end{array}$$(10a)

$$\begin{array}{cc}\hfill v_{\Vert }& =\mathit{v}·\mathit{b}(\mathit{r}),\hfill \end{array}$$(10b)

$$\begin{array}{cc}\hfill {\mathit{v}}_{\perp }& =\mathit{v}-v_{\Vert }\mathit{b}(\mathit{r})-{\mathit{v}}_E(\mathit{r}),\hfill \end{array}$$(10c)

$$\begin{array}{cc}\hfill \mu & =\frac{m{v}_{\perp }^2}{2B(\mathit{r})}·\hfill \end{array}$$(10d)

Here we have assumed that the perpendicular drift at r is dominated by v_E.

2.3.2. Re-sampling initial positions

In the context of the entire simulation domain, non-thermal electrons are rare. It is very unlikely that an electron sampled from the electron density distribution will become non-thermal. This will lead us to simulate a very large number of ‘uninteresting’ electrons. To deal with this problem, we used the results of an initial, unbiased test particle run as an indication of where acceleration occurs. Specifically, we used the relative energy gain as a function of the initial position to create a probability density distribution from which we drew the initial positions of a second run. This re-sampling will be biased towards sites that accelerated particles in the first ensemble, and consequently increase the number of accelerated particles in the second ensemble. To maintain the statistical representation of the plasma, the particles in the second run are weighed by their importance weights (see Appendix B for details on importance sampling). The importance weight is defined as

$$\begin{array}{c}\hfill w(\mathit{r})=\frac{f(\mathit{r})}{g(\mathit{r})},\end{array}$$(11)

where r is the particle initial position, f is the electron density distribution, and g is the probability density distribution from which we drew the initial positions in the re-sampling.

2.4. Solving the equations of motion

To solve Eqs. (3a) and (3b) we used the Julia package DifferentialEquations (Rackauckas & Nie 2017). Specifically, we use a second and third order L-stable Rosenbrock-W method, with adaptive time-stepping. All particles were solved for the same time span, but they were stopped if they exited the simulation domain, or if the GCA was considered invalid. The latter was determined by evaluating the criterion

$$\begin{array}{c}\hfill r_L/L_B<\mathrm{tolerance},\end{array}$$(12)

which we call the GCA tolerance. The criterion was evaluated at each integration time step. If the ratio is larger than the tolerance, the particle was stopped. To investigate the effects of the tolerance value, we ran the test particle ensemble multiple times with varying tolerance.

3. Results

3.1. The MHD snapshot

The MHD snapshot is shown in Fig. 1. Here, the plasma density gradients illuminate the separatrices of the magnetic configuration. The 2D fan and inner spine are merging with the outer spine at the null point, where two plasmoids are visible as high-density disks and with closed magnetic field lines. The upper-left plasmoid is on its way out into the upper-left reconnection exhaust, whereas the lower-right plasmoid is growing in size while moving towards the opposite exhaust.

Fig. 1. Plasma density of the MHD snapshot where we embedded test particles. The black streamlines show the magnetic field direction, and the two left panels shows the reconnection region in greater detail. The grid size is 81922, and the resolution is 3.9 km.

3.2. Test particle simulation

We first simulated an ensemble of 10 million electrons for a duration of 1 second. The initial positions of the particles were distributed proportional to the electron density in the whole domain of the MHD snapshot in Fig. 1. Figure 2 shows the energy distribution of all electrons before and after the simulation. The two distributions are hardly distinguishable.

Fig. 2. Normalised energy distribution of 10 million electrons before and after the first test particle simulation. The initial positions were sampled from the electron density of the MHD snapshot in Fig. 1. The initial velocities were sampled from a Maxwell-Boltzmann distribution at their initial temperature. The non-thermal electrons are too rare to be in the sample of this initial ensemble.

To increase the chances of sampling non-thermal electrons, we re-sampled initial conditions and ran a second ensemble of 3 million electrons. The new initial conditions were sampled from a proposal distribution constructed from the relative energy gain as a function of initial position of the first ensemble. This re-sampling procedure was repeated two more times, with an ensemble size of 5 million and 13 million, respectively. The third and final proposal distribution is displayed in Fig. 3. Figure 4 shows the initial and final electron energy distribution after simulating the third re-sampled electron ensemble. The final energy distribution is shown both with and without statistical weighing. The bins of the unweighted distribution have a height that corresponds to the number of electrons that fall within that bin, while the bins in the weighted distribution have a height that corresponds to the sum of the weights of the electrons that fall within that bin. The final energy distributions are also shown with and without a filter that excludes all particles not in the acceleration region, defined as the rectangular region displayed by the upper panel in Fig. 5. In contrast to Figs. 2 and 4 shows clear non-thermal power laws.

Fig. 3. Median relative energy gain as a function of initial position of the third electron ensemble. The function is normalised so that it can be used as a probability density distribution. The function is practically zero outside the reconnection region. Values above 15 all have the same colour.

Fig. 4. Normalised energy distribution after the third re-sampling of 13 million electrons. In this ensemble, the initial positions were sampled using Fig. 3 as a proposal distribution. The figure shows the final energy distribution with and without statistical weighting to account for the bias introduced by sampling from the proposal distribution, instead of the electron density. It also shows the distributions after applying a filter that excludes electrons outside the acceleration region, defined as the region displayed by the upper panel in Fig. 5. The distribution for energies lower than 10 eV is not shown.

Fig. 5. Magnetic field strength and field line segments of the current sheet overplotted with the trajectories of the 608 electrons with a final energy greater than 1 keV (excluding any particles with an initial energy above 1 keV). Their initial positions are shown as white dots, and the final positions as black dots.

From the final energy distribution, we selected the accelerated electrons with a final energy greater than 1 keV and plotted their trajectories in Fig. 5. The initial positions are mostly clustered into two groups, located around an O-type and an X-type magnetic null point associated with the upper-left plasmoid of Fig. 1. Figure 6 plots the trajectory of one of these electrons in detail, showing the characteristic path. Its energy evolution is highlighted in Fig. 7, which shows the energy evolution of all non-thermal electrons. The electrons have a gradual increase in energy and many of them have plateaued or even decreased in energy before 1 s. The GCA tolerance used to produce the results in Figs. 4–6, was 10⁻². We re-ran the electron ensemble with a lower (10⁻³) and a higher tolerance (10⁻¹), and their respective final energy distributions are shown in Fig. 8.

Fig. 6. Trajectory of a betatron-accelerated non-thermal electron. The electron is initially upstream of the O-type null point, E × B drifting towards it. When it arrives downstream of the O-point, it drifts away from it.

Fig. 7. Time evolution of the kinetic energy of the 608 non-thermal electrons in the final test particle ensemble. The highlighted line corresponds to the trajectory in Fig. 6.

Fig. 8. Initial and final energy distribution with a varying tolerance on the GCA’s validity. All distributions are statistically weighted, and electrons outside the acceleration region are excluded. The distribution for energies lower than 10 eV is not shown.

4. Discussion

4.1. Non-thermal power law formation

The test particle simulations show sufficient electron acceleration to form a power law. However, we needed careful statistical sampling of initial positions to find high energy electrons. Out of the ∼10⁷ electrons simulated, only ∼10³ became non-thermal. In addition, the non-thermal particles have statistical weights much smaller than one, so when we in Fig. 4 look at the energy distribution of all electrons (not filtered), the weighted distribution differs significantly from the unweighted distribution. On the other hand, if we filter out the electrons outside the region around the acceleration sites, we get a significant power law even in the weighted distribution. The difference in the filtereddistributions is lower because the non-thermal electrons are more common in the acceleration region, and because the sampling and electron densities are more uniform. The non-thermal tail of the weighted and filtered energy distribution carries an energy fraction of approximately 2 ⋅ 10⁻⁴, which increases to 6 ⋅ 10⁻⁴ for the corresponding unweighted distribution. For the unfiltered distributions, the total energy fractions are 10⁻⁷ and 10⁻⁵, respectively. Here we defined the non-thermal tail as the distribution above 1 keV. A higher number of non-thermal electrons would increase the precision of the power law estimate, but the estimate would not necessarily be more accurate, since the power law formation is also heavily dependent on the GCAtolerance.

The fact is that all the non-thermal electrons emerge from areas where the GCA is questionable, around magnetic null points. Figure 8 shows that the power law hardens with increasing GCA tolerance, from no obvious power law to a power law with exponent −4. This suggests that the GCA is not sufficient to describe the motion of all potentially non-thermal electrons in this MHD configuration. To obtain a converging non-thermal power law, we would need to solve Eq. (1) directly in the areas of high r_L/L_B, similar to the hybrid schemes employed by Birn et al. (2004), Stanier (2013), Bacchini et al. (2020), and Pallister et al. (2021). The need for full particle motion is emphasised by the positive correlation in Fig. 9, where we have binned the final electron energies against the maximum r_L/L_B value they experienced during their orbit. This figure indicates that the less valid the GCA is, the more likely the electron is to be accelerated. For the non-thermal electrons in Fig. 5, the GCA is least valid when they are closest to the null points. This is a result of r_L/L_B being proportional to B⁻².

Fig. 9. Correlation of an electron’s relative energy gain with its highest experienced value of rL/LB. In this ensemble, the GCA tolerance was 10−2, hence the abrupt cutoff at rL/LB = 10−2. The colour maps the counts on a logarithmic scale. The non-thermal electrons in the cluster around the O-type null point (of Fig. 5) are shown in blue. They all fall within the diagonal strand in the upper-right part of the histogram, where the relative energy gain is above 1 and rL/LB is above 10−4. The non-thermal electrons in the cluster around the X-type null point are shown in green and fall within the more diffuse diagonal strand below.

4.2. Acceleration mechanisms

The non-thermal electrons of Fig. 5 drift away from their initial positions close to the magnetic null points due to the E × B drift. At the same time they are mirrored back and forth due to parallel gradients in the magnetic field strength. This is similar to the bouncing particles in the turbulent coronal loop-top simulation by Bacchini et al. (2024) where the electrons are trapped by magnetic mirroring and the mirroring represents the dominating parallel acceleration. However, the total energy increase is caused by betatron acceleration as the magnetic field strengthens when the electrons drift away from the null point, making Eq. (9) positive. The mirroring acceleration itself does not alter the total energy. The electrons with a trajectory similar to the one in Fig. 6 all start upstream of the O-type null point of the largest plasmoid. Here, the electric field points out of the plane, so that the electrons initially drift towards the O-point. When the magnetic field strength is too weak for magnetic mirroring, they orbit to the other side of the null point where the electric field points into the plane. From this point, their trajectories proceed like the electrons initialised downstream of the O-type null point, drifting away and betatron-accelerated.

There are no parallel electric fields, so the first term of Eq. (3b) is always zero, but there is curvature in the magnetic field direction along the electron path, making the third term of Eq. (3b) (the first-order Fermi acceleration term) non-zero. The energy change from this curvature acceleration is positive for the upper-left electron cluster in the upper panel of Fig. 5 because here the curvature is parallel to v_E, making Eq. (5) positive. It is negative for the lower-right cluster, where the curvature is anti-parallel to v_E. Both clusters experience a gradual increase in total energy, indicating that betatron acceleration dominates over the first-order Fermi acceleration. The plateauing of the energy evolution that can be seen in Fig. 7 can be explained by the fact that v_E ⋅ ∇B gradually reduces and changes sign in the middle between the two null points. The reduction in v_E ⋅ ∇B is true for both clusters, but the sign change is only true for the lower-right cluster.

Standing out from the clusters, there is one non-thermal electron in Fig. 5 that maintains a circular orbit and gains energy mainly from first-order Fermi acceleration. However, the acceleration is not very efficient as the relative energy gain is only 8%. An efficient Fermi acceleration could be achieved closer to the O-type null point where the curvature is stronger. This requires that the electron has an initial position close to the null point, or a sufficient E × B drift towards it. On the other hand, orbits close to the null point – where the characteristic magnetic field length is short and where the orbit size is comparable to the snapshot resolution – are prone to numerical errors. The errors can produce slightly spiralling trajectories, and if the electron spirals towards the centre of the O-type null point, it experiences artificially strong first-order Fermi acceleration. For this reason, we filtered out the electrons that have their final position within a box of size 12.7 km × 4.5 km and 4.5 km × 4.5 km around the O-type null points associated with the upper and lower plasmoid, respectively. Of the 13 million electrons in the ensemble, 17 thousand end up in these boxes.

In general, all non-thermal electrons are mainly accelerated from initial positions in the current sheet region. Here there are large gradients in the magnetic field and a significant electric field, which are required for both betatron and first-order Fermi acceleration (see Eqs. (5) and (9)). Figure 3 shows that the change in energy is particularly correlated with the larger upper-left plasmoid, which has started to merge with the reconnection exhaust. This plasmoid has weaker gradients than the rest of the current sheet, but has a stronger electric field.

When it comes to the electron acceleration as a function of initial energy, we know from Eq. (8) that the energy change due to betatron acceleration is proportional to the magnetic moment of the particle, so that a high initial energy at low pitch angles increases the potential of strong betatron acceleration. And in Fig. 10, plotting the correlation between initial and final energy, we do see a tendency towards acceleration for high initial energies. All non-thermal electrons have an initial energy around or above 100 eV.

Fig. 10. Correlation between the initial (E0) and final energy (Ef) of the electrons in the final test particle ensemble.

In Fig. 11 we compare the non-thermal electron of Fig. 6 with its equivalent direct solution of the Lorentz force. The comparison indicates that the betatron acceleration is, in this case, dependent on the conservation of magnetic moment. For all of the electrons investigated, a direct solution fails to reproduce the total energy increase achieved with the GCA. Although the trajectories are similar, the directly solved electrons lose magnetic moment, which decreases the average pitch angle and consequently increases the mirror length. The small reduction in kinetic energy may be first-order Fermi de-acceleration due to the curvature being anti-parallel to the E × B drift. This type of evolution is observed with several of the betatron-accelerated, non-thermal electrons, independent of the value of r_L/L_B. The energy gain in the GCA seems to be correlated with the variance of magnetic moment seen in the equivalent direct solution. For electrons that experience low r_L/L_B, the magnetic moment is closer to being conserved, and the energy gain is closer to zero. The electrons that experience a high r_L/L_B gain more energy, but the magnetic moment is more variant. This correlation could explain the dependence of the non-thermal power law on the GCA tolerance (see Fig. 8).

Fig. 11. Comparison between the GCA and a direct solution to the Lorentz equation. The blue GCA trajectory is the same as in Fig. 6, and the orange trajectory is its equivalent full orbit. The lower panels show the pitch angle time evolution for both methods, as well as the time evolution of the magnetic moment (relative to the initial magnetic moment) and the kinetic energy for the direct solution.

5. Conclusion

We have developed a test particle code using the GCA and transported electrons in a static, 2D coronal, fan-spine configuration produced by the MHD code Bifrost. By using importance sampling in multiple stages, we were able to detect the formation of a non-thermal power law energy distribution from a statistically representative sample of electrons. Without the statistical weighing of the electrons, the energy fraction in the non-thermal power law increased by up to 2 orders of magnitude, indicating that the electron sampling technique has a significant impact when estimating non-thermal energy distributions using test particles.

Furthermore, the shape of the non-thermal power law heavily depends on our GCA tolerance. This is partially because the GCA breaks down when the electrons get too close to a magnetic null point. A higher tolerance allowed us to simulate electrons with an initial position closer to the null points, and since betatron acceleration is proportional to the magnetic moment, which is inversely proportional to the initial magnetic field strength, it had a large effect on the number of non-thermal electrons. We also find a positive correlation between the tolerance and electron energy gain in general. Comparisons of the guiding centre trajectories with direct solutions to the Lorentz force indicated that the invariance of the magnetic moment was key to the betatron acceleration and that the energy gain was stronger for less adiabatic orbits. The sensitivity of the non-thermal power law to GCA tolerance, and the computational infeasibility of simulating particle ensembles with direct integration of the Lorentz force, suggests that a hybrid scheme is necessary to get a converging estimate of the true non-thermal power law in this static 2D MHD configuration. In addition to these conclusions, the following points should be consideredin future work.

The particles were followed for 1 s. The MHD timescales in and around the current sheet are of the same order. Assuming the electric and magnetic fields are static for 1 s is therefore at the limit of what can be considered to be tolerable. A dynamic electromagnetic field may also help particles escape plasmoids, since the field lines are allowed to diffuse and reconnect over time.

The filtering of electrons close to O-type magnetic null points excludes potential non-thermal particles. These particles can be simulated more accurately using an interpolation scheme that ensures a divergence-free magnetic field, such as the one in Mackay et al. (2006). In addition, a hybrid scheme switching from GCA to full orbit integration is necessary when the particle gets too close to the magnetic null point. Non-adiabatic, O-type null point acceleration of the type discussed in Zhao et al. (2021) would then be possible, although only for very energetic electrons since for this type of acceleration, r_L should be larger than the range of motion in the xz plane, which means r_L should be larger than the MHD grid size of 3.9 km.

There is also the question of whether sustained, trapped, circular (or mirrored) motion around O-type null points will occur in more realistic 3D simulations. A third dimension will increase the complexity of the magnetic field, allowing the particles to escape through an additional dimension. A 3D version of the MHD simulation would also better facilitate the escape of accelerated electrons along the open field lines of the outer spine, and along the fan and inner spine towards the footpoints, as in Pallister et al. (2021). Moreover, 3D makes it easier for the particles to jump from one plasmoid or acceleration site to another, facilitating acceleration in consecutive stages.

Acknowledgments

This research is supported by the Research Council of Norway through its Centres of Excellence scheme, project number 262622. We thank Jan Karsten Trulsen for invaluable discussions on particle sampling. The simulation of test particle and data analysis were performed using the Julia programming language (Bezanson et al. 2017). The Figures were produced using the open-source Julia-package Makie (Danisch & Krumbiegel 2021). The authors thank the referee for valuable comments and suggestions.

References

Appendix A: Rejection sampling

To sample initial particle positions from the electron density, which is a 2D inhomogeneous distribution always greater than zero (here called the target distribution f(x)), we used rejection sampling (see e.g. Sect. 7.3.6 in Press et al. 2007). We drew our samples x_i from a uniform probability density distribution with limits identical to the target distribution, and accept the samples with a probability proportional to f(x_i). In this way, our accepted draws will reflect the shape of the target distribution, at the cost of rejected draws. The algorithm goes as follows:

1. Draw a random number u from the uniform probability distribution U(x_a, x_b). The limits x_a and x_b must coincide with the x-limits of the target distribution f = f(x, y).

2. Draw another random number v ∈ U(y_a, y_b), where the limits y_a and y_b coincide with the y-limits of f.

3. Draw a third random number w ∈ U(0, c), where the constant c is equal to or greater than the maximal value of f(x, y).

4. If w ≤ f(u, v) accept the sample (u, v). Otherwise, reject (u, v).

5. Repeat the above steps until you have the desired number of accepted samples.

The ratio of accepted draws to rejected draws strongly depend on the shape of the target distribution. For an optimal choice of c, the acceptance ratio approaches 1 as the target distribution approaches uniform. At the other extreme, the acceptance rate would be zero for a Dirac delta function. To approximate the target distribution from the drawn samples, we had to bin the samples in a histogram. The accuracy of each bin is inversely proportional to the square root of the number of samples in that bin.

Appendix B: Importance sampling

In importance sampling (see e.g. Sect. 7.9.1 in Press et al. 2007) we draw samples from a selected proposal distribution instead of the target distribution. To recreate the target distribution, we multiply the samples with importance weights. The more probable the sample is in the target distribution, the higher the weight, and the more probable the sample is in the proposal distribution, the lower the weight. The proposal distribution is required to cover the same domain as the target distribution, but is otherwise arbitrary, only limited by the feasibility of sampling from it. The benefit of importance sampling is two-folded. First, the proposal distribution may be easier or faster to sample from (e.g. a uniform distribution). Second, biasing the sampling towards values that otherwise would be rarely sampled allows us to increase the statistical accuracy of the target distribution at those locations. As an example, Fig. B.1 displays importance sampling of a Maxwell-Boltzmann probability density distribution. Direct sampling will estimate the high-value tail of the distribution with poor accuracy, since this part of the distribution is rare. A well-placed proposal distribution, like a Gaussian shifted towards the tail, will significantly increase the accuracy of the high-value estimates, at the cost of the accuracy of the lower value estimates.

Fig. B.1. Showcase of importance sampling. The upper panel shows a histogram of 2000 samples directly sampled from a Maxwell-Boltzmann probability density distribution, f(x). In the lower panel, the same amount of samples are sampled from a Gaussian probability density distribution, g(x), which is shifted towards higher x values compared to the Maxwell-Boltzmann distribution. This shift increases the fraction of samples with high x-values compared to the direct sampling, and if the samples are weighted with the ratio f(x)/g(x) when binned, the resulting histogram will estimate the Maxwellian-probability distribution, but this time with higher statistical accuracy for the higher values of x.

B.1. Calculating the importance weights

Consider the arbitrary function h(x) and the probability density distribution f(x). The integral

$$\begin{array}{c}\hfill \mathcal{I}={\displaystyle \int h}(x)f(x)\mathrm{d}x=E^f[h(\mathrm{x})]\end{array}$$(B.1)

is the expectation value of h(x) in f, where x is a random variable that is distributed according to f(x). The integral can be approximated by the sample mean,

$$\begin{array}{c}\hfill \mathcal{I}=E^f[h(\mathrm{x})]\approx \frac{1}{n}{\displaystyle \sum _i^n}h(x_i),\end{array}$$(B.2)

where x₁, …, x_n are realisations (samples) of x, and n is the number of these realisations. Next, we introduce the proposal probability density distribution g(x) into the integral,

$$\begin{array}{c}\hfill \mathcal{I}={\displaystyle \int h}(x)f(x)\frac{g(x)}{g(x)}\mathrm{d}x=E^g\left[h(\mathrm{x})\frac{f(\mathrm{x})}{g(\mathrm{x})}\right].\end{array}$$(B.3)

In this case, ℐ is formulated as the expectation value of h(x)f(x)/g(x) in g, and x is a random variable distributed according to g(x). We can then approximate the integral with a different sample mean, namely

$$\begin{array}{c}\hfill \mathcal{I}\approx \frac{1}{m}{\displaystyle \sum _j^m}h(x_j)\frac{f(x_j)}{g(x_j)},\end{array}$$(B.4)

where x₁, …, x_m are realisations of the random variable x, distributed according to g(x). The ratio f(x_j)/g(x_j) is the importance weight of the sample h(x_j). In the example in Fig. B.1, h(x) = x, f(x) is the Maxwell-Boltzmann target distribution, and g(x) is the Gaussian proposal distribution.

All Figures

	Fig. 1. Plasma density of the MHD snapshot where we embedded test particles. The black streamlines show the magnetic field direction, and the two left panels shows the reconnection region in greater detail. The grid size is 81922, and the resolution is 3.9 km.
In the text

	Fig. 2. Normalised energy distribution of 10 million electrons before and after the first test particle simulation. The initial positions were sampled from the electron density of the MHD snapshot in Fig. 1. The initial velocities were sampled from a Maxwell-Boltzmann distribution at their initial temperature. The non-thermal electrons are too rare to be in the sample of this initial ensemble.
In the text

	Fig. 3. Median relative energy gain as a function of initial position of the third electron ensemble. The function is normalised so that it can be used as a probability density distribution. The function is practically zero outside the reconnection region. Values above 15 all have the same colour.
In the text

Fig. 4. Normalised energy distribution after the third re-sampling of 13 million electrons. In this ensemble, the initial positions were sampled using Fig. 3 as a proposal distribution. The figure shows the final energy distribution with and without statistical weighting to account for the bias introduced by sampling from the proposal distribution, instead of the electron density. It also shows the distributions after applying a filter that excludes electrons outside the acceleration region, defined as the region displayed by the upper panel in Fig. 5. The distribution for energies lower than 10 eV is not shown.

In the text

	Fig. 5. Magnetic field strength and field line segments of the current sheet overplotted with the trajectories of the 608 electrons with a final energy greater than 1 keV (excluding any particles with an initial energy above 1 keV). Their initial positions are shown as white dots, and the final positions as black dots.
In the text

	Fig. 6. Trajectory of a betatron-accelerated non-thermal electron. The electron is initially upstream of the O-type null point, E × B drifting towards it. When it arrives downstream of the O-point, it drifts away from it.
In the text

	Fig. 7. Time evolution of the kinetic energy of the 608 non-thermal electrons in the final test particle ensemble. The highlighted line corresponds to the trajectory in Fig. 6.
In the text

	Fig. 8. Initial and final energy distribution with a varying tolerance on the GCA’s validity. All distributions are statistically weighted, and electrons outside the acceleration region are excluded. The distribution for energies lower than 10 eV is not shown.
In the text

Fig. 9. Correlation of an electron’s relative energy gain with its highest experienced value of rL/LB. In this ensemble, the GCA tolerance was 10−2, hence the abrupt cutoff at rL/LB = 10−2. The colour maps the counts on a logarithmic scale. The non-thermal electrons in the cluster around the O-type null point (of Fig. 5) are shown in blue. They all fall within the diagonal strand in the upper-right part of the histogram, where the relative energy gain is above 1 and rL/LB is above 10−4. The non-thermal electrons in the cluster around the X-type null point are shown in green and fall within the more diffuse diagonal strand below.

In the text

	Fig. 10. Correlation between the initial (E0) and final energy (Ef) of the electrons in the final test particle ensemble.
In the text

Fig. 11. Comparison between the GCA and a direct solution to the Lorentz equation. The blue GCA trajectory is the same as in Fig. 6, and the orange trajectory is its equivalent full orbit. The lower panels show the pitch angle time evolution for both methods, as well as the time evolution of the magnetic moment (relative to the initial magnetic moment) and the kinetic energy for the direct solution.

In the text

Fig. B.1. Showcase of importance sampling. The upper panel shows a histogram of 2000 samples directly sampled from a Maxwell-Boltzmann probability density distribution, f(x). In the lower panel, the same amount of samples are sampled from a Gaussian probability density distribution, g(x), which is shifted towards higher x values compared to the Maxwell-Boltzmann distribution. This shift increases the fraction of samples with high x-values compared to the direct sampling, and if the samples are weighted with the ratio f(x)/g(x) when binned, the resulting histogram will estimate the Maxwellian-probability distribution, but this time with higher statistical accuracy for the higher values of x.

In the text