© ESO, 2016

1. Introduction

Understanding the dynamics of charged particles in reconnecting solar processes, in particular the acceleration of particles in solar flares, is of fundamental importance for advancing our knowledge of plasma behaviour in general (for example, recent reviews on particle acceleration in flares include Vilmer et al. 2011; Fletcher et al. 2011; Zharkova et al. 2011; Cargill et al. 2012). There has been growing observational evidence that accelerated particle populations also occur in smaller scale processes such as micro-flares (see e.g. Hannah et al. 2011, and references therein) or nano-flares (Testa et al. 2014); these processes are also often associated with magnetic reconnection.

In the context of the solar atmosphere, it is now widely accepted that magnetic reconnection contributes to the maintenance of thermal (contributing to coronal heating as discussed in Parnell & De Moortel 2012; Parnell et al. 2015) and sporadic generation of non-thermal (high-energy) particle populations. Since magnetic reconnection is intrinsically associated with field aligned electric fields (e.g. Schindler et al. 1988, 1991; Hesse & Schindler 1988) and these electric fields are well-known accelerators of charged particles, the question of how the local particle dynamics on all scales are affected by these reconnection electric fields remains open and important. Additionally, magnetic reconnection can lead to rapid magnetic field changes that can also cause particle acceleration, such as collapsing magnetic traps (e.g. Grady & Neukirch 2009; Grady et al. 2012; Eradat Oskoui et al. 2014; Eradat Oskoui & Neukirch 2014).

Previous studies of particle behaviour typically use a test particle approach to study a wide variety of configurations, all of which are underpinned by or directly investigate magnetic reconnection. For instance, the acceleration of particles at specific topological features has been explored in studies at/near 2/2.5D null points (e.g. Bulanov & Sasorov 1976; Bruhwiler & Zweibel 1992; Kliem 1994; Litvinenko 1996; Browning & Vekstein 2001; Zharkova & Gordovskyy 2004, 2005; Wood & Neukirch 2005; Hannah & Fletcher 2006; Drake et al. 2006) and, more recently, 3D null points (e.g. Dalla & Browning 2005, 2006, 2008; Guo et al. 2010; Stanier et al. 2012), stressed magnetic fields/current sheets (e.g. Turkmani et al. 2005; Onofri et al. 2006; Gordovskyy et al. 2010) and 3D magnetic separators (Threlfall et al. 2015, 2016). Acceleration in more complex 3D models, such as twisted coronal loops (Gordovskyy & Browning 2011; Gordovskyy et al. 2014), have also been investigated. A significant step forward from the analysis of a localised reconnection site was taken by Baumann et al. (2013), who combined localised particle-in-cell (PIC) analysis in the vicinity of a reconnecting null-point with a magnetohydrodynamic (MHD) model of the surrounding active region (AR). This allowed the global impact sites of the particles to be identified.

Instead of considering an isolated reconnection process, whether embedded in a global field or not, we take a different approach and investigate the particle acceleration in an observationally driven MHD simulation of an AR. A decade ago, the first MHD simulations of downscaled ARs were produced that were based on extrapolated observed magnetograms, driven in a manner which approximated motions arising from convection cells on a range of scales (e.g. Gudiksen & Nordlund 2005a,b). Only recently have large-scale computations been performed that model a full AR without downscaling (Bourdin et al. 2013).

The MHD simulation we will use was designed to model the coronal response to photospheric driving by advection of field-lines resulting from horizontal granular motions in the photosphere, where the plasma beta is greater than one (Bourdin et al. 2013). The base of the simulation is prescribed by an observed AR magnetogram with the subsequent field-line braiding inducing currents in the corona caused by quasi-static changes propagating at the Alfvén speed along the field into the corona. In effect, this is one implementation of the field-line braiding mechanism of Parker (1972) and results in coronal loops with temperatures of around 1.5 MK. The location, apex height, and plasma flows along these loops all broadly match observations (Bourdin et al. 2013). The feasibility of this field-line braiding mechanism to heat the corona sufficiently by Ohmic dissipation of induced currents is demonstrated in Bourdin et al. (2015). This suggests that reconnection occurs regularly throughout the AR, and, hence, there must be widespread parallel electric fields.

In this paper, we investigate what, if any, particle acceleration occurs in this MHD simulation of an observed non-flaring AR (Bourdin et al. 2013). In this kind of model there are a wide range of locations (magnetic configurations) at which reconnection occurs, as well as a range of different reconnection regimes in operation. Thus, the MHD model we investigate here incorporates many of the particle acceleration scenarios that have been considered in isolation, as mentioned above. Thus, in this paper we consider the electromagnetic fields that have been generated in a non-flaring active-region simulation based on observations and determine the effects they have on charged particle dynamics.

Fig. 1 Panels, from L-R, showing a) a Hinode/SOT magnetogram of an AR core seen on the 14th Nov. 2007 saturated at ± 300 G; b) a Hinode/EIS image of the Fe xv emission (at 1.5 MK) of the corona above this region and c) the co-spatial and co-temporal simulated parallel electric field strength and orientation within this region (where pink/purple represents positive/negative E||, and the opacity represents the relative strength of E||− fully opaque regions identify ${\mathit{E}}_{\mathrm{\left|}\mathrm{\right|}}\mathrm{=}\mathrm{0.25}{\mathit{E}}_{\mathrm{\left|}\mathrm{\right|}}^{\mathrm{peak}}$ while transparent regions show ${\mathit{E}}_{\mathrm{\left|}\mathrm{\right|}}\mathrm{=}\mathrm{0.01}{\mathit{E}}_{\mathrm{\left|}\mathrm{\right|}}^{\mathrm{peak}}$, for a peak coronal electric field strength ${\mathit{E}}_{\mathrm{\left|}\mathrm{\right|}}^{\mathrm{peak}}\mathrm{~}\mathrm{3.5}\mathrm{V}\hspace{0.17em}{\mathrm{m}}^{-1}$ in the simulation domain). The right hand panel is generated by simulations described in Bourdin et al. (2013), which also describes field lines which cross the intensity maxima of several EUV-emissive coronal loops, which have been overlaid in each panel; “CL” indicates core loops, while “SL” indicates shorter loops. The circles indicate EUV-emissive loop footpoints in the model and the corresponding plasma motion (and observed Doppler shifts) at these locations; red represents draining (i.e. motion along the line-of-sight away from an observer) while blue represents upflows (motion towards an observer).

The primary objective of the present work is to survey the particle dynamics which result from MHD simulations of a non-flaring AR, and to determine the extent and locations to which particles may be accelerated by this type of model. The paper is organised as follows: in Sect. 2 we discuss the model itself, introducing both the MHD simulation domain (described in Sect. 2.1) and the equations which govern test particle motion (in Sect. 2.2). A survey of orbit behaviour throughout the AR core is outlined and discussed in Sect. 3, before examining specific test particle examples in detail in Sect. 4. A discussion of our findings is presented in Sect. 5 before conclusions and future areas of study are outlined in Sect. 6.

2. Model setup

Our model can be broadly split into two parts: a time-dependent MHD simulation (modelling the solar corona above an AR) on time-scales of hours into which we insert particles (using a single evolved snapshot), and the test particle motion itself on time scales of seconds. A brief overview of both parts are described in the following sections:

2.1. MHD active region model

