Dynamic coupling between particle-in-cell and atomistic simulations

We propose a method to directly couple molecular dynamics, finite element method and particle-in-cell techniques to simulate metal surface response to high electric fields. We use this method to simulate the evolution of a field emitting tip under thermal runaway by fully including the 3D space-charge effects. We also present a comparison of the runaway process between the two tip geometries of different widths. The results show with high statistical significance that in case of sufficiently narrow field emitters, the thermal runaway occurs in cycles where intensive neutral evaporation alternates with the cooling periods. The comparison with previous works shows, that the evaporation rate in the regime of intensive evaporation is sufficient to ignite a plasma arc above the simulated field emitters.


I. INTRODUCTION
Field emitting tips play a detrimental role in various vacuum devices that require high electric fields, such as vacuum interrupters [1], X-ray sources [2], fusion reactors [3] and particle accelerators [4]. When a field emitter is subjected to currents exceeding a certain threshold, the emission becomes unstable and swiftly results in a vacuum arc [5,6]. The latter is characterized by the ignition of a plasma in the vacuum gap, which drives a large current and converts the gap into a short-circuit. The ignition of an arc in vacuum (vacuum breakdown) can cause catastrophic failure in electron sources [5,7].
Although vacuum arcs have been studied more than half a century and it is well-known that they appear after intense field electron emission [6,8], the physical mechanisms that lead from field emission to plasma ignition without the presence of any ionizable gas are not understood yet. One hypothesis commonly used to explain the plasma build-up in vacuum after intense field emission is the "explosive emission" scenario [8][9][10]. According to it, when intensive field emission takes place, there is a critical current density beyond which heating is produced at a rate the emitter cannot dissipate, which leads to heat accumulation, sufficient to cause instant explosion and plasma formation. This phenomenological description does not give an insightful physical understanding of the underlying processes and cannot a priori * mihkel.veske@helsinki.fi predict the behavior of a specific system. Such empirical considerations have been common in the vacuum arcing theory since the 1950's, as the modern computational tools necessary to describe the physical processes of such a complex phenomenon as vacuum arcing were not available. Thus, there is a growing need for the development of theoretical concepts and computational tools that describe the vacuum arc ignition in a rigorous quantitative manner and provide a deep understanding in the physical processes that lead to it.
This has changed recently, since various computational methods have been employed to study vacuum arcs. In a recent work [11], plasma formation in Cu cathodes was studied by means of particle-in-cell (PIC) simulations. It was found that plasma can build up in the vicinity of an intensively field emitting cathode under certain assumptions regarding the interactions between the cathode surface and the plasma. One of them is that the cathode, apart from electrons, emits also neutral Cu atoms at a minimum rate of 0.015 neutrals per electron. However, the origin of the neutral atoms required to initiate a selfsustainable plasma is not yet clear.
Recent multiscale atomistic simulations on Cu nanotips [12] revealed a thermal runaway process which causes evaporation of Cu atoms and nanoclusters at a rate exceeding the above minimum required to ignite plasma. Thus, it gives a possible mechanism to explain the supply of neutral atoms required for the ignition of plasma. Nevertheless, this process is only indirectly connected to the plasma simulations via the mean evaporation rate. The preliminary results of [12] do not form a solid proof that this is the mechanism that initiates a arXiv:1906.08125v1 [cs.CE] 19 Jun 2019 vacuum arc. Furthermore, the interactions between the evaporated material, the metal surface and the emitted electrons that lead to plasma expansion cannot be understood without a direct coupling between the metal surface and plasma simulations.
In order to decipher the plasma ignition mechanism, it is necessary to consider the coupling between the processes unraveling on the metal surface under high electric field and the ones taking place in vacuum that lead to plasma formation. Until now, these two domains have been simulated independently. The evolution of the metal surface under high electric field has been simulated using hybrid molecular dynamics (MD), electrodynamics, electron emission and heat diffusion calculations [12][13][14][15]. On the other hand, the evolution of the vaporized and ionized material that leads to plasma ignition in the vacuum has been simulated with the PIC [11].
In this paper we develop a method that directly couples PIC simulations with the finite element method (FEM) and MD-based calculations of the metal surface response to the field. We use that method to calculate the evolution of a field emitting tip under thermal runaway. The utilization of the PIC method allows for a full 3D calculation of the space-charge effects, without imposing the semi-empirical corrected equivalent planar diode model used previously [12]. We also present a comparison of the runaway process between two tip geometries of different widths. That comparison takes us closer to obtaining an estimation for the runaway occurrence probability in realistic structures observed in experiments [16]. The incorporation of the PIC method into our simulations forms a decisive step towards the direct coupling between atomistic and plasma simulations, in order to fully reveal the processes that take place during the initiation of a vacuum arc.

