Multiscale particle simulation of the temporal evolution of heat flux into poloidal gaps of castellated divertor with edge-localized modes

Edge-localized modes (ELMs) produce an intense transient heat flux on castellated divertor mono-blocks (MBs). The parallel transport of electrons and ions from the ELM burst varies due to their different velocities, which significantly influences the power load on the MBs. This study shows that two heat load phases are caused by an ELM burst. In the first phase, the horizontal surface withstands the intense heat flux because the faster electrons generate a high sheath potential drop. The leading edge of the MBs suffers a great power load from energetic ions in the second phase. This study applies a multiscale hybrid numerical approach to investigate the temporal evolution of heat flux on the poloidal gap from the ELM source in the scrape-off layer to the castellated divertor surface. Our results indicate that the power load on both the MB surface and the leading edge should be considered in the future design of castellated divertors.


Introduction
Plasma-facing units (PFUs) made of tungsten mono-blocks (MBs) withstand the most intense particle and heat fluxes [1,2]. To minimize the thermomechanical stress, inter-PFU gaps can be introduced between neighboring MBs [3]. As a result, the particle and heat fluxes on the castellated divertors are nonuniform. The flux into the region between MBs is extremely high, leading to serious degradation of the PFUs, including melting [4,5], sputtering [6], macro-crack formation [7], and even deep cracking [8]. Experiments conducted in the COMPASS tokamak have investigated the nonuniform flux distribution on poloidal and toroidal castellated divertor gaps (CDGs) [9] and misaligned edges [10]. Particle-in-cell (PIC) simulations of plasma-material interactions have focused on the poloidal gaps between MBs in the high-heat-flux areas of the outer vertical target in the ITER tokamak [11].
Edge-localized modes (ELMs) can produce transient high heat loads on PFUs, resulting in serious heat flux in the MB gaps. Thus, studies have been performed to examine the transient power load on the leading edge and the misaligned edges. The efficiency of an optical approximation model for the power deposition on the misaligned castellated tungsten blocks has been investigated under the steady state and ELM heat flux in the Magnum-PSI and Pilot-PSI linear devices [12]. The asymmetric plasma deposition in the castellated gaps during ELM can be modeled by a self-consistent PIC simulation [13]. The heat load on the MBs produced by a simple ion orbit model showed that the magnetic shadowing provided by surface beveling in the toroidal direction protects the poloidal leading edges at the inter-PFU poloidal gaps from overheating, while partially mitigating the risk of ELM-induced melting [14]. The effect of toroidal and poloidal beveling of the MBs on the ITER divertor with ELM has also been studied [15]. Serious melting and droplet ejection occur on the ITERlike castellated divertors made of pure tungsten and lanthanum tungsten under ELM and disruption heat loads in a high-flux plasma gun [16]. During ELMs, the risk of melting of MB edge and full surface is influenced by the shapes/structures of MBs [17]. The effect of misaligned edges and the magnetic field orientation on plasma deposition into the gaps during ELMs has been studied for projected ITER parameters [18]. The transient heat load under an energetic ELM on the misaligned W lamella was studied by a series of numerical simulations and experiments [19], and the effect of impurities on the plasma and power loads onto and between divertor tiles was studied using the PIC method [20]. The kinetic model has been applied in the COMPASS tokamak scrape-off layer (SOL), allowing the parallel heat flux on the outer and inner target divertors to be compared [21].
Although the temporal evolution of the heat load during an ELM burst has been considered in previous studies, the density and temperature evolution of electrons and ions were not differentiated. In other words, the synchronous evolution of electrons and ions has been implicitly assumed, but this is not sufficient [22]. Dai et al found that the evolution of the impact energy of charged particles is significantly influenced by particle transport in the SOL due to the different velocities of electrons and ions. Thus, two heat load phases are caused by an ELM burst. The horizontal surface of the MBs suffers a serious heat load in the first phase because of the high electron temperature and high sheath potential drop. In the second phase, the leading edge withstands intense power due to the high ion temperature, as we will show in this paper.
In this work, we model the plasma transport in both the SOL and the CDG to establish a multiscale and self-consistent simulation for the heat flux distribution on the EAST CDGs. The remainder of this paper is organized as follows. The hybrid simulation method is introduced in section 2, and the simulation results are presented in section 3. Plasma transport in the SOL region is discussed as the basis of the temporal evolution of electrostatic potential on the CDG. The temporal evolution of the electrostatic potential significantly influences the heat flux distribution on the CDG. Finally, the conclusions to this study are summarized in section 4.

Method
Two-dimensional in space and three-dimensional in velocity (2D3V) PIC simulations are appropriate for studying the temporal evolution of the power flux into CDGs. However, these simulations are complicated by the great differences between the required accuracy and the spatial and temporal scales. In terms of the spatial scale, the width of the castellated gap is 0.5 mm. This means that the grid length must be significantly smaller than 0.5 mm to obtain sufficient accuracy to reflect the effect of the CDGs. The distance from the source of the ELM (which is near or at the midplane) to the divertor target is over 50 m. Regarding the temporal scale, the time step in PIC simulations for a tokamak environment is of the order of 10 −12 s. The duration of a typical ELM burst is about 0.3 ms. Such large differences between the simulation accuracy (0.5 mm and 1 ps) and the simulation dimensions (50 m and 0.3 ms) pose an unconquerable difficulty under present computational capabilities. As a result, we divide the simulation process into three steps.
The first step is to calculate the parallel transport of plasma from the ELM source to the PFU surface in the SOL region. Plasma transport is important for determining the power load.
The key factor that influences the transport is the so-called connection length, which is the distance from the ELM source to the target. A longer connection length creates a more significant difference between the arrival times of the electrons and ions. Thus, the main physical processes we focus on are (1) the particle dynamics and (2) the electrostatic potential evolution. We do not consider the full set of physics. The thickness of the plasma sheath (specifically, the magnetic pre-sheath) on the CDGs is of the order of ρ i − 1 mm, where ρ i is the Larmor radius of an ion. As the sheath thickness is much smaller than the dimension of parallel transport, the effect of CDGs on the parallel transport can be ignored. We also neglect the perpendicular transport of particles, because this scale is much greater than the gap width. As a result, a 1D PIC simulation is sufficient to track particles in the SOL and provide the necessary information. Coulomb collisions, as well as atom and molecular processes, are included in our simulations because the particle energy is influenced by these processes. The second step is to obtain the electrostatic sheath potential profile on the CDGs, for which a 2D3V PIC simulation is necessary. The plasma density and temperature at the sheath entrance are provided by the plasma parallel transport step as input parameters.
Based on the electrostatic potential calculated in the 2D PIC step, the trajectory of particles is tracked in the final step through a 2D3V particle-tracing (PT) method. The distribution of power flux on the CDG can be obtained from the trajectory tracking. The hybrid simulation method is described in the following subsections.

