CORONAL ELECTRON CONFINEMENT BY DOUBLE LAYERS

Coronal electrons can be accelerated to very high energies, over 100 keV (Lin et al. 2003; Krucker et al. 2010, 2007; Tomczak 2009), during flares and coronal mass ejections. Acceleration is believed to occur high in the corona (Fletcher et al. 2011). Hard and soft X-rays emitted by the accelerated electrons are observed in both the corona and chromosphere. Emission from the chromosphere is usually brighter because of the higher density, which allows for more efficient interaction with the X-ray-producing energetic electrons. Improved detection sensitivity, however, has led to more measurements of X-ray sources in the corona. These sources are situated at or distinctly above the top of the soft X-ray flaring loops. A comprehensive analysis suggests that above-the-looptop sources are a common feature of all flares (Petrosian et al. 2002).

The transport of energetic electrons from the corona to the chromosphere is crucial to understanding energy release in flares. Transport effects can modify the energy distribution of the propagating electrons and hence the observed X-ray spectra, affecting the interpretation of acceleration models. There is evidence for both free-propagation and confinement of energetic electrons. Evidence for free-propagation of energetic electrons, i.e., no interaction with the ambient plasma, was reported in previous time-of-flight measurements of hard X-ray emission (Aschwanden et al. 1995, 1996; Aschwanden & Schwartz 1995). A systematic time delay between lower-energy hard X-rays with respect to the higher energies was measured from a large sample (over 600) of solar flares. The time-of-flight delays indicate that lower-energy electrons arrive at the chromosphere after those with higher energy as they freely stream down the flare loop from the corona. On the other hand, observations of above-the-looptop sources reveal that the lifetime of the energetic electrons is two orders of magnitude longer than their free-streaming transit time across the source (Masuda et al. 1994; Krucker et al. 2010, 2007). This requires confinement of the electrons in the source region. Recently, a systematic study of solar flares with both looptop and footpoint emissions found that the number of electrons required to explain observations is two to eight times higher at the looptop than at the footpoint (Simões & Kontar 2013). This suggests electron accumulation at the looptop. Another systematic study reported that the difference in spectral index between non-thermal coronal and footpoint emission in some events is considerably greater than two, which is expected for the simple case that the same electron population, which produces looptop emission through thin target bremsstrahlung because of the lower density in the corona, free-streams to the footpoints to emit through thick target bremsstrahlung due to the higher density in the chromosphere (Battaglia & Benz 2006). This requires a filter effect in the propagation preferentially reducing the distribution at lower energies. Such a filter can be an electric field. An important question on electron transport is why some electrons appear to freely propagate while some are filtered or even trapped. The physics of this subject, however, remains poorly understood.

The transport of electron heat flux has previously been modeled as anomalous conduction (Manheimer 1977; Tsytovich 1971) in which anomalous resistivity arises from electron scattering by turbulent wave fields. Turbulence is excited by instabilities that result from the interaction between energetic electrons and ambient plasma. In a similar model, it was suggested that ion acoustic turbulence driven by a return current (RC) in response to hot electron streaming and an induced polarization electric field would drive a potential that would suppress electron transport (Levin & Melnikov 1993). An anomalous conduction front that moved along the flare loop at the head of an expanding hot electron source was considered as a means to confine hot electrons for the production of hard X-rays (Smith & Lilliequist 1979). The model, based on a one-dimensional (1D) one-fluid code, however, did not resolve ion inertial lengths, let alone capture processes occurring at electron scales. Later, 1D electrostatic particle-in-cell (PIC) simulation that resolved the shortest electron scale, the Debye length, did not reveal a conduction front (McKean et al. 1990). Recently, the existence of a thermal front, literally defined as a region that links plasmas in thermal nonequilibrium and sustains the temperature difference for longer than the electron free-streaming time, was studied in 1D electrostatic Vlasov simulations (Arber & Melnikov 2009). The formation of a temperature difference propagating at a speed comparable to the ion acoustic speed was observed and the behavior was identified as a conduction front. The physics of the responsible mechanism was not identified or investigated. We suggest, on the basis of the simulations and analysis in this paper, that the propagating front is an ion acoustic shock and the temperature difference to be a result of shock heating due to its extremely sharp transition. Other observed non-propagating temperature jumps that were associated with a potential jump were likely double layers (DLs), as mentioned by the authors.

More recently, the transport of coronal energetic electrons was studied by electromagnetic two-dimensional (2D) PIC simulations (Li et al. 2012, hereafter LDS). It was shown that transport suppression began as the electrons propagated away from the acceleration site. The suppression was caused by the formation of a DL and associated potential barrier that reflected electrons back to the source region. The DL reduced the electron heat flux by nearly 50%.

A DL is a localized region that sustains a potential drop in collisionless plasmas (Block 1978; Raadu & Rasmussen 1988). The potential drop comes from a large-amplitude electrostatic electric field sandwiched between two adjacent layers of equal and opposite charges. The structure is globally neutral, but quasi-neutrality is locally violated within it, which occurs at scales of ∼10 Debye lengths λ_De. An ideal DL is a monopolar electric field, but more generally, a DL can be bipolar. Therefore, instead of a monotonic drop in the potential ϕ, there can be dips or bumps at the low or high potential sides, but an overall potential drop ϕ_DL across the structure. ϕ_DL is the measure of the strength of a DL. A DL reflects particles if ϕ_DL is greater than their kinetic energy. In LDS, it is the reflection of energetic electrons by a DL that reduces their transport from the source.

