On the Wake Structure in Streaming Complex Plasmas

The theoretical description of complex (dusty) plasmas requires multiscale concepts that adequately incorporate the correlated interplay of streaming electrons and ions, neutrals, and dust grains. Knowing the effective dust-dust interaction, the multiscale problem can be effectively reduced to a one-component plasma model of the dust subsystem. The goal of the present publication is a systematic evaluation of the electrostatic potential distribution around a dust grain in the presence of a streaming plasma environment by means of two complementary approaches: (i) a high precision computation of the dynamically screened Coulomb potential from the dynamic dielectric function, and (ii) full 3D particle-in-cell simulations, which self-consistently include dynamical grain charging and non-linear effects. The applicability of these two approaches is addressed.


Introduction
Complex (dusty) plasmas have been proven to be an instructive reference system for strong coupling and correlation effects [1,2,3]. Thereby, complex plasma research inherently connects key issues from several fields including low temperature physics, warm dense matter physics, surface and solid-state physics, as well as material science [4,5,6,7,8,9].
The first principle description of complex plasmas is a theoretically challenging problem. A dusty plasma consists of electrons, positive ions, neutral atoms, and micrometer-sized dust grains, i.e., components with distinct mass asymmetry giving rise to dynamics on very different spatial and temporal scales. This multi-scale behavior, combined with the non-ideal character of a complex plasma, makes a full temporal resolution on all scales a numerically unfeasible task with current computer technology. However, there are two numerical approaches that can address the challenges associated with the complex plasma system: the One-Component Plasma (OCP) model and the Particle-In-Cell (PIC) method.
The OCP approach relies on the idea to reduce the quasi-neutral, multi-component complex plasma to the subsystem of dust grains which interact via a screened Coulomb potential. In particular, the statically screened (Yukawa-type) potential yields good agreement with the experiments for various specific setups, see for example Refs. [10,11]. Also, this isotropic potential is well suited for fundamental analytical and numerical investigations on cooperative effects such as self-organized structure formation [3,12], collective dynamical processes [13,14,15], spectral properties [16,17], or the melting behavior [18,19] in strongly correlated (screened) Coulomb systems. The multi-component and non-equilibrium nature of a complex plasma requires, however, a careful analysis of plasma streaming and dynamical grain charging effects, which were experimentally shown to become of utmost importance at least near the plasma sheath, see e.g. [20,21,22,23,24,25]. In a plasma with ions streaming at a uniform velocity, the dust potential becomes (strongly) anisotropic and takes the form of an oscillatory wake structure in the downstream direction, which was subject to various analytical and numerical studies e.g. [26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42] 1 . This wake can result in attractive, non-reciprocal forces between the equally charged dust grains and lead to remarkable structural and dynamical consequences [43,44,45,46,47,48,49,51,52]. Streaming effects can be incorporated in the OCP dynamical screening model by means of Linear Response (LR) theory, as presented in section 2.
The second powerful numerical ansatz to tackle the multi-scale problem of the different plasma components is the PIC approach, in which the trajectories of plasma particles are followed in electric fields self-consistently [55]. To obtain a self-consistent charge on the grains, the plasma dynamics need to be resolved in the vicinity of the grain on spatial and temporal scales which are smaller than the Debye length and electron plasma period, respectively. However, due to the large mass of grains (as compared to electrons and ions) as well as the large distances between the grains, the study of the dynamics of several dust grains requires simulation times that are much longer than the dust charging (which is on the ion plasma period time scale). Thus, a selfconsistent 3D PIC calculation, i.e., a dynamic evolution of the 3D system with more than one dust grain with a temporal resolution on the electron or ion plasma period time scale, required to ensure self-consistency, easily exceeds the capabilities of present high-performance computing systems [51,56].
Consequently, short-time, small-scale PIC simulations of dust charging or, alternatively, less demanding LR calculations of the dynamical plasma shielding have to be coupled with large-scale OCP simulations which incorporate the interaction between the grains on a more abstract level [46,47,48]. Such a coupled multi-scale numerical approach may ensure the description of the correlated system dynamics with proper charges on the grains as well as an accurate potential distribution in the vicinity of grains. With the goal to form a method-spanning picture, the present publication comprises a systematic analysis of the electrostatic potential distribution around a dust grain in the presence of a streaming plasma environment by (i) a high precision computation of the dynamically screened Coulomb potential from the dynamic dielectric function, and (ii) a critical assessment of these LR results, in particular in the view of non-linear effects and dynamical grain charging processes by means of self-consistent 3D PIC simulations.
The discussion of the LR results, section 2, extends over a broad range of plasma streaming velocities, different electron-to-ion temperature ratios and includes an evaluation of the influence of collisional damping. The corresponding 3D PIC simulations, section 3, are performed for the collisionless limit and apply, therefore, to the low pressure limit. Finally, it is vital to verify whether and to what degree the results obtained by a linear superposition of single grain potentials do coincide with the corresponding full 3D PIC simulations. Such a direct comparison is performed for the case of two grains in section 4.

