Dynamic coupling between particle-in-cell and atomistic simulations

TL;DR

This paper introduces a coupled simulation method combining molecular dynamics, finite element, and particle-in-cell techniques to study metal surface responses under high electric fields, revealing cyclic thermal runaway behaviors.

Contribution

It presents a novel integrated simulation approach for metal surface response to electric fields, including 3D space-charge effects and comparisons of different emitter geometries.

Findings

Thermal runaway occurs in cycles for narrow emitters.

Intensive evaporation can ignite plasma arcs.

The method accurately captures complex surface phenomena.

Abstract

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.

Tables

Table 1. Table 1: Results of the measurements.

	wide	thin 1	thin 2	ArcPIC Timko et al. (2015)
$t_{evap}$ [ps]	$\mathrm{229\pm 7}$	$\mathrm{194\pm 4}$	$\mathrm{613\pm 10}$	—
$t_{cool}$ [ps]	$\mathrm{176\pm 18}$	$\mathrm{175\pm 20}$	$-$	—
$r_{Cu}$ [ps${}_{}{}^{-1}]$	$\mathrm{22\pm 2}$	$\mathrm{110\pm 5}$	$\mathrm{238\pm 12}$	404
$r_e$ [fs${}_{}{}^{-1}]$	$\mathrm{6.7\pm 0.2}$	$\mathrm{8.2\pm 0.6}$	$\mathrm{3.7\pm 0.2}$	27.1
$R_{Cu}$ [nm-2 ns${}_{}{}^{-1}]$	$\mathrm{220\pm 20}$	$\mathrm{1100\pm 50}$	$\mathrm{2380\pm 120}$	40.4
$R_e$ [nm-2 ps${}_{}{}^{-1}]$	$\mathrm{67\pm 2}$	$\mathrm{82\pm 6}$	$\mathrm{37\pm 2}$	2.7
$r_{Cu/e}$ $[\times {10}^{-3}]$	$\mathrm{3.2\pm 0.3}$	$\mathrm{13.5\pm 1.6}$	$\mathrm{64\pm 6}$	15

Equations

$$\vec{\nabla }\cdot (ϵ\vec{\nabla }\Phi )=-\rho \text{,}$$

$$\vec{\nabla }\Phi \cdot \vec{n}=E_0$$

$$\vec{\nabla }\Phi \cdot \vec{n}=0$$

$$\Phi =0$$

$$\dot{\vec{r}}$$

$$\dot{\vec{v}}$$

$${\vec{v}}_{k+0.5}$$

$${\vec{r}}_{k+1}$$

$${\vec{v}}_{k+1}$$

$$\rho (\vec{r})=\sum _i{q}_i\delta (\vec{r}-{\vec{r}}_i)\text{,}$$

$$f_i=\frac{q{w}_{sp}}{ϵ}\sum _j{N}_i({\vec{r}}_j)+{\int }_{{\Gamma }_1}{N}_i{E}_{0}d\Gamma \text{,}$$

$$n_{sp}=\frac{J_{e}A\Delta t}{e{w}_{sp}}\text{,}$$

$$n_{sp}^{{}^{\prime }}$$

$${\vec{v}}_0=\frac{q}{m}\vec{E}\Delta t\left(\frac{1}{2}+R\right)$$

$$\Delta \vec{r}={\vec{v}}_0\Delta tR\text{.}$$

$${\vec{F}}_L(\vec{r})=\frac{1}{2}q\vec{E}(\vec{r})\text{.}$$

$${\vec{F}}_C({\vec{r}}_i)=\frac{1}{4\pi {ϵ}_0}\sum _{j\ne i}\frac{q_i{q}_j}{r_{ij}}{\hat{r}}_{ij}\exp (-\xi r_{ij})\text{,}$$

$${\oint }_{\Gamma }\text{vv}E\cdot \text{vv}dA=\frac{Q}{{ϵ}_0}\text{,}$$

$$Q_i={ϵ}_0{E}_i{A}_i\text{, }i=1,2,...,n_{faces}\text{.}$$

$$q_j=\sum _i{w}_{ij}{Q}_i\text{, }j=1,2,...,n_{atoms}\text{.}$$

