Fractional Laplacian spectral approach to turbulence in a dusty plasma monolayer

Chaos and nonlinearity are ubiquitous in the natural world, yet the current understanding of many related phenomena is limited by the lack of solid theoretical framework. An outstanding example is the origin of turbulence in charged media. Unlike conventional fluids, where turbulent behavior is typically dominated by inertial forces, the excitation of instabilities in charged fluids is further influenced by electromagnetic forces and collective effects.¹ In addition, dynamical instabilities, such as vortices, can lead to a redistribution of charge throughout the fluid, resulting in modification of the electric field. These effects are particularly prominent in systems where inertial (heavy) particles are suspended within the flow.^2,3 Such systems include atmospheric clouds,^4–6 fluidized bed reactors,^7–9 charged industrial sprays,¹⁰ and climate-related dust lifting processes.^11,12 Thus, knowledge of the fundamental physics mechanisms underlying the coupling between charge-related phenomena and the onset of instabilities is key to understanding turbulence across a wide range of physical systems.

Investigation of such mechanisms is especially interesting at low-to-medium Reynolds numbers (⁠$Re≲100$⁠), where the fluid dynamics are not solely dominated by inertial effects. The onset of turbulence in this regime is important in the study of highly viscous liquids or media with small-scale flows, including viscoelastic fluids,^13–16 porous materials,^17–19 and airfoil devices.^20–22 The ubiquitous presence of charge inhomogeneities and inertial particle suspensions in these systems requires knowledge of the coupling between charge-mediated interactions and small-scale excitations. Another interesting question is turbulence in active systems, such as insect flights^23,24 and stochastic self-propelled devices,^25,26 where energy driven at small scales can be transported to larger scales, causing a global instability. As these various systems are often characterized by small-scale, high-frequency dynamics, which violate the assumptions needed for a hydrodynamic approximation, investigation of the onset of turbulence in this case requires a kinetic treatment.

The complexity of phenomena affecting the onset of turbulence in charged media and/or active systems can be studied at the kinetic level employing complex (dusty) plasma as a model system. The field of complex plasmas investigates the dynamics of mesoscopic particles (or dust) suspended in weakly ionized low temperature gas. Dust grains immersed in plasma become negatively charged and are subject to both ion drag forces and collective interactions. Dusty plasmas can self-organize into strongly coupled fluids and crystalline structures,^27–30 which makes them ideal for the study of self-organization and stability, phase transitions, and transport phenomena. Experiments and numerical simulations using dusty liquids at low-to-medium Reynolds numbers have already been used to analyze self-excited^31–33 and externally driven^34,35 vortices, shear instabilities,^36,37 Kolmogorov flows,³⁸ and wave turbulence.^39–41 More importantly, the mesoscopic particles, which are the tracers of the turbulent dynamics in a dusty plasma liquid, are directly observable on easily accessible spatial and temporal scales, which allows an investigation at the kinetic level of the connection between global instabilities and individual particle statistics.

Here, we study the onset of turbulence in a dusty plasma monolayer, where energy is driven at the level of each individual particle and transported to larger scales, causing an “inverse” cascade, as opposed to the classical Kolmogorov cascade.⁴² As an inverse cascade is typical for two-dimensional (2D) turbulence, it has been observed in active matter monolayer systems, such as bacterial suspensions⁴³ and self-propelled camphor swimmers.⁴⁴ Dusty plasma monolayers, or quasi-2D structures, are intrinsically non-equilibrium systems, where energy is constantly driven and dissipated at small scales due to the interaction of individual dust particles with the neutral gas and plasma environment. Additionally, dust charging processes in these systems yield both nonlocal interactions and stochastic fluctuations in the dust particle spatial distribution. Thus, the investigation of self-excited instabilities in dusty plasma can provide insight into the connection between charge-driven processes and the inverse energy cascade in 2D turbulence.

This paper presents many-body (MB) simulations where dusty plasma monolayers are formed in conditions relevant to laboratory experiments. Specifically, we assume a distribution of dust particle sizes consistent with values reported by dust particle suppliers. Deviations from the mean particle size result in fluctuations of the corresponding particle charge, interparticle spacing, and interaction potential, all of which determine the dusty plasma dynamics and stability. To mimic conditions achievable in dusty plasma experiments in large electrode chambers,⁴⁵ we consider 10 000-particle monolayers with radial confinement acting most strongly at the cloud exterior, as set by a tenth-order polynomial confinement force. The kinetic energy in the system is regulated by the gas pressure: collisions between dust and neutral gas particles provide both a drag force, which damps particle motion, as well as a thermal bath, which provides “kick” excitations of various strength and direction to each dust grain.

Comparisons between stable and unstable monolayers are made by considering three pressure regimes, $p=5, 1, and 0.1 Pa$⁠, which allows for simulation of both crystalline and liquid-like monolayers. At $5 Pa$⁠, a stable monolayer is formed with distinct hexagonal symmetry and large crystalline domains. At $1 Pa$⁠, both the gas drag force and the strength of kicks from the thermal bath are decreased and the dust particles are mobilized, which results in enhancement of defect lines and shrinking of the crystalline regions. As the pressure drops to $0.1 Pa$⁠, the crystalline domains are destroyed, and the monolayer develops an instability. In each case, energy is imported at the individual particle level by the thermal bath and transported across larger spatial scales through electrostatic scattering events, which are enhanced at low pressure due to increased particle mobility. This allows for an inverse energy cascade from small to large scales for given system parameters. The direct Kolmogorovian energy cascade from large to small scales is also possible in dusty plasmas, for example, using electron beams or optical lasers to induce a large-scale shear flow in the structure.^46,47 This concept will be explored in our future work. The connection between 2D turbulence in dusty plasma monolayers and inertial turbulence is briefly discussed in Sec. III A.

For each set of conditions, the observed dust dynamics is compared against the predictions from the Fractional Laplacian Spectral (FLS) model, which determines the time-evolved energy state of the systems solely from knowledge of the underlying Hamiltonian structure. The spectral approach employed here (introduced in Refs. 48 and 49) identifies whether the energy spectrum of a given Hamiltonian operator includes an absolutely continuous part, which indicates the existence of transport in the form of scattering energy states. To explore the effect of random disorder and nonlocal interactions on the energy transport, the Hamiltonian used in this work is the random fractional discrete Schrödinger operator.⁵⁰ The potential energy term of this operator models random disorder, while nonlocality is allowed by employing a fractional Laplacian operator for the kinetic energy term. However, the spectral approach is a general mathematical tool, which can be applied to any Anderson-type Hamiltonian of the form defined in Ref. 51.

Scaling between the many-body (MB) simulation and the FLS technique is performed in the following manner. In the MB simulation, the variation of dust size (and, therefore, mass and charge) and properties of the radial confinement result in fluctuations of the mean interparticle separation and interparticle potential in space. The standard deviation from the mean particle spacing is used as a measure of random disorder, which determines the corresponding potential energy term in the Hamiltonian $H$ used in the FLS calculation. The kinetic energy term in $H$ is obtained from the characteristics of the dust particle diffusion observed in the MB simulation (the plot of mean squared displacement as a function of time and the probability distribution function of particle velocities). Once a Hamiltonian for a given set of MB simulation conditions is determined, the FLS technique is used to identify the existence of transport as a function of spatial scales, solely from the properties of the corresponding energy spectrum. Thus, an agreement between the (kinetic) many-body simulations and the (spectral) FLS calculations provides a powerful modeling tool for the study of energy transport in turbulent dynamics.

