On the wake structure in streaming complex plasmas

The plasma parameters are scaled to relevant base units, which are, in our case, the Debye length $\lambda _{{\mathrm { D}}_\alpha }=\sqrt {\varepsilon _0 k_{\mathrm {B}} T_\alpha /(n_\alpha q_\alpha ^2)}$ and the plasma period τ_α = 2π/ω_pα, with the corresponding plasma frequency $\omega _{\mathrm {p}_\alpha }=\sqrt {{n_\alpha q_\alpha ^2}/({\varepsilon _0 m_\alpha })}$, where α denotes electrons (e) or ions (i). T_α is the corresponding temperature, n_α density, m_α mass, q_α charge, k_B the Boltzmann constant and ε₀ the permittivity of vacuum. The relative streaming velocity ${\bf u}_{\mathrm { i}}$ of the ions with respect to the dust is characterized by the Mach number M ≡ u_i/c_s, where the ion sound (Bohm) speed is given by $c_{\mathrm {s}}\equiv \sqrt {{k_{\mathrm { B}} T_{\mathrm {e}}}/{m_{\mathrm {i}}}}$. Furthermore, the thermal velocity is $v_{T_\alpha }=\sqrt {{k_{\mathrm { B}} T_\alpha }/{m_\alpha }}$.

The dynamically screened dust potential and in particular the characteristics of the wake are studied with respect to three dimensionless plasma parameters:

• the ion drift velocity, expressed in terms of the Mach number M,
• the electron-to-ion temperature ratio T_r ≡ T_e/T_i,
• the ion–neutral scattering frequency (normalized to the ion plasma frequency) ν_in/ω_pi.

While the temperature ratio T_e/T_i controls the effect of Landau damping, the frequency ratio ν_in/ω_pi determines the effect of collisional damping and can be effectively controlled in experiments by changing the neutral gas pressure. A simple approximation for the pressure as a function of the ion–neutral collision frequency and ion temperature is given by

$$\begin{equation} p=\frac{\nu_{\mathrm{in}}/\omega_{\mathrm{p}_{\mathrm{i}}}\cdot e_0 \sqrt{k_{\mathrm{B}} T_{\mathrm{i}}\cdot n_{\mathrm{e}}/\varepsilon_0}}{\sigma_{\mathrm{in}}}, \end{equation} \tag{ 1 }$$

where e₀ = 1.602 × 10⁻¹⁹ C is the elementary charge, and σ_in ≃ 5.0 × 10⁻¹⁹ m² is the ion–neutral collision cross-section as defined in [57]⁶.

For the LR calculations, we consider argon, which is a typical working gas in complex plasma experiments [8]. The ions are considered to be singly charged (q_i = e₀). Consequently, ion and electron densities are equal, and are set to n_e = n_i = 2.0 × 10¹⁴ m⁻³ and are considered as spatially homogeneous distributions. The considered electron temperature is T_e = 2.585 eV (≃30 000 K). The corresponding electron Debye length is λ_{D e} = 845 μm. An atomic mass of argon m_n = m_i = 6.634 × 10⁻²⁶ kg leads to the Bohm speed of c_s = 2500 m s⁻¹. For the considered set of parameters, the ion plasma frequency takes the value of ω_pi = 3.0 × 10⁶ Hz. For the PIC simulations, we simulate ions with a reduced mass m_i/m_e = 120, and thus ω_pi and c_s are increased accordingly (for details see section 3.1). Note that the LR and PIC results presented in this paper do not involve any kind of free fit parameters.

The dust grains are charged by plasma currents and other currents, such as those induced by the photoelectric effect or secondary electron emission. For the charge equilibrium, the net current to the grain is zero, and such a grain is at floating potential Φ_fl with respect to the surrounding plasma. If charging is only due to plasma currents (i.e. electron and ion fluxes to the surface), then in electropositive plasmas, the charge Q_d on the grain and the floating potential will be negative, due to the high mobility of electrons. The floating potential on a spherical, finite-sized grain can be approximated with the orbit–motion–limited (OML) approach for both stationary and flowing collisionless plasmas [51]. We note that in the presence of collisions, the OML theory can give incorrect results, as has been noted by several authors [59–61]. The presence of trapped ions created by ion–neutral collisions can increase the ion current to the grain, and thus reduce the charge on the grain [62–64].

In a stationary plasma, the charging of a grain takes approximately one ion plasma period, τ_i, a time in which the ions can collectively respond to the perturbations in the electric field. A typical charging characteristic for a single, spherical conducting grain is shown in figure 1, where the data are taken from 3D PIC simulations for the stationary plasma with T_e/T_i = 30. In the static case, the charge on a small grain can be calculated with the capacitance model for a given floating potential [65]:

$$\begin{equation} Q_d=4\pi \epsilon_{0} a \left( 1+ \frac{a}{\lambda_{\mathrm{D}}}\right) \Phi_{\mathrm{fl}}, \end{equation} \tag{ 2 }$$

where a is the grain radius and λ_D the total Debye length defined in terms of the electron and ion Debye lengths λ⁻²_D ≡ λ⁻²_De + λ⁻²_Di for the static case, i.e. M = 0. When charging is only due to plasma currents, the floating potential in a collisionless plasma is typically of the order of Φ_fl ≈ −2k_BT_e/e₀. Thus, from the capacitance model, one can derive a formula for an approximate charge on a negligible small grain with a ≪ λ_D:

$$\begin{equation} Q_d/e_0=-1400 a_{[\mu \mathrm{m}]} T_{\mathrm{e}, [\mathrm{eV}]} , \end{equation} \tag{ 3 }$$

