Dust charging in dynamic ion wakes

Complex (or dusty) plasma systems are a special type of low-temperature plasma, where micrometer-sized particles (or dust) can levitate and self-organize into strongly coupled fluids, field-aligned chains, and crystalline solids.^1–5 Such structures are commonly formed in streaming plasmas, where the ion and electron particles flow along the direction of an externally applied electric field. When dust grains are immersed in such an environment, they become negatively charged, which makes them repel each other. However, the dynamical interaction of the dust particles with the ion flow often leads to attractive and non-reciprocal forces, which affect structure formation in those systems.^6–8 Therefore, the study of self-organization, dynamics, and stability of dusty plasmas requires a proper understanding of the dust interaction with a streaming ion flow.

As the ions pass near a negatively charged dust particle, depending on their velocity and proximity, their trajectories can be deflected, resulting in the formation of a wakefield structure downstream of the grain. One theoretical approach for modeling the highly anisotropic dust potential is linear response theory, where the effects of streaming ions on the dielectric response of the plasma are accounted for by assuming that the dust grains interact via a dynamically screened Coulomb potential. This model has been employed to demonstrate an oscillatory wake potential structure behind the dust grain^9–11 [Fig. 1(a)]. This oscillating wakefield potential has been confirmed by Particle-In-Cell (PIC) simulations, which have been employed to determine the structure of the ion wakefield downstream of a dust grain and calculate the resulting nonlinear grain-grain interactions.^12–15 Collisions between ions and neutrals tend to damp out the oscillations, however, leaving only one primary peak in the potential downstream of a grain.

A simplification of this wake structure is to represent the ion wakefield as a point-like region of positive space charge (called the wakefield focus) located downstream of the grain [Fig. 1(b)].^16–18 In this model, the wakefield focus can provide an attractive force for negatively charged dust particles located downstream, yielding field-aligned structure formation. Although the point-charge representation has been successful in demonstrating the non-reciprocal character of the grain–grain interaction, its application is often limited to strongly coupled dusty plasma structures where the downstream position and magnitude of the point charge can be assumed to be constant.

The features of the ion wakefield behind charged dust grains have also been explored using molecular dynamics (MD) simulations of the ions within the plasma flow.^19–21 The efficiency of this numerical technique was recently enhanced by the development of high-performance graphics processing units, capable of parallel processing with more than 1000 processors. Piel¹⁹ introduced a molecular asymmetric dynamics (MAD) code, where ions move in the “naked” Coulomb potential of dust grains, while interacting with each other through Yukawa force accounting for electron shielding. The MAD code employs superions with the same charge-to-mass ratio (and therefore dynamical behavior) as a single ion, which improves time efficiency and allows for large simulation volumes. Despite these approximations, the MAD technique can reproduce important features of the ion wake potential in reasonable agreement with previous PIC simulations.

The numerical techniques discussed so far commonly treat dust grains as static obstacles within the streaming plasma. However, investigation of phase transitions and fluid phenomena in dusty plasmas requires a model where the dust grains are both free to move and allow their charge to vary in the field of streaming ions. In this case, the ion wakes are modified by the presence of the downstream dust particles and the charging of the downstream grains is, in turn, affected by the location of the ion wakes. Therefore, the ion wake formation and grain charging are coupled and should be treated in a self-consistent manner. This paper introduces a molecular dynamics (MD) simulation, which resolves the motion of both dust grains and ions and simultaneously allows for charging of the dust particles by collisions with ions in a flowing plasma.

This paper is organized as follows. The details of the numerical model are presented in Sec. II, including the treatment of the ions, the dynamics of the dust, and the calculation of dust charging. The use of the numerical code is illustrated by modeling the behavior of two vertically aligned dust grains confined within a glass box placed on the lower electrode of a Gaseous Electronics Conference (GEC) rf reference cell. The lower particle is perturbed using a laser, causing it to oscillate within the wake of the upper particle. The results of these simulations are presented in Sec. III. Section IV presents a discussion of the results, with conclusions in Sec. V.

The numerical model DRIAD (Dynamic Response of Ions And Dust) is a molecular dynamics simulation designed to resolve the motion of both the ions and the dust grains on their individual time scales, while allowing the dust charge to vary in response to the changing ion density in the ion wakefield. Following the method described by Piel to develop the MAD (molecular asymmetric dynamics) code,¹⁹ the forces between pairs of ions and among ions and dust particles are treated in an asymmetric manner. The electric field $E→=−∇Φ$ is assumed to arise from contributions from charged dust, ions, and thermal electrons (with temperature $Te$⁠), which are treated as a fluid governed by the Boltzmann factor $exp(eΦ/kBTe)$⁠. Two limiting cases are identified. Far from the charged dust grains, $eΦ≪kB Te$⁠, and a linear approximation of the Poisson–Boltzmann equation yields a solution for shielded (Debye or Yukawa) potentials for the ions. Very close to a dust grain, $eΦ≫kB Te$⁠, and the small electron density allows them to be neglected, resulting in Coulomb potentials for the dust and ions. To address the continuum of transitions between these two cases, the ion–ion interactions are treated as Yukawa interactions, while the force on the ions from the dust arises from the dust Coulomb potential. This allows the nonlinear shielding of the dust grains by the ions to be properly addressed at the expense of an approximate treatment of the electron shielding. However, the force acting on the dust particles is calculated from the shielded ion charges. Comparisons with PIC simulations¹² showed that this asymmetric treatment of the dust-ion forces provides a reasonable agreement for the equilibrium potential distribution and interparticle forces.

Ions will typically reach an equilibrium distribution within one ion plasma period,

