Numerical modeling of tokamak breakdown phase driven by pure Ohmic heating under ideal conditions

We have simulated tokamak breakdown phase driven by pure Ohmic heating with implicit particle in cell/Monte Carlo collision (PIC/MCC) method. We have found two modes can be differentiated. When performing breakdown at low initial gas pressure, we find that it works at lower density and current, but higher temperature, and requires lower heating power, compared to when having a high initial pressure. Further, two stages can be distinguished during the avalanche process. One is the fast avalanche stage, in which the plasma is heated by induced toroidal electric field. The other is the slow avalanche stage, which begins when the plasma density reaches 1015 m−3. It has been shown that ions are mainly heated by ambipolar field and become stochastic in the velocity distribution. However, when the induced electric field is low, there exists a transition phase between the two stages. Our model simulates the breakdown and early hydrogen burn-through under ideal conditions during tokamak start-up. It adopted fewer assumptions, and can give an idealized range of operative parameters for Ohmic start-up. Qualitatively, the results agree well with certain experimental observations.


Introduction
Tokamak plasma, like many other fusion plasmas, is inherently hot. However, the fusion plasma must be initialized from dilute neutral hydrogen gas at room temperature, to fully ionized plasma. This process, often referred to as start-up [1], is driven mainly by a central solenoid that supplies magnetic flux and induces a toroidal electric field in tokamaks, with the existence of external toroidal magnetic field. This process is often called Ohmic heating. In most practical tokamaks, start-up is not just powered by Ohmic heating, but also with auxiliary heating methods, like electron cyclotron resonance heating (ECRH) [2,3], electron Bernstein wave [4] or coaxial helicity injections [1,5].
The start-up stage is a complex problem but of less concern historically in tokamak research: it gets attention only if there is a failure. Even if start-up failure is common in tokamak operation (for example, there were over 100 non-sustained breakdown shots on JET experiments in 2009 [6]) it is not a major problem as there are large parameter spaces, so that successful start-up is always possible for most tokamaks. However, this problem has drawn much attention because of the emerging technical challenges faced by ITER [7,8]. As a superconductive tokamak, the maximum electric field for start-up at ITER is E < 0.3 V m −1 , at the lower limit of successful breakdown observed with purely inductive start-up in existing tokamaks [9]. Without fine-tuned parameters, reliable start-up of ITER will be hard to achieve. In fact, ITER-like superconductive facilities like EAST [10] and KSTAR [11] have some problems leading to unsuccessful start-up at initial operations [1]. For several years now, nearly all the experimental work conducted on tokamaks all over the worldincluding DIII-D [12] in USA, JET [13,14], ASDEX [15] and FTU [16] in EU, JT60 [17,18] in Japan, KSTAR [11] and VEST [19] in Korea, SST-1 [20] in India, HL-2A [21], EAST [1] and J-Text [22] in China-has been to investigate the optimum operational parameters in the start-up phase. Understanding and optim ization of parameters in the start up phase is critical not only for the future successful ITER startup scenario [7], but also for saving the valuable volt-second [9], which prompts more rapid burn-through and ultimately extends the pulse length of the discharge.
Despite the intense experimental study of the start-up stage, this critical phase is far from understood either theoretically or numerically. Sometimes the start-up failure can be attributed to hardware issues but at most other times, the cause is not known. Part of the reason is that most present theoretical works depend on simplified analytical or empirical models, which have limited prediction abilities. For example, albeit with great success, an empirical formula by Lloyd [23,24] has been used for decades for the breakdown phase, while no explanation or simulation has ever been given for it. Even with over 40 years of research, the parameter ranges of the initial plasma for most tokamaks are still not determined by numerical simulations, but by trial-and-error methods. A better understanding and a predictive simulation tool for tokamak start-up, especially for ITER, is urgently needed.
On the other hand, to understand the start-up process in tokamaks, different codes [25][26][27] are used to simulate the plasma initiation, even for the full tokamak discharge [28]. Generally, there are two kinds of methods to simulate the start-up process; sometimes two models can be coupled [29] together. The first method is based on the 2D equilibrium model [30], which calculates all the plasma parameters based on the particle equilibrium with electric and magnetic field structure. The second is based on the global transport model, which was pioneered by Hawryluk [31] 40 years ago and Lloyd [23] 25 years ago. Hawryluk and Lloyd's work has been extended and widely used in a large number of subsequent works. In fact, the model for burn-through and ramp-up phase is a recently developed 0D code called DYON [13,14,32], which includes the wall-plasma interactions and impurities. DYON has been used to predict some behavior in JET and ITER with great success.
From the point of view of gas discharge simulation, if the global model can be used to simulate a breakdown or discharge, the particle in cell/Monte Carlo collision (PIC/MCC) method can solve the same problem with more information. In fact, a global model assumes a Maxwell distribution, which is not appropriate for all partially ionized plasmas. Such a model has been tried recently in the 2D case [33] based on drift-kinetic approximation, but the results are limited as the toroidal directionwas omitted. In this paper, we have developed a new model and simulated the tokamak start-up process from cold neutral gas to fully ionized plasma. The outline of the paper is as follows. In section 2, the physical and numerical model is briefly explained. In section 3, we will study the start up processes under different conditions. Finally, our concluding remarks are made in section 4.

