Electrostatic-elastoplastic simulations of copper surface under high electric fields

Maximizing the performance of modern linear accelerators working with high gradient electromagnetic fields depends to a large extent on ability to control breakdown rates near metal surfaces in the accelerating structures. Nanoscale voids, presumably forming in the surface layers of metals during the technological processing, can be responsible for the onset of the growth of a surface protrusion. We use finite element simulations to study the evolution of annealed copper, single crystal copper and stainless steel surfaces that contain a void under high electric fields. We use a fully coupled electrostatic-elastoplastic model in the steady state. Gradually increasing the value of an external electric field, we analyze the relationship of surface failure and depth of the void for the chosen materials with different elastoplastic properties. According to our results, the stainless steel and single crystal copper surfaces demonstrate the formation of well-defined protrusions, when the external electric field reaches a certain critical value. Among the three materials, annealed copper surface starts yielding at the lowest electric fields due to the lowest Young’s modulus and yield stress. However, it produces the smallest protrusions due to a significant strain hardening characteristic for this material.


I. INTRODUCTION
The Compact Linear Collider (CLIC)-a new particle accelerator designed at CERN [1,2]-aims to perform electron-positron collisions at the energies from 0.5 to 5 TeV, with optimal performance at 3 TeV. Such high energies can be achieved in room temperature accelerating structures, where the accelerating gradient can reach the values higher than 100 MV=m [3]. One of the significant problems arising at such kind of high electric fields and limiting the maximum accelerating gradient is repetitive electrical breakdowns at metal surfaces of the accelerating structures [4]. Previously it was established that the performance of accelerating structures depends to a large extent on the rate of local breakdowns: the higher the rate, the greater the power loss during the operation of the accelerator, and the stronger the damage on the high precision geometry structures [5].
The problem of vacuum breakdown is long known and has been extensively studied in many fields. However, while the complete understanding of the breakdown phenomenon is still lacking, field emission currents are believed to be either a direct cause of the breakdown process or have a significant role in it [6]. Since field emission currents can be controllably studied using dc electric fields, the experimental studies of dc electrical breakdowns are currently being carried out at CERN. This approach also allows one to focus on separate processes indistinguishable in more complex breakdown events taking place in rf-electromagnetic fields [6,7]. These experiments have shown that the breakdowns appeared at the cathode surfaces in the dc setup at values of the surface electric fields similar to the accelerating gradients in CLIC components, which are of the order of 100-150 MV=m, leading to local macroscopic surface fields higher than 200 MV=m. These values are 50-100 times lower than the breakdown electric fields expected from theoretical predictions and experimental observations of the onset of field emissions of atoms [6,8].
The field emission currents in high electric fields follow the Fowler-Nordheim equation or its modifications [9][10][11][12]. In these equations, the experimentally observed currents can be obtained if instead of the actually applied electric field E 0 , its magnified value, E ¼ βE 0 , is used. Here β is the field enhancement factor, which in copper can be in the range β ¼ 30-140 [6].
In literature, the field enhancement factor is usually associated with surface irregular features whose geometric aspect ratios are of the same order of magnitude as the observed enhancement of the applied electric fields [13]. If the field enhances sufficiently for the field emission, these features can become relatively stable field emitting tips. Impurities or rough features remaining on the surface after the manufacturing process in the presence of an external electric field may also act as field emitting tips [14]. Indeed, improvement of surface preparation before applying a high electric field leads to a considerable decrease of the breakdown rate; however, achieving values of the breakdown field close to the theoretical limit still remains beyond reach.
Although field emitting tips are usually assumed to preexist, the features able to enhance the external field locally up to hundreds of times have not been observed experimentally [6]. On the other hand, it was noticed that there exists a clear correlation between the value of average minimum electric fields which caused breakdowns at different metal surfaces with the lattice structures of these metals [15], linking the onset mechanisms of electrical breakdowns to the dislocation activity and plastic deformations of accelerating structure material. This hypothesis was further supported in [16] where an analytic model, capable of fitting breakdown rate and applied electric field, was presented.
The experimental study of the electrical breakdowns is complicated by the range of values of the related parameters differing by many orders of magnitude. The estimation from the Fowler-Nordheim equation [17] gives that a single field emission area covers typically only 10 -20 − 10 −16 m 2 of the surface, single emitter, capable of melting during rf experiments has diameter in range 17-25 nm [18], while a typical single damage area seen on the surface of the cathode is of the order of 10 −10 − 10 −8 m 2 [19]. The background currents 10-12 A rise up to 100 A [20] and the current densities at the field emitting tips rise quickly to 10 11 -10 12 A=m 2 [21,22] destroying the tip in a submicrosecond time interval [23,24]. Since this kind of spontaneously created nanoscale sharp stable features (tips) on metal surfaces with dimensions matching the observed significant enhancement of the external electric field were not observed experimentally, it can be suggested that a metal surface evolves under high electric fields in such a way that the features form dynamically, giving rise to field emission currents. Nevertheless, the exact mechanism which can explain the formation of these field emitting tips is not yet known.
One of the hypotheses, which can explain the dynamic formation of a tip on metal surfaces, is a surface protrusion promptly appearing in the presence of an external electric field. Previously several studies were addressing the possible scenarios of surface evolution in the presence of high gradient accelerating rf fields. In studies conducted by Norem et al., different mechanisms possibly leading to a breakdown event are simulated by using a molecular dynamics method [25,26]. For instance, it was demonstrated [25] that large clusters of atoms can evaporate from preexisting protrusions under these conditions. In other studies [27] the influence of Joule heating and Nottingham effect on heating of a microprotrusion was analyzed. In [28], the role of micrometer size molten particles was investigated in high gradient systems. However, in these works, preexisting protrusion or surface defect was assumed. On the other hand, in Refs. [29,30] by the finite element method it was shown that the field enhancement can also appear at microcracks, which are usually observed on surfaces due to fatigue after many operating cycles. The protrusion formation hypothesis was evaluated in [31], where the dislocation nucleation from a near surface void under static tensile stress due to an external dc electric field was simulated. The results showed that mass transport above the void may initiate the formation of a surface protrusion, possibly leading to the formation of the field emitting tips. Moreover, in [32] it was shown that stability of the field emitting tips depends on their dimensions and the Joule heating, due to strong field emission currents can significantly modify the shape of needlelike nanoscale protrusions. The dynamic formation of protrusions can be also motivated by the studies made in [33], where the authors saw the flattening of Cu nanoclusters placed on a copper surface to a single monolayer due to the effect of minimization of surface energy. In [31,32,34] a preexisting stress concentrator in the shape of a void was simulated (this was selected for the concreteness and simplicity of the model and serves as a source of a shear stress component). However, sufficient explanations for existence/appearance of field emitters is not yet provided. Understanding of mechanisms leading to the modifications of metal surfaces under high external electric fields, which may lead to vacuum breakdowns, is one of the essential factors in the successful design of future accelerators such as CLIC [1,2]. If these mechanisms are understood, new materials or more effective conditioning methods for accelerating structures can be developed. However, the complexity of the problem -the nanoscale phenomena appearing at random during experimental time scales up to hours and days-limits the available simulation methods. The atomistic simulation methods [31,32,34], such as molecular dynamics (MD) techniques, can give an atomic level resolution, but are limited to picosecond to nanosecond time scales. Moreover, for reaching reasonable computational times, the applied forces and electric field in the MD simulations need to be significantly exaggerated [31,32,34], leading to a difficult comparison between the experimental and simulated results. To overcome this problem, we use the finite element method as an intermediate tool between the experiment and MD simulations. Even at the nanoscale limit of its applicability, this method can give a general understanding of the process, e.g., the effect of different mechanical properties on the onset of plastic deformations. However, unlike molecular dynamics, the finite element method (FEM) does not consider individual atoms, but rather averages the material properties and calculates the quantities over finite size elements. This allows us to study longer time scales and larger systems than possible by MD, but introduces limits to the smallest available length and time scales. The range of the applicability of the method can be estimated from previously conducted works. In [35] the surface effects and effective mechanical properties of nanocomposites were estimated. It was found that the nanoscale size effects had significant influence to the material properties if the simulated nanovoids had r < 10 nm. In [36] the limits of continuum approach for describing the inelastic deformation of material (creep) were studied, requesting the size of the nanostructures to have radius >2 nm. In comparison with molecular dynamics methods [31,32,34], FEM has the advantage of being able to easily model macroscopic systems and long time scales to predicting the behavior of experimental macroscopic structures. It enables scanning of many different configurations to pinpoint the most intriguing ones for future studies.
In the current paper we aim to investigate the surface behavior of the metals, which are considered for building the CLIC accelerating structures [7]. The aim of the paper is fulfilled by studying plastic deformations due to applied high electric fields in three different model metals by using the finite element method (FEM). We investigate the effect of void location, by varying its depth with respect to the material surface, on the process of plastic deformation that lead to the field enhancing surface modifications in these materials.

