Low-frequency electrostatic modes in partially ionized complex plasmas: a kinetic approach

A kinetic model of partially ionized complex plasmas is employed for the numerical analysis of low-frequency longitudinal modes for typical laboratory plasmas. The approach self-consistently includes the effects of plasma particle absorption on dust, collisions with neutrals and electron impact ionization. In addition to the typical dust acoustic mode, the results reveal the existence of a novel long-wavelength mode, attributed to the interplay between the mechanisms of plasma production and loss. The main properties of mode dispersions are investigated through their dependence on plasma and dust parameters.


Introduction
The dust acoustic (DA) mode, a manifestation of the slow time scales of dust dynamics, as compared to ion time scales, has been extensively studied in both laboratory and micro-gravity conditions owing to the possibility of direct observations on a particle level (see, e.g., [1][2][3][4]).
The mode was predicted theoretically by Rao et al [5] through the analysis of the collisionless fluid equations for wave-like perturbations. Later works addressed this problem, including the effect of neutrals [6,7], dust charge variations [8][9][10] and particle drifts and ionization [11][12][13][14], by linearizing the collisional fluid equations. Such hydrodynamic models have been the most widely used descriptions for the interpretation of related experiments.
Kinetic treatment of DA waves has been addressed in the works [15,16] through the multicomponent approach, which implies fixed dust charge and omission of dust charge fluctuations and absorption of plasma fluxes on dust grains. In [17], the standard kinetic model was used, where dust charge fluctuations are also taken into account. Recently, the self-consistent 'full' kinetic model, originally developed in [18][19][20], was employed to investigate DA waves in fully ionized plasmas [21]. To implement the results of kinetic models for the analysis of experiments in typical low-temperature discharges [3], an extension to account for neutrals is essential.
In this paper, extending the work of [21] to include the effect of neutrals, the low-frequency longitudinal dusty plasma modes are studied employing a kinetic model of partially ionized complex plasmas [22]. The model considers the absorption of plasma fluxes on dust particles and the dust charge fluctuations self-consistently while taking into account collisions with neutrals via the Bhatnagar-Gross-Krook (BGK) approximation [23] and using a fluctuating plasma source to sustain the plasma. The total permittivity of the system is evaluated using the low-frequency responses that were derived in [24], encompassing the effect of neutrals on grain charging and with the assumption of Maxwellian distributions. 3 Numerical solutions of the kinetic dispersion equation are presented for typical laboratory plasma parameters. The results show the emergence of a novel long-wavelength (LW) mode in addition to the DA mode. The LW mode is present only at large wavelengths for strong ionization and is related to the interplay between absorption and ionization. The main characteristics of the LW mode are elucidated through the dependence of the dispersion relations on plasma and dust parameters. Based on the numerical results, the optimal conditions for experimental observation of the LW mode and for the observation of kinetic effects on the DA mode are discussed.

Theoretical aspects
In the framework of the linear fluctuation theory [25,26], the starting point is the Klimontovich equations of all species [27], which are decomposed in ensemble averaged and fluctuating parts to yield: the Boltzmann kinetic equations that contain the modified collision integrals and also the permittivity ( k,ω ) of the system as a function of the low-frequency responses. For the full kinetic model the low-frequency responses have been derived in [24] in a form ready for numerical evaluation with the additional assumption of Maxwellian distributions and taking into account the effect of neutrals on the ion capture cross-sections. In this work the dispersion equation k,ω = 0 is solved numerically and the ω r (k), ω i (k) relations for the low-frequency longitudinal modes are obtained.
In this section, before proceeding to numerical results, we discuss the conditions of validity of the kinetic model, the models for the description of non-binary collision processes and the analytical form of the permittivity, the method of numerical solution of the dispersion equation and the existence of multiple roots of the dispersion equation.
In this paper the following notation is used: e for the elementary charge and q for the dust charge, q eq = −eZ d for the equilibrium dust charge with Z d the dust charge number, a for the dust grain radius, m α , n α , T α for the species mass, density and temperature, respectively, v T α for the thermal velocities, ν n,α for the average collision frequencies with neutrals and λ Dα , ω pα for the species Debye lengths and plasma frequencies.

