PARTICLE-IN-CELL SIMULATION OF ELECTRON ACCELERATION IN SOLAR CORONAL JETS

Solar jets have been shown to be triggered by magnetic reconnection, similar to solar flares, while their released energy is much lower than what is released in a medium-sized flare event, and the timescales are usually shorter. Nevertheless, their high frequency of occurrence makes them a significant contributor to the solar ejecta, particularly the solar wind originating from coronal holes. Solar jets feature upflow velocities of more than 150 km s⁻¹ (Savcheva et al. 2007; Chifor et al. 2008) and are observed at EUV down to X-ray wavelengths primarily in coronal holes (Kamio et al. 2007), but also in active regions (Chifor et al. 2008). Using RHESSI data, Krucker et al. (2007, 2011) have further investigated the electron impact regions of solar jets as a result of interchange reconnection.

There have been several studies in the past employing fully three-dimensional kinetic models to study particle acceleration. Two approaches are most popular: particle-in-cell (PIC) simulations and test particle simulations. Their main difference is the back-reaction from the particles onto the fields, which is only taken into account in the former method, while it is assumed to be negligible in the latter one, justified by using a low number of test particles for such simulations. These kind of simulations provide much information on particle trajectories and favored locations of particle acceleration (Turkmani et al. 2005, 2006; Dalla & Browning 2005, 2006, 2008). Nevertheless, there are severe limitations to test particle simulations. These limitations and their consequences have been discussed in detail in Rosdahl & Galsgaard (2010). On the other hand, the approach using self-consistently evolving fields such as in PIC codes is subject to several resolution constraints, reducing the possible physical box size that can be simulated to far below length scales of solar jets. In order to bypass these limitations, we made use of modifications of the constants of nature. Thereby we are able to present results from fully three-dimensional PIC simulations of a solar jet, using essentially the same initial setup as in Rosdahl & Galsgaard (2010), but with self-consistently evolving fields. Such modifications have previously to some extent been used by Drake et al. (2006) and Siversky & Zharkova (2009). A comparison of simulations using different modifications exceeds the scope of this Letter. We therefore refer to Baumann et al. (2012) for a description and analysis of the physical consequences of modifications. For the present study we use a choice of modifications as close as we can get to reality with the currently available computational resources.

Section 2 provides an overview of the experiment and describes the MHD and PIC simulations as well as their intercoupling. In Section 3 the results are presented and discussed. Finally, in Section 4 conclusions are drawn.

The solar jet experiment at hand starts out with a fully three-dimensional resistive compressible MHD simulation of a twisted emerging flux rope, initially positioned 1.7 Mm below the photosphere, similar to the setup used by Moreno-Insertis et al. (2008). A constant magnetic field of 3.3 G is imposed on the entire computational box, inclined 65° in the yz-plane. The maximum magnetic field strength of the flux rope is slightly higher than 1000 G and hence much larger than the background magnetic field. The atmosphere is initially in hydrostatic equilibrium with a one-dimensional atmospheric profile similar to the one used in Archontis et al. (2005): the sub-photospheric temperature at the bottom is 5.5 × 10⁴ K, with a maximum mass density ρ of about 9 × 10⁻⁶ g cm⁻³ at a depth of 3.7 Mm below the surface. The "chromosphere" has a constant temperature of around 5600 K and the corona starts out with T = 2.2 × 10⁶ K and ρ = 6 × 10⁻¹⁶ g cm⁻³, as illustrated in Figure 1.

Zoom In Zoom Out Reset image size

The simulations are performed using the Stagger MHD code, as in Moreno-Insertis et al. (2008), assuming an ideal gas law and neither taking heat conduction nor radiative cooling into account. Viscosity and resistivity are locally defined, depending mainly on the velocity divergence, in order to provide a suitable (to maintain code stability), but minimal amount of dissipation. A mesh with 512³ cells covers a box with physical extents of 33.8 × 38.1 × 32.5 Mm. The mesh is non-uniform in all directions, with a minimal mesh spacing of (x_min, y_min, z_min) = (0.034, 0.034, 0.030) Mm around the reconnection region, and with a mesh spacing less than 10% larger than that in a region of size 12.5 × 10.7 × 5.1 Mm. The z-coordinate is the direction normal to the solar surface.