where $mi$ is the mass of an ion, $ni$ is the number density of the ions, and $e$ is the elementary electric charge. For the range of plasma densities characteristic of the sheath of a rf discharge, the plasma period is typically $τi≈1 μ$s. Similarly, the dust response time is on the order of the dust oscillation period, $τd=2π/ωd$⁠, which for typical experimental conditions is on the order of 100 ms. Fully resolving particle motion in numerical simulations usually requires that the time steps $Δt$ used to advance the simulation are on the order of $Δt=τ/100$⁠. As the ions in the simulation are accelerated in the vicinity of a dust grain, the DRIAD code resolves the motion of the streaming collisional ions on even smaller time scales, typically $Δti≤τi1000$ s, with the simulation advancing the ions for $NΔti=$ 100–200 ion time steps (the exact number used is determined by the plasma conditions and run-time). The ions are then frozen, and the dust is advanced one dust time step, $Δtd=τd1000=10−4$ s, with the appropriate parameters associated with the ions averaged over the elapsed ion time steps. Depending on the parameters chosen, the ion density and dust charge reach equilibrium within 10–15 dust time steps, during which time a dust particle has moved a distance less than the particle radius.

In the sheath region of a rf discharge, the electron density and electron flux are complex, in part due to the fact that they are dependent on the phase within the rf period. Simulations, which resolve the electron as well as the ion motion on the relevant spatiotemporal scales, have shown that the electrons can even develop convective patterns within a stationary dust cluster.²² The present model does not resolve the motion of individual electrons, instead focusing on the relevant time and spatial scales for the charging and dynamics of the dust grains. For the purposes of calculating the ion–ion interaction force and dust charging, in this work, the electrons are assumed to have a Boltzmann distribution with a background density $ne0=ni0$⁠.

Given the high number density of ions in the plasma, the numerical model computes the trajectories of superions, each of which represents a cloud of ions with the same charge-to-mass ratio (and hence the equation of motion) as a single ion. Allowing ∼100–200 ions per superion results in a reasonable CPU time (typically 12–24 computational hours to resolve one second of dust motion). The equation of motion for an ion with charge $qi$ is given by

where the electric field $E→i$ consists of contributions from the other ions in the simulation, the charged dust particles, and any electric fields present in the plasma (for example, the electric field in the sheath of a rf discharge) and $F→in$ is the ion-neutral collision force.

The ion–ion interactions are calculated assuming a Yukawa potential

where $qj$ is the charge on ion $j$⁠, $rij$ is the distance between ion $i$ and ion $j$⁠, and the Boltzmann electrons provide the screening with shielding length $λDe2=ϵ0kBTe/nee2$⁠. Here, $kB$ is the Boltzmann constant, $Te$ is the temperature of the electrons, and $ne$ is the number density of the electrons. As the region surrounding the negatively charged dust is depleted in electrons, the force acting on the ions from the dust is calculated using a Coulomb potential,

where $rid$ is the distance between the ion and dust grain with charge $Qd$⁠.

As we are particularly interested in the charging and dynamics of multiple dust particles within a flowing plasma, we define a cylindrical simulation region with the cylinder's axis aligned with the ion flow. Ions, which leave the simulation region or are absorbed by a dust grain, are reinjected on the boundary in a manner that is consistent with a shifted Maxwellian distribution with drift velocity $vdr$ and number density $ni0$⁠. The algorithm to determine the velocity and position of the incoming ions is adapted from the algorithm presented in Ref. 23 for insertion of ions on a spherical boundary.

The ions in the cylindrical region are subject to a confinement force from the assumed infinite homogeneous distribution of ions outside the simulation region. In the study by Piel,¹⁹ an analytical expression for the electric field inside of a spherical cavity in homogeneous Yukawa matter²⁴ was used to provide this boundary condition. A closed-form analytical expression does not exist for the electric field inside a cylindrical cavity. Instead, the electric field from these ions is determined by first numerically calculating the Yukawa potential of a homogenous distribution of ions of density $n0$ within the cylindrical simulation region. This potential is then subtracted from a constant uniform background potential, yielding the potential within a cylindrical cavity inside the homogeneous Yukawa material. The confining electric field is calculated from the negative gradient of this potential. This is done once at the beginning of the simulation on a sufficiently fine grid, and the confinement force acting on each ion at a given location is interpolated from the electric field calculated on the grid points.

At high pressures, ion–neutral collisions significantly reduce particle charge.^25–27 Ion–neutral collisions also create a drag force that balances the acceleration due to external electric fields (such as the electric field in the sheath of an rf plasma or the constant electric field in a DC plasma), leading to a constant drift velocity $vd$ in the direction of the applied field. In the present model, the ion-neutral collisions are calculated using the null-collision method.²⁸ The data for the Ar-Ar⁺ collision cross sections are taken from the Phelps database (hosted by the LxCat project).²⁹ 

As ions are accelerated as they approach the charged dust grains, the time step for ions close to a dust grain is refined by a factor of up to 256 to ensure that the distance an ion travels during a time step is smaller than the dust radius. This allows the ion orbits around the dust particles to be resolved and ensures that ion-dust collisions are detected. The trade-off is that each superion in the model, which has the same dynamics as a single ion with the same charge-to-mass ratio, represents a cloud of ions. A superion can orbit close enough to a dust grain that parts of the ion cloud intersect the grain surface, even though the center of the cloud, the location of the superion, does not. For computational expediency, the model treats each superion as a point particle and maintains the number of superions in the simulation, as well as the charge per superion, constant. This is one of the sources of extra noise in the charging of the dust grains.

The equation of motion for a dust grain with mass $md$ and charge $Qd$ in a typical laboratory experiment is given by

where $Fdd$ is the force between dust grains, $Fid$ is the force exerted by the ions on a dust grain, $mdg$ is the gravitational force, $E$ is the confining electric field within the region, $β$ is the neutral drag coefficient, and $ζrt$ is a thermal bath, with $rt$ being a random number uniformly distributed between −1 and 1. The force between pairs of dust grains $Fdd$ is assumed to be a Coulomb interaction since screening is mainly provided by the ions, and this is included explicitly in the model through the force $Fid$⁠. $Fid$ includes the forces from all ions in the simulation region calculated from the Yukawa interaction potential, averaged over the elapsed ion time steps, as well as the force from the homogeneously distributed ions outside the simulation region, as described above.