In Sec. II, we introduce the main features of the many-body simulation and discuss how realistic experimental conditions are scaled as input parameters for the simulation. Structural and statistical analyses of the resulting dusty plasma monolayers are presented in Sec. III, where we also discuss scaling between the observed properties of the monolayers and the corresponding Hamiltonian structure. In Sec. IV, the FLS method is introduced and used to compute the probability for energy transport as a function of scale for each set of conditions. Discussion of the observed similarities between the kinetic simulation and the spectral model is provided in Sec. V. In Sec. VI, we summarize the results and outline directions for future research. More details of the many-body simulation are discussed in  Appendix A. Detailed description of the FLS method and related mathematical proofs can be found in Ref. 50.  Appendix B includes further details on the velocity distribution functions obtained from the many-body simulations.

An outstanding challenge in the experimental study of global instabilities in complex plasmas on Earth is the need to generate and sustain extended structures, consisting of many dust particles. A main motivation for the present numerical study is the recent development of dusty plasma experiments where extended monolayers, consisting of more than 10 000 particles, have been achieved. One such experimental setup is cell 3 [Fig. 1(a)], located in Baylor University's CASPER lab. Cell 3 is based on the GEC RF reference cell chamber,⁵² with several modifications to the original design required for suspending dust particles in the discharge and visualizing them using different optical systems.^53–55 A unique feature of cell 3 is its large geometrical size, with an 8-inch (⁠$≈20.3$ cm) diameter lower electrode, which allows for the formation of very large dust clouds. Experiments with cell 3 have already produced dust crystals (monolayers and multi-layered structures), with radial spread of about 7.5 in (⁠$≈19$ cm). Figure 1(b) shows the outer region of such a dusty plasma monolayer suspended in argon plasma at gas pressure $1.33 Pa$⁠. At this pressure, the structure exhibits both crystalline domains and liquid-like disordered regions, which suggests the possibility of interesting energy transport throughout the structure.

To the best of our knowledge, dust structures of comparable diameter have been reported in only one other modified GEC RF cell within the dusty plasma community, as described by Meyer et al.⁴⁵ The structures discussed by Meyer et al. consist of $≈10 000$ dust particles, forming monolayers of diameter $≈27 cm$⁠. At $0.15 Pa$⁠, these extended monolayers were observed to exhibit wave instabilities and mode coupling, which suggests that this lower pressure is appropriate for modeling fluid dynamics phenomena, such as turbulence.

Motivated by the experimental realization of extended dusty plasma structures in large chambers, here we focus on modeling monolayers consisting of $10 000$ particles at pressures $1$ and $0.1 Pa$⁠; conditions for which a liquid-like dynamics have been observed. For completeness and comparison to more ordered structures, we also consider a case where the monolayer is generated at $5 Pa$⁠. At this higher pressure, neutral gas drag stabilizes the dust particle motion and the resulting monolayers are expected to exhibit large-scale crystalline domains and well-defined hexagonal symmetry (e.g., see the experiments reported in Ref. 30). Table I summarizes the experimentally relevant conditions used as inputs for the many-body simulations of dusty plasma monolayers at the three considered pressure cases.

In Table I, the dust charge $Qd$ is in units of elementary charge $e−$⁠, $λD$ is the Debye screening length due to the electron and ions in the plasma, and $Ωh$ is the radial confinement frequency (discussed in Sec. II B). For the low pressure cases, $p=1 , 0.1 Pa$⁠, the selected Debye length is $λD=1 mm$⁠, which has been reported as a typical value at lower pressures in Argon RF glow discharge.^45,56 In Ref. 45, the charge on the dust grains was estimated to $Qd=24 648 e−$ which accounts for the reduction of charge due to interactions with ion wakefields. At $p=5 Pa$⁠, the higher plasma density leads to a reduced Debye length of $λD=0.6 mm$⁠, as discussed in Ref. 30. Since dust charge reduction due to ion-neutral collisions is also expected, in this case the simulation employs a dust charge $Qd=18 570 e−$⁠, as suggested in Ref. 57. In Table I, the values of $Qd$ represent the mean charge on the dust particles. However, dust particles used in laboratory experiments are known to exhibit slight variations from the mean particle diameter, which can lead to observable differences in the monolayer dynamics, such as the bistability reported in Ref. 56. Thus, here the dust particle diameters are selected from a normal distribution with a mean of $9.19 μm$ and a standard deviation of $±0.9%$⁠. This also leads to variation in the corresponding mass and charge on each dust grain.

To simulate the dynamics of dusty plasma monolayers, we employ an $N$-body code optimized to include specified external potentials suitable for modeling the confinement forces and external electric fields characteristic of laboratory dusty plasma experiments. The code explicitly simulates the dynamics of the dust particles. The presence of other species in a dusty plasma (electrons, ions, and neutral atoms) is implicitly accounted for by including the assumed effect of these species on the dynamics of the dust. The code allows for specification of dust particle parameters, such as size, charge, and material density. In this work, the interaction between the dust grains and the surrounding plasma environment is simulated by choosing the appropriate Debye shielding length and form of the interaction potential. This code has been extensively used in modeling the dynamics of charged dust in astrophysical environments^58–62 and in GEC RF reference cells in Earth-based experiments.^63–66 In all cases, we use a simulation box size with $85 cm$ sides, which was selected to mimic the $85 cm$-diameter electrode chamber used for the experiments in Ref. 45.

In all simulations, we use a fourth-order Runge–Kutta adaptive time step with minimum time step $10−6 s$ and maximum time step $10−4 s$⁠. As a rule of thumb, the selected time step needs to be small enough to resolve the particle motion on a meaningful scale for the velocities involved. As discussed in Sec. III (see Fig. 2), for the 5 and 1 Pa cases, the horizontal and vertical velocities are $∼10−5 ms−1$⁠, while in the 0.1 Pa case, the vertical velocity for individual particles is as large as $∼10−2 ms−1$⁠. Therefore, in all simulations, the maximum particle motion is $∼10−4 s ×10−5 ms−1=10−9 m$ for the 0.1 and 5 Pa cases, and $∼10−4 s ×10−2 ms−1=10−6 m$ for the 0.1 Pa case. As the dust size is $≈10 μm=10−5 m$⁠, the selected time step can well resolve the dust motion.

In the many-body simulation, the equation of motion for a single dust particle with mass $md$ and charge $Qd$ is given by

where $Fdd$ is the force between pairs of dust grains, $mdg$ is the gravitational force, $QdE$ is the confining electric field force, and $βẋ$ is the neutral drag force. The force between pairs of dust grains $Fdd$ is assumed to be a shielded Coulomb (Yukawa) interaction, where the screening is mainly provided by the ions and determined from the Debye length $λD$⁠. The last term $ζrt$ is a thermal bath, which simulates diffusive motion of the dust grains due to random collisions with the neutrals in the gas. Here, $ζ$ is the maximum acceleration exerted from the neutrals to each dust particle and $rt$ is a random number selected from a normal distribution. Further details on the dust-dust interaction force, the axial (along the direction of gravity) electrostatic confinement force, neutral drag force, and thermal bath can be found in  Appendix A.

To model the radial electrostatic confinement in large-electrode cells, the radial potential acting on the $ith$ particle is a tenth-order polynomial of the form suggested in Ref. 45,

where $Ωh$ is the horizontal confinement frequency, $ρi$ is the radial position of the ith particle in the $xy$-plane, and $R$ is the approximate horizontal radius of the monolayer (radial extent). Following Meyer et al.,⁴⁵ we assume $Ωh=2π×0.12 s−1$ and $R=63 mm$⁠. Data for the particle positions, velocities, and accelerations are output every 0.001 s, which is selected to match a frame rate of $1000 fps$⁠, typically needed to resolve the motion of fast particles in laboratory dusty plasma experiments. Thus, in the following discussion, timesteps in the many-body simulation are sometimes referred to as frames.

