Anisotropic anomalous diffusion in microgravity dusty plasma. Part One: Nonequilibrium Statistical Analysis

In the absence of long-range interactions or correlations, the particle diffusion can be described by Brownian motion and the model differential equation is the well-known diffusion equation

where $\Delta=\frac{\partial^{2}}{\partial x^{2}}$ is the Laplacian operator, D is the diffusion constant, and p(x,t) is the distribution function. Equation 1 is also the linear Fokker-Planck equation with no drift. The one-dimensional solution to the diffusion equation has a Gaussiona distribution functional form given by

The mean squared displacement (MSD) is the second moment of this equation: $\langle(x-x_{0})^{2}\rangle\equiv\int(x-x_{0})^{2}p(x,t)dx=Dt$. As can be seen, the classical diffusion equation yields MSD that grows linearly with time. Assuming discrete particle positions, the MSD expression for N particles can be generalized in two or three dimensions to

where $d$ is the dimension and $\alpha$ is an exponent that quantifies deviations from a linear dependence in time (i.e., the presence of anomalous diffusion. When $\alpha=1$ (linear MSD), the diffusion is classical, while $\alpha\neq 1$ implies anomalous diffusion. While anomalous diffusion is ubiquitous in nature [40, 41, 42, 43, 14, 30, 44], its mathematical description is challenging. One way to study anomalous diffusion is with a nonlinear diffusion equation, the nonlinear Fokker-Planck differential equation, which has solutions including Porous Medium equations (see page 249 of [62]) and Tsallis q-Gaussian distributions [63] (also known as Kappa distribution or Student’s t distribution). Anomalous diffusion may also be studied by use of fractional spatial derivatives [27, 29, 64, 18] which may have solutions such as Lévy distributions. We note that it is not well understood if universal scaling relations exist that link the different proposed analytical formulations for anomalous diffusion. While there are some direct relations between the nonlinear Fokker-Planck equations and the q-Gaussian solutions, their relation to fractional differential equations is not straightforward even though they can both describe similar regimes of anomalous diffusion.

The nonlinear Fokker-Planck equation with no drift

has the solution

Here $\nu$ is related to the nonextensive parameter $q$ by $\nu\equiv 2-q$ [63]. $A_{q}$ is a normalization factor given by

The parameter $q$ can also be related to the distribution kurtosis $\kappa$, through the following relationship: $k=\frac{15-9q}{7-5q}$, for $q<7/5$ [65] where $k=Kurtosis=\frac{\langle x^{4}\rangle}{\langle x^{2}\rangle^{2}}$. This can be used to quantify the deviations from a Gaussian distribution and the related anomalous diffusion (e.g., as discussed in [66]). The primary interest of the present study is the range $1<q<3$ where the distributions exhibit leptokurtic or ‘fat-tail’ behavior. This regime is also where anomalous diffusion is superdiffusive ($\alpha>1$), though we will also discuss subdiffusive cases where $q<1$ or $\alpha<1$ (see Figure 2 for example distributions representing each regime). Anomalous diffusion is described by taking the second moment of Eq. 5, which yields a non-linear relation between the MSD and time

where $D_{q}$ is a new diffusion coefficient. Nonextensive statistics were proposed by Tsallis [67] and are described fully in the textbook on the subject [63]. The probability distribution function in this formulation is called a q-Gaussian, where the parameter $q$ quantifies the nonextensivity of the system. The nonextensivity feature of the formulation makes it independent of initial conditions and suitable for modeling many-body complex systems with long-range interactions, such as plasma. Systems with this nonextensive behavior no longer have additive entropies, $S_{q}(A+B)\neq S_{q}(A)+S_{q}(B)$, and the different microstates are not in equilibrium. Therefore, quantifying the nonextensivity is needed to make rigorous statistical claims and understand equilibrium properties of systems with long-range interactions and correlations. Thus, $q$ can be understood as a parameter quantifying the strength of correlations or non-local interactions causing the system to move away from equilibrium. The nonextensive parameter $q$ can be found by extremizing the Tsallis entropy $S_{q}$ with Lagrange constraints, such as normalization and momement of the distribution. Essentially, it can be thought of as the nonlinear form of the Maxwellian distribution, with normalization $A_{q}$, kinetic energies $E_{i}$, potential $U$, and inverse temperature $\beta_{q}$. This results in the following probability distribution function, known as the q-Gaussian

Note that in the limit where $q\rightarrow 1$, the Maxwellian distribution is recovered. This is accomplished by letting the potential energy $U_{q}=0$, $\beta_{q}=1/(v_{qth}^{2})$, where $v_{qth}$ is the q-Gaussian thermal velocity, and $E=\frac{1}{2}mv^{2}$, recalling the fact that $e^{x}=\lim_{n\to\infty}\left(1+\frac{x}{n}\right)^{n}$ and letting $n=\frac{1}{1-q}$

The relationship between the Maxwellian temperature $T_{M}$ and the q-Gaussian temperature $T_{q}$ was discussed in [56] and is given by

Eq. 10 relates the q-Gaussian variance fits (i.e. $\sigma^{2}_{q}$, $D_{q}$, $v_{qth}^{2}$, or $T_{q}$) to the corresponding Gaussian or Maxwellian value. This step is necessary as the kurtosis of a distribution affects the variance, which should be taken into account to find a more appropriate variance value.

Figure 2 shows different q-Gaussian distributions scaled to the same height to exemplify how the distribution shape deviates from a Gaussian ($q=1$) for different values of the nonextensive parameter $q$. Equation 5 and Eq. 8 have the same functional form, but have different variables, coefficients, and exponents as the former is in position space, while the latter is in velocity space. In this paper, we distinguish between the two by using $q_{p}$ and $q_{v}$, respectively. The variances of these distributions are finite for $q<5/3$ but diverge for $5/3\leq q<3$. This means that, as $N\rightarrow\infty$, the distribution approaches a Gaussian for $q<5/3$, while, for $5/3\leq q<3$, it approaches a Lévy distribution. For the position displacement distributions, $q_{p}>5/3$ the diffusion is a Lévy process which yields large ’jumps’ or Lévy flights. Since we utilize the notation $q_{v}$, we call a process ’energetic’ or ’suprathermal’ if $q_{v}>\frac{4+d}{2+d}$, where $d$ is the dimensionality. Thus, for $d=1$, $q>\frac{5}{3}$ there are ’large jumps’ in velocity space, as expected for a Lévy process. For $q<\frac{5}{3}$ (implying the mathematical variance is finite), there are no Lévy processes or Lévy flights. In part two of this paper, "Spectral Analysis" we will go into more detail on this topic and the use of the one-dimensional Fractional Laplacian to model the different sub-regimes of anomalous diffusion.

Here we consider PK-4 particle observation camera videos in which the laser sheet was fixed in the central region of the cloud for extended periods of time (about $20~{}s$). The camera frame rate used was $71~{}fps$ resulting in a $0.014~{}s$ time step. The region of interest in the videos was a rectangular section in the center of each cloud with size $14~{}mm$ by $2~{}mm$. The typical number of particles detected per frame was $~{}470$. This yields statistically significant datasets with $20k$ data points or more in each set. The particle positions in the xz-plane were obtained using the open source particle tracking Mosaic Suite plugin of Fiji, which is a distribution of ImageJ [68]. The dust positions were converted from pixels to $\mu$m assuming a pixel resolution of $14.20~{}\mathrm{\mu m}$ [31]. The particle trajectories imported from ImageJ include trajectories that do not start at the same time and do not last for the same duration. As shorter tracks can correspond to particles that move in and out of the plane illuminated by the laser sheet, we discarded tracks where the particles appeare in fewer than 10 frames (or $0.14~{}s$). The remaining trajectory data were used as inputs in the @msdanalyzer code to subtract drift, compute MSDs, and construct histograms of particle displacements and velocities. The nine pressure-current cases considered here use identical particle tracking datasets as those in [61], which investigated the anisotropic structure of the PK-4 clouds using pair correlation analysis. Thus, we expect that the results obtained here and in [61] are directly comparable.

The open-sourced code @msdanalyzer [40] was used to analyze the data. The developers of @msdanalyzer also developed the Mosaic Suite in Fiji, providing good compatibility between the code and the particle tracking data to construct the MSD (with options to provide the standard deviation and check for localization error) and to create the histograms of dust displacements and velocities. The @msdanalyzer code was also used to calculate the velocity autocorrelation, detect directed motion, and subtract drift. Those steps were necessary to ensure that the statistical analysis is performed on the diffusion portion of the dust motion, subtracting convective motion. In addition, motion along the direction of the external electric field $E$ (the x-axis or the $\parallel-$ direction) is considered separately from motion across the the direction of $E$ (the z-axis or the $\perp-$ direction). This was necessary due to the anisotropic coupling of the dust particles discussed in [61].

Each PK-4 dataset used in our analysis was extracted from a large field of view, encompassing a large portion of the dust cloud. As a result, we could break up the field of view in to smaller sections and calculate different drifts in different domains or sections of the cloud. This was necessary as nonhomogeneous drift was observed for several datasets. To address this issue, a method of nonhomogeneous drift subtraction (NHDS) was developed. In this method, the large field of view is separated into 12 smaller regions, as shown in Figure 3 A, and homogeneous drift subtraction is performed for each region. The number of smaller regions was selected to balance accuracy in drift subtraction and reasonable computation times. The NHDS method was validated against cases where the normal homogeneous drift subtraction worked, such as the case at 70 Pa pressure and 0.7 mA dc current. We found the two methods produced the same results, with similar velocity auto-correlation and the total drift plots.

Figure 3 shows the case where the most pronounced nonhomogeneous drift was observed for 70 Pa pressure and 0.35 mA dc current. As can be seen in Fig. 3 B, in the presence of nonhomogeneous drift, applying the drift subtraction to the large region of interest produces particle trajectories that exhibit some coordinated motion (instead of showing random diffusive motion). The NHDS method seems to work better as shown in 3 C, though some small areas still exhibit coordinated motion instead of the desired diffusive behavior. The drifts for each smaller domain region are shown in 3 D. While the total drift simply implied that the dust particles were linearly drifting to the left, the domain drifts show a spread of drift directions.

As a quantitative check of the NHDS application of drift correction for the 70 Pa 0.35 mA case, we calculated the velocity autocorrelation of all domains before and after, shown in Fig. 4. While the velocity autocorrelation before subtraction (red) showed a nonzero value, after the NHDS subtraction (blue), the total autocorrelation dropped to zero. Using these drift-subtracted particle trajectories, we proceeded to find the values of $\alpha$ from a fit to the MHDS plots and the values of $q_{p}$ and $q_{v}$ from fits to the position and velocity histograms. The values of $\alpha$ obtained from fits to the MSD plots at time delay of five seconds are very similar for two drift subtraction methods. However, the application of the NHDS method decreased the magnitude of th measured MSD (reflected in the vertical axis) by an order of magnitude.

In-house codes were developed to obtain nonlinear fits to the MSDs at different time delay intervals $\tau$ and distribution fits to the histograms. Fits to the displacement and velocity histograms were obtained using aMaxwellian (Eq. 9), a single q-Gaussian, and a Bi-q-Gaussian distributions of the form

Here $A,v_{th},v_{0}$, and $q$ are the fitting parameters for the velocity histograms. The fits to the position histograms used the same distributions where $v_{th}$ was replaced by $D$ and $(v_{r}-v_{r0})^{2}$ was replaced by $(r(\tau)-r_{0})^{2}$. To minimize fitting errors, the normalization in Eq. 6 was set to $1$. Dropping this normalization has a negligible effect. For all cases, it was found that the axial component of the velocities (displacements) is best approximated by a single q-Gaussian, while the radial component was best approximated by a Bi-q-Gaussian. The value of $q$ found from the position and velocity histogram fits are denoted as $q_{p}$ and $q_{v}$, respectively. Data obtained from particle displacements/velocities along the axial direction is labeled by a subscript $\parallel$, while cross-field direction data have a subscript $\perp$.

Figure 5 shows the MSD plots obtained from dust positions within a 2D plane for all cases after drift subtraction. A nonhomogeneous drift subtraction was necessary for the two of the cases: (70 Pa, 0.35 mA) and (46 Pa, 1 mA). For all three pressure cases, it seems like highest current (1 mA) leads to the largest slope of the MSD plots, suggesting an enhanced superdiffusive behavior (see dark red, dark blue, and dark green line on the plot). Normally, in the presence of linear drift in the dust motion, the MSD plot will resemblesuperdiffusivity. An increased dc current is expected to result in an increased electric field strength which would enhance the drift of all charged particles. However, in this case, the fast switching of the dc polarity should still prevent dust response. In addition, examination of the velocity autocorrelation suggests that drift has been successfully subtracted from the particles trajectories, leaving only random motion. Thus, it is likely that the enganced superdiffusive trend is related to the effect that dc current has on the ion stream velocity, charge on the dust, and the resulting ion wakefield structure surrounding each dust particle.

To better understand the role of directional anisotropy, we calculated MSDs from particle displacements along $\parallel$ and across $\perp$ the direction of the external electric field. Separate fits were performed at different time delays to assess the role of distinct physical processes. Only data corresponding to time delays smaller than $10s$ was used for the fits. At larger time delays, the standard deviation from the mean MSD increases considerably due to decreasing number of data points. Figure 5 shows representative plots of $MSD_{\parallel}$ and $MSD_{\perp}$ for the 70Pa, 1mA case. For all pressure-current conditions the $MSD_{\parallel}$ has a magnitude much greater than $MSD_{\perp}$ and the plots of $MSD_{\parallel}$ look almost identical to those of the combined MSDs in Fig. 5. The $MSD_{\parallel}$ for all conditions has a positive concavity at all time delays, suggesting superdiffusion. The one exception is the (30 Pa, 0.35 mA) case, which has a brief negative concavity at its start. This dataset had the most noisy trajectories and the least number of data points at long time delays. Thus, the deviation in observed behavior may be due to the poorer quality of the data. The $MSD_{\perp}$ for all cases exhibits a brief negative concavity for a delay period of about two seconds, followed by a positive concavity at lager time scales. This suggests that there may be a trapping mechanism causing sub-diffusive behavior at time scales smaller than $2s$. These effects are showcased in the plots in Fig 6.

Figure 6 is representative of the MSD curves in the parallel and perpendicular directions across all the pressure-current cases. Key features to notice are the larger y-axis scale for the parallel case, the distinct superdiffusive $\alpha>1$ curve for the parallel direction, and subdiffusive $0.1<\tau<~{}1$ to linear shape ($\tau>5$) of the perpendicular direction. These are consistent across all pressure-current cases.

The exponent $\alpha$ extracted from $MSD_{\parallel}$ and $MSD_{\perp}$ at time delay between $4~{}s$ and $5~{}s$ is plotted for all cases in Fig. 7. This time period was chosen as it has the most consistent slope fitting while not being too short of a time delay. Recall that Tsallis’ theory claims that the diffusion is a Lévy process when $q>5/3$. Using the scaling relation $MSD=\langle r^{2}\rangle\propto\tau^{\frac{2}{3-q_{p}}}$ [63] yields the criterion $\alpha>3/2$ for a Lévy process. Note that this scaling is valid only for $q_{p}$ obtained from the displacement distributions. A dashed line on Fig. 7 marks the location of $\alpha=3/2$, indicating that almost all distributions for the parallel cases are a Lévy processes, while the perpendicular direction are not. Later we compare this value of $\alpha$ to that found from $q_{p}$.

Representative histograms of the displacements are shown in Fig. 9 for 70 Pa 0.7 mA case. The histograms of displacements parallel to the electric field are best approximated by a single q-Gaussian distribution, while the histograms of cross-field displacements are best approximated by a Bi-q-Gaussian distribution in most cases. These are displacements of all particles at time delay $\tau=5$. The directional difference of the distributions type implies there is anisotropic diffusion. The nonextensive $q_{p}$ parameters extracted from the displacement histograms are plotted in Figure 8.

Notice that the coefficients for the parallel direction are again larger than for the perpendicular direction. Dashed line indicates a Lévy process. The coefficients extracted from the histogram of parallel displacements suggest that two of the 70Pa cases are consistent with a Lévy process. A crossover to a Lévy process is also seen for the cross-field direction in the 70 Pa, 1mA case if a single q-Gaussian is fitted to the histogram of cross-field displacements. However, the more accurate Bi-q-Gaussian fit to the cross-field displacement histograms suggest that the process is likely superdiffusive but not Lévy. These results are in qualitative agreement with the conclusions drown from the $\alpha$ coefficients extracted from the MSD fits at time delay $\tau=5$ shown in Figure 7. It seems that higher pressures create larger spread in the coefficients and exhibit a larger difference between the parallel and perpendicular directions suggesting a more pronounced anisotropic effect. The values for the diffusion $D\left(\frac{\mu m^{2}}{s}\right)$ calculated from Eq. 10, are shown in Fig. 10. Here we can see that, in the direction parallel to the field, an increase in pressure causes diffusion to increase. However, in the cross-field direction, the interpretation is less straightforward. $D_{\perp 1}$ and $D_{\perp}$ generally decrease with pressure with some exceptions. $D_{\perp 2}$ does not show a clear trend with pressure or current.

Next, we examined how the exponent $\alpha$ found directly from MSD fits compares to $\alpha$ calculated from the displacement histogram fits using the scaling in Eq. 7. Figure 11 shows the percent difference $\frac{\alpha_{MSD}-\alpha_{q_{p}}}{\alpha_{q_{p}}}$. We can see that on average the error is lower, indicated by the dashed line at 20%, for the cases where a Bi-q-Gaussian fit was used rather than a single q-Gaussian. This seems to be in agreement with findings from [30, 14].

Next we examine Fig. 12, which shows representative histograms and distribution fits for the velocity components parallel and perpendicular to the direction of the external electric field. These plots show that a Gaussian or a Maxwellian distribution (black line) does not approximate the ‘fat’ tails of the histograms, while a q-Gaussian (magenta line) and a Bi-q-Gaussian (green line) provide good fits to the $\parallel$ and the $\perp$ velocity histograms, respectively. It is again observed that the distribution features change substantially with direction, suggesting anisotropic dust diffusion caused by the external electric field. Similar to the displacement data, we observe that the velocity histograms in the cross-field direction are best approximated by a Bi-q-Gaussian fit, suggesting a superposition of two distinct diffusion processes. These trends in histogram shapes are observed for all pressure-current cases. This will be further discussed in Sec. 6.2.

Figure 13 shows the parallel and perpendicular nonextensive q parameters obtained from the velocity histograms for each pressure-current case. A slight negative slope is observed for increasing pressure for $q_{v\parallel}$ and none of the coefficients are above the Lévy line. The single q-Gaussian fits for the cross field velocity histograms had an $R^{2}=0.973$ and a normalized root mean square error $NRMSE=0.027$ while the Bi-q-Gaussian fit had $R^{2}=0.993$ and $NRMSE=0.013$. A larger value of $R^{2}$ and lower value of $NRMSE$ imply better fit thus the single q-Gaussian fits for the cross field velocity histograms are not as representative of the system as well as the Bi-q-Gaussian. With the Bi-q-Gaussian fits the $q_{\perp 1}$ values show little change with plasma conditions (most of them are around 1.2), while $q_{\perp}2$ values seem to decrease with increasing pressure, except for the 1 mA case. The parameter $q_{v}$ provides information on the equilibrium of the system, suggesting that in the direction parallel to the electric field, the dust ensemble is close to an equilibrium as seen in Figure13 a), despite the strong superdiffusion observed in Figure 9 a).The trends in the cross-field direction are opposite, suggesting that the identification of the two distinct diffusion patterns drive the system away from equilibrium. This will be further discussed in Sec. 6.3.