Applicability of the kinetic model
The kinetic model is a linear model; for multi-component plasmas this implies a first order expansion in 1/N Dα , where N Dα = n α λ 3 Dα is the number of particles of a species α in the Debye cube, meaning that all plasma components should be in the gaseous state [28]. Such a condition for dust, assuming that the mean intergrain distance is d n −1/3 d and that the effective screening length in small distances is λ λ Di , can be written as where is the coupling parameter. For the parameter regimes investigated here, this condition is approximately satisfied.
The kinetic model takes into account only fluctuating fields, which implies that the system must be homogeneous and that there is no external field. Hence, all species are assumed nondrifting. It is worth pointing out that it is not self-consistent to include drifts just by assuming shifted Maxwellians; the existence of an external field alters the fluctuating parts and therefore the expression of the permittivity.
The kinetic model considered is a low-frequency model, ω ku s , where u s is the ion sound speed; this is imposed not only by the assumption that electrons/ions are treated as continuous 4 Vlasov fluids and only the dust discreteness is taken into account, but also by the adiabatic description of electrons.
Furthermore, there is an ad hoc description of collisions with neutrals through the BGK collision term, which is a relaxation time approximation of the exact Boltzmann collision integral [29]. The collision operator uses a single relaxation time 1/ν n,α for a non-equilibrium distribution function to evolve into an equilibrium one. However, in partially ionized complex plasmas there are multiple collision processes with different relaxation times leading the system to equilibrium. As the effect of neutrals becomes stronger, the approximation becomes more exact.
The dust density parameter P(= n d Z d n e ) must be large enough for electron/ion binary collisions to be neglected with respect to collisions with dust, P Z d > 1 [18]. For the parameters studied here, this criterion is always satisfied.
Finally, the low-frequency responses are computed with the assumption of Maxwellian distributions instead of the solutions of the Boltzmann equations.

Description of collision processes
In the general formulation of the kinetic model the description of capture cross-sections, collisions with neutrals and electron impact ionization is left arbitrary.
For the dust-neutral collision frequency we have ν n,d = 8 √ 2π 3 δa 2 m n m d n n v T n [30] with δ = 1 + π 8 [31]. For the ion-neutral collision frequency ν n,i = n n v T i σ n,i , where σ n,i is the velocityindependent collision cross-section of the dominant resonant charge transfer collisions [32]; typical values are, e.g., σ n,i = 10 −14 cm 2 for neon gas and σ n,i = 2 × 10 −14 cm 2 for argon.
For a qualitative treatment of electron impact ionization, the classical Thompson model [32] of ionization cross-sections will suffice. The instantaneous ionization frequency is ν I (v) = n n σ i z (v)v and the average ν e = 1 where Ei(z) = ∞ z e −t t dt is the tabulated exponential integral. Here U i z is the first ionization energy of the neutral gas, which is 15.76 eV for argon and 21.56 eV for neon.
For the capture cross-sections of plasma particles on dust: in the weakly collisional regime [3,33], the electrons are considered to be collisionless and their capture cross-sections are given by the orbit motion limited (OML) model, σ e (q, v) = πa 2 (1 + 2qe am e v 2 ) [3], while for the ions, considering the effect of charge transfer collisions inside the perturbed radius R 0 , the cross-sections can be approximated by [24,34] where l n,i = 1 n n σ n,i is the mean free path in ion-neutral collisions. A self-consistent computation of R 0 (q) is provided in [35], where it is shown that R 0 and the effective screening length λ are functions of the parameter β = z τ a λ Di only, R 0 (β) = (−0.1 + 0.8β 1/2 )λ Di and λ(β) = (1 + 0.2β 1/2 )λ Di , with τ = T i T e being the plasma temperature ratio and z = Z d e 2 aT e the dimensionless charge.

5
Finally, the average absorption frequency of ions on dust (average capture frequency), defined by the average of the instantaneous ion capture frequency ν d,

The permittivity k,ω
The permittivity derived from the kinetic model of partially ionized complex plasmas, described in [22,24], has the following form: The above auxiliary responses are algebraic functions of the low-frequency responses of the medium and are defined by The expressions for the low-frequency responses, namely χ d,eq , were originally derived in [24] under the assumptions listed in section 2.1 and for the cross-sections given in section 2.2, in a form ready for numerical calculations. For completeness we provide their explicit forms in the appendix.