Using experimentally relevant conditions from Table I, simulations of extended dusty plasma monolayers were conducted for the three pressure cases, $p=5, 1, and, 0.1 Pa$⁠, which are expected to yield different structures and dynamics. Figures 2(a)–2(c) show the dust particle positions after 30s of simulation time, which is sufficient to allow for damping of instabilities due to initial conditions. Comparison of the three pressure cases clearly demonstrates that the monolayers form a stable structure at $5$ and $1 Pa$⁠, while the $0.1 Pa$ case yields a less ordered, liquid-like structure. In each simulation, the properties of the radial confinement in Eq. (2) are fixed, which yields the following common structural feature: the force due to the tenth-order polynomial potential almost vanishes in the central region of each monolayer, while beyond a radial distance $≈63 mm$ [the value of $R$ in (2)], the confinement force grows rapidly. As a result, particles with radial positions $≲63 mm$ are mostly confined by an external layer of particles located at radial positions $≳63 mm$⁠. As the internal structure in each case is more ordered and less dependent on the strong boundary conditions, in Secs. III A–III C and Sec. IV the diffusion and spectral analysis was performed only using the particles located within a region with radius $Rc=50 mm$⁠, indicated by red circles in Figs. 2(a)–2(c).

The stability of the system at each pressure is visible from the time evolution of the radial and axial velocities, as shown in Figs. 2(d)–2(f). For both the $5$ and $1 Pa$ cases, after 30 s of simulation time, the horizontal velocities limit to $Vx≈Vy∼10−5 m/s$⁠, while the vertical velocities approach $Vz∼10−6 m/s$⁠. The $0.1 Pa$ case exhibits distinct behavior as particle velocities experience a rapid increase to values $∼10−4 m/s$⁠, with large fluctuation in both the vertical and horizontal directions, suggesting the onset of an instability. As shown in Fig. 2(f), the vertical velocity increases steadily for the first 14 s, followed by a rapid growth which begins to saturate in the last 5 s. The growth in the horizontal velocities occurs $≈5 s$ after the beginning of the growth in $Vz$⁠. The onset of this instability results from the interplay between high dust mobility due to the low pressure and out-of-phase dust oscillations in the $z$-direction due to the assumed variations in dust particle size. A similar phenomenon has been observed experimentally by Gogia and Burton,⁵⁶ where a dusty plasma monolayer exhibited bistability at low pressure, which was attributed to dust mass variation of similar magnitude as the variation assumed in our simulation.

Dust grains in laboratory dusty plasma monolayers tend to self-organize into lattices of triangular symmetry, i.e., on average, each dust particle is located in the center of a hexagonal cell with six nearest neighbors located at each vertex. In Sec. II, we pointed out that the present simulations of dusty plasma monolayers assume particle diameters that can vary from the mean by $±0.9%$⁠, which is typical for laboratory dusty plasma experiments. This size variation leads to variation of both dust mass and charge, which in turn, affects the spatial distribution of particles and introduces lattice defects throughout the monolayer. Two-dimensional (2D) lattice defects can be defined as the fraction of particles with the number of nearest neighbors $NN$ different than six. Commonly, these are particles with five or seven nearest neighbors, located around defect lines. To quantify the amount of lattice defect in each considered monolayer, we employ a crystallinity code by Hartmann et al.,³⁰ which uses a Delaunay triangulation algorithm to calculate the number of nearest neighbors for each dust grain and the complex bond-order parameter

where $NNi$ is the number of nearest neighbors of the ith particle and $Θil$ is the angle of the lth nearest-neighbor bond measured with respect to the $X$-axis. When the code is applied to dust fluid structures with badly defined primitive cells, the Delaunay triangulation function returns an error due to the insufficient number of unique points or overlapping particles. In other words, in these structures, the function encounters numerous points lying on the same line, in which case the triangulation does not exist. This is used to distinguish between monolayers with a well-defined crystalline structure and monolayers with a liquid-like structure.

Figure 3 shows the complex argument $ArgG6$⁠, which measures the overall angular orientation of a given particle neighborhood. The three plots in Fig. 3 are obtained from dust particle positions at $t=5 s$⁠, at which time the $5$ and $1 Pa$ cases have reached equilibrium and the instability has not yet developed in the lowest pressure case (see Fig. 2). In each plot, dust particles with six nearest neighbors are marked by black dots, while dust particles with seven or more (five or less) are marked by circles (triangles). The plots of the complex argument show well-defined crystalline domains in the $5$ and $1 Pa$ cases, while a long-distance order is lost in the $0.1 Pa$ realization. In addition, the fraction of particles with $NN ≠6$ increases as pressure decreases, yielding high defect concentration at $0.1 Pa$⁠.

Since the monolayer at $0.1 Pa$ does not preserve well-defined hexagonal cells throughout the simulation, it is more appropriate to define a one-dimensional (1D) disorder of the form

where $σa$ is the standard deviation of the average interparticle separation $a$ computed for all dust particles within the region of interest (red circle of radius $Rc=50 mm$ shown in Fig. 2). We note that typical interparticle separations in this central region are $≈2 mm$ in each pressure case, which is larger than the separation of $≈1.7 mm$ obtained when all particles in the monolayer are accounted for. These interparticle separations are comparable to those measured for extended dusty plasma monolayers in Baylor University's cell 3 (Fig. 1). Table II lists the 1D disorder concentration values $c$ for dust particles in the region of interest for each pressure case. These values are used as input for the Fractional Laplacian Spectral (FLS) method discussed in Sec. IV.

In this work, although we consider passive dust particles, we note that meaningful parallels can be drawn between the features of turbulence observed in our simulations at pressure $0.1 Pa$ and turbulence in active matter. Active matter is commonly defined as a collection of active particles, each of which can convert the energy coming from their environment into directed motion, therefore driving the whole system far from equilibrium.^67,68 Active organisms tend to self-organize and develop coherent collective behavior/motions. Such systems share the feature of being intrinsically out of equilibrium as energy is constantly injected at the level of each individual entity. Active systems do not have to consist of living organisms (e.g., systems of self-propelled objects, such as camphor disks studied in Ref. 44) All of these features are characteristic of dusty plasmas, which are non-equilibrium, driven, and dissipative, and exhibit self-organization and long-distance interactions. The main difference from classical self-propelled particles is that typically dust particles in dusty plasma interact with each other via electrostatic forces/scattering events and “respond” to gradients in the electrostatic potential in the environment. We note that, beyond the parallels drawn in the present discussion, the dynamics of active systems can be directly studied in dusty plasmas, where the particles are self-propelled by “rocket force” due to material ablation,⁶⁹ or dusty plasmas with Janus active particles.⁷⁰ Simulations of such active dusty plasmas are beyond the scope of the present work but represent and interesting direction for future work.

A main characteristic of 2D turbulence and turbulence in active matter is the presence of large-scale correlations. In Ref. 44, the transition from uncorrelated to interacting regime was achieved by increasing the number of active swimmers in the system to $Np=20–60$⁠, while the system was found to freeze for $Np≳60$⁠. In our simulation of dusty plasmas, the charged dust grains are interacting at $p=0.1 Pa$⁠, while the system is observed to freeze at $p≳1 Pa$⁠. The transition from an uncorrelated to a correlated regime is visible in the mean squared displacement of the particles plotted as a function of time delay, which develops a linear (diffusive) region between the ballistic region (short time scales) and the region where the MSD oscillates until it approaches a constant value (long time scales). Such dynamics are characteristic of passive particles following Lagrangian dynamics and for particle tracers in inertial turbulence, as discussed in Ref. 44. This feature can be seen in Fig. 4(a), where quadratic (purple line) and linear (red line) functions were fitted to the $MSD(τ)$ plot obtained for particles in the $0.1 Pa$ simulation. The quadratic function yields the best fit to the data in the region $τ∈0, 3 s$⁠, while the linear function yields the best fit in the intermediate region $τ∈3, 10 s$⁠. For $τ≳10 s$⁠, the MSD plot deviates from a linear function fit and oscillates. In the present work, instead of dividing the MSD plot into individual regions and treating the dynamics in each region separately, we fit the entire curve with a power law [blue line in Fig. 4(a)], which allows us to treat the dynamics as anomalous diffusion, as further discussed in Sec. III C. However, the presence of these regions within the MSD plot is clearly visible in Fig. 4(a), which suggests transition from uncorrelated to correlated regime in the $0.1 Pa$ case for time delays $τ≳10 s$⁠.