Numerical simulations with a kinetic model have been used to study DL dynamics since they can both resolve kinetic scales down to λ_De, where DLs occur, and to probe nonlinear phenomena. Methods used to generate DLs include a current and an applied potential. In LDS and this work, the DL forms as a result of strong currents driven by an imposed field-aligned temperature jump in the initial state. Earlier 1D PIC simulations of systems driven by a subthermal electron current showed that ion trapping in the potential dip of a (bipolar) DL was associated with decay of a DL (Barnes et al. 1985). The trapping of ions slowed down the movement of the DL structure as the trapped ions added an inertial drag to it (Chanteur et al. 1983). It was argued that such slowing down led to the decay of the DL (Barnes et al. 1985). On the other hand, 1D PIC simulations with an applied potential reported that a DL decayed as a result of the propagation of solitons, which developed from the DL structure, across the DL width (Sato & Okuda 1981). Only very limited work exists on the saturation of DLs.

The degree of transport suppression depends on the maximum strength of the DL at saturation. In LDS, the generation mechanism of a DL was studied and identified as the Buneman instability involving the ions and RC electrons. LDS showed that the potential increased with increasing hot electron temperature, but the functional dependence on this temperature was not studied. In this work, we investigate the saturation mechanism of the DL with simulations and analytic modeling. We carry out a series of simulations with different initial values of the hot electron temperature and show that the DL potential at saturation scales linearly with the hot electron temperature. The DL saturates when its potential jump is large enough to accelerate ions above the local sound speed. The result is a parallel ion acoustic shock that stabilizes the Buneman instability and therefore saturates the DL. We demonstrate that the shock formation criterion predicts a maximum DL strength that is proportional to the hot electron temperature in agreement with simulations. At its maximum strength, the DL is observed to reflect, and hence contain, a significant fraction of electrons in the source.

In addition, we observe that an anisotropy with T_⊥ > T_∥ develops in the hot electrons that pass through a DL. Such an anisotropic (pancake) distribution favors further hot electron trapping in the presence of a magnetic mirror. Since the magnetic geometry of flaring loops resembles a mirror configuration, we explore combining a magnetic mirror with the DL. We demonstrate that the combination significantly enhances confinement.

An outline of the work is as follows: the setup of the simulations and parameters used for flare settings are described in Section 2; in Section 3, results from the simulations and the shock model are presented, and we provide evidence for the model and determine the saturation amplitude of a DL; in Section 4, we combine the DL with a magnetic mirror and show that further suppression of transport is obtained; further discussion of our results is given in Section 5; and we summarize in Section 6.

Two-dimensional electromagnetic PIC simulations are carried out using the p3d code (Zeiler et al. 2002). The initial setup here is the same as in LDS. A rectangular domain (Figure 1) represents a symmetric local segment of a flare loop centered at the looptop where a very hot electron source is located. The long direction, x, is along the loop axis and an initial uniform background magnetic field B₀ is applied in that direction. The short direction, y, is mainly for data averaging to produce smoother data in x. The initial density n₀ is uniform. A third of the electrons, centered in the domain, have higher initial parallel temperature T_{h0, ∥} than the rest. It represents a pre-accelerated hot looptop source. Higher T_{h0, ∥} is used because parallel transport is of more interest since it dominates over perpendicular transport in the presence of a strong magnetic field, which is true for a coronal flare loop. Electrons are initially modeled by bi-Maxwellian distributions and the ions by a Maxwellian distribution. We consider it to be more natural that the hot population does not have a preferred direction of propagation, so the simulated energetic electrons are not beamed in the initial state. We note that soft X-ray counterparts of hard X-rays are observed in the high corona (Masuda et al. 1994; Krucker et al. 2010). Soft X-ray spectra are well fit by Maxwellian distributions (Masuda et al. 2000; Tsuneta et al. 1997) and hard X-ray spectra are usually fit to a combination of a thermal core and a non-thermal tail distribution. The use of Maxwellian distributions in our simulation setup is consistent with soft X-rays and the thermal core of hard X-rays, but may not directly apply to the combined distribution that produces hard X-rays; particularly, the fraction of escaping hot electrons derived in Section 3.3 may not be directly applied to that case.

Zoom In Zoom Out Reset image size

We perform three simulations with increasing T_{h0, ∥} = 0.5, 1, and 2. Temperatures are normalized to $m_ic_A^2$, where m_i is ion mass and c_A = B₀/(4πm_in₀)^1/2 is the Alfvén speed. All hot electrons have unity plasma beta β (ratio of plasma pressure to magnetic pressure) in the parallel direction. Using unity β for hot electrons is consistent with recent coronal flare observations (Krucker et al. 2010). β_{h0, ∥} is maintained at the same value for different T_{h0, ∥} by varying B₀, which is not expected to affect the results since the phenomenon observed here is dominantly electrostatic. The magnetic field is not significantly perturbed in the simulations. Outside of the hot electron region, the ambient (cold) electron temperature T_c0 is 0.1, as are both the perpendicular temperatures throughout the domain and the ion temperature. The ratios of hot to cold (parallel) electron temperatures are, therefore, T_{h0, ∥}/T_{c0, ∥} = 5, 10, and 20, respectively, for the three runs. This increasing contrast in temperature allows us to study how the strength of a DL scales to the real system that has over three orders of magnitude in separation of scales. We will drop "∥" in T_{0, ∥} from here on for simplicity.

