Dynamical Modulation of Solar Flare Electron Acceleration due to Plasmoid-shock Interactions in the Looptop Region

In solar flares, magnetic reconnection is believed to play a crucial role in the explosive release of magnetic energy in the corona (Shibata & Magara 2011; Benz 2017). Observations have shown that an enormous number of particles are accelerated to high energies (up to tens of MeV for electrons) within tens of seconds to minutes and contain a significant fraction (as high as 10%–50%) of the released flare energy (Lin & Hudson 1976; Emslie et al. 2012; Aschwanden et al. 2017). Various mechanisms may contribute to particle energization in flares, including acceleration by parallel electric field, contracting/merging magnetic islands, large-scale compression in the reconnection layer (e.g., Drake et al. 2006; Oka et al. 2010; Zhou et al. 2015; Li et al. 2018a, 2018b), stochastic acceleration by magnetic turbulence or plasma waves (e.g., Miller et al. 1996; Petrosian & Liu 2004; Fu et al. 2020), acceleration by the shrinkage of reconnected magnetic field lines (e.g., Somov & Kosugi 1997; Karlický & Bárta 2006), and acceleration by a fast-mode shock, often referred to as the flare termination shock (TS), driven by the reconnection outflows (e.g., Tsuneta & Naito 1998; Mann et al. 2009; Warmuth et al. 2009; Guo & Giacalone 2012; Kong et al. 2013, 2019; Li et al. 2013; Nishizuka & Shibata 2013; Chen et al. 2015). However, it remains controversial which process plays a dominant role and can explain various observational signatures of particle energization in solar flares (Miller et al. 1997; Aschwanden 2002; Zharkova et al. 2011).

Hard X-ray (HXR) and radio observations provide primary diagnostics of the acceleration and transport of energetic electrons in solar flares. Nonthermal looptop sources suggest that particle acceleration takes place above the top of flare loops and the TS is one of the promising candidates as the acceleration mechanism (e.g., Masuda et al. 1994; Melnikov et al. 2002; Liu et al. 2008, 2013; Krucker et al. 2010; Krucker & Battaglia 2014; Oka et al. 2015; Gary et al. 2018). Although the TS has long been predicted in magnetohydrodynamic (MHD) simulations when fast reconnection outflows impinge upon the top of newly reconnected magnetic loops (e.g., Forbes 1986; Magara et al. 1996; Yokoyama & Shibata 1998), there is rarely solid observational evidence. One of the reasons is that the spatial size of the TS is expected to be very limited (Chen et al. 2019). In several studies, slow-drift radio emissions similar to Type II radio bursts (associated with shocks driven by coronal mass ejections) were interpreted as the radio signature of the TSs (Aurass et al. 2002; Mann et al. 2009; Warmuth et al. 2009). Recently, using the Karl G. Jansky Very Large Array, Chen et al. (2015) presented a high-cadence radio spectroscopic imaging observation of stochastic spike bursts at decimetric wavelengths during a solar flare, which had a morphology and dynamic evolution most likely representative of a flare TS as suggested by numerical simulations. The observation suggested that the TS is located at the ending points of plasma downflows and slightly above a coronal HXR source. Chen et al. (2019) further performed detailed analysis of the split-band feature in the same radio spike burst event, and found that the high-frequency band is located slightly below the low-frequency band which, in turn, supported the interpretation that the two bands were emitted in the downstream and upstream side of the TS, respectively.

In addition to efficient particle acceleration, in some flare events confinement of electrons in the looptop region is also required to account for the nonthermal looptop emissions (e.g., Simões & Kontar 2013). Several mechanisms have been suggested, including magnetic mirroring and turbulent pitch-angle scattering (e.g., Simões & Kontar 2013; Kontar et al. 2014; Musset et al. 2018; Ruan et al. 2020). In theoretical models of particle acceleration in solar flares, the TS is often considered as a planar standing shock (e.g., Tsuneta & Naito 1998; Mann et al. 2009; Nishizuka & Shibata 2013). However, recent MHD simulations have shown that the TS evolves dynamically and exhibits complex structures due to the impacts of jets/plasmoids (e.g., Takasao et al. 2015; Takahashi et al. 2017; Takasao & Shibata 2016; Shen et al. 2018; Cai et al. 2019; Zhao & Keppens 2020). Particularly, a concave-downward magnetic structure has been shown to be present below the TS in the looptop region (Takasao et al. 2015; Shen et al. 2018; Kong et al. 2019). Such a magnetic configuration is favorable for trapping electrons because it is more difficult for particles to travel transverse to the magnetic field than along it (Guo et al. 2010; Kong et al. 2015, 2016).