The confining force $QdE=QdErr̂+Ezẑ$ arises from the electric fields within the simulation. In the present paper, we consider a laboratory experiment where the dust is levitated against the force of gravity $mdg$ by the vertical electric field $Ez$ in the sheath above the lower electrode of a GEC rf reference cell. The sheath electric field is typically assumed to vary linearly with the distance from the electrode. The electric field profiles obtained from a fluid model of the plasma within CASPER's GEC cell show that this is a good approximation for the vertical region where the dust floats and that the electric field steepens with increasing power.³⁰ Here, we are interested in laboratory experiments with vertically aligned dust structures, where confinement in the horizontal direction is provided by the charged walls of a glass box placed on the lower electrode, which has been shown to be very nearly radial near the middle of the box.³¹ Accordingly, we parametrize the forces by the oscillation frequencies $ωr2=Qdmd∂Er/∂r$ and $ωz2=Qdmd∂Ez/∂z$⁠.

The damping factor $β$ depends on the neutral gas pressure and temperature, with

where $δ$ = 1.44 (measured for melamine formaldehyde dust in argon), a is the dust radius, n the gas number density, m_g the molecular mass of the gas, $Tg$ the gas temperature, and m_d the mass of the dust. A thermal bath is realized by subjecting the particles to random kicks, with the maximum acceleration imparted by a kick,

where $Δtd$ is the dust time step.

Typically, dust grains charge through collisions with electrons and ions in the plasma, which constitute currents to the grain's surface that depends upon the dust surface potential. In low-density plasmas, where the ions can be treated as collisionless, the grain charge is often calculated using orbital motion limited (OML) theory.³² However, at higher pressures, charge-exchange collisions between ions and neutral gas particles can generate ions that are trapped around the dust particle, causing OML theory to over-predict the particle charge.^25–27 Here, we use a combined MD and OML approach to determine the dust charge.

The positive charge accumulated each time step is the number of ions that cross the surface of the grain multiplied by the ion charge, $ΔQdi=Nicqi$⁠. The electrons, which are not explicitly modeled, are assumed to have a Boltzmann distribution, and the electron current from OML theory is used,

where $ne$ and $me$ are the electron density and mass, $a$ is the grain radius, and $Φd$ is the dust surface potential. The negative charge accumulated during the time step $Δti$ due to the electron current is $ΔQde=IeΔti$⁠. Thus, the total charge accumulated at each time step is $ΔQ=ΔQde+ΔQdi$⁠.

As an example, in Fig. 2(a), we show the charge calculated at each ion time step (⁠$Δti=10−9$ s) for a single isolated grain in a stationary (ion drift velocity $vd$= 0 M, where the Mach number is in units of the ion sound speed) and a flowing plasma (⁠$vd$ = 0.4 M; $Ez$ = 1300 V/m). The dust grain has a radius of 4.45 $μ$ m and is immersed in argon plasma with an electron temperature of $Te$ = 5.0 eV (58 000 K), an ion temperature of $Ti$ = 290 K, $ne=ni=2×1014$ m⁻³, and a neutral gas pressure of $P$ = 6.67 Pa. The ion plasma period for this condition is $τi=2.1 μ$s. The charge collected in a collisionless plasma (dashed lines) is compared with the reduced charge when ion-neutral collisions are included (solid lines). In the stationary plasma, ion-neutral collisions reduce the charge predicted by OML theory, $Q¯OML=33 400e−$⁠, to $Q¯=$ 12 300 $e−$ corresponding to a grain surface potential of $Φd=$−4.0 V. In the flowing plasma, the effect of ion-neutral collisions is to reduce the charge from $Q¯OML=50 400 e−$ to $Q¯=18 570 e−$ with a grain surface potential of $Φd=$ −5.9 V. Including the ion-neutral collisions in the model allows the dust grains to equilibrate more quickly. In the stationary plasma, the dust charge reaches its equilibrium value within one to two ion plasma periods, while in the flowing plasma, it takes ∼5 ion plasma periods to attain its equilibrium charge.

The orbit motion limited current for flowing ions is given by^33,34

where $ξ=vdr/2kbTi/mi$ (this reduces to the OML current in a stationary plasma in the limit $vdr→0$⁠). Numerical solutions of Eqs. (8) and (9) yield a predicted surface potential of $Φd$= −10.9 V in a stationary plasma and $Φd$= −16.3 V in the flowing plasma, and the results from the numerical model agree well with the predicted charges, shown by the dashed red lines in Fig. 2(a). As expected, due to ion-neutral collisions, the average dust charge is reduced by a factor of 2–3 from that predicted by OML theory for stationary electrons and flowing ions in a collisionless plasma^35–37 (⁠$Q¯OML=$ 2.7 $Q¯$ for these conditions).

As is common in MD simulations of grain charging,^14,19,22 the fluctuations about the average charge are quite large since each ion collected by the dust surface is a superion, with ∼200 elementary charges. The fluctuation in the charge on the particle shown in Fig. 2(b) (⁠$vdr=0.4 M,$ with collisions) has a standard deviation of $σ=$ 1250 e⁻. The theoretical charge variance can be estimated by³⁸ 

where $Ii$ and $Ie$ are the total ion and electron collection currents given in Eqs. (8) and (9) and primes indicate the derivatives with respect to $Q$ (all quantities are evaluated at the grain potential associated with the average charge $Q¯$⁠). Equation (10) predicts a charge standard deviation of $σ$ = 140 e⁻, which is considerably smaller than that of the numerical simulation.