II. METHOD
The model, which we present here, combines classical MD, FEM, PIC and electron emission calculations. The model allows to assess the effect of different physical processes on the dynamic evolution of the atomic structure of metal surfaces exposed to high electric fields. While the movement of the atoms is followed by an MD algorithm, the movement of the emitted electrons is tracked using the PIC method. FEM provides the tools for concurrent and self-consistent electric field and heat diffusion calculations, whose output affects the force and velocity calculations in MD and the force exerted on the electrons in PIC. Such a methodology forms the conceptual basis for the full coupling between metal surface and plasma simulations. A full schematic of the multiscale model is provided in appendix A.
The model is based upon the previously developed framework FEMOCS [15]. In FEMOCS, an unstructured flexible mesh is built, the density of which is controlled independently in all regions of the simulation domain. The mesh consists of linear hexahedral elements that are built into tetrahedra that allow achieving satisfactory accuracy while keeping the computational cost within reasonable limits. In the following sections we provide an overview of the methodology that was added to FEMOCS as compared to one reported earlier [15].

A. Electric field
A common way to take into account an electric field E within a computational model is to define it as the negative gradient of the electrostatic potential Φ that, in turn, can be calculated by solving the Poisson's equation where is the dielectric constant and ρ is the space-charge density. We solve eq. (1) iteratively with the boundary conditions (BC) defined in the vacuum (see fig. 5) as where E 0 is the long-range applied electric field. The boundary value problem (1)-(2) is solved in FEMOCS by means of FEM. The derivation of weak and algebraic formulation of eq. (1)-(2) is provided in appendix B.
To obtain the space-charge distribution ρ( r) in the vacuum above the metal surface, which enters eq. (1), we use the PIC method [17][18][19]. Within PIC, the spacecharge is formed by a set of particles that interact with each-other and with the electromagnetic field. Typically it is not feasible to track every real electron and ion in the system. As the interactions depend on the chargeto-mass-ratio instead of the charge or the mass of the particles, it is reasonable to group charges of the same type into superparticles (SP) that act within the model as single particles. Thus, the SPs represent a sampling of the continuous phase space occupied by the real particles like electrons or ions. The size of each sample (SP) is determined by its weight w sp , which corresponds to the (fractional) number of particles it represents. The SP weight is determined independently for each simulation, in order to achieve a sufficiently smooth charge distribution with a feasible number of SPs.
In general, the interactions of different types of particles in a system can be considered within PIC. Our current model, however, includes only the electrons, ignoring the presence of ions as well as electron-neutral interactions. Within this approximation, we can study the processes at the early stage of plasma formation, which is the focus of the current work. However, the model can be extended further if later stages of plasma evolution are of interest. In the following four sections, we give an overview of the various stages of the PIC model, that are summarized in fig. 1.

Particle mover
The propagation of SPs is performed by solving iteratively Newton's equations of motion for each particle. The acceleration of each SP by the electromagnetic field is found as˙ where B is the magnetic field at the location of the SP and q m is its charge-to-mass-ratio. Generally speaking, this equation can be solved numerically by means of the Boris method [20]. In our simulations, however, we apply only electrostatic fields and, since the magnetic field which is induced by the currents within the tip is negligible, we assume B = 0. With this approximation, we can employ a more efficient scheme by Verlet [21]. In this method, the positions and the velocities of SPs are calculated as v k+0. 5 where the subscript indicates the time step, Q = q m ∆t 2 and ∆t is the PIC time step. Note, that the calculation of v k+1 requires the value of the electric field at the time step k + 1. This is obtained in the field solver between steps (4b) and (4c), using as input the positions r k+1 . Since we have used periodic BCs for the electric field in the lateral directions, we apply the same type of symmetry for the SPs. This means that an SP that crosses the side boundaries Γ 2 , will reappear at the symmetric point at the opposite boundary of the box. On the other hand, particles crossing the top boundary Γ 1 and the material boundary Γ 3 (the latter happens very rarely for electrons) are removed from the simulation. Note that the top boundary with a Neumann BC does not represent a physical anode electrode, but rather the uniform far field behavior of the gap between two macroscopically flat electrodes. Thus the box is chosen in our simulation to be sufficiently high so that the electric field becomes uniform and the space charge negligible, which motivates the SP removal at that side.
On the cathode boundary, in practice, the absorbed particles may cause secondary effects, such as secondary and higher order emission, plasma recycling, impurity sputtering etc. Since these effect are beyond the scope of the current model and almost all electrons are removed at the top boundary, we leave their implementation for future work. Some attempts to take into account the above-mentioned effects can be found elsewhere [22][23][24][25][26].