Another feature, commonly reported in 2D fluid turbulence,^71–73 is that the PDFs of velocities in the 2D plane are observed to deviate from a Gaussian. Figure 4(b) shows the histogram of velocities in the $x$-direction (the histogram for $vy$ has a similar shape) for the $0.1 Pa$ simulation, which is found to deviate from a Gaussian, and is best described by a modified Maxwellian with Kappa distribution tails. The parameter $κ$ characterizes how far a system is from thermal equilibrium. When $κ→∞$⁠, a Kappa distribution approaches a Maxwellian, but when $κ$ is finite, the distribution function has high-energy tails, suggesting an increased number of high-velocity particles. For $κ<10$⁠, the Kappa distribution has a power law tail.⁷⁴ Here, the best fit parameter yields $κ≈11$ for the PDFs of both $vx$ and $vy$⁠, suggesting that the 0.1 Pa case is not at thermal equilibrium but neither does it exhibit power-law tails. Further discussion on the velocity distributions fits for each case, and the relation to dust temperature is included in  Appendix B.

As discussed at the beginning of Sec. II, an instability occurs in the dusty plasma monolayer at $p=0.1 Pa$ [Fig. 2(f)], which suggests the possibility of turbulent dynamics. As this instability develops at later simulation times (⁠$t≳14 s$⁠), it is not attributed to fluctuations due to initial conditions. It should be further noticed that the instability is self-excited since conditions remain unchanged during the examined time period. In the simulation, energy is imported at the individual particle level by the thermal bath and transported across larger spatial scales through dust-dust electrostatic scattering events. The increased particle mobility at low pressures allows for both enhanced dust oscillations in the vertical direction and enhanced transport in the 2D plane, which facilitates the electrostatic scattering events. When deflected from their equilibrium positions, dust particles in dusty plasmas convert electrostatic potential energy into kinetic energy, resulting in directed motion (often in the form of a restoring force). At high pressures, this kinetic energy is dissipated due to collisions with the neutral gas particles, and the dust grains oscillate around their equilibrium positions with small amplitudes. At low pressures, neutral gas damping is insufficient to dissipate the energy, resulting in increased kinetic energy of the dust particles. However, due to electrostatic confinement forces and Yukawa interactions among the dust particles, the monolayer does not become dilute or gas-like. Instead, the observed state has features similar to those of 2D turbulence observed in active systems.

The above discussion implies that the most probable direction of energy flow in the simulated dusty plasma monolayers is from smaller to larger scales, i.e., an inverse energy cascade, commonly observed in active matter.^43,44 Bourgoin et al.⁴⁴ showed that in the presence of turbulence in active matter, the second-order Eulerian structure function $SE2$ exhibits a Kolmogorovian scaling (⁠$SE2∼r2/3$⁠) over slightly less than one decade of scales $r$⁠. To determine if the instability in the $0.1 Pa$ case is turbulent, we computed the second-order Eulerian structure function given by

where $r$ is a bin of spatial scales, and the average $·$ is taken over all pairs $i,j$ of particles with separation $rij$ within that bin. The velocity difference $Vijt=Vit−Vjt$ between each particle pair $i,j$ was computed using the horizontal and vertical components of the velocity vectors. As can be seen from Eq. (5), the structure function calculates the dissimilarity among particles velocities as a function of their spatial separation. At the onset of the instability, large velocity changes occur at small scales as the particles approach each other at distances smaller than the Debye length, which results in electrostatic scattering events. Prior to the onset of the instability, the structure function is a flat line because the particles do not approach each other at a distance much smaller than the average interparticle separation. This can be seen in Figs. 5(a) and 5(b), which shows the structure function obtained at the tenth second (prior to the onset of the instability) and at the 15th second (after the onset of the instability).

Figure 5(c) shows a log –log plot of the structure function $SE2r$ obtained from position and velocity data for the $0.1 Pa$ case. The data were partitioned in bins of size $≈50 μm$⁠, which is a few times larger than the mean dust diameter (⁠$10 μm$⁠), but much smaller than the Debye length (⁠$1 mm$⁠). The structure function was calculated using all 10 000 particles and averaging over 10 frames, which yields the error bars on the plot. The resulting structure function [Fig. 5(c)] exhibits a power law dependence $SE2=arb$ in the interval $270 μm<r<870 mm$ (shaded in the plot) suggesting direct energy cascade at these small scales. A power law fit of the form $arb$ was performed using the points within the region $270−870 μm$⁠, which yields an exponent $b=0.6628$ and a constant $a=2.937×106$ with $R2=0.9$⁠. Since $2/3 ≈ 0.6667$⁠, the result is reasonably close to the expected $r2/3$ dependence.

At larger scales, the structure function tends to a constant asymptotic value $SE2≈2.4×108 μm2/s2$⁠, which is expected for uncorrelated particles. This shape of the structure function is very similar to the one reported in Ref. 44 [Fig. 4(a) in that paper] for self-propelled interfacial particles exhibiting active turbulence. The extent of the region which exhibits $S2 ∼ r2/3$ dependence suggests that the energy transfer occurs when a pair of particles is separated by less than half the interparticle separation (typically $2 mm)$⁠, which also coincide with the choice of Debye length $λD=1 mm$ for this simulation.

From the structure function calculation, we also established that no two particles approach each other at a distance smaller than $rij≈300 μm$⁠. This is different from what is expected for classical self-propelled particles that can collide with each other in the absence of charge. In a typical dusty plasma experiment, each negatively charged dust grain is surrounded by a cloud of positively charged ions, which provides effective screening of the dust particle charge, with characteristic scale of this screening given by the Debye length $λD$⁠. It is expected that for spatial scales $≫λD$⁠, dust grains behave as uncorrelated particles, while for scales $≲λD$⁠, correlations become significant. Therefore, it is also reasonable to observe that the energy cascade occurs at scales much bigger than the dust diameter but comparable or slightly smaller than the Debye length.

As discussed in Sec. III B, large-scale correlations determined from the Debye length in a dusty plasma can be related to the observed energy cascade driving the system's dynamics. One manifestation of correlation effects is the observation of anomalous particle diffusion, i.e., deviations from Brownian (uncorrelated) motion. A common method for assessing particle diffusion is computing the mean squared displacement (MSD) as a function of time delay $τ$⁠. The MSD for an individual particle is given by

Equation (6) can be computed and averaged over all particles in the ensemble and plotted as a function of $τ$⁠. In the normal diffusion regime, characteristic of uncorrelated particles, the mean square displacement (MSD) of the particle ensemble increases linearly with time, i.e., $x2∼tα$ where $α=1$⁠. Deviations from normal diffusion correspond to a nontrivial topology of the corresponding phase space and spatiotemporal correlations.⁷⁵ As a result, exponents $α≠1$ are also possible, yielding two distinct examples of anomalous transport: subdiffusion when $α<1$ and superdiffusion when $α>1$⁠. Anomalous diffusion has been experimentally observed in various strongly coupled fluids, including ultracold neutral plasma,⁷⁶ 2D and quasi-2D Yukawa liquids,^77–80 and dusty plasmas.^81–83 It has been predicted that the character of the observed diffusion is sensitive to the screening length, coupling strength, energy dissipation, and choice of measurement of time scales for strongly coupled charged liquid.⁸⁴ 