Physical model and parameters
As is well known, the tokamak start-up process consists of the plasma breakdown phase, the burn-through phase and the subsequent ramp-up process of plasma current until it arrives at a plateau state. The empirical formula by Lloyd [23,24] has been used for breakdown phase and been widely verified by many devices. The global transport model by Hawryluk [31] and Lloyd [23] has been adopted to study the burn-through phase, in which moderate plasma density and electron temperature are used as the initial conditions. Our model will treat the breakdown and early hydrogen burn-through phases in the same framework, under ideal conditions.
Simplification and approximation are the art of plasma physics. Actually, tokamak start-up is a complex 3D problem. In order to make the simulation possible, we have made several simplifications and approximations, following the work by Hawryluk [31], Lloyd [23] and Kim [13,14,32].
(1) We assume cylindrical tokamak and uniform plasma parameters. According to Hawryluk's [31] and Lloyd's [23] work, the space averaged or 0D model is a good approx imation for the start-up phase, because particle and energy transport is not the dominant loss mechanism. Therefore uniform plasma parameters can be used. The breakdown is determined along the axis in which the main magnetic and toroidal magnetic field point and the minor radius r; variations in magnetic field, error field, transport along r are all not included. However, to include the energy transfer between the ions and electrons self-consistently, we need a 1D model along the toroidal direction instead of 0D model. As we consider purely Ohmic heating, the plasma is homogeneous in the toroidal direction, therefore we have used a periodic boundary condition. In our model, the tokamak is treated as an infinitely long cylinder and only the field variation along the axis or toroidal direction is considered. In other words, our model is a short slice of the infinitely long axis of a cylindrical tokamak. As the field boundary condition is periodic, we adopt direct-finite-difference equation solutions [34]. We only consider the toroidal magnetic field B T and omit the poloidal magnetic field B P . We have found the results are the same if the poloidal magnetic field is included, because the poloidal transport process is omitted. In our simulations, B T is fixed at 2.3 T according to DIII-D [32]. Here B T is an important parameter, because it determines the velocity direction, especially after a collision. We also only consider the toroidal electric field, which contains two components: ambi . E ind is induced by the central solenoid and treated as a constant, and it accounts for the Ohmic heating, P oh . In fact, for any tokamak, E ind is always time-dependent as the impedance of the plasma changes, but adopting a constant value will make our results more universal. E ambi is the ambipolar diffusion electric field-the source of P equi , which represents the energy transferred from electrons to ions and is calculated self-consistently.
(2) We assume fully dissociated and pure H gas with argon cross sections. Therefore, we only follow the electrons and + H ions. The validity of adopting this fully dissociated H gas model has been discussed by Lloyd [24] and Kim [32]. This is not very accurate for the breakdown phase, but still does not significantly affect the results because the dissociation energy (∼1 eV) is small compared to the ionization and excitation energies (∼10 eV).
For simplicity and clarity, we used the argon cross sections, which include the electron elastic, excitation and ionization collisions, as well as the ion elastic and charge exchange collisions [35]. The cross sections are from the experimental data and have been widely verified in many works. There are two good reasons why we have used cross sections of argon instead of those of hydrogen: The first reason is that, in the context of gas breakdown and discharge, the effects of differences in cross section are small, because both the collision processes and Paschen's law are similar: the velocity directions of the incident particles are changed, part of the incident particles are lost by excitation, and new electron-ion pairs are created by ionization. According to our previous work and other similar simulations and theories for gas discharges [36,37], discharges of electropositive gases, like argon, helium, hydrogen, are always similar. The second reason is that, argon is the best studied monatomic gas. For gas discharges, argon or helium are always used as reference gases. The cross section data is the most accurate, easily interpreted and compared with other gases. Hydrogen, on the other hand, is a diatomic gas, which is much more complex to simulate and interpret. Our results can be a good reference for simulations of hydrogen gas in our following work. The electron excitation energy loss P iz , ionization energy loss P rad , as well as the ion elastic and charge-exchange energy loss, P CX and P elastic , are calculated self-consistently by the MCC method.
Before the tokamak begins to work, the hydrogen or deuterium gas is filled into the vessel in advance. As ionization proceeds, the neutral gas will be depleted and its temperature will increase because of collision processes.
In our model we consider the neutral depletion effect by counting the ionization in progress. The gas temperature is fixed at room temperature of 300 K. We have found the gas temperature does not change the results significantly. The validity of this assumption has also been discussed and verified in Lloyd's [23,24] and Kim's [32] work.
(3) We assume there is no particle loss and thus no plasmawall interactions (PWI), and consequently no impurities and no recycling, which are related to PWI.
Parallel and perpendicular transport of charged particles are always important factors which could influence the start-up process. In the breakdown phase, the parallel transport of electrons is dominant, and electrons escape more quickly than ions. This process is strongly affected by stray magnetic fields, which must be kept low so that the electrons can travel long distances to ensure sufficient collisions with neutrals before being lost to the wall. The stray field varies shot by shot, and can be characterized by a connection length, L f , which has the ideal value of infinity and the practical value of about 1000 m [23,24].
Here we take L f as infinite and omit all the particle loss. In the burn-through phase, the perpendicular transport of electrons and ions dominates and can be affected by × E B drift, ∇B drift and curvature drift, which will prompt charged particles to bombard the wall, leading to neutral recycling and release of impurities. Those processes we have omitted will be considered in future work; the reasons for their omission are as follows.
Firstly, it is preferable to have our results more universally valid for various facilities. Changing the gas pressure and induced electric field over ranges of values is always possible and easy for most devices, and the effects of the pressure and induced electric field depend only weakly on devices and on shot-to-shot variations. On the other hand, other factors, like initial profiles of stray magnetic fields, wall materials and conditions, strongly depend on devices and even vary significantly shot by shot. Clearly these factors will affect tokamak start-up significantly as well. However, if we consider these factors, it will make our results more reasonable and closer to the practical conditions, but will clearly lose universality. Therefore we have left them to our future works.
Secondly, all the processes we have omitted are energy loss mechanisms, therefore our model will give the most optimistic and bound limited estimations. It is critical to include the heating source correctly: for electrons it is the Ohmic heating and for ions it is mainly the power transfer from electrons; both are correctly included in the model here. When omitting parts of energy loss mechanisms, we will give optimistic estimations but the validity of the model is still maintained. Note the cross sections are the fundamental properties of the gas and are unchangeable, so we can always have accurate estimations of energy loss by collisions with the neutral gas, such as excitation or line radiation, ionization and charge-exchange.
Thirdly, some energy loss mechanisms, like transport loss and impurity-related terms, are less important in the early hydrogen burn-through phase before the neutrals are fully ionized. We can see clearly in [18] the ionization and radiation power loss through electron-hydrogen atom collision and the equipartition loss due to the power transferred are the dominant processes during the early burn-through phase, at least one order larger than other loss mechanisms in the first 20 ms. In addition, the electron temperature is very low, so the transport loss, i.e. thermal conduction loss, is small because the transport loss dominates only after the electron temperature rises, which will happen in the end of burn-through process. The reason is that at this early phase, the ions are just begin to bombard the wall, and the PWI and neutral recycling are just about to start. The fact that the electron transport loss dominates during the end of the plasma burn-through phase has been confirmed in [14]. Furthermore, in [13] Kim simulated a pure deuterium plasma and the results show that during the plasma burn-through phase the radiation and ionization loss are dominant in the total electron power loss, which we have considered in this work. Moreover, during the early burnthrough process deuterium radiation accounts for the majority of the radiation power [14] rather than the radiation power from the impurities. In short, the assumption that we neglect part of power loss and the effect of impurities is reasonable, especially in the breakdown and early burn-through phase that we consider here.
With the approximations above, our simulations effectively operate under ideal conditions, therefore our results are not quantitatively valid but qualitatively valid. The actual parameter space may be narrower than that we report here. But whatever the devices and shot parameters are, the actual parameters will always fall in the range that we will give here, as our results give the bound-limit parameters. The simulation results may not be sufficient to explain all start-up failures, especially in the burn-through phase, but can be used to explain some start-up failures under extreme conditions.
It is inspiring to analyse the start-up process after we have made such approximations based on Hawryluk's [31] and Lloyd's [23] models. Although their models are used to analyze the burn-through phase, we still need them in the breakdown phase in which there are also excitation, ionization and elastic scattering processes for electron energy balance and charge-exchange, elastic scattering for ion energy balance. The only power input is the Ohmic power. For a successful start-up, particle balance and energy balance criteria must be satisfied. For particle balance, the electron-pairs generation rate by electron-neutral ionization collision must exceed the loss rate of particle drift and recombination. In tokamaks, as the electron temperature is about 10 eV, the recombination plays a minor role and thus is neglected. So particle balance is always satisfied in this work.
For the electron energy balance equation, we have On the right-hand side of equation (2), P oh is the Ohmic power input, P iz and P rad are ionization and radiation loss, respectively, and P equi is the power transferred from electrons to ions through collisions. We have omitted all energy loss induced by transport, bremsstrahlung and impurities, as discussed above. The ion energy balance equation is as follows, where P CX and P elastic are the charge-exchange and elastic collision power loss respectively. Again, we have omitted all the energy loss from the transport and impurities, while we include the energy loss from the ion-neutral elastic collision. Note here we did not really solve the above equationswe calculated them self-consistently with the method we will describe below. Note in the above equation, P equi is the electron dragging force on the ions, so it is electron energy loss but ion heating mechanism. This term comes from the ambipolar electric field in the toroidal direction, due to the electrostatic waves excited by local charge separation.
We should emphasise that there is no arbitrary input param eter in our model, and all physical quantities are calculated from first principles. Only gas pressure, toroidal voltage/ electric field or stray magnetic field are controllable parameters in tokamak operation and they are variable input parameters. We treat B T and collision cross sections as fixed parameters, to which the results are not sensitive in our simulations. All plasma parameters shown below are given by self-consistent simulations.

