Simulation of a Solar Jet Formed from an Untwisting Flux Rope Interacting with a Null Point

Coronal jets are one of the most common phenomena in the solar atmosphere. They are typically collimated beams with helical motions in a transient process. Modern observations of jets have already been taken through a wide range of wavelengths from X-ray to extreme-ultraviolet (EUV; Shibata et al. 1992; Cirtain et al. 2007; Patsourakos et al. 2008). Observations with high spatiotemporal resolutions from the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) on board the Solar Dynamics Observatory (SDO; Pesnell et al. 2012) and simultaneous spectral observations from the Interface Region Imaging Spectrograph (IRIS; De Pontieu et al. 2014) provide an opportunity to study the origin and development of jets (Joshi et al. 2020a; Yang et al. 2020). Although the origins and morphology of jets have been investigated for decades, the results are still inconclusive (Liu et al. 2011a; Shibata et al. 1994; Raouafi et al. 2016; Schmieder 2022). Usually jets can be classified into standard jets, which fit the standard reconnection picture, and blowout jets, analogous to a miniature version of blowout explosion at their bases (Moore et al. 2010). Coronal jets extend in length from tens of megameters to about 400 Mm and possess a speed from about 100 to 600 km s⁻¹ (Nisticò et al. 2009; Paraschiv et al. 2010). They are not isolated phenomena and are believed to be closely related to other activities like coronal mass ejections, particle acceleration in the solar wind, and even coronal heating (Nitta et al. 2006; Yu et al. 2014; Joshi et al. 2020b). There are still many questions to be solved for solar coronal jets.

Magnetic reconnection is assumed to be responsible for explosive jet events in the corona. Fan-spine structure is composed of a fan-like dome with a null point connecting an outer spine to a remote region (Parnell et al. 1996; Priest & Titov 1996). The reconnection can easily occur at a null point and transfer the plasma to the far end of the spine. When reconnection takes place at a null point, energy stored in the magnetic field is continuously converted into the kinetic energy and thermal energy of plasma. And then triggered by a breakout structure, an intensive explosion or rapid motion of plasma can be observed in the corona. In terms of the magnetic topology, fan-spine structures (Liu et al. 2011b), eruption of magnetic flux ropes (Adams et al. 2014; Zhu et al. 2017), and bald patches with a continuous buildup of electric current (Guo et al. 2013; Schmieder et al. 2013) may easily trigger jet events. A bald patch is a region where the magnetic field lines are tangential to the photosphere and the direction of field lines points from negative polarity to positive polarity (Titov et al. 1993). The magnetic configuration of parasitic polarities situated in a large opposite polarity tends to generate such topologies. In some studies, successive jets in coronal holes exhibit a homologous behavior of recurrence, resulting from impulsive reconnection between closed and open field lines (Pariat et al. 2010).

Although spectral and EUV images have been widely used in the analysis of the jet formation and evolution, the lack of simulation is still a restriction on exploring the nature of jets (Schmieder 2022). Recently, state-of-the-art magnetohydrodynamic (MHD) models with high-performance computing techniques help us to understand the origin of jets. Fang et al. (2014) simulated the process of a twisted flux rope emerging in an open magnetic field and found that magnetic reconnection is responsible for the explosive events. González-Avilés et al. (2020) performed simulations of different parameters to probe the impacts on the origin and morphology of the jets. They found the significance of the resistivity and thermal conductivity in constructing the jets. Wyper et al. (2019) simulated a helical jet caused by a sigmoidal flux rope eruption, with artificial surface motions injecting free energy into the active region (AR). This model exhibits a coupled mechanism of the breakout reconnection and ideal MHD instability in the process of jet eruption. However, most of the simulations about jets mentioned above are ideal MHD models controlled by artificial settings. Simulations utilizing magnetograms from observations have not yet been widely studied. A force-free field is widely used to study the stable magnetic configuration before eruptions (Wiegelmann et al. 2006; Wiegelmann & Sakurai 2021). Such fields can be described by $\left({\rm{\nabla }}\times B\right)\times B\,=\,0$, and it becomes ${\rm{\nabla }}\times B=\alpha \left(r\right)B$ if considering the divergence-free condition. When the torsional parameter $\alpha \left(r\right)$ equals zero or is a constant, the field is a potential field or a linear force-free field, which can be solved in such particular situations. To further explore the actual physical conditions, a nonlinear force-free field (NLFFF) is solved where $\alpha \left(r\right)$ varies with position. Although no general analytic solutions exist because of the nonlinearity of the equation, various algorithms have been proposed and implemented to solve this issue (Schrijver et al. 2006; Metcalf et al. 2008). In this aspect, an NLFFF extrapolated from the observed magnetogram and the data-driven simulation based on such a field further help us to understand the magnetic structures of the AR and nature of flux rope eruptions (Guo et al. 2019), including some complex solar eruptions (Jiang et al. 2022). Jiang et al. (2016) simulated a coronal magnetic field that transformed from the pre-eruptive to eruptive state following a long-duration quasi-static evolution, and they found the establishment of a positive feedback between emerging flux and external magnetic reconnection. Therefore, a promising method is to perform the data-constrained or data-driven simulation that can reproduce the jet eruption with helical motions as observed.

