Particle-in-cell simulation of the effect of dust charge fluctuation on ion acoustic waves in a dusty plasma.

A new hybrid-particle-in-cell (PIC)-Monte Carlo Collision (h-PIC-MCC) algorithm is presented here. The code correctly simulates the damping of ion acoustic wave due to dust charge fluctuation in a dusty plasma along with other kinetic effects such as Landau damping. In the model, on event of a collision between a charged particle and a dust particle, a randomised probability determines whether the charged particle is absorbed by the dust with the collision cross section being determined dynamically by the overall interaction scenario. We find that this method is versatile enough as it can also include the size and mass distribution for the dust particles, in addition to the charged species dynamics. As such, it can be adopted to study numerous phenomena that occur in diverse dusty plasma environments. We believe that the damping of the ion acoustic wave through dust charge fluctuation is being demonstrated, for the first time, with a PIC code, in this work.

Introduction

Dusty plasma provides a broad environment for carrying out extensive research in the field of plasma physics. Examples of dusty plasma can range widely, from laboratory to space plasma environments[1-3]. Such a system is characterised by the presence of dust particles whose radius can vary from nanometers to micrometers in size. The fundamental question that arises in any dusty plasma environment is how these particles charge when immersed in plasma. Thus, the dust charging process has always been a subject of wide interest[4-6]. As a consequence of the dust charging process, numerous phenomena can occur in the plasma. A significant outcome among them is the exhibition of self-consistent dust charge fluctuation that results as a response to the variation in the background plasma currents into the dust particles[5,7-9]. This fluctuation can lead to the damping of the plasma waves as well as development of other instabilities[10-15]. We consider the effect of such dust charge fluctuation on the electrostatic ion acoustic wave as the primary subject of this investigation. Some notable early contributions in this field include the works of Jana et al.[10] and Ma et al.[12], which consider the damping of ion acoustic wave due to dust charge fluctuation, using fluid and kinetic theory, respectively. In this work, we perform a kinetic particle-in-cell (PIC) simulation to study the characteristics of the ion acoustic wave in presence of dust charge fluctuation in a dusty plasma. Following the general outline[16], we adopt a modified PIC algorithm to estimate the collisions of the dusts with the electrons and ions. We thus use our hybrid-PIC-Monte Carlo Collision (h-PIC-MCC) method to estimate the damping of the dust-ion-acoustic (DIA) wave, which we believe, has been carried out for the first time with a PIC simulation. We find that our numerical results corroborate well with the theoretical approximations.

In this work, the primary emphasis is on the effect of dust charge fluctuation on the ion acoustic wave rather than the detailed dust charging process. As we shall indicate in the relevant sections, our numerical scheme does compromise on the details of the dust charging process at the expense of gaining speed in the simulation. We show that despite this shortcoming, the simulation model correctly captures all the important signatures of the effect of dust charge and its dynamic behaviour on ion acoustic wave. Throughout our work, in many places we refer to the dust-ion-acoustic (DIA) wave as ion acoustic wave (IA) as the dust particles remain stationary and dynamically the dust dynamics does not enter into the ion acoustic oscillations. In “The theory of dust-ion-acoustic wave” section, we briefly describe the fluid and kinetic theories of DIA wave. In “Simulation model” section, we describe our simulation model in detail. The simulation results are described in “Simulation results” section. We conclude our work in “Conclusions” section. An “Appendix” outlining the computational details can be found at the end.

The theory of dust-ion-acoustic wave

We first consider the excitation of electrostatic ion acoustic wave in an e-i plasma without either resorting to plasma approximation or Boltzmannian electrons. The model is described by the electron and ion continuity and momentum equations,where for electrons and ions respectively, is the electronic charge, are thermal pressures and other symbols have their usual meanings. We have expressed the temperature in the energy unit. The equations are closed by the Poisson equationThe linear dispersion relation of this fluid model can be written as[17,18],with are the respective plasma frequencies and are the thermal velocities, where is the equilibrium plasma density. The above dispersion relation is being used to benchmark our simulation.