A. Materials of interest
We investigated the mechanism behind the growth of field enhancing tips in type 304 stainless steel [37], single crystal copper [38] and high conductivity oxygen-free copper (Cu-OFE REF. UNS C10100 Grade 1, 99.99% purity), annealed 2 h at 1000°C (annealed Cu) and used in the accelerating structures of CLIC [39]. However, all these materials were modeled by using the corresponding elastic parameters averaged over experimental sample volumes, i.e., grains are not modeled explicitly. The choice of model materials was motivated by the clear difference in their mechanical properties. The stiff mechanical properties of the steel and the soft mechanical properties of the annealed copper are compared to the single crystal copper, which was simulated previously by MD [31]. The elastoplastic properties of all the materials are presented in Table I. In the calculations, we use isotropic materials. The parameters of simulated annealed copper are phenomenological and are fitted to follow the experimental data to characterize the experimental stress-strain behavior [ Fig. 3(a), dashed line]. In single crystal copper simulations, the (001) crystal face is exposed to the external electric field.
The geometry of the simulated system consists of a sample containing a single void and vacuum, where the electric field is applied. The shape of the void is chosen spherical since this shape minimizes the surface energy. The cross section of the simulated geometry is presented in Fig. 1. The whole simulated cell, where the calculations are actually done, is much bigger than the region used for the visualization postprocessing, marked with the dashed outline. The geometrical parameters are specified in Table II. Since the void is of a spherical shape, the number of parameters of interest for the current study can be reduced to two, the radius and the depth of the void. Moreover, the spherical shape allows for the computational efficiency by utilizing the axial symmetry of the simulated box relative to the center of the box, as the tensorial formulation of the structural mechanics problem is coordinate system invariant. As a result, numerical implementation of the calculations is optimized so that the computations are conducted only in one symmetrical (2D) slice of the problem without loss of generality.