Another approach to identifying the presence of anomalous diffusion is to examine the probability distribution function (PDF) of the particle displacements or velocities. Normal diffusion of uncorrelated particles is characterized by a normally distributed PDF (e.g., Gaussian or Maxwellian), while anomalous diffusion is associated with non-Gaussian statistics. Using numerical techniques from Tarantino et al.,⁶⁶ we computed the MSD plots and velocity histograms for each pressure case obtained from position and velocity data for particles within the central region of interest, as shown in Fig. 6. The number of particle trajectories detected in the central region is $≈2000$ for the $5$ and $1 Pa$ cases and $≈4000$ for the $0.1 Pa$ case. The increased number at the lowest pressure is attributed to increased dust mobility, yielding particles leaving and entering the central region throughout the selected time period of $20 s$⁠, or $2000$ frames. Thus, $∼106$ data points are considered for the statistics in each case. As shown in Figs. 6(a)–6(c), the MSD in all three cases deviates from linear growth and is better approximated by a power law fit. Table II lists the extracted exponents from the best fits of the form $MSD∼τα$⁠.

At the highest pressure of $5 Pa$⁠, the MSD fit grows slightly faster than linearly, yielding an exponent $α≈1.19$⁠, which suggests enhanced probability for superdiffusion. However, note that the vertical axis scale is in units of $∼10−7 m2$ and the maximum difference between the linear fit and the power law fit (at $t=20 s)$ is on the order of $∼10−8 m2$⁠, which is two orders of magnitude smaller than the area of a typical elementary cell $∼10−6 m2$ (determined by the distance to the nearest neighbors, which is $a∼10−3 m$⁠). Thus, although $α>1$⁠, superdiffusion most likely will not be observable in a laboratory experiment. This conclusion is further supported by the velocity histograms for this case [Fig. 6(d)], which show a minimal spread of particle velocities around the zero value.

In the $1 Pa$ case, Fig. 6(b) shows that the MSD grows slower than linearly, and a power law fit to this plot yields exponent $α≈0.83$⁠, suggesting subdiffusive behavior. The velocity histograms [Fig. 6(e)] for this case show a small spread around the zero value, and the maximum difference from linear dependence in the MSD plot is $∼10−7 m2$ [at $t=20 s$⁠, Fig. 6(b)]. As this difference is again smaller than the area of a typical elementary cell within the structure, we expect that, although at $1 Pa$ the particles are overall more mobile, they do not exhibit vastly different dynamics than the $5 Pa$ case.

The diffusive dynamics obtained using the $0.1 Pa$ data are markedly different. The MSD growth [Fig. 6(c)] is slower than linear with a power law fit coefficient $α≈0.89$⁠, indicating the possibility of subdiffusion. The difference from linear growth at long times in this case is $∼10−4 m2$⁠, which exceeds the typical area of an elementary cell size by two orders of magnitude. This suggests the possibility for large displacements across the structure, which is supported by the observed tail- broadening of the velocity histograms [Fig. 6(f)].

The conclusions drawn from the MSD plots and velocity histograms are confirmed when plotting the trajectories of all particles over the simulation time period (⁠$30 s$ or $3000$ frames). Figures 7(a) and 7(b) show the trajectories for all particles in the central region of interest over the period of $30 s$ for pressures $5$ and $1 Pa$⁠. As can be seen, in the $5 Pa$ case, dust particles are arranged in a triangular lattice and the particle displacements barely deviate from their equilibrium positions. At $1 Pa$⁠, the dynamics are very similar to the $5 Pa$ case, but small excitations at the individual cell level can be observed. In the $0.1 Pa$ case, the dust particles can have large displacements and the particle trajectories completely fill in the region of interest. For clarity, Fig. 7(c) only shows the trajectories of 100 particles within the region of interest, which is sufficient to demonstrate the different scale of particle motion in this case. Figures 7(d)–7(f) show the “average” particle trajectory in each case, which is computed by averaging the displacements of all particles at each time step. The average trajectories for $5$ and $1 Pa$ seem typical for a Brownian-like motion, while the trajectory computed from the 0.1 $Pa$ data resembles features of Lévy flights.

The spectral approach is a mathematical technique that determines the existence of a continuous component in the energy spectrum of an Anderson-type Hamiltonian of the form defined in Ref. 51. The authors have previously applied the spectral approach to the study of transport in disordered 1D, 2D, and 3D lattices,^{48,59–62,88} where the amount of stochastic disorder $c$ was varied in the potential energy term of the examined Hamiltonian. The recently developed Fractional Laplacian Spectral (FLS) method^50,85 applies the spectral approach to a Hamiltonian with random disorder in the potential energy term and a kinetic energy term represented by the discrete fractional Laplacian, which allows for modeling disordered media with nonlocal interactions.

In this paper, we are interested in energy transport through a dusty plasma monolayer, which exhibits both disorder [discussed in Sec. III A, Eq. (4)] and anomalous diffusion due to nonlocal interactions (discussed in Sec. III C). To model energy transport in such system, we consider the random fractional discrete Schrödinger operator

where $−Δs, s∈0,2$ is the discrete fractional Laplacian, $δi$ is the $i$ th standard basis vector of the 1D integer space $Z$⁠, $·,·$ is the $ℓ2Z$ inner product, and $ϵi$ are independent variables, identically distributed according to a uniform (flat) distribution on the interval $−c/2, c/2$⁠, with $c>0$⁠. It has been shown^86,87 that the nonlocal characteristics of anomalous diffusion can be modeled using a fractional Laplacian operator $−Δs$⁠, where $s∈0, 2$⁠. It is expected that in the subdiffusive regime (⁠$s>1$⁠) transport is suppressed due to negative correlations, while for $s<1$⁠, transport is enhanced due to positive correlations. Padgett et al.⁵⁰ obtained the following 1D series representation of $−Δs$ for the interval $s∈0,2$⁠:

where

Here, $u$ is a discrete function on $Z$ and $Γ$ is the standard Gamma function. The series representation of the fractional Laplacian has been recently extended to an arbitrary order in Ref. 88. Numerical simulations by Padgett et al.⁵⁰ provided confirmation that the representation in (8) and (9) yields enhanced transport (superdiffusion) for $s∈0,1$ and enhanced localization (subdiffusion) for $s∈1,2$⁠, when compared to the classical case $s=1$⁠. The present work considers the Hamiltonian in (7) with a fractional Laplacian given by (8) and (9). The main inputs in these equations are the fractional power of the Laplacian $s$ and the disorder concentration $c$ (that determines the variables $ϵi$⁠). In Sec. III A, we obtained the value of the dimensionless disorder $c$ using Eq. (4). The value $s$ is known to asymptotically approach $∼1/α$⁠, where $α$ is the exponent extracted from the MSD plots, discussed in Sec. III C and listed in Table II.

Once the appropriate Hamiltonian $Hs,ϵ$ is defined, the time evolution of the initial state of the system is generated under the iterative application of $Hs,ϵ$⁠. Since this paper considers dusty plasma monolayers where energy is imported to the system at the individual particle level, we define the diameter of a single dust particle (⁠$≈10 μm$⁠) as the smallest relevant spatial scale. Larger scales are represented as multiples of this elementary scale. Physically, this means that the dust particles are viewed as rigid spheres that cannot interpenetrate and that need to move a distance of at least one dust diameter to be considered in a new spatial location. Smaller displacements are viewed as disorder, or slight deviations from the initial position. Mathematically, this amounts to viewing the standard basis vectors $δi$ of the 1D integer space $Z$ as discrete scales proportional to one dust diameter. We let the initial state of the system be given by $δ0$⁠, located at the origin of the 1D integer space $Z$⁠, and assign to this state a normalized initial energy equal to $1$⁠. The question of interest is how does this energy spreads across scales under the iterative application of the Hamiltonian?