The size of the simulations is L_x × L_y = 6553.6 × 25.6 $\lambda _{{\rm De}}^2$ with a grid size of 0.2 × 0.2 $\lambda _{{\rm De}}^2$ for T_h0 = 0.5. λ_De = v_th0/ω_pe is the Debye length, v_th0 = (2T_h0/m_e)^1/2 is the thermal speed based on the initial hot parallel electron temperature, and ω_pe = (4πn₀e²/m_e)^1/2 is the electron plasma frequency. For T_h0 = 2, λ_De is twice as big, so the grid size is 0.1 λ_De. In all cases, λ_De is well resolved. This resolution is for the stability of the code and is also important for capturing the structure of DLs that occur near the λ_De-scale. For T_h = 5 keV and n = 10⁹ cm⁻³, which are typical of coronal thermal X-ray sources, L_x ∼ 100m, which is of course far smaller than a realistic flare loop. Space and time are normalized to λ_De and $\omega _{pe}^{-1}$. A mass ratio m_e/m_i of 1/100 and speed of light c/c_A of 100 are used. The system is periodic in both directions. Because of the periodic boundaries, the simulations are evolved for less than the electron transit time of the domain at 1.5v_h0, so the majority of hot electrons will not reach a boundary during a run. The domain size is the same when measured in terms of electron inertial lengths d_e = c/ω_pe for all three T_h0. Hence, higher T_h0 runs with hotter electrons are evolved for shorter periods of time than the lowest T_h0 run. We note that our simulation domain is effectively 1D since the transverse direction is ⩽25 λ_De wide. The second dimension allows data averaging for noise reduction.

In addition, we perform an isotropic (T_{h0, ⊥} = T_{h0, ∥}) simulation using similar parameters as the highest temperature run to study how the DL changes the isotropy of the hot electron distribution as the electrons pass through the DL. This will be discussed in Section 4.

3.1. Scaling of DL Strength in Simulations

Figure 2 shows the time evolution of the DL strength in three simulations with increasing initial hot electron temperatures T_h0 = 0.5, 1, and 2. A similar cycle of DL growth followed by saturation is observed in all three runs. The growth mechanism was discussed in detail in LDS. As the contrast between hot and cold temperatures increases, the instability becomes stronger and develops more rapidly, leading to faster and further strengthening of the DL. For the same system size, i.e., for about the same distance traveled by the hot electrons, the maximum strength $e\phi _{{\rm DL}}^{{\rm max}}$ the DL reaches is higher in higher temperature runs. $e\phi _{{\rm DL}}^{{\rm max}}$ is plotted in Figure 3. A good linear fit (dashed line) of $e\phi _{{\rm DL}}^{{\rm max}}$ = 0.73 T_h0−0.11 with a coefficient of determination of 0.9998 is obtained. It indicates that the maximum strength of the DL scales linearly with the initial hot electron temperature. It also implies a threshold on the lowest possible T_h0 for DL formation. However, this threshold is not considered to be reliable as the regime of very low T_h0 is not explored and the behavior is unknown. We therefore take the fit to be approximately $\tilde{\phi }^{{\rm max}}\equiv e\phi _{{\rm DL}}^{{\rm max}}/T_{h0}\sim$ 0.73. Arber & Melnikov (2009) reported a proportionality of the potential drop to the final temperature drop across the thermal front from simulations based on two sets of T_h0's and T_c0's. The temperature drop could be close to T_h0 and hence, the proportionality to the temperature drop is similar to the linear scaling with T_h0 presented here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Shock Model

The physics occurring at the contact between the hot and cold electrons plays a crucial role in whether the hot electrons propagate out or are confined within the region. Initially, the hot electrons propagate by free-streaming along the magnetic field. Ions, being less mobile, lag behind. This creates a charge imbalance in the hot region, which draws in an electron beam from the cold electrons as an RC. The RC electrons drift relative to the ions, exciting an ion–electron streaming instability, identified as the Buneman instability (LDS). It is the driver of the DL. As the DL forms at the contact between the two populations, hot electrons with kinetic energy less than the DL potential eϕ_DL are reflected. The stronger the DL, the more reflected electrons and hence the fewer escaping hot electrons.

The maximum strength of a DL depends on its saturation mechanism. Saturation is observed to be linked to the formation of a shock in our simulations. In the following, we show data from the highest temperature run to demonstrate this. Figure 4 shows the time evolution of the electric field E_x from that run. Data are averaged in y. The bright white feature around ω_pet = 250–350 is the DL. It drifts to the left over time. A rightward propagating shock forms and starts to part from the DL at ω_pet ≳ 350, forming an intense negative electric field later in time. Around ω_pet ∼ 350 is when the DL starts to saturate (red curve plateauing in Figure 2). It is observed from all the simulations that whenever a shock forms, the DL stops strengthening (see, also, data from the lowest T_h0 run in Figure 3 from LDS). The shock is driven by ions accelerated to high velocities in the positive x direction, which will be further discussed later, so it propagates to the right and decouples from the DL.

