TEST-PARTICLE ACCELERATION IN A HIERARCHICAL THREE-DIMENSIONAL TURBULENCE MODEL

Understanding the mechanisms that accelerate charged particles in dynamical plasmas is an intense subject of research, with applications to astrophysical, space, and laboratory environments (Cargill et al. 2012). Charged particles are energized in solar flares, planetary magnetospheres, at interplanetary shocks, and in the interstellar medium. In the solar corona, electrons can be accelerated to tens of MeV, while ions gain energies up to several GeV, and high-energy particles may carry up to 30% of the energy budget of a flare (for a review, see Aschwanden 2002; Klein & MacKinnon 2007). The fast release of energy in solar flares is commonly attributed to magnetic reconnection events. In such complex natural systems, several processes may contribute to the observed particle energization, and these may be distinguished in part by the locations and scales at which they occur. In this way, the resulting acceleration might be viewed as a multistage process (see, e.g., de Jager 1986; Gordovskyy et al. 2005; Brown et al. 2009).

In particular, the role of scattering from magnetic irregularities is well known to be an important factor (Jokipii 1966): particles can be accelerated by the large-scale electric fields (Speiser 1965; Litvinenko 1996), as well as by the smaller scales of reconnecting fields (Drake et al. 2005) and/or discontinuities, including acceleration in shocks (Giacalone & Jokipii 1996; Giacalone 2005). Furthermore, they are also accelerated by resonances with structures at various scales, and acceleration mechanisms may be strongly influenced in various contexts by changes of magnetic connectivity (e.g., Rappazzo et al. 2012). Thus, essentially all length scales may be important in the investigation of the energization process, giving rise to a rather complex behavior in the case of turbulent fields, characterized by the presence of energy at all scales from large to small and of current and vorticity coherent structures that can cause resonant or nonresonant interactions (Matthaeus et al. 1984; Kobak & Ostrowski 2000). Indeed, although the smooth or average electromagnetic fields are associated with a variety of possible particle motions (Speiser 1965; Sonnerup 1971; Cowley 1978), magnetic fluctuations can strongly influence both the local electric field and particle motion. For example, strong resonant pitch angle scattering, as well as the formation of transient trapping centers such as secondary islands, are directly associated with the presence of electromagnetic fluctuations (Matthaeus & Lamkin 1986; Ruffolo et al. 2003; Drake et al. 2006).

Such effects are expected to be widespread since turbulent magnetic fluctuations are observed in space plasmas in practically all environments. Fully developed turbulence itself involves a hierarchical process in which many scales of motion are engaged in the energy transfer. As particles become energized, their interactions with the turbulence sample a correspondingly wide range of scales. From an observational point of view, different spectra and intensity variations have been detected from sources of nonthermal radiation in many flares (see, e.g., Takakura et al. 1995; Emslie et al. 2003; Fletcher et al. 2007), which is also hard to explain invoking only a single acceleration mechanism.

Because of the complexity of particle motion, a useful approach has been to resort to numerical experiments based on test particle simulations, in which no interaction among the particles and no back reaction to the imposed electromagnetic fields is considered (Litvinenko & Somov 1993; Turkmani et al. 2005; Dalla & Browning 2005; Onofri et al. 2006). However, processes that span a wide range of scales pose profound challenges because the different scales require different time and spatial resolutions, and qualitatively different phenomenologies may apply. For instance, in the description of a coronal loop, the typical scale of the system L measures thousands of kilometers. In contrast, the ion inertial scale d_i, the typical scale at which dissipative effects become important (Smith at al. 2001), is of the order of a meter and is comparable with the typical gyroradius, ρ = v_th/Ω, of a particle moving at the thermal speed v_th with a gyrofrequency, Ω.

Such conditions, spanning a wide range of length scales, are typical in space plasmas. For example, in the solar wind, it can be estimated that L/d_i ∼ 10⁵. In a coronal hole (at about one solar radius) and in magnetically confined regions of the solar corona, L/d_i ∼ 10⁶. Attempting to resolve this full range of scales requires an enormous cost in terms of computing power at present standards. In galactic astrophysics, even larger L/d_i may be relevant. To attack such problems, approximate models, such as the guiding center and gyrokinetics, have been employed to describe particle motion in a strong magnetic field (see, e.g., Lee 1983; Browning & Vekstein 2001; Giuliani et al. 2005; Gordovskyy et al. 2005; Lehe et al. 2009). These models are derived using magnetic moment μ conservation, so their validity is limited to specific circumstances. In fact, the equation of motion of the particle guiding center can be solved in the case of small gyroradius and large gyrofrequency, ρ ≪ l and Ω ≫ T⁻¹, where l and T should be interpreted as the typical length and timescales of the local electromagnetic field. As a consequence, even if particles are strongly magnetized globally, non-adiabatic behavior may develop at local scales (Brown et al. 2009; Dalena et al. 2012b), in particular in the presence of sharp field gradients, near reconnection sites, or in the presence of well-developed MHD turbulence. In these situations, the presence of resonant effects and/or coherent structures could violate magnetic moment μ conservation.

To represent turbulent fields at all relevant scales, one approach is to use gridless methods with prescribed analytical fields (see, e.g., Giacalone & Jokipii 1996; Qin et al. 2002; Bykov et al. 2008; Fraschetti & Melia 2008; Fraschetti & Giacalone 2012). In principle, this method can represent a wide range of scales by correspondingly increasing the number of modes retained in the representation. However, the direct gridless implementation becomes extremely onerous in that case. These models also most often employ random phases which prevent representation of coherent small scale structures, such as current sheets and vorticity concentrations, that can play a key role in particle acceleration.

Therefore, we have developed a new multi-box, multiscale test-particle technique to resolve the wide range of time and length scales usually involved in the description of space and astrophysical plasmas. The main idea is to study particles dynamic at different steps during their evolution, moving from several simulation boxes of different dimensions (for a schematic picture see Figure 1). The nonrelativistic equations of motion of a particle moving in an electromagnetic field are solved numerically, keeping approximately constant the resolution of the particle gyroradius ρ going from a smaller to a larger box. This has been accomplished by choosing a grid scale Δx for every box that satisfies the condition Δx < ρ. Different spatial and temporal resolutions in different regions allow us to treat all the scales involved in the system with the same accuracy. In this paper, we present this new numerical technique applied to the study of particle acceleration in a coronal loop. We report results that describe acceleration in direct simulations of reduced magnetohydrodynamic (RMHD) turbulence at different stages, over a broad range of time, length, and energy scales.

