Free Access
Issue
A&A
Volume 542, June 2012
Article Number A47
Number of page(s) 14
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201218857
Published online 04 June 2012

© ESO, 2012

1. Introduction

Observations of Hard X-ray (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 non-thermal accelerated particles. The dominant HXR sources are chromospheric foot-points of the flaring loops, at which there is continuum free-free 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 neutron-capture and nuclear de-excitation (see e.g. Vilmer et al. 2011, for a review), suggesting ions with energies up to  ~100 MeV/nucleon.

When the emission from the foot-points 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 non-thermal 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 non-thermal coronal HXR source. Super-Dreicer (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 X-points (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 frozen-in 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 flux-pileup 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 time-dependent 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 spine-fan 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  ~  107 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, steady-state, 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, steady-state, and incompressible MHD equations. These are normalised in the usual way by the characteristic magnetic field strength B0, at a typical length scale L0 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 L0 boundary is half the plasma beta pe = β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 ~ 1012 − 1014.

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 Bs 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 z-axis 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 z-direction QS = Z(x,y). For the fan model, it distorts the spine in the x-direction (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 co-ordinates (r,φ,z) is (9)(Craig et al. 1997) in terms of the confluent hypergeometric (Kummer) function M(a,b,ζ) (Abramowitz & Stegun 1972). The flux pile-up factor Bs 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 frozen-in 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 Pr = αr/2 is the radial part of the potential field.

This electric field is curl-free (as required for steady-state) 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 x-direction 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 z-direction.

thumbnail Fig. 1

a) Representative magnetic field lines for the spine model with parameters λ = 0.75, Bs = 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 null-point 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 non-adiabatic 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 p0 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 p0 are negative, as α < 0, so constraints must be put on the values of α and Bs 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 E0 = 1500 V/m on the r = 1, φ = π/2 boundary (or normalising by suitable solar coronal values, vAe = 6.5  ×  106 ms-1, B0 = 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 Bs ~ η-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 Bs ~ η-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η) ~ Bs, giving magnetic pressure there from Eq. (17) . To avoid negative thermal pressure in the model this requires the null point pressure p0 > (Z(rη))2 ~ η-2 which is unphysically large for the values of η considered.

Craig et al. (1997) showed that Bs 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 ≈ p0 − P2, where P(1) ~ α. We must constrain α ≤ Bs or this dynamic pressure will require the gas pressure at the null to be even larger. The maximum value we can take for p0 is the largest possible external hydromagnetic pressure pe,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, Be,max = 0.3 T. This gives a normalised saturation value Bs,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 α ≈  − Bs 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  | α |  = Bs. 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 Bs = 10, α =  −0.1, this limits the reconnection exhaust close to the spine current sheet to sub-Alfvé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, Xmax ≈ Bs. 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 z-dependence; 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 spine-fan 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 L0, effectively limiting the size of the current sheet.

thumbnail Fig. 2

Representative magnetic field lines for the fan solution with parameters λ = 0.75, Bs = α = 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 Bs ~ η − 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 Bs,max = 30 and we have α ≤ Bs 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) Bs, 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 α = Bs (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 x2 + y2 = 1, z = zη to the local Alfvén speed.

Using the asymptotic approximation (14) we find that the z-component of the electric drift that brings the particles to the fan plane is, for x,z ≫ η1/2, (26)for the optimised solution α = Bs. This gives electric drift inflow for positive xz (as Pz < 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 Variable-Step Variable-Order Adam’s method, where the step size is recalculated to properly resolve gyro-motion, 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)mc2. For each simulation we find that this is conserved up to 5 significant figures. Also, to check the code handles non-adiabatic 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 L0, the normalising length scale, to be L0 = 104 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 L0 = 106 m for simulations where particles are initially distributed within the current sheet, as velocities are typically much faster here. Note that a change in L0 also changes the value of η as given in Eq. (1).

All magnetic fields mentioned in Sect. 2 have dimensions of B0 = 0.01 T, typical for the solar corona. We set v0 = vAe = 6.5  ×  106 ms-1 (corresponding to a number density n0 = 1.126  ×  1015 m-3). All dimensionless times quoted are in terms of the gyro-period, Tω = 2πm/(qB0).

To examine single particle trajectories in both models we choose values of Bs = 10, λ = 0.75, η = 10-6. This η value is rather large, towards the highest possible anomalous resistivities (with L0 = 104 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 = 106 K (86 eV) in the spine model fields. The protons have positions from a uniform random distribution at a global distance x2 + y2 + z2 = 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 x-axis and β = 0 is the fan plane).

Figure 3 shows the final spatial distribution and energies of the particles at time t = 1.6  ×  106Tω,p ≈ 10 s, at which the energy spectrum becomes approximately steady-state. 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 Bs = 10, η = 10-6, α =  −0.1. This value of α limits the bulk flow exhaust speed to be sub-Alfvé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.

thumbnail Fig. 3

Angular distribution of protons from the null point for spine model at t = 1.6  ×  106Tω,p, at which the energy spectrum is steady state, for parameters λ = 0.75, Bs = 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.

thumbnail Fig. 4

Typical proton global trajectory from population “A” for parameters λ = 0.75, Bs = 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)vAe. 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 spine-axis, the solid line in the centre is the line Bz(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 (x0,y0,z0) = (−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  ×  105Tω,p ≈ 2 s (inset). At this point the proton becomes unmagnetised as the gyro-radius becomes comparable to the length-scale of magnetic field gradient ρ/LB > 1 (we typically find that gyro-motion starts to break down when ρ/LB ≳ 10-2).

thumbnail Fig. 5

Typical proton trajectory from population “B” in the many particle simulation, with initial position x = (−0.54,0.78,0.31)L0 and velocity v = (−0.004, −0.006,0.002)vAe. 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 dash-dotted 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 x-direction, parallel to the current at the spine, with du/dt ≈ qE0/m as it crosses r = 0. For small displacements in the y-direction a strong Lorentz force due to the Bz field returns it to y = 0 line. These oscillations are Speiser-like (Speiser 1965) with frequency approximately ω ∝ t1/2.

The Speiser-like motion finishes and the first gyrations start (not shown) when the proton reaches r ≈ 5rη at which ρ/LB ≲ 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 x-direction (note the y-axis 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 Speiser-like oscillations following the Bz(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 Bz is the dominant component of the magnetic field.

After the proton becomes re-magnetised at r ≈ 5rη it has weak electric drift, vE ≪ 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 z-direction, 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 =  −5L0 which we do consistently throughout these simulations. At the time of stopping it is actually losing energy as it re-crosses the same electric potential contours. However, some other protons from the many-particle 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 = 5L0 spherical boundary we use the energy at the instant of crossing (if this is not done some protons reach order  ~ 102L0 which becomes unrealistic as the background field increases without bound away from the null, also causing the time-step 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 ~ qExacc where E ≈ E0 ≈ ηBs/rη   [vAeB0 ]  and xacc ≈ 2rη   [L0 ]  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  | α |  < Bs should not have a large effect on this result.

thumbnail Fig. 6

Energy spectrum from the many particle simulation for protons in the spine model, with parameters λ = 0.75, Bs = 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 Bs = α = 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 α < Bs) as the external electric drift (Eq. (26)) scales more favourably with the resistivity. Protons that cross the R = 5L0 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 anti-clockwise 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 ≈ 107 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 ρ/LB ~ 10-3, due to very fast electric drifts close to the sheet. Here the first adiabatic invariant, the constancy of μ, is also violated.

thumbnail Fig. 7

a) Angular distributions of protons in fan model at t = 0, with initial temperature T = 86 eV. The x-axis is φ = 0 and the fan plane is β = 0. Protons are coloured by the final energy at t = 4000. Parameters used are λ = 0.75, Bs = α = 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 x-axis 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 y-axis 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 Speiser-like motion. They become slightly unmagnetised, with maximum ρ/LB ≈ 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.

thumbnail 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, Bs = α = 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 (x0,y0,z0) = (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 x-axis. The first adiabatic invariant is not violated, μ = const. and the maximum ρ/LB ~ 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)L0. The azimuthal electric drift is weak close to the x-axis and the proton drifts down to the fan current sheet. It enters the sheet at (x,y) = (−0.02,0.15) and becomes unmagnetised: ρ/LB > 1 and μ is not conserved. We observe Speiser-like oscillations as the proton is accelerated in the y-direction. 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 = qE0t2/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 Speiser-like throughout as ρ/LB < 10-2 when the proton reaches y = 1.5L0. The proton reaches a global distance in the y-direction and becomes magnetised by the background By 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.

thumbnail Fig. 9

Energy spectrum for the many particle fan simulation for protons with parameters λ = 0.75, Bs = α = 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 ≈ 105 eV and ϵk ≈ 107.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 non-uniform 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 107.5 − 8 eV in Fig. 9 becomes a “bump on tail” centred at 108 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 By(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 x-axis, 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 non-adiabatic to adiabatic motion.

Firstly, we place particles within the sheet so that they are initially unmagnetised by the By(y) component of the background field. They are magnetised only by the strong Bx(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, Bs = α = 5, η = 10-8. We increase the dimensional box length to L0 = 106 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 L0 in Eq. (1). We again artificially stop the particles as they cross the R = 5 spherical surface.

thumbnail Fig. 10

Proton positions after initial distribution within the fan current sheet, such that  | x0 |  < 1,  | y0 |  = 0, . Parameters used are λ = 0.75, Bs = α = 5, η = 10-8, L0 = 106 m. The dashed lines are representative magnetic field-lines inside the current sheet (note the difference in scale of the z-axis). The solid black line is the line (x1,0,z1) such that Bx(x1,0,z1) = 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 Bx(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, Ez 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 vEz ≈  − (η/z)   [vAe ] . It becomes unmagnetised at the fan plane, z ≈ 0, close to the null point and is directly accelerated in the y-direction. We compare this trajectory to the analytical WKB solutions of Litvinenko (2006). The ejection time for a non-relativistic 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 non-adiabatic region and the displacement magnetic field gradient is much stronger than the gradient from the background component, Bs/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 y-direction during this time, (33)which we compare with size of the non-adiabatic region from the null in this direction. The particle begins to be re-magnetised 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 gyro-frequency of a particle gyrating around By(y). We use v(y) = (2qEyy/m)1/2 from direct acceleration (if we use v(y) = Ey/By the value for y ∗  differs by 21/3), assuming that there was no initial y-velocity 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 gyro-turning 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 By(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 Bx(x,z) = 0, at which they become unmagnetised and accelerated. We do appear to see some gyro-turning for protons starting at  | x |  ≈ 1. This is probably due to the strong component of the Lorentz force, vyBz, that acts to turn the trajectory to the x-direction. For particles starting at x = 0 the Bz magnetic field switches sign during the z-oscillations, but at x = 1 the proton is unmagnetised below the fan plane and Bz stays positive. The proton is turned in the positive x-direction 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 field-line with a parallel component of the electric field E|| = E·P/ | P | .

We artificially stop the protons at R = 5 L0 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 α = Bs, 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 Bs. This expression (thin line) fits the energies of simulated particles (circles) as they cross R = 5 L0 very well.

thumbnail Fig. 11

Energy distribution of particles as they cross the R = 5 boundary, where initial position is within  | x0 |  < 1,  | y0 |  = 0, . Here, L0 = 106 m and the results for different values of η and Bs 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

thumbnail Fig. 12

Proton positions after initial distribution within the fan current sheet such that  | x0 | , | y0 |  < 1 and . Parameters used are λ = 0.75, Bs = α = 5, η = 10-8, L0 = 106 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” By(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 ρ/LB in the range 10-4 − 10-2. They are not ejected due to gyro-turning in the sense of Speiser (1965) as this requires non-adiabatic motion. The trajectories in the y − z plane seem to follow the magnetic field lines closely, although they have strong electric drift from the Ez component of the electric field. Within the current sheet,  | z |  < zη the truncated power series in Eq. (13) gives the z-component 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 x-direction for protons in the upper half of the sheet 0 < z < zη, and in the positive x-direction for  − zη < z < 0. This electric field also contributes to vEy 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, Bs = α = 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 non-adiabatic.

thumbnail 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 η, Bs and λ. Each data point is from a many particle simulation with initial Maxwellian distribution (T = 86 eV) of 5000 protons. The set-up is the same as that in Fig. 7. a) For different values of η with fixed Bs = α = 5, λ = 0.75. The solid line is a least squares fit to the points. b) For different values of Bs (with Bs = α) for fixed η = 10-8, λ = 0.75. The solid line is a least squares fit. c) For different λ with fixed Bs = α = 5, η = 10-8. No curve was fit.

Open with DEXTER

thumbnail Fig. 14

Scalings of steady-state energy spectra for global fan simulation. a) For different values of η with fixed Bs = α = 5 and λ = 0.75. The time taken to reach steady state was t = 8  ×  104, t = 2  ×  104 and t = 6.4  ×  103Tω,p for η = 10-10, η = 10-8 and η = 10-6 respectively. b) For different values of Bs = α with fixed η = 10-8 and λ = 0.75. Steady-state was reached at t = 2  ×  105, t = 2  ×  104 and t = 5  ×  103Tω,p for Bs = 1, Bs = 5 and Bs = 30 respectively. c) For different values of λ with fixed Bs = α = 5 and η = 10-8. Time to steady state is t = 105, t = 2  ×  104 and t = 8  ×  103Tω,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 Bs. 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 Bs 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 steady-state. 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 steady-state, 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 = L0 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 steady-state 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 λ, Bs and η as free parameters (for the optimised solution α = Bs, 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 vE(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, Bs = α = 10 and λ = 0.75 we find the energy spectrum from a distribution of protons starting with thermal energy at R = L0 from the null becomes power law at steady state for particles accelerated due to fast non-uniform drifts (index  ≈  − 1.5 between 105 eV and 107.5 eV). The electric drift in the external region is much quicker than the spine model, vE ~ η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 re-magnetised 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 non-thermal 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 self-consistent 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 time-dependence and compressibility, by using electromagnetic fields from MHD simulations, to see how particles behave in so called spine-fan 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

All Figures

thumbnail Fig. 1

a) Representative magnetic field lines for the spine model with parameters λ = 0.75, Bs = 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
thumbnail Fig. 2

Representative magnetic field lines for the fan solution with parameters λ = 0.75, Bs = α = 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
thumbnail Fig. 3

Angular distribution of protons from the null point for spine model at t = 1.6  ×  106Tω,p, at which the energy spectrum is steady state, for parameters λ = 0.75, Bs = 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
thumbnail Fig. 4

Typical proton global trajectory from population “A” for parameters λ = 0.75, Bs = 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)vAe. 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 spine-axis, the solid line in the centre is the line Bz(x,y,0.95) = 0.

Open with DEXTER
In the text
thumbnail Fig. 5

Typical proton trajectory from population “B” in the many particle simulation, with initial position x = (−0.54,0.78,0.31)L0 and velocity v = (−0.004, −0.006,0.002)vAe. 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 dash-dotted lines show contours of the electric potential, with the intersecting tick mark indicating lower potential to the right.

Open with DEXTER
In the text
thumbnail Fig. 6

Energy spectrum from the many particle simulation for protons in the spine model, with parameters λ = 0.75, Bs = 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
thumbnail Fig. 7

a) Angular distributions of protons in fan model at t = 0, with initial temperature T = 86 eV. The x-axis is φ = 0 and the fan plane is β = 0. Protons are coloured by the final energy at t = 4000. Parameters used are λ = 0.75, Bs = α = 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
thumbnail 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, Bs = α = 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
thumbnail Fig. 9

Energy spectrum for the many particle fan simulation for protons with parameters λ = 0.75, Bs = α = 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
thumbnail Fig. 10

Proton positions after initial distribution within the fan current sheet, such that  | x0 |  < 1,  | y0 |  = 0, . Parameters used are λ = 0.75, Bs = α = 5, η = 10-8, L0 = 106 m. The dashed lines are representative magnetic field-lines inside the current sheet (note the difference in scale of the z-axis). The solid black line is the line (x1,0,z1) such that Bx(x1,0,z1) = 0. The particles are stopped at R = 5 from the null point.

Open with DEXTER
In the text
thumbnail Fig. 11

Energy distribution of particles as they cross the R = 5 boundary, where initial position is within  | x0 |  < 1,  | y0 |  = 0, . Here, L0 = 106 m and the results for different values of η and Bs 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
thumbnail Fig. 12

Proton positions after initial distribution within the fan current sheet such that  | x0 | , | y0 |  < 1 and . Parameters used are λ = 0.75, Bs = α = 5, η = 10-8, L0 = 106 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
thumbnail 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 η, Bs and λ. Each data point is from a many particle simulation with initial Maxwellian distribution (T = 86 eV) of 5000 protons. The set-up is the same as that in Fig. 7. a) For different values of η with fixed Bs = α = 5, λ = 0.75. The solid line is a least squares fit to the points. b) For different values of Bs (with Bs = α) for fixed η = 10-8, λ = 0.75. The solid line is a least squares fit. c) For different λ with fixed Bs = α = 5, η = 10-8. No curve was fit.

Open with DEXTER
In the text
thumbnail Fig. 14

Scalings of steady-state energy spectra for global fan simulation. a) For different values of η with fixed Bs = α = 5 and λ = 0.75. The time taken to reach steady state was t = 8  ×  104, t = 2  ×  104 and t = 6.4  ×  103Tω,p for η = 10-10, η = 10-8 and η = 10-6 respectively. b) For different values of Bs = α with fixed η = 10-8 and λ = 0.75. Steady-state was reached at t = 2  ×  105, t = 2  ×  104 and t = 5  ×  103Tω,p for Bs = 1, Bs = 5 and Bs = 30 respectively. c) For different values of λ with fixed Bs = α = 5 and η = 10-8. Time to steady state is t = 105, t = 2  ×  104 and t = 8  ×  103Tω,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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.