Kinetic modelling of start-up runaway electrons in KSTAR and ITER

Understanding the formation of start-up runaway electrons (REs) is essential to ensure successful plasma initiation in ITER. The design of ITER start-up scenarios requires not only predictive simulations but also a validation of assumptions. The objective of this study is to strengthen the physical background required for predictive simulations aimed at ITER plasma start-up design, by validating the model assumptions. Through kinetic simulations, this study examines the validity of steady-state models for Dreicer generation under slowly-varying time scales relevant to plasma start-up and investigates the finite energy effect, commonly neglected, on the runaway avalanche growth rate. The research findings provide insights into situations where kinetic simulations are necessary. To secure a margin-of-control scheme without kinetic simulation, we suggest a strategy of scanning the Coulomb logarithm in fluid simulations as an alternative to predict runaway current takeover and avoid RE dominant scenarios. Ultimately, this paper seeks to offer a robust physical background, practically supporting the successful design of ITER start-up scenarios.

Original Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. less attention given to the generation of REs during the startup phase since an early start-up RE study [13]. A recent study on plasma initiation in ITER [14] indicated that ITER plasmas are at risk of producing REs due to the low prefill gas pressure required for the plasma burn-through, and consequently the low plasma density during plasma initiation. Relatively slow timescales for RE generation during the start-up allow for a control of the action [15], but this can be difficult to reconcile with maintaining the plasma at a specified current and temperature [16]. To ensure the plasma initiation in ITER, it is important to avoid runaway-dominant discharges where the RE current takes over the total current [15].
The ohmic start-up consists of three main stages: the breakdown, the burn-through and the plasma current ramp-up. The breakdown phase, characterized by the Townsend avalanche model, is where bound electrons begin to be dissociated from atoms [17]. A sustained increase in the plasma current demands the plasma burn-through-a phase where the main species and impurities are fully ionized. The burn-through criterion is that the ohmic heating power exceeds the total loss power consisting of radiation, ionization and others [18]. To satisfy the burn-through criterion with a low electric field, a low prefill gas pressure is required [14]. The final stage is the plasma current ramp-up, where the plasma current increases and the plasma equilibrium is controlled. Without auxiliary heating, strong fuelling after the burn-through phase can prevent the plasma temperature from rising due to the increased density, leading to current ramp-up failure [15]. On the other hand, delayed fuelling can make the plasma vulnerable to RE generation [15]. Due to this complexity, predictive simulations, as well as the validation of models, are necessary to design successful ITER start-up scenarios [14]. Recently, RE codes oriented on start-up applications have been proposed [16,19]. An optimal fuelling strategy needs to be developed to suppress significant RE generation and avoid failed burnthrough in ITER. Such an optimal fuelling strategy was studied for ITER [16]. The impact of impurity content and wall conditions on RE development was studied in the JT-60SA start-up RE study [19].
There are two RE generation mechanisms relevant for the start-up, the Dreicer mechanism [20] and the runaway avalanche mechanism [21]. Start-up RE modelling [15,16,19,22,23] has employed the fluid description where n RE is the RE density, S Dreicer is the Dreicer generation rate, γ ava is the runaway avalanche growth rate and τ RE is the RE confinement time. The Connor-Hastie (CH) formula [24] and Rosenbluth-Putvinski (RP) formula [25] can describe S Dreicer and γ ava , respectively. A current carried by REs can be estimated by approximating their velocity as the speed of light, implicitly assuming free acceleration of the REs.
where j RE is the RE current density, e is the charge of an electron and c is the speed of light. An expensive computational cost of the full electromagnetic plasma burn-through model [26] as well as a demand for flexible scenario development [14] has urged usage of the fluid description for start-up RE modelling. Models of the Dreicer generation rate were validated in dynamic scenarios with the plasma parameters fixed [27]. We aim to verify a method for extracting the Dreicer generation rate from the distribution function in dynamic scenarios where the plasma parameters are changing. Historically, the Dreicer mechanism and runaway avalanche mechanism have been studied by artificially decoupling each other. For example, Connor-Hastie solved the Fokker-Planck equation without a knock-on source [24]. On the other hand, Rosenbluth-Putvinski solved the Fokker-Planck equation in the zero-temperature limit [25]. It was later revealed in [28,29] that nonrunaway electrons produced by knock-on collisions barely contribute to the runaway avalanche growth rate, due to the scarcity of large angle collisions with the RP source [25] and the Chiu source [30,31]. However, in start-ups with moderate electric fields, the energy distribution of REs produced by the Dreicer mechanism can influence the runaway avalanche growth rate. It will be examined whether the steady state of the runaway avalanche growth rate is possible or not.
The paper is organized as follows. In section 2, we introduce a start-up RE model. In section 3, we confirm the adequacy of the analytic formulas (CH and RP) for the start-up RE modelling. We discuss certain circumstances where predictions from kinetic modelling may differ from those from fluid modelling. In section 4, we apply our analysis to KSTAR and ITER.

