Issue 
A&A
Volume 542, June 2012



Article Number  A47  
Number of page(s)  14  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201218857  
Published online  04 June 2012 
Solar particle acceleration at reconnecting 3D null points
^{1} Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK
email: stanier@jb.man.ac.uk
^{2} Jeremiah Horrocks Institute, University of Central Lancashire, Preston PR1 2HE, UK
Received: 20 January 2012
Accepted: 14 March 2012
Context. The strong electric fields associated with magnetic reconnection in solar flares are a plausible mechanism to accelerate populations of high energy, nonthermal particles. One such reconnection scenario, in a fully 3D geometry, occurs at a magnetic null point. Here, global plasma motion can give rise to strong currents in the spine axis or fan plane.
Aims. We aim to understand the mechanism of charged particle energy gain in both the external drift region and the diffusion region associated with 3D magnetic reconnection. In doing so we aim to evaluate the efficiency of resistive spine and fan models for particle acceleration, and find possible observables for each.
Methods. We used a full orbit test particle approach to study proton trajectories within electromagnetic fields that are exact solutions to the steady and incompressible magnetohydrodynamic equations. We studied the acceleration physics of single particle trajectories and found energy spectra from many particle simulations. The scaling properties of the accelerated particles with respect to field and plasma parameters was investigated.
Results. For fan reconnection, strong nonuniform electric drift streamlines can accelerate the bulk of the test particles. The highest energy gain is for particles that enter the current sheet, where an increasing “guide field” stabilises particles against ejection. The energy is only limited by the total electric potential energy difference across the fan current sheet. The spine model has both slow external electric drift speed and weak energy gain for particles reaching the current sheet.
Conclusions. The electromagnetic fields of fan reconnection can accelerate protons to the high energies observed in solar flares, gaining up to 0.1 GeV for anomalous values of resistivity. However, the spine model, which gave a harder energy spectrum in the ideal case, is not an efficient accelerator after pressure constraints in the resistive model are included.
Key words: Sun: flares / Sun: particle emission / Sun: Xrays,γrays / magnetic reconnection / acceleration of particles / Sun: corona
© ESO, 2012
1. Introduction
Observations of Hard Xray (HXR) and γray emission from solar flares by the RHESSI space telescope (Lin et al. 2002) suggest that a large proportion of magnetic energy is converted into kinetic energy of nonthermal accelerated particles. The dominant HXR sources are chromospheric footpoints of the flaring loops, at which there is continuum freefree emission from a beam of energetic electrons in collision with ambient plasma (Brown 1971). This continuum spectrum gives beam electron energies from around 10 keV up to almost 100 MeV (Lin 2006). There is line emission at the γray end of the spectrum from processes involving accelerated ions such as neutroncapture and nuclear deexcitation (see e.g. Vilmer et al. 2011, for a review), suggesting ions with energies up to ~100 MeV/nucleon.
When the emission from the footpoints is weak, or when they are occulted by the solar limb, a weaker HXR emission source is sometimes observed above the top of the flare loops (Masuda et al. 1994). Recent observations of two such flares indicate that this emission is nonthermal and that the source is actually the acceleration site for a significant number of energetic electrons (Krucker et al. 2010; Ishikawa et al. 2011). The estimated number of energetic electrons is a significant fraction of the emission site density, setting tough efficiency constraints on any proposed acceleration mechanism.
It is well accepted that magnetic reconnection plays the key role in the dissipation of magnetic energy during a flare and there is a growing body of observational signatures for the process (McKenzie 2011). The site of reconnection in the standard (CSHKP) flare model (e.g. Priest 2000) is above the thermal loops, not in disagreement with the site of the nonthermal coronal HXR source. SuperDreicer (Dreicer 1959) electric fields associated with reconnection are one plausible mechanism for particle acceleration and much theoretical work has been done to investigate the efficiency of this mechanism (for review, see Zharkova et al. 2011).
Early work on charged particle trajectories within a reconnecting current sheet concentrated on single particle motion and energy gain in idealised field configurations. Speiser (1965) found that charged particles in the simplest current sheet, of oppositely directed magnetic field and constant electric field, are trapped so the energy gain is limited only by the sheet length. With an additional small and constant magnetic field component normal to the current sheet plane, the particles are turned by 90° and ejected from the current sheet into the external drift region. Zhu & Parks (1993), Litvinenko & Somov (1993) and Litvinenko (1996) also included a third component of the magnetic field, a guide field parallel to the electric field. Above a critical guide field the trajectory is stabilised against ejection and the energy gain is once again only bounded by the sheet length (Litvinenko 1996).
This early analytical work was extended with 2D (or 2.5D where the fields are invariant in the third dimension) test particle simulations, many of which have been carried out using simple prescribed magnetic and electric fields that would be expected in a reconnection solution to the MHD equations. These simulations consider the effect of the guide field on trajectories and energy spectra within the Harris current sheet (Zharkova & Gordovskyy 2004, 2005; Wood & Neukirch 2005) and magnetic Xpoints (Vekstein & Browning 1997; Hannah & Fletcher 2006; Hamilton et al. 2005).
Heerikhuisen et al. (2002) and Craig & Litvinenko (2002) used magnetic and electric fields from the exact analytical solutions of Craig & Henton (1995) to the 2D incompressible, resistive MHD equations. Also, an approach combining numerical MHD simulations with a test particle code has been used to study 2D forced reconnection (Gordovskyy et al. 2010a,b). These simulations can include compressibility and time evolution, making them more realistic for coronal plasma. However, analytical solutions are essential to study acceleration due to reconnection in very large Lundquist number plasma at present.
The complexity of the coronal magnetic field in a flaring Active Region motivates the study of test particle motion in fully 3D reconnection geometries. Reconnection models in 3D are comparatively new, but it is clear that there are significant qualitative differences from the familiar 2D models (see e.g. Pontin 2011). Reconnection in 3D can occur both with and without magnetic null points. However, the simplest 3D magnetic configuration is based on the potential magnetic field about such a null point. Here, the solenoidal condition defines the magnetic topology of a 1D spine line and a 2D fan surface (called γline and Σsurface by Lau & Finn 1990) that separates different magnetic flux domains (a linear description of magnetic configurations at nulls is given by Parnell et al. 1996). Although these null points cannot be measured in the corona at present, there is some indirect evidence for their existence. Nulls are common features of magnetic topology models that reconstruct the magnetic field from photospheric magnetograms into the corona (see Longcope 2005, for a review). The application of this method at several flare sites suggests the importance of these nulls in certain flares (Des Jardins et al. 2009; Aulanier et al. 2000; Fletcher et al. 2001). Reconnection at magnetic null points is also thought to be important for the emergence of new flux from beneath the photosphere into the corona (Török et al. 2009; Maclean et al. 2009; Liu et al. 2011).
The type of reconnection that occurs at a 3D null depends upon the magnetic configuration and global plasma flow. Priest & Titov (1996) proposed two models of reconnection using a potential magnetic field and prescribed global flows that satisfy the ideal MHD equations. In ideal spine reconnection a shear flow across the fan plane causes frozenin flux inflow that converges on the spine axis. At the spine the field reconnects in the presence of singular electric field. For ideal fan reconnection the singular electric field occurs in the fan plane, driven by a shearing flow across the spine axis. Craig et al. (1995), Craig & Fabling (1996) and Craig et al. (1997) found exact solutions to the steady and incompressible resistive MHD equations at 3D null points by considering a fluxpileup disturbance field superposed with the background potential magnetic field. The disturbance field induces a current sheet in the spine axis and in the fan plane for the resistive spine and resistive fan models respectively. Craig & Fabling (1998) found corresponding timedependent solutions to these steady models, and the numerical simulations of Heerikhuisen & Craig (2004) found reconnection scalings in agreement with both steady state and time dependent models at peak reconnection rate. The 3D MHD simulations of Pontin et al. (2007a,b) also found a hybrid of the spine and fan models, named spinefan reconnection, when compressibility is included. Recent numerical and analytical study gives additional models for null reconnection when the global plasma motion is rotational rather than a shear flow (see for review Priest & Pontin 2009; Pontin et al. 2011).
It is not yet known if these null points are effective particle accelerators. Previous work by Dalla & Browning (2005, 2006, 2008) and Browning et al. (2010) used the ideal electromagnetic fields of Priest & Titov (1996) in a test particle code, finding the ideal spine reconnection model was effective to accelerate protons and electrons to high energies (max ~ 10^{7} eV) for solar coronal parameters. The ideal fan reconnection model was less effective for protons, partly as the geometry of the electric drift streamlines was less efficient at delivering particles to regions of high electric field. Guo et al. (2010) used null point magnetic and electric field configurations from MHD simulations, finding that strong electric fields due to convective plasma motion can be efficient at accelerating protons but less so for electrons. Litvinenko (2006) used the WKB method of Bulanov & Cap (1988) to show that single protons and electrons close to the null in the reconnecting fan current sheet can achieve the high energies observed in flares. However, this energy is limited as the particles become unstable in the sheet, due to the potential background field, and are ejected.
In this paper we examine test particle trajectories and energy spectra of protons in electromagnetic fields which are exact solutions to the 3D, steadystate, incompressible and resistive MHD equations at magnetic null points. These are the resistive spine and resistive fan solutions given in Craig et al. (1995); Craig & Fabling (1996) and Craig et al. (1997). We consider trajectories that start both in the outer ideal region, for comparison with particle acceleration results in the ideal models (Browning et al. 2010), and those that start directly inside the resistive fan current sheet, to compare with the analytical work of Litvinenko (2006).
The paper is organised as follows. In Sect. 2 we describe the model fields used, along with parameters chosen considering the pressure constraints and optimisations of Craig et al. (1997) and Craig & Watson (2000). We also derive the electric fields and potentials used in the code, and give approximate external drift velocity scalings. In Sect. 3 we give the results of test particle simulations for thermal distributions of protons starting in the drift region for each model. We choose typical high energy particles from these simulations and follow the single trajectories to understand the energy gain mechanism. We see how drift times and energy spectra scale with the different parameter choices in the fan model. In Sect. 4 we give a summary and conclusions on the efficiency of each model for accelerating protons.
2. Model fields and the test particle code
We consider two of the reconnection solutions at 3D null points found by Craig et al. (1995) and developed in Craig & Fabling (1996) and Craig et al. (1997). The fields satisfy the resistive, steadystate, and incompressible MHD equations. These are normalised in the usual way by the characteristic magnetic field strength B_{0}, at a typical length scale L_{0} from the null point, and by density ρ_{0}. This choice leads to the natural units for velocities in terms of the external Alfvén speed, . The thermal pressure is normalised by the dynamic pressure at Alfvén speed, , so that the dimensionless pressure on the L_{0} boundary is half the plasma beta p_{e} = β_{e}/2. The dimensionless resistivity, η, is given by (1)where η^{d} is the dimensional resistivity (Spitzer resistivity in the case of purely collisional plasma), μ_{0} is the magnetic permeability, and S is the Lunquist number which is typically very large in the solar corona S ~ 10^{12} − 10^{14}.
For completeness, the main properties of the solution are given here; see Craig & Fabling (1996) and Craig et al. (1997) for more detail. After normalisation the governing equations consist of the momentum equation, which in curled form is (2)and the induction equation, (3)with the solenoidal and incompressibility conditions, (4)Here J is the current density and ω is the vorticity in terms of the bulk plasma velocity u. In this normalised form they are (5)The three dimensional analytic solutions of Craig et al. (1995), Craig & Fabling (1996) and Craig et al. (1997) have magnetic and flow fields of the form where the scalar 0 ≤ λ < 1 gives the shear between the B and u fields. The vector field P(x,y,z) is a potential background field of strength α, and Q is a disturbance field of strength B_{s} which gives rise to current in the models.
For comparison with the particle acceleration results at 3D nulls in ideal MHD (Dalla & Browning 2005), we choose the zaxis to be aligned with the spine, and z = 0 as the fan plane. It must be noted that this choice of axis differs from that used by Craig et al. (1997). We study only the proper radial null (Priest & Titov 1996) where the background magnetic field lines in the fan plane lie in the radial direction. This background field is then written as (8)with α giving the sign and strength of the field. For the spine model, the displacement field distorts the fan plane in the zdirection Q_{S} = Z(x,y)ẑ. For the fan model, it distorts the spine in the xdirection (the more general fan case given in Craig et al. 1997, of has not been covered here).
2.1. Spine analytic fields
The disturbance field for the spine model in cylindrical coordinates (r,φ,z) is (9)(Craig et al. 1997) in terms of the confluent hypergeometric (Kummer) function M(a,b,ζ) (Abramowitz & Stegun 1972). The flux pileup factor B_{s} gives the approximate strength of the magnetic field at a dimensionless distance r_{η} from the spine axis, where r_{η} is defined as (10)It is the radius of a cylindrical region centred on the spine axis where resistive effects become significant (Craig et al. 1997).
The form of the displacement field in Eq. (9) is only a solution to the governing equations provided α < 0. This gives frozenin plasma inflow along the fan plane converging on the spine and outflow in the ± z directions away from the null point. The magnetic field in the outer (ideal) region is also directed inwards along the fan plane and outwards along the spine axis. Some representative magnetic field lines are shown in Fig. 1a, the displacement term shears the fan plane at φ = ± π/2 towards the spine axis while the fieldlines in φ = 0,π of the fan plane remain perpendicular to the spine.
To integrate particle trajectories using a test particle model we require the electric field. We calculate this from the uncurled form of Eq. (3) as (11)where P_{r} = αr/2 is the radial part of the potential field.
This electric field is curlfree (as required for steadystate) and so we can calculate the electric potential V to use as a check of energy conservation. This can be found by integrating E = − ∇V to get (12)where f(r) is the radial part of the displacement field in (9), Z(r,φ) = f(r)sinφ.
We can study the behaviour of these fields at small and large distances using the truncated power series and asymptotic formulae for the Kummer function respectively (Abramowitz & Stegun 1972). For all our cases the third argument in the Kummer function is negative. We have for 0 ≤ ξ ≪ 1, (13)and for ξ ≫ 1, (14)in terms of the Gamma function Γ(b). Near the spine axis r ≈ 0, the only contribution to the electric field is from current in the xdirection and (15)The full current distribution is plotted in Fig. 1b, it forms two cylindrical vortex structures that are localised with respect to the resistive region and invariant in the zdirection.
Fig. 1 a) Representative magnetic field lines for the spine model with parameters λ = 0.75, B_{s} = 3.4, α = −2, η = 3 × 10^{3}. The field lines are seeded from the top and base of the spine axis. b) Showing the direction and relative strength of the current in a plane of constant z for the same parameters. Here, is the size of the resistive region centred on the spine axis. 

