Atomistic modeling of metal surfaces under electric fields : Direct coupling of electric fields to a molecular dynamics algorithm

The effect of electric fields on metal surfaces is fairly well studied, resulting in numerous analytical models developed to understand the mechanisms of ionization of surface atoms observed at very high electric fields, as well as the general behavior of a metal surface in this condition. However, the derivation of analytical models does not include explicitly the structural properties of metals, missing the link between the instantaneous effects owing to the applied field and the consequent response observed in the metal surface as a result of an extended application of an electric field. In the present work, we have developed a concurrent electrodynamic– molecular dynamic model for the dynamical simulation of an electric-field effect and subsequent modification of a metal surface in the framework of an atomistic molecular dynamics (MD) approach. The partial charge induced on the surface atoms by the electric field is assessed by applying the classical Gauss law. The electric forces acting on the partially charged surface atoms (Lorentz and Coulomb) are then introduced in the MD algorithm to correct the atomic motion in response to the applied field. The enhancement factor at sharp features on the surface for the electric field and the assessment of atomic charges are discussed. The results obtained by the present model compare well with the experimental and density-functional theory results.


I. INTRODUCTION
Many different computational methods have been developed to investigate the effects of modification of materials caused by ion beam [1,2] or laser impacts [3][4][5] on any surface. Even though many experiments showed that significant modification of the material can be caused by applied electric fields [6][7][8][9], there is no plausible technique for dynamic simulation that could predict the possible modification of the structure when an electric field is applied. However, the need for such a technique is increasing with the miniaturization of electric devices, or increasing magnitude of the operated electric fields in large-scale machines such as linear accelerators with the sophisticated design of the accelerating structures [10]. A number of analytical models exist to explain the phenomenon of the field ion microscopy [11][12][13], but the gap between experimental observations and the theoretical predictions remains unbridged. It has also been proved that the probability of the vacuum arcs to occur near metal surfaces in ultrahigh vacuum conditions depends on the structural properties of the metals [14]. Thus knowing the actual dynamics of the motion of surface atoms of a certain metal under a high electric field can shed light on the triggering process of vacuum arcs.
The main reason for the atoms to behave differently under electric fields is the shift in electronic densities of surface atoms caused by an external electric field. The higher the field, the larger the shift is expected. Without discussing the nature of excess or depletion of an electronic density near surface atoms, we will use the term "atomic charge" to describe the state of the partial ionization of an atom. Note that the value of this charge is not an integer number of electrons and can well be less than 1.
There is a previous implementation of the surface charge into molecular dynamics (MD) algorithms [15], but the calculation of the charge on an atom is not done in correspondence with Gauss law. Instead, some atoms obtain a full integer charge, while others remain neutral although they are exposed to the field. Moreover, the given charges in that work do not change dynamically even if the surface shape changes strongly. However, the controllable operation of many tools designed for use in electric fields will be strongly dependent on the accuracy of the predictive model. The process of atom detachment from the surface under high electric fields is extensively debated in field evaporation theory. The models that exist to explain the experimental observation of the atom detachment under the high electric field [16][17][18][19] are based on the generalized approach to give a qualitative description of the interaction of a metal surface with a n-charged ion using a standard potential [13,20]. Such an approach practically cannot guarantee a reliable prediction on the surface breakage for a particular metal. Moreover, these models experience difficulties in distinguishing the two clearly different states of the atom: the strongly bound state of the atom, which belongs to the surface (with a maximum number of neighbors in its vicinity), and the weakly bound state of the atom, which is located above the surface (an adatom). These cases can be distinguished only if the real interatomic interactions are taken into account. The first-principles calculations of the external field effect on a metal surface [12,21,22] are able to account for the properties of the particular metal, however, they fail to encompass the dynamic evolution and complex features on the surface, which develop under the field. Nevertheless, these types of calculations can give a reliable estimation of the charge distribution owing to the rigorous account of the superposition of the external and internal electric fields.
In this paper we present a concurrent electrodynamicsmolecular dynamics (ED-MD) model for a real-time simulation of the development of arbitrarily shaped surfaces under electric fields; some results on charge distribution on surface atoms owing to the electric field will be discussed. In the presence of a sharp feature on the surface, electronic effects cannot be ignored. Even if the field applied to the metal surface is much lower than the threshold of the field-emission process, the field can grow to this value above the top of a tip. In this case the electric current caused by the field-emission phenomenon from the tip will raise its temperature by Joule heating, competing with the conduction of heat to the rest of the cell. These processes are also introduced in our model, but are discussed in greater detail elsewhere [23].