Particle weighting
An important part of any PIC simulation, both conceptually and computationally, is mapping of the electric field to the position of SPs. In our model, we first locate the finite element that surrounds the particle and then use the shape functions of that element to interpolate the electric field at the position of the SP. The search for the element is carried out in a similar manner as during the particle injection (see appendix E), but instead of the three barycentric coordinates, here we need four, i.e. one for each tetrahedral vertex.
To solve eq. (1) we have to perform an inverse operation since the SP charge is required in the right-hand-side (RHS) of the matrix equation (B7). This, however, cannot be done immediately, as the SP in PIC is assigned a discrete charge q, while in eq. (1) the charge appears in form of charge density ρ. To overcome this issue, we represent the charge of each SP by means of the delta function, so that where δ( r − r i ) is the Dirac delta function and the summation goes over the SPs that are located in the finite element where ρ is evaluated. By doing such a substitution, the RHS of eq. (B7) becomes where r j are the coordinates of the jth SP that is located in the element where f i is computed. It is important to note that to prevent self-forces, the same shape functions N i are used as during the mapping of the field to the position of the SPs and during the distribution of the charge inside an element.

Particle injection
In our simulations, the space-charge is built up due to intensive electron emission. Using the field emission tool GETELEC [14], we calculate the current density J e in the centroid of each quadrangle that is located at the metalvacuum boundary. J e determines the number of electron SPs, n sp , that will be injected from a given quadrangle with area A at given time step: where e is the elementary charge. In general, n sp is a real number that has integer and fractional parts. Since the number of injected SPs can be only integer, we use a uniformly distributed random number R ∈ [0, 1] to decide whether to round n sp up or down: The emitted SPs are distributed uniformly on the quadrangle surface. More details about it can be found in appendix E. After injection, the electric field E = E( r) is used to give them an initial velocity and displacement The random number R provides a fractional push by a random fraction of the previous time step, preventing the particles to form artificial bunches [11].

Electron superparticle collisions
In addition to interaction with the field, every electron experiences the effect of other charges present in the system. A naive approach to take into account the Coulomb forces would be to calculate exactly the interaction between each pair of particles within the Debye sphere. This, however, would result in O(n 2 ) complexity and is therefore not suitable for the large systems aimed in the present work. Instead, we applied a more reasonable Monte Carlo binary collision model [27], which can be described in the following steps: 1. Divide SPs into groups in such a way that all particles within the group are located in the same cell.
2. Pair the SPs in the group in a random way.
3. Collide each pair elastically by assuming a normal distribution of the scattering angles.
The pairing of the SPs is implemented here as described in [27]. More details on the statistical model used to describe the collisions of SPs, the last step in applying the binary collision model, is provided in appendix F.

B. Electrostatic forces
The charge induced by an electric field on the metal surface affects the dynamics of the atomic system due to additional electostatic interactions. To assess this effect, the value of the induced charge on the atoms must be estimated and introduced to the MD simulations where electrostatic interaction between the atoms are included as a force perturbation [13,15]. The main perturbation is due to the Lorentz force while a less significant contribution comes from the Coulomb interactions between the surface charges In the above equations, q i is the charge of ith atom, r ij the distance between the ith and jth atoms andr ij the unit vector in the direction of r ij . The exponent in eq. (12) describes the screening of ith charge by the conduction electrons in metals and the value of the screening parameter ξ depends on the material and its crystallographic orientation. In our case, ξ(Cu) = 0.6809Å -1 as defined in [13]. The surface charge can be estimated according to the Gauss's law [13]˛Γ where Q is the charge inside a closed surface Γ. Using eq. 13, one can discretize the continuous distribution of surface charges on a metal surface to assign a partial charge to each atom on the surface corresponding to the value of the electric field applied at the atom position. However, it is not trivial to decide how to identify the area which contributes to the charge on a given atom. In [12,13], a rectangular geometry of the grid points associated with each atom was suggested. Yet, this approach is limited by a high computational cost for large-scale simulations.
In the current work, we selected a different approach to calculate the surface charge, which allowed for higher computational efficiency while maintaining satisfactory accuracy. This approach consists of two steps. At first, we calculate the charge of the mesh faces in the metalto-vacuum boundary, hereafter referred as face-charge: Notice, that in eq. (14) it is enough to use the scalar field and area, since E A due to the Dirichlet BC. By using a weight function w ij ≡ w( r i − r j ), the face-charge is distributed between the atoms as The total charge must remain conserved, therefore the weight function must satisfy As the exact mathematical form of w ij can be chosen quite arbitrarily, we chose a form that is computationally cheap and prevents singularities for close points: where the cutoff distance r c determines the range where the charge Q i is distributed. The tests show that on flat regions the results are almost independent of the exact value of r c until it exceeds the maximum edge length of surface triangles. For this reason in our simulations we use an empirical r c value of 0.1 nm. The face-charge method allows calculating an accurate surface charge for atoms that are located on atomic planes. Charges in the regions with protruding features require further effort to achieve satisfactory accuracy. We improve the accuracy by building the Voronoi tessellation around the surface atoms in those regions as shown in fig. 2. The facets of the Voronoi cells are assumed to constitute the surface area for the given atom to estimate the charge on it. Before generation, however, additional support points need to be built above the surface to limit the extension of the Voronoi cells into the infinity. We build such support cloud by the following process: 1. copy surface atom as a new point p; 2. locate the triangle f where the projection of p, in direction of f normal n, lies within f ; 3. move p in direction of n for distance λ; 4. loop until each surface atom has a support point.
Empirical tests have shown that good results are obtained if the length of the shift vector λ equals to one lattice constant of the atomistic system. The coordinates of the surface atoms and the support points are input in the TetGen [28] to generate a Voronoi tesselation around these atoms. Note that the bulk atoms are in fig. 2 shown for visual purposes only, and are not considered during the Voronoi tesselation generation. For each surface atom, we determine the Voronoi facets that are exposed to the vacuum, which are used in eq. (13) to assess the partial charge associated with the atom. Assuming that the electric field in the location of the facets equals approximately to the field at the atom position, E i ≈ E j , we approximate eq. (13) for atoms in the regions with increased roughness as where the summation goes over the Voronoi facets of the jth atom that are exposed to vacuum.