$$\sum _i{w}_{ij}=1\forall j\text{.}$$

$$w_{ij}=\frac{\exp (-|{\vec{r}}_i-{\vec{r}}_j|/r_c)}{{\sum }_j\exp (-|{\vec{r}}_i-{\vec{r}}_j|/r_c)}\text{,}$$

$$q_j={ϵ}_0{\vec{E}}_j\cdot \sum _i{\vec{A}}_i\text{,}$$

$$P_J(T)=\sigma (T)(\nabla \Phi {)}^2\text{,}$$

$$\vec{\nabla }\cdot (\sigma \vec{\nabla }\Phi )=0\text{,}$$

$$\sigma \vec{\nabla }\Phi =\overrightarrow{J_e}$$

$$\sigma \vec{\nabla }\Phi =0$$

$$\Phi =0$$

$$\vec{\nabla }\cdot (\kappa \vec{\nabla }T)+P_J=C_v\frac{\partial T}{\partial t}\text{,}$$

$$T(t=t_0)=T_{amb}$$

$$\kappa \vec{\nabla }T\cdot \vec{n}=P_N$$

$$\kappa \vec{\nabla }T\cdot \vec{n}=0$$

$$T=T_{amb}$$

$$\kappa (T)=LT\sigma (T)\text{,}$$

$${\int }_{\Omega }w\vec{\nabla }\cdot (ϵ\vec{\nabla }\Phi )d\Omega +{\int }_{\Omega }w\rho d\Omega =0\text{.}$$

Full text

Keywords: vacuum arc, vacuum breakdown, electron emission, plasma initiation, thermal runaway, metal nanotip, multiscale-multiphysics simulation

I Introduction

Field-emitting tips play a detrimental role in various vacuum devices that require high electric fields, such as vacuum interrupters Slade (2007), x-ray sources KNAW (2019), fusion reactors McCracken (1980), and particle accelerators CERN (2018). When a field emitter is subjected to currents exceeding a certain threshold, the emission becomes unstable and swiftly results in a vacuum arc Dyke and Trolan (1953); Dyke et al. (1953). 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 Hartmann and Gundersen (1988); Dyke and Trolan (1953).

Although vacuum arcs have been studied more than half a century and it is well-known that they appear after intense field electron emission Dyke et al. (1953); Anders (2008), 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 (Anders, 2008; Mesyats, 1993, 2005). 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 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 of 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 Timko et al. (2015), plasma formation in Cu cathodes was studied by means of particle-in-cell simulations. It was found that plasma can build up in the vicinity of an intensively field-emitting cathode, assuming that not only electrons, but also neutral atoms are emitted from the cathode surface. Moreover, it was shown that it is sufficient to supply about 0.015 neutrals per electron to initiate the ionization avalanche leading to the formation of a stable plasma. Recent multiscale atomistic simulations on Cu nanotips Kyritsakis et al. (2018) revealed a thermal runaway process which causes field-assisted evaporation of Cu atoms and nanoclusters at a rate exceeding the threshold found in Timko et al. (2015). Thus, combining these two independently obtained results, one can consider intensively field-emitting tips as plausible sources of plasma, since they can supply both species, electrons and neutrals, at individual rates.

However, it is not clear yet whether a field-emitting tip is indeed able to produce sufficient densities to initiate the ionization avalanche, since there is no technique which can address this issue concurrently. To answer this question, a fully coupled molecular dynamics (MD) and particle-in-cell (PIC) simulation technique must be developed. The MD part will describe the behavior of atoms under given conditions, such as heat generated by field emission currents and stress exerted by the field, while the PIC part will follow the dynamics of the particles above the surface. In this way, the space-charge buildup above the field-emitting tip can also be found directly, without a need for simplifying assumptions.

In Kyritsakis et al. (2018), we attempted to take into account the space-charge effect by using a semiempirical approximate model, as proposed in Forbes (2008). This model treats the emitter as planar, using a universal correction factor to account for the three-dimensional (3D) nature of the geometry, which was an external adjustable parameter for our simulations, with its value roughly approximated based on previous PIC calculations Uimanov (2011). Although we obtained promising results with this model, the semiempirical correction does not capture the complex 3D distribution of the space-charge developed around the dynamically varying tip geometry. Moreover, the accuracy of the calculations of the space-charge affects directly the value of the field emission currents and all the subsequent processes developing in the tip due to the field emission process.