Modeling of parallel transport
The modeling of plasma parallel transport in the SOL region is performed by the PICS1+ code, which has been developed based on the PICS1 code [23]. PICS1 is a 1D collision-less PIC code for modeling plasma sheaths. We extend the utility of this code to the SOL region by updating and adding various algorithms. The physical processes of interest in the parallel transport and numerical methods are listed in table 1. Further details are given in the following sections.
2.1.1. Geometry of simulation box and magnetic field. The simulations use the magnetic field from the EAST experiment # 69634 discharge (upper-single null). Particles move along a selected magnetic field line whose normalized poloidal flux is ψ n = 1.0043. The separatrix surface and the selected field line pass through the coordinates (R, Z) = (2.324, 0) and (R, Z) = (2.325, 0) m, respectively. The purpose of the selection is to reflect the heat flux evolution near the strike point. The reason we did not choose the separatrix surface directly is that the length of magnetic field line will be numerically infinite because the poloidal magnetic field strength at the X point is zero.
The three-dimensional (3D) view of the magnetic field lines are shown in figure 1(a) in the R, Z, and Θ coordinates. Θ is the toroidal angle. The color scale is the magnetic field strength |B|. We track the magnetic field line in the R-Z-Θ coordinate space in both directions until the line reaches the outer and inner targets. The magnetic field data are from the EAST database. The field line is traced by the MATLAB-defined streamline calculation function. The total length of the magnetic field line (i.e. the connection length from outer target to inner target) is L c = 110 m. The magnetic field turns about four rounds near the upper X point. This is because the R and Z components of the magnetic field, B R and B Z , are small near the X point, and so the pitch angle θ = arctan the Θ component of the magnetic field). The small pitch angle is also the reason for the long connection length. The projection of the magnetic field line on the poloidal cross-section at Θ = 0 • is shown in figure 1(b). The crosspoint between the selected field line and the outer target is the location of the CDG we study, as denoted in this figure.
Figure 1(c) shows the magnetic field strength |B| as a function of the curve length s, which is defined as the distance to the outer target along the magnetic field line. For convenience, coordinates are described by the curve length s rather than the 2D coordinates (R, Z).
Energetic particles are generated during ELM burst from the middle of the selected magnetic field line (as denoted in figure 1(c)) and transport from the source region to outer and inner targets. The lower and higher boundaries of source region are s = L − and s = L + , respectively. The lower source boundary is L − ≡ 40 m, which locates at the midplane of the low field side. The selection of L − is based on previous studies [24]. The higher source boundary is L + ≡ L c − L − = 70 m. We have to concede that the value of L + is arbitrary because ELM particles are predominantly ejected from the low field side [24]. But the issue we focus on is the heat flux on the outer target, rather than the inner target. So, the high field side boundary will not significantly influence our results. The selected magnetic field line is discretized by an even grid with a fixed grid size of ∆s, as shown in figure 2. The outer and inner divertor targets are at s = 0 and s = L c , respectively.
Before the ELM, the plasma density and temperature at the source are n e0 = 10 19 m −3 and T 0 = 50 eV, respectively.  The subscript 0 refers to the source. The Debye length and plasma frequency in this condition are λ D = √ ε 0 T 0 /e 2 n e0 = 16.7 µm and ω pe = √ e 2 n e0 /ε 0 m e = 1.84 × 10 11 s −1 , respectively. For a classical PIC simulation, the grid size δs and time step δt are strictly constrained by δs < λ D and δt < 1/ω pe to keep the calculation stable. As the connection length L c ≫ λ D and the total simulation time τ total is at the order of milliseconds, we have that τ total ≫ 1/ω pe , making it almost impossible to use a standard PIC method for our simulations.
To balance the calculation time against the simulation accuracy, a speed-up technique is implemented [25] by artificially enhancing the value of the vacuum permittivity ε 0 to a higher value ε * 0 . A shortening factor is defined as k = √ ε * 0 /ε 0 . By doing this, the new Debye length λ * D = √ ε * 0 T 0 /e 2 n e = kλ D and the new plasma frequency ω * pe = √ e 2 n e /ε * 0 m e = ω pe /k.
Thus, we can choose a much higher-resolution grid size ∆s and time interval ∆t. In the present study, we choose k = 10 3 , δs = 10 µm, and δt = 5 ps. Thus, the actual grid size and time interval are ∆s = 10 mm and ∆t = 5 ns, respectively. The enhancement of ε 0 changes Debye length λ D . Thus, this treatment is not available on the studies closely related to Debye sheath, such as simulations on plasma-wall interactions. On the other hand, plasma transport in the SOL can be effectively simulated by the treatment due to the quasi-neutral condition [25]. The correct ε 0 is applied in the 2D PIC simulation for the heat load, as we will illustrate in the following part.
The maximum temperature in the source region is T 0,m = 1 keV. Similar to the ELM assumptions in the ITER simulations [17], the ELM pulse is assumed to have a triangular waveform. The temperature of plasma source T 0 (t) rises from 50 eV to 1 keV linearly from t = 0 to t = τ r and then drops linearly from 1 keV to 50 eV in τ d . After the ELM burst, T 0 returns to 50 eV. The plasma temperature T 0 (t) is synchronous in the whole source region (L − ⩽ s ⩽ L + ).
The averaged initial parallel velocity of ions during the ELM burst is u i,∥ = √ T 0,m /2m i ≈ 150 km s −1 . The factor of 1 2 gives the average of the ELM plasma temperature. The rise time of the ELM pulse is given by the flight time of ions from the mid-point L c /2 to the low field boundary L − , τ r ≈ (L c /2 − L − )/u i,∥ = 0.1 ms. We further assume that the decay time τ d = 2τ r . The factor of 2 is supported by experiments [26]. Therefore, the duration of the ELM burst considered in the present study is τ ELM = τ r + τ d = 0.3 ms. The ion transit time from the source to the outer target is τ i,∥ = L − /u i,∥ ≈ 0.3 ms. The total duration of the simulations is τ total = 4 ms > τ ELM + τ i,∥ . Therefore, the total numbers of grid cells and calculation steps are N = L c /∆s = 11 000 and N t = τ total /∆t = 800 000. Before the ELM burst, a pre-calculation is performed to obtain the initial particle distribution.
The particles are injected in the source region s ∈ (L − , L + ) with uniformly distributed initial positions. The initial velocities of the electrons and ions follow the Maxwell distribution. In the source region, the electrons and ions have average velocities of zero. We assume that T e0 (t) = T i0 (t) = T 0 (t) in the source for the whole simulation.
In each PIC step, the electrostatic field ⃗ E = −∇Φ is updated by solving Poisson's equation consistently coupled with charged particles by the successive overrelaxation method. By applying the speed-up technique, the discretized Poisson's equation is written as instead of the original form Φ i = 1 2 . Here, i is the cell index, ρ = ∑ α Z α n α is the net charge density, and n α is the density of species α.
The plasma parameters (particle number, velocity, temperature, and kinetic energy) are calculated by weighting the particles linearly on the nearby cells according to the particle positions as shown in figure 2. The weighting function is here i is the cell index and s i = (i + 0.5) ∆s is the position of the ith cell. s p is the position of the specific particle. The electric field acting on the particle is determined by here E(s i ) is the electrostatic field strength in the ith cell. The same interpretation method is applied for the magnetic field strength.
The magnetic field is assumed to be constant. In other words, the ELM-perturbed magnetic field is ignored and the pre-ELM magnetic field is applied in the whole simulation. The reason is that the perturbed magnetic field is significantly smaller than the background magnetic field [27].