Recently, statistical studies on the origin of coronal jets draw a conclusion that mini-filament eruptions play a critical role in triggering both standard and blowout jets (Baikie et al. 2022; Sterling et al. 2022). The rotational motion of a mini-filament along with its upward motion to the top of an arch system can result in magnetic reconnection. The reconnection between the inner and ambient field lines decreases the bound force imposed on the mini-filament and thus leads to its eruption. The jet followed by the mini-filament eruption shows a mixture of both cold and hot material and exhibits a helical motion, often seen in coronal mass ejections (Sheeley & Wang 2007) and AR flares (Zhang et al. 2022). These characteristics suggest that the relations between solar eruptions at different scales show quite similar features. Besides, ideal experiments of magnetic reconnection have also been performed. Farid et al. (2022) aimed to check the presence of a mini-filament. They used the flux rope insertion method to create an energized NLFFF model of coronal jets and identified the consequent eruption as an untwisting jet. Pariat et al. (2015) demonstrated that for getting a helical jet it is a primary requirement to store sufficient magnetic free energy and impulsively release it. They thought that the pre-eruption magnetic configuration might not be so important. Wyper et al. (2017), on the other hand, proposed that the magnetic breakout is a universal model for solar eruptions including coronal jets. However, current simulations cannot spontaneously produce a mini-filament structure directly from observed magnetograms.

In this paper, we analyze the motion of a jet through multiwavelength observations from AIA and IRIS. Then, we simulate the AR with a zero-β MHD simulation using the Message Passing Interface Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC; Keppens et al. 2003, 2012; Porth et al. 2014; Xia et al. 2018) and find the structure of a null point and a flux rope at the jet base. In Section 2, we show the observations of the jet and depict its kinematic evolution. The extrapolation and simulation methods are also described. Section 3 shows the results, and Section 4 presents a summary and discussion about our simulation. We show the importance of the twisting flux rope in the reconnection of this jet event.

2.1. Instruments

SDO/AIA provides regular full-disk EUV observations at several different wavelengths with a spatial resolution of 06 and a temporal cadence of 12 s. The AIA channels respond to a wide range of temperature (0.1–10 MK) and can be used to explore solar features over a wide range of scales.

IRIS is widely used to detect the transition region between the chromosphere and corona. It switches to different raster modes and obtains spectral data with different resolutions. In addition, it can provide high-resolution slit-jaw images (SJIs) in several wave bands, such as C ii λ1330, as a supplement to the SDO/AIA observations. The C ii λ1330 SJIs cover a specific region with a field of view of 60'' × 68'' and a pixel size of 016 × 032. We also adopt simultaneous spectroscopic observations right above the reconnection site to analyze the Doppler shift in this region with the Mg ii k lines, which form at chromospheric temperatures. The 1330 SJIs are used for the comparison of the structures of AIA.