Zoom In Zoom Out Reset image size

3.2.1. Shock Stabilizing the Instability

Shock formation saturates the DL growth by stabilizing the Buneman instability. In Figure 2, for example, the evolution of the DL strength eϕ_DL from the highest T_h0 run (red) indicates a growth phase before ω_pet = 350 and in general a saturation phase after that. DL saturation is evident from the decrease in E_x of the DL after ω_pet = 350 in Figure 4. In the following, we explain how the shock stabilizes the instability by reducing the speed of the RC beam that is the driver of the instability. Evidence is observed in the electron phase space. We show in Figure 5 the electron phase space at a time during the growth phase, at ω_pet = 250 (panel (a)), and during the saturation phase, at ω_pet = 600 (panel (b)). E_x (orange) is overlaid in panel (b) to show the location of the shock that has a large-amplitude negative peak in E_x. A vertical black line denotes the low potential side of the DL. To its left is the hot electron population (note the reflected hot electrons (cyan) in panel (a) with their velocities wrapping from positive to negative). To the right of the vertical black line is the RC beam. The speed of the RC beam during the saturation phase in panel (b) is |v_d| < 0.5v_th0, which is noticeably lower than that during the growth phase, which is |v_d| > 0.5v_th0 in panel (a). At ω_pet = 600, the shock is at x ∼ −2λ_De. Its large-amplitude negative leg significantly reduces the speed of the RC beam that is drifting in the negative x direction, and even reflects a significant number of RC electrons, i.e., from negative to positive velocities. The surviving beam that reaches the DL (vertical black line) has a much lower speed and is stable to the Buneman instability. Note that the reflection of the RC electrons at the DL due to trapping by the Buneman instability, which drives up the DL amplitude (see LDS), is absent at ω_pet = 600.

Zoom In Zoom Out Reset image size

3.2.2. Driver of the Shock and Maximum DL Strength

It is observed from the ion phase space that the shock is driven by ion acceleration to high velocities, beyond the ion sound speed, at the DL. Figure 6 shows the evolution of the ion phase space before and during the formation of the shock. The electric field E_x (orange) is overlaid in panel (a) to show the position of the DL, and the parallel ion pressure P_ix (magenta) in panels (d)–(f) to show the steepening of the shock. Velocities are normalized to the sound speed $c_{s0}\approx \sqrt{ T_{h0}/m_i}$ based on the initial hot electron temperature (recall that the simulations and the coronal energetic electron sources are in the limit where the electron temperature is much greater than the ion temperature). In panel (a), at ω_pet = 300, ions are strongly accelerated at the location of the DL, which is the large positive peak in E_x. Over time, they are further accelerated. Around ω_pet = 400 (panel (c)), some ions reach c_s0 and to their right, a jump in their velocity is developing. From ω_pet = 450 (panel (d)) and later, ion velocities at the DL surpass c_s0 with a velocity jump forming at the right of the "supersonic" ions. Velocities drop from supersonic to subsonic over a narrow transition less than the DL width (which is ∼10λ_De). P_ix also abruptly increases across the shock as expected. Thermalization of the ions is likely due to trapping in the large amplitude waves within the transition (Figure 6) as studied previously (Quest 1988). These observations lead us to conclude that the shock is driven by ion acceleration at the DL. The stronger the DL, the higher the ion velocity. Figure 7 shows the evolution of the maximum ion velocity measured from the ion phase space (solid), and ion velocity based on acceleration by the DL potential, v_{i, DL} = $\sqrt{2 e\phi _{{\rm DL}}/m_i}$ (dashed). v_{i, DL} resembles the shape of the red curve in Figure 2 since it is based on the DL potential. The solid and dashed curves are consistent with each other, and hence indicate that the accelerated ions are produced by the DL potential.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that by the time the DL saturates in the simulation, the local electron temperature T_h at the DL drops to ∼0.75 T_h0 (because a fraction of hot electrons escape from the source region), so the local sound speed c_s($\approx \sqrt{ T_h/m_i}$) is lower than c_s0. That is probably why the shock is generated even though the accelerated ion velocities only marginally exceed c_s0 in Figure 6(c).

Since the DL accelerates ions, which eventually drives a shock that in turn stabilizes the instability and saturates the DL growth, we can use the shock formation criterion, namely, ion acceleration to supersonic speed, to determine the maximum strength of the DL.

$$\begin{eqnarray} && v_i \sim \sqrt{\frac{2 e\phi _{{\rm DL}}}{m_i} } > c_s \sim \sqrt{\frac{T_h}{m_i} } \nonumber\\ && \quad \therefore \frac{e\phi _{{\rm DL}}}{T_{h}} > \frac{1}{2}. \end{eqnarray} \tag{ 1 }$$

Normalizing to T_h0 for comparison with simulations, we have:

$$\begin{equation} \frac{e\phi _{{\rm DL}}}{T_{h0}} > \frac{1}{2}\, \frac{T_h}{T_{h0}}. \end{equation} \tag{ 2 }$$

