First principles simulation of early stage plasma initiation process in ITER-scale tokamak

A first principles 6D kinetic model is developed to study the earliest times of unassisted plasma breakdown in an International Thermonuclear Experimental Reactor (ITER)-scale tokamak. This is then used for a comparative study of the predicted ionisation rate and the electron parallel velocity between the standard model for tokamak breakdown, assuming a zero-D (OD) Townsend avalanche, and the new kinetic model. The detailed model allows us to study the influence of the magnetic field configurations on the formation of plasma while explicitly resolving the electron trajectories. We introduce a ‘back-traced’ connection length L bt as a useful predictive tool for the spatial distribution of charged particles during the breakdown process. It is also found that the ionisation rate and the mean electron parallel velocity from the kinetic model generally exceed the 0D model predictions, demonstrating a growth in the total electron population from 103 to the order of 108 in approximately 1 ms. This implies that the 0D model can still serve as a conservative prediction for the first plasma campaign on ITER.


Introduction
Ohmic plasma initiation in a tokamak is achieved via the mechanism of Townsend breakdown ionisation of a pre-filled gas.During the breakdown phase of a tokamak plasma, free electrons are accelerated by a toroidal electric field induced from the rising poloidal magnetic field strength, causing additional electrons to be created through impact ionisation with neutral molecules that results in an exponential ionisation rate-the Townsend avalanche.In order to ensure an efficient start of International Thermonuclear Experimental Reactor (ITER) plasma operation, a predictive simulation of the breakdown process is helpful for navigating the narrow parameter space in which breakdown is possible and robust.
The theoretical model reviewed by Papoular is commonly used for simulating plasma breakdown in tokamaks [1].It describes various aspects of the breakdown process, such as loss as a function of the first Townsend coefficient α, the connection length L and the parallel drift velocity V de as seen in equation (1), = n e γ τ −1 ion = αV de ; τ −1 loss = V de /L. ( Upon closer scrutiny, the V de in the model is first assumed to be a constant function of electric field strength E and pre-fill pressure p, followed by fitting to experimental results.As a result, the formulated V de expression fails to consider the effects of magnetic field configurations, thus reducing the predictive ability for tokamaks that have different operating parameters.The connection length L is also an assumed fixed length depending on an averaged poloidal and toroidal magnetic field strength (B θ and B ϕ ), which neglects the spatial dependence of L across the poloidal plane.
Later works have attempted to address these issues.In order to improve on the simplistic ionisation rate description, a prescription of the rates of reaction types is introduced to determine the growth rate of certain ion species directly [4][5][6].While this is a definite improvement over Lloyd's model, the net ionisation rate is still not directly obtained via the emergent behaviour of the electron-neutral scattering events.Jiang et al developed a 1D numerical code that couples the 0D model with a Monte-Carlo method [7], removing the need for reaction rate definitions.However, the 1D simulation domain is unable to capture the main spatial dependency of the ionisation rate, e.g.via the field strength gradient along the radial direction.Even for models that consider the evolution of 2D poloidal flux surfaces through the Grad-Shafranov equation, the growth rate of charged particles is still essentially 0D [3,4].Any numerical model with lower spatial dimensions than 2D or three-dimensional (3D), requires several corrections: via drift terms to approximate particle orbits, approximation of the effects of pitch angle scattering, or approximate charged particle losses via the averaged stray fields.Such corrections, while physically motivated, essentially serve as refinements to numerical frameworks that still remain a simplification of the plasma breakdown in tokamak.
Additional uncertainty comes from the sheer scale of ITERlike tokamaks, which will have approximately eight times the plasma volume of the largest operational tokamak JET [14], which has a plasma volume of ∼100 m 3 [15].While the 0D model has been successfully validated against existing experiments, it remains to be seen whether the 0D model will be sufficiently predictive for tokamaks of such scale.The trend towards larger tokamaks will likely continue given the scaling of the Q factor as a function of the tokamak's major radius R [16,17], motivating the need to study the accuracy of the approximations of the 0D model.
Recently, de Vries et al revisited the topic of plasma initiation in general and utilised a fit of the elastic scattering cross section data of electron-H 2 by Tawara et al [18] as well as a force balance between the acceleration by background electric field and collisional drag to arrive at an updated V de in equation (2), This equation is then used in equation (1) for the calculation of electron population's rate of change.The derivation of equation (2) indirectly assumes a restricted degree of freedom in the velocity vector of electrons, imposing the condition that scattered electrons experience a drag force that is equal and opposite the accelerating force.The scattering angle distribution used in a 6D kinetic model is less harsh, allowing forward scattering which reduces the energy loss experienced by electrons.Therefore, the 0D model may overestimate collisional drag compared to a numerical model that simulates three dimensions in both space and velocity.The aforementioned treatment of L as a single isotropic quantity could also lead to difficulties in capturing the exact growth rate of the ionisation fraction.Altogether, there are a number of reasons to develop a numerical model that is capable of simulating the early state of tokamak plasma initiation, starting from a very low density charged medium (∼10 m −3 ) to study its rise in density over time due to the Townsend avalanche phenomenon in full three spatial dimension.Aside from providing a means to compare the 0D approximated growth rate as a function of L and V de , the detailed time development of the resulting plasma spatial distribution influenced by B θ as well as self-consistent electric fields is possible.This in turn permits an early prediction of the likely location of the first closed magnetic flux surface due to the growth of current densities.Finally, the developed numerical model also opens the possibility of explicitly modelling localised ECRH-assisted breakdown in future works.
A new first principles numerical solver is presented here for the purpose of simulating the earliest times of plasma breakdown, when the plasma current and subsequently generated self-consistent poloidal magnetic field is still negligible.The numerical study is conducted in a toroidal geometry, representing the breakdown region in an ITER-scale tokamak.This paper will briefly describe the numerical model in section 2 while section 3 focuses on the obtained results in comparison to the 0D equation in equation ( 1) as well as equation (2).

Numerical model
The numerical model developed for this work is a meshless particle kinetic code that resolves the electric and magnetic fields experienced by particles individually at particle's position.The simulated domain is bound by a simplified torus geometry which encapsulates the breakdown region centred around the magnetic null, shown in figure 1 along with the origin of the coordinate system.For future use-cases, the artificial toroidal boundary can easily be replaced by a torus with arbitrary aspect ratio.The parameters in table 1 are used throughout the simulations presented in the following sections.
All charged particles are individually simulated in our model, while the neutral molecules are not represented explicitly.At the beginning of the simulation, an initial population of free electrons with zero kinetic energy is seeded homogeneously within a smaller torus nested in the simulated domain.The electrostatic potential energy of electron at position x is determined in a straightforward manner via the electric potential Φ(x) obtained from equations (10) and (11).The smaller torus has a major and minor radius of 5.8 m and 1 m respectively.When a charged particle travels beyond the minor radius of 1.75 m, it is removed from the simulation and considered lost for further ionisation.A counter is implemented to record the number of lost charges throughout the simulation.
The numerical model simulates scattering and ionisation events between charged particles and neutrals via Monte Carlo method.The numerical treatment is described briefly in section 2.1.Newly created electron-ion pairs through electron-H 2 impact ionisation are added to the list of particles, to be simulated respectively as individual charged particle.All charged particles are subjected to the equation of motion described in section 2.2, while the experienced electromagnetic fields are described in sections 2.3 and 2.4 respectively.

Electron-neutral scattering
Our numerical study is conducted with H 2 as the neutral prefill gas with an initial electron number density of ∼10 m −3 .We employ the random acatter model developed in our prior work [19] to describe the charged particle trajectories through the pre-fill neutral gas medium, since it yields good agreement with a previous discharge experiment from Rose [20], in both the α/p value and electron parallel drift velocity at E/p ≈ 300 V m −1 Pa −1 .The E/p value is close to the operating conditions in ITER.
Cross-section data described in our prior work [19] together with additional cross-sections from the EIRENE database [22] for neutral dissociation of H 2 molecules are included to model the ionisation of the pre-fill gas.Since the charged particles will experience L > 2 km in the simulated ITER-sized tokamak, electrons that carry energies in excess of 1 keV will be present.The presented cross-sections in figure 2 are thus extrapolated towards the higher energy ranges assuming an exponential fall-off.The alternative of simply setting the crosssections in this range to zero is found to introduce an artificial hot, collisionless electron population.At every time step, the probability of collision P between electron and background neutral is calculated via the null collision method [23], expressed here as where ν T refers to the collision frequency that depends on the electron energy ε, while ∆t refers to the time step size.The expression for ν T is σ T refers to the total electron-H 2 scattering cross section, n denotes the number density of H 2 and ∥v∥ describes the velocity magnitude of the electron.At each time step, every electron will be checked for collisions by comparing a random number R 1 = [0, 1] to the computed probabilities, If equation (5) returns false, then the electron does not experience scattering in that time step.Otherwise, probability for all considered cross sections in figure 2 is calculated and compared to a second random number R 2 = [0, 1], the actual scattering outcome is then determined.Since the cross sections included in the simulations are not exhaustive, the sum of all cross sections σ sim,T in figure 2 is less than σ T for all ε.As such, the probability of an actual collision occurring computed with σ sim,T will always be lower than the probability calculated with σ T .The difference between σ T and σ sim,T is considered a null cross-section where the electron is treated as if there is no scattering event.

Charged particle equation of motion
The work presented here focuses on the plasma initiation up to the time of ∼1 ms, when the self-consistent poloidal magnetic field is negligible relative to the background stray field.Throughout the simulations, the charged particle motion is calculated using the gyrophase-error corrected Boris algorithm proposed by Zenitani and Umeda [24].The time discretised relativistic equation of motion is shown as where u = v/ √ 1 − ∥v∥ 2 /c 2 and vt = ( u t+∆t/2 + u t−∆t/2 ) /2. Equation ( 6) is not numerically solved as is, but computed through a series of equations shown here with A prominent feature of the Zenitani model is the analytical description of the charged particle's experienced θ angle rotation in equation (7).The obtained u t+ ∆t 2 is then used to update the particles' positions in the next time step via The ∆t used in the simulations is in the order of picoseconds to fully resolve the gyromotion of the electrons, dictated by the magnetic field strength in the order of 3 T to 5 T, for which the corresponding electron gyrofrequency is O(10 11 Hz).The chosen time resolution is much finer than what is needed to account for the expected plasma frequency of 3 MHz for electrons at a number density of 10 11 m −3 , predicted via equation ( 1) to be reached by ∼3 ms.This leaves the values E t and B t unresolved before the set of equation ( 7) can be solved numerically.E t and B t are the total electric and magnetic fields experienced by the charged particles at time t.Before the calculation in equation ( 7) is done, the spatial coordinate of the charged particles are used to first compute the experienced field through the description given in section 2.3 for E t and section 2.4 for B t .We considered implementing Littlejohn's Lagrangian guiding centre motion equations, specifically the extension for strong electric field cases by Morozov et al [25,26].However, the resulting guiding centre displacement is not exact in an inhomogeneous magnetic or electric field.Rather, truncated to the first order terms.We favour resolving the gyromotion of the electron in detail, in order to first create high fidelity results which can be compared against future approximating equations.Due to the high computational cost in each simulation, the numerical studies were conducted on the JURECA-DC and JURECA Booster supercomputers of Forschungszentrum Jülich [27].Over 10 million core-hours on these machines were committed to the simulations presented here.

Electric field calculation
The simulated low density plasma system will initially have localised self-consistent electric fields due to Debye lengths in excess of 6 m, rendering Debye shielding effects negligible for charged particles within the simulated breakdown region.This makes the evaluation of the self-consistent electric fields necessary throughout the simulation due to charged particles' spatial distribution.Therefore, we employ a particle-based tree code algorithm, the pretty efficient parallel Coulomb (PEPC) solver [28][29][30][31] to compute the electrostatic potentials experienced by all charged bodies, which also allows a seamless transition to a Coulomb-collision dominant plasma medium at higher charged particle density.The presence of a localised self-consistent electric field E int motivated the usage of a full 3D spatial domain, since the net electric field experienced by charged particles invalidates the toroidal symmetry assumed in a 2D simulated domain.
The PEPC solver is a highly parallelised implementation of the Barnes-Hut algorithm [32].An extensively simplified overview of the electrostatic potential calculation is presented here.Readers are encouraged to study the work of Barnes et al for the detailed accurate description of the algorithm.The calculation of the potential experienced by charged particle at x i is resolved via two formulations, depending on the spatial proximity of all the other simulated particles.The criterion that determines which of the two formulation is used in the computation of electrostatic potential is the multipole acceptance criterion (MAC).The implemented MAC in our model is with s denoting the spatial resolution of the particle grouping that contributes to the experienced potential at x i and d refers to the distance connecting x i to the centre of charge x coc of the aforementioned particle grouping.Charged particle entity that is within the grouping is indexed by j.In the case when equation ( 9) is met, the electrostatic potential contribution of a group of particles to point x i is calculated via where a, b = 1, 2, 3, r = x i − x coc , and d j = x j − x coc .The condition where equation ( 9) is not met, is when the particles are within the radius of approximately 1.75 m from x i .In such a case, the electrostatic potential is evaluated as a direct summation via with r j = x i − x j .Finally, the total electrostatic potential at point x i is a summation of each evaluated Φ n for each particle group, as well as the obtained Φ l .
In order to conduct a 1 ms duration simulation with a ∆t at picosecond order, approximately 10 9 timesteps are required.Evaluating E int alone consumes approximately 700 core-hours (or ∼90% of the time) per iteration on JURECA-DC system for 10 8 simulated charged particles.For this reason, the E int field is updated every ∼2 µs instead (update frequency of 0.5 MHz).The interval of 2 µs corresponds to the time taken by a 10 eV electron to complete a revolution around the torus geometry, at which the resulting E int is updated.The significance of 10 eV is suggested by Papoular as the energy barrier that electrons cannot overcome in an electronneutral collision dominated setting [1].The chosen interval saves computational costs while the E int is being updated at least six times more frequently than the predicted plasma frequency of ∼75 kHz.The predicted frequency originates from the time extrapolation of the electron number density based on equation (1), which should reach 7 × 10 7 m −3 at 2 ms.Such an evaluation frequency of E int is also more than sufficient to capture pitch angle scattering due to Coulomb collisions.For example, an electron population with assumed averaged energy of ∼100 eV and number density of ∼1.3 × 10 13 m −3 (Debye length of ∼6.6 mm and ω pe ≈ 200 × 10 6 rad s −1 ) has an electron-ion collision frequency ν ei ≈ 1.2 × 10 −2 Hz [33].
In contrast, the description of the background toroidal electric field E ϕ is straightforward.It is calculated via equation (12), while using a constant loop voltage V loop = 22 V for all simulation scenarios.As a result, the calculated toroidal electric field strength is ∥E ϕ ∥ ≈ 0.6 V m −1 for ρ = 5.8 m, through The value of ∥E ϕ ∥ and the prescribed pre-fill gas pressure p in the simulations are twice of ITER's eventual operation [34], in order to speed up the simulated plasma breakdown process while keeping the ratio in the range of 250 < E ϕ /p < 425 V m −1 Pa −1 with the given parameters in table 1.The increased ionisation rate is expected from the raised pressure p, following Townsend's proposed expression for the experimentally fitted α value with A = 3.83 m −1 Pa −1 and B = 93.6V m −1 Pa −1 the aforementioned fitting variables for H 2 gas.

Magnetic field calculations
The prescribed background magnetic field is assumed to be toroidally symmetric.As a result of this assumption, the poloidal field is treated as constant along the φ direction shown in figure 1.The calculations of the prescribed background toroidal and poloidal magnetic fields at the particle coordinate of (ρ, z) are done via In the above field calculations, we define the toroidal magnetic field B ϕ , which depends inversely on the radial distance ρ, and prescribed field strength B 0 explicitly.The poloidal field is decomposed into the field along the radial B ρ and vertical B z directions respectively [35].They are assumed to be generated from a set of N coil circular poloidal field coils, each parallel to the x-y plane and centred at x = 0, y = 0. Every coil has its defined radius R coil , vertical offset z coil from the x-y plane, and a predefined current I θ .K(k) and E(k) are the complete elliptic integrals of the first and second kind respectively.The coil arrangements are simple idealised setups designed to centre the null point at ρ = R 0 , z = 0: these parameters together with B 0 used in the simulation scenarios (runs 1 to 6) are tabulated in table 2. A positive sign for I θ implies a current flowing in the anti-clockwise direction along the poloidal field coils in the top-down view and the definition of B ϕ is such that its vector is clockwise in the same top-down view.The expressions of magnetic fields in equation ( 14) are divergence free in accordance with Gauss' law for magnetic fields.It is worth mentioning that the self-consistent magnetic field is not currently computed during the simulation, but can be easily evaluated a posteriori.The justification for neglecting it will be given later in section 3.5.2.
The first four scenarios (runs 1 to 4) consider a set of four poloidal field coils as opposed to the eight coils used in runs 5 and 6.Run 1 represents a scenario with expected stray field magnitude of 3 × 10 −5 T at the field null point.The other scenarios assume a near perfect null with B θ,null at 10 −8 T. Run 2 is identical to run 1 aside from the stated B θ,null , designed to study the impact of the depth of the magnetic null on γ and electron V ∥ distribution.Run 3 has twice the B ϕ than run 2,  The resulting magnitudes of the fields are summarised in table 3. B θ,null denotes the poloidal magnetic field strength B θ of the null located at ρ = 5.8 m, z = 0.0 m. ⟨B θ ⟩ describes the averaged B θ along the poloidal plane, ⟨B θ,edge ⟩ describes the averaged B θ along the blue circle in figure 3, and B ϕ,null is the toroidal magnetic field strength B ϕ at the null position.The blue circle corresponds to the torus minor radius of 1 m in table 1, inside which free electrons are initially seeded at the start of the numerical simulation.The free electrons are at rest and they are distributed homogeneously in space.Prior spatial distribution of the free electrons with higher concentration at the vicinity of the torus' minor axis is found to have negligible effect on the ionisation fraction growth rate.
For legibility, the parameters in table 2 are quoted in limited precision, and cannot be directly correlated to the values in table 3. Specifically, the positioning of the magnetic nulls in the case of runs 5 and 6 is highly sensitive to changes in R coil , z coil and I θ .For example, an increment of 80 A in I θ,5-8 in run 5 shifts the location of minimum B θ by up to 30 cm.Therefore, contact with the corresponding author to obtain the exact parameter set is recommended should there be an interest to reproduce the poloidal magnetic field.

Results and discussion
We now focus on the comparison of the values predicted by equation ( 1) and the results from the numerical simulations, i.e. ionisation rate γ sim and the V de,sim .Naturally, a direct comparison of 0D and 6D kinetic simulation results cannot be made: additional reduction of the simulation results is needed in order to make meaningful comparisons.However, the numerical study uniquely allows the effect of the B θ stray field structure and its impact on the spatial distribution of charged particles to be evaluated ab initio.Despite limiting the simulation time to a fraction of the completed breakdown time which is expected to be of order 8 ∼ 10 ms, this information already enables us to introduce a convenient method to predict the charged particle concentration leading up to this threshold.Finally, a method to determine the time duration where the reported results remain valid is presented.It is found that the spatial distributions of charge and current densities reach quasi-steady state (while continuing to grow exponentially) when V de,sim settles to a constant value after around ∼0.8 ms.This finding enables us to terminate the simulations at an earlier time, resulting in the simulation durations of approximately 1 ms and further allowing us to consider a range of scenarios in table 3.As a point of reference, the simulations for both run 1 and run 2 cost up to 3 million core-hours on the JURECA Booster module each.Run 3 consumed a higher amount due to the doubling of total time steps, as a result of the 2 × higher B ϕ .

Charged particle spatial distribution
The B θ field created from the definition in table 2 is found to play a major role in shaping the distribution of charged particles.Figure 4 shows the comparison of the charged particle density distribution between run 1 and run 5.The stark difference between the two scenarios originates from the poloidal magnetic vector field shown in figure 3, where charged particles end up following the stray fields and eventually exit from the simulation domain.Additionally, the combination of E ϕ orientation (anti/-clockwise) and the B θ stray field vector also influences the orientation of the charged particles distribution.In the example of run 1, E ϕ is in the clockwise direction when viewed from +z and electrons are accelerated in the opposite direction.The motion of the electrons confined by the overall B field yields a concentration along the topleft/bottom-right diagonal.By contrast, electrons will concentrate along the top-right to bottom-left diagonal if the E ϕ is in the opposite direction.It is also clear from the plot that there is a higher concentration of charged particles at the torus inner wall, where stronger local ∥E ϕ ∥ allows electrons to gain more energy, contributing to enhanced ionisation.This result is consistent with the reported outcome in Lloyd's work [2].
Since the charged particles are created via electron-neutral impact ionisation events, a localised concentration of charged particles results from the ionisation path of free electrons.In order to predict the concentration of charged particles during the plasma initiation process, one needs only to retrace the electron path.The resulting retraced electron path length is termed the 'back-traced' connection length L bt .Prior works have already discussed the influence of the connection length L on the confinement time of charged particles [1,[36][37][38][39].It is then a natural extension via equation (1) that a longer τ loss will improve the ionisation growth rate γ as well.While prior works usually treat L as a single measure that is spatially independent, retaining the spatial information in computing L can provide information on the spatial concentration of charged particles during the plasma initiation phase.It is worth noting here that the numerical integration is performed with the step length resolution of 1 µm, and the resulting error at the L bt of 10 km is in the order of tens of metres.
Since it was assumed that the B θ is constant along φ, a poloidal diagnostic plane is created at arbitrary ϕ angle, discretised according to a chosen resolution along ρ and ẑ direction.For every grid point of the discretised diagnostic plane, a streakline integral is performed along the prescribed background B field, in the direction that fulfils B • E ϕ > 0.
It is found that the L bt within ∼0.5 m radius of the null point along the poloidal plane can exceed 70 km for runs 5 and 6.This requires an excessively long numerical integration time to produce a full 2D map of the L bt .A more reasonable alternative is introduced that involves introducing a limit of 52.8 km at which the integration stops.The cap corresponds to the distance that is travelled unobstructed by an electron initially at rest, in ∼1 ms while experiencing continuous acceleration  from E ϕ .The duration of ∼1 ms is relevant since all scenarios are simulated up to approximately the stated time.It must also be mentioned that the scenarios labelled runs 1 to 4 have the maximum L bt below the introduced cap. Figure 5 shows the result of the computed L bt map for runs 1 and 5 respectively, which compares well with the results in figure 4. It must be noted, however, that the prediction for the spatial distribution via L bt is only valid in the time prior to the formation of localised closed magnetic flux surfaces.Should the numerical study extend beyond the end of the plasma breakdown phase, it is expected that the self consistent magnetic field resulting from the plasma current is sufficient to interact and alter the overall B θ field.Therefore, the duration of validity is determined later in section 3.6.Table 4 shows the averaged backtraced connection length ⟨L bt ⟩ for all scenarios.It can be observed that there are two distinct groups for the configurations with four poloidal coils (runs 1 and 2 and runs 3 and 4 respectively), while the configurations with eight coils stand on their own.Identical ⟨L bt ⟩ between runs 3 and 4 provides an indication that halving the B θ stray field magnitude and doubling the B ϕ strength (as tabulated in table 3) achieves the same effect, which should double the overall L bt value from runs 1 and 2. The trend is slightly different for runs 5 and 6, halving the averaged B θ for run 6 increased the ⟨L bt ⟩ by slightly less than twice compared to run 5.

Mean electron parallel drift velocity, V de
In order to make a meaningful comparison to the V de in equation ( 2), we define a comparable value for the parallel drift velocity in the simulation: V de,sim is taken to be the mean toroidal velocity component V ∥ of all electrons aligned with the background acceleration field vector Êϕ .A negative sign is given to V ∥ for electrons that are scattered in the direction opposing the acceleration vector.This definition approximates equation (2), in that it is a result of the force balance between acceleration via the background electric field and collisional drag.
The measured V de,sim for all scenarios after settling to a constant value are tabulated in table 5.An example for run 1 is shown in figure 6, where it has settled to V de,sim ≈ 2.74 × 10 6 m s −1 from 0.75 ms onward.Looking at figure 7, the median population of the electron parallel velocity magnitude is at <5 × 10 6 m s −1 , suggesting that electrons are impeded by scattering with background neutrals.Additionally, energy is  spent ionising neutrals during impact ionisation events which results in electrons with lower energy, further damping the growth of V de,sim even as the highly energetic electron population increases.Figure 6 also shows the calculated V de = 1.78 × 10 6 m s −1 from equation ( 2), at ∥E ϕ ∥ = 0.6 V m −1 and p = 2 mPa.This value is approximately 35.2% lower than V de,sim found in run 1.The cause for the lower V de arises from the assumption in the derivation of equation (2).Specifically, the collisional drag within the aforementioned force balance is based purely on the elastic scattering collision frequency, which is determined with the assumption that electron population has mean kinetic energy of less than 20 eV [37].This will be shown as a gross underestimation in section 3.3, where the assumed mean energy turns out to be at least a factor of four smaller than simulated values.It follows that an excessive collisional drag is retained in the force balance, thus resulting in an underestimation of V de compared to the result from the numerical model.
Comparisons of V de,sim with the other scenarios can be made via table 5, as well as the values for ⟨L bt ⟩ in table 4. V de,sim increases with ⟨L bt ⟩ and we find that the doubling of ⟨L bt ⟩ between runs 1 and 2 and runs 3 and 4 corresponds to an increased V de,sim by approximately 0.28 × 10 6 m s −1 .However, this trend does not apply to runs 5 and 6 despite the approximate factor of two difference in ⟨B θ ⟩.This suggests that the dominant influence on electron's V de,sim is the continuous addition of low energy electrons from ionisation events, rather than being dictated by small population of highly energetic electrons as a result of long L bt (shown in figure 8(b).As a side note, it is important to remember that direct comparison between configurations with four and eight poloidal coils is not well founded, because the single averaged ⟨L bt ⟩ does not reflect the much broader area with large values of L bt found only in configurations with eight coils.

Parallel electron velocity distribution
The change of the electron V ∥ distribution over time for run 1, starting from ∼0.78 ms when V de,sim has settled, is shown in figure 7(a).A key finding here is that the V ∥ distribution exhibits a time-independent, self-similar profile scaled by a single parameter, the electron population growth rate which will be discussed later.This is seen explicitly in figure 7(b), where the distributions have been normalised by their respective total electron populations, and are seen to overlap over a dynamic range of six decades.It is worth noting that this self-similarity is observed for all scenarios, provided that the respective V de,sim has settled.It is then possible to project the time evolution of the electron V ∥ distribution when the electron growth rate for a respective scenario is known, coupling the growth rate with the exhibited self-similar profile.
We can also study the influence of both B θ and B ϕ on the electrons' V ∥ distribution.To that end, each respective parallel velocity distribution f ( V ∥ ) is normalised by their respective total population of electrons.The resulting normalised distribution f is then plotted in figure 8(a) for comparison.
It is clear that the distributions fall into three distinct groups.Runs 1 and 2 coincide with each other, sharing the computed ⟨L bt ⟩ as shown in table 4 as well.A similar observation can be made for runs 3 and 4.However, it is shown a longer ⟨L bt ⟩ in run 6 compared to run 5 did not alter the distribution in a meaningful manner.A common trait seen in runs 1 to 4 is that a longer ⟨L bt ⟩ yielded a wider distribution, which can be explained by the longer confinement time of electrons, in turn recording the existence of electrons with ever higher energy.In the case of runs 5 and 6 which the V ∥ distribution closely approximates each other despite differing ⟨L bt ⟩, it is possible that there is a limit to ⟨L bt ⟩ above which other factors, such as collisional frequency between electrons and neutrals, play a dominant role in determining the overall electron velocity distribution.A more detailed parameter scan of ⟨L bt ⟩ would be required to reach a clearer conclusion.
The charged particle system at ∼1 ms still has not reached a density where Coulomb collisions dominate.As such a thermal equilibrium cannot be expected for electrons due to the lack of thermal diffusivity from electron-electron or electron-ion collisionality.This then hinders the fitting of the electron energy distribution with a Maxwellian distribution.As evidenced from figure 8(a), V ∥ is influenced by the static background E ϕ which causes a higher electron population in the +V ∥ direction.This feature is most notable on runs 5 and 6, which have longer L bt that improves the confinement of highly energetic electrons within the simulated domain.
Figure 8(b) shows the kinetic energy distribution for the electron population for each scenario.Recall from equation (1) that the assumed L is at 1 km.Considering the imposed E ϕ is approximately 0.6 V m −1 , the maximum electron energy derived from the 0D model would be approximately 600 eV assuming such electron experienced zero collisions throughout its lifetime.This is a stark difference from the plotted distribution, where the maximum kinetic energy is in the order of keV.This is a direct result of the higher local L bt for the respective prescribed B θ in this work.
Table 6 shows the mean electron kinetic energy for respective scenarios at the end of the simulation.Once again, the tabulated energy shows three groupings similar to figure 8. Specifically, runs 1 and 2 shares similar mean value, while runs 3 and 4 and runs 5 and 6 make up the remaining two groups.The grouping is consistent with the ⟨L bt ⟩ measure, indicating that an overall longer L bt raises the mean electron energy.Comparison between the simulation end time of runs 3 and 4 and the resulting mean electron energy suggests that   the mean electron energy also settles to a constant after V de,sim settled.This conjecture is made via a combination of observations, the first being the similarities in the L bt profile and ⟨L bt ⟩ measure between the two scenarios.The second is that an earlier termination of simulation time in run 3 compared to run 4, yielded a similar V ∥ profile in figure 8.The third and last is the observation that V ∥ has a self-similar profile over time, once the V de,sim settled.The self similarity suggests that the energy distribution of electrons remains the same over time and therefore the mean kinetic energy should also remain the same.Lastly, it is within expectation run 6 records a higher mean electron energy than run 5 due to the longer ⟨L bt ⟩, tabulated in table 4.

Ionisation growth rate
Figure 9 shows the rise of the total electron population throughout the simulation for respective scenarios.Once again, the growth rate shows three distinct groups, which corresponds to both V de,sim and the kinetic energy distribution plot.The tabulated growth rates in table 7 are obtained by numerically fitting the respective electron growth with an exponential function, i.e. the same form as the analytical solution of equation (1).Since the simulations also track the ion population from the ionisation events, its growth rate is recorded as γ i,sim while the growth rate of electrons is labelled γ e,sim respectively.Before making the comparison with the growth rate γ in the 0D model, it must be pointed out that the ionisation growth rate varies across the poloidal plane as seen in figure 4.However, the locality of the growth rate is not considered in equation (1).In order to make meaningful comparisons, the overall electron population growth rate is considered a proper approximation to the 0D model's γ definition as it removes the spatial dependence.
Given the operating condition of p = 2 mPa, ∥E ϕ ∥ = 0.6 V m −1 and L = 1 km, the resulting electron growth rate γ in equation ( 1) is calculated to be 8257.85s −1 .The growth of the electron population as predicted by the 0D model is also shown in figure 9.When compared to γ e,sim of runs 1 0D for given V e,sim (table 5) and γ e,sim (table 7).and 2, γ e,sim is approximately 28% higher than γ from the 0D model.This can be explained by a longer connection length experienced by electrons in the simulation, compared to the assumed 1 km in the 0D model.Aside from that, a 3D spatial domain allows more degrees of freedom in the velocity of electrons.This avoids introducing excessive collisional drag experienced by electrons since they are no longer forced to back-scatter and lose energy over time.It is now possible to compute the required 0D model connection length L 0D to recover the numerically fitted γ e,sim via equation (1), using V de,sim while retaining the α value.Figure 10 shows a comparison between the tabulated result in table 4 compared to the computed L 0D .It is immediately clear that L 0D is roughly constant throughout a range of γ e,sim with the respective V de,sim .This is a stark difference from the measured ⟨L bt ⟩, showing the 0D model's inability to correctly capture the relationship between connection length and the growth rate γ e,sim , and would require significant corrections in order to account for the the ionisation growth rate behaviour revealed by the 3D model.

Self consistent fields
Aside from the prescribed background fields B ϕ , B θ , and E ϕ which determine the trajectories of charged particles within the breakdown region, there are field contributions from the spatial distribution of charged particles and their motion.Since a background stray field B θ is ever present in all considered scenarios, a loss of charged particles is expected.Due to the mass ratio between electrons and much heavier ions, electrons are much more mobile.Thus, the created charged particle medium will initially acquire a net positive charge over time due to electron losses.Figure 11 shows the resulting E int at the time of ∼1.14 ms for run 1.The outward pointing vectors along the top-left to bottom-right diagonal indicates that the charged medium is net positively charged.Even though the magnitude of E int is still low compared to the externally applied fieldhere O(10 −3 V m −1 ), it is sufficient to influence the shape of the charged medium.Figure 11(b) shows the same case where the simulation setup is identical to run 1 but the calculation of E int is entirely excluded.The result is a higher density of charged particles near the simulation boundary, which is in stark contrast to figure 4(a).This shows the importance of resolving the self-generated electric fields E int even for a charged particle medium that has very low number density (∼10 6 m −3 ).

Magnetic field B int .
Unlike E int which is incorporated into the simulation runs, B int was not computed during the simulations.The main reason is the low current density built up to the time of ∼1 ms, creating a negligible B int while introducing additional algorithmic complexity and computational costs.Figure 12 is a post-processed calculation of the current density and B int taken at ∼1.14 ms for run 1.B int is computed via numerical integration of the Biot-Savart Law where ∆V denotes the grid volume of the uniform Cartesian grid and J(x j ) refers to the current density at grid point x j .Referring back to table 3 for the B θ,null parameter, all scenarios start with a minimum B θ magnitude of ∼10 −8 T. In order for B int to exert an appreciable influence on B θ , it should be of comparable magnitude.Therefore, at this stage of the simulation a growth of the current density of another five orders of magnitude would be required in order to trigger formation of localised closed magnetic flux surfaces.This also confirms the validity of the numerical studies of all scenarios up to the time scale of ∼1 ms.Even so, it is useful to estimate the time t crit.when the B int magnitude becomes comparable to the prescribed B θ , the prerequisite for the formation of closed magnetic flux surfaces, which in turn will dramatically increase the connection length, which could introduce an even faster ionisation growth rate onwards.This is discussed in the next section.

Validity range of the results
While the presented results are physically consistent and coherent, one of the major simplifications incorporated into this study is simply neglecting B int altogether.This limits the time regime for which the extrapolation of the presented results remains valid.More specifically, the dependence of V de,sim , γ e,sim , γ i,sim , and charged particle spatial distribution on L bt no longer holds when local closed magnetic flux surfaces begin to form.An estimate of the t crit time when this will happen is provided here.The main idea is to extrapolate the time evolution of the computed B int from t 0 = 0.8 ms onward.The choice of t 0 = 0.8 ms is based upon the observation that V de,sim has settled for all scenarios by then, implying that the γ e,sim growth factor has settled as well (from equation ( 1)).There are two observed features of B int that are shared across all scenarios.The first is that the B int is independent of the ϕ angle.The second is that the growth over time of the local ∥B int (ρ, z, t)∥ (in a poloidal plane discretised into 45 × 45 grid) resembles a simple exponential function.Making use of these observations, ∥B int (ρ, z, t)∥ can be extrapolated onto a single poloidal plane with: where C B (ρ, z) is a numerically fitted constant from equation ( 16), based on the time series starting from ∥B int (ρ, z, t 0 )∥.In order to extrapolate to the time when a closed magnetic flux surface is present, an additional assumption is introduced such that the B int unit vectors are time invariant.This means that a closed magnetic contour similar to figure 12(b) will emerge when is satisfied for all ρ and z.An example extrapolation (equations ( 16) and ( 17)) applied to run 1 is shown in figure 13.
It is obvious that equation ( 17) is satisfied for each (ρ, z) point at its own t crit.loc. .While it might take longer for all points to satisfy equation ( 17   The predicted t crit.also provides a limit up to which the presented results (L bt , γ sim , V de,sim ) remain valid, and will no longer apply to their respective scenarios at the time beyond t crit. .It is also worth pointing out that the t crit. is very early in a typical experimental tokamak pulse duration and that the predicted electron number density at t crit. is orders of magnitude lower than a typical low-density plasma (∼10 19 m −3 ) [40][41][42].

Conclusion
A 6D kinetic model study of the very early stage plasma initiation in ITER-scale tokamak from first principles was conducted.Simulations permitted quantitative comparisons of the numerically obtained mean electron parallel drift velocity V de,sim and ionisation growth rate γ sim to the 0D model.One of the major findings from numerical simulations of six magnetic field configuration variants is that both the γ sim and V de,sim are consistently higher than their counterparts in equation ( 1).This implies that the 0D model generally underestimates growth rates of charged particles for an ITER-scale tokamak.One explanation for the discrepancy is the inherent neglect of other degrees of freedom in the charged particle motion in equation ( 2), effectively exaggerating the collisional drag on accelerating electrons.This in turn diminishes the rate of ionisation in the 0D model, as well as suppressing V de .It is also found that the reported average electron kinetic energies in table 6 are higher than the assumed 20 eV in the derivation of equation (2).
Another significant difference is the connection length use in the 0D model.It is found that the resulting back-traced connection length L bt from the prescribed combination of B ϕ and B θ peaks well above the 1 km that is commonly used in the 0D model.This, in turn, yielded higher V de,sim , γ sim and lower t crit. .Having a magnetic configuration that yields higher L bt will naturally widen the V ∥ distribution as demonstrated in figure 8, due to longer confinement time of free electrons which in turn allows them to gain a higher amount of energy over time.This study also introduced the use of a poloidal map of the L bt as a new tool to predict the spatial concentration of charged particles over time.
A further important conclusion of this work is that the presented results can only be considered valid up to the time when closed magnetic flux surfaces are absent.Depending on the L bt map (as a function of B ϕ and B θ ), the influence of a self-consistent magnetic field B int needs to be accounted for at an earlier time.One other aspect that has not been studied in depth is the influence of wider areas with long L bt to the mentioned observable such as the V de , γ and V ∥ distributions.It is possible that a magnetic configuration with a wide area of long L bt , like in run 5 shown in figure 5, will lead to a faster ionisation rate and current growth compared to a run 1-like configuration of comparable ⟨B θ ⟩ and B ϕ .
The self-consistent electric field E int is also shown to have an appreciable influence on the charged particle's spatial distribution, even at an early time scale of up to 1 ms.In conclusion, this study has shown that while the 0D model is deficient in approximating the results obtained from the simulated ITER-scale tokamak, it can still serve as a conservative prediction model for the plasma initiation process.

Figure 1 .
Figure 1.LHS shows the sketch of torus geometry with the red dashed line being the torus' minor axis.The sketch on RHS is an exaggerated view of the blue circle that represents the simulation domain boundary along the x-z plane.The red circles represent the approximate projected location of the poloidal field coils in the configuration with four coils, while the red crosses are the additional coils for the configuration with eight coils.

Figure 3 .
Figure 3. (a) The poloidal magnetic field vector for run 1.(b) poloidal magnetic field vector for run 5.The magnitude of the poloidal field peaks in the order of mT.

Figure 4 .
Figure 4. (a) The electron number density ne distribution for run 1 at 1.1375 ms.(b) ne distribution for run 5 at 0.985 ms.

Figure 6 .
Figure6.V de,sim over time plot for run 1.At the last time step, V de,sim is calculated to be 2.74 × 10 6 m s −1 while V de = 1.78 × 10 6 m s −1 .

Figure 7 .
Figure 7. (a) The V ∥ distribution of electrons (f(V ∥ )) in run 1 at different times.(b) The corresponding velocity distributions scaled by the total electron population of respective times.

Figure 8 .
Figure 8.(a) Normalised parallel velocity distribution f(V ∥ ).(b) The normalised kinetic energy distribution of electrons for all scenarios, taken at the final time step of respective simulations.

Figure 9 .
Figure 9.Total net electrons over time for all considered scenarios.

Figure 11 .
Figure 11.(a) The E int due to charge imbalance between electrons and ions in run 1.(b) The electron density in the absence of E int in a setup identical to run 1.

Figure 12 .
Figure 12.(a) Calculated current density J. (b) The corresponding internal magnetic field B int at the time of 1.1375 ms for run 1.
), localised closed magnetic flux surfaces can form earlier.The closed surfaces will clearly alter the local confinement behaviour, drastically improving the local

Table 1 .
Simplified ITER-scale tokamak dimensions and operating conditions.

Table 2 .
Parameters for B calculation for all scenarios.

Table 4 .
Computed ⟨L bt ⟩ for all scenarios.

Table 5 .
V de,sim for all scenarios.

Table 6 .
Mean electron kinetic energy in eV at respective simulation end time for all scenarios.

Table 7 .
Numerically fitted γ e,sim and γ i,sim for all scenarios.

Table 8 .
Predicted t crit.andphysical parameters related to the maximum predicted local electron number density ne for all scenarios.Max.localne, m −33.4439 × 10 14 2.4927 × 10 14 1.4364 × 10 12 ωpe, rad s −1 1.0469 × 10 9 8.9068 × 10 8 6.7614 × 10 7 λ by charged particles.A conservative definition of t crit. is given in equation (18), to avoid overestimating the time at which L bt , γ sim , V de,sim are still considered valid t crit.= min{t crit.loc.}. (18) Applying the condition in equation (18), the resulting t crit. is approximately 2.75 ms.The exercise is repeated for all scenarios and the results are tabulated in table 8, along with other common plasma parameters at t crit.for reference.Three distinct groups are observed once again, with runs 5 and 6 expected to reach t crit.the fastest.The predicted t crit.reduces as ⟨L bt ⟩ from table 4 or γ sim values in table 7 increase.