# Ion wakes in two-dimensional complex plasma crystals

Lénaïc Couëdel, Benny Rodriguez, Diana Jimenez Marti, Lorin Matthews

► To cite this version:

Lénaïc Couëdel, Benny Rodriguez, Diana Jimenez Marti, Lorin Matthews. Ion wakes in two-dimensional complex plasma crystals. 2026. hal-05530362

HAL Id: hal-05530362

https://hal.science/hal-05530362v1

Preprint submitted on 27 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Ion wakes in two-dimensional complex plasma crystals

L. Couëdel,$^{1,2}$ B. Rodriguez Saenz,$^3$ D. Jimenez Marti,$^3$ and L. Matthews$^3$

$^{1)}$Aix-Marseille Université, CNRS, PIIM UMR 7345, 3013 Marseille, France

$^{2)}$Physics & Engineering Physics department, University of Saskatchewan, SK S7N 5E2, Saskatoon, Canada

$^{3)}$Center for Astrophysics, Space Physics, and Engineering Research (CASPER), Baylor University, Waco, TX, United States of America

(Dated: 27 February 2026)

Ion wakes in two-dimensional complex plasma crystals play a crucial role in determining the structural and dynamical properties of these strongly coupled systems. Using the DRIAD molecular dynamics code [L. S. Matthews, et al., Phys. Plasmas 27, 023703 (2020)], we investigate the formation and morphology of ion wakes beneath a hexagonal cluster of 19 microparticles under two limiting conditions: a collision-less ion flow and an electric-field-driven collisional flow. Our simulations reveal that, in both regimes, individual wakes form downstream of each microparticle, but their shape and intensity differ significantly. Collision-less flows produce elongated wakes with moderate ion-density enhancement, whereas collisional flows generate shorter, more concentrated wakes due to ion-neutral collisions. A striking feature observed in both cases is the emergence of a collective ion-density enhancement beneath the entire cluster, accompanied by a measurable increase in electric potential (up to ~ 0.5-0.6 V).

## I. INTRODUCTION

Dusty plasmas or complex plasmas are ionised gases containing micron-sized solid particles (dust) that acquire a (usually negative) high electric charge due collisions with the surrounding plasma ions and electrons. The charged dust interacts strongly with each other and the surrounding plasma environment, leading to rich collective behaviour$^1$. Beyond their relevance to basic plasma physics, dusty plasmas occur naturally in astrophysical environments such as planetary rings$^2$, cometary tails$^3$, and interstellar clouds$^{4,5}$, and they also play a role in technological contexts, including semiconductor manufacturing and plasma processing, where dust contamination and control are critical issues$^{6,7}$.

In a capacitively coupled radio-frequency discharge or in a dc discharge, it is possible to form a two-dimensional particle suspension by injecting calibrated spherical microparticles that are trapped in the sheath above a confining electrode where the electric force can balance the gravity force$^8$. Under specific conditions, the microparticle suspension can self-organise into ordered structures, the so-called two-dimensional (2D) complex plasma crystals, where the interparticle forces are dominated by screened Coulomb (Yukawa) interactions$^{8–11}$. These systems are of great interest because they provide a unique platform for studying at the kinetic (particle) level fundamental processes in strongly coupled systems$^{12}$, including phase transitions$^{13}$, wave propagation$^{9,10,14}$, and transport phenomena$^{15}$, at the kinetic level and in real time.

Because of the ion flow coming from the bulk plasma and directed toward the electrode, ion wakes are formed down-stream of each microparticle, creating positive potential structures downstream of negatively charged dust grains$^{16–19}$. The negatively charged particles' interaction is thus a superposition of Yukawa-repulsive interactions with an anisotropic, nonreciprocal attraction$^{16–21}$. These wakes can lead to vertical pairing and symmetry-

breaking configurations$^{22–24}$. Additionally, it has been shown experimentally that negatively charged dust particles levitating above a confining electrode modify the local sheath potential and self-consistently redistribute the plasma particles fluxes, as evidenced by an average increase of the ion drift velocity in the presence of a dust-particle layer$^{25}$.

In 2D complex plasma crystals under specific conditions, these nonreciprocal interactions can couple the longitudinal in-plane dust-lattice mode with the out-of-plane transverse lattice mode and trigger the mode-coupling instability (MCI) that can cause crystal melting$^{26–29}$: energy from the flowing ions is transferred to the microparticles and heats up the microparticle suspension. Experiments and simulations reveal distinct fingerprints of MCI due to ion wake effects, such as hybrid wave modes, mixed polarisation, and wake-induced bending of monolayers$^{27,28,30,31}$.

Understanding ion wakes is crucial for controlling the structural and dynamical properties of complex plasmas. Theoretical models range from simplified point-wake approaches$^{26,32,33}$ to self-consistent kinetic descriptions that incorporate ion-neutral collisions and sheath electric fields$^{19,34}$. Despite significant progress, open questions remain regarding the quantitative mapping of wake parameters to experimental observables and the universality of stability criteria for multilayer and fluid systems$^{35,36}$. Numerical studies emphasise the dependence of wake characteristics on gas pressure and ion flow velocity, while discharge power plays only a minor role$^{37,38}$.