The AR MHD simulation of Bourdin et al. (2013) that we use is observationally driven; a potential-field extrapolation from an observed line-of-sight (LOS) magnetogram on the base determines the initial magnetic field state within the numerical domain. During the simulation the magnetic field is self-consistently computed within the domain, while the lower boundary is evolved in line with the observed time-series of LOS magnetograms. These magnetograms were obtained from the Hinode solar observatory (Kosugi et al. 2007; Culhane et al. 2007; Golub et al. 2007). The region in question is located close to the centre of the solar disk on the 14th November 2007; while not assigned an NOAA AR number, the observed region displayed many typical AR characteristics (for example, a system of closed loop structures, visible in both extreme ultraviolet (EUV) and X-ray wavelengths, linking a pair of strong magnetic field patches of opposite polarity; see Fig. 1).

While data from the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI, Lin et al. 2002) was unavailable for this date, the Geostationary Operational Environmental Satellite (GOES) X-ray data suggests that solar activity at this time was particularly low. The most significant GOES event on the 14th November (and indeed for the preceding fortnight) was an A1-class event, while the emergence of a second AR (NOAA 10974) coincided with a B1-class event on the 17th November. From the available data, we characterise the chosen region as a relatively isolated, stable, and slowly evolving AR which (crucially) does not produce a flare within one day before and two days after the observation. Hence, it was an ideal choice for the original aim of the simulation, namely to study coronal heating.

The computational domain is periodic in the horizontal directions. To sufficiently isolate the main AR polarities from their periodic images, the AR observation is surrounded by a region of quiet Sun. The full horizontal extent of 235 × 235 Mm² is covered by 1024² grid points; approximately a quarter of this area covers the AR core (see Fig. 1, which shows this core region alone). Thus the grid cells in the xy directions are square with lengths of 230 km on either side.

The initial model atmosphere contains vertical stratification and is initialised smoothly to avoid spurious oscillations (see Bourdin et al. 2014b, for more details). Thus the domain has a stretched grid in the vertical direction that reaches up to 156 Mm above the photosphere. The vertical extent of each cell covers between 100 km, near the photosphere, up to a maximum of 800 km in the corona. For the top boundary, a potential-field extrapolation allows closed and “open” field lines to relax into a force-free state.

The compressible resistive MHD simulation was performed using the Pencil Code (Brandenburg & Dobler 2002) and incorporates gravity, viscosity, and an isotropic magnetic diffusivity η(= 1 /μσ) = 10¹⁰ m²/ s (for electrical conductivity σ and permeability of free space μ; more details may be found in Bourdin et al. 2013). This value achieves diffusion of the generated current density structures at the grid scale, and is equivalent to setting the magnetic Reynolds number around unity for length scales on the order of the simulation grid scale (which is necessary to achieve numerical stability).More details of the implementation and selection of finite resistive effects for this (and related) model(s) may be found in e.g. Galsgaard & Nordlund (1996), Bingert & Peter (2011), Peter (2015). We address implications of the chosen parameter regime for the recovered particle behaviour in Sect. 5. An image highlighting the distribution of the parallel electric field in a single snapshot of the simulation domain can be seen in Fig. 1. Parallel electric fields are a necessary and sufficient condition for magnetic reconnection, which dissipates electric currents via Ohmic heating.

The initial magnetic field is evolved according to the advection of photospheric field-lines using a time-series of LOS magnetograms and caused by a driver on the base which simulates horizontal granular motions in the photosphere. In this region the plasma beta is greater than one, but falls off in height as the density falls, such that the plasma beta is less than one in the corona. During the simulation run the magnetic field evolves self-consistently according to compressible resistive MHD, and subject to an energy equation that incorporates both heat conduction and radiative losses.

In particular, the energy input in the model results from the Poynting flux across the base, which provides approximately 10⁷ erg^-1 cm^-2 (vertical signed flux average) that passes through the transition-region layer and reaches the corona. This balances the expected coronal energy requirement of an AR core region by (partial) conversion of magnetic energy into heat, which forms EUV-emissive coronal loops (Bourdin et al. 2014a). A Spitzer-type heat conduction term is included in the energy equation to model energy losses in the corona due to fast electrons (Spitzer 1962). In addition, X-ray and EUV emission is implemented through a radiative loss function which is parametrised by the local plasma density and temperature (Cook et al. 1989).

The lack of flaring activity and the reported agreement with observations makes this simulation an ideal candidate to study particle behaviour in a realistic solar configuration.

2.2. Relativistic particle dynamics

Having established the global environment which we will study, all that remains is to outline the equations which will govern particle behaviour. Despite the lack of flaring activity, we anticipate particle velocities may achieve values of a significant fraction of the speed of light (c). We therefore make use of the full relativistic set of guiding-centre-motion equations, outlined in Northrop (1963), based on the treatment of Vandervoort (1960), presented here in normalised form: $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{u}}_{\mathrm{\parallel }}}{\mathrm{d}\mathit{t}}& \mathrm{=}& \frac{\mathrm{d}}{\mathrm{d}\mathit{t}}\left(\mathit{\gamma }{\mathit{v}}_{\mathrm{\parallel }}\right)\mathrm{=}\mathit{\gamma }{{u}}_{\mathit{E}}\mathrm{·}\frac{\mathrm{d}{b}}{\mathrm{d}\mathit{t}}\mathrm{+}{\mathrm{\Omega }}_{\mathit{s}\mathrm{cl}}{\mathit{t}}_{\mathrm{scl}}{\mathit{E}}_{\mathrm{\parallel }}\mathrm{-}\frac{{\mathit{\mu }}_{\mathit{r}}}{\mathit{\gamma }}\frac{\mathit{\partial }{\mathit{B}}^{\mathit{\star }}}{\mathit{\partial s}}\mathit{,}\\ {{Ṙ}}_{\mathrm{\perp }}& \mathrm{=}& {{u}}_{\mathit{E}}\mathrm{+}\frac{{b}}{{\mathit{B}}^{\mathit{\star }\mathit{\star }}}\mathrm{\times }\{\frac{\mathrm{1}}{{\mathrm{\Omega }}_{\mathit{s}\mathrm{cl}}{\mathit{t}}_{\mathrm{scl}}}[\frac{{\mathit{\mu }}_{\mathit{r}}}{\mathit{\gamma }}\left({\nabla }{\mathit{B}}^{\mathit{\star }}\mathrm{+}\frac{{{\mathit{v}}_{\mathrm{scl}}}^{\mathrm{2}}}{{\mathit{c}}^{\mathrm{2}}}{{u}}_{\mathit{E}}\frac{\mathit{\partial }{\mathit{B}}^{\mathit{\star }}}{\mathit{\partial t}}\right)\\ & & \mathrm{+}{\mathit{u}}_{\mathrm{\parallel }}\frac{\mathrm{d}{b}}{\mathrm{d}\mathit{t}}\mathrm{+}\mathit{\gamma }\frac{\mathrm{d}{{u}}_{\mathit{E}}}{\mathrm{d}\mathit{t}}]\mathrm{+}\frac{{{\mathit{v}}_{\mathrm{scl}}}^{\mathrm{2}}}{{\mathit{c}}^{\mathrm{2}}}\frac{{\mathit{u}}_{\mathrm{\parallel }}}{\mathit{\gamma }}{\mathit{E}}_{\mathrm{\parallel }}{{u}}_{\mathit{E}}\}\mathit{,}\\ \frac{\mathrm{d}\mathit{\gamma }}{\mathrm{d}\mathit{t}}& \mathrm{=}\frac{{{\mathit{v}}_{\mathrm{scl}}}^{\mathrm{2}}}{{\mathit{c}}^{\mathrm{2}}}\left[{\mathrm{\Omega }}_{\mathit{s}\mathrm{cl}}{\mathit{t}}_{\mathrm{scl}}\left({{Ṙ}}_{\mathrm{\perp }}\mathrm{+}\frac{{\mathit{u}}_{\mathrm{\parallel }}}{\mathit{\gamma }}{b}\right)\mathrm{·}{E}\mathrm{+}\frac{{\mathit{\mu }}_{\mathit{r}}}{\mathit{\gamma }}\frac{\mathit{\partial }{\mathit{B}}^{\mathit{\star }}}{\mathit{\partial t}}\right]\mathit{,}& \\ {\mathit{\mu }}_{\mathit{r}}& \mathrm{=}& \frac{{\mathit{\gamma }}^{\mathrm{2}}{\mathit{v}}_{\mathrm{\perp }}^{\mathrm{2}}}{\mathit{B}}\mathrm{·}\end{array}$Here μ_r is the relativistic magnetic moment, for a particle with rest-mass m₀ and charge q, whose guiding centre is located at R, subject to an electric field E and a magnetic field B (with magnitude B(= | B |) and unit vector b(= B/B)). Local conditions will dictate aspects of the orbit behaviour, particularly through guiding centre drifts; the largest in magnitude is typically the E × B drift, which has a velocity u_E(= E × b/B). The component of velocity parallel to the magnetic field is ${\mathit{v}}_{\mathrm{\parallel }}\mathrm{(}\mathrm{=}{b}\mathrm{·}\dot{{R}}\mathrm{)}$, while E_∥(= b·E) is the magnitude of the electric field parallel to the local magnetic field, $\dot{{R}}{}_{\mathrm{\perp }}\mathrm{(}\mathrm{=}\dot{{R}}\mathrm{-}{\mathit{v}}_{\mathrm{\parallel }}{b}\mathrm{)}$ is the component of velocity perpendicular to b, and s is a line element parallel to b. Finally, γ is the Lorentz factor ${}^{\mathrm{(}}\mathit{\gamma }^{\mathrm{2}}\mathrm{=}\mathrm{1}\mathit{/}{\mathrm{(}}^{}\mathrm{1}\mathrm{-}{\mathit{v}}^{\mathrm{2}}\mathit{/}{\mathit{c}}^{\mathrm{2}}{\mathrm{)}}^{}\mathrm{=}{\mathit{c}}^{\mathrm{2}}\mathit{/}{{}^{\mathrm{(}}\mathit{c}^{\mathrm{2}}\mathrm{-}{\mathit{v}}^{\mathrm{2}}{\mathrm{)}}^{}}^{\mathrm{)}}$. Using this factor, we define a relativistic parallel velocity u_∥(= γv_∥) for simplicity of notation.

