A non-Maxwellian kinetic approach for charging of dust particles in discharge plasmas

Nanoparticle charging in a capacitively coupled radio frequency discharge in argon is studied using a particle in cell Monte Carlo collisions method. The plasma parameters and dust potential were calculated self-consistently for different unmovable dust profiles. A new method for definition of the dust floating potential is proposed, based on the information about electron and ion energy distribution functions, obtained during the kinetic simulations. This approach provides an accurate balance of the electron and ion currents on the dust particle surface and allows us to precisely calculate the dust floating potential. A comparison of the obtained floating potentials with the results of the traditional orbital motion limit (OML) theory shows that in the presence of the ion resonant charge exchange collisions, even when the OML approximation is valid, its results are correct only in the region of a weak electric field, where the ion drift velocity is much smaller than the thermal one. With increasing ion drift velocity, the absolute value of the calculated dust potential becomes significantly smaller than the theory predicts. This is explained by a non-Maxwellian shape of the ion energy distribution function for the case of fast ion drift.


Introduction
The charging of dust particles in a discharge plasma is crucial in the description of dust-plasma interaction. Usually the floating potential of the dust is calculated from the balance of electron and ion currents on the dust particle surface using a Maxwellian distribution function for electrons and ions (see, for example [1]). However, it is known from experiments that in a capacitively coupled radio frequency (ccrf) discharge without dust the electron energy distribution function (EEDF) can be far from Maxwellian even in a bulk plasma. The EEDFs were measured in [2] in a ccrf discharge in helium and argon to illustrate the difference in the electron kinetics in non-Ramsauer and Ramsaur gases. It was shown that in helium the shape of the EEDF remains Maxwellian for all pressures. In argon the EEDF shape changes from Druyvesteyn to bi-Maxwellian with decreasing gas pressure. The dust particle charge was investigated in [3] in a low pressure argon discharge plasma using a Boltzmann equation for the electrons and fluid equations for the ions. It was found that the presence of the dust particles modifies the shape of the EEDF from Druyvesteyn to Maxwellian due to the fact that the dust grains collect fast electrons.
Unfortunately the modification of the ion energy distribution function (IEDF) with the presence of dust is not completely understood. The shifted Maxwellian distribution is the commonly used assumption for the IEDF in a non-zero external electrical field. Particle in cell Monte Carlo collision (PIC-MCC) simulations of charging of unmovable dust were performed in [4]. However, the shapes of the EEDF and IEDF, as well as the procedure of equilibration of the electron and ion currents on the dust particle surface were not discussed. In [4] the unmovable dust was taken to be uniformly distributed over the discharge gap including the electrode sheaths and this assumption makes the results less reliable, because the heavy negatively charged dust is not able to penetrate in the electrode sheath area. Recently the formation of the steady state nanosize dust profile in the rf discharge and the influence of dust size on the plasma parameters were studied in argon in [5] and in a mixture of argon and acetylene in [6].
In this paper, we study the charging of the nanoparticles in a ccrf argon discharge by using the EEDF and IEDF from PIC-MCC simulations. The discharge parameters with the static dust distribution are found from a self-consistent PIC-MCC simulation. Special attention is paid to compare the PIC-MCC EEDF and IEDF with the Maxwellian ones.
This paper is organized as follows. In section 2, we present the model of a discharge plasma with movable nanoparticles. Section 3 is devoted to the description of our new method of calculation of the dust particle charge. In section 4, we present the results of our PIC-MCC 3 calculations and compare them with the analytic model results. Our conclusions are given in section 5.

