The 3D magnetic topology and plasma dynamics in open stochastic magnetic field lines

The plasma disruption in tokamaks is one of the most critical problems that must be overcome for a successful fusion reactor, because it can cause excessive heat and current loads on the walls resulting in irreversible damage. In general, there are two distinct phases in the plasma disruption: the thermal quench (TQ) and the current quench (CQ). The TQ is a sudden collapse of the plasma thermal energy in a device, and its mechanism can differ depending on the type of TQ. In the case of TQ triggered by locked modes,^1,2 magnetic islands growing in a non-rotating plasma can overlap each other creating a stochastic magnetic field over a very wide radial range, which could lead to TQ. Since the stochastic field lines are connected to the wall boundary in this case, the so-called open stochastic magnetic field plays an important role in enhancing electron heat transport and thus decreasing electron temperature rapidly. The underlying mechanism of TQ triggered by locked modes has not been clearly revealed yet due to the complexity of the three-dimensional plasma dynamics in the open stochastic magnetic field lines.

Plasma transport in stochastic magnetic fields has been a long-standing problem of particular interest because stochastic magnetic fields are often present in various plasma experiments (e.g., reversed field pinch³ and stellarator^4,5) and various conditions [e.g., resonant magnetic perturbation (RMP)^6,7] The behavior of complex trajectories of stochastic field lines is usually described in terms of a diffusive process.^8,9 In this context, plasma transport can be characterized by the stochastic field line diffusion and particle collision-induced cross field-line decorrelation, assuming the dominant parallel plasma dynamics along each field line. Based on this simple physical picture, plasma transport in the stochastic field was characterized in various parametric regimes. In the case of toroidal devices, due to the non-uniform magnetic field, some electrons get trapped by magnetic mirror forces and do not exhibit stochastic behavior.^10,11 Considering the consistent coupling between electron and ion dynamics, the ambipolar potential must develop in the stochastic layer in order to maintain the quasi-neutrality of the plasma. A simplified one-dimensional model was developed to take into account the effects of the ambipolar potential on the plasma transport.¹² Plasma turbulent transport in the stochastic magnetic field has also been an interesting topic of recent theoretical studies.^13,14

Although various models have been developed to characterize plasma transport in the stochastic field, these models cannot be directly applied to the situation of TQ. This is because the typical analytic model does not take into account the open stochastic magnetic field line connected to the wall within the finite length that is important in TQ. The behavior of open stochastic field lines cannot be simply characterized by the concept of diffusion because field line trajectories could hit the wall before exhibiting sufficient diffusive nature. The finite length of the open field line, the so-called connection length, has a complex three-dimensional structure depending on the characteristics of the magnetic perturbations that produce the stochastic field. It is worth noting that the 3D structure of the connection length is also observed in RMP experiments.^15,16 However, the stochastic layer due to RMP is much thinner and localized in the edge pedestal region, unlike a global stochastic layer presented in the TQ situation. Due to particle loss along the 3D open field lines with different lengths, the ambipolar potential must have a 3D structure to balance the electron and ion losses along each field line. In the RMP experiments, the existence of the 3D ambipolar potential associated with the edge stochastic layer and remnant islands has been verified and studied in various devices.^17–20 

The 3D ambipolar potential is more important in the TQ situation because the stochastic layer of TQ is much wider than that of RMP. As a very wide stochastic layer needs the 3D ambipolar potential on a global scale for quasi-neutrality, we must consider its macroscopic effects on the overall plasma transport. However, in many thermal quench studies, the 3D magnetic topology and associated ambipolar potential have been often ignored or oversimplified due to the complexity. It is also generally assumed that parallel plasma transport is dominant in the system, resulting in a very high parallel conductivity $χ ∥$ and a very low perpendicular conductivity $χ ⊥$⁠, i.e., $χ ∥ ≫ χ ⊥$⁠. In our previous study,²¹ we proposed a comprehensive physical picture of the collisionless electron dynamics in the open stochastic field, taking into account the 3D magnetic topology and the self-consistent ambipolar potential. It was found that open stochastic field lines act as 3D magnetic mirrors, and perpendicular E × B transport and mixing effects induced by the 3D ambipolar potential are significant. In particular, the 3D magnetic mirror plays a crucial role in determining electron temperature in the system. Considering the 3D magnetic mirror and the E × B mixing effects, we discovered a novel mechanism that effectively decreases the electron temperature due to the selective loss of high-energy electrons.

In this study, we newly elucidate ion dynamics in open stochastic magnetic fields, which have not been studied in depth previously. In a similar manner to understanding the electron dynamics, the combination of 1D parallel dynamics and 3D cross field effects provides insights into ion kinetic physics. As the initial plasma is hot and collisionless at the onset of the thermal quench, this work aims to understand collisionless mechanisms of plasma transport in stochastic magnetic fields for the earlier stage of the thermal quench. In particular, we investigate in detail how the 3D magnetic topology relates to electron and ion temperature decrease rates. For this purpose, we examine two different magnetic configurations to study how the electron and ion dynamics change with the 3D magnetic topology. The 3D kinetic simulations based on a global gyrokinetic model^22,23 have been performed on two magnetic configurations. We emphasize the importance of the 3D magnetic mirror and the ambipolar potential that determine electron and ion temperature decrease rates, by analyzing and comparing 3D magnetic topology and simulation results.

Section II describes the 3D topology of open stochastic field lines with two key concepts, the connection length L_c and the effective magnetic mirror ratio $M eff$⁠. Section III illustrates the comprehensive physical picture of the collisionless electron dynamics in the open stochastic field lines step-by-step. Section IV explains the ion dynamics by focusing on the roles of magnetic drift motion and ambipolar potential. Section V discusses the difference between electron and ion temperature evolution. Section VI summarizes and concludes the study.

The magnetic islands can be produced at resonant flux surfaces when resonant magnetic perturbations are caused by internal instabilities or external sources. If the radial widths of magnetic islands are large enough to overlap with adjacent islands, the stochastic magnetic field is generated in the overlapping region.^24,25 If that region extends beyond the last closed flux surface (LCFS), the stochastic field lines wander inside the plasma with complex 3D trajectories and eventually hit the wall boundary as shown in Fig. 1(a). The opening of the stochastic magnetic field lines toward the wall can be observed in RMP experiments or the thermal quench events triggered by locked modes. Whereas a typical RMP generates open stochastic field lines within a very thin layer at the edge region, even the partial thermal quench can have a very wide global stochastic layer, from half the minor radius to the LCFS, as observed in recent DIII-D experiments.¹ It is important to understand the 3D topology of open stochastic field lines since it determines the overall plasma dynamics of those systems. For this purpose, we analyze the 3D magnetic topology in terms of two key concepts: the connection length L_c and the effective magnetic mirror ratio $M eff$⁠. In particular, we will highlight the importance of $M eff$ because the 3D magnetic mirror strongly influences electron heat transport and temperature in the system.