Zoom In Zoom Out Reset image size

The paper is organized as follows. In Section 2, we describe the equations and properties of the MHD fields and of the particles. In Section 3, the multi-box technique is introduced. In Section 4, results are given for energization of protons obeying the nonrelativistic equations of motion. Section 5 contains our conclusions.

The procedure consists of an RMHD simulation, appropriate for the low frequency dynamics of a plasma in a strong externally supported magnetic field, B₀ = B₀e_z, followed by integration of test particle orbits. The RMHD equations are solved in a Cartesian box with an aspect ratio of l_∥/l_⊥ = 10 (ratio of characteristic lengths parallel and perpendicular to B₀), and a numerical grid of 256 × 256 × 120 points. For more details on the numerical methods and code, see Rappazzo et al. (2007, 2008).

Particle trajectories are integrated in a single snapshot of the RMHD field, taken at t = t_nl = 10t_A, where t_nl is the nonlinear time and t_A is the Alfvén crossing time. The nonlinear stage is characterized by the presence of current sheets elongated along the axial direction and represents a statistical steady state in which energies fluctuate around a mean value, and on average, total dissipation and Poynting flux are in balance. The quasi-static approach in which magnetic and electric fields do not change with time is justified by assuming that the particle acceleration time is shorter than the characteristic time of MHD system evolution.

Another possible way to attack the multiscale problem is to study the acceleration of particles moving in an analytic electromagnetic field (Tu & Marsch 1993; Gray et al. 1996; Bieber et al. 1996), given as a Fourier series. In the limit of an infinite number of wave modes, the turbulence can approach a pure symmetry and may be, for example, isotropic and spatially homogeneous (Batchelor 1960). Methods that sum a sparse set of plane waves with randomly distributed propagation direction, and with random polarizations and phases, may be referred to as gridless methods (Giacalone & Jokipii 1999). These have been employed in studies of the transport of cosmic rays performing numerical simulations of test-particles interacting with a time stationary turbulent magnetic field. These methods typically use a smaller number of Fourier modes due to the requirement of direct summation of the series at many points. These methods also lack coherent structures, such as current sheets.

While it has its own limitations, the multi-box technique we introduce here does allow for fields that have coherent structures and, given the large number of modes it permits, it allows for development of resonances and related diffusion processes that require a high density of modes. Furthermore, by rescaling the fluctuations across multiple boxes, the multi-box technique permits the study of particle acceleration including these important effects: coherent structures and resonances across all relevant scales during the multistage acceleration process.

2.1. Reduced MHD Simulation

A coronal loop is a closed magnetic structure threaded by a strong axial field, with the footpoints rooted in the photosphere. This makes it a strongly anisotropic system, as measured by the relative magnitude of the Alfvén velocity, $v_{_{\!A}} = 2000$ km s⁻¹ (associated to the strong axial field B₀), compared to the typical photospheric velocity, v_ph ∼ 1 km s⁻¹. These relatively slow motions shuffle the footpoints of the axially elongated magnetic field lines giving rise to a coronal magnetic field b that is small compared to the guide field δb/B₀ ≪ 1, and the plasma dynamics are then expected to develop in the reduced MHD regime (Strauss 1976) (for a more detailed analysis see Einaudi et al. 1996; Dmitruk et al. 1997). The loop dynamics may be studied in a simplified geometry, neglecting any curvature effect, as a "straightened out" Cartesian box, with an orthogonal square cross section of size L = l_⊥, and an axial length l_z embedded in an axial homogeneous uniform magnetic field, B₀ = B₀e_z.

The dynamics are integrated with the (nondimensional) equations of RMHD (Kadomtsev & Pogutse 1974; Strauss 1976; Montgomery 1982), well suited for a plasma embedded in a strong axial magnetic field. In dimensionless form, they are given by:

$$\begin{equation} \partial _t \mathbf {u} + \mathbf {u} \cdot \nabla \mathbf {u} = - \nabla P + \mathbf {b} \cdot \nabla \mathbf {b} + c_{_{\!A}} \partial _z \mathbf {b} + \frac{ 1 }{Re} \nabla ^{2} \mathbf {u}, \end{equation} \tag{ 1 }$$

$$\begin{equation} \partial _t \mathbf {b} + \mathbf {u} \cdot \nabla \mathbf {b} = \mathbf {b} \cdot \nabla \mathbf {u} + c_{_{\!A}} \partial _z \mathbf {u} + \frac{ 1 }{Re_m} \nabla ^{2} \mathbf {b}, \end{equation} \tag{ 2 }$$

$$\begin{equation} \nabla \cdot \mathbf {u} = \nabla \cdot \mathbf {b} = 0. \end{equation} \tag{ 3 }$$

Here, gradient and Laplacian operators have only transverse (x–y) components as do velocity and magnetic field vectors (u_z = b_z = 0), while P is the total (plasma plus magnetic) pressure. Incompressibility in RMHD equations follows from the large value of the axial magnetic field (Strauss 1976), and they also remain valid for low plasma β systems (Zank & Matthaeus 1992; Bhattacharjee et al. 1998) such as the corona, where the plasma β ∼ 0.01 in closed regions. The plasma is then assumed to have a uniform constant density, ρ₀.

To render the equations nondimensional, we have first expressed the magnetic field as an Alfvén velocity $[ b \rightarrow b/\sqrt{4\pi \rho _0} ]$, and then all velocities have been normalized to the velocity u* = 1 km s⁻¹, the order of magnitude of photospheric convective motions. Lengths and times are expressed in units of the perpendicular length of the computational box, ℓ* = ℓ_⊥, and its related crossing time, t* = ℓ*/u*. As a result, the linear terms ∝∂_z are multiplied by the dimensionless Alfvén velocity c_A = v_A/u*, where $v_A = B_0 / \sqrt{4\pi \rho _0}$ is the Alfvén velocity associated with the axial magnetic field.

