Molecular dynamics simulation of cathode crater formation in the cathode spot of vacuum arcs

A three-dimensional model based on molecular dynamics has been developed to describe the formation of a single cathode spot in vacuum arcs. The formation of the cathode spot is assumed to be controlled by the plasma ions, the effect of which is simulated in LAMMPS through the process of ion bombardment. The cathode is represented by structured copper atoms, while the ions are continuously injected into the domain with a certain velocity towards the cathode surface. Ion bombardment leads to the appearance of a crater, which is caused by the accumulation of pressure effect against the relaxation of substrate temperature. The size of the crater is found to be determined by the spatial distribution of the injected ions. The formation of the cathode spot is also scrutinised by electron emission from the cathode surface with variable surface temperature during the cathode spot development process. In addition, the evaporated atoms forming the metal vapour are observed. This study provides a description of the formation of the cathode spot at microscale, which shall be helpful to further studies of the arc rooting and arc contact (electrode) erosion in vacuum environment.


Introduction
The dielectric strength of vacuum gaps has become the focus of research and development activities in the global effort to * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. extend the application of vacuum switching technology into high and extra high voltage levels [1]. Arc induced contact erosion adversely affects the dielectric strength of a vacuum gap as a direct consequence of damage to the contact surface from macroscopic to microscopic scales. Arc contact erosion, especially for the cathode, is closely associated with the existence of the cathode spots and their dynamic behaviour [2,3]. Understanding of the formation dynamics of the cathode spots, which is a rather complex process, is vital for arc control, minimisation of contact damage, and improvement of dielectric performance in the development of high voltage vacuum circuit breakers.
A cathode spot is a small region on the cathode surface conducting a portion of the arc current [1]. Interacting with the plasma particles such as ions, the cathode surface is locally heated up and damaged. The high local surface temperature and electric field generate intense surface electron emission, which further ionizes the evaporated atoms [4]. When the surface temperature is insufficient to maintain surface emission of electrons, the cathode spot stops discharging, and a new cathode spot will need to be formed for current continuity.
Cathode spots have been the subject of experimental investigation over many decades. Different types and features of the cathode spots were proposed, such as Type-1 cathode spots appearing on contaminated surfaces and Type-2 cathode spots detected on clean surfaces [5][6][7][8][9][10]. The observed sizes of the cathode spots were basically in the range of micrometres [9]. The current that each cathode spot can carry depends on the contact material [11].
To explain the characteristics of cathode spots various theories have been proposed. In the early stage, the behaviours of ions, electrons, and atoms in and near the cathode spot were analysed and described by a surface layer model under an assumed stationary cathode surface [12,13]. Meanwhile, the transient process on cathode spot formation was explained by the 'ECTON' model, where a micro-explosion induced by intense Joule heating ignites a cathode spot [14,15]. This model was further developed in [16] to link the explosion to thermal runaway of a thin neck of current conduction formed under surface deformation. Thereafter, surface deformation under plasma influence becomes the key factor of this transient process of spot formation. Recent simulation of cathode spots commonly adopted the fluid assumption of the cathode to exhibit surface deformation by hydrodynamics (HD). For example, in [17][18][19], the bulk cathode was assumed as a fluid entity with distributed temperature-dependent viscosity. The influences of plasma on cathode surface are simplified into the boundary conditions including the energy flux density, current density, with the formation of surface crater observed as the result. Furthermore, in [20], the effects of self-generated plasma on the crater formation are also considered to reach a self-consistent model. However, the 'ECTON' theory, which is by far the most accepted explanation of the transient process of cathode spot, is still a hypothesis. The cathode spot simulations based on the HD method have not reached consensus in results, hence could not validate the 'ECTON' theory.
This paper presents a cathode spot study based on the molecular dynamics (MD) method. It is based on fundamental principles, i.e. the activity of individual atoms within a system of atoms according to Newton's law. The MD method has been successfully used to simulate the dynamic change of a metallic structure with a size in nanometres such as the deformation of a nanotip [21]. Also [22], firstly adopted the MD method to simulate the surface modification in the plasma-surface interactions, presenting the craterlike damage on the surface by the fluxes of energetic ions. Compared to the HD method where the cathode body is assumed as a masscontinuous fluid, the MD method regards the cathode body as discrete atoms. Therefore, the detailed physical processes concerning the cathode spot can be analysed by modelling the microscopic behaviour of individual particles.
In this work, the formation of a single cathode spot is simulated in large-scale atomic/molecular massively parallel simulator (LAMMPS), aiming to analyse the surface deformation and the evolution of surface temperature happened on a cathode spot [23]. The plasma effects are represented by ion bombardment, and the surface morphology is reflected by the evolution of atomic structure. First, the mechanism of surface deformation under ion bombardment is modelled and analysed, leading to the conclusion that the size of the cathode crater is directly related to the spatial distribution of the plasma ion moving towards the cathode. Accordingly, a scaled-down cathode spot model in nano scale is adopted under the limitation of the spatial scale in the MD method, much smaller than the size of the practical cathode spots (micrometres). The thermal effects of plasma ion heating, thermal conduction, and material evaporation are all considered, without coupling the Nottingham heat nor the influence of ionization of the metal vapour. As a result, the formation of the cathode spot is assessed by the surface electron emission derived from the surface temperature. This research helps understand the molecular process of cathode spot formation and is promising in providing useful information needed for future vacuum arc research.