We set up two different magnetic configurations that have very wide stochastic layers to mimic the thermal quench situation triggered by locked modes. For this purpose, we apply prescribed static magnetic perturbations to the “cyclone base case”²⁶ equilibrium with an artificial circular wall boundary. The prescribed magnetic perturbation δB of multiple harmonics is defined using the perturbed parallel vector potential,^27–29 

where $δ A ∥$ is the perturbed parallel vector potential, $B 0$ is the equilibrium magnetic field, α is the perturbation parameter, $ψ t$ is the square root of toroidal flux as the radial location, $Γ ( m , n )$ is the radial shape function for each mode, ζ and θ are toroidal and poloidal angles in the PEST coordinate system,^30,n and m are toroidal and poloidal mode numbers, and $η ( m , n )$ is the phase angle for each mode. We chose multiple mode numbers resonant to the safety factor profile of the cyclone base equilibrium such as $( m , n ) = [ ( 2 , 1 ) , ( 3 , 2 ) , ( 4 , 2 ) , ( 5 , 2 ) , ( 5 , 3 ) , ( 6 , 3 ) , ( 7 , 3 ) , ( 8 , 3 ) ]$⁠. For simplicity, every mode has the same shape function $Γ ( m , n ) ( ψ t ) = ( 6 × 10 − 4 ) ψ t exp ( − | ( ψ t − 0.7 ) / 0.35 | 6 )$⁠. To focus on understanding plasma dynamics in open stochastic field lines, we apply strong magnetic perturbations that create a very wide stochastic layer from almost half the minor radius to the wall boundary without remnant magnetic islands as shown in Fig. 1(b). The two magnetic configurations have a different set of phase angles $η m , n$⁠. The first configuration, the same as studied in Ref. 21, has the uniform phase angle for all modes, i.e., $η ( m , n ) = π ∀ ( m , n )$⁠, so that the O-points of the magnetic islands are perfectly aligned in the outboard midplane at ζ = 0 poloidal plane. Conversely, in the second configuration, the phase angle is randomly chosen to represent a situation where the O-points are not aligned: $η ( m , n ) = [ 1.5 , 6.2 , 0.2 , 3.4 , 2.8 , 4.3 , 5.5 , 1.9 ]$ where the order of values corresponds to that of the mode numbers described above. Global gyrokinetic particle simulations were conducted using the GTS^22,23 code to study plasma transport in these systems. More details about the numerical models can be found in our previous work.²¹ 

It is worth noting that there is a key difference in the topological understanding of open stochastic magnetic field lines as opposed to an internal stochastic layer generated by overlapping internal magnetic islands. In the internal stochastic layer, the stochastic magnetic field lines wander infinitely without termination. Their excursions in the radial direction and associated particle transport can be characterized by considering the radial diffusivity of the stochastic fields^8,9,12 or the subdiffusive nature of chaos due to a bounded radial domain.³¹ On the other hand, an open stochastic magnetic field line has a finite length, the so-called connection length L_c, due to the presence of the wall boundary that terminates the excursion of the stochastic field line. The connection length of the open stochastic field line is directly related to the particle loss toward the wall and hence plays important roles in plasma dynamics, which will be described in Sec. III. For open magnetic field lines, the concept of diffusivity is ambiguous because the calculation of diffusivity is less likely to converge while tracking the set of open magnetic field lines for finite length L_c. Moreover, the connection length of each open field line can have completely different values in a very wide range, from very short (⁠ $L c ≲ 10 m$⁠) to much longer (⁠ $L c ≳ 1000 m$⁠), depending on its 3D trajectory. For example, Fig. 1(a) shows that two field lines passing through two very close points 5 cm apart on ζ = 0 poloidal plane (denoted by glowing small spheres) have about 12 times different connection lengths, such as $L c ∼ 45 m$ (blue) and $L c ∼ 550 m$ (yellow), respectively. This bifurcation nature of the connection length cannot be described by the typical concept of the average diffusivity of stochastic field lines. Therefore, we have to consider the individual 3D trajectories of stochastic field lines in order to properly understand the plasma dynamics in such a system.

For this purpose, we construct and analyze the three-dimensional structure of the connection lengths by calculating the volumetric average of L_c of multiple field lines based on spatial 3D grid nodes with high resolution. Figure 1(c) shows the 2D structure of the connection length on the poloidal plane ζ = 0. It is worth noting that the shorter and longer L_c regions are intertwined in complex ways to form specific patterns. The shorter L_c regions in blue color show the fingerlike structure extending from the wall boundary to the inner stochastic region. The shape and radial length of fingerlike structures depend on the mode numbers, amplitudes, and phase angles of magnetic perturbations that are resonant at the edge region. For example, Figs. 3(a) and 3(c) show the connection length of two magnetic configurations at ζ = 0 poloidal plane, in which the uniform $η ( m , n )$ configuration shows longer connection lengths than the random $η ( m , n )$ configuration. Field lines passing through the shorter L_c region (blue) hit the wall very quickly before spreading enough to show chaotic trajectories. On the other hand, field lines passing through the longer L_c region (yellow) wander inside the stochastic layer with more chaotic trajectories and eventually hit the wall. It is worth noting that this kind of pattern is also observed in the RMP experiments in which the L_c structure is often distinguished into laminar and ergodic zones.^32,33

Since the connection length has a complicated structure with sharp gradients between distinct L_c regions, the 3D magnetic topology cannot be simply characterized by some one-dimensional quantities such as the radial profile of average connection length. Instead, we directly deal with the 3D magnetic topology to understand associated plasma dynamics in the system. In particular, the sharp gradients between different L_c regions play important roles in the evolution of the ambipolar potential and the generation of E × B mixing effects, which will be described in Sec. III.