The magnetogram used in the magnetic field extrapolation and numerical simulation comes from the Spectro-Polarimeter (SP; Lites et al. 2013) of the Solar Optical Telescope (SOT; Ichimoto et al. 2008; Shimizu et al. 2008; Suematsu et al. 2008; Tsuneta et al. 2008) on board Hinode (Kosugi et al. 2007). It was obtained by scanning this region step by step along the east–west direction in a fixed band of wavelengths centered on the Zeeman-sensitive Fe i lines at 6302 Å. The observation started at 01:02:05 UT and ended at 01:59:30 UT, right before the jet event, which started at around 02:03 UT. The whole magnetogram covers a field of view of 270'' × 164'' with a pixel size of 032 along the slit and a scanning step size of 03. To compensate the stretch caused by the long-time scanning, we scale its spatial resolution with the assistance of the magnetogram from the Helioseismic and Magnetic Imager (HMI; Scherrer et al. 2012; Schou et al. 2012) on board SDO by comparing the main common features in the magnetograms, which have already been corrected for projection effects. The pixel size of the HMI magnetogram is 05. The position information between the instruments we used is shown in Figure 1.

Zoom In Zoom Out Reset image size

2.2. Observations

The jet event started at around 02:03 UT on 2019 March 22 in NOAA AR 12736 and lasted for several minutes. This event has been intensively studied from an observational point of view (Yang et al. 2020; Joshi et al. 2020a, 2021a, 2021b; Schmieder 2022; Schmieder et al. 2022). Multi–wave-band observations from SDO/AIA recorded this intense eruption in all its wave bands. We select four EUV images and one UV image provided by SDO/AIA as shown in Figure 2, which reflect different temperatures in the chromosphere and corona. At the EUV brightness peak time around 02:06 UT, the helical motion of the jet around its axis can be recognized in most of the wave bands through the bright and dark patterns. Meanwhile, IRIS observed the AR, and one of its slit positions was just above the reconnection site of the jet base at the early eruption phase. Figure 2 shows some observational features, such as the shape of the jet base. An X-shaped reconnection site and a flux rope can also be inferred by the heated brightening structures as shown in Figure 2. SOT/SP on board Hinode obtained magnetic field on the photosphere from 01:02 UT to 01:59 UT by scanning this region with a spectrograph. A potential field can be extrapolated when the magnetic field is in a stable or quasi-stable state. Usually it is valid doing so before eruptions occur. As the magnetogram shows, this magnetic field possesses a major positive polarity at northeast and a major negative polarity at southwest, as well as several parasitic negative polarities around the positive one. A fan-spine structure may easily form in such a magnetic configuration around those parasitic polarities. In this case, the fan lies above main polarities in the east, and the outer spine connects to the far side in the west region.

Download figure:

Video Standard image High-resolution image

We focus on the small field of view of the jet base using the SDO/AIA full-disk observations. To study the propagation of the twisting jet, we select a slit along the direction of the jet indicated as the solid white line in Figure 3(a) and present the time–distance map of the slit in Figure 3(b). The figure shows the evolution of the jet spire, a bulk motion of bright plasma with a velocity of ∼300 km s⁻¹ followed by a motion of dark structure at around 02:04 UT. This suggests that the hot outflows set in through magnetic reconnection at first, followed by the eruption of the cold material, which may come from the flux rope in the jet base, and move along the reconnected field lines. The curved pattern in Figure 3(b) also shows the helical motion of the jet, indicating that a twist structure, probably a flux rope, may exist previously at the jet base. To further look into the helical motion of the jet, we draw five equi-distanced cuts perpendicular to the slit in Figure 3(a). Figure 4 displays the time–distance maps of the five cuts. The eruptive motion of the bright and dark patterns, shown by slices in two perpendicular directions in both Figure 3(b) and Figure 4, and the recurrent features, labeled by dashed white lines in Figure 4, exhibit the helical motion of the jet. The panels of Figure 4 in sequence also show the separation between these two patterns. The dark wave feature in the last panel of Figure 4, which is sliced from the transverse direction, clearly shows an oscillatory motion in the eruption.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

IRIS performed spectral observations with a four-step raster for the region covering the reconnection site during the whole process of the jet event. At the time of the jet eruption, the profile of the Mg ii lines shows a broad extension in both the red and blue wings, which indicates complex material motions in the transition region. The presence of complex motions has been separately verified by Joshi et al. (2020a). We adopt the moment method to calculate the average Doppler velocity regardless of the irregular profile. We mainly analyze the Mg ii k lines of the first slit position, which is located just above the reconnection site during the peak time. We thus obtain the Doppler velocity from observations around the peak time and then compare it with the numerical values calculated from the data-constrained simulation in Section 3.

