Sheath parameters for non-Debye plasmas: simulations and arc damage

This paper describes the surface environment of the dense plasma arcs that damage rf accelerators, tokamaks and other high gradient structures. We simulate the dense, non-ideal plasma sheath near a metallic surface using Molecular Dynamics (MD) to evaluate sheaths in the non-Debye region for high density, low temperature plasmas. We use direct two-component MD simulations where the interactions between all electrons and ions are computed explicitly. We find that the non-Debye sheath can be extrapolated from the Debye sheath parameters with small corrections. We find that these parameters are roughly consistent with previous PIC code estimates, pointing to densities in the range $10^{24} - 10^{25}\mathrm{m}^{-3}$. The high surface fields implied by these results could produce field emission that would short the sheath and cause an instability in the time evolution of the arc, and this mechanism could limit the maximum density and surface field in the arc. These results also provide a way of understanding how the"burn voltage"of an arc is generated, and the relation between self sputtering and the burn voltage, while not well understood, seems to be closely correlated. Using these results, and equating surface tension and plasma pressure, it is possible to infer a range of plasma densities and sheath potentials from SEM images of arc damage. We find that the high density plasma these results imply and the level of plasma pressure they would produce is consistent with arc damage on a scale 100 nm or less, in examples where the liquid metal would cool before this structure would be lost. We find that the sub-micron component of arc damage, the burn voltage, and fluctuations in the visible light production of arcs may be the most direct indicators of the parameters of the dense plasma arc, and the most useful diagnostics of the mechanisms limiting gradients in accelerators.


I. INTRODUCTION
Comparatively little is known about the vacuum arcs and gradient limits that are important in determining the cost and overall parameters of large linear accelerator facilities. Vacuum arcs are involved in many fields, from particle accelerators, plasma devices, high power switching, surface coating and a variety of laboratory and commercial applications, and these arcs have been under study for many years. Nevertheless, the properties of these dense plasmas are not well understood, although the general behavior of these arcs has been known and under study since 1901 [1][2][3][4][5] and these plasmas seem to limit the performance of both major accelerator and tokamak projects and facilities [6][7][8]. In part, the reason for this situation is that the arcs are small, and many parameters (which are individually hard to measure) evolve very rapidly over a very wide range. Theory and modeling are complicated by the large number of mechanisms that seem to be involved in arc evolution and high density plasmas, which requires a non-Debye analysis of basic properties.
While arc damage has been measured and catalogued for over a hundred years, there has not been any clear correspondence between specific types of arc damage and the past or subsequent behavior of the arc. We argue in this paper that the causes of arc damage are due to the high density, high surface electric field plasma that produces a high plasma pressure in the liquid metal surface underneath the plasma arc. This pressure produces very small scale structures, at or below a few hundred nm. We have developed a self consistent model of arc evolution and show that these structures may or may not survive the subsequent cool-down of the liquid surface, however the cool down itself also seems to produce cracks with small scale structures [9][10][11][12]. This damage can produce future breakdown events.
Recent work has shown that the development of the arc can be explained by two mechanisms: 1) mechanical failure of the solid surface due to Coulomb explosions caused by high surface fields [12], and, 2) the development of unipolar arcs [13], that can act as virtual cathodes and produce currents that can short the driving potential. Once an arc starts, the surface electric field, and field emission increase, increasing ionization of neutrals, causing an increase in the plasma density. The density increase decreases the Debye length and causes an increase in the surface electric field, thus both the electric field and the density increase exponentially with time, roughly described by the arrow in Fig. 1. PIC simulations of the unipolar arc model for vacuum arcs relevant to rf cavity breakdown [9][10][11][12] show that the density of plasma formed above the field emitting asperities can be as high as 10 26 m −3 . The temperature of such plasma is low, in the range of 1 − 10 eV.
These high densities can make the Debye screening length, become smaller than the mean inter-particle distance, or the number of particles in the Debye sphere, to become less that unity. This implies the failure of the ideal plasma approximation, as well as most of the assumptions used in simple calculations. Processes in such a dense plasma can be affected by three body particle collisions so that the Particle In Cell (PIC) method which relies on a simple collisional model, with two body collisions, becomes inappropriate, as shown in Fig 1, where the arrow shows the approximate range of parameters for evolving arcs from Ref [11,12], as well as the approximate region of validity for PIC and Molecular Dynamics (MD) codes.
In this paper we calculate the parameters of the surface environment underneath the plasma sheath for the high density plasma conditions using direct two-component MD simulations where the interactions between all electrons and ions are computed explicitly. Although MD simulations have limited space and time scales their results can be considered as the lower level output for the multiscale approach.
Equilibrium and non-equilibrium nonideal plasmas have been studied extensively by MD in the past several decades [14][15][16][17][18]. Nevertheless there are few studies of the spatially inhomogeneous systems such as electric double layers or plasma sheath. In this paper we report on the first results for MD simulations of the nonideal plasma sheath near a metallic surface.
We can compare modeling with experimental measurements, but both the modeling and experiments are somewhat indirect. An understanding of the surface fields, combined with plasma density give estimates of the plasma pressure that can be compared with experimental data. It is possible to make indirect experimental measurements of the plasma pressure by comparing the linear dimensions of structures seen in the surface with estimates obtained from comparing the surface tension of the liquid metal with the plasma properties. We describe how the scale of structures frozen in damage as the arc cools, can be used to set limits on the plasma properties.

