Thermal runaway of metal nano-tips during intense electron emission

Thermal runaway of metal nano-tips during intense electron emission

A. Kyritsakis Deparment of Physics and Helsinki Institute of Physics, University of Helsinki, PO Box 43 (Pietari Kalmin katu 2), 00014 Helsinki, Finland    M. Veske Deparment of Physics and Helsinki Institute of Physics, University of Helsinki, PO Box 43 (Pietari Kalmin katu 2), 00014 Helsinki, Finland    K. Eimre Intelligent Materials and Systems Lab, Institute of Technology, University of Tartu, Nooruse 1, 50411 Tartu, Estonia    V. Zadin Intelligent Materials and Systems Lab, Institute of Technology, University of Tartu, Nooruse 1, 50411 Tartu, Estonia    F. Djurabekova Deparment of Physics and Helsinki Institute of Physics, University of Helsinki, PO Box 43 (Pietari Kalmin katu 2), 00014 Helsinki, Finland
July 12, 2019

When an electron emitting tip is subjected to very high electric fields, plasma forms even under ultra high vacuum conditions. This phenomenon, known as vacuum arc, causes catastrophic surface modifications and constitutes a major limiting factor not only for modern electron sources, but also for many large-scale applications such as particle accelerators, fusion reactors etc. Although vacuum arcs have been studied thoroughly, the physical mechanisms that lead from intense electron emission to plasma ignition are still unclear. In this article, we give insights to the atomic scale processes taking place in metal nanotips under intense field emission conditions. We use multi-scale atomistic simulations that concurrently include field-induced forces, electron emission with finite-size and space-charge effects, Nottingham and Joule heating. We find that when a sufficiently high electric field is applied to the tip, the emission-generated heat partially melts it and the field-induced force elongates and sharpens it. This initiates a positive feedback thermal runaway process, which eventually causes evaporation of large fractions of the tip. The reported mechanism can explain the origin of neutral atoms necessary to initiate plasma, a missing key process required to explain the ignition of a vacuum arc. Our simulations provide a quantitative description of in the conditions leading to runaway, which shall be valuable for both field emission applications and vacuum arc studies.

I Introduction

Field electron emission plays a crucial role in numerous modern technological applications, from electron microscopy to flat displays Egorov and Sheshin (2017). However, if the applied field at an emitting cathode is increased beyond a limit, violent discharges in the form of arcs appear in the vacuum. Vacuum arcs, also known as vacuum breakdowns, play a significant role in various technological applications by being either exploitable or highly undesirable. For example, they can be used in ion sources Brown et al. (1985) or in physical vapor deposition Anders (2008). On the other hand, they cause catastrophic emitter failure in electron sources Hartmann and Gundersen (1988); Dyke and Trolan (1953); Descoeudres et al. (2009a). Moreover, they hinder the function and limit the performance of various vacuum devices that require high electric fields, such as fusion reactors McCracken (1980), vacuum interrupters Slade (2007) and powerful linear accelerators to be used in particle colliders for new insightful experiments at CERN (CLIC) (Aicheler et al., 2012).

Over decades, vacuum arcing strongly attracted the attention of researchers from different fields Dyke and Trolan (1953); Dyke et al. (1953); Dolan et al. (1953); Hartmann and Gundersen (1988); Anders et al. (1993); Latham (1995); Anders (2008); Descoeudres et al. (2009b), since various physical phenomena are involved. However, the mechanisms of plasma onset are still under debate. During an arc, plasma is ignited in the vacuum and burns until the available energy from the power source is consumed. It is well-known since the 1950’s (Dyke et al., 1953; Dolan et al., 1953; Descoeudres et al., 2009b; Anders, 2008) that vacuum arcs appear after intense field electron emission, but the physical mechanism that leads from the latter to plasma ignition is not yet fully understood.

Recent experimental studies show that the improvement of surface and vacuum quality cannot diminish the probability of vacuum arcing Degiovanni et al. (2016). Moreover, arcs were found to appear on single metal tip cathodes when their current density exceeded a critical value even for well-controlled emission conditions Dyke and Trolan (1953). These indicate the existence of an inherent mechanism of surface response that initiates breakdown. One hypothesis commonly used to explain the plasma build-up is the ”explosive emission”(Anders, 2008; Mesyats, 1995, 1993, 2005), which assumes that an instant explosion of the emitting tip leads to plasma formation. However, this hypothesis is mainly based on phenomenological considerations, since the atomic level insight in such a complex phenomenon as vacuum arcing was not available.