Solution of the dispersion equation
The longitudinal dispersion equation, k,ω = 0, is to be solved in its general form { (ω r , ω i , k)} = 0, { (ω r , ω i , k)} = 0 when allowing for temporal attenuation of waves. This 6 is a coupled system of two nonlinear equations and is a formidable numerical task. However, considering the physical mechanisms of energy input and dissipation allows us to devise approximate yet accurate methods of solution.
The damping mechanisms in the low-frequency partially ionized plasma system are as follows. (i) Landau damping due to wave resonance with the dust and ion species. Landau damping on dust is always very weak, unless the dust temperature greatly exceeds the ion temperature by several orders of magnitude and the wavenumbers of interest are very large; only then does the ratio v ph /v Td reach unity, with v ph = ω/k the phase velocity. Landau damping on ions is also very weak, since the ions are assumed to be non-drifting. (ii) Damping due to plasma absorption and charge fluctuations. It is important mostly for very small wavenumbers that are comparable to the mean free path in plasma capture collisions and very small frequencies (so that the charge varies in the characteristic time scale of the phenomenon). (iii) Damping due to collisions of the heavy species with neutrals. Because of the use of the BGK description, temporal damping due to this effect is equal to the hydrodynamic result ω i = − ν n,d 2 and is wavenumber independent.
On the other hand, the only source of free energy is electron impact ionization of neutrals, which is not very high in low-temperature plasmas (T e 5 eV). Hence, it is expected that the low-frequency waves are always damped and ω i < 0.
The strength of dissipation by neutrals dictates the method of solution of the dispersion equation; the dust neutral collision frequency is a function of the gas pressure (P n ) and the radius of the dust grains.
• Very low pressures (below 1 Pa for a = 1 µm): the total dissipation is low, |ω i | ω r , and the weak growth rate approximation is applicable, i.e.  k). Therefore, we can confine the values of |ω i | for which { k,ω } = 0 can have a solution to a small interval {|ω 1 |, |ω 2 |}, then pick any ω * i in this interval and solve the equation { (ω r , ω * i , k)} = 0 independently with negligible error.

The dust acoustic mode and the long-wavelength (LW) mode
Preliminary results of [24] have shown that { k,ω } can have two roots. The root corresponding to smaller wavelengths is identified as the DA mode, since as in multi-component kinetic models, the sign of { k,ω } changes from positive to negative over this zero crossing. The root 7 corresponding to larger wavelengths is a new LW mode; over this root the sign of { k,ω } changes from negative to positive, which implies that { 0,ω } can be negative for partially ionized complex plasmas. The LW mode is absent for non-fluctuating plasma sources and also for multi-component kinetic models. Therefore, it can be attributed to the competition between electron impact ionization and plasma absorption on dust.
Indeed for neon gas, due to the high ionization energy (21.56 eV) the electron impact ionization frequency is negligible and the LW mode does not exist for a plasma discharge with the typical T e = 3 eV. Meanwhile, for argon gas, whose ionization energy is much lower (15.76 eV) and thus the energetic tail of the Maxwellian electrons leads to an ionization frequency of the order of the plasma capture frequency, the LW mode emerges.
It should be noted that the LW mode is not a pure kinetic mode; similar modes have been predicted by Tsytovich and Watanabe [14]. In that work, a self-consistent hydrodynamic model is used and the new modes are attributed to the perturbations of the complex plasma equilibrium state, which should be defined by quasineutrality, plasma flux balance on the grain and plasma generation-loss balance; hence the necessity to include a plasma source depending on the physical scenario. In that work, as in our model, the source is electron impact ionization of neutrals.

Numerical results
In the numerical calculations, n i = 10 9 cm −3 , T i T d T n 0.03 eV and the mass density of the spherical dust is ρ m = 1 g cm −3 .