Whereas the concept of connection length has been very widely adopted and used to study plasma transport in the open stochastic field, the magnetic mirror force, $b · ( − μ ∇ B )$ where b is the unit vector of the magnetic field and μ is the magnetic moment of the charged particle, has been often neglected or oversimplified. However, in our previous study,²¹ we demonstrated the importance of three-dimensional magnetic mirror effects that strongly influence electron heat transport in open stochastic fields. First, we need to understand how open stochastic field lines act as the three-dimensional magnetic mirror. The small magnetic perturbation $δ b$ barely affects the magnitude of local mirror force, i.e., $b · ( − μ ∇ B ) ≈ ( b 0 + δ b ) · ( − μ ∇ B 0 )$ where the subscript 0 stands for the equilibrium magnetic field. Interestingly, the small change in the direction of the magnetic mirror force results in spatially varying magnetic potential (⁠ $μ B 0$⁠) along 3D trajectory of the stochastic field line between two wall endpoints as shown in Fig. 2. This implies that the charged particle cannot flow freely along the field line and its parallel motion is strongly influenced by the magnetic potential structure. For example, the charged particle can be magnetically trapped between two bouncing points along the 3D trajectory of the stochastic field line. Therefore, it is important to understand the key role of the 3D magnetic mirror in determining the electron confinement and heat transport, which will be described in Sec. III.

Although the detailed structure of the 3D magnetic mirror along the stochastic field line can be complicated as shown in Fig. 2(b), its overall configuration can be simplified by a few key variables defined at each position x along the field line

where $x w ± ζ$ is the position of wall endpoints in $± ζ$ directions and l is a parametric variable along the field line. It should be noted that the values of $B max + ζ ( x )$ are not constant along the field line and depend on the measuring position x even within the same field line. Based on these variables, a novel concept of effective magnetic mirror ratio $M eff$ at position x is defined as follows:²¹ 

The $M eff$ is a very useful concept to understand the 3D magnetic mirror effects in a concise manner. Figures 2(b) and 2(c) show that the overall structure of magnetic potential along the field line can be distinguished into either of two configurations depending on the value of $M eff$⁠: a magnetic well (⁠ $M eff$> 1) configuration and a magnetic hill configuration (⁠ $M eff < 1$⁠). In the case of the magnetic well configuration, some of the charged particles can be magnetically trapped between two bouncing points due to the magnetic mirror force. The value of $M eff$ greater than one indicates the capability of trapping particles; the higher $M eff$ the deeper well and more trapping as shown in Fig. 2(c). On the other hand, in the case of the magnetic hill configuration, the magnetic mirror force accelerates particles toward the wall boundary so that no particles are magnetically trapped. The value of $M eff$ less than one describes how strong the magnetic mirror force accelerates the particle toward the wall; the smaller $M eff$ the stronger acceleration. This acceleration toward the wall by the magnetic hill plays a crucial role in enhancing electron heat transport in the presence of the ambipolar potential, which will be described in Sec. III.

It is worth noting that multiple magnetic hills can exist within a single field line as shown in Fig. 2(b). In general, the presence and location of magnetic hills are highly dependent on the 3D trajectory of each field line and the measuring position within it. To understand the topology of $M eff$⁠, it is useful to construct the 3D map of $M eff$ in a similar way as the connection length analysis. Figure 2(a) shows the presence of magnetic hills (blue) with a complex pattern at the edge region. The structure of magnetic hills strongly depends on the characteristics of the magnetic perturbations resonant to the edge region. Figures 3(b) and 3(d) show the 2D structure of $M eff$ on the ζ = 0 poloidal plane for the uniform and random $η ( m , n )$ configurations, respectively. The uniform $η ( m , n )$ configuration has wider and stronger magnetic hills, whereas the random $η ( m , n )$ configuration thinner and weaker magnetic hills. This difference in the structure of magnetic hills results in a significant difference in the electron heat transport and the structure of electron temperature, which will be described in Sec. III C.

In this section, we will describe a comprehensive picture of plasma dynamics by considering the evolution of the 3D ambipolar potential in the open stochastic field line. Subsections III A–III C illustrate the comprehensive picture step-by-step: one-dimensional dynamics along each field line, three-dimensional effects across the field line, and a novel mechanism that efficiently decreases electron temperature by considering both 1D and 3D dynamics. As initial conditions for the GTS simulation, we assume the uniform plasma with a density of $1.6 × 10 19 m − 3$ and a temperature of $5 keV$ to focus on how quickly the electron temperature can drop in the presence of the open field line mimicking the thermal quench. The mean-free-path of the high-temperature electrons is about of the order of $10 km$⁠, which is much longer than the connection length $L c ≲ 1 km$⁠. Therefore, in this study, we focus on a collisionless mechanism of plasma transport in the complex 3D open stochastic field lines. The reference electrostatic potential of the system is zero at the wall boundary, which is the grounded conductor.

In this subsection, we illustrate the one-dimensional parallel plasma dynamics along each field line by assuming that each field line is independent of the other. First, we will briefly describe the basic dynamics of test particles in 3D open stochastic fields without considering any coupling between electrons and ions. Second, we will explain how the ambipolar potential evolves inside the stochastic layer to maintain the quasi-neutrality of the plasma.

The test particle simulation is a numerical experiment that generates some test particles within some region of interest and to traces their trajectories without considering any coupling between the test particles. This method is typically used to estimate the diffusivity or other statistical characteristics of particle trajectories which could be different from that of stochastic field lines due to the particle's drift motion and the Larmor radius effect. Unlike typical studies, we do not estimate the diffusivity but directly investigate the behavior of the test particles associated with the 3D magnetic topology such as the connection length L_c and the effective magnetic mirror ratio $M eff$⁠. It is noteworthy that the L_c and $M eff$ can rigorously illustrate the spatial-temporal evolution of the test electron population including its complex three-dimensional structures.²¹ The $M eff$ precisely determines whether the test electron at any 3D position x will be trapped by the magnetic mirror force according to the following trapping condition $| v ∥ ( x ) / v ⊥ ( x ) | < M eff ( x ) − 1$ where $v ∥$ and $v ⊥$ are the parallel and perpendicular velocities of the electron, respectively. The population of passing electrons decays rapidly along each field line with the characteristic confinement time $τ pass ∼ L c / ( 2 v e th )$ where L_c is the connection length of each field line and $v e th$ is thermal electron velocity. On the other hand, the trapped electrons can be well confined within the stochastic layer for a very long time unless the toroidal precession or collisions cause detrapping depending on the structure of $M eff$⁠. Although the test electron dynamics can be clearly understood by the 3D magnetic topology analysis, it should be noted that it does not reflect real-world physics since it does not consider any coupling between electrons and ions.