In this paper, we present an approach of the coupled MD-PIC technique for a more accurate estimation of processes contributing to plasma onset. In particular, we develop a method that directly couples PIC and MD simulations of the metal surface response to the field aided by the use of the finite element method (FEM). We use the model to simulate similar process as in Kyritsakis et al. (2018), but including accurately the 3D space-charge distribution that builds up above the highly curved nanometric field-emitting tip. There are some attempts in the literature to simulate the evolution of the space-charge using an MD-like method Torfason et al. (2015, 2016) where all interelectron interactions are considered. For specific applications, such models may be sufficient, but for computational efficiency we prefer to use the PIC method. The latter provides a sufficiently reliable kinetic description of space-charge or plasma evolution and allows incorporating the complex nature of the atomic processes and plasma-surface interactions, while maintaining a feasible computational load as compared to the MD-like methods Tskhakaya et al. (2007); Timko et al. (2015). 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 structures observed in experiments Timko (2011). 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

We present here a model that combines classical MD, FEM, PIC and electron emission calculations (fig. 1). The model allows us 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 model is provided in appendix A.

The model is based upon the previously developed framework FEMOCS Veske et al. (2018). In FEMOCS, an unstructured mesh is dynamically built around surface atoms, 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. For more details about meshing, local solution extraction, and optimization features in FEMOCS, see Veske et al. (2018). In the following sections we provide more details about the newly developed methodology for calculating electric field, electrostatic forces, and velocity perturbation.

II.1 Electric field

A common way to take into account an electric field within a computational model is to define it as the negative gradient of the electrostatic potential $\Phi$ that, in turn, can be calculated by solving the Poisson’s equation in the vacuum (domain $\Omega_{1}$ in fig. 2)

[TABLE]

where $\epsilon$ is the dielectric constant and $\rho$ is the space-charge density. We solve eq. (1) iteratively with the boundary conditions (BC) defined in the vacuum (see fig. 2) as

[TABLE]

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.

II.2 PIC method

To obtain the space-charge distribution $\rho(\vec{r})$ in the vacuum above the metal surface, which enters eq. (1), we use the PIC method Birdsall and Langdon (1985); Buneman (1959); Dawson (1962). Within PIC, the space-charge is formed by a set of superparticles (SPs), which is a common approach for space-charge and plasma simulations. Each SP is a group of charges of the same type, electrons or ions, that act within the model as a single particle. In other words, the SPs represent a sampling of the continuous phase space occupied by the real particles like electrons or ions. Such a grouping is justified by the fact that the electrostatic and electrodynamic interactions depend on the charge-to-mass ratio instead of the charge or the mass of the particles. The size of each 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.

The movement of SPs is affected by the global electric field, which includes the long-range electrostatic interactions between particles. The short-range interactions with nearby SPs are taken into account by means of binary Monte Carlo collisions (see appendix F for details). This makes it possible to significantly enhance the computational efficiency and simulate large systems for long times. Furthermore, the binary collision technique can be used to consider the interactions of different types of particle species in a plasma, such as impact ionizations. 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, which are summarized in fig. 3.

II.2.1 Particle mover

The propagation of SPs is performed by solving iteratively Newton’s equations of motion for each SP:

[TABLE]

where $\frac{q}{m}$ is the charge-to-mass ratio of SPs. In eq. (3), SPs are affected only by the electrostatic field, since the magnetic field which is induced by the currents within the tip is negligible. We solve eq. (3) numerically by using the leapfrog integration scheme Birdsall and Langdon (1985). In this method, the positions and the velocities of SPs are calculated as

[TABLE]