Equations (11a) have been expressed in dimensionless form, using a field strength b_scl, lengthscale l_scl and timescale t_scl. Dimensionless quantities are related to their dimensional equivalent through: $\begin{array}{ccc}\overline{{B}}\mathrm{=}{\mathit{b}}_{\mathrm{scl}}\hspace{0.17em}{B}\mathit{,} \mathit{x̅}\mathrm{=}{\mathit{l}}_{\mathrm{scl}}\hspace{0.17em}\mathit{x,} \mathit{t̅}\mathrm{=}{\mathit{t}}_{\mathrm{scl}}\hspace{0.17em}\mathit{t,}& & \end{array}$where the barred quantities represent dimensional counterparts of the variables seen in Eqs. (11a). This choice of quantities fixes the remaining normalising constants; for example, velocities in the model are scaled by v_scl(= l_sclt_scl^-1), energies by KE_scl(= 0.5mv_scl²) and (assessing the dimensions of Faraday’s Law) electric fields are scaled by e_scl(= b_scll_sclt_scl^-1 = b_sclv_scl). This investigation is motivated by the behaviour of particles in a solar coronal environment. Equations (11a) have been normalised accordingly; unless otherwise stated, all experiments take b_scl = 0.001 T, l_scl = 100 km and t_scl = 10 s.

To further simplify Eqs. (11a), only electrons or protons are considered here; this fixes the rest mass m₀ = m_e = 9.1 × 10^-31 kg and charge q = e = −1.6022 × 10^-19 C for electrons, or m₀ = m_p = 1.67 × 10^-27 kg and q = | e | = 1.6022 × 10^-19 C for protons. In this way, several normalising constants are expressed in terms of a normalising electron/proton gyro-frequency, Ω_scl(= qb_sclm₀^-1). The factor of Ω_sclt_scl controls the scales at which certain guiding centre drifts become important.

Finally, several terms in Eqs. (11a) now also depend on the ratio of perpendicular electric field (E_⊥) to the magnitude of the magnetic field (B); for a given magnetic field strength B, B^⋆ and B^⋆⋆ are defined as $\begin{array}{ccc}{\mathit{B}}^{\mathit{\star }}\mathrm{=}\mathit{B}{\left(\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{{\mathit{c}}^{\mathrm{2}}}\frac{{{\mathit{E}}_{\mathrm{\perp }}}^{\mathrm{2}}}{{\mathit{B}}^{\mathrm{2}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathit{,} {\mathit{B}}^{\mathit{\star }\mathit{\star }}\mathrm{=}\mathit{B}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{{\mathit{c}}^{\mathrm{2}}}\frac{{{\mathit{E}}_{\mathrm{\perp }}}^{\mathrm{2}}}{{\mathit{B}}^{\mathrm{2}}}\right)\mathrm{·}& & \end{array}$(NB. B^⋆ and B^⋆⋆ retain the dimensions of B).

We evolve each of Eqs. (11a) in time using a 4th order Runge-Kutta scheme with a variable timestep, subject to the numerically derived values of electric and magnetic fields found in a single snapshot of the MHD simulations of Bourdin et al. (2013; discussed in Sect. 2.1). A similar approach has been used by, for example Gordovskyy & Browning (2011), Gordovskyy et al. (2014). Each orbit is also terminated after the normalising timescale is reached, i.e. after 10 s of real time. The MHD simulation timescale is determined by the cadence of the observed magnetogram timeseries used to drive the photospheric boundary, which (for the results reported in Bourdin et al. 2013) is 90 s. A comparison of the time between MHD snapshots and the timescale used in the guiding centre approximation code implies that these scales are well separated. We also assume that the spatial scales of the gyro-motion and the global MHD environment are similarly separated, and will test this assumption for all orbits (see Sect. 3).

Fig. 2 Illustration of initial particle grid and electromagnetic environment in our study. a) shows the grid of initial particle positions (shown as black orbs) located in the coronal region of MHD simulations of a stable AR (Bourdin et al. 2013). Both a) and b) show interpolated magnetic field lines (grey lines) and observed magnetogram at the base of the simulations, and isosurfaces of parallel electric field E||; a) shows opaque pink and purple surfaces of positive and negative E|| respectively, encompassing regions of current at 25% of the peak value of E|| in the corona, while b) also includes transparent surfaces of E|| at 1% of the coronal peak value (~3.5 V m-1). For reference, three green orbs labelled “A”, “B” and “C” indicate the initial positions of orbits studied in detail in Sect. 4.

Fig. 3 Survey of the behaviour of 2000 a) electron and b) proton orbits, for initial position grid seen in Fig. 2a. Final particle positions are colour-coded by peak energy achieved during the orbit calculation lasting 10 s (or upon leaving the MHD simulation domain) with interpolated magnetic field lines (thin grey lines) included for reference.

