EDP Sciences
Open Access
Issue
A&A
Volume 561, January 2014
Article Number A107
Number of page(s) 13
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201322199
Published online 16 January 2014

© ESO, 2014

1. Introduction

A solar flare is believed to occur via the release of energy stored in the Sun’s magnetic field when a magnetic reconnection event occurs in the solar corona. Some of this energy is released in the form of accelerated particles that transport energy throughout the solar atmosphere and produce X- and gamma-radiation (Aschwanden 2002). The exact mechanism of particle acceleration is still debated but some particles will certainly be accelerated by the electric field in the reconnection region. Martens (1988) suggested that electric field acceleration in a collisionless current sheet could release the majority of the flare energy in the form of fast particles, primarily ions (see Speiser 1965). Since then many authors have studied such a situation by calculating test particle trajectories in model electromagnetic fields (as a random sample: Petkaki & MacKinnon 1997; Wood & Neukirch 2005; Hamilton et al. 2005; Hannah & Fletcher 2006; Browning & Gordovskyy 2010; Petkaki & MacKinnon 2011; Burge et al. 2012; see also the review of Zharkova et al. 2011).

An unexpectedly bright, coronal hard X-ray (HXR) source was first discovered using the HXT instrument on Yohkoh (Masuda et al. 1994) and several more examples have been found in RHESSI data (Krucker et al. 2008a). These sources have attracted great interest particularly because they may represent the location of the electron acceleration region (Krucker et al. 2008a). In at least one case the number of accelerated electrons appears to be comparable to the total present in the emitting volume, adding weight to such an interpretation (Krucker et al. 2010). Fletcher & Martens (1998) point out that electrons accelerated in a reconnection region, where very low magnetic fields would be expected, would mostly be trapped in the corona by magnetic mirroring (see also Petkaki 1996). A localised, coronal hard X-ray source would result.

Masuda et al. (1994) suggested that a coronal HXR source could be created by a very high-temperature plasma at the top of a flaring loop. However, as Fletcher (1995) pointed out, a HXR source created by heating should be seen to increase in size as the plasma expands. That this is not seen would then require some kind of plasma confinement at the top of the loop. If instead the loop top source results from non-thermal particles contained by transport effects, no source of thermal emission is required.

Direct, imaging localisation of the acceleration region would be an exciting step forward, opening up new possibilities for probing the physical processes that dominate there. It also places new demands on modelling, however. In particular, we know that electrons must undergo collisions in order to emit bremsstrahlung HXRs. Therefore self-consistent models for coronal HXR sources must include collisional energy loss and scattering alongside the accelerating agent. Most test particle modelling is collisionless, but HXR radiation would not be detectable if collisions were truly negligible.

Hamilton et al. (2003) developed a method of including collisions when following particle trajectories, but their method includes only energy change, not pitch angle scattering. Accompanying energy loss and occurring on the same timescale (e.g. Spitzer 1962; Trubnikov 1965), pitch angle scattering via collisions must be expected to change the outcomes for test particles. It will lead at least some electrons to re-encounter the dissipation region more often, possibly offsetting collisional energy losses. The idea that scattering in angle, even if accompanied by energy loss, can lead to enhanced acceleration was suggested in a different, ionospheric context by Gurevich et al. (1985). Any such scattering effects here would be compounded with mirroring, because particles encounter greater magnetic field strengths away from the null. Thus the net effect of collisions on particle acceleration needs more detailed consideration.

Gordovskyy et al. (2012) studied test particle evolution in a detailed, MHD model of a twisted lop undergoing reconnection, but the guiding centre approximation was appropriate in their work. Here the chaotic character of particle orbits near the null point (Martin 1986; Chen 1992) rules out a guiding centre description, as well as implying that small collisional effects may have substantial consequences for particle evolution.

In this paper we describe work aimed at including the effects of binary collisions in reconnection test particle calculations. In Sect. 2 we write down the model equations including both Lorentz force terms from large-scale electromagnetic fields, and the stochastic (Wiener process) terms that describe collisional energy loss and scattering. We also briefly describe the simple electromagnetic fields we invoke to model a reconnection region. In Sect. 3 we describe the numerical method used to follow test particle trajectories, based on Honeycutt’s (1992) extension of the RK4 method to stochastic differential equations. In an Appendix we test this method on the 1D problem described by MacKinnon & Craig (1991), thereby verifying that it reproduces analytical results for the distribution function, at least as well as simpler numerical methods. Section 4 looks at the changes in single particle behaviour resulting from collisions, while Sect. 5 studies the distributions that emerge. In Sect. 4 we study the modifications to electron acceleration near null points that result from binary collisions. Section 7 discusses the implications of our results, and in particular their possible application interpretation of coronal HXR sources.

2. Model for test particle evolution

2.1. Governing equations

Test particle calculations start from the equations of motion in prescribed electric and magnetic fields: \noindentHere E and B are the electric and magnetic field vectors, v is the test particle velocity, q and m are the particle charge and mass respectively, and c is the speed of light. We use c.g.s. units. External considerations lead to forms of E and B that describe the (comparatively) large-scale fields met in the reconnection region. As in previous work (Petkaki & MacKinnon 1997, 2007, 2011) we consider a 2D system with translational invariance in the third (z) direction, and concentrate on particles near a 2D X-type neutral point, where strong electric fields can develop and particles may be freely accelerated. As noted in the introduction, the presence of a null rules out a guiding centre description. To focus attention on effects due to collisions we make the simplest possible assumption, that steady reconnection results in a static, uniform electric field. Then we have where hats denote unit vectors, and d is the size of our acceleration region.

