Observation of Kolmogorov turbulence due to multiscale vortices in dusty plasma experiments

We conducted multiple campaigns to study flows and turbulence to ensure the reproducibility and consistency of observed dust rotation phenomena. A total of six independent experimental campaigns were conducted at the same working pressure of 0.12 mbar, as listed in Table 1. All of these campaigns have been performed using polydispersive Kaolin microparticles. The schematic in Fig. 2 shows tentative dust rotation locations from different campaigns. These campaigns successfully captured the X-Z plane of the dust cloud, ensuring that the laser sheet remained at the center of the glass flange at location $y=0$. The hemispherical, shaded region represents the dust void. In each campaign, we observed dust rotation with variations in flow characteristics. Specifically, I, II, III, and VI form in reasonably different plasma environments. We anticipate an interplay of different sets of forces Kaur et al. (2015b); Chai and Bellan (2016); Vaulina et al. (2003a); Choudhary et al. (2018); Samarian et al. (2002) to be responsible for individual vortices in the extended diffused region. It is interesting to note that in all these conditions, the nature of turbulence almost follows the standard forced Kolmogorov theory.

Table 2 enlists parameters for the background plasmas and the dust cloud. The electron temperature, $T_{e}$, and the plasma density, $n_{e,i}$ (in quasineutral plasma), were measured using the single Langmuir probe above the electrode and in the diffused plasma region in the absence of dust Bose et al. (2017). Over the cathode, $T_{e}$ ranges $1.5-4$ eV from the center to the edge. In the diffused region where dusty plasma typically forms, $T_{e}$ is relatively high, ranging from 3 - 7 eV at different locations within the domain span. The plasma density varies by up to two orders of magnitude from the electrode to the diffused region of interest, as shown in the table.

We used a microscope to obtain the kaolin particles’ size distribution. Most particles ranged from 5 to 8 $\mu$m, although there were a few outliers with sizes as large as 50 $\mu$m. Based on the known mass density of kaolin used in the experiments, the typical mass of dust particles ranged from $1.6\times 10^{-13}$ to $6.7\times 10^{-13}$ kg. Out of the six campaigns, we could track particle positions with reasonable accuracy using the ImageJ software Schindelin et al. (2012) for two campaigns. We used the particle data from Fig. 3(b) for 145 successive frames and a particle tracking algorithm to gather particle velocities (in mm/s). Determining the accurate average interparticle separation, $(a)$, for dust clouds composed of particles of varying sizes and experiencing different spatial plasma environments is challenging. However, we calculated “$a$” from regions near vortex motion with visually no flow but only thermal motion. We obtained an estimated value of $a\approx 224~{}\mu$m and a typical dust density, $n_{d}\approx 2.3\times 10^{4}$ cm^-3 in the diffused plasma region. Figure 3(b) depicts one such particle arrangement (from Campaign V) we used to estimate $a$.

Dust particles typically have a charge between [2 - 5] $\times 10^{4}$$e$. The floating potential on a spherical dust particle is the basis for the collisionless orbital motion limited model, which is used to compute the charge. We fitted the theoretical Gaussian velocity distribution to the Probability Distribution Function (PDF) of the particle velocity data (averaged PDF from 145 frames), which we plotted in Fig. 3(a). A theoretical fitting in comparison with experimental data yields dust temperature ranges of 0.6 - 2.2 eV. We see that in a three-dimensional situation, the dust temperature will vary since it is computed using a two-dimensional velocity field. For two datasets (II and VI), we additionally computed the 2D radial distribution function (RDF). The dynamics are thermal, and there is no bulk flow in the area used for this computation. With a tiny initial peak, the RDF in Fig. 4 illustrates the equilibrium dust dynamics in a somewhat liquid regime.

Image processing and PIV are crucial for accurately extracting fluid velocity in the dust cloud. We detail the processes adopted from recorded videos to ascertain velocities over a square grid. We processed video data from six campaigns (I to VI) to study fluid velocities and turbulent mixing. Campaigns I and II were recorded at 25 fps, III and IV at 101 fps, and V and VI at 125 fps. All original images were pre-processed using the open-source ImageJ software by converting them to grayscale. Thresholding was then applied to convert all pixels with an intensity range of 50-220 to maximum brightness, while other pixel intensities were set to zero. A convolution filter sharpened the particles for better clarity, especially given the significant particle-size asymmetry. This data was then ready to export for PIV studies.