Approximating T_h ∼ 0.75 T_h0, it gives eϕ_DL/T_h0 ≳ 0.4. As soon as the DL reaches a strength of 0.4 T_h0, a shock forms, leading to DL saturation. This equation shows that the saturated DL strength scales linearly with the hot electron temperature, which is consistent with the result in Figure 3. The maximum DL strength based on this shock model is $\tilde{\phi }^{{\rm max}}= e\phi _{{\rm DL}}^{{\rm max}}/T_{h0}\sim$ 0.4. This value agrees with the measurement from our simulations within a factor of two.

We note that this prediction is independent of the ion-to-electron mass ratio. It is therefore applicable to realistic systems with real mass ratio.

3.3. Escaping Electron Density

To quantify the degree of confinement, we can compare either the reflected electron density in the hot source or the opposite, the escaping electron density, to the total density. In the following, we calculate the escaping hot electron density n_esp after passing through a potential barrier eϕ_DL. Given the initial hot electron distribution function f_e0, which is a Maxwellian distribution, and eϕ_DL, n_esp can be expressed analytically by integrating f_e0 over velocities of the escaping electrons. f_e0 is a function of both parallel and perpendicular velocities, so a three-dimensional (3D) volume integral ∫d³v needs to be performed. However, the DL is a parallel electric field and only reduces electron velocities in the parallel direction. The integral ∫f_e0d³v is therefore separable in v_∥ and v_⊥. The perpendicular contribution gives unity, so n_esp is reduced to a 1D integral over only the parallel velocities. For simplicity, we drop "∥" in "v_∥" in the calculation below. Note that it suffices to consider either side of the contact since the system under consideration is symmetric. We choose the right side where the escaping hot electrons move to the right and thus have positive velocities. We then have:

$$\begin{eqnarray} n_{{\rm esp}} &=& \int ^{\infty }_{b} \! f_{e0} \, dv; \quad b\equiv \sqrt{\frac{2e\phi _{{\rm DL}}}{m_e}} \nonumber \\ &=& \frac{n_0}{\sqrt{\pi }\,v_{th0}} \int ^{\infty }_{b} \! e^{-v^2/v_{th0}^2} \, dv \nonumber \\ &=& \frac{n_0}{2}\, {\rm erfc}(\tilde{\phi }^{1/2}); \quad \tilde{\phi }^{1/2}=\frac{b}{v_{th0}} , \end{eqnarray} \tag{ 3 }$$

where ${\rm erfc}(z) = ({2}/{\sqrt{\pi }}) \int ^{\infty }_{z}\,e^{-t^2}\,dt$ is the complementary error function. The total density of hot electrons moving to the right is n_tot = $\int ^{\infty }_{0} \! f_{e0} \, dv$ = n₀/2, so the ratio of escaping to total density is

$$\begin{equation} \tilde{n}_{{\rm esp}}\equiv \frac{n_{{\rm esp}}}{n_{{\rm tot}}}={\rm erfc}(\tilde{\phi }^{1/2}). \end{equation} \tag{ 4 }$$

As $\tilde{\phi }$ increases, $\tilde{n}_{{\rm esp}}$ decreases, i.e., fewer hot electrons escape. This tendency is also observed in the simulations. For example, at lower values of $\tilde{\phi }$ before the DL saturates, $\tilde{n}_{{\rm esp}}$ is measured to be higher than the value at saturation. Using $\tilde{\phi }\sim$ 0.7 obtained in Section 3.1, $\tilde{n}_{{\rm esp}}$ is predicted to be ∼0.24 with Equation (4). In the simulations, $\tilde{n}_{{\rm esp}}\sim$ 0.25n₀/n_tot = 0.5, which is within a factor of two from the prediction. This value of escaping density implies that about 50% of the total density is trapped in the source region. Thus, a substantial number of electrons are confined.

In the escaping hot electrons, the parallel velocities are reduced by the DL while the perpendicular velocities are not affected. This causes an anisotropy downstream from the DL. We indeed observe a pancake distribution, as discussed below.

It is easier to identify an anisotropy due to the DL when starting with an isotropic distribution, so a simulation with the same parameters as the highest temperature run (T_{h0, ∥} = 2) described in Section 2, except for using an isotropic distribution for the hot electrons, i.e., T_{h0, ⊥}=T_{h0, ∥}, and L_x being half of the original and L_y being twice the original, is performed. Figure 8 shows the electron distributions versus velocities parallel (v_ex) and perpendicular (v_ey) to the background magnetic field, near the end of the run, at ω_pet = 350. They are from locations (panel (a)) upstream, i.e., before the hot electrons encounter the DL, and (panel (b)) downstream, i.e, after some of them escape from the DL. They are sampled over a space of 50 λ_De in x, several times the width of the DL. A lower-velocity region is masked out since it also contains the cold RC electron beam, which has a much higher peak value than the hot electrons and will therefore overshadow the latter. In panel (a), the distribution takes on a circular shape. Hot electrons are nearly isotropic upstream. After passing through the potential barrier of the DL, the parallel velocity is reduced while the perpendicular velocity is unchanged. This gives rise to an anisotropy downstream from the DL in panel (b). The parallel temperature is now lower than the perpendicular one, which results in a pancake distribution. A cutoff at higher v_ex develops as high velocity electrons leave the regions upstream and downstream of the DL by the end of the simulation.

