Multi-diagnostic experimental validation of 1d3v PIC/MCC simulations of low pressure capacitive RF plasmas operated in argon

The particle-in-cell/Monte Carlo collisions (PIC/MCC) simulation approach has become a standard and well-established tool in studies of capacitively coupled radio frequency (RF) plasmas. While code-to-code benchmarks have been performed in some cases, systematic experimental validations of such simulations are rare. In this work, a multi-diagnostic experimental validation of 1d3v electrostatic PIC/MCC simulation results is performed in argon gas at pressures ranging from 1 Pa to 100 Pa and at RF (13.56 MHz) voltage amplitudes between 150 V and 350 V using a custom built geometrically symmetric reference reactor. The gas temperature, the electron density, the spatio-temporal electron impact excitation dynamics, and the ion flux-energy distribution at the grounded electrode are measured. In the simulations, the gas temperature and the electrode surface coefficients for secondary electron emission and electron reflection are input parameters. Experimentally, the gas temperature is found to increase significantly beyond room temperature as a function of pressure, whereas constant values for the gas temperature are typically assumed in simulations. The computational results are found to be sensitive to the gas temperature and to the choice of surface coefficients, especially at low pressures, at which non-local kinetic effects are prominent. By adjusting these input parameters to specific values, a good quantitative agreement between all measured and computationally obtained plasma parameters is achieved. If the gas temperature is known, surface coefficients for different electrode materials can be determined in this way by computationally assisted diagnostics. The results show, that PIC/MCC simulations can describe experiments correctly, if appropriate values for the gas temperature and surface coefficients are used. Otherwise significant deviations can occur.