Modeling of structural mechanics
We simulate the deformation of the materials by using COMSOL MULTIPHYSICS 4.3B and NONLINEAR STRUCTURAL MECHANICS TOOLBOX [42]. The model includes fully coupled elastoplastic-electrostatic interactions where the geometrical nonlinearity, arising during the deformation of the sample, is handled by utilizing the deformed meshes and large strains. Such a model enables for direct investigation of effect of electric fields which is self-consistently recalculated following the surface topography. The mechanical stress and deformations are calculated by using the elastoplastic model of isotropic materials with nonlinear material hardening and large strains, using multiplicative decomposition of the deformation gradient [43][44][45]. The details of the approach are specified in Appendixes.
The mechanical stress applied to the surface of each sample due to the electric field is calculated by using the Maxwell stress tensor for the electrostatic field [31], where p is pressure, ε 0 the vacuum permittivity and E n the electric field at the surface, which is normal to the metal surface in the dc case.
In this manner the applied electric field can cause deformation in the simulated samples, which can be calculated by FEM. We modeled this deformation by using arbitrary Lagrangian-Eulerian formulation (ALE) [46]. The ALE method allows for flexible deformation of material boundaries without remeshing the geometry. In the material region, the mesh is deformed naturally by using the calculated material deformation from the elastoplastic deformation calculations. In the vacuum region, however, the precalculated deformation is not available. Here, the geometry is altered by deforming the vacuum region boundary corresponding to the sample deformation. The subsequent mesh relaxation is performed by the Laplace smoothing with the following constraints: (i) Deformation of the sample boundary is calculated by the Nonlinear solid mechanics toolbox. (ii) Normal directional deformation at the sides of the vacuum region is not allowed. (iii) Top side of the vacuum region is fixed.

Validation of material model and material parameters
Experimental tensile strength tests are reproduced by simulations to ensure the expected behavior of material model at large plastic strains. For validation we used annealed copper parameters as it is currently the most favored material for CLIC components. Also, annealed copper demonstrates highly nonlinear behavior during tensile tests, which is fairly challenging for computational test simulations. In tensile strength tests only mechanical stress (no electric field) is applied to the samples, thus, we calculate only elastoplastic deformation under the constant external stress. The simulated test sample is a single metal bar, with 50 mm height and 10 mm diameter. The sample is deformed by incrementally increasing displacements in a parametric simulation by applying displacement boundary condition at the bottom and the top of the sample, until the ultimate tensile strength is reached. The maximum von Mises stress is interpreted to be equivalent to the experimental stress values, since only normal stress is present in the calculations.
In a uniaxial tensile test, only one stress component is present (σ 11 ) and the material fails, if σ xx ¼ σ crititcal . The von Mises stress (Eq. (B4) [47]) can be considered as a uniaxial equivalent to the multiaxial stress; it combines the stress tensor components into a single quantity and is often used for material failure criteria [47]. The plastic deformation starts, when the components of the stress tensor or their combination exceed a critical value. If the von Mises stress becomes larger than yield strength of the material, the material deforms plastically (details of the material yielding are specified in Appendix).
During the validation, cylindrical symmetry of the geometry in conjunction with the isotropic material assumption was utilized, making it possible to conduct calculations characterizing the whole 3D sample by 2D simulations. In the calculations, ∼2100 quadratic quad mesh elements were used.
To evaluate the validity of the applied model on the nanoscale as well, we performed also the nanoindentation simulation and compared the result to the available experimental nanoindentation data [38]. For these simulations we used the single crystal stress-strain curve from Ref. [38] in order to estimate the initial yield stress and hardening function of the single crystal copper. To ensure the comparability of our results and the experimental data we reproduced the main details of the nanoindentation experiment [38] in our simulation: the pressure to the sample of 15 μm in diameter and 10 μm in height was applied by a 3.4 μm radius needle. The radius of the sample, thus, was sufficiently large compared to the needle and the contact area. In our calculations, the needle was considered rigid, penetrating into the surface parametrically from 0 to 150 nm. The applied force was calculated from the material response by integrating the pressure over the top surface of the sample. All calculations in the parametric depth sweep were performed in the steady state and the interaction between the sample and the needle was simulated by using a frictionless contact pair. Since the problem was highly nonlinear, we utilized adaptive meshing [48] to obtain a mesh independent solution. In the final mesh, ∼7000 triangular linear elements were used.

Effect of electric field
The electric field is calculated by solving the Laplace equation: At the material surface and at the top of the simulation box we keep the value of the potential fixed, but allow for ramping of the voltage incrementally from zero until final deformation. The size of the vacuum (h V ) is chosen so that even the highest protrusion seen in the simulations, compared to the vacuum region height, is very small, as even at the largest deformations, h protrusion =h V ¼ 0.0014. This guarantees that the field is constant at the top of the vacuum (Fig. 1, point 1) within an insignificantly small numerical error. The applied (not enhanced) electric field during the simulation is measured at point 1 as shown in Fig. 1. Using this type of boundary condition for calculation of the electric field in the chosen geometry, we provide additional numerical robustness and increase the computational efficiency of the simulations while emulating the experimental dc setup, where the metal anode is used. Finally, at the sides of the vacuum box (in the lateral direction) we applied the n · ∇φ ¼ 0 boundary condition (where n is unit normal vector at the boundary), keeping the normal electric field zero, allowing us to emulate the continuation of the material infinitely far in the lateral directions.