where the grain radius a is given in micrometers and temperature in eV. In the case of no plasma flow relative to the dust grain, the negatively charged grain is shielded by a stationary, positively charged ion cloud. The electric potential around a small (non-absorbing) charged grain takes the form of a statically screened Coulomb potential, which is typically referred to as the Debye–Hückel or Yukawa potential

$$\begin{equation} \Phi_{\mathrm{Yuk}}(r)=\frac{1}{4\pi \varepsilon_0} \frac{Q_d}{r} \mathrm{e}^{-\frac{r}{\lambda_{\mathrm{D}}}}. \end{equation} \tag{ 4 }$$

Plasma absorption, ion–neutral collisions and plasma flow can modify the potential distribution. Notably, the long-range asymptote of the plasma-mediated potential is often not screened exponentially, but has a power-law decay instead; see, for example, [33, 53, 54, 67, 68].

Zoom In Zoom Out Reset image size

Plasma flow breaks the symmetry in the grain charging. In the mesothermal velocity regime, when the flow velocity u_i is larger than the ion thermal velocity v_Ti but less than the electron thermal velocity v_Te (v_Ti < u_i < v_Te), the grain will receive a larger ion flux to the surface on the upstream than on the downstream side, while the electron flux will be similar on both sides. Anisotropic charging can lead to an electric dipole on the grain [30, 66]. The electric dipole moment depends on the material of the grain and is more pronounced for supersonic flows [37]. For small conducting grains the electric dipole can often be neglected, as the charge is redistributed on the surface, but for larger conducting grains, the potential distribution in the vicinity of a grain can modify the charge distribution on the grain surface.

A streaming plasma leads to a wake in the density and potential distributions around the grain. The main mechanisms for the wake formation are plasma absorption on the grain surface and the influence of the electric field on plasma particle trajectories in the vicinity of the grain [41]. Plasma absorption leads to density depletion on the downstream side, while electric fields can scatter ions into the wake region forming the ion focus region. Both effects are more pronounced for cold ions, as well as for supersonic flows, when also a Mach cone forms.

The charging time in a streaming plasma is longer than in a stationary plasma and depends on the material of a dust grain, being longer for insulating grains, for which it can take several ion plasma periods to reach stationary conditions.

Considering a weak (linear) response of the plasma to the presence of the dust grain, the electric potential can be computed by a 3D Fourier transformation

$$\begin{equation} \Phi({\bf r}) = \int\!\mathrm{d}^{3}{\bf k} \frac{Q_d}{2 \pi^2 k^2 \epsilon({\bf k}, {\bf k} \cdot {\bf u}_{\mathrm{i}})} \mathrm{e}^{\mathrm{i} {\bf k} {\bf r}} , \end{equation} \tag{ 5 }$$

where ${\bf r}$ denotes the distance from the grain's center of mass. The (shifted) Maxwellian plasma is represented by the following longitudinal dielectric function, which includes Bhatnagar–Gross–Krook (BGK)-type ion-neutral collisions [69, 70]:

$$\begin{equation} \epsilon^l({\bf k},\omega)= 1 + \frac{1}{k^2 \lambda_{\mathrm{D}_{\mathrm{e}}}^2} + \frac{1}{k^2 \lambda_{\mathrm{D}_{\mathrm{i}}}^2} \left[ \frac{1+\zeta_{\mathrm{i}} Z(\zeta_{\mathrm{i}})}{1+\frac{{\rm i}\nu_{\mathrm{in}}}{\sqrt{2} k v_{T_{\mathrm{i}}}}Z(\zeta_{\mathrm{i}})} \right] , \end{equation} \tag{ 6 }$$

where we use the substitution

$$\begin{equation} \zeta_{\mathrm{i}}=\frac{{\bf k}({\bf v}_{d}-{\bf u}_{\mathrm{i}})+\mathrm{i}\nu_{\mathrm{in}}}{\sqrt{2} k v_{T_{\mathrm{i}}}}. \end{equation} \tag{ 7 }$$

The screening by electrons is considered as static (Yukawa-type), since the electron flow is of no relevance (u_e ≪ v_Te). The plasma dispersion function is defined as

$$\begin{equation} Z(\zeta)= \frac{1}{\sqrt{\pi}}\int_{-\infty}^\infty \frac{\exp(-t^2)}{t-\zeta} \mathrm{d}t = \mathrm{i}\sqrt{\pi} \exp(-\zeta^2) \mathrm{Erfc}(-\mathrm{i} \zeta) . \end{equation} \tag{ 8 }$$

The function Z(ζ) and the product ζ·Z(ζ) are plotted in figure 2. In the limiting case of M = 0, that is, $({\bf v}_d-{\bf u}_{\mathrm { i}}) \rightarrow 0$, equation (6) reduces to static electron and ion shielding, which corresponds to the spherically symmetric Yukawa potential equation (4). For $|{\bf v}_d-{\bf u}_{\mathrm { i}}| \gg v_{T_{\mathrm {i}}}$ the product ζ_i·Z(ζ_i) approaches −1 such that the ion-screening contribution vanishes and only electron screening remains in equation (6); see the right panel of figure 2. The finite electron contribution ensures the numerical convergence of the Fourier integral equation (5).

Zoom In Zoom Out Reset image size

The accuracy of the dynamical screening model depends crucially on the 3D discrete Fourier transform (DFT), equation (5). A sufficient sampling rate in k- and r-space at the same time requires a large number of sample points, which are, however, in 3D typically limited by computer resources to N_x = N_y = N_z = 256. Even with the optimal choice of the sampling interval, corresponding to the grid spacing in the reciprocal space, the result for many of the relevant plasma parameters is often not satisfying and lacks accuracy. In particular, the cutoff error in the streaming direction may lead to pseudoperiodicity effects from the DFT, such as artificial oscillations in the upstream direction.