There has been a lack of studies that couple a kinetic energetic-particle model with a realistic MHD simulation of solar flare region. Recently, Kong et al. (2019) presented a model by numerically solving the Parker transport equation (Parker 1965) with the plasma velocity and magnetic field obtained from MHD simulations. They found that the low-energy electron spectrum resembles a power law, which agreed well with the prediction of diffusive shock acceleration (DSA) theory. The accelerated electrons are concentrated in the looptop region due to the acceleration at the TS and confinement by the magnetic trap structure, in agreement with HXR and microwave observations. Therefore, the model in Kong et al. (2019) has strong implications regarding the acceleration and confinement of electrons in the nonthermal looptop sources.

During solar flares, fast temporal variations of flare emissions, including the so-called quasi-periodic pulsations (QPPs), are commonly observed in multiple wavelengths including radio and HXRs (Nakariakov & Melnikov 2009; McLaughlin et al. 2018). In addition, radio and HXR imaging occasionally shows multiple sources at or above the flare looptops (e.g., Tomczak 2001; Petrosian et al. 2002; Sui & Holman 2003; Liu et al. 2008, 2013; Gary et al. 2018; Chen et al. 2020; Yu et al. 2020), suggesting more complicated dynamics during electron acceleration and transport. Kong et al. (2019) explored electron acceleration and distribution in a quasi-steady and nearly symmetric phase of the TS. Recent theoretical work and numerical simulations have shown that a large-scale current sheet can break into numerous plasmoids (e.g., Shibata & Tanuma 2001; Loureiro et al. 2007; Bhattacharjee et al. 2009). As shown in Shen et al. (2018), plasmoids are produced intermittently in the reconnection current sheet and interact dynamically with the TS. As a result, key properties of the TS, such as the compression ratio, Mach number, and shock oblique angle, can significantly vary with time.

In this study, we investigate the dynamical modulation of electron acceleration and transport in the looptop region due to plasmoid-shock interactions. This Letter is structured as follows. We briefly introduce our numerical methods in Section 2 and present detailed analysis of the simulation results in Section 3. We summarize and discuss the implications of this work in Section 4.

The numerical methods used in this study are nearly identical to those in Kong et al. (2019). Here we provide salient details for completeness. We first simulate the magnetic reconnection-driven TS in a classic two-ribbon solar flare geometry with a 2.5D resistive MHD method (Chen et al. 2015; Shen et al. 2018), then we model the acceleration and transport of electrons in the looptop region using the MHD fields in a post-processing manner.

We perform a MHD simulation of a solar flare in the x–y plane. The resistive MHD equations are solved with the Athena code (Stone et al. 2008). The initial configuration is a vertical Harris-type current sheet. The magnetic field lines are fixed to the bottom boundary by using a line-tied condition to model the two-ribbon flare. We use a uniform resistivity that gives the Lundquist number S = 10⁵. The normalization units are given as: the length L₀ = 75 Mm, the magnetic field B₀ = 40 G, the velocity V₀ = 810 km s⁻¹, and the time t₀ = 92 s. Interested readers are referred to Shen et al. (2018) for details of the MHD setup and simulation results.

While Kong et al. (2019) has investigated electron acceleration and transport during a relatively steady period between 96.5–97.5 t₀, here we will pay special attention to a more dynamical evolution period during 97.4–98.5 t₀. To see the influence of plasmoid-shock interactions more clearly, we first model the acceleration of electrons for 3 t₀ (between 94.4 and 97.4 t₀ in Figure 2(a)) by using the fixed MHD frame at 97.4 t₀, which provides a quasi-steady distribution of energetic electrons before the plasmoids arrive. After 97.4 t₀, the particle simulation is performed with regard to the time-dependent MHD frames when the plasmoids interact with the flare TS. We select the flare looptop region given by x = [−0.15, 0.15] and y = [0.4, 0.7] for particle simulation. The temporal cadence between adjacent MHD frames is 0.002 t₀, compared to 0.01 t₀ in Kong et al. (2019), which better captures the dynamical evolution of the TS at shorter timescales.