II. SIMULATION TECHNIQUE
The two-component fully ionized electron-ion plasma is considered. Neutral atoms are not taken into account which can affect relaxation times at relatively low plasma densities when the density of neutrals is high enough. It should not, however, affect the stationary distribution of charges. In the present work the simulations are restricted to the singly ionized plasma with Z = 1.
The electron-electron and ion-ion interactions are given by the Coulomb potential. For electrons and ions it is modified at short distances to account for quantum effects. The equation below assumes a Gaussian wave function for an electron where the r 0 parameter that equals to 0.21 nm in out case to match the ionization energy for copper at r = 0: U (0) = −7.73 eV (see Fig. 2). The similar interaction model was used e.g. in [18][19][20] for simulations of ionized metallic clusters. More accurate electron-ion and electron-electron interaction models are discussed e.g. in [21,22] although they seem to be redundant for this particular case. In fact the results are weakly dependent on the short distance part of the potential as the change of the U (0) value from 7.73 eV to 5.1 eV does not change the results within simulation accuracy. The Leap-Frog integration scheme is used to solve the classical equations of motion for electrons and ions [23]. The method takes into account the conservation of the total energy of the finite system, as long as there is no external potential. To follow the electron dynamics, time steps of 0.001 − 0.01 fs were taken to calculate the time evolution.
The general simulation scheme follows the method described in [24] and shown in Fig. 3. First an equilibrium MD trajectory is calculated for the system at given density and temperature using the nearest image method (periodic boundary conditions) for all dimensions. The simulation box size and other parameters are summarized in Table I. The Langevin thermostat [25] is used initially to bring the system to an equilibrium state while it is switched off for a production run. Then the system becomes adiabatic which ensures that all thermodynamic quantities are conserved in average. The ion mass is set to be equal to the electron mass for better mixing of ionic trajectories at this phase. The nonideality parameter, Γ, is the ratio of the average Coulomb potential energy and the average kinetic energy per electron [26].  T is the initial electron temperature, ne is the initial number density of electrons (or ions), Lz is the transversal simulation box size, Lx is the longitudinal simulation box size, Ni is the number of ions which is equal to the initial number of electrons, Γ = e 2 (4πne/3) 1/3 /(4π 0kB T ) is the nonideality parameter, Θ = 2 (3π 2 ne) 2/3 /(2mekBT ) is the degeneracy parameter, λD is the Debye length.
At the second phase the particle positions and velocities at particular time moments are taken from the equilibrium trajectory to be used as the initial states for nonequilibrium simulations of the plasma sheath. The interval between those points should be large enough for the initial states to be statistically independent from the microscopical point of view. However, all these states correspond to the same macroscopical conditions as they are taken from a single equilibrium trajectory. Then a bunch of trajectories is computed starting from the given ensemble of initial states and the results are averaged over the ensemble.
In order to study the plasma sheath, the XY plane at z = 0 axis is considered as a metallic surface whereas a reflecting wall is introduced on the other side of the box at z = L z . The periodic boundary conditions are still applied for transverse dimensions x and y. When an electron crosses the surface it is always meant to be absorbed. Therefore it is removed from the system and the overall surface charge is incremented by its charge q surf ← q surf − e.
A non-zero surface charge produces an electrostatic field which influence the particles. where σ = q surf /(L x L y ) is the surface charge density and L x , L y are the box sizes in the transverse dimensions. Assuming L x = L y = 2a one can obtain (see e.g. [27]) where E z is the longitudinal component of the electric field. It can be shown that Eq. (4) tends to the field expression of a uniformly charged plane E z = σ/(2 0 ) as z → 0 and to the Coulomb field E z = σa 2 /(π 0 )/z as z → ∞.
It is important to use Eq. (4) in the simulation with the given boundaries instead of the field of a charged plane E z = σ/(2 0 ) as the later cannot be screened by plasma particles at a large distance. As the surface field grows it starts to repel electrons from the surface until the stationary state is reached. We do not compute dynamics of ions at this phase as the ions are too heavy to contribute to the simulation results at the electron time scale. At the same time the ions are movable at the equilibrium trajectory that is used for generation of the initial states. Thus the averaging over an ensemble means the averaging over different configurations of ions. In a real system the number of ions will vary with time, due to ions entering and leaving the plasma from their thermal motion and self sputtering. Because the ion velocity is 340 times smaller than the electron velocity, this process is very slow, and we have neglected these effects. This is equivalent to assuming that the self sputtering coefficient for copper ions is near unity.
We have checked that the final results are independent of the simulation box size. If the box is doubled the deviation of the results are in within the statistical errors.
The thermodynamics parameters was maintained in the course of simulation. It was found that the overall electron temperature deviates in the range of 1 − 10% due to absorption of the most energetic electrons to the surface.

