{'fraction_non_alphanumeric': 0.03145012998541627, 'fraction_numerical': 0.018197958277851752, 'mean_word_length': 4.399349537829511, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics of double layers, ion acceleration and heat flux suppression during solar flares

Observations of flare-heated electrons in the corona typically suggest confinement of electrons. The confinement mechanism, however, remains unclear. The transport of coronal hot electrons into ambient plasma was recently investigated by particle-in-cell (PIC) simulations. Electron transport was significantly suppressed by the formation of a highly localized, nonlinear electrostatic electric potential in the form of a double layer (DL). In this work large-scale PIC simulations are performed to explore the dynamics of DLs in larger systems where, instead of a single DL, multiple DLs are generated. The primary DL accelerates return current electrons, resulting in high velocity electron beams that interact with ambient ions. This forms a Buneman unstable system that spawns more DLs. Trapping of heated return current electrons between multiple DLs strongly suppresses electron transport. DLs also accelerate ambient ions and produce strong ion flows over an extended region. This clarifies the mechanism by which hot electrons in the corona couple to and accelerate ions to form the solar wind. These new dynamics in larger systems reveal a more likely picture of DL development and their impact on the ambient plasma in the solar corona. They are applicable to the preparation for in-situ coronal space missions like the Solar Probe Plus.


INTRODUCTION
Novel advances in imaging and spectroscopy have made possible observations of solar eruptions such as flares and coronal mass ejections (CMEs) with unprecedented resolution (Lin et al. 2002). This has led to better detection of radiation from energetic electrons in the tenuous corona, which, due to its lower plasma density, does not usually emit detectable radiation compared to the generally brighter emissions from the denser chromosphere. The corona is also proposed to be the electron acceleration site by flare models. Electrons can be accelerated to over 100 keV (Lin et al. 2003;Krucker et al. 2010Krucker et al. , 2007Tomczak 2009), three orders of magnitude higher than the ambient coronal temperature of ∼100 eV, during flares and CMEs. X-rays are emitted from the interaction between the accelerated electrons and ambient plasma. As energetic electrons propagate from the corona to the chromosphere, transport effects can modify the energy distribution of the propagating electrons and hence the observed Xray spectra, affecting the interpretation of acceleration models. The transport of electrons from a coronal acceleration site is therefore crucial to understanding energy release in flares.
Evidence for both confinement and free-propagation of energetic electrons is present in observations. A recent systematic study of solar flares with wellobserved looptop (coronal) and footpoint (chromospheric) emissions found that the number of electrons required to explain observations is 2-8 times higher at the looptop than that at the footpoint (Simões & Kontar 2013). This implies electron accumulation at the looptop. Another systematic study of the relation between coronal and footpoint X-ray sources indicated that the simple scenario of electron free-streaming from the corona towards the footpoints could not explain the observed difference in photon spectral indices between the coronal and footpoint emissions (Battaglia & Benz 2006). A filter effect in the propagation preferentially reducing the distribution at lower energies was required. Such a filter can be an electric field. The observations of abovethe-looptop sources reveal that the decay time of hard X-rays (HXR), which is a measurement of the lifetime of the HXR-producing electrons at the source, was over two orders of magnitude longer than the free-streaming transit time across the source (Masuda et al. 1994;Krucker et al. 2007Krucker et al. , 2010. This suggests that the energized electrons are trapped at the looptop. A statistical survey of 55 partially occulted flares in which the bright footpoints were blocked by the solar limb showed that many looptop X-ray observations share similar decay times and sizes (Krucker & Lin 2008), implying electron confinement in the corona is a common phenomenon in solar flares.
On the other hand, there is evidence for freepropagation of energetic electrons, i.e., no interaction with the ambient plasma, in time-offlight measurements of HXR emission (Aschwanden et al. , 1996. A systematic time delay between lower (25-50 keV) and higher (50-100 keV) energy HXR emission was reported from a statistical correlation study of over 600 flares . This demonstrates that electrons with lower energies arrive at the chromosphere after those with higher energies as they freely stream down the flare loop from the corona. The measurements supported the free-streaming scenario, apparently contradicting the confinement scenario. An important component of solar flare models is the transport of energetic electrons from the corona to the chromosphere and from the corona to the solar wind. This, however, remains poorly understood.
The transport of electron energy has previously been modeled as both classical conduction (Spitzer 1962) and anomalous conduction (Manheimer 1977;Tsytovich 1971;Smith & Lilliequist 1979;Levin & Melnikov 1993). Recent analytical and numerical studies of the transport of superhot electrons with energies greater than 10 keV showed that the classical conduction model produced a heat flux significantly higher than the real energy fluxes reported from multi-wavelength observations of solar flares (Oreshina & Somov 2011). This means that processes other than classical Coulomb collisions (between energetic electrons and ambient plasma) need to be considered as well.
In the anomalous conduction model, transport is limited by anomalous resistivity due to electron scattering in turbulent wave fields. Turbulence is produced by instabilities that result from the interaction between energetic electrons and ambient plasma. An anomalous conduction front at the head of an expanding hot electron source was considered as a way to confine hot electrons for the production of HXRs (Smith & Lilliequist 1979). Since it was based on a one-fluid description, the model did not resolve the ion inertial length, let alone capture processes developing at electron scales. Later, 1D electrostatic particle-incell (PIC) simulations 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 conduction front was studied in 1D electrostatic Vlasov simulations (Arber & Melnikov 2009). They identified a propagating temperature jump as a conduction front. It was, however, suggested by recent PIC simulations and analysis (Li et al. 2013) that the propagating front was an ion acoustic shock and that the temperature jump was a result of shock heating at its extremely sharp transition. Other observed nonpropagating temperature jumps associated with a potential jump were likely double layers (DLs), as discussed in Arber & Melnikov (2009).
More recently, the transport of coronal energetic electrons was investigated by electromagnetic PIC simulations (Li et al. 2013). It was shown that transport of flare-heated electrons from the source region was significantly suppressed. The suppression was due to the formation of a DL. Its associated electric potential reflected electrons back to the source region. A substantial fraction, about 50%, of the total hot electron density was confined by the DL. There was also a population of escaping electrons, which may represent the free-streaming population discussed above. This model is in qualitative agreement with the observational evidence of confinement and free-propagation of energetic electrons.
A DL is a localized region that sustains a large potential drop in a collisionless plasma (Block 1978;Raadu & Rasmussen 1988). The potential drop comes from a large-amplitude electrostatic electric field sandwiched between two adjacent layers of opposite charge. As a whole, the structure is neutral, but quasi-neutrality is locally violated within the layer, which occurs at scales of ∼10 Debye lengths λ De . An ideal DL is a monopolar electric field, corresponding to a monotonic drop in the potential φ. In general, a DL can be bipolar. There can be dips or bumps at the low or high potential sides, but an overall potential drop φ DL develops across the structure. The strength of a DL is measured by φ DL . Particles are reflected by a DL if φ DL is greater than their kinetic energy. For a thermal distribution, a DL with φ DL equal to the temperature of the distribution can reflect the entire thermal bulk of the distribution, which contains particles with velocities lower than the thermal speed. In Li et al. (2013), it is the reflection of energetic electrons by a DL that confines them in the source region.
Although DLs have not yet been observed in the solar corona, they are widely seen in various space plasmas, including the solar wind (Mangeney et al. 1999), the Earth's magnetosphere (Mozer et al. 1985;Ergun et al. 2001Ergun et al. , 2009) and Earth's radiation belt (Mozer et al. 2013). They are also inferred to exist in the magnetospheres of Jupiter (Hess et al. 2009) and Saturn (Gurnett & Pryor 2012). In the Earth's plasma sheet, a large number of DLs were detected by the THEMIS (Angelopoulos 2008) spacecraft during periods of high magnetic activity, implying that DLs may frequently occur in such situations (Ergun et al. 2009). In the solar atmosphere, many DLs are likely to develop during flares. Li et al. (2013) reported electron confinement by a single DL in PIC simulations. An important question is whether much larger, more realistic systems develop multiple DLs or a single dominant DL. If multiple DLs develop in larger scale systems, the transport properties of electrons are likely to be very different from that in a single DL system. Understanding this subject is important for future in-situ observations of the corona from missions like NASA's Solar Probe Plus (Guo 2010).
In this work, we investigate DL dynamics in simulation domains larger than those in earlier work. While only a single dominant DL is present in smaller domains, multiple DLs develop in larger domains. This gives rise to new dynamics not observable in a single-DL system. These dynamics are studied here for the first time. They include the following: the primary DL accelerates return current electrons to high velocities, producing electron beams that interact with ambient ions to drive the system Buneman unstable, resulting in the generation of many secondary DLs; trapping of electrons occurs between DLs, suppressing electron heat transport across these regions; DLs also accelerate ambient ions to produce strong ion flows across an extended region, which has important implications for understanding the formation of the solar wind.
In the following, we describe the setup of the simulations and the parameters suitable for solar flare settings (Section 2). The time evolution and key features of the simulations are summarized in Section 3.1. The generation of multiple DLs in larger domains are presented in Section 3.2. In Sections 3.3 and 3.4, new dynamics involved in the production of secondary DLs and the interaction of these DLs with the ambient plasma are delineated. In Section 4, we briefly discuss the implications of the present results for understanding the mechanism for driving the solar wind, compare them with previous relevant studies and finally, we summarize our results.

SIMULATION SETUP
We perform two-dimensional electromagnetic PIC simulations using the p3d code (Zeiler et al. 2002). The initial setup is the same as in Li et al. (2013). We use a rectangular domain (Figure 1) to represent a symmetric local segment of a flare loop centered at the looptop where a hot electron source is initialized. An initial uniform background magnetic field B 0 is applied in the long direction, x, which is along the loop axis. The initial density n 0 is uniform. Although the simulation domain is two-dimensional, no significant structure develops in the direction transverse to B 0 . Hence, we average over that coordinate to reduce noise. A third of the electrons, centered in the domain, have higher initial parallel temperature T h0, than the rest. It represents a preheated hot looptop source. Higher T 0, is used because parallel transport dominates over perpendicular transport in the presence of a strong magnetic field, which is typical for a coronal loop. At the initial state, the parallel temperature transitions from a higher to a lower value over a length scale of a few tens of λ De 's. The transition scale is not expected to affect the results because the temperature gradient broadens by more than an order of magnitude before DLs form. The broadening of the transition can be seen in Figure 2 of Li et al. (2012). In late time of the present simulations DLs form at various spatial locations far away from the the initial temperature gradient.
Electrons are initially modeled by bi-Maxwellian distributions and ions by a Maxwellian distribution. It is more natural that the hot population does not have a preferred direction of propagation. The simulated hot electrons are therefore not beamed in the initial state. We note that soft X-ray (SXR) counterparts of HXRs are observed in the high corona (Masuda et al. 1994;Krucker et al. 2010). SXR spectra can be well fit by a Maxwellian distribution (Masuda et al. 2000;Tsuneta et al. 1997) and HXR spectra are usually fit to a combination of a thermal core and a nonthermal tail distribution. Although we use Maxwellian distributions in our simulation setup, we do not expect our results to be sensitive to the specific form of the initial electron distribution as long as there is a sufficient difference in the energy between the hot and cold regions. We use unity plasma beta β (ratio of plasma pressure to magnetic pressure) in the parallel direction for hot electrons, corresponding to T h0, =1. Using unity β for hot electrons is consistent with recent coronal flare observations (Krucker et al. 2010).
Temperatures are normalized to m i c 2 A , where m i is ion mass and c A =B 0 /(4πm i n 0 ) 1/2 is the Alfvén speed. 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.
We perform four simulations with increasing domain sizes in the parallel direction, L x /L x0 =1, 2, 4 and 8. L x0 is the length of the smallest run, which has a size of L x0 × L y = 4634.1 × 18.1 λ 2 De . The grid size is 0.14×0.14 λ 2 De . There are 400 particles per cell. λ 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 0 e 2 /m e ) 1/2 is the electron plasma frequency. Using typical parameters for coronal thermal X-ray sources, T h =5 keV, n=10 9 cm −3 , we obtain λ De ∼ 2 cm. The largest simulation thus has L x ∼ 1 km, which is very small compared to the scale of a realistic flare loop. Space and time are normalized to λ De and ω −1 pe . A mass ratio m e /m i of 1/100 and speed of light c/v th0 of 7 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. With increasing domain size, a simulation can be evolved for longer times, which allows us to study the long-time behavior of the system.

The Impact of System Size on the Evolution of DLs
We present in figure 2 the time evolution of the total DL strength eφ DL in simulations with various system domain sizes. The largest domains contain more than one DL; for those the total DL strength is measured across all of them. eφ DL is measured from the electric field data, which is averaged † to eliminate initial fluctuations at the contact of the hot and cold electron regions. The fluctuations fade out after the first ∼1200 ω −1 pe . eφ DL at early times in the two larger runs are not shown † The electric field data is averaged over one electron plasma period in the two smaller runs and over five periods in the two larger runs. Data is not sampled as frequently in the latter as in the former, leading to a longer averaging time and hence less effective subtraction of the fluctuations in the larger runs. eφ DL of the smaller runs are therefore smoother while that of the largest run (red) appears to oscillate (about an overall growth trend) before ωpet=1200, but in general becomes smoother after that time. No averaging is needed after 1200ω −1 pe since the averaged and unaveraged data overlap. because of the presence of oscillatory fluctuations that are not sufficiently averaged out due to the less frequent sampling of data, used for these runs to improve computational efficiency in larger systems.
The two smaller runs with L x /L x0 =1 and 2 have similar evolutionary trends. The DL(s) in both cases grow, saturate and then grow again. The growth of DLs is driven by the Buneman instability due to the interaction of return current electrons and ambient ions (Li et al. 2012). The saturation mechanism is the formation of a sound wave shock that develops from ions being accelerated by the DL potential to supersonic speeds (Li et al. 2013). The shock stabilizes the Buneman instability and saturates further amplification of DLs.
The larger simulation domains in the two largest runs allow us to follow the dynamics well beyond the growth phase of the DLs. At late time, eφ DL /m i c 2 A in the two largest simulations settles to an average value of ∼ 0.6. The DL strength is sustained throughout the entire course of the simulations without showing any signs of decay. In earlier 1D particle simulations in which DLs were driven by a strong applied potential, DLs decayed as many solitons, in the form of spiky wave trains, propagated toward the high potential side of the DLs (Sato & Okuda 1981). They were observed to have a lifetime (of about 500 ω −1 pe ) comparable to the transit time of the solitons across the DL width. In our case, the DLs remain for a much longer time and exhibit no sign of decay.

Formation of Mulitple DLs in Larger Systems
In the smallest run, a single dominant DL is observed (not shown). In the second run (L x /L x0 =2), a weaker second DL emerges. We show in Figure 3 the time evolution of the electric field E x from that run. The bright white feature indicated by a cyan arrow (ω pe t∼1400) is the dominant DL. It sustains throughout the simulation. Another white feature indicated by a green arrow (ω pe t∼1100) is the weaker second DL. A DL is signified by an overall jump in the electric potential eφ, like those in Figure 5(d) or Figure 6. A shock wave (blue arrow) is generated at ω pe t∼300 (ion heating across the shock front in the form of ion trapping can be seen in Figure  6 of Li et al. (2013)). It leaves the DL at that time, causing saturation of the DL strength from ω pe =300 to 500. The plateauing of the orange curve in Figure 2 indicates saturation during that period. We note that there is no significant jump in the electric potential at a shock after it leaves a DL although a DL and a shock appear similar in Figure 3. This will be demonstrated in the next figure.
In the two larger runs (L x /L x0 =4, 8), many DLs are generated. In the following, we focus on data from the largest run to illustrate the dynamics in a multiple-DL system. To simplify notation, we drop the units of x and t from here on. The evolution of E x from the largest simulation is shown in Figure 4 . The first DL (cyan arrow) is the bright feature that emerges at (x, t) ∼ (0, 100), and is self-sustaining until t ∼ 3500 at x ∼ -140. It strengthens over time and propagates to the left. A shock associated with it forms at (x, t) ∼ (-14, 400) and separates from it thereafter. Note that at (x, t) ∼ (-15, 600), the shock, seen as a negative peak in E x in Figure 5 (d), does not cause a significant jump in eφ like a DL does. A second DL (green arrow) to the left of the first one emerges at (x, t) ∼ (-45, 400). A strong shock with an intense negative electric field is produced by this DL at t ∼ 2600 and separates from the DL. This will be further discussed later.

Generation Dynamics in Multiple-DL Systems
The fundamental generation mechanism for multiple DLs is the same as that for a single DL. As noted, it is due to the Buneman instability driven by the relative drift between a return cur-rent (RC) electron beam and the ambient ions (Li et al. 2012). The RC beam comes from the cold electron region and is a response to the streaming hot electrons coming from the hot electron region. However, in a larger system, there is more time for the system to evolve and significantly different dynamics and interactions are observed. The original RC beam (that drives the first DL) is accelerated by the first DL and the resulting high-velocity beam is responsible for the generation of the second DL. Figure 5 illustrates the early time evolution of the electron phase space near the interface between the hot and (ambient) cold electrons. The initial (t=0) electron temperature (in green) is overlaid in (a) to show the hot and cold regions. To the right of x ∼0 is the cold electron region. To the left is the hot electron region. A RC beam that is drawn from the cold electrons into the hot region is the beam that excites the Buneman instability. On top of initial fluctuations in E x at t=100 (black), the first DL emerges among the unstable waves near x ∼ -10.  The second DL is not yet present at this time. In (b), the first DL is developing at x ∼ -10 and the small second DL emerges at the position located by a vertical black line. The RC beam that is accelerated to v ex ≈ -0.5 by the first DL is the driver of the second DL. The beam density (cyan) is being dragged from negative to positive velocities near the vertical black line, indicating partial reflection of the beam by the second DL. The electrons reflected to the right of the second DL, at x ∼ -45 in (c), increase from t=250 (b) to 500 (c), transferring energy to the second DL. As a result, the second DL strengthens noticeably from t=500 to 600 (d). This can be seen from the increase in the potential (magenta) drop across the second DL. Therefore, the first DL accelerates the RC beam and drives the formation of the second DL. After a shock forms and departs at t ∼2600, the second DL weakens (Figure 4) because the shock reflects and hence holds back the RC electrons that can otherwise sustain the second DL (Li et al. 2013). As the shock diminishes after t ∼5000, the second DL revives. In Figure 4, one sees a gradual brightening since (x, t) ∼ (-270, 5000). There is also increasing wave activity to the right of the second DL. Around t ∼ 5300, a third DL, indicated by a blue arrow in Figure 4, develops at x ∼ -145 and remains until the end of the simulation. It appears to propagate slowly to the right, opposite to what is usually observed. This is because the ions are accelerated by the second DL, producing an ambient flow to the right in the region x > -270. As a result, the third DL is carried by the ambient plasma flow to the right.
In late time, the system reaches a state where many DLs arise. We show in Figure 6 the electron phase space at t=5600, some time after the emergence of the third DL. In addition to the second DL (at x ∼ -300) and the third DL (x ∼ -140), where significant potential jumps occur, a couple of weaker DLs are generated and overall there is substantial wave activity. Magenta arrows indicate the associated potential jumps. Note that the RC beam is reflected by the negative E x of the second DL and also by the weaker DLs (for instance, around x ∼ -235).
In larger domains that allow the system to evolve for a much longer time, the behavior at late time is much more dynamic than what was seen in smaller domains (Li et al. 2013, Sato & Okuda 1981. Multiple DLs are generated downstream of the RC beam from a primary DL (the first DL in this case). The growth of these downstream (to the left of the primary DL) DLs produces an ambient ion flow that pushes the upstream DL further to the right, opening up a space (in between) for more DLs to form. A chain of many DLs results.

Interaction with Ambient Ions
As discussed, DLs interact with ambient ions and create strong ion flows in a multi-DL system. We present in Figure 7 the space-time evolution of the ion flow velocity. When compared to the electric field evolution in Figure 4, one can see that the ions are strongly accelerated at the locations of DLs. The change in ion flow therefore resembles the pattern of the spatial and temporal development of the electric field. Following the path of the first DL, significant flow velocities of ∼0.4 c A ∼1 v ti (= 2T i /m i ) are produced. After being accelerated by the DLs to velocities as high as ∼2 v ti at the red patches in Figure 7, ions flow to the right (due to a drop in eφ) and the flow extends over a wider scale than the acceleration site, which is Fig. 6.-Electron phase space at t=5600, a time when many DLs form. Arrows indicate potential drops associated with a couple of small DLs between two larger DLs at x∼ -300 and -143. The same format is used as in Figure 5.
localized at the DLs. This is essentially the ambient plasma flow that carries the third DL to the right as mentioned in Section 3.3. Over time, as DLs spawn across the system, away from the proximity of the initial contact between the hot and cold electron regions, strong ion flows spread out in space as well. We discuss the implication of this process in Section 4.

Electron Trapping between DLs
Another new dynamics observed in multi-DL systems is the trapping of electrons between DLs. We show in Figure 8 the long time evolution of the electron phase space around two DLs. The electric field of both DLs grows over time. The second DL (on the left in (b) at x ∼ -75) is weaker than the first (on the right at x ∼ -45), i.e., it has a smaller drop in the potential than the first. Apparently, it is much more bipolar than an ideal monopolar DL as it develops a large negative E x (in (d)), which corresponds to a deep dip in the potential. This dip together with the potential drop at the first DL takes the form of a "cap" (an inverted potential well). Such a potential is capable of trapping electrons. The cap enlarges as the dip deepens from (a) to (d). With the following calculation, it can be demonstrated that the same (cyan) population of electrons stay within the cap during this entire period of ∼1000 ω −1 pe . These electrons are centered at approximately zero velocity with a velocity range of ∆v/v th0 ∼ ±0.3 (in (d)). That means that their typical velocity is of order ∼0.3. Given the separation between the two DLs being ∆x ∼ 20 λ De , the transit time of the electrons across the separation is τ transit = ∆x/∆v = 20/0.3 ∼ 70 ω −1 pe . This is much shorter than the 1000 ω −1 pe time scale. Therefore, they are trapped between the two DLs. In order to trap these electrons, the potential cap needs to have a depth of ∆φ/T h0 ∼ ∆v 2 /v 2 th0 ∼ 0.3 2 ∼ 0.1. In (d), taking the shallower left side, the cap has a depth of eφ cap /T h0 ∼ 0.2 (note T h0 =1), enough to trap these electrons. The trapped electrons are essentially a peak in phase space because of an enhancement in phase space density there. It has a spatial extent much wider than the typical DL width (of ∼ 10 λ De ) because the two DLs are widely separated. Given that electrons are trapped between DLs, there will not be any transport of energy across these regions. The more DLs develop, the

DISCUSSION AND SUMMARY
The impact of DLs on the ambient ions has important implications for understanding how hot electrons in the corona drive the solar wind. Hot electrons produced in the corona try to escape along the ambient magnetic field both toward the chromosphere and outward from the sun. In the process they drive vast numbers of DLs, which in turn accelerate the ambient ions, creating strong ion flows over large regions. Along open field lines, such ion flows move away from the the upper corona into the interplanetary space where they make up the solar wind. In the exospheric models of (fast) solar wind generation (Lamy et al. 2003), an assumed interplanetary potential is required to accelerate ions to the solar wind speed. Such a potential is likely to be manifest as a chain of DLs. Indeed, based on DL measurements by the WIND spacecraft near 1 AU, it was estimated that a succession of DLs from the corona to the Earth could produce the total electric potential required by the exospheric models (Lacombe et al. 2002).
We briefly compare some basic characteristics observed in the multi-DL systems in this work and those reported in previous 1D electrostatic PIC simulations driven by a strong applied potential (Sato & Okuda 1981). In this earlier work it was found that the number of DLs formed increased with the system length. We also observe a similar tendency here. It was also shown that the total strength of the DLs became stronger with the system length as more DLs were formed, but was far weaker than the externally applied potential. In our case, the total DL strength eφ DL increases moderately between the two smaller systems, indicated by the maxima of the black and orange curves in Figure 2. As the domain size further increases, eφ DL , however, appears to saturate (blue and red curves) and stays roughly at a constant value at late time. That value is comparable to the initial hot electron temperature. In Sato & Okuda (1981), since eφ DL was far less than the applied potential, which could be considered as the maximum available energy in the system, there was room for the DLs to grow stronger. In our case, however, eφ DL already comes close to the initial temperature, which can similarly be understood as an upper limit on the potential. Further strengthening beyond this value is unlikely. This may be why more DLs in larger systems do not result in a stronger total potential jump here as they did in Sato & Okuda (1981).
In conclusion, using large-scale PIC simulations, we reveal the self-consistent generation of multiple DLs during the transport of flare-heated electrons. Dynamics only observable in multi-DL systems are investigated for the first time. The primary DL enhances the return current, which is responsible for driving DLs, and results in subsequent breeding of secondary DLs. The chain of DLs occurs over an extended region. Electrons are trapped between DLs, suppressing the transport of electron heat in the trapped regions. The chain of DLs accelerates ambient ions, creating strong ion flows across the system. The simulations therefore reveal the mechanism by which hot electrons heated in the corona couple to ions to drive the strong outflows that make up the solar wind. This in turn separates adjacent DLs, expanding the DL breeding territory. These dynamics provide a more realistic picture of DL occurrence in the solar corona that is important for understanding energetic electron transport and in preparing for future space missions such as Solar Probe Plus.