Numerical model and parameters
As the plasma is strongly magnetized, most numerical simulations based on particle method for tokamak plasma are based on gyrokinetics, either delta-f or full-f. However, it is not applicable to use gyrokinetics methods for tokamak start-up. The first reason is that the closed magnetic surface has not formed and the magnetic field configuration is not in steady state yet. The other reason is that during most of the start-up time the collision processes are dominated by electron-neutral and ionneutral collisions. It is practically impossible to include the collision operators for these kinds of collision processes, as they are strongly nonlinear. However, the MCC model will be a better and more straightforward choice to include these collision processes.
The PIC-MCC model is a particle-mesh based method in which macro particles with finite size in a Lagrangian frame interact with fields defined on grids in an Euler frame. Here the particles and fields are updated in laboratory coordinates, instead of toroidal coordinates, as we use a 1D model along the toroidal direction with a periodic boundary condition. The PIC-MCC model is widely used for low temperature and space plasma, but is only occasionally used for tokamak plasma simulations. In our work, we adopt an implicit PIC-MCC method which allows us to use much larger space and time steps while keeping the same accuracy as an explicit one. Our model is based on some earlier work and the code has been tested in a lot of work [38,39] for low temperature plasmas, but has also been significantly modified and employed in tokamak simulations.
The implicit PIC method in this work is based on [40] and we apply weighted particles [38] as well as a uniform grid space, which means super particles are assigned to an identical number of physics particles and the number of macro particles per cell is identical, so the density in every cell is the same. However, with the discharge proceeding, the number of physical particles will increase to a quite large number to affect the efficiency of simulation, so we make a judgement: when the number exceeds a value, the number of particles will be halved and at the same time the weighting will be doubled.
For the MCC method, we have used Vahedi's model [41], and a standard MC procedure consists of the following steps:(1) choose particles stochastically participating in collision; (2) calculate energies and the corresponding cross sections of particles chosen; (3) choose the reactive type based on the cross sections and velocities of the particles; (4) calculate the scattering angle and azimuth angle; (5) calculate the velocity after collision. Since plasma density is relatively small, we consider only electron-atom and ion-atom collisions; short range ( λ < D ) coulomb collisions are negligible, but long range coulomb collisions are included self-consistently by the PIC model.
The physical simulation length is 6.4 cm, which is divided into 32 cells uniformly for all different input parameters, so the spatial step is 2 mm. The physical simulation time is 25 ms and the time-step is × − 4.0 10 10 s. Therefore one simulation will run for × 6.25 10 7 time steps, which will take about 2-7 d on a modern PC. Since we can run several programs on one PC at the same time, we did not parallelize the code. To valid ate the accuracy of our code, we have tested the effects of space and time steps, as well as other numerical parameters over a very broad range. We have found that resolution of the gyro-frequency does not affect the results, as we do not include particle transport. In addition, to reduce the noise, all the results below except the electric field are averaged across the whole simulation domain and 1000 time steps.  Figure 1 shows the results at × − 6.65 10 4 Pa, and figure 2 shows the results at × − 6.65 10 3 Pa, both with E = 1 V m −1 . The time evolution of the neutral gas, electron and ion densities is shown in figure 1(a) for = × − P 6.65 10 4 Pa and figure 2(a) for = × − P 6.65 10 3 Pa. There are common and different features for the two cases. At the beginning the electron and ion densities increase exponentially with the ionization, which means Townsend avalanche occurs. This phase lasts about 2 ms for both cases until the ionization rate becomes about 0.01 for = × − P 6.65 10 4 Pa and 0.001 for = × − P 6.65 10 3 Pa. Afterward, the plasma density still increases nearly exponentially, but with much smaller rate. This second phase is longer for = × − P 6.65 10 4 Pa, which means it will take a longer time to reach full ionization. The neutrals are nearly fully ionized at the end of the second phase, while the density of neutral gas decreases abruptly in a short time. The plasma density increases much faster for = × − P 6.65 10 3 Pa, which implies that the ionization is more effective for = × − P 6.65 10 3 Pa. So we can conclude that for high pressure operations the avalanche is shorter, which is in good agreement with Lloyd's theoretical derivation [23] about the avalanche duration and experimental results from JET [6]. After the second phase, the ion and electron densities barely change, since no recycling is included and thus no further ionization will occur. The density is also higher for = × − P 6.65 10 3 Pa by one order, which means for both cases the neutrals are fully ionized. Figure 1(b) shows the ' α D ' emission. It is not really ' α D ', but the excitation rate of the neutrals; however, it has the same implications as the real ' α D ' line and we can see that it is in good agreement with the experimental observations. ' α D ' comes from the de-excitation of the electrons, and the life time of the excitation state (<1 ns) is much shorter than the simulation time. So we calculate the electron excitation rate in short intervals and use it as the ' α D ' emission line. We can see that it appears as a density peak, similar to the experimental observations. This is because at the start the electron density is low, the α D emission is low, and the excitation increases with increasing density until the density, as well as the α D emission, reaches a maximum. Subsequently, the α D emission decreases due to the ionization of neutrals and the consequent depletion of electrons. We can also see that ' α D ' peaks at the same time as the second phase described above, where the peak for = × − P 6.65 10 3 Pa is much sharper than that for = × − P 6.65 10 4 Pa. This is because in the latter case the second phase is much longer.