Kinetic dispersion relations
To elucidate the effect of neutrals on the dispersion, it is instructive to use neon gas and low electron temperatures, because in this case only the DA mode exists. For a = 1 µm, P = 1, T e = 3 eV and varying pressure from 0 to 15 Pa in figure 1, we notice small differences in the dispersion relation, while for varying pressure from 20 to 60 Pa, we notice the introduction of wavenumber cut-offs at zero frequency, with a value increasing with the increase of P n . Such cut-offs can also be predicted from a simple hydrodynamic model taking into account only collisions with neutrals, but not charge fluctuations, absorption on dust and ionization. After the standard analysis of the hydrodynamic equations for wave-like perturbations [3], the dispersion relation for the DA mode will then read, for k > k crit = ν n,d 2ω pd λ D , where λ −2 D = λ −2 Di + λ −2 De , while for k < k crit it will be purely imaginary (aperiodic). However, the hydrodynamic critical wavenumbers will be much larger than the kinetic theory ones, obviously due to the effect of the neglected processes.
By increasing the electron temperature still considering neon gas, the electron impact ionization of neutrals becomes stronger and the LW mode appears. In figure 2, the dispersion relations for P = 1, a = 1 µm, P n = 20 Pa and varying electron temperature from 3 to 10 eV are plotted. At 5 eV and above, the LW mode appears and as the ionization frequency rises it gets less and less steep. It is characterized by large phase velocities (50 cm s −1 and higher)  compared to the DA phase velocities (less than 10 cm s −1 ) and also by a negative group velocity v ph /v g < 0 with v ph = ω r kk and v g = ∂ω r (k) ∂kk (backward wave), implying that the energy and disturbance propagate in opposite directions. It is present only for very low wavenumbers and has an abrupt frequency cut-off. From figure 2, one can see that the DA mode shares the same frequency cut-off; that is, below a frequency that is mostly T e , , U i z dependent neither mode can propagate. Moreover, as the ionization frequency increases the common cut-off frequency rises, and the LW mode can propagate for larger wavenumbers. It appears that both dispersion relations could be joined in a single continuous function ω r (k); in this sense the LW mode could be thought of as a second branch of the DA mode. However, the different physical processes related with each wave and the different behavior of the permittivity in the neighborhood of the zero crossings hinder such a statement.
One should note that for very small wavenumbers the mode eigenfrequencies tend to reach ion characteristic time scales. For the description of the mode in such regimes, k → 0, ion dynamics also have to be taken into account. Here, to ensure operation in the low-frequency regime, the dispersion relations are always truncated in frequencies lower than the dust plasma frequency ω pd .
In order to evaluate the effect of neutrals on both modes, we plot the dispersion relations for T e = 5 eV, a = 1 µm, P = 1 and varying P n from 20 to 60 Pa in figure 3. One notices that for increasing pressure: for the DA mode, the wavenumber cut-off is still increasing and also the mode eigenfrequencies for constant wavenumber are monotonically suppressed, while for the LW mode, there appear larger wavenumbers where the mode can exist and also the steepness of ω r (k) is decreasing.
In figure 4 we investigate the effect of the dust radius on the dispersion relations for T e = 10 eV, P = 1, P n = 20 Pa and varying dust radii from 0.5 to 7.5 µm. The electron temperature was chosen to be sufficiently high in order to ensure that the LW mode exists even for the largest grain radius (e.g. for 5 eV the LW mode disappears already for a = 5 µm). Also note that P is kept constant and not the dust density n d , which decreases with increasing grain radius; such a choice was made to ensure the smallness of . We conclude that an increase of radius leads to lower frequency cut-offs of the LW mode, but slightly affects its behavior for very low wavenumbers, while in the DA mode it greatly suppresses the eigenfrequencies for constant wavenumbers. The latter is to be expected, since an increase of the radius decreases the dust plasma frequency that scales as ω pd = 4π n d Z 2 d e 2 m d ∼ a −1/2 when roughly assuming that z is radius independent in the weakly collisional regime for ions. The dust plasma frequency sets the upper limit in the frequency of electrostatic oscillations the massive grains can respond to and is also the small-wavelength limit of the DA mode, i.e lim kλ Di 1 ω r (k) = ω pd . In figure 5, the dispersion relations are plotted for varying dust density parameter P from 1 to 0.05, P n = 20 Pa, T e = 5 eV and a = 1 µm. The LW mode exists only for P = 1, 0.5; for smaller dust densities it stops propagating and the critical wavenumbers k crit of the DA mode reappear at zero frequencies (as in figure 1). As expected from the hydrodynamic model k crit increases as P decreases, since roughly k crit ∼ P −1/2 . Hence, apart from the condition of large ionization frequency, the LW mode requires a sufficient amount of dust to propagate, since the average plasma capture frequency on dust is proportional to the dust density.
Finally, in figure 6 we investigate the dependence of the dispersion relations on the gas type, which is mainly due to the difference in the first ionization potential but is also due to the different resonant charge transfer cross-sections. The parameters are a = 1 µm, T e = 5 eV, P = 1 and P n = 20 Pa. One can notice that as the ionization potential decreases, the LW mode can propagate for dramatically larger wavenumbers. Already for argon gas the length scale of the oscillation is much lower than 1 cm, which means that the waves can be easily observed in typical laboratory discharges. Moreover, the cut-off frequency also increases dramatically (but is still lower than the dust plasma frequency, ω pd 900 s −1 , for all operating gases). Similar effects can be noted for the DA mode.
We note that even for the largest ionization frequencies considered in the above results, the ionization characteristic length l iz = v T e ν e is much larger than the perturbed radius R 0 around the grain and therefore the charging process is not affected by ionization of neutrals. We should also report that the results are sensitive to the dust charging model (e.g. a change in σ n,i by a factor of two leads to significant deviations).