Preparation for simulation
The formation of the cathode spot is the result of plasmasurface interaction. Particles that are involved in this interaction include plasma ions, plasma electrons, surface emitted electrons, surface evaporated atoms, secondary emitted electrons, back ions that are derived from the evaporated atoms and return to the cathode surface, etc [24,25]. According to the results of previous cathode spot simulations in HD [17][18][19][20], among all the species, the ions from leftover plasma dominate the contribution (over 80%) to the energy flux density and surface pressure during the formation of the cathode spot [18]. Additionally, the contribution of Joule heating was proved considerably lower than the other effects [20]. Therefore, the assumption adopted in this paper is that the cathode spot is formed under the stimuli of the leftover plasma ions, while the contribution of other species, Joule heating, and radiation are ignored. In the vacuum arc, the plasma ions are derived from the ionisation of metal vapour evaporated from the electrode surface. In this work, the cathode material is set as pure copper. Therefore, the type of plasma ions should be copper ions.
According to the surface layer model [12,13,26], all the plasma ions travelling through the sheath fall can be assumed with the same energy q ip : where the first term represents the Bohm velocity when ions are passing the interface of pre-sheath and sheath, the second term represents the enthalpy of each ion when approaching the cathode surface, and the third term represents the kinetic energy gained during the acceleration by the electric field in the sheath. z i is the average charge number of the ions, U stands for the cathode sheath voltage, k is the Boltzmann factor, T e represents the electron temperature in plasma, while T i is the surface temperature. Right before the incident copper ion bombard the surface, it is Auger-neutralised to be a copper atom by an electron excited from the surface [27]. Since the plasma ions bombard the cathode bulk after being neutralised, the incident particles can be perceived in the form of atoms. The energy q a of the incident atom is calculated as: where V i represents the total ionisation energy of the incident ions, and φ eff represents the effective work function required to extract an electron from the surface for the neutralisation. The condensation energy of the atom is not considered as a part of this external energy as it is involved internal the process of ion bombardment in MD.
Since the drift velocity obtained from sheath acceleration is much greater than the thermal velocity of ions, the directional randomness of the thermal velocity is ignored. Therefore, the direction of the plasma ion current is assumed normal to the surface. As the ion current density J ip is uniform in the area of current path S, the plasma ions current I ip can be expressed as: where n i stands for the volume density of plasma ions, and v i represents the velocity of plasma ions. Further, the equivalent interval of every ion bombardment on the cathode surface T ib in the current path is determined by the density and velocity of ions as: With ions velocity fixed by the calculation of ions energy in (1) and (2), a greater ions density leads to a shorter interval between each ion bombardment, which means higher ion bombardment frequency.