Here we present multi-physics atomistic simulations that reveal a thermal runaway process on emitting nano-tips. The latter is based on the gradual deformation of an emitting tip, due to the field-induced forces. The simulations show that both high current density and sufficiently large tip size are needed to initiate the melting at its apex. After that, the field-induced forces gradually deform the tip, elongating and sharpening its apex in a process similar to the Taylor cone formation in liquid metal ion sources Krohn and Ringo (1975); Wagner (1982); Swanson (1983). This process, in combination with the decrease of the electrical and thermal conductivities at high temperatures, leads to a positive feedback thermal runaway mechanism. Eventually, large fractions of the emitter evaporate in the form of neutral atoms, but also as charged nano-clusters. The total number of evaporated atoms (both isolated and clustered) is compatible with the minimum evaporation rate required to ignite plasma, as found by recent Particle-In-Cell (PIC) calculations Timko et al. (2011, 2015).

Ii Method

In order to study the atomic level processes leading to the initiation of plasma, we perform Molecular Dynamics (MD) simulations. However, the thermal runaway process involves various physical phenomena, which classical MD cannot describe implicitly due to its atomistic nature. The electric field interacts with the material, inducing charges and forces on the surface atoms and generating electron emission. The latter heats the nano-tip due to the Nottingham Nottingham (1936); Charbonnier et al. (1964) and Joule effects. The combination of the above causes significant changes in the structure and geometry of the material and hence they have to be quantified and included in the calculations.

In this article we propose a multi-physics model, which extends the capability of the classical MD method by solving the electrostatic, electron emission and heat diffusion equations concurrently for a dynamically evolving nano-emitter shape. In the next three sections we explain how these effects are incorporated in the MD simulations. In appendix A a flowchart of the whole simulation algorithm is given. Figure 1 illustrates a schematic representation of the tip geometry that along with a compilation of the partial differential equations with their corresponding boundary conditions used to obtain the electric field and temperature distributions.

Figure 1: Schematic of the partial differential equations and their corresponding boundary conditions.

ii.1 Electric field

In order to introduce the electrostatic phenomena, we have updated our HELMOD model Djurabekova et al. (2011) that couples MD simulations with electric field calculations, and includes the charge-induced forces in the MD interactions. Although HELMOD provided many interesting insights into the interaction of high electric fields with metal surfaces, the rigid grid used to calculate the electric field imposed critical limitations for large systems, dynamically evolving at high temperatures. To overcome this limitation we used our recently developed model FEMOCS Veske et al. (2016a, 2017) to calculate the electric field distribution around the tip on a flexible mesh and integrate this information into HELMOD.

As HELMOD, FEMOCS also solves the Laplace equation


but uses the Finite Element Method (FEM) on an unstructured mesh. This mesh is generated automatically from the atomic positions in order to follow the varying nano-tip shape. A Dirichlet boundary condition


is then applied on the nano-tip surface and a Neumann


at the top of the simulation box, where is the applied macroscopic field. Once the electrostatic potential and field


are obtained, we calculate the charges induced in each atom due to the field (details on how this charge is calculated are given in appendix B). Given the charge on each atom, we calculate and add to the MD interaction the Lorentz force


The term stands for the fact that the atoms are exposed to the electric field only from the one side of the material-vacuum surface. The inter-atomic Coulomb forces are also calculated and included in the MD interaction by the method described in ref. Djurabekova et al. (2011).

The electric field effects are calculated concurrently with MD, so that the field distribution evolves according to the emitter shape. When the emitter shape changes significantly, the electric field distribution is recalculated. At every timestep, the position of the atoms is compared to the last full-calculation step. If the root-mean-square (RMS) average of the atomic displacement (compared to the last field calculation step) is smaller than 0.38Å, the last calculated field solution is reused. This technique is used in order to keep the computational time of our simulation feasible. Choosing the value for the threshold is a trade-off between computational efficiency and approximation accuracy. The aforementioned threshold was used because it reduced the CPU time by about two orders of magnitude, while the level of accuracy of the field calculations deteriorated only by 2% as compared to the accuracy o A full analysis of the CPU time and field calculation error as a function of the threshold can be found in ref. Veske et al. (2017).

ii.2 Electron emission

Our previous calculations Parviainen et al. (2011); Eimre et al. (2015) showed that the thermal effects caused by electron emission play an important role in the nanotip evolution under high fields. The classical electron emission equations Schottky (1923); Fowler and Nordheim (1928); Murphy and Good (1956) are inadequate to describe combined thermal-field emission Jensen and Cahay (2006) from nanometrically sharp emitters He et al. (1991); Kyritsakis and Xanthakis (2015, 2016). However, our recently developed computational tool GETELEC Kyritsakis and Djurabekova (2017) provides with the means to consistently and efficiently calculate the electron emission current and Nottingham heat distributions from sharp nano-tips, even at temperatures beyond the melting point. We also note that the local fields and the current densities calculated here are close to the Space Charge (SC) limit; thus the SC effect Child (1911) is also included in our model.