III. SIMULATION RESULTS AND FIT FORMULAS
Typically the relaxation of the electric field is observed for about 1 ps (see Fig. 4). The development of the electron profile is shown in Fig. 5. The stationary density profiles obtained after the relaxation are shown in Fig. 6. As the ions does not move, their distribution mimics the uniform distribution obtained from the equilibrium trajectory with full periodic boundary conditions. On the contrary, electrons form the well pronounced layer of plasma near the surface with a positive charge which we consider as the plasma sheath.
The plasma charge density profile is given by the difference between the ion and electron densities as presented in Fig. 7 in the semilogarithmic scale. It is seen that starting from the surface the charge density σ(z) decays exponentially which agrees with the Debye approximation. At high densities, however, the exponential decay is preceded by a non-exponential area. This regions makes difference between calculation of the sheath length from the slope λ exp of the exponential decay σ(z) ∼ e −z/λexp and from the distance λ at which the charge density decreases at the value of e = 2.71 (σ(λ) = σ(0)/e) as illustrated in Fig. 7c.
Both quantities λ exp and λ are presented in Fig. 8 depending on the plasma density and temperature. It is seen that λ exp follows the Debye-like dependence on density λ exp ∼ n −1/2 e whereas the real sheath length λ scales with a slightly different exponent.
The fits for MD data are T e = 10eV : Fig. 9 shows how the ratios λ exp /λ D and λ/λ D , depend on plasma temperature, electron density and the nonideality parameter. While Fig. 9a shows the dependence of the screening length on plasma temperature, Fig. 9b shows that the values of λ/λ D for different temperatures are close to each other when plotted versus the parameter Γ. This implies that the nonexponential charge density decay and the difference between λ and λ D are primarily a function of the plasma nonideality, defined as the ratio of average electrostatic potential energy divided by average kinetic energy.
The values of the electric field at the metal surface are presented in Fig. 10 depending on both temperature and density. The solid line corresponds to Eq. (2) from [13] where M i and m e are the masses of electron and ion. If the Debye radius in Eq. (7) is substituted by the MD result (5) or (6) Fig. 11 shows the plasma potential calculated using the simple relation of φ = E/λ where both the surface electric field E and the sheath length λ are obtained from MD simulations. A more rigorous result can be found by integration of the electric field distribution in plasma but it requires a more accurate evaluation of the space charge away from the sheath area and will be the subject of future work.