Particle pushing.
To obtain high-accuracy particle positions and velocities, the 4th-order Runge-Kutta method is applied to the particle pushing. The position and velocity are advanced four times in a Runge-Kutta step. The electrostatic potential is updated in every sub-step to obtain the electric force acting on the charged particles. The guiding-center approximation is used in each sub-step. The position s and parallel velocity u of particles of species α are described by the guiding-center approximation. q α and m α are the charge and mass for species α, respectively. The equations of movement are as follows [25,28]: The magnetic momentum µ = m α v 2 ⊥ /2|B| remains constant during collision-less motion. For 2D or 3D simulations, there are additional drift terms v E×B and v ∇|B| for E × B and curvature ∇|B| drifts. For 1D simulations, both drifts can be ignored because they are perpendicular to the magnetic field. The polarization drift is also neglected.

Coulomb collision.
A four-step random-pairing binary collision algorithm is applied to describe the Coulomb interactions between charged particles (electron-electron e ↔ e, electron-ion e ↔ i, and ion-ion i ↔ i). The steps are (1) particle grouping, (2) random sorting, (3) particle pairing, and (4) velocity changing. The particles are grouped according to the coordinates in each cell. For the particles in a cell, the index of the particles is randomly re-ordered. The random sorting step ensures that the particles are selected randomly. To guarantee that each particle collides at least once in each collision step, particle pairing is executed as follows. For the self-collision cases (e ↔ e, i ↔ i), if the number of particles in a cell N is even, particles with indexes 2j ↔ 2j + 1 collide, where j runs from 0 to N/2 − 1. When N is odd, the first three particles collide with each other, namely, 0 ↔ 1, 1 ↔ 2, 0 ↔ 2. The particles are then paired as 2j + 1 ↔ 2j + 2, where j runs from 1 to (N − 1)/2 − 1. In the electron-ion collision, the particles are paired directly in order when the particle numbers are equal, N e = N i . When the particle numbers are not equal, N e > N i , for example, some ions suffer multiple collisions with different electrons. The particle pair j e ↔ j i has the relation j e = kN i + j i , where k = 0, 1, · · · is an integer.
The Takizuka-Abe (TA) binary collision model [29] is applied in the velocity changing step. In the time interval ∆t, the electrons and ions in a cell are randomly selected for binary collision.
The velocity of the guiding-center approximation is set The guiding-center velocity components should be transformed to the Cartesian system v = (v x , v y , v z ).

In the magnetic field
where ξ is a random number for the gyration phase that is uniformly distributed between zero and 2π and b xy = is a uniformly distributed random number. The scattering angle θ, which represents the angle between ⃗ u(t + ∆t) and ⃗ u(t), is obtained in two ways according to the variable D = e 4 nL ln Λ 8π ϵ 2 0 m 2 12 u 3 ∆t, here e, n L , ln Λ and m 12 are electron charge, the lower density between n 1 and n 2 , Coulomb logarithm and reduced mass m 1 m 2 /(m 1 + m 2 ), respectively. When D < 1, θ is calculated from the relation θ = 2 arctan δ. δ is a Gaussian distributed random number with the mean ⟨δ⟩ = 0 and the mean square (Here, ⟨·⟩ denotes the average value.) When D ⩾ 1, θ is randomly selected between zero and π uniformly.
The updated velocities of particles 1 and 2 are The TA model obeys the energy and momentum conservation laws. Obviously, This means that the TA model satisfies the momentum conservation law. It is easy to know that (⃗ u + ∆⃗ u) 2 = ⃗ u 2 , which implies the conservation of total kinetic energy.
Moulton et al [30] indicates that the velocity distribution is affected by Coulomb collisions. Although the collision frequency may not be high under the condition n e ∼ 10 19 m −3 and T e ∼ 50 eV, the main consideration under Coulomb collisions is that the collisions occur over the whole SOL transport process, not just near the target. The connection length from source to target is over 40 m and the transit times of electrons and ions are ∼ 4 µs and 0.3 ms. The effect of Coulomb collisions over such a long time may be significant. As the effect of Coulomb collisions is not the main concern of this paper, we will not discuss this topic further, although we include the Coulomb collision model in the simulations.
2.1.5. Atom and molecular processes. The A&M processes (ionization, charge exchange, and excitation) are treated by a simple direct Monte Carlo method. All particles, including electrons, ions, and molecules, are treated in the same way. The collision pairs are selected from the grouping in the TA boundary collision. As we are concerned with the heat flux on CDGs, a divertor attachment state is assumed. Thus, a relatively low neutral gas pressure is applied, and we set n H2 = 2 × 10 19 m −3 with a gas temperature of 300 K. Fifteen reactions are considered in our simulations, including eleven electron-impacting reactions and four proton-impacting reactions, as listed in table 2. The right column (Index) refers to the reaction index in the HYDHEL database [31].
The A&M interaction processes are simulated for an electron/ion and neutral hydrogen molecule pair. For an electron or ion, the total probability of an A&M interaction is P = 1 − e −∆t ∑ j nH 2 σju , where j is the reaction index, u is the relative velocity, and n H2 is the density of H 2 . The collision cross- ) is a function of the impact energy E = (mu 2 )/2. The coefficients a i are taken from the HYDHEL database. When κ < P, one A&M interaction happens. Here, κ is a random number between zero and unity. An A&M interaction is then selected and randomly weighted by the counterpart cross-section σ. If no new particles are generated (excitation, for example), the impacting particle loses ∆E energy, as listed in table 2. The excited molecule is quickly deexcited and produces a photon. As the process of de-excitation is very fast, the influence of excited A&M processes on ionization can be ignored. In other words, the effect of excitation in our simulation is energy loss. If new particles are generated, besides the energy ∆E, extra energy is lost from the impacting particle and transported to the new particles. The value of the extra energy is randomly chosen from 0 to E − ∆E. The velocity directions of the new particles are randomly selected. The positions of the new particles are the same as the colliding particles. The recombination mainly occurs in a dense and cold plasma, i.e. detachment state. Our study assumes that the divertor plasma is in the attachment state with low neutral gas density and hot plasma, so the recombination process can be ignored.