Despite the slightly denser corona compared to the previously performed simulations by Moreno-Insertis et al. (2008), the overall evolution of the experiment is the same: the twisted flux tube is made buoyantly unstable by applying a density perturbation at its center, which causes it to rise up into the much rarer corona, against the Lorentz force from the bending tube, while expanding as described in Archontis et al. (2004). Above the photosphere, the expansions of the magnetic flux rope, now due to the high magnetic pressure, continue rapidly as more and more flux reaches coronal heights. At the same time, the corona is locally pushed upward where plasma as well as magnetic flux emerges from below (cf. Archontis et al. 2005). The interaction between the two magnetic field domains defined by the initial coronal magnetic field and the flux rope leads to a destabilization of the field configuration and causes the formation of a thin dome-shaped current sheet where the magnetic field lines are most inclined relative to each other. As in Moreno-Insertis et al. (2008), most magnetic field lines of the two domains end up being almost anti-parallel at their first encounter, which makes their interaction maximally powerful. The sheet is subject to ohmic dissipation, causing it to reach temperatures as high as about 8 × 10⁶ K. Reconnection gradually occurs between various field lines along the twisted tube and the ambient coronal magnetic field. Due to this restructuring of the magnetic field, several distinct flux domains form in the corona, as shown in Moreno-Insertis et al. (2008, see their Figure 2). A hot plasma jet pair emerges from the high temperature and gas pressure region of the reconnection region, propagating sideward together with the expelled reconnected magnetic field lines. The plasma in the jets is fed to the reconnection site from both sides of the current sheet; the region below supplies dense and cold photospheric plasma, while the plasma coming from above is much hotter and rarer coronal plasma.

After a large fraction of the dense plasma enclosed in the flux rope has been reconnected and ejected in the form of plasma jets, as well as drained off along the magnetic field lines due to gravity acting on the heavy sub-photospheric gas, we initialize the PIC experiment with a cut-out of size 22 × 22 × 22 Mm from the MHD data set. The reconnection process continues in the PIC simulation expelling coronal low-density plasma and reconnecting field lines from both connectivity domains.

The PIC simulations are performed using the Photon-Plasma code (Haugbølle 2005; Hededal 2005), which solves the Maxwell equations together with the relativistic equation of motion for charged particles. We fix the magnetic fields at the boundaries to the values given by the MHD data set, and leave the boundaries open for particles to exit and enter (T. Haugbølle et al. 2012, in preparation). To initialize the electric field in the PIC simulation, only the advective electric field ($-\textbf {u} \times \textbf {B}$, where u is the bulk speed) is passed on. The particles are initially given a random thermal velocity drawn from a Maxwellian distribution, plus the bulk velocity from the MHD simulation. Electron velocities consist of additionally the velocity due to the initial electric current:

$$\begin{equation} v_{{{J}}} = \frac{1}{\mu _0 q n} (\nabla \times \mathbf {B}), \end{equation} \tag{ 1 }$$

where the magnetic field B and the density n are provided by the MHD snapshot data set. The jet velocity of this specific MHD data set at that point in time is on the order of 400–800 km s⁻¹. These jets are dominated by the thermal motion in this high-temperature and low magnetic field region.

We conducted several PIC runs with grid dimensions of 400³ and 800³ cells on a uniform grid with cell sizes of 55 km and 27.5 km respectively, in each case with 20 particles per species (protons and electrons) per cell, simulating up to 7.5 solar seconds.

To minimize computational constraints, the MHD snapshot is cut at 1.1 Mm below the bottom of the corona, hence in the transition region, shown as a vertical dashed line in Figure 1. This limits the density span to a factor of 4 × 10⁴, small compared to the span of about 2 × 10¹⁰ covered by the MHD simulation, but still large for a PIC simulation. Additionally, the temperature is reduced by a factor of four in order to avoid drowning the high-energy tail in the Maxwellian distribution.