Calculating fully the charge density distribution induced by the emitted electrons in three dimensions (3D) is a complicated and computationally expensive calculation Jensen (1999), which is usually done with the PIC method Chen et al. (2009); Uimanov (2011). Nevertheless, Forbes Forbes (2008) claimed that the standard 1D SC model for field emission (Barbour et al., 1953) can be used to obtain the reduction of the surface field if a correction factor is introduced to account for the non-planar nature of the emitter.

According to the standard 1D SC model, the local field on the emitter (uniform in 1D) is found as , where is the dimensionless field-lowering factor and is the field found by the Laplace equation, i.e by ignoring the SC effects. is a function of the emitted current density , the local field and the total applied voltage , which is here an external parameter.

However, for a 3D nanometric emitter, and are not well-defined. In a quasi-planar emitter, the values at the apex can be used, since the SC is localized around it in a region much smaller than the radius of the emitter, where these quantities are practically constant. On the contrary, in a nanometric emitter the whole tip surface contributes to the SC, which means that the apex values are no longer representative.

Here we will define the following representative values for the field and the current density correspondingly


In eq. (II.2) the surface integrals are performed on the full-width-half-maximum emission surface , i.e the surface where . On this surface, is the mean value of and is the weighted mean value of , with the weight being the emitted current density.

The field suppression factor can then be obtained by solving the equation Forbes (2008)




where is the electron mass, the elementary charge and the vacuum permittivity.

The whole potential and field distributions on the emitter are then multiplied by and the current density distribution along with are recalculated according to the new field distribution. This procedure is repeated iteratively until self-consistency is reached, i.e. the - suppressed field and the corresponding current density distributions give that satisfy eq. (7), (8). Note that the obtained self-consistent is used to multiply the whole electric field distribution, thus affecting not only the electron emission, but the charge and force calculations described in section II.1.

In order to validate our 1D space charge model, we compare its results to full PIC simulations performed by Uimanov Uimanov (2011). We reproduced the solution of the Laplace equation for the geometries presented in ref. Uimanov (2011), and then used our model to obtain the field-lowering factor , using a corrected applied voltage . is obtained for various values of the applied voltage and for two of the geometries simulated by Uimanov: ”FE4” and ”FE5”. ”FE4” is an emitter with 10nm radius and ”FE5” with 1nm, thus being the most relevant to our simulated tips that have a radius nm (see section II.4).

In figure 2 we see the (values at the apex) curves reported in ref. Uimanov (2011) (markers) along with the results of our model (lines) for the same geometries. The corresponding fitted values for the correction factor are shown in the legends. We see that our model is in very good agreement with the full 3D PIC calculation with an error of less than 1%, at a range of up to about 35GV/m.

Figure 2: Local apex field including the SC versus the one ignoring SC for the two different emitter geometries simulated by Uimanov Uimanov (2011). Markers correspond to the results of ref. Uimanov (2011) and the lines to the results obtained by our 1D SC model.

In view of the above, the 1D SC model described above is a good approximation for nanometric size emitters, given that a valid correction factor is used. Although we cannot obtain here its value for the specific geometry we simulate, the fitted values obtained above for similar hemisphere-on-a-cone geometries indicate that the used value is reasonable for our nano-tip. This value was chosen for nm, by linearly interpolating between and , with respect to the corresponding .

ii.3 Heating

In order to obtain the temperature distribution on the emitter, we solve the heat diffusion equation Bejan and Kraus (2003)


where is the volumetric heat capacity, is the heat conductivity and is the local deposited volumetric heating power density (in ).

The geometrical structures which we simulate, have one dominant dimension, i.e. their height is much longer than the lateral dimensions. Hence eq. (9) can be reduced to one dimension


where is the cross-sectional area of the tip at a given height and is the mean deposited heat on the corresponding cross-section. Details on the exact definition of and the derivation of (10) can be found in appendix C.

The above equation is solved using the Finite Difference Method (FDM) with an explicit Euler scheme. We divide the tip in finite slices of width , that coincide with the slices of the rectangular grid used for the charge calculation (see appendix B). On each surface grid cell, the current density and the deposited Nottingham heat surface density are calculated as described in section II.2. Then the total current and total Nottingham heating power contributions of each cell ( and correspondingly) are calculated by multiplying the densities with the cell face areas, similarly to the charge calculation described by equations (14) and (15). We consider that the emission quantities have the same direction as the electric field of the cell.

The total current flowing through the cross-section of the -th ( at the apex) slice is cumulative from the apex to the base as dictated by the continuity equation, i.e


where denotes the cells of the -th slice that are exposed to vacuum. The total heat deposited on the -th slice is the sum of the Joule and the Nottingham heat components, i.e


For the Joule heating, we need to obtain the electric conductivity which depends on the local temperature and the size of the nano-tip. For a given temperature , we obtain by linearly interpolating tabulated values found in the literature Matula (1979); Gathers (1983) for the resistivity .