Physical model
The PIC-MCC model of a high frequency discharge plasma with nanoparticles is based on the electron and ion kinetic equations, including their interaction with dust and the Poisson equation. In this model, the velocity distribution functions of electrons, f e (t, x, v), and ions, f i (t, x, v) (three-dimensional over velocity and one-dimensional in space) are obtained by solving the Boltzmann equations numerically, where v e , v i , n e , n i , m e and m i are the electron and ion velocities, densities and masses, respectively, E is the electric field, J e , J i are the collisional integrals for electrons and ions which include elastic and inelastic collisions with background gas and negatively charged dust particles.
The electron and ion kinetics in the PIC-MCC model are simulated with test particles with weight W , which means the number of electrons/ions represented by this test particle. Below we will call the test particles electrons and ions. The equations of motion are integrated for electrons and ions in the discharge electric field accounting for the collisions with neutral atoms and with negatively charged dust. The probabilities of various events (elastic scattering, excitation of neutrals, ionization and absorption on the dust surface) during an electron/ion time step are set by the cross sections depending on the test particle's energy and simulated by using the null-collision method [7] and a random number generator. If an ionization event takes place, a new electron and ion with the same coordinate and weight are introduced. Electrons/ions leaving the calculation area through the electrodes or captured by dust are excluded from the simulation. The interelectrode gap is covered with a calculation grid (one-dimensional in this case). The macroscopic characteristics of the discharge such as the plasma densities, ion and electron average energy, the ion/electron flux are calculated from the EEDF and IEDF over the calculation grid. The electrical field distribution is calculated at each time step from the Poisson equation: where n e and n i are obtained from the PIC-MCC calculation. The boundary conditions are the sinusoidal applied voltage on one electrode and zero voltage on the second grounded electrode. In this model, the space charge of the dust particles with the density n d (x) and the mean particle charge Z d (x) is included. The particle charge is found self-consistently during the simulation, as described below. The calculation starts from chosen initial distributions for the electrons, ions and dust corresponding to some analytical approximation for the discharge plasma. The iterative procedure of solving equations (1)-(3) proceeds until a steady state solution is achieved. Usually this takes about 1000 simulated rf cycles for discharge with unmovable dust. The number of particles representing the electrons and ions in our PIC-MCC simulations is usually 10 4 -10 5 .

4
The interelectrode distance is divided into a calculation grid with 101-201 nodes for definition of the macroscopic plasma parameters. The time step is τ e = (0.3 − 1) × 10 −11 s for electrons and τ i = (1 − 3) × 10 −10 s for ions. For electron scattering in argon the cross sections are taken from [8]. For ions we take into account the resonant charge exchange collisions. Additionally our model incorporates the cross sections of electron and ion absorption by negatively charged dust particles, taken from the orbital motion limit (OML) theory [9]: here ϕ d is the floating potential of the dust particle surface gained in the discharge plasma, ε e,i is the kinetic energy of the electron/ion, r d is the radius of the dust particle, e is the elementary charge. Since the calculations were performed for low gas pressure and for nanometre-size particles, the conditions of validity of the OML approach l e,i λ D r d (where l e,i is the mean free paths for electrons and ions, λ D is the Debye screening length) are fulfilled here. Note that our model is not limited by the OML approach, and it is possible to use any other absorption cross sections in analytical or tabular form. In this work, we restrict our study to static dust and concentrate on the determination of the dust floating potential.