Along with the macroscopic fields represented by E and B, electrons feel a continually fluctuating electric field due to all the other discrete charges in the plasma. These result in both systematic slowing down and diffusion in velocity space (Chandrasekhar 1943; Trubnikov 1965). The evolution of the particle distribution function f(r,v,t) is then governed by a Fokker-Planck equation with drift and diffusion matrices A and D: (5)\noindentwhere we use the summation convention for indices. The calculation of A and D for test particles undergoing binary collisions is detailed in Trubnikov (1965) or Rosenbluth et al. (1957).

The distribution function f governed by Eq. (5) gives the envelope of the trajectories followed by test particles, each of which evolves according to a Langevin equation: (6)Here W(t) is a vector of independent Wiener processes, i.e. random quantities with mean 0 and variance equal to t. Also D1/2 denotes the matrix whose square is equal to D: D1/2.D1/2 = D. As long as D is diagonalizable, D1/2 may clearly be constructed thus: (7)where V is a matrix whose columns are the unit eigenvectors of D, and V-1 is its inverse. Q is the diagonal matrix whose diagonal values are the square roots of the eigenvalues of D. Below, however, we adopt a coordinate system in which D is diagonal so we only need take the square roots of the individual diffusion coefficients to construct D1/2.

The two descriptions (5) and (6) are equivalent (e.g. Gardiner 1985; MacKinnon & Craig 1991): the distribution function f may be generated via solution of (5), or by solving the Eqs. (6) and (2) for very many test particles. Electron collisional evolution is dominated by the myriad slightly fluctuating electrostatic fields from the most distant particles within the Debye sphere. Thus we may include the effects of collisions in test particle calculations by straightforwardly combining the right hand sides of (1) and (6): (8)

2.2. Drift and diffusion coefficients; coordinate system

Above we wrote the test particle governing equations quite generally, with only the fields (4) specified in a particular coordinate system. In what follows we assume that test particle energies are ≫thermal energies (although this may not necessarily be the case in the corona). Then diffusion in speed is negligible and we have (Spitzer 1962; Trubnikov 1965; Emslie 1978) (9)where Λ is the Coulomb logarithm, usually taken to be 25 in the solar corona. There are significant simplifications if we write v in spherical polar coordinates, using specifically (v,φ,θ) where and the angle θ is measured in the x-y plane, anticlockwise from . Then the diagonal elements of D become \noindentand all off-diagonal elements are identically 0.

The equations of motion are made dimensionless in the same manner as the equations of motion in Petkaki & MacKinnon (1997, 2007). Distances are normalised to a scale de and times to te, the electron gyroperiod at distance de from the null). Since we are interested in >keV energy electrons, speeds are normalised to c (and thus particle kinetic energies to rest mass energies). We avoid the appearance of any further dimensionless number in the problem if distances are normalised to , where B0 is the magnetic field strength gradient, a constant across the whole system. de is the radial distance at which the local electron gyroradius equals the light travel time to the null. If we take B = 100G and L to be a typical coronal length scale of 109cm then de = 1.3 × 105cm. The normalising timescale te then follows from these two quantities: te = de/c = 4.33 × 10-6 s. The electric field and the magnetic field are both normalised to B0de. Throughout the rest of this paper we use dimensionless units unless explicitly stated otherwise, and assume the values of de and te quoted here in converting results to dimensional units.

The full, dimensionless equations of motion (1) and (2) now become where Wφ and Wθ are Wiener processes (given subscripts φ and θ only to emphasise their independence). The dimensionless parameter controls the importance of collisions and is given in a fully ionised hydrogen plasma by (19)where ne is the electron number density, Λ is the Coulomb logarithm and te is our normalising time for electrons. Taking Λ = 25, where n10 = ne/1010 and B07 = 107B0. The smallness of does not mean that collisions are unimportant. First, since the (non-adiabatic) particle orbits near the null are chaotic (Martin 1986; Chen 1992), extra terms that are small in magnitude may nonetheless have a strong influence on test particle outcomes. Second, in the 2D configuration studied here particles can easily be trapped near the null for times comparable to their total collisional lifetimes so collisions become cumulatively important. We return to these points below.

2.3. Electric field

A variety of observational techniques have yielded estimates of peak electric field strengths in flares and erupting prominences (Somov et al. 2008; Foukal et al. 1987; Qiu et al. 2002). These have ranged from a few tens of V/m to more than 1 kV/m, i.e. roughly 0.005−0.3 in our dimensionless units (assuming the typical coronal parameters mentioned above).

A check on our adopted values is provided by estimating the associated inflow velocity vin, In what follows we will mostly adopt E = 0.001 or smaller values (equivalent to 3.9 Vm-1). With a 100 G strength field on the boundary of our system, this implies an inflow velocity of 0.39 km s-1, slightly lower than the values observed in flares which are typically on the order of a few kilometers per second (see e.g. Yokoyama et al. 2001).

3. Numerical integration method

We wish to generate a numerical approximation to the solution of the stochastic differential Eqs. (13)−(18). In the absence of the stochastic, collisional terms, we have the same sort of test particle problem already discussed by previous authors using a variety of numerical methods to integrate the equations of motion. Here we prefer the fourth-order Runge-Kutta (RK4) numerical method (e.g. Press et al. 1992). With a numerical timestep of Δt, RK4 generates a numerical estimate of the solution to a precision of (Δt)5. We see below that our RK4 integrations yield distributions with the same properties as more sophisticated methods (Hannah & Fletcher 2006; Petkaki & MacKinnon 2011) − even if individual particle trajectories are less precisely described. For our purposes RK4 is the obvious choice because there is a rigorously founded extension which deals with stochastic force terms (Honeycutt 1992). As is well known, trajectories of particles subject to stochastic forces, e.g. Brownian motion particles, are continuous everywhere but differentiable nowhere, so that methods for numerical integration must be reconsidered carefully. Stochastic differential equations are often solved numerically using the Euler method or variants (MacKinnon & Craig 1991; Higham 2001). These methods would give a poor description of the Lorentz force terms in (13)−(18), however.

