Self-consistent calculation of the optical emission spectrum of an argon capacitively coupled plasma based on the coupling of particle simulation with a collisional-radiative model

We report the development of a computational framework for the calculation of the optical emission spectrum of a low-pressure argon capacitively coupled plasma (CCP), which is based on the coupling of a particle-in-cell/Monte Carlo collision simulation code with a diffusion-reaction-radiation code for Ar I excited levels. In this framework, the particle simulation provides the rates of the direct and stepwise electron-impact excitation and electron-impact de-excitation for 30 excited levels, as well as the rates of electron-impact direct and stepwise ionization. These rates are used in the solutions of the diffusion equations of the excited species in the second code, along with the radiative rates for a high number of Ar-I


Introduction
Capacitively coupled plasmas (CCPs) are highly important in plasma processing technologies, see, e.g.[1][2][3][4][5] and exhibit also a manifold of physics effects [6] making them attractive subjects of experimental and modeling / simulation studies.From the computational side, CCPs can be described by fluid calculations [7,8], particle based simulations [9,10], and the combination of these, termed as hybrid methods [11,12].Among these approaches, particle simulations and hybrid methods comprising a particle simulation component are applicable in the kinetic regime and eliminate the need for any assumptions concerning the electron energy distribution function (EEDF).Particle simulations are mostly accomplished using the particle-in-cell method complemented with Monte Carlo treatment of collision processes (known as 'PIC/MCC'), see, e.g.[9,13,14].
PIC/MCC simulations have tremendously aided the exploration of the physics of low-pressure plasma sources.This approach has frequently been used in studies of charged particle dynamics and particle energy distribution functions in various gases and their mixtures, over a wide range of conditions, like electrode configurations, pressure, driving voltage waveform and amplitude, as well as electronegativity [15][16][17][18][19][20][21].Most of such simulations of electropositive (e.g.noble gas) discharges include only two charged species: electrons and singly-charged atomic ions of the parent gas [22].At higher pressures and/or in molecular gases this simple approach cannot be kept, in such settings a number of ionic species need to be accounted for in the calculations and an increasing number of elementary processes need to be considered [23][24][25][26][27][28][29][30].
The effects of the presence of excited atoms / molecules in an appreciable number in the plasma on the discharge characteristics has already been recognized [31], and addressed in a number of studies, see e.g.[32][33][34].In argon discharges, e.g. both in dc [35] and rf [32,33] operation the stepwise ionization, e − + Ar * → Ar + + 2e − and the pooling ionization, Ar * + Ar * → Ar + Ar + + e − processes were identified to contribute significantly to the overall ionization rate at pressures in the range of 100 Pa.At such conditions, some of the Ar * excited levels, in particular the metastable (1s 5 and 1s 3 ) and the resonant (1s 4 and 1s 2 ) levels, acquire significant density.This makes, e.g.stepwise ionization efficient and has a strong effect on the EEDF.
Recently, a number of works have investigated the importance and the effects of the excited levels on the characteristics of argon CCPs [33,37,38].In these works, a significant influence of the presence of these excited levels on the plasma density was reported as well as an emerging dominance of the stepwise and pooling ionization processes when the pressure approaches the ≈100 Pa range.
The basis of the extensions of models with excited levels are comprehensive sets of stepwise cross sections [39][40][41][42][43][44].These data allow accounting for the re-distribution of the populations of the atoms in the excited levels due to electron-impact collisions (stepwise excitation and de-excitation).Beside these collisional mechanisms, optical transitions between certain levels influence the populations of the excited levels, as well.Models that include both these effects are called collisionalradiative models (CRM) [45][46][47][48][49].They are used primarily for diagnostic purposes by relating the experimental spectra to the plasma parameters such as the electron density and temperature.The models account for the specificity of the investigated discharge.For example, in [45], the calculations were presented for a DC discharge operated at 1000 V voltage, 1 Torr Ar pressure and 2 mA current.Measurements of the electron density and electron temperature based on a CRM were reported for a DC discharge as well in [46].Schulze et al [47] presented a method for the measurement of metastable and resonant atom densities from emission spectra in argon and argon-diluted low pressure inductively coupled plasmas (ICP).The application of a CRM for the determination of the same characteristics in non-equilibrium plasmas was addressed in [48].In [49], the reliability of electron temperature and density diagnostics in an Ar CCP and ICP by optical emission spectroscopy and a collisional radiative modeling was investigated via a comparison with a Langmuir probe diagnostic.
CRMs usually require as input data the electron density and the EEDF for the calculations of the densities of the excited levels and the intensities of the spectral lines that connect the levels considered [36,48,50,51].It has been shown that from the optical emission spectrum one can determine the features of the EEDF [52] and obtain the densities of excited atoms in the metastable and resonant levels [53,54].Measurements of the intensity ratios of certain pairs of spectral lines may facilitate determination of plasma parameters, like the electron density and the electron temperature [55,56], while avoiding the need for absolute calibration of the optical system.
In [31], the possibility of connecting a PIC/MCC simulation with a CRM was examined.In this work, the EEDF and the electron density calculated from the PIC/MCC simulation have been used as input data of the CRM.The PIC/MCC simulation did not include calculation of the density of the excited levels and the CRM did not provide a feedback to the PIC/MCC.This 'uni-directional', i.e.PIC/MCC → CRM coupling was shown to yield reasonably correct spectral line intensities at pressures below ∼ 20 Pa, in the domain where the excited level densities are low to moderate.However, this initial work hinted also at deficiencies of this 'uni-directional' approach when the population density of the excited levels increases.To yield improved results a PIC/MCC ↔ CRM coupling should be realised, which is the subject of our present work.In [57], we have already reported the extension of the PIC/MCC simulation with a set of 1s and 2p excited levels and developed a diffusion-reaction-radiation (DRR) model for the atoms in the excited levels to compute their spatial density distributions taking into account the rates of all electron-impact (direct and stepwise) processes, various ionization channels, quenching reactions, etc.The aim of that work was to derive the spatial density distribution of the Ar 1s 5 metastable atoms and benchmark this with measurements conducted with a laser absorption method.
The present work uses essentially the same models but focuses on the calculation of the optical emission spectrum of the plasma.For that some adjustments have been made to the calculations, such as the introduction of a second escape factor when computing the final spectra in order to provide more realistic emission intensities.The results are then compared to experimental data.To our knowledge, this work is the first to explore the populations of the excited levels and the corresponding spectra in a CCPs while accounting both for kinetic effects and the mutual influence of the excited levels on the electron population.
The paper is structured as follows: In section 2 we describe the elements of the computational framework (the extended PIC/MCC code in section 2.1 and the DRR model in section 2.2).Section 3 gives a brief description of the experimental methods, which were already presented in our previous work [31] in some details.The results of the calculations are presented and compared with experimental results in section 4, while a summary of our work is given in section 5.