To make the plasma micro-scales marginally resolvable, the charge per particle is reduced, and to ease the time step constraint from the propagation of electromagnetic waves the speed of light is reduced, as explained in Baumann et al. (2012). The electron skin depth in the current sheet is resolved with about 5–10 grid cells. The electron gyroradius varies between a fraction of a cell in the flux tube interior to many cells inside the current sheet—because of the near cancellation of oppositely directed magnetic fields there are several dynamically evolving null points inside the current sheet.

In addition to these stratified atmosphere simulation runs, control runs with a constant density of 3 × 10 − 15 g cm⁻³ were also performed. Both types of runs show similar overall results.

During the MHD flux emergence simulation a diffusive electric field $\mathbf {E_{{\rm res}}} = \eta \mathbf {j}$, where η is the resistivity and j the electric current density, builds up in the reconnection region, approximately cospatial with the current sheet. The parallel (in relation to the magnetic field) electric field E_∥ is part of the diffusive (non-ideal) electric field, and provides information on the rate of reconnection as well as on favored regions for particle acceleration. The diffusive component of the electric field is, unlike the advective electric field ($\mathbf {-u} \times \mathbf {B}$), not inherited by the PIC code, but builds up self-consistently. Figure 2 compares the location of the diffusive electric field component for the chosen snapshot of the MHD simulations with the E_∥ field of the PIC simulation 0.5 s after start. The parallel electric fields reach in general much higher magnitudes in the MHD simulations compared to the PIC simulations.

Zoom In Zoom Out Reset image size

E_∥ is the most efficient particle accelerator, since its force acts on the particles without being affected by the perpendicular particle gyromotion. Its maximum in the PIC simulations is located inside the current sheet, equivalent to the diffusive electric field E_res in the MHD simulation. Accelerated electrons (see Figure 3) are located in the plasma outflow parts of the reconnection region. The electron bulk velocity in the jet is on the order of 2000 km s⁻¹. The proton bulk jet flow on the other hand is only about 270 km s⁻¹, which difference to the electron bulk speed defines the electric current required by the magnetic field configuration. The lower plane of the PIC simulation visualization in Figure 2 (right) shows the electric current density. At the bottom center of this figure resides the flux rope, whose twisted field lines are indicated by the strongest electric current pattern. Additionally, to the left of the flux rope signature, the current sheet features a turbulent structure. This is the result of the plasma transport caused by reconnection. In addition, we observe fast up- and down-flowing plasma, as can be seen in Figure 3, in which upward moving electrons are shown in purple/light gray and downward moving electrons are shown in green/gray.

Zoom In Zoom Out Reset image size

By tracing particles that win energy over a time period of a second, it becomes clear (see Figure 4) that the acceleration mechanism is a systematic DC electric field, particularly present in the proximity of the current sheet, and mainly directed parallel to the magnetic field. However, since the magnetic field is very weak in the current sheet area, particles are almost decoupled from the magnetic field lines, leading to pitch angle fluctuations seen in the third panel to the left in Figure 4. We only write out every 25th particle data, hence in Figures 3 and 4 one particle represents in reality 25 of its type. Note that in Figure 4 the electric field fluctuations in the second panel on the left are to a large extent Monte Carlo noise, due to low numbers of particles per cell. This has been verified with a control experiment using twice the number of particles per cell. Further, the background electric current density image is not exactly at the particle positions, which results in a slight projection offset.

Zoom In Zoom Out Reset image size

Only a limited number of particles happen to be close enough to the current sheet to be efficiently accelerated by E_∥, thus being able to contribute to the electric current required by the magnetic field topology. After experiencing a certain amount of acceleration, they get expelled from the current sheet, as they follow magnetic field lines which leave the current sheet region. The particles lost from the current sheet are replaced by new particles, which again need to be accelerated. Since the magnetic field in the reconnection region is very low, due to the chosen orientation of the initial background magnetic field geometry, it is easy for electrons to get misguided, as they are no longer tightly attached to their magnetic field lines and therefore may encounter different electric field structures on their large gyroradius trajectories.