II. METHOD
In principle, the effects of an external electric field on the electrons and nuclei of a metal surface should be calculated by a quantum-mechanical solution of the electronic states of the system. Such calculations have been done until now for the static surface configuration and small atomic systems only [12,21,22]. However, to enable the dynamic simulation of morphology evolution on extended surfaces, we introduce an atomic-level interpolation approach to the application of the classical Gauss law [24]. The effect of the electric field on the surface atoms that we consider is through the effect on the charges induced on these atoms (a result of the depletion or excess of electronic densities according to Gauss law).
We developed our ED-MD model on the basis of the classical MD code PARCAS [25]. The main idea is to introduce, into the Newton's equations of atomic motion, additional forces of an electric nature that act on the charges induced by the applied electric field. Within the Born-Oppenheimer approximation (which is conventionally used in MD simulations [26]), we presume that the relaxation of the electronic system in response to the electric field and the subsequent Coulombic attraction of nuclei by electrons are momentaneous compared to the thermal motion of the atomic nuclei. Thus, the atom is considered as a whole entity with respect to the external electric field.
We consider further that the Gauss law stated for a continuous metal surface (given by the "pillbox" technique, which is a surface charge per area, see Sec. II C [24]) is applicable on an atomic level. We do not account for electron tunneling effects responsible for the ionization of surface atoms at the moment of their detachment from the surface. The process of detachment is considered if a surface atom (or an adatom) is provided sufficient energy from the electric field to overcome the surface barrier. When an atom resides at a metal surface, it can be partially charged, i.e., have an electron density that is higher or lower than that of an average atom inside the (uncharged) bulk in a way that fulfills Gauss law. On the other hand, when an atom is leaving the surface, the electron density must at some distance from the surface collapse into a state that corresponds to an integer number of electrons on the detached atom. We also do not account for this loss or gain of energy owing to the electron wave collapsing during the detachment process; this is why there is no additional energy or energy penalty for these processes accounted for in the model.
Thus, the degree of ionization of the atom cannot be considered reliably predicted within the present approach. Note also that the present approach has been developed for the electrostatic field; no magnetic components or time dependences of a rf field as well as its characteristically significant skin depth are taken into account. The electric field is not allowed to penetrate deeper than the first atomic layer, which correspond to the skin depth ∼2Å, as estimated for metal surfaces in the classical electrodynamics textbooks [24], although the boundary between the internal field of the surface atoms and the external electric field is considered to be sharp.
The atoms leaving the surface in the ionized state, or those that experienced a postionization, as well as the strong current of electrons escaping from the surface, affect the field around the tip. Strictly speaking, to find the field distribution around the tip at the fields of interest of this paper, one must solve the Poisson equation with the space-charge component in it [24]. Somewhat simplified approaches (the case of planar electrodes in Ref. [27]) can help to take into account the space-charge screening effect by scaling the value of an electric field. However, a field screening effect owing to a space charge will not be included in our approach owing to the reasons specified in the next section.