Positive ambipolar potential must develop inside the stochastic layer to balance the electron and ion losses for maintaining the quasi-neutrality of the plasma. Otherwise, the electron loss is much faster than ion loss due to the faster electron thermal velocity, resulting in nonphysical consequences of violating the quasi-neutrality of the plasma in the stochastic layer. Figure 4(a) describes how the positive ambipolar potential develops from the wall boundary toward the inner region along the field line. The propagation speed of the potential development is of the order of electron thermal speed since small charge perturbations causing the potential generation is carried by electrons. After a short time of the order of $τ ∼ L c / v e th$⁠, the positive ambipolar potential fully develops inside stochastic layer. The level of the ambipolar potential is determined dynamically to satisfy the ambipolar transport by accelerating ions toward the wall and decelerating or reflecting electrons.

As the plasma density and temperature near wall boundary decrease due to the parallel loss toward the wall boundary, the level of required ambipolar potential near the wall boundary also reduces correspondingly. The ambipolar potential drop propagates toward the plasma interior along the field line at the speed of the order of ion sound speed c_s, which is the characteristic speed of the ambipolar transport. After some time $Δ t$⁠, the ambipolar potential drop propagates along each field line to a distance of about $c s Δ t$ from the wall boundary. For shorter L_c lines, the propagation can travel the entire length of the field line during time $Δ t$⁠, i.e., $c s Δ t ≳ L c$⁠, thus significantly reducing the overall level of the ambipolar potential along the field line. On the other hand, for longer L_c lines, since the propagation has still only traveled a small portion of the field line from both wall boundaries, i.e., $c s Δ t ≪ L c$⁠, the most part of the field line can maintain the high ambipolar potential. As a result, the ambipolar potential becomes low at shorter L_c regions while it is still high at longer L_c regions as shown in Fig. 4(b). Therefore, the ambipolar potential will have a complex 3D structure that is strongly correlated with the L_c structure. The 3D effects by sharp gradients of the ambipolar potential across distinct L_c regions will be described in Subsection III B.

In the remainder of this subsection, we will focus on understanding the electron confinement along each field line by considering both the magnetic mirror and ambipolar potential effects. To understand the electron confinement, it is essential to identify whether the electron is passing or trapped along the 3D trajectory, because passing and trapped electrons have completely different confinement times. For simplicity, we assume that the ambipolar potential has a smooth and slowly varying profile along the field line like Fig. 4(b) such that there are no strong local fluctuations. This is a reasonable assumption to consider the average effect of the ambipolar potential when determining the passing-trapped condition for electrons. The basic role of the positive ambipolar potential in the stochastic layer is straightforward. Considering the zero reference potential at the wall (⁠ $Φ wall = 0$⁠), the positive ambipolar potential inside the stochastic layer acts as an electrical barrier that inhibits the escape of electrons. In the simplest case without considering magnetic mirror force, the electron will be reflected before hitting the wall boundary if the parallel kinetic energy of the electron at position x, $ε ∥ ( x )$⁠, is less than the height of the electrostatic potential barrier, i.e., $ε ∥ ( x ) < q e ( Φ wall − Φ ( x ) ) = e Φ ( x )$⁠, where $q e$ is the charge of the electron. Figure 5(c) shows schematic diagrams of the electron trapping by the ambipolar potential barrier near the wall boundary and the corresponding passing-trapped regions in the velocity space $( v ∥ , v ⊥ )$⁠.

In addition to the electric force, we must consider the magnetic mirror force together to get the correct physical picture. In general, we can define the total potential energy for the charged particle at the position x by summing the electrostatic and magnetic mirror potential energies

where q_s is the charge of the particle. Considering the total potential energy along the field line, the electron at the position x will be trapped if any of the following conditions is satisfied:²¹ 

where $B w eff ≡ min ( B w + ζ , B w − ζ )$ is an effective wall magnetic field and $Δ V$ is an effective difference in total potential energy between the current position x and the wall endpoints $x w ± ζ$⁠. Equation (7) is the pure magnetic trapping condition characterized by $M eff$⁠, which is still accurate for electrons with high magnetic moments or undergoing small ambipolar potential changes. Equation (8) is a global trapping condition between two wall endpoints without considering local changes in total potential energy along the field line. The $Δ V$ in Eq. (8) can be interpreted as the height of an effective total potential barrier that the electron must overcome in order to escape [see Figs. 5(b) and 5(c)].

It is important to note that the nature of the trapping conditions is highly dependent on the electron's magnetic moment and the magnetic potential configuration. In the case of the low-μ electron, regardless of the magnetic configuration, its trapping is likely determined by Eq. (8) with the dominant electrostatic potential, i.e., $ε ∥ < Δ V ≈ e Φ ( x )$⁠. This is similar to the case of the magnetic flat (⁠ $M eff ∼ 1$⁠) where the magnetic potential change is negligible as shown in Fig. 5(c). For higher-μ electrons, the magnetic mirror force plays an opposite role in the magnetic well and hill, respectively. In the magnetic well (⁠ $M eff > 1$⁠), higher-μ electrons are more likely trapped by the magnetic mirror force as shown in Fig. 5(b), since the higher magnetic moment increases the RHS term of Eq. (7). As a result, the magnetic well has a hyperbola-like passing-trapped boundary in the velocity space as shown in Fig. 5(b). On the other hand, in the magnetic hill (⁠ $M eff < 1$⁠), it never satisfies Eq. (7) regardless of μ because the RHS term is always negative, i.e., $( M eff − 1 ) < 0$⁠. It is worth noting that the height of the potential barrier $Δ V$ in Eq. (8) is reduced in the magnetic hill due to the negative second term of $Δ V$⁠, i.e., $B w eff − B ( x ) < 0$⁠. Since the second term of $Δ V$ is proportional to μ, the higher-μ electron is less likely trapped. Figure 5(d) shows that high μ electrons can become passing electrons with the help of high magnetic hills, and the passing-trapped boundary in the velocity space is the ellipse shape. This opposite nature of high-μ electrons in magnetic wells and hills results in a novel mechanism that is critical for determining the electron heat transport in complex 3D stochastic layer, which will be described in Sec. III C.