Modelling of cathode spot
Following the stated assumptions, the formation of the cathode spot is realized by the copper substrate bombarded by copper atoms. As shown in figure 1, the model is composed of a copper bulk and a headroom for the distributed incident atoms to be injected continuously. The thermostat boundary is set for ±x, ±y, and −z direction by having multiple atom layers maintained around 300 K. For z direction non-periodic boundary is set to allow evaporated atoms to leave this box. Therefore, several fixed atom layers at the bottom are necessary to prevent the atoms leaving the domain from the bottom. In this model, the interactions between the copper atoms are governed by the embedded atom method potential describing the stable metallic structure, combined with the Ziegler-Biersacke-Littmard repulsive potential when the interatomic distance is too short. This potential has been justified and verified in the simulation of the ion bombardment in radiation [28]. The potential energy is evaluated for the atoms with a distance between 0.02 and 6 Å. The crystallographic orientations of the copper bulk are . The standard Velocity-Verlet algorithm is used to integrate Newton's equations of motion.
According to [1,9], the effect of plasma ions on the cathode is determined by the sheath above the cathode surface instead of the whole vacuum gap. The kinetic velocity gained by the ion during the sheath acceleration is much higher than their thermal velocity in the vacuum gap. Considering that the thickness of the sheath is several nanometres [9], and the calculated velocity of the ion is in the degree of 10 nm ps −1 , in this work the timestep is set as 1 fs, and the duration of simulation is set as hundreds of picoseconds.
This simulation commences with a relaxation period of 10 ps for the copper substrate to be fully relaxed at 300 K under canonical ensemble. After this relaxation, the period of ion bombardment begins. The incident ions are injected in the headroom continuously, with the incident velocity of each ion and the interval of ion injection controlled as the input, while the injection position is random. From the discussion in section 2.1, these two parameters v i and T ib can be derived from the basic parameters describing the plasma ions environment, including z i , U, T e , and n i . Their values have been measured from experiments and are shown varied over a wide range [9,15]. In this work, T e = 2 eV, z i = 2, U = 15 V have been assumed. The corresponding ion energy is 62 eV, and the velocity of the ions is 137 Å ps −1 . The peak ions density n i0 = 1 × 10 26 m −3 is selected as the reference, with the distribution of ions density discussed later. During this period the whole system is under microcanonical ensemble with the total energy of the system updated after new plasma ions insertion.

Mechanism of cathode crater formation under ion bombardment
In previous cathode spot simulations, the cathode spot is observed by the formation of the cathode crater, which is the result of surface deformation [17][18][19][20]. In this section, different strategies of ion bombardment are compared to analyse the effect of the two input parameters v i and T ib on the surface deformation. Afterwards, the mechanism of crater formation under ion bombardment is concluded.

Single ion bombardment: the effect of v i
The formation of a cathode crater and the generation of atom sputtering usually appear after ion bombardment, as analysed in the case of irradiation [29,30]. However, the ion energy in the applications of irradiation is generally 500 eV or even higher, much greater than the ion energy in the case of cathode spot that is 62 eV. Therefore, different phenomena are expected to be induced by the different ion energies. Two sets of single ion bombardment simulations are compared here with the ion energy 62 eV and 500 eV, respectively. The model shown in figure 1 is adopted with a size of 7.23 × 7.23 × 7.23 nm. The initially relaxed substrate temperature is 300 K. Each simulation set is repeated 100 times, and the typical result is presented as follows.
In the case of 500 eV, the incident ion may be injected into a deep position of the substrate. Correspondingly, a thin but deep crater is formed temporarily, and a sputtered atom is observed. While in the case of 62 eV, the incident ion always stays at the surface of the substrate, without the crater formed nor the atom sputtered out. Figure 2 shows the evolution of the position and velocity of the incident ion during the bombardment process. The motion of the sputtered atom in the case of 500 eV is also attached.
The surface deformation of ion bombardment lies in the process of the ion moving deeper and squeezing the substrate atoms away. Therefore, this pressure effect can be reflected by the value of injection depth, which is defined as the maximum depth the incident ion reaches as the injection depth. In the case of 62 eV, the injection depth is always nearly zero, which means that the ion energy is too low to produce an obvious pressure effect on a cold cathode.