The velocity distribution fits were used to extract dust temperatures for the directions parallel and perpendicular to the electric field. Figures 14 and 15 shows the temperatures found using a Maxwellian fit, and a q-Gaussian fit, respectively.

Kinetic temperature is typically defined as the variance of a Maxwellian distribution. Since the velocity histograms for the dust in the PK-4 experiment are best approximated by a non-Maxwellian, the notion of "temperature" here is not well defined. The range of kinetic temperatures found from the Maxwellian fits (Fig. 14) is $2000K-6500K$, or about $0.2eV-0.6eV$, suggesting that the majority of the dust particles have very high kinetic energy even though the experiment is at room temperature and the pressure is low. Anomalously high dust temperature is a known phenomena that is caused by electrostatic fluctuations [69]. In Fig. 15, we used Eq. 10, $T_{q}(\frac{5q-3}{2})=T_{Mq}$ (from [56]), since it provides a relation between the temperature or thermal velocity found from a q-Gaussian and to a ‘Maxwellian-like’ temperature.

As shown in Figure 12 b), the perpendicular velocity histograms are best approximated by a Bi-q-Gaussian fit, where the the bulk distribution (the one with the larger variance and smaller $q_{v}$ value) more closely resembles a Gaussian and, thus, can be used to extract the thermal velocity. The second distribution in the Bi-q-Gaussian accounts for the "fat tails". Thus, in Fig. 15, $T_{\perp 1}$ is extracted from a Gaussian-like sub-population distribution, while $T_{\perp 2}$ is obtained from fitting to a ’tail-halo’ distribution. It may seem unintuitive why the Gaussian sub-population would have higher temperature values, while the ’halo-tails’ have lower temperature values. This will be discussed in 6.2. Essentially we convert a variance from a distribution which has some ’tailed-ness’ to a variance which does not have any, so that we can use the same definition of temperature. We argue that this is not a different treatment of temperature in dusty plasma, but an attempt at calculating temperature more carefully with nonextensive statistics. This allows us to derive more accurate dust temperatures for the PK-4 dusty plasma with values in the range $2500K-5500K$, or about $0.2eV-0.5eV$, which is smaller than the typically calculated values of $~{}10-300eV$ [69].