IV. EXPERIMENTAL ESTIMATES OF PLASMA PARAMETERS
The internal parameters of the arc, the surface field and the metal under it, are not directly accessible, but experimental measurements can provide some indirect evidence of the internal structure and the active mechanisms. Three phenomena are, in principle, sensitive to these numbers: 1) the mechanism that limits the exponential increase of the density and electric field with time, 2) the properties of the metal surface (melting point, self sputtering yield, etc.) may determine the burn voltage through the mechanism of self sputtering, and, 3) the scale of surface damage frozen into the surface should be sensitive to the plasma pressure that created it. We discuss these briefly.

A. Field Emission Driven Plasma Instabilities
The surface fields and plasma densities described above can be very high and it has been shown that these fields and densities increase exponentially with time during the evolution of the discharge [10][11][12]. The magnitude of these fields suggests that field emission over the entire active area could short out the sheath, causing an instability or oscillation in the plasma limiting this exponential increase.
As the surface field increases above 2 GV/m, it becomes possible for field emission currents to short out (or significantly reduce) the sheath field in times on the scale of nanoseconds. Assuming this occurs, we would expect that the sheath would rapidly reestablish itself due to the short collision time, the comparatively large plasma volume, and the high plasma density, and this behavior would produce fluctuations in the visible emission of the arc and fluctuations in the thickness of the sheath. This phenomenon could be described as the plasma "bouncing" on the metal surface. It is known that arcs are unstable, and fluctuations of this sort have been described by Jüttner [3] and Anders, see Fig 3.22 of Ref. [2]. In rf accelerator cavities, we see oscillations in visible light detected by a phototube with a frequency of approximately 200 MHz that could be due to this effect, see Fig 12 and Fig 13. We can understand the parameter range involved by estimating the current density required to short the sheath in a given time.
where ∆t = 1 ns implies currents of about 20 MA/m 2 and σ is the surface charge density. We can approximate Fowler-Nordheim field emission expression for current density at low fields [28] with, where the current density i is expressed in A/m 2 and the electric field is in V/m. Thus, when electric fields of 3 GV/m appear on the surface they will produce field emission currents capable of shorting the sheath in 1 ns. Assuming the plasma takes a few ns to return to the original density this would imply an instability with a time constant of a few ns. Instabilities in arc evolution are a well studied phenomenon. These strong oscillations may be related to the "ecton" model Mesyats has described, where he assumes micro-explosions with a timescale of 10 −8 s [4]. The characteristic, discontinuous "chicken track" traces on the interior of tokamak cavities could also be driven by these instabilities [2,5]. These instabilities are the subject of a further study, and will be be reported elsewhere.
We note that the field emission current densities discussed here, multiplied by the areas of melted copper in 805 MHz cavity damage spots, on the order of 2.5 ×10 −7 m 2 , would produce currents on the order of 4 A, roughly equal to the shorting current expected in 805 MHz breakdown events, see Fig 1 of Ref [12].