We used MATLAB-based PIVlab software for the PIV analysis of the pre-processed image data Thielicke and Sonntag (2021). After importing the images, we selected the square region of interest (RoI) to study vortex and turbulence flows. Unwanted regions (blank spaces with no dust) were removed by masking. For the PIV analysis, we used a $32\times 32$ pixel interrogation area with a 75% overlap of windows. The FFT-based cross-correlation algorithm calculated velocities (pixels/frame) for all windows within the RoI. The velocity data was then calibrated based on actual experimental dimensions, using pre-marked references in the images for length and fps for timescales. Unwanted high-velocity data significantly higher than the bulk flow speed were removed using velocity histogram cropping. Finally, we obtained N-1 sets of velocity field data from N images and used them to profile velocity fields, vorticity, and energy spectra.

The dust cloud in the diffused plasma region showed a variety of collective dynamics as we varied the discharge voltage between 230 and 600 volts. As listed in Table 1, the self-sustained dust rotation has been observed for values at both edges of the voltage range. Figures 5 and 6 show time-elapsed images for different campaigns where dust rotation is visible from particle trajectories. The inset in Fig. 5(I) shows the tracking of selected particles in the box region. The particle tracking was performed using the ImageJ software’s tracking algorithm. The exact location of each dust cloud with reference to the center of the cathode (in the Z-direction) and the surface of the glass (in the X-direction) is mentioned as axis labels. All other dust rotations are observed typically at the same location and similar discharge parameters except for campaigns I and II. For each campaign, the vortex motion is quite stable, and dust continues to rotate until we change the plasma creation conditions. We have provided dust rotation movies for each campaign as supplementary material mov .

We understand that a typical dust rotation in an unmagnetized dusty plasma may arise due to counteracting forces at spatially different locations (Bose et al., 2019). While gravity is always one force acting downwards, the radially outwards ion drag force in the void, the spatially varying force due to the sheath electric field are few possibilities in the counter direction, or at least their one component is opposite to the gravity.

A comprehensive analysis to establish the cause of dust rotation is beyond the scope of this paper. However, we ruled out some of the possibilities. A temperature gradient in neutral gas can give rise to thermophoretic force, which can affect dust dynamics (Shukla and Mamun, 2015). Typically, in glow discharge, the temperature of the cathode can increase reasonably due to ion bombardment after a few hours of operation (Kaur et al., 2015b). In our experiment, we observed the dust vortex soon after reaching the operating discharge voltage range before the cathode was reasonably heated. This suggests that thermophoretic force may not influence the vortex dynamics significantly. We also ruled out the possibility of a neutral flow in the vacuum chamber affecting the dust dynamics. We repeated the vortex experiment by simultaneously closing the pump and the Ar gas inlet to minimize any kind of neutral flow in the vacuum chamber. In this condition, too, we observed the dust vortex, thus ruling out any role of neutral flow.

In earlier dusty plasma experiments with vortex observation, the causes involved variations in microparticle sizes (Yokota and Honda, 1996; Rubin-Zuzic et al., 2007), charge gradients (Vaulina et al., 2000, 2003b), or a non-zero curl of the forces exerted by the plasma on the microparticles (Goedheer and Akdim, 2003; Akdim and Goedheer, 2003; Kaur et al., 2015c). In the future, we plan to study the force balance comprehensively in our system. The goal of this paper is limited to demonstrating turbulence in a dusty plasma medium.

Figures 7 and 8 provide flow velocity vectors and vorticity plots for all campaigns. The vorticity plots were obtained after taking the mean of velocity field profiles from 100 frames, which profoundly helped to visualize the large-scale rotations. The averaging cleans up the fluctuations and reflects the large-scale vortex structures. However, these fluctuations are important when evaluating the energy spectrum in the next section. The vorticity profiles of all campaigns reflect that the dust rotation is clockwise for campaigns III-VI as shown in Fig. 8, and the same is counterclockwise for campaigns I and II (see Fig. 7). This is observable from the velocity flow vectors and the colors of the vorticity profiles in different plots.

An important feature of turbulent flows is that they exhibit multi-scale energy cascades, which are typically demonstrated through the spatial energy spectrum (Kolmogorov, 1991a; Frisch and Kolmogorov, ). Figure 9 presents subplots displaying the energy spectra for each campaign, which demonstrate the distribution of kinetic energy across different scales or wavenumbers. The spectra are calculated in the following way: For each recorded dataset, we take the PIV and get multiple velocity-flow profiles for each campaign. Once we have a coarse-grained flow velocity ${\bf u(x},t)$ for a given snapshot, we perform a Fourier transform on ${\bf u(x},t)$ to obtain ${\bf u(k},t)$, from which we calculate the modal energy $E({\bf k},t)=|{\bf u(k},t)|^{2}/2$. After that, we find the one-dimensional shell spectrum $E(k,t)$ by adding up the modal energies within a certain range $k-1\leq k^{\prime}\leq k$, as described by Frisch and Kolmogorov .