A. Calculation of an electric field above the metal surface of an arbitrary shape
We find the distribution of the electric field over the metal surface of an arbitrarily rough shape (not flat) by solving the Laplace equation where is the electrostatic potential [24]. The motivation for the employment of the Laplace equation rather than the general case of the Poisson equation [24] for electric-field distribution, accounting for a space charge, is that that the primary interest of the present research is focused on the processes within the metal surface. The charges induced on a conductive surface are maintained by the voltage from the external source. These define the electric field between electrodes that will be distorted in the space owing to the roughness on the surface of the electrodes, or a space charge growing with the electron current from a surface emitter. These emitters, in the shape of surface protrusions, will distort the field efficiently, causing the redistribution of the charge on the surface atoms of the protrusions and in their vicinity. We aim to find this redistribution by applying the Laplace equation. The field enhancement on a protrusion will be screened eventually by the space charge over the top of the protrusion, but in the first approximation we can neglect this effect, which can lower the value of a charge by ∼ 20%, as estimated in Ref. [27], which the current approach can tolerate because of the other uncertainties caused by the assumptions discussed in the following sections.
For the solution of Eq. (1) we employ the fine threedimensional (3D) grid and apply the finite-difference method using Gauss-Seidel iterations [28]. The size of a grid point is the same order of magnitude as the size of a lattice atom, therefore one atom belongs in one grid point in such a way that there is at least one grid point between two neighbor atoms if the atoms occupy the lattice sites according to the perfect crystal structure. This algorithm was developed for {100} planes in fcc lattice structures and can be also applied to {110} planes of bcc lattice structures. The grid points of the Laplace solver are not required to be strict cubes, but they have measurements proportional to the lattice constants of a metal of interest in all three dimensions.
We solve Eq. (1) with the mixed boundary conditions; the Dirichlet type at the bottom of the grid [ (x,y,0) = const on a metal surface and the constant can be assumed zero for certainty] and the Neumann type at the top of the grid ( ∇ = − F ex ) [28]. In infinity the field is aligned parallel, thus only the z component of the potential has a gradient; this is equal to the value of the external electric field. The upper boundary is written then as ∂ ∂ẑ = − F ex . The bottom boundary has a complex shape as it is fully defined by the position of surface atoms. Every atom that belongs to the surface (an atom with at least one bond to the rest of the cell) is considered as a bottom boundary and so are the empty grid points between the atoms. If the surface has a complex shape with sharp features on it, must be considered as a function of the radius vectors of the surface atoms positions, in other words, the bottom boundary condition for the Laplace solver is (x,y,z bottom ) = 0, where z bottom is the z coordinate of the bottom atoms. Because the shape of the surface is dynamically changing during the MD simulation, the time dependence is implicitly included in (x,y,z bottom ) function and is taken into account via the replacement of the bottom boundary shape every time the Laplace solver is called to find the electric-field distribution. In this manner the bottom boundary is dynamically changing with time, depending on the momentary positions of surface atoms.
For the best performance of the Laplace solver, we require the convergence criterion x,y,z ( (x,y,z) i ) − (x,y,z) i−1 ) 10 −12 V, where the indices i and i − 1 denote the current and previous iteration steps, respectively, for calculating the potential (x,y,z). The efficiency of the finitedifference algorithm is, however, fairly low. The iterations over all the grid points at every MD step demand extensive computational time, especially if the number of atoms that form the surface is large. While solving the Laplace equation at every MD step, the efficiency can be significantly improved by exploiting the current solution as an initial guess for the next step. In this way, if the positions of the atoms did not change so much that the atoms remain within the same grid points during one MD step, the iterations will not be needed as the convergence will be reached already at the very first step. The other trick is to increase the efficiency of a Laplace solver to apply more advanced iteration techniques, which can save significantly the computational time, such as a multigrid method.

B. The multigrid method
Solving Laplace's equation using the Gauss-Seidel iteration requires calculating n 3 finite differences for every iteration, assuming a cubic grid with n points in each direction. If a large accuracy is desired, the iteration must be performed over thousands of steps which can be very time consuming for large systems.
However, by using the multigrid method, the process of solving Laplace's equation is accelerated significantly. A comprehensive explanation of the application of the multigrid method can be found in Ref. [28]. In practice, we solve the equation at first approximately by using a few iterations of the Gauss-Seidel algorithm, and calculating a correction term, which improves upon the initial solution. This correction term is solved on a coarser grid with half the number of grid points in each direction, i.e., 1 8 of the total number of points. We use the full-weighting restriction and linear interpolation to transfer values between grids of different coarseness.
The determination of the correction term is recursive, i.e., a correction term to the correction term is calculated by using a coarser grid, and this is repeated until no coarser grid can be formed. At the last solution on the coarsest grid, we apply the plain Gauss-Seidel iteration. For increased efficiency, movement occurs back and forth between finer and coarser levels in a W shape, i.e., the solver returns to a coarser level after temporarily moving back up to a finer grid level (known as a W cycle).
We also note that the efficiency of the multigrid method strongly depends on the number of grid points in one direction on the original grid. If we write the equation for the number of grid points as n = a2 b , where a and b are integers, it is clear that a points will remain on the coarsest grid, because the number of points is halved between grids. Thus, it is expected that no acceleration can be obtained if the size of the system is odd because the number of points will remain unchanged, while a size which is a power of 2 will lead to a maximum calculation speed.
The efficiency of the multigrid method implemented in our model is illustrated in Fig. 1. Here we compare the performance of the Laplace solver with various system sizes using a pure Gauss-Seidel iteration to that utilizing the multigrid method, while the accuracy of the calculation remains the same (<10 −12 V). The multigrid algorithm is orders of magnitude faster than the plain version, and the importance of grid size in the multigrid case is also clearly seen. Somewhat paradoxically, increasing the number of grid points can lead to a faster solution, e.g., the solution of a 124 × 124 × 124 system takes ten times longer to calculate than for a 128 × 128 × 128 system. In our simulations we thus always choose systems with dimensions close to a power of 2.