Modeling of the plasma sheath on CDGs
The geometry of the poloidal CDG is shown in figure 3(a) (in 3D) and (b) (in 2D). The poloidal gap is parallel with the poloidal direction. The direction of heat flux is approximately perpendicular to the poloidal CDG because the poloidal magnetic field The CDG that we are studying locates at (R CDG , Z CDG ) = (1.747, 1.041), which is the cross-point between the selected magnetic field line and the EAST boundary (denoted in figure 1(b)). The strike point is denoted by the black dashed line, which is quite near to the line of heat load position of heat flux, as shown in figure 3(a). The reason we did not choose the strike point as the load position is illustrated above (section 2.1.1). The distance between the strike point and flux position is 2 mm, which is significantly smaller than the heat width on the target (∼30 mm [32]), so our results can reflect the heat flux evolution at the strike point.
In this step, a 2D3V PIC simulation is performed to study the electrostatic field on the PFC surface and CDGs. The PICS2 code is applied. This code is designed to perform simulations of magnetized plasma in contact with solid objects, and has been successfully used for the study of plasma deposition in the vicinity of castellated PFC [33]. Although 3D geometry of castellated gaps has been employed in tokamaks [9], we applied 2D PIC simulation in present study for EAST CDGs. There are two reasons. First, the geometry of CDGs between EAST MBs is simple. In other words, there is no complex structure that needs 3D modeling (such as misaligned gaps). Second, electrons can influence the heat load distribution by moving through the crossings between poloidal and toroidal gaps [34], which causes the necessity of 3D modeling near the crossings. However, the length of CDGs (i.e. size of MBs) is several centimeters and the width of CDGs is 0.5 mm.
The electrons, which are constrained by magnetic field, are hard to influence the entire heat flux distribution of CDGs. Thus, 2D modeling is sufficiently for EAST CDGs. The Cartesian coordinate system is used to describe the geometry for the 2D simulations, as shown in figure 3(b). The x and z directions represent the coordinates parallel and perpendicular to the PFC surface, respectively. Because the difference between the Debye length and the simulation spatial scale is quite large, the 2D PIC calculations are very timeconsuming. As a result, we reduce the size of the PIC simulation box and only retain the parts necessary for calculating the electrostatic potential. The horizontal range is from x = 0.0-1.0 mm. The gap is from x = 0.0-0.5 mm, i.e. the width of the gap is x gap = 0.5 mm. The PFC surface runs from x = 0.5-1.0 mm, a width of x surf = 0.5 mm. The total width is x PIC = x gap + x surf = 1.0 mm. The periodical boundary condition is applied in the x direction. The depth of the CDGs is z gap = 1 mm and the height above the PFC surface is z h = 4.0 mm. The geometry of the 2D PIC simulations is shown in the inner box (denoted by the dotted line) in figure 3(b).
As the simulation dimensions are very small, the magnetic field is assumed to be constant. The magnetic field vector is ⃗ B = (−B 0 cos θ, 0, −B 0 sin θ), with B 0 = 2.524 T. The glancing angle between the magnetic field and the PFC surface is θ = 5 • . The direction of incident flux is the same as ⃗ B. The plasma density and temperature along the x direction are constant at the sheath entrance z = z h ≡ 4 mm (denoted by the blue line in figure 3(b)) where electrons and ions are injected.
The collisionless motion of a particle of the species α (mass m α and charge q α ) is described by The 4th-order Runge-Kutta algorithm is applied to increase the accuracy of particle position and velocity. Due to the electrostatic field, the particle trajectory cannot be written in an analytical form. Thus, the differential form is applied to describe the particle trajectory T are the coordinate and velocity, respectively. The subscript n denotes the nth step. The trajectory is obtained by iteration of equation (7). The sheath entrance is determined by finding the position s s where n e ̸ = n i . The plasma density at the sheath entrance n s = n(s s ) comes directly from the 1D PIC simulations. We assume that electrons/ions are composed of hot particles from the ELM and cold background particles. Thus, the velocity distribution at the sheath entrance is obtained by fitting the double-shifted Gaussian function where represent the velocity distribution of cold and hot particles, respectively. The coefficient A ∈ (0, 1) denotes the weight of the two components. v c and v h are the shifted velocities of the cold and hot distributions, respectively. The drift velocity in the perpendicular direction is zero, i.e. v c,⊥ = v h,⊥ = 0. A, v h , T h , v c,∥ , and v h,∥ are obtained by the least-square fitting of the counterpart velocity distribution near the sheath entrance from the PIC simulations. The fitted temperature satisfies the relation T s ≈ AT c + (1 − A)T h , where T s is the electron/ion temperature at the sheath entrance from the 1D PIC step. The bias is introduced because the particles in several cells near the sheath entrance provide sufficient velocity data for the velocity fitting, while T s is the temperature of the sheath entrance cell.
An implicit assumption is that the parallel and perpendicular temperature and coefficient are the same, namely, T ∥ = T ⊥ and A ∥ = A ⊥ . In the collision-less case, this does not hold because T ∥ changes with the parallel velocity, which is influenced by the electric field, while T ⊥ is constant. The reason for assuming that T ∥ = T ⊥ is that Coulomb collisions between like-particles produce an exchange of energy in the parallel and perpendicular directions and make T ∥ ≈ T ⊥ . We concede that this treatment is not sufficiently accurate, but the fitting result shows that such an approximation is acceptable.
The modeling of plasma-surface interactions has been simplified in the present research by making two assumptions. The first is that all particles are deposited on the PFC surface and the reflection rate is zero. The second assumption is that the secondary and thermal electron emissions can be neglected. There are two reasons for these simplifications. First, it is very hard to estimate the reflection and electron secondary and thermal emission rates. Second, the present research focuses on the effect of the ELM, and other factors are beyond the scope of this paper. When the temperature of MB is close to the melting point of tungsten, the thermionic emission can collapse the sheath potential and our result is not available in such case.