C. Heating
We have shown previously [12,29,30] that thermal effects caused by field emission play a significant role in the evolution of nanotips under high field. To take these effects into account, we calculate the electron emission current density J e and the Nottingham heat P N [31,32] on the emitting surface by using our field emission tool GETELEC [14]. To calculate the volumetric resistive heating power density, we use Joule's law as where σ(T ) is the electric conductivity and Φ is the electric potential obtained by solving the continuity equation in a metal which comes with the following boundary conditions: The resistive and Nottingham heat cause a nonuniform temperature distribution inside the tip. This distribution is found by solving the time-dependent heat equation in a metal where C v is the volumetric heat capacity. The initial and boundary conditions for eq. (22) are where T amb is the ambient temperature and P N is the boundary heat flux due to the Nottingham heating. In eq. (22) and (23), κ = κ(T ) is the heat conductivity as given by the Wiedemann-Franz law [33] κ(T ) = LT σ(T ), (24) where L is the Lorentz number (see next paragraph). PDEs (20) and (22) with the corresponding boundary conditions are solved in FEMOCS by means of FEM. The details about how we discretized these equations can be found in appendices C and D.
Since we aim to simulate the dynamic evolution of nanotips, the mean free path of electrons inside such systems becomes comparable with the dimensions of the tip itself. Similarly to our previous work [12], we took into account the finite-size effects (FSE) by modifying the bulk value of σ(T ) [34] by a correction factor ν = ν(T, d), as calculated by a simulation method developed by Yarimbiyik et al. [35]. The characteristic size of the nanotip d equals to the mean diameter along the tip for the initial geometry. Similarly, for the Lorentz number in eq. (24) we use an FSE-reduced value of L = 2.0 · 10 −8 W Ω K -2 , as reported for a Cu film of 40 nm thickness [36].

A. Model validation
Prior to using the proposed model for actual simulations of physical systems, we validate our space-charge and surface charge calculation methods. The field solver and emission module were verified earlier in [15] and [14].
We validate our PIC model by calculating the currentvoltage curve for a planar cathode-anode system, for which a 1D semi-analytical solution of Poisson's equation is available [37][38][39]. For this, we apply a voltage V between two parallel plane electrodes of area A = 135 nm 2 that are separated by a gap distance of d =18.2 nm and run PIC simulations for a few tens of fs until the total emission current I converges to a steady-state value.
For this geometry, the current density J can be accurately calculated by utilizing a stationary-point iteration to obtain the self-consistent values of the cathode electric field as obtained by the analytical solution of Poisson's equation and the current density J as found by the Fowler-Nordheim equation [12]. For high voltages, this curve asymptotically converges to the Child-Langmuir law [37,40,41].
The results of the PIC simulation, together with the semi-analytical and the Child-Langmuir curves are shown in fig. 3. We observe that our PIC model is in a very good agreement with the semi-analytical curve with an RMS error of 1.9%.
To validate the surface charge calculation method, we compare it against our earlier model HELMOD [13]. For this, we calculate the charge distribution on the apex of an R = 3 nm, h = 11.5 nm nanotip that is placed on a W = 93.2 nm substrate. Both the tip and substrate are cut out from a single-crystalline 100 Cu block. We measure the surface charge both with HELMOD and with the model presented here and calculate the relative difference between them. The variation of along the surface of {100} and {110} slice of the apex is presented in fig. 4. The data show that the current model tends to provide slightly smaller charges than HELMOD with a mean difference of 6.5±1%. The largest difference occurs below the atomic steps and on {110} facets.

