Modelling of two differently sized dust species in plasmas under micro-gravity

A self-consistent two-dimensional hydrodynamic model for a dusty argon plasma has been developed to model more than one dust species. Results are presented for situations where dust particles with two different diameters have been included. The final steady state solution is achieved after three injection steps of the dust particles. At every injection phase dust particles of only one size are let in and the simulation is continued until the steady state solution is achieved. Results show that the differently sized dust particles form crystals at different positions. These dust clouds have an influence on each other by means of positive space charge layers created due to the argon ions, which cannot match the steep dust crystal boundaries. The screened Coulomb interaction between the two differently sized dust species is neglected. The electric potential, ion density, electron density and electron energy show significant changes after each injection phase, even at an amount of dust that is small compared to that studied during the micro-gravity experiments.

. Schematic diagram of the PKE reactor.

Introduction
Dusty plasma experiments performed under micro-gravity conditions have shown interesting phenomena such as voids surrounded by crystalline regions and rotating dust clouds which usually appear outside the crystalline regions. In their PKE (Plasmakristall-Experiment) chamber (figure 1) Morfill et al [1,2] usually injected uniform size dust particles under micro-gravity. In an earlier paper [3], we used our dust-argon hydrodynamic model to study these plasma crystal experiments. Our modelling results showed good agreement with the experimental results. This offers a sound basis to use our model as a predictive tool.
Theoretical and numerical studies up to now have basically followed dust particles in the electric field and particle fluxes of an undisturbed discharge. An important aspect not covered is the influence of the dust particles on the discharge due to recombination on their surface [4,5]. In discharges containing a considerable amount of dust, as in the PKE experiments, this approximation is not correct. For this a fully self-consistent model is needed. We have developed such a model for a dust containing radio-frequency (rf) discharge in argon and used it to investigate the behaviour of dust particles. The model contains a dust-argon fluid part which has been described in a previous paper [3]. We have extended the single-dust fluid model to a multi-dust fluid and studied the behaviour of plasmas containing dust particles of different sizes. Our model accounts for the influence of the different dust fluids on the plasma and for their transport as fluids. This makes it a sophisticated tool for studying dusty discharges which contain dust particles of different sizes. In this paper we describe the multi-dust fluid model and present the results obtained from it.