To further complicate matter, large fluctuations (>1 $σ$⁠) can persist over time scales longer than one ion plasma period.^38–40 As the charge fluctuations occur on a fast time scale, the dust is unable to respond to these fluctuations and the dust dynamics should be governed by the average charge. Figure 2(b) shows the charging history of the grain in the flowing plasma $vdr$ = 0.4 M as it is spread over the dust time steps with $NΔti$= 135 per dust time step (a detail of this graph is shown in the inset). Even though the charge is averaged over the previous N elapsed ion time steps to obtain $Qavg$ before calculating the dust dynamics (shown in red), the large charge fluctuations persist on the dust time scale. In this case, the time averaged charge is $Qavg$ = 18 700 e⁻, with a standard deviation of $σ=$ 670 e⁻. To smooth out the large fluctuations on the dust time step, which is several hundred times longer than $135Δti$⁠, the dynamic dust charge is calculated from a weighted average of the dust charge at the previous dust time step and the average charge at the current dust time step,

The coefficients in Eq. (11) are chosen to minimize the charge variance without causing excessive lag in the particle charge. In Fig. 2(b), $Qd$ is represented by the dark blue dashed line. The time averaged dynamic charge is $Qd$ = 18 700 e⁻ with a standard deviation of $σ=$ 370 e⁻. Thus, the charge fluctuations are reduced to ∼2% of the average charge. As shown in the inset, after the dust reaches its equilibrium charge, the dynamic dust charge lags the charge calculated on the ion time step by less than 1 ms, which is relatively fast compared to the dust period of $τd≈$ 100 ms.

In this section, we illustrate the connection between dust particle charging and the ion wakefield by modeling the interaction between two dust grains confined within the sheath above the lower electrode of a GEC reference cell. The simulation parameters are chosen to reproduce a common experiment designed to investigate the transition of a two-particle system from a horizontal to vertical configuration.^41–43 With the dust grains in a vertical pair, a laser pulse can be used to perturb one of the particles, allowing the particle interaction to be probed^8,44,45 and wakefield to be mapped.⁴⁶ 

To examine this experiment numerically, we model the charging and dynamics of a two-particle system with three different ion drift velocities, $vdr$ = 0.4, 0.6, and 1.0 M, where the Mach number M is in units of the ion sound speed $cs=kBTe/mi$⁠. The sheath electric fields that correspond to these drift velocities are $Ez$ = 1300, 2600, and 5500 V/m. The simulation parameters assume spherical melamine formaldehyde dust particles with a radius of 4.45 $μ$m. The plasma conditions for the argon plasma are based on Langmuir probe measurements taken in a GEC cell under similar operating conditions with a gas pressure of 6.67 Pa.⁴⁷ The plasma bulk density is set to $n0$ = 2 × 10¹⁴ m⁻³, the electron temperature is $Te$ = 5 eV (58 000 K), and the ion temperature is $Ti$ = 290 K. The electron Debye length for this set of parameters is $λDe$ = 1.2 mm and the ion Debye length $λDi=$ 83 $μ$m. The cylindrical simulation region has a radius of 1.5 $λDe$ and a height of 5 $λDe$⁠. The dust confinement in the radial and vertical directions is set so that the ratio of confining forces used is $(ωz/ωr) ∼ 1$⁠, slightly less than the instability criterion $ωz/ωrcr=Ω ∼ 3$ for transition from horizontal to vertical pairing for similar plasma conditions.⁴³ In the simulation, the gradient of the electric field is set, and the equilibrium dust charge determines the values of $ωx$ and $ωz$⁠. The parameters employed in the simulation for the three drift velocities are summarized in Table I.

Initial short runs with a single dust particle were made to determine the equilibrium charge and position of the grain. Subsequently, two dust particles were placed near their expected equilibrium positions, with the initial charge on the lower grain reduced to 80%–85% of that on the upper grain. As the lower grain was offset from its equilibrium position, it oscillated vertically within the wake of the upper particle, allowing the charging within the wake to be probed. At time t = 0.15 s, the lower dust particle is given a horizontal acceleration of 0.5 m s⁻² for a time Δt_laser = 0.05 s to simulate a laser pushing the particle outside of the wake. Figure 3(Multimedia view) shows the path of the lower particle P2, relative to the upper particle P1, for three separate runs at the three ion drift velocities. The progression of time is shown as the path changes from light to dark. In the latter two cases, P2 and P1 switch places briefly. For each condition, we examine the characteristics of the ion wake field as the particles interact, the charge as a function of the particle separation, and the acceleration caused by the ion-mediated dust–dust interaction.

The ion density maps for the various drift velocities are shown in Fig. 4(Multimedia view) for three different particle separations. At low ion drift velocities (left column), each dust particle has a distinct ion cloud surrounding it, except when the two particles are very close together, $Δz=0.3 λDe$⁠, as shown in the top row. As the ion velocity increases, the ion density cloud is elongated in the direction of ion flow, and the focusing region moves downward. Thus, for the highest ion velocities, only a single high-density cloud is formed beneath the lower particle.

Maps of the combined ion dust potential for the same conditions are shown in Fig. 5(Multimedia view). The total potential at each point $p=x,z$ is calculated from

where the sum i runs over all the ions, $ΦY$ is the Yukawa potential given in Eq. (3), the sum j runs over the dust grains that interact through the Coulomb potential $ΦC$ given in Eq. (4), and $Φout$ is the potential of the uniformly distributed ions outside the simulation region. As seen in PIC simulations of the ion wake at subsonic velocities with ion–neutral collisions, there is only a single positive potential in the ion wake downstream of the dust, with no oscillations in the potential.¹² At low drift velocities, a positive potential region forms between the two dust grains when the particles are far apart. The region then disappears as the particles approach each other (left column). The positive potential region between the two particles becomes less pronounced as the drift velocity increases and disappears for the highest drift speed $vdr$ = 1.0 M.