Results and discussion
Comparing the electron and ion currents in figures 1(d) and (g), we can see that the electron current is always positive, while the ion current is positive initially, then becomes negative and increases linearly. The reason is similar to that of the density: The density is low, thus the effective E ambi is much smaller than E ind . Both electrons and ions are heated solely by the E ind . But their velocities are opposite as they have different charge. Thus the currents are both positive. When the plasma density is high enough, E ambi becomes dominant. Electrons move faster than ions, so the electrons drag the ions to move together, therefore the ion current density becomes negative. The current is mainly carried by electrons.
Comparing cases for the two different pressures of = × − P 6.65 10 4 Pa and = × − P 6.65 10 3 Pa, we can see that the plasma current increases much faster, and the current is one order larger for = × − P 6.65 10 3 Pa than for = × − P 6.65 10 4 Pa. This means that a larger electric field is required, hence more volt-seconds, to achieve such breakdowns. There are also many small peaks in the electron current for = × − P 6.65 10 4 Pa. These peaks have a duration of ∼ 2 5 ms, which can be understood as follows. The pressure is low, so that the collision effect is small, and most electrons are heated along the toroidal direction until they are heated to an energy above the excitation and ionization thresholds. After that, the electron energy will reduce rapidly due to excitation and ionization collisions, so the current decreases significantly and therefore a current peak will occur. No peaks appear for = × − P 6.65 10 3 Pa case, in which the collisional process is dominated.
From figures 1(e) and 2(e), we can see that the electron temperature growth is different for different pressures. At = × − P 6.65 10 4 Pa, E ind is dominant over the collisions, and the electron temperature can increase to above 100 eV at first in this phase, then decreases to about 9 eV, finally barely changing after that. At higher pressure = × − P 6.65 10 3 Pa, collisions dominate over E ind , so that the electron temperature increases to only about 9 eV and then decreases to 5 eV. The temperature increases to about 7 eV again and barely changes after that.
According to figures 1(h) and 2(h), the ion temperature keeps the trend of rising. At first, the ion temperature is low, because the E ind is not sufficient to heat the slow moving ions. The ion temperature significantly increases as E ambi increases and becomes dominant. So ion temperature increases as fast as the electron density and E ambi increase. After the plasma is fully ionized, the difference between the ion and electron temperatures becomes smaller, so the effective E ambi becomes smaller. The E ind drives the plasma current and heats the plasma again. So the electron and ion temperature increase slowly. , the ionization power loss P iz , the line radiation power loss P rad and the total collision power loss + P P iz rad . Initially, the Ohmic heating power is low, as the impedance of the plasma column is high, so the main current is the eddy current on the wall, not the plasma current. As the plasma density increases, the resistance of the plasma column decreases rapidly, so the heating power also increases significantly.
We also show the electron and ion total heating rate < ⋅ > E J T e and < ⋅ > E J T i , see figures 1( f ), (i) and 2( f ), (i).
In figure 1( f ) there are several steps which correspond to the spikes in figure 1(d). There is no step for = × − P 6.65 10 3 Pa. It is also clear to see that the ion heating rate is very stochastic, showing that E ambi and P equi are both stochastic. Ions are accelerated and decelerated randomly by E ambi , but ions are still heated as a whole, so the ion temperature keeps the increasing trend, see figure 2(h).
In order to figure out the kinetic behaviours of tokamak plasma in breakdown phase, we show the electron energy probability function (EEPF) and parallel velocity distribution function ( ∥ V DF) which is parallel to the torodial magnetic field in figure 3 at different times corresponding to the point A-D in figure 1 for = × − P 6.65 10 4 Pa. The Boltzmann equation describes the evolution of the distribution function, so the EEPF can be computed by solving it. When the collision frequency is constant, the EEPF is a Maxwell distribution. For a Maxwell distribution, EEPF can be expressed by where ε represents the electron energy, so g ln p is linear with ε. If the collision cross section is constant, the EEPF is a Druyvesteyn distribution [42] in which the most energetic electrons lose energy more rapidly. Here A is the fast avalanche, B is the slow avalanche, C and D are the slow rampup phases. The corresponding ion energy probability function (IEPF) and parallel velocity distribution function velocity distribution function are plotted in figure 4.
As can be seen, the electron energy probability function is close to a Maxwell distribution at point A. Since A is close to the end of the fast avalanche phase, electrons are heated by E ind , so that the velocity distribution is two-temperature drifting Maxwell distribution, with the significantly fast drifting component and a smaller slow but opposite drifting component. This slow but opposite drifting component may come from the E ambi heating, which makes the ions drag the electrons. In figure 4, we can see that the ion energy distribution is also close to Maxwellian, but a small fraction (<10 −4 ) of high energy ions can be observed, due to the E ambi heating. The ion velocity distribution is a one-temperature drifting Maxwellian distribution, but with positive drifting velocity. This implies the ions are still heated by E ind .
Point B is in the slow avalanche phase. The electron energy probability function is not a Maxwell but Druyvesteyn distribution, which also implies that E ambi becomes the dominant heating or cooling source. The ion energy probability function is close to a Maxwell distribution, but slightly skewed towards to lower energy, which comes from the ion charge exchange collision. The ion velocity distribution is a one-temperature drifting Maxwell distribution, but spikes begin to appear, induced by the E ambi or P equi . From figures 1(i) and 2(i), we can also see that the heating of ions by E ambi is stochastic. From the E ambi at different times, we can also conclude that the amplitude of E ambi is about 100 V m −1 , its wave length is several mm, and it is changing fast and stochastically.
Points C and D are the slow ramp-up phases, where the electron energy probability function comes closer to a Maxwell distribution. The velocity distribution is still a twotemper ature drifting Maxwell distribution, but the drifting velocity is increasing and drifting velocity differences between two electron groups become smaller. The ions are still in Maxwell distribution, but their velocity distribution becomes very stochastic. We can see in figure 6 that the ion drifting velocity increases significantly from point C to point D. Figures 5 and 6 show the EEPF and IEPF, respectively, at different times corresponding to points A-D in figure 2 for = × − P 6.65 10 3 Pa. Comparing figure 3(a) with figure 5(a), we can clearly see that in the fast avalanche phase the latter has larger energy probability. This can be explained by the observation that at the lower pressure, the collision rate is low, so the induced electric field is dominant rather than collisions, as we have mentioned above. In addition, we can see that at point A the EEPF is a Druyvesteyn distribution. Nevertheless, it becomes close to a Maxwell distribution at point B-D, which is rather similar to that in figure 5-while the electron ∥ V DF is quite different from that in figure 5. In this case the parallel velocity distribution of electrons changes from a one-temperature drifting Maxwell distribution at point A to a two-temperature drifting Maxwell distribution at points B-D. Besides, the drifting velocity is increasing and the difference between two electron groups becomes larger. Figure 3 also shows the averaged ambipolar diffusion field at different times corresponding to the point A-D in figure 2, and in this     It is well known that 0.3 V m −1 is the upper limit in ITER start-up due to the limitations on the superconducting poloidal coils, and this is also close to the lower limit of successful start-up by pure Ohmic heating. So it is meaningful to simulate the start-up phase under such conditions. In order to easily compare the effect of induced electric field, we also choose the same gas pressure of × − 6.65 10 4 Pa and × − 6.65 10 3 Pa. Figures 7 and 8 show the kinetic plasma behaviours and the corresponding pressures are × − 6.65 10 4 Pa and × − 6.65 10 3 Pa respectively. We can clearly see that both in figure 7(a) and 8(a) there are three different increasing rates, which is a big difference from figures 1(a) and 2(a). This can be explained by the observation that when the induced electric field is low, E ambi increases slowly and at the second phase it is still smaller than E ind . On the other hand we can say E ambi is closely related to E ind . Furthermore, the increasing rate during the third phase is much larger even than that during the first phase in figure 8, which is different from the other three cases. The reason may be that during the second phase, the collision frequency is low because the energy of most electrons does not reach ionization energy, which we can obtain from figure 8(e). Thus electrons are cooled at the start and then heated till the third phase. In the third phase, most electrons happen to collide with neutrals, in either excitation or ionization collisions, which makes the plasma density increase fastest. Besides, the breakdown is greatly delayed as the pressure is increased comparing figure 7(a) with figure 8(a) and as the induced electric field is decreased comparing figure 2(a) with figure 8(a). We can also obtain these results by comparing the occurrence of the ' α D ' peak, which is the typical definition of breakdown time [23].
In addition, from figures 7(c) and 8(c) we see that in both cases the total electron collision loss overtakes the Ohmic heating rate at some time: this is because in these two cases the induced electric field is small, which means the Ohmic heating is small, so the breakdown is not as efficient as the cases in which the induced electric field is 1 V m −1 . Furthermore, the total electron collision loss in figure 8(c) surpasses the Ohmic heating rate by much more than that in figure 7(c), which is because of high collision rate resulting from the high pressure. Accordingly, the maximums of electron and ion temperature are lower than those in figures 7(e) and (h).
However, the evolution of ion temperature in figure 8 is much different from the other three cases. In figure 8(e) we can see that the electron temperature decreases at the second phase and this trend is sustained for a long time while the ion temperature begins to increase slowly and then increases very fast-nearly instantaneously. Here at the second phase electrons are cooled, which means the right-hand side of equation (2) is less than zero and ions are heated greatly by P equi , which can also be confirmed from figure 8(c) where the total electron collision loss is much more than the Ohmic heating rate. Comparing the electron current density in figures 7 and 8, we can see that in the former case there are some peaks, which is similar to that in figure 1. Accordingly, the electron heating rate in figure 7 has several steps, similar to that in the figure 1. This is because the trends of electron and ion heating rate in figures 7 and 8 are similar to that in figures 1 and 2, except the magnitude.
In order to better understand the plasma characteristics, the EEPF, IEPF, electron and ion ∥ V DF, averaged ambipolar diffusion field at different times corresponding to the point A-D are shown in figures 7 and 8. In comparison, electrons can be heated much more efficiently during the first breakdown phase when the pressure is lower, which we can obtain from the EEPF in figures 9(a) and 10(a), so the first breakdown phase is quickly finished. Subsequently, electrons are cooled owing to the collisions with neutrals and then the temperature rises slowly as the breakdown proceeds (see figures 9(b)-(d) and 10(b)-(d)). In addition, the ambipolar diffusion field is much smaller in figures 9 and 10 than that in figures 3 and 5 during the first breakdown phase, which proves that the induced electric field can affect the ambipolar diffusion field.
From the IEPF, ion ∥ V DF in figures 11 and 12, whatever the gas pressure is, the ion temperature is rising all the time, which means ions are heated in the whole process. Besides, we can see that the ion heating is more effective under lower pressure when the induced electric field is the same. In figure 12(a), we cannot see the IEPF, because of the low induced electric field. The Ohmic heating is small and the ambipolar diffusion field which is the heating source for ions is quite small (see figure 10(i)) at the beginning, so ions are not heated effectively. Subsequently, electrons are heated by the induced electric field to have higher energy than ionization and excitation energy, and collisions happen, which makes the P equi increase so that ions are heated.