Finally, for velocity q-distributions, Jiulin and Haining [70, 71] have derived the following expression for $q_{v}$ in terms of a thermophoretic force and Lorentz force

$T_{s}$ is the temperature for each species: electron, ion, or dust. In the PK-4, there is no magnetic field, thus, $\vec{v}\times\vec{B}=0$. In the absence of a temperature gradient, $q_{v}=1$ and the system is in equilibrium, i.e., the velocity histograms should be best described by a Maxwell-Boltzmann distribution. Equation 13 has has not yet been verified against any experimental measurements, though we will give a qualitative argument for its possible credibility.

Here we provide physical arguments why the Bi-q-Gaussian fit for the perpendicular histograms is needed. First, we consider the domain separation on Fig. 3 A) which was made for the NHDS method and reconstruct the velocity histogram in each domain. Figure 16 shows the velocity distribution for domain 11 (dark blue in Fig. 3 A) in both the parallel and perpendicular directions for from the 30 Pa 0.7 mA case.

Although the Maxwellian fit still underpredicts the tailsof the distribution, the single q-Gaussian fits the perpendicular histogram data much more closely than in Figure 12, though the Bi-q-Gaussian is also a reasonable fit. This phenomenon is apparent in all of the other domains and for all pressure cases with current $0.7mA$. Other current cases were not analyzed. We attribute these differences to the transition between local and global dynamics as the system size increases. The nonextensive Tsallis parameter can be used to assess the nonequilibrium state of the system with equilibrium represented by $q=1$ and nonequilibrium by $q\neq 1$. The different local domains in a system may be more or less in equilibrium depending on the properties of the corresponding local distribution functions. As the sample of velocities used to create the histogram is increased to include multiple domains with varying equilibrium properties, the conclusions about global diffusion and thermodynamics may be significantly altered by the averaging process. A smaller domain size is more likely to provide accurate information on the nonextensive parameter $q_{v}$, whether the particles in the small region are in equilibrium or not. Of course, this requires that the number of data points in the smaller domain is sufficient to yield statistically significant results. Table 3 provides the values of all the nonequilibrium $q_{v}$ parameters in each domain with each table entry representing the domain location in the field of view of the particle tracks, (Figure 3 A). This yields a spatial map of the equilibrium properties at different pressure of the PK-4 dust cloud. Table 3 shows that regions near the outer part of the dust cloud are typically further away from equilibrium than regions within the central region.