Recent molecular dynamics simulations have shown that, while a single microparticle produces a positive potential region downstream associated with increased plasma density, the addition of multiple grains leads to complex ion wakes and electric fields that cannot be predicted by analytical models$^{39,40}$. Building on this, Watson et al.$^{41}$ introduced an innovative point-charge ap-
proach in which each wake charge is dynamically influ-
enced by neighbouring dust grains, enabling more ac-
curate modeling of large dusty plasma crystals. Never-
theless, no studies have focused on the wake structure
downstream of large two-dimensional clusters, and it re-
mains unclear whether current wake models capture the
essential physics of ion wakes.

This article is organised as follows. In Section II, we describe the simulation methodology, including the DRIAD molecular dynamics code and the initialisation of the dust cluster. Section III describes the steady-state ion velocity distribution for both collision-less and collisional regimes. Section IV presents the simulation results, focusing on ion density and electric potential distributions for the two cases. Finally, Section V discusses the implications of these findings and concludes with perspectives for future work.

II. SIMULATION

A. DRIAD

The DRIAD (Dynamic Response of Ions And Dust) molecular dynamics code was used to calculate the ion wakes downstream of a cluster of microparticles. The code models ions and dust dynamics while self-consistently computing dust charge and was extensively described in other articles³⁷,³⁸,⁴⁰,⁴²⁻⁴⁵. In DRIAD, electrons are approximated by a Boltzmann distribution when the dust potential is sufficiently negative to repel most electrons in the vicinity of the grain⁴⁶. This leads to electron depletion near the dust grain so that ions interact with the grain via a Coulomb potential, and with each other at long range through a screened Yukawa potential³⁹. Since the ions provide the shielding for the negative dust grains, and all ions are explicitly included in the simulation, this shielding is naturally taken into account.

In this work, the simulation domain is a cylinder of
height *H* and radius *R*. The dust microparticles are fixed
at prescribed positions and only the ion motion is calcu-
lated. The equation of motion for each ion *i* is given
by:

$$
m_i \ddot{\vec{r}} = \vec{F}_{ij} + \vec{F}_{id} + \vec{F}_E(z) + \vec{F}(r, z) + \vec{F}_{in}, \quad (1)
$$

where the electrostatic force between the ions *i* and *j*
separated by a distance *r*<sub>*ij*</sub> is:

$$
\vec{F}_{ij} = \sum_{i \neq j} \frac{q_i q_j}{4\pi \epsilon_0 r_{ij}^3} \left(1 + \frac{r_{ij}}{\lambda_{\text{De}}}\right) \exp\left(-\frac{r_{ij}}{\lambda_{\text{De}}}\right) \vec{r}_{ij}, \quad (2)
$$

the electron Debye length is $\lambda_{D_e} = \sqrt{\epsilon_0 k_B T_e / n_e e^2}$, $n_e$ is the electron number density, $T_e$ is the electron temperature, $e$ is the elementary charge, $q_i$ is the charge of an ion, $k_B$ is the Boltzmann constant and $\epsilon_0$ is the vacuum permittivity. The force due to the dust particles $\vec{F}_{id}$

arises from Coulomb interactions between the ions and
the dust particles³⁹:

$$
\vec{F}_{id} = \sum_{d} \frac{q_i Q_d}{4\pi \epsilon_0 r_{id}^3} \vec{r}_{id}, \quad (3)
$$

where the sum is over the *d* microparticles present in the
system and $r_{id}$ is the distance between an ion *i* and a
microparticle. The force from the vertical electric field
in the sheath is given by:

$$
\vec{F}_E(z) = q_i E(z) \hat{z}, \quad (4)
$$

where $\hat{z}$ is the unit vector in the vertical ($z$) direction and $E(z)$ is the sheath electric field at the ion position. In the current research, only a constant electric field is investigated ($E(z) = E_0$).

The force of the external ions (boundary conditions assuming an infinite region of homogeneously distributed ions of density $n_{i0}$ outside the simulation domain) is given by:

$$
\vec{F}(r, z) = q_i \vec{E}_B(r, z) \quad (5)
$$

where $\vec{E}_B(r, z)$ is the electric field within the cylindrical simulation domain when filled with homogeneous Yukawa matter. The corresponding potential is calculated numerically and subtracted from a constant potential to obtain the potential $\Phi_B(r, z)$ within the cylindrical cavity, from which the resulting electric field can be calculated as $\vec{E}_B(r, z) = -\nabla \Phi_B(r, z)$.

DRIAD allows for the inclusion or exclusion of ion-neutral collision effects, represented as $\vec{F}_{in}$. These collisions can play a crucial role in complex plasma environments, particularly in the charging process, when the ion-neutral mean free path becomes shorter than the screening length. To accurately model this phenomenon, ion-neutral collision cross-sections are sourced from the Phelps database⁴⁷ hosted by the LXCat project for both isotropic and back-scattering collisions. Collisions between ions and neutrals are calculated using the null collision method⁴⁸,⁴⁹.