where the subscript indicates the time step, $Q=\frac{q}{m}\frac{\Delta t}{2}$, and $\Delta t$ is the PIC time step. Note that the calculation of $\vec{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 $\vec{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 $\Gamma_{2}$, will reappear at the symmetric point at the opposite boundary of the box. Such a mapping, however, occurs extremely rarely, as the simulation box is chosen wide enough to prevent particles reaching sides before hitting the top boundary $\Gamma_{1}$ where they are absorbed. For that reason, the contribution of mirror images of the simulation box to the total SP flux is negligible. 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.

Next to the top boundary $\Gamma_{1}$, the cathode boundary $\Gamma_{3}$ also acts as a particle absorber. On $\Gamma_{3}$, the absorbed particles may cause secondary effects, such as secondary and higher order emission, plasma recycling, impurity sputtering, etc. Since these effects 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 Brooks and Naujoks (2000); Taccogna et al. (2004a, 2007, b, 2005).

II.2.2 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) by solving a matrix equation (31), an inverse operation needs to be performed since the SP charge is required in the right-hand side (RHS) $f_{i}$ of eq. (31). 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 a form of charge density $\rho$. To overcome this issue, we represent the charge of each SP by means of the $\delta$ function, so that

[TABLE]

where $\delta(\vec{r}-\vec{r}_{i})$ is the Dirac $\delta$ function and the summation goes over the SPs that are located in the finite element where $\rho$ is evaluated. By doing such a substitution, the RHS of eq. (31) becomes

[TABLE]

where $\vec{r}_{j}$ are the coordinates of the $j$th 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.

II.2.3 Particle injection

In our simulations, the space-charge is built up due to intensive electron emission. Using the field emission tool GETELEC Kyritsakis and Djurabekova (2017), we calculate the current density $J_{e}$ in the centroid of each quadrangle that is located at the metal-vacuum 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:

[TABLE]

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\in[0,1]$ to decide whether to round $n_{sp}$ up or down:

[TABLE]

This discretization scheme guarantees that on average, the number of emitted SPs remains the same as given by eq. (7). 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 $\vec{E}=\vec{E}(\vec{r})$ is used to give them an initial velocity

[TABLE]

and displacement

[TABLE]

We add a term $\frac{1}{2}$ in eq. (9) to accelerate the particle by half a time step and initialize the velocity for proper leapfrog stepping. It is needed, because to enhance computational efficiency, we inject SPs between particle mover steps of (4a) and (4b) that were described in sec. II.2.1. The additional random number $R$ provides a fractional push by a small random fraction of the previous time step, preventing the particles from forming artificial bunches Timko et al. (2015). This small initial push enhances computational stability, without significantly affecting the results, as shown in sec. III.1.

II.2.4 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 $\mathcal{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 Takizuka and Abe (1977). More details about it can be found from appendix F.

II.3 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 is included as a force perturbation. In this way the total force on an atom becomes $\vec{F}=\vec{F}_{EAM}+\vec{F}_{L}+\vec{F}_{C}$, where $\vec{F}_{EAM}$ is the force obtained from the interatomic potential. The main perturbation is due to the Lorentz force

[TABLE]

In this equation, the term $\frac{1}{2}$ originates from the Maxwell stress tensor Landau and Lifshitz (1984) and takes into account that only half of the charge is exposed to the field; the rest is shadowed by the material domain. A less significant contribution to the total force $\vec{F}$ comes from the Coulomb interactions between the surface charges

[TABLE]

where $q_{i}$ is the charge of $i$th atom, $r_{ij}$ the distance between the $i$th and $j$th atoms, and $\hat{r}_{ij}$ the unit vector in the direction of $\vec{r}_{ij}$. The exponent in eq. (12) describes the screening of the $i$th charge by the conduction electrons in metals and the value of the screening parameter $\xi$ depends on the material and its crystallographic orientation. In our case, $\xi(\text{Cu})=0.6809$ Å-1 as defined in Djurabekova et al. (2011).

The surface charge can be estimated according to Gauss’s law Djurabekova et al. (2011)

[TABLE]

where $Q$ is the charge inside a closed surface $\Gamma$. 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 Djurabekova et al. (2011); Kyritsakis et al. (2018), 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. First, we calculate the charge of the mesh faces in the metal-to-vacuum boundary, hereafter referred to as face-charge:

[TABLE]

Notice, that in eq. (14) it is enough to use the scalar field and area, since $\vec{E}\parallel\vec{A}$ due to the Dirichlet BC. By using a weight function $w_{ij}\equiv w(\vec{r}_{i}-\vec{r}_{j})$ (not to be confused with the PIC superparticle weight $w_{sp}$), the face-charge is distributed between the atoms as

[TABLE]

The total charge must remain conserved, therefore the weight function must satisfy

[TABLE]

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:

[TABLE]

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. 4. 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:

Copy surface atom as a new point $\vec{p}$; 2. 2.

Locate the triangle $f$ where the projection of $\vec{p}$, in direction of $f$ normal $\vec{n}$, lies within $f$; 3. 3.

Move $\vec{p}$ in the direction of $\vec{n}$ for distance $\lambda$; 4. 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 $\lambda$ equals one lattice constant of the atomistic system. The coordinates of the surface atoms and the support points are input in the TetGen Si (2015) to generate a Voronoi tesselation around these atoms. Note that the bulk atoms are in fig. 4 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 the field at the atom position, $\vec{E}_{i}\approx\vec{E}_{j}$, we approximate eq. (13) for atoms in the regions with increased roughness as

[TABLE]

where the summation goes over the Voronoi facets of the $j$th atom that are exposed to vacuum.

II.4 Heating

We have shown previously Parviainen et al. (2011); Eimre et al. (2015); Kyritsakis et al. (2018) that thermal effects caused by field emission play a significant role in the evolution of nanotips under a high field. To take these effects into account, we calculate the electron emission current density $J_{e}$ and the Nottingham heat $P_{N}$ Charbonnier et al. (1964); Paulini et al. (1993) on the emitting surface by using our field emission tool GETELEC Kyritsakis and Djurabekova (2017). To calculate the volumetric resistive heating power density, we use Joule’s law as

[TABLE]

where $\sigma(T)$ is the electric conductivity and $\Phi$ is the electric potential obtained by solving the continuity equation in a metal (domain $\Omega_{2}$ in fig. 2)

[TABLE]

which comes with the following boundary conditions:

[TABLE]

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 (domain $\Omega_{2}$)

[TABLE]

where $C_{v}$ is the volumetric heat capacity. The initial and boundary conditions for eq. (22) are

[TABLE]

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), $\kappa=\kappa(T)$ is the heat conductivity as given by the Wiedemann-Franz law Franz and Wiedemann (1853)

[TABLE]

where $L$ is the Lorentz number (see the next paragraph). Equations (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 Kyritsakis et al. (2018), we took into account the finite-size effects (FSEs) by modifying the bulk value of $\sigma(T)$ Schuster et al. (2001) by a correction factor $\nu=\nu(T,d)$, as calculated by a simulation method developed by Yarimbiyik et al. Yarimbiyik et al. (2006). The characteristic size of the nanotip $d$ equals 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\cdot 10^{-8}$ W $\Omega$ K-2, as reported for a Cu film of 40 nm thickness Nath and Chopra (1974).

The resulting temperature distribution is used to adjust the atomistic velocities in MD. For this, the atoms are first grouped according to mesh tetrahedra that surround them. After this, a Berendsen thermostat Berendsen et al. (1984) is applied separately to each such group. The target temperature $T$ of the thermostat is obtained by averaging nodal temperatures of a tetrahedron, while the macroscopic temperature of the atomistic domain $T_{0}$ equals the average of microscopic temperatures of all the atoms of the same group. For the control, we use a time constant of $\tau=1.5$ ps. Such $\tau$ has shown previously Kyritsakis et al. (2018) to prevent the simulation artifacts in MD, while being significantly lower than the relaxation time of the heat equation.

III Results and discussion

III.1 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 Veske et al. (2018) and Kyritsakis and Djurabekova (2017).

III.1.1 Validation of PIC

We validate our PIC model by calculating the current-voltage curve for a planar cathode-anode system, for which a 1D semianalytical solution of Poisson’s equation is available Child (1911); Barbour et al. (1953); Forbes (2008). For this, we apply a voltage $V$ between two parallel plane electrodes of area $A=135$ nm2 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 Kyritsakis et al. (2018). For high voltages, this curve asymptotically converges to the Child-Langmuir law Child (1911); Langmuir (1913, 1923).

The results of the PIC simulation, together with the semianalytical and the Child-Langmuir curves, are shown in fig. 5. We observe that our PIC model is in very good agreement with the semianalytical curve with an RMS error of 1.9%.

III.1.2 Validation of surface charge model

To validate the surface charge calculation method, we compare it against our earlier model HELMOD Djurabekova et al. (2011). 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 $\left<100\right>$ Cu block.

We measure the surface charge both with HELMOD and with the model presented here and calculate the relative difference $\epsilon$ between them. The variation of $\epsilon$ along the surface of the $\{100\}$ and $\{110\}$ slice of the apex is presented in fig. 6. The data show that the current model tends to provide slightly smaller charges than HELMOD with a mean difference of $6.5\pm 1$%. The largest deviation occurs on sites that correspond to atoms below atomic steps and kinks. This indicates that Voronoi cells that are built around those sites tend to be poorly exposed to vacuum, and for that reason, a smaller surface area is assigned for those sites as compared with HELMOD. From the $\{110\}$ slice it is also visible that $\{110\}$ surfaces are a source of systematic difference. This occurs because the Voronoi cell-based model tends to assign charge also for atoms on the second layer of the $\{110\}$ surface, while HELMOD considers those atoms as totally shadowed.

III.2 Simulation setup

We used the model described in sec. II to simulate the behavior of a Cu nanotip under a 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, which 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 that defines the vacuum domain $\Omega_{1}$ and the metal domain $\Omega_{2}$ is illustrated in fig. 2. In $\Omega_{1}$, PIC and electric field calculations are performed, while MD and the heat solver operate in $\Omega_{2}$. Notice, that MD uses only the upper part of $\Omega_{2}$ (enclosed with dashed lines in fig. 2), while FEM covers all of it. A conical tip with an aperture angle $\gamma=3^{\circ}$, 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 atomistic region is cut out from a monocrystalline $\left<100\right>$ Cu box with a lattice constant of 3.64 Å. 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 of $r_{thin}=17$ nm and $r_{wide}=54$ nm.

The simulation geometry is motivated from the experimental observations that suggest the existence of field emitters with an aspect ratio of 20-100 Kildemo (2004) on flat Cu surfaces. The size of the tip is chosen to be large enough to prevent instability due to excess surface energy Veske et al. (2016); Jansson et al. (2016) and sufficiently small to result in a reasonable computational cost. In order to obtain results that are directly comparable to our previous work Kyritsakis et al. (2018), we chose the dimensions of the thin system to be identical to the ones used in Kyritsakis et al. (2018). In order to obtain sufficient mesh density, the characteristic distance between mesh nodes in the apex region is chosen to be one half of lattice constant of the initial atomistic system.

The MD simulations were carried out by means of the classical MD code PARCAS Nordlund (1994); Nordlund et al. (1998); Ghaly et al. (1999). The FEM calculations in FEMOCS are based on open-source C++ library Deal.II Bangerth et al. (2007). For the MD simulations we used the interatomic EAM potential developed by Mishin et al. Mishin et al. (2001). This potential has been successfully used in our previous works Kyritsakis et al. (2018); Veske et al. (2016); Parviainen et al. (2011) and has shown accurate reproduction of nonequilibrium system energetics Mishin et al. (2001). The stochastic nature of the thermal effects were taken into account by running 50 independent simulations with different random 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.

III.3 Course of the simulation

We simulated the runaway process in a nanotip similarly to the one in Kyritsakis et al. (2018). Despite the usage of a more simplistic 1D space-charge model in Kyritsakis et al. (2018), 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 duration of the entire process until the tip stops emitting neutrals.

In fig. 7, 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 of the tip. Normally this leads to an abrupt change in the tip height, but due to averaging this change is smoothed 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. 7, are provided on the right of fig. 7.

In fig. 8, 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 space-charge, 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. 7a). After the melting point was reached, the tips start growing, causing a step-by-step increase in the local field, current, and temperature (fig. 7b). The 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. 7d). In some cases, a smaller cluster of the size up to a few atoms might also be detached (fig. 7c). 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. 7e), 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. 7h).

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. 9, 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% Schuster et al. (2001). 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. 7-9 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. 7 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. 10. 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 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}=\frac{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. 10). 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 evaporation 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 1.

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. 8. 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 Timko et al. (2015) 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 Timko et al. (2015), the neutral evaporation rate $r_{Cu}$ was considered proportional to the electron emission rate $r_{e}$, with a minimum ratio of $r_{Cu/e}\equiv r_{Cu}/r_{e}$ of 0.015 required to lead to plasma formation. In our previous work Kyritsakis et al. (2018), this ratio was found to have an average value of 0.025$\pm$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. 11 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. In both 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. 11), the slope of which gives the evaporation-emission rate. The calculated parameters are summarized in table 1.

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 four to five 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.