Zoom In Zoom Out Reset image size

Since a pancake distribution is preferentially trapped in a magnetic mirror and the geometry of a flare loop resembles a magnetic mirror, we next explore combining the DL with a mirror for the possibility of enhanced confinement.

4.1. Combining with a Magnetic Mirror

In the presence of a magnetic mirror, particles outside of a loss cone will be trapped. For the 2D velocity space sketched in Figure 9, the loss cone is bounded by the two dashed lines. When including a potential (provided by a DL here), the loss cone becomes a hyperboloid defined by the following equation:

$$\begin{equation} v_\parallel ^2 - v_\perp ^2(r-1) = \frac{2e\phi _{{\rm DL}}}{m_e} , \end{equation} \tag{ 5 }$$

where r is the mirror ratio and ϕ_DL is the DL potential. See Appendix A for a derivation. In the limit ϕ_DL → 0, this reduces to a cone. Also shown in Figure 9 are two contours of constant total energy. When going through a DL, an electron labeled as a blue blob will move to lower parallel velocity while its perpendicular velocity remains unchanged (following the blue arrow). This will create an anisotropy, helping to move electrons out of the loss cone.

Zoom In Zoom Out Reset image size

The hyperboloid intersects with the v_∥ axis at a value b = $\sqrt{2e\phi _{{\rm DL}}/m_e}$ (labeled in red). For a stronger DL, i.e., larger ϕ_DL, b has a greater value. This moves the loss hyperboloid to the right, resulting in a reduced overlap with the electron distribution. Therefore, fewer electrons will be lost. The same is true for a larger mirror ratio for which the hyperboloid closes closer to the v_∥ axis, making it thinner and hence the phase volume smaller.

4.2. Escaping Electron Density with the Addition of a Mirror

We now quantify the degree of electron confinement in the combined system. A similar calculation to that in Section 3.3 is carried out to determine the escaping electron density. We integrate the initial distribution over the loss hyperboloid. A derivation is given in Appendix B. The escaping density normalized to total density is given by

$$\begin{eqnarray} \tilde{n}_{{\rm esp}} &=&{\rm erfc}(\tilde{\phi }^{1/2}) - \exp [\tilde{\phi }/A(r-1)] \, {\rm erfc}([ \tilde{\phi }X ]^{1/2})/X^{1/2}, \nonumber \\ A & \equiv & \frac{T_\perp }{T_\parallel }, \, X\equiv 1+\frac{1}{A(r-1)}. \end{eqnarray} \tag{ 6 }$$

$\tilde{n}_{{\rm esp}}$ is plotted for several mirror ratios r in Figure 10, assuming isotropic initial electron distributions (A = 1). The r = 1 case, corresponding to the presence of a DL alone, is taken from Equation (4) and plotted for comparison with cases with a mirror. With neither a DL nor a mirror, i.e., for $\tilde{\phi }$ = 0 and r = 1, all electrons escape, so $\tilde{n}_{{\rm esp}}$ = 1. $\tilde{n}_{{\rm esp}}$ decreases with increasing $\tilde{\phi }$ and r. A typical mirror ratio of 2–2.5 from looptop to footpoint can be inferred from magnetohydrodynamic simulations of solar flares (Yokoyama & Shibata 1998; Birn et al. 2009). The mirror ratio for a region near the looptop is therefore expected to be moderate. Using $\tilde{\phi }\sim$ 0.7 as determined from simulations, $\tilde{n}_{{\rm esp}}$ goes down from ∼0.24 when no mirror is present to ∼0.14 for a moderate mirror ratio of 1.5. Therefore, combining a DL with a moderate mirror substantially enhances the confinement of hot electrons.

Zoom In Zoom Out Reset image size

We have demonstrated linearity in the scaling of the maximum DL strength with the hot electron temperature using PIC simulations and a simple theoretical model. In the latter, a single DL is assumed as the dynamics involving only one DL is considered. This is consistent with the simulations in which a single dominant DL is observed. However, in larger and more realistic systems, multiple DLs might possibly emerge. Formation of multiple DLs was reported in earlier 1D PIC simulations driven by a strong applied potential across longer systems that evolved for longer time than our simulations (Sato & Okuda 1981). The strengths of those DLs added up, constituting a greater total potential jump. If several DLs are generated at widely separated locations (as were observed in Sato & Okuda, 1981), it is unlikely that only one shock will form. Any shock that forms could affect the RC beam at its associated DL and also neighboring DLs, complicating the dynamics of the entire system. Whether the maximum strength for a single DL determined with this model can be straightforwardly applied to obtain the total maximum strength of multiple DLs remains to be studied.