The charge on a microparticle *Q* is calculated in a
self-consistent way by calculating the electron and ion
currents to that grain. Details are given in previous
publications³⁷,³⁸,⁴⁰,⁴²⁻⁴⁵

To calculate the ion density at each ion time step, the
simulation region is discretised into a grid with spacing
(1/8 × 1/8 × 1/6.4)λ<sub>D<sub>e</sub></sub>. For each grid point (x<sub>i</sub>, y<sub>i</sub>, z<sub>i</sub>),
the total number of super-ions within a sphere of radius
R = λ<sub>D<sub>e</sub></sub>/6 centred at that point is stored. The spheres
centred at each grid point overlap slightly, ensuring that
all ions in the simulation are counted and providing a
smoothed ion density. The ion number density at the
grid point is calculated by dividing the total number of
physical ions contained in the super ions by the volume
of the sphere. The same grid is used to calculate the ion
potential. The resulting ion potential at a grid point is
<figure><img src="image_2.png" /><figcaption>FIG. 1. Cluster of 19 particles. The blue circle represents the edge of the cylindrical simulation domain.</figcaption></figure>

obtained by summing the Yukawa-like potential contributions from all individual super-ions in the simulation region. Before plotting, the background ion-density $n_0$ and potential $\Phi_0$ profiles in the absence of dust were also calculated and subtracted from the results obtained with dust under identical boundary conditions in order to isolate the ion wakes.

## B. Initialisation of the dust cluster

The aim of the paper is to demonstrate that ion wakes beneath a monolayer are more complex than a simple superposition of individual wakes and the exact representation of the plasma sheath is beyond its scope. Instead, we focus on two limiting cases: (i) collision-less ion flow in the presence of an external electric field, approximating a collision-less sheath; and (ii) electric-field-driven ion flow with ion-neutral collisions, approximating a collisional sheath. For this reason, a hexagonal cluster with 19 dust grains is implemented in the simulation. Each grain has a radius $r = 4.55$ µm and an initial charge $q_d = -17900$ e (the charging is then evolved self-consistently during the run). To generate the perfectly hexagonal cluster, integer axial indices $(m, n)$ are generated on a 2D stencil and a Boolean mask $\{|m| \le N, |n| \le N, |m+n| \le N\}$, with $N = 2$ for a 19-grain hexagonal cluster. The selected axial coordinates are then mapped to Cartesian positions using the conventional hexagonal lattice conversion,

$$x = \frac{\sqrt{3}}{2} m \lambda_{D_e}, \qquad (6)$$

$$y = \left(\frac{m}{2} + n\right) \lambda_{D_e}, \qquad (7)$$

thereby setting the in-plane interparticle (lattice) spacing to one $\lambda_{D_e}$, which is typical in two-dimensional complex plasma crystals$^{50}$. The cylindrical simulation region has a radius $R = 4.0\lambda_{D_e}$ and height $H = 20\lambda_{D_e}$. The grains are placed at a fixed height $z = 2.5\lambda_{D_e}$. A two-dimensional view of the cluster is shown in Fig. 1. Note that although the dust grains are held at fixed positions, simulations in which the dust grains to oscillate about

their equilibrium positions due to thermal noise show no significant change in the wake structure.

## III. ION VELOCITY DISTRIBUTION

Strongly coupled complex plasmas typically consist of levitating monodisperse spherical microparticles within the sheath above the powered electrode of a capacitively coupled RF argon discharge$^{42,50}$. The pressure range for such experiments spans from the sub-Pascal regime$^{10,50}$ to several tens of Pascals$^{42,51}$, significantly altering the collisionality of the RF sheath. We therefore consider two limiting cases: (i) a collision-less ion flow with a constant electric field, approximating a collision-less sheath; and (ii) an electric-field-driven ion flow, approximating a collisional RF sheath.

In this work, we assume an argon plasma with electron temperature $T_e = 2.15$ eV and an undisturbed density $n_\infty = n_{i\infty} = 8 \cdot 10^8$ cm$^{-3}$.

At the start of the simulations, 131072 macro-ions are injected with random positions and velocities inside the domain (351 real ions per macro-ions). When an ion exits the domain or is absorbed by a dust grain, it is re-injected at the boundary to keep the total number of macro-ions constant. Re-injection uses random positions and velocities sampled from a distribution based on the differential flux of particles across the boundary, derived from a shifted Maxwellian distribution$^{46}$. To reduce edge effects, the injection boundary is set to be at least two Debye lengths away from the dust layer. Further details on ion (re-)injection are given in Refs.$^{39,43}$.

### A. Collision-less flow

To model a collision-less ion flow, the ion velocity component $v_z$ along the electric field $\vec{E} = E_0\hat{z} = -1850\hat{z}$ V/m is sampled from a shifted Maxwellian distribution:

$$f_M(v_z) = \frac{1}{\sqrt{2\pi}v_T} \exp\left(-\frac{(v_z - v_{fl})^2}{2v_T^2}\right), \qquad (8)$$

where the thermal velocity is $v_{T_i} = \sqrt{k_B T_i / m_i}$, with $T_i$ the ion temperature (assumed equal to the neutral gas temperature $T_n = 300$ K) and $m_i$ the ion mass. The shift corresponds to a Mach number $M = v_{fl} / v_B \simeq 2.2$, where $v_B = \sqrt{k_B T_e / m_i}$ is the Bohm velocity, and the flow is directed downward along the electric field. The $x$ and $y$ velocity components are sampled from a non-shifted Maxwellian at the same temperature $T_i$. Ions that leave the simulation region or are absorbed by the dust are re-inserted on the boundaries following the algorithm described in$^{46}$.

Figure 2 shows the steady-state ion velocity distribution function (IVDF) after 15 µs of physical time. Panel (a) depicts the radial IVDF, which remains close to a
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) Two-dimensional histogram of radial ion velocity versus vertical position z, normalised to $v_{T_i}$ for collisionless flow. (b) Histogram of vertical ion velocity versus z, normalised to $v_B$. The dashed red line marks the position of the dust cluster. The insets are the velocity histograms at the dust cluster position.</figcaption></figure>

Maxwellian. Panel (b) shows the IVDF along z: due to the electric field and microparticles, the most probable velocity is near $3v_B$, with a secondary peak around $2v_B$ most probably due to ion re-injection. No significant IVDF changes occur above or below the dust cluster.

## B. Electric field driven collisional flow

Ion injection follows the same procedure as in the collision-less case, but ion-neutral collisions are now enabled. As a full model of the RF sheath is beyond the scope of this study, we assume a constant axial electric field over the simulation region ($\vec{E} = E_0\hat{z} = -1850\hat{z}$ V/m). Some codes, such as OpenDust, handle ion-neutral collisions assuming a constant collision frequency$^{52}$. Within the sheath region, charge transfer is the dominant collision outcome and the ion-neutral cross section is characterised by a weak velocity dependence$^{19}$. Assuming a constant collision frequency yields a velocity distribution qualitatively similar to that given by a constant-cross-section model$^{19}$.

In DRIAD, ion-neutral collisions are handled through the null-collision method$^{48,49,53,54}$ with collision cross sections taken from the LxCAT database. In this case, a subset of ions is selected for potential collisions based

on a reasonable maximum collision frequency. The collisions outcomes for the selected ions are then determined by a Monte Carlo method using a velocity-based collision probability$^{49,53}$. Note that the use of an equilibrium ion velocity distribution deduced from mobility data measured at constant applied electric field is not relevant for the sheath region of low-pressure discharges where the electric field is non-constant and the number of ion-neutral collisions is too small to reach equilibrium$^{55}$.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Same as Fig. 2 for an electric-field-driven collisional flow.</figcaption></figure>

Figure 3 shows the steady-state ion velocity distribution function (IVDF) after 15 µs of physical time. Panel (a) depicts the radial IVDF, which remains close to a Maxwellian, as in the collision-less case. Panel (b) shows the IVDF along z: the most probable velocity is near $2.2v_B$. The main difference is the presence of many slow ions produced by charge-exchange collisions. The vertical IVDF differs significantly from the equilibrium IVDF for an electric-field-driven collisional ion flow$^{19}$. This discrepancy is likely due to (i) the limited domain size, which reduces the number of collisions needed for equilibrium, and (ii) re-injected ions being sampled from a shifted Maxwellian. Again, no significant IVDF changes occur above or below the dust cluster.
<figure><img src="image_5.png" /><figcaption>FIG. 4. Ion–density and electric–potential axial profiles in the y–z plane for the collision-less case. (a,b) Single microparticle; (c,d) superposition of wakes from 19 isolated microparticles arranged as in the cluster; (e,f) full 19–particle cluster. The cluster generates a substantially stronger and more coherent wake structure than the simple linear summation of isolated wakes.</figcaption></figure>

## IV. RESULTS

### A. Cluster of microparticles in a collision-less flow

In Fig. 4, the axial ion–density profiles and the axial electric–potential profiles in the y–z plane are shown for three configurations in a collision-less flow with Mach number $M = 2.2$, a background plasma density $n_\infty = n_e = 8 \times 10^8$ cm$^{-3}$, an electron temperature $T_e = 2.15$ eV, and an applied electric field $E_0 = 1850$ V/m. Panels (a,b) display the profiles for a single grain. Panels (c,d) show the profiles obtained by summing the wakes of 19 isolated grains placed at the cluster positions. Panels (e,f) show the profiles for the full 19–grain cluster.

As seen in the figure, the ion–density and potential enhancements for 19 grains (whether summed or fully simulated) are markedly larger than for a single grain. However, the summed wakes do not reproduce the cluster results: the ion–density enhancement is stronger for the cluster, and the potential maximum is both larger and shifted downwards relative to the superposed case. In addition, the wakes of grains located at the cluster periphery are bent inwards towards the axis of the simulation domain.