C. Calculation of a charge induced by the electric field on a surface atom
Nevertheless, even by having the fastest Laplace solver, we still encounter the inevitable problem that the discrete information obtained in such a Laplace solver must be transferred to the continuous atomic motion in the MD algorithm. This is one of the challenging parts while developing a hybrid ED-MD model. In our approach we require that two atoms cannot occupy a single grid point, but they are allowed to occupy as close as two neighboring grid points. With respect to the Laplace solver, the atoms are considered to be the same cuboid shape as a grid point of the Laplace solver; the center of the atom is always forced to the center of the grid point of the Laplace solver, where the atom was found. The real position of an atom (obtained in the MD algorithm) can be different, but this difference always remains within the uncertainty ±0.5a x,y,z , where a x , a y , and a z are the sizes of a grid point in all directions. In this fashion the charge of the surface atom as a function of its position with regard to the neighbors is slightly discrete, and still monotonically corresponds to the change of the field around the atom. However, the charge obtained in the discrete grid of the Laplace solver is eventually assigned to the atom in the real position according to the MD algorithm. This compromise allows for the agreeable conjunction of the discrete technique of the Laplace solver and the continuous atomic motion calculated by solving the system of Newton's equations [26].
By the charge of surface atoms owing to an applied electric field, we understand the following. If the field is applied between two planar electrodes, it will polarize the surface atoms, inducing either depletion (on an anode) or excess (on a cathode) of electron density around the atoms. The process of partial charging of surface atoms (polarization) can be understood in the framework of the jellium model of a metal, which has been highly exploited to calculate the metal properties in analytical and density-functional theory (DFT) approaches [11,13,21]. Nevertheless, in the present approach we do not discuss the nature of polarization of metal atoms, because we are not doing rigorous calculations of electronic states of the atoms. Presently we employ the Gauss law, which gives the linear dependence of surface charge density σ on a value of an applied electric field F in SI units [24] as Here, ε 0 = 8.85 . . . × 10 −12 F/m is the vacuum permittivity and | F | is the electric field directly over the metal surface. We also assume that the Gauss law is valid at the atomic level; then the charge induced on a surface atom will be where A at is an area of the surface atom exposed to the field. The topography of atom location in the {100} face of a fcc structure (Cu{100}) or in the {110} face of a bcc structure (W{110}) enables the representation of an atomic layer as a chessboard with filled (occupied by an atom) and empty (hold no atom) cuboids (Fig. 2). In other structures the surface pattern must be reconsidered, although the main algorithm will be still valid. In this figure the dark cuboids represent atoms that are forced to occupy the same position as the grid point of the Laplace solver where the center of the atom belongs (while in reality the atom can be only partly there). Thus we approximate every atom to a cuboid shape in such a way that every surface atom has four neighbors on a perfect lattice. We consider three Cartesian components of the electric field F i x,y,z in which the grid point i possesses from the solution of the Laplace equation (1). Every component F i x ,F i y ,F i z has an opportunity to contribute to the entire value of the atom charge, where A j yz,xz,xy is the area of the side of the atom j perpendicular to the component F i x,y,z , respectively. The total charge of the atom j will be calculated as a sum of all the contributions q i j from the grid points surrounding the atom j : where N j is the number of grid points with F i obtained as the Laplace solution. As can be seen from Fig. 2, in the case of a flat surface there is only one grid point, which is located directly above the atom. If we calculate the charge of the atoms taking into account only these grid points, the total charge of the surface will be underestimated by a factor of 2. At the same time, the grid points, which are not located directly on the top of the atom, cannot be neglected because the value of the field in those grid points is also significant. Thus weencounter another "continuous↔discrete conjunction" problem. In this case, the continuous distribution of the field should be applied FIG. 3. Two-dimensional projection of the atomic cell model with the example of the field direction due to the distortion on the surface feature. In the model for side atom 1 the contribution of the F x component is taken into account (F y is ignored) and for atom 2 the situation is inverse.
to the discrete number of surface atoms, which are approximated to have a geometrically finite shape. The conjunction between the continuous field and the discrete atoms is done as follows.
Every grid point is considered in the middle of eight other grid points that together form a large complete cuboid. All eight grid point cuboids are divided into three groups: The first nearest neighbors (nn) are on the faces of the large cuboid, the second nn's are on its edges, and the third nn's are in its corners. Only those grid points that potentially host the atoms (not parts of the Laplace solver) are considered. If the atoms are found, the calculation of the charge contribution to each atom is made according to the priorities; the field of the given grid point does contribute in the charge of each atom around, but not in full strength. It is divided into two halves; one half is for 1 nn atoms and the other half is evenly divided for second and third nn atoms. If the atoms are missing, the strength is distributed between the rest of the atoms with the same rule of priorities. The equations for the fractions that are accounted for for each atom found in first, second, and third nn's are written as follows Here N inn is the number of atoms found in the inn shell; g 0 inn and N 0 inn are the fraction of a grid point given to the atoms of the inn shell and the maximum number of atoms, which can be found in this shell in a perfect structure, respectively (g 0 1nn = 0.5, N 0 1nn = 1; g 0 2nn = 0.25, N 0 2nn = 4; g 0 3nn = 0.25, N 0 3nn = 4). In this fashion the sum of the charges assigned to the discrete atoms that belong to a flat surface is equal to the charge distributed on this surface according to Eq. (2).
Note also that the distortion of the field on sharp features on the surface can lead in this algorithm to the presence of a field component parallel to the metal surface, which contradicts the Gauss law (Fig. 3). This is an undesired inherent feature of the discrete grid of the Laplace solver. The finer the grid, the smaller the probability for this artifact to appear. Presently we ignore the parallel component of F i and take into consideration only that which is perpendicular to the metal surface.
The charge obtained in this way is used in our MD simulations to calculate the force exerted on the surface atom by the electric field and owing to the Coulombic interactions. Then, the total force on an atom in MD algorithm becomes Here, f EAM is the force obtained from the interatomic potential and f i,L = q i F i isthe Lorentz force exerted on the atom i by the electric field F i , which is the field obtained as a vector sum of all the grid points surrounding the atom i. The Coulombic force f i,C is calculated as follows: Here, the screening factor e −r ij ξ describes the screening of a bare charge by conduction electrons in a metal [29]. The value of ξ can be determined by comparing the classical electrodynamics result of the energy of a flat surface with the summation of the screened Coulomb potential for the atomic arrangement of a particular surface. The calculation is described in detail in Appendix A, and gave a value of ξ = 0.6809Å −1 for the Cu (100) surface. However, we note that for the high electric fields employed here, the Lorentz force generally dominates, and hence the results are not sensitive to the value of ξ as long as it remains in a reasonable range ∼0.1-1.0Å −1 .
The value of ξ is not easy to determine for a charged metal surface, however, we found that for large electric fields and reasonable choices of ξ ∼ 1/(1Å)-1/(10Å), the electric forces dominates over the Coulomb force and hence the choice of ξ does not affect the results significantly.
The sharp features assume the highest partial charge, which is, for soft metals, still a fraction of an electron, but for the hard melting metals (with s high cohesive energy of the surface atoms) the value of the partial charge on the surface atoms can exceed significantly a value of 1e before the surface breaks.