Introduction
Among the variety of low temperature plasma sources available, low pressure radio frequency capacitively coupled plasmas (RF CCPs) are of high technological relevance [1][2][3]. They are used for several applications, e.g. etching and deposition processes. Despite their widespread use, some fundamentals of their operation are not understood and, therefore, represent an area of active research and debate. In particular, power absorption mechanisms and plasma-surface interactions are not understood in many cases [4][5][6][7][8][9].
To reveal such fundamentals of RF CCP operation, a synergistic combination of experiments and simulations is typically required, since the simulations provide access to plasma parameters that are difficult or impossible to measure with a high spatial and temporal resolution. If validated experimentally, computational simulations can be used for plasma source design and process development. In this way significant time and costs can be saved that would otherwise occur to build prototypes and perform a huge number of experimental tests [10].
While being self-consistent, such simulations typically require a number of input parameters such as surface coefficients that describe the material, particle and incident energy dependent probability of secondary electron emission from boundary surfaces per incident particle [22][23][24][25][26][27][28], and the probability for electron reflection at such surfaces [29][30][31]. If complex reactive plasmas are investigated computationally, absorption, reaction, and etch/sputter probabilities for different particle species at boundary surfaces are also required [25,[32][33][34]. Finally, the gas temperature is used as an input parameter as well in many cases, although self-consistent implementations exist [28,35].
A constant, uniform value for the gas temperature T g is usually adopted over the whole range of operation parameters (typically 300 K T g < 400 K [11,26,36,37]).
In the simulations, based on the ideal gas law, the gas density n g depends on the pressure p (also an input parameter) and the gas temperature according to n g = p/(k B T g ), where k B is the Boltzmann constant. Thus, gas heating depletes the gas density. This, in turn, affects the mean free path of electrons and ions according to λ e,i = 1/(n g σ e,i ), where σ e,i is the crosssection for collisions of electrons and ions, respectively. In this way the gas temperature can affect the spatio-temporal electron impact ionization dynamics, the mode of discharge operation, and the shape of ion flux-energy distribution functions at boundary surfaces.
To obtain realistic simulation results, depending on the discharge conditions, the required input parameters include heavy particle (ion, neutral) and electron induced secondary electron emission coefficients (SEECs) as well as probabilities for electron reflection at boundary surfaces. However, such coefficients are often unknown or suffer from large uncertainties [38]. In many cases, the plasma facing materials and their conditions are unknown, due to unknown effects of the plasma on the surface and a lack of in-situ surface diagnostics. Previous works have demonstrated that some of these parameters, for example SEECs, can be determined under plasma exposure by computationally assisted diagnostics [39] and combinations of current measurements and modeling [40,41]. Alternatively, they can be calculated based on ab initio models [42][43][44][45]. Typically, however, such surface coefficients are either neglected, guessed, or taken from beam-experiments performed under ultra-high vacuum conditions without plasmas that typically modify the surface. Depending on the discharge conditions, this way of handling surface coefficients and the gas temperature represents a major source of uncertainty with respect to the simulation results. Typical values used for the surface coefficients are 0.2 for the electron reflection probability [26,46] and γ 0.1 for the ion induced SEEC for clean metallic electrodes [2,[47][48][49][50]. Other surface processes/coefficients are often not considered.
Systematic comparisons between the results of different simulation codes and between simulation and experimental results are required to verify the computational implementations and to validate the discharge models used. One of the first comparisons between measurements performed in an RF CCP and PIC/MCC simulation results was done by Vahedi et al in argon gas [51]. The authors found good agreement between electron energy distribution functions obtained from their code and those measured by Godyak et al [52], but the plasma densities differed by a factor of about two at pressures of about 10 Pa. Rakhimova et al [47] compared hybrid simulation results (PIC/MCC approach for electrons combined with a fluid model for the ions) and experimental results for an RF CCP driven at 13.56 MHz and 81 MHz. The results obtained from the different approaches show deviations of almost one order of magnitude in plasma density at 13.56 MHz, depending on the modelling technique and the choice of the ion induced SEEC γ. With γ = 0.1 a good agreement between PIC simulation and experimental results was found. Braginsky et al [53] compared PIC/MCC simulation results of an RF CCP operated in argon at a low driving frequency of 1.76 MHz with measurements of the plasma density, the RF current, and the spatio-temporally resolved electron dynamics at selected values of the gas pressure and RF power. They concluded that an accurate implementation of plasma-surface interactions in the simulation is essential to obtain good agreement of the computational results with measurements. Derzsi et al [54] have carried out experimental and simulation studies of RF oxygen discharges driven by tailored voltage waveforms [55] and have compared the discharge power, the self-bias voltage, the ion flux, as well as the ion flux-energy distribution function over a wide domain of the gas pressure and the driving voltage. These studies have found a reasonable agreement between the measured and computed characteristics provided that a proper value for the quenching rate of singlet delta oxygen molecules at the electrode surfaces was adopted. In a subsequent study [56] experimental observations of the electron power absorption modes and mode transitions were also found to be correctly reproduced by this simulation code.
Turner et al performed a computational benchmark study for capacitive RF discharges operated in helium [57] by comparing the results of five different, independently developed PIC/MCC codes to each other under selected discharge conditions. While a very good agreement between the results of the different implementations of the PIC/MCC scheme was obtained, the authors stated that these results deviate considerably from those available in the literature. These differences were attributed to the simple physical model (i.e. no emission of particles from surfaces, no effects of excited states) used for the code benchmarking.
Recently, RF phase resolved electron density measurements using a hairpin probe were presented [58]. These data are expected to be an excellent basis for an experimental validation of PIC/MCC simulations.
In this work, we provide results of a systematic multidiagnostic experimental validation of 1d3v electrostatic PIC/MCC simulation results for RF CCPs operated in argon gas at 13.56 MHz. For this, we perform measurements in a specifically designed reference reactor, which is geometrically symmetric and, thus, suitable for comparison with 1d3v PIC/MCC simulation results for plane parallel RF CCPs [47,[59][60][61][62]. The simulations include a careful treatment of plasma-surface interactions. In the experiment, stainless steel electrodes are used and multiple diagnostics are applied to measure the driving voltage waveform at the powered electrode, the central plasma density, the ion flux-energy distribution function at the grounded electrode, the spatio-temporally resolved electron impact excitation dynamics, and the gas temperature as a function of the gas pressure and driving voltage amplitude. We find that the simulations can describe the experiments quantitatively correctly, i.e. good agreement is found for all plasma parameters considered for all gas pressures and driving voltage amplitudes investigated, if the correct gas temperature and appropriate surface coefficients are used in the simulations. At high gas pressures, the gas temperature is found to increase considerably beyond room temperature. From a systematic comparison of experimental and computational results the effective electron reflection probability at the electrodes is found.
To the best of our knowledge this work corresponds to one of the first successful systematic and multi-diagnostic experimental validations of 1d3v PIC/MCC simulations of CCPs for the simplest scenario of a single frequency, geometrically symmetric argon discharge. While code-to-code benchmarks and qualitative comparisons between simulation and experimental results exist, such quantitative validation efforts are rare and important. Good agreement between results of different codes and the ability of simulations to qualitatively reproduce experimentally observed parameter trends is often not sufficient to rely on such simulations for plasma process development and understanding fundamental phenomena. The key outcome of this work is, thus, the quantitative experimental validation of the simulation results itself as well as its sensitivity on distinct input parameters rather than any plasma physical effects. In fact, such validation studies are required for a variety of discharge conditions and plasma sources. This paper might trigger further investigations of this type.
The manuscript is structured in the following way: in section 2, the methods of our investigations are introduced including the experimental setup, experimental diagnostics, and the simulation method. The results of our study are presented and discussed in section 3. Finally, the work is summarized and conclusions are drawn in section 4.

Experimental setup
The experimental validation of 1d3v PIC/MCC simulations of plane parallel RF CCPs requires a plasma source that resembles the geometry assumed in the simulation. In particular, care has to be taken to ensure that a symmetric discharge establishes in the source. This requires having equal powered and grounded surface areas around the plasma, i.e. the design should be geometrically symmetrical. Constructing a geometrically symmetric RF CCP is a challenge. Essentially all commercial and research reactors are geometrically asymmetric, since the chamber walls are typically grounded and, thus, the grounded surface area is larger than the powered surface area. Even if a dielectric plasma confinement is used, there will typically be a capacitive coupling of these dielectric surfaces with external grounded surfaces, which will make the discharge asymmetric. In asymmetric reactors, a DC selfbias voltage builds up. In a truly symmetric system this voltage should approach zero. This requirement can be viewed as a test condition for the success of the reactor design. The other requirement that must be fulfilled for reasonable comparisons of the experimental results with 1D simulation data is that edge effects in the experimental system should be negligible. This can be ensured by having an electrode diameter that significantly exceeds the electrode separation.
For these reasons and similar to previous work [59,63], we constructed a cylindrical reactor whose lateral walls are made of borosilicate glass (see figures 1 and 2). It consists of a glass cylinder of 17 cm height. DN120 flat glass flanges are located at the top and bottom of the cylindrical reactor and two KF40 glass flanges, that oppose each other, are connected to the main cylinder from the side according to figure 2. The KF40 flanges are directly melted onto the cylinder wall and connected to the inside of the cylinder by 15 mm diameter holes. The top and bottom DN120 flanges are closed by stainless steel lids and FKM O-rings. A KF40 flange is welded onto the bottom lid and is used to connect pressure gauges and the retarding field energy analyzer (RFEA) electronics. Two planar stainless steel electrodes of 12 cm diameter are placed inside the reactor according to figure 1. The inner diameter of the glass cylinder is 120 mm ± 0.5 mm so that the gap between the electrodes and the glass wall varies between 0 mm and 1 mm depending on the azimuthal position. Each electrode is connected to its adjacent lid via a copper rod, which holds the electrode in place and provides an electrical connection to the respective lid. Adjusting the length of these rods allows to control the electrode gap. One of the electrodes is driven by the RF voltage, while the other  electrode is grounded. In this way a geometrically symmetric RF CCP is realized, in which essentially no DC self-bias is generated (the self-bias voltage is less than 5% of the driving voltage amplitude).
The reactor is evacuated by a two-stage pumping system (Leybold Divac 2.5E and Turbovac SL80) via one of the two KF40 side flanges. The base pressure of the system is 5 × 10 −5 Pa. The operating pressure is set by controlling the pumping cross section via a gate valve (VAT HV-Schieber Reihe 14.0), which can be operated manually or by using a computer controlled stepper motor. The pressure in the system is monitored by three manometers mounted at the bottom KF40 flange, i.e. a cold cathode gauge (Pfeiffer Vacuum IKR 261) and two capacitance gauges with different ranges (MKS Baratron 627B 10 Pa, MKS Baratron 626A 1000 Pa). A mass flow controller (Bronkhorst F-201CV) regulates the flow of 6.0 purity argon gas into the chamber via the opposing KF40 side flange (flow: 0.5 sccm-2 sccm).
The driving voltage is provided by an RF generator (Advanced Energy Cesar 136) and applied to the top electrode via an impedance matchbox (Advanced Energy VarioMatch 1000). The voltage at the top lid (including the DC self-bias) is measured by a high voltage probe (PMK PHVS 662-L-RO) connected to an oscilloscope (LeCroy LT364). There is no noticeable voltage drop between the position where the voltage probe is attached and the powered electrode itself. This was checked by measuring the voltage at both points with two probes simultaneously in the vented chamber.
The plasma is operated in argon at gas pressures of 1 Pa to 100 Pa with driving voltage amplitudes ranging from 150 V to 350 V at 13.56 MHz. The electrode gap is fixed at 40 mm. This pressure range is chosen to validate the simulation results in the low pressure non-local as well as in a more local regime at higher pressure. This electrode gap and range of driving voltage amplitudes are selected to ensure a stable discharge operation within this pressure range. No parasitic plasma was observed on the backside of the electrodes for the conditions chosen.
Five experimental diagnostics are used: (i) a high voltage probe to measure the time resolved voltage drop across the discharge, (ii) phase resolved optical emission spectroscopy (PROES) to measure the spatio-temporally resolved electron impact excitation dynamics, (iii) a Langmuir probe to quantify the plasma density in the center of the reactor, (iv) a RFEA to obtain the ion flux-energy distribution at the grounded electrode, and (v) a tunable diode laser absorption spectroscopy (TDLAS) system to measure the gas temperature.
For PROES, we use an ultra fast gated ICCD camera (Stanford Computer Optics 4Picos) in combination with a telecentric lens (Thorlabs MVTC23013) and an interference filter (Andover Corporation 010FC37-25/750.4-D, FWHM = 0.914 nm). The camera is aimed at the center of the reactor chamber. It is triggered by the RF Generator. A delay generator (Stanford Research Systems DG535) is used to reduce the frequency of the trigger signal to 3/4 of the maximum camera input trigger frequency (200 kHz → 150 kHz) and to shift the trigger in time with respect to the driving voltage waveform. The combination of the camera and the lens results in a spatial resolution of more than 6 pixels per millimeter. A camera gate time of 1 ns is used. In this way space and time resolved measurements of a selected emission line are performed. We use the Ar-I 750.39 nm emission line and control the camera and the delay generator with a custom LabVIEW program, so that the acquisition of the PROES images is fully automated. The electron impact excitation rate from the ground state into the Ar 2p 1 state (lifetime: 22.2 ns [64]) is calculated from the measured plasma emission using a model described in [65]. As the threshold energy for this electron impact excitation process is 13.5 eV, the space and time resolved dynamics of energetic electrons above this energy threshold is measured.
To measure the plasma density, a modified version of the Langmuir probe system described in [66] is used. To keep the influence of the probe on the plasma as small as possible, we extended the tip of the probe by a ceramic capillary (length: approx. 60 mm, outer diameter: 1 mm, inner diameter: 0.5 mm). The probe wire is located inside this capillary. The auxiliary electrode of the original probe tip was connected via a wire to a piece of copper tape wrapped around the tip of the ceramic capillary. The capillary and the wire between the original probe tip and the copper tape are covered by Kapton tape. In this way, only the capillary is inserted into the plasma, while the bigger main probe body stays inside the side flange (see figure 1). We used a tungsten wire with a diameter of 50 μm which extends 3.2 mm into the plasma. The raw data collected by the probe system are evaluated using another custom LabView program. The electron density is determined via the 'i 2 -method' proposed first in [67] and refined in [68]. To apply this method no information on the electron temperature is needed for the determination of the plasma density. Small fluctuations of the potential produce only minimal changes of the probe current and hence of the slope of an i 2 vs V p plot, where i is the measured current and V p is the potential applied to the probe [67]. The noise level of the measured data did not allow the determination of the EEDF.
To measure the ion flux-energy distribution function at the grounded electrode a modified RFEA system (Impedans Semion) is used. In the original system, the RFEA sensor is placed in the middle of a holder, which is then placed on the electrode in the reactor. This was not possible in our reactor, because the diameter of the connector is larger than the diameter of the holes in the reactor walls. We also wanted to avoid using the holder, since it has a diameter of 70 mm, while the electrodes have diameters of 120 mm. This means that the holder would cover about 34% of the area of one of the electrodes and thereby change the plasma-facing surface material and reduce the electrode gap at the position of the holder. Thus, we designed an electrode in which the RFEA sensor is embedded in the center. In this way the ion flux-energy distribution function is measured at the electrode surface and most of the plasma-facing material remains stainless steel. Moreover, the electrode gap remains unchanged and laterally uniform. The original Impedans Semion data acquisition and control system was used.
TDLAS is used to determine the gas temperature. This is done by measuring the absorption profile of the transition Ar(1s 5 → 2p 6 ) at 772.376 nm. We use a Toptica fiber coupled laser head combined with its dedicated controller (head: DFB pro 100 mW, 772 nm + Fiberdock; controller: DLC pro). The laser head has a temperature and a current scan option. With the temperature change, a wavelength range of about 2 nm can be scanned mode hop free. The current scan spans a wavelength range of 0.05 nm and the laser line width Δν is less than 2 MHz. For all measurements, we set the temperature to 30 • C and use the current scanning mode. The experimental setup for TDLAS is shown schematically in figure 2. The laser beam is coupled into a single mode optical fiber with integrated beam splitter (Thorlabs TN785R2A1). The fiber coupling makes the system a lot more versatile compared to a system in which the laser has to be guided through open air. It also reduces the complexity of adjustment procedures and reduces the noise caused by vibration of the system. 90% of the laser power is transferred to a Fabry-Perot interferometer (FPI, Toptica FPI 100-750-3V0, 1 GHz free spectral range). Based on the FPI the wavelength step of the scanning laser is monitored. Moreover, it allows to verify that the laser wavelength is shifted linearly as a function of the current. For all the measurements, we used a laser scanning frequency of 10 Hz. The remaining 10% of the laser power is guided through a second fiber with a graded index collimator at the end (Thorlabs 50-780-APC) towards the plasma reactor. Before entering the reactor, the Gaussian shaped beam is again split according to a ratio of 9/1, so that the total laser power entering the discharge is of the order of 150 μW. In this way, absorption saturation effects are avoided. The beam passes the reactor in the middle and is detected by a photodetector (Thorlabs DET10N2) on the opposite side. Both the photodetector and the FPI are connected to a PC-based oscilloscope (PicoScope 6402C) which is controlled via LabVIEW to process and store the measured data on a computer. The oscilloscope is triggered by a signal of the laser controller with a frequency equal to the scanning frequency (here: 10 Hz). For each measurement, four data sets are recorded: (1) plasma and laser on, (2) plasma on and laser off, (3) plasma off and laser on, (4) plasma off and laser off. The sets are averaged up to 100 times. The data evaluation is done semi-automatically based on a custom LabVIEW program. Both the photodiode and the FPI data are loaded for the evaluation. The peaks in the FPI data are detected automatically and their positions on the horizontal axis are then used to convert the original unit of the axis (time) to wavelength. It is checked for each measurement, if the distance between all interference peaks is uniform. This shows that the wavelength changes linearly with the laser current and no mode hopping occurs. With the wavelength information from the FPI data, the measured absorption profile is calculated and plotted. A Gaussian fit to the data is performed automatically and yields the metastable density as well as the gas temperature. We observe no relevant changes in the values when we fit a Voigt instead of a Gauss function. The observed line width, therefore, is mainly a consequence of Doppler broadening, and pressure broadening can be neglected. The evaluation procedure is more extensively described in [69]. Information about the principles of TDLAS can be found in [70][71][72].

PIC/MCC simulation
Our simulation studies are based on a 1d3v electrostatic PIC/MCC code [20]. The code traces electrons and Ar + ions in a homogeneously distributed background gas of which the density is defined by the pressure (p) and the gas temperature (T g ). The electron impact cross section set (that comprises the elastic momentum transfer, a lumped excitation, and an ionization cross section) is adopted from [22], while for the Ar + ions the data from [73] are used. In the case of the electrons, all collision processes are characterised by isotropic scattering. In ionization processes the scattered and the ejected electrons are assumed to share the remaining kinetic energy equally. For the Ar + ions we take into account elastic collisions only and account for the anisotropic scattering by approximating the differential cross section with an isotropic and a backward channel as advised in [73].
Particles arriving at the electrodes are treated in the following way. (i) Electrons can be either absorbed or reflected. The latter process is characterised by an effective reflection coefficient, r eff , i.e. we do not distinguish between the various possible microscopic processes (elastic and inelastic reflection, and emission of secondary electrons, see e.g. [8,25,27,28,74]). This simplification is done to limit the number of input parameters of the simulation to keep the simulations as simple as possible and to allow finding combinations of input parameters that yield computational results in quantitative agreement with the experiment more easily. (ii) Upon the impact of ions at the electrodes we consider the emission of secondary electrons. At the ion energies relevant to this study, the potential ejection mechanism is foreseen to prevail, therefore we set the value of the ion induced SEEC to a fixed value, γ = 0.07, based on [22,39].
The simulations provide a number of plasma characteristics, e.g. the spatio-temporal distributions of particle densities and fluxes, and the flux-energy distributions of the charged particles impinging the electrode surfaces. To facilitate comparison with the experimental PROES measurements as well, the excitation rate from the ground state into the Ar 2p 1 state is also computed via accumulating the contributions to this rate by individual electrons along their trajectories (based on the electron impact cross section of this specific state [75]).
The parameters of the simulations are set in a way to comply with the stability and accuracy criteria of the PIC/MCC technique [76,77]. We use 12 000 time steps per RF period and 1800 grid points within the electrode gap. The number of superparticles (per species) is chosen to be ≈3 × 10 5 in order to achieve accurate results [78].

Results
Below, we present the results of our studies: the gas temperature, the central plasma density, the ion flux-energy distribution function at the grounded electrode, as well as the spatio-temporally resolved electron impact excitation dynamics are measured as a function of the gas pressure and the driving voltage amplitude at 13.56 MHz with an electrode gap of 40 mm. The measured gas temperature is used as input for our 1d3v PIC/MCC simulations. The other measured quantities are compared to the corresponding simulation results to perform a multi-diagnostic experimental validation of the computational results. In the following, each of these plasma parameters is discussed separately.

Gas temperature
As the gas temperature can affect the simulation results significantly, an accurate value must be known as input to the simulation as a basis for a meaningful experimental validation. Therefore, T g was measured in the center of the reactor by TDLAS as a function of the gas pressure and the driving voltage amplitude. Following any change of the pressure and the voltage amplitude, enough time was given to the system to stabilize. The measurements were repeated several times and excellent reproducibility was obtained within the error bars shown in figure 3, which shows the results of the temperature measurements. The gas temperature is found to depend weakly on the driving voltage amplitude, but to increase strongly as a function of the gas pressure from about 300 K at 1 Pa to up to 400 K at 100 Pa.
The dependence of the gas temperature on the driving voltage amplitude is more pronounced at high compared to low pressures. Increasing the voltage leads to a higher ion flux and mean ion energy at the electrodes at constant pressure. At high pressure, the ion mean free path is much shorter than the maximum sheath width and energetic ions, that are accelerated towards the adjacent electrode by the sheath electric field, collide frequently with the background gas and transfer energy to the gas. At higher voltage, such collisions happen more frequently and more energy is transferred to the gas [37]. At the lowest pressure of 1 Pa, the ion mean free path is of the order of 1 cm and, thus, comparable to the maximum sheath width (see figure 7). Therefore, very few collisions of energetic ions with the neutral gas happen and the effect of the voltage on the gas temperature is weaker at low pressures. Figure 4 shows the electron density in the center of the discharge obtained from Langmuir probe measurements and simulations performed with different input parameters for the surface coefficients and the gas temperature as a function of the gas pressure (1 Pa-20 Pa) at different driving voltage amplitudes. The values of the temperature include two constant values of 300 K and 350 K (often adopted in discharge simulations) and the measured gas temperature indicated as T TDLAS . At pressures exceeding 20 Pa the plasma confinement inside the cylindrical glass reactor was reduced due to the fact that the discharge penetrates into the side flange with the Langmuir probe inside. Also, an additional discharge around the Langmuir probe tip started to appear at pressures above 20 Pa. Such effects are not included in the simulation and, thus, no comparison of the measured and computationally obtained plasma densities can be performed at such high pressures.

Plasma density
Generally, the plasma density is found to increase as a function of the gas pressure. For each driving voltage amplitude, three different simulation results are compared to the measurements. These are obtained based on different input parameters for the effective electron reflection probability (r eff ), the ion induced SEEC (γ), and the gas temperature (T g ). The choices of the input values for the three simulation cases shown in figure 4 are: i. r eff = 0.2 γ = 0.07 T g = 350 K, 'base case', red triangles. ii. r eff = 0.2 γ = 0.07 T g = 300 K, blue triangles. iii. r eff = 0.7, γ = 0.07, T g measured), 'best case', green squares.
The choices of input parameters for the 'base case' are based on previous work [26,46] (r eff ), [22,39] (γ). For this 'base case' (red triangles) good agreement of the computational results of the plasma density with the measurements is only found at higher pressures. At low pressures of 1 Pa and 2 Pa strong deviations of up to one order of magnitude are found. The second set of computational results (blue triangles) is obtained based on the same effective electron reflection coefficient and ion induced SEEC, but using a lower gas temperature of 300 K, which agrees better with the measured gas temperature at low pressures (see figure 3). At low pressures, T g = 300 K results in a significant increase of the plasma density in the simulation under otherwise fixed input parameters, because the gas density and, thus, the collisionality are increased. Correspondingly, energetic electrons generated by sheath expansion heating at the electrodes ionize the background gas more likely before they leave the discharge by hitting the opposite electrode (see also section 3.4). However, this increase of the plasma density obtained by lowering the gas temperature is not enough to reproduce the measured values at low pressures. For the 'best case' simulation results (green squares) the effective electron reflection probability is increased to 0.7 and the measured values of the gas temperature as well as an ion induced SEEC of 0.07 are used. This scenario yields good agreement with the measured plasma densities at all pressures and driving voltage amplitudes. In the 250 V and 350 V cases, the measured electron density drops from 10 Pa to 20 Pa, while it continues to increase in the simulation. The experimentally observed decrease might be the result of an imperfect RF compensation of the probe at the highest pressure studied in this work, i.e. it is due to an experimental error.
We also tested if the same good agreement can be obtained based on other combinations of the two surface coefficients, γ and r eff , but this turned out not to be possible. Increasing γ leads to much higher plasma densities at higher pressures, but has only a weak effect at lower pressures. This is because ioninduced secondary electrons play an important role only under the conditions of a significant electron multiplication inside the sheath regions that is possible at high pressures only. At low pressures there is no such avalanche effect and the effect of γ on the plasma density is low. Increasing γ, therefore, does not allow to bring the results of the simulation 'base case' closer to the measured values of the plasma density. At low pressures, however, the plasma density is sensitive to r eff , which is not the case at higher pressures. This is caused by the long mean free path of electrons at low pressures, which can be even larger than the electrode gap. Thus, energetic beam electrons generated by sheath expansion heating at one electrode will propagate collisionlessly through the bulk and hit the opposite electrode during the local sheath collapse, where they are lost according to 1 − r eff (see section 3.4). Thus, under such conditions increasing r eff leads to a much better confinement of energetic electrons in the plasma and a higher plasma density. At higher pressures increasing this surface coefficients has almost no effect due to the short electron mean free path, i.e. such energetic electron beams cannot reach the opposite electrode [65]. Overall, the low pressure simulation results are not sensitive to γ and the high pressure results are not sensitive to r eff . For a known gas temperature, this behaviour allows to find r eff and γ for different materials and operating conditions within the accuracy limits of the simulation method and based on the comparison of simulated and measured data, since the simulation results depend uniquely on only one of these coefficients in a given pressure range.
This comparison of measured plasma densities to those obtained computationally shows that good quantitative agreement is found for the 'best case' set of simulation input parameters, i.e. for a distinct choice of surface coefficients and gas temperature. Figure 5 shows the ion flux-energy distribution function at the grounded electrode obtained experimentally and computationally based on the 'best case' scenario of input parameters (r eff = 0.7, γ = 0.07 and T g taken from the TDLAS measurements) for gas pressures ranging from 1 Pa to 10 Pa at an exemplary driving voltage amplitude of 350 V. For these input parameters good agreement between the measured and computationally obtained plasma densities was found. Simulation results smoothed to resemble the energy resolution of the RFEA (10 eV) are also shown. We contemplate that this relatively poor energy resolution is caused by the fact that we use our own cables to connect the RFEA sensor to its electronics. These cables pick up RF noise from the plasma, and, thus, the potential of the grids inside the RFEA sensor is not perfectly constant, but modulated, leading to a reduction of the RFEA energy resolution.

Ion flux-energy distribution function
Good agreement between the smoothed experimental data and the simulation results is found with a small, ≈5 eV shift between the two data sets. This is caused by the marginal  DC self-bias, which is generated in the experiment, since the reactor is not perfectly symmetric. Regardless of this small energy shift, the shape of the distribution function including the number and width of the lower energy peaks are reproduced correctly by the simulation.
In both, experiment and simulation, a high energy peak of the ion distribution function is present at low pressures of 1 Pa and 2 Pa. This peak disappears with increasing pressure. This happens, because the mean free path of the ions decreases with increasing pressure. Figure 6 shows how adjusting the input parameters of the simulation and smoothing the computationally obtained ion distribution function according to the experimental energy resolution improves the agreement between experiment and simulation at the lowest pressure of 1 Pa. The 'base case' scenario in the simulation (r eff = 0.2, γ = 0.07 and T g = 350 K) yields an ion distribution function that is shifted towards higher energies compared to the experimental results. Changing the simulation input parameters to the 'best case' scenario (r eff = 0.7, γ = 0.07 and T g taken from the TDLAS measurements) causes a shift of about 10 eV towards lower energies with respect to the 'base case' result. According to the simulation results this predominantly happens, since the mean electron energy is reduced at the higher value of r eff , because the electron confinement is improved and the mean electron energy is the result of a particle balance for the electrons. Consequently, the RF floating potential and the sheath voltage are also reduced in the case of the higher electron reflection probability. This explains the observed shift of the ion distribution function towards lower energies. Smoothing the 'best case' computational result yields a distribution function that is in good agreement with the experiment except for a small remaining shift along the energy axis due to the low DC self bias present in the experiment, but not in the simulations.
These results show that the choice of input parameters for the simulation (surface coefficients and gas temperature) that yielded optimum agreement with experiments for the central plasma density (see figure 4) also yields excellent agreement of the computationally obtained ion flux-energy distribution function at the grounded electrode with experimental results. While figure 6 only shows simulation results for selected combinations of simulation input parameters (r eff , γ, T g ), each of the input parameters was changed independently, while leaving the others constant. In this way the effects of individual input parameters on the ion distribution function were studied and the 'best case' scenario was found to yield optimum results. Figure 7 shows spatio-temporal plots of the electron impact excitation rate from the ground state into the Ar2p 1 state for a fixed driving voltage amplitude of 250 V and for different gas pressures of 2 Pa, 10 Pa, and 50 Pa. The electron energy threshold for this process is 13.5 eV. Thus, these distributions sample the density of high-energy electrons. Simulation results for the 'base case' (left column, r eff = 0.2, γ = 0.07 and T g = 350 K) and for the 'best case' scenario (right column, r eff = 0.7, γ = 0.07 and T g taken from the TDLAS measurements) as well as experimental results (middle column) are shown. Under all discharge conditions the plasma is operated in the α-mode [65,79], i.e. excitation occurs as a consequence of the energy gain of the electrons in the vicinity of the edges of the expanding sheaths at both electrodes. At high pressure, the electron mean free path is so short that the energetic beam electrons generated in this way cannot propagate into the plasma bulk, while this is possible at lower pressures. This mode of operation is observed experimentally and by the simulation in all cases. Except for the lowest pressure of 2 Pa, details of the experimentally observed spatio-temporal excitation rate plots are well reproduced by the simulation for the base and the 'best case' scenario. At 2 Pa, however, the maximum sheath width is clearly larger in the 'base case' simulation scenario compared to the experiment. In agreement with the results discussed in previous sections, this indicates that the plasma density is underestimated by the simulation in the 'base case' scenario. The 'best case' simulation scenario, which uses a higher electron reflection probability and a lower gas temperature, yields a maximum sheath width that is very similar to the one observed experimentally. This confirms that a realistic plasma density is obtained. In agreement with the results of previous sections, this shows that the simulation results are sensitive to these input parameters at low pressure in the non-local regime. Figure 8 shows the spatio-temporally resolved electron impact excitation dynamics at 1 Pa and 350 V as a result of the simulation 'base case' (left column), the experiment (middle column), and the simulation 'best case' (right column). In contrast to the higher pressure cases (see figure 7), there are even bigger differences between the result of the simulation 'base case' and the measurement, while the simulation 'best case' reproduces the experiment well. In particular the experimentally observed maximum sheath width is only reproduced by the simulation 'best case'. In the simulation, the generation of multiple electron beams during a single sheath expansion phase is observed in agreement with the results of Wilczek et al and Berger et al [12,80]. Due to the limited temporal resolution of the PROES method, this could not be observed in the experiment.

Spatio-temporal electron impact excitation dynamics
Similar to the results obtained for other plasma parameters, which were shown in the previous sections, good agreement between experiment and simulation at low pressures is again only found for the 'best case' simulation scenario, i.e. based on the selection of distinct surface coefficients and gas temperatures as input for the simulation.

Conclusion
Results of a systematic multi-diagnostic experimental validation of 1d3v electrostatic PIC/MCC simulations of RF CCPs operated at 13.56 MHz in argon at pressures between 1 Pa and 100 Pa as well as driving voltage amplitudes from 150 V to 350 V were presented. Measurements were performed in a custom built geometrically symmetric and cylindrical reactor to resemble the 1D geometry assumed in the simulation. The lateral walls of this reactor are made of glass and stainless steel electrodes of identical surface areas were used. Experimentally, the voltage waveform at the powered electrode, the gas temperature, the central plasma density, the ion flux-energy distribution function at the grounded electrode, and the spatiotemporally resolved electron impact excitation dynamics were measured as a function of pressure and voltage. These experimental results were compared for the same conditions with those obtained from PIC/MCC simulations, in which (i) the gas temperature (T g ), (ii) the effective electron reflection probability at the electrodes (r eff ), and (iii) the ion induced SEEC (γ), were variable input parameters.
Good quantitative agreement between the experimental and computational results was found for all pressures and driving voltage amplitudes, if the measured gas temperature and distinct choices of surfaces coefficients were used as input for the simulations, i.e. r eff = 0.7 and γ = 0.07. Based on extensive variations of the coefficients in the simulations, we estimate the accuracy of r eff , determined in this way, to be ±0.1. The literature value of γ for metal surfaces and ion impact energies similar to those present in our experiments [22] yields good results and is, thus, confirmed by our work. The gas temperature was found to increase as a function of pressure from 300 K at 1 Pa to about 400 K at 100 Pa. In case of using other input parameters in the simulation, significant deviations from experimental results were found. For instance, based on a gas temperature of 350 K and an electron reflection probability of 20%, which represent standard choices of these parameters, the plasma density obtained from the simulation was found to be one order of magnitude lower compared to the experiment at a low pressure of 1 Pa.
The simulation results were found to be sensitive to the gas temperature and the electron reflection probability, but insensitive to the ion induced SEEC (γ) at low pressure. At high pressure, the simulation results are sensitive to the ion induced SEEC, but not to the electron reflection probability and the gas temperature. The sensitivity of the computational results to γ at high pressures is caused by a collisional multiplication of secondary electrons inside the sheaths, which affects the plasma density and other plasma parameters. At low pressure, beams of energetic electrons generated by sheath expansion heating at each electrode propagate collisionlessly through the bulk and arrive at the opposite electrode during the local sheath collapse. Thus, they hit the opposite electrode and their confinement is determined by the choice of the electron reflection probability in the simulation. If they are confined well, the ionization and, thus, the plasma density are enhanced. Via gas depletion the gas temperature affects the mean free path of electrons and, therefore, the confinement of energetic electrons at low pressures. Correspondingly a realistic gas temperature must be used in the simulation to obtain realistic results. For a known gas temperature, the unique sensitivity of the simulation results to r eff at low pressure and to γ at high pressure provides an opportunity to determine these surface coefficients via computationally assisted plasma diagnostics. For different surface materials, distinct plasma parameters could be measured as a function of external control parameters and the experimental results could be compared to results of simulations, where the surface coefficients are changed until agreement with the experiment is found.
This work demonstrated that PIC/MCC simulations can yield realistic results that are in quantitative agreement with experiments for a variety of plasma parameters and over a wide range of discharge conditions. This however, necessitates using realistic gas temperatures and surface coefficients in the simulation. Otherwise significant deviations occur and unrealistic computational results can be obtained.
Clearly, additional experimental validations of plasma simulations should be performed for other gases, plasma sources, and other types of simulation codes.
Our experimental investigations can also serve as the basis of checks of surface coefficients obtained from theoretical models of the plasma-solid interface.