Figure 6 shows the charge of the two particles (upper panels) as well as their relative vertical separation (lower panels) as a function of time. Initially, the particles are vertically aligned and the lower particle oscillates up and down within the wake of the upper particle. The charge on the lower particle is significantly reduced as it approaches the upper particle, while the charge on the upper particle remains relatively constant. Figure 7 shows the charge on the two vertically aligned particles during this time (before the laser push is applied) in units of $Q0$⁠, the average charge on the upstream particle. The decharging of the downstream particle (shown in blue) is almost linearly proportional to the vertical separation between the two particles. Interestingly, there appears to be a hysteresis in the charge, depending on whether the downstream particle is approaching or receding from the upstream particle [the hysteresis is evident even though the charge has been corrected for the 1 ms lag caused by smoothing the charge using Eq. (11)]. The velocity of the particles is less than a few cm/s, and so this effect is not due to a difference in the relative streaming velocity of the ions with respect to the dust. Rather, it is likely due to the fact that the lower dust grain is either approaching or receding from the region of high ion density in the wakefield. Note that at the slowest ion drift speed, the charge on the upstream particle is also affected by the location of the downstream particle, with the charge on the upstream particle increasing slightly as the particle separation decreases.

As the lower particle moves out of the upstream particle's wake field, its charge increases and fluctuates about $Q0$⁠. Figure 8 shows the charge on P2 over the entire trajectory. The region where the ion wake decharges the lower grain is clearly distinguishable, with the horizontal extent of the wakefield region decreasing as the ion drift velocity increases. This is more clearly illustrated in Fig. 9, where the relative charge on the downstream particle is shown as a function of relative horizontal separation $Δx$ for several different relative vertical separations, $Δz$⁠. The data shown here are for multiple simulation runs, with data taken over the entire trajectory, in which case the particles may change positions. In this case, the charge on the downstream particle is designated $Qdn$ and the charge on the upstream particle is $Qup$⁠. The smallest drift velocity, $vdr$ = 0.4 M, produces a wake that broadens in the horizontal direction as the vertical separation decreases, and the maximum decharging (minimum value of $Qdn/Qup$⁠) occurs at the smallest particle separation [Fig. 9(a)]. At the highest drift velocity, $vdr$ = 1.0 M [Fig. 9(c)], the horizontal extent of the wake is nearly constant with vertical separation $Δz$⁠, and the charge ratio $Qdn/Qup$ reaches its minimum value at approximately $Δz=0.5 λDe$⁠.

As two vertically aligned particles approach each other in the z-direction, the size of the wake changes and the location of the maximum ion density shifts, eventually forming downstream of the second particle. Figure 10 shows profiles of the average ion density downstream of the two particles for varying locations of P2, comparing the profiles with those of a single grain, which is indicated by the red line. When the particle separation $Δz$ > 0.4 $λDe$ (indicated by the dashed lines), there is a local maximum in the density profile downstream of each particle, the magnitude of which is relatively constant. As the particle separation decreases below 0.4 $λDe$⁠, the two wakes merge, with a single ion focus (local maximum in the ion density) downstream of P2. For subsonic ion flow velocities [Figs. 10(a) and 10(b)], the maximum density in the single combined wake increases markedly.

An estimate of the excess positive charge in the ion cloud associated with each grain is calculated from the charge contained within a region downstream of a dust grain where the ion density $ni>1.6 n0$⁠, as indicated by the contour lines in Figs. 4(g)–4(i). The radial and axial extents of the wake regions as a function of the particle separation are shown in Fig. 11. The filled blue points indicate the wakes where only a single local maximum in the ion density is detected downstream of P2. As noted above, for the subsonic flow velocities, this primarily occurs for $Δz$ < 0.4 $λDe$⁠, and the radial and axial extent of the combined wake expands. For the case where $vd$ = 1.0 M, the ion enhancement downstream of P1 is very weak and barely exceeds the density threshold set to detect the wake, leading to cases where only a single wake is detected downstream of P2 over the range of interparticle distances.

The spatial extent of the wakes (assumed to be azimuthally symmetric) is used to calculate the total excess charge in the wake downstream of each particle,

(where j is the index of the grid points with $ni>1.6 no$ and $Δr$ and $Δz$ are the dimensions of each grid cell). As shown in Figs. 12(a) and 12(b), at subsonic drift velocity, the excess charge below P1 decreases slightly as the particles approach each other [even though the maximum charge density remains relatively constant, as shown in Figs. 10(a) and 10(b)], while the excess charge below P2 increases. Below a distance of $∼0.4λDe$⁠, a single wake is created, with a charge that generally exceeds the wake charge below P2 alone. However, for $vdr$ = 1.0 M, almost the entire wake charge is located downstream of P2 [Fig. 12(c)]. The limiting value of the total charge contained within the wake downstream of P1 and P2 at $Δz=λDe$ is given in Table II, showing that as the ion drift speed increases, the wake shifts below the second particle. At the subsonic speeds investigated, the charge in the lower wake is relatively constant with ion flow speed.

The changing location and magnitude of the focused ion wake has implications for the particle dynamics. Many numerical models of coupled dust motion make use of the “point-charge model,” where the positive ion cloud downstream of a dust particle is assumed to have a fixed charge $q0$ and fixed location $ΔL$ beneath each particle.^17,48–50 This may be a fair approximation if particles are in a horizontal plane. As shown in Figs. 13(a)–13(c), the ion densities downstream of two particles with approximate horizontal alignment (small $Δz$⁠) are nearly the same, though for the small horizontal spacing of 0.4–0.5 $λDe$ shown here, there is a single combined positive potential region downstream of the two grains [Figs. 14(d)–14(f)].

For vertically aligned particles, the result is more complicated as the location and magnitude of the wake charge are affected by the downstream particle. Here, we compare our simulation results, $Φsim$ [calculated from Eq. (12) on the symmetry axis connecting the two particles] with (1) a Coulomb potential, (2) a spherical point charge model, and (3) an ellipsoid point charge model. The Coulomb potential considered for case (1) is given by

where $j$ is the index of the two dust particles. For case (2), the positive point charge is represented by a sphere of charge $qw$ and radius $rw$⁠, as shown in Figs. 12 and 13, centered at the maximum in the ion density, as shown in Fig. 10. The resulting potential as a function of position $z$ along the symmetry axis is given by