Fig. 4 Surveys of electron and proton orbits, for initial position grid seen in Fig. 2a, subject to electric field reductions of a factor 10 (a)−b)) and 100 (c)−d)) from the original electric field (Eorig) determined by the MHD active region simulations. Final particle positions are colour-coded by peak energy (using the same colour bar as in Fig. 3 to enable direct comparisons) achieved during the orbit calculation lasting 10 s for electrons or 1000s for protons (or upon leaving the MHD simulation domain) with interpolated magnetic field lines (thin grey lines) included for reference.

3. Characteristics of global behaviour

We begin by studying the particle orbit behaviour for a uniform grid of initial positions spread throughout the computational domain in the horizontal direction at several different heights. Our aim is to obtain an overview of the global orbit behaviour one might typically expect to find throughout the AR.

For this investigation, an identical number of particles for each species were studied; 2000 electron orbits and 2000 proton orbits were given an initial pitch angle of 45° and an initial kinetic energy of 20 eV. These choices are motivated by the findings of previous work (Threlfall et al. 2015), where it was shown that (if present) the contributions from any local parallel electric fields typically outweigh the choice of initial kinetic energy and pitch angle in impacting particle behaviour. Furthermore, while the initial kinetic energy chosen here is relatively low compared to the equivalent thermal energy (20 eV ≈ 0.23 MK), this simply serves to aid illustration of particle orbit behaviour. We have also checked that a tenfold increase in initial particle energy does not significantly alter the overall findings reported here. Our particles are distributed in a 20 × 20 × 5 grid in (x, y, z); the grid is evenly spaced and extends from x,y ∈ [ 10,200 ] Mm in five z −planes, at z = [ 60,65,70,75,80 ] Mm above the photosphere. A visual representation of this grid is shown in Fig. 2a. This distribution of initial conditions of particle orbits was chosen to establish particle behaviour in a region encompassing the main AR features (including examples of both open and closed field structures), at large enough heights that both chromospheric effects (such as collisionality and partial ionisation) and numerical effects (primarily the inclusion of a stretched grid in z close to the photospheric boundary) could be minimised. Figure 3 shows the results of both the electron and proton orbit calculations for this distribution of initial conditions.

We will now highlight some key findings illustrated by Fig. 3. Firstly, all orbits (both protons and electrons) are accelerated to some degree, almost always because of the local parallel electric field. The majority of orbits achieve at least keV energies; the vast majority are accelerated to energies of 50−500 keV. By comparing Fig. 3a with Fig. 3b, it is clear that a similar number of electron and proton orbits achieve ≳ MeV energies. All of the orbits which achieve the highest energy levels in the survey have trajectories which head towards the lower solar atmosphere. Additionally, more electron than proton orbits may be found at the lowest energy levels. It is also apparent that several orbits, both proton and electron, appear to be continually reflected between mirror points; we will label these as “trapped” orbits.

To estimate the number of trapped orbits, we count any orbit which contains two or more mirror points (locations where the parallel velocity reverses sign) as being “trapped”. We estimate that, out of the 2000 orbits initialised in the manner described earlier, 1061 electron and 228 proton orbits are mirrored twice or more. To present a fair comparison (with protons taking longer to reach two or more mirror points), we extended the simulation time of proton orbits from 10 s to 1000 s. While this is much longer than the 90 s cadence of the MHD simulations, our objective here is simply to assess under what conditions similar number of protons and electrons may be trapped. With an extended time-range, the number of trapped proton orbits rises to 840. Thus, this trapping occurs for a significant fraction of the orbits considered in these surveys, but the true number of particles per species trapped by this mechanism is related to the time over which the MHD fields evolve, how these fields evolve, and the species mass and charge. Individual particle orbit examples, including trapped orbits, are studied in greater detail in Sect. 4.

The maximum kinetic energy of any orbit was achieved by an electron, whose kinetic energy reached 41.5 MeV (the highest energy reached by a proton was 28.7 MeV). The minimum kinetic energy was achieved by an electron, which gained 0.09 eV (the lowest energy gained by a proton was 0.5 eV). Regarding the question of scale separation, the maximum electron gyro-radius recorded was 0.41 m, while the maximum proton gyro-radius achieved during any orbit in the survey was 16.2 m. The x/y grid-point resolution of the MHD simulation is always 230 km; the resolution in the z-direction varies between 100 and 800 km (caused by a non-equidistant grid). By comparison with the largest gyro-radii recorded, we can state that, at least in the case of the initial survey, the length scales of the MHD simulation and the particle orbits are sufficiently separate that our use of the guiding centre approximation is clearly justified.

The amount of particle acceleration recovered in these orbit calculations is extreme, both in the number of particles which achieve modest/high energy gains, and the sizes of the gains themselves. The acceleration is a direct result of the electric field configuration used in the MHD active region model; the magnitudes of the electric fields are a consequence of the parameters used which were chosen such that the emission derived from the model matches the observed emission. This electric field is super-Dreicer over 40% of the volume, particularly low down in the atmosphere. Super-Dreicer electric fields are known to lead to runaway particle acceleration, where collisional effects can no longer compensate for increasing particle velocities (see e.g. Holman 1985). In these circumstances, it is perhaps unsurprising that we observe that many orbits are strongly accelerated, to highly non-thermal energies. What is surprising, though, is that not all particles which encounter strong electric fields immediately reach one of the boundaries. Many particles orbits which achieve moderate to high energies are retained at various heights in the atmosphere (resulting from particle trapping, as mentioned earlier).

To investigate the impact of the original MHD electric field upon particle behaviour, and in light of the amount of acceleration illustrated in Fig. 3, we repeat the orbit calculations, subject to the same initial conditions outlined earlier, but for electric fields which have been reduced by a certain factor. Figure 4 illustrates how the peak orbit energies and final positions change when the electric field strength of the MHD simulation (which we will henceforth label as E_orig) is reduced by a factor of 10 (for electron orbits in Fig. 4a, and proton orbits in Fig. 4b), and 100 (electrons/protons in Figs. 4c–d). These figures show that each decrease in electric field strength brings a corresponding reduction in the peak energy gains achieved by the orbits, while having minimal impact on the final particle positions.

The maximum energy gained in each experiment falls linearly with the electric field strength; peak energies in the experiment with a tenth of the original electric field strength (i.e. E = 0.1 E_orig) were found to be 4.1 MeV/2.87 MeV, exactly one tenth of the peak energies recorded in the original orbit surveys for electrons/protons respectively. This pattern continues for a hundred-fold electric field strength reduction (E = 0.01E_orig), where the peak energies fall to 0.41 MeV/0.287 MeV, again exactly one hundredth of the original peak energy values. Finally, we have once again counted the number of orbits which exhibit two or more mirror points. The trapped particle populations remain relatively consistent; out of 2000 orbits initialised in these reduced field strength cases, 821 electron orbits remaining trapped when E = 0.1 E_orig, and 843 when E = 0.01 E_orig. For protons (which were followed for 1000 s), we recover 759 and 850 orbits in each respective field strength case which remain trapped within this electro-magnetic environment.

To investigate the reasons behind the range of energies recovered and particle behaviours seen in this single AR system, we now focus on several examples of individual particle orbit behaviour.

4. Examples of orbit behaviour encountered

A fairly broad range of behaviour was recovered in our initial survey of orbit behaviour, detailed in Sect. 3. We will now study specific examples of particle orbits whose initial positions are highlighted in Fig. 2; these examples illustrate different types of characteristic particle behaviour observed in the simulations, which are discussed in detail below. It should be noted that, while each example highlights a single specific type of generic behaviour, the other particle behaviour may be present in each example, to a greater or lesser extent. For instance, directly accelerated particles also exhibit guiding centre drift and vice versa.