The time evolution of the initial state $δ0$ under the action of the Hamiltonian $Hs,ϵ$ is given by the sequence $δ0,Hs,ϵδ0,Hs,ϵ2δ0,…, Hs,ϵNδ0$⁠, where $N∈N$ is the number of timesteps (equivalently, $N$ is the number of iterations of $Hs,ϵ$⁠). Let $φ0′, φ1′,φ2′,…,φN′$ be the sequence of $ℓ2Z$ orthogonal vector states obtained from Gram–Schmidt orthogonalization of the original sequence $δ0,Hs,ϵδ0,Hs,ϵ2δ0,…, Hs,ϵNδ0$⁠. This step allows for a proper definition of a mathematical distance between subspaces in the Hilbert space. Then, for any nontrivial vector $ν≠δ0$ (representing any spatial scale of interest), we define the distance parameter (mathematical distance) as

where $φk′$ is the $k$ th term of the sequence $φ0′, φ1′,φ2′,…,φN′$⁠. Here, $·,·$ is the $ℓ2Z$ inner product and $·2=·,·$⁠. Equation (10) was originally derived in Liaw,⁴⁸ where results from spectral theory were used to verify the following conjecture:

For an Anderson-type Hamiltonian, if one can find a nontrivial vector $ν$⁠, for which the limit of the distance parameter approaches a positive value as time approaches infinity, then the spectrum of the Hamiltonian includes an absolutely continuous (ac) part, which indicates the existence of extended energy states.

The spectrum of a Hamiltonian $H$ consists of: (i) an absolutely continuous part, corresponding to extended states and (ii) a singular part, which includes discrete eigenvalues and poorly behaved transitional states (called singular-continuous part of the spectrum). If the spectrum of $H$ coincides with the singular part (i), transport in the examined problem is localized. In the presence of non-vanishing absolutely continuous part of the spectrum, de-localization occurs in the form of extended states (by the RAGE theorem, see e.g., Sec. I B of⁸⁹) In other words, if one demonstrates with positive probability that

then the time-evolved transport behavior of the system under the action of the examined Hamiltonian includes extended energy states. In Ref. 50, it was shown that the extended state conjecture also holds for the random fractional discrete Schrödinger operator in Eq. (7).

In Sec. IV B, we use the criterion (11) to determine if energy transport is likely to occur at a given spatial scale. An important distinction has to be made between energy transport and energy dissipation across scales. Here, energy transport represents the existence of extended (scattering) energy states, which is characteristic of ordered systems, such as crystalline lattices. In contrast, energy dissipation across scales is a process, which decreases the probability of energy transport at any particular scale. Thus, in the following discussion, we interpret $limN→∞Ds,ϵN>0$ as high probability for energy transport and the existence of ordered structures at the examined spatial scale, while $limN→∞Ds,ϵN=0$ as high probability for energy dissipation and the existence of disordered (dissipative) structure at the examined scale.

As discussed in Sec. IV A, the two important inputs in the FLS method are the disorder concentration $c$ and the fractional power of the Laplacian $s$⁠, which we determined from structural and diffusion analysis of the dusty plasma monolayers at the three simulated pressure cases. The values listed in Table II were used to define a Hamiltonian appropriate for each case. Next, we calculated the distance parameter from Eq. (10) for increasing values of the spatial scale, represented by the vector $v$ in Eq. (10). The final input needed for this calculation is the physically relevant $range$ of nonlocal interactions in the dusty plasma monolayer, which allows us to determine how many basis vectors $δi$ (neighboring states) are included in the calculation at each time step.

In a typical laboratory experiment, the negatively charged dust particles cause the formation of ion clouds or locally correlated Debye spheres. It has been argued in Ref. 74 that the Debye length assigns a type of large-scale uncertainty in position $Δr2∼λD$ and the positional variance can be expressed as

where $d$ is the dimension of the space. Thus, the $range$ of nonlocal interactions can be viewed as the number of standard deviations needed to quantify the expected positional variation. For a 1D case (⁠$d=1$⁠), this can be expressed as $range=nσr≅nλD$⁠, where $n$ is the number of standard deviations. We further normalize this value by the size a single dust diameter $dD≈10 μm$ (the smallest relevant spatial scale), which yields the dimensionless value $range*≅nλD/dD$⁠. For the $5 Pa$ simulation, the assumed Debye length is $λD=0.6 mm$⁠, yielding $range*≈n60$⁠. In the $1 Pa$ and $0.1 Pa$ cases, $λD=1 mm$⁠, resulting in $range*≈n100$⁠.

The value $range*$ is a numerical truncation, which determines how many neighboring states are considered when computing the action of the fractional Laplacian at each time step. We note that the application of a fractional Laplacian $−Δs$ results in interactions that, in principle, extend to infinity. Therefore, the application of the series representation in (8) and (9) is exact if an infinite number of neighbors are accounted for at each time step. However, in a dusty plasma monolayer, the effective range of nonlocal interactions has a characteristic length scale given by $λD$⁠. Thus, the numerical cutoff is not unphysical. The remaining question is whether (8) and (9) provide an accurate representation of a fractional Laplacian after the cutoff. For the three cases considered in this work, the smallest values $range*=60, 100$ are obtained if only one standard deviation $σ$ is considered. The corresponding remainder removed by the truncation is $R∼1/range2 ∼10−4$⁠. As this remainder is small, we expect it to yield a negligible contribution to the energy transport calculation. Based on the diffusion analysis in Sec. III C (see Fig. 6), although a single $σr$ may be sufficient to describe the positional spread observed for the $5 Pa$ case, the development of tails in the velocity distributions at lower pressures suggests that a choice of $2σ$ or $3σ$ is more appropriate.

Figure 8 shows the time evolution of the distance parameter, calculated for $N=3000$ timesteps, which matches the simulation time in the many-body simulations. For each set of conditions, the resulting distance plot is an average of $10$ realizations of the same numerical experiment, which minimizes fluctuations due to the random distribution of disorder values $ϵi∈−c/2, c/2$ in (7). For each case, the input values $c$⁠, $s$⁠, and $range*$ are listed in Table II. The reference vector in Eq. (10) is a linear combination of $L$ number of basis vectors in the Hilbert space, with equal weights, and has the form $v=1/L∑j=1Lδj$⁠. In this representation, a single basis vector corresponds to the minimum relevant scale (approximately equal to one dust diameter, or $10 μm$⁠). Thus, increasing the number $L$ amounts to considering larger spatial scales. Here, the examined dimensionless scales are $L=10, 1000$⁠, which corresponds to spatial scales $100 μm−1 cm$⁠.

In Fig. 8(a), corresponding to conditions obtained from the $5 Pa$ simulation, the distance values remain close to $1$ for all examined spatial scales. This corresponds to high probability for energy transport across all scales, which is typical for highly ordered crystalline lattices. At $1 Pa$⁠, the probability for transport is slightly decreased at scales $L∈200, 450$ or $2 mm−4.5 mm$⁠, which are proportional to the range of nonlocal interaction $range*=200$ or $2 mm$ for this case. Nevertheless, at the final time step considered $N=3000$⁠, the distance values computed for all scales in the $1 Pa$ case never drop to zero. Instead, the minimum value, which occurs at $L=200, 250$⁠, is $Ds,ϵ3000≈0.6$⁠. Thus, in the $1 Pa$ case, there is strong evidence that the probability for transport is still nontrivial at all scales, which suggests a well-defined crystalline structure.