B. Simulation setup
We used the aforementioned model to simulate the behavior of a Cu nanotip under high electric field. The expected outcome of the simulation is that the upper part of the tip will experience extensive deformations, while the lower section remains practically constant. Such a scenario helps us to optimize the computational cost by simulating in MD only the upper half of the initial system. The bottom region, that completes the geometry for calculating the physical quantities that affect the atom motion, is built separately and remains unchanged for the whole simulation. We choose different widths for this extension in order to investigate the impact of the heat dissipation on the processes at the apex of the tip.
The geometry of the tip is illustrated in fig. 5. A conical tip with an aperture angle γ = 3 • , initial height h = 93 nm and a hemispherical cap of radius R = 3 nm is placed into a simulation box of height H = 10h (H is adjusted with the change in nanotip height). The tip lies on a substrate of width W = 620 nm and height H b = 7.3 nm. The MD region is cut out from a monocrystalline 100 Cu box. The bottom layer of the atoms is fixed in place to ensure a smooth transition to the extended part. Two different extension geometries were used, with characteristic radii r = 17 nm (thin system) and r = 54 nm (wide one). The simulation geometry is motivated from the experimental observations that hint the existence of field emitters with aspect ratio of 20-100 [42] on flat Cu surfaces. The size of the tip is chosen to be large enough to prevent instability due to excess surface energy [43,44] and sufficiently small to result in a reasonable computational cost. In order to obtain results that are directly comparable to our previous work [12], we chose the dimensions of the thin system to be identical to the ones used in [12].
The MD simulations were carried out by means of the classical MD code PARCAS [45][46][47]. The FEM calculations in FEMOCS are based on open-source C++ library Deal.II [48]. For the MD simulations we used the interatomic EAM potential developed by Mishin et al. [49]. This potential has been successfully used in our previous works [12,29,43] and has shown accurate reproduction of non-equilibrium system energetics [49]. The stochastic nature of the thermal effects were taken into account by running 50 independent simulations with different ran-dom seeds. In all cases, the initial velocities were sampled from the Maxwellian distribution with T = 300 K. No periodic boundaries were applied for the MD cell, while in FEM and PIC calculations, the periodic boundaries were applied in lateral dimensions. A constant time step of 4.05 fs was used for MD, 0.51 fs for PIC and 40 fs for the heat equation. An electron SP weight of w sp = 0.01 in PIC gave a sufficiently dense sampling of the phase space, in order to ensure a smooth space-charge distribution.
We applied a long-range electric field E 0 = 0.6 GV/m for the thin tip. The lower field enhancement factor of the wide tip as compared to the thin one is compensated by adjusting the applied field to E 0 = 0.608 GV/m to ensure the same field emission current from both tips.