Electromechanical coupling
Since in our study we focus on the effect of electric field, we follow in the simulations the ramping of the field, which is done parametrically from 0 to 5000 MV=m for the annealed copper and up to 10 000 MV=m for all other materials. The annealed copper demonstrated the most ductile and stainless steel the stiffest properties from the simulated materials, requiring corresponding scaling of the maximum applied electric field. All materials are simulated for different void depths, while controlling the aspect ratio of the void. Due to the scale invariance, we do not modify the void radius, but the aspect ratio, defined as a ratio of depth to void radius (h=r), in the range from 0.2 to 2. The scale invariance has been confirmed by obtaining the same field enhancement factor, electric field and stress distribution while scaling the system with the preserved aspect ratio for all the tested voids.
The parametric electric field ramping was implemented by increasing the applied electric field value in each parametric step at most by 5 MV=m. For each parametric, electric field value, steady state equations for electric field, elastoplastic deformation and mesh deformation/ smoothing were solved. For numerical solution of the equations, we used the damped Newton method as the nonlinear solver [49] in conjunction with MUMPS [50,51].
The calculation of the large deformations caused by the high electric field requires a series of stopping conditions to be specified. We consider the simulations to be ended, if any of the following conditions are met: (i) the maximum electric field is reached; (ii) mesh quality becomes lower than 0.3 (1 for perfectly regular mesh element and 0 for degenerated mesh element); (iii) explosive protrusion growth-when the simulated electric field increment became lower than 50 V=m. This corresponds to the mechanical failure of the material.
In the calculations, we use the quadratic triangular finite elements. The mesh is divided in two regions with the different element densities. The visualization area (Fig. 1) has ∼12000 elements; in the rest of the geometry there are ∼3000 elements. The same type of mesh is used for all materials for the corresponding cases of identical aspect ratios. The mesh convergence was tested for all simulations.

A. Validation of material model
In Fig. 2 we show the comparison of the tensile strength test results for annealed copper, obtained experimentally in [39] and calculated in the present simulations as described in Sec. II A1-A2. The simulated stress-strain curve shows linear behavior until the initial yield stress is reached. Then, strain hardening takes place exhibiting the nonlinear behavior of the stress-strain curve, which was also taken into account in the simulations (Appendix C). The experimental and simulated results agree very well until the ultimate tensile strength is reached. At this point, necking of the sample starts, restraining the continuation of calculations. The model manages to capture effectively the necking effect, as it is seen in Fig. 2(a); however, as the neck develops, it concentrates the stress and further material deformation takes place in this area, decreasing progressively its diameter during the test and changing significantly the geometry of the sample. The simulations of the tensile stress after the ultimate tensile strength is reached, requires experimental data of the corresponding neck radii. As currently the latter is not available, the accuracy of the simulation results is given by the ultimate tensile strength values.
The force required for a nanoindenter (a probe) to be moved to a certain depth as the function of this depth is given in Fig. 2(b) for the experimental [38] and the simulated results. The experimental nanoindentation data [38] is mostly linear, with a small nonlinear part at the beginning of the test. The simulated force-depth curve follows the experimental results closely, slightly overestimating the force required for the deformation up to the probe depth ∼50 nm. Then, the value of the simulated force decreases until it becomes equal to the measured value at the probe depth ∼85 nm. The decrease continues until the force is underestimated by ∼20%, compared to the experiment, at the probe depth 150 nm. This kind of behavior can be explained by the insufficient experimental data: single crystal copper strain hardening data [38] is available only up to strains 0.1. When the effective plastic strain ε pe ¼ 0.1 in the simulated nanoindentation system is reached at the probe depth ∼50 nm, the further tests must rely on the extrapolation of strain hardening data. The maximum strain 0.15, at which the ultimate tensile strength is reached in annealed copper, corresponds in the nanoindentation test to the probe depth of 70-80 nm. In this range, the simulated nanoindentation results coincide almost completely with the experimental data, providing accurate description of the material deformation and making a continuum model suitable also for nanoscale and microscale simulations, if suitable material data is available. However, the model has several limitations, like the assumption of elastically isotropic material and practical requirements for medium strains (at very large strains the deformation of the mesh becomes unreasonable). Clearly the model fails to capture the crystal structure anisotropy at larger strains. For the current purpose of a study of plastic deformation in metals under high electric fields, however, these limitations are acceptable.