Nonequilibrium $q_{v\perp}$ in each domain

We use the newly-founded q-Gaussian fits to calculate temperature in each domain using Eq. 10. Figure 17 a) $T_{\|}$ and b) $T_{\perp}$ show temperatures in units of Kelvin for 20 domains calculated for the 30 Pa 0.7 mA case. Each domain contained about 500 dust tracks, a minimum of 10 time frames per trajectory, which yields statistically significant amount of data for the velocity fits. The data indicates that temperatures gradients exist in the PK4 system. As seen from Eq. 13, the temperature gradient and and the gradient in potential are related through the $q_{v}$ parameter. MD simulations of ions and dust in the PK-4, with conditions closely resembling the 70 Pa 0.7 mA case, reveal gradients in the electric potential surrounding each dust particles, as shown in Fig. 18 [38]. Examination of Figures 17 and 18 suggests that the dust clouds in PK-4 exhibit both temperature and electric potential gradients, in addition to q-Gaussian velocity distribution functions. This provides a qualitative argument for applying Eq. 12, which is dependent on the ratio of the thermophoretic force over the electric field force. Providing a quantitative validation for Eq. 13 is beyond the scope of the present paper and will be explored in future work.

When analyzing the MSD plots fitting (Fig. 7), it seems like both $\alpha_{\parallel}$ and $\alpha_{\perp}$ increase with pressure for dc currents 0.35 mA and 0.7 mA. The trends are less obvious at 1 mA. These fits also indicate that none of the extracted $\alpha_{\perp}$ exponents correspond to a Lévy process, while all data for $\alpha_{\parallel}$ suggest a Lévy processes. Meanwhile, the $q_{p}$ coefficients extracted from the position histogram fits (Fig. 9) also suggest that $q_{p\|}$ increases with pressure (i.e., enhanced superdiffusion), in agreement with the $\alpha_{\parallel}$ trends. However, only the highest pressure-currents case (70 Pa 1mA) in Fig. 9 indicates a Lévy processes (i.e., $q>5/3$). The $q_{p\perp}$ values extracted from single q-Gaussian fits show similar trends for the same currents as $\alpha_{\perp}$ but the 70Pa 1mA case suggests a Lévy process, unlike $\alpha_{\perp}$. The Bi-q-Gaussian fits yield $q_{p\perp 2}$ values that exhibit closer similarity to the $\alpha_{\perp}$ trends. We can also see from Figure 11 that the fits for $q_{p\perp 2}$ and $q_{p\perp 1}$ yield lower percent error between the predicted and measured $\alpha_{\perp}$, strengthening the argument that the Bi-q-Gaussian is a more appropriate fit. The calculated diffusion constants (Fig. 10) show an overall increase in parallel diffusion with pressure, while the perpendicular diffusion seems to decrease with pressure when a single q-Gaussian fit is used. These trends are not as clear for $D_{\perp 1}$ and $D_{\perp 2}$ which were obtained from a Bi-q-Gaussian fit. Both the MSD fits and the displacement histogram fits indicate a pronounced difference in the diffusion regime when comparing the direction along the electric field versus the cross-field direction, which can be attributed to an anisotropy in the ion wake-mediated dust-dust interaction potential.