III.4 Discussion

The 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 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}\approx 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}\leq 0.2$ GV/m Descoeudres et al. (2009).

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 Barbour et al. (1953); Uimanov (2011), which assumes a uniform current density in space. While simulating field emitters, however, the emission current in space is highly nonuniform due to the nonplanar 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 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 1, the ratio $r_{Cu/e}$ for the second runaway event is significantly higher than the threshold value of 0.015 that was determined in Timko et al. (2015). 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 $\left<r_{Cu/e}\right>=0.026\pm 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 Timko et al. (2015). However, it is premature to conclude that thick field emitters are not able to produce enough neutral atoms to ignite plasma. First, in the simulations performed by using the ArcPIC model in Timko et al. (2015), 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 Timko et al. (2015), a much wider electron emission and neutral injection area ($10^{4}$ nm2) was assumed. In the present simulations, we can safely consider that all emission and evaporation originate from an area smaller than 100 nm2 (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 1 we see that the values for $R_{Cu}$ and $R_{e}$ obtained from the present results significantly exceed the ones assumed in Timko et al. (2015), 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 Timko et al. (2015).

IV Conclusions

We have investigated the dynamic evolution of Cu field emitters of different widths 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 a lower neutral evaporation rate, which leads to the conclusion that very sharp field emitters may be necessary to enable plasma formation.

Acknowledgements

The current study was supported by the Academy of Finland project AMELIS (Grant No. 1269696), CERN CLIC K-contract (No. 47207461), Estonian Research Council Grant PUT 1372, and the Estonian national scholarship program Kristjan Jaak (No. 16-4.2/653), 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

Appendix B Discretization of Poisson’s equation

Multiplying eq. (1) with a weight function $w$ and integrating over the whole domain gives

[TABLE]

Using the properties of differentials, the previous expression could be rewritten as

[TABLE]

The usage of Green’s theorem on the first term removes double differentiation of $\Phi$:

[TABLE]

By plugging BCs (2) into leftmost term, we obtain the weak form of eq. (1):

[TABLE]

To transfer that formula into FEM format, we use the Galerkin approach and equalize the weight function with shape functions and discretize the potential as a linear combination of them:

[TABLE]

where $n$ is the number of degrees of freedom of the finite element. By also combining $\epsilon$ and $\rho$ together, (28) becomes

[TABLE]

To solve that system of equations numerically, it is handy to express it in matrix form

[TABLE]

where

[TABLE]

Appendix C Discretization of continuity equation

Notice the similarity between eq. (20) and (1). Therefore, by replacing dielectric constant $\epsilon$ with conductivity $\sigma$ and by taking $\rho=0$, we obtain the representative form of eq. (27) for (20):

[TABLE]

The application of BCs (21) into the second term gives the weak form of eq. (20):

[TABLE]

To transfer that formula into algebraic form, we use the same strategy as in appendix B and equalize weight function with shape functions and expand the potential as a linear combination of them. In that way the eq. (33) becomes

[TABLE]

which looks in matrix form as

[TABLE]

where

[TABLE]

Appendix D Discretization of heat equation

To discretize the heat equation (22), we use similar strategy as in appendix B and multiply (22) with a weight function $w$ and integrate over the whole domain:

[TABLE]

The properties of differentials allow expressing (36) as

[TABLE]

By means of Green’s theorem, we remove the double differentiation of $T$ from the first term of (37),

[TABLE]

After plugging BCs (23) into the previous formula and rearranging terms, we obtain the weak form of eq. (22)

[TABLE]

To express the weak form as a matrix equation, we follow the procedure from appendix B and equalize the weight function with shape functions and expand the temperature as a linear combination of them. In that way we separate spatial and temporal parts in $T$, and eq. (39) becomes

[TABLE]

which looks in matrix form as

[TABLE]

where

[TABLE]

To solve semidiscrete eq. (41), the time also must be discretized. It can be done by first approximating

[TABLE]

and then introducing parameter $\Theta\in[0,1]$, so that

[TABLE]

After combining (42) and (43) together, we get a $\Theta$-scheme for solving eq. (41):

[TABLE]

To handle the nonlinearity, we assume in this formula, that thermal and electric conductivities are weakly temperature dependent, i.e., $\kappa(t_{n+1})\approx\kappa(t_{n})$ and $\sigma(t_{n+1})\approx\sigma(t_{n})$. In the case of $\Theta=0$, $\Theta=0.5$, and $\Theta=1$, we get the explicit Euler, Crank-Nicolson and implicit Euler scheme, respectively. In our simulations we use the 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 the mesh rebuild. Although the implicit Euler scheme shows a lower order of convergence than the Crank-Nicolson method, i.e., smaller $\Delta t$ needs to be used Lewis et al. (2004), 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. 12a). The coordinates of random point inside a parallelogram $ADEF$ are given by