B. Evolution of a protrusion in the simulated materials 1. Effect of material properties
Snapshots of the simulation results of all materials at the end of the simulations are presented in the Fig. 3, where we show the electric field and corresponding von Mises stress distributions. The isocontours show the effective plastic strain values at 0, 0.01 and 0.05 levels. The isocontour of the zero effective plastic strain separates plastically and elastically deformed regions, while isocontours with higher effective plastic strain values border the areas where most of the irreversible deformation is taking place. All the presented materials, stainless steel, single crystal copper and annealed copper demonstrate a fairly similar shape of protrusions. The small shape variations are explained by the slightly different end points of each calculation which depend on the stopping conditions. However, the major difference appears in the sizes and shapes of the plastically deformed regions. In the stainless steel, plastic deformation is distributed only near the void, illustrating clearly the stress concentration properties of the void. In case of the annealed copper, the material is almost fully plastically deformed, as the stress caused by the applied electric field exceeds the initial yield stress of the sample. Only a small localized region beneath the void is deformed elastically.
At higher effective plastic strain values, at ε pe ¼ 0.01, the spatial distribution of plastic deformation in the single crystal copper and stainless steel is already fairly similar, while the annealed copper demonstrates the widest distribution of the effective plastic strain exceeding that of the other materials. The most significant material deformation area is enclosed by the 0.05 effective plastic strain isocontour, where the majority of the irreversible  [38] for the nanoidentation test in single crystal copper. In the simulation the experimental geometry is fully reproduced to ensure the comparability of the results (see text for details).
deformation is taking place. The 0.05 isocontour is located in the area close to the void surface, characterizing its stress concentrating properties and its interaction with the nearby surface of the sample, as the maximum plastic strain develops between the side of the void and the close-by material surface. Since the effective plastic strain in the rest of the material, except in close proximity of the void, is small, the influence of the bulk material to the material failure is limited and dominated by the localized interaction of the void and surface.
In the distributions of the von Mises stress we observe only small differences, although the stress values are considerably different. In stainless steel, the maximum stress is concentrated near the top of the void, while in annealed copper the maximum stress is concentrated at the void sides.
In the following we will mostly focus on the results for the stainless steel and the annealed copper since the mechanical properties of the single-crystal copper lie between these two extreme cases.
The simulation results of the stainless steel with h=r ¼ 0.5 are presented in Fig. 4(a). Here, the series of snapshots demonstrates the evolution of the von Mises stress and the electric field distributions during the ramping of the applied electric field. The first snapshot shows the surface response at the external electric field 5000 MV=m; at this moment the material starts to yield and the plastic deformation appears, as the dislocations nucleate and propagate. A further increase of the electric field widens the plastically deformed volume and at the field strength 5000-6240 MV=m, the plastic deformation spreads all over the upper part of the void. Within the deformed volume, as the electric field is ramped and the stress in the material is increased, the interactions of the dislocations lead to the strain hardening. As a result, the yield stress increases, which prevents appearance of a protrusion. The protrusion on the surface becomes visible later on, when a sufficiently large volume of the metal around the void is plastically deformed. This happens at the field ∼6650 MV=m. We observe somewhat similar behavior in the annealed copper. The results presented in Fig. 4(b) show that a protrusion appears only after the plastic deformation fills the volume around the void entirely. Similarly to the stainless steel, we observe the first plastic deformation at the sides of the void since the shear stress in this area is maximal, but at the notably lower field-already at 500-600 MV=m, due to the considerably lower yield stress. The plastically deformed volume widens fast during the ramping of the electric field and at 1500 MV=m, the void is completely surrounded by the plastically deformed material. We see a small protrusion growth taking place already at the field strength 3500 MV=m. The significant deformation of the material above the void, though, happens at 4500 MV=m when the strength of electric field becomes sufficient for the large plastic deformations around the void.
This result can be explained by the mechanical properties of the stainless steel and annealed copper (see Table I). The Young's modulus of the copper is almost 2 times smaller than that of the stainless steel. The yield stresses of both materials also differ more than 50 times, making the annealed copper to yield at the lower stress. However, extremely soft characteristics of the annealed copper are compensated by significant strain hardening, decreasing the ductility of the material.

Effect to the void depth
In Fig. 5, the formation of a protrusion above the voids placed at the different depths under the surface in the annealed copper is shown. The colors represent von Mises stress and electric field strength. The isocontour represents the zero effective plastic strain, outlining the plastically deformed volume of material. Keeping the radius of the void the same, we change the aspect ratio of the void by varying its depth: h=r ¼ 1-0.2, with the increment 0.2.
We see that in all the cases the protrusion growth takes place only in the plastically deformed regions of the material. It is interesting to note that the location of the material failure depends on the h=r ratio. When the h=r is small (h=r ¼ 0.2), the deformation takes place mainly between the top of the void and the surface. The increase of the h=r ratio moves the location of the maximum von Mises stress towards the sides of the void. Thus, we conclude that the failure of material is caused by different mechanisms depending on the depth of the void.
Comparing the first and the last snapshots of protrusions in Fig. 5, h=r ¼ 0.2 and h=r ¼ 1, respectively, we see that the protrusion caused by the shallow void is well defined and leads to the significant material deformation, while the protrusion caused by the deepest void is small. The field around the top of these protrusions is enhanced 1.7 and 1.3 times, respectively. When the void is located even deeper in the material, the field enhancement is even smaller. Further, analyzing the distribution of the von Mises stress near the void we see that the deeper void results in lesser values of the von Mises stress at the surface of the material. The distribution of the von Mises stress near the even deeper void (h=r > 1) decays to zero at the surface, preventing the dislocation-mediated formation of a protrusion on the surface. Hence, another possible mechanism must explain the protrusion growth from a void positioned deep beneath the surface.