B. Material Dependence
We expect some dependence of the sheath potential and sheath parameters on the properties of the surface material. It has been shown that the evolution of the plasma is primarily driven by field emitted electron beams at high electric field and self-sputtering of surface material driven by ions falling through the sheath potential [10][11][12]. Self-sputtering produces a flux of neutral atoms that can raise the plasma density and also fuel the plasma, permitting long plasma lifetimes [5]. Numerical simulations of sputtering yields show that this mechanism is very sensitive to the sheath potential, ion charge distribution, surface (melting) temperature [29], and even grain orientation [30].
Anders, in Section 3.7 of Ref [2] explains the "burn voltage" of arcs, i.e. the voltage drop at the cathode, in terms of the Cohesive Energy Rule, where the cohesive energy of the cathode material is essentially the binding energy of the surface atoms. The larger the cohesive energy, the larger the burn voltage (proportional to the sheath potential) required to maintain a plasma. (The burn voltage is related to the sheath potential, but not equal to it, since electrons emitted from the plasma are not necessarily emitted at the plasma potential.) We believe that the mechanism that correlates the cohesive energy to the burn voltage is likely to be self-sputtering, which is determined by interatomic bonding. This data suggests that the sheath potential should be related to the atomic bonding energy, since we have shown that the sheath potential is primarily a function of the plasma temperature and weakly dependent on the density (see Fig. 11), we assume that this the plasma temperature is primarily involved. Measurements relating the burn voltage to plasma temperature (and perhaps crystal orientation) in different materials might explore this correlation.

C. Plasma Pressure and Surface Damage
We believe that the nature of the surface damage can provide information on the parameters of the sheath and the arc. The plasma ion flux hitting the surface should rapidly melt the top few layers of the surface. The plasma pressure pushing on the liquid metal surface can generate a Tonks-Frankel like instability [31], and uneven surfaces produced by this instability will be opposed by the surface tension force, which will tend to flatten the surface. As the dimensions of this instability become smaller, the surface tension force becomes more dominant, producing a correlation between the plasma pressure and the spatial scale of damage. The experimental problem is that surface tension will tend to smooth over the whole melted area when the liquid surface cools, making the melted areas polished and smooth. Our approach is to look for the smallest scale structure visible in arc damage, and assume that cooling has been rapid enough to preserve some evidence of the plasma pressure. For the scale of damage we observe ( ∼ 0.2 × 10 −6 m) thermal decay times are on the order of, where D is the thermal diffusivity constant, approximately 1.1 × 10 −4 m 2 /sec, implying times on the order a few ns or less, depending very strongly on the geometry of both the material and the details of the heat pulsing. We find experimentally that arcs moving in a transverse magnetic field produce the most fine structure. The pressure exerted on the surface by the plasma is due to: a) the plasma, and b) the electric field. The thermal plasma pressure is due to ion impacts, p i = nkT , and the electric field pressure is defined by, p E = 0 E 2 /2 [27], the total is then, p = p i − p E , since the ions push on the surface and the electric field, generally much smaller, pulls. In the limit of small E and T , the ion pressure can be a function primarily of the sheath potential, p i = neφ. If, due to a variety of reasons, the pressure is unevenly applied, it will produce a deformation in the liquid surface that is opposed by surface tension, see Fig 14. The approximate scale of these effects is set by the equilibrium radius, r, where the radius where the surface tension is balanced by the plasma pressure can obtained by equating the surface tension force around the circumference with the pressure over the whole area of a hemispherical bubble [32], with γ equal to the surface tension constant, approximately 1 N/m, depending on temperature, giving r ∼ 2γ/p [33]. For small structures it has been shown by Tolman that this expression should be corrected by a factor δ, using the expression, where δ is the Tolman length [34]. The Tolman length is generally evaluated using Molecular Dynamics, and estimates vary from tenths of molecular dimensions to hundredths of atomic dimensions. For radii, r, on the order of 100 nm this correction is not significant.
There are many types of arc damage that have been seen SEM images [2,5,35]. The damage from a single event is generally circular, in the range of 5 -200 µm in diameter, and frequently craterlike with a raised rim. The damage usually shows signs of melting. If the surface has absorbed significant energy, fine structure from the arc can be lost as the metal solidifies, However, if the arc deposits little energy to the surface or cools quickly, for example in SEM images of damage, Fig 16, a) from 201 MHz rf coupler, and b) from arc damage from Castano [35] and images from laser damage [13], we find complex structures on the scale of 100 -300 nm, which are not seen in arc damage where large amounts of energy (∼1 J) were present. We assume that if large amounts of energy are transmitted through an arc crater there is less small scale structure, consistent with high stored heat keeping the metal liquid until the surface tension smoothed off the surface. Classic unipolar arc tracks [5] (where magnetic fields move the arc in rapid retrograde motion) are associated with more fine structure, consistent with faster liquid cool down preserving this fine structure.
Simulations of unipolar arcs using PIC codes [12] have shown that the plasma potential seems to stay approximately 50 to 75 V during the development of the arc, thus the variation in plasma pressure is primarily due to the plasma ion density. Schwirzke showed that unipolar arcs could produce holes 5 times deeper than their diameter ( 0.7 µm) [13]. If we assume that these structures grew from craters with r ≤ 0.2 -0.35 µm, and the plasma potential, φ was 50 V, this would imply that the density of the plasma had to be at least 1 − 4 × 10 24 m −3 , see Fig 15. This is consistent with estimates made from the PIC code, which would not be expected to be reliable at these high densities.
The primary arc damage that results in further high enhancement factors and further breakdown events is likely due to this sub-micron damage, coming either from the plasma pressure itself producing a turbulent surface if it can quickly cool, or cracks produced when the large molten area beneath the arcs cools from the melting point of copper to room temperature leaving a network of surface cracks. The production of high enhancement factors in surface cracks has been demonstrated in Ref [10,11]. The sub-micron component of arc damage thus appears to be both the most direct indicator of the internal parameters of the arc plasma, and (when cracks and crack junctions are considered) the most likely to produce further breakdown events due to high enhancement factors.