Comparisons of experiments and theoretic analysis
Theories and simulations have their merits only after they have been verified by experimental observation. Unfortunately, like many previous works, only very limited results here can be compared with experiments directly. There are two reasons for this: the first is that most of the diagnostic tools in tokamaks are designed for much denser and hotter plasma, which are unavailable in the early phase; the second is that we have made many assumptions under ideal conditions with the kinetic method for the first time. For a reasonable comparison, specific devices and working conditions will be needed, which will be considered in our future work. However, we can still make some comparisons qualitatively.
(1) The breakdown phase can be finished during several milliseconds and breakdown delay is also observed with the increasing of gas pressure and with decreasing of electric field(see figures 7(a) and 8(a)), which agree with the experimental data in DIII-D [23] and JET [6]. (2) The fact that ' α D ' peak occurs during the breakdown phase is also consistent with experimental observation in [23].
(3) Our results show that both ion and electron temperatures are less than 10eV at the end of the breakdown phase due to the thermal ionization barrier, which is in good    agreement with theoretic analysis and experimental observations [43]. In addition, even during the early burn-through phase, the electron temperature saturates to a constant level below 10 eV for a pure deuterium plasma [14], which is also consistent with our result that the electron temperature remains unchanged during the early burn-through phase. (4) The plasma current is mainly carried by electrons. When ' α D ' peaks, we can see that the current amplitude is comparable to the typical experimental data [44].
As for the electron and ion heating rate, because we have considered the ambipolar diffusion field self-consistently, which has never been considered in other work, so it cannot be compared. However we can clearly see the state of electron and ion heating and the from the result it is reasonable to conclude that after the breakdown phase the collisions decrease greatly, then the heating rate reaches its maximum. In short, to a certain extent, our simulation results qualitatively agree with the experimental observations and theoretic analysis.