From Fig. 8(c), we see that for $0.1 Pa$⁠, the distance values visibly drop from $1$ for scales $L∈100, 500$ or $1−5 mm$ and approach $0$ for scales $L∈200, 400$ or $2−4 mm$ [dark blue color in Fig. 8(c)]. For vanishing values of the distance, the spectral method cannot conclude the existence of transport at the considered vector scale $v$⁠. Thus, our interpretation of these calculations is that there is low probability for transport at the corresponding spatial scales. Instead, closer examination of the distance plots corresponding to $L∈10, 450$ [Fig. 8(d)] suggests that in this interval, the probability for transport rapidly decreases as scales increase, which can be interpreted as energy dissipation across scales. This suggests the presence of an instability in the corresponding dusty plasma monolayer, which is consistent with the observations from the many-body simulations.

Since in Fig. 8 the examined dimensionless reference scales $L∈10, 1000$ were normalized by the dust particle diameter $dD≈10 μm$⁠, the corresponding spatial scales vary in the range $100 μm$–$1 cm$⁠. Comparison to the cell orientation maps of Fig. 3(a) shows that for the $5 Pa$ case, scales in the range $100 μm−1 cm$ are smaller than the typical size of the observed crystalline domains (with characteristic 1D scales $≳2 cm$⁠). This explains the observed high probability for transport at all scales. Similarly, examination of Fig. 3(b) suggests that at $1 Pa$⁠, small sub-structures are formed within the larger crystalline domains. Those consist of individual cells with different spatial orientation or small clusters of “defect” cells that have five or seven nearest neighbors. The 1D scale of an individual cell is approximately equal to twice the mean interparticle separation, or $≈2×2 mm$ for the central regions of all the examined monolayers. In Fig. 8(b), the FLS method predicts decreased probability for transport at dimensionless scales $L∈200, 450$⁠, corresponding to the spatial range of $2−4.5 mm$⁠, which is in agreement with the observation of substructure formation within the larger crystalline domains in the $1 Pa$ simulations.

For the $0.1 Pa$ case, Fig. 3(c) clearly shows the formation of sub-cell structure as indicated by the variation of colors across portions of individual cells. Additionally, we note that this cell orientation map could only be generated in the initial timesteps of the simulation before the onset of the instability. At later timesteps, the dust particle monolayer loses crystallinity, and dust positions are observed to substantially vary from the average interparticle separation. This is evident both from Fig. 7 and from the high value of the corresponding dimensionless disorder $c=1.3×10−3$ listed in Table II. Figures 8(c) and 8(d) suggest that in the interval $L∈200, 400$⁠, corresponding to spatial scales of $2−4 mm$⁠, transport is not likely to occur. As these spatial scales coincide with $1−2$ interparticle distances, this result confirms the breaking of dust ordering into well-defined cells, observed in the simulation. The higher probability for energy transport at larger scales $L∈500, 1000$⁠, or $5−10 mm$⁠, can be explained by the presence of global liquid-like behavior of the monolayer. In other words, although dust ordering in individual cells is lost, the structure is not chaotic and still preserves the characteristics of a strongly coupled monolayer.

In this article, we presented a theoretical study of energy transport in dusty plasma monolayers. Specifically, we considered three sets of experimentally relevant conditions, which are expected to result in a highly ordered crystalline structure at $5 Pa$⁠, a disordered crystalline structure at $1 Pa$⁠, and a liquid-like structure at $0.1 Pa$⁠. The latter case is observed to develop a self-excited instability, having the characteristics of active turbulence, where energy cascades from smaller to larger scales. We modeled these three structures using a many-body simulation of 10 000 spherical dust particles, whose dust diameter was allowed to vary by $±0.9%$⁠. The variation in the particle size introduces variation of dust particle mass and charge, which affects the spatial distribution of particles within the monolayer. To mimic the confinement forces expected in large-electrode dusty plasma cells, we assumed a tenth-order radial confinement potential, which has a flat profile in the center of the simulation box and steep slope close to the exterior. The resulting dusty plasma structures form central regions of particles, where radial confinement is provided indirectly through an external later of particles. In this manner, focusing the analysis on these central regions of dust particles allows to minimize the effect of boundary conditions.

For each examined pressure case, we performed structural and diffusion analysis of the simulated dusty plasmas, which confirmed that for $5$ and 1 $Pa$⁠, the obtained monolayers are crystalline, while the $0.1 Pa$ condition yields a liquid-like structure that develops a self-excited turbulent instability. For each case, we computed a dimensionless disorder $c$ from the standard deviation from the mean interparticle separation, which was found to increase inversely with pressure. The dust particle diffusion regime was identified from analysis of the mean squared displacement (MSD) as a function of time and from velocity histograms. While the particles in all three cases are found to move according to anomalous diffusion, only at $0.1 Pa$ does the observed deviation from classical diffusion result in a significant change of the individual particle trajectories. This is also confirmed by the velocity histograms, which exhibit tail broadening at this pressure. For each case, a power law fit to the plot of $MSDτ$ was used to extract an exponent $α$⁠, which quantifies the deviation from classical diffusive behavior. The observation of anomalous diffusion, or $α≠1$⁠, suggests the presence of nonlocal interactions, which is expected for dusty particles interacting through a shielded Coulomb (Yukawa) potential. The characteristic range of such nonlocal interactions was determined from the Debye length $λD$ selected for each simulation.

To study energy transport across scales in each monolayer, we used a Fractional Laplacian Spectral (FLS) method, where the time evolution of the initial energy state of the system is obtained from the iterative application of the corresponding Hamiltonian. To account for both random disorder and nonlocal interactions, we employed a random fractional discrete Schrödinger operator as the Hamiltonian. The fractional power $s$ of the Laplacian is known to asymptotically approach $1/α$⁠, where $α$ is the exponent characterizing time evolution of the MSD plots. This fraction and the range of nonlocal interactions determined from the many-body simulations were used to obtain the fractional Laplacian for each case. The dimensionless disorder calculated for each monolayer was used to define the potential term in the corresponding Hamiltonian. Using the FLS method, we computed the distance parameter in Eq. (10), which quantifies the probability for transport (in the form of extended energy states) as a function of spatial scale. The results showed that at $5 Pa$⁠, energy transport is likely to occur at all scales, in agreement with the fact that the corresponding structure is highly ordered with large crystalline domains. At $1 Pa$⁠, the probability for energy transport is slightly decreased at energy scales proportional to the size of individual elementary cells, which coincides with the observed formation of substructures within the large crystalline domains in that case. Finally, at $0.1 Pa$⁠, the FLS calculations demonstrated that energy transport is suppressed for spatial scales proportional to the size of individual cells but preserved at larger spatial scales. These results were found to coincide with the observed breaking of elementary cells as the monolayer transitioned to unstable liquid-like structure.

The possible relationship between classical inertial turbulence and 2D turbulence has been studied in monolayers of active particles,⁴⁴ but still presents a relatively unexplored question in plasma physics. Here we showed that, for a dusty plasma monolayer at $0.1 Pa$⁠, the time evolution of the MSD plots exhibits a transition from uncorrelated to correlated regime at intermediate time scales. It was further demonstrated that the velocity histograms deviate from Maxwellian distributions by developing “fat” tails, which suggests the anomalous diffusion of the dust particles. While both these features are expected for 2D turbulence, the presence of $r2/3$ dependence in the corresponding structure function plot is a feature typical of 3D inertial turbulence. The coexistence of these features at small scales is possible due to the quasi-2D nature of dusty plasma monolayers, where random fluctuations in the $z$-displacements of the particles can result in the onset of a global instability in the monolayer plane. Additionally, electrostatic scattering among charged particles in a dusty plasma introduce a mechanism for energy transfer at scales on the order of the Debye length (much bigger than scales proportional to the particle size, where energy is transferred through collisions). Thus, dusty plasma monolayers present an appropriate platform for studying the fundamental physical mechanisms driving both direct and inverse energy cascades in complex systems, where both random disorders and nonlocal interactions may be present.