is capped at the lowest available value of 3500K, i.e. . Furthermore, the mean free path of the electrons in the material is decreased due to the nanometric size of the nano-tip. To correct our values for this finite-size effect, we use the simulation method and tool of Yarimbiyik et. al. Yarimbiyik et al. (2006). The latter can calculate the mean free path reduction for a given nano-wire diameter. The value we use for the latter is the mean diameter along the tip for the initial geometry.

The heat conductivity is calculated from according to the Wiedemann – Franz law . The Lorentz number is known to be reduced for nanometric size structures Nath and Chopra (1974). We used its smallest value reported in ref. Nath and Chopra (1974) for a Cu film of thickness. We note that for a nano-wire of smaller size such as the ones we simulate, might decrease even more, which will promote the thermal runaway even for smaller fields or emitter sizes.

The Euler scheme for equation (10) is . Discretizing equation (10) yields


where the subscript denotes the slice number for all quantities. A constant Dirichlet boundary condition is applied at the bottom of the tip, i.e K. Note that in (13) values for are needed to obtain . Thus we assume a virtual point for which and . The latter corresponds to a Neumann boundary condition of zero heat flux at the apex.

The forward difference Euler scheme of eq. (13) is numerically unstable if and are not sufficiently small. Here we used and . is determined by the rectangular grid which might change as the MD simulation evolves and is approximately half a lattice constant. Our tests showed that these values are small enough to safely avoid numerical instabilities. Besides, the computational cost of iteratively updating eq. (13) is negligible comparing to the rest of the simulation, thus there is no need to optimize their values as long as they are small enough.

The conductivities are updated at each heating time step according to the current . On the other hand, and are recalculated and updated when either the field distribution is recalculated (meaning that the atoms have been displaced more than a limit) or has an RMS difference greater than 1% compared to the last step when and where obtained.

Finally, the resulting temperature distribution is used to control the temperature of the MD system. The velocities of the atoms residing inside the th slice are scaled according to a Berendsen control scheme Berendsen et al. (1984) with control temperature and relaxation time . This value is much smaller than the relaxation time of the heat equation, but also big enough to avoid artefacts appearing in MD when intense velocity scaling is applied.

The above 1D model has the advantage of its simplicity and computational efficiency. In order to validate the model, we compare its results with our previous 3D Finite Element Method (FEM) model Eimre et al. (2015) (see appendix D for details). This comparison shows that the 1D heat model described above produces results that are in excellent agreement with ones from 3D FEM.

ii.4 MD simulation set-up

We simulate a conical Cu nanotip whose major axis aligned with the crystallographic direction (see figure 3). The cone is terminated with a hemispherical cap of a radius . The tip has a total height and a full aperture angle of . In order to keep the computational time feasible, only the upper half of the tip is fully simulated with MD. According to our estimations, only this part is significantly heated, providing high kinetic energies to the atoms. The bottom half maintains practically a constant shape and we consider it fixed; the field and temperature calculations are extended to the continuous limit in that region. Thus the total number of simulated atoms is 206000 and the MD simulation box has a size of nm (shown with grey lines in figure 3).

Figure 3: Simulated tip shape illustrated by the MD atoms (in the box) and the constant ”virtual” atoms (out of the box). The color coding represents the steady-state temperature distribution for constant shape. The coordinate reference system is placed in the middle of the tip at the bottom of the MD domain. The inset illustrates the atom charges by the color coding and, the corresponding forces by the black-line arrows.

However, the electrostatic and heat equations are solved on the whole tip domain, which includes an extension of the MD system shown by the ”virtual” particles lying out of the box in figure 3. These points are used to generate the FEM mesh on a nm box. This size is sufficient to assume periodic boundary conditions at the sides and a Neumann boundary condition at the top. Only a small part of this box is shown in the figure. A macroscopic electric field with a total voltage kV is applied instantly after the MD system is run to relax for 0.5ps. These values for are relevant to the surface fields used in modern high-gradient accelerating structures Degiovanni et al. (2016); Wu et al. (2017).

For the MD simulations we used two different inter-atomic potentials with a constant timestep fs. The Sabochick and Lam (SL) EAM potential Sabochick and Lam (1991) and the one from Mishin et. al. Mishin et al. (2001). In order to take into account the stochastic nature of the process, we repeated the same simulation 7 times for each potential. Each repetition was initialized with different atomistic velocities (randomly sampled from the Maxwellian distribution with ) by giving different seeds to the corresponding random number generator.

Since our tip is isolated and continued to the bottom, no periodic boundaries were applied. The two bottom atomic layers of the MD systems were fixed to obtain a smooth connection to the rest of the tip and avoid translational motion of the whole structure. Finally, the ten bottom layers above the two fixed were controlled to linearly increasing temperatures from 0 to in order to avoid possible artefacts caused by extreme velocity mismatch. The temperature on the remaining of the tip was controlled according to the calculated .