2.3. Magnetic Field Extrapolation and Data-constrained Simulation

Based on the vector magnetic field obtained from SOT/SP, we make an NLFFF extrapolation to study the magnetic topology structure in AR 12736. To further investigate the relationship between the magnetic structure and observational morphology, we scale the spatial resolution of the SOT/SP magnetogram by comparing it with the SDO/HMI magnetogram that has already been corrected for the projection effects. Meanwhile, we resize the original data to reduce the spatial resolution by twice for relieving the computational burden. The NLFFF is reconstructed from the magneto-frictional module implemented by Guo et al. (2016a, 2016b) in MPI-AMRVAC. It can relax the potential field together with the bottom boundary to NLFFF in the computation domain. The method can be applied to multiple situations, such as Cartesian or spherical coordinates with uniform grids or adaptive mesh refinement grids. We note that it is hard to attain a completely force-free field relaxed from an observed magnetogram because the Lorentz force deduced from observations is not nicely balanced in the photosphere owing to the role of gravity (Zhu et al. 2016). Hence, the NLFFF might not be in a complete equilibrium state initially.

Figures 5(a) and (b) show the vertical magnetic field and the distribution of current density magnitude in the extrapolated field 1.2 Mm above the photosphere. We zoom in on the region of interest showing a current density isosurface of 0.03 A m⁻² and the distribution of J_z in Figures 5(c) and (d), respectively. The concentrated current reveals a strong magnetic flux rope structure in the vicinity of the small parasitic negative polarities in the north. It is the same structure detected in the HMI vector magnetograms by Joshi et al. (2020a). And the elongated double-peaked J_z pattern with hooks confirms the existence of the flux rope (Aulanier & Dudík 2019; Joshi et al. 2020a). Figure 6(a) displays two magnetic structures in the extrapolated field, a flux rope and a null point. Figure 6(b) shows the flux rope lying above a bald patch region low in the photosphere. By tracing the magnetic field lines of the flux rope, we find that they are actually a part of the bald patch field lines. Figure 6(c) shows the null point in the south above the flux rope. The X-shaped reconnection site is reproduced with the spine connecting to a remote region in the west.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To probe the nature of the jet, we conduct a data-constrained MHD simulation based on the result of the NLFFF extrapolation. The model we adopt is the same as that in Guo et al. (2019), which solves the zero-β MHD equations. For the boundary conditions, we select the data-constrained case, which fixes the bottom boundary, rather than the data-driven case, because of the limitation of the scanned magnetogram from SOT/SP. The MHD equations are solved in the dimensionless forms in MPI-AMRVAC. In this model, the gravity and gas pressure gradient are omitted in the momentum equation, and the energy equation including the heat conduction term and others is dropped entirely. In other words, we only solve the continuity equation, momentum equation, and magnetic induction equation to get the density, velocity, and magnetic field.

The initial condition is set to be the relaxed NLFFF, and the initial density distribution is calculated by assuming a temperature distribution of a stepwise function of height. We set the computation grid to be 288 × 200 × 200, namely 166 × 85 × 85 Mm³ in real scale with resolution of 0.58 × 0.42 × 0.42 Mm³. In addition, an anomalous resistivity is added in this model to facilitate the occurrence of magnetic reconnection. When the current density J of a local point is higher than the critical value J_c , the resistivity η is set to be $\eta ={\eta }_{0}\times {(J-{J}_{c})}^{2}/{{J}_{c}}^{2}$ and otherwise to be zero. Here we set η₀ and J_c to be 5 × 10⁻⁴ and 1000 in normalized units, respectively. The whole simulation time is about half an hour. We output the snapshots of the simulation results every 24 s for further comparisons with the AIA and IRIS observations.

In the NLFFF extrapolation, we recognize a flux rope low in the photosphere related to the north parasitic polarities and an X-shaped current sheet high in the corona in a south region between the two main polarities, which is also the reconnection site as shown in Figure 6. At the early time of the simulation, the eastern footpoint of the flux rope slowly drifts close to the positive polarity in its southeast. Meanwhile, the main body of the flux rope rotates around its axis counterclockwise and approaches the reconnection site located high in the south region. The rotation slowly changes the magnetic configuration, and the kinetic energy of plasma is accumulated through the continuous reconnection. When the field lines of the flux rope encounter the X-shaped current sheet, they join the reconnection process with other field lines within the sheet. Then, the explosive burst takes place and the jet with helical motion comes out of the null point along the spine connecting to the faraway side in the major negative polarity. At last, the field lines shrink to a low magnetic arch structure after the burst. Thus, we can detect the hot reconnection outflows and later the cold material from the flux rope, corresponding to the bright and dark patterns in AIA observations, respectively. The general evolution is revealed by four typical moments shown in Figure 7.