Fig. 5 Position A: example of electron which becomes ”trapped”. a) displays a zoomed view of the electron orbit (thick black line), which begins at the location shown by a green orb, and ends at the location shown by the blue orb, while black trapezoids indicate locations where the parallel velocity of the particle changes sign. The thin grey line indicates the locally interpolated B-field line based on the initial position (but almost exactly matches the particle orbit path). The red isosurface denotes where the parallel electric field changes sign (E|| = 0). Orbit properties for this electron are displayed in b) and c); b) shows the change in kinetic energy (KE) and parallel electric field component (E||) as a function of time, while c) shows the normalised parallel velocity (v||/ max(v||)) and magnetic field strength (| B |) for this orbit.

4.1. Particle A: trapped orbits

We begin with the particle labelled “A” in Fig. 2b. This is one of many orbits which are accelerated to high (but not extreme) energies over the course of the simulation. This orbit also remains within the coronal region of the AR core throughout the orbit lifetime, as with many of the particles seen in Fig. 3. For reference, the specific electron orbit we discuss is initially located at the centre of the grid, at position (x,y,z) = (150,150,75) Mm, and remains trapped in the immediate vicinity of this position for the duration of the calculation (10 s), achieving a maximum kinetic energy of 14 keV in that time.

In Fig. 5, we present a close-up view of this specific orbit during this time, together with a study of the local conditions (determined by the MHD simulation) encountered by the particle over the course of its orbit. From Fig. 5a, it is clear that the particle is trapped between mirror points (locations where the local parallel particle velocity reverses sign). One common cause of particle trapping is the magnetic mirror effect (where regions of increasing magnetic field strength cause the particle to reverse direction because of the adiabatic invariance of the magnetic moment, μ). To establish the role played by the magnetic mirror effect in this particular case, we study both the local parallel velocity and encountered magnetic field strength, in Fig. 5c, which suggests that changes in sign of v_|| are not strongly associated with locations of peaks of | B |. To further clarify the role of the magnetic mirror effect, we calculate the loss cone angle for this position. Particles with pitch angle θ are naturally trapped by the magnetic mirror effect (in a so-called “magnetic bottle”) if $\begin{array}{ccc}\sin \mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{>}\sqrt{\frac{\mathit{B}}{{\mathit{B}}_{\mathit{M}}}}\mathit{,}& & \end{array}$where B_M is the maximum magnetic field strength achieved at any location along a given field line; particles which escape this magnetic trap (i.e. those with pitch angles which do not satisfy the above inequality) are said to be in the “loss cone”. The field line orbited by this electron (represented in Fig. 5b by a thin grey line) is actually a closed field line, anchored in the photosphere at both ends. We estimate the peak magnetic field strength at one foot-point to be 137 G, and 558 G at the other; these are the locations of greatest magnetic field strength at any point along this single field line (which are of course of opposite sign in the vertical direction at the photosphere). Thus, to be in the loss cone for this particular field line, a particle must have a pitch angle θ< 8.68° in one direction and θ< 4.30° in the other. We repeated the orbit calculation for an electron with the same initial position, but with a reduced initial pitch angle of θ = 4.0°; this now would ordinarily be sufficient for the particle to escape the natural magnetic trap caused by increasing magnetic field strength near the field-line foot-points. The orbit yielded by this second calculation was almost identical to that with θ = 45°; the particle remains trapped regardless of initial pitch angle. We can therefore exclude magnetic mirroring as the cause of this type of electron trapping. This naturally leads to a question; if the electrons are not magnetically trapped, then how are they trapped?

The answer lies in the variation of the local electric field. For this electron, the orbit path encounters a change in the sign of the local parallel electric field, E_||, as shown by the red curve in Fig. 5b. An iso-surface of the contour for which E_|| = 0 can be seen in Fig. 5a; the electron repeatedly crosses this contour over the course of its orbit.

Upon crossing the contour of E_|| = 0, the electron begins to decelerate. At the time when the local parallel velocity is approximately zero, the parallel electric field strength peaks; the orbit reverses direction (parallel to the magnetic field), re-crosses the E_|| = 0 contour and re-enters the region of oppositely signed E_||, leading to further (and repeated) mirroring. This pattern repeats for the entire duration of the orbit simulation. This is not an isolated example; the vast majority of electrons which do not encounter the simulation boundaries are typically “trapped” this way. This trap is particularly effective for electrons. Even relatively weak regions of oppositely signed parallel electric field strength are sufficient to repeatedly accelerate electrons this way.

By comparison, protons are much heavier; for a proton to be trapped in this manner requires both extreme and prolonged changes in the local value of E_|| in order for the proton to be decelerated sufficiently to reverse direction. Figure 3b shows that a significant number (840) of protons are also trapped, and hence are retained close to the centre of the AR core. This is unsurprising, from the size and extent of the E_|| regions near the footpoints of field-lines which thread the AR core shown in Fig. 2b. An example of a proton which is trapped in this manner will be studied in more detail in Sect. 4.3.

This trapping mechanism is clearly highly dependent on the (orientation of the) local electric field. However, the electric fields used in the simulation might be considered to be too strong for an active region which does not produce a flare, as they lead to highly accelerated particles. We therefore investigate how this trap is affected if the simulated electric field strength is reduced. In Fig. 6, we illustrate the parallel electric field strength and parallel velocity of three electron orbits for electrons that all have the same initial conditions and start from the same initial position, (110,110,80) Mm, but are affected by electric fields whose strengths have been scaled to different values. Beginning with the original electric field strength, E = E_orig, we see from Fig. 6b that the orbit repeatedly mirrors, at the instants when the electric field strength peaks in Fig. 6a. Reducing the electric field strength by a factor of 10 (E = 0.1 E_orig, shown in purple) shows that the orbit still repeatedly mirrors, but takes longer to reach each mirror point compared to the original field strength case. Further reducing the electric field strength so E = 0.01E_orig (shown in orange) demonstrates that this trap is robust; once again the orbit repeatedly mirrors until the orbit calculation is terminated, even when the orbit calculations are allowed to continue for longer timescales.

Fig. 6 Impact of decreasing electric field strength upon electric field trapping mechanism. An electron orbit (which initially is trapped between regions of strong, oppositely oriented electric field) continues to repeatedly mirror when the electric field strength is reduced by a factor of 10 (purple) and 100 (orange). a) shows the parallel component of electric field, E||, experienced in each field strength case (see legend) while b) shows the corresponding changes in parallel velocity, v||.

4.2. Particle B: proton drifts

Turning our attention to protons, we now focus on another type of behaviour seen in the initial survey. Several proton orbits are also seen to gain only a relatively small amount of energy during the 10 s calculation. One such orbit (labelled “B” in Fig. 2b) was initialised at (x,y,z) = (125,125,80) Mm and gained (at its peak) 361 eV in kinetic energy. Once again we present a close-up view of the orbit path, together with local plasma conditions in Fig. 7.

It is clear from Fig. 7a that this proton experiences a relatively constant (anti-)parallel electric field strength throughout its lifetime; the orbit path remains distant from the contour representing E_|| = 0 (also shown by the value of E_|| in Fig. 7b). However, despite the lack of a change in sign of E_|| or a clear peak in | B | (see Fig. 7a) which would potentially signal the magnetic mirror effect, a mirror point is clearly visible close to the point of initialisation (green orb), and in-between the initial and final particle positions (cf. Fig. 5a, where the mirror points enclose the initial and final particle positions!). The mirror point location may be understood by again considering local plasma conditions.