The statistical analysis suggests that at high pressure-current conditions, a transition from superdiffusion to a Lévy process is expected in the axial direction. In the cross-field direction, the Bi-q-Gaussian fits suggest that the observed behavior is a superposition of a classical diffusion and superdiffusion, but no crossover to Lévy process is expected. These observations can be explained when considering the physical mechanisms causing the dust particles in the PK-4 experiments to organize into field-aligned filamentary structures. The alignment and the strong coupling of dust grains along the direction of the electric field is caused by enhancement of the ion wakefield focusing and elongation of the wakefield structure surrounding the dust grains (as discussed in [38]). The formation of positive space charge due to ion focusing will cause a negatively charged dust grain to drift locally towards a nearby concentration of ions. This local drift is a nondeterministic microscopic effect, which is why it affects the microscopic motion by causing pronounced superdiffusion. The elongated ion wakefield structure surrounding the dust grains also causes a restoring force in the cross-field direction that keeps the individual dust particles aligned within the filament. This can explain the subdiffusive (or trapping) trends observed at small time scales in the cross-field MSD plots. Finally, if a dust particle escapes the confining potential that otherwise keeps it within a crystal-like filamentary structure, it cannot easily find force balance in between filaments. Instead, the dust particle makes large-scale jumps across the cloud until force balance is found within another filament. These jumps are frequently observed in video data from PK-4.