This work was supported by the NSF Grant Nos. 1903450 (EGK, JLP, CDL, and LSM), 1707215 (LSM and TWH), and 1740203 (TWH and LSM), NSF-DMS Grant No. 1802682 (CDL), NASA Grant No. 1571701 (TWH and LSM), and DOE Grant No. SC0021284 (EGK).

Since August 2020, C. D. L. has been serving as a Program Director in the Division of Mathematical Sciences at the National Science Foundation (NSF), USA, and as a component of this position, she received support from NSF for research, which included work on this paper. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

In Sec. II B, we briefly outlined some features of the many-body simulations, which were most relevant to the discussion. Here, we elaborate on each term in the force equation (repeated for convenience), which was used to obtain the dust particle dynamics

Below we discuss the dust-dust interaction force $Fdd$⁠, the interplay between gravity $mdg$ and confinement forces $QdE$⁠, the drag force $βẋ$⁠, and the thermal bath $ζrt$ used in the many-body simulations of the present study.

In typical complex plasma experiments, the primary dust-charging mechanism is the collection of electrons and ions from the environment. Due to the electrons' higher thermal speed, dust grains generally become negatively charged and are surrounded by a region of positive space charge due to ion shielding. The resulting local interparticle interaction is commonly assumed to be of the Yukawa (screened Coulomb) form

where $r$ is the distance from a particle with charge $Qd$ and $λD$ is the Debye shielding length (the scale length over which a charged grain is shielded by the plasma). It has been reported⁹⁰ that in strongly coupled liquids with Yukawa interactions, the laminar-to-turbulent transition at a fixed Reynolds number $Re$ can be induced by increasing the strength of the interparticle potential alone. This behavior has been verified across different ranges of Reynolds numbers (low $Re<0.1$⁠, intermediate $2≤Re≤35$⁠, and high $Re>50$⁠). Thus, in the present study of active turbulence, it was considered appropriate to derive the dust-dust interaction force $Fdd$ from a Yukawa potential. In a future study, it will be interesting to consider how modifications on the Debye spheres and the corresponding dust–dust interaction influence the observed dust diffusion and the onset of turbulent instabilities.

The confinement forces in Eqs. (1)/(A1) are electrostatic in nature and have the form

where the radial confinement is derived from the tenth order polynomial potential

which was discussed in Sec. II B. The vertical confinement in a dusty plasma experiment is typically provided by the electrodes in the cell. Here, we model such confinement by a fifth order polynomial electrostatic field, which was obtained from a fluid model of the plasma in a GEC cell^91,92

where $z$ is the height above the lower electrode, as shown in Fig. 9.

Note that $ẑ$ is positive in the direction from the lower to the upper electrode, which is opposite to the direction of gravity. In the sheath region of the plasma, the direction of the electric field coincides with the direction of gravity, i.e., $−ẑ$⁠. Since the dust particles are negatively charged, the resulting electric field force is in the $+ẑ$ direction, balancing the force due to gravity, which allows dust levitation at the exact location where forces balance. The fifth order polynomial was obtained from fits to a fluid model code and experimental data, discussed in Ref. 53. This polynomial fit was obtained for particular plasma conditions. In general, as the power delivered to the plasma increases, the confinement force becomes stronger with a steeper $Ez$⁠. In the region where the dust is levitating in the sheath, the electric field is often approximated as linear. The bottom electrode is assumed to be located at $z=0 m$⁠, and the top at $z=0.0254 m$⁠, after which the electric field takes a constant value of $8083 V/m$⁠.

As a result, a dust grain of mass $m$ and negative charge of magnitude $Qd$ acquires acceleration

This electric force is exerted against gravity, so that the electro-gravitational (eg) acceleration in the $z$ direction is given in the simulation by

The gas drag force in the many-body simulation is given by $Fdrag=−βẋ$⁠, where the damping factor $β$ depends on the neutral gas pressure and temperature, with

Here, $δ=1.44$ (measured for melamine formaldehyde dust in argon), $rd$ is the dust radius, $n$ is the gas number density, $mg$ is the molecular mass of the gas, $Tg$ is the gas temperature (assumed to be $300 K$⁠), and $md$ is the mass of the dust. The Epstein damping times for the simulations presented here are 0.1497 s (at 5 Pa), 0.7484 s (at 1 Pa), and 7.4838 s (at 0.1 Pa). Note that even for the longest damping time, $≈7.5 s$⁠, the 30 s simulation time covers four damping periods.

A thermal bath is realized by subjecting the particles to random force “kicks,” $Fth=ζrt$ with the maximum acceleration imparted by the amplitude

where $Δtd$ is the dust time step. Notice that each kick represents a cumulative effect from neutral collisions with the dust particle, which yield a measurable displacement at the characteristic timescale of the dust. The random number $rt$ is selected from a Gaussian distribution with a unit normal distribution (zero mean and unit variance). Gaussian distribution of random variables should yield a classical diffusion of the dust grains (also called, Brownian motion or a random walk). The dust particles are expected to exhibit classical diffusion in cases where all interactions are local, which is true for collisions with neutral gas particles from the environment. A Gaussian distribution with a narrow width leads to small random numbers $rt$⁠, which is appropriate for modeling strongly coupled dust crystals at high gas pressure. As the variance of the Gaussian is increased, bigger random numbers $rt$ are more likely, which leads to less stable dust grains. Thus, in the present study, the thermal bath is responsible for mobilization of dust grains at lower pressures. However, the observed anomalous dust diffusion results from nonlocal interactions due to charge effects and not from the thermal bath.

Extracting temperature from Gaussian fits is not straightforward for interacting particles in general and is particularly difficult for dusty plasma monolayers, which are known to exhibit non-Maxwellian velocity PDFs. Additionally, the kinetic temperature measured from fits to the velocity histograms is known to be highly anisotropic and larger than room temperature at low pressures (as discussed in Refs. 93 and 94). In Fig. 10, we show different fits to the velocity PDFs for each pressure case. For each pressure, the $vz$ components exhibit strongly non-Maxwellian distributions and are best described by a combination of Maxwellian distribution and Cauchy distribution tails. The $vx$ and $vy$ components in each case can be described by a Maxwellian fits, which indicate that the temperature increases substantially as pressure decreases: $≈130–135 K$ (at 5 Pa), $≈650–700 K$ (at 1 Pa), and $≈4300 K$ (at 0.1 Pa). Although the kinetic temperature extracted for the 0.1 Pa case seems very high, temperatures on the order of eV are consistent with temperatures reported by Williams and Thomas in Refs. 93 and 94.

Alternatively, one can fit a Maxwellian distribution with Kappa tails to the $vx$ and $vy$ PDFs in the 0.1 Pa case, which yields a much lower temperature $≈400 K$⁠. This is consistent with the temperature obtained from the fit to the $vz$ PDF. The parameter $κ$ characterizes how far the system is from thermal equilibrium. When $κ→∞$⁠, the Kappa distribution approaches a Maxwellian, but when $κ$ is finite, the distribution function has high-energy tails, with more high velocity particles than for a Maxwellian. For $κ<10$⁠, the Kappa distribution has a power law tail. Here, $κ≈11$ for both $vx$ and $vy$ PDFs, suggesting that the 0.1 Pa case is away from thermal equilibrium but does not exhibit power-law tails. The modified Maxwellian distribution with Kappa tails has been used to describe the velocity PDFs of a weakly coupled 3D dusty plasma cloud, suspended in microgravity Neon plasma, powered by a RF coil.⁹⁵ Those studies found $κ=1.72$ and $T=161 K$⁠. Note that these experiments did not include external heating or perturbation.

The data that support the findings of this study are available from the corresponding author upon reasonable request.