Based on the plasma velocity and magnetic field from the MHD simulation, we model the acceleration and transport of electrons by numerically solving the Parker transport equation (Parker 1965). It is achieved by integrating stochastic differential equations corresponding to the Fokker–Planck form of the transport equation using a large number of pseudo-particles (e.g., Zhang 1999; Giacalone & Neugebauer 2008; Guo et al. 2010; Kong et al. 2017; Li et al. 2018b). The transport of charged particles in the magnetic field is described by the spatial diffusion coefficient. The diffusion coefficient parallel to the magnetic field is calculated from the quasi-linear theory (Jokipii 1971; Giacalone & Jokipii 1999), with an energy dependence ${\kappa }_{\parallel }={\kappa }_{\parallel 0}{(E/{E}_{0})}^{2/3}$. Following Kong et al. (2019), we take κ_∥0 = 0.005 κ₀, where κ₀ = L₀V₀. We also consider perpendicular diffusion similar to results of test-particle simulations in synthetic turbulence (Giacalone & Jokipii 1999) and set κ_⊥/κ_∥ = 0.1. In MHD simulations the TS can be resolved in several cells, which means that the shock width is on the order of one grid cell, Δx ∼ 0.001 L₀. Therefore, in nearly all regions the characteristic diffusion length at the lowest energies is larger than the grid cell, κ_nn/V₀ > Δx, where κ_nn is the diffusion coefficient in the shock normal direction. We set the time step Δt = 10⁻⁵ t₀ at the injection energy to ensure that pseudo-particles can "see" the shock transition. Kong et al. (2019) has shown that with these parameters the electron spectrum in the low-energy range resembles a power law, close to that predicted by the DSA theory. Because the shock shape is very complex and evolves dynamically, here we inject particles uniformly in the particle simulation domain. The initial energy of electrons is fixed to be E₀ = 0.5 keV (corresponding to the typical energy of ∼6 MK plasma in the flaring region) and a total of 10⁸ pseudo-particles are injected at a constant rate.

In the MHD simulation, plasmoids form intermittently at the primary X-point at y ∼ 1.5 L₀ and move upward/downward surfing in the reconnection outflows. We examine three downward-moving plasmoids as they crash into the looptop region during 97.4–98.5 t₀ and interact with the TS. Figures 1(a)–(b) show the maps of magnetic field strength B and plasma velocity divergence ∇ · V at 97.4 t₀. The magnetic field is stronger in the plasmoids than that of the ambient, likely due to compression. The TS locates at y ∼ 0.6 L₀, as well manifested by negative ∇ · V regions. Shen et al. (2018) showed that the TS is a sharp transition layer, across which the flow speed decreases and the Mach number quickly drops from ∼2 to less than 1. Figure 1(c) shows the time–distance plot of ∇ · V across the TS. We find that the height of the TS varies dynamically. The average shock compression ratio X ranges from 1.7 to 2.5, and shock angle θ_Bn ranges from 30° to 75° (see Figures 2(d)–(e)). As shown below, the TS can be deformed and restored, and multiple shocks/compression regions can be generated during the impacts of plasmoids. Kong et al. (2019) has shown that electrons can be accelerated efficiently by the TS, leading to a concentrated electron population in the looptop region. We expect that the dynamical plasmoid-shock interactions will affect both the shock properties and the acceleration of electrons.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 2(a) shows temporal variations of the numbers of electrons at different energies, 5–10 keV, 15–25 keV, 25–50 keV, and >50 keV, respectively. The electron numbers are integrated over the looptop region, x = [–0.1, 0.1] and y = [0.5, 0.7], and normalized to their respective values at 97.4 t₀. As noted in Section 2, the MHD background before 97.4 t₀ is fixed. The number of energetic electrons at different energies rise gradually. Collectively, they develop a power-law energy spectrum, and reach a quasi-steady state prior to 97.4 t₀. After 97.4 t₀, when the energetic-particle simulations are performed with regard to the time-dependent MHD evolution, the three plasmoids collide with the TS successively with similar time intervals. The electron number profiles differ strongly at different energies. Overall, after the series of plasmoid collisions, compared to the pre-collision steady-state values, the total number of low-energy (5–10 keV and 15–25 keV) electrons increases, but that of high-energy (>25 keV) electrons decreases.