In Subsection III A, we assumed that the parallel plasma transport along each stochastic field line is the dominant mechanism and is independent of that of the other field lines. However, as the ambipolar potential evolves in the stochastic layer, we must consider the perpendicular E × B transport induced by the 3D structure of the ambipolar potential. Figure 4(c) shows that a sharp gradient of the ambipolar potential across shorter and longer L_c regions produces a strong electric field $E ⊥$⁠. The $E ⊥$ generates the E × B drift motion across the field lines. While the average part of $E ⊥$ drives a mean E × B flow mainly along the equipotential surface, the fluctuating part of $E ⊥$ induces turbulent E × B mixing between shorter and longer L_c regions. As the charged particles are actively interchanged between different field lines by the E × B mixing, the ambipolar potential structure along the field lines can fluctuate significantly, as shown in Fig. 4(c). In the 3D space, the mixing between higher $Φ$ at longer L_c regions and lower $Φ$ at shorter L_c region results in higher harmonic modes such as radial eddies. The overall structure of the ambipolar potential is still correlated with the L_c structure, since the one-dimensional parallel dynamics still operate along each field line on average. For example, Fig. 6 shows the GTS simulation results in the case of the random $η ( m , n )$ configuration. Figures 6(a) and 6(b) show the ambipolar potential and electron density structures. The overall structure is correlated with the L_c structure [see Fig. 3(c)] such that the ambipolar potential and electron density are higher at longer L_c and lower at shorter L_c regions. The radial eddies with higher poloidal harmonics are mainly generated between shorter and longer L_c regions. It is worth noting that the higher harmonic modes, or radial eddies, are created in the process of relieving sharp gradients of the 3D ambipolar potential on the macroscopic scale, and are not a kind of micro-instability that grows slowly from initial equilibrium.

In general, the E × B drift can enhance radial particle and heat transport toward the wall because the charged particles can be directly delivered into the radial direction by radial eddies across the field lines, often referred to as convective cells. In addition to the typical effects of convective cells, we found that the E × B mixing between magnetic wells and hills has a strong effect on decreasing electron temperature when compared to only the 1D parallel dynamics along field lines. The confinement of electrons along field lines is highly dependent on the magnetic configuration, especially for high-μ electrons as described in Sec. III A. The high-μ electron can be deeply trapped at the magnetic well, whereas in the magnetic hill, it can easily escape to the wall along the field line. Considering only these 1D parallel dynamics, the electron heat transport of the system will stagnate because a large amount of energy of high-μ electrons in the magnetic well regions is trapped and not emitted toward the wall for a very long time.²¹ Therefore, in order to properly understand electron heat transport in the system, three-dimensional effects must be considered together.

The E × B mixing between magnetic wells and hills provides an efficient three-dimensional mechanism for enhancing electron heat transport and decreasing electron temperature. The E × B mixing can transfer some high-μ electrons trapped in magnetic well regions to magnetic hill regions, where the high-μ electrons get detrapped and then have a chance to escape toward the wall along field lines. Although some high-μ electrons may return to the magnetic well regions by the E × B mixing before escaping, they can repeatedly enter the magnetic hill regions and have a chance to escape. Therefore, more and more high-μ electrons exit toward the wall over time. Figure 5 shows a schematic diagram of the collisionless detrapping. The high-μ electron initially trapped in magnetic well (⁠ $M eff > 1$⁠; red region) flows across different $M eff$ regions by the mean E × B and the E × B mixing. While the electron goes through a magnetic flat (⁠ $M eff ∼ 1$⁠; white region) and a weaker magnetic hill (⁠ $M eff ≲ 1$⁠; light blue region) regions, it is still trapped due to the dominant ambipolar potential barrier there. When it gets into a stronger magnetic hill (⁠ $M eff < 1$⁠; blue region) region at $x H$⁠, the electron gets detrapped and escapes toward the wall along the field line as shown in Fig. 5(d). It is worth noting that the collisionless detrapping only occurs for high-μ electrons and not for low-μ electrons. The collisionless detrapping mechanism decreases electron temperature more efficiently than convective cells because it selectively releases the thermal energy of high-μ electrons while trapping the thermal energy of low-μ electrons. In this sense, the collisionless detrapping mechanism could be the major cause of the rapid electron temperature drop in a strong and broad stochastic layer during the thermal quench.

The electron temperature structure is strongly correlated with the $M eff$ structure because the collisionless detrapping in magnetic hills is the main mechanism for lowering electron temperature. For instance, Fig. 6(c) shows that the electron temperature of the magnetic hill region is much lower than that of the magnetic well regions [see Fig. 3(d)]. It is worth noting that the low electron temperature in the magnetic hill leads to a decrease in the nearby electron temperature due to the E × B mixing effects. Therefore, the effectiveness of the collisionless detrapping and consequent electron temperature drop is highly dependent on the structure of $M eff$⁠, specifically how wide and long the magnetic hills are. Comparing two magnetic configurations, the uniform $η ( m , n )$ and random $η ( m , n )$⁠, clearly shows how effectively the collisionless detrapping decreases electron temperature depending on the structure of $M eff$⁠. As shown in Fig. 3, the uniform $η ( m , n )$ configuration has longer and wider magnetic hills (blue) than the random $η ( m , n )$ configuration. Therefore, the collisionless detrapping occurs much more frequently in the uniform $η ( m , n )$ case due to the deeper and longer magnetic hills. Figure 7 shows a large difference between two cases for the electron temperature decrease rate. Both cases show a sharp electron temperature drop during a very initial phase (⁠ $t ≲ 20 μ s$⁠) due to the loss of initial passing electrons. The uniform $η ( m , n )$ case shows a faster temperature drop to a lower value of about $4000 eV$ because it initially has more passing electrons due to the deeper and longer magnetic hills than the random $η ( m , n )$ case. As the ambipolar potential develops in the stochastic layer to maintain the quasi-neutrality, the decrease rate of electron temperature is greatly reduced. Since the majority of electrons get trapped by the ambipolar potential and magnetic wells, the collisionless detrapping becomes the dominant mechanism to decrease the electron temperature by releasing the thermal energy of the trapped high-μ electrons. The uniform $η ( m , n )$ case shows a steady temperature decrease at a rate of $− 540 eV / ms$⁠, comparable to the rate observed in the experiment.¹ Interestingly, the decrease rate of the uniform $η ( m , n )$ case is faster than the random $η ( m , n )$ case by more than a factor of 2 because the collisionless detrapping occurs much more frequently due to the deeper and longer magnetic hills. It is worth noting that the random $η ( m , n )$ case shows a slower temperature drop even though its connection length is shorter than that of the uniform $η ( m , n )$ case. This indicates that the effective magnetic mirror ratio could be more important than the connection length in determining electron heat transport and temperature change.