In case (3), the positive point charge is represented by an ellipsoid of charge $qw$⁠, with semimajor axes $a=b=rw$ and $c=Δzw2$⁠, centered at the maximum in the ion density. The potential on the symmetry axis as a function of $z$ is given by

where $φs=a2+sb2+sc2+s$⁠, and the lower integration bound is 0 if a point is inside the ellipse or $λ$ if a point is outside the ellipse; $λ$ is the greatest root of $fs=z2c2+s−1=0.$⁵¹ The simulation results for the potential upstream and downstream of the upper particle, $Φsim$ [Eq. (12)], are shown in Fig. 14 by the solid blue lines, for varying particle separations. Within the simulation region, the potential from the ions is smaller at the bottom of the cylinder than at the top (especially for the higher drift velocities), with a slope that ranges from −3.5 V/m for $vd=0.4 M$ to 35 V/m for $vd=1.0 M$⁠. This slope is subtracted from the total potential, for comparison with the Yukawa and point charge potentials, which assume a uniform background potential. For comparison purposes, the zero of each potential is set by its value at a distance $z=λDe$ and then normalized by $Vo=−14πϵo QupλDe$⁠, the magnitude of the Coulomb potential of the upper particle at a distance of $z=λDe$⁠.

At the greatest ion flow speed $vdr$ = 1.0 M [Fig. 14(c)], the potential downstream of P1 is well-represented by the point charge model, as all the excess ion density tends to be downstream of P1. This also means that the upstream potential is matched fairly well by the Coulomb potential. The spherical point charge model is less accurate as the ion drift speed decreases, as it does not fully capture the increased ion density that extends upstream of each particle. In this case, the ellipsoidal point charge model yields the best match at small interparticle spacing and does a better job of capturing the overall shape of the potential for larger particle separations. However, the cutoff used to calculate the total charge in the ion wake (⁠$ni>1.6 nio$⁠) does not capture the full effect of the enhanced ion density in the region of the particles.

The ion wake field has long been successfully employed to explain the non-reciprocal interaction between upstream and downstream particles, as well as the attractive horizontal force that the upstream particle exerts on the downstream particle.^{6,18,43,50,52} The extent of this effect can be examined by calculating the net electric force exerted by the ions on a dust grain from the gradient of the ion potential at the location of the particle, $Fid=Qd∇Φi$⁠.

Figure 15 shows the horizontal acceleration calculated for the grains as a function of their relative position for the three different ion flow speeds. The maps were generated by overlaying the region where the data were collected with a grid with a spacing of 100 $μ$m. The average acceleration $Fdd,x+Fid,xmd$ was calculated for all data points within a radius of 70.7 $μ$m of each grid point. The effect of the ion wake field downstream of the particle (marked by the black dot at the origin) is readily apparent in Fig. 15, where positive values indicate acceleration toward the right. There is a wedge-shaped region below the upstream grain over which the attractive force from the wake potential exceeds the repulsive particle interaction force. The location of this region moves downward and becomes narrower as the ion flow speed increases.

In the vertical direction, the electric force exerted by the ions on the dust grains is always downwards, reflecting the positive potential region formed downstream of the dust. However, there is an asymmetry in the drag, which is exerted on each particle, as shown in Fig. 16, which is in the center-of-mass (COM) system. Since the particles have the same mass, in the COM system, $x1, z1COM=−x2,z2COM$⁠. An upstream particle (⁠$zCOM>0$⁠) always experiences a larger force from the ions than a downstream particle. As two particles approach the COM in the vertical direction, the downward force on the upstream particle increases, as does the upward force on the downstream particle. The net effect of the ion-dust force in the vertical direction is to push the two particles together, in effect reducing the repulsion between the grains. This asymmetry lessens as the ion drift speed increases.

The DRIAD simulation allows for the self-consistent calculation of particle charging as the dust grains move in a streaming ion environment. Here, we examined ion wakefield characteristics, dust charging processes, and the resulting effects on the dust dynamics for three choices of ion velocity: $vdr$ = 0.4 M, 0.6 M, and 1.0 M, where the Mach number M is in units of the ion sound speed. Specifically, we considered the case where a downstream grain is perturbed by a simulated laser pulse, allowing it to move both within and outside the wake of the upper particle, as shown in Fig. 3.

A similar experiment was modeled by Miloch et al.¹⁴ using a PIC code (DiP3D) for the flowing ions and electrons, while keeping the dust configurations static. Although Miloch et al. examined supersonic ion velocities in the range (1.1–2.5) M, they assumed lower electron temperature (3 eV vs 5 eV used here), yielding a smaller value of the ion sound speed. Thus, the sonic drift speed $vdr$ = 1.0 M in our simulation is close to their reference ion drift speed of 1.2 M. The large grain sizes used in the DiP3D simulation, $r$ = 0.1 $λDe ≈$ 400 μm, result in larger grain charges and increased ion focusing, which were compensated to some extent by using a lower plasma density. Nevertheless, a comparison between DRIAD and DiP3D simulations yields several results, which are qualitatively in agreement. Miloch et al. showed that for $vdr$ = 1.2 M, the stable position of the dust pair occurs at particle separation $z≈1.1λDe$⁠, which corresponds to the maximum in the potential distribution in the wake of a single grain. In our simulation, the downstream particle oscillates about a position approximately 0.6 $λDe$ (⁠$≈$ 0.72 mm) downstream of the upper particle, with the maximum wakefield potential for the sonic case located at $z≈1λDe$ (⁠$≈$ 1.2 mm) downstream of the second particle [see Figs. 5(c), 5(f), and 5(i)]. Although the exact location (in mm) of the wakefield maximum does not coincide, we confirm the qualitative characteristics of the sonic/supersonic flow wakefield, including its elongated narrow structure and the tendency to form a single wake below both grains for small vertical separation of the dust grains [see Figs. 4(c), 4(f), and 4(i)].