A. Laplace solution for protrusions on a metal surface
Prior to collecting results from the proposed model we tested the Laplace solver exploited in our approach. As it was mentioned in Sec. II A, to solve the Laplace equation we employ the finite-difference method with the mixed boundary condition. The accuracy of the numerical integration depends on the level of the grating of the basic grid of the solver, which in turn is restricted by the computational efficiency of the technique employed for the dynamic simulations. Although the finite-difference method is well established and widely used for the integration of the differential equations, the level of the grating of the basic grid for the solver can play a crucial role for the reliability of the results. Because within the present model we are dealing with a perfectly conductive material, some known solutions for the field distortion on hemispherecapped cylindrical protrusions can be tested for verification purposes. The solver can be verified easily by the testing the dependence of field enhancement | F max |/| F ex | on the aspect ratio of a surface protrusion. Here we use β to denote the field enhancement factor (in line with the common practice of using this parameter in the Fowler-Nordheim equation for the field-emission current density) and γ will denote the aspect ratio as in Refs. [30][31][32].
The generally accepted concept of the electric field linearly enhancing with the aspect ratio of a metal protrusion on a metal surface [33] was verified in many articles [30][31][32]34]. Although the different calculations give slightly different fitting functions to describe the dependence, the main conclusion of a significant growth of the field on a protrusion seems to be a strong tendency in all models. Although the slope of the linear dependence is high, it is less than the unity as given in Ref. [33]: The analytical model by Kokkorakis et al. [34] gives the simplified formula| F max |/| F ex | = a 0 + 0.72r tip /h tip , and those based on finite-element calculations Edgcombe and Valdrè [30] derive a similar but sublinear dependence. These formulas give approximately the same slope of 0.72, which can be a good probe testing for our Laplace solver. For this purpose we used the same solver as implemented in our model, but without an atomic system, to be able to check the effect of the level of the grating on the results. The cubic-shaped grid used in this series of calculations was 384 × 384 × 384. The Dirichlet boundaries were applied to the bottom of the cube, where in the middle the hemisphere-capped cylinders with the aspect ratios between 1 and 10 were placed to imitate an atomic protrusion. The shape of the field distribution calculated for a protrusion with an aspect ratio of 10 is demonstrated in Fig. 4. Here the contour lines show the equipotential lines of the field distorted around a sharp protrusion.
The protrusion with an aspect ratio of 1 (a hemisphere) was included in the calculations as a reference to the value of a conductive sphere, which can be found analytically. The comparison of the solutions for a conductive sphere in an external electric field obtained analytically and by the Laplace finite-difference solver can be found in Appendix B.
In Fig. 5, the result of the dependence| F max |/| F ex | on the aspect ratio (the radius of 12 grid points for the hemispheres was the same for all the protrusions) is shown. As one can see, the data can be fit well with the linear dependence with a slope of ∼0.73, which is in good agreement with the slopes from Refs. [34] and [30], taking into account the discreteness of our Laplace solver. Here also the enhancement factor of the hemisphere (an aspect ratio of 1) is close to 3 (the result for a conductive sphere). The ratio of the level of grid grating and the calculation efficiency has been checked for a hemisphere  [30] and [34]. The convergence of the value obtained for β hemisphere toward 3, the analytical solution for a conductive sphere in the external field, with an increase of the grating level, is shown in the inset. The upper limit was defined by the computer capacity limitations. and for a protrusion with an aspect ratio of 10. In both cases, the tendency to underestimate the field enhancement was observed with the grid growing more crude (see the inset of Fig. 5). However, the finer the grid, the better the agreement obtained with the analytical estimations of the field enhancement can be found, which confirms the adequacy of the chosen technique for the solution of the Laplace equation, considering the simplicity and efficiency required for the dynamic simulations.