To compare the simulation with the observations of SDO/AIA and IRIS, we project the simulation to the line of sight and superpose it on the AIA observations in Figure 8. Since the magnetogram from SOT/SP is scanned through a period of time, we could not give an exact time that the simulation starts from. Thus, we use T₀ to represent the start time and focus on the time evolution of simulations and observations. The simulation matches the observations quite well in the perspective of the triangular shape at the base. The field lines of the flux rope fit with the anemone-like arcade, and the field lines at the null point display a triangular shape like the jet. Note that only magnetic field lines are displayed in Figure 8, and the plasma is assumed to move along the field lines. We find that the time-sequential images match roughly synchronously with the real-time evolution revealed in the observations. All these results indicate that the simulation reproduces the observations to a satisfactory degree.

We further extract the velocity and density information from the simulation snapshots. Then, we analyze the Doppler velocity along the line of sight and compare it with that from observations. The results are shown in Figure 9. The Doppler velocity from the simulation is derived from the Gaussian fitting based on the distribution function of density squared versus velocity along the line of sight as shown by each white line in Figure 9(b). Note that every white line is perpendicular to the IRIS slit and can be projected to one point on the slit. Then, we can obtain the velocity distribution from simulation at these points along the slit. We discard the extremely large values compared with the maximum absolute value from observations and further make a spatial interpolation to the interval of observations. This is intended for calculating the cross-correlation between Doppler velocities from observations and simulations. Because of the lack of a lot of information about specific matters in the zero-β simulation, we take such an averaged method in analyzing the Doppler velocities. For the observational result at a selected time, we compare it with each simulation snapshot. And the best result with a correlation coefficient of 0.71 is shown in Figure 9(e). The time difference between the observation and simulation is about 6.5 minutes for this best-matched case, which is acceptable considering the simplified MHD equations used in the simulation.

Zoom In Zoom Out Reset image size

This jet event has already been studied by two different groups, as was mentioned in Section 2. Different opinions were proposed regarding the trigger mechanism of the jet event. Yang et al. (2020) did an NLFFF extrapolation to derive the main magnetic structure of this event. They confirmed a fan-spine magnetic topology and found two flux ropes in the extrapolated results. In addition, they found a null point responsible for the reconnection site, which is the same as what we find in this work. In their opinion, the breakout-type reconnection caused by one of the flux ropes leads to the eruption of the jet. By contrast, Joshi et al. (2020a) checked the magnetogram evolution prior to the occurrence of the event and analyzed the vector magnetic field map in comparison with an ideal MHD simulation. They found a flux rope along the polarity inversion line in the north, and the flux rope migrates toward the south owing to the photospheric motions, leading to the formation of a small parasitic bipole. The rotational motion of the flux rope was verified by them through analysis on the IRIS spectra in their Figure 8(b). Though not so obvious, the rotational characteristics can also be exhibited in our Figure 9(d). A strong extended bidirectional shifts at around y = 75, and from y = 75 blueshifts fade along the y-axis, while redshifts fade against the y-axis. They found that the bipole is also a bald patch region, above which a current sheet forms. Thus, the twist of the flux rope is transferred into the jet through reconnection without direct eruptions. These two studies come into a consensus that the helical motion of the jet comes from the preexisting twist in the flux rope. Nevertheless, the critical difference of these two studies lies in how the jet erupts, namely, how is the twist transferred into the jet. Such a difference could be ascribed to the different methods used in these two studies, since the study by Yang et al. (2020) is based on the magnetic field extrapolation together with a high-resolution observation from NVST, while that of Joshi et al. (2020a) focuses on the evolution of the magnetogram and comparisons with ideal MHD simulations. Neither of them performed dynamic simulations about this AR. Thus, we make a zero-β MHD simulation in this work to probe how the jet is triggered. The Doppler velocity of our simulation is in accordance with IRIS observations from the aspect of correlation analysis. Specifically, our simulation shows that the flux rope lies low at the bald patch and the null point is located high in the corona. However, they are not far from each other when observed from the viewing angle of the spacecraft. In that sense, it is hard to distinguish their distinct features only according to the spectrum that is actually an integration along the line of sight. Therefore, pure analysis of observations might not lead to a solid conclusion. The combination of MHD simulations and observations is needed for learning the physical nature of solar eruptive events.