The three–dimensional profiles for the cluster, shown in Fig. 5 for an extended domain, reveal that the max–

<figure><img src="image_6.png" /><figcaption>FIG. 5. Collision-less case: (a) three-dimensional ion–density distribution and (b) corresponding electric–potential distribution. The coherent enhancement observed downstream of the cluster develops into a well-defined structure distinct from the individual wakes.</figcaption></figure>

ima of the ion–density and electric–potential distributions coincide and are located at approximately $10\lambda_{D_e}$ downstream of the cluster. Unlike the individual wakes directly beneath each particle, this strong ion–density enhancement exhibits no distinct contributions from separate grains. Weak traces of individual wakes remain
visible, especially for the central grain and its near-
est neighbour, where small potential increases of order
~ 0.01 V appear and correlate with local ion-density en-
hancements (see Fig. 4f). These electric-potential inho-
mogeneities, produced by ion wakes, can strongly influ-
ence microparticle dynamics, though this effect is beyond
the scope of the present study.

The most prominent feature is the relatively large po-
tential increase of approximately 0.6 V in the region of
enhanced ion density located around 10$\lambda_{D_e}$ below the
cluster. This enhancement may explain observed modifi-
cations of sheath potentials in the presence of levitating
microparticle monolayers$^{25}$, and possibly the reported
bending of complex plasma monolayers$^{31}$. Further exper-
imental studies of the potential structure beneath mono-
layers would therefore be highly valuable.

<figure><img src="image_7.png" /></figure>

FIG. 6. Electric charge of each grain in the collision-less case. Colour indicates the magnitude of the time-averaged charge, and circle size represents the standard deviation. The near-uniformity of charge across the cluster confirms the negligible role of electron depletion.

Fig. 6 shows the time-averaged grain charge $Q$ and the charge standard deviation $\sigma_Q$ in the collision-less case. The central grain carries approximately $25500e^{-}\pm400e^{-}$, while grains at the periphery carry about $25900e^{-}\pm480e^{-}$ (surface potential of $\approx -8$ V). Thus, the reduction in charge for the central grain is small and consistent with enhanced ion density due to the focusing of the ions. The electron depletion around an individual grain due to the charging of the dust is negligible. For an intergrain spacing of $1\lambda_{D_e}$, the charge density within one Debye sphere is $Q/V_{\text{sphere}}$, where $V_{\text{sphere}} = (4/3)\pi\lambda_{D_e}^3$. For a typical grain charge of $2.5 \times 10^4 e^{-}$ and $\lambda_{D_e} = 385 \ \mu\text{m}$, this corresponds to a charge density of $\sim 1 \times 10^8 \ \text{cm}^{-3}$, eight times smaller than the background electron density $n_{\infty}$. Electron depletion could, however, become significant in much denser clusters.

B. Cluster of microparticles in a collisional flow

In Fig. 7, the ion-density profiles and the axial elec-
tric-potential profiles in the y-z plane within the simu-
lation domain are shown for three cases in an electric-

field-driven collisional ion flow in argon at a pressure
$p = 0.75$ Pa. The applied electric field is $E_0 = 1850$ V/m,
the neutral gas temperature is $T_n = 300$ K, the back-
ground plasma density is $n_\infty = n_e = 8 \times 10^8$ cm$^{-3}$,
and the electron temperature is $T_e = 2.15$ eV. These parame-
ters correspond to a flow Mach number of approximately
$M \approx 2.2$. Panels (a,b) show the profiles for a single grain.
Panels (c,d) show the profiles obtained by summing the
wakes of 19 isolated grains arranged as in the cluster.
Panels (e,f) show the profiles for the full 19-grain clus-
ter.

As shown in Figs. 7(a,c,e), the collisional flow produces enhanced ion density downstream of each grain. For the single-grain case, this behaviour agrees with earlier studies of electric-field-driven collisional flows$^{37,38,52}$. As in the collision-less case, the ion-density and potential enhancements for 19 grains (summed or clustered) are substantially larger than for a single grain (Figs. 7(b,d,f)). However, the summed wakes do not reproduce the full-cluster profiles: the ion-density enhancement is significantly stronger in the cluster, and the potential maximum, though of similar magnitude, is shifted upwards towards the grains. The wakes of grains at the cluster edge bend towards the simulation axis and are noticeably shorter than the central wake.

The density enhancement within each individual wake
is also much stronger than in the collision-less case, while
the wakes themselves are less elongated. This is consis-
tent with the known weakening of ion focusing due to
ion-neutral collisions⁵⁶. A striking feature of the collis-
sional regime is the pronounced cone-shaped region of
increased ion density that forms immediately below the
cluster and extends to approximately 6$\lambda_{D_e}$. This collec-
tive structure is absent when simply summing isolated
wakes. Although the individual wakes are still visible
near each grain, the dominant enhancement arises from
the cluster as a whole.