Kinetic model
The evolution of the electron distribution function is governed by the kinetic equation In this equation, F(p, µ) = 2π p 2 f(p, µ) represents the normalized distribution function. The momentum, p, is normalized by m e c, where m e is the electron mass. The variable µ = cos θ represents the cosine of the pitch angle θ. The variable τ = t/τ c represents the time, t, normalized to the relativistic collisional time, τ c = 4π ε 2 0 m 2 e c 3 e 4 n free lnΛnorm . Here, ε 0 represents the permittivity of free space, n free represents the free electron density, and lnΛ is the Coulomb logarithm. The subscript 'norm' indicates that it is the logarithm used for normalization. The parameterĒ = E/E c is the resistive electric field, E, normalized to the critical electric field E c = mec eτc . C{F} is the Coulomb collision term which represents the collision operator suitable for the startup RE modelling, comprising the Fokker-Planck operator, C FP {F}, and the knock-on source operator, S{F}. Other effects such as the synchrotron and Bremsstrahlung radiation were excluded because they have minor effects on RE development during the start-up. In certain regimes, e.g., near the critical regime [32], their role could be significant. They are, however, out of the operational domain, which will be handled in section 4.1. Near-threshold regimes will be investigated in section 4.2, but we will only discuss the Dreicer generation there.

Fokker-Planck collisional model
The Fokker-Planck model is appropriate if small angle collisions are more probable than large angle collisions [33]. However, in a weakly ionized system at low temperature, such as during the breakdown phase, inelastic hard collisions can dominate. This stage is not included in our kinetic study, since the confinement time is very short, and we assume all particles are quickly lost to the wall. The Fokker-Planck operator for electron-electron collisions, especially free electrons, can be linearized [34]. The linearized operator is suitable for the start-up RE after the closed flux surface formation. This is because the condi-tion´p ≫p th dpdµδF ≪´dpdµF M is met, where F M is the Maxwellian distribution, δF is the distribution of the perturb- 2T e is the thermal momentum normalized by m e c andT e = T e /m e c 2 is the electron temperature, T e , normalized to the electron rest energy. The linearized Fokker-Planck operator has a test particle part and a field particle part. The field particle part is often negligible due to its exp(−p 2 /2T e ) factor [35]. Of course, the field particle part in the large angle collision operator plays a crucial role in the RE amplification process.
The Fokker-Planck collision operator that describes collisions with free and bound electrons and scattering on partially stripped ions is given by It has three characteristic frequencies: ν ee S is the slowing down frequency, ν ee ∥ is the velocity space diffusion frequency, ν ee D is the pitch angle scattering frequency. The upper script ee represents electron-electron collisions. ν ei,I D is the pitch angle scattering frequency for electron-ion collisions or electronneutral collisions and the upper script I represents species of ions and neutrals. n I bound is the electron density bounded by species I.
Collisional frequencies are normalized by relativistic collisional frequency, 1/τ c ) n free lnΛ norm , (9) where β = v/c is the normalized velocity, γ = 1/ √ 1 − β 2 is the relativistic factor, X G = v/v th is the velocity normalized by the thermal velocity v th ≈ p th c, G(X G ) is the Chandraseckhar function, Φ(X G ) is the error function, Z I is the net charge, Z I 0 is the nucleus charge, I I 1 (y) is the first screening coefficient and I I 2 (y) is the second screening coefficient [36].
The ion scattering part of the collisional operator is described using the Lorenz scattering operator with the collisional frequency accounting for the screening effect of partially ionized impurities according to the Thomas-Fermi model [36].
Collisions with free electrons are represented with the collisional integral in the form given in [37], which interpolates between two important regions: the relativistic superthermal limit and non-relativistic collisions with the thermal Maxwellian population.
Collisions with bound electrons may play a role during the plasma start-up. In the Fokker-Planck part, we introduce these collisions in a form of supra-thermal relativistic collisional integral [34], with the collisional frequencies taken to be consistent with the Bethe stopping power [38] at high energy.
The dominance of small angle collisions is represented by the Coulomb logarithm. The free electron part, lnΛ free , takes the Coulomb logarithms in the form of [36] as a high p asymptote, which are asymptotically matched with the Coulomb logarithm [39] as a low p asymptote. From inelastic collisions, the bound electron part of the slowing-down and pitchangle scattering frequencies, lnΛ I bound , took the Coulomb logarithm from the Bethe model [36,38] as a high p asymptote taking into account the modified mean excitation energy [40] due to increased temperature. It is subsequently asymptotically matched with 0 as a low p asymptote. Note that the asymptotically matched logarithm in the bound electron part, rigorously speaking, is not exactly 0, but approaches 0 as p goes to 0. Free electrons cannot lose their energy through inelastic collisions with bound electrons if the incident electron's energy is below the threshold energy of excitation. However, using this logarithm may lead to an incorrect representation of the physics, especially in a low temperature system with many neutrals, where the collisional effect of neutrals on test particle dynamics becomes comparable to that of free electrons.

