{'fraction_non_alphanumeric': 0.031644987820011954, 'fraction_numerical': 0.016132738888633542, 'mean_word_length': 4.732446318008168, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Advanced fluid modelling and PIC/MCC simulations of low-pressure ccrf discharges

Comparative studies of capacitively coupled radio-frequency discharges in helium and argon at pressures between 10 and 80 Pa are presented applying two different fluid modelling approaches as well as two independently developed particle-in-cell/Monte Carlo collision (PIC/MCC) codes. The focus is on the analysis of the range of applicability of a recently proposed fluid model including an improved drift-diffusion approximation for the electron component as well as its comparison with fluid modelling results using the classical drift-diffusion approximation and benchmark results obtained by PIC/MCC simulations. Main features of this time- and space-dependent fluid model are given. It is found that the novel approach shows generally quite good agreement with the macroscopic properties derived by the kinetic simulations and is largely able to characterize qualitatively and quantitatively the discharge behaviour even at conditions when the classical fluid modelling approach fails. Furthermore, the excellent agreement between the two PIC/MCC simulation codes using the velocity Verlet method for the integration of the equations of motion verifies their accuracy and applicability.


Introduction
Capacitively coupled radio-frequency (ccrf) discharge plasmas are widely used in plasma processing technologies. Typical examples for the application of ccrf discharges are plasma-enhanced chemical vapor deposition and plasma etching [1]. Besides experimental diagnostics, numerical modelling and simulation of ccrf discharges provide established tools to get detailed insights into the discharge physics [2][3][4].
Particle-in-cell/Monte Carlo collision (PIC/MCC) simulation is the most recognised method for the theoretical description of ccrf discharges [3,[5][6][7][8][9][10][11]. In the PIC/MCC method, a collection of particles is followed in space and time taking into account particle-particle and particle-wall interactions as well as the effect of the self-consistently determined space charge field [2,12]. Fluid (or continuum) models, which are based on a hydrodynamic description of the plasma, provide an alternative approach for the analysis of ccrf discharges [4,[13][14][15][16][17]. Compared to PIC/MCC simulations, the numerical solution of fluid models is computationally less demanding and hence more attractive for practical applications. A disadvantage of the fluid approach is its limited application range when low-pressure ccrf discharges are considered. Here, the fluid description is applicable as long as the charged particles' mean free path is much smaller than the characteristic dimension of the discharge [18]. Different modelling approaches have been developed that aim at a combination of the advantages of kinetic simulations and fluid models in so-called hybrid methods [19][20][21][22][23] or at an improvement of the accuracy of the classical fluid description [18,[24][25][26][27][28][29][30][31]. Since the latter is in the focus of the present paper, an overview of recent approaches is given here.
In non-thermal plasmas, the energy is mainly delivered through the electrons. Therefore, an accurate description of the electron component is crucial for the reliability of an integral plasma model [32] and hence researchers have spent large effort to improve the description of the electron component. Earlier works aimed at a customized description of fast electrons by means of so-called beam models, see, e.g., [18,33,34] and references therein. Robson et al. [24] have introduced a physically based fluid model for electrons in low-temperature plasmas based on the so-called heat flux ansatz. This approach has been applied in Ref. [35] to describe periodic electron structures in a constant electric field by means of a fluid model and it has been generalised in Ref. [25] to electron swarms under the influence of nonuniform electric fields. In another work of Robson et al. [26], the accuracy of fluid models for light particles has been improved by a direct substitution technique that uses swarm transport data instead of cross sections for the calculation of the collision terms. Rafatov et al. [27] have proposed a fluid model that includes a nonlocal ionization source term in order to overcome fundamental shortcomings of the classical fluid description. A non-local collisionless electron heat flux has been considered in Ref. [28] to get an improved description of the collisionless electron heating effect in low-pressure, high-frequency ccrf discharges. Furthermore, fluid models for electrons in non-thermal plasmas based on a four-moment description have been developed [29,30]. Note that in the frame of the hydrodynamic plasma description two moment equations (particle density and flux) are commonly considered for heavy particles, while three moment equations (particle density, flux and mean energy density) are usually taken into account for the electron component [4,36]. In Refs. [37,38] it has been shown that the four-moment model proposed in Ref. [29] is more accurate than conventional fluid approaches for the theoretical description of negative streamer fronts in nitrogen and neon, while the classical fluid model using the local field approximation (LFA) [39,40] gives reasonably good results under these conditions, too [38]. In contrast, Grubert et al. [41] have shown that the LFA is not applicable for the investigation of low-pressure gas discharge plasmas and the local mean energy approximation (LMEA) [15] should be used instead. A novel LMEA based drift-diffusion approximation for electrons has been derived in Ref. [31] from the four-moment model proposed in Ref. [30]. This model has been found to be more accurate than the classical LMEA based fluid description for dc glow discharge plasmas at low and atmospheric pressure [31,42]. The main advantage of this drift-diffusion approach established in Ref. [31] compared to most of the other modified fluid descriptions is that it is not limited to specific discharge conditions and it does not increase the computational cost [42].
Comparisons between fluid and particle methods for the simulation of low-pressure ccrf discharges have been carried out before, e.g., in Refs. [8,[43][44][45][46]. From the rigorous comparison of different PIC/MCC simulation and fluid modelling results for helium ccrf discharges in the pressure range from 4 to 40 Pa in Surendra [46] it can be concluded that the degree of agreement between fluid and PIC/MCC simulation methods strongly depends on the gas pressure and is also divergent for different plasma properties. For the lowest pressure considered in [46], larger differences between PIC/MCC and fluid results of about 50 to 60 % have been found for the predicted rf voltage and plasma density while generally good agreement for other important properties, such as the ion flux to the electrodes (error of 20 to 30 %), has been reported. For the same benchmark situation as in Ref. [46], increasing differences in the plasma density obtained by the classical LMEA based fluid model and by PIC/MCC simulations with rising gas pressure (20 % at 4 Pa to 35 % at 40 Pa) have been found in [8].
In Ref. [43] moments of the Boltzmann equation have been calculated using PIC/MCC simulations. It has been shown that the convective term in the momentum balance equation for ions is of particular importance if the collisionality of the sheaths is low, while it can be neglected for electrons if secondary electron emission is of minor importance. Furthermore, it has been pointed out that the spatial and temporal variation of the electron energy plays a predominant role in low-pressure ccrf discharges.
Different fluid models for low-pressure ccrf discharges have previously been compared, e.g., by Young et al. [47] and Chen et al. [28]. It has been shown that the reliability of fluid models for ccrf discharges can be enhanced by an adequate description of electron momentum and energy transport [28,47] and, particularly, by the introduction of an integral form for the electron heat flux which provides an improved prediction of electron heating mechanisms [28].
In Refs. [30,31,42,48,49] it has been pointed out that an improved prediction of electron heating and energy transport mechanisms can also be achieved by the recently developed drift-diffusion model [31]. In the present paper, this fluid modelling approach is described and applied to the analysis of low-pressure ccrf discharges. The fluid modelling results are compared to results obtained by PIC/MCC simulations in order to evaluate the applicability range of the improved fluid modelling approach and to demonstrate its advantages against the classical fluid description for low-pressure ccrf discharges. Following the strategy of Turner et al. [8] and Surendra [46], well defined benchmark conditions for ccrf discharges in helium and argon are considered for this purpose.

Methods
The evaluation of the applicability range and accuracy of two different fluid modeling approaches for the analysis of low-pressure ccrf discharges is performed considering simple, spatially one-dimensional discharge situations in helium and argon. The configuration and main discharge features are sketched in figure 1. The plasma between Particle density Electric potential Figure 1. Sketch of the spatially one-dimensional discharge configuration with time t, voltage amplitude V 0 , frequency f , electrode separation d and maximum plasma density n max . Dashed grey lines: electric potential at two different times t; dashed blue lines: electron density at two different times t; solid red line: ion density.
plane electrodes separated by the distance d is driven by a sinusoidal voltage with amplitude V 0 and frequency f . The specific discharge parameters used for four different cases in helium and three different cases in argon are listed in table 1. For the conditions considered in the present work, ions cannot follow the oscillating electric field while the electron density strongly changes with time t in the sheath regions.
Advanced fluid modelling and PIC/MCC simulations of low-pressure ccrf discharges 5 The maximum plasma density n max typically occurs at the centre of the discharge region. As in Ref. [8], it is assumed that the plasma is composed only of electrons and positive ions in the background gas helium or argon. Collision processes are limited to interactions between these charge carriers and the neutral background gas.
The fluid description of the electron component is performed by means of two different drift-diffusion approaches: the novel drift-diffusion model introduced in Ref. [31] and abbreviated by DDAn (cf. section 2.1) and the commonly used classical drift-diffusion model, named DDA53, using simplified electron energy transport coefficients [50] (cf. section 2.2). Even though elaborations of both approaches have been published elsewhere, specific details of the present implementations are given in sections 2.1 and 2.2. In addition, a time-dependent two-moment model for ions is applied. It takes into account the continuity equation and the momentum balance equation and is described in see section 2.3. The balance equations for electrons and ions are complemented by the Poisson equation for the electric potential φ. Here, e 0 , ε 0 , n i and n e denote the elementary charge, the vacuum permittivity and the ion and electron densities, respectively. Details about the boundary conditions are given in section 2.4 and the numerical solution of the coupled set of fluid model equations is described in section 2.5. In addition to the fluid modelling approaches, two PIC/MCC simulation codes were developed independently and are applied for mutual verification and benchmarking of the different fluid models for the parameter range considered. Details of the PIC/MCC simulation procedures are given in sections 2.6.

Novel drift-diffusion approximation for electrons
The particle density n e and the energy density w e of electrons with mass m e are determined by the solution of the coupled balance equations ∂ ∂t n e (z, t) A consistent drift-diffusion approximation for the particle flux Γ e and the energy flux Q e of the electrons has been deduced by an expansion of the electron velocity distribution function (EVDF) in Legendre polynomials and the derivation of the first four moment equations from the electron Boltzmann equation [30,31]. It reads − e 0 m eνe and includes the momentum and energy flux dissipation frequencies ν e andν e as well as the transport coefficients ξ 0 , ξ 2 ,ξ 0 andξ 2 . These properties are given as integrals of the isotropic part f 0 and the first two contributions f 1 and f 2 to the anisotropy of the EVDF over the kinetic energy U of the electrons, respectively, according to Here, λ e is the electron-energy-dependent mean free path of the electrons [41]. Their particle density n e , energy density w e , particle flux Γ e and energy flux Q e are given by The coefficients (6)- (11) are determined within the framework of the common LMEA as functions of the mean electron energy U e = w e /n e by solving the stationary, spatially homogeneous Boltzmann equation in multiterm approximation for different values of the electric field [41,51], see section 3.1.
The source term S in equation (2) represents the gain of electrons due to ionization of neutral gas atoms in electron-neutral collisions while P in equation (3) describes the loss of electron energy in elastic, exciting and ionizing collisions of electrons with neutral gas atoms. They are given by S(z, t) = k io n gas n e (z, t) , Here, n gas = p/(k B T gas ) is the density of the background gas with pressure p, temperature T gas and mass M , k B denotes the Boltzmann constant, k ex j and k io are the respective rate coefficients for excitation and ionization processes with energy thresholds U ex j and U io , respectively. The energy rate coefficient for energy dissipation in elastic collisions is denoted byk el . The rate coefficients are given by where Q ex j , Q io and Q m are the electron-neutral collision cross sections for excitation, ionization and momentum transfer in elastic electron-neutral collisions, respectively. As for the transport coefficients, the rate coefficients are incorporated into the fluid model as functions of the mean electron energy in the framework of the local mean energy approximation.

Classical drift-diffusion approximation for electrons
The conventional drift-diffusion model for electrons in nonthermal plasmas comprises the balance equations (2) and (3) with the expressions for the particle and energy fluxes [41,51]. The particle and energy diffusion coefficients D e andD e as well as the mobilities b e andb e are given by integrals of the EVDF over energy space [41,51]. Because numerical problems arise in many situations when the consistent expressions for the energy transport coefficientsD e andb e are used [31,52], the simplified expressionsD e = 5/3 D e andb e = 5/3 b e are usually applied [27,51,[53][54][55]. This classical drift-diffusion approximation with the transport coefficients is used here for comparative studies. Note that the simplified energy transport coefficients are valid in case of a Maxwellian EVDF, only.

Two-moment model for ions
For the description of low-pressure ccrf discharges the ion inertia must be taken into account by considering an effective electric field [17] or by solving the time-dependent momentum balance equation [23]. The latter approach is chosen here and the system of moment equations is solved to determine the density n i and the particle flux Γ i of ions with mean velocity v i = Γ i /n i and mass m i . The ion pressure p i is given by the ideal gas law p i = n i k B T i where it is assumed that the heating of ions is negligible, i.e., T i = T gas . Although the assumption T i ≈ T gas is generally not valid in the sheath regions of ccrf discharges, it does not have any significant impact since the pressure term in equation (26) is generally unimportant [43]. The ion momentum dissipation frequency ν i is obtained from measured ion mobilities b i according to ν i = e 0 /(m i b i ), see section 3.1.

Boundary conditions and initial values
The Poisson equation (1) is supplemented with the boundary conditions φ(0, t) = V 0 sin(2πf t) and φ(d, t) = 0, according to the setup shown in figure 1. Note that the desired amplitude J 0 of the discharge current density is used as input for the fluid modelling and the amplitude V 0 of the rf voltage applied to the powered electrode is automatically adapted in each period according to until a periodic state is reached. Here, J calc is the actual amplitude of the discharge current density. This procedure ensures that J calc = J 0 when periodic state is reached. In order to exclude any uncertainties regarding the implementation of boundary effects in the different modelling approaches, the boundary conditions for the particles are set as simple as possible. It is assumed that neither reflection of particles nor emission of secondary electrons occur at the electrodes located at z = 0 and z = d. The expression applied at the boundaries for electrons (j = e) and ions (j = i) reads where ν = −1 at z = 0 and ν = 1 at z = d. The drift velocity v d j is deduced from the respective expressions in equations (4), (21) and (26) and v th j = 8k B T j /(πm j ) denotes the thermal velocity of species j. For electrons, the "temperature" T e = 2 U e /(3 k B ) is used here. Similarly, the boundary condition for the electron energy balance equation (3) is given by ξ 0 E in case of the novel drift-diffusion approximation (5) and v d e = −5/3 b e E for the classical drift-diffusion approximation. The choice of boundary conditions for the hyperbolic system (25)-(26) requires special care because the number of required boundary conditions at each boundary depends on the direction of the characteristics as well as the ion sound speed c i = k B T i /m i as described, e.g., in Refs. [56,57]. Table 2 lists the number of boundary conditions that is applied here. The condition ∂n i /∂z = 0 is used in addition to equation (29) if two boundary conditions are to be set. A logarithmic extrapolation of ion properties is performed to complete the number of physical boundary conditions where required. Table 2. Number of boundary conditions to be set for the hyperbolic system (25)-(26) depending on the mean ion velocity v i and the ion sound speed c i .
In addition to boundary conditions, initial values are to be set for the particle, momentum and energy balance equations (2), (3), (25) and (26). Here, a quasi-neutral state with a homogeneous initial density n j (z, 0) = n init of charged particles and mean electron energy U e (z, 0) = U init e is assumed. The mean ion velocity is set to zero at t = 0. Since we do not allow for the emission of secondary electrons from the electrodes, the initial charge-carrier density n init must be large enough to sustain the plasma. However, it was verified that the periodic discharge behaviour is not influenced by the specific values used.

Numerical solution of fluid models
For the numerical solution of the system of nonlinearly coupled equations a finitedifferences discretization in space is performed and a semi-implicit time-stepping scheme is applied. The number of spatial grid points and the size of the time step used for the numerical solution are specified below. Both parameters are chosen such that converged results are obtained.
For Poisson's equation (1) the standard central difference quotient of second order [58] is used while the parabolic drift-diffusion equations (2) and (3) for the particle and energy density of electrons with fluxes (4) and (5) or (21) and (22), respectively, are discretized by means of the exponentially fitted Scharfetter-Gummel finite-difference scheme [59] as described in Ref. [60]. The hyperbolic system of moment equations (25) and (26) for the density and the flux of ions is discretized in space by the standard first-order upwind scheme [58] and a predictor-corrector approach is chosen for the time-coupling of the equations as detailed in Ref. [42]. Here, the number of boundary conditions affects the solution method by evaluating the ion sound speed and the mean ion velocity for each time step and setting the required physical and numerical boundary conditions according to table 2. The correctness of the computer code has previously been verified by the comparison with other methods and computer codes for a wide range of test problems and discharge conditions [32,42,61].

PIC/MCC simulation procedure
Two independent PIC/MCC simulation codes for low-pressure ccrf plasmas were developed at the Institute of Theoretical Physics and Astrophysics (ITAP), University of Kiel, Germany (named PIC(ITAP)) and at the Leibniz Institute for Plasma Science and Technology (INP) Greifswald, Germany (named PIC(INP)) for mutual verification and for validation of the fluid models. Details of the PIC(INP) method, which has been extended from a previous model for streamer simulations [62], are given in Ref. [11]. Although the general procedure of both PIC/MCC codes is the same a brief description is given here.
Compared to the fluid models discussed above, PIC/MCC simulations are particlebased, i.e., they track the trajectories of so-called superparticles under the influence of the self-consistent electric field determined by solving equation (1). The two PIC/MCC simulation codes used here resolve one space dimension and trace all three velocity components, usually referred to as 1d3v. Particles are represented in the cloud-incell scheme, where the weight factor for the charge distribution on the grid decreases linearly from the particle position toward the grid points. The weight of superparticles is constant in all simulations, i.e., the adaptive particle management available in the PIC(INP) code is disabled to avoid possible uncertainties. Both PIC/MCC codes are parallelized using the Message Passing Interface (MPI).
For the time integration of the equations of motion the velocity Verlet algorithm is used as opposed to the leapfrog method usually applied in many other codes [8]. It has been shown in Ref. [11] that the velocity Verlet method converges faster than the leapfrog scheme with respect to the size of the time step used in PIC/MCC simulations of low-pressure ccrf discharges.
Particles that reach the electrode surface at z = 0 or z = d are fully absorbed and are removed from the simulation. This is consistent with the boundary conditions applied for the fluid models (see section 2.4). Except for backscattering in elastic ionneutral collisions, the particles are scattered isotropically in the collision events and the remaining energy in ionizing collisions is shared in equal parts between the two electrons. Instead of calculating the collision probability for all particles individually, the null-collision method is used. Further details on the PIC/MCC algorithm can be found, e.g., in Refs. [11,12,63]. Table 3 summarizes the applied fluid and PIC/MCC simulation methods used for the comparative investigations. 3. Input data

Physical data
Since the aim of this work is to verify the reliability of two different fluid approaches in comparison with PIC/MCC simulations, the sources of uncertainties are reduced by making the physical discharge model of the considered low-pressure ccrf discharges in helium and argon as simple as possible. Table 4 lists the considered collision processes with their energy thresholds where applicable and the sources from which the corresponding cross section data sets are taken. For helium, the same set of collision processes and cross sections as in Ref. [8] is used which include elastic electron-neutral collisions, excitation of the triplet and singlet helium states, direct ionization of helium in its ground state and elastic ionneutral collisions. The input data used for the present calculations are depicted in figure 2. Note that the ion cross-sections are used for the PIC/MCC simulations only and are not considered for the fluid modelling. There, the collisional impact of ionneutral collisions is taken into account by the momentum dissipation frequency ν i as obtained from measured ion mobilities (cf. section 2.3).
A similar set of processes is considered for argon. It comprises elastic electronneutral and ion-neutral collisions as well as the total electron impact excitation with an energy loss of 11.55 eV and ionization of argon. As, e.g., in Refs. [68,69] the cross section data set from the JILA database of Phelps [65] is used here. More specifically,  the fit formulas given in Refs. [66,67] with a lower limit of 0.1 eV for the ion energy are applied for the generation of the cross section input data shown in figure 3.
The electron-neutral collision cross sections represented in figures 2a and 3a are not only used as input for the PIC/MCC simulations but are additionally utilized for calculation of the momentum and energy flux dissipation frequencies (6) and (7), the transport coefficients (8)- (11), (23) and (24) as well as the rate coefficients (18)- (20). This is done by solving the stationary, spatially homogeneous Boltzmann equation in eight-term approximation for the required range of electric field strengths. Here, the method described in Ref. [70], which has been modified to take into account  nonconservative electron collision processes correctly, has been employed. For utilization in the fluid models, the obtained transport and rate coefficients are tabulated as functions of the mean electron energy which is also obtained by the solution of the electron Boltzmann equation.
The ion momentum dissipation frequency appearing in equation (26) is determined from measured ion mobilities b i depending on the reduced electric field E/n gas according to ν i = e 0 /(m i b i ). For helium, the He + ion mobility is taken from the measurements of Patterson [71] and the measured data of Ellis et al. [72] is used for Ar + .

Numerical parameters and initial values
In order to ensure the general validity of the conclusion to be drawn from the present investigations, the considered pressure range for the ccrf discharges in helium and argon is extended as much as possible. For helium, four different pressures (cases C1-C4) in the range from 10 to 80 Pa are considered (see table 1). At lower pressures the validity of the drift-diffusion approximation used for the fluid description of electrons becomes inadequate. On the other hand, at pressures above 80 Pa the PIC/MCC simulations require long computing times and hence they become less appropriate.
The pressure range used for argon spreads from 20 to 80 Pa (see table 1). Here, no converged PIC/MCC simulation results could be obtained for 10 Pa because of a strong sensitivity of the simulation results on the initial superparticle number N sp . As pointed out by Turner [73,74], this effect is caused by velocity space diffusion and cannot easily be circumvented in certain situations.
The numerical input parameters and initial values used for the fluid modelling and PIC/MCC simulations for the four discharge cases in helium and three cases in argon are given in table 5. Note that special care was placed on choosing the numerical parameters adequately to obtain converged results. Nevertheless, as discussed in Ref. [8], the accuracy of PIC/MCC simulation results remains difficult to assess as Table 5. Numerical input parameters and initial value used for the simulation of ccrf discharges in helium and argon at a pressure of 10 (C1), 20 (C2), 40 (C3) and 80 Pa (C4) and a gas temperature of 300 K. it depends on several numerical parameters and their combination. The development of a generally accepted effective procedure for successively refining the results would be highly desirable. It is also worth mentioning that the initial plasma density n init and energies U init e,i might affect the number of rf periods required to obtain periodic results but do not influence the periodic behaviour investigated in the following section. About 1000 (C1) to 4000 (C4) rf periods were required to reach the periodic state for the conditions given in table 5. Prior to the evaluation of the two different fluid modelling approaches DDAn and DDA53 in section 4.2, the discharge behaviour predicted by the PIC/MCC simulation procedures is discussed in the following section.

Main discharge features
The temporal change of the applied voltage and the discharge current density at the powered electrode as well as the spatiotemporal behaviour of the electron and ion densities at 20 and 80 Pa in helium and argon are depicted in figure 4. Here, normalized  For all considered conditions ions cannot follow the electric field and hence their density is stationary in time. In contrast, electrons respond almost instantaneously to the change of the electric field in the sheath regions. Hence, the electron current at the momentary anode is maximal close the instant of the largest applied voltage as shown in figures 4a, b. These figures also show that the phase shift between current and voltage decreases with increasing pressure both in helium and in argon. Because the collision frequency increases with raising pressure, the width of the sheath regions is smaller and Here, electrons are considered to be "fast" if their energy is larger than e 0 V 0 /4, where V 0 is the voltage amplitude given in table 1. In both gases the maximum density of fast electrons occurs in front of the momentary cathode just before the applied voltage reaches its maximum at t × f = 0.25 (positive applied voltage, cathode at z = d) and at t × f = 0.75 (negative applied voltage, cathode at z = 0), respectively. However, in argon the profile of the fast electron density is much more localized than in helium which indicates that in argon the electrons loose their energy more rapidly when propagating towards the anode. The corresponding time and space averaged isotropic part of the EVDF f 0 (U )/n e normalized by the electron density n e is shown in figure 5. Obviously, it is increasingly influenced by inelastic electron-neutral collisions for larger pressures.
Advanced fluid modelling and PIC/MCC simulations of low-pressure ccrf discharges 17 The electron impact excitation and ionization processes lead to a marked depletion of the electron population above the lowest threshold energy for exciting collisions, i.e. 19.82 eV for helium and 11.55 eV for argon. The agreement of the two independently developed different PIC/MCC simulation codes for all considered discharge conditions illustrated in figure 5 and also found for all macroscopic properties mutually verifies their correctness.

Comparison of fluid modelling and PIC/MCC simulation results
In order to evaluate the accuracy of the novel fluid model DDAn and the classical fluid model DDA53 for the considered discharge conditions, their results are compared to macroscopic quantities derived from the kinetic PIC/MCC simulations. For this purpose, the amplitude of the current density J 0 obtained by the PIC/MCC simulations as indicated in table 6 is used as input and the amplitude V 0 of the rf voltage is automatically adapted during the calculations according to equation (28). This is in accord with the procedure of Turner et al. [8] who also used a fixed discharge current for benchmarking. Figure 6 table 1). It is found that the rf voltage required by the fluid models to reach the same current density J 0 is generally 10 to 30 % larger than that of the PIC/MCC simulations. This holds for both the fluid models. Only the deviation of the rf voltage required by DDA53 to sustain the prescribed current density in argon increases monotonically with increasing pressure to 40 % at 80 Pa. Figure 7 shows the results obtained by the two PIC/MCC simulation codes and the fluid models for the time averaged ion density at 20 Pa in helium (figure 7a) and argon (figure 7b). As for the time and space averaged isotropic distribution f 0 /n e (see figure 5), the predictions of the different PIC/MCC simulation codes PIC(ITAP) and relation to that of the PIC/MCC simulations n PIC i . It is found that the large deviation of the maximum ion density obtained by DDA53 at 20 Pa in argon reduces to 30 % at a pressure of 80 Pa. At the same time, the differences between DDAn and PIC/MCC for argon reduce from 20 to 10 % if the pressure is increased from 20 to 80 Pa. By contrast, the smallest deviation between PIC/MCC and both fluid modelling results for the maximum of the ion density in helium is obtained at the lowest pressure of 10 Pa. This might be explained by the fact that the EVDF in helium at 10 Pa is almost Maxwellian (cf. figure 5a). In such situations fluid approaches are generally more adequate.
The temporal variation of the electron and ion fluxes at the powered electrode (z = 0) are presented in figure 9 for helium and argon at 20 Pa. Due to the symmetry of the discharge configuration, the same behaviour can be observed at the grounded electrode (z = d) with a time shift of t × f = 0.5. Again, the results of both fluid modelling approaches are in qualitative agreement with the PIC/MCC simulations for helium. However, only the novel drift-diffusion approximation DDAn is in conformity with the PIC/MCC simulations for argon while larger deviations are obtained when using the DDA53 fluid model. The differences between the PIC/MCC simulation results and the fluid results for the particle fluxes at the electrodes are quantified in figure 10 by means of the time averaged ion fluxes Γ i at the powered electrode. Note that the accurate knowledge of this parameter is of particular importance for many applications [3]. Figure 10 shows that for helium the fluid models DDA53 and DDAn underestimate the average ion flux by about 10 to 20 % compared to the PIC/MCC simulations. Here, the deviations for DDA53 are slightly smaller than for DDAn. In argon, DDAn underestimates the average ion flux by about 20 to 30 %. By contrast, DDA53 largely overestimates the average ion flux in argon. The differences to the PIC/MCC simulation results increase from 50 % at 20 Pa to 60 % at 80 Pa.
In summary, the plasma density and the average ion flux at the electrodes as "global" plasma parameters are predicted by the novel drift-diffusion approximation DDAn with an uncertainty of less than 30 % compared to PIC/MCC simulations. Much larger errors of up to 60 % are to be expected for these parameters for the classical fluid model DDA53 in argon. For helium, smaller deviations of less than 40 % are observed. Similar differences between classical fluid models and PIC/MCC simulations have previously been reported for ccrf discharges in helium [46] and argon [36] at comparable conditions.
In order to get deeper insights into the differences of the considered fluid descriptions, the spatial variations of macroscopic quantities derived from fluid modelling and PIC/MCC simulations at the instant t × f = 0.25 are compared in the following. At this time the voltage applied at z = 0 reaches its maximum, V 0 , and the momentary cathode is at z = d. Figure 11 shows the spatial variation of the first four moments of the EVDF, namely density, mean energy, particle flux and energy flux of electrons together with the electric field, the ionization rate and the electron heating rate, as obtained by the fluid models DDAn and DDA53 as well as the PIC/MCC simulations. The density of fast electrons with energies higher than e 0 V 0 /4 obtained by PIC/MCC simulations is additionally depicted in figure 11b. In general, both fluid models are able to reproduce most discharge features qualitatively and quantitatively fairly well. Larger differences can particularly be observed in the mean electron energy (figure 11c). Both fluid approaches overestimate the mean electron energy in the plasma bulk region by about 20-30 %. A large discrepancy can also be observed in the mean electron energy obtained by DDAn and DDA53 in the cathode region close to z = d. However, this is of minor importance for the present situation because no emission of secondary electrons is considered and hence the electron density in the sheath region is extremely low (n e < 1 m −3 ). Note that this is also the reason why no PIC/MCC data for the mean electron energy U e = w e /n e is available for that region. The overestimation of the mean electron energy in the bulk plasma by the fluid models is caused by the fact that, particularly at lower pressures, highly energetic electrons significantly contribute to  Figure 11. Results for helium at 20 Pa: Spatial variation of electric field E (a), electron density n e (b), mean electron energy U e (c), electron heating rate −e 0 EΓ e (d), ionization rate S (e), particle flux Γ e (f) and energy flux Q e (g) of electrons obtained by the different modelling approaches at time t × f = 0.25. the ionization budget. At the same time, the density of these fast electrons (figure 11b) is too low to have any effect on the mean electron energy. Since the ionization rate coefficient used in the fluid description depends on the mean electron energy only, (see section 3.1), a higher mean energy is enforced by the need to deliver the predefined current. This lack of a precise description of electrons contributing to the electron generation but not to the mean electron energy also causes a more localized profile of the ionization rate determined by the fluid models compared to the PIC/MCC results (figure 11e).
Another important difference between the fluid modelling results and the PIC/MCC solution concerns the obtained sheath width. Particularly the spatial profiles of the electric field (figure 11a) and the fluxes (figures 11f, g) show that both fluid models overestimate the sheath width by about 5 %. This difference causes the spatial shift between the PIC/MCC and fluid results, e.g., for the maximum of the heating rate in the sheath/plasma transition region (figure 11d). Apart from this, the results of DDAn and DDA53 for the electron heating rate as well as the particle and energy fluxes are in good agreement with the PIC/MCC results. It is worth mentioning that the fluxes Γ e and Q e are highly transient quantities which makes the direct comparison of results obtained by different methods at a certain time difficult.
The corresponding spatial behaviour of the electric field, electron density, mean energy, heating and ionization rate, electron particle and energy fluxes for argon at 20 Pa and the instant t × f = 0.25 is shown in figure 12. The main difference to the  Figure 12. Results for argon at 20 Pa: Spatial variation of electric field E (a), electron density n e (b), mean electron energy U e (c), electron heating rate −e 0 EΓ e (d), ionization rate S (e), particle flux Γ e (f) and energy flux Q e (g) of electrons obtained by the different modelling approaches at time t × f = 0. 25. situation for helium is that here only the novel fluid description for electrons, DDAn, is able to predict most of the macroscopic quantities with the same accuracy as in helium. Most likely, the presence of the Ramsauer minimum in argon causes the occurrence of nonlocal transport effects which can be captured by the DDAn approach but not by the classical model DDA53 due to the coinciding particle and energy transport of the electrons induced by the assumption of a Maxwellian EVDF [31]. The divergences between the results of DDA53 and PIC/MCC particularly affect the electron heating rate (figure 12d) and the energy flux (figure 12g) where large differences in the spatial profiles can be observed. In contrast, the fluid model DDAn including a consistent description of electron energy transport, provides a comparatively good prediction of these quantities when compared with the PIC/MCC simulation results. However, the comparison of the PIC/MCC results for the spatial profile of the fast electron density (figure 12b) with the spatial profile of the ionization rate (figure 12e) indicates that highly energetic electrons predominantly determine the ionization rate, similar to the behaviour found in helium. Hence, the marked differences between the results of DDAn and the PIC/MCC simulations for the spatial profiles of the mean electron energy (figure 12c) and the ionization rate (figure 12e) are again caused by the improper consideration of ionization in the plasma bulk induced by fast electrons.
The approach introduced by Rafatov et al. [27] for low-pressure dc glow discharges aiming at an enhanced description of nonlocal ionization by adding an additional source term was found to be not applicable for the modelling of ccrf discharges. The separate description of highly energetic electrons by the Monte Carlo collision method (see, e.g., [19,75]) could be a more promising extension of the present fluid model description to overcome the remaining shortcomings. Note that the discrepancy between fluid and PIC/MCC methods regarding the spatial distribution of the ionization rate drops with raising pressure due to the increase of the collisionality and the associated decrease of the impact of fast electrons.

Conclusions
In the present work the applicability and the accuracy of two different fluid approaches for the analysis of low-pressure ccrf discharges were investigated by benchmarking them against PIC/MCC simulations. The considered fluid methods comprise time-dependent particle and momentum balance equations for ions as well as a novel drift-diffusion approximation (DDAn) and the classical drift-diffusion approximation with simplified energy transport coefficients (DDA53) for electrons, respectively. In order to assure the general validity of the findings and to provide a test bed for future studies, simple ccrf discharge configurations in helium and argon at pressures ranging from 10 to 80 Pa were considered. Main findings of the comparative studies are the following: • Results of the novel as well as the classical fluid model are in good qualitative and quantitative agreement with macroscopic quantities derived from PIC/MCC simulations for ccrf discharges in helium.
Here, the novel drift-diffusion approximation provides a slightly better prediction of the plasma density than the classical drift-diffusion approximation with deviations of less than 30 % for the DDAn model and 40 % for the DDA53 approach compared to PIC/MCC results.
• For argon, the classical fluid model fails to reproduce the discharge features predicted by PIC/MCC simulations. In contrast, the novel drift-diffusion approximation maintains its applicability and reliability and provides a prediction of relevant plasma parameters with deviations of less than 30 % compared to the PIC/MCC simulation results.
• At lower pressures both fluid models fail to correctly reproduce the spatial profile of the ionization rate. This is caused by the lack of an adequate description of highly energetic electrons which contribute to the ionization budget but not to the mean energy of electrons.
Future studies will address the questions how the influence of highly energetic electrons can be better included in the novel fluid modelling framework presented here and how plasma-boundary interactions can be described correctly.