C. Course of the simulation
We simulated the runaway process in a nanotip similarly to the one in [12]. Despite the usage of more simplistic 1D space-charge model in [12], the overall picture of the runaway process was confirmed in current simulations. In addition, here we provide the statistical analysis of the nanotip behavior and estimate the time duration of the entire process until the tip stops emitting neutrals.
In fig. 6, we show the evolution of the height of the thin and wide nanotips that is averaged over the parallel runs. The error bars in this graph indicate the standard error of the mean height and reflect the differences in evolution of individual tips. On the graph one can observe a tendency in height variations over time in the form of a wave. Due to the high temperatures and strong forces the droplets of molten Cu once in a while leave the apex the tip. Normally this leads to an abrupt change in the tip height, but due to averaging this change is smoothened producing the shape of a wave. A rise of the tip height after its fall is seen in the majority of the simulations indicating that the runaway is able to revive. This tendency is stronger for thinner tips, since the efficient heat conduction in the wider tips cools the shortened tip faster. To illustrate the process of height change in the tips evolving under the applied electric field, some characteristic simulation snapshots designated by the corresponding letters on the left of fig. 6, are provided on the right of fig. 6.
In fig. 7, we analyze the correlation between the averaged total emitted current and the apex temperature. In this graph, we see that due to the build-up of spacecharge, the field emission current drops suddenly within the first few fs of the simulation. In the following approximately 20-25 ps, there is another slight but consistent drop in the current. This feature appears due to the shape modification of the apex region after melting, which leads to a reduction of the emitting area.
As long as the temperature in the tips does not exceed the melting point, the height of the tips does not change ( fig. 6a). After the melting point was reached, the tips start growing, causing a step-by-step increase in the local field, current and temperature (fig. 6b)  increase lasts until the thermal energy of the apex atoms exceeds their binding energy, which causes a gradual disintegration of the tip, either atom-by-atom, or in a form of clusters. The latter is often preceded by the formation of a neck below the apex ( fig. 6d). In some cases, a smaller cluster of the size up to a few atoms might also be detached ( fig. 6c). Within a necking region, the high current density causes intensive Joule heating leading to a self-amplifying increase in local temperature until the top of the tip is detached. After the break-up, the tips cool down ( fig. 6e), but within a few tens of ps, the temperature increases again and the whole process reappears. The cycle lasts until the enhanced local field is not sufficiently high to generate the field emission current capable to heat the tip above the melting point. In this case, the tips stabilize and the runaway process ceases ( fig. 6h).
Although initially the growth rate is similar in thin and wide tips, the temperature rise in the wide tip soon starts lagging behind the one in the thin tip. The lag occurs despite the same amount of heat generated initially in both systems. The more efficient cooling in the wider tip prevents the fast expansion of the region of molten Cu at the top of the tip, hence reducing also the growth rate. This is seen in fig. 8, which shows the relative height of the molten region in a tip. In this figure we see that in the thin tips, the intensive evaporation starts when about 30% in the height of the tip is molten, remaining higher than this value for the whole evaporation period. In the wide tip, on the other hand, that ratio barely reaches 15%. The molten region between solid Cu and the apex effectively acts as a thermo-electric insulator, since beyond the melting point, the thermo-electrical resistivity of Cu increases abruptly more than 60% [34]. The higher percentage of the molten region in the thin tip reduces efficiently its heat conduction, inducing a stronger accumulation of the heat at the apex of the tip. We note that figs. 6-8 show the averaged data of multiple runs. Due to averaging, we do not observe in these plots any abrupt transitions. In these graphs, we observe a smooth modulation of the data. We will refer to the well-pronounced peaks in fig. 6 as tip growth phases. During each such growth phase the tip first sharpens and after an evaporation event or a sequence of them, it blunts down until the emission currents re-heat the tip, restarting the growth phase.
However, by plotting the height evolution of individual tips, we see that the height does not change smoothly, but rather abruptly due to large evaporation events, as demonstrated in fig. 9. Here, we see that the first growth phase does not have a single well defined peak. Instead, while the tip is still in the growing phase, smaller peaks indicating detachment of small clusters of atoms appear, although these do not affect dramatically the dynamics of the tip evolution. These events are stochastic and the time interval between them may vary from a few up to hundreds of ps. Since evaporation starts already from at the beginning of the growth phase (approximately when the first small peak occurs) and continues until the large piece of molten copper detaches from the surface, the time interval t evap until the major evaporation event is not well defined. To describe quantitatively the dynamics of the field-emitting tip and enable the comparison between different cases, we measured the time intervals until the start and the end of the nanotip growth phase, t s and t e , and defined t evap = 1 2 (t s + t e ). The time intervals t s and t e are found at boundaries of half maximum of the height curve around the first intensive evaporation period (see fig. 9). During the second tip growth phase the tip was growing until a large droplet of the molten copper was detached. Since this was always a single event, we quantified the time interval until the second major evopration event by determining the time when the tip reached its second maximum in its height. We also analyzed the cooling time by determining the interval between the time instances when the tips reached their maximum and minimum heights and calculated t cool = t min − t evap . The results of these measurements are summarized in table I. As already mentioned, the two tip growth phases observed in thin tips resulted in evaporation processes with different characteristic features. The first growth phase was accompanied by several evaporation events before the detachment of the large droplet, while during the second growth phase only one large disintegration occurred. During the first growth phase, the heat is rapidly accumulating at the apex of the tip, increasing the kinetic energy of the atoms. The pulling effect of the strong local electric field increases the probability of the energetic surface atoms to fly off the surface as single atoms or as small clusters of atoms. After a few hundreds of ps, the temperature within the tip reaches an equilibrium, with a large fraction of the tip beneath the very top being molten. Such a molten tip is flexible and responds gradually as a whole to the tensile stress exerted by the field. This leads to stretching of the tip into the vacuum and a consequent necking. The initial necking leads to the increase of the current density in the narrower regions, increasing the Joule heating and thus the local temperature in this region. Eventually, a droplet forming above the neck flies off the surface leading to a large evaporation event, which completes the first growth phase.
During the second growth phase the apex temperature does not reach the same value as during the first one, see fig. 7. For that reason, during the second growth phase the tip grows gradually and no small modulations of the height similar to the first growth phase are observed. However, a larger fraction of the tip is already molten since the beginning of the second phase, rendering it flexible and responsive to the tensile stress exerted by the field. This results in a necking under the top of the tip, forming a droplet since the beginning of the growth phase. Hence, we observe only a single large evaporation event which abruptly reduces the height of the tip.
It has been previously shown [11] that a certain evaporation rate of neutral atoms is required in order to build up the vapor density necessary to trigger an ionization avalanche, which eventually leads to plasma formation. In [11], the neutral evaporation rate r Cu was considered proportional to the electron emission rate r e , with a minimum ratio of r Cu/e ≡ r Cu /r e of 0.015 required to lead to plasma formation. In our previous work [12], this ratio was found to have an average value of 0.025±0.003 over the whole evaporation process. Here we shall examine the evolution of the evaporation rate in greater detail.
For this purpose, we plotted in fig. 10 the average cumulative number of evaporated atoms and emitted electrons. We observe that the runaway process consists of two alternating regimes, characterized by intensive evaporation and a metastable stage with no or very few evaporation events. Both in the thin and wide tips, the violent evaporation at constant rate remains to the intensive evaporation regime. During this regime we observe large variations in height, current, field and temperature. In the metastable regime, on the other hand, the evaporation is practically not observed, while field emission continues. In order to quantify the mean evaporation rate within the limits of each evaporation regime, we fit a straight line to each curve (bold dashed lines in fig. 10), the slope of which gives the evaporation-emission rate. The calculated parameters are summarized in table I.
The data show, that thin tips start evaporating several tens of ps earlier than wide ones, but there is no significant difference in their cooling periods. During the first growth phase, the thin tip is able to provide neutral atoms at a rate 4-5 times higher than the wide one, while the total currents differ only by 10-30%. From these results, we deduce that thin field emitters have higher probability to trigger plasma build-up than the wide tips.