The location of the wake maximum (in mm) predicted in the present simulation is in agreement with the experimental study by Jung et al.⁵³ They estimated that the wake is localized between 0.7 mm and 2 mm downstream of the upper particle, while the estimated position using DRIAD is 1.2 mm for the sonic case examined here. For vertical separation smaller than 0.7 mm, Jung et al. also established that a single ion focus is formed downstream of the lower particle, which is well reproduced by the present simulation. We found that a single wakefield below the particle pair forms for vertical separation smaller than 0.4 $λDe$ (0.48 mm) in the subsonic case [Figs. 10(a) and 10(b) and 11(a)–11(e)], while at sonic flow speed, the same effect can be observed for all examined vertical separations [up to 0.8 $λDe$ = 0.96 mm, as shown in Figs. 10(c) and 11(c) and 11(f)]. Although we did not consider larger vertical separations, in Fig. 13(c), we show that the vertical extent of ion wakefields in the sonic regime can reach $≈$2 mm, suggesting that wake-induced decharging can be observed at such large vertical separation, as demonstrated by Jung et al.

Miloch et al.¹⁴ found that the total positive charge contained in the wake of a grain was about 15% of the charge on the upstream grain. Here, we find that at large interparticle separations $≥λDe$⁠, the charge in the wake downstream of P2 is ∼30% of the charge on the upstream grain. While the wake charge is relatively constant for sonic drift velocity $vd$ = 1.0 M, at subsonic velocities as the particles approach each other and the wakefields merge, the charge in the combined wake increases substantially and becomes comparable to the charge on the upstream grain. Another feature that is clearly seen in the DRIAD simulation, as in the DiP3D simulations, is the non-linear addition of the wakes for small interparticle distances (Fig. 10), where the additional negative charge of the downstream grain enhances the ion focus.¹⁴ At subsonic ion drift speeds, the nonlinear superposition is not only evident when the two grains are aligned with the vertical flow, but also the shape of the combined wake is substantially altered when the particles are at oblique angles to the direction of ion flow.

The location and existence of a combined wakefield directly affect charging processes as the downstream particle moves toward and away from the upstream grain. The charge on the downstream grain at the closest approach is found to be decreased by 20%–30% as compared to the charge of a single isolated grain. This value is comparable to the results from the study by Miloch et al., who found 14%–36% decharging, and Jung et al., who estimated a charge reduction of 16% at larger particle separations. Miloch et al. found that horizontally aligned grains have independent wakes and charges, but that when perturbed, the lower grain begins decharging as soon as it enters the Mach cone originating from the ion focus below the upper grain. At the subsonic and sonic velocities modeled here, the simulation results show that decharging is related to the spatial extent of the ion wake. As the ion drift speed increases, the vertical extent of the region where the grain is decharged increases, even as the radial extent of this region narrows. We also find that horizontally aligned grains have independent wakes and charge [Figs. 13(a)–13(c)].

Jung et al. reported that decharging is reduced and almost disappears for particle separations smaller than 0.5 mm (Fig. 5 of that paper). Our simulations can reproduce this result only for cases where the perturbed particle (P2) moves above the unperturbed one (P2) or when, in addition to small vertical separation, there is a horizontal separation of at least 0.3 mm–0.5 mm (as shown in Fig. 8). In other words, depending on the ion flow velocity, the decharging effect disappears at small separations only if P2 is outside the Mach cone created by P1. In cases where the wakefield from P1 is located between the two grains (typical for the subsonic velocities), decharging of the lower grain increases as the vertical separation is reduced since P2 is entering the wake of P1. This can be seen in Fig. 9(a), where the charge ratio $Qdn/Qup$ reaches its minimum when the vertical separation is about 0.3 $λDe$ = 0.36 mm (light blue open circles in that figure). In contrast, when a combined wake forms downstream of the dust pair (characteristic of the sonic case), $Qdn/Qup$ reaches its minimum at about 0.5 $λDe$ = 0.6 mm and decharging is reduced for the smaller separation [Fig. 9(c)]. Thus, we conclude that the experimental observations examined by Jung et al. probably correspond to sonic/supersonic ion velocity, which is possible for their experimental conditions (dust particles in the sheath of the rf plasma close to the lower electrode).

It should be noted that in the DRIAD simulation, the electric field accelerating the ions is constant, resulting in a constant flow velocity, rather than the (almost) linear electric field that exists in the sheath. The linear sheath electric field accelerates the ions, and the number density of both the electrons and the ions decreases toward the lower electrode, with the electron density decreasing at a faster rate. Thus, the charge on the dust will be a function of the location within the sheath as well as the location in the wakefield, changing the dynamics of the particles. Future simulations need to include the effects of the changing electron and ion density to better recreate the conditions experienced by particles in the sheath.

As Piel reported for the initial results of the MAD code,¹⁹ the horizontal restoring force diminishes as the ion drift velocity decreases, which is confirmed by the present simulation (as shown in Fig. 15). The acceleration maps in Fig. 15 reveal that the maximum in the horizontal restoring force also moves outward from the alignment axis as the ion flow velocity decreases. For this reason, the restoring force does not vanish for small drift velocities, such as $vd$ = 0.4, as the non-linear superposition of the wakes increases the ion density between the two particles and not just downstream of each particle. As also reported by Piel,¹⁹ there is no net vertical attraction on the lower particle from the wake of the upper particles. Piel concluded that the lower particle feels a stronger attraction from its own wake charge than from the wake charge of the upper particle. However, in the MAD code, the charge on each particle was held fixed at the same value and did not include the decharging of the lower grain. Here, we find that the difference in the downward ion drag acting on each particle effectively pushes the two particles together. The upstream particle experiences a stronger drag force, and the reduced drag on the downstream particle allows it to draft in the upstream wake.

Miloch et al. noted that in experiments the perturbed lower grain moves to a metastable position slightly downstream and to the side of the upper grain, with the lower grain slowly subsiding from this position over several seconds.¹⁴ While we do not reproduce the experimental conditions that lead to this metastable position, its existence is hinted at in the trajectories of the grains for $vdr$ = 0.6 and 1.0 M [Figs. 3(b) and 3(c)], where the perturbed grain can be seen to oscillate about a position to the side of the upper grain. As the lower particle reaches the edge of the Mach cone, it decharges and moves below the upper grain (Fig. 5), as also observed in the experiments.