Calculation of heat flux in CDGs
We extend the 2D PIC simulation domain in the PT step because the width in the second PIC step (1.0 mm) is not sufficient to reflect the flux distribution on the MBs. The range of the PT step is from x = −2.0 mm to x = 2.5 mm. The periodical boundary condition is applied in the x direction. The extended electrostatic potential Φ ext in the PT step is obtained from the potential Φ in the PIC step: As Φ ext is fixed, the restriction of dt < 0.1/ω ce disappears because the accurate electron trajectory is not needed. The restriction of the time step for electrons is widened to dt e < dx/v te = dz/v te , where dx and dz are the grid sizes in the x and z directions, which are the same as the grid sizes in the 2D PIC step. Another advantage of the PT method is that different time steps can be used for electrons and ions. Because the Larmor radius of ions ρ i = (m i v ti )/qB ∼ 0.4 mm ≈ x gap , the full-orbit of ions must be obtained. The movement of ions is described by equations (6) and (7). The restriction from the gyration motion is dt i < 0.1/ω ci = 0.8 ms. Moreover, dt i < dx/v ti = 300 ps should be satisfied. Thus, dt i = 300 ps. The ions are tracked by the 4th-order Runge-Kutta algorithm. ⃗ r(t) and ⃗ r(t) + ∆⃗ r are the positions before and after particle pushing, respectively. ⃗ r(t) + ξ ∆⃗ r is the cross-point between the particle trajectory and mono-block surface. L is the curve length from A to the cross-point.
The guiding-center algorithm of electrons can be applied: The vectors ⃗ e x and ⃗ e z are the normal vectors in the x and z directions, respectively. Electrons and ions are injected at the upper boundary z = z h . The initial velocity distribution is the same as that in section 2.2. A particle is deleted when it hits the surface of the MBs or CDGs, or escapes from the simulation box at the upper boundary.
The positions before and after particle-pushing of a particle are ⃗ r(t) and ⃗ r(t + ∆t), respectively. The position change is ∆⃗ r =⃗ r(t + ∆t) −⃗ r. If ⃗ r(t + ∆t) is in the MB (as shown in figure 4), the particle is treated as landing on the surface and its kinetic energy is added to the power flux. By calculating the cross point between the particle trajectory and target boundary, the landing point (which can be expressed as ⃗ r + ξ ∆⃗ r) and the interpretation factor ξ can be obtained. The velocity of the landing particle hitting the material is Because the time interval of particle pushing is small enough and the trajectory between ⃗ r(t) and ⃗ r(t + ∆t) is approximately straight, the linear interpretation can provide the accurate value of velocity.
The points A-E are defined to describe heat flux distribution in figure 3  n is obtained by weighted summation of the kinetic energy E j = 1 2 m α |⃗ v h | 2 of the jth particle whose landing point is L j : The recycled flux (i.e. reflected particles) and the surface recombination process are not considered. The reason is that the major kinetic energy deposits on the surface are independent of whether a particle is reflected or recombined.  5(c)), T 0 reaches 1 keV. ⟨T e ⟩ s is higher than that in figure 5(b) and ⟨T e ⟩ s is close to 0.5 keV. By contrast, the ion temperature T i has not significantly increased. T i near the source is 0.1 keV, while T i at the target is still around 20 eV. In this period, the electron temperature is significantly greater than the ion temperature. The period for which T e > T i is from 0 ms to about 0.3 ms, which is named the electron dominating phase (EDP).

Temporal evolution of plasma parameters in SOL
T i is significantly greater than T e when t = 0.25 ms, as shown in figure 5(d). There is a significant peak of T i at s peak ≈ 25 m. The peak is generated by the peak of T 0 , because the location is consistent with the distance of the ion transport from the source region. The distance moved by the energetic ions from τ r to t is v i,∥ (t − τ r ) ≈ L − − s peak . Moreover, T e becomes flat and ⟨T e ⟩ drops to around 0.14 keV.
At t = 0.75 ms (figure 5(e)), the distance moved by the ELM ions is v i,∥ (t − τ r ) > L − . Thus, there is no significant peak in the profile and T i becomes flat. ⟨T i ⟩ s is around 0.1 keV, while ⟨T e ⟩ s is about 50 eV. In this period, the ion temperature is greater than the electron temperature. The reason is that most of the energetic electrons arrive at the target due to their higher velocity, whereas the ELM ions are still in the SOL. This period is called the ion dominating phase (IDP). The IDP lasts from 0.5 ms to about 2.5 ms.
After the IDP, T e ≈ T i (figure 5(f ), t = 2.5 ms), which is defined as the post-ELM phase (PEP).