We note that the total computational time needed for the simulation is quite high. Every simulation runs on a single core and it needed approximately 10 days to run the presented 350ps simulation time either on a standard workstation or cluster. The memory consumption of the simulation does not exceed 2GB.

Iii Results

After the electric field is applied, the maximum local field at the apex as calculated by the Laplace equation is . However, such a high field induces an emitted current density with an average value at the bottom of the MD simulation domain (middle of the conical tip) . This will suppress the fields due to the SC effect by a factor of and converge to a final current density of . The inset of figure 3 shows a closer view on the tip apex with the color coding corresponding to the calculated field-induced charges for each atom and the black arrows corresponding to the forces . All the numerical results reported throughout the next two sections in the format ”” or with plot error bars correspond to the mean values and standard errors as obtained from the various simulation runs with different random seeds.

Figure 4: Time evolution of the maximum temperature on the tip of a constant shape (dashed blue line) and of a variable shape as calculated by the multi-scale atomistic method (solid red line).

Figure 4 demonstrates the evolution of the maximum temperature on the nano-tip for one of the simulations with the SL potential. In the initial geometrical configuration, the deposited heat power density is very high (reaches at the apex) and is dominated by the Nottingham effect. If the heat equation is solved under constant geometry, the heat is rapidly dissipated towards the bulk due to the high thermal conductivity of Cu and the temperature distribution converges to a steady state, with a maximum value of about K at the apex.

Figure 7: (a) Evolution of the height of the tip (red line, left axis) and the Cu atoms evaporation rate (blue dashed line, right axis). The designated points (a-g) correspond to the frames depicted in figure (b). The color coding of (b) corresponds to the local temperature. All evaporated atoms and nano-clusters have been removed. The full animation is available in the supplementary material.

On the other hand, if we let the tip shape evolve according to the multi-scale method described in section II, we obtain a completely different picture (solid red line). The flexible shape of the molten material on the top of the tip significantly affects the field and temperature distributions; this causes a thermal runaway process leading to extremely high temperatures.

Figure 10: Time evolution of: (a) the total deposited Nottingham heating power (solid blue line), Joule heating power (dashed red line) and (b) the representative electric field (solid blue line, left axis) and the total emitted current (dashed red line, right axis). The inset of (a) shows the same data in their full range.

This shape-related runaway process is illustrated in figures 7 and 10 for the same simulation. The solid red line of figure (a)a shows the time evolution of the height of the tip. Figure (b)b demonstrates the snapshots (a-g) of the tip evolution at the time steps designated on the graph respectively. Figure (a)a shows the corresponding evolution of the deposited heat components (dashed red for Joule and solid blue for Nottingham) and figure (b)b gives the evolution of the representative field of eq. (II.2) (solid blue line) and the total emitted current (dashed red line - right axis) for the same simulation.

The tip shape initially stays approximately constant (), while its temperature gradually increases towards its steady-state value of 1820K. In this stage the Nottingham heating power is dominant (about an order of magnitude higher than Joule), but as the temperature increases it slowly reduces. The steady-state temperature exceeds the melting point at the apex region and the corresponding atoms become mobile. Then the field-induced forces can pull them upwards, thus sharpening and elongating the tip. This leads to a slow increase of the local fields, as shown in the figure (b)b (ps). The increased field induces higher forces that elongate and sharpen the tip even further.

At this stage the total emitted current stays roughly constant (see figure (b)b) because the emission is deeply in the space-charge limited regime. Nevertheless, the increase of the tip hight in combination with the increased resistivity hinders the heat conduction and causes an increase in the Joule heat and temperature. After about ps, the above processes have formed a positive feedback loop leading to thermal runaway. The system height, field, temperature and Joule heat components constantly increase feeding the increase of each other.

Shortly after ps the deformed tip has exceeded 3000K (frame b in figure 7). At this temperature the Nottingham effect is inverted, producing cooling instead of heating due to the increased thermionic component of the emission (blue curve in figure (a)a becomes negative). However, the high temperature also results in significantly increased resistivity. Furthermore, the field-induced forces cause neck-thinning which becomes more pronounced as the simulation evolves (frame c). This thinning confines the total emitted current in a small cross-sectional area, thus producing a high local current density. The combination of the above produces high Joule heating, which exceeds the Nottingham cooling and leads to a rapid increase of the local temperature exceeding 10000K for a very short time. The emitted current at this point increases abruptly because of the thermionic component of the emission. We see that the corresponding field decreases rapidly at the same point due to the space charge limitation, but this does not affect the emission due to its thermionic nature at that point.