Computational methods
The computational framework, whose the scheme is shown in figure 1, consists of two main parts.The first one is the PIC/MCC module that implements the simulation of the kinetics of the charged species (electrons and Ar + ions) in the plasma, while the second one is the DRR module that accounts for the balance and the diffusion of the Ar atoms in the excited levels and for the calculation of the rates of radiative transitions [57].These two modules are executed iteratively.The PIC/MCC module delivers spatially resolved information to the DRR module about the rates of the collision processes (electron-impact processes (that include direct and stepwise excitation and ionization, as well as de-excitation), from which in the DRR module the net source functions for the Ar atoms in the excited levels are computed and used for the solution of the diffusion equations for the excited levels.The DRR module, in turn, transfers to the PIC/MCC module the spatially resolved information about the excited levels.In this study, we take into account the lowest 30 excited levels of the Ar atoms [58].The gas-phase elementary processes considered in the discharge model are summarized in table 1.In the following, the two modules are discussed in some details.

The PIC/MCC module
The PIC/MCC module implements a 'classical' 1d3v (onedimensional in real space and three-dimensional in velocity space) electrostatic simulation of the kinetics of electrons and Ar + ions in a discharge chamber equipped with two planeparallel electrodes.One of the electrodes is driven by a harmonic RF voltage with V pp = 300 V and f = 13.56MHz.Unlike most simulations of Ar plasmas, besides Ar atoms in the ground state (GS), Ar atoms in 30 excited levels are also considered as targets in electron-neutral collisions.The spatial distribution of the densities of the Ar atoms in the distinct excited levels is provided by the calculations in the DRR module.We use a comprehensive set of excitation cross sections that includes excitation from the GS to 30 excited levels and 435 stepwise excitation channels between the excited levels, taken from the BSR database [59] of LxCat [67,68] (see also table 1).The inverse of these processes, i.e. 435 electronimpact de-excitation reactions to lower-lying excited levels and 30 de-excitation processes to the GS are also considered from the excited levels.The cross sections of the de-excitation processes are obtained from the corresponding stepwise excitation cross sections using the principle of detailed balance, see, e.g.[42].Ionization is accounted for via electron-impact from the GS, as well as via stepwise electron-impact ionization from the 1s and 2p levels (14 levels altogether), and pooling ionization (between atoms in any of the 1s levels).The cross sections / rates of the stepwise and pooling ionization processes are taken from [60,61], respectively.For the Ar + ions, we consider elastic collisions only, adopting an isotropic channel and a backward scattering channel in the Ar + +Ar collisions as recommended by Phelps [62].
The GS Ar atoms are supposed to be uniformly distributed within the electrode gap, with a density n gas corresponding to the prescribed gas pressure, p, and temperature, T g , the latter being determined experimentally [31].The depletion of the GS Ar atom density due to the creation of excited level populations is marginal and is, therefore, neglected.
In the PIC/MCC simulation, a uniform spatial grid with N g = 600 − 800 points is used for the calculation of the charged particle densities and the electric potential.The RF cycle is split into N t time steps for the integration of the equations of motion of the electrons (via the leapfrog scheme) with a time step of ∆t e = T RF /N t .N t ranges between 4000 (at the lowest pressure of 2 Pa) and 12 000 (at the highest pressure of 100 Pa).For the ions, subcycling is used with a time step of ∆t i = 20∆t e .These settings respect the stability and accuracy criteria of the PIC/MCC scheme [22,69,70].The surface model in the PIC/MCC code includes ion-induced electron emission characterized by an effective secondary electron yield of γ = 0.07 and electron reflection, characterized by an effective value of r = 0.7 [71,72].
During the execution of a given number of RF cycles in the PIC/MCC simulation, spatially resolved data for the rates of   [49,66] the electron-impact processes (including direct and stepwise excitation and ionization, as well as de-excitation) are calculated by computing the collision frequencies of the respective processes along the trajectories of individual electrons traced in the simulation.The data are interpolated to a grid consisting of N f = 60 points, which provides high enough spatial resolution.Following proper normalization, these (time-averaged) rates, which are the primary output of the PIC/MCC simulation, are written to a data file, which is subsequently read by the DRR module.

The DRR module
The core of the DRR module is a set of spatially onedimensional time-dependent differential equations for the diffusion of the Ar atoms in the 30 excited levels (see, e.g.[33]): where n k is the density and D k is the diffusion coefficient of species k, and S k (x) is the source function that includes the rates of all creation and loss channels for the given species.The diffusion coefficient of the excited atoms is taken to be D • n gas = 1.7 × 10 20 m −1 s −1 [66], with a temperature dependence given in [49].We note that appreciable diffusion takes place for the Ar atoms in the 1s levels only.Atoms in the higher excited levels have very short life times due to radiative channels.
The set of the above equations is solved with an explicit, finite difference forward-time-centered-space method [73], as in [33], with a boundary condition where γ k is the recombination coefficient and v k is the mean velocity of species k at the surface.For γ k a value of 0.5 is adopted [33].As mentioned in section 2.1, the contributions of the electron-impact processes to the temporally averaged (over several RF periods) source functions appearing in equation ( 2) are provided by the PIC/MCC module.The diffusion equations are solved for the steady state.Beside the rates of electron-impact processes accounted for in the PIC/MCC module, the source functions S k (x) also include the rates of quenching processes of the excited levels by neutrals and the radiative transition rates.For the metastable 1s 5 and 1s 3 levels both 2-body and 3-body quenching reactions are included based on [63].For the other excited levels 2-body reactions by neutrals are taken into account [64,65].Although the discharge model is one-dimensional, in order to approximate better the azimuthally symmetric geometry of the experimental system, the radial diffusion losses for the excited species are taken into account as well, via a lifetime defined by their diffusion coefficients and the chamber radius.
The set of 143 radiative transitions is adopted from the work of Siepa [36] along with the corresponding Einstein coefficients.The transitions originating from the resonant 1s 4 and 1s 2 levels (and from some higher levels directly connected radiatively to the GS) are heavily absorbed at the conditions of this study.Therefore for a proper calculation of the populations of these levels the inverse of the radiation processes, i.e. the re-absorption of the radiation has to be considered as well.This is accomplished via adopting an escape factor, which allows a simplified treatment of the radiation trapping.
The escape factor for a given spectral line is usually given in the literature as a functional form, with an argument being the optical depth τ , defined as the product of a characteristic dimension of the system and the absorption coefficient at the center of the given line emitted via spontaneous transition between levels j and i.The absorption coefficient (see, e.g.[49]) is given as where n i is the population of (upper) level i, g i , g j are the statistical weights of the respective levels, A j →i the Einstein coefficient, M the mass of an Ar atom, and k B the Boltzmann constant.
In our model, based on the geometry of the plane-parallel electrode arrangement, we adopt an escape factor that is appropriate for a slab geometry, based on [74,75].(More precisely, we use equations (6a) and (6b) of [75], which originate from [74], but in [75] typographical errors of the corresponding equations of [74] have been corrected.)Actually, a number of the Ar 2p→1s transitions is weakly/moderately trapped as well and, therefore, the escape factor is computed for each transition and is taken into account in the density balance of each excited level.The inclusion of a set of excited levels and the collisional and radiative transitions between these 'automatically' allows the computation of the intensity of the spectral lines of which the upper and lower levels are within the set of excited levels considered.The goal of this work to provide such synthetic spectra of an Ar CCP in the wide pressure range of 2 Pa−100 Pa and their comparison with experimental data measured and published earlier [31].As the optical emission intensity of the spectral lines have been observed from a direction that is perpendicular to the discharge axis, a separate escape factor is used for the computation of the intensity of the radiation that escapes the plasma in the radial direction.We use escape factors proposed by Mewe [76], Apruzese [77], and by Bhatia and Kastner [78].Their values agree well with each other in the asymptotic limits (i.e. in the cases of very strong and very weak trapping), however, depending on the assumptions adopted they differ in the intermediate regime of moderate trapping.As it will be shown later, the calculated optical emission spectrum indeed depends to some extent on the choice of the specific form for the escape factor.
In the experiments reported in [31] the intensities of only a subset of the spectral lines, considered in the current model were measured.A list of these lines, along with some important parameters of these transitions is provided in table 2.

Coupling of the codes
As already mentioned at the beginning of section 2, the PIC/MCC and the DRR codes are executed iteratively to obtain a converged solution for a given set of discharge conditions.This convergence is usually achieved after a few thousands of RF cycles, similarly to 'standard' PIC/MCC simulations.In the iterative solution of our model, the PIC/MCC code is executed for 2-5 RF cycles, which is followed by executing the DRR module, and then this cycle is repeated hundreds to thousands of times.The rates of elementary processes are updated after each run of the PIC/MCC code using a moving average.
We note that the metastable atom density distributions obtained from this computational framework were earlier benchmarked with experimental data obtained with spatially resolved TDLAS measurements [57].This comparison has shown an agreement within a factor of two in the density of 1s 5 metastables within a wide pressure range.Here, we focus on the kinetics of the other excited species and the spectral intensities closely related to the density of these species.

Experimental
The computational results presented in this work are compared to experimental data pertaining (i) the optical emission intensities of selected spectral lines belonging to the set of Ar 2p→1s transitions and (ii) the density of the Ar 1s 5 metastable atoms at the midplane of the discharge.Moreover, the computations use as input data gas temperature values determined experimentally.Since the details of the experimental setup and methods, as well as the data evaluation procedures have been described in [31] to full details, here only a brief reminder about these is given.
The experiments have been conducted using a geometrically highly symmetric plane-parallel electrode CCP source (Budapest v.3 Cell) [79].This source is equipped with two Table 2. Ar 2p→1s radiative transitions whose intensities were recorded in the experiments.The closely separated transitions marked with * are not resolved in the experiments; for these the sum of the intensities is measured and computed.Note that this is a subset of the transitions considered in the DRR module, for a full list of the 143 transitions see [36].A is the Einstein coefficient for the given transition and g i , g j are the statistical weights of the excited levels.disk electrodes, which have 14.2 cm diameter and are made of 304 grade stainless steel.The electrodes are enclosed within a quartz cylinder, their distance is set to L = 4 cm.The plasma was created in Ar with a f = 13.56MHz Tokyo HY-Power RF-150 generator and using a Tokyo HY-Power MB-300 matching network.The RF voltage was measured in the vicinity of the electrode by a Solayl Vigilant Probe.The peakto-peak value of the RF voltage was fixed at V pp =300 V and the gas pressure was varied between p = 2 Pa and 100 Pa.The spectral line intensity measurements were based on a Carl Zeiss Jena PGS-2 spectrometer equipped with an APHALAS CCD-S3600-D-UV detector.Light from the central, ≈1 cm-diameter region of the plasma was captured in these measurements with an optical fiber oriented perpendicularly to the principal axis of the discharge.This spectrometer allowed the separation of the closely spaced pairs of lines at 750.4 nm and 751.5 nm, as well as 800.6 nm and 801.5 nm.The lines at 772.38 nm and 772.42 nm, remained, however, unresolved.The measured intensity values for these latter lines represents the sum of their intensities.Relative sensitivity calibration of the system was carried out using an RS-15 Total Flux Calibration Light Source having a certified calibration report (that specifies the radiant flux of the lamp as a function of the wavelength in the range between 300 nm and 1100 nm) provided by Gamma Scientific.This calibration procedure was based on the measurement of the intensity of the radiation emitted by this lamp (at the conditions 12 V, 8.333 A), with the same optical components (including the quartz cylinder, fiber collimator lens, and fiber optic cable) as in the plasma emission measurements.
The Ar 1s 5 metastable atom density and the gas temperature were determined using Tunable Diode Laser Absorption Spectroscopy [80][81][82][83].The approach relies on the measurement of the absorption on a selected spectral line over which the laser wavelength is scanned through.In [31], the transition Ar(1s 5 → 2p 6 ) at a wavelength of 772.376 nm was chosen for this purpose.In the experimental setup we used a laser diode of type Toptica LD-0773-0075-DFB-1 driven by a control unit Toptica DLC DFB PRO L. The laser light passed through the plasma along its diameter, at the middle of the electrode gap and was detected at the other side of the chamber by a photodiode.To perform proper background subtraction, detector signals were recorded with and without discharge, both with laser on and off states.Assuming dominating Doppler broadening, the amplitude of the absorption provides information about the line-integrated metastable atom density, while the width of the line conveys information about the gas temperature.For a more detailed description, the reader is referred to [31,84].

Results
The discharge conditions used in the computational studies correspond to those in the experimental investigations presented in [31]; the CCP is driven by an RF voltage with a peak-topeak amplitude of V pp = 300 V, a frequency of f = 13.56MHz, and the Ar pressure ranges from p = Pa to 100 Pa.The electrode gap is L = 4 cm.For each set of discharge conditions, measured gas temperature values [31] are used in the computations.The temperature grows from K at 2 Pa to 348 K at 100 Pa gas pressure.
The plasma is sustained by various ionization mechanisms: electron-impact ionization of the GS Ar atoms, i.e. 'direct ionization', electron-impact stepwise ionization (e − + Ar * → Ar + + 2e − ) and pooling ionization (Ar * + Ar * → Ar + Ar + + e − ).The rates of these processes are shown in figure 2, as a function of Ar pressure, temporally averaged  over an RF cycle.The spatially averaged (see panel (a)) direct ionization rate shows no significant dependence on the pressure.However, at the center of the discharge, this rate decays significantly with increasing p as the electrons that gain energy near the edges of the expanding sheaths deposit their energy over a shorter distance due to the higher neutral density and the shrinking energy relaxation length.The spatially averaged rates of the stepwise ionization and pooling ionization processes are nearly equal and grow remarkably (about two orders of magnitude) as the pressure increases.In the center of the discharge (see panel (b)), the stepwise process becomes dominant at p = 30 Pa and even the rate of pooling ionization exceeds that of the ionization of GS atoms at pressures above 40 Pa.This behavior is due to the significant population of the excited levels and the decrease in the number of energetic electrons able to ionize the atoms in the GS.
The computed densities of the Ar atoms in the 1s excited levels in the center of the discharge, are shown in figure 3, as a function of the gas pressure.The densities of the metastable levels 1s 5 and 1s 3 exhibit a peak at 10-20 Pa Ar pressure, while the densities of the 1s 4 and 1s 2 resonant levels increase monotonically with increasing gas pressure.Among these levels, 1s 5 has the highest density.For this level, the corresponding experimental data, measured at identical conditions in [31], are also included in figure 3 for comparison.The computed density matches the experimental data fairly well, considering the complexity of the discharge model, the uncertainties of its input data, the experimental errors and the fact that the measurements capture the line averaged metastable density along the diameter of the plasma, whereas the radial variation of the density is not accounted for in the spatially one-dimensional discharge model [31,57].The increase at pressures below 30 Pa is mostly due to the increase in the electron and neutral gas density, whereas at higher pressures the depletion of the EEPF (discussed later, figure 4) leads to stronger decrease in their production rate constant and hence of the 1s 5 density.In turn, the variation of the 1s 5 metastable atom density with the pressure explains the trends in the pooling ionization rate, observed in figure 2(b): as the 1s 5 level exhibits the highest density, the pooling process is most efficient between atoms in this level, and the rate is defined nearly by the square of the 1s 5 density.Note, that the stepwise excitation rate from this level does not decrease as the electron density grows with increasing pressure, compensating for the decrease of the density of Ar atoms in the 1s 5 level.
At pressures above ≈50 Pa the populations of the resonant 1s 4 and 1s 2 levels exceed that of the metastable 1s 3 level.The monotonic increase of the 1s 4 and 1s 2 density can be explained by the enhancement of the radiation trapping with increasing pressure: the loss rate of the resonant levels due to radiation is −nAη, where the escape factor, η decreases from 7.3 × 10 −4 at 2 Pa to 2.3 × 10 −5 at 100 Pa for the 1s 2 level and from 2.7 × 10 −3 to 8.9 × 10 −5 for the 1s 4 level.At the same time, i.e. with increasing pressure, the loss of the metastables increases due to stepwise excitation from these levels to 2p levels, as will be analyzed below in more details.We note, that a similar relation between the metastable and one of the resonant level densities was found in an earlier work on a direct current hollow cathode discharge [85].
Figure 3 also shows the electron density in the center of the discharge.Over the pressure range covered, n e grows about an order of magnitude, from ≈ 2 × 10 9 cm −3 to ≈ 2 × 10 10 cm −3 .At the lowest pressures, the density values match well those obtained from PIC/MCC simulations, which did not include the excited level populations [31].However, at the highest pressure of 100 Pa, the present simulations predict a factor of two higher density as compared to the previous, more simple simulations, in agreement with the findings of [57].As CRMs require the electron density and the EEDF or the electron energy probability function (EEPF) for the calculation of the spectral line intensities, the accurate computation of these plasma characteristics is essential for reliable spectral properties.As mentioned above, the simulations that include excited level populations provide more accurate values for the electron density.This is also the case for the EEDF/EEPF, the comparison of the different simulation approaches in [57] indicated significant changes of the EEPF at pressures above 20 Pa.The EEPFs obtained from the present calculations in the central part of the plasma are shown in figure 4. The distributions exhibits the well-known 'two-temperature' shape at low pressures characteristic for the α-mode of operation [86][87][88].With increasing pressure, the high-energy tail of the distribution gets gradually depleted due to the more frequent collisions.Due to the stepwise excitation/ ionization processes this depletion is enhanced as compared to cases when such processes are not considered, see, e.g.[57].
The populating and de-populating processes of the 1s levels are shown in figure 5 in the center of the discharge.The data in panel (a), corresponding to the 1s 5 level, indicate that at low pressures radiation from higher levels (data marked as 'Radiative source' in figure 5) and electron-impact (EI) excitation from the GS (data marked as 'EI direct') are the major sources.These are balanced by radial diffusion and stepwise excitation to 2p levels.At the highest pressures, radiative transitions from higher levels and electron-impact stepwise 1s 5 → 2p transitions remain the two processes with the highest rates.The radial diffusion loses its importance, while stepwise excitation to other closely separated 1s levels becomes appreciable.Note, that the two processes with the highest rates (i.e.stepwise excitation to the 2p levels and radiation from these levels to the 1s 5 level) basically cancel each other out.This essentially means that the energy the electrons lose in inelastic collisions with the atoms in the 1s levels is converted into light.Further, the two processes, i.e.EI excitation from the 1s levels and radiative transitions from the 2p levels closely link the two groups of levels but does not change their net population.The density balance of the 1s 5 level is more influenced by lower-rate processes.Specifically, the 1s levels are populated by direct excitation from the GS and are lost by the resonant radiation from the 1s 4 and 1s 2 levels.The loss of the metastable levels involves an intermediate step of collisional transfer to one of these resonant levels.Then, overall, the production and loss of the metastable atoms involves electron collisions, leading to relative independence of their density on the electron density.This allows the decrease in their density seen in figure 3 despite the steady increase of the electron density with the pressure.In other words, the effect is caused by changes in the collisional rates that reflect the changes in the EEDF, especially in the region of higher electron energies instead of changes in the densities of the participating species.
Between 10 Pa and 20 Pa, the axial diffusion changes from a loss process to a source process of the 1s 5 atoms.At low pressures, their source is concentrated within the central domain Rates involved in the balance of the 1s excited levels, including electron-impact (EI) rates from ground state, to/from other 1s and 2p levels, the rate of radiative source and loss, the loss rate due to pooling ionization, as well as axial and radial radial diffusion source/losses.The data correspond to the grid points nearest to the midplane of the discharge.Discharge conditions: Vpp = 300 V, f = 13.5 MHz, L = 4 cm. of the plasma and the metastables diffuse 'outwards', while at higher pressures their creation by electron-impact is confined to the vicinity of the positions of the maximum sheath width and therefore 'inwards' diffusion towards the discharge center takes place.(As the logarithmic scale of figure 5 does not allow displaying a function that changes sign, the positive and negative contributions are shown separately as 'source' and 'loss' functions.)The 1s 3 level behaves in a quite similar way, but the rates of important processes are about an order of magnitude lower as compared to those for the 1s 5 level.
For the 1s 4 level, electron-impact excitation from GS and radiation from higher levels are the two dominant populating processes over the whole pressure range, with a contribution of stepwise excitation and de-excitation from other 1s levels.The dominant loss is due to radiation, with some contribution of stepwise processes to other 1s and 2p levels at the highest pressures.The behavior of the 1s 2 level is also similar to the other resonant level, 1s 4 , with comparable rates of the important processes.Stepwise ionization and pooling ionization, and electron-impact de-excitation from 2p levels are less important as loss/source channels for the 1s levels.The diffusion is only important for the metastable levels at low pressures and is unimportant for the resonant levels due to their shorter effective lifetime as compared to those of the metastable levels.
The 2p levels develop considerably lower populations as compared to the 1s levels, as shown in figure 6.As a result, these excited levels have a negligible contribution to the ionization balance, compared to the contribution from the levels in the 1s block.The densities of the 2p levels are in the range of 10 5 cm −3 -10 7 cm −3 , the values increase monotonically with increasing gas pressure for most of the levels.An exception among the higher density levels is the 2p 9 level, whose density is maximum at p ≈ 20-30 Pa, which is followed by a slight decrease.The explanation of this behavior is that this level is dominantly excited by a stepwise process from the 1s 5 level and, therefore, its density closely follows the pressure dependence of the density of the 1s 5 level (see later).Another exception among the lower density levels is the 2p 1 level, which exhibits a density maximum at p ≈ 10 Pa.This level is known to be dominantly populated by electron-impact excitation of GS Ar atoms.The decay of the density of atoms in this excited level at the center of the discharge at higher pressures can be explained by the decrease of the density of high-energy electrons in the discharge center (most of the excitation occurs near the sheath boundary that moves towards the electrode as the pressure is increased).
Next, we analyse the populating and depopulating processes of these two excited levels, Ar 2p 1 and 2p 9 .The rates of these processes in the center of the discharge are shown in figure 7, as a function of the pressure.For the 2p 1 level (panel (a)), at low pressure the electron-impact excitation of GS atoms is the dominant populating process, which is balanced by the radiation from this level.Cascade radiation ('Radiative source') from higher levels contributes in the order of 10% to the populating processes.The kinetics of this level is then close to the corona limit.The importance of stepwise excitation from 1s levels is negligible at low pressures.However, this process becomes dominant even for this level at the highest pressure of 100 Pa.Electron-impact excitation from other 2p levels and from higher levels as well as electron-impact de-excitation to GS and to other 2p levels and 1s levels is insignificant at all conditions studied.For this level and for any other 2p levels the axial and radial diffusion play a minor role even at the lowest pressures (and thus it is not shown in this figure).Quenching by neutrals is found unimportant for the 2p 1 level, too.The 2p 9 level exhibits a different behavior as compared to the 2p 1 level.At the lowest pressure of 2 Pa, electron-impact excitation from the GS, stepwise excitation from the 1s levels, and radiation from higher (>2p) levels populate this level with approximately the same rate, while the dominant loss is represented by radiation from this level.Regarding the populating processes: (i) electron-impact excitation from the GS to the 2p 9 level is efficient since at low pressure the EEDF contains a lot of high-energy electrons (figure 4), (ii) for the same reason also higher excited levels (>2p) are efficiently populated by direct electron-impact excitation of ground-state Ar atoms and due to the specific nature of the system of radiative transitions, these processes also provide a substantial populating source for the 2p 9 level, and (iii) the 2p 9 level is known to have a high cross section for stepwise excitation from the highly populated 1s 5 level.Towards the higher pressures, the electron-impact excitation from the GS loses its importance, also the radiative source becomes less important.The major populating process remains the stepwise electronimpact excitation from the 1s levels.All of these observations can be linked to the depletion of energetic electrons in the EEDF with the increase of the pressure (figure 4).Neutral quenching of this level is more important, as compared to the 2p 1 level.The observations for the other, lower-rate processes in the case of the 2p 1 level remain valid for the 2p 9 level, too.
While it could be argued that the changes of the excited level populations as a function of pressure are ultimately due to the changes of the EEDF, it should also be kept in mind that the EEDF is also influenced by the populations of the excited levels as the rates of the stepwise excitation and de-excitation processes depend on these populations and the change of the electron energy after each collision leads to a modification of the EEDF.The present computational approach ensures that all these processes are accounted for in a self-consistent manner.Consequently, the shape of the EEDF and the rates of the elementary processes with various threshold energies (which 'define' the densities of the Ar atoms in the excited levels) are highly correlated.At low pressures, e.g. when many energetic electrons are present, high-threshold-energy processes, like electron-impact excitation of ground-state Ar atoms are important and play an essential role in building up the populations of the 1s and 2p levels (as shown in figures 5 and 7, respectively).With increasing pressure, the high-energy tail of the EEDF gets depleted and we observe that while the source of the 1s levels by direct electron-impact excitation remains substantial, the source of 2p levels starts to decay above 10 Pa.These different behaviors are caused by the differences in the excitation thresholds (which are 11.55 eV for the lowest-lying 1s level and 13.07 eV for the lowest-lying 2p level) as well as by the fact that for all levels above the 1s block an alternative source with lower threshold, namely the 1s levels, exists.The increasing number of electrons in the intermediate (∼ few eV) energy range at high pressures promotes the stepwise excitation processes, which results in Rates involved in the balance of the 2p 1 (a) and 2p 9 (b) excited levels, including electron-impact (EI) rates from/to the ground state, 1s, other 2p, and higher levels, the rates of radiative source and loss, the quenching rate with neutrals, and the loss rate due to stepwise ionization.The data correspond to the grid points nearest to the midplane of the discharge.Discharge conditions: Vpp = 300 V, f = 13.5 MHz, L = 4 cm.substantial increase of the 2p level populations from the 1s levels as shown in figure 7.
The main goal of this work, as mentioned earlier, is the computation of the optical emission spectrum of the plasma.As discussed above, the model yields the populations of Ar atoms in the excited levels (n j ) in the discharge.From this information, the relative intensities of the spectral lines can be calculated in a straightforward way, as these can be associated with the radiative losses of the specific 2p levels towards the 1s levels, i.e.
using the Einstein coefficients (A j →i ) given in table 2, which lists the spectral lines of interest here.Here we used the notation η * for the escape factor, to emphasize that, as the optical emission intensity of the spectral lines has been observed from a direction that is perpendicular to the discharge axis, a separate escape factor (appropriate for a cylindrical geometry) is used for the computation.This escape factor is different from the one used in the model for the plane-parallel electrode, 'slab' geometry.This is a notable change in comparison to the model in [57].
Figure 8 shows the computed spectral line intensities in comparison with the experimental results.The data are presented in terms of normalized (dimensionless) intensities: as in the experiments only relative intensities were measured, normalization was carried out both for the experimental and computed intensities, to ensure that the sum of the line intensities gives 1, i.e. and In the calculation of the intensities, Mewe's escape factor was used in the DRR module.There is an overall good agreement between the experimental and calculated data.We find a good agreement at 2 Pa and 4 Pa, larger differences (up to about 30%) show up at pressures of 7-10 Pa.At higher pressures the good agreement is restored.Selected data are also re-plotted in figure 9 using a logarithmic scale for the relative intensity.While this representation suppresses the small deviations between the experimental and computed data, it allows checking the level of the agreement between these data sets over a wider dynamic range, i.e. for the low-intensity lines, as well.This figure confirms that the computational reproduction of the spectrum is quite good over three orders of magnitude of the spectral line intensities.We note, that the spectral line intensities actually depend to some extent on the choice of the 'radial' escape factor used.This is illustrated in figure 10, where the intensities calculated with escape factors proposed by Mewe [76], Apruzese [77], and Bhatia and Kastner [78] are displayed.These data were also normalized as defined above.One can observe significant differences between the spectral line intensities obtained with these various escape factors.This finding actually points out an inherent weakness of adopting an escape factor for the computations.The reason of this sensitivity is that most of the lines considered are partially trapped at the conditions of the present study, even at the lowest pressure.
It is notable that for most lines, the escape factor of Bhatia and Kastner results in higher intensities while the formula of Mewe gives lower values than the one provided by Apruzese.However, for the strongly absorbed lines (with escape factor from Mewe of less than 0.4), the trend is exactly the opposite.This is consistent with the dependence of the escape factors on the optical depth, illustrated in figure 11.For low optical depths τ < 1, the escape factor of Bhatia and Kastner is larger, i.e. it predicts less self-absorption than the estimations by the other authors, while the expression of Mewe gives the smallest escape factor and, consequently, the largest self-absorption.At large optical depths (τ > 1), the situation is reversed.Naturally, all three estimations have the same asymptotic values at negligible (η → 1 for τ ≪ 1) and at infinite optical depths (η → ∞ for τ ≫ 1).Another interesting property of the escape factor expressions is that for τ = 1 they all give a value close to 1/2.The escape factor of Bhatia and Kastner then gives a quick transition between negligible self-absorption to fully self-absorbed lines within a narrow range of variation of the optical depth.In contrast, the expression of Mewe provides a more gradual change.As figure 11 shows, this means that even at an optical depth of τ = 10 −2 , the relation of Mewe predicts some self-absorption.This is illustrated in figure 12, where the escape factor values for each of the spectral lines are shown as calculated from Mewe's formula [76].Only the radiation on the lines for which η * ≈ 1 can freely escape the plasma, while lower values indicate gradually more important trapping.At the lowest pressure of 2 Pa the density of the excited levels is the lowest and the self-absorption should be low.Still, the plasma is optically thin only for about half of the lines, as indicated by the respective escape factors η * ≈ 1.For two of the lines: 2p 6 → 1s 5 (λ = 763.5 nm) and 2p 9 → 1s 5 (λ = 811.5 nm), however, we find η * ≪ 1 already at this low pressure.Note that these lines share the 1s 5 lower level and the decreased escape factor is indicative of a significant population of this metastable level even at 2 Pa.Due to this strong absorption by the 1s 5 level, these lines are a good choice for laser absorption measurements of the density of this metastable level even when it is not strongly populated.Note, however, that a significant 1s 5 population does not result 'automatically' in radiation trapping.For example, for some other lines, like the 2p 4 → 1s 5 (λ = 714.7 nm) line, because of the relatively low Einstein coefficient or, equivalently, oscillator strength of this transition (A = 6.25 × 10 5 s −1 , respectively f = 2.87 × 10 −3 ) compared to those of the 2p 6,9 → 1s 5 transitions (for which A = 2.45 × 10 7 s −1 and A = 3.3 × 10 7 s −1 , respectively, f = 2.14 × 10 −1 and f = 4.6 × 10 −1 ), see table 2. The λ = 714.7 nm transition is not trapped even at 100 Pa pressure despite the significant population of the 1s 5 lower level.For a comparison, the escape factors for the resonant lines originating from the 1s 4 and    from the GS under wide range of plasma conditions, as illustrated also by figure 7(a).This is unlike any of the other 2p levels, where significant, often dominant, contribution to the population comes from excitation from the levels in the 1s block.Since the population in these levels usually depends on the electron density, this means that ratios involving the 750.4 nm line would be more sensitive to the electron density.Indeed, the ratio I 696.5 /I 750.4 has been recommended by [55] as being sensitive to the electron density, based on the CRM developed in that work.In [56] it was shown experimentally that the ratio I 763.5 /I 750.4 is even more sensitive to this parameter.
The other two line ratios involving the 750.4 nm line, i.e.I 811.5 /I 750.4 and I 763.5 /I 750.4 , discussed in [49] match well the ratios found in the experiment, similarly to the I 696.5 /I 750.4 intensity ratio.The ratios I 763.5 /I 750.4 , and I 696.5 /I 750.4 keep the good agreement throughout the entire pressure range investigated.The third ratio involving the 811.5 nm line (2p 9 → 1s 5 ) Figure 13.Ratios of intensities of selected pairs of spectral lines as a function of the pressure, as obtained from the experiment and simulations.The calculation is based on Mewe's escape factor [76]. Discharge conditions: Vpp = 300 V, f = 13.5 MHz, L = 4 cm.
shows deviations from the experimental values likely due to the peculiarities of the 2p 9 level.For this level the 811.5 nm line is the only radiative transition.Hence, this line is usually one of the most intensive ones (figure 8).For this reason it usually experiences also a significant self-absorption (figure 12), thus making it particularly sensitive to the radiation trapping model used (figure 10).Due to the uncertainties about the reliability of the different escape factor models outlined above, the use of this line ratio is not recommended.
The last line ratio investigated, I 763.5 /I 738.4 , shows notable differences at low pressures.For this pair of lines, the agreement improves with increasing pressure and reaches very good agreement above 30 Pa.While the 763.5 nm line is always notably self-absorbed, the transition at 738.4 nm experiences some self-absorption only at higher pressures (figure 12).The latter fact is related to the lower level of the transition being the resonant 1s 4 level, whose density attains appreciable values only at the higher pressures investigated (figure 3).This means that the 763.5 nm line is always affected by the escape factor model while the predicted intensity of the 738.4 nm line is influenced by it only at higher pressures.When both intensities become self-absorbed, the uncertainty due to the escape factor model probably largely compensate each other and the line ratio becomes less influenced by this aspect, thus improving the agreement with the experimental observations.The differences observed in the line intensity ratios, and between the measured and computed spectra, in general, point out the limitations of the model, which are briefly discussed below.We can identify the following points which may decrease the accuracy of the modeling calculations: (a) While the experimental system has a two-dimensional (cylindrical) symmetry, the discharge model is onedimensional, i.e. it cannot take into account the radial variations of the densities of the species.The somewhat higher computed 1s 5 metastable atom density (observed at most pressure values) with respect to the measured data can well be attributed to this: as the density of the metastable atoms can be supposed to decay towards the radially limiting surface of the cell, the simulations overestimate the metastable density.Due to the line averaged nature of the measurement, the experimental value is likely somewhat underestimated.(b) Treating the radiation trapping with an escape factor is itself a simplification.Moreover, we found that using escape factors from various literature sources changes the computed spectral line intensities to some extent.(c) The model considers a high number (30) of excited levels.
However, the computations find the excitation rate to even higher levels (termed as 'Rydberg excitation' in the cross section set adopted) significant.This is illustrated in figure 14, where the rates of the most important processes that determine the balance of the populations of the 2p 1 and 2p 9 levels are reproduced from figure 7, along with the excitation rate of the 'Rydberg' levels.In the present model, the (presumably) radiative decay of these higher lying levels, which is supposed to populate the 2p levels, is not accounted for.Figure 14(a) shows that at the lowest pressures the Rydberg excitation rate is about an order of magnitude higher than the direct electron-impact excitation rates of the 2p 1 and 2p 9 levels.The importance of radiative transitions originating from high-lying levels was already pointed out in several previous studies, see, e.g.[45,54,[89][90][91][92].To estimate the uncertainty in neglecting the contribution of the 'Rydberg' levels, the excitation rate in figure 14 can be assumed to be distributed equally among the ten 2p levels.The result will be an additional source, that is comparable with the sources accounted for by the present model.This leads to an uncertainty in the intensities of 50% at the lower pressures.At the higher pressures, studied here, this source of error becomes negligible.

Summary
In this work, we reported the development of a computational framework for the calculation of (a part of) the optical emission spectrum of a low-pressure argon CCP, focusing on a Rates of the most important processes involved in the balance of the 2p 1 (a) and 2p 9 (b) levels, as in figure 7, but showing only the processes of major importance.Additional data for the excitation rate of higher-lying levels included in the electron-impact cross section set as a single entry (referred to as 'Rydberg' excitation [59]), are also shown.
subset of 2p→1s transitions.The calculations have been based on a particle-in-cell/Monte Carlo collisions (PIC/MCC) simulation code and a DRR code, which were executed iteratively.A basic PIC/MCC code was extended to include 30 excited levels of the Ar atoms.Electron collisions with these excited atoms could result in stepwise excitation and ionization, and in de-excitation processes.The rates of these processes were computed in the PIC/MCC simulation and were transferred to the DRR code, which solved the diffusion equations of the four 1s and ten 2p excited levels, also taking into account pooling ionization, quenching reactions, as well as radiative transitions between the 30 excited levels and between some of these and the GS of Ar atoms.The computations were carried out for a RF driving voltage with a peak-to-peak value of V pp = 300 V, a frequency of f = 13.56MHz, Ar pressures in the range of p = 2 Pa to 100 Pa and an electrode gap of L = 4 cm, matching the conditions of the experimental investigations reported in [31].
The iterative solution of the PIC/MCC and DRR modules provided the rates of all reactions and the densities of Ar atoms in fourteen (1s and 2p) excited levels.The correctness of the calculations was validated earlier in [57] in a comparison of the computed spatial density distributions of Ar 1s 5 metastable atoms at various discharge conditions with corresponding experimental data obtained via laser absorption spectroscopy.The present computational framework also yields the intensities of the spectral lines.Unlike a 'stand-alone' collisional radiative model, the present approach includes the self-consistent calculation of the electron density and the EEDF and is, thus devoid of any assumptions regarding the shape of this function.The present approach represents significant improvement compared to our earlier work [31].In [31], the EEDF and the electron density obtained in a simple PIC/MCC simulation (that did not include the populations of the excited levels) was coupled into a CRM, which yielded the spectral line intensities.In that work, the approach was identified to be applicable only at pressures not exceeding about 20 Pa, where the modification of the EEDF by the excited atoms is insignificant.The 'bi-directional' coupling (i.e. the iterative execution of the PIC/MCC and DRR parts of the calculations (where the DRR contains the CRM)) adopted here, can fully account for these changes of the EEDF and therefore it provides much better spectral line intensity results, especially at pressures above 10 Pa.
We have found that beside the resonant transitions, which are very heavily trapped, most of the 2p→1s lines are also partially trapped.Therefore, an escape factor was adopted for all the transitions upon the calculations of the radiative transition rates.For the self-consistent calculations of the excited atom densities an escape factor appropriate for a slab geometry (that approximates best the experimental apparatus) was used.For the calculation of the spectral line intensities observed perpendicularly to the discharge axis, as in the experiments, another form of the escape factor, corresponding to a cylindrical geometry was adopted.We have found, that using the functional form of the escape factor proposed in [76] a fairly good agreement can be obtained between the measured spectra [31] and the computed ones over the whole pressure range of interest (2-100 Pa).Non-negligible differences were, however, found between the spectra computed with different escape factors available in the literature.This indicates the need for a dedicated study of the accuracy and the range of validity of the various escape factor models available in the literature.
We have also analyzed the limitations of the present approach, which were found to be (i) the one-dimensional nature of the simulation which cannot capture the twodimensional density distributions that exist in the experimental system, (ii) the sensitivity of the synthesized spectra on the choice of the escape factor, and (iii) omission of the cascade populating processes from the high-lying excited levels (termed as 'Rydberg levels' in the present cross section set).Solutions to these issues would require (i) developing a twodimensional version of the simulation package, (ii) using a more sophisticated approach for the escape factor to account for radial effects, or the implementation of a Monte Carlo treatment of the radiation as in [93,94], and (iii) inclusion of a more detailed cross section set for the higher-lying excited Ar-I levels to account for the radiative channels (cascades) that populate the excited levels that are already considered in the present model.Overcoming these limitations is considered as future work.

Figure 1 .
Figure 1.Scheme of the computational framework.

Figure 2 .
Figure 2. Spatially averaged (a) and central (b) ionization rates as a function of the gas pressure.The data are time-averaged, the discharge conditions are: Vpp = 300 V, f = 13.5 MHz, L = 4 cm.

Figure 3 .
Figure 3.Computed populations of argon atoms in the 1s excited levels and electrons in the midplane of the discharge, as a function of the pressure.The curve labelled as '1s 5 Exp'.shows the experimental data for the 1s 5 level density, from TDLAS measurements presented in[31].Discharge conditions: Vpp = 300 V, f = 13.5 MHz, L = 4 cm.

Figure 4 .
Figure 4. Time-averaged EEPF in the central part of the discharge (in the region 0.45 ⩽ x/L ⩽ 0.55), as a function of the gas pressure.Discharge conditions: Vpp = 300 V, f = 13.5 MHz, L = 4 cm.

Figure 5 .
Figure 5.Rates involved in the balance of the 1s excited levels, including electron-impact (EI) rates from ground state, to/from other 1s and 2p levels, the rate of radiative source and loss, the loss rate due to pooling ionization, as well as axial and radial radial diffusion source/losses.The data correspond to the grid points nearest to the midplane of the discharge.Discharge conditions: Vpp = 300 V, f = 13.5 MHz, L = 4 cm.

Figure 6 .
Figure 6.Computed populations of the 2p excited levels at the midplane of the discharge, as a function of pressure.Discharge conditions: Vpp = 300 V, f = 13.5 MHz, L = 4 cm.

Figure 7 .
Figure 7.Rates involved in the balance of the 2p 1 (a) and 2p 9 (b) excited levels, including electron-impact (EI) rates from/to the ground state, 1s, other 2p, and higher levels, the rates of radiative source and loss, the quenching rate with neutrals, and the loss rate due to stepwise ionization.The data correspond to the grid points nearest to the midplane of the discharge.Discharge conditions: Vpp = 300 V, f = 13.5 MHz, L = 4 cm.

Figure 9 .
Figure 9.Comparison of the normalized spectral line intensities: measured values vs. computed values (based on Mewe's escape factor), using a logarithmic scale to illustrate the behavior over a wider dynamic range.Discharge conditions: Vpp = 300 V, f = 13.5 MHz, L = 4 cm.

Figure 10 .
Figure 10.Comparison of the computed normalized spectral line intensities obtained with different escape factors.Discharge conditions: Vpp = 300 V, f = 13.5 MHz, L = 4 cm.

Figure 13
Figure 13  examines the ratios of the intensities of selected pairs of spectral lines, often discussed.Three of these pairs contain the 2p 1 → 1s 2 line at 750.4 nm.The 2p 1 level is known to be predominantly populated by direct excitation

Figure 11 .
Figure 11.Escape factors for radiative losses in the radial direction in cylindrical geometry as a function of the optical depth τ , given by different authors.

Figure 12 .
Figure 12.Escape factors for the spectral lines of interest, for selected values of the gas pressure.The values are calculated according to the formula by Mewe [76].Discharge conditions: Vpp = 300 V, f = 13.5 MHz, L = 4 cm.

Figure 14 .
Figure 14.Rates of the most important processes involved in the balance of the 2p 1 (a) and 2p 9 (b) levels, as in figure7, but showing only the processes of major importance.Additional data for the excitation rate of higher-lying levels included in the electron-impact cross section set as a single entry (referred to as 'Rydberg' excitation[59]), are also shown.

Table 1 .
Gas-phase elementary processes considered in the model.Here, Ar * and Ar * * denote excited levels, with ε(Ar * ) < ε(Ar * * ).Ar r and Ar m represent the lowest resonant (1s 4 and 1s 2 ) and metastable (1s 5 and 1s 3 ) levels, respectively.GS stands for the ground state.For the processes involving radiation, '×2' indicates that both spontaneous emission and re-absorption processes are considered.