Having initialised the proton orbit with a positive parallel velocity (i.e. a velocity parallel to B), the proton decelerates from its slow initialisation speed (caused by the sign of the local parallel electric field). This deceleration takes place over a long timescale compared to the equivalent deceleration of an electron, resulting from the difference in species mass. During the phase where parallel orbit velocities are small, particle drifts (primarily the E × B drift) begin to play a significant role in the proton trajectory.

Figure 7a includes interpolated magnetic field lines based on the changing position of the proton guiding centre over time. From this, it is clear that the proton no longer follows the same field line for all time, but instead drifts across many field lines¹, before finally beginning to accelerate in the opposite direction to the initial direction of travel. For this particular orbit, after 10 s, the perpendicular displacement of the proton from its original field line is approximately 77 km; by comparison, the final parallel displacement is 1.03 Mm (with a peak positive displacement along the field of approximately 31 km, before being accelerated in the opposite direction over a distance of 1.06 Mm). Over the course of the orbit, the proton experiences a near-constant E × B drift velocity of approximately 7.81 km s^-1; integrating this value over the 10 s orbit lifetime yields a perpendicular displacement which agrees with that observed in Fig. 7a. Furthermore, the perpendicular displacement is well-aligned with the local E × B vector. While the parallel displacement is likely to further increase with time and typically will dominate over perpendicular drifts (particularly for large values of E_|| or long timescales), this case demonstrates that proton orbits may in fact drift over relatively vast distances, particularly upon encountering regions of deceleration.

Fig. 7 Position B: example of proton which exhibits a slow guiding centre drift. a) displays a zoomed view of the proton orbit (black), which begins at the location shown by a green orb, and ends at the location shown by a blue orb, with black trapezoids indicating the location of a change in sign of parallel velocity. The thin grey lines indicate locally interpolated B-field lines based on several orbit positions. The red isosurface denotes where the parallel electric field changes sign (E|| = 0). Orbit properties for this proton are displayed in a) and c); b) shows the change in kinetic energy (KE) and parallel electric field component (E||) as a function of time, while c) shows the normalised parallel velocity (v||/ max(v||)) and magnetic field strength (| B |) for this orbit.

Fig. 8 Position C: example of strong proton and electron acceleration. a) displays a zoomed view of the proton (blue) and electron (black) orbits, which both begin at the location shown by a green orb, ending at a location(s) shown by a blue orb, and with changes in sign of v|| indicated by black trapezoids (noting that, on this scale, several symbols occur close together; these symbols are shown in the dark green box). The thin grey lines indicate locally interpolated B-field lines based on the initial position (which, for the most part, closely match orbit paths), while contours of strong positive/negative E||, and the photospheric magnetogram of vertical magnetic field strength are included for context (see Fig. 2). Orbit properties for the electron and proton are displayed in b) and c) as diamonds (electron) and plus-signs (proton); b) shows the change in kinetic energy (KE) and parallel electric field component (E||) as a function of time, while c) shows the normalised parallel velocity (v||/ max(v||)) and magnetic field strength (| B |) for each orbit.

4.3. Particle C: Accelerated electrons/trapped protons

Many examples of both the electron “trapping” and proton drifting can be observed in our survey. However, whether or not an orbit becomes trapped is determined by the local electromagnetic configuration, the timescale over which the orbit is simulated, and the particle species. Some orbits in the initial survey are strongly accelerated in only a single direction, until they leave the simulation domain. To understand this category of behaviour, we illustrate the orbits of both an electron and a proton which are both initialised in the same location (labelled “C” in Fig. 2b), and which (initially) encounter (the same) strong parallel electric fields. As with previous examples, we will illustrate a close-up view of the orbit path and the local plasma parameters encountered, in Fig. 8. For reference, the initial position used in this example is (x,y,z) = (100,200,80) Mm. To further emphasise the differences in behaviour between species, we extend the simulation time from 10 s to 100 s for this case. This extends the particle simulations beyond the cadence of the MHD simulations; in extending the simulation time in this manner, we simply aim to illustrate particle behaviour over a longer term.

In Fig. 8a, the initial position is located near the apex of a closed field line which links the foot-points of the AR studied. The proton (blue orbit) and electron (black orbit) are accelerated towards opposing foot-points (as one might expect from the charge difference between the species). Furthermore, while the electron is rapidly accelerated almost all the way to the photospheric boundary, the proton never leaves the (upper) corona, but is ultimately reflected and returns along the initial trajectory. In Fig. 8c, we see that the electron (diamonds) reaches the photosphere after approximately 0.5 s. In this time, the proton (plus signs) has barely moved. The proton achieves maximum parallel speed after approximately 4 s and mirrors approximately seven times in 100 s, with the first mirroring taking place after approximately 10 s.

To assess the reason for the proton becoming trapped in this way, we once again calculate the loss cone (see discussion in Sect. 4.1 for more details). Along the field line in question, we estimate that initial pitch angles less than 5.81° and 6.61° would able particles to escape the trap caused by increasing magnetic field strength. As with the case shown in Sect. 4.1, repeating this experiment with an initial pitch angle of 4.0° shows negligible differences in proton/electron behaviour; both the acceleration and trapping behaviour displayed in Fig. 8 are caused by the electric field. In extending the simulation time from 10 s to 100 s, we demonstrate that the electric field continues to trap the proton in the upper atmosphere for all time, while accelerating the electron until it reaches the loop foot-point at the photospheric boundary.

All particle orbit calculations are automatically terminated upon reaching a 5 Mm layer above the lower simulation boundary (which results from the limited applicability of our test-particle approach in this region). At termination, the electron has travelled approximately 118 Mm in less than 0.5 s. By comparison, the distance between mirror points for the proton is approximately 38 Mm. The reason for the differences between the proton and electron results are twofold; from Fig. 8c, we can see that the proton achieves its peak velocity after approximately 4 s, taking much longer to accelerate. Upon encountering a region of oppositely oriented E_||, the proton is unable to overcome the resulting deceleration, and mirrors. By comparison, the electron is continually accelerated as it plummets towards the photosphere. The sign of the electric field strength is positive throughout the lifetime of the electron orbit (see Fig. 8b), hence the electric field trapping mechanism plays no role in the electron orbit. Furthermore, the magnetic mirror effect also plays no role in this orbit, as the strong, uni-directional electric field encountered by this electron results in a high parallel velocity (1875 km s^-1 or 0.006c).

5. Discussion

We have conducted a survey of test-particle behaviour in an environment which aims to describe the behaviour of a solar AR. From this survey (detailed in Sect. 3), it is clear that this environment (at the electric field strengths derived by the MHD model) would cause a significant amount of particle acceleration to non-thermal energies. To avoid confusion, we must first emphasise that our results are based on particle orbits and not particle fluxes. In our surveys, we have used a relatively coarse grid of initial positions to describe particle behaviour at specific locations, for a single value of initial kinetic energy and pitch angle. Such choices mean that a direct comparison with observed particle spectra from this region is not possible; we can only broadly compare common types of particle orbit behaviour.