D. Discussion
The newly implemented PIC module within the hybrid MD-FEM model allowed us to achieve a more accurate description of the emission currents, as compared to that given by the previously used simplified 1D space-charge model. The new model resulted in an increased emission for the same geometry and applied field, enabling us to observe the thermal runaway process at lower applied electric field E 0 ≈ 0.6 GV/m, compared to 0.8 GV/m according to earlier estimations. This result brings the theoretical estimations closer to the experimental values, where vacuum breakdowns are commonly observed at external electric fields E 0 ≤ 0.2 GV/m [50].
The observed reduction in the field and the growth in the current values are natural. The simplified model is based on the equivalent planar diode approximation [51,52], which assumes a uniform current density in space. While simulating field emitters, however, the emission current in space is highly non-uniform due to the non-planar geometry, which is reducing the suppression of the local field due to the space-charge.
As stated earlier, we selected the applied field in our simulations, so that the total emission current at the beginning of the simulation in both thin and wide tips is the same. We took a special care for this to ensure that both tips develop initially under similar conditions. However, we see that the evaporation rates from the thin and wide tips can differ by an order of magnitude. Also, in the thin tips, we observed two growth phases which resulted in two runaway events. As it can be seen in table I, the ratio r Cu/e for the second runaway event is significantly higher than the threshold value of 0.015 that was determined in [11]. During the first runaway process the ratio r Cu/e is quite close to this value. If we consider the average rates since the first evaporation event until the end of the simulation, we obtain a value r Cu/e = 0.026±0.008, which clearly exceeds the lower limit for triggering an avalanche ionization process.
On the other hand, for the wide tip, the ratio r Cu/e is significantly lower than the one assumed in [11]. However, it is premature to conclude that thick field emitters are not able to produce enough neutral atoms to ignite plasma. First of all, in the simulations performed by using the ArcPIC model in [11], a linear relationship between atom evaporation and electron emission rates was assumed. Such a simplified assumption of the behavior of surface atoms in the pre-breakdown condition may affect the estimated threshold values severely. It is evident from our results that the linear relationship between atom evaporation and electron injection rate is valid only as a rough qualitative approximation and the actual dependency is much more complex, with intermittent high evaporation events and metastable periods, for which electron emission continues. Furthermore, in [11], a much wider electron emission and neutral injection area (10 4 nm 2 ) was assumed. In the present simulations, we can safely consider that all emission and evaporation originate from an area smaller than 100 nm 2 (a hemisphere of 4 nm radius). Given such a significant mismatch of the space scales, it is more reasonable to compare the surface evaporation and emission fluxes R Cu and R e , i.e. the rates divided by the emission areas. This is because it is the vapor and electron densities that determine whether the avalanche ionization reaction, needed to ignite plasma, can take place. These densities are not dependent on either the absolute overall evaporation or emission rates or their ratio. They are rather connected more to the fluxes R Cu and R e , the increase of which promotes plasma ignition. From table I we see that the values for R Cu and R e obtained from the present results significantly exceed the ones assumed in [11], meaning that even wide tips might provide an electron and neutral density sufficient to ignite plasma.
Although the runaway process and its connection to larger scale plasma onset simulations via the approximation of the r Cu/e coefficient or the evaporation flux R Cu offers a plausible explanation of the vacuum arc initiation, the exact mechanisms that lead to plasma formation are not yet clear. Our current PIC model includes only electron SPs and omits the interaction of the evaporated material with the electrons and the appearance of positive ions that reduce the space-charge in the vicinity of cathode, increase field emission current and trigger surface sputtering and heating. However, the present developments form the basis for future expansion of the methodology to include more particle species in the PIC domain and properly handle the plasma-metal interactions. This aspires to shed light to the detailed processes that lead to vacuum arc ignition, without the need for the rough qualitative approximations applied in [11].