In what follows, we briefly discuss the kinetic theory behind the dust charge fluctuation dynamics and related kinetic effects[9,12]. We consider the 1-D linear ion acoustic wave with a characteristic frequency lying in the range , where is the dust plasma frequency. The kinetic theory of ion acoustic wave with dust charge dynamics is due to the Boltzmann equationwhere are the ion and electron distribution functions, is the electrostatic potential, and is the collision operator which represents the rate of electron and ion capture by the dust particles thus incorporating the dust charge dynamics. In the linear theory, for negatively charged dust particles, which is the case of this present work, the oscillating current to the dust grain is given bywhere is the effective collision cross section of the respective species with the dust particles and is the charge relaxation rate are the dust charges, and C is the capacitance of a dust particle. The lower limit of the integration , respectively for electrons and ions, is the minimum velocity with which the electron or the ion must reach the surface of the dust particle. In the above expression and in the expressions hereafter, the subscript ‘0’ and ‘1’ refer to the equilibrium and perturbed values of the respective quantities. The expression for the collision cross section iswhere is the average dust radius. The equations are closed by the linearised Poisson equationThe linear dispersion relation for ion acoustic wave can be written aswherewhere are the ratios of equilibrium dust density to that of electron and ion densities. In the above expressions, are the dielectric susceptibilities and represent the coupling of dust charge fluctuation with electron and ion density perturbations. Eq. (10) predicts damping of the ion acoustic wave, which is due to the dust charge fluctuation. Besides, also contains the Landau damping term for the wave.

We also note that, even without dust charge fluctuation, the ion acoustic wave is subjected to Landau damping. By neglecting the dust charge fluctuation in the above analysis, the Landau damping rate for ion acoustic wave can be approximated as[9]where is the ratio of equilibrium electron to ion densities, , and

Simulation model

Our primary simulation model comprises of a 1-D electrostatic e-i dusty plasma with periodic boundary conditions. The dynamics is controlled with an electrostatic particle-in-cell (PIC) code. Both electrons and ions are considered as thermal particles, distributed with their respective distributions, while the dust particles are assumed to be cold. The simulation parameters correspond to a typical laboratory situation with plasma density and dust density of . The electron and ion temperatures are respectively and , and an e-i mass ratio of . In a simulation box length of , both electrons and ions are represented with macro-particles, while the dust particles are represented with 2500 macro-particles. We carry out our simulation with 600 equally spaced 1-D cells. The electron Debye length of the plasma is and the evolution time step is . With these simulation parameters, we are able to have full spatial resolution with a temporal resolution of the order of the electron time scale. In Fig. 1 (left panel), we have shown the spectral plot of the ion acoustic dispersion relation in absence of dust particles as determined from the simulation, superimposed with the theoretically calculated curve. The oscillations in Fig. 1 (right panel) show the electron and ion density oscillations. While the solid curve represents oscillations in , the curve with fine-scale oscillations is that in , where the fine-scale oscillations are due to Langmuir (electron plasma) oscillations modulated by the ion acoustic oscillations. It can be seen that the simulation results are in good agreement with the theory. It is also seen from these curves that the ion density is closely followed by electron density thereby validating the use of plasma approximation (not used here though).

Figure 1