Field enhancement effect due to protrusion
We also analyze the field enhancement factor, a universal quantity depending only on geometric parameters of surface protrusions and not on the mechanical properties of the surface itself. Moreover, the field enhancement characterizes field emission current, which can lead to development of an electrical breakdown spot.
The evolution of this quantity, relating closely to the protrusion size, does depend on the surface mechanical properties. The field enhancement factor for the stainless steel as a function of the external electric field with the different h=r ratios of voids is shown in Fig. 6. The field enhancement is almost constant and close to 1 (no enhancement), until a certain critical value is achieved. Then, the field enhancement factor starts increasing in the catastrophic manner, following the sudden growth of a protrusion as seen in Fig. 6: the sudden protrusion growth on the stainless steel surface occurs after plastic deformation begins. The deeper the void, the stronger the external field needs to be applied to activate the process of formation of a protrusion on the surface. The catastrophic growth of the field enhancement is observed up to h=r ¼ 1; after that the growth is fast, but less dramatic.
We observe somewhat similar behavior for the annealed copper [ Fig. 6(b)]. Unlike the case of the stainless steel, the field enhancement factor for the annealed copper grows continuously and smoothly during the ramping of the electric field for all the simulated h=r ratios. This corresponds to the slow but continuous deformation of the material surface near the void. Since the hardening behavior in the tensile test of both, the annealed copper and the stainless steel, is fairly similar (Fig. 2) and the initial yield stress affects only the electric field value which starts the plastic deformation, the continuous deformation of the annealed copper can be explained solely by its low Young's modulus.
The simulation results show that the stainless steel, while leading to a protrusion growth at larger field strength than other simulated materials, undergoes explosive and very fast surface deformation and corresponding field enhancement increase. Thus, if sufficiently high electric field is applied to the stainless steel, even a small variation in the field strength might lead to the development of an electrical breakdown event.
In the current simulations, although using FEM, the macroscopic electric field is still up scaled, as we consider the idealized material surface (no naturally preexisting surface roughening) and a single bulk material defect (a spherical void). The absence of surface defects and the mechanisms activated by electric currents and temperature rise can explain the difference between the simulated value of the field, when we observed the first deformation in metals, and the average lowest electric field before a breakdown event takes place, measured in the experiments: 1000 MV=m in our simulations versus ∼200 MV=m in dc experiments [18]. Similarly, we observe a field enhancement β lower than the values typically reported in the experimental studies (β ¼ 30-140) [6]. In the case of real materials, some residual surface roughness would always remain. Moreover, if an electric breakdown occurs, it will lead to melting of breakdown area and introduce additional roughness in range of tens of micrometers [52], which act as field enhancing features and lead to local enhancement of applied electric field. The field enhancement factors of stacking surface defects multiply due to Schottky conjecture and increase local surface electric field near the actual surface significantly, to a higher level than applied field [53]. For example, even a semispherical surface defect will cause the field enhancement ∼3 times [54]. However, our simulations clearly demonstrate that a material surface is modified under the tensile stress due to the field, when subsurface material defects are present, resulting in a clear local field enhancement. In the simulations we use the elastic parameters of materials close to the experimental ones. Nevertheless, the simulated systems remain idealized; already nonspherical voids may cause higher local stress concentrations inside the material, which would enhance the surface feature growth, lead to more complex distribution of plastic deformation and hence lower the field at which breakdown occurs. However, these results suggest that combination of several, possibly nonspherical voids or other stress concentrating defects would finally be needed to trigger the breakdown event. Local surface roughness above the void or formation of sharp edges around the protrusion due to severe plastic deformations will also lower the value of the breakdown field. Moreover, in a more complex realistic picture of a vacuum breakdown event, there are many other mechanisms (such as electromigration and dislocation reactions) which can interplay with the near-surface void, enhancing its effect. Unfortunately these mechanisms, which are responsible for the modification of presurface regions under high electric fields, are not accessible by FEM methods, but will be subject to future studies by different techniques. Currently we analyze the effect of plastic response of materials with different elastic properties on the tensile stress due to the field.