Conclusion
In conclusion, we have simulated the breakdown process using a one-dimensional implicit PIC-MCC method. The behaviors of electron, ion and neutral gas densities, α D emission, radiation loss, ionization loss, total electron collision loss and Ohmic heating rate, electron and ion current density, averaged electron and ion heating rate, electron and ion temperature are obtained for different electric fields and pressures. For the first time, the simulation work gives much information to explain the physics of tokamak start-up, thus might open up a new paradigm in simulating the tokamak start-up process.
According to the variation trend of plasma parameters shown above, two modes of tokamak start-up can be distinguished. For specific tokamaks, the actual mode depends on the size of the facility and the power of the Ohmic heating coil.
(1) Low pressure mode; If the pressure is low, E ind is dominant over the collision process. The electron avalanche and density increase are slower. There are small peaks during the electron current increase. The electron temperature can increase to high values and then decrease to less than 10eV. After full ionization, plasma density and current are relatively low, which will require smaller heating power. Small tokamaks may run on this mode. (2) High pressure mode; If the pressure is high, E ind can be balanced with the col lision process. The electron avalanche and density increase are much faster. There are no peaks during the electron current increase. The electron temperature can increase to several eV, then decreases and finally increases again to several eV. After full ionization, plasma density and cur rent are an order higher than those at low pressure, which will require much greater heating power. Large tokamaks like JET will operate on this mode.
For the whole simulated process, three stages can be further distinguished; between the first and third one there can be a transition phase: (1) Stage I is the fast avalanche stage. This stage is analogous to the the breakdown phase. In this stage, the plasma density is low but increases exponentially and fast. The density is low and thus effective E ambi is much smaller than E ind . Both electrons and ions are heated solely by the E ind . Their velocities will be opposite as they have different charge, so that the currents are both positive. The electron temperature increases rapidly from 300 K to high values, while the ions are hardly heated. The reason is that the density is so low that the particles are solely heated by E ind . The duration of this stage strongly depends on the electric field, while only weakly on the pressure. (2) Stage II is a transition stage.
This stage is analogous to the end of the breakdown phase and the beginning of the burn-through phase. It is prominent only when the induced electric field is low at 0.3 V m −1 , otherwise the duration is too short to be recognizable. Because the effective E ambi is closely related to E ind , when E ind is small, E ambi increases slowly and needs much more time to exceed E ind . In this stage, the plasma density is low and still increases exponentially but slower. The electron temperature decreases rapidly at first, implying the electrons are cooled, or the heating rate is smaller than the cooling rate. The ion temperature increases slowly. In this stage, the E ambi is still small and effectively comparable to E ind . (3) Stage III is the slow avalanche stage.
This stage is analogous to the early hydrogen/deuterium burn-through phase, without influence from the impurities or later burn-through phase. In this stage, the plasma density is high, and still increases exponentially. The effective E ambi dominates over E ind . Electrons are heated by E ind and cooled by E ambi , while ions are heated by E ambi . The electrons drag ions to move in the same direction, so their currents are opposite. The electron temperature increases slowly or even decreases for some time. In our simulations, this stage will start whenever density increases to about 10 15 m −3 , regardless of the pressure and electric field, and end after the plasma is fully ionized. The duration of this stage strongly depends on both the electric field and the pressure. Higher pressure will lead to shorter duration of this stage.
Finally, it is worth discussing the applicable range of the model here. Particle losses and impurities are the dominant loss mechanisms in the later burn-through phase, therefore the model must include these processes to simulate the whole burn-through phase. To simulate the ramp-up phase, the magn etic field produced by the plasma current should also be considered, which will require an electromagnetic instead of electrostatic model. In future work, we will perfect our physics model step by step, considering the particle losses and impurities which are also important in the start-up process.