Even with the constraint that we cannot directly compare our findings with observed spectra, it is clear that typical orbit energy gains found in our work are much greater than would be expected from a non-flaring active region. MHD models (like the one considered here) have been shown to accurately reproduce specific observational aspects of the solar feature they have been designed to model. Our orbit calculation results show that kinetic processes also have implications for the parameter regimes chosen in such models; in our case, our orbit calculations imply that the extent of regions of strong electric field created within the MHD model would lead to significant amounts of runaway particle acceleration, which (crucially) was not observed for this specific active region. Thus particle orbit calculations could and should be used in future to inform in the design of MHD models of specific structures, with the goal of being able to closely represent the physical processes taking place within those structures, and accurately reflect observations of both heating and particle behaviour.

The effects of the strength and extent of the electric fields determined in this model are clearly illustrated in the orbit results seen in Fig. 2. We estimate the peak parallel electric field in the simulation domain to be ${\mathit{E}}_{\mathrm{\left|}\mathrm{\right|}}^{\mathrm{peak}}\mathrm{=}\mathrm{295}$ V m^-1, while in the corona (i.e. above the lowest 10 gridpoints), ${\mathit{E}}_{\mathrm{\left|}\mathrm{\right|}}^{\mathrm{peak}}\mathrm{\approx }\mathrm{5}$ V m^-1. For comparison, we estimate the Dreicer electric field to approach 6 × 10^-2 V m^-1 (for n = 1 × 10¹⁶ m^-3, T = 1 MK). Regions of super-Dreicer electric field cover extensive areas of the domain (resulting in the orbit energies seen in Fig. 2). In Threlfall et al. (2015), a magnetic reconnection experiment with ${\mathit{E}}_{\mathrm{\left|}\mathrm{\right|}}^{\mathrm{peak}}\mathrm{\approx }\mathrm{0.1}$ V m^-1 was sufficient to rapidly accelerate many particles to non-thermal energies, despite the apparent strength appearing too small to be considered relevant for solar flares. In fact, both the electric field strength and the extent of the reconnection region (which, for some cases in Threlfall et al. 2015, reached 20 Mm) determine the reconnection rate and (ultimately) the energy gained by the particle through direct acceleration. Thus both the strength and extent of the electric fields in the work of Bourdin et al. (2013) are responsible for the non-thermal particle populations recovered here.

We have also shown that reductions of the electric field strength from the MHD simulations (by factors of both ten and one hundred) reduce the peak energies gained by the particles by the same amount. Once again, with the energy gains directly resulting from the presence of parallel electric field, E_||, reductions in the strength of E reduce the voltage through which each particle may be accelerated, lowering the ultimate amount of energy they gain.

While many of the particles experience direct acceleration, not all are accelerated out of the domain. A significant number remain trapped in the upper regions of the model atmosphere instead. In Sect. 4.1, we illustrate an electron which is “trapped” between weak regions of oppositely oriented parallel electric field, continuously re-accelerating to a kinetic energy of approximately 10 keV. In Sect. 4.3, we show an example of a proton which is similarly trapped; notably protons take much longer to be repeatedly re-accelerated, but can still achieve high kinetic energy gains (100 keV in the case of the proton in Sect. 4.3) in between mirror points. The strength and extent of the two regions of oppositely oriented electric field are the primary factor that controls the peak energy which may be repeatedly gained and lost by particles of either species in an electric field trap.

The presence of this type of trap has been noted in studies of magnetic reconnection at stochastic current sheets (Turkmani et al. 2006) and 3D magnetic null-points (Stanier 2013), but neither paper investigates the cause of the trapping. In our work, we have studied this trap in detail, and find that the number of particles trapped by this mechanism may remain relatively unchanged despite reductions in electric field strength in the MHD model.

For a general coronal structure, one might expect that particles are likely to encounter regions where the parallel electric fields change sign (for example between different coronal loop footpoints, or at fragmented current sheets). Our results show that this potentially may occur at many locations throughout a specific active region, and trap energised particles at heights above the photosphere, particularly in the corona. A likely way for particles to escape such a trap would be from changes in the global electro-magnetic environment; such changes would take place over much longer timescales than those considered here.

Whether or not this type of trapping would yield an observational signature (e.g. in radio wavelengths) is at present unknown, and may be worthy of further investigation. While trapped, these particles would also likely lose energy via collisional effects, particularly lower down in the atmosphere, or at low speeds (e.g. during deceleration/reflection).

Other, more familiar types of particle behaviour are also present in our orbit surveys. Particularly in cases which involve proton deceleration (in Sect. 4.2), particle drifts may play a significant role in particle trajectories, enabling a proton to drift upto (potentially) tens/hundreds of kilometres from its original field line. Such drifts may contribute to the lack of alignment in X-ray and γ-ray features observed during a single solar flare event (e.g. Hurford et al. 2003).

Finally, it is worth noting that, of the characteristic types of behaviour discussed in Sect. 4, some types are more affected by the choice of numerical resistivity used in the simulation than others. Reducing the resistivity would weaken the strength of the electric fields generated, which (in turn) reduces the magnitude of the E × B drift velocity, and the voltage differences through which particles are accelerated. However, we have shown in Sect. 3 that the number of trapped orbits remains relatively constant even when decreasing electric field strengths by factors of up to 100.

Further investigation (Sect. 4.1) reveals that electric field strength reductions simply allow orbits to travel further prior to mirroring. This results from the orientation of the electric field, which is unaffected by the scaling so that in weaker electric fields orbits have further to travel before encountering a sufficiently strong electric field to be reflected, causing the majority of particle orbits to remain trapped. Indeed, with broader regions of no parallel electric field, more particles find themselves mirroring between smaller weaker regions of oppositely oriented parallel electric field. Only by encountering regions of strong E_|| (and thus achieving high velocities) or by-passing the electric field regions all-together (resulting from long-term changes in the MHD environment in which the orbits are performed) are such orbits capable of avoiding the trap set by these weak electric field regions.

6. Conclusions and future work

We have studied the relativistic guiding centre behaviour of test particles in an observationally driven large-scale 3D magnetic reconnection experiment (Bourdin et al. 2013).

We find that the majority of particle motion is controlled by the strength, orientation and extent of the local parallel electric field; in a significant number of cases, the electric field causes strong direct acceleration of particles to non-thermal energies. On the one hand, the strength of the electric field (which reaches super-Dreicer levels throughout a significant fraction of the domain) means that collisions would play no role in the particle behaviour. On the other hand, the lack of evidence for a flare-like event within this AR combined with our results imply that the resistivity and current structures created in the MHD simulation generate too many orbits which are ultimately accelerated to unacceptably high kinetic energies. Reducing the electric field strength of the simulations leads to less strongly accelerated particles. However, even weak electric fields stretched over long distances have the capability to accelerate particles to very high energies.

As part of this investigation, we have identified several behavioural characteristics of particles within this environment. Direct acceleration and drifting of particles were anticipated and are weakened by simulations containing weaker/shorter electric field regions. However, we also observe a type of trapping of particles between regions of oppositely oriented parallel electric field which has previously only been seen in isolated, idealised studies of magnetic reconnection (Turkmani et al. 2006; Stanier 2013). This type of trap results from a change in sign of the parallel electric fields, and is likely to be present in any test particle calculations based on compressible resistive MHD simulations of AR behaviour. Reducing the electric field strengths in the model does not appear to affect the efficiency of this trap. Whether or not this type of trap would produce an observational signature is unknown, but if it does it could allow for further insights into electric field structures within the solar atmosphere.