The three-dimensional profiles in Fig. 8 confirm the
formation of this cone-shaped density enhancement. In
Fig. 8(b), together with Fig. 7(f), it is clear that the
individual wakes are no longer distinguishable in an ex-
tended domain. Instead, a single collective wake–poten-
tial structure forms, with a maximum enhancement of
approximately 0.5 V along the central axis extending to
about 4$\lambda_{D_e}$. The potential enhancement is non-uniform
and peaks along the axis. This constitutes a potentially
measurable modification of the electric potential beneath
the monolayer.

Fig. 9 shows the time-averaged charge $Q$ and the stan-
dard deviation $\sigma_Q$ in the collisional case. The central
grain carries approximately $17200e^{-} \pm 240e^{-}$, whereas
peripheral grains carry much larger charges of about
$19000e^{-} \pm 200e^{-}$ (surface potential of $\approx -6$ V). Intermedi-
ate grains have charges around $18000e^{-}$, larger than that
of the central grain. The strong charge reduction towards
the cluster centre is fully consistent with the collective
ion-focusing effect generated by the cluster.
<figure><img src="image_8.png" /><figcaption>FIG. 7. Ion-density and electric-potential axial profiles in the y-z plane for the collisional case. (a,b) Single microparticle; (c,d) sum of wakes from 19 isolated microparticles; (e,f) full 19-particle cluster. The cluster produces a markedly stronger and more collective wake structure than the simple superposition of isolated wakes.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. Three-dimensional collisional case. (a) Ion-density distribution showing the extended cone-shaped enhancement below the cluster. (b) Corresponding electric-potential distribution, exhibiting a collective wake-potential structure rather than signatures of individual wakes.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 9. Electric charge of each particle in the collisional case. Colours represent the magnitude of the time-averaged charge, and circle radii denote the corresponding standard deviation. Strong charge reduction occurs towards the cluster centre, consistent with the collective ion-focusing effect.</figcaption></figure>

V. DISCUSSION AND CONCLUSION

In this work, we have examined ion-wake formation beneath two-dimensional microparticle clusters using the DRIAD molecular-dynamics code. Our results

demonstrate that the wake structures generated by a two-dimensional multi-particle clusters cannot be understood as a simple linear superposition of individual single-grain wakes. Instead, both the collision-less and collisional regimes exhibit a pronounced collective ion-density enhancement beneath the cluster, accompanied by a measurable increase in electric potential of up
to ~ 0.5 V.

The simulations also highlight that the morphol-
ogy, spatial extent, and intensity of the wakefield are
highly sensitive to the ion-velocity distribution func-
tion (IVDF). Collision-less flows produce elongated
wakes with maxima located more than 10λ<sub>De</sub> down-
stream, whereas collisional flows result in more com-
pact, cone-shaped structures extending to roughly 6λ<sub>De</sub>.
These differences underscore the importance of modelling
ion flows appropriate to specific plasma sheath conditions
and motivate further efforts aimed at implementing more
realistic, experimentally benchmarked IVDFs.

The potential increase reported here is small com-
pared with typical surface potentials of dust grains, yet
it arises in a region where the direct electrostatic influ-
ence of the grains themselves is negligible. This suggests
that collective wakefields may contribute to phenomena
such as sheath modification²⁵ and monolayer bending³¹,
as observed experimentally. The differing charge distri-
butions between central and peripheral grains, and in
particular the significant charge reduction towards the
centre of the cluster in the collisional regime, are con-
sistent with collective ion focusing and may influence
mode-coupling processes, including the mode-coupling
instability in two-dimensional crystals. These findings in-
dicate that wake models assuming point-like or isolated
wakes are inadequate for describing interparticle forces
in extended monolayers, and that collective wakefields
should be incorporated into future kinetic and fluid-level
descriptions of complex-plasma stability.

At the same time, the present study highlights several
areas where improved modelling is needed. The strong
dependence on the ion-velocity distribution and the sim-
plifying assumption of a constant electric field in the
sheath motivate more advanced treatments that incor-
porate spatially varying electric fields, realistic boundary
conditions, and the known RF-induced modulation of the
IVDF⁵⁷. Similarly, while the asymmetric MD approach
employed here is physically justified within a dust-cluster
context, systematic benchmarking against PIC simula-
tions in large domains would strengthen confidence in the
validity of the method for long-range structures. Further
investigations of domain-size convergence, charging mod-
els in dense clusters, and the influence of finite electron
depletion would also help to assess the generality of the
conclusions.

Overall, the results presented here provide clear evi-
dence that collective ion wakes are an intrinsic feature
of dust monolayers and extended clusters in flowing plas-
mas. Future work should combine improved kinetic mod-
elling with targeted laser-induced fluorescence measure-
ments capable of resolving ion-density and potential per-
turbations beneath microparticle layers. Such studies
will help determine the extent to which collective wake-
fields contribute to crystal stability, non-reciprocal forces,
and mode-coupling phenomena in two-dimensional com-
plex plasma crystals.

ACKNOWLEDGMENTS