Simões & Kontar (2013) presented evidence for electron accumulation at the looptop from a systematic study of solar flares with both looptop and footpoint hard X-ray sources. Magnetic mirroring was considered as a means to trap electrons along the flare loops. To explain the observation, mirror ratios of 2.5–5 are needed for isotropic distributions in a more realistic scenario that took into account Compton backscattering at the chromosphere and a neutral target at the footpoint (as opposed to zero backscattering and a fully ionized chromospheric target). Such mirror ratios are higher than estimates of less than 2.1 from investigations of large flare samples of 40–80 flares (Aschwanden et al. 1999; Tomczak & Ciborski 2007). This suggests that magnetic mirroring alone is not sufficient to trap the electrons at the looptop. As demonstrated in Section 4.2, DL formation, together with a magnetic mirror, can enhance electron trapping. For instance, from Figure 10, for a moderate mirror ratio of r = 1.5, having a DL with a strength of $\tilde{\phi }\sim$ 0.7 traps (1- $\tilde{n}_{{\rm esp}}$) = (1–0.14) = 0.86 of the total electron density while having the mirror alone traps a much lower density of (1–0.4) = 0.6. DL formation can, therefore, provide an explanation for electron accumulation at the looptop reported in Simões & Kontar (2013).

From recent PIC simulations of a pre-accelerated hot electron source, chosen as a basic model of flare-heated coronal sources, it is observed that the transport of the hot electrons is significantly suppressed by the formation of a DL (Li et al. 2012). The degree of suppression depends on the strength of the DL. In this work, a series of PIC simulations are performed to obtain a scaling of the DL strength with the hot electron temperature. A linear scaling relation is observed. The amplitude of the DL is limited by the formation of sound wave shocks produced by ions accelerated through the DL potential. An analytic calculation based on this model yields a linear scaling with the hot electron temperature, a result consistent with simulations. This study shows that a substantial fraction of electrons is confined by the DL. Thus, DLs can produce the electron confinement required by the prolonged decay time of X-ray looptop emissions in the solar corona.

The computations were performed at the National Energy Research Scientific Computing Center. This work has been supported by NSF grants AGS-1202330 and ATN-0903964.

The total energy of an electron in an electric potential ϕ(x) and subjected to a magnetic field B(x) is given by

$$\begin{equation} W = m_e v_\parallel ^2/2 - e\phi (x) + \mu B(x) , \end{equation} \tag{ A1 }$$

where μ = $m_ev_\perp ^2/2B(x)$ is the magnetic moment. Without other external forces, W is conserved. The potential jump of a DL can be approximated as a step function so that

$$\begin{equation} \phi (x)= \left\lbrace \begin{array}{@{} ll}\phi _{{\rm DL}} & x<x_{{\rm DL}} \\ 0 & x>x_{{\rm DL}} \end{array}\right., \end{equation} \tag{ A2 }$$

where x_DL is the DL position. Before the electron passes through the DL, i.e., x < x_DL, we have

$$\begin{equation} W= \frac{1}{2}m_e v_\parallel ^2 - e\phi _{{\rm DL}} + \mu B_1 , \end{equation} \tag{ A3 }$$

where B₁ is the value of the magnetic field at some location x < x_DL. After it passes through the DL, it continues to travel until its v_∥ reaches zero at some point with a stronger magnetic field B₂, so

$$\begin{equation} W= \mu B_2. \end{equation} \tag{ A4 }$$

Since W is a constant, we equate Equations (A3) and (A4), obtaining

$$\begin{eqnarray} \frac{1}{2}m_e v_\parallel ^2 - e\phi _{{\rm DL}} &=& \mu (B_2 - B_1) \nonumber \\ &=& \mu B_1(r-1)\,;\, r\equiv B_2/B_1 \nonumber \\ &=& \frac{1}{2}m_e v_\perp ^2 (r-1) \nonumber \\ \therefore \, v_\parallel ^2 -v_\perp ^2 (r-1) &=& \frac{2e\phi _{{\rm DL}}}{m_e}. \end{eqnarray} \tag{ A5 }$$

This is a hyperbolic equation in (v_∥, v_⊥) space.

In this case, the integral is not separable in v_∥ and v_⊥ as it is for n_esp (Equation (3)). In the following, we carry out the full 3D integral over (v_∥,v_⊥,ψ),where ψ is the azimuthal angle in the 2D plane spanned by the two perpendicular velocity components. Visualizing the loss hyperboloid in Figure 9 in 3D, the other perpendicular velocity axis is directed out of the page. v_⊥ of the loss hyperboloid is simply the cross-sectional radius of a cut of the hyperboloid taken at different values of v_∥ along the v_∥ axis. This radius is given by $v_{\perp l}^2$ = $({v_\parallel ^2-b^2})/({r-1})$ (from Equation (5)), where $b\equiv \sqrt{2\,e\phi _{{\rm DL}}/m_e}$. We then calculate the escaping electron density:

$$\begin{eqnarray} n_{{\rm esp}} &=& \int ^{2\pi }_0 \! d\psi \int ^{\infty }_{b} \! dv_\parallel \int ^{v_{\perp l}}_{0}\! v_\perp dv_\perp \, f_{e0}(v_\parallel ,v_\perp) \nonumber \\ &=& 2\pi \int ^{\infty }_{b} \! dv_\parallel \int ^{v_{\perp l}}_{0} \! v_\perp dv_\perp \frac{n_0}{\pi ^{3/2} v_{t\parallel } v_{t\perp }^2} \nonumber\\ &&\times \exp \left[-\left(\frac{v_\parallel ^2}{v_{t\parallel }^2} + \frac{v_\perp ^2}{v_{t\perp }^2}\right)\right] \nonumber \\ &=& \frac{n_0}{\sqrt{\pi } v_{t\parallel }}\int ^{\infty }_{b} \! dv_\parallel \, e^{-\frac{v_\parallel ^2}{v_{t\parallel }^2}} \, \left(1-\exp \left[- \frac{v_\parallel ^2 -b^2}{(r-1)v_{t\perp }^2} \right] \right). \nonumber\\ \end{eqnarray} \tag{ B1 }$$