Open with DEXTER 
At large distances from the spine, the electric field goes as (16)This has the same functional form as the ideal spine solution of Priest & Titov (1996), the subject of previous work on 3D nullpoint particle acceleration by Dalla & Browning (2005, 2006) and Browning et al. (2010). Indeed, by simple choice of parameters we could set the magnetic and electric fields to asymptotically match those of the ideal case. Some care must be taken here as Dalla & Browning (2005) studied positive nulls, where α > 0, and with E(0 < φ < π) > 0. We have opposite sign for both electric and magnetic fields giving the same electric drift inflow quadrants but different sign for the convective electric field. Particles that become nonadiabatic in the external region r ≫ r_{η}, and gain energy parallel the electric field, will rotate about the spine in the opposite direction to those in Dalla & Browning (2005, 2006). In this paper we will only qualitatively compare particle trajectories in the ideal and resistive spine models as an asymptotic match will give rise to unphysical hydromagnetic pressures on the edge of the resistive region r ≈ r_{η} that were absent in the simplified ideal model (see below).
The thermal pressure profile for the spine model can be found from integrating the uncurled form of Eq. (2). It is given in Craig et al. (1997) to be (17)where p_{0} is the gas pressure at the null point, the first term inside the brackets is due to dynamic pressure from the background flow and the other two terms are from balance with magnetic pressure. All terms except for p_{0} are negative, as α < 0, so constraints must be put on the values of α and B_{s} in order to avoid unphysical negative pressures as discussed in Litvinenko et al. (1996); Litvinenko & Craig (1999); Craig et al. (1997) and Craig & Watson (2000). We give some of the arguments here for the sake of completeness (see above references for more detail).
The strong electric field (fast electric drift) simulations for the ideal spine model studied by Dalla & Browning (2005) were characterised by the dimensional value of the electric field E_{0} = 1500 V/m on the r = 1, φ = π/2 boundary (or normalising by suitable solar coronal values, v_{Ae} = 6.5 × 10^{6} ms^{1}, B_{0} = 0.01 T, gives E ≈ 1/40). This value can be equated with Eq. (16). Crucially, to match the external electric field in the resistive model to the fixed amplitude electric field in the ideal spine reconnection model requires the scaling B_{s} ~ η^{1} as η is reduced to suitable solar coronal values (Craig et al. 1997, showed that if we require displacement field at the boundary Z(1,π/2) ~ 1 this also gives B_{s} ~ η^{1}). However, this scaling gives rise to large magnetic pressure on the sheet edge. The maximum of the displacement field occurs at r ≈ r_{η}, where Z(r_{η}) ~ B_{s}, giving magnetic pressure there from Eq. (17) . To avoid negative thermal pressure in the model this requires the null point pressure p_{0} > (Z(r_{η}))^{2} ~ η^{2} which is unphysically large for the values of η considered.
Craig et al. (1997) showed that B_{s} must be limited to a saturation value on r = r_{η}, giving weak electric fields and small amplitude displacement field on the boundary: Z(1) ≪ 1. Also, at r = 1, z ≪ 1, we have dynamic pressure due to bulk fluid inflow p ≈ p_{0} − P^{2}, where P(1) ~ α. We must constrain α ≤ B_{s} or this dynamic pressure will require the gas pressure at the null to be even larger. The maximum value we can take for p_{0} is the largest possible external hydromagnetic pressure p_{e,max} available to drive the reconnection. We follow Craig et al. (1997) and take where the maximum external magnetic field is that of a sunspot at the photosphere, B_{e,max} = 0.3 T. This gives a normalised saturation value B_{s,max} = 30.
So far we do not know the value α should take, but expect that the bulk fluid exhaust from the reconnection region is of the order of the local Alfvén speed. The exhaust on the edge of the current sheet at a global distance from the null, r = r_{η}, φ = π/2, z = 1, is given by where the local Alfvén speed is for our choice of normalisation. As we are not interested in the case where λ = 1 (where there is no shear between the velocity and magnetic fields) we have α ≈ − B_{s} for Alfvénic exhaust. This is the largest magnitude of α we can take without having problems due to dynamic pressure. It also leads to the thinnest current sheet and thus maximises the current density in the resistive region. However, as Craig & Watson (2000) show, the dissipation rate is (18)which has no α dependence as the increase in current density due to resistive region thinning is cancelled by the dependence of the total dissipation volume.
The electric drift velocity in the external region is given by (19)which is very slow when  α  = B_{s}. It is thus necessary to limit the magnitude of α so that results can be obtained with reasonable integration times. For the simulations in Sect. 3 we use B_{s} = 10, α = −0.1, this limits the reconnection exhaust close to the spine current sheet to subAlfvénic speeds.
2.2. Fan analytic fields
The displacement field for the fan model is (20)(Craig et al. 1997). We define z_{η} as (21)the approximate height at which X takes the maximum value, X_{max} ≈ B_{s}. It is a measure of the height of a resistive region centred on the fan plane, z = 0. This form of solution is only valid for α > 0 which gives a positive null point, the field is washed in from the global boundaries at z = ± 1 and it exits the simulation box radially along the fan plane. Some representative magnetic field lines are shown in Fig. 2; the displacement field shears the spine axis as it approaches the fan plane giving rise to strong current inside the resistive region.
The electric field is (22)and the electric potential is (23)where J(z = 0) = X′(0)ŷ is current density at the centre of the sheet. The current only has zdependence; it is infinite in extent in the x and y directions. This is clearly unrealistic, although resistive MHD simulations by Pontin et al. (2007b) find that spinefan reconnecting current sheets formed due to shear flows around a null point spread out along the fan plane in the incompressible limit. Note that analytic multiple null solutions found by Craig et al. (1999) have finite current sheets, avoiding this problem. In our simulations below we consider particle acceleration only within a restricted range of 5 L_{0}, effectively limiting the size of the current sheet.
Fig. 2 Representative magnetic field lines for the fan solution with parameters λ = 0.75, B_{s} = α = 10, η = 10^{6}. The lines are again seeded from the top (solid lines) and base (dashed lines) of the spine. 

