A fluid model for colloidal plasmas under microgravity conditions

A numerical model is presented to simulate steady states of complex (`dusty') plasmas under microgravity conditions. The model uses a fluid description for the plasma as well as for the dust microparticles. This is achieved by using an appropriate equation of state for the crystalline phase of the dust. The only forces assumed to act on the dust particles are the electric force, ion drag and pressure gradients. The model is used to study the formation of the stable `void' present in recent microgravity experiments. The structure of the dust clouds when the void is present is examined and explained.


Introduction
A dusty plasma is an ionized gas containing small particles of solid matter which acquire a large electric charge by collecting electrons and ions from the plasma. Dust particles can be efficiently confined in radio-frequency (rf) discharge plasmas. Because of the electrostatic interaction potential between the dust particles and of their confinement by the discharge electric field, these dusty plasmas may undergo phase transitions and condense to form 'liquid plasmas' and 'plasma crystals' [1,2].
In many experimental cases the dust particles have large enough mass to be affected by gravity which essentially confines the dust particles to a narrow 'plane' with a thickness of a few particle layers. In order to counter this effect experiments have been performed in microgravity conditions (see, e.g., [3,4]). In these experiments three-dimensional structures were observed but in addition a stable dust-free region ('void') was created in the centre of the discharge. The origin of the void and the nature of the forces that sustain it has been the subject of much debate recently; Goree et al [5] have studied through an analytic model the conditions for a stable void to appear for the case of collisionless ions. A similar approach, but for the collisional ions case, was taken by Tsytovich et al in [6]. However, their assumptions on the ionization profile in the discharge lacks self-consistency. A numerical approach was followed by Akdim and Goedheer [7] who used a fluid model similar to the one presented here and found that the dust void could only be created if the ion drag force was increased by at least a factor of five.
In this paper, we present a fluid model which we believe is an improvement over previous models in as far as the treatment of the dust particles is concerned. Furthermore, we show that the ion drag is sufficient to create the dust void and no other force is necessary.

Model description
The numerical model is based on the fluid approach and both plasma particles (ions and electrons) and dust particles are described in terms of the first moments of the Boltzmann equation. The semi-implicit technique of Ventzek [8,9] is used to couple Poisson's equation to the transport equations. The electron energy equation is interpreted using the implicit method of Hagelaar et al [10] whereas the boundary conditions are as in [11].

Plasma species
For each of the plasma species s the continuity equation can be written ∂n s ∂t + ∂j s ∂x = k i n e n n − β s n s n d , where n s and j s are the species density and flux, respectively. The first term on the right-hand side describes production with n e being the electron density, n n the neutral density and k i (¯ ) the ionization rate coefficient depending on the mean electron energy¯ . The second term describes losses to the dust particles with n d being the dust density and β s the collection rate coefficient defined as where V f = Q d /r d is the floating potential on the dust, Q d is the charge per dust particle, r d is the particle radius, k B denotes the Boltzmann constant and e the elementary charge. The particle flux is written as j s = ±µ s n s E − D s ∂n s /∂x, where µ s is the species mobility, D s = k B T s µ/e is the diffusion coefficient, E is the electric field and T s is the species temperature. An 'effective' electric field is used for the ions to account for inertial effects [12,13].
The ion temperature is assumed constant and the local field approximation is used so that the mobility is a function of the local electric field. For the electrons, however, this is not sufficient and the electron mobility as well as the ionization frequency are assumed to be a function of the local mean electron energy¯ [9]. This is obtained from the solution of the electron energy transport equation ∂n e¯ ∂t where the last term on the right-hand side represents electron energy losses due to collisions with neutrals and K l is an effective energy loss coefficient dependent on¯ . Energy losses due to absorption by the dust particles are not included since they have been found to be negligible.

Dust particles
The dust is also modelled as a fluid. It is assumed that there are no sources or sinks of dust and the initial dust population is maintained at all times. The forces assumed to act on the dust are the electric force and the ion drag. Neither gravity nor thermophoresis is taken into account. Equation (1) (without the source and sink terms) is used to solve for the dust density. The dust flux is written as where F ion is the ion drag force and D d , µ d ∝ ν −1 m are the effective dust diffusion and mobility coefficients. The dust collision frequency, ν m , can be arbitrarily chosen-we are only interested in steady state solutions (where j d = 0), and ν m essentially sets the timescale for the dust motion. This is the reason why we can also neglect the influence of the dust viscosity. The mobility is taken to be µ d = Q d /(m d ν m ). In order to obtain the charge on each dust particle we use the orbital-motion-limited (OML) theory which is usually valid if the particle radius is smaller than the linearized Debye length, where λ e and λ i are the electron and ion Debye lengths, respectively. The charge on the dust is changing dynamically according to the collection rate coefficients given by equations (2) and (3). The determination of the effective diffusion coefficient D d is more complicated and will be treated in the next section.

Equation of state (EOS) for the crystalline phase and determination of the diffusion coefficient of the dust
The effective dust diffusion coefficient D d cannot be obtained from the dust mobility using Einstein's relationship for an ideal gas. This is because dust particles in microgravity and laboratory experiments are normally in strongly coupled states. In principle, the determination of the diffusion coefficient requires knowledge of the equation of state (EOS) for the dust system in different phases. In this paper, it is assumed that the ideal EOS is valid for the gas state of the dust and a first attempt to solve for EOS for the crystal case is provided. The EOS for the liquid phase can be derived using, e.g., the model of hard spheres (see, for example, [14]).
We will now proceed with a calculation of the EOS for the dust in crystalline structures. Given a system of N d particles with temperature T d , its free energy is given by [15] where ω i = ω i (n d ) is the frequency of the ith normal mode of the crystal. Considering that the vibrational part of the free energy is small compared to the coupling part N d E 0 the pressure is given by is the dust density. The coupling energy per particle in a crystal can be well approximated by the Yukawa coupling energy between the particle and its nearest neighbours. That is, is the mean inter-particle distance, N nn is the number of nearest neighbours (N nn = 8 and α = 3 1/2 /2 2/3 1.09 for bcc lattices, N nn = 12 and α = 2 1/6 1.12 for fcc and hcp lattices [16]), and λ D the Debye length. Finally we obtain for the crystalline pressure, where ≡ Q 2 d /4π 0 k B T d is the coupling parameter, κ = /λ D is the lattice parameter and P g = n d k B T d is the gas pressure. In the above we have assumed that ∂Q d /∂n d and ∂λ D /∂n d are negligible which is reasonable when the Havnes parameter [17] does not exceed unity.
In most cases particles form crystalline structures in the steady state. The simplest approach to derive the appropriate EOS for our simulations is to interpolate it from the gaseous state to the crystalline state. The width of the transition region is a function of the coupling parameter . In the weakly-coupled case ( < 1) the gas state equation P = P g is used. When > cr , where cr (κ) is the crystal-liquid curve, we apply P = P cr from equation (4). Finally, for 1 < < cr a linear interpolation is performed so that P is equal to P l = P g (1 − x) + P cr x, where x = ( − 1)/( cr − 1). The crystal-liquid curve can be very well approximated by [18] cr (κ) = 106 exp(κ)/(1 + κ + κ 2 /2). With the above approximate but complete EOS one can 'define' the effective diffusion coefficient to be

Ion drag force
An accurate determination of the ion drag force is important for the study of the void structure.Expressions for the ion drag have been provided by several authors [19]- [22]. However, Khrapak et al [23] have recently shown that the calculation of the Coulomb logarithm for the case of ion-dust collisions is not carried out properly. In their approach the Coulomb logarithm is given by where ρ C (v) = eQ d /4π 0 m i v 2 is the Coulomb radius (v is the ion velocity). Given (5) the ion drag can be calculated by evaluating the following integral over the ion velocities where σ s (v) = 4πρ 2 C (v) C (v) is the orbital cross section, u = j i /n i is the (average) ion velocity and v T = √ k B T i /m i is the ion thermal velocity. The difference between equation (6) and usual expression for the ion drag force (e.g. [19]) is the following: the Coulomb radius ρ C ('length of interaction') for the ion scattering can be very large, because of the high value of the dust charge Q d . In typical laboratory dusty plasmas, ρ C exceeds the Debye length λ D . Then the usual Coulomb scattering approach with the cut-off at λ D (which is used in [19]) is no longer applicable-most of the momentum transfer is associated with ions scattering outside the Debye sphere. Equations (5) and (6) take these 'far collisions' into account.

Model results
Since we are not interested in the dynamics of the initial void creation itself (which might be a kinetic effect that cannot be present in a fluid representation) but rather in whether the conditions exist to sustain a stable void or not, our initial dust configuration is such that a void already exists (see figure 1). Then, the model will gradually converge and the void will close or converge to a steady state void.
A typical case where no void is present at steady state is shown in figure 2. The discharge parameters for that run are a gap of 6 cm between electrodes, a driving voltage amplitude of 80 V at 13.56 MHz and argon gas at a neutral pressure of 0.5 Torr. The particle radius is taken to be r d = 2 µm and the kinetic temperature of the particles is assumed to be room temperature. At steady state the charge per dust particle is approximately equal to 1.2 × 10 4 elementary charges.  If we lower the neutral pressure the ion drag can overcome the electric force that acts on the particles and a stable void can be maintained. This happens because the ion flux in the void pointing towards the cloud gets higher as the pressure goes down, whereas the electric field in the void does not change significantly. Such a situation is shown in figure 3 where the conditions are the same as before, except that neutral pressure has been set to 0.4 Torr.
An interesting feature that appears is the sharp void-dust interface; the transition from the void to the dust cloud is abrupt, the dust density rapidly increasing near the boundary and decreasing again inside the dust cloud. This peak is not an artifact of the model-the same is often observed in experiments, where at the void-dust interface the dust appears a few times denser in a very narrow layer of thickness equal to about 2-3 inter-particle distances [3,4]. The reason for that is that ions get trapped in a potential well that exists next to the void-dust interface (see figure 4) and get absorbed by the dust. This inverts the direction of the ion flux (see figure 5)  in the vicinity of density peak so that it points to the peak. Thus, the ion drag, being superior to the electric force (which also changes sign), compresses from both sides the dust population near the interface and causes the jump in the density. This is illustrated in figure 6. Note also that under the conditions presented here the dust liquifies only in the peak region and remains in gas phase everywhere else. Thus, at least for these conditions, the only effect that the EOS has on the model is in limiting the dust compression at the interface and plays no significant role for the rest of the dust cloud.
Furthermore, let us comment that because in the vicinity of the density peak the ion flux points towards it, it is implied that there has to be an increased ionization rate in the region of the increased dust density to balance the ion losses to the dust. This indeed appears in the model, and is the result of electron heating by the electric field close to the density peak which is enough to compensate the reduced electron density in the region and increase the ionization. This result questions the validity of the assumption in [5,6] that ionization is null in dust cloud and instead takes place in the void region.
Finally, if the pressure is further decreased, increasing the ratio of the ion drag to the electric force, a situation is reached where the dust is concentrated in a narrow region adjacent to thesheaths and the density profile forms a single peak.

Conclusions
A self-consistent numerical model has been developed in order to simulate a colloidal ('dusty') plasma in microgravity conditions. The model uses a fluid description for the plasma as well as for the dust microparticles. This has been made possible by deriving an approximation to the equation of state of a 'liquid' dust cloud. The forces acting on the dust which are accounted for are the electric force, ion drag and the force due to pressure gradients. The fact that the void can be reproduced by the model indicates that ion drag alone is sufficient to sustain it and no other force (such as thermophoresis) is needed.
The structure of the dust cloud obtained by the model has been examined and has been found in good agreement with microgravity experiments [3]. Contrary to what is commonly stated, the ion drag does not necessarily always point towards the electrodes, because of ion trapping inside the dust cloud. This also explains why at the void-dust interface there is a very narrow region of increased dust density.
Despite the improvement over previously presented fluid models [7] and the success of the model in interpreting some of the features of the dust cloud, there are some limitations to the present approach; purely kinetic effects which might be important for the further understanding of dust dynamics inside the plasma cannot be included in a fluid representation without some difficulty. In addition, a detailed study of the phenomena occurring at the void-dust boundary requires models capable of dealing with distances smaller or of the order of the average dust inter-particle distance. To this goal, we believe that research using kinetic models (such as models based on the particle-in-cell method) will be useful in the future.