B. Charge induced by an electric field on surface atoms
Because there is no well-established recipe for the calculation of charge distribution on a metal surface, we intend at first to verify the proposed technique for a single, static snapshot of the surface morphology evolution.
The direct comparison of the charge induced on a surface atom or adatom is not possible, as there is no information in the literature giving the precise value from the experiment or other calculations. Moreover, the atomic charge is not a uniquely defined concept, and quantum-mechanical codes employ a wide variety of charge calculation methods [35]. One possibility of verifying the accuracy of the estimation of a charge is to use the Gauss law for the continuous flat surface Eq. (2)], and calculate precisely the charge per surface atom. However, this knowledge is very limited and fails to predict the value of a charge on a single atom for nonflat surfaces because of the field distortion on the surface features.
Nevertheless, the indirect verification of a charge induced by an electric field on an adatom is possible if the value of the critical electric field for field evaporation F c ev for a certain metal surface is known. For instance, a DFT study of the field evaporation of a tungsten surface [22] reported a value of F c ev close to the ones observed in the experiment [16]. This gives us a plausible opportunity to compare the results obtained with the present model with the results from Ref. [22], which were calculated independently by a different simulation technique, which employs a deeper description of the electronic structure.
To fulfill the aimed comparison, we constructed a tungsten (011) surface with an adatom added to the next layer in the energetically most favorable position, as shown in Fig. 6. This adatom was also considered in Ref. [22] as the most stable one. The MD code PARCAS was used to relax the surface with the adatom at 0 K during 10 ps. In these calculations we used three different empirical interatomic potentials: the Finnis-Sinclair EAM (embedded atom method) potential [36] referred to as "FS" hereafter, the Tersoff-like potential [37], hereafter referred to as "JUS," and the Dudarev-Björkas EAM potential [38], hereafter referred to as "DB." Three different interatomic potentials were used only for ensuring that a common trend exists, irrespective of the choice of the potential. Applying a drag method [39], we produced a series of static calculations by using the proposed ED-MD model to estimate the resulting forces acting on the charged adatom from the internal interatomic interactions and the external field given by Eq. (3). These forces then were integrated over the distance along the dragging path and compared to the potential energy values from Ref. [22]. The result of this comparison can be found in Fig. 7.
An analysis of Fig. 7 reveals the following. The interatomic interactions given by the FS potential are described very closely by the DFT calculations, which is seen from the comparison of the potential energy curves in Fig. 7(a) calculated for the adatom in the absence of an external electric field. When the curves pass the potential well, their behavior on the positive side of the z coordinates slightly diverge, showing that the interactions in the DFT approach are stronger and of shorter range compared to the empirical FS potential. The same conclusion is deduced from the comparison of the surface binding energies given by both the DFT and our simulation,  [36] (FS) potential and the results of DFT calculations from Ref. [22] (DFT). One curve obtained with the Dudarev-Björkas EAM potential [38] (DB) at 65 V/nm is shown for the comparison. (b) Potential energy of a charged adatom according to the DB potential and the Tersoff-like potential [37] (JUS). The same curve for the case of no external electric field from Ref. [22] is also added. which are 7.88 eV from Ref. [22] and 6.47 eV given by the FS potential. This explains the main result of our simulation, which gives the value of the critical field for W (011) surface at ∼40 V/nm, while the DFT calculation obtains no barrier if the applied field exceeds 60 V/nm. This result seems fascinatingly close to the experimental data, however, it was shown previously [37] that DFT calculations can overestimate the cohesive energy of W by ∼1.5 eV; this is why the agreement should not be overestimated. Applying another potential, for instance, the DB EAM potential [38], the critical field can be even higher, as seen in Fig. 7(b). This indicates the sensitivity of the results to the choice of the model potential. However, here we do not discuss the advantages of different models to describe the cohesive interactions between metal atoms, but rather verify the suggested technique. The reasonable agreement between our results and the results from the calculation by the DFT methods confirms the adequacy of the present approach, which was developed independently with no data calibrated against each other. The main advantage of the ED-MD model is that it allows for the simulation of dynamic change of the charging of surface atoms depending on their position regarding the neighbor positions under an electric field. Combined with the model of electronic effects described in Ref. [23], the present ED-MD model is an efficient and fairly accurate tool for the dynamic simulation of extensive metal surfaces exposed to high electric fields. The results do somewhat depend on the model potential [cf. Fig. 7(b)], thus the latter must be tested well prior to employing it within the present model.
Applying the proposed technique, we simulated the evolution of the shape of an asperity on the Cu {100} surface.
The initial shape of a tip was a simple cylinder (∼1 nm × 1 nm in height and diameter). The tip was placed on the Cu {100} surface 11.71 × 11.71 nm and ∼3 nm thick (32 × 32 × 8 unit cells). Three bottom layers were fixed to prevent a possible shift of the original cell position under the electric field. The tip was allowed to relax during 20 ps at the given temperature (T = 600 K) before the external field was applied. The Berendsen control of temperature [40] with a time constant of 20 fs was applied to prevent artificial heating of the system. The time step for the simulation was chosen as 4.06 fs for Cu and 6.9 fs for W, previously found to give good energy conservation in near-equilibrium simulations [25]. The field F ex = 8 V/nm was exerted in a ramping mode (linearly increasing) over 5 ps. A combination of values of the field and the temperature were chosen ( F ex < F b = 11 V/nm, where F b is the field at which the flat surface starts to break, and T room < T < T m , where T room and T m are the room temperature and the melting point of Cu, respectively) to achieve the efficient field-enhanced evaporation of atoms, but at a reasonable evaporation rate. The snapshots of the simulation can be found in Fig. 8(a). Here the first image (upper right-hand side) shows the snapshot when the applied field has already been ramped to the desired magnitude ( F ex = 8 V/nm) and the last one (lower right-hand side) shows the snapshot after a dramatic change of the shape of the surface tip. Two intermediate snapshots are shown to complete the imaging of the dynamics of the process. It is clear that the thermal fluctuations of the atoms elevated at 600 K enhance the atom detachment, which can be understood as field-enhanced atom evaporation. Because in the present approach we do not consider the atom ionization, the detachment is a purely thermal process, where the effect of an electric field can be considered as a permanent tensile stress applied to the surface, but following the shape of the surface morphology. If an atom has changed its position owing to the thermal fluctuations and hopped higher than the position of its neighboring atoms, it attains a higher partial charge, increasing the probability of its detachment.
The solution of the Laplace equation that is obtained at every MD step, when the surface has been changed, is shown in Fig. 8(b). The images correspond to the initial and final snapshots shown in Fig. 8(a).