Open with DEXTER 
The thermal pressure profile for the fan model is (Craig et al. 1997), (24)However, in this case a displacement field of order unity on the z = 1 boundary, X(1) ~ 1, gives the scaling B_{s} ~ η^{ − 1/4} (Craig et al. 1997). This gives much weaker hydromagnetic pressure on the current sheet edge compared to the spine model, but it is still too large for the values of η considered. Again we saturate B_{s,max} = 30 and we have α ≤ B_{s} to avoid problems from dynamic pressure.
Craig & Watson (2000) show that the Ohmic dissipation rate from the fan model is (25)and so in this case, for fixed (saturated) B_{s}, the maximum dissipation occurs with the thinnest current sheet (the so called optimised solution). The thinnest sheet we can have subject to the dynamic pressure constraint is when α = B_{s} (for any fixed value of λ). Also, as this choice gives the largest current density, it maximises the resistive electric field within the sheet which is interesting for particle acceleration. As above, this choice of α sets the bulk fluid exhaust at x^{2} + y^{2} = 1, z = z_{η} to the local Alfvén speed.
Using the asymptotic approximation (14) we find that the zcomponent of the electric drift that brings the particles to the fan plane is, for x, z ≫ η^{1/2}, (26)for the optimised solution α = B_{s}. This gives electric drift inflow for positive x, z (as P_{z} < 0), and outflow for positive z and negative x. It is much faster than the spine case due to the more favourable scaling with resistivity. There are also fast drift streamlines in the x − y plane that that can be found from the numerical (or approximate analytical) solution of (27)we numerically plot these streamlines on top of the single particle trajectory results.
2.3. Test particle code and parameter choice
We modify the test particle code of Dalla & Browning (2005, 2006, 2008) and Browning et al. (2010) to use the electromagnetic fields given above (from the solutions of Craig et al. 1997). A VariableStep VariableOrder Adam’s method, where the step size is recalculated to properly resolve gyromotion, is used to integrate the relativistic Lorentz equation. (28)where p is the momentum of the particle, γ is the Lorentz factor, q and m are the charge and rest mass and E and B are the analytic expressions for the electric and magnetic fields for each model.
We use the expressions for the electric potential V, given in Eqs. (12) and (23), to calculate the electric potential energy at each time step. With this we verify that the total energy (29)is conserved where ϵ_{k} = (γ − 1)mc^{2}. For each simulation we find that this is conserved up to 5 significant figures. Also, to check the code handles nonadiabatic motion in strong magnetic field gradients of a current sheet we reproduce the results of Speiser (1965), including the ejection time for the case with background field.
We choose L_{0}, the normalising length scale, to be L_{0} = 10^{4} m for global simulations to keep integration times short. This size of simulation box can be considered as a local region around the null at which the linear background field and flow in Eq. (8) are good approximations. We use a larger value of L_{0} = 10^{6} m for simulations where particles are initially distributed within the current sheet, as velocities are typically much faster here. Note that a change in L_{0} also changes the value of η as given in Eq. (1).
All magnetic fields mentioned in Sect. 2 have dimensions of B_{0} = 0.01 T, typical for the solar corona. We set v_{0} = v_{Ae} = 6.5 × 10^{6} ms^{1} (corresponding to a number density n_{0} = 1.126 × 10^{15} m^{3}). All dimensionless times quoted are in terms of the gyroperiod, T_{ω} = 2πm/(qB_{0}).
To examine single particle trajectories in both models we choose values of B_{s} = 10, λ = 0.75, η = 10^{6}. This η value is rather large, towards the highest possible anomalous resistivities (with L_{0} = 10^{4} m), but is useful to observe particles entering the current sheet. We vary all three of these parameters to produce scalings of energy gain and drift times, where we use values as low as those expected in purely collisional plasma i.e. η = 10^{12}.
3. Results
3.1. Spine global trajectories
Initially, we place a distribution of 5000 protons with Maxwellian velocities of temperature T = 10^{6} K (86 eV) in the spine model fields. The protons have positions from a uniform random distribution at a global distance x^{2} + y^{2} + z^{2} = 1 from the null point. We only discuss here protons that start in the upper right inflow region of longitude 0 < φ < 180° and latitude 0 < β < 90° (here φ = 0 is the xaxis and β = 0 is the fan plane).
Figure 3 shows the final spatial distribution and energies of the particles at time t = 1.6 × 10^{6} T_{ω,p} ≈ 10 s, at which the energy spectrum becomes approximately steadystate. Those protons starting in the lower left inflow region have final distributions as in Fig. 3 after reflections in both φ = 0 and β = 0, apart from statistical differences. The parameters used here are B_{s} = 10, η = 10^{6}, α = −0.1. This value of α limits the bulk flow exhaust speed to be subAlfvénic, but it increases the electric drift speed in the external region (see Eq. (19)) due to weaker magnetic field on r = 1. This gives reasonable simulation times, but there are still some particles in the upper right inflow quadrant at the end of the simulation.
Fig. 3 Angular distribution of protons from the null point for spine model at t = 1.6 × 10^{6} T_{ω,p}, at which the energy spectrum is steady state, for parameters λ = 0.75, B_{s} = 10, α = −0.1, η = 10^{6}. The initial distribution was Maxwellian at T = 86 eV in the upper right inflow region. 

Open with DEXTER 
There are two main populations of accelerated particles. The population labelled “A” in Fig. 3 is close to the fan plane,  β  ≲ 10°, with energy ϵ_{k} ≳ 1 keV, and with longitude − 90° ≲ φ ≲ 90° comprising of about 8% of the total proton number. The maximum particle energy of this population is about 15 keV. Note that the current in the spine axis is aligned with φ = 0: through the centre of this population. There are also some high energy protons scattered at large positive latitudes for φ ≲ 0, and at large negative latitudes for φ ≳ 0. To look more closely at what is happening here we will choose a typical proton from this population and follow its trajectory below.
For those particles that have crossed the fan plane, β = 0, into the lower right outflow quadrant, the spatial and energy distribution looks similar to the ideal spine case in Dalla & Browning (2006). The accelerated population which has ϵ_{k} ≳ 1 keV and β ≲ − 85° is labelled “B”. This population is about 6% of the total protons in the simulation and the maximum kinetic energy in this population is ϵ_{k,max} ≈ 12 keV. The angular distribution differs slightly with the ideal case in that there are few particles found between the latitudes − 70° < β < − 85°; particles appear to be closer to the negative spine axis in the resistive case.
Fig. 4 Typical proton global trajectory from population “A” for parameters λ = 0.75, B_{s} = 10, α = −0.1, η = 10^{6}. The particle is taken from the many particle simulation, having initial position x = (−0.52,0.80,0.29) and velocity v = (−0.0044,0.0013, −0.0088)v_{Ae}. The magnetic field lines (thin dashed) are a projection of the field from the plane of the trajectory φ ≈ 120°. Inset shows the 3D trajectory of the proton as it crosses the spineaxis, the solid line in the centre is the line B_{z}(x,y,0.95) = 0. 