Temporal evolution of plasma parameters at the sheath entrance
We are especially interested in the plasma density and temperature at the sheath entrance on the CDG because they affect the electrostatic potential profile. The simulation results between t 1 = 0.0 ms and t 9 = 3.0 ms are shown in figure 6. The duration of the simulation time is 4.0 ms, but the plasma parameters from 3.0 to 4.0 ms remain the same and are not shown. The ELM burst begins at t 1 . Before t 1 , a 400 000-step PIC simulation is performed to obtain the steady state of the pre-ELM plasma, which is not shown.
The source temperature is shown in figure 6(a). The red solid line denotes the T 0 evolution. The shifted temperature profiles are also shown in this figure; these are helpful for explaining the tendency of the plasma density and temperature. The temperature profile shifted by the electron transit time τ e,∥ ≈ 4 µs (green dotted line), T e,shift = T 0 (t − τ e,∥ ), almost overlaps the source temperature profile because τ e,∥ ≪ τ ELM . By contrast, the temperature profile shifted by the ion transit time τ i,∥ ≈ 0.3 ms (blue broken line), T i,shift = T 0 (t − τ i,∥ ), completely separates from the source profile because τ i,∥ ∼ τ ELM = 0.3 ms. T i,shift reaches its peak value at t 5 = τ r + τ i,∥ .
The evolution of the normalized plasma density n s /n 0 at the sheath entrance is shown in figure 6(b). The plasma density of 0.4 remains almost constant between t 1 and 0.1 ms, when the sheath has not been significantly influenced by the ELM. After that, n s /n 0 begins to increase at about 0.1 ms and reaches its first peak at about t n1 ≡ t 2 = 0.16 ms. The value of the first peak is 0.5. The plasma density then drops quickly to the initial level at t 3 = 0.2 ms. The duration of the first peak is very short (∼t 3 − t 2 = 0.04 ms). At t 4 = 0.3 ms, n s begins to increase again and takes about 0.1 ms to reach its second peak at t n2 ≡ t 5 . The height of the second peak is 0.6, greater than the value of the first peak. The duration of the second peak is also much longer than that of the first peak. n s returns to its normal value at about t = t 8 , indicating that the second peak lasts about 1.7 ms.
The temporal evolution of T es and T is is shown in figure 6(c). The solid and broken lines represent T es and T is , respectively. T es increases from 20 eV to its maximum value T es,peak ∼ 0.4 keV at t Te ≡ t 2 . After t 2 , T es begins to decrease monotonously and drops to 20 eV at 1.5 ms. T is remains almost constant from 0.0 ms to t = 0.12 ms, and then begins to increase. T es is greater than T is before t 4 (EDP phase). At t 4 , T es = T is , whereafter T is > T es . This is the start of the IDP. T is reaches its maximum value of 0.2 keV at t Ti ≡ t 5 , and then decreases. At t 8 , T is drops to 50 eV, at which point the PEP starts. In the PEP, T e and T i decrease slowly and monotonously to 20 eV.
The evolution of the plasma can be qualitatively analyzed by a two-component assumption. By assuming that each plasma species (namely, electrons and ions) is composed of hot ELM particles and cold background particles, n s can be written as n s = n c + n h , where c and h represent cold and hot particles, respectively. The sheath temperature is written as T s = T h ⟨n h ⟩ + T c ⟨n c ⟩. T h and T c are the temperature of hot and cold particles, respectively, and T h > T c . ⟨n h ⟩ = n h /n s , ⟨n c ⟩ = n c /n s , and ⟨n h ⟩ + ⟨n c ⟩ = 1.
The evolution of n s is influenced by electron and ion transport. As τ e,∥ ≈ 4 µs ≪ τ i,∥ ≈ 0.3 ms, the ELM electrons reach the target much earlier than the ions, which leads to a significant charge difference. The charge difference generates a strong electrostatic field, which attracts background ions and causes n s to increase. Namely, n c is enhanced by the fast ELM electrons, which is the mechanism of the first peak in n s . The evidence for this assertion is that the time of the first peak is between t 1 and t 4 (i.e. in the range of T e,shift ). In this time, n h ∼ 0 because τ i,∥ > t n1 . The second peak of n s is at t n2 = t Ti = τ i,∥ + τ r . Thus, the second peak is caused by the ELM ions from the source region. In other words, the tendency of n s to change after t = t 3 is caused by the evolution of n h . The evidence of this speculation is that the time of the second peak is consistent with T i,shift .
The present result is basically consistent with [22]. However, in the present simulations, the peaks of n s are completely separate from each other. The gap between the two peaks is about t 4 − t 3 ≈ 0.1 ms. The reason is that the selected magnetic field line is close to the separatrix surface and the connection length is long, as mentioned in section 2.1.1. The long connection length prolongs the difference τ i,∥ − τ e,∥ = Both T is and T es increase monotonously to their peak values and then decrease. The times of the peak values and the slopes either side are different, as determined by the transit time. The reason for t Te < t Ti is that τ e,∥ ≪ τ i,∥ . Different values of τ ∥ lead to T es,peak > T is,peak because hot electrons arrive at the target in a shorter time, which leads to a higher ⟨n e,h ⟩.
The consistency of t n1 = t Te agrees with our previous analysis that the first peak of n s is induced by high-energy electrons, because higher T es results in a higher sheath potential drop. However, t n1 is greater than τ r + τ e,∥ = 0.105 ms (the time of the maximum value of T e,shift ). The reason is that the electric field generated by fast ELM electrons accelerates ions and decelerates electrons, which slows down the electrons and delays the time at which T es,peak is reached.
For ions, the T is profile is significantly wider than T i,shift . There are two reasons. The first is that the length of the source region is not negligible (L + − L − = 30 m). Second, the initial velocity is not unique. The influence of these two facts is more significant for ions due to the higher transit time. This also explains why T is starts to increase at t = 0.12 ms, which is before t 4 . Some hot ions with very high velocities reach the target very early, but the ratio ⟨n i,h ⟩ is small compared with ⟨n i,c ⟩. As more hot ions arrive at the target, ⟨n i,h ⟩ increases and enhances T is . When ⟨n i,h ⟩ reaches its maximum value, T is also reaches a peak. This is why t Ti = t n2 . Thereafter, the temperature begins to drop monotonously to 20 eV because the ratio ⟨n i,h ⟩ drops. The simulation includes Coulomb collisions and the A&M process meanwhile the two-component assumption does not. Thus, the two-component approximation cannot quantitatively explain the evolution of n s and T s due to the absence of Coulomb collisions and the A&M process. However, this approximation still reveals the effect of the different parallel transport of electrons and ions in the SOL, which is the primary factor behind the tendencies of n s and T s . Therefore, the simplified model is still helpful for understanding the evolution of n s and T s .