Both the position and velocity histograms obtained for the $\perp$ direction (Figure 8, 12) were best fit by a bimodal distribution, i.e., the complex shape of the distribution can be viewed as a superposition of two simpler distributions. The single q-Gaussian which yielded the best results and fewest errors in the parallel direction, does not fit well in the perpendicular direction. We concluded that a Bi-q-Gaussian provided the best fit. Liu and Goree [15] showed that the dust velocity histogram in the PK-4 experiments can be described by what they called a Maxwellian ’core’ and a Kappa distribution ’halo’,

Equation 14 is very similar to the Bi-q-Gaussian as the Maxwellian ’core’ is recovered by setting $q=1$ for one of the q-Gaussian distributions, while the Kappa ’halo’ is identical to the q-distribution when one uses the known scaling relation $q=1+1/\kappa$ in the other q-Gaussian. In addition, Liu and Goree concluded that the cloud exhibited liquid-like coupling as evidenced by the calculated pair correlation function. The dust cloud in these experiments was confined in rf plasma produced by an rf coil and the observed cloud structure was homogeneous. In contrast, the present experiments were conducted in a pure dc discharge and the observed cloud structure was found to exhibit anisotropic coupling, suggesting both liquid and crystalline properties [61]. In both cases, the distributions consist of two terms, which implies a superposition of two distinct diffusion processes. Figure 19 shows a decomposition of the Bi-q-Gaussian into individual q-Gaussians.