A straightforward approach to reduce the numerical complexity of the Fourier integral equation (5) is to make use of the axial symmetry in the streaming direction by introducing cylindrical coordinates. However, the numerical evaluation of the resulting radial Hankel transform is very time-consuming and converges poorly due to an oscillating Bessel function in the integrand. Instead, an alternative approach was used here to improve the memory and time efficiency by exploiting the radial symmetry. The 3D integration area was sliced along the k_z-direction and N_z 2D-DFTs were performed on a k_x,y = 256 × 256 grid. Notably, for each k_z-value only the resulting N_x/2 = 128 values in one radial direction (we used the positive half of the x-axis) have to be stored and used for the subsequent k_z-integration. Hence, instead of N_x·N_y = 256², only N_x/2 = 128 DFTs in the k_z-direction have to be computed. Moreover, using this scheme the integration area can be easily adjusted to asymmetric spatial dimensions. In particular, for M ⩾ 0.5 in the case of weak Landau and collisional damping, we used N_z = 512 instead of N_z = 256 grid points, which considerably reduced the pseudoperiodicity effects and ensured convergence of $\Phi ({\bf r})$.

For the sake of completeness, we should mention that because of different length scales of the Debye potential and the wake modulations we did not transform $\Phi ({\bf k})$ directly, but the difference $\Delta \Phi ({\bf k}) = \Phi ({\bf k}) - \Phi _{\mathrm { Yuk}}({\bf k})$ [46], where $\Phi _{\mathrm { Yuk}}({\bf k})$ denotes the collisionless static case. Subsequently, utilizing the linearity of the Fourier transform, the Yukawa potential equation (4) is added in real space again. In real space, the DFT yields a grid resolution of λ_{D e}/8 for the wake contribution $\Delta \Phi ({\bf r})$. The data are post processed with a spline interpolation.

In order to ensure the convergence of the Fourier integration, sufficient collisional and/or Landau damping is required. More specifically, for the collisionless case (ν_in = 0) the temperature ratio is limited to T_e/T_i ⩽ 25. For reliable results over the full range of M and temperature ratios up to T_e/T_i = 100, weak collisional damping ν_in = 0.1ω_pi is required. The Fourier integration performs better for small values of M.

The LR calculations discussed here use the plasma parameters (electron temperature, density, ion mass, etc) as introduced in section 1.1. Neglecting the dust drift, i.e. $|{\bf v}_d|=0$, the only dust parameter which enters the calculations is the grain charge, which is set to Q_d = −10⁴e₀ = −1.602 × 10⁻¹⁵ C. Since the wake potential of a point-like grain (a/λ_De ≃ 0) scales linearly with Q_d, a pre-computed potential $\Phi ({\bf r})$ can easily be adapted to any other grain charge required; see equation (5).

Figure 3 shows the shape of the dynamically screened (wake) potential $\Phi ({\bf r})$ in real space for a range of typical plasma parameters taken from the experiments [8, 23]. To capture the general trends, we consider three characteristic electron-to-ion temperature ratios T_r = T_e/T_i = 10,30,100 (top to bottom row) and drift velocities in the subsonic regime (M = 0.5, left), at the Bohm speed (M = 1.0, middle column) and supersonic ion flow (M = 1.5, right). The ion–neutral scattering frequency is fixed at a value of ν_in/ω_pi = 0.1. According to equation (1) this corresponds to pressures of p = 31,18 and 10 Pa for T_r = 10,30,100, respectively. The plasma flow gives rise to an oscillating wake structure downstream of the grain. On account of that, the grain potential $\Phi ({\bf r})$ essentially deviates from the static Yukawa potential in all considered cases. Typically, the potential's anisotropy, i.e. the symmetry breaking, increases with M. The wakefield is most pronounced for large values of T_r.

Zoom In Zoom Out Reset image size

For M = 0.5 and T_e/T_i = 10 (p = 31 Pa), the range and height of the wakefield oscillations are strongly reduced by the overlapping effect of Landau and collisional damping. Only a relatively weakly pronounced primary potential maximum, i.e. a single node directly behind the grain, emerges in the subsonic regime.⁷ For M = 1.0, and even more pronounced for M = 1.5, the wave front becomes cone-shaped. In the supersonic case, the peak is lower than for M = 1.0. The same is observed for the second (third) potential peak for T_r = 30 (T_r = 100). Further maxima and minima besides the first downstream potential maximum start to appear for temperature ratios T_r > 10. Considering T_r = 100, there are in total three (significant) positive space charge regions, which may result in an attractive (non-reciprocal) force between the equally charged dust grains. For M = 1.0 the form of the wake extends far in the cross-streaming direction. With growing M the wake becomes increasingly stretched in the streaming direction while the peak height goes down. A plane Φ = 0 wave front is observed for M ≈ 0.875 (T_r = 30).