Honeycutt (1992) developed an extension of Runge-Kutta methods to stochastic differential equations. We adopt this extended RK4 method here to retain a suitably precise description of gyro-motion while including stochastic effects of binary collisions. Then we can compare results with and without collisions, confident that the comparison has not been contaminated by artefacts resulting from two different numerical methods.

We used a fixed step size, selected by calculating the energy distribution of electrons at a fixed time for a variety of step sizes, and noting the largest step size below which the electron energy distribution does not change. Since we were able to find such a timestep a shadowing theorem presumably applies. After several initial experiments of this sort we mostly adopted a step size Δt = 0.03, probably too large for useful calculation of individual particle outcomes, but yielding robust distributions with acceptable running times.

For illustration, consider the 1D equation (20)With a timestep Δt the RK4 method generates a solution to (20) via repeated use of the algorithm (21)where Now, following Honeycutt (1992), consider the one-variable additive noise equation (22)which adds a stochastic forcing term to (20). The factor of 2 is included so that the corresponding diffusion term in the associated Fokker-Planck equation is just . Honeycutt shows that a numerically generated solution will reproduce the 4th order accuracy of RK4 and reproduce the statistical properties of the exact solution to the corresponding order in the Taylor expansion if it is calculated using (21) but with extra Wiener type terms added to each of the iterates for the next timestep: (23)where Note that precise reproduction of the statistical properties of the solution requires all of the numerical noise terms to have the same form; the terms in F2 and F3 may no longer be thought of as estimates of the solution at half a timestep, Δt/2, as they are in conventional RK4.

3.1. The test problem

In order to develop and test the stochastic RK4 algorithm, we consider a problem which already has a known solution. The problem used was the problem considered in MacKinnon & Craig (1991), which dealt with pitch-angle scattering of particles in a non-magnetised medium. The spatially integrated FP equation for this problem also has a known analytical solution (for the spatially homogeneous case) which is given in terms of the Legendre polynomials. Stochastic RK4 treatment of many test particles yields a distribution in very good agreement with this analytical solution, discussed in more detail in the appendix.

4. Effect of collisions on single particle behaviour

4.1. Magnitude of collisional terms

From (16), the distance over which collisions will make a significant difference to particle energy and direction is of order (24)whereas the size of the non-adiabatic region around the null within which particles can be energised is of order (Petkaki & MacKinnon 2007) (25)Collisional energy loss will have an important influence on the energy gained in the reconnection region if i.e. (26)\noindentEquation (26) will be satisfied only for electron energies well below even coronal thermal energies, depending only very weakly on ne and B0. However we see next that even small levels of pitch-angle scattering can make a big difference to particle evolution.

4.2. Single particle trajectories

Figure 1 shows two trajectories followed by electrons near the null. Both particles start with identical initial conditions (see figure caption) and encounter the same electromagnetic fields (E = 0.001 in our dimensionless units, equivalent to 3.9 V/m; 2D X-type B as above, ) but one trajectory is for a vacuum while the other experiences collisions in a medium such that K = (equivalent to an ambient density of 1010 cm-3 with the typical coronal parameters mentioned above). These trajectories have not been chosen randomly but to illustrate how pronounced the consequences of collisions can be.

thumbnail Fig. 1

Trajectory of one electron, integrated with and without the addition of pitch angle scattering. Both electrons had identical starting conditions. All electrons start with an energy of 1.3 keV. The electric field has magnitude 0.001 in our dimensionless units, which is 3.9 V/m for our typical, coronal parameters. The particle started at x = −0.079, y = 0.046 and z = 0, and with angles φ = 1.57 and θ = 2.61.

Open with DEXTER

The general character of particle orbits in this configuration is well understood (e.g. Bulanov & Sasorov 1976; Burkhart et al. 1990). Near the null, particles are not closely tied to field lines and their trajectories become chaotic (Martin 1986; Chen 1992). Beyond a distance given roughly by rAD above, their gyroradii become smaller than the field strength variation length scale and they behave in the usual adiabatic way, conserving magnetic moment ,showing curvature, gradient and E × B drifts, etc. In this 2D configuration E is perpendicular to B everywhere so particle energy gains can only be significant in the non-adiabatic region. The energy gained by a particle depends primarily on the length of time it spends in the dissipation region, which in turn will depend on its initial position (as is seen in the approximate, analytical discussion of Bulanov & Sasorov 1976).

The trajectory shown in black in Fig. 1 is for an electron in vacuum. It starts close to the null, where it can gain energy, moving away until it starts to move adiabatically. It mirrors at a fairly small value of r, returns to cross the vicinity of the null and gain energy again, and finally moves away beyond the boundary of the figure. The pitch angle it attains as it starts to move adiabatically is a consequence of its chaotic motion nearer the null, not predictable in any simple way except by this sort of numerical calculation. In this instance it finds itself moving with a small pitch angle so it is able to travel far from the null and does not mirror again within the time of the simulation (equivalent to 0.015 s, i.e. its stopping time in the absence of an electric field).

The trajectory shown in red is for the electron experiencing collisional scattering. Initially it closely resembles the vacuum case but slight differences have more and more pronounced consequences as the integration proceeds. Its movement as it travels back towards the null is more complex, it spends more time in the region where it can gain energy and it subsequently attains a large pitch angle and so moves adiabatically at a smaller distance from the null. For the remainder of the integration it bounces left to right along field lines, at the same time E × B drifting upwards, in a manner often seen for such test particles (e.g. Hannah & Fletcher 2006) but with added irregularities that result from continual, small, stochastic changes in pitch angle. Its energy is close to constant during this final period. The outcome for the particle is completely different from the vacuum case because of the sensitivity of the orbit near the null to even slight changes in velocity. In this particular case the particle experiencing collisions actually gains more energy than the vacuum case.

Table 1