Temporal evolution of electrostatic field with ELM
The electrostatic fields at t 1 -t 8 are obtained by the 2D3V PIC simulations, revealing the field evolution on the CDGs during the ELM burst. The 2D PIC simulations are performed using the PICS2 code. Several assumptions are made for the PIC simulations. First, the 2D PIC simulations are quasi-static. That is, n i , T es , and T is are constant at the sheath boundary in each case. Second, the electromagnetic perturbation induced by the evolution of the electron and ion flux is ignored. The reason for this treatment is that the time required to obtain the static state is τ static ∼ (ρ i /v ti )/ sin θ ∼ 80 ns, much smaller than the timescale of plasma evolution f/(df/dt) ∼ 0.01 ms (here, f is n s , T es , or T is ). The parameters are listed in table 3.
The time step in the 2D PIC step is set to dt = 0.2 ps to satisfy the condition dt < 0.1/ω ce ; dt also satisfies the condition that dt < T es /v te . According to our experience, the total simulation duration should be significantly greater than τ static because the velocity of the ions is not uniform. Thus, a onemillion-step PIC simulation is applied for each case. The Debye length λ D is listed to determine the grid size dx. The grid sizes in the x and z directions are equal, namely, dx = dz. In Case 1, T es = T is = 20 eV is applied for the pre-ELM phase and PEP. Cases 2 and 3 study the effect of different electron temperatures in the EDP. Case 4 considers the EDP→IDP transition. Cases 5-8 study the potential profile in the IDP.
The electrostatic field is shown by the vector arrays in figure 7. The length of the vectors is scaled by where E x and E z are the x and z components of the electrostatic field ⃗ E. There is a strong downward electric field on the MB top (z = 0 and 0.5 mm < x < 1.0 mm). This is reasonable because the Debye sheath forms and the electrostatic potential drops dramatically on the MB top. The electric field at the gap entrance (GE, located at z ∼ 0 and 0.0 < x < 0.5 mm) is also strong. The electric field at the GE points to the vertical gap surface, which drags ions toward it.
The most important difference between the EDP and IDP is the direction of the electric field at the GE. The electric field is downward in the EDP, while ⃗ E points upward in the IDP. The contour plots of E z are also shown in figure 7. In the EDP case, E z changes gradually (Cases 2-4) and is negative. There is a significant protrusion at the GE in the pre-ELM (Case 1) and IDPs (Cases 5-8). The value of E z (x = 0.25, z = 0) are given in each figure. Apparently, positive E z repels ions and influences the heat flux distribution on the gap surface.
The reason for this difference is the so-called potential peak at the GE [35]. Electrons are shadowed by the vertical gap wall due to the shorter Larmor radius ρ e , which leads to lower electron density at the GE. By contrast, ions with a larger Larmor radius of ρ i can reach the GE. The difference in the electron and ion density produces a positively charged region at the GE, generating a higher electric potential. The height of the potential peak is closely related to λ D , as [35] describes. The potential peak can be about 1.4T e with λ D ∼ 10 µm, forming a transport barrier inside the gap. For large λ D ⩾ 30 µm, the potential peak is almost invisible.
Our results are consistent with [35]. E z is negative at the GE in the EDP1, EPD2, and EDP→IDP stages, as the corresponding Debye length λ D = 66.6, 71.5, and 44.4 µm are significantly greater than 30 µm. In the IDP cases, λ D is between 20 and 30 µm, and E z at the GE is positive.
The positive E z (or potential peak) can be weakened when toroidal gaps exist. Simulations using a full 3D3V PIC approach have shown that electrons are transported to poloidal gaps via gap crossing, neutralizing the positively charged poloidal gap [34]. This effect is not included in present study. The reason is that electrons, which are constrained by magnetic field, are not likely to influence the heat flux distribution of the entire CDGs because the CDGs length (∼several centimeters) is significantly greater than the width of CDGs (0.5 mm).

Power loads on CDGs
The distribution of power flux P as a function of curve length L (section 2.3) is shown in figure 8. The MB top is the horizontal surface from A to B. The vertical surface from B to C facing plasma is named the wetted side. The bottom of the CDG runs from C to D. The shadowed side is from D to E.
The step interval of the flux distribution is δL = 0.045 mm. The flux on the MB top, P B − , is almost constant for all cases; B − is a short segment of the surface from B − δL to B. Similarly, B + is from B to B + δL. There is a sudden jump from P B − to P B + . P wet = P L∈(B,C) is the flux on the wetted surface.
The flux distribution in the pre-ELM (Case 1) is shown in figure 8(a). P B − and P B + are 2 and 4.2 MW m −2 , respectively. The ratio P B + /P B − = 2.1 is significantly smaller than