Knock-on collisional model
The linearized Fokker-Planck operator does not accurately describe the large angle collisions, which are essential for the avalanche multiplication of RE. To address this, we employ the Chiu operator [30,31] as the field particle source operator in our model.
While the RP source [25] assumes that incident electrons have infinite energy, the Chiu operator takes into account the finite energy of incident electrons (only approximating their pitch-angle to be zero). Both operators assume, however, that collisions occur with a stationary target, leading to an implicit constraint on the velocity of both the secondary electron and scattered electron after a collision to be much higher than the thermal velocity [29]. This constraint necessitates the introduction of a certain momentum threshold, p cut-off , above which close collisions are evaluated. We will discuss this in section 3.2. We replace n free as n free + n bound , which may result in errors due to a finite ionization potential at most n bound /n free .
In binary collisions, electrons are indistinguishable. A naming convention is adopted in which the electron with higher energy after the collision is referred to as the primary electron and the other is referred to as the secondary electron. Accordingly, the energy of a secondary electron cannot exceed half of the sum of the energy of an incident electron and the rest energy of an electron. Subsequently, the minimum energy of test incident particles required for large angle collisions can be obtained from the condition [29,32] where γ 0 is the relativistic factor of the test incident particle. Then, for the knock-on source operator we have where p 0 is the normalized momentum of the test incident particle, p 0 min is the normalized momentum of the particle with γ 0 = 2γ − 1, H is the Heaviside step function, δ is the Dirac delta function and F 0 (p 0 ) =´dµ 0 F(p 0 , µ 0 ) is the isotropic part of the normalized distribution. Σ(γ 0 , γ) is the Moller cross section [41] normalized by 2π r 2 e where re is the classical electron radius, . (12) In this work, the large angle collision operator only comprises the field particle part of the Chiu operator. The test particle part of the Chiu operator accounts for the energy loss of the test particles from the large angle collisions, which are approximately incorporated into the Coulomb logarithm in the Fokker-Planck operator. The dominance of the small angle collisions appears to justify ignoring the test particle term. However, rigorous justification without doubly counted collisions demands including the test particle term, as well as modifying the Coulomb logarithm [29], lnΛ = lnΛ − ln(cot (θm/2)), where θm is the centre-of-mass scattering angle corresponding to the cut-off, which is over the scope of this work. Thus, for simplicity, we neglect the test particle term of the Chiu operator.