Open with DEXTER 
A typical proton trajectory from population “A” is shown in Fig. 4. The proton which starts at (x_{0},y_{0},z_{0}) = (−0.52,0.80,0.29) in the upper right hand inflow quadrant initially moves away from the null, but mirror bounces and travels back towards the spine along the fan plane. The electric drift speed increases towards the spine causing the proton to enter the resistive region, which has radius r_{η} ≈ 0.01, about the spine axis. It enters at (x,y,z) ≈ (−0.01,0,0.95) after t = 3 × 10^{5} T_{ω,p} ≈ 2 s (inset). At this point the proton becomes unmagnetised as the gyroradius becomes comparable to the lengthscale of magnetic field gradient ρ/L_{∇B} > 1 (we typically find that gyromotion starts to break down when ρ/L_{∇B} ≳ 10^{2}).
Fig. 5 Typical proton trajectory from population “B” in the many particle simulation, with initial position x = (−0.54,0.78,0.31)L_{0} and velocity v = (−0.004, −0.006,0.002)v_{Ae}. The dashed lines show the projection of the magnetic field from the plane of motion, φ ≈ 110°, onto the y − z plane. Inset shows the motion in the x − y plane close to the spine axis. The purple arrows show the direction and relative magnitude of the gradient drift velocity and the dashdotted lines show contours of the electric potential, with the intersecting tick mark indicating lower potential to the right. 

Open with DEXTER 
The proton is then directly accelerated in the xdirection, parallel to the current at the spine, with du/dt ≈ qE_{0}/m as it crosses r = 0. For small displacements in the ydirection a strong Lorentz force due to the B_{z} field returns it to y = 0 line. These oscillations are Speiserlike (Speiser 1965) with frequency approximately ω ∝ t^{1/2}.
The Speiserlike motion finishes and the first gyrations start (not shown) when the proton reaches r ≈ 5r_{η} at which ρ/L_{∇B} ≲ 1. However, the energy gain of ϵ_{k} ≈ 11 keV is localised to within x ≈ 2r_{η}, during which the trajectory does not deviate much from the xdirection (note the yaxis scale in the inset of Fig. 4). In effect, the proton has left the localised current sheet while unmagnetised but before it can be ejected by the background field components, in contrast to 2D current sheet configurations with weak guide field (e.g. Speiser 1965; Litvinenko 1996). Figure 4 may give the impression that the particle is being ejected, however, this is just the centre of the Speiserlike oscillations following the B_{z}(x,y,z = const.) = 0 line (which here is not straight as in the usual 2D configurations). This behaviour is evident considering the F = qv × B force for the unmagnetised proton if B_{z} is the dominant component of the magnetic field.
After the proton becomes remagnetised at r ≈ 5r_{η} it has weak electric drift, v_{E} ≪ v_{ω}. It follows the fieldlines closely and mirror bounces, travelling back towards the spine: there the proton is taken up to high latitude before it bounces again. This mirror bouncing is the reason for the “scattered” accelerated protons in Fig. 3, some of which are at large latitudes.
A typical particle trajectory chosen from population “B” is shown in Fig. 5. The proton starts at ( − 0.54,0.78,0.31) and drifts towards the spine but bounces and crosses the fan plane instead. It exits the simulation box down the base of the spine axis, reaching an energy ϵ_{k} = 6.72 keV as it crosses z = −5. As there is no electric field in the zdirection, the energy gain must occur due to motion in the x − y plane, which is also shown in Fig. 5. The proton enters the region close to the spine axis parallel to a contour of the electric potential, but then drifts across the contour due to strong gradient drift. While the proton gains energy, the gradient drift is larger than the electric drift by a factor of 2 with the latter directed inwards towards the current sheet. The proton is stopped as it reaches z = −5L_{0} which we do consistently throughout these simulations. At the time of stopping it is actually losing energy as it recrosses the same electric potential contours. However, some other protons from the manyparticle simulation in Fig. 3 reach the current sheet at low latitudes, gaining higher energy.
The energy spectrum for the spine simulation is shown in Fig. 6. If protons cross the R = 5L_{0} spherical boundary we use the energy at the instant of crossing (if this is not done some protons reach order ~ 10^{2}L_{0} which becomes unrealistic as the background field increases without bound away from the null, also causing the timestep to decrease and simulation time to increase). The initial Maxwellian spectrum hardens to a broken power law with maximum energy of about ϵ_{k} ≈ 15 keV. This maximum energy can be understood as the difference in potential energy across the spine current sheet, ϵ_{k} ~ qEx_{acc} where E ≈ E_{0} ≈ ηB_{s}/r_{η} [v_{Ae}B_{0} ] and x_{acc} ≈ 2r_{η} [L_{0} ] is the acceleration distance (from − r_{η} ≲ x ≲ r_{η}), as the reconnection electric field drops off quickly for  x  > r_{η}. For the parameters used, this gives ϵ_{k} ≈ 13 keV. This approximate expression has no dependence upon the parameter α, so the limiting of  α  < B_{s} should not have a large effect on this result.
Fig. 6 Energy spectrum from the many particle simulation for protons in the spine model, with parameters λ = 0.75, B_{s} = 10, α = −0.1, η = 10^{6}. For particles leaving the R = 5 sphere, the energy at the time of crossing is used. 

Open with DEXTER 
3.2. Fan global trajectories
The many particle simulation for the fan model is shown in Fig. 7 for the optimised solution B_{s} = α = 10, with η = 10^{6}, λ = 0.75. The initial distribution has thermal energy ϵ_{k} = 86 eV with uniform random position in the upper inflow quadrant −90° < φ < 90° and 0 < β < 90°. The final angular distribution is taken from when the proton distribution reaches a steady state in energy at t = 4000 T_{ω,p} ≈ 0.025 s. This is more than two orders of magnitude faster than the spine model for similar parameters (even after the spine drift was increased by limiting α < B_{s}) as the external electric drift (Eq. (26)) scales more favourably with the resistivity. Protons that cross the R = 5L_{0} spherical boundary from the null point before this time are stopped and the energy and angular position at time of crossing is used. The t = 0 angular distribution has some structure in terms of final energy gain as the initial random thermal velocities are dominated by the strong electric drift.
Within this structure there is some asymmetry in φ. Indeed, we do not expect symmetry between particles drifting clockwise and anticlockwise about the null as in the ideal case (Dalla & Browning 2006) now that there is a current in the y direction. Those protons with ϵ_{k} ≈ 10^{7} eV at φ ≈ − 20° (the yellow vertical band to the left of the green vertical band in Fig. 7a) do not enter the current sheet but gain high energy, as they are unmagnetised slightly with ρ/L_{∇B} ~ 10^{3}, due to very fast electric drifts close to the sheet. Here the first adiabatic invariant, the constancy of μ, is also violated.
Fig. 7 a) Angular distributions of protons in fan model at t = 0, with initial temperature T = 86 eV. The xaxis is φ = 0 and the fan plane is β = 0. Protons are coloured by the final energy at t = 4000. Parameters used are λ = 0.75, B_{s} = α = 10, η = 10^{6}. b) Angular distributions at time t = 4000T_{ω,p} when the energy spectrum has reached steady state. 

Open with DEXTER 
Typically, the high energy protons of Fig. 7 start either close to the xaxis at low to mid latitudes (about 7% of the total number at latitude β ≳ 1° with final energy ϵ_{k,fin} ≳ 10 MeV), or they start at very low latitude close to the fan plane ( < 1% of total at β ≲ 1° and ϵ_{k,fin} ≳ 10 MeV). At t = 4000 T_{ω,p} these energetic protons are found at β ≈ 0 either side of φ = 90°; the yaxis in the direction of the fan current.
At t = 4000 T_{ω,p} there are a small number of high energy protons scattered at high latitudes (about 0.1% with ϵ_{k} > 10 MeV). These enter the current sheet temporarily at negative longitude far from the null point, but exit again without any Speiserlike motion. They become slightly unmagnetised, with maximum ρ/L_{∇B} ≈ 10^{2}, following complicated trajectories. As they are not typical they are not investigated further in the external region, but the behaviour within the current sheet is discussed below.
Fig. 8 Two typical proton trajectories from the many particle simulations in the fan model. Proton “1” is represented by a thin line with initial position (0.86,0.41,0.30), and proton “2” by a thick line with initial position (0.8,0.003,0.6). The parameters are λ = 0.75, B_{s} = α = 10, η = 10^{6}. The solid lines are representative magnetic field lines (seeded from the top of the spine axis and projected into the 2D planes) and the arrows show the direction and relative magnitude of the electric drift velocity. a) In the x − y plane, where the electric drift arrows are from the edge of the current sheet z = z_{η}. b) In the x − z plane close to the current sheet, where the electric drift arrows are plotted on y = 0. The initial positions are not shown in this plane. 