C. Material failure mechanisms
To investigate the material failure mechanisms in the simulated systems, we investigate the maximum von Mises stress at the material and void surfaces. These results, as a function of the external electric field for several aspect ratios of voids, are presented in Fig. 7. The plotted data represent voids, having most important material failure characteristics. Figure 7(a) shows the results for the stainless steel and Fig. 7(b) for the annealed copper samples. The von Mises stresses at the surfaces of the material (solid lines) and of the void (dashed lines) are shown separately; the applied electric field for all presented aspect ratios starts from zero while each grid segment in the figure corresponds to 2000 MV=m.
The curves in Fig. 7(a) have three clearly distinguished regions. In the first region, the material deforms elastically until the initial yield stress is reached. In the following (the second region) the material yields keeping the von Mises stress at the constant value until a critical field strength is reached. In this region of the von Mises stress plateau the material around the void hardens and the plastically deformed volume around the void grows. Finally, if the sufficient volume around the void is plastically deformed, the third region in the curve is reached and the protrusion at the surface starts to grow fast.
Comparing the maximum von Mises stress at the material and the void surfaces we observe the general trend. The material failure can start at either surface, depending on the aspect ratio of the void, if the stress is sufficiently large. For the small h=r ratios, the stress at the sample surface develops faster than the stress at the void surface, as demonstrated for h=r ¼ 0.2. It is interesting to note that the aspect ratio h=r ¼ 0.3 equalizes the probability of material failure at both surfaces since the stress there develops simultaneously (both curves coincide fully at h=r ¼ 0.3). For deeper voids the yield stress is reached first at the void surface, and then, it yields at the material surface. Thus we conclude that the mechanism of material failure in the surface, which contains a void, depends on the depth of the void. For the shallow voids, the depth of the void is about one-fifth of its radius the material fails at its surface, while after the void depth becomes deeper that a third of its radius, the failure appears at the void surface at first. The deeper the void the bigger the difference between the electric field values at which the material fails first at the void surface and then, at the surface of the material, see Fig. 7.
We observe a very different behavior of the maximum von Mises stress for the annealed copper. Here, at first the material behaves elastically until the yield stress is reached. Since the yield strength of the annealed copper is very low, the elastic limit is reached very fast. Instead of a clearly visible plateau seen at this point, as in the stainless steel stress-electric field curves, there is only a very small region, where the behavior of the material changes. To explain this kind of elastoplastic response of the annealed copper on the stress generated by high electric fields, we consider three factors defining the response of the material in the present material-Young's modulus, yield stress, and strain hardening behavior. The yield stress determines the height of the plateau, the stress at which the plastic deformation starts. The stress-strain behavior of the materials is approximately the same for both materials in the region where the linear approximation is valid [ Fig. 2(a)], due to combination of elastic and plastic properties of the materials. However, a significant difference lies in the hardening behavior of the materials, which characterizes combined interactions of dislocations in the sample. When the stainless steel shows linear hardening during the whole range of applied fields (and resulting stresses), the hardening of the annealed copper is nonlinear, decreases slowly and approaches to perfectly plastic behavior (and possibly softening) at the end of the simulations ( Fig. 9 in Appendix), resulting in a more monotonic increase of the von Mises stress compared to the stiff stainless steel surfaces. The bottleneck of the material failure shifts from the surface of the material to the surface of the void in the similar manner as in the stainless steel. However, the soft properties of the copper lead to the surface of the material being susceptible to the earlier failure compared to the void surface of the slightly deeper voids than in the stainless steel. For the void positioned at 1=3 of its radius below the surface, as can be seen in Fig. 8(b), the gradual increase of the electric field causes the material surface failure before the yield of the void surface occurs.
For completeness of the present study, we also performed the same analysis of the maximum of the von Mises stress in the single crystal copper surface. The behavior of the stress-electric field curve in this surface was very similar to the one in the stainless steel. The initial, elastic response was followed by the von Mises stress plateau at the value of the yield strength followed by the last stage of the fast increase of the stress. Since the single crystal copper was softer than stainless steel (but stronger than annealed copper), the electrical field required to deform the sample was considerably lower. The void with the aspect ratio h=r ¼ 0.5 caused initial yielding of the material at 2450 MV=m. Fast protrusion growth, similar to the stainless steel in Fig. 7(a) started at ∼3500 MV=m.
The anomaly at the end of the von Mises stress-electric field curves in Fig. 7(b) for the annealed copper corresponds to the ultimate tensile strength level. From this value on, no further strain hardening takes place and material is approximated to behave as perfectly plastic due to the lack of experimental data for this regime of the von Mises stress. However, this kind of approximation has a very small influence on the overall simulation results, since the limits of the experimental hardening data are hit only at the very end of the simulations, when significant deformation has already occurred. Moreover, when we compare the experimental and the simulated stress-strain responses in Fig. 2, we also see that the experimental annealed copper actually undergoes some softening after the ultimate tensile stress is reached. Since the simulated material response in this regime is perfectly plastic it requires higher stress for the deformation to continue. Thus, in our simulations of the annealed copper we slightly overestimate the required electric field for the final breakdown of the material.
We also analyze the external electric field at the material yielding as a function of the aspect ratio of the voids in all three materials. The results are present in Fig. 8. Here again, the solid and dashed lines show the electric fields at which the stresses at the material and the void surfaces, respectively, reach the yield stress. In all three cases the curves saturate at some values of h=r ratio. This indicates that at some depth (at h=r ∼ 1) the interaction between the void surface and the material surface through the strain field caused by the external electric field disappears, although the yielding of the material still occurs. Note that in all cases in the first region, the electric field needed to yield the material grows almost linearly until the critical depth is reached. The yielding at the void surface does not show a strong dependence on the depth of the void and happens at the electric field strengths much lower than those required to yield the material surface. This tendency is seen for all the voids at the depth greater than a 1=3 of their radius. Figure 8 allows for identification of three separate mechanisms of material failure due to the presence of near surface void. For h=r < 0.3 the material yields at the surface above the top of the void. If h=r ¼ 0.3, the yielding starts equally at the both surfaces-the material and the sides of the void. If h=r > 0.3, the stress is concentrated on the sides of the void, where the material yields first. When the critical depth of the void is reached around h=r ∼ 1, the location of the void becomes too deep in the material and the deformation of the sample in this case is caused by the changes in the effective bulk properties of the material-a too deep void does not interact with sample surface, but it acts as a stress concentrator by decreasing the cross section area of the material.