Kinetic simulation in dynamic scenarios
We adopted the self-consistent kinetic solver [2] using the FiPy Finite Volume partial differential equation Solver [42] and implemented the kinetic model introduced in section 2. In this section, we examine the details of RE formation under start-up conditions. In particular, we will verify the adequacy of analytical Dreicer and avalanche theories for the given parameters.
The linearity of the governing equation (3) allows for the representation of the distribution function as a superposition of the Dreicer component, F Dreicer , and the avalanche component, Fava, creating the primary REs and secondary REs, respectively: Equation (3) should be decomposed in a manner that possess the essential physics of both mechanisms. The evolution of ∂τ F Dreicer should appropriately describe a speed-up process of nonrunaway electrons. Meanwhiles, ∂τ Fava should not omit secondary REs created from knock-on collisions.
If we determine F Dreicer as satisfying the Fokker-Planck equation without the Chiu operator and initialize it as the Maxwell distribution, the evolution of F Dreicer can be decoupled from Fava mathematically and in terms of solutions, Although the avalanche component makes a contribution to the Dreicer generation, as we will demonstrate in section 3.2, this contribution is negligible. As a result, only the evolution of Fava is influenced by changes in F Dreicer , Given that two mechanisms have been separately studied, it makes sense to write E and C FP parts in the above equations.
To ease the interpretation, without analytically studying two mechanisms with the momentum dependent Coulomb logarithm, we set lnΛ free = 20 in this section. The kinetic simulation can take into consideration the dependence of the Coulomb logarithm on momentum and other kinetic effects, such as the reduced screening effect. The steady-state computation from the kinetic simulation will incorporate these effects in section 4. Figure 1 displays the time evolutions of the plasma parameters. The plasma parameters change linearly over time, the electron density ne, just denoting n free from this section, increasing from 10 18 m −3 to 1.5 × 10 18 m −3 , Te rising from 150 eV to 300 eV, the effective charge Z eff increasing from 1 to 1.5 and E/Ec decreasing from 35 to 30 for 0.15 s. Their evolution is reversed from 0.15 s. The evolutions are relevant for those of the KSTAR ohmic start-up. The time scales for these evolutions are slower than the critical collisional time τ crit = τc · (p crit /γ crit ) 1.5 ≈ τc · (p crit ) 1.5 where p crit = 1/ √ E/Ec − 1 and γ crit = √ p 2 crit + 1.  The time scales are calculated as 1/T A = 1/A · dA/dt where A is the plasma parameter and T A is the time scale for its evolution. We solve equation (14) without the Chiu operator by initializing F Dreicer with the steady-state distribution function at t = 0 s. Figure 2 verifies the adequacy of the CH formula in dynamic scenarios. The final RE density predicted by the CH formula is in agreement with that predicted by the kinetic simulation at t = 0.3 s. However, there is a discrepancy before 0.2 s when the RE density is evaluated using p crit , which is often referred to as a runaway region boundary. This discrepancy can be attributed to the method used for extracting the Dreicer generation rate from the distribution function in the kinetic simulation. In the Connor-Hastie work [24] solutions in five regions of the phase space are found and matched asymptotically. An associated steady state particle flux is found. We borrow the nomenclature of [24] and roughly denote regions IV and V as phase spaces with p near p crit and above p crit , respectively. It is typically interpreted that the RE region boundary can be chosen arbitrarily. p crit is often taken to be such a boundary. REs in the upper part of region IV, where p > p crit is satisfied, can be considered runaways, but the leading order terms in region IV are pitch angle scattering and velocity space diffusion, not slowing down and electric field acceleration. Particles in region IV are not freely accelerated. Therefore, a more precise determination of region V, where electrons can be freely accelerated, is desirable in order to extract the Dreicer generation rate from the distribution function in a dynamic scenario.

