Diffusion coefficients of electrorheological complex (dusty) plasmas

Equilibrium molecular dynamics (EMD) simulations have been executed to investigate the parallel (D║) and perpendicular (D┴) diffusion coefficients for three-dimensional (3D) strongly coupled (SC) electrorheological complex (dusty) plasmas (ERCPs). The effects of uniaxial (z-axis) AC electric field (MT) on dust grains have been investigated along with various combinations of plasma parameters (Γ, κ). The new outcomes obtained by mean squared displacement of Einstein relation show diffusion coefficients for low-intermediate to high plasma couplings (Γ) for varying MT. The D║ and D┴ at MT = 0.01 agree well with earlier available data obtained from the Green–Kubo and Einstein relation for 3D SC-Yukawa systems. The simulation data show that D║ increased with an increase in moderate MT strength and that D┴ decreased for the intermediate to large MT strength. Both (D║ and D┴) remained nearly constant for low MT values. The investigations show that the current EMD scheme is more efficient for nonideal gas-like, liquid-like, and solid-like states of SC-ERCPs. It has been demonstrated that present simulation outcomes extended the MT range up to 0.01 ≤ MT ≤ 10 to understand the diffusive and rheological behaviors of dusty plasma systems.


Introduction
Ionized gases containing electrons, ions, and neutral and micron-sized charged dust grains are called dusty or complex plasmas (CPs). Plasma species, including dust grains, show a collective behavior that characterizes a complex system. Micron-sized dust grains are charged under typical laboratory conditions. The amount and types of charge depend on charging phenomena. Mostly, a dust grain has a negative charge due to the high mobility of electrons [1][2][3][4][5][6]. The minimal charge-to-mass ratio of dust grains establishes that these particles are coupled and arrange themselves the same as an atom in the liquid and solid-like phases [1,7]. These dust grains significantly influence the constituent's composition and dust-free plasma equilibrium, resulting in a substantial alteration in thermophysical properties. CPs became an interesting research field after observing dust Coulomb crystals in a laboratory experiment [7]. CPs have been well investigated concerning fundamental aspects in space and industrial applications. Dust grains are naturally present in the Earth's atmosphere, space, and laboratories. CPs in astrophysical conditions include comets, interstellar clouds, and interplanetary liquid in the white dwarf, neutron stars, supernova cores, etc. DC discharge and radiofrequency discharge plasmas are typical examples of laboratory low-temperature CPs. Dust grains can be insert manually or grown in a laboratory through sputtering and chemical reactions. CPs have several applications in science and engineering, such as semiconductors processing, plasma display, filtration process, paper and petroleum industries, biomechanics, insulation systems, and ceramic process. They also play an important role in nuclear fusion devices [1,5,8].
Conventional electrorheological (ER) fluids contain nonconducting particles, drastically changing their properties when an external electric field is applied [9][10][11][12]. ER fluids have various applications in industries, such as vibration control in smart materials [10], ultrasonic transmissions [11], and sound transmission with small losses [12]. In the presence of the external DC or AC modes of an electric field, CPs also show electrorheological behaviors, so-called ER-CPs (ERCPs) [13][14][15][16][17][18][19]. Under the action of the external electric field, the attractive forces between charged dust grains are introduced through ion streaming (wake potential). Due to attractive forces, dust grains are polarized and arranged as strings, sheets, or crystalline structures same as ER fluids. In the case of the DC electric field, the system is non-Hamiltonian and induced wake potential is asymmetric. The AC electric field (frequency of field is much smaller than ion frequency and larger than dust grain frequency in plasma) system is Hamiltonian, and induced wake potential is symmetric. In the presence of an external electric field, the phase transition phenomena occurs in CPs [14][15][16][17][18][19]. During the phase transition, phenomenon transport properties of CPs should change. ERCPs were discovered by Ivlevet al [14] through microgravity experiments (PK-3 DC discharge) on the board of the international space station (ISS) and molecular dynamics (MD) simulations where they observed a phase transition with increasing AC electric field [14,15]. ER and magnetorheological fluids have significant applications on an industrial scale, such as in medical engineering, hydraulics, or photonics. In the future, for precise tailoring, new smart materials may be developed with the help of ERCPs. ERCPs can play an essential role in the manufacturing of smart materials [15].
The CP research has critical issues for understanding the effects of external fields such as electric or magnetic on charged dust grains in a typical plasma background. The external electric field significantly alters different properties of CPs, such as static and dynamic, structure, shape of Debye screening (Debye length), waves and instabilities, and phase transition [13][14][15][16][17][18][19]. The external electric field in CP research directions has opened new dimensions of work. The effect of the external electric field has widely been investigated in various physical systems [39] and ionic liquids [40,41]. It is understood that the diffusion of dust grains are connected with the phase transition process under the influence of external electric and magnetic fields in CPs. Very little literature is available on the effects of an external magnetic field [22][23][24] and an electric field [25,26] on diffusive characteristics in CPs. Parallel and perpendicular diffusion coefficients were calculated using the EF-MD simulation by Ott and Bonitz [22] and Dutta et al., [23], and GK by Begum and Das [24] in magnetized plasmas. Recently Shakoori et al. used GK-MD to investigate the effect of an external electric field on diffusion in CPs [25,26].
However, the literature demonstrates that the self-diffusion coefficients of CPs were extensively investigated using both formulas (EF and GK) by employing theoretical techniques, computer simulation methods, and experiments without an external electric field [20,21,[27][28][29][30][31][32][33][34][35][36][37]. To the best of our knowledge, no data reported previously related self-diffusion coefficients under the influence of an external electric field except in two references where they used a uniaxial AC electric field [25] and a homogeneous external electric field [26] with the GK-EMD method and reported a preliminary data. In the present study, we investigated parallel and perpendicular self-diffusion coefficients to the external uniaxial AC electric field using the EF-EMD technique for a wide range of plasma Coulomb coupling (Γ) and Debye screening (κ) parameters and electric field strength.
The major goal of the current work is to investigate the normalized parallel and perpendicular diffusion coefficients in SC-ERCPs by using the EF-EMD technique under varying uniaxial AC electric fields (M T ) for various plasma states (Γ, κ). This work provides accurate data of diffusion coefficients through an EF-EMD method with computational time and cost-effectiveness in SC-ERCPs and can be used in other related physical systems. Recently, we have reported some preliminary results for SC-ERCPs [25], where we have published results with varying M T using GK-EMD for small N and selected a range of intermediate to high Γ. This paper facilitates a review of the present picture of diffusion coefficients in SC-CPs and new simulation results of parallel and perpendicular diffusion coefficients for SC-ERCPs for the first time. Moreover, the simulation outcomes demonstrate that ERCPs under the action of an electric field behave the same as conventional ER fluids regarding the diffusion of dust particles, indicating phase transitions. The presented EF-EMD scheme can help to understand the microscopic diffusive phenomenon and linear behavior of ERCPs. This EMD technique is used to study the algorithm's performance and compare the outcomes to those obtained through the MD method in 3D CPs and OCPs.

Numerical model for MD simulation
Nowadays, the art of MD simulations has become a common computational tool employed to investigate physical properties in various areas of science and technology. MD simulations follow the equation of motion for spherical particle systems [2][3][4][5][6]. We have used a homemade code based on MD simulations [42] to calculate self-diffusion coefficients (D) with the Einstein relation. Equilibrium MD (EMD) simulations are used in the microcanonical ensemble to examine the diffusion coefficients in SC-ERCPs. The Yukawa (screened Coulomb) potential is used to calculate the interaction (isotropic) between dust grains that are produced by the shielding effects of free ions and electrons in the plasma background. It is good approximation that makes it possible to describe many experimental results. The Yukawa potential is helpful for qualitative analytical and numerical estimations of structural and dynamical properties of CPs and self-organization of dust particles [7]. The Yukawa potential is commonly used as a pair potential for CPs and other physical systems such as chemical, polymers, physics, medicine and biology, astrophysics, and environmental sciences. Most scholars used the Yukawa pair potential for investigations of transport properties of CPs such as thermal conductivity [2][3][4][5], shear viscosity [38], and diffusion [27][28][29][30][31][32][33][34][35][36][37]. In the present study, we have applied a uniaxial ac electric field as additional interaction (anisotropic) with Yukawa pair interactions. This idea has been taken from PK-3 Plus laboratory under microgravity conditions on board of ISS and MD simulations [14,15], theoretical approach [13], and experiments during a parabolic flight on board of ISS [19]. For ERCPs, charged dust grains interact with each other through the Yukawa pair potential and quadrupole (anisotropic) interactions: the equation of an anisotropic interparticle potential is given as where the first term is the Yukawa potential l [2][3][4][5], λ D is Debye length, Q d is the charge on dust grains, r is the interparticle distance, and ϵ 0 is the permittivity of free space. The second term shows the interaction between the charge of one grain and the quadrupole part of the wake produced by another grain.
The charge-quadrupole interaction is the same as the two equal and parallel dipoles with a magnitude ≡ 0.65 M T Q d λ D : thus, the second term can be referred to as "dipole." The θ is the angle between the electric field (E) and vector r connecting two dust grains are shown in Fig. 1, and M T is the "thermal" Mach number. The square of "thermal" Mach number (M T ) can be written as T and is normalized with the thermal velocity of dust grains; here, v 2 T = T n ∕m d (equal to that of neutrals), where T n and m d are the temperature of neutral and mass of dust grains [14][15][16]. The external applied AC electric field (z-direction) caused the oscillation of dust grains with the velocity amplitude u d (drift velocity) and made the potential symmetric. Thus the "thermal" Mach number becomes proportional to the electric field (M T α E). At a large distance, the potential can be extended into a series over a Mach number with coefficients that are proportional to corresponding multipoles such as charge, dipole, and quadrupole. The interactions between charged dust grains at a low strength of the uniaxial electric field in ERCPs are similar to the dipolar interactions in conventional ER fluids; further details can be found in References. [14-16, 18, 19, 25]. It should be mentioned here that Eq. (1) has already been used in the CPs system for theoretical studies [18], experiments [14,15,19], and MD simulations [14,15,25]. Two important parameters are fully characterized CP systems. The first parameter is the Coulomb coupling where a is the Wigner-Seitz radius written as (3∕4 n) 1 3 , n is equilibrium dust density, k B is Boltzmann constant, and T is the system temperature. The second essential parameter is the Debye screening l e n g t h , = a∕ D . T h e p l a s m a d u s t f r e qu e n c y Sketch illustrating the interaction of dipoles that are induced through an external uniaxial AC electric field (E). For all moments in the direction of the field (z-direction), the interaction is computed by the angle ɵ and distance r adopted from [15] pd = √ Q 2 d ∕3 0 ma 3 describes the time scale of a CP system, where m is the dust grain's mass [3,6].

Einstein relation for diffusion coefficients
The Einstein formula has been used to calculate the diffusion coefficients of CP liquids by employing MD simulations. Using the Einstein method, diffusion motions are analyzed from the slope of computed mean squared displacement in a long time limit and can obtain diffusion coefficients [20-23, 30-32, 35]. Here, we split the EF into parallel and perpendicular parts of self-diffusion coefficients corresponding to field direction given by where r(t) represents the initial position, r(t 0 ) is the final position of j th particles, and ⟨… ⟩ is the average of the entire particle ensemble [22,[33][34][35]. The phase transition of fluids is connected with the diffusion of particles. The external electric and magnetic fields affect the motion of particles in parallel and perpendicular directions separately [22][23][24]40]. ERCPs change their phases as the external electric field varies [14][15][16]. We employed Eqs. (2) and (3) to compute parallel and perpendicular self-diffusion coefficients in ERCPs (anisotropic plasmas). It should be mentioned here that diffusion coefficients were already estimated using GK [24] and the Einstein relation [22,23] for an anisotropic plasma due to an external magnetic field in both parallel and perpendicular directions.

Execution step of MD simulations
The simulations start from face-centered cubic (FCC) lattice configurations of Yukawa identical particles having randomly initial velocities sampled from the Maxwellian distribution function concerning the required initial temperature value. N = 2048 particles are used for the whole calculation at different κ and Γ. Three-dimensional periodic boundary conditions for each direction are imposed on the simulation box to simulate the infinite system. The dimension of a cubic simulation box with an equal edge length is L∕a = (4 N∕3) 1∕3 [2,3,6]. The standard Ewald summation method was used to calculate the log-range interaction between dust particles. In the Ewald summation technique, the total interaction potential is divided into two parts, one is real space and the second is reciprocal space, and both are rapidly convergence. The Ewald convergence parameter value (γ = 5.6/L) is used. With higher κ values, the real space term alone shows enough accuracy and performance [43]. The interaction and forces are computed into two parts: the repulsive interaction through the typical Yukawa potential given in the first part of Eq. (1) and the anisotropic interaction using an electric field was given in the second part of Eq. (1) [17]. The equation of motion is integrated using the standard predictor-corrector integration method [5,33,34]. When an electric field is switched on, many physical phenomena occur; these novel regimes are not limited to fluid-solid transition but also exhibit chain, ring, sheet, vortex, etc., depending on the nature of the electric field [17]. For intermediate to high strength of M T , isotropic interactions between particles are ineffective; only anisotropic interactions dominate and cause the phase transition. Phase transition phenomena are connected with the diffusion of particles. The nature of diffusion of particles in parallel and perpendicular to the electric field direction is investigated [15][16][17].
Appropriate input parameters are chosen for the consistency and accuracy of the employed numerical scheme. These input parameters such as system temperature (T = 1/Γ), screening strength (κ), uniaxial electric field (M T ), number of particles (N), simulation time step, and total run time are selected to provide SC-ERCP diffusion coefficients. MD simulations are executed 2 × 10 5 /ω pd t time unit to record the simulation outcomes [6]. In this paper, the EMD method is used for SC-ERCP systems to estimate self-diffusion coefficients over a wide range of plasma coupling (Γ = 1-1000), screening (1 ≤ κ ≤ 3) and the number of particles (N = 2048).

Simulation results and discussion
In this section, the diffusion measurements obtained for SC-ERCP systems that we have computed for systematic EMD simulations using Eqs. (2) and (3), respectively, are presented at different values of Coulomb coupling (Γ), Debye screening strength (κ), and normalized electric field strength (M T ). The diffusion coefficient parallel D ║ (Γ, κ, M T = 0.01) and perpendicular D ┴ (Γ, κ, M T = 0.01) to the electric field are normalized by plasma frequency (ω p ) with applied electric field strength and are discussed and compared with the previously known theoretical and simulation data at nearly the same sets of plasma state points (κ, Γ) of SC Yukawa systems in the nonideal gas-like, liquid-like and solid/ crystal-like states. The diffusion coefficients are reported here with appropriate normalizations for a wide range of coupling (Γ ≥ 1), screening length (κ ≥ 1), and varying M T strengths (0.01 ≤ M T ≤ 10). The normalization for the diffusion coefficients has been used in previous studies on Yukawa (CPs) systems [35][36][37]. The D ║ and D ┴ for the SC-ERCP system may be normalized by the plasma frequency D 0 ║ = D ║ /ω p a 2 and D 0 ┴ = D ┴ /ω p a 2 . Interestingly, the current EF-EMD scheme provides good estimations of diffusion coefficients over an extensive range of plasma states (Γ, κ) than those earlier employed different numerical schemes [24,35] for SC Yukawa systems. It is generally experienced in the present EMD scheme that there are various restrictions that can be controlled in order to find well-organized results. These restrictions include the κ, the system temperature, M T , ω, dt, and total simulation length are changed to determine how these restrictions may influence the diffusion of SC-ERCP systems. It is stated here that notable results have been computed for the diffusion of ER systems employing EMD simulations at nearly the same sets of plasma points (Γ, κ) with appropriate system size N as those used in earlier simulation methods to demonstrate that the present simulation results of diffusion are well identified with earlier available simulation results.
The employed appropriate varying uniaxial AC electric field strength gives near equilibrium diffusion D ║ (Γ, κ, M T = 0.01) and D ┴ (Γ, κ, M T = 0.01) measurements, which are suitable for the whole range of plasma states of Γ ≡ (1, 1000) and κ ≡ (1, 3). The obtained data have matched well with the earlier simulation data of 3D CPs at electric field strength M T = 0.01. The wide spectrum of these electric field strengths is useful for calculating the steady-state values of diffusion D 0 ║ and D 0 ┴ [equilibrium values of diffusion coefficients], where the signal-to-noise ratio is appropriate for mentioned plasma state points (Γ, κ). In the present case, it is exercised that the varying electric field strength 0.01 ≤ M T ≤ 10 have diffusion coefficients within statistical error for the computed plasma parameters.
In this article, we have employed an already established relation for the anisotropic interparticle potential [14,15,19] account as ERCPs and EMD technique [16][17][18][19][20][21][22][23][24]. It should be mentioned that EF-EMD has already been used in earlier studies to compute parallel and perpendicular diffusion coefficients in a magnetized plasma [22,23].  Figure 2 shows our major results and compare former simulation trends from EMD measurements of SC Yukawa systems by Ohta and Hamaguchi [35] at a small N, and Daligault [36,37] at different screening and coupling parameters. It should be mentioned here that the data of diffusion coefficients of Ohta and Hamaguchi [35] are rescaled from Einstein frequency (D * = D/ω E a 2 ) to plasma frequency (D 0 = D/ω p a 2 ) for proper comparison with the present and earlier results. Moreover, the present EMD investigations of diffusion D ║ (Γ, κ, M T = 0.01) and D ┴ (Γ, κ, M T = 0.01) are also discussed with those estimated from different computational methods for Yukawa systems of Dzhumagulova et al. [27,32], Vaulina et al. [20,21,30,31], Ott and Bonitz [22], and Begum and Das [24]. It is examined that the present results are in practically excellent agreement with earlier numerical estimations for SC Yukawa systems (SCYSs), and it is shown that the current data obtained for ER systems and former results for SCYSs have a comparable efficiency, all generating close trends for diffusion D 0 ║ and D 0 ┴ coefficients. Furthermore, it is noted that the present EMD outcomes are found to be in a reasonable agreement with experimental results at a low Coulomb coupling range (Γ = 0.3, 4) for the mentioned screening strength (κ = 0.1, 0.58) [29] as well as the Enskog kinetic theory for SC-CPs [28].
The present scheme is tested first by calculating the D ║ (M T = 0.01) and D ┴ (M T = 0.01) for SC-ERCPs with constant uniaxial AC electric field strengths (≡ 0.01) to verify the accuracy and consistency of the employed numerical method. It is analyzed that the improved EMD scheme derived from the Einstein relation is more appropriate for diffusion coefficients. The current numerical scheme gives an alternative to the more usual method of computing the suitable diffusion coefficients of ERCPs. It is shown that the ER and Yukawa systems have a comparable performance and yield diffusion results through the GK relation and Einstein relation with the EMD technique, which is in satisfactory agreement with earlier known simulation results.  (3) different simulations corresponding to the mentioned Γ (a total of 68 simulation datasets for the three cases), covering from the non-ideal phase (Γ = 1) to the strongly coupled phase (Γ = 1000) depending on κ. These panels show the main results and are compared with earlier numerical results obtained from the EMD simulations of SCYSs by Ohta and Hamaguchi (rescaled by plasma frequency) [35] and Daligault [36,37] of classical one-component SCYSs for almost the same and close state points of (Γ, κ). Moreover, the present EF-EMD simulations of D 0 ║ (Γ, κ,M T = 0.01) and D 0 ┴ ( Γ, κ,M T = 0.01) have fair agreements with the computed 3D EMD estimations of Begum [24], Dzhumagulova et al. [27], and Vaulina et al. [20,21,30,31]. It is found that the present results of diffusion coefficients for the CP system are in generally excellent agreement with earlier estimations of SCYSs, and it is shown that the SC-CPs system and SCYSs have a comparable efficiency, all producing the same trends for D 0 ║ ( Γ, κ,M T = 0.01) and D 0 ┴ ( Γ, κ,M T = 0.01). The panels in Fig. 2 show that the present results of diffusion coefficients are matched well to the former MD results calculated from the Einstein relation (D E ) and GK relation (D G ) by Ohta and Hamaguchi [35] and D G by Daligault [37]. and c κ = 3 parameters. The present results are compared with earlier EMD investigations reprinted by Ohta and Hamaguchi [35] and Daligault [36,37] for strongly coupled Yukawa systems We have also calculated the plasma D ║ (M T = 0.01) and D ┴ ( M T = 0.01) for temperature scaling law at normalized constant M T = 0.01. The temperature scaling expression shows normalized plasma diffusion coefficients by Einstein frequency (D * ║ = D ║ /ω E a 2 and D * ┴ = D ┴ /ω E a 2 ) as a function of normalized scaled melting temperature written as T * (≡ T/ T m = Γ m /Γ), where Γ m is a liquid-solid critical point depending on plasma Debye screening length (κ). Further details associated with melting temperature for the 3D Yukawa system and scaling law for diffusion coefficients in earlier literature [35] are given as where A, B, and C are the unknown fitting parameters that can be obtained using Eq.    [25] and SC-CPs in a homogeneous electric field [26] and MD simulation data for ionic liquids [40]. Now consider the cross-field self-diffusion coefficients  ┴ decreased with increase in the M T strength (higher than the mentioned regime). It should be said here that the D 0 ┴ diminished at higher Γ (≡ 300, 500) and M T (> 3). However, it is noted that the diffusion coefficients decrease with increase in Γ, and it slightly increases with increase in κ. The behavior of D 0 ┴ (M T , Γ, κ) in electrorheological is consistent with magnetized plasma [22,23]. These results are within the statistical uncertainty range, confirming earlier conventional ER fluids and ionic liquids [40].
The evaluation of calculated parallel (D 0 || ) and perpendicular (D 0 ┴ ) diffusion coefficients illustrates that the difference between the various sets of simulation results is reported at a wide range of plasma Coulomb coupling. The EMD data obtained in the current work clearly describes the effects of uniaxial AC electric field at diffusion coefficients of SC-ERCPs. First, the EMD method needs simulations for various values of M T to establish the linear regime. For possible reasons for increasing behaviors of parallel diffusion coefficients with increase in M T , one is that the electric field produces the plasma flow. The external M T on plasma affects the interactions of dust particles that it causes to "tune" from isotropic to anisotropic interactions between dust grains. Furthermore, CPs apply an electric field that polarizes grains, inducing additional dipole-dipole coupling. The drift velocity of the dust particle also increased with increase in M T [14-16, 18, 19]. The effects of M T can change interactions between dust particles and cause an additional force. Second, in the usual EMD scheme considering the electric field, the breakdown of the ion cage structure needs lower activation energy for the free region [14]. The increasing trend of diffusion coefficients indicates that the system may go from subdiffusion to superdiffusion and disordered to ordered states in the presence of the AC mode of an electric field. The current simulation outcomes of D 0 || are found to be consistent with CPs [26], ERCPs [25], and D 0 ┴ with magnetized plasmas [22,23].

Conclusions
In conclusion, we have reported equilibrium molecular dynamics (EMD) simulation results based on the Einstein formula (EF) to illuminate the diffusive characteristics of 3D SC-ERCPs under the influence of an external uniaxial (z-axis) AC electric field (M T ). The parallel (D ║ ) and perpendicular (D ┴ ) diffusion coefficients to M T were investigated for a wide range of coupling (Γ) and screening length (κ) parameters. The results obtained through the EF-EMD simulations for SC-ERCPs and those measured through EF and GK-EMD investigations for Yukawa systems are found to be in reasonable agreement, generally overpredicting within statistical limits. The accuracy of the modified EF-EMD scheme for ER systems is excellent in its wide range of fields 0.01 ≤ M T ≤ 10. It is accomplished that the present scheme for measuring D ║ and D ┴ through the EF-EMD method yields reliable results, and this method is quite accurate for ERCPs. Three regimes of electric-field dependence of D ║ and D ┴ have been identified at strong coupling. D ║ is independent on the low strength of M T , has an increasing behavior for intermediate M T strength, and again constant for large M T strength. D ┴ also is independent on M T for small strength, faster decay at large, and low Γ (= 5, 10). We have also observed an intermediate decay of D ┴ at intermediate to higher coupling parameters (≡ 20, 500) and M T (≡ 1, 10). These are the unique features of SC-ERCPs induced by an external electric field. Our investigations are expected to be directly relevant to electrorheological materials in a compact of smart, precise tailoring materials. Moreover, this numerical scheme can be used to explore other transport properties and phase transition phenomena of ERCPs.