IV. CONCLUSIONS
In the current study, we investigate three different materials, considered for building the CLIC accelerating structures. We studied the metal surface behavior under high external electric field to investigate the mechanisms leading to the changes in surface morphology, which can be viewed experimentally as localized spots of high field emission currents-precursors of breakdown events. We show that the presence of a near surface void in metals indeed can provoke a local field enhancing protrusion growth if the metal surface is exposed to high electric fields. We restrict ourselves to three representative cases, the stainless steel, the single crystal copper and the annealed copper, to study the different regimes of the material failure in such an extreme FIG. 8. External electric field at which the stress at the material (solid lines) and the void (dashed lines) surfaces reach the yield stress as a function of the void aspect ratio for the stainless steel, the single copper and the annealed copper surfaces containing voids. condition by simulating elastoplastic deformation of the material using continuum mechanics methods. Although, applying FEM in the simulations reduces seriously exaggerated electric fields necessary in MD studies, the dynamic processes are not accessible by FEM. Moreover, the smallest size accessible by FEM is 10 nm, which is the border line of applicability of this method.
We see that all simulated materials behave similarly with respect to the depth of void location: very shallow voids cause first the failure of material surface, and at the depths greater than a third of their radius the failure always starts at the void surface. All three materials show the failure at the void surface at lower values of electric fields than the corresponding failure at the material surface. While the simulations demonstrate material surface deformation due to the applied electric field, field enhancement factor in the range 50-100 is not obtained leading to very high surface field values at the end of the simulation. Higher field enhancement factors could be achieved, if the simulations would continue until the material fracture and opening of the void or if simulations would incorporate preexisting surface imperfections. However, the simulations clearly demonstrate the possibility of field enhancing material surface modification due to the subsurface defects.
The main difference in the elastoplastic response for all materials is mostly explained by plastic behavior of the materials. In the annealed copper the von Mises stress grows monotonically even after the yield stress is reached showing the strong intermixing of elastic and plastic deformations. Although the strength of the electric field is different at the yielding of the investigated materials, we clearly see that only shallow voids (with aspect ratio h=r < 1) allow for the volume between the void and material surfaces to be sufficiently deformed for nucleation of a surface protrusion. However, the failure at the void surfaces does not depend significantly on the void depth and takes place at electric fields significantly lower than those needed to yield the material surfaces.

ACKNOWLEDGMENTS
This work was supported by Estonian Research Council Grant No. PUT 57, ERDF programme "Research internationalisation", Mecbeth project granted by Academy of Finland and Magnus Ehnrooth foundation for travel support. The authors would like to thank Walter Wuensch and Markus Aicheler from CERN for stimulating discussions and insights.

APPENDIX A: SOLID MECHANICS
During the simulations, large deformations are assumed. Thus, to simulate accurately the deformation of the sample, Green-Lagrange strains (ε) and second Piola-Kirchoff stresses (S) are used. Since the system has no initial stress or strain, the material model is expressed as [45] S ¼ C∶ε el : Since large plastic strains are assumed, the deformation gradient tensor F is multiplicatively decomposed to elastic and plastic parts [43,55,56]: The plastic deformation is removed from the total deformation gradient as F el ¼ FF −1 p and the elastic Green-Lagrange strain tensor is Similarly, the plastic Green-Lagrange strain tensor is calculated as The elastic and plastic Green-Lagrange strain tensors are related as

APPENDIX B: YIELD FUNCTION AND ISOTROPIC PLASTICITY
If associate flow rule is applied and the plastic potential coincides with yield surface, the plastic flow for the large strains, is written as [43,55] with the Kuhn-Tucker condition for the plastic multiplier λ and the yield function F y : λ ≥ 0; F y ≤ 0 and λF y ¼ 0: In current formulation, the elastic left Cauchy-Green tensor is written in terms of the deformation gradient and right Cauchy-Green tensor is B el ¼ FC −1 p F while the plastic flow rule is solved for the inverse of the plastic deformation gradient F −1 p . Finally, the plastic Cauchy-Green rate and elastic left Cauchy-Green tensor are expressed as In the case of associate flow rule, the yield function must be continuously differentiable with respect to the stress: When the von Mises yield condition is used, the effective stress ϕðσÞ is defined as φðσÞ ¼ σ mises ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 ½ðσ xx − σ yy Þ 2 þ ðσ yy − σ zz Þ 2 þ ðσ zz − σ xx Þ 2 þ 3ðτ 2 xy þ τ 2 yz þ τ 2 zx Þ r ; ðB4Þ where σ ij are the normal stress components and τ ij the shear stress components (i; j ¼ x; y; z) When the material hardening takes place, σ ys ¼ σ ys ðε pe Þ, where ε pe is effective plastic strain and σ ys ðε pe Þ ¼ σ ys;0 þ σ h ðε pe Þ: ðB5Þ

APPENDIX C: MATERIAL HARDENING MODEL
In the present study we use isotropic nonlinear material hardening to capture highly nonlinear stressstrain behavior in the simulated materials. If the material undergoes plastid deformation, it becomes gradually more difficult to deform, until ultimate tensile strength of the sample is obtained. To capture this behavior, the strain hardening of the material, the initial yield stress is allowed to increase as a function of plastic deformation according to Eq. (B5). The plastic deformation dependent part of the yield strength is presented in Fig. 9, where the red curve corresponds to the annealed copper, the green curve to the single crystal copper and the magenta curve to the stainless steel. The material hardening data is extracted from the experimental stressstrain data by using a MATLAB script, assuming the offset yield point at 0.2% of the strain and taking into account corresponding Young's modulus and initial yield stress (Table I).