The production of energetic electrons due to DSA depends on a number of factors such as the compression ratio, the shock angle, and magnetic field configuration in the looptop region. Here we attempt to explore how the variation of number of energetic electrons depends on the properties of the TS. From the Parker equation we can find that the rate of particle acceleration is related to the magnitude of plasma compression, 〈dp/dt 〉 ∼〈−p∇ · V〉, where p is the particle momentum and V is the flow velocity (Jokipii 2012). Figures 2(b)–(c) plot the temporal variation of the total number of cells with compression in the looptop region. The compression area is illustrated by the number of cells where the velocity divergence ∇ · V < −50 $[{t}_{0}^{-1}$], and we use a criterion of ∇ · V < −150 $[{t}_{0}^{-1}]$ for strong compression (mainly located at the TS). The number of low-energy (5–10 keV) electrons correlates well with the total number of compression cells, with a correlation coefficient cc = 0.88, and high-energy (>50 keV) electrons show good correlation with the total number of strong compression cells, with cc = 0.70. Figures 2(d)–(e) display the temporal variations of the shock compression ratio and shock angle averaged over the looptop region. There may be weak correlations with the >50 keV electrons (cc = 0.38 and 0.49, respectively). Note that here we only consider the primary shock at the front of downward reconnection outflows (i.e., the TS). As shown in Figures 4 and 5, multiple shocks/compression regions can form in the looptop region and contribute to particle energization. We can only find moderately good correlations for high-energy electrons because they demand cruder conditions on various factors, such as the shock compression ratio and shock angle, and the background magnetic field configuration, and the acceleration timescale is longer than the variation timescales of these factors.

We find that the number of energetic electrons exhibit rapid variations at short (<0.1 t₀ or <9 s) timescales. However, the variations during each plasmoid-shock interaction are not exactly the same. This is because various factors, such as the size of plasmoids, the TS properties, and the magnetic field configuration, can affect the acceleration and transport of electrons. In particular, the first plasmoid, denoted "P1," is relatively small in size, which merges into the looptop region immediately after crossing the TS front. The size of the plasmoid is about 0.02 L₀ = 1.5 Mm, or about 2'' on the Sun. Chen et al. (2015) showed that a temporary disruption of the TS by a fast plasma downflow with about the same size coincides with the reduction of the HXR and radio flux. In this study, we mainly focus on examining the first plasmoid-shock interaction, which is comparable to the situation in Chen et al. (2015).

The first plasmoid-shock interaction occurs between 97.5–97.85 t₀. As shown in Figures 2(b)–(e), at low energy of 5–10 keV, the electron number first decreases slightly and then increases gradually. At higher energies, the numbers decrease continually. The amount of decrease is more significant at higher energies. Eventually, the number of 15–25 keV electrons is reduced by ∼20% and nearly 40% for 25–50 keV. Both the trends and magnitudes are generally consistent with the evolution of X-ray flux observed by the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI) and Fermi/Gamma-ray Burst Monitor (GBM) in Chen et al. (2015). Figure 3(a) shows the energy spectra of accelerated electrons at three times as marked by black arrows in Figure 2(a). The energy spectra below 10 keV are close to a power law with a spectral index of δ = 2.5. The spectral index is consistent with the DSA prediction: if we take X ≈ 2, a typical value for the flare TS as predicted in the MHD simulation (Shen et al. 2018) and inferred from the observed split-band feature (Chen et al. 2019), the DSA predicts a power-law distribution with the spectral index δ = (X + 2)/[2(X − 1)] ≈ 2 in the non-relativistic limit. We note, however, in RHESSI observations, the X-ray emission from the low-energy part of the nonthermal electron spectrum (at  ≲10–20 keV) is usually "buried" under that from flare-heated thermal plasma (Holman et al. 2011). The energy spectra above 20 keV deviate from a simple power law and get softer with time. Figure 3(b) shows a more detailed view of the spectra at >10 keV. At the energy range of several tens of keV, the degree of softening is ∼1. It is similar to that reported in Chen et al. (2015), but we note that in their event the nonthermal "tail" is only observed (i.e., above the background) up to ∼25 keV, and they fitted the RHESSI X-ray spectra in the 10–25 keV range with a single power law.

Zoom In Zoom Out Reset image size