One-site continuous ion bombardment: the effect of T ib
Since no surface deformation is detected after a single ion bombardment on the cold cathode, the cathode crater formation should owe to the temperature rise in the substrate. In this part, the simulation sets with various ion insertion intervals T ib are compared. The incident ions of 62 eV are continuously inserted at the same position, with the T ib set as 200 fs, 400 fs, 800 fs, and 2000 fs, respectively. After the continuous ion bombardment for 100 ps, the cathode crater is observed, which is bigger in the case with smaller T ib . In the last set, the cathode crater is absent. The evolution of the injection depth of the first ten incident ions is shown in figure 3(a). The curves of the first three simulation sets commonly exhibit a trend that the next incident ion has a bigger injection depth, and this change is more obvious under shorter T ib . On the contrary, in the case of 2000 fs, the incident ions always stay at the surface. By this comparison it can be seen that the formation of the cathode crater is related to the increase of the injection depth.
The corresponding evolutions of the local temperature of the bombardment site are compared. Here the local temperature is calculated by the average thermal velocity of the atoms in a small region representing the bombardment site. Figure 3(b) presents the comparison between 400 fs and 800 fs for example. After each ion bombardment, the bombardment site goes through a transient period first with the local temperature rising sharply, then a relaxation period with the local temperature decreasing by the effect of heat conduction. Energy dissipation through heat conduction is determined by the temperature gradient. Therefore, during the relaxation period, the temperature decreases very fast in the beginning and gradually gets slower. From figure 3(b), with the ion energy similar, a longer relaxation period leads to lower local temperature when the next ion bombardment happens. Name this temperature as the substrate temperature.
The evolutions of substrate temperature during the first ten ion bombardments are shown in figure 3(c). This substrate temperature grows faster with shorter bombardment intervals. Higher bombardment frequency uprises the substrate temperature faster, thus generating a more prominent pressure effect and producing a severe cathode crater. In addition, the curve  of 200 fs shows that this stable substrate temperature is limited by the boiling point of substrate material, which is 2650 K for copper in this model. After reaching the boiling point the substrate temperature starts to fluctuate. Further, the curve of 2000 fs shows that even though the incident ion always stays at the surface, the energy input still heats the substrate continuously. Although no cathode crater is observed within 100 ps, given a sufficiently long time of ion bombardment, the substrate temperature will be high enough, with the atomic structure of the bombardment site less compact. As a result, an obvious pressure effect can be generated.
In conclusion, the pressure effect of ion bombardment depends on the ion energy and bombardment frequency simultaneously. The mechanism of cathode crater formation can be concluded as the accumulation of pressure effect against the relaxation of substrate temperature.