Open with DEXTER 
Figure 8 shows the trajectory of two typical protons taken from the simulation. Proton “1” starts at (x_{0},y_{0},z_{0}) = (0.86,0.41,0.30) and drifts around the null point due to the strong azimuthal electric drift. Although it drifts down towards the current sheet, it reaches a minimum height of z ≈ 15 z_{η} before it flows into the outflow quadrant, not entering the sheet. The main velocity contribution is electric drift as it moves around the null point, but v_{ ∥ } becomes dominant as the particle exits the simulation box parallel to the negative xaxis. The first adiabatic invariant is not violated, μ = const. and the maximum ρ/L_{∇B} ~ 10^{4} at closest point of approach to the sheet. The proton is strongly magnetised throughout. Despite not reaching the current sheet the energy gain is still considerable, reaching 0.5 MeV as it crosses the R = 5 sphere.
Particle 2 starts at (0.8,0.003,0.6)L_{0}. The azimuthal electric drift is weak close to the xaxis and the proton drifts down to the fan current sheet. It enters the sheet at (x,y) = (−0.02,0.15) and becomes unmagnetised: ρ/L_{∇B} > 1 and μ is not conserved. We observe Speiserlike oscillations as the proton is accelerated in the ydirection. At t = 0.846 ms after entering the sheet, it passes out of the simulation box at R = 5. Here, the particle is still within the sheet with v_{ ∥ } = 0.36c and ϵ_{k} = 67 MeV. Using this time period in the direct acceleration formula, y = qE_{0}t^{2}/2m, with the electric field on z = 0, from (22), gives y ≈ 5. Thus the proton is directly accelerated in the current sheet for the entire length of the simulation box. However, this motion is not Speiserlike throughout as ρ/L_{∇B} < 10^{2} when the proton reaches y = 1.5L_{0}. The proton reaches a global distance in the ydirection and becomes magnetised by the background B_{y} component of the magnetic field, which acts as a kind of guide field. When the simulation is run without stopping the particle at R = 5, the proton is not ejected from the current sheet throughout the whole simulation time t = 4000 T_{ω,p}.
This particle enters the current sheet at a distance R ≈ 0.15 from the null point; however, this distance is not typical for the many particle simulation in Fig. 7. In the simulation 9.3% of the total particles reach the current sheet, after a mean time of about 800 T_{ω,p}. The average distance from the null point of particles entering the sheet is R ≈ 2.2; some remain magnetised by the background magnetic field inside the sheet.
Fig. 9 Energy spectrum for the many particle fan simulation for protons with parameters λ = 0.75, B_{s} = α = 10, η = 10^{6}. For particles leaving the R = 5 sphere, the energy at the time of crossing is used. 

Open with DEXTER 
The energy spectrum for the fan simulation is shown in Fig. 9. Almost all the protons are accelerated as the fast electric drift speed in the fan model is typically much larger than initial thermal velocity. The energy spectrum appears steady state by t = 4000T_{ω,p}. Between ϵ_{k} ≈ 10^{5} eV and ϵ_{k} ≈ 10^{7.5} eV the spectrum is approximately power law shape: f(E) ∝ E^{ − γ} with γ ≈ 1.5. It is interesting to note that this is the same spectral index that is expected for direct acceleration without a guide field in a 2D current sheet (see for e.g. Heerikhuisen et al. 2002). However, for the case of Fig. 9, the protons in this energy range have not entered the current sheet. Instead the energy gain occurs in the external region due to nonuniform drifts parallel to the electric field. The protons that do enter the current sheet gain the highest energies and can be seen as a flat tail at the hard end of the spectrum. We do not give a spectral index for this population as it has not reached steady state: as these protons are trapped in the sheet, the energy gain by direct acceleration depends upon the size of the simulation box. As a test we repeat the simulation but stop the protons at a spherical surface of radius R = 10 from the null point, instead of R = 5 that has been used consistently throughout these simulations. Now the “flat tail” at 10^{7.5 − 8} eV in Fig. 9 becomes a “bump on tail” centred at 10^{8} eV (not shown) disconnected from the main distribution. The population of protons that is trapped in the sheet as it crosses R = 5 due to the strong “guide field” remains trapped at R = 10 where B_{y}(y) has doubled in strength. We do note that for the rest of the protons, which have not entered the current sheet, the energies remain approximately the same when stopped at R = 10. This can be understood as follows: the protons that do not enter the current sheet typically leave the simulation box parallel to the negative xaxis, so they do not cross contours of the electric potential, V = V(y,z) in Eq. (23), after they have drifted around the null point and are in outflow.
3.3. Fan current sheet trajectories
The simulations considered thus far concern proton trajectories starting from the external region, at a distance R = 1 from the null point. However, most of the protons entering the current sheet do so far from the null point. In the following, protons are initially distributed within the fan current sheet close to the null, to study the transition from nonadiabatic to adiabatic motion.
Firstly, we place particles within the sheet so that they are initially unmagnetised by the B_{y}(y) component of the background field. They are magnetised only by the strong B_{x}(x,z). The protons are uniformly distributed in the area  x  < 1; y = 0;  z  < z_{η} with initial thermal energy T = 86 eV. Figure 10 shows the position of 2000 protons at t = 2500 T_{ω,p} (a), and t = 17 500 T_{ω,p} (b), during this simulation for the parameters λ = 0.75, B_{s} = α = 5, η = 10^{8}. We increase the dimensional box length to L_{0} = 10^{6} m as velocities in the current sheet are typically fast, giving reasonable integration times. This makes our results more comparable to the approximate analytic solutions of Litvinenko (2006). Note that η decreases due to the increase in L_{0} in Eq. (1). We again artificially stop the particles as they cross the R = 5 spherical surface.
Fig. 10 Proton positions after initial distribution within the fan current sheet, such that  x_{0}  < 1,  y_{0}  = 0, . Parameters used are λ = 0.75, B_{s} = α = 5, η = 10^{8}, L_{0} = 10^{6} m. The dashed lines are representative magnetic fieldlines inside the current sheet (note the difference in scale of the zaxis). The solid black line is the line (x_{1},0,z_{1}) such that B_{x}(x_{1},0,z_{1}) = 0. The particles are stopped at R = 5 from the null point. 

Open with DEXTER 
At t = 2500 T_{ωp} most of the protons are strongly magnetised by the B_{x}(x,z) magnetic field. Inside the current sheet,  z  < z_{η} we can use Eq. (13) to get approximate expressions for the electric and magnetic fields, E_{z} is small except at global distance in y (see below).
For a proton starting at x = 0, y = 0, z = z_{η}, on the edge of the current sheet, the background components of the magnetic field are negligible. The proton drifts towards the vertical centre of the sheet v_{Ez} ≈ − (η/z) ẑ [v_{Ae} ] . It becomes unmagnetised at the fan plane, z ≈ 0, close to the null point and is directly accelerated in the ydirection. We compare this trajectory to the analytical WKB solutions of Litvinenko (2006). The ejection time for a nonrelativistic proton that is unmagnetised close the null point in the fan current sheet, x ≈ 0, z ≈ 0, for our parameters is (32)(Litvinenko 2006) provided that the proton remains within the nonadiabatic region and the displacement magnetic field gradient is much stronger than the gradient from the background component, B_{s}/z_{η} ≫ λα. The second assumption is valid for our simulation; however, we do not observe proton energy gain limited by ejection in these simulations. To understand this, we consider the distance travelled in the ydirection during this time, (33)which we compare with size of the nonadiabatic region from the null in this direction. The particle begins to be remagnetised by the background field at a global distance y^{ ∗ } such that (34)where v(y) is a typical proton velocity and ω_{By(y)} is the gyrofrequency of a particle gyrating around B_{y}(y). We use v(y) = ^{(}2qE_{y}y/m^{)}^{1/2} from direct acceleration (if we use v(y) = E_{y}/B_{y} the value for y^{ ∗ } differs by 2^{1/3}), assuming that there was no initial yvelocity and the particle entered the sheet at y ≈ 0. We recover the result of Litvinenko (2006), that in dimensional form (35)The ratio of these two distances is (36)where we have ignored factors of order unity. The ratio of the two timescales is the square root of this. There is little gyroturning for protons starting close to the null point as this ratio is necessarily small for the fan current sheet. The proton is magnetised by the B_{y}(y) “guide field” and trapped in the sheet; the energy gain is only bounded by the length of the sheet.
Figure 10 also shows the more general case of protons starting at y = 0,  z  < z_{η} and at a global distance in x. These protons drift vertically until they reach the diagonal line where B_{x}(x,z) = 0, at which they become unmagnetised and accelerated. We do appear to see some gyroturning for protons starting at  x  ≈ 1. This is probably due to the strong component of the Lorentz force, v_{y}B_{z}, that acts to turn the trajectory to the xdirection. For particles starting at x = 0 the B_{z} magnetic field switches sign during the zoscillations, but at x = 1 the proton is unmagnetised below the fan plane and B_{z} stays positive. The proton is turned in the positive xdirection but is quickly magnetised by the guide field when it reaches a distance of about y^{ ∗ } (see Eq. (35)). In Fig. 10 it can be seen that particles are accelerated radially outwards from the null. They continue to gain energy as they become magnetised about the background field P, on a fieldline with a parallel component of the electric field E_{} = E·P/  P  .
We artificially stop the protons at R = 5 L_{0} from the null point. At t = 50 000 T_{ω,p} all of the protons in the simulation have reached this distance without being ejected and we fit the energies of the particles by the expression (37)for the optimised solution α = B_{s}, where φ is the azimuthal angle (φ = 90° is parallel to the current). Figure 11 shows the energies of 5000 protons in three simulations with identical setup to Fig. 10 but with different values of η and B_{s}. This expression (thin line) fits the energies of simulated particles (circles) as they cross R = 5 L_{0} very well.
Fig. 11 Energy distribution of particles as they cross the R = 5 boundary, where initial position is within  x_{0}  < 1,  y_{0}  = 0, . Here, L_{0} = 10^{6} m and the results for different values of η and B_{s} are plotted. The solid points are protons from the three simulations and the thin lines show the sinφ relationship in Eq. (37). 