(left panel) The spectral plot of the ion acoustic dispersion (without any dust) relation from the simulation. The solid curve is the theoretically calculated one. (right panel) The corresponding ion acoustic wave in electron and ion densities . The solid line is the ion density and the other one is electron density fluctuation. (Generated using Gnuplot, Version 5.2 patchlevel 7, URL : http://www.gnuplot.info).

A few clarifications about our simulation model is in order here. As mentioned at the beginning, the goal of this work is to demonstrate the proof of principle about the effect of dust charge fluctuation on ion acoustic waves in an e-i plasma, implemented in a PIC code. Toward this, we have used the simplest possible plasma model, namely a fully ionised e-i hydrogen plasma with stationary dust particles. Various plasma parameters are such that which can be readily implemented in a laboratory[19-22]. We note that though a hydrogen plasma is not the very usual choice for laboratory plasmas, there are laboratory experiments[19,20] which use hydrogen plasmas. Computationally, as we work in the ion acoustic domain, simulating an ion acoustic wave with heavier ions such as in argon or potassium plasmas and hence longer wavelengths, becomes more challenging. Nevertheless, our computational model, in principle, can be extended to any such laboratory plasma situations. A brief description of the computational model is given in the “Appendix”. In our model, we have neglected the neutrals in the plasma. The presence of neutrals in an e-i plasma affects the plasma dynamics through neutral collisions. However, as we are not considering any neutral collision here and focus our attention mainly on electrostatic interactions, the neutrals can be neglected as far as the ion acoustic wave is concerned.

Dust charging with PIC simulation

We now describe our dust-PIC model, which is used to study the dust charging mechanism and its effect in a dusty plasma. It is worth noting that the classical PIC model does not account for collisions[16]. Besides every particle in a PIC model is actually a macro-particle comprising of a number of real-life particles. As such, any collision implemented under the PIC formalism will actually account for collisions en masse. However, there are ways to include collisions among participating particles in PIC, one of which is the popular PIC-Monte Carlo Collision (PIC-MCC) algorithm[23], usually used for neutral collisions. It consists of using a randomised probability to account for the collisions based on the theoretical estimation of the collision cross sections. Multistep Monte Carlo Collision is another algorithm adopted by Gatsonis et al.[24] which is based on an algorithm similar to the above but considers dust collisions with the charged species. We however use here a different method to estimate for the collisions of dust particles with electrons and ions, which can be described as the hybrid-PIC-MCC (h-PIC-MCC) method.

The hybrid-PIC-MCC scheme

Our hybrid-PIC-MCC (h-PIC-MCC) scheme is based on the Monte Carlo (MC) scheme described by Birdsall[23]. Consider the collision frequency between a dust and a charged particle ‘s’ (either an electron or an ion) ,where is the local density of dust particles in the particular computational cell. The absorption cross section is in the event of a collision between a dust particle and a charged particle ‘s’ and is the relative velocity between the charged particle and the dust particle. However, as in our case the dust particles are stationary, is the velocity of the incoming particle. The absorption cross section is given by,where is the potential difference between the dust particle and the local plasma and is the radius of the dust particle, and and are the charge and mass of the colliding charged particle. So, the probability of collision of a charged particle ‘s’ within a time interval is given by[23],Note that the above relation for a single species of physical particle can be equivalently written for a computational particleas with sufficient number of collisions, the resultant velocity distribution will resemble that of a single physical particle[23]. We now introduce a uniform random number [23,24] and assume that if , an absorption occurs, otherwise not. It should be noted that this scheme allows only one collision per time step for a charged particle where either the charged particle is absorbed or remains in the bulk plasma, unlike the multistep MCC process[24]. In a multistep MCC process, the dust particle is treated as the incoming particle which collides with a target particle—a charged particle. Once there is a collision, it checks for more collisions within the time interval over which the system evolves. In both the cases, when an absorption occurs, the computational particle is removed from the bulk plasma and an equivalent charge number is updated for the dust particle accordingly.

However unlike typical PIC-MCC process[23], we estimate through the Coulomb interactions at the grid, which accumulates the charges at the grid points from actual positions of the charged particles through an interpolation process (see “Appendix”). So, instead of estimating the actual absorption cross section for a single dust particle, the process estimates a combined cross section for all the dust particles which are there in a particular cell. We then determine which of the dust particles in a cell would actually absorb the charged particle, depending on the final position of the charged particle and position of the dust particle in that cell. This process is repeated for all the computational particles in the system. This way, we can considerably speed up the calculation of cross section in each time step at the cost of loosing certain charging information for each dust particle which are there in one cell. However, on average the process correctly estimates the dust charge accumulation averaged over a computational cell. In order to account for the dust charge fluctuation in our simulation, we assume that whenever a collision of the dust particle with an electron occurs, it contributes to an increase of negative charge on the surface of the dust particle. This is effected by decreasing the number of electrons in the simulation domain accompanied by an equivalent increase of electron dust charge number . The accumulation of positive charge on the surface of a dust particle is, however, modeled by assuming that whenever an ion-collision occurs, an electron is ejected from the dust particle which causes to decrease (equivalently charging the dust particle positively) and an increase of a plasma electron in the simulation box[8]. Hence, the total ion number in the simulation domain remains constant while the total electron number and dust charge fluctuate depending on the type of collisions.

We note that all the above simplifications in our h-PIC-MCC makes it impossible for this scheme to exactly determine the charging history of individual dust particles while capturing the charging and charge fluctuation dynamics averaged over a computational cell in an ion acoustic time scale. As a result, we can expect the effect of dust charge fluctuation on ion acoustic wave (or dust-ion-acoustic) to be correctly estimated while compromising on the charge history of individual dust particles. This can be considered as a bargain as our emphasis is primarily on the effect of dust charge fluctuation on ion acoustic wave, rather than the charging dynamics of the dust particles.

Dust particle versus dust site

As the dust particles are considered stationary, the only way they contribute to the dynamics of the e-i plasma is through electrostatic interaction. As such the mass of the dust particle as well as their size do not involve in the dynamics of the plasma except for estimating the absorption cross section, which is measured in terms of . Besides, as far as the electrostatic interaction of the plasma particles with dust particles are concerned, each dust particle behaves like a collection of electrons (as they are negatively charged) with zero temperature (hence zero velocity). Which is why, we would like to refer them as dust sites rather than dust particles. So a dust site with number of charges on it behaves as number of electrons, occupying the same position in the phase space. In the PIC framework, this is perfectly admissible, as each computational particle can be considered as point particle (compared to the scale length of the plasma) and they can overlap each other[25].

Simulation results

Dust charging

As mentioned in the previous section, we consider the dust component of the plasma to be cold and stationary. The plasma quasineutrality condition[9] is maintained at all time,where are the instantaneous physical densities of the ions, electrons, and dust particles. Computationally, the above relation can be represented aswhere is the instantaneous dust charge number averaged over all the computational particles. The quantities for each species are called weights, where L is the simulation domain length, is the initial physical density for the species ‘s’ (at the beginning of the simulation), and is the total number of computational particles for the species ‘s’ at any instant of time t. The weight is a dimensionless quantity which measures number of physical (actual) particles of a species that are being represented in a single macro-particle (or computational particle) of the species. So, is the total number of actual physical particles present in the simulation domain at a time t. Note that as per our scheme, and are the number of dust sites (or macro-particles) which are fixed in space and time. We use 2500 dust sites, distributed uniformly in the domain with zero net charge. For calculation of the absorption cross section , the dust radius is fixed at . The number of dust sites are chosen so as to have no overlapping of the sites as they occupy the simulation domain. The dust sites can be considered to be equivalent to dust macro-particles, each having radius . We note here that each computational particle (or macro-particle) consists of a large number of physical particles determined by the weight factor w and each computational particle has its shape function (see the “Appendix”), which in our case is . Or in other words, the computational particles can still be considered as point particles, so long as the scale length of interactions among themselves is much larger than their size. In our case, this is easily satisfied as our interaction scale length is and . Physically, each dust site behaves as a collection of stationary electrons (assuming the dust particles to be negatively charged). We can calculate the equivalent dust charge number by considering the charging current for the dust sites. An orbit motion limited (OML) estimate of the charging current in a dusty plasma is given by[1]where is the charging time and is the average dust charge. Computationally we can only find out the depletion of electron macro-particles in each time step as they get absorbed in the dust sites, which provides the total current to the dust particles over the entire computational domain. As an example, assume that number of electron macro-particles are depleted in time . The total charge accumulated in the dust sites over the entire computational domain is given byso that . The charging current for an individual dust particle can be written asfrom which, using OML approximation as given by Eq. (23), we can find out the average dust charge number asWith the plasma parameters in simulation as mentioned in the beginning of the “Simulation model” section, we get a dust-charging current at saturation of the dust charge with an average dust charge . In Fig. 2, the average charging of single dust particle is shown, calculated from the total count of electron depletion in the simulation domain. As can be seen from the figure, on average, a single dust particle attains an equilibrium net charge of about . A theoretical estimate for can be written as[1]which is only about half the value calculated from the simulation and can be considered to be in good agreement. According to the OML approximation[9,26] also, the initial charging current can be written aswhich is of the same order of magnitude as the simulation value. Note that initially the electron and ion plasma density is equal . In the above expression, the factor is the potential difference between a dust grain and the bulk plasma, which is initially. For the parameters used in the simulation, a theoretical estimate of the initial charging time (till a dust particle attains an average equilibrium charge)[1]which in our simulation is calculated as .

Figure 2

The charging of a single dust particle. is the average dust charge number for a single dust particle. The time is in the unit of ion acoustic time scale . (Generated using Gnuplot, Version 5.2 patchlevel 7, URL : http://www.gnuplot.info).

Constant dust charge

We now present the simulation results for the effect of constant dust charge on ion acoustic wave (or dust-ion-acoustic wave). It is well known that the presence of negatively charged dust particles causes the phase velocity of the ion acoustic wave to increase which is also associated with the Landau damping of the wave becoming less severe due to the dust particles. A linear dispersion relation of the dust-ion-acoustic (DIA) wave can be derived from the fluid theory as[9,27]where is the dust plasma frequency. For stationary dust particles, this relation reduces towhich can be further simplified for asFor , with our given parameters, while Eq. (30) predicts an increase of the phase velocity by an increase of about as compared to the phase velocity when there are no dust particles for , for Eq. (31) the value is . For varying k, we need to solve the ion acoustic dispersion relation numerically as given by Eq. (10). Our choice of k corresponds to mode 6 of perturbation in the periodic domain, with .

In Fig. 3, we show this effect as simulated through our dust charging model. Both the panels in the figure show the theoretically calculated -k curves from solutions of Eq. (10) for (blue colour, no dust) and 0.5 (red colour, fifty percent depletion of electrons on dust particles) without any dust charge fluctuation. Superimposed in both the panels, are the reconstructed -k spectra through the simulation—for (left panel) and (right panel). As we can see, both the simulated spectra closely follow the theoretical curves. The associated effect corresponding to this increase of phase velocity is the overall decrease of severity of Landau damping of the ion acoustic wave[9,28-30]. There are several experimental works[28,29] which investigated this effect in the past. We note that though the parameter regimes over which these experimental works are carried out favours decrease of the Landau damping rate , an actual numerical evaluation of the same via Eq. (16) may also show increase of the damping rate[31], depending on the values of and k. We also note here that in absence of any dust dynamics, the effect of dust charge on Landau damping of ion acoustic wave, as per Eq. (16), is limited only to the value of . In left panel of Fig. 4, we show the estimated Landau damping rate from simulation as one varies from 0 to 1. The solid line in the figure represents the theoretical estimate. We find that the simulation results on average mimic the theoretical values. The whole computation takes about 1 h in a 4-core CPU machine for .

Figure 3

The simulated -k spectra shown alongside the theoretical curves. The solid curves are same in both the panels, representing (blue) and 0.5 (red). The corresponding spectra are shown in the left and right panels, respectively. (Generated using Gnuplot, Version 5.2 patchlevel 7, URL : http://www.gnuplot.info).

Figure 4

Linear Landau damping rate for an e-i plasma with constant dust charge. (Generated using Gnuplot, Version 5.2 patchlevel 7, URL : http://www.gnuplot.info).

Dust charge fluctuation

The dust charge fluctuates due to varying electron and ion currents to the dust particles. This charge fluctuation primarily leads to a noticeable damping of the ion acoustic wave[9,10,12]. As discussed above, even without dust charge fluctuation, the ion acoustic wave suffers Landau damping, which is shown in Fig. 4. For the plasma parameters used in this simulation and for , the Landau damping rate (see Fig. 4). The additional linear damping due to dust charge fluctuation can be calculated numerically from the full dispersion relation Eq. (10), which for the simulation parameters, comes out to be . This effect is very clearly visible as the damping rate of the ion acoustic wave in presence of dust charge fluctuation is almost one order of magnitude larger than the damping when the fluctuation was absent for the same dust charge .

In the left panel of Fig. 5, the dots represent the average amplitude of the perturbation and the solid line is the best fit line of the linear damping rate , which is estimated to be about , which is quite close to the theoretically estimated value. The panel on the right in Fig. 5 shows the actual simulation data, superimposed with the theoretically calculated damping oscillations (solid red colour line). In the simulation, we find out the damping rate from the electrostatic oscillations, measured in terms of the square root of potential energy .

Figure 5

The damping of the ion acoustic wave due to dust charge fluctuation. The panel on the left shows the amplitude of the ion acoustic wave with time fitted with the best fit damping rate . The panel on the right shows the actual simulation data, superimposed with the theoretically calculated damping curve with . (Generated using Gnuplot, Version 5.2 patchlevel 7, URL : http://www.gnuplot.info).

Conclusions

To conclude, we have investigated the effect of dust charge fluctuation on ion acoustic wave using a 1-D hybrid-PIC-MCC algorithm. Our model is able to predict the damping of ion acoustic wave due to dust charge fluctuation, with the damping rate agreeing quite well with the corresponding theoretical estimate, which is demonstrated for the first time in this work, using a PIC code. This damping is a result of the coupling between the dust charge fluctuation and the ion acoustic wave through the plasma density variations. Our numerical model could reproduce all the benchmark results for ion acoustic wave in presence of negatively charged dust particles. These results include increase of phase velocity of ion acoustic wave in presence of negatively charged dust particles as well as decrease in Landau damping rate of the wave as dust charge increases. All these results agree quite well with the corresponding theoretical estimates. On a future outlook, this method can very well be extended for two and three dimensional geometry.

Nevertheless, we would like to emphasise that our numerical scheme does not represent the dust charging dynamics in details as far a single dust particle is concerned but it does result compatible effects of the dust charging dynamics averaged over the computational cell. This compromise leads to a considerable speed up in the calculation while capturing the effect of the dust charge dynamics on ion acoustic wave. Finally, we would also like to mention that we have used this particular numerical concept as proof of principle to demonstrate its efficacy. Though we have used a number of simplified assumptions in respect of dust charging dynamics, these assumptions can very well be relaxed and the scheme should still work. For example the constancy of ion density during dust charging, can be relaxed and the same method which we have adopted for electron-dust collision, can also be applied for ion-dust collision. We have taken care to see that our simulation parameters represent a typical laboratory dusty plasma, which can easily be used to study other dusty plasma environments such as space and astrophysical plasma as well.

The CIC charge interpolation scheme. (Generated using PGF/TikZ, Version 3.1.4a, URL : https://github.com/pgf-tikz/pgf).