Fractional Laplacian Spectral Approach to Turbulence in a Dusty Plasma Monolayer
Evdokiya G Kostadinova, Rahul Banka, Joshua L Padgett, Constanze D Liaw, Lorin S Matthews, Truell W Hyde
11 Active turbulence in a dusty plasma monolayer
E G Kostadinova ¹, R Banka ¹, J L Padgett ¹ ٫ ², C D Liaw ¹ ٫ ³, L S Matthews ¹, & T W Hyde ¹ ¹ CASPER and Department of Physics, Baylor University, Waco, TX 76706, USA ² Department of Mathematical Sciences, University of Arkansas, Fayetteville, Arkansas 72701, USA ³ Department of Mathematical Sciences, University of Delaware, Newark, DE 19716, USA E-mail: [email protected], [email protected], [email protected], [email protected], [email protected], [email protected] Abstract.
This work presents a theoretical investigation of active turbulence in a dusty plasma monolayer, where energy is injected at the individual particle level and transported to larger scales, leading to the spontaneous formation of spatially disordered flow patterns. Many-body simulations of 10,000-particle dusty plasma monolayers are used to demonstrate how the global dynamics depend on the statistical properties of the dust assembly for realistic laboratory conditions. We find that spatial defects due to variations in the dust size distribution and charge-driven nonlocal interactions resulting in anomalous dust diffusion, are key factors for the onset of instabilities. The resulting dynamics exhibits features of inertial turbulence over slightly less than one decade of scales, down to scales on the order of the particle diameter. These processes are examined analytically using a recently developed Fractional Laplacian Spectral (FLS) technique, which identifies the active energy channels as a function of scale, disorder concentration, and nonlocality. The predictions from the theoretical (spectral) analysis demonstrate agreement with the results from the many-body (kinetic) simulations, thus providing a powerful tool for the study of active turbulence.
Keywords: complex (dusty) plasma, active turbulence, fractional Laplacian, spectral approach
I. INTRODUCTION 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 e xample 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 [1]. 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 [10], 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 (
𝑅𝑒 ≲ 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], stochastic self-propelled devices [20], [21], airfoil devices [22]–[24], and insect flights [25], [26]. 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. However, such fluids are often characterized by small-scale, high-frequency dynamics, violating the assumptions needed for a hydrodynamic approximation. Thus, 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 can be studied at the kinetic level in the low-to-medium Reynolds number regime 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 [38], 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. This offers a unique opportunity to investigate at the kinetic level the fundamental charge-related processes responsible for the onset of turbulence. Here we study active 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 [42]. As an inverse cascade is typical for two-dimensional (2D) turbulence, it has been observed in active matter monolayer systems, such as bacterial suspensions [43] and self-propelled camphor swimmers [44]. 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 each dust particle with the gas and plasma environment. Additionally, the dust-plasma charging processes yield both non local interactions and stochastic fluctuations in the dust particle spatial distribution. Thus, the investigation of self -excited instabilities in such structures can provide insight into the connection between charge-driven processes and the inverse energy cascade in active turbulence. Here, we conduct many-body 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 [45], 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 ‘kicks’ of various strength and direction to each dust grain.
Comparisons between stable and unstable monolayers are made by considering three pressure regimes, 𝑝 = [5, 1, 0.1]𝑃𝑎 , which allow for simulation of both crystalline and liquid -like monolayers. At , a stable monolayer is formed with distinct hexagonal symmetry and large crystalline domains. At , 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 crystalline regions. As the pressure drops to , the crystalline domains are destroyed, and an instability develops in the dynamics of the monolayer. In each case, energy is only imported in the structure at the individual particle level by the thermal bath, and subsequently distributed across the monolayer in the form of excitations of various scales. This allows us to investigate 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. 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 dynamical state solely from knowledge of the underlying Hamiltonian structure. The spectral approach employed here (introduced in [48], [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 [50]. 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 [16]. In the many-body 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. The kinetic energy term is obtained from the characteristics of the dust particle diffusion observed in the simulation, including 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 simulation conditions is determined, the FLS technique is used to identify the existence of transport as a function of scale, solely from the properties of the corresponding energy spectrum. Thus, an agreement between the (kinetic) may-body simulations and the (spectral) FLS code predictions provides a powerful modeling tool for studying energy cascade in turbulent dynamics. In Sec. II, we introduce the main features of the many-body simulation and discuss scaling between simulation input parameters and relevant experimental conditions. Structural and statistical analyses of the resulting dusty plasma monolayers are presented in Sec. III, which also provides scaling between the dynamical features 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 and spectral models 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 [50].
II. SIMULATION OF DUSTY PLASMA MONOLAYERS A. Scaling of Experimental Parameters An outstanding challenge in the experimental study of global instabilities in complex plasmas 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 1a), located in Baylor University’s CASPER lab. Cell 3 is based on the GEC RF reference cell chamber [51], with several modifications to the original design required for suspending dust particles in the discharge and visualizing them using different optical systems [52]–[54]. 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 1b shows the outer region of such a dusty plasma monolayer suspended in argon plasma at gas pressure . 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 in Meyer et al. [45]. The structures discussed in Meyer et al. consist of ≈ 10,000 dust particles, forming monolayers of diameter ≈ 27𝑐𝑚 . At , these extended monolayers are observed to exhibit wave instabilities and mode coupling, which suggests that this lower pressure is appropriate for modeling fluid dynamics phenomena, such as turbulence. Thus, to investigate energy transport through extended dusty plasmas, here we focus on modeling monolayers consisting of 10,000 particles at pressures and ; 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 . 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 [30].) Table I summarizes the experimentally relevant conditions used as inputs for the numerical simulations of dusty plasma monolayers at the three considered pressure cases.
Fig. 1. a) Cell 3: a large-electrode RF cell, located at Baylor University’s CASPER lab. b) Outer region of an extended dusty plasma monolayer generated in Cell 3, argon plasma, pressure . At this pressure, the monolayer exhibits both crystalline domains and disordered regions.
Disordered region
Crystalline region a) b) Table I. Conditions used in the many-body simulations of dusty plasma monolayers Dust Parameters Gas / Plasma Environment Type MF spheres discharge RF glow, Argon Diameter [𝜇𝑚]
Pressure [𝑃𝑎] Mass [× 10 −13 𝑘𝑔] 𝜆 𝐷 [𝑚𝑚] < 𝑄 𝑑 > [ 𝑒 − ] 𝛺 ℎ [𝑠 −1 ]
2𝜋 × 0.12
Particle number
Frame rate [𝑓𝑝𝑠]
In Table I, the dust charge 𝑄 𝑑 is in units of elementary charge 𝑒 − , 𝜆 𝐷 is the Debye screening length due to the electron and ions in the plasma, and 𝛺 ℎ is the radial confinement frequency (discussed in the next section) . For the cases where 𝑝 = 1𝑃𝑎, 0.1𝑃𝑎 , the selected Debye length is 𝜆 𝐷 = 1𝑚𝑚 , which has been reported as a typical value at lower pressures in Argon RF glow discharge [45], [55]. In [45], the charge on the dust grains was estimated to 𝑄 𝑑 = 24,648𝑒 − which accounts for the reduction of charge due to interactions with ion wakefields. At 𝑝 = 5𝑃𝑎 , the higher plasma density leads to a reduced Debye length of 𝜆 𝐷 = 0.6𝑚𝑚 , as discussed in [30]. Since dust charge reduction due to ion-neutral collisions is also expected, the simulation employs a dust charge 𝑄 𝑑 =18,570𝑒 − , as suggested in [56]. In Table I, the values of 𝑄 𝑑 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 [55]. Thus, here the dust particle diameters are selected from a normal distribution with a mean and a standard deviation of ±0.9% . This also leads to variation in the corresponding mass and charge on each dust grain. B. Main Features of the Many-Body Simulation To simulate the dynamics of dusty plasma monolayers, we employ a self-consistent 𝑁 -body code optimized to include specified external potentials suitable for modeling the confinement forces and external electric fields characteristic of laboratory complex 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 [57]–[61] and in GEC RF reference cells in Earth-based experiments [62]–[65]. In the many-body simulation, the equation of motion for a single dust particle with mass 𝑚 𝑑 and charge 𝑄 𝑑 is given by 𝑚 𝑑 𝒙̈ = 𝑭 𝑑𝑑 + 𝑚 𝑑 𝒈 + 𝑄 𝑑 𝑬 − 𝛽𝒙̇ + 𝜁𝒓(𝑡), (1) where 𝑭 𝑑𝑑 is the force between pairs of dust grains, 𝑚 𝑑 𝒈 is the gravitational force, 𝑄 𝑑 𝑬 is the confining electric field force, and 𝛽𝒙̇ is the neutral drag force. The force between pairs of dust grains 𝑭 𝑑𝑑 is assumed to be a shielded Coulomb (Yukawa) interaction, where the screening is mainly provided by the ions and determined from the Debye length 𝜆 𝐷 . The last term 𝜁𝒓(𝑡) 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 r(t) 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 𝑖 𝑡ℎ particle is a tenth-order polynomial of the form suggested in [45] 𝑉 𝑖 = 0.5𝑚 𝑑 (𝛺 ℎ2 𝜌 𝑖10 𝑅 ), (2) where 𝛺 ℎ is the horizontal confinement frequency, 𝜌 𝑖 is the radial position of the 𝑖 𝑡ℎ particle in the xy-plane, and 𝑅 is the approximate horizontal radius of the monolayer (radial extent). Following Meyer et al. [45], we assume 𝛺 ℎ = 2𝜋 × 0.12𝑠 −1 and 𝑅 = 63𝑚𝑚 . Data for the particle positions, velocities, and accelerations is output every 0.001 s. This output frequency is selected to match a frame rate of , typically used in laboratory dusty plasma experiments. Thus, in the following discussion, timesteps in the many-body simulation are sometimes referred to as frames. III. STRUCTURE AND DYNAMICS OF SIMULATED MONOLAYERS Using experimentally relevant conditions from Table 1, simulations of extended dusty plasma monolayers were conducted for the three pressure cases, 𝑝 = 5𝑃𝑎, 1𝑃𝑎, and, 0.1𝑃𝑎 , which are expected to yield different structure and dynamics. Figure 2abc show the dust particles 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 and , while the case yields a less ordered, liquid-like structure. In each simulation, the properties of the radial confinement in (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𝑚𝑚 (the value of 𝑅 in (2)), the confinement force grows rapidly. As a result, particles with radial positions ≲ 63𝑚𝑚 are mostly confined by an external layer of particles located at radial positions ≳ 63𝑚𝑚 . As the internal structure in each case is more ordered and less dependent on the strong boundary conditions, in the following sections we focus the analysis only on the particles located within a region with radius 𝑅 𝑐 = , indicated by red circles in Fig. 2abc. The stability of the system at each pressure is demonstrated by the time evolution of the radial and axial velocities, as shown in Fig. 2def. For the case, after 30s of simulation time, the horizontal velocities limit to 𝑉 𝑥 ≈ 𝑉 𝑦 ≈ 2 𝑐𝑚 𝑠⁄ , while the vertical velocities approach 𝑉 𝑧 ≈ 0.3 𝑚𝑚 𝑠⁄ . Similar vertical velocities are observed at , while the limiting horizontal velocities approach greater values of 𝑉 𝑥 ≈ 𝑉 𝑦 ≈ 11 𝑐𝑚 𝑠⁄ . The case exhibits distinct behavior as the particle velocities experience a rapid increase with large fluctuation in both the vertical and horizontal directions, suggesting the onset of an instability. As shown in Fig. 2f the vertical velocity increases steadily for the first 14s, followed by a rapid growth in the last 16s. The growth in the horizontal velocities occur ≈ 5𝑠 after the beginning of growth in 𝑉 𝑧 . 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 𝑧 -direction due to the assumed variations in dust particle size. This phenomenon has been observed experimentally by Gogia and Burton [55]. A. Active Turbulence As discussed in the previous section, an instability occurs in the dusty plasma monolayer at 𝑝 =0.1𝑃𝑎 (Fig. 2f), which suggests the possibility of turbulent dynamics. As this instability develops at later simulation times, 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 simulation period. The only source of energy input in the simulation occurs at the individual particle level, where random ‘kicks’ from the thermal bath are used to model collisions with neutral gas atoms. This implies that the only possible direction of energy flow in the structures is from smaller to larger scales, i.e., an inverse energy cascade, commonly observed in active matter [43], [44]. Bourgoin et al. [44] suggested that in the presence of active turbulence characterized by such an inverse cascade, the second-order Eulerian structure function 𝑆 𝐸2 exhibits a Kolmogorovian scaling ( 𝑆 𝐸2 ~𝑟 ) over slightly less than one decade of scales 𝑟 . To determine if the instability in the case is turbulent, we computed the second-order Eulerian structure function given by Fig. 2. Time-evolved radial positions of dust particles for a) , b) , and c) . The red circle in each plot indicates the central region ( 𝑅 𝑐 = 50𝑚𝑚 ) of particles used for the diffusion and spectral analysis. Plots of the horizontal and vertical velocities as a function of simulation time are shown for d) , e) , and f) . Shaded areas indicate regions where large velocity changes are observed. a) b) c) d) e) f) 𝑆 𝐸2 (𝑟) = 〈|(𝑽 𝒊𝒋 (𝑡) ∙ 𝒓 𝒊𝒋 ) 𝑟 𝑖𝑗 ⁄ | 〉, (3) where 𝑟 is a bin of spatial scales, and the average 〈∙〉 is taken over all pairs (𝑖, 𝑗) of particles with separation 𝑟 𝑖𝑗 within that bin. The velocity difference 𝑽 𝒊𝒋 (𝑡) = 𝑽 𝒊 (𝑡) − 𝑽 𝒋 (𝑡) between each particle pair (𝑖, 𝑗) was computed using the horizontal and vertical components of the velocity vectors. The structure function 𝑆 𝐸2 was computed for all particles in a given frame and averaged over the last 20 frames. Figure 3 shows a log-log plot of the structure function 𝑆 𝐸2 (𝑟) obtained from position and velocity data for the case, together with a curve fit of 𝑆 𝐸2 ~𝑟 . The data was partitioned in bins of size ≈ 10𝜇𝑚 , which is close to the mean dust diameter. The shaded region in the plot highlights the distances in the range [30,90] 𝜇𝑚 where the structure function exhibits a dependence 𝑆 𝐸2 ~𝑟 , suggesting turbulent dynamics at small scales. At larger scales, the structure function tends to a constant asymptotic value 𝑆 𝐸2 ≈ 2.7 × 10 𝜇𝑚 𝑠 ⁄ , which is expected for uncorrelated particles. This shape of the structure functions is similar to the one reported in [44] for self-propelled interfacial particles exhibiting active turbulence. An important difference here is the higher values of 𝑆 𝐸2 at the two smallest scales, corresponding to distances in the range ≈ (10 − 20)𝜇𝑚 , which was much smaller for the self -propelled swimmers in [44] (see Fig. 4a of that paper). Large values of 𝑆 𝐸2 indicate a large differences in velocities of particles separated by distances 𝑟 𝑖𝑗 ≈ (10 − 20)𝜇𝑚 , which is to be expected in our system, where the particles are all negatively charged and strongly repel each other at small separations. B. Crystallinity and Disorder Concentration 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 particle s and introduces lattice defects throughout the monolayer. Two-dimensional (2D) lattice defects can be defined as Fig. 3.Second-order Eulerian structure function 𝑆 𝐸2 from particle positions and velocities obtained in the simulation. The dashed line represents the fit 𝑆 𝐸 ~𝑟 , while the shaded region highlights the scales at which 𝑆 𝐸 from the simulation follows a 𝑆 𝐸 ~𝑟 trend. the fraction of particles with number of nearest neighbors 𝑁𝑁 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. [30], which uses a Delaunay triangulation algorithm to calculate the number of nearest neighbors for each dust grain and the complex bond-order parameter 𝐺 (𝑖) = 16 ∑ 𝑒𝑥𝑝(𝑖6Θ 𝑖 (𝑙)) 𝑁𝑁 𝑖 𝑙=1 , (4) where, 𝑁𝑁 𝑖 is the number of nearest neighbors of the 𝑖 th particle and Θ 𝑖 (𝑙) is the angle of the 𝑙 th nearest-neighbor bond measured with respect to the 𝑋 -axis. When the code is applied to dust fluid structures with badly defined primitive cells, the Delaunay triangulation function returns an error due to 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 well-defined crystalline structure and monolayers with liquid-like dynamics. Figure 4 shows the complex argument 𝐴𝑟𝑔(𝐺 ) , which measures the overall angular orientation of a given particle neighborhood. The three plots in Fig. 4 are obtained from dust particle positions at 𝑡 = 5𝑠 , at which time the and 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 and cases, while long-distance order is lost in the realization. In addition, the fraction of particles with 𝑁𝑁 ≠ 6 increases as pressure decreases, yielding high defect concentration at . Fig. 4. Complex phase angle of the bond- order parameter 𝐺 , calculated from dusty plasma simulations for pressure a) , b) , and c) . In each plot, black dots mark dust particles with six nearest neighbors, circles mark dust particles with seven nearest neighbors or more, and triangles mark dust particles with five nearest neighbors or less. a) b) c) Since the monolayer at does not preserve well-defined hexagonal cells throughout the simulation, it is more appropriate to define a one-dimensional (1D) disorder of the form 𝑐 ≈ 𝜎 𝑎 𝑎 , (5) where 𝜎 𝑎 is the standard deviation of the average interparticle separation 𝑎 computed for all dust particles within the region of interest (red circle of radius 𝑅 𝑐 = shown in Fig. 2). We note that typical interparticle separations in this central region are ≈ 2𝑚𝑚 in each pressure case, which is larger than the separation of ≈ 1.7𝑚𝑚 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 𝑐 for dust particles in the region of interest for each pressure case. Those values are used as input for the Fractional Laplacian Spectral (FLS) method discussed in Sec. IV. C. Dust Particle Diffusion As discussed in Sec. III A, when assemblies of charged particles are considered, the small-scale behavior of the second-order Eulerian structure parameter deviates from that which is expected for active turbulence 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. The characteristic scale of this screening is the Debye length 𝜆 𝐷 , which represents the distance from a dust particle at which the electrostatic potential is reduced by , where 𝑒 is Euler’s number. Thus, it is expected that for spatial scales ≫ 𝜆 𝐷 , dust grains behave as uncorrelated particles, while for scales ≲ 𝜆 𝐷 , correlations become significant. One manifestation of correlation effects is the observation of anomalous particle diffusion, i.e., deviations from Brownian 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 𝑀𝑆𝐷 𝑖 (𝜏) = ∑ [𝑟 𝑖 (𝑡 + 𝜏) − 𝑟 𝑖 (𝑡)] . (6) 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., 〈𝑥 〉~𝑡 𝛼 where 𝛼 = 1 . Deviations from normal diffusion correspond to a nontrivial topology of the corresponding phase space and spatiotemporal correlations [66]. 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 [67], 2D and quasi-2D Yukawa liquids [68]–[71], and dusty plasmas [72]–[74]. 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 [75]. 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 . [66], 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. 5. The number of particle trajectories detected in the central region is ≈ 2,000 for the and cases and ≈ 4,000 for the 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 , or frames. Thus, ∼ 10 data points are considered for the statistics in each case. As shown in Figs. 5abc, 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 𝑀𝑆𝐷 ∼ 𝜏 𝛼 . Table II. Parameters measured from structural and diffusion analysis of the dusty plasma monolayers. Highlighted in red are parameters used as inputs in the FLS analysis in Sec. IV.
Pressure [Pa]
5 1 0.1 𝑹 𝒄 a [ × 𝟏𝟎 −𝟑 m] 𝝈 𝒂 [ × 𝟏𝟎 −𝟕 m] 𝜶 c [ × 𝟏𝟎 −𝟒 ] 𝒔~ 𝟏 𝜶⁄ 𝒓𝒂𝒏𝒈𝒆
200 300 300
At the highest pressure of , the MSD fit grows slightly faster than linearly, yielding an exponent 𝛼 ≈ 1.19 , which suggests enhanced probability for superdiffusion. However, notice that the vertical axis scale is ~10 −7 𝑚 and the deviation itself is on the order of ~10 −8 𝑚 , which is two orders of magnitude greater than the area of a typical elementary cell ~10 −6 𝑚 (determined by the distance to the nearest neighbors, which is 𝑎~10 −3 𝑚 ). Thus, although we obtained 𝛼 > 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. 5d), which show a minimal spread of particle velocities around the zero value. In the case, Fig. 5b 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. 5e) for this case show a small spread around the zero value and the deviation from linear dependence in the MSD plot is ~10 −7 𝑚 (Fig. 5d). As this spread s again smaller than the area of a typical cell within the structure, we expect that, although at the particles are overall more mobile, they do not exhibit vastly different dynamics than the case. The diffusive dynamics obtained using the data are markedly different. The MSD growth (Fig. 5c) is slower than linear with a power law fit coefficient 𝛼 ≈ 0.89 , indicating the possibility of subdiffusion. The deviation from linear growth in this case is ~10 −4 𝑚 , which exceeds the typical area of a cell size by two orders of magnitude. This suggests the possibility for large displacements across the structure, which is supported by the o bserved tail broadening of the velocity histograms (Fig. 5f). The conclusions drawn from the MSD plots and velocity histograms are confirmed when plotting the trajectories of all particles over the simulation time period ( or frames). Figure 6ab shows the trajectories for all particles in the central region of interest over the period of for pressures and . As can be seen, in the case, dust particles are arranged in a triangular lattice and the particle displacements barely deviate from their initial positions. At , the dynamics are very similar to the case, but small excitations at the individual cell level can be observed. In the case, the dust particles can have large displacements and the particle trajectories completely fill in the region of interest. For clarity, Fig. 6c only shows the trajectories of 100 particles within the region of interest, which clearly demonstrates the different scale of particle motions in this case. Figure 6def shows the ‘average’ particle trajectory in each case, which is computed by averaging the displacements of all particles at each timestep. The av erage trajectories for and seem typical for a Brownian-like motion, while the trajectory computed from the 0.1 𝑃𝑎 data resembles a Lévy flight. Fig. 5. Mean squared displacement (MSD) as a function of time delay for a) , b) , and c)
The gray shaded region shows the range of MSD values computed for all of the particles, with the average MSD shown by the black line. In each plot, a linear fit to the data is indicated by a red line, while a power law fit to the data is represented by a blue line. The corresponding velocity histograms are shown for d) , e) , and f) . a) b) c) d) e) f) IV. ENERGY TRANSPORT A. Fractional Laplacian Spectral (FLS) Approach 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 [76].
The authors have previously applied the spectral approach to the study of transport in disordered
Fig. 6. Dust particle trajectories during the entire simulation time ( ) for a) , b) , and c) . Average trajectories generated from the average displacements of all particles in each frame for d) , e) , and f) . a) b) c) d) e) f) 𝑐 was varied in the potential energy term of the examined Hamiltonian. The recently developed Fractional Laplacian Spectral (FLS) method [50], [77] 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 B, Eq. (5)) and anomalous diffusion due to nonlocal interactions (discussed in Sec. III C). To model energy transport in such sys tem, we consider the random fractional discrete Schrödinger operator 𝐻 𝑠,𝜖 ∶= (−Δ) 𝑠 + ∑ 𝜖 𝑖 〈∙, 𝛿 𝑖 〉𝛿 𝑖𝑖∈ℤ , (7) where (−Δ) 𝑠 , 𝑠 ∈ (0,2) is the discrete fractional Laplacian, 𝛿 𝑖 is the 𝑖 th standard basis vector of the 1D integer space ℤ , 〈∙,∙〉 is the ℓ (ℤ) inner product, and 𝜖 𝑖 are independent variables, identically distributed according to a uniform (flat) distribution on the interval [−𝑐 2, 𝑐 2 ⁄⁄ ] , with 𝑐 > 0 . It has been shown [78], [79] that the nonlocal characteristics of anomalous diffusion can be modeled using a fractional Laplacian operator (−Δ) 𝑠 , where 𝑠 ∈ (0, 2) . It is expected that in the subdiffusive regime ( 𝑠 > 1 ) transport is suppressed due to negative correlations, while for 𝑠 < 1 , transport is enhanced due to positive correlations. Padgett et al. [50] obtained the following 1D series representation of (−Δ) 𝑠 for the interval 𝑠 ∈ (0,2) (−Δ) 𝑠 𝑢(𝑛) = ∑ (𝑢(𝑛) − 𝑢(𝑚))𝐾 𝑠 (𝑛 − 𝑚) 𝑚∈ℤ;𝑚≠𝑛 , (8) where 𝐾 𝑠 (𝑚) = {4 𝑠 Γ(1 2⁄ + 𝑠)√𝜋|−Γ(−𝑠)| ∙ Γ(|𝑚| − 𝑠)Γ(|𝑚| + 1 + 𝑠) , 𝑚 ∈ ℤ\{0} ,0, 𝑚 = 0 . (9) Here 𝑢 is a discrete function on ℤ and Γ is the standard Gamma function. Numerical simulations by Padgett et al . [50] provided confirmation that the representation in (8) and (9) yields enhanced transport (superdiffusion) for 𝑠 ∈ (0,1) and enhanced localization (subdiffusion) for 𝑠 ∈ (1,2) , when compared to the classical case 𝑠 = 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 𝑠 and the disorder concentration 𝑐 (that determines the variables 𝜖 𝑖 ). In Sec. III B, we obtained the value of the dimensionless disorder 𝑐 using equation (5). The value 𝑠 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 𝐻 𝑠,𝜖 is defined, the time evolution of the initial state of the system is generated under the iterative application of 𝐻 𝑠,𝜖 . 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 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 need to move a distance of at least one dust diameter The series representation of the fractional Laplacian has been recently extended to arbitrary order in [80]. 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 {𝛿 𝑖 } of the 1D integer space ℤ as discrete scales proportional to one dust diameter. We let the initial state of the system be given by 𝛿 , located at the origin of the 1D integer space ℤ , and assign to this state a normalized initial energy equal to . The question of interest is how does this energy spread across scales under the iterative application of the Hamiltonian . The time evolution of the initial state 𝛿 under the action of the Hamiltonian 𝐻 𝑠,𝜖 is given by the sequence {𝛿 ,𝐻 𝑠,𝜖 𝛿 , 𝐻 𝑠,𝜖2 𝛿 , … , 𝐻 𝑠,𝜖𝑁 𝛿 } , where 𝑁 ∈ ℕ is the number of timesteps (equivalently, 𝑁 is the number of iterations of 𝐻 𝑠,𝜖 ). Let {𝜑 , 𝜑 , 𝜑 , … , 𝜑 𝑁′ } be the sequence of ℓ (ℤ) orthogonal vector states obtained from Gram-Schmidt orthogonalization of the original sequ ence {𝛿 , 𝐻 𝑠,𝜖 𝛿 , 𝐻 𝑠,𝜖2 𝛿 , … , 𝐻 𝑠,𝜖𝑁 𝛿 } . This step allows for a proper definition of a mathematical distance between subspaces in the Hilbert space. Then, for any nontrivial vector 𝜈 ≠ 𝛿 (representing any spatial scale of interest), we define the distance parameter (mathematical distance) as 𝐷 𝑠,𝜖𝑁 ∶= √1 − ∑ ( 〈𝜈 , 𝜑 𝑘′ 〉‖𝜈‖‖𝜑 𝑘′ ‖) , (10) where 𝜑 𝑘′ is the 𝑘 th term of the sequence {𝜑 , 𝜑 , 𝜑 , … , 𝜑 𝑁′ } . Here, 〈∙,∙〉 is the ℓ (ℤ) inner product and ‖∙‖ = 〈∙,∙〉 . Equation (10) was originally derived in Liaw [81], where results from spectral theory were used to verify the following conjecture: Extended States 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. In other words, if one demonstrates with positive probability that lim
𝑁→∞ 𝐷 𝑠,𝜖𝑁 > 0 , 𝐷 𝑠,𝜖𝑁 ∈ [0,1] (11) then the time-evolved transport behavior of the system under the action of the examined Hamiltonian includes extended energy states. In [50], it was shown that the extended states conjecture also holds for the random fractional discrete Schrödinger operator in Eq. (7). In the following section 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 lim 𝑁→∞ 𝐷 𝑠,𝜖𝑁 > 0 The spectrum of a Hamiltonian 𝐻 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 𝐻 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. 1.2 of [82]). as high probability for energy transport and the existence of ordered structure s at the examined spatial scale, while lim 𝑁→∞ 𝐷 𝑠,𝜖𝑁 = 0 as high probability f or energy dissipation and the existence of disordered (dissipative) structure at the examined scale. B. Energy Across Scales As discussed in the previous section, the two important inputs in the FLS method are the disorder concentration 𝑐 and the fractional power of the Laplacian 𝑠 , which we determined from structural and diffusion analysis of the dusty plasma monolayers at the three pressures of interest. The values listed in Table II were used to define a Hamiltonian appropriate for each case. Next, we calculated the distance parameter from equation (10) for increasing values of the spatial scale, represented by the vector 𝑣 in equation (10). Finally, we need to determine the physically relevant 𝑟𝑎𝑛𝑔𝑒 of nonlocal interactions in the dusty plasma monolayer. 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 [83] that the Debye length assigns a type of large-scale uncertainty in position √〈∆𝑟 〉~𝜆 𝐷 and the positional variance can be expressed as 𝜎 𝑟2 = 〈∆𝑟 〉 = 〈𝑟 〉 − 〈𝑟〉 ≅ 𝑑 + 12 𝜆 𝐷2 , (12) where 𝑑 is the dimension of the space. Thus, the 𝑟𝑎𝑛𝑔𝑒 of nonlocal interactions can be viewed as the number of standard deviations needed to quantify the expected positional variation and can be expressed as 𝑟𝑎𝑛𝑔𝑒 = 𝓃𝜎 𝑟 ≅ 𝓃𝜆 𝐷 , where 𝓃 is the number of standard deviations and we assumed 𝑑 = 1 . We further normalize this value by the size a single dust diameter 𝑑 𝐷 ≈ 10𝜇𝑚 (the smallest relevant spatial scale), which gives the dimensionless value 𝑟𝑎𝑛𝑔𝑒 ∗ ≅ (𝓃𝜆 𝐷 ) 𝑑 𝐷 ⁄ . For the simulation, the assumed Debye length is 𝜆 𝐷 = 0.6𝑚𝑚 , yielding 𝑟𝑎𝑛𝑔𝑒 ∗ ≈ 𝓃60 . In the and cases, 𝜆 𝐷 = 1𝑚𝑚 , resulting in 𝑟𝑎𝑛𝑔𝑒 ∗ ≈ 𝓃100 . The value 𝑟𝑎𝑛𝑔𝑒 ∗ is a numerical truncation, which determines how many neighboring states are considered when computing the action of the fractional Laplacian at each timestep. We note that the application of a fractional Laplacian (−Δ) 𝑠 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 timestep. However, in a dusty plasma monolayer, the effective range of nonlocal interactions has a characteristic length scale given by 𝜆 𝐷 . 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 𝑟𝑎𝑛𝑔𝑒 ∗ = 60, 100 are obtained if only one standard deviation 𝜎 is considered. The corresponding remainder removed by the truncation is 𝑅~ 1 (𝑟𝑎𝑛𝑔𝑒) ~10 −4 ⁄ . As this remainder is small, we expect it to yield a negligible contribution to the ene rgy transport calculation. Based on the diffusion analysis in Sec III C (see Fig 5), although a single 𝜎 𝑟 may be sufficient to describe the positional spread observed for the case, the development of tails in the velocity distributions at lower pressures suggests that a choice of or is more appropriate. Figure 7 shows the time evolution of the distance parameter, calculated for 𝑁 = 3,000 timesteps, which matches the simulation time in the many-body simulations. For each set of conditions, the resulting distance plot is an average of realizations of the same numerical experiment, which minimizes fluctuations due to the random distribution of disorder values 𝜖 𝑖 ∈ [−𝑐 2, 𝑐 2 ⁄⁄ ] in (7). The normalized range of nonlocal interaction in each case is 𝑟𝑎𝑛𝑔𝑒 ∗ = 3𝜎 𝑟 . The reference vector in equation (10) is a linear combination of 𝐿 number of basis vectors in the Hilbert space, with equal weights, and has the form 𝑣 = (1 𝐿⁄ ) ∑ 𝛿 𝑗𝐿𝑗=1 . In this representation, a single base vector corresponds to the minimum relevant scale (approximately equal to one dust diameter). Thus, increasing the number 𝐿 amounts to considering larger spatial scales. In Fig. 7a, corresponding to conditions obtained from the simulation, the distance values remain close to 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 , the probability for transport is slightly decreased at spatial scales 𝐿 ∈ [100, 450] , which are proportional to the range of nonlocal interaction 𝑟𝑎𝑛𝑔𝑒 ∗ = 200 for this case. Nevertheless, at the final timestep considered 𝑁 = 3,000 , the distance values computed for all scales never drop to zero. Instead, the minimum value, which occurs at
𝐿 = [200, 250] , is 𝐷 𝑠,𝜖3,000 ≈ 0.6 . Thus, in the 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. 7c, we see that for , the distance values visibly drop from for scales 𝐿 ∈[100, 500] and approach for scales 𝐿 ∈ [200, 400] (dark blue color in Fig. 7c). For vanishing values of the distance, the spectral method cannot conclude the existenc e of transport at the considered vector scale 𝑣 . 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 𝐿 ∈ [10, 450] (Fig. 7d) 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. V. COMPARISON OF KINETIC AND SPECTRAL APPROACH Since in Fig. 7 the examined dimensionless reference scales
𝐿 ∈ [10, 1000] were normalized by the dust particle diameter 𝑑 𝐷 ≈ 10𝜇𝑚 , the corresponding spatial scales vary in the range to . Comparison to the cell orientation maps of Fig. 4a shows that for the case, scales in the range are smaller than the typical size of the observed crystalline domains (with characteristic 1D scales ≳ 2𝑐𝑚 ). This explains the observed high probability for transport at all scales. Similarly, examination of Fig. 4b suggests that at , 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𝑚𝑚 for the central regions of all the examined monolayers. In Fig. 7b, the FLS method predicts decreased probability for transport at dimensionless scales 𝐿 ∈ [100, 450] , corresponding to spatial range of , which is in agreement with the observation of substructure formation within the larger crystalline domains in the simulations. For the case, Fig. 4c 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. 3 and from the high value of the corresponding dimensionless disorder 𝑐 =
Fig. 7. Time evolution of the distance parameter from equation (10), calculated for increasing reference vector scales for the three pressure cases a) , b) , and c) . The range of nonlocal interactions in a) is , while the range in b) and c) is . Other parameters used for each case can be found in Table II. All distance calculations are averages of realizations. Plots in d) correspond to the region of scales in c) marked by the white dashed rectangle. a) b) c) d) −3 . Fig. 7cd suggest that in the interval 𝐿 ∈ [200, 400] , corresponding to spatial range of , transport is not likely to occurs. As these spatial scales coincide with 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
𝐿 ∈[500, 1000] , or , 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. VI. CONCLUSION
In this article, we present 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 , a disordered crystalline structure at , and a liquid-like structure at . 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 is allowed to vary by ±0.9% . The variation in 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 and 1 𝑃𝑎 , the obtained monolayers are crystalline, while the condition yields a liquid-like structure that develops a self -excited turbulent instability. For each case, we computed a dimensionless disorder 𝑐 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 velocity histograms. While the particles in all three cases are found to move according to anomalous diffusion, only at the observed deviation from classical diffusion results in 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 𝑀𝑆𝐷(𝜏) 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 du sty particles interacting through a shielded Coulomb (Yukawa) potential. The characteristic range of such nonlocal interactions was determined from the Debye length 𝜆 𝐷 selected for each simulation. To study energy transport across 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 𝑠 of the Laplacian is known to asymptotically approach , 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 equation (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 , 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 , 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 , 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. ACKNOWLEDGMENTS This work was supported by the NSF grant numbers 1903450 (EGK, JLP, CDL, and LSM), 1707215 (LSM and TWH), and 1740203 (TWH and LSM), NSF-DMS grant number 1802682 (CDL), and NASA grant number 1571701 (TWH and LSM), and DOE grant number SC0021284 (EGK). Since August 2020, CDL 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. DATA AVAILABILITY The data that support the findings of this study are available from the corresponding author upon reasonable request. Bibliography [1]
J. Lu and R. A. Shaw, “Charged particle dynamics i n turbulence: Theory and direct numerical simulations,”
Phys. Fluids , vol. 27, no. 6, p. 065111, Jun. 2015, doi: 10.1063/1.4922645. [2]
J. K. Eaton and J. R. Fessler, “Preferential concentration of particles by turbulence,”
Int. J. Multiph. Flow , vol. 20, pp. 169 – J. R. Fessler, J. D. Kulick, and J. K. Eaton, “Preferential concentration of heavy particles in a turbulent channel flow,”
Phys. Fluids , vol. 6, no. 11, pp. 3742 – G. E. Thomas and J. Olivero, “Noctilucent clouds as possible indicators of global change in the mesosphere,”
Adv. Space Res. , vol. 28, no. 7, pp. 937 – ohy, “Aircraft measurements of high average charges on cloud drops in layer clouds,” Geophys. Res. Lett. , vol. 31, no. 14, Jul. 2004, doi: 10.1029/2004GL020465. 1 [6]
L. Zhou and B. A. Tinsley, “Production of space charge at the boundaries of layer clouds,”
J. Geophys. Res. Atmospheres , vol. 112, no. D11, Jun. 2007, doi: 10.1029/2006JD007998. [7]
G. Hendrickson, “Electrostatics and gas phase fluidized bed polymerization reactor wall sheeting,”
Chem. Eng. Sci. , vol. 61, no. 4, pp. 1041 – R. G. Rokkam, A. Sowinski, R. O. Fox, P. Mehrani, and M. E. Muhle, “Computational and experimental study of electrostatics in gas –solid polymerization fluidized beds,”
Chem. Eng. Sci. , vol. 92, pp. 146 – R. G. Rokkam, R. O. Fox, and M. E. Muhle, “Computational fluid dynamics and electrostatic modeling of polymerization fluidized- bed reactors,”
Powder Technol. , vol. 203, no. 2, pp. 109 – “Characterisation of charged hydrocarbon sprays for application in combustion systems | SpringerLink.” https://link.springer.com/article/10.1007/s003480050310 (accessed Sep. 24, 2018). [11] F. Esposito, R. Molinaro, C. I. Popa, C.E.S.A.R.E. Molfese, F. Cozzolino, L. Marty, K. Taj‐Eddine , G. Di Achille, G. Franzese, S. Silvestro, and G. G. Ori, “The role o f the atmospheric electric field in the dust- lifting process,”
Geophys. Res. Lett. , vol. 43, no. 10, pp. 5501 – J. F. Kok and N. O. Renno, “The effects of electric forces on dust lifting: Preliminary studies with a numerical model,”
J. Phys. Conf. Ser. , vol. 142, no. 1, p. 012047, 2008, doi: 10.1088/1742-6596/142/1/012047. [13]
A. V. Malm and T. A. Waigh, “Elastic t urbulence in entangled semi-dilute DNA solutions measured with optical coherence tomography velocimetry,”
Sci. Rep. , vol. 7, no. 1, p. 1186, Apr. 2017, doi: 10.1038/s41598-017-01303-4. [14]
J. Mitchell, K. Lyons, A. M. Howe, and A. Clarke, “Viscoelastic po lymer flows and elastic turbulence in three- dimensional porous structures,”
Soft Matter , vol. 12, no. 2, pp. 460 – A. Groisman and V. Steinberg, “Elastic turbulence in a polymer solution flow,”
Nature , vol. 405, no. 6782, pp. 53 –
55, May 2000, doi: 10.1038/35011019. [16]
A. Groisman and V. Steinberg, “Elastic turbulence in curvilinear flows of polymer solutions,”
New J. Phys. , vol. 6, no. 1, p. 29, 2004, doi: 10.1088/1367-2630/6/1/029. [17] S. Bu, J. Yang, Q. Dong, and Q. Wang, “Experimental study of flow transitions in structured packed beds of spheres with electrochemical technique,”
Exp. Therm. Fluid Sci. , vol. 60, pp. 106 – c, “Onset of turbulence in a regular porous medium: An experimental study,” Phys. Fluids , vol. 21, no. 4, p. 045104, Apr. 2009, doi: 10.1063/1.3091944. [19]
J. W. Fox, “Onset of Turbulent Flow in certain Arrays of Particles,”
Proc. Phys. Soc. Sect. B , vol. 62, no. 12, p. 829, 1949, doi: 10.1088/0370-1301/62/12/309. [20]
R. Golestanian and A. Ajdari, “Stochastic low Reynolds number swimmers,”
J. Phys. Condens. Matter , vol. 21, no. 20, p. 204104, 2009, doi: 10.1088/0953-8984/21/20/204104. [21] A. Najafi and R . Golestanian, “Coherent hydrodynamic coupling for stochastic swimmers,”
EPL Europhys. Lett. , vol. 90, no. 6, p. 68003, 2010, doi: 10.1209/0295-5075/90/68003. [22]
J. C. M. Lin and L. L. Pauley, “Low -Reynolds- number separation on an airfoil,”
AIAA J. , vol. 34, no. 8, pp. 1570 – M. S. Selig and J. J. Guglielmo, “High - Lift Low Reynolds Number Airfoil Design,”
J. Aircr. , vol. 34, no. 1, pp. 72 –
79, 1997, doi: 10.2514/2.2137. [24] W. Shyy, Y. Lian, J. Tang, D. Viieru, and H. Liu,
Aerodynamics of Low Reynolds Number Flyers . Cambridge University Press, 2007. [25]
Z. Jane Wang, “Two Dimensional Mechanism for Insect Hovering,”
Phys. Rev. Lett. , vol. 85, no. 10, pp. 2216 – Ellington, “The novel aerodynamics of insect flight: applications to micro - air vehicles,” J. Exp. Biol. , vol. 202, no. 23, pp. 3439 – J. H. Chu and L. I, “Direct observation of Coulomb crystals and liquids in strongly coupled rf dusty plasmas,”
Phys. Rev. Lett. , vol. 72, no. 25, pp. 4009 – H. Thomas, G. E. Morfill, V. Demmel, J. Goree, B. Feuerbacher, and D. Möhlmann, “Plasma Crystal: Coulomb Crystallization in a Dusty Plasma,”
Phys. Rev. Lett. , vol. 73, no. 5, pp. 652 – A. P. Nefedov, O. F. Petrov, V. I. Molotkov, and V. E. Fortov, “Formation of liquidlike and crystalline structures in dusty plasmas,”
J. Exp. Theor. Phys. Lett. , vol. 72, no. 4, pp. 218 – “Crystallization Dynamics of a Single Layer Complex Plasma,” Phys. Rev. Lett. , vol. 105, p. 115004, 2010, doi: 10.1103/PhysRevLett.105.115004. [31] M. Schwabe, S. Zhdanov , C. Räth, D. B. Graves, H. M. Thomas, and G. E. Morfill, “Collective Effects in Vortex Movements in Complex Plasmas,”
Phys. Rev. Lett. , vol. 112, no. 11, p. 115002, Mar. 2014, doi: 10.1103/PhysRevLett.112.115002. [32] G. E. Morfill, M. Rubin-Zuzic, H. Rothermel, A. V. Ivlev, B. A. Klumov, H. M. Thomas, U. Konopka, and V. Steinberg , “Highly Resolved Fluid Flows: ``Liquid Plasmas’’ at the Kinetic Level,”
Phys. Rev. Lett. , vol. 92, no. 17, p. 175004, Apr. 2004, doi: 10.1103/PhysRevLett.92.175004. [33] M. Rubin-
Zuzic, H. M. Thomas, S. K. Zhdanov, and G. E. Morfill, “Circulation’ dynamo in complex plasm a,”
New J. Phys. , vol. 9, no. 2, p. 39, 2007, doi: 10.1088/1367-2630/9/2/039. [34]
M. Klindworth, A. Melzer, A. Piel, and V. A. Schweigert, “Laser -excited intershell rotation of finite
Coulomb clusters in a dusty plasma,”
Phys. Rev. B , vol. 61, no. 12, pp. 8404 – G. Uchida, S. Iizuka, T. Kamimura, and N. Sato, “Generation of two -dimensional dust vortex flows in a direct current discharge plasma,”
Phys. Plasmas , vol. 16, no. 5, p. 053707, May 2009, doi: 10.1063/1.3139252. [36]
V. Nosenko and J. Goree, “Shear Flows and Shear Viscosity in a Two -Dimensional Yukawa System (Dusty Plasma),”
Phys. Rev. Lett. , vol. 93, no. 15, p. 155004, Oct. 2004, doi: 10.1103/PhysRevLett.93.155004. [37] R. Heidemann, S. Zhdanov, K.
R. Sütterlin, H. M. Thomas, and G. E. Morfill, “Shear flow instability at the interface among two streams of a highly dissipative complex plasma,”
EPL Europhys. Lett. , vol. 96, no. 1, p. 15001, 2011, doi: 10.1209/0295-5075/96/15001. [38] A. Gupta, R. Gane sh, and A. Joy, “Kolmogorov flow in two dimensional strongly coupled dusty plasma,”
Phys. Plasmas , vol. 21, no. 7, p. 073707, Jul. 2014, doi: 10.1063/1.4890488. [39]
S. Zhdanov, M. Schwabe, C. Räth, H. M. Thomas, and G. E. Morfill, “Wave turbulence observe d in an auto- oscillating complex (dusty) plasma,”
EPL Europhys. Lett. , vol. 110, no. 3, p. 35001, 2015, doi: 10.1209/0295-5075/110/35001. [40] Y.-Y. Tsai, M.-
C. Chang, and L. I, “Observation of multifractal intermittent dust -acoustic-wave turbulence,”
Phys. Rev. E , vol. 86, no. 4, p. 045402, Oct. 2012, doi: 10.1103/PhysRevE.86.045402. [41]
J. Pramanik, B. M. Veeresha, G. Prasad, A. Sen, and P. K. Kaw, “Experimental observation of dust - acoustic wave turbulence,” Phys. Lett. A , vol. 312, no. 1, pp. 84 –
90, Jun. 2003, doi: 10.1016/S0375-9601(03)00614-5. [42] “The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers,”
Proc. R. Soc. Lond. Ser. Math. Phys. Sci. , Jul. 1991, doi: 10.1098/rspa.1991.0075. 3 [43] H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, “Meso - scale turbulence in living fluids,” Proc. Natl. Acad. Sci. , vol. 109, no. 36, pp. 14308 – Bizonne, F. Raynal, R. Volk, and C. Ybert, “Kolmogorovian Active Turbulence of a Sparse Assembly of Interacting Marangoni Surfers,”
Phys. Rev. X , vol. 10, no. 2, p. 021065, Jun. 2020, doi: 10.1103/PhysRevX.10.021065. [45]
J. K. Meyer, I. Laut, S. K. Zhdanov, V. Nosenko, and H. M. Thomas, “Coupling of Noncrossing Wave
Modes in a Two-
Dimensional Plasma Crystal,”
Phys. Rev. Lett. , vol. 119, no. 25, p. 255001, Dec. 2017, doi: 10.1103/PhysRevLett.119.255001. [46]
C. M. Ticoş, D. Ticoş, and J. D. Williams, “Kinetic effects in a plasma crystal induced by an external electron beam,”
Phys. Plasmas , vol. 26, no. 4, p. 043702, Apr. 2019, doi: 10.1063/1.5092749. [47] V. Nosenko, M. Pustylnik, M. Rubin-Zuzic, A. M. Lipaev, A. V. Zobnin, A. D. Usachev, H. M.Thomas, M. H. Thoma, V. E. Fortov, O. Kononenko, and A. Ovchinin , “Shear flow in a three -dimensional complex plasma in microgravity conditions,”
Phys. Rev. Res. , vol. 2, no. 3, p. 033404, Sep. 2020, doi: 10.1103/PhysRevResearch.2.033404. [48]
C. Liaw, “Approach to the Extended States Conjecture,”
J. Stat. Phys. , vol. 153, no. 6, pp. 1022 – W. King, R. C. Kirby, and C. Liaw, “Delocalization for the 3 -D discrete random Schroedinger operator at weak disorder,”
J. Phys. Math. Theor. , vol. 47, no. 30, p. 305202, Aug. 2014, doi: 10.1088/1751-8113/47/30/305202. [50]
J. L. Padgett, E. G. Kostadinova, C. D. Liaw, K. Busse, L. S. Matthews, and T. W. Hyde, “Anomalous diffusion in one- dimensional disordered systems: a discrete fractional Laplacian method,”
J. Phys. Math. Theor. , vol. 53, no. 13, p. 135205, Mar. 2020, doi: 10.1088/1751-8121/ab7499. [51]
J. K. Olthoff and K. E. Greenberg, “The Gaseous Electronics Conference RF Reference Cell— An Introduction,”
J. Res. Natl. Inst. Stand. Technol. , vol. 100, no. 4, pp. 327 – V. Land, E. Shen, B. Smith, L. Matthews, and T. Hyde, “Experimental and computational characterization of a modified GEC cell for dusty plasma experiments,”
New J. Phys. , vol. 11, no. 6, p. 063024, Jun. 2009, doi: 10.1088/1367-2630/11/6/063024. [53]
H. M. Anderson and S. B. Radovanov, “Dusty Plasma Studies in the Gaseous Electronics Conference Reference Cell,”
J. Res. Natl. Inst. Stand. Technol. , vol. 100, no. 4, pp. 449 – J. B. Pieper, J. Goree, and R. A. Quinn, “Experimental studies of two‐dimensional and three‐dimensional structure in a crystallized dusty plasma,”
J. Vac. Sci. Technol. A , vol. 14, no. 2, pp. 519 – G. Gogia and J. C. Burton, “Emergent Bistability and Switching in a Nonequilibrium Crystal,”
Phys. Rev. Lett. , vol. 119, no. 17, p. 178004, Oct. 2017, doi: 10.1103/PhysRevLett.119.178004. [56] L. S. Matthews, D. L. Sanford, E. G. Kostadinova, K.
S. Ashrafi, E. Guay, and T. W. Hyde, “Dust charging in dynamic ion wakes,”
Phys. Plasmas , vol. 27, no. 2, p. 023703, Feb. 2020, doi: 10.1063/1.5124246. [57]
L. S. Matthews and T. W. Hyde, “Effect of dipole–dipole charge interactions on dust coagulation,”
New J. Phys. , vol. 11, no. 6, p. 063030, 2009, doi: 10.1088/1367-2630/11/6/063030. [58]
L. S. Matthews and T. W. Hyde, “Effects of the charge -dipole interaction on the coagulation of fractal aggregates,”
IEEE Trans. Plasma Sci. , vol. 32, no. 2, pp. 586 – L. S. Matthews and T. W. Hyde, “Charged grains in Saturn’s F - Ring: interaction with Saturn’s magnetic field,”
Adv. Space Res. , vol. 33, no. 12, pp. 2292 – S. Matthews and T. W. Hyde, “Charging and Growth of Fractal Dust Grains,”
IEEE Trans. Plasma Sci. , vol. 36, no. 1, pp. 310 – L. S. Matthews and T. W. Hyde, “Gravitoelectrodynamics in Saturn’s F ring: encoun ters with
Prometheus and Pandora,”
J. Phys. Math. Gen. , vol. 36, no. 22, p. 6207, 2003, doi: 10.1088/0305-4470/36/22/349. [62]
M. Sun, L. S. Matthews, and T. W. Hyde, “Effect of multi -sized dust distribution on local plasma sheath potentials,”
Adv. Space Res. , vol. 38, no. 11, pp. 2575 – L. S. Matthews, K. Qiao, and T. W. Hyde, “Dynamics of a dust crystal with two different size dust species,”
Adv. Space Res. , vol. 38, no. 11, pp. 2564 – K. Qiao, J. Kong, Z. Zhang, L. S. Matthews, and T. W. Hyde, “Mode Couplings and Conversions for Horizontal Dust Particle Pairs in Complex Plasmas,”
IEEE Trans. Plasma Sci. , vol. 41, no. 4, pp. 745 – K. Qiao, J. Kong, E. V. Oeveren, L. S. Matthews, and T. W. Hyde, “Mode couplings and resonance instabilities in dust clusters,”
Phys. Rev. E , vol. 88, no. 4, p. 043103, Oct. 2013, doi: 10.1103/PhysRevE.88.043103. [66] “Phys. Rev. E 48, 168 similar transport in incomplete chaos.” https://journals.aps.org/pre/abstract/10.1103/PhysRevE.48.1683 (accessed Sep. 14, 2018). [67]
T. S. Strickler, T. K. Langin, P. McQuillen, J. Daligault, and T. C. Killian, “Experimental Measurement of Self-
Diffusion in a Strongly Coupled Plasma,”
Phys. Rev. X , vol. 6, no. 2, p. 021021, May 2016, doi: 10.1103/PhysRevX.6.021021. [68]
T. Ott, M. Bonitz, Z. Donkó, and P. Hartmann, “Superdiffusion in quasi -two-dimensional Yukawa liquids,”
Phys. Rev. E , vol. 78, no. 2, p. 026409, Aug. 2008, doi: 10.1103/PhysRevE.78.026409. [69]
B. Liu and J. Goree, “Superdiffusion and Non -Gaussian Statistics in a Driven-Dissipative 2D Dusty
Plasma,”
Phys. Rev. Lett. , vol. 100, no. 5, p. 055003, Feb. 2008, doi: 10.1103/PhysRevLett.100.055003. [70]
B. Liu and J. Goree, “Superdiffusion in two - dimensional Yukawa liquids,” Phys. Rev. E , vol. 75, no. 1, p. 016405, Jan. 2007, doi: 10.1103/PhysRevE.75.016405. [71] L.-
J. Hou, A. Piel, and P. K. Shukla, “Self -Diffusion in 2D Dusty-Plasma Liquids: Numerical-Simulation
Results,”
Phys. Rev. Lett. , vol. 102, no. 8, p. 085002, Feb. 2009, doi: 10.1103/PhysRevLett.102.085002. [72] S. Nunomura, D. Samsonov, S. Zhdanov, and G.
Morfill, “Self -Diffusion in a Liquid Complex
Plasma,”
Phys. Rev. Lett. , vol. 96, no. 1, p. 015003, Jan. 2006, doi: 10.1103/PhysRevLett.96.015003. [73]
T. Ott and M. Bonitz, “Anomalous and Fickian Diffusion in Two - Dimensional Dusty Plasmas,”
Contrib. Plasma Phys. , vol. 49, no. 10, pp. 760 – O. S. Vaulina and S. V. Vladimirov, “Diffusion and dynamics of macro -particles in a complex plasma,”
Phys. Plasmas , vol. 9, no. 3, pp. 835 – T. Ott and M. Bonitz, “Is Diffusion Anomalous in Two - Dimensional Yukawa Liquids?,”
Phys. Rev. Lett. , vol. 103, no. 19, p. 195001, Nov. 2009, doi: 10.1103/PhysRevLett.103.195001. [76]
V. Jakšić and Y. Last, “Simplicity of singular spectrum in An derson- type Hamiltonians,”
Duke Math. J. , vol. 133, no. 1, pp. 185 – E. G. Kostadinova, J. L. Padgett, C. D. Liaw, L. S. Matthews, and T. W. Hyde, “Numerical study of anomalous diffusion of light in sem icrystalline polymer structures,”
Phys. Rev. Res. , vol. 2, no. 4, p. 043375, Dec. 2020, doi: 10.1103/PhysRevResearch.2.043375. [78]
V. V. Uchaikin, “Self -similar anomalous diffusion and Levy- stable laws,”
Phys.-Uspekhi , vol. 46, no. 8, p. 821, 2003, doi: 10.1070/PU2003v046n08ABEH001324. [79]
R. Metzler and J. Klafter, “The random walk’s guide to anomalous diffusion: a fractional dynamics approach,”
Phys. Rep. , vol. 339, no. 1, pp. 1 –
77, Dec. 2000, doi: 10.1016/S0370-1573(00)00070-3. 5 [80] T. F. Jones, E. G.
Kostadinova, J. L. Padgett, and Q. Sheng, “A series representation of the discrete fractional Laplace operator of arbitrary order,”
ArXiv210103629 Cs Math , Jan. 2021, Accessed: Feb. 11, 2021. [Online]. Available: http://arxiv.org/abs/2101.03629. [81] C. Li aw, “Approach to the Extended States Conjecture,”
J Stat Phys , vol. 153, pp. 1022 – D. Hundertmark, “A short introduction to Anderson localization,” in
Analysis and stochastics of growth processes and interface models , Oxford: Oxford University Press, 2008, pp. 194 – G. Livadiotis, “Chapter 5 - Basic Plasma Parameters Described by Kappa Distributions,” in
Kappa Distributions , G. Livadiotis, Ed. Elsevier, 2017, pp. 249 – esh, “Molecular dynamics study of flow past an obstacle in strongly coupled Yukawa liquids,” Phys. Plasmas , vol. 23, no. 12, p. 123703, Dec. 2016, doi: 10.1063/1.4971449.
Appendix A. MANY-BODY SIMULATION OF DUSTY PLASMA 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 𝑚 𝑑 𝒙̈ = 𝑭 𝑑𝑑 + 𝑚 𝑑 𝒈 + 𝑄 𝑑 𝑬 − 𝛽𝒙̇ + 𝜁𝒓(𝑡). (13) Below we discuss the dust-dust interaction force 𝑭 𝑑𝑑 , the interplay between gravity 𝑚 𝑑 𝒈 and confinement forces 𝑄 𝑑 𝑬 , the drag force 𝛽𝒙̇ , and the thermal bath 𝜁𝒓(𝑡) used in the many-body simulations of the present study.
1. Interparticle potential
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 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
𝑉(𝑟) = ( 𝑄 𝑑 𝑟) 𝑒 − 𝑟𝜆 𝐷 , (14) where 𝑟 is the distance from a particle with charge 𝑄 𝑑 and 𝜆 𝐷 is the Debye shielding length (the scale length over which a charged grain is shielded by the plasma). It has been reported [84] that in strongly coupled liquids with Yukawa interactions, the laminar-to-turbulent transition at a fixed Reynolds number 𝑅𝑒 can be induced by increasing the strength of the interparticle potential alone. This behavior has been verified across different ranges of Reynolds numbers (low 𝑅𝑒 < 0.1 , intermediate [2 ≤ 𝑅𝑒 ≤ 35] , and high
𝑅𝑒 > 50 ). Thus, in the present study of active turbulence, it was considered appropriate to derive the dust-dust interaction force 𝐹 𝑑𝑑 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. 2. Confinement forces and gravity The confinement forces in equation (1) / (13) are electrostatic in nature and have the form 𝑭 𝑐𝑜𝑛𝑓 = 𝑚 𝑑 𝒂 𝒄𝒐𝒏𝒇 = 𝑄 𝑑 𝑬 ⟹ 𝑎 = 𝑄 𝑑 𝑚 𝑑 𝑬 , (15) where the radial confinement is derived from the 10 th order polynomial potential 𝑉 𝑖 = 0.5𝑚 𝑑 (𝛺 ℎ2 𝜌 𝑖10 𝑅 ), (16) 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 5 th order polynomial electrostatic field 𝑬 𝑧,𝑐𝑜𝑛𝑓 (𝑧) = [−8083 + 553373𝑧 + (2.0 × 10 )𝑧 − (3.017 × 10 )𝑧 + (1.471 × 10 )𝑧 − (2.306 × 10 )𝑧 ]𝒛̂ , (17) where z is the height above the lower electrode, as shown in Figure 8. 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 5 th order polynomial was obtained from fits to a fluid model code and experimental data, discussed in [52]. 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 𝐸 𝑧 . 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 𝑧 = 0𝑚 , and the top at 𝑧 = 0.0254𝑚 , after which the electric field takes a constant value of . As a result, a dust grain of mass m and negative charge of magnitude 𝑄 𝑑 acquires acceleration 𝒂 𝑧,𝑐𝑜𝑛𝑓 (𝑧) = − 𝑄 𝑑 𝑚 𝑑 𝑬 𝑧,𝑐𝑜𝑛𝑓 (𝑧) = − 𝑄 𝑑 𝑚 𝑑 𝐸 𝑧,𝑐𝑜𝑛𝑓 (𝑧)(−𝒛̂) = 𝑄 𝑑 𝑚 𝑑 𝐸 𝑧,𝑐𝑜𝑛𝑓 (𝑧)𝒛 ̂ . (18) Fig. 8. Fifth order polynomial used for the vertical electric field. The red shaded region is the sheath above the lower electrode, where the electric field is approximately linear. The region shaded blue is the plasma bulk, where the electric field is essentially zero. The region between the bulk and the sheath is the pre-sheath. This electric force is exerted against gravity, so that the electro -gravitational (eg) acceleration in the 𝑧 direction is given in the simulation by 𝒂 𝑧,𝑒𝑔 (𝑧) = 𝑄 𝑑 𝑚 𝑑 𝐸 𝑧,𝑐𝑜𝑛𝑓 (𝑧)𝒛̂ − 9.81𝒛̂ . (19) 3. Gas drag and thermal bath The gas drag force in the many-body simulation is given by 𝑭 𝑑𝑟𝑎𝑔 = −𝛽𝒙̇ , where the damping factor 𝛽 depends on the neutral gas pressure and temperature, with 𝛽 = 𝛿 4𝜋3 (𝑟 𝑑 ) 𝑛 𝑚 𝑔 𝑚 𝑑 √8𝑘 𝐵 𝑇 𝑔 𝜋𝑚 𝑔 . (20) Here 𝛿 = 1.44 (measured for melamine formaldehyde dust in argon), 𝑟 𝑑 is the dust radius, 𝑛 the gas number density, 𝑚 𝑔 the molecular mass of the gas, 𝑇 𝑔 the gas temperature, and 𝑚 𝑑 is the mass of the dust. A thermal bath is realized by subjecting the particles to random force “kicks” , 𝐹 𝑡ℎ = 𝜁𝒓(𝑡) with the maximum acceleration imparted by the amplitude 𝜁 = √2𝛽𝑘 𝐵 𝑇 𝑔 𝑚 𝑑 𝛥𝑡 𝑑 , (21) where 𝛥𝑡 𝑑 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 𝒓(𝑡) 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, Brow nian 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 𝒓(𝑡) , which is appropriate for modeling strongly coupled dust crystals at high gas pressure. As the variance of the Gaussian is increased, bigger random numbers 𝒓(𝑡)𝒓(𝑡)