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.


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 j = e, i for electrons and ions respectively, q i,e = ±e is the electronic charge, p j = γ j n j T j 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 equation The linear dispersion relation of this fluid model can be written as 17,18 , with ω pj = n 0 e 2 /(ǫ 0 m j ) are the respective plasma frequencies and v thj = T j /m j are the thermal velocities, where n 0 = n i0 ≃ n e0 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 ω pd ≪ ω ≪ kv the , where ω pd is the dust plasma frequency. The kinetic theory of ion acoustic wave with dust charge dynamics is due to the Boltzmann equation where f j ≡ f i,e are the ion and electron distribution functions, φ is the electrostatic potential, and F coll 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 by where σ j is the effective collision cross section of the respective species with the dust particles and η is the charge relaxation rate Q d0,d1 are the dust charges, and C is the capacitance of a dust particle. The lower limit of the integration v mj = v m , 0 , 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 is where r d is the average dust radius. The equations are closed by the linearised Poisson equation The linear dispersion relation for ion acoustic wave can be written as where (1) ∂n j ∂t + ∇.(n j v j ) = 0, (2) m j n j dv j dt =q j n j E − ∇p j , (3) ǫ 0 ∇ · E = e(n i − n e ).
pe ω 2 − γ e k 2 v 2 the = 0, ǫ 0 ∇ 2 φ 1 = en e1 − en i1 − Q d0 n d1 − Q d1 n d0 . where δ e,i are the ratios of equilibrium dust density to that of electron and ion densities. In the above expressions, χ e,i are the dielectric susceptibilities and η e,i 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, η i 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 γ L for ion acoustic wave can be approximated as 9 where δ n = n e0 /n i0 is the ratio of equilibrium electron to ion densities, c s = ω pi De , 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 n 0 ∼ 10 16 m −3 and dust density of ∼ 10 11 m −3 . The electron and ion temperatures are respectively T e ∼ 1 eV and T i ∼ 0.01 eV , and an e-i mass ratio of (m e /m i ) −1 ∼ 1835. 16 . In a simulation box length of 0.004 m , both electrons and ions are represented with 10 5 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 De ∼ 7.4 × 10 −5 m and the evolution time step is ∼ 10 −11 s . 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 (n i,e ) oscillations. While the solid curve represents oscillations in n i , the curve with fine-scale oscillations is that in n e , 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 n i ∼ n e (not used here though).
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][20][21][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 (12) 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) ν s , where n d is the local density of dust particles in the particular computational cell. The absorption cross section is σ ab in the event of a collision between a dust particle and a charged particle 's' and v s is the relative velocity between the charged particle and the dust particle. However, as in our case the dust particles are stationary, v s 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 r d is the radius of the dust particle, and m s and q s are the charge and mass of the colliding charged particle. So, the probability of collision of a charged particle 's' within a time interval ∆t is given by 23 , Note that the above relation for a single species of physical particle can be equivalently written for a computational particle as 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 0 ≤ R ≤ 1 23,24 and assume that if P p > R , 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 ν p 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 (18)  www.nature.com/scientificreports/ 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 z d . 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 z d 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 r d . 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 z d number of charges on it behaves as z d 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 n i,e,d are the instantaneous physical densities of the ions, electrons, and dust particles. Computationally, the above relation can be represented as where z d (t) is the instantaneous dust charge number averaged over all the N d computational particles. The quantities w s = L(n s0 /N s ) for each species are called weights, where L is the simulation domain length, n s0 is the initial physical density for the species 's' (at the beginning of the simulation), and N s (t) is the total number of computational particles for the species 's' at any instant of time t. The weight w s 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, N s (t)w s is the total number of actual physical particles present in the simulation domain at a time t. Note that as per our scheme, N i (t) ≡ N i = const. and N d 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 σ ab , the dust radius is fixed at r d ∼ 10 −2 De . 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 r d . 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 S x (see the "Appendix"), which in our case is S x (x − x p ) = δ(x − x p ) . 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 ∼ De and r d ≪ De . 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 z d 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 Q d = ez d 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 I tot to the dust particles over the entire computational domain. As an example, assume that ∆N e (22) n i = n e + z d n d , With the plasma parameters in simulation as mentioned in the beginning of the "Simulation model" section, we get a dust-charging current I d ∼ −9.33 × 10 −3 µA at saturation of the dust charge with an average dust charge z d ∼ 670 . 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 Q d ∼ 670e . A theoretical estimate for I d 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 as which is of the same order of magnitude as the simulation value. Note that initially the electron and ion plasma density is equal n e0 ∼ n i0 ≡ n 0 . In the above expression, the factor ∆φ is the potential difference between a dust grain and the bulk plasma, which is ∼ 0 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 ∼ 0.01 µs.

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 (ω/k) 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-ionacoustic (DIA) wave can be derived from the fluid theory as 9,27  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 δ n = 1 (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 δ n = 1 (left panel) and δ n = 0.5 (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 γ L , an actual numerical evaluation of the same via Eq. (16) may also show increase of the damping rate 31 , depending on the values of δ n 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 δ n . In left panel of Fig. 4, we show the estimated Landau damping rate γ L from simulation as one varies δ n 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 t ≃ 20τ ia .
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 1 > δ n > 0.4 , the Landau damping rate γ L ≃ 6 × 10 5 s −1 (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 γ D ∼ 4 × 10 6 s −1 . 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 z d .
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 e −γ D t , which is estimated to be about γ D ∼ 2 × 10 6 s −1 , 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 |eφ| 1/2 .

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 n i 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  www.nature.com/scientificreports/ a typical laboratory dusty plasma, which can easily be used to study other dusty plasma environments such as space and astrophysical plasma as well.

Appendix: Computational details
For the sake of clarity and completeness, in this "Appendix", we very briefly describe our hybrid-particle-in-cell (PIC)-Monte Carlo Collision (h-PIC-MCC) electrostatic simulation model, that we have used to investigate the physical problem. For a detailed description of the PIC model, the reader can refer to some excellent monographs and review articles available elsewhere 25,32 . The following description is also detailed in one of our earlier papers 33 .
It is well known that the PIC algorithm uses computational macro-particles (also known as super-particles) instead of real particles, so that the distribution function f s for each species is a superposition of several computational macro-particles 33-35 , where the superscript p denotes the attribute of the computational particles. The particle distribution function f p can be defined as, where S x,v are the shape functions of the computational particles and N p is the number of computational particles. In our case , so that all the physical particles within a computational particle have the same velocity. Besides S x (x − x p ) = δ(x − x p ) , which means that they are point particles. It should be noted that the shape functions must satisfy the unitarity condition 25,34 where ξ stands for any phase space coordinate. As all the dynamics are treated in terms of the computational particles, the Vlasov equation for species 's' in 1-D can be written as, The equations of motions for the macro-particles can be obtained by taking the zeroth order, first order moments with respect to x, and first order moment with respect to v in the phase space of Eq. where φ is the electrostatic potential. As mentioned in "Simulation results" section, the quasi-neutrality of the plasma is maintained by demanding (assuming negatively charged dust particles) at every time step, where z d (t) is the instantaneous dust charge number averaged over all the N d computational particles. The weights w s = L(n s0 /N s ) for each species, where L is the simulation domain length, n s0 is the initial physical density for the species 's' (at the beginning of the simulation), and N s (t) is the total number of computational particle for the species 's' at any instant of time t. As mentioned before, the weight w s is a dimensionless quantity which www.nature.com/scientificreports/ measures number of physical (actual) particles of a species that are being represented in a single macro-particle of the species. So, N s (t)w s is the total number of actual physical particles present in the simulation domain at a time t. Note that as per our scheme, N i (t) ≡ N i = const. (see "Dust charging with PIC simulation" section). It is also worth nothing that as this is a 1-D simulation domain, the actual volume of the simulation box is L which can be thought to be equivalent to a three dimensional box of length L with a unit breadth and width-L ≡ L × 1 × 1.

Charge interpolation.
We adopt the cloud-in-cell (CIC) linear weighting scheme in order to interpolate the charge from the particle positions on to the grid points 23,33 . A schematic representation of the scheme is shown in Fig. 6, where the 1-D grid points are marked as j − 1, j, j + 1, j + 2 and the macro-particle with charge q is between the jth and (j + 1) th grid points, marked by the arrow. The assignment of charge onto the jth and (j + 1) th grid points are simply the fraction of charge q weighted by fraction of the length from the particle to the (j + 1) th and jth grid points, where l is the length of the cell and h is the distance of the charged particle from the jth grid point. Fig. 7. The cycle begins by distributing the macro-particles with their charge and respective velocities. This is also the step where any kind of perturbation is launched. In our case, we have given a velocity perturbation to the ions which is of the form a sinusoidal wave with amplitude ∼ v thi , which ensures the linearity of the perturbation. The accumulation of charges from each cell is used to calculate the charge density ρ j from the neighbouring particle positions through a standard interpolation function F CIC as described above, where the summation extends over the particles in the neighbouring cells and total number of different species. The Poisson equation then can be discretized as, where the subscript 'j' of the field quantities denotes the value at the middle of the computational cell. Apparently, the cell vertices are at x j±1/2 . The Poisson equation is solved with the help of a successive-over-relaxation method, which provides the electric field, defined as, Figure 6. The CIC charge interpolation scheme. (Generated using PGF/TikZ, Version 3.1.4a, URL : https :// githu b.com/pgf-tikz/pgf).

Discretisation. The full computational cycle is schematically shown in
Inititial Distribution Absorption / Emission at Dust Sites which provides the new positions and velocities for the macro-particles. However before the next computational cycle begins, the absorption (or emission) of the electrons at the dust sites is determined for each maro-particle depending upon their closest approach to the dust sites (as already described).
Benchmarking. The benchmarking results of the code can be seen in Fig. 8, where the electrostatic potential, kinetic, and total energy, evaluated as the average energy (see below) per particle, plotted against time in absence of any external perturbation. All the initial conditions are same as described before. The left panel of Fig. 8 shows the energy history of a pure e-i plasma without any dust, while the right panel for dusty plasma where 50% of the electrons are depleted in the dust particles. The quantity δ n is the ratio of equilibrium electron density to the equilibrium ion density, δ n = n e0 /n i0 . The dust case is without any dust charge fluctuation. The energies plotted in the figure are the average energies of a particle, where the summation is over all the computational particles (macro-particles) denoted by N j , j = e, i, d for ions, electrons, and dust particles. The quantity Q j , M j are the respective charge and mass of the macro-particles. The <> denotes the average over all the computational particles N = j=e,i,d N j . As one can see that the energy fluctuations which are measured in terms of electron thermal energy can very well be termed as random fluctuations, as desired.