The potential profile can be studied in more detail in figure 4, where cuts through the potential surface are plotted for M = 0.375,0.6875,1.0 and T_r = 30. For the considered grain charge of Q_d = −10⁴e₀ the potential height takes a value around 10⁻²k_BT_e/e₀ ≈ 25 mV, which is consistent with  [22, 42]. Evidently, an oscillatory wake is not only formed when the speed of the ion flow exceeds the ion-acoustic velocity [54, 71]. Rather, a significant ion focus with a potential peak height of 21.7 mV is already found for M = 0.375, i.e. well below the ion sound speed. In fact, the highest amplitude of the wakefield is reached at M = 0.6875; see also figure 5. In the upstream direction, the range of Φ_para is already significantly enhanced for M = 0.375, while in the cross-stream direction Φ_ortho is still relatively close to the corresponding Yukawa potential for that value. When M is increased, the (predominantly repulsive) potential profiles in the upstream and cross-stream directions approach each other, as discussed in the context of equation (9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the following, we will discuss the functional characteristics of the wakefield downstream from the grain. The potential peak height and its position as a function of M are plotted in figure 5. The positions of the maxima and minima of the wake potential are found to shift linearly away from the grain when M is increased. Averaging over the three temperature ratios T_e/T_i = 10,30,50 the best linear fit a_iM + b_i for the ith peak is obtained for the parameters a_i = (0.848,3.91,7.19,11.1,15.8) mm and b_i = (0.0553,−0.327,−0.745,−1.54,−2.78) mm; see dashed gray lines in figure 5 (top). Corresponding PIC results for the wake maxima (see table 1 in section 3) are shown with diamond symbols. Despite the finite and relatively large grain size considered in the PIC simulations, there is good agreement between both methods. On the basis of that result, a longitudinal wavelength of the wake equal to 2πλ_DeM, discussed in  [42], holds in the relevant range of M as a reasonable approximation.

In figure 5, we also show the magnitude of the potential peak heights as a function of M. As discussed in the context of figure 4, the wake amplitude reveals a clear maximum when M is varied. Interestingly, this maximum shifts to the subsonic regime when the ratio Te/Ti is increased (in figure 5, an arrow indicates the spline interpolated maximum). According to the OML theory [72], the dust charge QTr OML increases monotonically with M for the considered velocities, as shown for three temperature ratios in the bottom left panel of figure 5 (see the dashed black curves). Taking the charge increase into account, the maximum value of the peak height shifts to slightly higher values of M as displayed for the first potential peak (black solid line).

The left panel of figure 6 shows the temperature dependence in more detail. For high ion temperatures, i.e. T_r ⩽ 5, the largest peak amplitude is found in the supersonic regime, i.e. M ⩾ 1.5. An increase of T_r leads to a monotonic rise of the potential height (the effect of Landau damping is reduced). The greatest slope of the peak height with respect to T_r is observed for M = 0.5. So that in the limit of a large temperature ratio (a low ion temperature), the largest wake amplitude is found in the subsonic regime, that is, M = 0.5, whereas for the largest Mach number, M = 1.5, the amplitude of the wake oscillation is considerably smaller. The maximum value of the first peak's height as a function of M and its position are denoted by open black circles for T_r = 10,30,50,100. The spline through these data points reveals that the strongest first peak in the wake (the primary potential maximum) over the full range of Mach numbers shifts with an increase of T_r closer to the dust grain (top left panel). We note that a similar trend as that observed for the fixed grain charge (black dashed line), is found when one includes the effect of the charge increase according to OML theory (not shown). The right panel of figure 6 shows that an increase of ion–neutral scattering ν_in leads to an increased damping of the wake oscillation with the peak positions being changed only slightly.

Zoom In Zoom Out Reset image size

So far, we have mainly focused on the potential's wake. The form of the potential upstream and laterally from the grain can be captured by an effective screening length [33, 36, 70, 73, 74]. Figure 7 displays the Debye length as a function of M for different temperature ratios T_r. For low ion drift, u_i ≪ v_Ti, both electrons and ions provide shielding, but ion screening typically dominates since T_i < T_e. For zero flow, M = 0, the isotropic Yukawa potential, equation (4), applies with the total Debye length being λ_D = 0.30,0.14 in units of λ_De for T_r = 10,50, respectively. In the static case, the negatively charged (Coulomb) potential of the dust grain is predominantly screened by a positively charged ion cloud, which isotropically forms around the grain. In contrast, in the supersonic limit M ≫ 1, (isotropic) Debye screening is entirely contributed by the electrons⁸. For finite M, the characteristic screening length is found to be different in the upstream and cross-stream directions, with the screening in the upstream direction being larger; see figure 7. It needs to be mentioned that the characteristic screening lengths were obtained by fitting the analytic Yukawa potential to potential cuts through the dust grain as shown in figure 4. However, these potential cuts do not exactly match the functional form of the Yukawa potential. For instance, for M = 0.375 (left panel in figure 4) it is found that the orthogonal slice slightly overshoots Φ = 0, i.e. the potential possesses a weakly attractive branch side ways from the grain⁹. Therefore, the screening lengths in figure 7 cover the form of the potential close to the grain, but have an approximate character only. Also, a systematic error is observed in the limit M → ∞, where the value of the effective screening length in the upstream direction converges toward a value that is about 5% larger than the appropriate screening length, that is, the electron Debye length λ_De. In general, we find that the Yukawa potential matches the potential profile perpendicular to the flow essentially better than in the upstream direction, which is in agreement with [73].

Zoom In Zoom Out Reset image size

Finally, we wish to point out that the characteristic screening length can be roughly approximated by the interpolation formula

$$\begin{equation} \lambda_{\mathrm{s}}= \lambda_{\mathrm{D}_{\mathrm{e}}} \left(\frac{k_{\mathrm{B}} T_{\mathrm{i}} +\eta m_{\mathrm{i}} u_{\mathrm{i}}^2}{k_{\mathrm{B}} T_{\mathrm{e}} + k_{\mathrm{B}} T_{\mathrm{i}} + m_{\mathrm{i}} u_{\mathrm{i}}^2}\right)^{1/2} = \left(\frac{1}{f_{\mathrm{s}}(\eta) \lambda_{\mathrm{D}_{\mathrm{i}}}^2} + \frac{f_{\mathrm{s}}(1)}{f_{\mathrm{s}}(\eta) \lambda_{\mathrm{D}_{\mathrm{e}}}^2}\right)^{-1/2}, \end{equation} \tag{ 9 }$$

where $f_{\mathrm {s}}(\eta )=1+\eta M^2 \frac {T_{\mathrm { e}}}{T_{\mathrm {i}}}$ and η is a heuristic fitting parameter which we introduce here. Considering η = 1, λ_s reproduces both screening limits correctly: (i) λ_s = λ_D for M = 0 and (ii) λ_s = λ_De for M → ∞ (see the dotted lines in figure 7) [36]. However, λ_s is found to converge much more slowly toward λ_De than the curves obtained by LR theory (figure 7) and PIC simulations (as will be shown in section 3.3, figure 11). Therefore, it is suitable to introduce the fitting parameter η = 1.27, which allows for an essentially better agreement with the present results in the relevant range $0\leqslant M \lesssim 2$.¹⁰ The concept of an effective screening length in a plasma with streaming ions has been also discussed by Khrapak et al ([74] and p 45 of  [5]). Here, the analytical form of the effective screening length is approximated by the formula

$$\begin{equation} \lambda_k = \left(f_{1(2)}(M_T) \lambda_{\mathrm{D}_{\mathrm{i}}}^{-2} + \lambda_{\mathrm{D}_{\mathrm{e}}}^{-2}\right)^{-1/2}, \end{equation} \tag{ 10 }$$

where the two fitting functions $f_{1}=\exp (-M_T^2/2)$ and f₂ = (1 + M²_T)⁻¹ are defined in terms of the thermal Mach number M_T ≡ u_i/v_Ti = Mc_s/v_Ti. Evidently, with the fitting function f₂, equation (10) is identical to equation (9) with η = 1 and is therefore not further considered in the following. For T_r = 10, equation (10) yields in conjunction with f₁ a very reasonable approximation (especially for the screening in the upstream direction) without any heuristic parameter; see dashed gray line in figure 7. Furthermore, λ_k shows (in contrast to λ_s) the correct analytic trend of a rapidly increasing slope with increasing temperature ratio T_r, which is evident in the intersection of curves at different temperatures. However, λ_k overestimates the decay of the ion contribution to the screening as shown for T_r = 50 in figure 7 and for T_r = 30 in figure 11.

The theoretical finding that the effective screening length converges toward the electron Debye length λ_De as M increases is supported by experimental results [75, 76]. In those experiments the effective screening length for the grains suspended in the sheath was found by analyzing collisions between the dust grains. In both cases it was shown that supersonic ions are unable to screen the dust particles and therefore the electron Debye length is the appropriate screening length for dust grains in supersonic plasma flows.

Analytical plasma models that also account for nonlinear effects are difficult to develop [77]. Therefore, one often employs numerical kinetic particle simulations, in which plasma particle trajectories are followed in self-consistent force fields without many approximations, allowing for a self-consistent evolution of the system. As plasma particle trajectories are the characteristics of the Vlasov equation, the PIC simulations can be considered as an alternative way of solving the Vlasov equation for arbitrary particle distributions [78].

Plasma simulations, in which forces are evaluated between all pairs of simulation particles, are typically very expensive in terms of computational time and memory. A simulation of n particles, with such direct force calculations, has complexity $\mathcal {O}(n^2)$, as approximately n² equations need to be solved to find the forces on all the particles. $\mathcal {O}$ relates to the computer resources usage time. As the number of simulated particles needs to be large, especially in 3D simulations, such an algorithm would be very inefficient.

Introducing a grid within the PIC scheme can significantly reduce the complexity of the algorithm [79]. In the PIC method, the physical quantities of each simulation particle are weighted to neighboring grid points to build appropriate density fields on the grid. In the electrostatic approximation, the electric potential is found from the charge density on the grid points by solving the Poisson equation. Thus, the complexity of the PIC algorithm is $\mathcal {O}(n)+\mathcal {O}(n_{\mathrm {g}} \log (n_{\mathrm {g}}))$, where n_g denotes the number of grid points with n_g ≪ n. $\mathcal {O}(n_{\mathrm {g}} \log (n_{\mathrm {g}}))$ is the complexity of solving the field equations.

Integration of particle trajectories puts limits on PIC simulations. With both electrons and ions simulated, the time resolution should be a fraction of the electron plasma period. The charging time is usually of the order of the ion plasma period and thus requires a large number of time steps. Simulations can be accelerated by assuming Boltzmann-distributed electrons, which leads to a hybrid PIC–fluid code [80]. This approximation is, however, not always valid, as trapped electrons or electron sources (e.g., due to photoemission) can give rise to local deviations from the Boltzmann distribution [81]. Another way of speeding up the evolution of the system is by using a reduced ion mass m_i, as the ion plasma frequency increases with a decrease in m_i. Ion mass values as low as m_i = 30m_e have been used in the literature [82]. In [37], it has been demonstrated that mass ratios m_i/m_e > 100 give reliable results for dust charging. While the reduction of ion mass leads to some quantitative differences, the results are qualitatively correct and scalable and can be presented in normalized units. Note that with the ion mass being reduced, the sound speed c_s is increased and the floating potential is slightly reduced. In turn, with a relatively large grain, the discreteness noise and charge fluctuations on the grain are diminished.

In the present study, we use the DiP3D (dust in plasma 3D) PIC code, which has been developed to study the charging of dust grains and other objects in various plasma environments [41, 51]. The code simulates a collisionless plasma in the electrostatic approximation, with electrons and ions represented as individual plasma particles. To weight particle charges to the regular grid and to project forces from the grid to the particles, first order linear weighting is used. The particle trajectories are advanced with the leap frog method characterized by a staggered time mesh [79]:

$$\begin{eqnarray} &&\begin{array}{l} x_{\mathrm{i}}(t+\Delta t) = x_{\mathrm{i}}(t)+v_{\mathrm{i}}(t+\Delta t/2)\Delta t,\\ v_{\mathrm{i}}(t+\Delta t/2) = v_{\mathrm{i}}(t- \Delta t/2) + f_{\mathrm{i}}(t) \Delta t /m_{\mathrm{i}}, \end{array} \end{eqnarray} \tag{ 11 }$$

where i refers to an electron or ion, f_i is the force projected on the ith particle from the nearest grid points and Δt is the computational time step. The plasma particles are initially uniformly distributed with Maxwellian velocity distributions in the simulation box of size L³ = (24λ_De)³. The particles are injected into the simulation box at each time step according to the prescribed particle fluxes. To simulate an open plasma system, the particles can leave freely through the boundaries of the simulation box. The code uses the Dirichlet boundary conditions for the potential. The potential at the boundaries is set according to the plasma potential Φ_pl = 0 V.

One or more spherical dust grains are placed inside the simulation domain far away from the boundaries. We assume a conducting grain, and each plasma particle that hits the surface is removed and contributes to the current to the grain. The charge of the removed particle is distributed on the grain surface to cancel internal electric fields. The grain is initially uncharged and becomes charged self-consistently during the simulation.

The DiP3D code is run on a computer cluster, with typically n = 10⁷ or more simulation particles distributed on several nodes, with a typical grid size of n_g ⩾ 128³. The message-passing-interface (MPI) is used for parallelization. The code can be stopped and then restarted with modified dust grain configurations. This feature is used to reduce the computation time in the case of systematic parametric studies. The code allows also for movement of the simulated object and for the inclusion of various charging processes such as photoemission. Recently, the code was upgraded to include external magnetic fields and ion–neutral and electron–neutral collisions [83]. With all these features, the code can be placed among other cutting-edge particle plasma codes for object–plasma interactions [36, 56, 84]. In this work, we make use of only the basic features of the code.

It should be noted that the dimensionality of the code is an important aspect of any simulation. For the ultimate goal of direct application of the numerical results to support laboratory studies, 3D simulations should be used. Many important effects, such as ion focusing, are still present when the dimensionality is reduced. However, in a 2D system, the charge of the plasma particle corresponds to the charge of an infinite rod, and a simulated circular grain represents an infinite cylinder. The plasma frequency ω_p, Debye length λ_D and charge density do not change with the dimensionality of the system. Reducing the dimensionality can give faster algorithms, but while the results of such simulations are qualitatively correct, they do not need to represent a 3D system quantitatively. Thus, the results from codes with reduced dimensionality should be compared with the corresponding theoretical models. In this paper, we present the results of 3D simulations.

With the DiP3D code, we performed simulations of a single grain in a collisionless plasma. The grain radius is a = 0.185λ_De = 0.1563 mm. The finite size of the grain leads to enhanced potential variations in the wake as compared to the previous LR results, because the self-consistent grain charge is of the order of Q_d ≈ −10⁵e₀. All other plasma parameters are as described in section 1.1. The grain charging and wake formation have been investigated for different electron-to-ion temperature ratios, $T_{\mathrm { r}} \in \left \{10,30,100\right \}$, and Mach numbers, $M \in \left \{0.75,1,1.5\right \}$, as shown in the contour plots of the electric potential in figure 8. In all cases, we observe a positive maximum in the potential distribution in the wake at a distance of approximately 1–2 λ_De downstream from the dust grain. This potential maximum is related to the ion focus in the wake of the negatively charged grain. The maximum in the potential is located further downstream from the grain than the enhancement in the ion density in the focus region, which is consistent with the Poisson equation ∇²Φ = −ρ/₀, where ρ is the charge density. Ion focusing is due to the deflection of ion trajectories by electric fields in the vicinity of a negatively charged grain [28, 49].

Zoom In Zoom Out Reset image size

Similarly to the LR results in section 2, in the PIC simulations the potential enhancements in the wake become more pronounced for larger M and larger electron-to-ion temperature ratios T_r, i.e. for colder ions. The first local potential maximum is followed by other potential extrema further downstream. For T_r = 100 and M = 0.75, the pattern in the wake shows a diverse structure with several minima and maxima being also located away from the flow symmetry axis. Note that this could also be an artefact due to the Poisson solver on the finite grid. In particular, for very low M and T_r, the possibility for instabilities can become an important issue for the convergence of the code [85]. For larger M and T_r, the extrema in the potential distribution are along the direction of the flow. With an increase of the ion temperature, the Landau damping becomes important, and thus ion focusing and wake patterns are less pronounced. For T_r = 10 and M = 0.75, only a weak, single maximum is observed. The considerable structural consequences of even a weak ion focus were the subject of recent investigations [52].

In figure 9, we present potential cuts along the flow direction z at y = 0 and x∈[0,2λ_De] for T_r = 30. In the figure, we also plot the potential cut in the cross-stream direction, as well as the potential cut for the grain in stationary plasma. In contrast with the LR results, the potential maximum monotonically increases with the flow speed, and so does its extension in the perpendicular direction whilst the Mach cone forms. However, there is an important difference between the LR and the PIC results. In the LR analysis, the charge on the grain is constant, whereas in the PIC simulations it changes with the flow. The floating potential on a dust grain in a stationary plasma is Φ_fl ∼ − 2.0 V and becomes larger with increasing M (see the discussion below). While the charge and potential on the grain increase with the flow, the amplitude of peaks in the wake potential is being enhanced.

Zoom In Zoom Out Reset image size

In order to compare the PIC results with LR calculations, in the left panel of figure 10, we present potential cuts through the grain along the z-direction and in the cross-stream direction for M = 1.0 and T_e/T_i = 30 as obtained from collisionless 3D PIC simulations and from LR calculations with collisional damping included (ν_in/ω_pi = 0.1). The Yukawa potential, equation (4), for a grain in stationary plasma is also shown. The grain charge of the LR result is adjusted to match the potential of the first peak in the PIC simulation (the same charge is used for the Yukawa potential). It can be inferred from the plot that the agreement between the PIC and LR results is good for the first two extrema. It is important to note that in contrast with the PIC simulations, the LR results do not take into account the finite particle radius a = 0.185λ_De. For this reason, the LR curves are slightly shifted toward zero. When the finite size of the grain is taken into account, the agreement between the LR and PIC results for the potential profiles is improved.

Zoom In Zoom Out Reset image size

As was discussed in the context of figure 5, we observe good overall agreement in the longitudinal wavelength between the PIC and LR results¹¹. An essentially larger discrepancy in the wavelengths was recently reported by Hutchinson [42]. Therefore, in the right panel of figure 10 a comparison for T_e/T_i = 50 and M = 1.0 is shown between our LR results, the COPTIC results for a point-like particle (figure 11 of [42]) and the linear calculation by Lampe et al (figure 3 of [33]). The main difference between Lampe et al's and our LR results lies in the fact that we include finite collisional damping (which in our case is necessary to avoid pseudoperiodicity effects for the given temperature ratio T_e/T_i = 50, as stated in section 2.2). The wavelength obtained from COPTIC is in perfect agreement with our LR calculation¹². The stronger decay of the wake oscillations can be attributed to finite collisional damping which is included in the LR calculations (ν_in/ω_pi = 0.1). However, in the left panel of figure 5, we observe a much better agreement in the wake amplitude at the first subsidiary valley. The stronger damping in our PIC simulation (left panel) must be due to the large grain size (and charge) giving rise to much stronger nonlinear effects than the point-like grain in the COPTIC simulation. As discussed in [42], a major effect of nonlinearity is the local increase of the effective ion temperature in conjunction with an enhanced Landau damping of the wake amplitude. The presence of (strong) nonlinearities in the present PIC (DiP3D) simulations becomes apparent from the similarity of the shape of the primary potential maximum for M = 1 and the temperature ratios T_r = 30 and especially T_r = 100 (figure 8) and the nonlinear wake structure presented in figure 5(a) in [42]. Therefore, we assume that the good agreement between our collisionless PIC simulation and the LR calculation with collisions included—not only in the wavelength of the wake oscillations but also in the decay of the wake amplitude—is due to a nonlinearity-induced damping mechanism that reduces the amplitude of the wake potential.

In the downstream direction, the primary repulsive branch of the potential is shifted closer toward the grain at finite (subsonic) flows, which is similar to the LR results (cf figures 4 and 9). As ions contribute more to the grain screening in stationary plasma than in flowing plasma, grain screening is strongest in the static (Yukawa) case and weakens for finite M; see again figure 9. A quantitative comparison of the PIC and LR results has been made by calculating the effective screening length from the potential profiles for different M for T_r = 30 in the cross-stream and upstream directions; see figure 11. The general trend of the PIC and LR results is that the effective screening length λ increases with the flow velocity M. The effective screening for the stationary plasma in the PIC simulations is λ = 0.23λ_De. This value is larger than the the results from the LR calculations and also larger than the total Debye length, which is λ_D = 0.18λ_De. The difference is due to a finite grid resolution in the simulation, where the grid spacing is Δx ≈ λ_D. For flow velocities M < 0.5, the cross-stream and upstream shielding lengths are less than 0.4λ_De. The values are roughly 0.1λ_De larger (smaller) than the cross-stream (upstream) results from LR calculations. For M > 0.5, the effective screening length from the PIC simulations strongly increases with M, with a steeper increase for the upstream direction. This is the same behavior as for the LR results. At M ⩽ 1, for both the LR and the PIC results, the screening length λ is larger on the upstream side than on the downstream side. The increase of the effective screening length with the flow is in agreement with previous studies [33, 70]; however, details of the trend for the cross-stream and upstream directions are different. In  [33] the cross-stream shielding was systematically larger than the upstream one. On the basis of our results from simulations with two different numerical methods, the upstream shielding length is larger than the cross-stream one for M ⩽ 1. For higher M, the trend continues for the LR calculations, whereas for PIC simulations, the shielding becomes larger on the cross-stream side than on the upstream side. The effective screening length converges toward the electron Debye length for supersonic flows, as was discussed in the previous section for the LR results.

Zoom In Zoom Out Reset image size

Finally, in table 1, we summarize the floating potential on the grain as well as the positions and heights of the maxima for different M values. Since the ion flux to the grain surface changes with the ion flow speed, the floating potential Φ_fl on the grain becomes more negative with increasing M. Changing T_r has little influence on the floating potential, as Φ_fl is mainly controlled by the electron temperature T_e, which is kept constant in our simulations. The heights of the potential maxima increase with the flow speed, and so does the distance between the peak and the center of the grain.

The dynamical screening potentials obtained from the LR approach are topologically similar to the corresponding results from the PIC simulations; see again figures 3 and 8. The PIC results are more noisy, as they account for the dynamics of plasma particles and fluctuations in the system. Since the PIC simulations consider a finite-sized grain that is self-consistently charged, the potential on the grain changes with the flow speed; see table 1. A comparison of the peak positions of the positive ion wake charges, figure 5, shows good agreement. Note that the exact position of weak peaks with a low signal-to-noise ratio (high peak numbers at small T_r and M values) can barely be resolved; cf figure 3.

The amplitude of the peak in the wake calculated with LR achieves a maximum at a certain ion flow velocity (figure 5). This maximum is at lower flow velocities for larger temperature ratios. We do not see a maximum in the peak amplitude for the PIC results summarized in table 1. However, in section 2, we considered the charge on the grain to be constant, whereas in the PIC simulations the potential on the grain changes with the flow. When scaling linearly all PIC results to a fixed potential, we obtain a maximum in the height of the peak as a function of flow velocity. This maximum is at lower velocities for higher temperature ratios, similarly to the LR results. Moreover, adjusting the correct grain charge, the potential pattern of the wakes obtained from the LR calculations agree well with the PIC results; see again figure 10.

Both LR theory and PIC simulations confirm that the effective screening of the grain changes with the flow. While the effective screening length calculations with the PIC method can be considered to be less accurate due to the potential representation on a finite grid, the trend for the effective screening length agrees with the LR results. For subsonic ion flow velocities, the screening length on the upstream side is usually larger than on the cross-stream side.

An important question for the OCP model is to determine to what extent the wake potential around two or more grains can be approximated by the linear superposition of single-grain potentials. We will thus study the impact of nonlinear effects in the multiple grain wakes by comparing the results from PIC simulations of two grains with the results obtained by a linear combination of data from the corresponding simulation of a single grain.

In figure 12(a), we show the results from the PIC simulations of two grains for T_e/T_i = 30 and M = 1.25. The downstream grain is located off the symmetry axis at a distance 0.93λ_De from the upstream grain. The charge of a downstream grain is affected only slightly by the upstream grain, as the downstream grain is located out of the ion focus region [51]. In the linear combination of two wakes originating from a single-grain calculation, the grain charges (and thus the grain potentials) are adjusted in such a way that they match the corresponding values from the simulations of two grains; see figure 12(b). The constructed wake captures the general features of the results from the simulation of two grains, i.e. the position and approximate height of the potential maxima and minima. However, the small differences between the two results indicate that nonlinear wake effects are present in the simulations of two grains. In figure 12(c), the potential difference ΔΦ = Φ_linear − Φ between the linear combination of wakes and the original result is plotted. In the system of two grains, the downstream grain influences the ion dynamics in the wake, and therefore leads to an asymmetric ion focus region, which cannot be fully represented by a linear combination. The potential resulting from the linear combination overestimates the maximum in the wake by up to 20%. Also the symmetry of the wake potential is broken.

Zoom In Zoom Out Reset image size

The corresponding results for two grains aligned with the flow that are separated by 0.74λ_De are shown in figure 13. The charge on the downstream grain is strongly affected by ion focusing, and it is only about 30% of the charge on the upstream grain. Thus in the linear combination of the wake, the contribution for the downstream grain has been adjusted accordingly. The constructed wake again represents well the topology of the wake from the simulations of two grains. In this case, the relative potential difference is up to 30% with the strongest discrepancies associated with the first maximum.

Zoom In Zoom Out Reset image size

The agreement between the constructed and the simulated wake is worse for subsonic flows. Figures 14 and 15 show the corresponding results for the subsonic case M = 0.75 for the downstream grain located off axis and aligned with the flow, respectively. In both cases, the differences ΔΦ are up to 50% in the closest vicinity of the grains.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In all considered cases the general features of the wake, such as positions of the maxima, are well reproduced. In the constructed wakes, the heights of the first maximum are overestimated, and for asymmetric grain arrangements, the symmetry of the wake further downstream from the grains is also modified. This is due to a single asymmetric ion focus region forming in a system of several grains with small separation, which are not aligned with the flow. Thus, constructing the wake from the linear combination of two wakes of single grains should be done with care. To construct such a wake it is necessary to know a priori the charge or potential on a downstream grain, which can be strongly influenced by the wake effects [50, 51]. While the main features of the wake can be well represented by the linear combination of single-grain wakes, there may be inaccuracy due to the nonlinear effects associated with ion focusing.

In order to study the wake structure behind a grain in streaming plasmas, we have employed two complementary methods: (i) the numerical evaluation of the linearized electrostatic potential that is created by a point-like charge and (ii) self-consistent, full 3D PIC simulations. These methods have the best performance in different regimes. The LR method has the best convergence for small T_r and large ν_in (and of minor importance small M), i.e. large Landau and/or collisional dampings are optimal for the performance of the code. The PIC simulations are easiest for higher M and small T_r, as the simulations for small M and cold ions are often favorable conditions for instabilities [85].

In the full PIC method, the ions often have reduced mass. The results with reduced ion mass are reliable when the analysis is carried out in dimensionless units [37]. The PIC method, which includes self-consistent charging of the grain, accounts for a finite grain size, while the LR assumes a point-like dust grain with a given charge. The considered particle radius a = 0.185λ_De in the PIC simulations corresponds to the nonlinear regime [42].

The two approaches are complementary and together allow for studies of a grain potential over a wide range of parameters. Moreover, for single-grain simulations, there is good agreement in the functional form of the wake (i.e. the number, shape, height and relative positions of the wake's maxima and minima). The main issue of the LR approach is the correct estimate of the grain charge as a function of M, and also for a grain that is shadowed by another grain located upstream. A possible solution to overcome this issue could be the use of pre-computed tables obtained from self-consistent PIC simulations or making use of data from experiments [24]. It is also worth mentioning that possible deviations from a shifted Maxwellian ion velocity distribution, as reported in [42, 86, 87], need to be considered in the development of advanced models of the dust plasma interactions applicable for a wide range of plasma parameters.