The cut-off frequencies and critical wavenumbers
In this subsection, we investigate the simultaneous cut-off of both modes by plotting the real part of the permittivity, { k,ω }, for frequencies near the cut-off frequency. As seen in figure 2, for T e = 5 eV, the cut-off frequency is about 90 s −1 . The existence of two modes in the same frequency range means that for any given frequency, { k,ω r }, when viewed as function of k, has two roots. Taking into account that the function is continuous and negative in the k → 0 limit, we conclude that { k,ω r } must have a maximum between the two roots and that it must be positive. Hence, the modes can only be in simultaneous cut-off, when this maximum becomes negative. This is confirmed by figure 7; for ω r = 90 s −1 the maximum is still positive and both modes propagate, while for ω r = 88 s −1 the maximum is negative and both modes are in cut-off. The maximum in { k,ω } is shifted towards lower values as the frequency decreases; therefore the mode spectrum is always continuous. This type of cut-off is qualitatively different from the cut-off of the DA wave which appears in the absence of strong ionization and due to the effect of neutrals (see figure 1). In that case, the DA wave is in cut-off for a critical wavenumber that corresponds to zero frequency and { k,0 } is always positive.

Properties of the LW mode
From the above presented results, the following characteristics of the LW mode emerge.
• It is determined by the interplay between electron impact ionization and plasma absorption.
• It does not exist for: (i) low electron temperatures compared to the gas first ionization energy, (ii) relatively low pressures and (iii) relatively low dust densities.
• Its dispersion relation appears to be adequately fitted by ω r (k) = A 1 + A 2 , where the coefficients A 1 , A 2 , A 3 depend on the discharge, dust grain and plasma parameters, while A 1 determines the cut-off frequency of the wave. Such dispersion relations are typical for ionization waves in discharge plasmas [36].
• It is present for large wavelengths. However, as the ionization frequency/pressure increases it can propagate at larger wavenumbers.  • For very small wavenumbers (k → 0) it reaches time scales characteristic of ion dynamics and the low-frequency kinetic model considered here does not suffice for a self-consistent solution.
• The phase velocities are an order of magnitude larger than the DA phase velocities.
• The group velocity is negative, v ph v g < 0, such backward waves are also typical of certain classes of ionization waves [36].
• The presence of the mode in the large-wavelength part of the Fourier spectrum leads to a cut-off of the DA mode in this range. The cut-off frequency/wavenumber of the DA mode is the same with the ones of the LW mode.