In this section, we describe the ion dynamics in the presence of open stochastic magnetic field lines. The ion dynamics has two key aspects that differ from the electron dynamics: (i) slower parallel motion and larger excursions by the magnetic drift motion and (ii) the opposite role of the ambipolar potential regarding the particle confinement along the field line. We will first explain the effect of the magnetic potential on ion dynamics without considering the electric field and then incorporate ambipolar potential effects.

Ions do not perfectly follow the field line due to their $∇ B$ and curvature drift motion that causes large excursions. Nevertheless, the concept of one-dimensional parallel dynamics along each field line is still useful for understanding ion parallel transport in open stochastic magnetic fields. The magnetic mirror force plays the same role in determining the ion confinement along the field line, as in the electron parallel dynamics. In the magnetic well region (⁠ $M eff > 1$⁠), some ions can be magnetically trapped by the reflecting magnetic mirror force as shown in Figs. 8(a) and 8(b). Conversely, in the magnetic hill region (⁠ $M eff < 1$⁠), no ions are magnetically trapped, so they can escape slowly toward the wall along the open magnetic field line. Unlike electrons, ions can easily move to other field lines due to their $∇ B$ and curvature drift motion. The ion trapped in the magnetic well region can be detrapped without collisions when it gets into the magnetic hill region by the magnetic drift motion. Figure 8(c) shows that some ions near the magnetically passing-trapped boundary are lost due to the collisionless detrapping induced by the magnetic drift motion. This collisionless detrapping by the magnetic drift motion is observed more strongly in ions than in electrons. Figure 9 shows the test particle simulation result of temporal evolution of magnetically passing and trapped ion densities at $ψ t = 0.7$⁠, without considering the electric field effects. The characteristic decay time or confinement time of the passing ion density is proportional to the connection length and inversely proportional to the thermal speed, i.e., $τ i , pass ∥ ∼ 0.5 L c / v i , th$⁠, which is typically hundreds of microseconds in this study. The trapped ion density decreases steadily due to the collisionless detrapping by the magnetic drift motion. The decrease rate of the trapped ion density is slower than that of the passing ion density, but it is faster than that of the trapped electron density [see Fig. 7(a) in Ref. 21].

Figure 10 shows the two-dimensional effect of the magnetic drift motion on ion transport. The magnetic drift motion (⁠ $v ∇ B + curv$⁠) is upward in the case of positive toroidal magnetic field as shown in Fig. 10(a). Figure 10(b) shows the magnetic drift motion in the PEST flux coordinate, which is directed radially inward at $0 < θ < π$ and radially outward at $− π < θ < 0$⁠. The opposite radial direction of the magnetic drift motion may cause some asymmetry in the ion transport and the associated plasma structures. The radially inward magnetic drift motion at $0 < θ < π$ accelerates the penetration of the negative ion perturbation from the wall boundary into the plasma interior. On the other hand, the radially outward magnetic drift motion at $− π < θ < 0$ impedes the penetration of the negative ion perturbation. As a result, the negative ion perturbation or decrease in ion density and temperature can penetrate slightly faster and deeper at $0 < θ < π$ than at $− π < θ < 0$⁠. For example, Figs. 10(c) and 10(d) show the test particle simulation results of the ion density and temperature at ζ = 0 poloidal plane for the uniform $η ( m , n )$ magnetic configuration. Although the uniform $η ( m , n )$ magnetic configuration has symmetric L_c and $M eff$ structures at the ζ = 0 poloidal plane, the ion density and temperature in the test particle simulation has slightly asymmetric structures. Both ion density and temperature are slightly lower at $0 < θ < π$ than at $− π < θ < 0$ because of the magnetic drift motion effect. This effect may be noticeable in the early stages, but it becomes less pronounced later as other effects related to complex 3D magnetic topology and ambipolar potential prevail.

In addition to the magnetic drift motion, the ambipolar potential is a key factor in understanding the ion transport and kinetic physics. As we described in Sec. III A and Fig. 4, the ambipolar potential must develop and propagate from the wall boundary to plasma interior in order to maintain the quasi-neutrality of the plasma in the stochastic magnetic fields. The ion and electron dynamics are strongly coupled through the ambipolar electric field. The ambipolar electric field accelerates and deaccelerates the propagation of ion and electron density perturbations, respectively. As a result, the plasma density perturbation propagates along the magnetic field line approximately at the ion sound speed. Figures 11(a) and 11(c) show the ion density evolution at $ψ t = 0.7$ for three different simulations in the case of the uniform and random $η ( m , n )$ configurations, respectively. The simulation considering only $E ∥$ (red) shows a faster decrease in the ion density than the test particle simulation (blue) due to the acceleration effect by the ambipolar potential. It is worth noting that the random $η ( m , n )$ configuration shows longer confinement times even though it has shorter connection lengths than the uniform $η ( m , n )$ configuration. This is because the random $η ( m , n )$ configuration has higher magnetic mirror ratios, so more ions are trapped, resulting in a slower confinement time.

The ambipolar potential strongly influences whether the ion is confined before escaping the wall boundary along the open magnetic field line. It is noteworthy that the role of the ambipolar potential for the ion confinement is opposite to that for the electron confinement. In the case of electrons, the ambipolar potential acts as an electrical barrier that traps electrons as described in Sec. III. On the contrary, the ambipolar potential accelerates ions toward the wall, helping them escape as shown in Fig. 8(d). Considering the 1D parallel dynamics, the ion confinement along the field line can be analyzed in a similar way to determining the passing-trapped conditions for electrons [Eqs. (7) and (8)]. However, it should be noted that determining the ion confinement is not straightforward as the electron confinement. This is because ions do not perfectly follow the magnetic field line due to the perpendicular drift motion, and the ambipolar potential can change during the long ion bounce time due to the slow parallel speed of the ions. Nevertheless, with some assumptions, the concept of the 1D parallel dynamics is still useful for obtaining insight into ion kinetic physics. Assuming the static ambipolar potential that smoothly decreases from the position x toward the wall along the field line, the ion at position x will be likely trapped if the following condition is satisfied:

where $x pk$ is the position of the nearest peak of the magnetic field from the position x (see Fig. 8) and $Δ V i$ is the height of total potential energy barrier for ions. Note that the second term of the right-hand-side, the contribution of the ambipolar potential, is negative so that $Δ V i$ decreases compared to the pure magnetic mirror case. In the magnetic well (⁠ $M eff > 1$⁠), more ions can be passing because the ambipolar potential lowers the total potential barrier as shown in Fig. 8(d). Figure 8(e) shows a hyperbolic passing-trapped boundary by considering both magnetic and electrostatic potential. Unlike the pure magnetic mirror case, ions with small magnetic moments can be passing ions regardless of the pitch angle. Figure 8(f) shows the only $E ∥$ simulation result of ion loss or negative ion perturbation (⁠ $δ f i < 0$⁠). The hyperbolic passing-trapped boundaries (black lines) are calculated from the simulation assuming that the electrostatic potential $Φ$ is zero at the point of the nearest peak of the magnetic field for simplicity, i.e., $Φ ( x pk ) ≈ 0$⁠. Unlike the test particle simulation [Fig. 8(c)], the significant ion loss is observed in the region of low magnetic moments and high pitch angles (blue contours under the black line) due to the ambipolar potential effect. Since more low-energy ions are lost compared to the pure magnetic mirror case, the ion temperature in the only $E ∥$ simulation becomes higher than the test particle simulation as shown in Figs. 11(b) and 11(d).

The E × B drift motion associated with 3D ambipolar potential affects the ion transport and kinetic physics. As explained in Sec. III B, the E × B drift motion enhances the average radial transport and mixes charged particles between the magnetic well and hill regions. In the case of electrons, the E × B mixing and the resulting collisionless detrapping in the magnetic hill region are the important mechanisms that can decrease the electron temperature steadily. In the case of ions, the E × B drift motion effect is more complex than the effect on electrons. Due to the slow parallel speed of ions, even the passing ions have the long confinement time of hundreds of microseconds to be lost to the wall boundary. The E × B mixing timescale $τ mix$ could be much shorter than the confinement times of both passing (⁠ $τ i pass$⁠) and trapped (⁠ $τ i trap$⁠) ions, i.e., $τ mix ≪ τ i pass < τ i trap$⁠. In these situations, the E × B mixing can cause a considerable collisionless trapping of ions. For example, some passing ions in the magnetic hill slowly escaping toward the wall may suddenly be transferred to the magnetic well by the faster E × B mixing, and the ions get trapped in the magnetic well. Since the E × B induces both collisionless trapping and detrapping, competing with each other, their effects on the ion transport can change depending on the situation. Figure 11 shows the effects of the E × B mixing on the ion density and temperature evolution for the uniform and random $η ( m , n )$ configurations, respectively. In the case of the uniform $η ( m , n )$ configuration, the $E ∥ + E ⊥$ simulation (green) considering the E × B drift motion shows slower decrease in ion density and temperature than the only $E ∥$ simulation (red) without considering the E × B drift motion. This implies that the collisionless trapping prevails over the collisionless detrapping in the case of the uniform $η ( m , n )$ configuration. On the other hand, the random $η ( m , n )$ configuration shows that the ion density and temperature decrease faster in the case of $E ∥ + E ⊥$ simulation, implying that the collisionless detrapping effect is more dominant. The detailed physics regarding the competition between collisionless trapping and detrapping of ions will be further investigated in future work.

It is worth discussing the difference between the ion temperature decrease rate with that of electrons. Figure 7 shows a comparison between the electron and ion temperature evolution in the case of $E ∥ + E ⊥$ simulations. At the initial stage, both electron and ion temperatures drop sharply due to the loss of initial fast passing particles. As the trapped particle population becomes dominant over time, the temperature decreases more linearly. Interestingly, for both magnetic configurations, the linear decrease rates of the ion temperature are much faster than those of the electron temperature. In the case of the uniform $η ( m , n )$ configuration, the ion temperature decrease rate of $− 1000 eV / ms$ is about two times faster than the electron temperature decrease rate of $− 540 eV / ms$⁠. For the random $η ( m , n )$ configuration, the ion temperature decrease rate of $− 1200 eV / ms$ is about five times faster than the electron temperature rate of $− 230 eV / ms$⁠. The reason for the faster ion temperature drop is that ions have greater radial and poloidal excursions due to the magnetic drift motion than electrons. In the case of electrons, the E × B drift motion is the dominant mechanism that carries electrons into magnetic hills, leading to the collisionless detrapping. For ions, in addition to the same E × B drift motion, the larger radial and poloidal excursions of ions allow a wider range of ions to enter magnetic hills. Hence, the ions have a higher chance of collisionless detrappings, which results in a faster temperature decrease rate. For example, Figs. 6(c) and 6(d) show different structures of the electron and ion temperature in the case of the random $η ( m , n )$ configuration. The electron temperature structure is strongly correlated with the magnetic mirror ratio structure. The electron temperature is the lowest in the magnetic hill, and the temperature in the surrounding area is lowered due to the E × B mixing effects. However, the electron temperature at the outboard midplane (⁠ $θ ≈ 0$⁠) still remains high since the E × B mixing by radial eddies is difficult to transfer electrons to magnetic hills far from the outboard midplane. Unlike the electron temperature, the ion temperature at the outboard midplane decreased because large excursions of ions by magnetic drift can cause ions to enter magnetic hills far from the midplane, leading to collisionless detrapping. As a result, the ion temperature drop is observed over a wider range, not limited to the magnetic hill regions. This wider range of temperature drop results in a faster temperature decrease on average, especially in the inner plasma region where there is no magnetic hill.

In this study, we proposed a novel understanding of the 3D magnetic topology and associated plasma dynamics in the presence of open stochastic magnetic field lines connected to the wall boundary. The complex 3D structure of stochastic field lines can be characterized in terms of two key concepts, the connection length L_c and the effective magnetic mirror ratio $M eff$⁠. The L_c determines the characteristic length and time scales of the 1D parallel plasma dynamics along each field line. The 3D structure of L_c, in which the longer and shorter L_c regions are intertwined, strongly affects the structures of the plasma density and the ambipolar potential. The new concept of $M eff$ describes the topology of the 3D magnetic potential, consisting of magnetic wells (⁠ $M eff > 1$⁠) and hills (⁠ $M eff < 1$⁠), which has a significant impact on electron and ion thermal energy confinements and their temperature structures. In particular, the existence of magnetic hills is a very important key to understanding the electron and ion temperature drop in the stochastic layer.