LC acknowledges support from the Natural Sciences
and Engineering Research Council of Canada (NSERC),
Grant No. RGPIN-2025-05757, from the ANR –
FRANCE (French National Research Agency) for its
financial support of the SARDINE project ANR-23-
CE30-0025 and from the french government under the
France 2030 investment plan, as part of the Initiative
d'Excellence d'Aix-Marseille Université – A*MIDEX ”
AMX-22-RE-AB-034. LSM, DJM, and BRS gratefully
acknowledge support from the US Department of Energy,
Office of Science, Office of Fusion Energy Sciences under
award number DE-SC-0021334 and the National Science
Foundation grant PHY 2308743.

¹P. K. Shukla, Physics of Plasmas **8**, 1791 (2001).

²O. Havnes and G. E. Morfill, Adv. Space Res. **4**, 85 (1984).

³M. Horanyi and D. Mendis, The Astrophysical Journal **294**, 357 (1985).

⁴T. W. Hartquist, W. Pilipp, and O. Havnes, Astrophys. Space Sci. **246**, 243 (1996).

5J. Beckers, J. Berndt, D. Block, M. Bonitz, P. J. Brugge-
man, L. Couëdel, G. L. Delzanno, Y. Feng, R. Gopalakrishnan,
F. Greiner, P. Hartmann, M. Horányi, H. Kersten, C. A. Knapek,
U. Konopka, U. Kortshagen, E. G. Kostadinova, E. Kovačević,
S. I. Krasheninnikov, I. Mann, D. Mariotti, L. S. Matthews,
A. Melzer, M. Mikikian, V. Nosenko, M. Y. Pustylnik, S. Ratyn-
skaia, R. M. Sankaran, V. Schneider, E. J. Thimsen, E. Thomas,
H. M. Thomas, P. Tolias, and M. van de Kerkhof, Physics of
Plasmas **30** (2023), 10.1063/5.0168088.

⁶G. S. Selwyn, J. Singh, and R. S. Bennett, J. Vac. Sci. Technol. A **7**, 2758 (1989).

⁷R. L. Merlino and J. A. Goree, PHYSICS TODAY. **57**, 32 (2004).

⁸D. Samsonov, J. Goree, H. M. Thomas, and G. E. Morfill, Phys. Rev. E **61**, 5557 (2000).

⁹V. Nosenko, S. Zhdanov, A. V. Ivlev, G. Morfill, J. Goree, and A. Piel, Phys. Rev. Lett. **100**, 025003 (2008).

¹⁰L. Couëdel, V. Nosenko, S. K. Zhdanov, A. V. Ivlev, H. M. Thomas, and G. E. Morfill, Phys. Rev. Lett. **103**, 215001 (2009).

¹¹A. V. Zampetaki, H. Huang, C.-R. Du, H. Löwen, and A. V. Ivlev, Phys. Rev. E **102**, 043204 (2020).

¹²G. E. Morfill and A. V. Ivlev, Rev. Mod. Phys. **81**, 1353 (2009).

¹³H. M. Thomas and G. E. Morfill, Nature **379**, 806 (1996).

¹⁴S. Nunomura, D. Samsonov, and J. Goree, Phys. Rev. Lett. **84**, 5141 (2000).

¹⁵J. Goree, B. Liu, and Y. Feng, Plasma Phys. Control. Fus. **55**, 124004 (2013).

¹⁶F. Melandsø and J. Goree, Phys. Rev. E **52**, 5312 (1995).

¹⁷O. Ishihara and S. V. Vladimirov, Phys. Plasmas **4**, 69 (1997).

¹⁸M. Lampe, G. Joyce, G. Ganguli, and V. Gavrishchaka, Phys. Plasmas **7**, 3851 (2000).

¹⁹R. Kompaneets, G. E. Morfill, and A. V. Ivlev, Phys. Rev. E **93**, 063201 (2016).

²⁰L.-J. Hou, Y.-N. Wang, and Z. L. Mišković, Phys. Rev. E **64**, 046406 (2001).

²¹R. Kompaneets, U. Konopka, A. V. Ivlev, V. Tsytovich, and G. Morfill, Phys. Plasmas **14**, 052108 (2007).

²²V. Steinberg, R. Sütterlin, A. V. Ivlev, and G. Morfill, Phys. Rev. Lett. **86**, 4540 (2001).

²³A. Samarian, S. Vladimirov, and B. James, JETP Lett. **82**, 758 (2005).

²⁴A. Melzer, V. A. Schweigert, and A. Piel, Physica Scripta **61**, 494 (2000).

²⁵C. Arnas, M. Mikikian, G. Bachet, and F. Doveil, Phys. Plasmas **7**, 4418 (2000).
²⁶S. K. Zhdanov, A. V. Ivlev, and G. E. Morfill, Phys. Plasmas **16**, 083706 (2009).

²⁷L. Couëdel, V. Nosenko, A. V. Ivlev, S. K. Zhdanov, H. M. Thomas, and G. E. Morfill, Phys. Rev. Lett. **104**, 195001 (2010).