New method for calculation of the dust particle charge
In the steady state of discharge glow, the floating potential of the dust particle ϕ d can be derived from the equality of the electron and ion currents onto the dust particle surface. The volume densities of both currents j e and j i onto the dust surface (in other words, the charge carrier absorption rate per unit volume) at the position x and time t are given by the following basic equations (see, for example [10]): where ε is the kinetic energy, V e,i are the magnitudes of electron and ion velocity vectors, respectively. Averaging of equations (5) and (6) over a discharge cycle leads to an integral equation for ϕ d at each position x: In the traditional OML approach, using the cross sections (4) and assuming the Maxwellian velocity distribution for electrons and ions, equations (5) and (6) give the analytical expressions for the currents onto a single dust particle [10]: 5 for electron current, and, taking into account that ions may have drift velocity u of the same order of magnitude as the thermal one, hence they are assumed to have shifted on u Maxwellian distribution: for ion current.
Here n e,i and T e,i are the electron/ion density and temperature, k is the Boltzmann constant, g = u/ √ kT i /m i is the ratio of ion drift and thermal velocities. Note also that expression (8) has two well-known limit cases: for the case g 1, and The above shown expressions are widely used in OML theory. If the values of n e,i (x) and T e,i (x) are known in every position x, the floating potential ϕ d may be calculated from equations (7) and (8) if I e = I i . To a first approximation this allows us to complete the model of rf discharge with nanoparticles, assuming the dust particle charge Z d (x) = r d ϕ d (x), which is correct for the particles with size much smaller than the Debye screening length.
However, since it is known that the energy distribution in a discharge plasma may be far from Maxwellian, the aim of the developed algorithm is to avoid the assumption about the Maxwellian shape of the energy distribution. The PIC-MCC simulation of the complex plasma provides two other possible ways for obtaining of the dust particle charge. The simplest way for the direct definition of the currents is to count the number of events of electron and ion absorption by the dust particles with the cross sections (4) during the test particles motion simulation, with adjustment of the floating potential to minimize the difference in the electron and ion currents. The problem of such an approach is the statistical noise of the currents, that requires averaging of their magnitudes over large numbers of rf cycles. We offer a more accurate and less time consuming method, namely to calculate the currents by the equations (5) and (6), using the EEDF and IEDF f e (t, x, ε) and f i (t, x, ε), obtained from the PIC-MCC simulation of the discharge plasma with dust.
In the PIC-MCC simulation, the integration over energy and the averaging over a discharge cycle in equations (5) and (6) is equivalent to a summation over the test particles ensemble: where k is a summation over electrons, and where the summation is over the ions. The above expressions are the direct analogue of equations (5) and (6) for the case of a discrete ensemble of test particles. Here, W k,l is the 6 weight of the kth electron/lth ion, V k,l and ε k,l are its velocity magnitude and kinetic energy, respectively. The summation in equations (9) and (10) is carried out during one rf cycle and only over the particles corresponding to the calculation grid node with coordinate x (the terms W k,l V k,l n d (x)σ id,ed are taken at each grid node using the standard PIC procedure). For the selfconsistent definition of the floating potential ϕ d , we apply the iterative method, in which after each rf cycle we add a correction ϕ to the value of ϕ d . To accelerate the convergence of the iterations, we take into account that the dependence of j e on ϕ d is exponential, whereas for j i it is linear. The correction is calculated at each iteration step in all nodes of the calculation grid from the following equation: where ϕ d is the value from the previous iteration and T e is the mean electron temperature taken from the PIC-MCC simulation (here it is used as the iteration parameter). Equation (11) is solved by the Newton method. Then the summation in equations (9) and (10) is performed again at the next rf cycle with the corrected value for the potential ϕ d + ϕ and equation (11) is solved again. After 3-4 iterations, the currents j PIC e and j PIC i become practically equal for all x. The proposed method for the calculation of the dust floating potential ϕ d is restricted neither by the OML assumptions (the cross sections σ ed and σ id may be given in another form) nor by the assumption about the Maxwellian distribution for the electrons and ions. Our method provides an accurate dust potential determination during the rf discharge simulation and can be used to test the validity of the traditional OML approach (equations (7) and (8)). It should be mentioned that the additional computational time required for the summation in equations (9) and (10) and solving equation (11) is very small in relation to the computation time needed for the solution of equations (1-3).