In Figure 4, the first column plots the maps of ∇ · V, and the other two columns plot the spatial distributions of accelerated electrons at 5–10 keV and 25–50 keV, respectively. Right before the collision, in panels (a)–(c), both low- and high-energy electrons are concentrated downstream of the TS, and the size at low energy is relatively larger, which is consistent with the results in our previous study (Kong et al. 2019). Around 97.64 t₀, in panels (d)–(f), the centroid of plasmoid P1 is moving across the TS. The low-energy electrons are more uniformly distributed, and the number of 25–50 keV electrons has been largely reduced. After the collision, in panels (g)–(i), the TS is distorted and separated into three parts, including one horizontal and two oblique shocks. Meanwhile, multiple shocks/compression regions form in the looptop region, possibly due to plasma backflow. Similar shocks have been discussed in previous MHD simulations (e.g., Takasao & Shibata 2016; Takahashi et al. 2017; Zhao & Keppens 2020) and detailed analysis is beyond the scope of this work. As discussed above, the number of 5–10 keV electrons increases, in accordance with the increase of compression cells. But the number of 25–50 keV electrons continually decreases. Chen et al. (2015) showed that nonthermal electron population is reduced but not eliminated during the disruption of the flare TS that coincided with the arrival of a plasma downflow. Our simulation results are consistent with their observations.

Zoom In Zoom Out Reset image size

We now analyze the effects of the other two plasmoids, P2 and P3. In Figure 5, the first two rows show the distributions of accelerated electrons during the impact of plasmoid P2. In the course of P2 moving across the TS, a very compact source appears both at low and high energies. When P2 reaches the downstream region, an intense source is distributed around the TS, similar to the case at 97.56 t₀ in Figure 4. As shown in Figures 2(d) and (e), during 97.9–98.0 t₀, the TS has a nearly perpendicular shock geometry and a high compression ratio, both of which are favorable for efficient particle acceleration. This explains why the electron numbers increase when P2 moves across the TS. We can also find another smaller source at y ∼ 0.55 L₀. The distance between the two sources is ∼0.06 L₀ ∼5 Mm. Later, during 98.0–98.1 t₀, the horizontal shock is deformed. As shown in Figure 2, both the shock compression ratio and shock angle decrease. As a result, the growth rate of electron number becomes slower, and the number of >50 keV electrons even is reduced. In the bottom row, multiple sources appear after the collision of plasmoid P3. The upper source is located near the TS. The lower source is located at ∼0.6 L₀, which coincides with a compression region existing for nearly 0.1 t₀, as shown in Figure 1(c).

Zoom In Zoom Out Reset image size

In this study, we perform a numerical simulation of electron acceleration and transport in the flare looptop region. We explore the impacts of plasmoids that intermittently crash into the TS. As a result, the TS evolves dynamically and can be distorted and restored. We find that the energetic electron population shows a rapid variation in both time and space, with distinct differences at low and high energies. The number of low-energy (5–10 keV) electrons correlates strongly with the total area of compression regions, while the high-energy (>50 keV) electrons show good correlation with strong compression cells but moderate correlations with the shock compression ratio and shock angle.

We have focused on studying the impacts of the first plasmoid-shock interaction, when the shock front first experiences a transition from a quasi-steady state to a highly dynamic period. The simulation results during this interaction compare favorably with the observations reported in Chen et al. (2015) regarding the evolution of the TS front (outlined by the radio spike centroids) as well as the associated X-ray and radio emissions. The plasmoid merges into the looptop region immediately after crossing the TS. The number of 15–25 keV electrons is reduced by ∼20% and nearly 40% for 25–50 keV, while the number of 5–10 keV electrons increases. The spatial distributions of accelerated electrons also show that the number of 25–50 keV electrons has been greatly reduced throughout the looptop region. In addition, in accordance with the observations of the X-ray spectra reported in Chen et al. (2015), the energy spectra above 10 keV evolves softer in time immediately following the arrival of plasmoid P1. During the subsequent two plasmoid-shock interactions, both the electron number profiles and spatial distributions display a more complicated temporal and spatial evolution, with multiple sources appear when the plasmoids interact and move across the TS front.