We estimated the energy spectra for each flow profile of the recorded datasets separately. Then, we plotted the average energy spectra over several time frames with error bars, as seen in Fig 9 for all the campaigns. The wavevectors are normalized with $\mathrm{k_{max}}$, where $\mathrm{k_{max}}$ is the maximum wavenumber associated with the system. The values of $\mathrm{k_{max}}$ for each campaign are $6.334$, $5.69$, $5.69$, $5.69$, $4.91$, and $4.91$ $\mathrm{mm^{-1}}$, respectively. Note that the units of $\mathrm{k_{max}}$ are $\mathrm{mm^{-1}}$, and the units of the energy spectra are $\mathrm{(mm/s)^{2}}$. We perform normalization to present the spectra on the same scale without losing generality.

We observe that the energy fed at the large scales cascades to the smaller scales following the standard -5/3 scaling for all campaigns, which is consistent with Kolmogorov’s theory (Kolmogorov, 1991b). The corresponding scaling in the real space is observed by analyzing the second-order structure function discussed in section III.3. However, note that the wavevector ranges vary slightly for different campaigns. Also, the total energy range in the spectra varies up to an order of magnitude as we look for the spectra for I-VI. It is based on the vortex domain size and finite flow velocities. For example, in Fig. 7, the vortex size of I is smaller than that of II; hence, its energy spectrum has a larger magnitude range.

The flow dynamics is three-dimensional, but the present analysis is carried out only in the x-z plane due to the limitation of recording only one 2D plane using a laser sheet and a camera. However, we have individual movies of the y-z and x-z planes (asynchronous) for campaigns V and VI. Here, we provide in Fig. 14 (see the Appendix), the energy spectrum and the vorticity plots for campaign V in the y-z plane. These plots demonstrate the 3-D characteristics of turbulent vortex flow in the cloud at the specific position of $x=5$ mm. The least-square fitted scaling (the slope) obtained is $-1.595\pm 0.099$, consistent with our claims of Kolmogorov scaling within experimental uncertainties.

In this study, we calculate the second-order structure function for all data sets by determining the difference in velocity between a reference point $r$ and varying separation distances $l$, as shown in Fig. 10. To maintain generality, we computed the structural functions using a time average of 90 velocity frames. The fitted slopes yield the scaling exponent close to $0.66$ with standard deviation, which is almost consistent with Kolmogorov theory (Kraichnan, 1974). We also see that the $S_{2}(l)$ values for campaigns III and IV are two orders of magnitude lower. This is because the velocity order is lower in these vortex regimes, which is also shown in the energy spectra, Fig 9.

Using analytical arguments, we also estimated the Kolmogorov dissipation length and time scales. These scales were calculated based on the mean energy transfer rate $\epsilon\sim v_{l}^{3}/l$, where $v_{l}=2~{}\text{mm/s}$ is the average velocity and $l=a$ is the average interparticle separation. The local viscosity is $\nu\sim v_{l}l$. With $\epsilon$ and $\nu$, we computed the Kolmogorov length scale as $l_{d}=(\nu^{3}/\epsilon)^{1/4}=110~{}\mu\text{m}$ and the velocity scale as $v_{d}=(\nu\epsilon)^{1/4}=2~{}\text{mm/s}$. Consequently, the dissipative time scale is $\tau_{d}=l_{d}/v_{d}=0.0478~{}\text{s}$, leading to a dissipation frequency of approximately $21~{}\text{Hz}$. This is also clear from Fig. 11, which shows the dissipation frequency at almost $21~{}\mathrm{Hz}$.

Since the turbulent velocity fields exhibit significant fluctuations, it’s essential to comprehend their statistical properties. Consequently, researchers have focused on studying velocity fields’ PDFs and their derivatives. Velocity PDFs provide insights into the distribution of flow velocities at different points in space (Pandit et al., 2009; Andrés and Banerjee, 2019). The PDF of velocity fields should follow a Gaussian distribution in fully developed turbulent flows, but their derivatives often deviate from the Gaussian distribution, especially at the tails (Gulitski et al., 2007).

Figure 12 presents the PDFs of velocity gradients, $\delta v_{z}$ and $\delta v_{x}$, which deviate from the Gaussian distribution (shown by the black dashed lines) at the tails, indicating intermittency. Note that the PDFs shown are averaged over 90 datasets for each campaign. The fitted plots are obtained by providing $\delta v_{x}$ and $\delta v_{z}$ dataset to the standard MATLAB function “fitdist” for “normal” distribution. The distribution returns a mean value and the associated standard deviation in each case. Table S2 detailing mean velocity, standard deviation, skewness, and kurtosis is included as supplementary material.