Maps of the ion wakefield for an interacting dust pair were made in an experiment performed in the CASPER lab where the horizontal confinement for the particles was provided by the charged walls of a glass box placed on the lower powered electrode.⁴⁶ Careful adjustment of the system power and pressure increases the horizontal confinement, allowing the particles to arrange themselves into a vertical pair. A laser pulse was then used to perturb the lower particle, and the motion of the two particles was then recorded using a side-mounted high-speed camera. Analysis of the motion of isolated dust particles was also used to reconstruct a map of the electric field within the region, which was then subtracted from the motion of the paired particles, allowing the ion-mediated dust interaction to be measured. The radial and vertical extent of the region where the ion wakefield provides a restoring force most closely matches the simulation results for the lowest ion drift velocity $vdr=0.4 M$⁠, as shown in Figs. 15(a) and 16(a). Although the magnitude of the vertical “attractive” force is similar, the horizontal acceleration measured in the experiment is lower by an order of magnitude than the accelerations found in the simulation. This may indicate that the ion drift velocity in the experiment is lower and the particle charges smaller than those in the simulation. These conditions can result from plasma depletion in the region due to ion/electron collection on the walls of the glass box.

The importance of simultaneously modeling the grain charge and the dust dynamics becomes evident in the appearance of a hysteresis in the grain charging, correlated with the grain motion in the wake. As shown in Fig. 7, the charge on the lower grain varies by $±3%$ depending on whether the downstream grain is approaching or receding from the upstream particle. As this hysteresis is evident even though the charge has been corrected for the 1 ms lag caused by smoothing the charge using Eq. (11), we expect it to result from an actual charging delay due to the finite charging time. A mechanism for such hysteresis was described by Nunomura et al.,⁵⁴ who proposed that below critical thresholds for gas pressure and electron density, the delay in charging as a dust grain moves through the plasma sheath can lead to energy gain and the onset of self-excited vertical oscillations of dust. According to this model, the charge on dust has a memory of its history of motion, which leads to non-Markovian kinetics. It is predicted that the particle will gain (lose) energy if $E→⋅∇Qd−eq<0$ (⁠$E→⋅∇Qd−eq>0$⁠), where $E→$ is the electrostatic field and $∇Qd−eq$ is the gradient of the equilibrium charge on the dust particle. In other words, the particle gains energy if the charge $Qd−down$ while falling down in an electrostatic potential is more negative than $Qd−up$⁠, the charge while the particle climbs up.

Nunomura et al. found this condition to be satisfied for grains levitating just below the sheath boundary, leading to a spontaneous oscillation of the particles in a dust crystal as the gas pressure was decreased. In our study, this condition is satisfied for the lower grain oscillating within the wake of the upper grain, with the largest gradient in the charge apparent for the lowest ion drift velocity (Fig. 6). As seen in Fig. 3(a), near the end of the simulation, the lower grain returns to its position beneath the upper grain, but the oscillation amplitude increases in time [Fig. 6(d)]. Such instabilities have also been observed in laboratory experiments where dust particles in a given vertical configuration exhibit a short period of rapid motion as a critical threshold is reached (in either plasma power or pressure) before transitioning into a new stable structure.^4,55 Using DRIAD to match the conditions where particular stable configurations exist or to predict a transition through an instability will allow investigation of the role the ion wakefield-dust interaction plays in structural stability.

A numerical model is presented, which simultaneously tracks the dynamics of ions and dust in a plasma environment with flowing ions. The dust charge and the ion wakefield are self-consistently calculated from the ion dynamics, allowing for detailed analysis of the wakefield-mediated interaction as the structural configuration of the dust grains changes.

As one particle approaches another, the downstream particle is decharged within the high-ion density region of the ion wake. As the ion drift speed approaches $vdr=$ 1.0 M, the magnitude of the decharging and the horizontal extent over which it occurs are relatively constant within the elongated wake. An interesting feature revealed by the dynamic charge calculation is the hysteresis in the charging/decharging of the lower grain as it moves up and down within the wake, especially evident at subsonic ion flow speeds.

The spatial extent of the ion wake characteristics as well as the location and magnitude of the maximum ion density changes with ion flow velocity as well as the relative dust position. As the particles approach each other, less charge is contained in the wake of the upstream particle, while the charge in the wake of the downstream particle grows. For subsonic ion flow, the wakes behind the two particles merge as they approach closer than 0.5 $λDe$⁠, creating a single maximum positive potential downstream of the second particle. For ion flow with $vdr=1.0$ M, however, the extent and maximum density of the wake are relatively independent of the particle separation, with the wake almost entirely downstream of P2 for all particle separations.

The non-reciprocal interaction force between the two dust grains is readily apparent. In the horizontal direction, the downstream grain experiences an attractive force, which increases in magnitude as the ion drift speed increases. Although there is no net attractive force in the vertical direction, there is an asymmetry in the interaction force due to the difference in the force exerted by the ions on the dust. The upstream particle always experiences a larger downward force compared to that of the downstream particle, an effect which decreases with the ion drift speed.

The simulated conditions shown here were for two dust grains confined within a glass box placed on the lower electrode of a rf cell. However, the boundary conditions in the numerical model can easily be adapted for different experiments such as the DC glow discharge in the PK-4 experiment onboard the International Space Station.⁵⁶ In this experiment, the polarity of the DC electrodes can be rapidly switched in order to trap the dust grains in the field of view of the cameras. This can create an ion wake both upstream and downstream of the dust grains, resulting in a homogeneous-to-string structural transition of the dust cloud, an interesting condition for further research.

Thanks to Peter Hartmann and Marlene Rosenberg for useful discussion about this work. Support from NSF Grant Nos. 1707215 and 1740203 and NASA Grant No. 1571701 is gratefully acknowledged.