Defining:

$$\begin{eqnarray} \,\,\,\, a&\equiv & \frac{1}{v_{t\parallel }^2} \,\,, \, \, c\equiv \frac{1}{(r-1) v_{t\perp }^2} \nonumber \\ n_{{\rm esp}} &=& n_0 \sqrt{\frac{a}{\pi }}\int ^{\infty }_{b} \! dx \left(e^{-ax^2} - e^{cb^2}\, e^{-(a+c)x^2} \right) ;\,x \equiv v_\parallel \nonumber \\ &=& n_0 \sqrt{\frac{a}{\pi }} \left(\int ^{\infty }_{b\sqrt{a}} \! \frac{dt}{\sqrt{a}} e^{-t^2} \, - \, e^{cb^2} \int ^{\infty }_{b\sqrt{a+c}} \!\frac{ds}{\sqrt{a+c}} e^{-s^2}\ \right); \nonumber\\ t &\equiv& x\sqrt{a}, \, s\equiv x\sqrt{a+c} \nonumber \\ &=&\frac{n_0}{2} \left[ {\rm erfc}(b\sqrt{a}) \, - \, e^{cb^2}\sqrt{\frac{a}{a+c}}\,{\rm erfc}(b\sqrt{a+c}) \right] , \end{eqnarray} \tag{ B2 }$$

where ${\rm erfc}(z) = ({2}/{\sqrt{\pi }}) \int ^{\infty }_{z}\,e^{-t^2}\,dt$ is the complementary error function. The total density of hot electrons moving to the right (i.e., the maximum n_esp if all hot electrons escape) is

$$\begin{equation} n_{{\rm tot}}=\int ^{\infty }_{0} \! d^3v \, f_{e0}(v_\parallel ,v_\perp)=\frac{n_0}{2} \end{equation} \tag{ B3 }$$

as expected. Therefore, the fraction of escaping to total density is

$$\begin{equation} \tilde{n}_{{\rm esp}}=\frac{n_{{\rm esp}}}{n_{{\rm tot}}} = {\rm erfc}(b\sqrt{a}) \, - \, e^{cb^2}\sqrt{\frac{a}{a+c}}\,{\rm erfc}(b\sqrt{a+c}). \end{equation} \tag{ B4 }$$

Expressing a, b, and c in terms of eϕ_DL, r and temperatures,

$$\begin{eqnarray} b\sqrt{a} &=& \frac{b}{v_{t\parallel }}=\sqrt{\frac{e\phi _{{\rm DL}}}{T_\parallel }}= \tilde{\phi }^{1/2} \nonumber \\ cb^2 &=& \frac{b^2}{(r-1) v_{t\perp }^2}=\frac{e\phi _{{\rm DL}}}{T_\perp (r-1)} \nonumber \\ &=&\frac{e\phi _{{\rm DL}}}{T_\parallel }\frac{T_\parallel }{T_\perp (r-1)}= \frac{\tilde{\phi }}{A(r-1)}; \quad A\equiv \frac{T_\perp }{T_\parallel } \nonumber \\ \frac{a}{a+c} &=& \frac{\frac{1}{v_{t\parallel }^2}}{\frac{1}{v_{t\parallel }^2}+\frac{1}{(r-1) v_{t\perp }^2}}=\frac{1}{1+\frac{v_{t\parallel }^2}{(r-1) v_{t\perp }^2} } \nonumber \\ &=&\frac{1}{1+\frac{T_\parallel }{T_\perp (r-1)}}=\frac{1}{1+\frac{1}{A(r-1)}}\nonumber \\ b\sqrt{a+c} &=& \sqrt{b^2\left(\frac{1}{v_{t\parallel }^2} + \frac{1}{(r-1) v_{t\perp }^2} \right) } =\sqrt{ \frac{e\phi _{{\rm DL}}}{T_\parallel }+ \frac{e\phi _{{\rm DL}}}{T_\perp (r-1)} } \nonumber \\ &=&\left[\tilde{\phi }\left(1+\frac{1}{A(r-1)} \right) \right]^{1/2}. \end{eqnarray} \tag{ B5 }$$

The final form of $\tilde{n}_{{\rm esp}}$ is therefore

$$\begin{eqnarray} \tilde{n}_{{\rm esp}} &=& {\rm erfc}(\tilde{\phi }^{1/2})-\exp \left[\frac{\tilde{\phi }}{A(r-1)} \right] \frac{1}{\sqrt{1+\frac{1}{A(r-1)}} }\nonumber\\ &&\times {\rm erfc}\left(\left[\tilde{\phi }\left(1+\frac{1}{A(r-1)} \right) \right]^{1/2} \right) , \end{eqnarray} \tag{ B6 }$$

where A ≡ (T_⊥/T_∥) is a measure of the electron anisotropy.