We adopt a standard simplified diffusion model in which both the magnetic resistivity η_m and dynamic viscosity ν are constant and uniform. The kinetic and magnetic Reynolds numbers are then given by

$$\begin{equation} Re = \frac{\rho _0\, \ell ^\ast \, u^\ast }{\nu }, \quad Re_m = \frac{4 \pi \, \ell ^\ast \, u^\ast }{\eta _m c^2}, \end{equation} \tag{ 4 }$$

where c is the speed of light, and numerically, they are given the same value Re = Re_m = 400.

We numerically solve Equations (1)–(3) written in terms of the potentials in Fourier space, i.e., we advance the Fourier components in the x- and y-directions of the magnetic field and velocity scalar potentials (for a more detailed description of the numerical code and methods, see Rappazzo et al. 2007, 2008). Along the z-direction, no Fourier transform is performed so that we can impose nonperiodic boundary conditions, and a central second-order finite-difference scheme is used. In the x–y plane, a Fourier pseudospectral method is implemented. Time is discretized with a third-order Runge–Kutta method. We use a computational box with an aspect ratio of 10, which spans 0 ⩽ x, y ⩽ 1, 0 ⩽ z ⩽ 10, and implement a numerical grid with 256 × 256 × 120 points.

As boundary conditions at the bottom and top plates (z = 0 and z = 10, the photospheric-mimicking surfaces), we impose two independent velocity patterns intended to mimic photospheric motions, made up of large spatial scale projected convection cell flow patterns constant in time (e.g., see Figures 1 and 2 in Rappazzo et al. 2008). These are essentially large-scale disordered vortices that stir the footpoints of the line-tied axially elongated magnetic field lines. Similar to standard forced MHD turbulence simulations, this boundary forcing injects energy into the system, giving rise to an MHD turbulence cascade dominated by magnetic energy (Einaudi et al. 1996; Dmitruk et al. 1997, 2003; Rappazzo et al. 2008). In physical space, current sheets form on very fast ideal timescales (Rappazzo & Parker 2013) and are aligned to the strong axial magnetic field (e.g., see Figure 18 in Rappazzo et al. 2008).

Given its quasi-2D properties (perpendicular magnetic and velocity fluctuations, current sheets elongated in the axial direction, and the presence of a strong guide field), RMHD turbulence is an adequate model in which to study particle acceleration, serving as a first simplified model with three-dimensional (3D) characteristics, in between the more simple and more extensively studied two-dimensional (2D) case and the more complex and challenging fully 3D setup. Furthermore, RMHD represents a good model for the low corona, largely characterized by the presence of a strong guide field, in both open and closed regions.

2.2. Particle Integration

The behavior of a test particle in electromagnetic fields obtained from RMHD is described by its time dependent position r(t) and 3D velocity v(t), advanced according to the nondimensional equations:

$$\begin{equation} \frac{d{\bf r}}{dt} = {\bf v}, \end{equation} \tag{ 5 }$$

$$\begin{equation} \qquad\frac{d{\bf v}}{dt} = \beta ({\bf E} + {\bf v} \times {\bf B}). \end{equation} \tag{ 6 }$$

The nondimensional electric field E is obtained through Ohm's law (normalized with the field E₀ = v_AB₀/c), from the RMHD simulation magnetic field B, plasma velocity u, and current density j:

$$\begin{equation} {\bf E} = -\frac{1}{c} {\bf u} \times {\bf B} + \eta {\bf j}. \end{equation} \tag{ 7 }$$

The nondimensional resistivity η = 1/Re_m is equal to 2.5 × 10⁻³. To render the equations nondimensional, we use as characteristic quantities the turbulence correlation length λ_c = l_⊥/4 = 1000 km, the Alfvén velocity $v_{_{\!A}} = 200$ km  s⁻¹, and the Alfvén crossing time $t_{_{\!A}} = \lambda _c/v_{_{\!A}}$ (Dalena et al. 2012a).

The dimensionless parameter $\beta =\Omega t_{_{\!A}}$ in Equation (6) couples particle and field spatial and temporal scales and provides a particularly useful means to relate numerical experiments to space and astrophysical plasmas. Indeed it can be rewritten as (Ambrosiano et al. 1988; Dalena et al. 2012b):

$$\begin{equation} \beta = \frac{q}{m} \frac{\sqrt{4\pi \rho _0}}{c} \lambda _c = \frac{\lambda _c}{d_{i}}, \end{equation} \tag{ 8 }$$

where $d_i = c/ \omega _i =(c^2 m_i / 4 \pi n_i q_i^2)^{1/2}$ is the ion inertial scale. Thus, the β parameter measures the range of scales involved in a particular system. In general, in a turbulent collisionless plasma, the bandwidth of the inertial range fluctuations may extend from large fluctuations near the correlation scale λ_c to small fluctuations near the ion inertial scale d_i. In this case, one expects that β ≫ 1. For the typical values in our simulations, d_i ∼ 2 m, and β = λ_c/d_i ≃ 5 × 10⁵.

Considering that the average value of the coronal temperature ∼10⁶ K corresponds to a velocity of ∼200 km  s⁻¹ for protons, the nonrelativistic approximations used in Equations (5) and (6) is adequate to investigate the acceleration processes up to values of $\gamma = 1/\sqrt{1-(v/c)^2} \sim 1.1$, corresponding to energies of ∼0.1 GeV for protons. The simulations presented here consider durations that allow protons to gain at most energies of the order of 0.1 GeV. In order to attain longer integration times, the relativistic analog of Equations (5) and (6) should be used, and this is left to future investigations.

The appearance of the large number β in Equation (6) makes it apparent why it is challenging to integrate test particle evolution with parameters appropriate to space and astrophysical plasmas. Equations (5) and (6) are very stiff (involving greatly differing time scales), while, as some particles may accelerate until gyroradii become comparable to L, there is a huge range of spatial scales involved. Therefore, it is very difficult to generate numerically relevant magnetic fields that have structures across the required range of scales.