Dreicer generation
A reasonable indicator for region V is a certain characteristic of the isotropic part of the distribution, or the 0th Legendre mode. After a straightforward calculation (see appendix), we obtain the form The isotropic part of the distribution function in region V is approximately constant along the momentum grid. Exact determination of the boundary between regions IV and V is difficult, but we can roughly estimate it as the perimeter where the isotropic part of the distribution function starts to vary. We therefore introduce a region-based critical momentum, p V , for the boundary between the regions IV and V. We systematically determine p V by fitting´dµF in the regions IV and V as lines and identifying the point where the two lines intersect as p V . Figure 3 illustrates snapshots of the 1D distribution function and the location of p V in the momentum space. This illustration is just an example and differs from the situation in figure 1. The In figure 2, the red dashed line denotes the RE density evaluated using p V (t). The RE density extracted from the distribution function with p V agrees well with the CH formula before 0.2 s. p V also supports the assumption in equation (2) that REs will be freely accelerated. Particles above p crit can be considered as REs. However, the decline in the RE density evaluated with p crit after t = 0.15 s indicates that many of them are not ultimately able to join the RE population and are decelerated when the Dreicer rate goes down. Figure 4 compares the corresponding Dreicer generation rates: ( dnRE dt ) p crit is higher than the steady state rate due to particle deposition in the upper part of region IV. ( dnRE dt ) pV has a better agreement with steady state results. Figure 5 illustrates the evolutions of F 0 in the dynamic scenario (dashed lines) and in the steady state (solid lines). F 0 in region IV evolves towards a steady state, enhancing ( dnRE dt ) p crit , as also seen in figure 4. F 0 in region V differs from the steady state. The decreasing isotropic part of the distribution function in the region V results from the increasing Dreicer generation rate. It is evident that a full kinetic simulation may be required to predict the runawaydominant start-up scenario.

Runaway avalanche
It is known that the choice of p cut-off in the avalanche simulations can affect the results. Figure 6 shows the sensitivity study for the runaway avalanche growth rate to the value of the p cut-off . The Te = 1 eV case corresponds to the Te → 0 limit,  and the Te = 300 eV case to a higher temperature relevant to startup. The normalized runaway avalanche growth rate,γava, was calculated using equations (15) and (19), We find thatγava is independent of the choice of p cut-off in both cases as long as it is lower than p crit . In the Te = 1 eV case, the absence of the diffusive process leads to a simple force balance determining the acceleration of the particles. If secondary electrons fail to exceed p crit they will continuously slow down. If they surpass p crit , they will become REs and contribute to RE amplification. If p cut-off is set too high, the missed population causes a decrease inγava.
In the Te = 300 eV case, particle diffusion near the p crit boundary may play a role. Despite appearing independent of p cut-off below p crit ,γava is still affected by the diffusive process which drives secondary electrons upward. However, this effect appears to be small. When S Dreicer ⩾ γavan RE , the evolution of F Dreicer can account for the Dreicer generation. If nonrunaway electrons produced from the knock-on collisions have a similar impact on the Dreicer generation as test particles, it means S Dreicer ≪ γavan RE and the Dreicer generation is negligible.
The observed trends from the sensitivity study, featured as using the equations (15) and (19) to evaluateγava, are consistent with previous research [28,29]. To use the Chiu operator, we approximate the Maxwell distribution as the Delta distribution, demanding p cut-off ≫ p th . We select p cut-off to meet the conditions p cut-off ⩽ p crit and p cut-off ≫ p th .
It is interesting to note thatγava depends on Te even though the difference is negligible. The Chiu operator accounts for the finite energy of the primary population so that the production rate of secondary REs from the knock-on collisions depends on the primary RE distribution function. The Fokker-Planck operator has a slight dependence on Te in the supra-thermal regime, leading to a small difference in the distribution function and thus a small difference inγava. Figure 7 shows the evolutions of γava in the dynamic scenario with p cut-off = 0.5 p crit . The evolutions of the plasma parameters are the same as figure 1. It's seen that γava is different from the RP formula due to the use of the Chiu operator in this startup scenario.
The entire population density n cap RE can be divided into two subpopulations: n cap,inf RE and n cap,fin RE based on whether their energy can be approximated to be infinite when calculating their knock-ons production rate, or not, correspondingly.
Avalanche growth Γ(p 0 ) is a product of the cross section and particles velocity β 0 = p 0 /γ 0 [32] . (22) Note that Γ(∞) = 1/(γ crit − 1) is for the infinite RE energy and is consistent with [25]. We consider RE belonging to n cap,inf RE if Γ ≈ Γ(∞), otherwise to n cap,fin RE . A boundary of two subpopulations is p∞. Figure 8 shows the averaged RE kinetic energy and ratios of the RE density. In stage I, the condition n incap RE ≪ n cap RE is not satisfied, due to the low averaged RE kinetic energy. Stage I occurs during seed RE formation due to increasing F Dreicer in equation (15). The evaluation of γava using p V results in a higher value than using p crit since the incapable population was less counted in the RE density (see figure 7).
In  Figure 9(a) illustrates that the knockon source replenishes particles in regions IV and V, building a higher level of plateau structure than F Dreicer at 0.3 s (see the black dashed line in figure 5). A moderate electric field during the start-up cannot accelerate particles quickly and prevents n cap,inf RE from being dominant. We introduce the total RE source rate created by primary REs with energy between p 0 crit and p as S 0 to estimate the contribution of n cap,fin RE to RE amplification, Figure 9(b) also suggests that n cap,fin RE participates in 40% of RE amplification. The finite energy effect has an impact on the increase in γava, which can persist until n cap,inf RE dominates n cap RE , though it is not critical. Therefore, γava can differ from the RP formula during the tokamak start-up.
Therefore, in startup cases where avalanches cannot be neglected, it may be necessary to consider the runaway avalanche process during stages I and II using the kinetic modelling.
A large initial peak in γava in figure 7, averaged RE kinetic energy, n cap RE /n RE and n cap,inf RE /n cap RE in figure 8 has just been determined by the initial distributions of F Dreicer and Fava. The effect of the initial distributions on the evolution of γava is blurred as the newly generated RE takes over the total RE population. In addition, the apparent breaks in γava in figure 7 and n cap RE /n RE in figure 8 take place just due to the abrupt reversal of parameter evolution shown in figure 1.