IV. CONCLUSIONS
In conclusion, an atomistic model of a metal surface held under an applied electric field has been developed. The model accounts concurrently for both electrodynamic effects and interatomic interactions between the atoms. The Laplace equation with the dynamically changing Neumann boundary condition on the surface is solved at every MD step, and the following charge redistribution is allowed to depend on the position of surface atoms. The Laplace solution has been shown to give the electric-field enhancement around a hemispherically capped cylinder that is in very good agreement with previous analytical and finite-element calculations. The model allows for the dynamic simulation of surface modification under the electric-field effect, including the lattice expansion under the field and the migration of charged adatoms. However, this sum clearly diverges, because the number of atoms for a circular shell of thickness dr is proportional to its area 2πr dr. To solve this dilemma, one can consider that bare charges in a conducting material are screened by rearrangement of the conduction electrons in a metal [29], which can, in the Debye-Hückel or Fermi-Thomas approximations, be described by the screened Coulomb potential FIG. 9. Solution of the Laplace equation for a conducting sphere in an external field F ex with the finite-difference numerical method, compared to the exact analytical solution, for the field at the surface. The field was determined along the azimuthal angle θ off the z direction. The number of grid points in each direction was 200, and the radius of the sphere was 1/12 of the total size of the solution region. and hence u a = 1 a 3 /4 By requiring that u a = u c , one can fix the value of ξ to give a u a that matches u c .  (A7) Using numerical calculation of Eq. (A6) with increasing distance from the central atom until convergence, we determined that ξ = 0.6809 1/Å fulfills u c = u a to four digits of accuracy. Because both u a and u c are ∝q 2 , this result is independent of the choice of q. This value is very close to a simpler estimate obtained from Thomas-Fermi theory of 0.55 1/Å [Eq. (17.55) and Table 2.1 in Ref. [29]] and hence can be considered to be of quite reasonable magnitude.