Summary of results of simulations, showing decrease in the number of particles left in the simulation at t = 0.015 s for decreasing super-Dreicer fields.

4.3. The Dreicer field

Resulting from the velocity-dependence of collisional scattering, the idea of the critical field (Dreicer 1959) is also important for understanding what happens to individual particles. If electrons have speed less than the thermal velocity, collisions happen with almost constant frequency, increasing in number as the thermal velocity is approached. If an electron is moving faster than the thermal velocity, the collision frequency scales as 1/v2 (cf. Eq. (16)), so collisions become less frequent as the electron’s speed increases (see e.g. Rozelot et al. 2000; Trubnikov 1965).

Since electron energy loss rate decreases with energy, there is a critical electron energy for which energy gain from electric field is greater than energy loss from collisions. Electrons above this critical energy can be freely accelerated out of the thermal distribution by the electric field. The Dreicer field is the strength of electric field for which this critical energy equals the thermal energy, i.e. all electrons in the plasma can be freely accelerated. The speed at which collisions become less important as known as the runaway speed and is given by (e.g. Holman 1985): (27)vTe is the thermal speed of the electrons, which is given by: (28)The electric field strength ED is the Dreicer field, given in statvolt cm-1 by (e.g. Holman 1985) : (29)where Λ is the Coulomb logarithm, λD is the Debye length, and T is the plasma temperature. For the plasma we will consider (T = 1.4 × 107 K, n = 1010 cm-3,Λ = 25, the Dreicer field is 1.8 × 10-7 statvolt cm-1, which is 5.4 × 10-3 V/m. The electric field we apply in our simulations is generally 10-3 in our dimensionless units, equivalent to 3.9 Vm-1 for the typical coronal parameters mentioned above.

For the plasma we will consider, vTe is 1.5 × 109 cm s-1. This gives a runaway speed of 5.6 × 107 cm s-1 for E = 0.001, which in our units is a speed of 1.9 × 10-3. This means that all of the electrons in our distribution are “runaway” electrons, and can be accelerated out of a thermal distribution.

5. Effect of collisions on accelerated electron distributions

5.1. Reference case, = 3.2 × 10-8

We now concentrate on the effects of collisions on electron energy distributions, starting by comparing the distributions that result in the presence of an electric field E = 0.001 in vacuum, and in a medium with i.e., with B0 = 10-7, a density of 1010 cm-3, reasonable for the corona. Figure 2 presents the distributions assembled from the final states of 104 test electrons, followed up to a time equivalent to 0.015 s, their collisional stopping time in a medium of this density. All particles start with an energy of 1.3 keV, randomly chosen directions, and positions distributed randomly within r < 1 . For these simulations, a particle is considered to have stopped if its energy is less than 5.11 × 10-2 eV. This value was chosen as the simulation was found to become unstable if the particle energy fell below 10-7 in dimensionless units; in any case this is far below coronal thermal energies. 8298 of the original 10 000 electrons stayed above this energy until the end of the run and their distribution is shown in Fig. 2. The number of surviving electrons at the end of the simulation is always quoted in the box at upper left of the energy distribution figures.

The energy distribution in the vacuum case reproduces results seen in Petkaki & MacKinnon (2007), who used the more sophisticated Bulirsch-Stoer integration method, reassuring us that our RK4 results correctly describe particle distributions. The code used is also capable of reproducing the results of other X-type neutral point simulations seen in e.g. Hannah & Fletcher (2006).

thumbnail Fig. 2

Energies at t = 0.015 s for 10 000 electrons whose trajectories have been integrated in a vacuum (black), and in a medium with ne = 1010 cm-3 (red). All electrons start with an energy of 1.3 keV. The electric field has magnitude 0.001 in our dimensionless units, which is 3.9 V/m for our typical, coronal parameters.

Open with DEXTER

In the vacuum case, particles are accelerated by the electric field but to varying degrees depending on their initial conditions. In the denser medium acceleration by the electric field now acts in competition with collisional slowing down. Particles that would have been less accelerated anyway by the electric field are less able to withstand collisional energy loss. The maximum energy achieved is the same in both cases but more particles achieve this energy in the absence of collisions. The dense medium distribution develops a long tail to low energies and 17% of the particles stop completely. Of course, in the absence of the electric field all of them would have stopped by the end of the calculation.

We also investigated the energy distributions for electrons if E = 0.0001 and E = 10-5 (two smaller, but still super-Dreicer fields). As the electric field decreases, lower energies are achieved, both with and without collisions. However, once again, the electrons undergoing collisions should have lost all of their energy in this time. We do not plot the distributions for E = 0.0001 and E = 10-5, but the results are summarised in Table 1.

It can be seen that the fraction of particles left after one stopping time decreases with decreasing (super-Dreicer) electric field. What happens if a sub-Dreicer field is applied to the electrons? The Dreicer field for an electron density of 1010 cm-3 and temperature 1.4 × 107 K is 5.4 × 10-3 V/m. We imposed an electric field E = 10-7, equivalent to 3.9 × 10-4 V/m. Electrons in a field of this magnitude have a runaway speed of 5.7 × 107 m s-1, which is 0.18 in dimensionless units, meaning our electrons are initially travelling below the runaway speed, and collisions will be more important. The effect of the sub-Dreicer field on electrons which both do and do not experience collisions can be seen in Fig. 3.

thumbnail Fig. 3

Energies at t = 0.015 s for 10 000 electrons whose trajectories have been integrated in a vacuum (black) and with collisions (red; with ). All electrons start with an energy of 1.3 keV. The electric field has magnitude 10-7 in our dimensionless units, which is 0.00039 V/m for our typical, coronal parameters.

Open with DEXTER

In the absence of collisions, the energy of electrons in such a low field effectively does not change over the time of the integration. Once again, 80% of the particles experiencing collisions do not stop completely, even in a stopping time, although their energies are not interesting from the point of view of X-radiation and a complete description would need to take account of diffusion in energy close to the ambient thermal speed (cf. Galloway et al. 2005). In fact, fewer particles stop completely with this sub-Dreicer field after the expected stopping time than are found for a small but super-Dreicer field. Again they display a spread in energy reflecting the variety of histories near the null. A similar spread of energies actually occurs in the vacuum case but is not discernible compared to the particles’ initial energies.

thumbnail Fig. 4

Energies at t = 0.015 s for 10 000 electrons whose trajectories have been integrated considering only collisional energy loss (black) and with pitch angle scattering included (red). All electrons start with an energy of 1.3 keV. The dimensionless electric field has magnitude 10-3.

Open with DEXTER

thumbnail Fig. 5

Energies at t = 0.015 s for 10 000 electrons whose trajectories have been integrated considering only collisional energy loss (black) and with pitch angle scattering included (red). All electrons start with an energy of 1.3 keV. The dimensionless electric field has magnitude 10-7.

Open with DEXTER

5.2. Role of pitch-angle scattering

In Sect. 1 we speculated that pitch-angle scattering might help some particles to retain high energies for longer, via more frequent returns to the dissipation region. We studied this possibility by artificially suppressing the collisional terms in the angular rates of change, Eqs. (17) and (18), while retaining a nonzero () energy loss term in (16). Figures 4 and 5 show the energy distributions at t = 0.015 s for electrons which have undergone pitch angle scattering and collisional energy loss, as well as energy loss unaccompanied by scattering, for two values of the electric field. The results of these simulations are summarised in Table 2.

thumbnail Fig. 6

Distance from the neutral point (r) for each particle at each saved timestep, considering only collisional energy loss (top) and with pitch angle scattering included (bottom). The electric field has magnitude 0.001 in our dimensionless units, which is 3.9 V/m for electrons. This simulation was done for a sample of 100 electrons, and their positions were recorded every 100 timesteps.

Open with DEXTER

Table 2

Summary of results of simulations.

We see that pitch-angle scattering does indeed have the effect of allowing some particles to spend more time near the null and thus gain higher energies − in spite of collisional energy loss. These particles are few in number, however, only about 1% of the total, so this does not appear to be a strong effect, at least for the parameter range of these initial simulations.

thumbnail Fig. 7

Energies at t = 0.003 s for 10 000 electrons which have been moving in a plasma of density ne = 5 × 1010 cm-3, for a variety of electric fields. All electrons start with an energy of 1.3 keV. The initial positions of the particles were randomly distributed within the region r < 1, and the initial particle pitch angles were randomly chosen to be between 0 and 2π. For this simulation, .

Open with DEXTER

Figure 6 investigates this question further by showing the distance (r) of 100 electrons from the central null, whose positions were recorded every 100 timesteps. We show the example where E = 0.001. It can be seen that adding collisions causes particles to travel much less far from the neutral point (the maximum values of r is approximately halved), and that more timesteps (approximately 10% more) are spent at very small r (r < 1).

5.3. Effects of varying density

We further investigate the combined effects of collisions and acceleration near the null by following sets of test particles for a range of densities and electric field strengths. In the reference case of Sect. 5.1, the plasma density was set to ne = 1 × 1010 cm-3 (i.e. ). We now simulate particle behaviour for a variety of different densities: ne = 5 × 1010,1011,5 × 1111 and 1012 cm-3. Four different electric field strengths were used: E = 10-3,10-4,10-5 and 10-7, spanning both sub- and super-Dreicer fields. The results can be seen in Figs. 7 to 10. Particles were followed up until their stopping times so we grouped them by and each set of four figures is drawn for the same physical time.

thumbnail Fig. 8

Energies at t = 0.0015 s for 10 000 electrons which have been moving in a plasma of density ne = 1011 cm-3, for a variety of electric fields. All electrons start with an energy of 1.3 keV. The initial positions of the particles were randomly distributed within the region r < 1, and the initial particle pitch angles were randomly chosen to be between 0 and 2π. For this simulation, .

Open with DEXTER

thumbnail Fig. 9

Energies at t = 0.0003 s for 10 000 electrons which have been moving in a plasma of density ne = 5 × 1011 cm-3, for a variety of electric fields. All electrons start with an energy of 1.3 keV. The initial positions of the particles were randomly distributed within the region r < 1, and the initial particle pitch angles were randomly chosen to be between 0 and 2π. For this simulation, .

Open with DEXTER

thumbnail Fig. 10

Energies at t = 0.00015 s for 10 000 electrons which have been moving in a plasma of density ne = 1012 cm-3, for a variety of electric fields. All electrons start with an energy of 1.3 keV. The initial positions of the particles were randomly distributed within the region r < 1, and the initial particle pitch angles were randomly chosen to be between 0 and 2π. For this simulation, .

Open with DEXTER

Clear trends, mostly unsurprising, are evident. For fixed electric field strength E, increasing density causes lower particle energies to be attained. The spread of energies also decreases. We can also see a trend in the number of particles which remain energised at the end of the simulation. This trend is plotted in Fig. 11, which shows that there appears to be an optimum density at which more particles are able to remain energised at the end of the simulation. This could be due to the fact that we see an increased energy loss rate for denser plasmas, but denser plasmas will also scatter particles more in angle leading them to spend more time near the null where they may be energised, either via mirroring or via greater complexity of orbit, as seen in Fig. 1, causing them to return to the region of the magnetic field where they can be energised by the electric field. The competition between these two effects means that there is an optimum density at which the energy lost by particles is less than the energy they gain from the electric field. The value of this optimum density changes depending upon the electric field strength.

For sub-Dreicer fields, the number of particles which remain energised simply increases until all particles remain energised at the end of the simulation. This could be due to the fact that for sub-Dreicer fields particles collisions are more important; therefore as collisions (which may send particles back to regions where they can gain energy) increase in frequency, the number of particles which remain in the simulation also increases. Recall that the pitch-angle scattering rate in (17) and (18) increases as particle energy decreases.

5.4. Long times

Can the imposed electric field maintain particles at high energies for times much longer than a stopping time (as might help with interpreting coronal HXR sources Krucker et al. 2008b)? To test this, we ran a simulation with ne = 1010 cm-3, E = 0.001 (i.e. 1.8 V/m), for a running time equivalent to 0.15 s, which is 10 times the collisional stopping time for this density. The results can be seen in Fig. 12. Figure 12 shows that particles are unable to achieve an energy greater than around 50 keV for this electric field strength and simulation time. This is because particles are also losing energy due to collisions. In the ne = 1010 cm-3, case, about 80% of the particles stop altogether, and the rest form a population which have energies of around a few tens of kev. The particles in the less dense plasma do not achieve higher energies, but there are more of them (almost all of the initial population), and there is also a low energy tail.

Crucially, for the higher density, we see a narrower distribution of “surviving” particles with energies of a few tens of keV, energies in the HXR range (see e.g. Holman et al. 2003). We also attempted to run the simulation for longer times for E = 0.001 for chromospheric densities (ne = 1012 cm-3), but no particles survived for significantly longer than their collisional stopping times.

6. Effect of guide field

In all simulations previously discussed in this paper, we have considered a magnetic field configuration with Bz = 0. We now study the effects of collisions when the magnetic field has a small z-component, Bz. We take Bz = 1 in dimensionless units. This value of Bz was chosen so that it would be important for particle trajectories near the null but negligible at the system boundary. Such a guide field may keep particles near the dissipation region for longer, enabling them to gain higher energies (Litvinenko 1996; Hamilton et al. 2005).

The equations of motion in the presence of nonzero Bz are unchanged, with the exception of the rate of change of θ, whose deterministic part now becomes (30)We repeated our previous variation of system parameters (E and ne) and explored the same parameter space as we did for Bz = 0. A sample of the results can be seen in Figs. 13 and 14. In Fig. 13 we see that adding a guide field leads to more efficient acceleration and narrower final energy distributions, with no low energy “tail”. In the presence of a guide field, electron distributions with and without collisions look much more similar than they do when Bz = 0.

thumbnail Fig. 11

Variation in number of particles remaining in the X-type neutral point simulations discussed in Sects. 4 and 5. Measurements were taken after one stopping time for a variety of different values of electric field, E.

Open with DEXTER

thumbnail Fig. 12

Final energies of particles at t = 0.15 s for E = 0.001 and two different densities. The simulation time was ten times the stopping time for ne = 1010 cm-3, and equivalent to the stopping time for ne = 109 cm-3.

Open with DEXTER

Illustrative particle trajectories give further insight. The top two panels of Fig. 14 show the trajectory and energy gain of a single sample particle with the same initial conditions in each case. It can be seen that the addition of a guide field substantially alters the particle’s trajectory and energy gain profile. With the addition of a guide field, particles no longer gain energy relatively smoothly, but instead gain and lose energy repeatedly as they travel in the direction of the electric field and then in the opposing direction, although most particles end with a net increase in energy. This was true for all densities and electric fields studied.

The third panel of Fig. 14 shows the variation in the cosine of the particle’s pitch angle with time. It can be seen that adding collisions has a strong effect on particle pitch angle if there is no guide field present. There are many changes in pitch angle, corresponding with the many changes in particle direction as it is scattered back and forth. If there is a guide field present, adding collisions does have a small effect in terms of when the particle pitch angle varies, but does not substantially change the character of that variation. The guide field means that all particles move adiabatically. In the corona, with , collisions away from nulls result in only small perturbations to particle orbits.

The bottom panel of Fig. 14 shows the x-locations at every 100th timestep for 100 particles with an identical set of starting conditions in the presence and absence of a guide field, and with and without collisions. It can be seen that the single factor which most affects the amount of time particles spend near the null is the presence or absence of a guide field. The addition of a guide field causes particles to spend more time near the null, meaning that they can be accelerated to higher energies. The addition of collisions reduces this effect somewhat, as particles may be scattered away from this region, and may also lose energy due to collisions.

Gordovskyy et al. (2012) studied test particle evolution in an MHD model for energy release in twisted flux tube. Electric fields in their model develop parallel to B and effects near null points play a minor role. Collisions in the corona play a minor role in their results, as they do here when the guide field is added.

thumbnail Fig. 13

Final energies of 10 000 electrons at t = 0.015 s, with and without collisions, and with and without a guide field. For each panel, E = 0.001 and ne = 1010 cm-3.

Open with DEXTER

thumbnail Fig. 14

Top panel: trajectory in x-y plane of a particle with identical starting conditions in the presence and absence of a guide field, and with and without collisions. In each case, the particle started the simulation at x = −0.079, y = 0.046, z = 0. and with an initial energy of 1.26 keV. These are the same initial conditions as for the particle whose trajectory is plotted in Fig. 1. Second panel: energy gain vs time for a particle with identical starting conditions in the presence and absence of a guide field, and with and without collisions. Third panel: cosine of particle pitch angle vs. time for a particle with identical starting conditions in the presence and absence of a guide field, and with and without collisions. Bottom panel: x position at every 100th timestep for 100 particles with an identical set of starting conditions in the presence and absence of a guide field, and with and without collisions. For each panel, E = 0.001 and ne = 1010 cm-3, i.e. . The legend shown in the bottom panel applies to all panels.

Open with DEXTER

7. Conclusion and discussion

In this paper, we have developed a method for including collisional scattering and energy loss in the calculation of test particle trajectories. This method was developed by extending the stochastic RK2 method of Honeycutt (1992) to an RK4 method, which was then tested using the problem of pitch angle scattering in an unmagnetised plasma, as studied in MacKinnon & Craig (1991).

The stochastic integrator was used to add collisions to particle trajectories at an X-type neutral point. The addition of collisions causes the particles to lose energy, but because the particles are scattered in pitch angle, some of them may return to the neutral point to be repeatedly energised by the electric field (although, as we have seen, the effect is small).

To gain understanding we have studied a highly idealised system so we note some of the limitations in our calculations. First, we neglected the diffusion in speed that becomes important as particles approach the thermal speed (e.g. Galloway et al. 2005). They would be correct for low temperature plasma and are instructive nonetheless. Our results in respect of acceleration to X-ray emitting energies are not affected by this omission.

We considered an initial value problem. In reality new, unaccelerated particles would continually advect into the diffusion region and the distributions attained would result from a convolution of the results given here with the time profile for the arrival of new material into the reconnection region (this point is also discussed by Hannah & Fletcher 2006).

We have chosen to concentrate on position and angle-integrated energy distributions, the most important quantity for calculating HXR emission. Test particle calculations of this sort can in principle yield more detailed information on the spatial and angular distributions of fast electrons. Spatial distribution in particular could be compared with spatially resolved HXR observations, e.g. from RHESSI although our highly simplified, 2D geometry will limit the value of such an exercise. Also much larger numbers of test particles would be necessary to give statistically useful results with phase space broken into smaller and smaller bins.

The proportion of particles left with nonzero energy decreases with decreasing electric field. However, if the electric field is less than the Dreicer field, collisions become more important, and more particles survive for a stopping time − though not necessarily with HXR-emitting energies − than are found for a small super-Dreicer field. Increasing the plasma density causes the acceleration process to become more efficient, up to a critical value at which the efficiency of the energisation decreases due to collisional energy losses. This critical density changes depending upon the electric field strength used in the simulation. A more detailed discussion of this finding would need to include the diffusion in velocity space that becomes more important as particles approach thermal energies (Galloway et al. 2005).

We also investigated the effect of adding a small guide field. It was found that adding a guide field caused more efficient acceleration of particles as the guide field caused particles to remain closer to the null. When collisions were also added, this additional energy gain was reduced as particles were scattered away from the null and also lost energy due to collisions.

Except for the case with guide field, collisions play a more important role in our results than in the work of Gordovskyy et al. (2012). There are two reasons for this. As already discussed, the presence of a null and the resulting non-adiabatic behaviour allow collisions to have a disproportionate influence on particle evolution. Also our accelerated particles are trapped near the null for times long enough for collisions to become cumulatively important. The geometry of Gordovskyy et al. (2012) does not allow such long-term trapping to happen.

Of course, in a more realistic model, the electric field is unlikely to be constant. Likewise, the plasma density will also change with position within the corona and may also be subject to temporal fluctuations. The variations in both electric field and density will cause changes in the efficiency of particle energisation. Whether collisional slowing down or energisation from the electric field is the dominant effect depends upon the exact values of the electric field and density.

The calculations here may be relevant to understanding coronal hard X-ray sources. As noted by Zharkova et al. (2011), however, the small volume associated with the reconnection region means that too few electrons would be accelerated to account for typical flare HXR fluxes, and a single region of the sort studied here can only be one ingredient of a model for these sources. More complex situations involving many nulls (see Parnell et al. 2010) would have to be involved. We note that such a model would add the coronal hard X-ray sources to the class of “Continuous Reacceleration Thick Target Models” described in a general way by Brown et al. (2009), in which a comparatively small number of electrons produces a large HXR fluence because the accelerator continues to act on them while they radiate.

Appendix A

In order to develop and test an RK4 algorithm, we considered a problem which already has a known solution. The problem used was that considered in MacKinnon & Craig (1991), which dealt with pitch-angle scattering of particles in a non-magnetised medium. The Fokker-Planck (FP) equation for this problem also has a known analytical solution (for the spatially homogeneous case) which is given in terms of the Legendre polynomials, which acts as a further check for the stochastic RK4 solution. MacKinnon & Craig (1991) developed a stochastic system for calculating the variation in particle pitch angle that makes use of the Ito form of a stochastic differential equation (s.d.e.). The FP equation can be replaced by a system of stochastic differential equations. As shown by MacKinnon & Craig (1991) this general equivalence in this particular case means that μ evolves according to the s.d.e. (A.1)where r(t) is a Gaussian random noise process. The initial distribution is monoenergetic. Speeds are normalised to the initial speed (v0), distances are normalised to and times are normalised to , where ne is the electron density of the plasma and Λ is the Coulomb logarithm. It should be noted here that v is also evolving with time; the particles are slowing down monotonically. This can be integrated using the Euler method or by using stochastic RK4 with a noise term (A.2)We carried out a comparison of the 2 methods. A particular example is shown in Fig. 1, at t = 0.06 (the stopping time for these particles is t = 1/3). This shows that all three solutions are in close agreement. Similarly good agreement is found for later times. This apparently simple process, of adding a Wiener noise term to each of the RK4 iterates, is justified in detail by Honeycutt (1992).

thumbnail Fig. 1

Comparison of stochastic RK4, Euler integration, and the exact solution evaluated using Legendre polynomials, t = 0.06.

Open with DEXTER

Acknowledgments

Solar physics research in Glasgow is supported by STFC Rolling Grant ST/I001808/1. C.A.B. thanks STFC for the support of a CASE Studentship held jointly between University of Glasgow and the British Antarctic Survey. We thank Mervyn Freeman for support and useful discussions. Some of A.L.M’s contribution was carried out while on study leave at CRAAM, Mackenzie Presbyterian University, São Paulo and he thanks P. Kaufmann and colleagues for support and hospitality and FAPESP for financial support. P.P. acknowledges support from Dr Helen Mason and the Atomic Astrophysics Group (DAMTP, Univ. Cambridge). We also thank the referee for many useful comments and questions which helped us to improve the paper.

References

All Tables

Table 1

Summary of results of simulations, showing decrease in the number of particles left in the simulation at t = 0.015 s for decreasing super-Dreicer fields.

Table 2

Summary of results of simulations.

All Figures

thumbnail Fig. 1

Trajectory of one electron, integrated with and without the addition of pitch angle scattering. Both electrons had identical starting conditions. All electrons start with an energy of 1.3 keV. The electric field has magnitude 0.001 in our dimensionless units, which is 3.9 V/m for our typical, coronal parameters. The particle started at x = −0.079, y = 0.046 and z = 0, and with angles φ = 1.57 and θ = 2.61.

Open with DEXTER
In the text
thumbnail Fig. 2

Energies at t = 0.015 s for 10 000 electrons whose trajectories have been integrated in a vacuum (black), and in a medium with ne = 1010 cm-3 (red). All electrons start with an energy of 1.3 keV. The electric field has magnitude 0.001 in our dimensionless units, which is 3.9 V/m for our typical, coronal parameters.

Open with DEXTER
In the text
thumbnail Fig. 3

Energies at t = 0.015 s for 10 000 electrons whose trajectories have been integrated in a vacuum (black) and with collisions (red; with ). All electrons start with an energy of 1.3 keV. The electric field has magnitude 10-7 in our dimensionless units, which is 0.00039 V/m for our typical, coronal parameters.

Open with DEXTER
In the text
thumbnail Fig. 4

Energies at t = 0.015 s for 10 000 electrons whose trajectories have been integrated considering only collisional energy loss (black) and with pitch angle scattering included (red). All electrons start with an energy of 1.3 keV. The dimensionless electric field has magnitude 10-3.

Open with DEXTER
In the text
thumbnail Fig. 5

Energies at t = 0.015 s for 10 000 electrons whose trajectories have been integrated considering only collisional energy loss (black) and with pitch angle scattering included (red). All electrons start with an energy of 1.3 keV. The dimensionless electric field has magnitude 10-7.

Open with DEXTER
In the text
thumbnail Fig. 6

Distance from the neutral point (r) for each particle at each saved timestep, considering only collisional energy loss (top) and with pitch angle scattering included (bottom). The electric field has magnitude 0.001 in our dimensionless units, which is 3.9 V/m for electrons. This simulation was done for a sample of 100 electrons, and their positions were recorded every 100 timesteps.

Open with DEXTER
In the text
thumbnail Fig. 7

Energies at t = 0.003 s for 10 000 electrons which have been moving in a plasma of density ne = 5 × 1010 cm-3, for a variety of electric fields. All electrons start with an energy of 1.3 keV. The initial positions of the particles were randomly distributed within the region r < 1, and the initial particle pitch angles were randomly chosen to be between 0 and 2π. For this simulation, .

Open with DEXTER
In the text
thumbnail Fig. 8

Energies at t = 0.0015 s for 10 000 electrons which have been moving in a plasma of density ne = 1011 cm-3, for a variety of electric fields. All electrons start with an energy of 1.3 keV. The initial positions of the particles were randomly distributed within the region r < 1, and the initial particle pitch angles were randomly chosen to be between 0 and 2π. For this simulation, .

Open with DEXTER
In the text
thumbnail Fig. 9

Energies at t = 0.0003 s for 10 000 electrons which have been moving in a plasma of density ne = 5 × 1011 cm-3, for a variety of electric fields. All electrons start with an energy of 1.3 keV. The initial positions of the particles were randomly distributed within the region r < 1, and the initial particle pitch angles were randomly chosen to be between 0 and 2π. For this simulation, .

Open with DEXTER
In the text
thumbnail Fig. 10

Energies at t = 0.00015 s for 10 000 electrons which have been moving in a plasma of density ne = 1012 cm-3, for a variety of electric fields. All electrons start with an energy of 1.3 keV. The initial positions of the particles were randomly distributed within the region r < 1, and the initial particle pitch angles were randomly chosen to be between 0 and 2π. For this simulation, .

Open with DEXTER
In the text
thumbnail Fig. 11

Variation in number of particles remaining in the X-type neutral point simulations discussed in Sects. 4 and 5. Measurements were taken after one stopping time for a variety of different values of electric field, E.

Open with DEXTER
In the text
thumbnail Fig. 12

Final energies of particles at t = 0.15 s for E = 0.001 and two different densities. The simulation time was ten times the stopping time for ne = 1010 cm-3, and equivalent to the stopping time for ne = 109 cm-3.

Open with DEXTER
In the text
thumbnail Fig. 13

Final energies of 10 000 electrons at t = 0.015 s, with and without collisions, and with and without a guide field. For each panel, E = 0.001 and ne = 1010 cm-3.

Open with DEXTER
In the text
thumbnail Fig. 14

Top panel: trajectory in x-y plane of a particle with identical starting conditions in the presence and absence of a guide field, and with and without collisions. In each case, the particle started the simulation at x = −0.079, y = 0.046, z = 0. and with an initial energy of 1.26 keV. These are the same initial conditions as for the particle whose trajectory is plotted in Fig. 1. Second panel: energy gain vs time for a particle with identical starting conditions in the presence and absence of a guide field, and with and without collisions. Third panel: cosine of particle pitch angle vs. time for a particle with identical starting conditions in the presence and absence of a guide field, and with and without collisions. Bottom panel: x position at every 100th timestep for 100 particles with an identical set of starting conditions in the presence and absence of a guide field, and with and without collisions. For each panel, E = 0.001 and ne = 1010 cm-3, i.e. . The legend shown in the bottom panel applies to all panels.

Open with DEXTER
In the text
thumbnail Fig. 1

Comparison of stochastic RK4, Euler integration, and the exact solution evaluated using Legendre polynomials, t = 0.06.

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.