The superposition of the two q-Gaussian distributions shown in Fig. 19 d) shows that one diffusion process (in red) is characterized by a small peak but large variance, while the other one (in blue) has both a large peak and large tails. The variance (and the corresponding temperature) of the blue distribution, however, is much smaller. Thus, the observed diffusion in the cross-field direction can be viewed as a superposition of two processes - one with higher temperature close to equilibrium and a second one with a lower temperature but farther from equilibrium.

This distribution shape is also called a ’knee and ankle’ distribution, shown in Figure 20.

The ’knee and ankle’ distribution has also been observed in solar wind ion and electron velocities and magnetic field fluctuations $\delta B$ [72, 73, 74] as well as in fluxes of cosmic rays [75]. In [74], a bi-kappa distribution was proposed to explain the scale-dependent changes in the solar wind PDFs obtained from spacecraft measurements. It was pointed out that small scale PDFs are highly non-Maxwellian, while Maxwellian is recovered for large scales. In PK-4, the dust particles seem to exhibit similar disparity between small and large scales as evident from the distribution fits. As the dust diffusion is highly dependent on ion-wake-mediated interaction potential, small time fluctuations in the ion wakefield can be the reason for the highly non-Maxwellian features in the PDFs. Such fluctuations are expected to result from high-frequency ionization waves that were recently discovered in PK-4 [37].

In Space Physics, it is common to describe the non-Maxwellian velocity distributions using the kappa-distribution which is related to the q-Gaussian by $\kappa=\frac{1}{q-1}$. In [73], three physical mechanisms for generating a Kappa distribution were proposed: (i) Debye shielding, (ii) magnetic field biding, and (iii) temperature fluctuations. It was proposed that $\kappa$ has a negative correlation with the Debye number $N_{D}$ (number of particles in the Debye sphere). As $\kappa$ is inversely related to $q$, we expect a positive correlation between $N_{D}$ and $q$. In other words, inclreasing number of ions in the Debye spheres surrounding the dust should lead to enhancement of the superdiffusion. It was also proposed in [73] that $\kappa$ is positively correlated with magnetic field fluctuations $\delta B$, which would imply negative correlation between $q$ and $\delta B$. While there is no magnetic field in PK-4, the charged dust particles are aligned to the external electric field. Thus, we conjecture that increase in electric field fluctuations should randomize the dust motion, thus, leading to classical diffusion. Finally, the negative correlation between $\kappa$ and the magnitude temperature fluctuations $\delta T^{2}$ proposed in [73] suggests a positive correlation between $q$ and $\delta T^{2}$. The role of temperature and temperature flucttuations in PK-4 will be further discussed int he next section. According to [76], the magnetic field fluctuations, $\delta B$, nearer the Sun around $6.9$–$9.7$ AU, or even $43.6$–$47.2$ AU, remain in a non-Gaussian, meta-equilibrium state. However, farther out between $80$–$87$ AU, where the system has had more time ($\approx 1$ year) to relax, the magnetic field fluctuations reach a Gaussian distribution suggesting an equilibrium in magnetic field fluctuations. In other words, at these distances the solar-wind relaxes back to an equilibrium after going through a transition in the diffusion regime from $q>5/3$ to $q<5/3$ between $\sim 47$ AU and $\sim 80$ AU. Similar nonequilibrium-to-equilibrium regime transitions can be explored more easily and cost-effectively in laboratory settings with a dusty plasma experiment.