The systematic parallel electric field building up in the current sheet area is capable of accelerating particles to non-thermal velocities. Figure 5 presents the energy histogram of electrons located in a cut-out of 19.0 × 13.2 × 6.5 Mm around the reconnection region. The initial energy distribution is shown by the dashed line. It is primarily a superposition of the two different plasma inflow domains of the reconnection region passed to the PIC code through the MHD bulk velocity: one of coronal origin, hence at higher temperatures and lower densities, and another one from plasma emerging from the flux rope, hence at lower temperatures and higher densities. To this the drift speed from particles in the current sheet is added, as defined by Equation (1). The distribution of probability for a particle to stay in the current sheet defines a power law in the high-energy tail of the Maxwell–Boltzmann distribution, typically featuring an index of around −0.5 on a logarithmic scale. This is the case for electrons and protons, as the acceleration mechanism inside the current sheet is for both species the same. A series of different simulations indicates that the influence of the physical resolution on the power-law index as well as on the upper energy limit is minuscule (see Figure 5), since the current per unit width of the current sheet is set by the magnetic field geometry. This independence also demonstrates numerical convergence. The tail slopes converge very rapidly toward these power-law indices, since the electron acceleration happens impulsively. Part of the accelerated electrons escape along open magnetic field lines (see Figure 3) into interplanetary space and hence contribute to the solar wind. Downward moving energetic electrons impact the much denser chromosphere, which could lead to bremsstrahlung emission on the Sun. In the case of high-energy electrons, this may cause hard X-ray signatures of the type observed by RHESSI. Figure 5 shows that most energy is deposited at the magnetic field line footpoints of the coronal loops. Observational hard X-ray signatures of solar jets described in Krucker et al. (2011) show a three-footpoint pattern (see, e.g., event 7 in Figure 3). Due to the choice of our flux rope being 90° inclined relative to the uniform background magnetic field, we rather get a three-dimensional structure of footpoints, represented by three lines of electron impact regions in Figure 5. Additionally, because of the magnetic field line orientation and twist, an impact region at the intersection of the flux tube and the photosphere develops. (The impact line to the right is slightly cut because of the choice of the computational box size.)

Zoom In Zoom Out Reset image size

On the basis of an MHD jet experiment, similar to the one conducted by Moreno-Insertis et al. (2008), but using stretched meshes to obtain higher spatial resolution, we have used PIC simulations to study the acceleration of charged particles in the three-dimensional reconnection region of a solar coronal jet. This is the first fully three-dimensional kinetic model employing self-consistent fields to investigate particle acceleration in the context of solar jets. It uses a new concept of combining macroscopic with microscopic simulations.

A strong correlation is found between a slowly evolving DC electric field located inside the current sheet and the location of the accelerated particles (electrons and protons). The magnetic field is weak and chaotic inside the current sheet, from where most of the particles that are accelerated are quickly lost, only to be replaced by new particles, which again need to be accelerated. This represents the dissipation mechanism at work. The systematic electric field required to constantly accelerate new particles is, in effect, a dissipative ("resistive," non-ideal) electric field, sustained even though the plasma particles in this experiment are collisionless. This physical process leads to an energy probability distribution of accelerated particles, featuring a power-law tail at the high-energy side of the initial Maxwell–Boltzmann distribution, with a power-law index of about −1.5. The power-law index and the invariant energy cut-off show numerical convergence of the simulations.

The impact region of energetic electrons on the lower boundary of the box is comparable to observations of similar events and magnetic field geometries, described by Krucker et al. (2011).

In the future, with increasingly available computational resources, we hope to be able to resolve turbulent regions sufficiently, to enable a study of stochastic particle acceleration, and expect to be able to accelerate particles to much higher energies.

Special thanks go to Klaus Galsgaard, Jacob Trier Frederiksen, and Troels Haugbølle for very valuable discussions. The work of G.B. was supported by the Niels Bohr International Academy, the SOLAIRE Research Training Network (MRTN-CT-2006-035484), and the SWIFF FP7 project (no. 263340, http://www.swiff.eu) of the European Commission. We acknowledge that the results have been achieved using the PRACE and John von Neumann Institute for Computing Research Infrastructure resource JUGENE/JUROPA based in Germany at the Jülich Supercomputing Centre. Furthermore, we acknowledge a DECI grant and grants from the Danish Center for Scientific Computing.