Eventually, the high temperature in combination with the field-induced forces cause detachment of the upper part of the tip, creating an evaporated charged nano-cluster. The latter is charged due to the partially charged atoms on its surface. In the one presented in figure (b)b it consists of 1784 particles that bear a total charge of -23.8e.

Although the full behaviour of evaporated nano-clusters requires further investigation (interaction with the field, collision with the evaporated atoms, emitted electrons, etc.), for the purpose of this work we consider them instantly removed and continue the simulation with the remaining tip. This is justified by the fact that they will be rapidly accelerated by the field towards the anode due to their negative charge.

All evaporated atoms and clusters are marked and removed by a cluster analysis algorithm Ester et al. (1996) implemented in FEMOCS Veske et al. (2017), with a distance cut-off equal to the one of the MD potential ( for SL and for Mishin). Figure (b)b shows only the non-evaporated atoms. Note that we have not explicitly forced any evaporation event. They appear naturally in the MD simulation due to the high temperature and the field-induced forces.

The dashed blue line of figure (a)a (right axis) demonstrates the evolution of the evaporation rate of Cu atoms. The peaks in that line correspond to the evaporation of nano-clusters. Shortly after the evaporation of the first cluster, the tip is very hot and many atoms evaporate either as isolated or clustered. The peaks on frames (d) and (e) in the height and evaporation rate graphs correspond to two relatively large clusters being detached.

After ps the tip height and temperature gradually decrease, and the system returns to a state similar to the one before ps. However, at about ps the runaway process is reinitialized in the remaining tip until a new large cluster is evaporated at frame (f). After that, the process finally stops as the tip blunts and cools down (g). The full animation of the evolution of the tip can be found in the supplementary material.

Although the runaway and the evaporation appeared at different times and the tips took various shapes in the 14 simulation runs with different random seeds, the steps described above remained qualitatively similar. After a certain time a runaway process leads to a large cluster evaporation (corresponding to frame (c) of figure 7) which is followed by smaller runaway–evaporation events (frames (d), (e) and (f)).

The statistical behaviour of the events is demonstrated in figure (a)a, where we plot the time when each event appeared along with the size of the corresponding evaporated nanocluster. The behaviour of these events was slightly different for the two potentials, hence we plotted their statistics separately. The SL potential gives the first runaway significantly earlier than the Mishin one, while the size of the corresponding nano-cluster is smaller. However, after the first event, the SL simulations produce more and larger events, because the corresponding height has not been reduced as much as for the Mishin. Thus, the mean number of events is for the SL potential and for the the Mishin (there is no simulation with 5 events for Mishin).

We see that the events shown in figure (a)a are quite localized in time (cf. small corresponding error bars), which implies a relatively small stochasticity of the runaway-evaporation process. This means that the latter substantially differs from the usual thermal evaporation, that would appear as purely random. Thus, although the runaway is a thermally-activated process, it is mainly driven by the electric field.

Figure 13: (a) Time (left axis - blue solid lines) that each evaporation event appears, along with the number of atoms of the corresponding evaporated nanocluster (right axis - red dashed lines). Round markers correspond to calculations with the Mishin potential square ones to the SL. Events #3 and #4 appeared only for one of the runs with the Mishin potential, hence the absence of error bars. (b) Time–averaged evaporation rate (bars). The width of the bar corresponds to the time interval over which the evaporation has been averaged. The error bars correspond to the standard error of all the runs for both potentials.

The events described above, along with the atomic thermal evaporation create a statistical distribution for the evaporation rate, which is shown in figure (b)b. Due to the sharp peaks of the various events, the mean rate is very noisy. For this reason we plot its time-average over 5000 MD time-steps (10.13 ps). The noisiness of this quantity is also the reason for the large error bars in the figure. The two peaks in the rate correspond to the times that the first evaporation event occur for the two different potentials.

Iv Discussion

An arising question is whether the above process is enough to initiate self-sustaining plasma. In a recent work Timko et al. (2011, 2015), PIC simulations showed that plasma can build up in the vicinity of an intensively field emitting spot. This happens if the cathode, emits not only electrons, but also neutral Cu atoms at a rate of at least 0.015 neutrals per emitted electron. However, the physical processes that can lead to the emission of neutrals along with electrons was not fully understood.

We define the average evaporation rate as the total number of atoms detached from the tip (isolated or clustered) over the total time interval between the first and the last evaporation events. Our calculations give atoms/fs with the total number of evaporated atoms being . If we divide by the average emitted electron current during the same period e/fs, we obtain the mean evaporation rate per emitted electron atoms/e. This rate is strikingly close to the values estimated in Timko et al. (2015) and even exceeds the reported minimum of required to ignite plasma.