Results and discussion
The most typical simulation results are shown for the case of ccrf discharge in argon at 50 mTorr pressure, with the interelectrode distance d = 8 cm, the voltage frequency 13.6 MHz and the voltage amplitude 120 V. The dust particle radius was chosen as r d = 100 nm, the dust density distribution was considered as constant. Two different dust density profiles were taken in our simulation: (a) one with a maximum in the center of the discharge gap, n d (x) = n max sin(π(x − l s )/(d − 2l s )), and (b) with two maxima near the electrode sheaths and a void in the gap center, which is often observed in experiments, n d (x) = n max sin 2 (2π(x − l s )/(d − 2l s )). The maximum of the density n max = 10 6 cm −3 , l s is the width of electrode sheath, where dust is absent (n d (x) = 0 for x < l s and x > d − l s ). The electrode sheath l s = 1 cm means the maximal sheath width in the cathodic phase of the discharge rf cycle. Together with a self-consistent determination of ϕ d (x) by equations (9)-(11), the traditional analytical estimation by the OML theory, obtained by solving the system of equations (7) and (8) for I e = I i , and called below ϕ OML , was also made for comparison. Figure 1 shows the densities of the charge carriers in the discharge plasma for both types of the dust density profile. It can be seen that the ion density is about equal to the sum of the densities of the electrons and the dust charge everywhere in the gap, except the electrode sheaths. This means that the central region of the discharge plasma remains quasineutral, and the negative charge of the particles is effectively screened by the ions. The shape of the electric potential of the plasma is similar for both dust profiles. The average electron energy in the gap center shown in figure 2 is 4.5 eV, and gradually decreases to 2.5 eV in the electrode sheaths. The ion energy in the quasineutral region is close to the thermal one of the neutral gas (kT i = 0.026 eV), but significantly increases toward the sheaths (see figure 2). The mean free path for ions l i under these conditions is near 1 mm, the Debye screening length λ D in the region containing dust is 0.02-0.06 mm. Figure 3 shows the volume density of the electron and ion currents onto the dust surface j e and j i . These values are the time-averaged amounts of the charge carriers captured by dust in 1 cm 3 per 1 s obtained by the direct count of MCCs (as described in the previous section) and from the equations (9) and (10). In the first case, the events of absorption with cross section (4) were counted during the simulation of the test particles motion for 300 rf cycles. Since these events are relatively seldom, the statistical noise is large despite the long averaging time. In the second case, the sums (9) and (10) are plotted, averaged for 10 rf cycles. The values of both sums practically coincide for the self-consistent determination of ϕ d . It is seen that the direct count of absorbtion events during the PIC-MCC simulation and equations (9) and (10) give equivalent current values, but the statistical noise for the second case is strongly suppressed. In figure 4, the self-consistent potential ϕ d is shown and compared with the analytical evaluation of ϕ OML . As the particle charge and potential are negative, the absolute value is plotted in figure 4 for convenience. The ϕ OML was obtained for the parameters shown in figures 1 and 2, with T e,i = 2ε e,i /3k. The disagreement between the self-consistent potential and ϕ OML is mainly due to the fact that the energy distribution functions f e,i are different from Maxwellian. In the central area of the discharge, where ions have thermal energy and Maxwellian energy distribution, the disagreement is minimal and does not exceed 0.3 V. The small local maximum on both curves in the gap center for the dust profile in figure 4(b) is explained by the minimum of the ion density where the dust has a void (see figure 1(b)) and coincides for both potentials. But near the electrode sheaths, where the electric field rises and the ion drift velocity u directed to the electrode increases (see figure 2), the self-consistent ϕ d and ϕ OML begin to diverge and at the edge of dust cloud the difference reaches 1.7 V or 25%.
The most probable reason for this disagreement is the increase of ion current to the dust in the case of fast ion drift (or large electrical field) in relation to the OML-based expression (8). This may be explained by the shape of the IEDF, which must be different from the shifted Maxwellian in the presence of ion resonant charge exchange. Since the resonant exchange produces ions with thermal energy, the abundance of low-energy ions having v u is much larger than the shifted Maxwellian IEDF gives. In order to illustrate this fact, in figure 5(a) the PIC-MCC and shifted Maxwellian IEDFs are shown for x = 1.1 cm. At this point, the ion drift velocity u = 500 m s −1 is two times larger than the thermal one √ kT n /m i = 250 m s −1 , where T n = 300 K is the neutral gas temperature (see figure 2). The shifted Maxwellian IEDF is plotted for the same u and T n . Both functions are normalized as f (E) dE = 1. For the PIC-MCC case the distribution: (i) is shifted to lower energy, and (ii) has a larger population of high energy ions. The combination of fast ion drift and continuous thermalization of the ions, due to the resonant charge exchange collisions, increases the ion current on the dust surface in comparison to the case of ion drift with the shifted Maxwellian distribution. For the equal ion flux toward the dust particle, in the first case the probability of ion absorbtion is much larger. Hence, the absolute value of ϕ d reduces toward the electrode sheaths in comparison with the behavior of ϕ OML .
It should be mentioned that the EEDF in this simulation was close to Maxwellian (see figure 5(b) and changed weakly over the region containing dust. So the electron current on the dust is adequately described by equation (7).
The incorrectness of the analytical evaluation of ϕ OML in the presence of ion resonant charge exchange was noticed in [11]. There the ion capturing by charged dust particles was simulated in the absence of an external electric field. It was also found that despite the condition of applicability of the OML approach being formally fulfilled (l e,i λ D r d ), the analytical evaluation of ϕ OML does not always provide an accurate value for the floating potential. In other work [12], where a direct current stratified glow discharge was considered, the dust potential was found using equation (8) for the ion current and equation (5) for the electron component,  with the electron distribution function obtained by solving the Boltzmann equation. Like in the present work, the non-Maxwellian approach gives up to 2 V difference from the OML results. Note that our model does not consider the detailed ion trajectories near the particles. The presence of dust in the discharge is described by giving the probability of electrons and ions capturing during their motion through the cross section σ id (4), which depends on the local dust concentration n d (x), potential ϕ d (x) and the instant electron/ion kinetic energy. For small dust particles r d λ D the capture cross section σ id does not depend on the details of the ion motion near the dust particle [10]. According to equations (5) and (6), the most important thing is the shape of the IEDF (assumed to be Maxwellian in the OML approximation).
Let us notice that in the central region of the discharge gap, where charge resonant exchange also occurs, the OML approach gives good agreement with our self-consistent calculation (see figure 4). This is explained by the fact that the ion drift in this region is slow (u √ kT i /m i ) and the IEDF is close to Maxwellian with T i = T n (neutral gas temperature), hence the resonant exchange does not influence the ion energy spectrum. But, when the ion drift velocity increases and u becomes comparable to the thermal velocity, the IEDF loses the Maxwellian shape and the disagreement between the PIC-MCC ϕ d and OML model results becomes visible.

Conclusion
A new method for calculation of dust charging in a ccrf discharge was developed based on the kinetic treatment of the electron and ion motion using the PIC-MCC method. The floating potential of the dust particles is calculated self-consistently with discharge dynamics through the integral balance of ion and electron currents onto the dust surface, taking into account the energy distribution functions of both the carriers. The PIC-MCC simulations show that the OML-approximation can give incorrect results even when the conditions of its applicability 12 (l e,i λ D r d ) are fulfilled. In the presence of ion resonant charge exchange collisions, the OML-approximation gives correct results only in a weak electric field, where the ion drift velocity is much smaller than the thermal one. With increasing ion drift velocity, the OMLapproximation gives up to 25% larger absolute value of the dust potential, than that obtained with the PIC-MCC calculations. This disagreement may be explained by the shape of the IEDF, which is different from Maxwellian or shifted Maxwellian. For the case of fast ion drift with resonant charge exchange, the abundance of low energy ions in the energy distribution shape is larger than that expected from a Maxwellian approximation, so the ion current on dust is enhanced in relation to the traditional OML-based expression. Our numerical approach allows to take into account the shape of the real energy distribution and its effect on the charging of the dust particles.