Moreover, when dealing with numerical simulation of test-particles in the electromagnetic field, one of the most important requirements is the resolution of the particle gyroradius ρ during the integration. Indeed, ρ determines the minimum scale of a magnetic field structure capable of affecting the dynamics of the particle. If the value of ρ is smaller than the grid scale value Δx, the test-particle will not be able to sample the different structures of the field, and it will move as if it were under the action of a constant field. The parameter β is related to the value of the particle gyroradius by the parameter $\epsilon = \rho /\lambda _c = v/\Omega = (v/v_{_{\!A}})/\beta$ (Mace et al. 2000), sometimes called the dimensionless particle rigidity.

The typical velocity of a particle traveling at the thermal speed in the solar corona is usually a fraction of the characteristic Alfvén speed $v_{_{\!A}}$, thus ρ ⩽ λ_c/β. The wider the range of length scales involved, the smaller the gyroradius will be in comparison with the largest scale present in the system. For instance, with the choice of our parameters, ∼ 9 × 10⁻⁷. On the other hand, the grid resolution scale Δx = L/N of the simulation box is determined by the number of grid points N and the box length L. The resolution of the particle gyroradius in this case would require a number of grid points of the order of 10⁷ in the x, y, and z directions, not achievable at the present standard.

Another subtle point in the description of test-particle motion is the potential for resonances that can contribute significantly to acceleration. If particles are injected with a gyroradius much smaller than the outer scale of the system λ_c, it becomes possible that the simulation will lack resonances. The dimensionless particle rigidity can be related to the bend-over wavenumber of the turbulent energy spectrum, k_bo = 1/λ_c, and the minimum resonant wavenumber, $k^{r}_{\rm min} = 1/\rho$, as $\epsilon = k_{\rm bo}/k^{r}_{\rm min}$. When ρ ≫ λ_c, particles experience all possible k-modes in few gyroperiods, resonating with the energy containing scale ($k^{r}_{\rm min} \ll k_{\rm bo}$). For lower energies, the test particles resonate in the inertial range ($k^{r}_{\rm min} \sim k_{\rm bo}$). However, when ρ ≪ λ_c, particles will require higher and higher wavenumbers to resonate with but it is possible that the numerical representation may not have sufficient power in the high wavenumber magnetic energy spectrum to produce these resonant effects in the numerical experiment (see, e.g., Figure 6 of Dalena et al. 2012b).

To deal with the numerical issues discussed above, we have developed a multi-box, multiscale test-particle approach to numerically integrate the nonrelativistic equations of motion of a particle in electromagnetic fields. Particle dynamics is studied at different stages during the energization process, moving between simulation boxes of different dimensions and spatial resolutions, each satisfying the condition Δx < ρ. With this technique, we have been able to investigate particle motion and acceleration at different length scales, retaining the resolution of the particle gyroradius simultaneously at each step. The method allows us to treat lengths of the order of few meters (comparable with the ion inertial scale in the system) as well as lengths up to several thousand kilometers (the size of a coronal loop) with the same accuracy. This new technique is useful to identify the mechanisms that, acting on different scales, are responsible for acceleration to high energies of a small fraction of the particles in the coronal plasma.

Numerical simulations in the different boxes are performed using the same realization of the turbulent magnetic, velocity, and electric fields obtained from the direct solution of RMHD equations. However, moving from the smallest to the biggest box, it is important to establish a physically reasonable scaling of these fields to avoid anomalous behavior. With this intention, we assume that the magnetic field B and the plasma velocity u are statistically similar in the different boxes, i.e., their spectra and their moments of the distribution function are the same in the different boxes. This essentially means that the statistics are scale-invariant in the inertial range. The only difference is that the same fields are rescaled and mapped, i.e., squeezed or expanded, into the new box with different dimensions. Indeed, going from one box to another, the quantity that is actually changing is the size L of the simulation box. Thus, we can use the same realization of the plasma velocity field of the RMHD simulation that is not dependent on the length L. Analogously, the magnetic field, written in terms of the Alfvén velocity, remains invariant with changing the characteristic length L. A different discussion must be done for the Ohmic component of the electric field E = ηj. From the dimensional analysis of this term, both the resistivity [η] = [VL] and the current density [j] = [B]/[L] are dependent on L. At smaller scales, the current density j becomes stronger, while the resistivity decreases. Thus, we need to rescale them separately so that their product, i.e., the Ohmic electric field, remains invariant (assumes the same range of values and has the same spectra) in the different simulations. We proceed as follows. First, we assume that in each box the grid length Δx is of the same order of the Kolmogorov microscale or dissipation scale l_d, determined by the condition that the dissipation rate equals the nonlinear energy transfer rate, τ_diss(l_d) = τ_NL(l_d) (Frisch 1995). For a given box, with characteristic length L and velocity u, this condition reads l_d/η² = L/u(l_d/L)^2/3, so at the length l_d, the resistivity is given by $\eta = (l_d^{4/3}/L^{1/3})u$, and η decreases with increasing L. The current j evaluated at the length l_d, where the strongest current sheets are present, is of the order j = B/l_d ∼ B/Δx. Considering the biggest and the smallest box in the models of grid size Δx_B and Δx_S, respectively (the subscripts B and S are associated with the biggest and the smallest box), the condition Δx_B ≫ Δx_S ensures that j_S ≫ j_B. Thus, at smaller scales, the current is stronger as it is expected to be. The requirement that the Ohmic electric field E = ηj spans the same range of values in the different boxes, i.e., E_S = E_B, leads to the condition

$$\begin{equation} \left(\frac{l_d}{L} \right) \sim \left(\frac{\Delta x}{L} \right) = {\rm const}. \end{equation} \tag{ 9 }$$

This is equivalent to assuming a constant magnetic Reynolds number R_m = uL/ηc² in the different boxes.

As the parameter β = λ_c/d_i in Equation (6) is the ratio between the correlation length and the ion inertial length (the largest and smallest physical scales in the system), we ensure that the particle in each box has knowledge of the real large-scale correlation length of the fields keeping the value of parameter β constant in all computational boxes and equal to its realistic value 5 × 10⁵. In all of the simulations, we maintain the same aspect ratio between the perpendicular and the parallel direction l_∥ = 10 l_⊥, and periodic particle reinjection boundary conditions are used in all directions. One of the conditions respected by the multi-box method is the requirement that the particle gyroradii are not so large in the rescaled boxes that a particle can gyrate at a box size in the perpendicular directions. Thus, avoidance of box-size effects due to periodicity is highly desirable. The RMHD fields in the parallel direction are not periodic while the particles that cross the z-boundary are reinjected at the other boundary, thus they experience a small discontinuity in the physical fields. However, this approximation is not expected to have a significant impact on the statistics of the particle dynamics. In fact, the current sheets are strongly elongated along z, and, following the reduced MHD scaling, the magnetic and velocity fields also change their values more slowly along this direction. Thus, the fields and structures experienced by the particle during its gyromotion do not change dramatically during the crossing of the box boundary in the parallel direction.