Application
In a real situation, the Dreicer generation should occur first to provide RE seeds. In our application, the increased temperature makes the Dreicer generation dominant.
In KSTAR, the trace amount of REs does not significantly affect plasma dynamics. It's the runaway-free discharge where the RE current is not zero but negligible. The plasma parameters may be deduced from the plasma burn-through simulations (without RE model). In section 4.1, we demonstrate the consistency between the RE dynamics described by the timedependent kinetic simulation and those computed by the analytical models in the runaway-free discharge. We also explore the possible role of impurities or main-species in RE production by running the kinetic simulation with the inferred parameters.
For the runaway-dominant discharge, self-consistent RE start-up simulations are necessary to infer the plasma parameters. In section 4.2, we will examine the nonstandard ohmic scenario, where the requirements on the prefill gas pressure imposed by the burn-through condition in ITER result in the plasma density being too low to hinder RE formation. We also discuss the intrinsic sensitivity of the Dreicer generation rate to the driving electric field and collisionality.

Standard ohmic discharge in KSTAR
We choose KSTAR experiment #12328 as an example standard ohmic case for our analysis. DYON code [18] was used to interpret the evolution of the plasma parameters.
The updated DYON code includes an improved plasma surface interaction model for carbon sputtering. The model considers the effective carbon sputtering coefficient during radial movement triggered by a robust positional control [43]. For example, the plasmas are limited to the outboard limiter at the initial phase, and then they touch the inboard limiter after the inward movement. It is noted that while the outboard limiter is composed of only three parts in a toroidal direction, the inboard limiter consists of sixteen parts fully covering the touched toroidal area. The carbon sputtering coefficient is set  as a function of time, while the recycling coefficient is set to 1.04. In KSTAR, the release of the main species can exceed retention, since the ohmic discharge follows the He glow discharge for wall conditioning [43]. Figure 10 shows the time evolution of the plasma parameters calculated using the DYON code [18] and the runaway current estimated from the time-dependent kinetic simulation solving equation (3). The timescales of plasma parameter evolutions are slower than τ crit when RE creation is crucial. RE losses are neglected since our interest is toward the phases after the closed flux surface formation. Figure 11 demonstrates the agreement between the RE density and current density predicted by the time-dependent kinetic simulation and those calculated using the steady-state kinetic simulation and the analytic formulas. The analytic formulas do not consider the reduced screening effect and momentum dependency of the Coulomb logarithm while these factors are included in the kinetic simulation. The final RE density and current density are relatively insensitive to the choice of p V . In fact, scanning p V from 0.2 to 0.3 resulted in errors of at most 8% and 2% in the final RE density and current density, respectively. The time-dependent p V can match the RE density before t = 0.4 s. However, this is not true for the RE current density because some REs have velocities less than the speed of light. Such a kinetic effect may result in a significant difference between the kinetic and fluid models when RE current dominates the total current dynamics.
Before 0.2 s, the plasma heating process increases Te and E/E D = E/Ec · Te/mec 2 , leading to an increase in the Dreicer generation rate. Near 0.2 s, the sputtering coefficient rises due to the inward movement of the plasmas, which are limited to the inboard wall. The enhanced radiative power due to the increased density of the partially stripped carbon species hinders the plasma heating and leads to an increase in the electric field. This, in turn, accelerates the Dreicer generation. After 0.27 s, the increased electron density due to the mainspecies recycling suppresses the Dreicer generation.
Penetration of low Z impurities into the plasmas can provide bound electrons, but they are likely stripped out after the closed flux surface formation, owing to the high temperature and low ionization potential compared to high Z impurities, such as Ne, Ar, Fe and W. This suggests that the reduced screening effect during the start-up after the closed flux surface formation is not as critical as during the disruption with high impurity contents. Consistency between the RE density predicted from the analytic formulas without the reduced screening effect and that calculated from the kinetic model supports this interpretation (see figure 11(a)). However, this does not necessarily mean that this is true for high Z impurities or during the breakdown and plasma-burn through phases.