Regional continuous ion bombardment: the effect of the ion density distribution
In the previous cathode spot simulation without referring to the current transfer, an assumed spatial distribution of ions density is commonly assumed [17][18][19][20]. In order to connect the effects of v i and T ib with the ion density distribution, in this part, the simulation of regional continuous ion bombardment is analysed.
According to section 2.1, the relationships between the parameters including ion current, ion current path, ion current density, ion energy, and bombardment frequency are analysed. With the ion energy at always 62 eV, four simulation sets are conducted as given in table 1. In these four sets the ion current density is assumed even distributed in the bombardment area. The simulation time for continuous ion bombardment is also 100 ps. The picture of the crater will be shown later, and the crater radius given in table 1 is the radius of the hole observed in the plane of z = 0.
In simulation set 1, 2, and 3, the bombardment area is the same while the ion current density is different. The result is similar to the comparison in section 3.2, as the ion current density, ion density, and bombardment frequency are directly related when the ion energy is fixed. In the case with higher current density, the crater is bigger and deeper. In simulation set 3, no crater is observed after 100 ps. Because the bombardment frequency decides the rate-of-rise of substrate temperature, from the perspective of cathode crater formation, the ion density decides the speed of cathode crater formation. During a certain input window, when the ion density is too low, the rise of substrate temperature will be insufficient to produce an obvious pressure effect. Therefore, a threshold of ion density can be perceived to exist for the formation of the cathode crater. The precise value of this threshold is hard to obtain but is roughly in the range of 1 × 10 10 ∼ 1 × 10 11 A m −2 for the input window of 100 ps and ion energy of 62 eV.
In simulation set 1, 4, and 5, the total ion current is the same, while the bombardment area and ion current density are changed reversely. As a result, the cathode craters with different size features are observed. According to the information of craters recorded at 30 ps, the case with greater ion current density has deeper crater depth, which is consistent with the previous analysis. However, the crater radius, which is already stable at this timestep, is shown nearly consistent with the size of the bombardment area. Therefore, the case with higher ion current density and smaller bombardment area has a small and deep crater, while the case with lower ion current density and bigger bombardment area shows a big and shallow crater.
From this comparison, with the same ion energy, the value of ion density is found critical to the intensity of cathode spot formation, while the stable radius of the crater is determined by the area of the ion density distribution [31]. The various features of the cathode crater observed in experiments may be explained by the coordination of these two parameters.
An assumption of the spatial distribution f (r) of plasma ions was proposed by [32] and commonly adopted by [17][18][19]: where n ip is the ion density, n i0 is the peak value of the ion density, and r is the radial position on the surface. The point of r 0 is the position where the ion density is half of the peak value, while k decides the rate of the decrease of n ip radially outwards.
Applying the findings in the simulation with uniform ion distribution to the case with this type of distribution, it can be concluded that the value of n i0 determines the intensity of cathode crater formation. Under a given input window, when n i0 is smaller than the threshold of ion density required, the cathode crater will not appear. Besides, the position where ion density equals to the threshold value is mutually decided by r 0 and k. With k usually set as 2, the size of the cathode crater is nearly proportional to the value of r 0 .
In conclusion, under the current simplified simulation scheme that the cathode spot is formed under the stimuli of plasma ions distribution, the developing speed of the cathode crater is decided by the ion current density, and the stable radius of the cathode crater is solely decided by the spatial distribution of plasma ions.

Simulation result of the cathode spot formation
The computational burden required in MD simulation is determined by the number of atoms involved, which limits the spatial size of the MD model. Compared to the physical size of the cathode spot in the degree of micrometres, the spatial size of the MD model built in LAMMPS is only in the degree of nanometres. However, according to section 3, through the ion bombardment, the size feature of the cathode crater is controlled by the size feature of ion distribution. Therefore, a scaled-down model can be adopted to present and analyse the process of cathode spot formation.
In this section, an example of cathode spot formation is presented. A model as shown in figure 1 with the size of copper substrate 20 × 20 × 10 nm is adopted. The spatial distribution of plasma ions distribution is assumed as (5). The ion energy and peak ion density are still 62 eV and 1 × 10 26 m −3 according to the model described in section 2, with r 0 adjusted to 0.6 nm for the equivalence. To reflect the ion distribution the bombardment region is divided into several ring-shaped areas with the same width. For each ring different interval of ion insertion T ib is set correspondingly. The continuous ion bombardment simulation for 800 ps is conducted. As a result, the evolution of the structure deformation and temperature distribution in the substrate is first presented. Then, the deformed surface is tracked, the temperature of which is obtained for the calculation of surface electron emission. According to the distribution of electron current density, the formation of the cathode spot is assessed.