As shown in Fig. 15, temperature decreases with increasing pressure for $T_{\|}$ and $T_{\perp}1$, aligning with the conclusions presented in [61], where it was seen that dust ’cooled off’ due to the increased role in dust neutral collision, which was calculated by use of pair-correlation functions. The temperatures $T_{\perp}2$, however, seem to show the opposite trend of increasing with increased pressure. The Tsallis coefficients $q_{v\parallel}$, never exhibit values greater than $5/3$, which suggests that there is no Lévy processes for this direction. We also notice that $q_{v\parallel}$ tends to decrease with increasing pressure (thus, improving equilibrium), while increasing current slightly increases $q_{v\parallel}$ (thus, driving the system away from equilibrium). This is reasonable as an increase in pressure causes an increase in neutral collisions, meaning that the random collision process brings the dust motion to equilibrium. Meanwhile, the current provides an electric field which affects the ion streaming and focusing surrounding the dust, which can cause local attractive interactions among dust particles and macroscopic regions of positive space charge.

While the coefficients $q_{v\perp}$ and $q_{v\perp 2}$ suggest Lévy processes, the trends are not as simple to interpret. To visualize this, we look at Figure 19 a), which shows a hot population, more so in equilibrium, while the cold, shown in blue, has much lower kinetic energy. This is most likely due to the trapping potential that creates the observed filamentary structures within the dust cloud. However, the blue population is far from equilibrium. This population is confined in the cross-field direction and much of its kinetic energy is in the parallel direction. When this energy transitions to the cross-field direction, it may do so suddenly, resulting in the larger tails of the distribution, or big jumps of the dust particles in space. In other words, Figure 19 illustrates a one population that is highly peaked with significant tails (kurtosis), giving a small variance, i.e., temperature, and another distinct population that exhibits high variance but less kurtosis. Based on this interpretation, increasing pressure appears to drive the first population $q_{\perp 1}$ away from equilibrium, while improving equilibrium for the second population $q_{\perp 2}$ for all but the 1 mA case. The confining forces keep the dust particles in a crystal-like structure, mostly aligned in the parallel direction, but these forces do not always succeed in maintaining this configuration. Similar to slipping out between tight tweezers, particles that escape these confining forces become energetic in the perpendicular direction. Once the particles escape, their movement is primarily dictated by the ion wakefield caused by the electric field, leading to high parallel diffusion. In other words, the forces confining the dust must act in the parallel direction; otherwise, the particles would accelerate due to the electric field.

An examination of the domain velocity distributions at 0.7 mA reveals that $q_{v\perp}$ tends to be higher on the edges of the cloud and smaller in the center at low pressure, but larger on the top and smaller on the bottom for 46 Pa, 70 Pa. We also see that for 30 Pa at 0.7 mA there exist clear gradients in both the $T_{\parallel}$ and $T_{\perp}$. (We remind the reader that there is only one $T_{\perp}$ here since a single q-Gaussian fits well the individual domain data in the cross-field direction.) Unexpected high dust temperature has been a point of contention in the dusty plasma community because the experiments are always in a room temperature setting, however this is better understood when considering all forms of energy in strongly coupled systems. We propose that, in addition to the kinetic temperature of the dust, the additional apparent energy comes from variation of the electrostatic potential. Typically one calculates temperature as $E=k_{B}T=1/2mv_{th}^{2}$. However, we propose $E+U=k_{B}T+Q_{dust}\langle\Delta\phi_{float}^{2}\rangle$, meaning that room temperature of the dust particles can be recovered by employing the following expression

An order of magnitude calculation using this equation reveals that variation in the floating electric potential on the order of $~{}10^{-5}$ Volts could explain temperature discrepancy. This conclusion is explained more rigorously in [69].