Nonstandard ohmic discharge in ITER
The nonstandard ohmic scenario studied in this section is just an example in which the low plasma density, resulting from the absence of fuelling and insufficient main-species recycling, leads to RE current taking over the total current. Figure 12 shows the nonstandard ohmic scenario investigated using the SCENPLINT code [44]. REs are modelled selfconsistently, using equations (1) and (2) with S Dreicer and γava taken from [24,32] with lnΛ = 21, respectively. Until t = 3 s, Te increases due to ohmic heating, similar to the KSTAR case. However, due to a much higher Te in the ITER case, RE generation is significantly triggered. At t = 3 s, the increase rate of I RE took over that in Ip, indicating that the applied voltage is consumed by the REs. The applied voltage also heats the thermal electrons, but it is not enough to overcome the total energy loss. As a result, the RE current becomes dominant in current dynamics. Unlike the KSTAR case, the main-species recycling in this scenario does not raise the electron density or suppress the Dreicer generation. The ITER nonstandard scenario does not show a noteworthy role of impurity in RE creation.
The nonstandard ohmic scenario in ITER is characterized by the dominant Dreicer generation and slower time scales of plasma parameter evolutions. This implies that the timedependent kinetic simulation can provide the equivalent RE density to that modelled by the analytic formulas. However, the RE current density differs because of REs with velocities less than the speed of light. Kinetic modelling of the runaway-dominant start-up scenario requires a direct coupling of the plasma burn-through simulation and the RE kinetic simulation. This is beyond the scope of this work. However, the timing of the RE current takeover predicted by the selfconsistent kinetic modelling may be within a window of the timing obtained by scanning the Coulomb logarithm in the self-consistent fluid modelling. Increasing the Coulomb logarithm reduces the RE current density due to enhanced collisionality. We can also set the reference of the Dreicer generation rate using the steady-state kinetic simulation.
The sensitivity of the Dreicer generation rate to the Coulomb logarithm is shown in figure 13. The plasma parameters are prescribed by the nonstandard ohmic scenario over time. The reduced screening effect has little impact on the Dreicer generation rate. However, the Dreicer generation is highly sensitive to the Coulomb logarithms. Such sensitivity was previously predicted by Gurevich [45], as discussed in [24], due to the dominant exponential factor in the Dreicer generation rate. The timing of the RE current takeover can be as different as nearly 0.5 s depending on the Coulomb logarithm. Collisional cross sections have a logarithmically divergent nature in their integral and the cut-off of the maximum impact parameter must be chosen [46], in order to estimate the net effect of collisions. The choice of this cut-off produces ambiguity in the Coulomb logarithm. For example, the maximum impact parameter for thermal plasmas are often selected to be the Debye length [46], while for fast electrons it can be extended to the electron skin depth [36,47]. The minimum impact parameter is often selected as the closest distance of two colliding particles [33], but for fast electrons it can be replaced as the de Broglie wavelength in the center-ofmass frame [47]. Moreover, as discussed in section 2.2, separating large angle scattering components from the Fokker-Planck operator can influence the Coulomb logarithm. Owing to the nature of the Coulomb logarithm, none of the possible choices can be rigorously justified, and, therefore, the ambiguity in the results should be accepted. That is, the Dreicer generation rate is exponentially sensitive to the value of the Coulomb logarithm used in the collisional frequency, and thus it is also sensitive to the choice of the limiting impact parameters.
In addition, the Dreicer generation is also sensitive to the plasma parameters, indicating that a small adjustment to the plasma parameters can offset an error arising from the different Coulomb logarithms in self-consistent start-up RE modelling. Once the RE current takes over the current ramp-up rate, the plasma parameters can regulate themselves to keep the RE current rate close to the current ramp-up rate. Therefore, it can be determined whether scanning the Coulomb logarithm for scenario development is necessary or not by investigating the RE current takeover condition with a self-consistent start-up RE code, which is beyond the range of this work. Another issue is the profile effect on the Dreicer generation, as studied for the start-up and disruption [19,48], respectively, which is again beyond the range of this work.  Evolutions of the Dreicer generation rate in the steady state from the full kinetic model (black square marker), without the reduced screening effect (blue circular marker). A shaded area represents the Dreicer generation described by the CH formula with the Coulomb logarithm from 18 to 22. The green triangular marker shows the CH formula with lnΛ = 21, as implemented in SCENPLINT. The red dotted line is a takeover condition of the RE current over the current ramp-up rate due to the Dreicer generation, dIp/dt ≈ 0.9 MA s −1 .