Evolution of bulk deformation and bulk temperature
By OVITO, which is a post-processing software for MD simulation, the evolution of atomic structure and atom temperature can both be observed. The atoms in the copper substrate keep vibrating, because of which the atom temperature or the atom velocity exhibits a kind of randomness. For better observation of the temperature, the copper substrate is meshed into small chunks with a size of 0.732 × 0.732 × 0.732 nm, which is two times the lattice constant of copper. The mean velocity of all the atoms in each chunk is calculated and regarded as the temperature of this chunk.
An example of the evolution of bulk temperature and bulk deformation is shown in figure 4. The colour bar in the pictures of the upper half represents the chunk temperature, with the corresponding surface mesh at the lower half presenting the deformation of the substrate. The cross-section of the substrate is shown while the boundary is not. In the beginning, the temperature of the central region increases first by the most frequent ion bombardment. The heat conduction outward and the less frequent ion bombardment in neighbour regions lead to the slower temperature increase of neighbour regions. Theoretically, the volume of the region with higher temperature should swell due to the larger inter-atomic distance. However, the pressure effect of incident ions leads to a sunken area in the centre. Consequently, the sunken centre and swelled vicinity together exhibit the central crater. Along with the heat conduction, the incident energy is spread further, with the crater expanding radially. Such crater growth is obvious during 0 ∼ 400 ps, as the comparison between figures 4(a)-(d). However, the difference between figures 4(d) and (e) is small, showing that the growth of cathode crater is already stable.
As for the assessment of the molten pool, the liquid area and solid area can be easily distinguished by the local lattice structure. The solid area is defined by the regions which remain the stable FCC metallic structure, while the rest are perceived as liquid. The size of the molten pool is plotted in figure 5, showing that the molten fool is quickly formed and expanded in the beginning but balanced afterwards. This trend is similar to the growth of the crater, both reflecting a thermal balance is already built after 400 ps.

Evolution of surface deformation and surface temperature
Surface temperature is an important parameter in the cathode spot theory for calculating surface emission. Also, surface deformation is the focus in recent simulations to involve its influence on the plasma-surface interactions. Therefore, the calculation of the surface state is essential to discuss the formation of a cathode spot. The method for tracking the deformed surface in LAMMPS is explained as follows.
Firstly, the evaporated atoms are distinguished by coordination analysis and not counted in the real-time substrate. Secondly, the substrate atoms are meshed in x and y directions and the highest position of each x-y chunk is determined and recorded. Thirdly, the atoms whose axial positions lie in a corresponding height range are grouped as the surface atoms. Finally, the surface atoms are meshed again in the x and y directions, with the centre-of-mass in the z-direction and the mean temperature of the chunks calculated.
The atomic structure and the temperature distribution of the tracked surface are shown in figure 6. As shown the deformed surface is successfully tracked dynamically in this model, and the evolution of the surface is representative of the evolution of the copper substrate. The centre of the surface is quickly melted and deformed when ion bombardment begins. After around 400 ∼ 500 ps the balance is reached. The maximum surface temperature observed during this simulation is around 4500 K. Worth mentioning that during the single ion bombardment, some substrate atoms will possess ultra-high drift velocity and transient temperature for several timesteps. The temperature presented here tried to avoid the influence of such a transient phenomenon and reach a steady-state temperature.
From figure 6 the surface temperature evolution is also observed to be nearly axisymmetric. Therefore, the surface chunks can be grouped according to their radial position to present the radial distribution of surface temperature. Here the surface chunks are divided into 16 radial layers with the width same as the chunk size. Further, the temperature of each chunk is time-averaged within 10 ps to weaken the influence of the randomness of atom temperature. Figure 7(a) presents the radial distribution of temperature with respect to time, and figure 7(b) shows the time evolution of temperature in each layer. From figure 5 the central region within the radius of 2 nm is quickly heated up over the melting point and even surpasses the boiling point in the first several tens of picoseconds. Then after 100 ps the temperature of the central layer vibrates around 3500 K, with the temperature of outer layers increasing slowly. There is nearly no difference between the curve of 400 ps and that of 800 ps. Therefore, a stable temperature gradient is established after 400 ps. At the balanced state, the inner layer possesses a higher temperature. As well, the inner layer reaches its balance temperature faster.