The correlation length of photospheric convective motions is ∼1000 km, the biggest box in the perpendicular direction is extended for 4 correlation lengths L_B = l_⊥ = 4000 km, with 256 × 256 grid points in the x–y plane and Δx = 3.9 × 10⁻³l_⊥ ≫ d_i (Figure 1). The perpendicular length of the smallest box measures L_S = L₁ = 86d_i = 4.87 × 10⁻⁵l_⊥, with a grid of 256 × 256 in the x–y plane and Δx₁ = d_i/3 = 1.9 × 10⁻⁷l_⊥, smaller than the minimum gyroradius. Considering the initial condition for the particles, five boxes are needed for particles to travel the entire length of the coronal loop.

Ten thousand particles are injected in the smallest simulation box (hereafter BOX1) with random position and the same thermal velocity $v_{\rm th} = 0.6\, v_{_{\!A}}$, corresponding to the coronal temperature T = 10⁶ K. In spherical coordinates, with the polar axis along the z-direction parallel to the mean magnetic field of strength B₀, the particle velocity components are:

$$\begin{equation} v_x = v\sin \theta \cos \phi \quad v_y = v\sin \theta \sin \phi \quad v_z = v\cos \theta. \end{equation} \tag{ 10 }$$

Particles' initial velocities are randomly distributed in the gyrophase ϕ and pitch angle θ. Time is advanced with a fourth-order Runge–Kutta integration method with an adaptive time-step (Press et al. 1992). The values of the fields at each particle position are obtained by linear interpolation in space from the grid of the RMHD simulation. The trajectories are calculated using a quasi-static approach in which magnetic and electric fields do not change with time.

The simulation in BOX1 is stopped after $5\times 10^5 \tau _g = 0.15 t_{_{\!A}}$, where τ_g = 2π/Ω is the mean particle gyroperiod (computed in the mean field). At this time, we select the high-energy (H-E) particle population, about ∼10% of the total BOX1 population, and inject them randomly in the next larger simulation box, BOX2. The gyroradius of the particles that are transferred is >10⁻⁶l_⊥ at the end of the BOX1 simulation. BOX2 has a grid scale of the order of their gyroradius (Δx₂ = 10⁻⁶l_⊥) and a perpendicular length L₂ = 2.56 × 10⁻⁴l_⊥. Analogously, the simulation in BOX2 is stopped at $t = t_{_{\!A}} \sim 3\times 10^6 \tau _g$, and a new high-energy population (with ρ ⩾ 10⁻⁵l_⊥) is selected and randomly injected in the third box (BOX3), which has a grid scale Δx₃ = 10⁻⁵l_⊥ and a perpendicular length L₃ = 2.56 × 10⁻³l_⊥. The BOX3 simulation is stopped at $t = 6 t_{_{\!A}} \sim 2\times 10^7 \tau _g$. In each box, the condition (9) is satisfied. The characteristic parameters used in the different simulations (box length L, grid size Δx, runtime t_F, average gyroradius $\bar{\rho }$) are summarized in Table 1. The grid scale in the second box is about five times larger than in the first box, while maintaining the same number of grid points. Furthermore, the grid scale in the third box is ten times larger than in the second box. This essentially means that larger scale structures are present in larger boxes. These structures are required for particles with a larger gyroradius to have a natural interaction with the fluctuations, including resonances.

Table 1. Values of the Parameters Used in the Three Simulations

BOX	L[l⊥]	Δx[l⊥]	$t_F[ t_{_{\!A}}]$	$\bar{\rho }[ l_\perp ]$
1	L1 = 4.87 × 10−5	1.9 × 10−7	0.15	10−7
2	L2 = 2.56 × 10−4	10−6	1	10−6
3	L3 = 2.56 × 10−3	10−5	6	10−5

Note. Box orthogonal length L, grid size Δx = L/256, runtime t_F, and average gyroradius $\bar{\rho }$, in units of the perpendicular cross section length scale of the physical loop l_⊥ = 4000 km and Alfvén crossing time t_A = l_∥/c_A (l_∥ = 10l_⊥).

Download table as:  ASCIITypeset image

At the end of the third simulation, particles have reached energies of 0.1 GeV and a nonrelativistic approach could not be suitable for higher energy. Moreover, the absence in our model of an escape mechanism and the lack of self-consistency intrinsic in test-particle simulations could become an important limitation too, especially in the presence of very high energy particles. Thus, we decided to stop this study with the simulation performed in BOX3 (please also note that L_B ≠ L₃).

In the next section, for each box, we show the results corresponding to the total particle distribution and separately for the associated high-energy (H-E) population.

We begin by showing diagnostics that characterize the energization of the particle populations in time. Figure 2 shows the time history of the average particle energy E in eV, for the entire simulation spanning BOX1, BOX2, and BOX3. Figure 3 shows the probability distribution functions (PDFs) of protons gyroradii, $\rho = v_\perp ^2/|{\bf B}|$, at different stages of the evolution. The PDF in BOX2 and in BOX3 are normalized to $1/N_1^*$ and $1/N_2^*$, where $N_1^*$ and $N_2^*$ are the percentage of particles selected at the end of the previous simulation and reinjected randomly in BOX2 and BOX3, respectively. We want to caution the reader at this point that while we present results in physical units, our goal is to study the physical mechanisms and not to reproduce any real data sets. There are numerous reasons that our simulation data would be expected to differ from space and astrophysical data, for example, the lack of consideration of mechanisms of escape, and the lack of self-consistency (test-particles), as well as other features dictated by the numerics.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The initial energy distribution is sharply peaked (delta function) at 89 eV. At the end of the simulation in BOX1, after $0.13 t_{_{\!A}}$, the energy of the H-E particles (dashed black line) is more than one order of magnitude larger than the energy of the total population. When simulation starts, ten thousand particles are injected in the box; when it is stopped, only 10% of the overall population gain sufficient energy for their gyroradii to be an order of magnitude bigger than the grid scale in BOX1, and therefore comparable with the grid size in BOX2, so that the particle can sample the structure of the fields in the second box. Indeed the H-E particles, increasing their energy, gyrate with a larger gyroradius and are allowed to sample and interact with different regions of the field. In contrast, particles that do not gain a sufficiently large amount of energy continue to gyrate with a gyroradius that is not large enough to explore different field structures, as if they were "confined" to the same range of scales and the same region of space.