Open with DEXTER 
Fig. 12 Proton positions after initial distribution within the fan current sheet such that  x_{0}  ,  y_{0}  < 1 and . Parameters used are λ = 0.75, B_{s} = α = 5, η = 10^{8}, L_{0} = 10^{6} m. Particles circled in black are those that start in y < 0 and cross z = z_{η} ≈ 9.6 × 10^{5} before t = 19 000 T_{ω,p}. 

Open with DEXTER 
In Fig. 12 we place 5000 protons in the fan current sheet with initial position in  z  < z_{η}, − 1 < x,y < 1 so that they are initially magnetised by the “guide field” B_{y}(y). This is the more general case, as protons reaching the current sheet from the external region will not typically do so at y ≈ 0. The protons that do not start close to y = 0 are directly accelerated without the initial drift phase. By t = 19 000 T_{ω,p} all of the protons have left the simulation box; either through the R = 5 boundary, or through the edge of the current sheet  z  = z_{η}. The particles that cross  z  = z_{η} in y > 0 start close to the edge and leave due to initial thermal velocity. However, those starting with y < 0 are ejected from well within the current sheet. These protons (19.7% of total number) are circled in Fig. 12. Typically they remain magnetised, with ρ/L_{∇B} in the range 10^{4} − 10^{2}. They are not ejected due to gyroturning in the sense of Speiser (1965) as this requires nonadiabatic motion. The trajectories in the y − z plane seem to follow the magnetic field lines closely, although they have strong electric drift from the E_{z} component of the electric field. Within the current sheet,  z  < z_{η} the truncated power series in Eq. (13) gives the zcomponent of the electric field from Eq. (22) as (38)which is stronger than the current electric field (30) for global y and z ≠ 0. However, it only contributes to strong electric drift (not shown) in negative xdirection for protons in the upper half of the sheet 0 < z < z_{η}, and in the positive xdirection for − z_{η} < z < 0. This electric field also contributes to v_{Ey} but this is dominated by the direct electric field acceleration.
The protons that are not ejected from the current sheet have an approximate sinusoidal dependence in kinetic energy gain, given by Eq. (37) (there is a thicker spread of points about the predicted lines than in Fig. 11 due to differences in initial potential energy).
3.4. Scalings
In the fan global simulation of Fig. 7 we chose the optimised parameters η = 10^{6}, B_{s} = α = 10 and λ = 0.75. With consideration to the large variation in both scale and behaviour in a given distribution of flares, it is interesting to see how the results of this simulation scale when the simulation parameters are varied.
The current sheet within the fan model is the most effective way to accelerate the particles. Thus, it is interesting to study the effect of varying parameters on the fraction of protons that enter the current sheet from the external region, and the average time taken to drift there from an initial position on the R = 1 sphere. These scalings are shown in Fig. 13. They are from simulations of 5000 protons at T = 1 MK starting at the upper inflow region at R = 1. We define the current sheet as for the fan model, although we note that not all of the protons reaching this height become nonadiabatic.
Fig. 13 Percentage of total particles (+) reaching current sheet at z = z_{η} from the external region (R = 1 in the inflow quadrant), and the mean time taken (*), for different values of η, B_{s} and λ. Each data point is from a many particle simulation with initial Maxwellian distribution (T = 86 eV) of 5000 protons. The setup is the same as that in Fig. 7. a) For different values of η with fixed B_{s} = α = 5, λ = 0.75. The solid line is a least squares fit to the points. b) For different values of B_{s} (with B_{s} = α) for fixed η = 10^{8}, λ = 0.75. The solid line is a least squares fit. c) For different λ with fixed B_{s} = α = 5, η = 10^{8}. No curve was fit. 

Open with DEXTER 
Fig. 14 Scalings of steadystate energy spectra for global fan simulation. a) For different values of η with fixed B_{s} = α = 5 and λ = 0.75. The time taken to reach steady state was t = 8 × 10^{4}, t = 2 × 10^{4} and t = 6.4 × 10^{3} T_{ω,p} for η = 10^{10}, η = 10^{8} and η = 10^{6} respectively. b) For different values of B_{s} = α with fixed η = 10^{8} and λ = 0.75. Steadystate was reached at t = 2 × 10^{5}, t = 2 × 10^{4} and t = 5 × 10^{3} T_{ω,p} for B_{s} = 1, B_{s} = 5 and B_{s} = 30 respectively. c) For different values of λ with fixed B_{s} = α = 5 and η = 10^{8}. Time to steady state is t = 10^{5}, t = 2 × 10^{4} and t = 8 × 10^{3} T_{ω,p} for λ = 0.9, λ = 0.75 and λ = 0.3 respectively. 