²⁸L. Couëdel, S. K. Zhdanov, A. V. Ivlev, V. Nosenko, H. M. Thomas, and G. E. Morfill, Phys. Plasmas **18**, 083707 (2011).

²⁹A. V. Ivlev, S. K. Zhdanov, M. Lampe, and G. E. Morfill, Phys. Rev. Lett. **113**, 135002 (2014).

³⁰T. B. Röcker, A. V. Ivlev, S. K. Zhdanov, and G. E. Morfill, Phys. Rev. E **89**, 013104 (2014).

³¹T. B. Röcker, A. V. Ivlev, S. K. Zhdanov, L. Couëdel, and G. E. Morfill, Phys. Plasmas **21**, 073711 (2014).

³²A. V. Ivlev and G. Morfill, Phys. Rev. E **63**, 016409 (2000).

³³V. V. Yaroshenko, A. V. Ivlev, and G. E. Morfill, Phys. Rev. E **71**, 046405 (2005).

³⁴T. B. Röcker, A. V. Ivlev, R. Kompaneets, and G. E. Morfill, Phys. Plasmas **19**, 033708 (2012).

³⁵A. Melzer, Phys. Rev. E **90**, 053103 (2014).

³⁶A. V. Ivlev and R. Kompaneets, Phys. Rev. E **95**, 053202 (2017).

³⁷L. S. Matthews, D. L. Sanford, E. G. Kostadinova, K. S. Ashrafi, E. Guay, and T. W. Hyde, Physics of Plasmas **27**, 023703 (2020).

³⁸R. Banka, K. Vermillion, L. Matthews, T. Hyde, and L. Couëdel, Plasma Physics and Controlled Fusion **65**, 044006 (2023).

³⁹A. Piel, Phys. Plasmas **24**, 033712 (2017).

⁴⁰K. Vermillion, R. Banka, A. Mendoza, B. Wyatt, L. Matthews, and T. Hyde, Physics of Plasmas **31** (2024), 10.1063/5.0203902.

⁴¹Z. Watson, P. Adamson, J. Martinez-Ortiz, K. Vermillion, C. Carmichael, S. Garcia-Rodriguez, L. Matthews, T. Hyde, and B. Wyatt, Fundamental Plasma Physics **13**, 100082 (2025).

⁴²K. S. Ashrafi, R. Yousefi, M. Chen, L. S. Matthews, and T. W. Hyde, Phys. Rev. E **102**, 043210 (2020).

⁴³K. S. Ashrafi, *Investigations of the Plasma Conditions in the Sheath Using Dust Grain as Probes*, Ph.D. thesis, Baylor University (2020).

⁴⁴K. Vermillion, D. Sanford, L. Matthews, P. Hartmann, M. Rosenberg, E. Kostadinova, J. Carmona-Reyes, T. Hyde, A. M. Lipaev, A. D. Usachev, A. V. Zobnin, O. F. Petrov, M. H. Thoma, M. Y. Pustynik, H. M. Thomas, and A. Ovchinin, Physics of Plasmas **29**, 023701 (2022).

⁴⁵L. Matthews, K. Vermillion, P. Hartmann, M. Rosenberg, S. Rostami, E. Kostadinova, T. Hyde, M. Pustynik, A. Lipaev, A. Usachev, A. Zobnin, M. Thoma, O. Petrov, H. Thomas, and O. Novitskiy, Journal of Plasma Physics **87** (2021), 10.1017/S0022377821001215.

⁴⁶I. H. Hutchinson, Plasma Physics and Controlled Fusion **44**, 1953 (2002).

⁴⁷A. V. Phelps, J. Appl. Phys. **76**, 747 (1994).

⁴⁸H. R. Skullerud, Journal of Physics D: Applied Physics **1**, 1567 (1968).

⁴⁹Z. Donkó, Plasma Sources Science and Technology **20**, 024001 (2011).

⁵⁰L. Couëdel and V. Nosenko, Phys. Rev. E **105**, 015210 (2022).

⁵¹J. Kong, T. W. Hyde, L. Matthews, K. Qiao, Z. Zhang, and A. Douglass, Physical Review E **84**, 016411 (2011).

⁵²D. Kolotinskii and A. Timofeev, Computer Physics Communications **288**, 108746 (2023).

⁵³Z. Donko, A. Derzsi, M. Vass, B. Horváth, S. Wilczek, B. Hartmann, and P. Hartmann, Plasma Sources Science and Technology **30**, 095017 (2021).

⁵⁴K. Nanbu, IEEE Transactions on Plasma Science **28**, 971 (2000).

⁵⁵J.-H. Mun, M. Muraglia, O. Agullo, C. Arnas, and L. Couëdel, Physics of Plasmas **31**, 073906 (2024).

⁵⁶A. Piel, *Plasma physics: an introduction to laboratory, space, and fusion plasmas*, edited by Springer (Springer, 2017).

⁵⁷A. D. Kuypers and H. J. Hopman, Journal of Applied Physics **63**, 1894 (1988).