With the simultaneous observations from AIA and IRIS, we analyze the evolution of the jet in NOAA AR 12736 on 2019 March 22. We draw a time-slice map for the slit along the jet axis and find that the jet is twisting around its axis as it propagates outward at a speed of ∼300 km s⁻¹. The jet erupts with hot matter followed by cold matter. We suppose that they correspond to the hot reconnection outflows through a null point and the cold material from the flux rope. The helical motion comes from the untwisting of the flux rope. Furthermore, we utilize an NLFFF to check the topological structure of the magnetic field with the photospheric vector magnetic field from SOT/SP. A null point structure in the south and a flux rope in the north of a bald patch region are found to exist in the AR. Then, we adopt a zero-β MHD simulation based on the results of the NLFFF extrapolation. It reveals that the jet may be produced by the magnetic reconnection at the null point. With the untwisting motion of the flux rope and slipping motion of its footpoint, the inner field lines of the flux rope encounter and then reconnect with the field lines originally situated at the null point. The reconnection outflows propagate along the spine line connected to the major negative polarity in the west region, which correspond to the bright and dark patterns in the AIA images.

We experiment with both the SDO/HMI magnetogram and SOT/SP magnetogram and the reconstruction of the NLFFF for further MHD simulations. An X-shaped current sheet was found in both cases; however, the flux rope can only be found in the SOT/SP case. Although we adopt the same magnetogram as that used in Yang et al. (2020), the magnetic structures obtained are different in the two cases. We attribute the different results to the usage of different extrapolation methods. Moreover, in our trials without extra anomalous resistivity, the rotational movement of the flux rope near the photosphere cannot be reproduced in the simulation. The resistivity is artificially added in the regions where the current density is large and thus the flux rope reconnects under high resistivity and reveals untwisting motions. This implies that the rotational untwisting motion of the flux rope at the beginning comes from the reconnection in the flux rope. Furthermore, the ideal MHD simulations conducted by Wyper et al. (2018) show that minority-polarity intrusions may add free energy into the field and can lead to both intermittent low-level reconnection and explosive, high-energy-release reconnection above these regions. Though the moving of the parasitic polarity is not considered in our data-constrained simulation, this may partly account for the impacts of the resistivity added in our experiments.

Our MHD simulation does not contain heat conduction and radiation process. Therefore, the physical parameters we obtain in the numerical simulations are not quantitatively accurate. Recently, Zhou et al. (2022) conducted an MHD simulation of a sheared arcade configuration, in which a magnetic flux rope is formed and erupts through reconnection. The flux rope in their simulation rotates against its twist and transfers the twist into the surrounding field lines, which resembles our results. However, further studies are needed for clarifying the mechanism of the rotation of the flux rope leading to the reconnection between the flux rope and the spine structure connected to the remote side of the AR. In particular, we should focus on the role of the anomalous resistivity in driving the reconnection in the inner part of the flux rope and leading to its rotational motion around its axis and its upward motion in the chromosphere. We should also improve the data-constrained simulations, such as adding some specific magnetic structures and photospheric motions in order to fit the observations more closely. On the other hand, we also aim to carry out simulations with full MHD equations including heat conduction and radiative cooling to learn the physics more accurately in such jet events.

J.H.Z., Y.G., and M.D.D. are supported by the National Key Research and Development Program of China (2022YFF0503004, 2021YFA1600504, and 2020YFC2201201) and NSFC (11773016 and 11733003). We thank the two anonymous referees very much for their constructive suggestions. We deeply appreciate the free access to the magnetogram and observation data provided by SDO/AIA, SDO/HMI, SOT/SP, and IRIS science teams. J.H.Z. thanks Ye Qiu and Ze Zhong for their help in the work. Computational resources were provided by the High Performance Computing Center (HPCC) at Nanjing University.