Pt cos θ
Pt sin θ = 1 tan θ ∼ 10, here P t = √ P 2 B + + P 2 B − is the total flux. The reason is that the Larmor radius of ions is close to gap width. On the wetted side, P wet decays exponentially from B + and vanishes at λ w = 0.5 mm. The wetted length λ w describes the depth at which the flux becomes zero.
The flux distribution in the EDP (Cases 2-4) is shown in figure 8(b). The most significant difference between the EDP and pre-ELM phase is that P B + is smaller than P B − , rather than greater. The power load P drops from P B − = 138 MW m −2 to P B + = 90 MW m −2 in EDP1. Additionally, the flux profile is flatter than that of the pre-ELM. The flux P C − = 44 MW m −2 is greater than zero. In other words, the flux distribution on the plasma wetted surface is lower and flatter. In the EDP2 case, P drops from P B − = 99 MW m −2 to P B + = 91 MW m −2 , and becomes 2 MW m −2 at C − . The wetted length λ w = 1.0 mm because the flux is zero on the plasma shadowed side. The flux distribution during the EDP→IDP transition is similar to that of the pre-ELM-P B − = 46 MW m −2 , P B + = 82 MW m −2 , and λ w = 0.7 mm.
The reason for P B + < P B − in EDP1 and EDP2 is the change of impact angle ϖ = arccos (⃗ n ·⃗ v/|⃗ v|) of the ions. ⃗ n is the normal vector to a surface and ⃗ v = (v x , v y , v z ) is the velocity vector, with x, y, and z components of v x ,v y , and v z , respectively. The normal vector is ⃗ n = (0, 0, 1) and (1, 0, 0) for the B − and B + surfaces, respectively. Thus, the impact angles on the B − and B + surfaces are ϖ B − = arccos (v z /|⃗ v|) and ϖ B + = arccos (v x /|⃗ v|), respectively. The heat flux of ions is proportional to P ∝ ∑ j v 2 j cos ϖ j ≈ N ⟨ v 2 ⟩ ⟨cos ϖ⟩, where j is the particle index and N is the number of particles. ⟨ v 2 ⟩ and ⟨cos ϖ ⟩ are the averaged total squared velocity and the cosine of the impact angle, respectively. The In the EDP, the  strong electrostatic field induced by the high T es accelerates v z , which reduces P B + /P B − . The strong electric field also enhances the wetted length λ w . The reason is explained in figure 9, which shows the trajectories of sampled particles with and without electric field. Figure 9(a) compares the different trajectories of ions in the electrostatic field ⃗ E = 0 (blue broken line) and ⃗ E = ⃗ E(EDP1)(red solid line). The initial positions and velocities are the same. Due to the downward electrostatic field, ions are dragged to the deeper part and the landing points move lower on the wetted surface, which enhances λ w .
By contrast, the electric field in the IDP enhances the ratio P B + /P B − because of the positive E z at the GE. P B − = 65, 24, 9, and 2.2 MW m −2 , P B + = 210, 74, 23, and 5 MW m −2 in IDP1, 2, 3, 4, respectively. The ratio P B + /P B − = 3.3, 3.2, 2.6, and 2.2 in IDP1, 2, 3, 4, respectively. This means that P B + /P B − gradually decreases as the peak of E z decreases. The reason is explained by figure 9(b), which shows the trajectory of an ion in the electric field E(IDP1). The ions at the GE are repelled by the positive E z and the landing points move upward, which leads to λ w is reduced. The shrinking of λ w makes the flux near the B points more concentrated.
Thus, although the heat flux is enhanced in both the EDP and IDP, the main mechanism is different. In the EDP, the high T es leads a strong electrostatic field, which accelerates the ions and causes a high heat flux on the PFC surface. In other words, the main mechanism of high flux in the EDP is energetic electrons. By contrast, T es < T is in the IDP and the effect of the electrostatic field is not significant. Therefore, the main mechanism of flux enhancement is energetic ions.
The distributions of heat fluxes carried by electrons are similar in different phases, as shown in figure 10. The ratio of the heat flux on edge over that on the MB top is P B + /P B − ≈ 1/ tan θ for all phases. The reason is that the Larmor radius of  electrons is much less than gap width. The ratio of heat power between ion and electron at the edge are P B + (i)/P B + (e) = 72, 12, 18, 70, 300 and 164 in the pre-ELM, EDP1, EDP2, EDP→IDP, IDP1 and IDP2 phases, respectively. The reason that the values of P B + (i)/P B + (e) in EDP1 and EDP2 are significantly smaller than those in other phases is that the electron temperature in EDP1 and EDP2 is significantly higher than ion temperature. Nevertheless, the heat flux of electrons is significantly less than that of ions. It is because that the hightemperature electrons produce strong electrostatic field with high potential drop, which accelerates ions. In other words, the kinetic energy of electrons translate into the kinetic energy of ions, as previous study indicated [30]. Thus, we ignore the heat flux carried by electrons in the following discussion.
The temporal and spatial distributions of heat flux P are shown in figure 11. The horizontal and vertical axes represent the time and distance, respectively. The flux direction is from top to bottom. The MB top, wetted surface, gap, and shadowed surface are labeled on the right. The flux distribution is obtained by the cubic interpretation of all the cases. The timescale is linear to clarify the duration of the phases. The EDP, IDP, and PEP are also labeled.
There are two peaks. The first appears at t = t 2 , marked by 1 and 2 for the B + and B − surfaces, respectively. The second peak is at t = t 5 and is marked by 3 and 4 for the B + and B − surfaces, respectively. The power values P 1 = P B − (t 2 ) = 89 MW m −2 , P 2 = P B − (t 2 ) = 140 MW m −2 , P 3 = P B − (t 5 ) = 192 MW m −2 , and P 4 = P B − (t 5 ) = 46 MW m −2 are labeled on the color bar. The first peak is introduced by the electrons due to the high T es in the EDP. The MB top suffers the most intense heat flux and the flux distribution extends deeper towards the bottom of the CDGs. The second peak is caused by the energetic ions in the IDP. A very high flux becomes concentrated in a very narrow region at the gap edge in the IDP. Figure 12 shows the heat flux evolution of B − and B + . Points 1-4 are the same as those in figure 11. The tendency of the flux on P B − and P B + is different. After reaching the  first peak P 2 , P B − drops dramatically to a low level (P B − = 46 MW m −2 ) at t = t 4 , and then increases slowly to the second peak P 4 . Thus, P 2 > P 4 . By contrast, P B + decreases slowly to 82 MW m −2 at t = t 4 after the first peak P 1 . Thereafter, P B + increases dramatically to P 3 , which is much greater than P 1 . The reason for this difference is that the first peak is induced by the high T es , which leads to a low P B + /P B − value. The second peak is induced by the high T is . This causes a high P B + /P B − .
The power quantity S (in MW) evolution as a function of time is shown in figure 13 to compare the magnitude of deposited energy for EDP, IDP, and PEP phases. S is obtained by the integration of the power load P from point A (L(A) = 0) to E (L(E) = 3.5 mm) along the whole toroidal direction P(t, L)dL (12) here R CDG is the R coordinate of the CDG and R CDG = 1.747 m. The ratio of the toroidal circumference to the toroidal extent 2π RCDG 1.5 mm is multiplied to the integration to estimate the power around the entire machine. The value of 1.5 mm is the x range from A to E ((L(B) − L(A)) + (L(D) − L(C))). The power width 2 mm is the distance from the strike point to the flux position on the divertor (in figure 3(a)). The total energy on CDGs in the EDP, IDP, and PEP phases is also given by The time points t 1 , t 4 , t 8 , and t 9 are denoted in figure 6 whose values are 0, 0.3, 2, and 3 ms, respectively. The total energy load during EDP is approximately equal to that in the IDP W(EDP) ≈ W(IDP). The reason is that the temperature of electron and ion in the source is equal, which implies electrons and ions carry the same quantity of energy. The energy load in the PEP phase is quite minor due to the lower plasma temperature in the phase. Experimentally, this ratio is closely related to the ELM frequency. As we simulate a single ELM burst, we cannot build the relation between F B + (total)/F B − (total) and the ELM frequency.

Conclusion
Multiscale three-step simulations were performed to study the temporal evolution of heat flux on poloidal CDGs. The simulations covered (1) the burst of the ELM and plasma transport in the SOL, (2) electrostatic sheath evolution on the gap, and (3) particle and heat deposition on the gap. The results reveal the relation between the ELM burst and heat load. The plasma transport, sheath formation, and PT were modeled by 1D PIC simulations with a speed-up technique, 2D PIC simulations, and 2D PT methods, respectively.
The 1D PIC simulations with the modeling of the ELM burst revealed the plasma temperature and density profile as a function of time. The results indicate that there are two special phases during the ELM burst. As electrons are faster than ions, electrons from the ELM burst reach the target earlier than the ions. In this phase, the electron temperature T e near the target is significantly greater than the ion temperature T i , which is defined as the EDP. After that, ELM electrons are absorbed by the target, which leads to T i > T e . This is the IDP. The mechanism for the difference between T e and T i is caused by the long connection length of the magnetic field, which has not been elaborated in heat load studies of ELM bursts.
Detailed 2D PIC simulations based on the plasma transport results indicate that there are significant differences in the sheath potential on the MB top in the EDP and IDP. First, the high T e leads to greater potential drop on the material surface. Second, the electric field direction on the gap is different in the EDP and IDP.
The different electric field profile changes the heat flux distribution, as illustrated by the 2D PT method. E z changes the velocity direction of the charged particles. In the EDP, the flux on the MB horizontal surface is enhanced by the high potential drop, and the flux on the plasma wetted region is widened. By contrast, the heat flux becomes concentrated within a very narrow region on the plasma wetted surface in the IDP. In other words, a significant difference has been found in the heat distribution during the EDP and IDP. The horizontal MB top suffers a more intense heat load in the EDP, while the leading edge faces the very high heat flux in the IDP.
In summary, there is a significant temporal evolution of the heat flux during ELM bursts. Regarding the design of the leading edge, the flux on both the MB top and edge should be taken into consideration.