This idea is reinforced examining the trajectories performed by the particles in BOX1, shown in the plot on the left of Figure 4. Three different trajectories are highlighted in the picture, with changing colors corresponding to increasing proton speed. Particles are observed to accelerate most when their trajectories are aligned to the mean axial field, B₀. For instance, two of them gain energy traveling a short distance in the axial direction and then start to rattle around trapped at constant z. On the other hand, the third particle with the trajectory more elongated along z rapidly ascends, gaining energy. This elevator-like path is associated with sheet-like structures of electric current density that are extended in the axial direction and provide channels for particle motion. Because at this early stage, the gyroradius is comparable or smaller than the ion inertial scale d_i (ρ_max = 0.6d_i at the injection, see black dotted line in Figure 3), when a proton encounters a current channel, it is likely to remain inside the channel over several gyroperiods.

Zoom In Zoom Out Reset image size

In this scenario, protons that find strong currents will reach larger energy mostly in the parallel direction, while protons that encounter only average or small currents will only acquire a moderate energy. It is also important to notice at this point that in BOX1, the largest difference between the averaged energy of the total distribution and of the H-E population is observed. This suggests that, although the energies reached at the end of the first stage are not extremely high, the H-E population is selected at the beginning of the evolution, namely for $t \ll t_{_{\!A}}$, by the effect of the Ohmic electric field acting in the direction parallel to the guide field B₀, along which the current density structures are preferentially aligned. Note that near-alignment of parallel current channels with a strong imposed magnetic field is expected for MHD flows that are mainly incompressive in nature (see, e.g., Oughton et al. 1994).

In the second box, we inject 1000 particles and follow their dynamics. These particles are taken from the high-energy tail of the distribution at the end of the simulation in BOX1, as shown by the overlapping of the gyroradius PDFs at the beginning of the simulation in BOX 2 (Figure 3). The energy difference between the H-E population and the total distribution (Figure 2) is substantially reduced, reinforcing the idea that the effects acting at the smallest scale are responsible for the selection of the high-energy particles.

The analysis of the particle trajectories in BOX2 (not shown) in comparison to BOX1 indicates that much fewer particles travel almost undisturbed along the axial direction. Indeed, if at the injection, few particles gyrate with a gyroradius smaller or comparable with the current sheet thickness, then at later times, (Figure 3) the probability distribution function of ρ moves toward larger values, i.e., ρ ≫ d_i. As a consequence, protons do not remain in the current channels, and they do not see a coherent parallel current, which could produce a net energization effect in the parallel direction, as happened previously in BOX1. Nonetheless, a net increase of the average energy is observed.

The same considerations enter the analysis of the dynamics in BOX3. It appears that all of the particles undergo the same acceleration mechanism here as in BOX2, reaching almost 0.1 GeV (green line in Figure 2). The motion in BOX3, shown in Figure 4, cannot be associated with the presence of strong current channels. Indeed, at the end of the third simulation, all of the particles have a gyroradius ρ ∼ 10⁻⁴ l_⊥ (Figure 3), greater than the thickness of the strong current sheets. Moreover, it is not possible to differentiate between two separate populations of particles.

A closer look at the H-E particles is required for a better understanding of their dynamics. In Figure 5, we show the position of the H-E particles in BOX1 and BOX2, together with contour plots of the absolute value of the current averaged in z, 〈|j|〉_z. This picture shows qualitatively that two different acceleration mechanisms are acting during different times of the particles' evolution. While the positions of the H-E particles are more or less equally distributed within the simulation BOX2, (with however some "voids" in evidence), the clusters of particles close to the current sheets structures in BOX1 indicates energization due to the Ohmic electric field within the z-aligned current channels. Eventually, some particles pitch angle scatter out of the channel, or, after attaining higher energy, leave the current channel because their gyroradius increases. This explains the presence of energetic particles in BOX1 further from the strongest current sheets.

Zoom In Zoom Out Reset image size

Another way to study particle motions is to examine pitch angle distributions. In Figure 6, the PDFs of the cosine of pitch angle α = v_∥/|v| in the three different boxes are shown. In BOX1 and in BOX2, distributions of α are shown for both the H-E population (dark blue) and the total population (light blue). As pointed out before, in the third box, it is not possible to pick out a high-energy distribution because all of the particles reach more or less the same final energy. Thus, the PDF in BOX3 (right plot) has been computed considering all the particles in the simulation. The H-E particles in BOX1 are mainly aligned to the guide field, α = ±1, reflecting the fact that the accelerating field is in the z direction. However, in BOX2, the PDF of pitch angle of the H-E distribution moves away from α = ±1 and is peaked around α = 0, indicating that the acceleration process is taking place mainly in the direction perpendicular to B₀. The pitch angle distributions support the picture of a two-stage acceleration process, which we further examine below.

Zoom In Zoom Out Reset image size

4.1. On the Nature of the Perpendicular Electric Field

We now analyze the different terms in the equation of motion to qualitatively explain the observed dynamics of the protons. In this analysis, we follow the procedure described in Dmitruk et al. (2004). For a general MHD representation of the plasma flows and fields, the equation for the parallel (v_∥) and the transverse (v_⊥) proton velocities are given by:

$$\begin{equation} \frac{d{\bf v}_\parallel }{dt} = \beta [ ({\bf v}_\perp - {\bf u}_\perp) \times {\bf b}_\perp + \eta {\bf j}_\parallel ], \end{equation} \tag{ 11 }$$