IV. CONCLUSIONS
We have investigated the dynamic evolution of Cu field emitters of different width during the intensive field emission in the pre-breakdown condition. For the study we developed a concurrent multiscale-multiphysics model that combines classical molecular dynamics, finite element method and particle-in-cell techniques. The high efficiency and robustness of the model allowed us to perform extensive statistical analysis of a highly stochastic thermal runaway process. The simulations show that the thermal runaway can start in a h = 93 nm Cu nanotip under macroscopic electric fields of E 0 = 0.6 GV/m. This value exceeds only three times the value repeatedly reported in experiments with flat copper electrodes. In a thin field emitter, the thermal runaway is a cyclic process of alternating regimes of intensive atom evaporation and a subsequent cooling process of a tip. The amount of evaporated atoms that are emitted from thin emitters is sufficient to ignite self-sustainable plasma. Increasing the width of the emitter lowers the atom evaporation rate and decreases the probability for the occurrence of more than one runaway cycles. Wide emitters show also lower neutral evaporation rate, which leads to the conclusion that very sharp field emitters may be necessary to enable plasma formation.
After plugging BCs (23) into previous formula and rearranging terms, we obtain the weak form of eq. (22) To express the weak form as a matrix equation, we follow the procedure from appendix B and equalize weight function with shape functions and expand temperature as a linear combination of them. In that ways we separate spatial and temporal parts in T and eq. (D4) becomeŝ which looks in matrix form as where To solve semi-discrete eq. (D6), also the time must be discretized. It can be done by first approximating ∂T n ∂t ≈ T n+1 − T n ∆t (D7) and then introducing parameter Θ ∈ [0, 1], so that After combining (D7) and (D8) together, we get a Θscheme for solving eq. (D6): To handle the non-linearity, we assume in this formula, that thermal and electric conductivities are weakly temperature dependent, i.e. κ(t n+1 ) ≈ κ(t n ) and σ(t n+1 ) ≈ σ(t n ). In case of Θ = 0, Θ = 0.5 and Θ = 1, we get explicit Euler, Crank-Nicolson and implicit Euler scheme, respectively. In our simulations we use implicit Euler scheme while solving eq. (22). This choice helps to quickly diffuse high frequency noise in the temperature distribution that is introduced while mapping temperatures during mesh rebuild. Although implicit Euler scheme shows lower order of convergence than Crank-Nicolson method, i.e. smaller ∆t needs to be used [53], in our simulations the small time step is determined by MD and therefore this shortcoming can be safely ignored in favor of more preferential diffusive properties.

Appendix E: Distributing injected superparticles
The injected superparticles must be distributed randomly and uniformly on the quadrangle. This can be done in a straightforward way, by noting that the mesh quadrangles are built by connecting the centroids of the mesh edges and the centroids of triangles (see fig. E1a). The coordinates of random point inside a parallelogram ADEF are given by where R 1 , R 2 ∈ [0, 1] are uniform random numbers. The point r, however, might not lay inside the quadrangle ADOF , as its area is smaller than the area of the parallelogram, S ADOF = 1 3 S ABC < S ADEF = 1 2 S ABC . Therefore its location in the quadrangle must be checked and if r turns to be out of ADEF , a new point should be generated and the check procedure repeated. The location of the point can be determined with barycentric coordinates [54]. Denoting λ 1 , λ 2 and λ 3 as the barycentric coordinates of point r with respect to the 1st, 2nd and 3rd node of a triangle ( fig. E1b), r is located in the left region of the triangle (vertical hatching) if λ 1 ≥ λ 2 and in the bottom diagonal region (diagonal hatching), if λ 1 ≥ λ 3 . charge q before and after the collision, u 1,2 and v 1,2 , are written as v 1,2 = u 1,2 ± 0.5∆ v. (F1) The velocity change ∆ v is calculated from the scattering of the relative velocity, where the rotation matrix M is expressed as where a = sin θ sin φ v v ⊥ , b x,y = sin θ cos φ vx,y v ⊥ and v ⊥ = v 2 x + v 2 y . The Monte Carlo technique comes into play while calculating the azimuth and scattering angle φ and θ. To calculate the first, we generate a uniform random number in the range of [0, 2π). Scattering angles θ, however, need to be distributed normally and their contribution to matrix (F3) are found from sin θ = 2δ 1 + δ 2 , where δ ≡ tan θ is a random number generated according to the Gauss distribution. The mean of that distribution δ = 0 and the variance for a system of SPs with identical charge and mass where n is the number of particles in the cell, V cell the cell volume, 0 the vacuum permittivity and Λ = 13.0 the Landau logarithm [55]. In case there is odd number of SPs in the cell, the first 3 pairs in that cell should be collided with a variance of δ 2 = 0.5 δ 2 [27].