Plasma Parameters
The plasma parameters are scaled to relevant base units, which are, in our case, the Debye length λ Dα = ε 0 k B T α /(n α q 2 α ) and the plasma period τ α = 2π/ω pα , with the corresponding plasma frequency ω pα = n α q 2 α /(ε 0 m α ), where α denotes electrons (e) or ions (i). T α is the corresponding temperature, n α density, m α mass, q α charge, k B the Boltzmann constant, and ε 0 is the permittivity of vacuum. The relative streaming velocity u 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 s ≡ k B T e /m i . Furthermore, the thermal velocity is v Tα = k B T α /m α . 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 ratio of electron-to-ion temperature T r ≡ T e /T i , • the ion-neutral scattering frequency (normalized to the ion plasma frequency) ν in /ω p i .
While the temperature ratio T e /T i controls the effect of Landau damping, the frequency ratio ν in /ω p i 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 function of the ion-neutral collision frequency and ion temperature is given by where e 0 = 1.602 · 10 −19 C is the elementary charge, and σ in 5.0 · 10 −19 m 2 is the ion-neutral collision cross-section as defined in Ref. [57]. 2 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 0 ). Consequently, ion and electron densities are equal, and are set to n e = n i = 2.0 · 10 14 m −3 and are considered as spatially homogeneous distributions. The considered electron temperature is T e = 2.585 eV ( 30000 K). The corresponding electron Debye length is λ De = 845 µm. An atomic mass of Argon m n = m i = 6.634 · 10 −26 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 ω p i = 3.0 · 10 6 Hz. For the PIC simulations, we simulate ions with a reduced mass m i /m e = 120, and thus ω p i and c s are increased accordingly (see section 3.1 for details). The LR and PIC results presented in this paper do not involve any kind of free fit parameters.