$$\begin{equation} \frac{d{\bf v}_\perp }{dt} = \beta [ ({\bf v}_\perp - {\bf u}_\perp) \times {\bf B}_0 + ({\bf v}_\parallel - {\bf u}_\parallel) \times {\bf b}_\perp + \eta {\bf j}_\perp], \end{equation} \tag{ 12 }$$

where u_⊥ and b_⊥ are the components of the plasma velocity and magnetic field perpendicular to the mean field B₀, and u_∥ is the parallel plasma velocity. If we specialize these equations to RMHD fields, such as those used in the present study, several terms in Equations (11) and (12) vanish identically because u_∥ ≡ 0 and j_⊥ ≡ 0 in RMHD. Moreover, given the ordering of RMHD fields, b_⊥ is small (b_⊥ < B₀).

The terms in Equation (11) all remain in RMHD, but in this highly anisotropic turbulence regime the inductive term, u_⊥ × b_⊥, will be small compared to ηj_∥ near the current sheets, and it will be negligible elsewhere as in the low corona velocity fluctuations concentrated around current sheets (Rappazzo & Velli 2011). Thus, the Ohmic electric field is the dominant term in the equation for the parallel velocity, and its efficiency is predominant at smaller scales, as we showed above.

After setting u_∥ = j_⊥ ≡ 0 in Equation (12), of the remaining terms, it is the first on the right-hand side that can effectively accelerate particles. It gives rise to a drift motion with velocity u_⊥ plus a gyromotion around the magnetic field B₀. For particles temporarily confined within a current channel, the parallel Ohmic field ηj_∥ will be dominant. However, as soon as their gyroradii become comparable with the current sheet thickness d_i, protons will see a transverse MHD-induced electric field E_⊥ = −(u_⊥ × B₀e_z)/c, that varies on scales comparable to their gyroradius. If kicks of this electric field are in phase with the proton velocity, this results in a net increase of the perpendicular energy. This may be viewed as a betatron-like resonance process (for a description of the betatron mechanism see Swann 1933). Thus, variations of the transverse electric field are produced by variations in the plasma velocity u_⊥. The plasma velocity is concentrated around the current sheets that are indeed embedded in a vorticity quadruple. Therefore, particle acceleration associated with the perpendicular electric field component will act more strongly in the region near current sheets, including the interior of flux tubes bordering the current sheets.

In order to verify the effectiveness of this mechanism, we computed the work done by the electric field on particles moving with velocity v, namely, (v · E). This calculation is done in BOX3 where the betatron-resonant process is thought to dominate. Figure 7(a) shows the power spectrum of (v · E) averaged over 100 particles. The parallel and the perpendicular terms in the product, v_∥E_∥ and v_⊥E_⊥, are both shown, and the frequencies are expressed in term of the characteristic ion gyrofrequency Ω_i. For this analysis, particle trajectories have been integrated over 10⁴ gyroperiods, saving the required data every Δt = 0.1τ_g. It is clear that the process involves a substantial degree of resonance, as strong peaks are found for ω = Ω_i and for the successive harmonics, from n = 2 to n = 5. These features indicate that a strong coupling between the plasma and the particles is occurring at the particle gyrofrequency.

Zoom In Zoom Out Reset image size

To verify that protons have been effectively accelerated in the direction perpendicular to B₀, we computed the PDF of q = (v · E)/|v| together with the PDFs of the parallel q_∥ = (v_∥E_∥)/|v| and perpendicular q_⊥ = (v_⊥E_⊥)/|v| terms. These quantities are shown in Figure 7(b). The PDFs have been computed considering all the values along the trajectories of 100 particles. The PDF of the total q and the PDF of the perpendicular term q_⊥ are nearly indistinguishable. The PDF of q_∥ is narrower, showing that the energization is larger in the perpendicular direction. As expected, q, q_∥ and q_⊥ display both positive and negative values. Indeed, throughout their evolution, particles can be accelerated but they can slow down as well. However, the average values of the three distributions are positive, confirming that a net acceleration is occurring. The average values and their standard deviations are ${\bar{q}} = 1.2 \times 10^{-6}$ and σ = 2.7 × 10⁻³, ${\bar{q}_\parallel } = 8.3 \times 10^{-7}$ and σ_∥ = 3.8 × 10⁻⁴, ${\bar{q}_\perp } = 4.7 \times 10^{-7}$ and σ_⊥ = 2.7 × 10⁻³.

Moreover, in Figure 7, the lower frequency part of the total power spectrum (panel (a)), i.e., for frequencies lower than 0.1Ω_i, is dominated by the parallel terms, while for higher frequencies, the perpendicular component dictates the shape of the spectrum. The ∼1/ω² behavior of the parallel frequency spectrum suggests second order Fermi acceleration of a standard type, in which scattering from counter-propagating fluctuations causes a momentum diffusion. However, it is also clear that this is not the dominant acceleration mechanism at this stage. Perpendicular acceleration accounts for about two orders of magnitude more power (panel (d)) and consists of systematically greater magnitudes of acceleration (panel (b)). We may then conclude that the system remains complex, with several mechanisms involved at this stage, but according to several measures shown here, the betatron mechanism may be viewed as the strongest mechanism that we have identified. It is also clear that not only a separation of length scales is involved, but also a separation of times scales.

As we pointed out above, the energy increase in the perpendicular direction arises if variations of the plasma transverse velocity field u_⊥ are in phase with the gyromotion of the particle. From this perspective, the process can be thought of as resonant. However, it is important to emphasize that the mechanism described should be considered as a first-order process being caused by a convective electric field E_⊥ = −(u_⊥ × B₀e_z)/c. It is also important to point out that, although it represents the main acceleration mechanism in the second and third box, the coupling between the convective electric field and the particles can take place as soon as the particle gyroradius becomes bigger than the ion inertial scale d_i.