Open with DEXTER 
The average time taken for the particles to reach the current sheet gives a measure of the external electric drift speed. The approximate drift scaling of Eq. (26), , is in reasonable agreement with these drift times.
The fraction of particles reaching the sheet increases typically with increasing η and decreasing B_{s}. Note that the size of the current sheet which we use to produce the scalings has the dependence , although this does not fully explain the result as the trends are not simple power laws. Figure 14 shows how the energy spectra vary with these parameters. As might be expected the spectra shift to the right for an increase in both B_{s} and η. Both the convective electric field (and so external electric drift) and the direct electric field within the sheet increase with larger values of these parameters.
Up until now we have used λ = 0.75 as a constant in all of the simulations. This parameter has been typically left as constant in the calculation of MHD energy dissipation scalings (Craig et al. 1997; Craig & Watson 2000) for the fan model as it can only be varied within an order one range. However, it has a large effect on the efficiency of the fan model for particle acceleration (see Figs. 13c and 14c). Varying λ within 0 ≤ λ < 1 has a comparable effect on the fraction of particles reaching the current sheet as varying η by six orders of magnitude. Also, as shown in Fig. 14, decreasing λ shifts the energy spectrum to higher energy and decreases the time taken to reach steadystate. These effects can be explained somewhat by an increase in external drift speed, although varying λ also has an effect on other quantities such as the current sheet height.
4. Summary and discussion
We investigated test particle motion in the electromagnetic fields of Craig et al. (1995); Craig & Fabling (1996) and Craig et al. (1997), that are solutions to the steadystate, incompressible and resistive MHD equations at a 3D null point. The study was carried out by modifying the code of Dalla & Browning (2005, 2006, 2008) and Browning et al. (2010). We considered initially Maxwellian (T = 86 eV) distributions of protons starting at a global distance R = L_{0} from the null point in resistive spine reconnection, where the electric current is within a thin cylinder about the spine axis, and resistive fan reconnection, with a current sheet in the fan plane. When the energy spectrum from the simulations reached steadystate we find the final angular position of the particles from the null and their energy distribution. We identified different populations of accelerated particles and, to understand the acceleration mechanism, examined a typical single particle trajectory in each case. For the fan model we ran additional simulations with the particles initially distributed within the fan current sheet, to study the effect of the null point on directly accelerated particles. We consider two cases, where particles are firstly unmagnetised and secondly magnetised in their initial position by a background “guide field”. Finally, we show how the external drift times and energy spectra for the fan model scale when treating λ, B_{s} and η as free parameters (for the optimised solution α = B_{s}, corresponding to the thinnest current sheet, Craig & Watson 2000).
We found that the spine model, which gave promising acceleration results in the ideal case (Dalla & Browning 2008), is much less effective when resistive effects are included (at least for this specific resistive model). The electric drift in the external region is weak, scaling with resistivity as v_{E}(r ≫ r_{η}) ~ η, giving very long drift times for protons to reach the spine axis. We find that there are two populations of accelerated particles. One of these escapes the simulation box down the base of the spine axis, similar to the proton jet found in the ideal model (Dalla & Browning 2006), and the other is close to the fan plane, where particles have crossed the current sheet in the spine axis. The energy gain for particles that reach the current sheet is low. The main limiting factor is the small electric potential energy difference across the current sheet, due to localisation of the reconnection electric field to a small cylinder about the spine axis. The apparent contrast with the ideal results arises from parameter constraints. In the ideal regime, the magnitude of the electric field for the spine and fan models was set to equal strength in the external region, at a global distance from the null. The electric field falls steeply as 1/r in the external region of the spine model which gives very strong acceleration close to the spine in the ideal case. However, when the pressure constraints in the resistive model are taken into account, namely the limiting of the displacement field on the edge of the current sheet to avoid unphysical magnetic pressures (Craig et al. 1997), this 1/r dependence gives very weak electric field, and slow drift, in the external region.
We found much higher proton energies in the resistive fan model for similar parameters. For η = 10^{6}, B_{s} = α = 10 and λ = 0.75 we find the energy spectrum from a distribution of protons starting with thermal energy at R = L_{0} from the null becomes power law at steady state for particles accelerated due to fast nonuniform drifts (index ≈ − 1.5 between 10^{5} eV and 10^{7.5} eV). The electric drift in the external region is much quicker than the spine model, v_{E} ~ η^{1/4}. It accelerates all of the particles in the simulation as it is faster than the initial random velocities associated with thermal motion at T = 86 eV. We find that the population with the highest energy gain, up to 0.1 GeV, corresponds to protons that have entered the fan current sheet. The energy gain for these protons is not limited by ejection due to unstable motion as they are remagnetised within the current sheet by a “guide field”. The upper bound in energy gain is only limited by the electric potential energy, associated with the length of the current sheet. However, we find that a number of protons that enter the current sheet upstream of the null point can be ejected while remaining magnetised. This is due to the geometry of the background field lines, namely that they diverge at the null. We will study this effect in the future when we consider electrons. Browning et al. (2010), and Guo et al. (2010) show that electrons remain magnetised at a closer distance to the null point, which may give a difference between the number of electrons and protons ejected in this manner.
We find that the parameter λ, which gives the degree of shear between the magnetic and velocity fields (such that 0 ≤ λ < 1) has a large effect on the final energy spectrum of protons in the fan model. In the limit of λ = 0 the magnetic field in the fan model is annihilated. As we expect magnetic field to still exist in the reconnection site after a topological change it would be more likely that λ ≈ 1.
In these simulations we have neglected the electromagnetic effects of the nonthermal particles onto the background fields. This is typical of the test particle approach, where it is assumed that the number of particles in the current sheet is a small fraction of the total number. For a large range of parameters in the fan model (see Fig. 13) this fraction is typically less than 5% of the total number of particles starting in the inflow region. To estimate the strength of the magnetic fields from these particles, and the polarisation electric field from any charge seperation, it is necessary to also consider electrons which we will do in future work. A fully selfconsistent approach, e.g. using Particle In Cell simulations, is computationally expensive at present, particularly in fully 3D geometries due to the large dynamic range of spatial scales. We have also neglected compressibility, a simplification used to get the analytic solutions (Craig et al. 1997). It would be interesting in future to include both timedependence and compressibility, by using electromagnetic fields from MHD simulations, to see how particles behave in so called spinefan reconnection (Pontin et al. 2007a,b). It would also be interesting to study the effect of fan plane curvature on acceleration within the fan current sheet, by using the analytic solutions of Watson & Craig (2002) that are based in cylindrical geometry rather than cartesian.
Acknowledgments
We are grateful to the anonymous referee for the kind comments and A.S. thanks STFC for studentship support.
References
 Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions (Dover Publications) [Google Scholar]
 Aulanier, G., DeLuca, E. E., Antiochos, S. K., McMullen, R. A., & Golub, L. 2000, ApJ, 540, 1126 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, J. C. 1971, Sol. Phys., 18, 489 [NASA ADS] [CrossRef] [Google Scholar]
 Browning, P. K., Dalla, S., Peters, D., & Smith, J. 2010, A&A, 520, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bulanov, S. V., & Cap, F. 1988, Soviet Ast., 32, 436 [NASA ADS] [Google Scholar]
 Craig, I. J. D., & Fabling, R. B. 1996, ApJ, 462, 969 [NASA ADS] [CrossRef] [Google Scholar]
 Craig, I. J. D., & Fabling, R. B. 1998, Phys. Plasmas, 5, 635 [NASA ADS] [CrossRef] [Google Scholar]
 Craig, I. J. D., & Henton, S. M. 1995, ApJ, 450, 280 [NASA ADS] [CrossRef] [Google Scholar]
 Craig, I. J. D., & Litvinenko, Y. E. 2002, ApJ, 570, 387 [NASA ADS] [CrossRef] [Google Scholar]
 Craig, I. J. D., & Watson, P. G. 2000, Sol. Phys., 191, 359 [NASA ADS] [CrossRef] [Google Scholar]
 Craig, I. J. D., Fabling, R. B., Henton, S. M., & Rickard, G. J. 1995, ApJ, 455, L197 [NASA ADS] [CrossRef] [Google Scholar]
 Craig, I. J. D., Fabling, R. B., & Watson, P. G. 1997, ApJ, 485, 383 [NASA ADS] [CrossRef] [Google Scholar]
 Craig, I. J. D., Fabling, R. B., Heerikhuisen, J., & Watson, P. G. 1999, ApJ, 523, 838 [NASA ADS] [CrossRef] [Google Scholar]
 Dalla, S., & Browning, P. K. 2005, A&A, 436, 1103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dalla, S., & Browning, P. K. 2006, ApJ, 640, L99 [NASA ADS] [CrossRef] [Google Scholar]
 Dalla, S., & Browning, P. K. 2008, A&A, 491, 289 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Des Jardins, A., Canfield, R., Longcope, D., Fordyce, C., & Waitukaitis, S. 2009, ApJ, 693, 1628 [NASA ADS] [CrossRef] [Google Scholar]
 Dreicer, H. 1959, Phys. Rev., 115, 238 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Fletcher, L., Metcalf, T. R., Alexander, D., Brown, D. S., & Ryder, L. A. 2001, ApJ, 554, 451 [NASA ADS] [CrossRef] [Google Scholar]
 Gordovskyy, M., Browning, P. K., & Vekstein, G. E. 2010a, A&A, 519, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gordovskyy, M., Browning, P. K., & Vekstein, G. E. 2010b, ApJ, 720, 1603 [NASA ADS] [CrossRef] [Google Scholar]
 Guo, J., Büchner, J., Otto, A., et al. 2010, A&A, 513, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamilton, B., Fletcher, L., McClements, K. G., & Thyagaraja, A. 2005, ApJ, 625, 496 [NASA ADS] [CrossRef] [Google Scholar]
 Hannah, I. G., & Fletcher, L. 2006, Sol. Phys., 236, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Heerikhuisen, J., & Craig, I. J. D. 2004, Sol. Phys., 222, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Heerikhuisen, J., Litvinenko, Y. E., & Craig, I. J. D. 2002, ApJ, 566, 512 [NASA ADS] [CrossRef] [Google Scholar]
 Ishikawa, S., Krucker, S., Takahashi, T., & Lin, R. P. 2011, ApJ, 737, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Krucker, S., Hudson, H. S., Glesener, L., et al. 2010, ApJ, 714, 1108 [NASA ADS] [CrossRef] [Google Scholar]
 Lau, Y., & Finn, J. M. 1990, ApJ, 350, 672 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, R. P. 2006, Space Sci. Rev., 124, 233 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, R. P., Dennis, B. R., Hurford, G. J., et al. 2002, Sol. Phys., 210, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Litvinenko, Y. E. 1996, ApJ, 462, 997 [NASA ADS] [CrossRef] [Google Scholar]
 Litvinenko, Y. E. 2006, A&A, 452, 1069 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Litvinenko, Y. E., & Craig, I. J. D. 1999, Sol. Phys., 189, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Litvinenko, Y. E., & Somov, B. V. 1993, Sol. Phys., 146, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Litvinenko, Y. E., Forbes, T. G., & Priest, E. R. 1996, Sol. Phys., 167, 445 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, W., Berger, T. E., Title, A. M., Tarbell, T. D., & Low, B. C. 2011, ApJ, 728, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Longcope, D. W. 2005, Liv. Rev. Sol. Phys., 2, 7 [Google Scholar]
 Maclean, R. C., Parnell, C. E., & Galsgaard, K. 2009, Sol. Phys., 260, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Masuda, S., Kosugi, T., Hara, H., Tsuneta, S., & Ogawara, Y. 1994, Nature, 371, 495 [NASA ADS] [CrossRef] [Google Scholar]
 McKenzie, D. E. 2011, Phys. Plasmas, 18, 111205 [NASA ADS] [CrossRef] [Google Scholar]
 Parnell, C. E., Smith, J. M., Neukirch, T., & Priest, E. R. 1996, Phys. Plasmas, 3, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Pontin, D. I. 2011, Adv. Space Res., 47, 1508 [NASA ADS] [CrossRef] [Google Scholar]
 Pontin, D. I., Bhattacharjee, A., & Galsgaard, K. 2007a, Phys. Plasmas, 14, 052106 [NASA ADS] [CrossRef] [Google Scholar]
 Pontin, D. I., Bhattacharjee, A., & Galsgaard, K. 2007b, Phys. Plasmas, 14, 052109 [NASA ADS] [CrossRef] [Google Scholar]
 Pontin, D. I., AlHachami, A. K., & Galsgaard, K. 2011, A&A, 533, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Priest, E. R. 2000, in High Energy Solar Physics Workshop – Anticipating Hess!, ed. R. Ramaty, & N. Mandzhavidze, ASP Conf. Ser., 206, 13 [Google Scholar]
 Priest, E. R., & Pontin, D. I. 2009, Phys. Plasmas, 16, 122101 [NASA ADS] [CrossRef] [Google Scholar]
 Priest, E. R., & Titov, V. S. 1996, Roy. Soc. London Philos. Trans. Ser. A, 354, 2951 [NASA ADS] [CrossRef] [Google Scholar]
 Speiser, T. W. 1965, J. Geophys. Res., 70, 4219 [NASA ADS] [CrossRef] [Google Scholar]
 Török, T., Aulanier, G., Schmieder, B., Reeves, K. K., & Golub, L. 2009, ApJ, 704, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Vekstein, G. E., & Browning, P. K. 1997, Phys. Plasmas, 4, 2261 [NASA ADS] [CrossRef] [Google Scholar]
 Vilmer, N., MacKinnon, A. L., & Hurford, G. J. 2011, Space Sci. Rev., 159, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Watson, P. G., & Craig, I. J. D. 2002, Sol. Phys., 207, 337 [NASA ADS] [CrossRef] [Google Scholar]
 Wood, P., & Neukirch, T. 2005, Sol. Phys., 226, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Zharkova, V. V., & Gordovskyy, M. 2004, ApJ, 604, 884 [NASA ADS] [CrossRef] [Google Scholar]
 Zharkova, V. V., & Gordovskyy, M. 2005, MNRAS, 356, 1107 [NASA ADS] [CrossRef] [Google Scholar]
 Zharkova, V. V., Arzner, K., Benz, A. O., et al. 2011, Space Sci. Rev., 159, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., & Parks, G. 1993, J. Geophys. Res., 98, 7603 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1 a) Representative magnetic field lines for the spine model with parameters λ = 0.75, B_{s} = 3.4, α = −2, η = 3 × 10^{3}. The field lines are seeded from the top and base of the spine axis. b) Showing the direction and relative strength of the current in a plane of constant z for the same parameters. Here, is the size of the resistive region centred on the spine axis. 