Grain charging and wake formation
The dust grains are charged by plasma currents and other currents, such as 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 Φ f l 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 it has been noted by several authors [59,60,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,63,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 Fig. 1, where the data is taken from 3D PIC simulations for the stationary plasma with T e /T i = 30. In the Typical charging characteristics for a spherical conducting grain in stationary plasma. Data is from 3D PIC simulations of a stationary plasma with T e /T i = 30, time t is normalized with the ion plasma period τ i . static case, the charge on a small grain can be calculated with the capacitance model for a given floating potential [65]: where a is the grain radius and λ D the total Debye length defined in terms of electron and ion Debye lengths λ −2 D ≡ λ −2 De + λ −2 D i for a static case, i.e., M = 0. When charging is only due to plasma currents, the floating potential in a collisionless plasma is typically on the order of Φ fl ≈ −2k B T e /e 0 . Thus, from the capacitance model one can derive a formula for an approximate charge on a negligible small grain with a λ D : 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 Debye-Hückel or Yukawa potential Plasma absorption, ion-neutral collisions, and plasma flow can modify the potential distribution.
In particular, the long-range asymptote of the plasma-mediated potential is often not screened exponentially, but has a power-law decay instead, see e.g. [33,53,54,67,68]. The plasma flow breaks the symmetry in the grain charging. In the mezothermal velocity regime, when the flow velocity u i is larger than the ion thermal velocity v T i and lower than the electron thermal velocity v Te (v T i < 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 to 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, 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.

Linear Response Approach
According to the OCP concept, the multi-component dusty plasma can be reduced to the subsystem of dust grains which interact via a dynamically screened Coulomb potential. This dynamical potential takes into account the effect of plasma streaming on the dielectric response of the plasma. Depending on the quality of dielectric function (k, ω), the considered dynamical screening model comprises an accurate representation of the most essential plasma properties including screening, wakefield oscillations, ion and electron thermal effects, Landau damping, as well as collisional damping [33,48].

Theoretical Approach
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 where 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] l (k, ω) = 1 + 1 where we use the substitution 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 The function Z(ζ) and the product ζ · Z(ζ) are plotted in Fig. 2. In the limiting case of M = 0, that is (v d − u i ) → 0, Eq. (6) reduces to static electron and ion shielding, which corresponds to the spherically symmetric Yukawa potential Eq. (4). For |v d − u i | v T i the product ζ i · Z(ζ i ) approaches −1 such that the ion screening contribution vanishes and only electron screening remains in Eq.

Numerical implementation
The accuracy of the Dynamical Screening Model depends crucially on the 3D Discrete Fourier Transform (DFT), Eq. (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 not satisfying and lacks accuracy. In particular, the cut off error in the streaming direction leads to pseudoperiodicity effects from the DFT, such as artificial oscillations in upstream direction.
A straightforward approach to reduce the numerical complexity of the Fourier integral Eq. (5) is to make use of the axial symmetry in 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 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 2 only N x /2 = 128 DFTs in 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/or collisional) damping, we used N z = 512 instead of N z = 256 grid points which considerably reduced pseudoperiodicity effects and ensured convergence of Φ(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 Φ(k) directly, but the difference [46], where Φ Yuk (k) denotes the collisionless static case. Subsequently, utilizing the linearity of the Fourier transform, the Yukawa potential Eq. (4) is added in real space again. In real space, the DFT yields a grid resolution of λ De /8 for the wake contribution ∆Φ(r). The data is 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ω p i is required. The Fourier integration performs better for small values of M .

Results
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. |v d | = 0, the only dust parameter which enters the calculations is the grain charge, which is set to Q d = −10 4 e 0 = −1.602 · 10 −15 C. Since the wake potential of a point-like grain (a/λ De 0) scales linearly with Q d , a pre-computed potential Φ(r) can easily be adapted to any other grain charge required, see Eq. (5). Fig. 3 shows the shape of the dynamically screened (wake) potential Φ(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 Bohm speed (M = 1.0, middle column), and supersonic ion flow (M = 1.5, right). The ionneutral scattering frequency is fixed at a value of ν in /ω p i = 0.1. According to Eq. (1) this corresponds to pressures of p = 31 Pa, 18 Pa and 10 Pa for T r = 10, 30, 100, respectively. The plasma flow gives rise to an oscillating wake structure downstream the grain. On that account, the grain potential Φ(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 .
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. 3 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 Fig. 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 4 e 0 the potential height takes a value around 10 −2 · k B T e /e 0 ≈ 25mV , which is consistent with Refs. [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 ion sound speed. In fact, the highest amplitude of the wakefield is reached at M = 0.6875, see also Fig. 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 profile in the upstream and cross-stream direction approach each other, as discussed in the context of Eq. (9).
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 Fig. 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   Table 1 in section 3) are marked with diamond symbols. Despite the finite, and relatively large grain size considered in the PIC simulations, there is a good agreement of both methods. On the basis of that result, a longitudinal wavelength of the wake equal to 2πλ De M , discussed in Ref. [42],  holds in the relevant range of M as a reasonable approximation. In Fig. 5, we also show the magnitude of the potential peak heights as a function of M . As it is discussed in the context of Fig. 4, the wake amplitude reveals a clear maximum when M is varied. Interestingly, this maximum shifts to the subsonic regime when the ratio T e /T i is increased (in Fig. 5, an arrow indicates the spline interpolated maximum). According to the OML theory [72], the dust charge Q Tr OM L increases monotonically with M as shown for three temperature ratios in the bottom left panel in Fig. 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 Fig. 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. 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 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 increase of T r closer to the dust grain (top left panel). We note that a similar trend, as 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 Fig. 6 displays 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.
So far, we 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]. Fig. 7 displays the Debye length as function of M for different temperature ratios T r . For low ion drift, u i v T i , both electron and ions provide shielding, but ion screening typically dominates since T i < T e . For zero flow, M = 0, the isotropic Yukawa potential, Eq. (4), applies with the total Debye length being λ D = 0.30, 0.18, 0.14 in units of λ De for T r = 10, 30, 50, respectively. In this static case, the negatively charged dust grain Coulomb potential is predominantly screened by a positively charged ion cloud which isotropically forms around the grain. In contrast, in the supersonic limit M 1, (spherically symmetric) Debye screening is entirely contributed by the electrons. 4 For finite M , the characteristic screening length is found to be different in the upstream and cross-stream direction, with the screening in the upstream direction being larger, see Fig. 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 such as shown in Fig. 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 Fig. 4) it is found that the orthogonal slice slightly overshoots Φ = 0, i.e. the potential possesses a weakly attractive branch sideways from the grain 5 . Therefore, the screening lengths in Fig. 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 upstream direction converges towards a value that is about five percent 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 accordance with [73].
Finally, we like to point out that the characteristic screening length can be roughly approximated by the interpolation formula where f s (η) = 1 + ηM 2 Te T 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 Fig. 7) [36]. However, λ s is found to converge much slower towards λ De than the curves obtained by LR theory (Fig. 7) and PIC simulations (as will be shown in section 3.3, Fig. 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 ≤ M 2. 6 The concept of an effective screening length in a plasma with streaming ions was also discussed by S.A. Khrapak et al., see [74] and on page 45 in Ref. [5]. Here, the analytical form of the effective screening length is approximated by the formula where the two fitting functions f 1 = exp(−M 2 T /2) and with the fitting function f 2 , Eq. (10) is identical to Eq. (9) with η = 1 and is therefore not further considered in the following. For T r = 10, Eq. (10) yields in conjunction with f 1 a very reasonable approximation (especially for the screening in the upstream direction) without any heuristic parameter, see dashed gray line in Fig. 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 Fig. 7 and for T r = 30 in Fig. 11.
Theoretical findings 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.

3D Particle-In-Cell Simulations
The LR theory, employed in the previous section, neglects some important processes associated with the dust grain charging and wake formation. In particular, in the regime of strong plasma flows and high grain charges, non-linear effects can be significant. Particle kinetic simulations, such as the PIC method, allow for self-consistent studies of the grain charging, and account also for non-linear phenomena.

Theoretical Approach
Analytical plasma models that account also for non-linear effects are difficult to develop. 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 [77].
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 O(n 2 ), as approximately n 2 equations need to be solved to find the forces on all the particles. 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 [78]. 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 O(n) + O(n g log(n g )), where n g denotes the number of grid points with n g n. O(n g log(n 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 it requires a large number of time steps. Simulations can be accelerated by assuming Boltzmann distributed electrons, which leads to a hybrid PIC-fluid code [79]. 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 [80]. 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 decreasing m i . Ion mass values as low as m i = 30m e have been used in the literature [81]. In [37] it has been demonstrated that mass ratios m i /m e > 100 give reliable results for the dust charging. While the reduction of ion mass leads to some quantitative differences, the results are qualitatively correct, 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.

Numerical implementation
In the present study, we use the Dust in Plasma 3D (DiP3D) particle-in-cell 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 [78]: where i refers to an electron or ion, f i is the force projected on the i-th 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 3 = (24λ De ) 3 . The particles are injected into the simulation box at each time step according to prescribed particle fluxes. To simulate an open plasma system, the particles can leave freely through the boundaries of the simulation box. The code uses Dirichlet boundary conditions for the potential. The potential at the boundaries is set according to the plasma potential Φ pl = 0 V.
One or several 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 7 or more simulation particles distributed on several nodes, with a typical grid size of n g ≥ 128 3 . The Message-Passing-Interface (MPI) is used for parallelization. The code can be stopped and then restarted with modified dust grain configurations, which can reduce the computation time in case of systematic parametric studies. The code also allows for a movement of the simulated object, and for the inclusion of various charging processes such as photoemission. Recently, the code has been upgraded to include external magnetic fields and ion-neutral and electron-neutral collisions [82]. With all the features, the code can be placed among other cutting-edge particle plasma codes for object-plasma interactions [36,56,83]. In the present work, we make only use of 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, results from codes with reduced dimensionality should be compared with corresponding theoretical models. In the present paper, we present results from 3D simulations.

Results
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 on the order of Q d ≈ −10 5 e 0 . All other plasma parameters are as described in section 1.1. The grain charging and wake formation have been investigated for different electronto-ion temperature ratios, T r ∈ {10, 30, 100}, and Mach numbers, M ∈ {0.75, 1, 1.5}, as shown in the contour plots of the electric potential in Fig. 8. In all cases, we observe a positive maximum in the potential distribution in the wake at a distance of approximately one to two λ 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 ∇ 2 Φ = −ρ/ 0 , 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 [49,28].
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., 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 located also 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 of instabilities can become an important issue for the convergence of the code [84]. 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 subject of recent investigations [52].
In Fig. 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 to 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 PIC results. In the LR analysis, the charge on the grain is constant, while in the PIC simulations it changes with the flow. The floating potential on a dust grain in a stationary plasma is Φ fl = −2.02 V and it becomes larger with increasing M (see 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.
In order to compare the PIC results with LR calculations, in the left panel of Fig. 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 /ω p i = 0.1). The Yukawa potential, Eq. (4),  Figure 10.
Left: Potential cuts through the grain along the z-direction (solid lines) and in cross-stream direction (dash-dotted lines) for M = 1.0 and T e /T i = 30 as obtained from collisionless 3D PIC (red) and LR calculations with ν in /ω p i = 0.1 (blue). The grain charge of the LR result and the corresponding Yukawa potential are adjusted to match the potential of the PIC simulation, which includes the grain charging process self-consistently. Note that in contrast to PIC, the LR results do not take into account the finite particle radius a = 0.185λ De , which explains the offset of the curves. Right: Potential cuts along the z-direction for M = 1.0 and T e /T i = 50. Shown are our LR calculations with finite damping included ν in /ω p i = 0.1 (blue solid line), COPTIC data by Hutchinson for a point-like grain in a collisionless plasma (black solid line) [42], and linear response results for the collisionless case by Lampe et al. (red dotted line) [33].
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 to 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 towards 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.
As discussed in the context of Fig. 5, we find an overall good agreement in the longitudinal wavelength of the PIC and LR results. 7 A more inconsistent discrepancy in the wavelengths was recently reported by Hutchinson [42]. Therefore, in the right panel of Fig. 10 a comparison for T e /T i = 50 and M = 1.0 is shown between our LR results, COPTIC results for a point-like particle (Fig. 11 in [42]), and the linear calculation by Lampe et al. (Fig. 3 in [33]). The main difference between Lampe et al.'s and our LR results is due to 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 full accordance with our LR calculation. 8 The stronger decay of the wake oscillations can be attributed to finite collisional damping which is included in the LR calculations (ν in /ω p i = 0.1). However, in the left panel of Fig. 5, we find 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 non-linear effects than the point-like grain in the COPTIC simulation. As discussed in [42], a major effect of non-linearity is the local increase of the effective ion temperature in conjunction with an enhanced Landau damping of the wake amplitude. The presence of (strong) non-linearities in the present PIC (DiP3D) simulations becomes apparent by 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 (Fig. 8) and the non-linear wake structure presented in Fig.5(a) in Ref. [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 towards the grain at finite (subsonic) flows, which is similar to the LR results (cf. Figs. 4 and 9). As ions contribute more to the grain screening in stationary than in flowing plasma, the grain screening is strongest in the static (Yukawa) case and weakens for finite M , see again Fig. 9. The quantitative comparison of the PIC and LR results has been done 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 Fig. 11. The general trend for 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, as well as 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 the flow velocities M < 0.5, the cross-stream and upstream shielding lengths are lower 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 steeper increase for the upstream direction. This is the same behavior as for the LR results. At M ≤ 1, for both LR and PIC results, the screening length λ is larger on the upstream than on the downstream side. The increase of the effective screening length with the flow is in agreement with previous studies [33,70], however, the details of the trend for cross-stream and upstream directions are different. In Ref. [33] the cross-stream shielding was systematically larger than the upstream one. Based on 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, while for PIC simulations, the shielding becomes larger on the cross-stream than on the upstream side. The effective screening length converges towards the electron Debye length for supersonic flows, as it was discussed in the previous section for the LR results. 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 Φ f l 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 center of the grain.

Method-Spanning Comparison
We now summarize the correspondence of the two methods and discuss their limitations.

Comparison of the potentials
The dynamical screening potentials obtained from the LR approach are topologically similar to the corresponding results from the PIC simulations, see again Figs. 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. The comparison of the peak positions of the positive ion wake charges, Fig. 5, shows a 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 just barely be resolved, cf. Fig. 3.
The amplitude of the peak in the wake calculated with LR achieves a maximum at a certain ion flow velocity (Fig. 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, while 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  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 Fig. 12(a), we show 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 only affected a little 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 Fig. 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 non-linear wake effects are present in the simulations of two grains. In Fig. 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 can not 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.
The corresponding results for two grains aligned with the flow that are separated by 0.74λ De are shown in Fig. 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.
The agreement between constructed and simulated wake is worse for subsonic flows. Figs. 14 and 15 show 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.
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 taken 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 a single-grain wakes, there may be inaccuracy due to non-linear effects associated with ion focusing.

Range of applicability of both simulation models
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 particle-in-cell 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 or collisional damping 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 condtitions for instabilities [84].
In the full PIC method, the ions have often 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 non-linear regime [42].
The two approaches are complementary, and together allow for studies of a grain potential over a wide regime of parameters. Moreover, for single grain simulations, there is good agreement in the functional form of the wake (i.e., number and relative positions of the wake potential 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,85,86], need to be considered in developement of advanced models of the dust plasma interactions applicable for a wide range of plasma parameters.

Summary and Outlook
In the present work, we have investigated the electrostatic potential distribution around a spherical dust grain in stationary and flowing plasmas by means of linear response calculations and particle-in-cell simulations. The streaming plasma leads to strong deviations from the statically screened Coulomb potential, giving rise to a distinct oscillatory wake structure behind the grain. For both sub-and supersonic flow velocities, positive ion wake charges downstream from the grain are present that can give rise to non-reciprocal dust grain interactions. We have presented a systematic overview of the LR results for different flow velocities, electron-toion temperature ratios, and ion neutral collision frequencies. The LR results have been tested against self-consistent particle-in-cell simulations. The results from the PIC simulations for the wake pattern agree qualitatively with the LR data when the finite size of the grain is accounted for.
Downstream from the grain, a pronounced oscillatory wakefield with several maxima develops with increasing ion flow. Independent of the temperature ratio T e /T i , the peak positions are shifted linearly away from the grain when M is increased. Upstream and to the side from the grain, the potential decay can be approximated by a Yukawa potential. The effective screening length increases in both directions with the flow, but in contrast to earlier LR results [33], it is found to be larger in the upstream direction than laterally from the grain.
In a system of several grains in a plasma flow, the charge and potential on the downstream grains can be significantly modified by wake effects. The results from PIC simulations show that the wake behind the two grains can include non-linear effects, and that the charge on the downstream grains can be significantly reduced. While the general wake pattern can be represented by a linear combination of wakes behind a single grain to a rather good accuracy, local deviations up to 50% of the potential are observed for subsonic flows in the closest vicinity of the grains. Thus, special care needs to be taken when applying LR results in OCP simulations, and at least the charge reduction on the downstream grains needs to be considered for the model to be credible.
The influence of wake effects on collective many-particle behaviour in 3D plasma crystals has become a very timely question and is currently under ongoing experimental and theoretical investigation [43,44,45,52]. The presented wake potentials are currently employed in Dynamical Dust Simulations which accompany the experimental exploration of string formation and externally triggered phase transitions in 3D confined dusty plasma clouds [44,45]. OCP simulations that include many particles require, however, an adequate adjustment of the charges of the individual grains, which are shadowed by the upstream grains. Such data may be obtained from the first principle PIC simulations or from experiments [24]. The development of such a grain-charging model is the subject of ongoing work.
Due to the high scalability of Coulomb systems, we expect these results to be of interest also to other types of multi-component plasmas. In particular, the generalization of the dynamically screened dust approach will open the way to a description of non-equilibrium quantum systems, such as warm dense matter, on the microscopic level. There, the classical dielectric function must be replaced by the corresponding quantum (Mermin) dielectric function [54,87,88]. This scheme then allows for the computation of collective many-body properties of strongly correlated ions embedded into a partially ionized quantum plasma of electrons by first principle molecular dynamics simulations.