V. CONCLUSIONS
We used molecular dynamics simulations to study the non-ideal plasma sheath at a metal surface for conditions we believe are appropriate to those found in accelerator cavities or unipolar arcs. The simulations started from the uniform equilibrium plasma state. Then the relaxation of the electron density profile with formation of the plasma sheath was observed. The relaxation time was found to be of the order of ∼ 100 fs. It was shown that the plasma sheath length depends on the electron number density in a slightly different way than the usual expression for the Debye radius due to a non-exponential charge profile at short distances. The values of the sheath length and the surface field were obtained for two values of temperatures and a wide density range with the non-ideality parameter Γ = 0.1 − 4. We compare the MD results with the contemporary theoretical models and with experimental data from damage. When we compare the plasma conditions that would result from these sheaths with data we find damage consistent with the high plasma pressures implied by the MD and PIC results.
We find that the high density plasma these results imply and the level of plasma pressure they would produce is consistent with the spatial scale of arc damage in rf cavities, in examples where the arc would cool before this structure would be lost. It appears that the sub-micron component of arc damage is both the most direct indicator of the internal parameters of the arc plasma, and, in the case of cracks, the most likely source of further breakdown events due to high enhancement factors. The high surface fields implied by these results could produce field emission that would short the sheath and cause an instability in the time structure of the arc. The relation between self sputtering and the burn voltage is not well understood but the two seem to be closely correlated. We find that the sub-micron component of arc damage, the burn voltage, and fluctuations in the visible light production of arcs may be the most direct indicators of the sheath parameters of the dense plasma.    [13]. Crosses correspond to the Eq. (2) from [13] where the Debye length is replaced by the sheath lengths λ obtained from MD (Fig. 8).