Conclusion
The main results of this work strengthen the physical background for start-up RE modelling by employing the kinetic description of REs.
The steady state Dreicer generation models are clearly valid if the timescales of plasma parameter evolutions are slow. However, a correct interpretation is needed to infer the Dreicer generation rate from the distribution function in a dynamic scenario. We introduce p V and demonstrate that the Dreicer generation rate from the time-dependent kinetic simulation is consistent with analytic models in dynamic scenarios. Nevertheless, a kinetic simulation may still be necessary to model the runaway-dominant start-up where the approximation in equation (2) may affect the RE current takeover timing.
We emphasize two stages where γava can deviate from the RP formula due to the finite energy effect. During stage I the avalanche rate is reduced significantly compared to the RP model due to a very low average RE population energy. This stage is relevant to seed RE formation, especially when the Dreicer generation rate increases exponentially over time. During stage II, the finite energy of the knock-on capable population somewhat increases the avalanche rate. Kinetic modelling may be necessary when Dreicer generation and runaway avalanches are so significant that RE current rise can be comparable to thermal current rise. However, this is an exceptional case. Dreicer generation is significant in KSTAR and ITER applications.
The standard ohmic discharge in KSTAR is the runawayfree discharge. It is verified that the final RE density calculated from the time-dependent kinetic simulation is consistent with that predicted by the analytic formulas. Moreover, the final RE current density can be modelled with equation (2). The conclusion can be extended to other runaway-free discharges as well.
We also considered a nonstandard runaway-dominant scenario in ITER. The sensitivity study of the Dreicer generation rate to the Coulomb logarithm suggests that the choice of Coulomb logarithm can influence the timing of RE current takeover by nearly 0.5 s. Scanning the Coulomb logarithm can be useful to secure a margin of control scheme.
In this work, we only consider the phase after the closed flux surface formation, assuming that all of the REs are lost to the wall. However, this assumption can be relaxed by taking local confinement into account. During this phase, there can be weakly ionized plasmas with low temperature. In future work, we plan to model this phase while considering the finite ionization potential.