[TABLE]

where $R_{1},R_{2}\in[0,1]$ are uniform random numbers. The point $\vec{r}$, however, might not lie inside the quadrangle $ADOF$, as its area is smaller than the area of the parallelogram, $S_{ADOF}=\frac{1}{3}S_{ABC}<S_{ADEF}=\frac{1}{2}S_{ABC}$. Therefore its location in the quadrangle must be checked, and if $\vec{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 Vince (2014). Denoting $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ as the barycentric coordinates of point $\vec{r}$ with respect to the first, second and third node of a triangle (fig. 12b), $\vec{r}$ is located in the left region of the triangle (vertical hatching) if $\lambda_{1}\geq\lambda_{2}$ and in the bottom diagonal region (diagonal hatching), if $\lambda_{1}\geq\lambda_{3}$.

Appendix F Collisions between superparticles of the same type

The Monte Carlo binary collision model can be described in the following steps:

Divide SPs into groups in such a way that all particles within the group are located in the same cell. 2. 2.

Pair the SPs in the group in a random way. 3. 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 Takizuka and Abe (1977). To collide elastically SP pairs, we use the statistical model as follows. By assuming energy and momentum conservation, the velocities of paired particles of identical mass $m$ and charge $q$ before and after the collision, $\vec{u}_{1,2}$ and $\vec{v}_{1,2}$, are written as

[TABLE]

The velocity change $\Delta\vec{v}$ is calculated from the scattering of the relative velocity,

[TABLE]

where the rotation matrix $\mathbf{M}$ is expressed as

[TABLE]

where $a=\sin\theta\sin\phi\frac{v}{v_{\perp}}$, $b_{x,y}=\sin\theta\cos\phi\frac{v_{x,y}}{v_{\perp}}$ and $v_{\perp}=\sqrt{v_{x}^{2}+v_{y}^{2}}$. The Monte Carlo technique comes into play while calculating the azimuth and scattering angle $\phi$ and $\theta$. To calculate the first, we generate a uniform random number in the range of $[0,2\pi)$. Scattering angles $\theta$, however, need to be distributed normally and their contribution to matrix (48) are found from

[TABLE]

where $\delta\equiv\tan\theta$ is a random number generated according to the Gauss distribution. The mean of that distribution $\langle\delta\rangle=0$ and the variance for a system of SPs with identical charge and mass

[TABLE]

where $n$ is the number of particles in the cell, $V_{cell}$ the cell volume, $\epsilon_{0}$ the vacuum permittivity and $\Lambda=13.0$ the Landau logarithm Krall and Trivelpiece (1973). In the case that there is odd number of SPs in the cell, the first three pairs in that cell should be collided with a variance of $\langle\delta^{2}\rangle^{{}^{\prime}}=0.5\langle\delta^{2}\rangle$ Takizuka and Abe (1977).