Surface emission
According to the T-F emission theory [3], surface temperature decides the energy of free electrons, and the surface electric field influences on the tunnel effect. With surface temperature and electric field intensity known, the current density of emitted electrons j em can be calculated as:  where λ R and A G are both constant coefficients for a given type of metal. T i represents the surface temperature, and W and ∆W are the work function and Schottky correction, respectively. The Schottky correction is calculated by: exhibiting the influence of surface electric field. The surface electric field is dependent on the motion of incident ions in this model and is calculated as: where j ip represents the incident ion current density and ε 0 is the vacuum permittivity. The energy flux density of the cooling effect by surface emission is: As the ion current density is set as constant in the model, the surface emission is determined by the surface temperature. Figure 8(a) shows the electron current density under different surface temperatures. Compared with the current density of plasma ion which is in the range of 3 ∼ 6 × 10 10 Am −2 in the central area, the electron current density is only comparable when the surface temperature is above 3000 K. According to the radial distribution of the surface temperature as shown in figure 7, only the central region with the radius of 1.5 nm meets the demand. Consequently, in this model, the surface emission is only intense enough in this central region. The electron current density in this region is calculated based on the surface temperature obtained from the simulation and its evolution is presented in figure 8(b) with the error bar. The peak current density of the plasma ions is also shown for comparison. From the perspective of current continuity, the transition from it being maintained by plasma ion to it being conserved by the surface emitted electron is realised at 100 ps. Therefore, it can be perceived that a stable cathode spot is established in this central region.
Correspondingly, the energy flux densities of surface emission under different surface temperatures are plotted in figure 8(c). This picture shows that this cooling effect is more prominent in regions with higher temperatures. Under the simulation conditions, the energy flux density of electron emission is only 8.6 × 10 10 Wm −2 at 3500 K, which is smaller than the peak energy flux density of incident ions of 1.2 × 10 12 Wm −2 , but higher than the cooling effect of evaporation as calculated above. Therefore, conclusions can be drawn that the cooling effect of surface emission is most intense in the central region, where its contribution to the thermal balance is non-negligible. The cooling process of surface emission needs to be incorporated in the future model.

Evaporation
In this model, atom evaporation from the surface is observed as shown in figure 1. The distribution of the evaporated atom numbers within 400 ps is shown in figure 9(a). The phenomenon that the most evaporated atoms are from the central region is reasonable because this region possesses the high temperature above the boiling point. Figure 9(b) presents the flux density of evaporation in the central square area of 8 × 8 nm, which fluctuates around 1 × 10 28 m −2 s −1 against time. The corresponding energy flux density of this cooling effect is around 5.12 × 10 9 Wm −2 . Theoretically, the evaporation flux density Γ vap is calculated as: where p v represents the saturated vapour pressure, and T i stands for the surface temperature [12]. m i means the ion mass, and k is the Boltzmann constant. The relationship between the theoretical evaporation flux and the surface temperature is plotted in figure 9(c). According to the surface temperature distribution with the thermostat boundary presented in figure 7, the central temperature is 3500 ∼ 4000 K. The corresponding theoretical evaporation flux density is nearly 1 × 10 28 m −2 s −1 as shown in the circle, which is comparable with the simulated evaporation flux density in the central region of 8 × 8 nm. This comparison validates the evaporation results of this simulation.