The plasma dynamics considering the consistent coupling between electrons and ions can be understood by a combination of 1D parallel dynamics along each field line and 3D effects associated with the magnetic topology. The positive ambipolar potential develops in the stochastic layer to balance the electron and ion losses along the field line for maintaining the quasi-neutrality of the plasma. As the plasma density falls faster in the shorter L_c region than in the longer L_c region, the positive ambipolar potential becomes lower in the shorter L_c region than in the longer L_c region. As a result, the overall structures of the plasma density and the ambipolar potential are strongly correlated with the L_c structure. Meanwhile, sharp gradients of the ambipolar potential across different L_c regions induce strong E × B transport and mixing effects. The radial eddies are created in the process of relaxing the sharp potential gradients between different L_c regions.

To understand electron and ion kinetic physics in open stochastic magnetic fields, the positive ambipolar potential and the 3D magnetic potential must be considered together. Due to the opposite charge of electrons and ions, the ambipolar potential plays the opposite role for electron and ion confinement. The positive ambipolar potential acts as a barrier to trap electrons before escaping toward the wall, whereas it helps ions escape along the field line by accelerating the ions. On the other hand, the magnetic potential plays the same role for electrons and ions as it does not depend on charge but on the magnetic moment μ. While magnetic wells (⁠ $M eff > 1$⁠) trap charged particles of high-pitch angles before escaping toward the wall, the magnetic hills (⁠ $M eff < 1$⁠) help charged particles to escape by accelerating them toward the wall. Combining both electric and magnetic potential effects, the passing-trapped conditions for electrons and ions in magnetic wells and hills can be derived, respectively. The passing-trapped conditions provide physical insight into electron and ion kinetic physics.

In the case of electrons, low-μ electrons are trapped in both magnetic wells and hills. High-μ electrons are trapped in magnetic wells but they can be passing in magnetic hills. The selective loss of high-μ electrons with high thermal energy efficiently cools the electron temperature in magnetic hills. As a result, the electron temperature in magnetic hill regions becomes much lower than the magnetic well regions where the majority of electrons are trapped. The high-μ electron trapped in the magnetic well region can enter the nearby magnetic hill region by the E × B mixing and have a chance to get detrapped and escape the wall along the field line. The collisionless detrapping of high-μ electrons occurs steadily by the E × B mixing between magnetic wells and hills, so that the electron temperature near the magnetic hill can be lowered.

In the case of ions, all ions in magnetic hills are passing ions. In magnetic wells, low-μ ions can be passing ions with the help of the positive ambipolar potential, but high-μ ions are trapped by dominant magnetic mirror force. Unlike electrons, ions have large radial and poloidal excursions due to $∇ B$ and curvature drift motion. The ions far from the magnetic hills can be transferred to magnetic hills by the large excursion and get detrapped. This collisionless detrapping of ions results in ion temperature drop in a wider range, and the temperature decrease rate is even faster than that of electrons. Meanwhile, the E × B mixing between magnetic wells and hills causes both collisionless detrapping and trapping of ions due to the slow parallel motion. Since the collisionless trapping and detrapping mechanisms compete with each other, the E × B mixing may accelerate or deaccelerate the ion temperature decrease rate depending on the situation.

It is important to note that the efficiency of the collisionless detrapping mechanism is highly dependent on the existence and structure of the magnetic hills in the stochastic layer. Global gyrokinetic GTS simulations were performed for two different magnetic configurations, the uniform and random $η ( m , n )$ cases, to compare the effects of the connection length L_c and the magnetic mirror ratio $M eff$ on the electron and ion temperature drop. It was found that, surprisingly, the uniform $η ( m , n )$ case that has longer connection length shows a twice faster electron temperature decrease rate than the random $η ( m , n )$ case that has shorter connection length. This is because the uniform $η ( m , n )$ configuration has more magnetic hills deeply penetrating into the plasma interior, leading to more collisionless detrappings by E × B mixing and faster electron temperature drop. This indicates that the electron temperature and heat transport are not simply determined by the connection length, but are strongly influenced by the existence of magnetic hills and the associated collisionless detrapping mechanism. Therefore, the collisionless detrapping mechanism could be a major cause of the rapid electron temperature drop in the presence of open stochastic magnetic field lines connected to the wall boundary. In particular, this mechanism could be more efficient in the thermal quench situation because very strong magnetic perturbations are more likely to create a wide stochastic layer with magnetic hills deeply penetrating into the plasma interior. The ion temperature shows a similar level of temperature decrease rate for both magnetic configurations. This is because the large excursion of ions by the magnetic drift allows ions can more easily access magnetic hills far away. Therefore, the dependence of the ion temperature decrease on the magnetic hill structure is less sensitive than the electron.

There are other effects that could play important roles in the TQ process, such as the Coulomb collision and particle recycling from the wall, which were not considered in this study. Since these effects are likely small in the early stages of TQ, where the plasma temperature is still sufficiently high, the collisionless mechanism related to the 3D magnetic topology may still dominate. However, as electron temperature decreases and a significant amount of lost particles interacts with the wall, the Coulomb collision and particle recycling must be considered for a more appropriate understanding. In this study, we assumed prescribed static magnetic perturbations for simplicity. In order to understand the entire thermal quench process, it is necessary to consider the temporal evolution of magnetic perturbations according to time-varying plasma density and current, which is an interesting but challenging work. In future work, we will further study the thermal quench physics on more realistic magnetic configurations, taking into account collision effects and time-dependent magnetic perturbations.

This work was supported by U.S. DOE Contract No. DE-AC02–09CH11466 and SciDAC Tokamak Disruption Simulation project. Numerical simulations were conducted on Traverse at PPPL.

The authors have no conflicts to disclose.

Min-Gu Yoo: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing - original draft (equal). W. X Wang: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Supervision (equal); Writing - review and editing (equal). E. Startsev: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing - review and editing (equal). C. H. Ma: Methodology (equal); Software (equal). S. Ethier: Software (equal); Writing - review and editing (equal). J. Chen: Software (equal). X. Z. Tang: Funding acquisition (equal); Project administration (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.