The consistency of our results with the previous independent simulations that used a different method (PIC) indicates that there is an intrinsic mechanism able to supply neutral atoms sufficient to ignite a vacuum arc near the cathode surface. Our calculations based on the state-of-the-art understanding of the electric field–material interactions suggest that an external source of neutral atoms might not be necessary to ignite the plasma; the latter can also be fed through self-evaporation of the surface atoms, in a process much less violent than the ”explosive emission” conventionally assumed Anders (2008); Mesyats (1995, 2005).

Another important question is what are the prerequisites to initiate thermal runaway. The latter is a complicated process that depends on various initial configuration parameters, namely material, geometry, applied field and voltage. A full analysis of all these parameters is out of the scope of this work and will be given in a forthcoming publication. However, as a general comment, two are the main contributing factors to the initiation of the thermal runaway: melting and force.

This means that the tip height and the current density have to be sufficient to cause melting at the apex region. Here we have to stress out that the height of the tip itself plays a significant role for the heating of the tip, not only because it increases the field enhancement and thus the emitted current, but also because it decreases the heat dissipation, thus producing higher temperatures at the apex for the same deposited heat.

height [nm] [K]
63.1 1045
73.1 1067
83.1 1099
93.1 1135
103.1 1152
Table 1: Maximum steady-state temperature for tips of various heights. The applied field is adjusted so that the average current density is .

This is demonstrated in table 1 where we show the maximum steady-state temperature (constant shape – no MD simulation) for tips of various heights. The shape of the upper part of the emitter is the same as in figure 3. The lower ”extension” part is adjusted to the given total height. The applied field is adjusted so that the mean current density is kept constant at . We see that although the current density is the same, the tip height itself affects the maximum temperature.

In general, the current density required to produce sufficient heating to reach melting temperatures is of the order of , in agreement with the experimental results of Dyke et. al Dyke et al. (1953); Dyke and Trolan (1953). Moreover, the importance of tip melting prior to an arc has been observed also experimentally by Batrakov et. al. Batrakov et al. (1999).

On the other hand, the balance between the field-induced forces and the surface stress plays also a significant role. For the simulations we presented here, a minimum local field at the apex of about 10-12GV/m is necessary. In view of the above, it would be safe to assume that tips with a height of at least several tens of nm are required in the case of Cu. The existence of nano-tips of such size, similar to the one we simulated here, is motivated by experimental observations. Field emission measurements on flat Cu surfaces, have exhibited high enhancement factors of the order 20-100 Kildemo (2004). Our previous calculations Veske et al. (2016b); Jansson et al. (2016) has shown that Cu nano-tips that are smaller than 1-2nm in radius cannot be stable. Thus, assuming a minimum tip radius of 1-2nm, the tips involved in arc initiation should have a height of at least several tens of nm, which is the main motivation for the geometry chosen to be simulated here.

V Conclusions

In conclusion, we have simulated the thermal and shape evolution of intensively electron emitting Cu nano-tips by means of a multi-scale atomistic model. Our results reveal a thermal runaway process initiated by the field-induced forces acting on the molten apex of a nano-tip. The thermal runaway leads to evaporation of metal fractions in the form of either atoms or nano-clusters, at a rate that exceeds the minimum needed to ignite plasma. This is a self-sufficient process, which does note depend on an external source of neutral atoms. Thus we show that the onset of a vacuum arc in ultra high vacuum is an intrinsic response of a metal surface to the applied high electric field.


The current study was supported by the Academy of Finland project AMELIS (grant No. 1269696), CERN CLIC K-contract (No. 47207461), Estonian Research Council Grants PUT 57 and PUT 1372 and the national scholarship program Kristjan Jaak, which is funded and managed by the Archimedes Foundation in collaboration with the Ministry of Education and Research of Estonia. We also acknowledge grants of computer capacity from the Finnish Grid and Cloud Infrastructure (persistent identifier urn:nbn:fi:research-infras-2016072533).

Appendix A Simulation flowchart

Figure A.1: Flowchart of the simulation algorithm.

Appendix B Field-induced charge calculation

In order to calculate the charges and forces on the surface atoms, we use an extension of our previously developed molecular dynamics – electrodynamics method HELMOD Djurabekova et al. (2011). We build a rectangular grid on the Molecular Dynamics (MD) simulation box so that every atom belongs to a certain grid cell. Since the atoms are in an FCC crystal (at least when the temperature is still low), there are numerous grid cells (roughly half) that contain no atoms. Furthermore, when the temperatures approach or exceed the melting point there might be some cells (a small percentage) containing more than one atoms, but this does not affect the functionality of our method.

Figure B.1: Atoms on the tip apex (maroon spheres) after a relaxation of 0.5ps along with the corresponding material-designated cubic grid cells. The color coding on the cells corresponds to the charge assigned to each cell. The inset demonstrates a closer view in the box area, with the grey small cubes being vacuum-designated cells.