More investigations have been conducted to better characterize the energization that occurs during the different stages. Figure 8 shows the time history of the ratio of average perpendicular energy and average parallel energy E_⊥/E_∥ in the different boxes for the total population and for the H-E particles. It is clear from the evolution of E_⊥/E_∥ in BOX1 that the H-E particles are preferentially accelerated in the direction parallel to B₀. When particle energy is sufficiently high, pitch angle scattering can transfer momentum from the parallel to the perpendicular direction, thus explaining the increase of the ratio at the end of the BOX1 simulation. Subsequently, the ratio E_⊥/E_∥ rapidly attains values near 2, i.e., the distribution tends to isotropize because of pitch angle scattering, as shown above by the PDFs of α for the H-E population in BOX3 and BOX2. However, in BOX2, it is still possible to observe strong peaks of E_⊥/E_∥ > 2, testifying to a strong acceleration in the perpendicular direction.

Zoom In Zoom Out Reset image size

In this paper, we investigated particle acceleration in turbulent electromagnetic fields embedded in a strong guide magnetic field, an important issue in astrophysical, space, and laboratory plasmas. We performed test particle simulations for protons in turbulent fields obtained from direct numerical solutions of RMHD turbulence (for a detailed description of the RMHD simulation, see Rappazzo et al. 2008). Although the test particle approach is not self-consistent, it provides a useful bridge between the large-scale macroscopic MHD and a more realistic particle kinetic physics description. The geometry of the system and the typical parameters used in the simulations are well-suited for the description of particle acceleration in coronal loops.

A common problem in the numerical description of plasma and/or particle dynamics in space plasmas is the involvement of a huge range of length and time scales. In the case of coronal loops, the large scale motion occurs on a distance comparable with the loop perpendicular length (thousands kilometers) but the acceleration process takes place at different stages and/or locations during the evolution of the system and may require the participation of several mechanisms. Therefore, it can be considered as a multistage process, and the numerical study may require that the various length and time scales should be followed with different time and spatial resolutions. In the study of this kind of system, the acceleration may be associated with strong and sudden changes, including changes of magnetic topology near reconnection sites as well as strong discontinuities. Consequently, the use of the guiding center equations may be inappropriate because non-adiabatic behavior may develop on local scales.

Motivated by these reasons, and to maintain the possibility of gyro-resonant interactions at all phases of the acceleration process, we have developed a multi-box, multiscale technique to solve the nonrelativistic equation of motion of a particle moving by the action of electromagnetic fields (see Section 3). This technique allowed us to resolve an extended range of scales present in the system, namely, from the ion inertial scale d_i, of the order of 1 m, up to macroscopic scales of the order of 10 km (1/100th of the correlation scale, the photospheric convection scale ∼10³ km), resolving the particle gyroradius at the same time. An important feature of the multi-box method is that the scaling of the MHD fields is controlled so that the magnetic field, plasma velocity, and electric field take on similar values in the rescaled boxes as their size varies from the initial smallest box to the largest final box. This amounts to changing length scales and resistivities in the sequence of boxes to maintain the required statistical continuity.

Employing the multi-box strategy, we performed a numerical experiment that shows how acceleration takes place typically in a two-stage process. First, a small fraction of initially loaded low-energy particles are entrained in elongated current sheets long enough to gain substantial energy in response to the parallel electric field in the current sheet. High parallel speeds are attained in this stage. Gradually, these particles escape the current channels due to pitch angle scattering and due to their larger gyroradius after energization. In the second stage, particles encounter nonuniform in-plane electric fields that enable a resonant betatron acceleration for those particles that remain entrained in these regions for a suitable amount of time. In this case, the acceleration is mainly perpendicular. However, once again, pitch angle scattering tends to isotropize the particles, which also then facilitates their escape from the accelerating region.

A point of interest is that the two-stage acceleration process described here is fully consistent with important features of test particle acceleration in turbulence that have previously been reported. For example, the initial parallel acceleration of test particles in current channels associated with 2D reconnection was reported some time ago by Ambrosiano et al. (1988). There too, it was observed that the energization is limited by the time of temporary trapping of particles near the acceleration region. The role of pitch angle scattering in disrupting the parallel acceleration process was also described. Furthermore, it is likely that secondary islands near or in current sheets may enhance the ability of sheet-like structures to entrain and accelerate particles (Ambrosiano et al. 1988; Drake et al. 2006). We have not attempted to describe the effect of secondary islands here.

It is also noteworthy that the parallel and perpendicular stages of acceleration that we described were anticipated, although as separate effects, in the study by Dmitruk et al. (2004). In that case, using a single simulation box, it was shown that at early times, protons experience perpendicular acceleration while electrons experience parallel acceleration. The consistency of that picture with the present one is established by noting that in that study, particles were initialized with zero speed but quickly acquired speeds comparable to the Alfvén speed. However, the current sheets have a thickness comparable to the ion inertial scale, so the protons almost immediately encounter the edges of the current sheets and therefore enter the perpendicular "stage two" process described here. Given their much lower inertia, electrons have very small gyroradii when moving at the Alfvén speed, and so can experience what we presently call the stage one parallel acceleration process, gaining substantial energy in a short time. The multi-box method allows the computation to follow both stages one and two for a single particle type, while maintaining resonant couplings that are central to spatial transport and therefore contributing to entrainment and escape processes.

For the particular example appropriate to coronal loops that we showed here, by the end of the third box in the multi-box sequence, particles attained energies up to 0.1 Gev. In order to cover a realistic large parallel coronal loop length, another two simulations in larger boxes would need to be performed. However, at the end of the third simulations, particles have reached energies of 0.1 GeV and a nonrelativistic approach could not be not suitable for higher energy. Continuing in this way, higher energies are expected at the end of the largest scale simulation, and an "escape" mechanism mimicking the collision of particles with the denser chromospheric layer (e.g., Gontikakis et al. 2013) should be included. Furthermore, it must be pointed out that the lack of self-consistency intrinsic in test-particle simulations could at some point become an important limitation, especially for the highest energy particles.

Although a refinement of the multi-box technique may also be required, the results shown in this paper indicate a promising way of understanding energetic particles in astrophysical, space, and laboratory environments and studying some basic issues in the complex topic of bridging the MHD and kinetic descriptions of a plasma.

This research is supported in part by NASA Heliophysics Theory program NNX11AJ44G, NSF Solar Terrestrial and SHINE programs AGS-1063439 & AGS-1156094, NASA Magnetosphere Multiscale mission through the Theory and Modeling Team, the Solar Probe Plus Project through the ISIS Theory team, and by EU Marie Curie project "Turboplasmas" at the University of Calabria.