Comparison with experiments and other simulations
The observation of transient phenomena in an individual cathode spot is limited by the existing experimental techniques. Therefore, the comparison between cathode spot simulations and experiments is represented by the evolution of the cathode crater. Figure 10 presents an example of the surface craters formed in our simulation, in a HD simulation [18], and observed in the experiments [8][9][10]. In figures 10(a) and (c), the cathode craters, only from the view of morphologies, are all composed of a sunken crater in the centre and a rim surrounding the crater. The spatial scale in figure 10(a) is several nanometres, while the spatial scale in the other two are several micrometres. The similar leftover plasma ion input with different spatial scales leads to similar crater profiles as expected. The dependence of the spot size on the leftover plasma ion input has been discussed in section 3. In addition, the observed radius of the cathode crater in experiments varied in a wide range from hundreds of nanometres to tens of micrometres [9], which might be explained by our finding that the size of the cathode spot is highly dependent on the density distribution of plasma ions.
Comparing the results between this model and other simulation works, the evolutions in morphology of the obtained surface crater are also similar, with the similar assumptions of the input leftover plasma ions, e.g. figures 10(a) and (b). In addition, the similarity is reflected by several key parameters obtained from the simulation results. In this paper and other simulation works, the peak surface temperatures of the Cu cathodes are all around 4000 K, and the current density of surface emitted electron are in the same order of 10 11 Am −2 [17][18][19][20].
Compared to other simulation methods, MD method has its uniqueness and advantages. First, MD method allows the observation of the deformation of the discrete atomic structure, instead of assuming the surface as an entity with masscontinuity. Therefore, this model provides a way close-toreality to analyse the surface deformation. Second, the effect of the leftover plasma ions is directly realised in this model through ion bombardment, without the requirement to isolate its thermal effect and pressure effect on the cathode surface, which is intrinsically coupled in MD simulation. Third, the density distribution of evaporated atoms can be observed directly through MD method, which is not available in HD method.
Compared to the many physical processes involved in a cathode spot, this simulation work has certain limitation, as only focusing on leftover plasma ions. The other factors such as the Joule heating, the back ions as well as the thermal effect of surface emission are not taken into consideration. As explained previously in section 1, these effects do exist in the cathode spot processes, however, they are rather secondary in the formation of cathode spot including the crater formation and the evolution of surface temperature. Nevertheless, this work applied the MD method to analyse crater formation of cathode spot and its follow-up work intends to eliminate these limitations-and reach a more self-consistent model by incorporating other physical processes.

Summary
A 3D model in nano scale has been developed based on the MD method to simulate the formation and development of a single cathode spot in vacuum environment. The formation of the cathode spot is represented by the structure deformation and evolution of surface temperature distribution. The mechanism of cathode crater formation is found to be the accumulation of pressure effect against the relaxation of local surface temperature. A threshold of plasma ion density is identified for the formation of the cathode crater within a certain effective duration of the leftover plasma. Also in this simulation scheme, the size of the cathode spot is decided by the spatial distribution of plasma ions flux.
The thermal processes during the formation of cathode spots are discussed, with contributions of thermal conduction and evaporation considered. The central region is quickly heated up by intense ion bombardment and stabilises around 3500-4000 K, sufficient to generate surface emission electrons with a current density higher than that of the plasma ions. By the thermal conduction and less frequent ion bombardment, the outer region is slowly heated and reaches lower temperatures. Eventually, a stable temperature gradient is built up, and the cathode crater and the molten pool stop growing. When the input of leftover plasma is controlled within, a cathode crater with a similar size is obtained, and a cathode spot with intense surface electron emission over 6 × 10 10 Am −2 is detected in the central region of around 30% of the input radius. The relationships among the sizes and the current densities are consistent with previous simulation studies.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).