Open with DEXTER  
In the text 
Fig. 2 Representative magnetic field lines for the fan solution with parameters λ = 0.75, B_{s} = α = 10, η = 10^{6}. The lines are again seeded from the top (solid lines) and base (dashed lines) of the spine. 

Open with DEXTER  
In the text 
Fig. 3 Angular distribution of protons from the null point for spine model at t = 1.6 × 10^{6} T_{ω,p}, at which the energy spectrum is steady state, for parameters λ = 0.75, B_{s} = 10, α = −0.1, η = 10^{6}. The initial distribution was Maxwellian at T = 86 eV in the upper right inflow region. 

Open with DEXTER  
In the text 
Fig. 4 Typical proton global trajectory from population “A” for parameters λ = 0.75, B_{s} = 10, α = −0.1, η = 10^{6}. The particle is taken from the many particle simulation, having initial position x = (−0.52,0.80,0.29) and velocity v = (−0.0044,0.0013, −0.0088)v_{Ae}. The magnetic field lines (thin dashed) are a projection of the field from the plane of the trajectory φ ≈ 120°. Inset shows the 3D trajectory of the proton as it crosses the spineaxis, the solid line in the centre is the line B_{z}(x,y,0.95) = 0. 

Open with DEXTER  
In the text 
Fig. 5 Typical proton trajectory from population “B” in the many particle simulation, with initial position x = (−0.54,0.78,0.31)L_{0} and velocity v = (−0.004, −0.006,0.002)v_{Ae}. The dashed lines show the projection of the magnetic field from the plane of motion, φ ≈ 110°, onto the y − z plane. Inset shows the motion in the x − y plane close to the spine axis. The purple arrows show the direction and relative magnitude of the gradient drift velocity and the dashdotted lines show contours of the electric potential, with the intersecting tick mark indicating lower potential to the right. 

Open with DEXTER  
In the text 
Fig. 6 Energy spectrum from the many particle simulation for protons in the spine model, with parameters λ = 0.75, B_{s} = 10, α = −0.1, η = 10^{6}. For particles leaving the R = 5 sphere, the energy at the time of crossing is used. 

Open with DEXTER  
In the text 
Fig. 7 a) Angular distributions of protons in fan model at t = 0, with initial temperature T = 86 eV. The xaxis is φ = 0 and the fan plane is β = 0. Protons are coloured by the final energy at t = 4000. Parameters used are λ = 0.75, B_{s} = α = 10, η = 10^{6}. b) Angular distributions at time t = 4000T_{ω,p} when the energy spectrum has reached steady state. 

Open with DEXTER  
In the text 
Fig. 8 Two typical proton trajectories from the many particle simulations in the fan model. Proton “1” is represented by a thin line with initial position (0.86,0.41,0.30), and proton “2” by a thick line with initial position (0.8,0.003,0.6). The parameters are λ = 0.75, B_{s} = α = 10, η = 10^{6}. The solid lines are representative magnetic field lines (seeded from the top of the spine axis and projected into the 2D planes) and the arrows show the direction and relative magnitude of the electric drift velocity. a) In the x − y plane, where the electric drift arrows are from the edge of the current sheet z = z_{η}. b) In the x − z plane close to the current sheet, where the electric drift arrows are plotted on y = 0. The initial positions are not shown in this plane. 

Open with DEXTER  
In the text 
Fig. 9 Energy spectrum for the many particle fan simulation for protons with parameters λ = 0.75, B_{s} = α = 10, η = 10^{6}. For particles leaving the R = 5 sphere, the energy at the time of crossing is used. 

Open with DEXTER  
In the text 
Fig. 10 Proton positions after initial distribution within the fan current sheet, such that  x_{0}  < 1,  y_{0}  = 0, . Parameters used are λ = 0.75, B_{s} = α = 5, η = 10^{8}, L_{0} = 10^{6} m. The dashed lines are representative magnetic fieldlines inside the current sheet (note the difference in scale of the zaxis). The solid black line is the line (x_{1},0,z_{1}) such that B_{x}(x_{1},0,z_{1}) = 0. The particles are stopped at R = 5 from the null point. 

Open with DEXTER  
In the text 
Fig. 11 Energy distribution of particles as they cross the R = 5 boundary, where initial position is within  x_{0}  < 1,  y_{0}  = 0, . Here, L_{0} = 10^{6} m and the results for different values of η and B_{s} are plotted. The solid points are protons from the three simulations and the thin lines show the sinφ relationship in Eq. (37). 

Open with DEXTER  
In the text 
Fig. 12 Proton positions after initial distribution within the fan current sheet such that  x_{0}  ,  y_{0}  < 1 and . Parameters used are λ = 0.75, B_{s} = α = 5, η = 10^{8}, L_{0} = 10^{6} m. Particles circled in black are those that start in y < 0 and cross z = z_{η} ≈ 9.6 × 10^{5} before t = 19 000 T_{ω,p}. 

Open with DEXTER  
In the text 
Fig. 13 Percentage of total particles (+) reaching current sheet at z = z_{η} from the external region (R = 1 in the inflow quadrant), and the mean time taken (*), for different values of η, B_{s} and λ. Each data point is from a many particle simulation with initial Maxwellian distribution (T = 86 eV) of 5000 protons. The setup is the same as that in Fig. 7. a) For different values of η with fixed B_{s} = α = 5, λ = 0.75. The solid line is a least squares fit to the points. b) For different values of B_{s} (with B_{s} = α) for fixed η = 10^{8}, λ = 0.75. The solid line is a least squares fit. c) For different λ with fixed B_{s} = α = 5, η = 10^{8}. No curve was fit. 

Open with DEXTER  
In the text 
Fig. 14 Scalings of steadystate energy spectra for global fan simulation. a) For different values of η with fixed B_{s} = α = 5 and λ = 0.75. The time taken to reach steady state was t = 8 × 10^{4}, t = 2 × 10^{4} and t = 6.4 × 10^{3} T_{ω,p} for η = 10^{10}, η = 10^{8} and η = 10^{6} respectively. b) For different values of B_{s} = α with fixed η = 10^{8} and λ = 0.75. Steadystate was reached at t = 2 × 10^{5}, t = 2 × 10^{4} and t = 5 × 10^{3} T_{ω,p} for B_{s} = 1, B_{s} = 5 and B_{s} = 30 respectively. c) For different values of λ with fixed B_{s} = α = 5 and η = 10^{8}. Time to steady state is t = 10^{5}, t = 2 × 10^{4} and t = 8 × 10^{3} T_{ω,p} for λ = 0.9, λ = 0.75 and λ = 0.3 respectively. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.