Dusty plasma fluid model
In the fluid part of our two-dimensional model the particle balances, the electron energy balance and the Poisson equation are solved, including the transport of the dust fluids. Further details about the algorithms and the set of equations solved in the argon fluid part can be found in [6]. Drift-diffusion expressions are used for the fluxes and in the electron energy balance only Joule heating and energy in inelastic collisions are taken into account. Problems related to the huge difference in the timescale of the dust motion (1-10 s) and the rf period (70 ns) have been solved by a time-splitting technique and an iterative procedure. The charge on a dust particle is calculated from the orbital-motion-limited (OML) probe theory. The (constant) charge on 20.3 the dust particle is obtained from the balance of the local (OML) electron and ion currents collected by the particle. Recombination of ions and electrons on the dust particle surface is also taken into account. Ion-neutral collisions have been included to simulate a possible gas heating mechanism. For this we have used a simple approximation by assuming that the energy taken up from the electric field by the ions is dissipated locally in collisions with the gas [7]. The heating of the dust particle surface is also taken into account and more details about the expressions used can be found in the following papers [8,9]. The dust particle (surface) temperature can affect the gas temperature which, in turn, could affect the other elementary processes in the discharge. The power released by the ions and the dust is the source in the temperature balance of the gas [9]. To obtain expressions for the dust particle fluxes, several forces acting on the dust particle are taken into account, e.g. the ion drag force, the neutral drag force, the electric force and the thermophoretic force.
A dust particle in a discharge is charged up to the floating potential relative to the surrounding plasma. This potential depends on the local ion and electron density and energy distribution. For a spherical dust particle with a radius r d , the OML theory [10] predicts a positive ion and electron current: Here, n e is the electron density, n i is the positive ion density, e is the elementary charge, k B is Boltzmann's constant, T i is the positive ion temperature, T e is the electron temperature, m i is the ion mass, m e is the electron mass and V f is the floating potential. All species are assumed to have a Maxwellian energy distribution. When the ions enter the plasma sheaths near the electrodes, they get a directed velocity v i due to the electric field. Therefore, we have replaced k B T i in the expression for the ion current by the mean energy E i , which is Equation (3) is obtained by using the mean speed expression of Barnes et al [11] given by By calculating 1 2 m i v 2 s , equation (3) is obtained. The directed velocity v i is calculated by solving the drift-diffusion equation for the ions.
In the model the charge Q d = 4π 0 r d V f on the dust is calculated by equating these currents. The floating potential of the dust is assumed to be constant during an rf cycle. This assumption is justified by the fact that the currents toward the dust surface are too small to change the charge significantly during an rf cycle.
When a dust particle becomes negatively charged, it will attract positive ions; these will recombine with an electron that has to be replaced again by an electron from the discharge to maintain the floating potential. As a result the equilibrium fluxes of positive ions and electrons 20.4 arriving at the dust surface will recombine and the energy released is used to heat up the dust particle surface [8,9]. The electron flux (equation (2)) results in a recombination rate In a plasma dust particles undergo a wide variety of forces. Assuming that a dust particle is a perfect sphere the gravitational force can be written as where r d is the dust particle radius, ρ d is the mass density and g is the gravitational acceleration.
For the often used melamine-formaldehyde dust particle ρ d is approximately 1.51 × 10 3 kg m −3 . When a dust particle gains a certain velocity relative to the neutral gas, it will experience a drag force due to momentum transfer from/to the gas. This neutral drag force has been discussed in detail by Graves et al [12]. It can be approximated by the following equation: where n n is the density of the neutral with mass m n , v d is the drift velocity of the dust particle and v th is the average thermal velocity of the gas. Because advection of the neutral gas is not included in the model, this force will only be present as a damping force on the velocity of the dust particles. The assumption is made that this damping force will be in equilibrium with the sum of all the other forces acting on the dust particle.
Another force caused by momentum transfer is ion drag. This force results from the positive ion current that is driven by the electric field. It consists of two components. The collection force represents the momentum transfer of all the ions that are collected by the dust particle and is given by where v s is the mean speed of the ions, v i is the ion drift velocity and b c is the collection impact parameter. The second component is the orbit force given by with b π/2 the impact parameter that corresponds to a deflection angle π/2 and the Coulomb logarithm: where λ L = ((1/λ e ) 2 + (1/λ i ) 2 ) −1/2 is the linearized Debye length, which is a combination of the electron Debye length, λ e , and the ion Debye length, λ i . The ion drag is discussed in more detail by Barnes et al [11]. Due to their charge, dust particles will experience an electric force. Daugherty et al [13] derived the following expression: 20.5 where Q d is the charge on the dust particle, E is the electric field and κ = 1/λ L . In a discharge the dust particle radius is much smaller than the linearized Debye length. Therefore the term between the brackets is approximately 1 and the electric force is given by This expression holds for situations where the dust particles are not shielded from the plasma by positive ions trapped in orbitals around the dust particle [14]. In that case the particle will behave as some kind of dipole. When a temperature gradient is present in a discharge, for instance due to cooling or heating of the electrodes a third force driven by momentum transfer will occur. This force is called the thermophoretic force. Atoms impinging from the hot side have more momentum than their companions on the cold side: this can result in a force pointing in the direction − ∇T .
For large Knudsen numbers Talbot et al [15] derived the following expression: where v th = [8k B T gas /(πm)] 1/2 is the average thermal velocity of the gas and κ T is the translation part of the thermal conductivity. α, the thermal accommodation coefficient of the gas, is taken equal to 1. For a stationary solution the five forces calculated in our model (the neutral drag, the electric force, gravitational force, thermophoretic force and the ion drag) are in equilibrium. With the introduction of a momentum loss frequency and a mobility and diffusion coefficient for the dust particles given by where p tot is the static pressure and m d is the dust particle's mass: it is possible to define a 'drift-diffusion' expression for the flux of the dust particles and treat them with the same numerical procedures as the other charged particles in the fluid model. Hence, in this case, we neglect the inertia of the dust particles. Because of the low mobility of the dust particles the effective field E eff can be approximated by the time-averaged rf field. With the use of the Einstein relation, the diffusion can be defined as a function of the mobility (equation (16)).
The drift velocity and the diffusion coefficient of the dust fluids are much smaller than those of the ions and electrons. Therefore it would take rather a lot of computational effort to achieve a steady state solution for the dust when it is followed during an rf cycle. We therefore have To solve this problem, we correct the artificial space charge by adapting the positive ion density distributions prior to the next series of rf cycles. With this method the required speed-up is established. Both for the plasma species and the dust fluids the transport equations are solved using the Sharfetter-Gummel exponential scheme [6]. The internal pressure of the crystal due to the inter-particle interaction has been included by means of a density dependence of the diffusion coefficient for the dust. The diffusion coefficient of the dust is increased by a factor exp(N d /N c ), where the reference density N c is chosen such that the dust density saturates at a value N crys . This models the incompressibility of the crystal. Actually, the (yet unknown) equation of state of the dust crystal should be used to account for the internal pressure. Since we were not primarily interested in the precise structure of the crystalline regions, we have chosen the simple and computationally robust exponential increase of D d . Plasma crystal experiments [1] have shown an inter-particle distance of about 300 µm for a dust particle with diameter of 15 µm. This results in an average 'crystal' density N crys of 3.7 × 10 10 m −3 . For the smaller dust particles we used a reference density that is a factor of 200 higher than for the bigger dust particles.