In solar flares, QPPs are ubiquitously observed in almost all wavelengths including HXR and radio emissions (Nakariakov & Melnikov 2009; McLaughlin et al. 2018). The typical periods range from a fraction of a second to several minutes (e.g., Tan et al. 2010; Yuan et al. 2019; Li et al. 2020), therefore possibly related to processes in the MHD regime. In our simulations, in response to the arrival of the three plasmoids, the total electron number at different energies shows a quasi-periodic temporal variation with a time interval of tens of seconds. Although the limited number of plasmoid interactions in our simulation renders it difficult to directly associate the temporal evolution of the energetic electrons (and the nonthermal emissions) to the so-called QPPs, our results suggest that the spatial and energy distribution of the electrons in the looptop region can indeed exhibit rapid temporal evolution in response to the arrival of the plasmoids at the TS front. If the formation of the plasmoids exhibits a quasi-periodic behavior, as was demonstrated in a number of numerical studies (e.g., Shen et al. 2011), our results may provide a viable explanation for the formation of the QPPs. We note that this scenario was previously suggested by Takasao & Shibata (2016), although they did not have particle simulations to explicitly show the temporal and spatial evolution of energetic electrons in the vicinity of the flare TS. The periodicity in our simulation is highly correlated with the production rate of plasmoids in the flare current sheet. Both in linear theory (Loureiro et al. 2007) and 2D MHD simulations (e.g., Cassak et al. 2009; Samtaney et al. 2009; Huang & Bhattacharjee 2010), it has been shown that the number of plasmoids formed in the reconnection current sheet depends on the value of Lundquist number. If the Lundquist number is much larger, as in the realistic solar corona, many more but smaller plasmoids can be generated. The size (magnetic flux) distribution of plasmoids has also been investigated (e.g., Fermo et al. 2010; Uzdensky et al. 2010; Huang & Bhattacharjee 2012; Shen et al. 2013; Ni et al. 2015). Therefore, the value of the Lundquist number could have an important effect on the properties of the TS, the acceleration of energetic electrons, and the periodicity in flare emissions. This will be explored in detail in our future work. In addition, the structure in the current layer has been shown to be even more complex in 3D numerical simulations (e.g., Guo et al. 2015; Yang et al. 2020). However, it remains unknown how the 3D turbulence and physics influence the physical properties and dynamical evolution of the flare TS.

X-ray and microwave emissions in the corona occasionally show multiple sources (e.g., Tomczak 2001; Petrosian et al. 2002; Sui & Holman 2003; Liu et al. 2008, 2013; Yu et al. 2020). Particularly, the double coronal X-ray sources often display energy-dependent feature, i.e., the centroid locations becoming closer at higher energies. The distance between the upper and lower sources range from a few arcseconds to tens of arcseconds and the sources can be visible for several minutes and even longer. The double X-ray sources are attributed to particle acceleration taking place both in the upward and downward reconnection outflow regions (Liu et al. 2013). MHD simulations have shown that TSs can be generated both at the flare looptop and the bottom of the flux rope (e.g., Takahashi et al. 2017; Zhao & Keppens 2020). In our work, we have only focused on the lower TS in the looptop region. The distance between the upper and lower sources is a few Mm and only lasts for a few seconds. Therefore, it cannot explain the double coronal X-ray sources as shown in Liu et al. (2013), but predicts that multiple sources can occur sporadically in the looptop region due to the modulations in electron acceleration and transport when jets/plasmoids interact with the TS. We note that a recent study by Yu et al. (2020) does show such a double X-ray (and microwave) source in the looptop region during the post-impulsive phase of the 2017 September 10 X8.2 flare.

There has been a lack of studies that couple a kinetic energetic-particle model with a realistic MHD simulation of solar flare region. We present a model by numerically solving the Parker transport equation with the plasma velocity and magnetic field from MHD simulations. Our simulations can reproduce the energy spectrum and spatial distributions of energetic electrons necessary for explaining the nonthermal looptop sources. In our future work, we will combine our MHD-particle model with radiation models to investigate the dynamical evolution in HXR and microwave emissions during solar flares.

The authors thank Xiaocan Li and Sijie Yu for helpful discussions. This work was supported by the National Natural Science Foundation of China under grants 11873036, 42074203 and 11790303 (11790300), the Young Elite Scientists Sponsorship Program by China Association for Science and Technology, and the Young Scholars Program of Shandong University, Weihai. F.G. is supported by NSF grant AST-1735414 and DOE grant DE-SC0018240. B.C. acknowledges support by NSF grant AST-1735405 to NJIT. C.S. acknowledges support by NSF grants AST-1735525, AGS-1723313, AGS-1723425 to SAO. J.G. acknowledges support from NSF grant 1735422. The work was carried out at National Supercomputer Center in Guangzhou (TianHe-2).