APPENDIX B: TEST OF ACCURACY OF LAPLACE SOLVER
The finite-difference method for solving the Laplace equation is exact in the limit of the grid point size going to zero [28]. However, computer capacity limitations prevent the use of a very large number of grid points, and hence it is important to determine the relation between grid point size and accuracy of the solution. We tested the accuracy of the finite-different Laplace solver by simulating the electric field around an uncharged conducting sphere in vacuum subject to an external electric field F = F exẑ . This case is a suitable test for our solver because it has an exact analytical solution [24] that is valid everywhere in space. Moreover, the solution at the top of a tip capped with a hemisphere (considered in the main text) is approximately the same as for the sphere in the direction of the field.
Our Laplace solver was applied to solve the electric field in three dimensions around a conducting sphere region. The sphere was implemented by setting a number of grid points in the center of the solution region to be conducting boundary grid points. Similar to the use of the solver described in the main text, the potential at the surface of the conducting sphere was fixed to zero, while at the bottom and top of the cubic solution region, the gradient was set fixed to −∇φ = F exẑ (gradient boundary condition). Periodic boundary conditions were used at the side x and y boundaries. We tested the accuracy of the method systematically as a function of total number of grid points in each dimension N = N x = N y = N z and the fraction of the sphere radius r compared to the cell size f r = r/N. The results are illustrated for the maximum field at the surface of the sphere in Fig. 9 and for cross sections of the solution region along the x and z axes in Fig. 10. In general, we found that for f r 0.2 and N 100 the solution is accurate to ∼5% or better everywhere in the solution region. The multigrid and single grid (Jacobi) solution methods were found to give the same results.