Results and discussion
We have modelled the PKE chamber used by Morfill et al [1]. The reactor is cylindrically symmetric. The simulation starts with zero dust density profiles. During the simulation the dust is injected from both electrodes by adding source terms in the dust particle balances for the first grid points below/above the electrodes: the injection rate is about 250 000 particles s −1 . Both electrodes are driven by a rf power source at a frequency of 13.56 MHz. The peak-to-peak voltage is 70 V, which results in a power dissipation of about 0.04 W. The pressure is 40 Pa. The equation of motion for the dust fluids (equation (17)) is solved for the time-averaged electric field, plasma densities and fluxes. First the simulation is started for a dust-free discharge and is run until it reaches steady state. After that, the total amount of dust is injected in three steps. These steps are as follows: in step 1, the injection of dust particles that have a diameter of 15 µm is started until an amount of 350 000 particles is reached; after that, the simulation is run to achieve a steady state solution. In step 2, 350 000 dust particles of 2 µm diameter are injected into the above steady state solution. This gives a mixture of small and large dust particles [16,17]. This simulation is also run to the steady state. Eventually, in step 3, a further 350 000 dust particles of 2 µm diameter are injected. This results in a total amount of 1.05 million dust particles. After that the sources are switched off and the simulation is continued until it reaches the final steady state. Figures 2(a)-(c) show the dust densities of both dust species (large and small) after each injection step. After injection step 1, the larger dust particles arrange themselves in a crystal around a central void [1]. The crystal is formed due to the balance of forces acting on the dust particles. When the smaller dust particles are injected in step 2, the crystal consisting of the 15 µm dust particles moves radially outwards towards the reactor wall. Also the crystalline regions between the electrodes become narrower due to an increase of the ion drag force and electric force acting on the larger dust particles. The smaller dust particles are pushing the larger particles to move radially towards the reactor wall. After step 3, this effect becomes more obvious. The shape of the crystal of the larger particles seems to be dictated by the smaller dust particles' cloud and the dust crystal is pushed more in the radial direction towards the wall. This movement can be explained by the presence of the smaller dust particles and the space charge layers which are formed by the presence of the dust crystals. Figures 3(a) and (b) show the space charge and dust density profiles along the axial symmetry axis (z = 0.027 m). It can be seen from figure 3(b) that, between the two dust density peaks, a space charge layer is created. The line integrated number of dust particles of 2 µm size along the axial symmetry axis is about a factor of 40 larger than for the 15 µm size particles. The difference in the charge on the dust particles (figures 5(a) and (b)), however, is approximately a factor of 6. To maintain quasi-neutrality the charge of the 2 µm dust particles must be compensated by at least a factor of 6 times more ions. The increase in ion density is shown in figure 4  acting on the crystal of larger particles results in the dust crystal moving radially towards the reactor wall. Figures 6(a)-(c) show the time-averaged potential distribution V (r, z) for three different situations. In the three cases the potential has its maximum in the bulk of the plasma between the electrodes. Comparing the potential distributions, a significant change in the plasma potential can be observed. The plasma potential decreases in the centre of the discharge after injection step 1 ( figure 6(b)). Also a decrease of the axial electric field can be observed. After step 3, the electric potential increases and its maximum becomes almost equal to that for the dust-free 20.11   formed by the larger dust particles ( figure 2(a)). This drop in electron density is caused by the charging of the dust particles and recombination of ions and electrons on the dust surface. After step 3 the electron density increases, which seems quite strange because the amount of dust has increased. The explanation for this phenomenon can be found from figure 2(b). After injection 20.13 step 3 the dust crystal formed by the 15 µm particles is pushed radially towards the reactor wall. As a result the dust crystal between the electrodes becomes thinner and moves to a position where there is less plasma. This results in less recombination of ions and electrons on the dust particle surfaces. The same effect is observed for the ion density: after step 3 ( figure 8(b)) the ion density shows a significant increase compared to the ion density after step 1 ( figure 8(a)). The ion density profiles have peculiar shapes. After injection step 1 the ions around the midplane in between the electrodes are pushed towards the centre. After step 3 ( figure 8(b)) the opposite effect can be observed. These outer contour lines in the ion density profiles can be explained by examining figures 2, 3 and 9. These figures show the space charge regions in the plasma and the dust density profiles. After injection step 1, a dust crystal is formed by the larger particles. To maintain quasi-neutrality the ions tend to compensate the charge carried by the dust particles. This results in an ion density profile which has the peculiar shape as shown in figure 8(a). The steep boundaries of the dust crystal give rise to the creation of space charge layers [9] (figure 9). After injection step 3, the space charge layer appears ( figure 3(b)). This space charge layer accelerates the ions in the direction of the reactor wall. This acceleration can be observed in the outer contour line of figure 8(b). Figures 10(a) and (b) show the time-averaged electron energy. A comparison of the figures shows that the electron energy after injection step 3 has decreased. This effect can be explained by an combination of the increase of the electron density (figures 7(a) and (b)) after injection step 3 and the constraint of constant applied driving voltage: more electrons mean less energy per electron. Therefore the time-averaged electron energy decreases.

Conclusions
The numerical simulation results show that self-consistent modelling of the plasma parameters in a dusty argon discharge is important. We have shown that injection of dust particles with two different diameters results in the formation of dust clouds at different positions, with the cloud of smaller dust particles inside the central void of the crystal of the larger particles. Neglecting the screened Coulomb force between two dust clouds, these clouds can still interact with each other through the formation of space charge layers, which are formed by the steep boundaries of the dust crystals.