Plausible physical scenarios connected to the LW mode
The LW mode depends on the trade-off between plasma generation and loss, yet appearing almost entirely in the low-frequency regime, the inertia has to be provided by dust motion. A possible physical mechanism leading to this mode could be the following: in the homogeneous complex plasma system, electrons/ions are continuously absorbed on the dust grains while new ones are created in the regions between the grains by electron impact ionization. However, the mobile electrons having less friction in their interactions with dust/neutrals are removed much faster from the ionization region, which creates a positive polarization charge, which is the source of an electrostatic field. Furthermore, the dependence of the ionization frequency on the electron density will lead to a polarization charge that is oscillating in space; thus, under the effect of the produced field, the grains will start oscillating, moving though alternating regions of attraction/repulsion with their neighboring grains. Such a physical process is reminiscent of the mechanism of long-range collective grain attraction [35,[37][38][39][40]. In particular, away from the region of strong dust-ion coupling and the position of the attractive well in the pair interaction potential, it has been shown [41] that for non-fluctuating sources the potential energy is exponentially decaying (in our case the LW mode does not exist for such sources), whereas for proportional sources it exhibits an oscillatory behavior (in our case the LW mode exists). To complete this analogy, also note that the corresponding space scales of both phenomena are the same. However, in [41] the steady-state system is solved; dust is therefore assumed to be immobile and the potential is only spatially varying.
It is important to point out again that all waves are damped due to collisions with neutrals. Naturally, the DA mode having a monotonically increasing ω r (k) dispersion relation as a function of k with the property lim kλ Di 1 ω r (k) ω pd can be completely described in the lowfrequency regime. On the other hand, the LW mode for very low wavenumbers (kλ Di 10 −3 ) does not saturate in frequencies and its eigenfrequencies reach ion time scales of response. For the description of this region, the low-frequency kinetic model used here clearly does not suffice, since in the computation of the ion kinetic responses the frequency is dropped. Hence, an instability cannot be excluded for very large wavelengths. In fact, another argument leading to such a possibility is that a self-consistent hydrodynamic description of the system in the case of dust-ion acoustic waves [42,43] has shown that an instability appears at such length scales (k → 0).

Experimental realizations
As is evident from the numerical results of section 3, the dispersion relations depend crucially on the ionization frequency. Two scenarios emerge: (I) for configurations with low ionization frequency, only the DA mode appears; (II) for configurations with high ionization frequency, both-DA and LW-modes appear. The corresponding wave properties can be summarized as follows.
(I) For relatively low pressures (see figure 1), the deviations from the fully ionized case are very small; therefore the results of [21] are valid, predicting large deviations from the multi-component description in long wavelength ranges for large grains, high plasma densities and relatively high dust densities. For higher pressures (figure 1) the DA waves can propagate at larger wavelengths than the ones the hydrodynamic model predicts. (II) Frequency cut-offs are introduced and the DA mode does not propagate at long wavelengths (figures 2-6), which are the wavelengths where kinetic effects (dust charging, plasma absorption on dust) are most pronounced [21]. The frequency cut-offs are common for both modes (see figure 7).
Hence, we conclude that an experiment for studying dust charging effects on DA waves would benefit from large grain radii, large wavelengths and small ionization frequency. Small ionization frequency can be achieved by operation with a gas like neon in typical dc or rf discharges, or by operating in configurations with low neutral gas pressures, e.g. magnetic cusp devices [44].
On the other hand, an experiment aiming to investigate the LW mode should operate with gases of low ionization potential (e.g. argon and krypton) and would benefit from plasma regimes with elevated electron temperatures and large pressures. Quantitatively speaking, for argon gas, P n = 20 Pa, and for T e = 5 eV, the frequency cut-off appears at k = 20 cm −1 (see figure 6), implying wavelengths below 1 cm, which could be feasible to observe provided they are not damped. Plasmakristall-4 (PK-4) [45,46], scheduled to fly on the International Space Station (ISS) in 2013, is ideal for such regimes: (i) the observation of long wavelengths is viable due to the 35 cm glass tube; (ii) the electron temperature is a slowly decreasing function of gas pressure, with T e ranging from 6 to 7 eV; (iii) when implemented on the ISS, grains with a > 1 µm can levitate in the bulk plasma under microgravity conditions and therefore the radius dependence can be explored; and (iv) experiments are scheduled for large pressures and neon/argon gas, under both laboratory and microgravity conditions.