Following on from the present investigation, several opportunities for further work present themselves. The work of Bourdin et al. (2013) merits repeating, particularly to see how the results presented here change for lower magnetic diffusivities. The use of observations to drive MHD simulations is important, and worthy of repeating for a more recent active-region system (including regions that produce or do not produce a flare) which can be observed by a full suite of modern instruments, e.g. RHESSI, AIA/SDO, STEREO and ground-based telescopes. This might also allow for a comparison of synthetic hard X-ray data (generated by test-particle simulations) and a specific flare-like event, in a similar manner to that used in Gordovskyy et al. (2014). Furthermore, development of the test-particle model (for example, by studying the effect of collisions, or through the use of a distribution function which better reflects the initial state of particles in a coronal environment) would also be beneficial.

Acknowledgments

The authors gratefully acknowledge the support of the UK Science and Technology Facilities Council [Consolidated Grant ST/K000950/1]. The research leading to these results has received funding from the European Commission’s Seventh Framework Programme FP7 under the grant agreement SHOCK (project number 284515). This work was supported by the International Max-Planck Research School (IMPRS) on Solar System Physics. The results of this research have been achieved using the PRACE Research Infrastructure resource Curie based in France at TGCC, as well as JuRoPA hosted by the Jülich Supercomputing Centre in Germany. Preparatory work has been executed at the Kiepenheuer-Institut für Sonnenphysik in Freiburg, as well as on the bwGRiD facility located at the Universität Freiburg, Germany. We thank Suguru Kamio for his help finding active region observations. Hinode is a Japanese mission developed, launched, and operated by ISAS/JAXA, in partnership with NAOJ, NASA, and STFC (UK). Additional operational support is provided by ESA and NSC (Norway).

References

All Figures

Fig. 1 Panels, from L-R, showing a) a Hinode/SOT magnetogram of an AR core seen on the 14th Nov. 2007 saturated at ± 300 G; b) a Hinode/EIS image of the Fe xv emission (at 1.5 MK) of the corona above this region and c) the co-spatial and co-temporal simulated parallel electric field strength and orientation within this region (where pink/purple represents positive/negative E||, and the opacity represents the relative strength of E||− fully opaque regions identify ${\mathit{E}}_{\mathrm{\left|}\mathrm{\right|}}\mathrm{=}\mathrm{0.25}{\mathit{E}}_{\mathrm{\left|}\mathrm{\right|}}^{\mathrm{peak}}$ while transparent regions show ${\mathit{E}}_{\mathrm{\left|}\mathrm{\right|}}\mathrm{=}\mathrm{0.01}{\mathit{E}}_{\mathrm{\left|}\mathrm{\right|}}^{\mathrm{peak}}$, for a peak coronal electric field strength ${\mathit{E}}_{\mathrm{\left|}\mathrm{\right|}}^{\mathrm{peak}}\mathrm{~}\mathrm{3.5}\mathrm{V}\hspace{0.17em}{\mathrm{m}}^{-1}$ in the simulation domain). The right hand panel is generated by simulations described in Bourdin et al. (2013), which also describes field lines which cross the intensity maxima of several EUV-emissive coronal loops, which have been overlaid in each panel; “CL” indicates core loops, while “SL” indicates shorter loops. The circles indicate EUV-emissive loop footpoints in the model and the corresponding plasma motion (and observed Doppler shifts) at these locations; red represents draining (i.e. motion along the line-of-sight away from an observer) while blue represents upflows (motion towards an observer).

In the text

Fig. 2 Illustration of initial particle grid and electromagnetic environment in our study. a) shows the grid of initial particle positions (shown as black orbs) located in the coronal region of MHD simulations of a stable AR (Bourdin et al. 2013). Both a) and b) show interpolated magnetic field lines (grey lines) and observed magnetogram at the base of the simulations, and isosurfaces of parallel electric field E||; a) shows opaque pink and purple surfaces of positive and negative E|| respectively, encompassing regions of current at 25% of the peak value of E|| in the corona, while b) also includes transparent surfaces of E|| at 1% of the coronal peak value (~3.5 V m-1). For reference, three green orbs labelled “A”, “B” and “C” indicate the initial positions of orbits studied in detail in Sect. 4.

In the text

	Fig. 3 Survey of the behaviour of 2000 a) electron and b) proton orbits, for initial position grid seen in Fig. 2a. Final particle positions are colour-coded by peak energy achieved during the orbit calculation lasting 10 s (or upon leaving the MHD simulation domain) with interpolated magnetic field lines (thin grey lines) included for reference.
In the text

Fig. 4 Surveys of electron and proton orbits, for initial position grid seen in Fig. 2a, subject to electric field reductions of a factor 10 (a)−b)) and 100 (c)−d)) from the original electric field (Eorig) determined by the MHD active region simulations. Final particle positions are colour-coded by peak energy (using the same colour bar as in Fig. 3 to enable direct comparisons) achieved during the orbit calculation lasting 10 s for electrons or 1000s for protons (or upon leaving the MHD simulation domain) with interpolated magnetic field lines (thin grey lines) included for reference.

In the text

Fig. 5 Position A: example of electron which becomes ”trapped”. a) displays a zoomed view of the electron orbit (thick black line), which begins at the location shown by a green orb, and ends at the location shown by the blue orb, while black trapezoids indicate locations where the parallel velocity of the particle changes sign. The thin grey line indicates the locally interpolated B-field line based on the initial position (but almost exactly matches the particle orbit path). The red isosurface denotes where the parallel electric field changes sign (E|| = 0). Orbit properties for this electron are displayed in b) and c); b) shows the change in kinetic energy (KE) and parallel electric field component (E||) as a function of time, while c) shows the normalised parallel velocity (v||/ max(v||)) and magnetic field strength (| B |) for this orbit.

In the text

Fig. 6 Impact of decreasing electric field strength upon electric field trapping mechanism. An electron orbit (which initially is trapped between regions of strong, oppositely oriented electric field) continues to repeatedly mirror when the electric field strength is reduced by a factor of 10 (purple) and 100 (orange). a) shows the parallel component of electric field, E||, experienced in each field strength case (see legend) while b) shows the corresponding changes in parallel velocity, v||.

In the text

Fig. 7 Position B: example of proton which exhibits a slow guiding centre drift. a) displays a zoomed view of the proton orbit (black), which begins at the location shown by a green orb, and ends at the location shown by a blue orb, with black trapezoids indicating the location of a change in sign of parallel velocity. The thin grey lines indicate locally interpolated B-field lines based on several orbit positions. The red isosurface denotes where the parallel electric field changes sign (E|| = 0). Orbit properties for this proton are displayed in a) and c); b) shows the change in kinetic energy (KE) and parallel electric field component (E||) as a function of time, while c) shows the normalised parallel velocity (v||/ max(v||)) and magnetic field strength (| B |) for this orbit.

In the text

Fig. 8 Position C: example of strong proton and electron acceleration. a) displays a zoomed view of the proton (blue) and electron (black) orbits, which both begin at the location shown by a green orb, ending at a location(s) shown by a blue orb, and with changes in sign of v|| indicated by black trapezoids (noting that, on this scale, several symbols occur close together; these symbols are shown in the dark green box). The thin grey lines indicate locally interpolated B-field lines based on the initial position (which, for the most part, closely match orbit paths), while contours of strong positive/negative E||, and the photospheric magnetogram of vertical magnetic field strength are included for context (see Fig. 2). Orbit properties for the electron and proton are displayed in b) and c) as diamonds (electron) and plus-signs (proton); b) shows the change in kinetic energy (KE) and parallel electric field component (E||) as a function of time, while c) shows the normalised parallel velocity (v||/ max(v||)) and magnetic field strength (| B |) for each orbit.

In the text