When each atom is assigned to a certain cell, the cells are separated into material and vacuum domains. Figure B.1 illustrates this separation. The cubes correspond to the cells of the material domain, for the initial configuration of the simulated tip after a relaxation of 0.5ps. The maroon spheres correspond to the atomic positions. In the inset a closer view of the area in the box is demonstrated, with the smaller grey cubes corresponding to some of the cells of the vacuum domain.

The separation is done in the following way. If a cell contains an atom or has a number of first nearest neighbours containing an atom , then it is designated as material (e.g. cell in the inset). If it does not contain an atom and it has it is designated as vacuum (e.g. cell in the inset).

Then the surface charges are calculated in every cell face separating a vacuum cell from a material cell according to the Gauss law


where is the surface charge on each cubic face, is the local electric field in the center of the face, and is the unit vector normal to the cubic surface. The electric field distribution has already been calculated by solving the Laplace equation using the recently developed Finite Element Method (FEM) tool FEMOCS Veske et al. (2017). Then the charge on each surface cubic cell is calculated as


where the summation is done over all the faces of the cell that are in contact with a vacuum one. The charges are shown in the color coding of the cells in figure B.1. Note that inner cubic cells that have not a common face with vacuum cells have .

Since we want to calculate the inter-atomic Coulomb interactions and the field-induced Lorentz forces, the charge on each atom needs to be calculated. For this, all the grid-point charges have to be assigned to atoms while the total charge is conserved. This is done in the following way. If a grid point contains one or more atoms, then all its charge is assigned to those atoms (divided equally). In the vast majority of the cases there is only one atom in each grid cell. On the other hand, if a grid cell does not contain any atoms, its charge is distributed equally to the atoms belonging to its first-neighbouring cells.

Appendix C One-dimensional heat equation

Since the length of the simulated structures is dominant, let us integrate eq. (9) on an infinitesimal ”slice” , between with width , cross-sectional area and volume


The integral of the left hand side equals to where denotes the mean temperature on . The second term in the right hand side equals to the total heat deposited on the slice. Using the above notation and Gauss’s divergence theorem we can rewrite eq. (16) as


where denotes the unit vector perpendicular to the surface around the slice . This surface has three components. The two cross-sections and and the ring around them. For the contribution of the ring is negligible compared to . Moreover, since the cross-sectional surfaces are along the plane, we have and the surface integral on the cross-section can be expressed as . Then eq. (17) yields


Dividing eq. (18) by and considering that the temperature does not vary significantly in the lateral directions, i.e. , yields eq. (10), the 1D heat equation for a variable cross-section and conductivity wire.

The total heat is given by adding the Joule and the Nottingham heat components. The Nottingham heat is expressed as surface density (in ), but it can also be expressed as a volumetric quantity, localized in the surface region with the help of the Dirac function. The total deposited heat at a given point is


where is the local electric conductivity, and is the Nottingham heat at a surface point . Then the mean heat on a cross-section is given by integrating (19) over the cross-section A(z). The surface integral of the second term turns into a line integral due to the function and it yields:


In the above equation is the total current flowing on the cross-section , and stands for the closed line surrounding the cross-section at .

Appendix D Validation of the one-dimensional heat equation

The 3D FEM model of ref. Eimre et al. (2015) is used to solve the steady-state heat equation on the same FEM mesh produced by FEMOCS. The emitted current and Nottingham heat are in that case is calculated using the standard General Thermal-Field (GTF) Jensen and Cahay (2006). Only for this comparison, we will ignore the nanometric emitter size effects in electron emission and the space charge effect. Note that originally the FEM model ignored the Nottingham effect. Nevertheless, it was here implemented by adding a Neumann boundary condition to the heat equation, similar to equation (4) of reference Eimre et al. (2015). The boundary condition is


where denotes the surface density of the deposited heating power and is the unit vector normal to the material surface.

Figure D.1: Steady state temperature distribution as calculated by the 1D FDM (blue lines) and the 3D FEM (red markers) models for two different tip geometries and field values. For the regular conical geometry ”regular” (solid lines) the applied field was and for the oblique geometry ”irregular” (dashed lines) .

Figure D.1 shows the comparison between the two models. The lines correspond to the final steady-state temperature distribution along the simulated tip, as calculated with the 1D model described above. The markers correspond to the temperature along the vertical axis of symmetry of the tip, as calculated by the 3D FEM model.

The calculation was performed for two different tip geometries, a ”regular” and an ”irregular”. The ”regular tip” geometry (solid lines and dots) corresponds to the original conical geometry (shown in figure 3) and the one for ”irregular tip” one(dashed lines and diamonds) to the deformed shape of frame (b) in figure (b)b. We see that the simplified 1D model is in excellent agreement with the 3D FEM calculation, even for the irregular geometry. The markers correspond to the mean value of the temperature for a given slice. The deviation of the temperature around that value is smaller than the size of the marker in the figure.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description