**2. MD simulation algorithm and parameters**

The Computer simulation provides a linkage between theoretical and experimental research work. In thermal fluids, science and engineering, the installations of experimental setups are very complex and high cost. In the age of modern technology, the cost-effective and low time consumption high computational power is most prominent. Nowadays, it is a trend that before starting the experimentations, first test with computational tools than verified with experiments. There are various computer algorithms and techniques are used for the calculation of various properties in various materials. We perform computer simulations to test a theoretical model, verify experimental data, and also for comparison purposes [38]. Different computational techniques can be designed for extreme conditions, such as for very low temperature and high density. It also acts as a bridge between microscopic length, time scale, and macroscopic worlds. MD and Monte Carlo (MC), and Langevin dynamics simulations are the main methods used to compute CDPs' physical properties. Different software's are also available to compute various physical properties and build a new model by following one simulation scheme.

MD simulations have become a prominent computational tool to investigate various properties such as thermophysical and dynamical properties in different types of fluids and materials. The MD simulations consist of the numerical solution of the Newton equation of motion for the spherical particles system [9]. This section explained the MD simulation algorithm to calculate the self-diffusion

coefficient (*D*). Equilibrium MD simulations are used to investigate the *D* in ER-CDPs. Yukawa potential (screened Coulomb) is the most commonly used potential for CDPs systems, while other physical systems such as physics of chemical and polymer, medicine and biology systems, astrophysics, and environmental research are adopted. Most scholars used Yukawa potential to screen charged dust particles in the dusty plasma results in isotropic interactions. In the present research used applied a uniaxial ac electric field for additional dipole–dipole-like interactions. This idea was taken from Pk-3 plus experiment under the microgravity conditions [20], theoretical approach [39], MD simulation study [24], and *Pk*-4 experiments [23]. For ER-CDPs, the charged dust particles interact with each other through Yukawa potential and Quadrupole interactions due to external applied ac electric field: the equation for ER-CDPs becomes as

$$W(r,\theta) = Q\_d^2 \left[ \frac{1}{4\pi\varepsilon\_o} \frac{\exp^{-r/\lambda\_D}}{r} - 0.43 \frac{M\_T^2 \lambda\_D^2}{r^3} \left( 3\cos^2\theta - 1 \right) \right] \tag{1}$$

First-term in the above equation presents pair-wise Yukawa potential in the absence of an electric field, and the second term shows the Quadrupole interaction between particles due to uniaxial ac external electric field [22]. Where *λ<sup>D</sup>* is Debye length, *Qd* is the charge on dust particles, *r* distance between interacting particles, and *ϵ<sup>o</sup>* permittivity of space. The *ɵ* is the angle between the electric field (*E*) and interacting dust particles. In the present study, we take *ɵ* = 0o for uniaxial ac electric field for anisotropic interactions. The *r* is the distance between the interacting particles, *r* = *r*<sup>j</sup> - *r*i. The *MT* is thermal Mach number normalized with the thermal velocity of charged dust particles *M*<sup>2</sup> *<sup>T</sup>* <sup>¼</sup> *<sup>u</sup>*<sup>2</sup> *d* � �*=v*<sup>2</sup> *<sup>T</sup>* where *v*<sup>2</sup> *<sup>T</sup>* ¼ *Td=md* the drift velocity is proportional to an electric field (*ud* α *E*), an electric field can be measured in the unit of *MT.* For small values of *MT,* in ER-CDPs, the interactions of particles are the same as dipolar interaction in conventional ERFs. In the prior study, strings (chain) of dust particles were observed in typical conditions when an external ac electric field was applied.

There are two central (dimensionless) parameters, which are fully-characterized complex plasma systems. The first parameter is the plasma coupling parameter (define same as Coulomb systems) *<sup>Γ</sup>* <sup>¼</sup> *<sup>Q</sup>*<sup>2</sup> *<sup>d</sup>=*4*πϵoaw:skBT*, *a*ws Wigner-Seitz radius defined as 3ð Þ *<sup>=</sup>*4*π<sup>n</sup>* <sup>1</sup> <sup>3</sup>πn �1 <sup>2</sup> Where *n* is the equilibrium dust number density, *kB* the Boltzmann constant, and *T* is the system's absolute temperature. The plasma coupling (*Γ*) parameter is the ratio of average potential energy to the average kinetic energy of interacting dust particles. The SCCDPs associated with *Γ* > 1 inversely account as weakly coupled dusty plasma (*Γ* < 1). Another essential parameter is the Debye screening parameter, *κ* ¼ *aw:<sup>s</sup>=λD*. The plasma dust frequency *ωpd* ¼ *Q*2 *d=*3*πϵ*0*ma*<sup>3</sup> � �<sup>2</sup> describes the time scale of a complex plasma system; here, *m* is the dust particle's mass.

### **2.1 Green-Kubo relation of diffusion coefficient**

Green-Kubo integral formula were used to calculate the self-diffusion coefficient (*D*) for 3D CDPs [30, 33], one component plasma [31] and ionic mixtures [32]. Here we have used the same numerical models for the investigations of *D* in 3D ER-CDPs; the Green-Kubo relation is given as following.

$$D = \frac{1}{3N} \int\_0^\infty Z(t)dt\tag{2}$$

*Studies of Self Diffusion Coefficient in Electrorheological Complex Plasmas… DOI: http://dx.doi.org/10.5772/intechopen.98854*

In Eq. (2), *Z*(*t*) is known as the velocity autocorrelation function (VACF), 3 shows the 3D system, and *N* represents the number of simulated dust particles. The *D* is the integral of VACF calculated over all the particle velocity product's ensemble average segments at a time (*t*) and initial time (*t*0). If Z(*t*) decays too slowly, Yukawa particles' motion is described as anomalous diffusion for the integral equation one converges. The transport coefficient's existence, autocorrelation function must rapidly decay enough for integral to convergence. Different types of diffusion motion were analyzed through VACF [12], defined as

$$Z(\mathbf{t}) = \left< v\_j(\mathbf{t}).v\_j(\mathbf{0}) \right> \tag{3}$$

Where *v*<sup>j</sup> (*t*) denote the *j th* particle's velocity at simulation time (*t*) and (0), the brackets h i *:* … represent the canonical ensemble average all over particles. *Z*(*t*) demonstrates the decay that is characterized by plasma frequency. The existence of the transport coefficient, autocorrelation function must rapidly decay enough for integral convergence [29, 33, 34].

In this method, classical molecular dynamics simulations are used to map the trajectories of *N* = 500 and calculated the diffusion coefficient. Periodic boundary conditions (PBCs) were used to be constrained enforced during the *N*-particles simulation under statistical microcanonical ensemble (*NVE*). The edge length of the 3D cubic simulation box is *L* with 12.79*a*. Particles are placed randomly in the simulation cubic cell for the beginning of the simulations [15]. An appropriate time step (*dt* = 0.001*ωpd*) was selected to calculate particle position and velocity at each step. The leapfrog integration method is used to integrate the equation of motion and calculate each particle's force, acceleration, and position for every time step. The whole scheme in MD dynamically simulates the equation of motion of *N*-dust particles with interaction through Yukawa potential, and Quadrupole interaction was given in Eq. (1). The Ewald summation method was used for calculations of Yukawa and dipolar potential energies and forces. Here we select the Ewald convergence parameter (*γ* = 5.6/*L*) for the Yukawa and anisotropic interactions. Suitable input parameters are selected to get the accuracy and consistency of the model. The appropriate system temperature (ð Þ *T* ¼ 1*=Γ* Yukawa Ewald-summation method, screening length, Mach number, time step, number of particles, simulation run time, etc., are selected to provide 3D ER-CDPs diffusion motion. Negative derivatives of Eq. (1) calculate forces on each interacting charged dusty particle in the Yukawa system. MD simulations are performed 2 x 10<sup>5</sup> /*ωpd* time unit to record the SDC (*Γ*, *κ*, *MT*). The 12-core-processor takes about three hours to perform each simulation. Total 156 simulations are performed and calculated the VACF and *D* for a wide range of plasma parameters and external applied uniaxial electric field. For the present chapter, the EMD method is used of ER-CDPs and calculated the *D* over a wide range of *Γ* = 2–500, 1 < *κ* ≥ 3, and *N* = 500. In addition, the conservation of total energies and momentum is checked and verify that system has selfconsistency.

## **3. Simulation results and discussions**

This book chapter used the equilibrium MD simulation method using Yukawa (screened Coulomb) potential and dipolar interactions. Without external forces (electric and magnetic), Yukawa potential was used and calculated thermophysical properties. In this work, we add dipolar interactions with the same numerical schemes and analyses the diffusion motion through VACF and *D.* First of all; we compute the VACF for Screening parameters (1 < *κ* ≥ 3), plasma coupling

parameters (1 < *Γ* ≥ 500), uniaxial electric field strength (0.0 ≤ *MT* ≥ 1.50) for the same number of particles (*N* = 500). After that, *D* is calculated for nearly the same parameters used for the VACF.

The single-particle motions are computed for 3D ER-CDPs through the equilibrium MD simulation using VACF for diffusion by employing Eq. (3). VACF are analyzed and discussed for plasma screening *κ* = 1.4, 2 and 3, plasma coupling (*Γ* = 2– 500), number of spherical charged dust particles (*N* = 500) with the variation of the external uniaxial electric field. The VACF has been widely used for 3D CDPs over a wide combination range of plasma parameters (*κ*, *Γ*) [12, 30, 36]. The VACF was computed at *MT* = 0.0 for comparison purposes with earlier MD results. It has been shown that MD outcomes have comparable fair agreements with previous results. Here we defined the critical strength of the applied external electric field (*MCT*). The *MCT* is directly proportional to plasma coupling strengths or inversely proportional to system temperature. We have been found that above the *MCT* system goes out of thermal equilibration. We have also got noisy results above the *MCT* strength of the external applied uniaxial electric field. The critical strength of the uniaxial electric field is different for the different combinations of CDPs parameters.

The simulation outcomes of VACF in 3D ER-CDPs as a function of the simulation time scale for four different plasma coupling values (*Γ* 20, 50, 100, and 250) at constant *κ* = 1.40 and *N* = 500 are displayed in **Figure 1**. The effects of variation uniaxial electric field (*ΜΤ* 0.01–1.30) on VACF are computed to analyze the diffusion motion of charged dust particles. It has been shown that *MT* significantly depends on plasma temperature and the screening of charged particles. At Higher plasma temperature (*Γ* = 20), one can easily observe the rapid decay of VACF with very weak particle oscillations. This regime *ΜΤ* does not have significant effects on the dynamics of dust particles. **Figure 1(b)** shows the oscillatory damped motion and slightly decreasing amplitude of VACF up to uniaxial electric field strength (*ΜΤ* = 0.90). For higher plasma coupling strength (*Γ* 100,250), results of VACF show the higher oscillations with slightly damping phenomena. **Figure 2** demonstrated the VACF at *κ* = 2.0 for four different regimes of strongly coupled CDPs (*Γ* = 20, 50, 200, 400) *N* = 500 with a variation of uniaxial electric field (0.0 < *MT* ≥ 1.50). It is observed that the effect of *MT* on VACF is nearly the same as was observed in **Figure 1**. It has also been shown that the *MCT* increases with increases in the screening strength. Furthermore, we increased the screening strength *κ* = 3.0 and *Γ* = 30, 75, 250, and 500 are shown in **Figure 6**. The trend of damping and oscillations of dust particles was observed the same as in **Figures 1** and **2**. Therefore, the VACF below the *MCT* simulation outcomes is founded in good agreement with previous data [13, 30, 35, 36].

It is observed that at higher temperatures (low coupling, *Γ*), the *MT* does not affect the dynamics and oscillation of the particles. So what we can say that at higher system temperature, the radius between interacting charged dust particles is large. Due to the large distance, Eq. (1) remains invalid for higher uniaxial electric field strength. The motion of a particle in this regime indicates only thermal motion. Complete damped oscillations of dust particles were observed for a long time at intermediate values of *Γ*. The magnitude of oscillations slightly decreased with increasing *MT*, increasing the diffusion at a low rate of thermal motion. In ER-CDPs, the effect of *MT* on VACF can be easily observed from intermediate to higher plasma coupling (*Γ*) strength. The simulation time length scales, maximum at intermediate values of *Γ* and further increasing *Γ*, the length scale slightly decreased. It was observed that the oscillation was not fully damped at higher *Γ*. It is concluded that a uniaxial electric field can affect diffusion motion in ER-CDPs for the sufficient strength of *MT*., so what can say the presented numerical method gives us precise and reliable data for comparative study.

*Studies of Self Diffusion Coefficient in Electrorheological Complex Plasmas… DOI: http://dx.doi.org/10.5772/intechopen.98854*

### **Figure 1.**

*VACF as a function of molecular dynamics simulation time scale for constant plasma screening parameter* (κ *= 1.40). MD simulations results obtained for* N *= 500 simulated dust particles with the variation of uniaxial external electric field (0.0 <* MT ≤ *1.30) and four different plasma regimes (a)* Γ *= 20, (b)* Γ *= 50, (c)* Γ *= 100 and (d)* Γ *= 250.*

Now we focused on the primary outcomes through the equilibrium MD method. The simulation data of this work are that the self-diffusion coefficient (*D*) of ER-CDPs can be obtained with minimum statistical error at *MT* = 0.0 by equilibrium MD. Moreover, it is shown that the MD with GKR has a fair good agreement with previous simulation outcomes of Ohta and Hamaguchi [30], Daligault [31, 32], and Begum and Das [33, 37] at *MT* = 0.0. **Figures 3**–**5** display the primary outcomes of *D* measured from the equilibrium MD method for ER-CDPs at *κ* = 1.4, 2, and 3, respectively. *D* of ER-CDPs was calculated for different uniaxial electric field strengths (0 < *MT* < 1.5). Here we have explained the *D* for a wide range of combinations of plasma parameters (*Γ*, *κ*) over acceptable *MT* ranges below the *MCT*. This section focused on the variation of plasma coupling (1/*T*), electric field, Debye screening length, and constant *N* = 500 simulated particles.

The effect of *MT* on *D* in 3D ER-CDPs is calculated for plasma coupling (*Γ* 20, 50, 100, and 200), Debye screening length (*κ* 1.4), uniaxial electric field (0.0 ≤ *MT* ≤ 1.30), and *N* = 500 number of particles shown in **Figure 3**. We have performed several MD simulations at constant *κ* = 1.4 and *N* = 500 to analyze the effect of *MT* on *D*. The reported data of *D* shows good agreements for small *MT* values with previous results [30–32]. With respective of *Γ* different regimes are under consideration. It is found that the *MCT* values are the same as in **Figures 1, 2, 6**

### **Figure 2.**

*VACF as a function of molecular dynamics simulation time scale for constant plasma screening parameter (*κ *= 2.0). MD simulations results obtained for* N *= 500 simulated dust particles with the variation of uniaxial external electric field (0.0 <* MT ≤ *1.30) and four different plasma regimes (a)* Γ *= 20, (b)* Γ *= 50, (c)* Γ *= 200 and (d)* Γ *= 400.*

for the same combinations of plasma parameters. It is noted that for plasma coupling (*Γ* = 20), the critical values (*MCT >* 0.42) for *κ* = 1.4 display in **Figure 3**. We have found the increasing trend of *MT* with increases the plasma coupling values or decreasing the system temperature such as *Γ* = 50, 100 and 200 critical *MCT* ≥ 0.88, 1.12, 1.25 respectively. The possible reason for the low strength of *MT* in high system temperature is the large interparticle distance show with the term (1/**r** 3 ) in Eq. (1). The increasing values of *D* with respective of *Γ* at *κ* = 1.4 can be easily observed from four panels of **Figure 3**. Upon further increasing *MT*, the MD simulation gives an error in the form of out of thermal equilibrium. At lower plasma coupling strength (*Γ* = 20), the *D* does not significantly increase under the applied external electric field's influence; only 0.03% observed the integral values of VACF. The increasing of diffusion was observed for *Γ* = 50 are 20%, *Γ* = 100 are 75% and *Γ* = 200 are 160%. The electric field effect on *D* in ER-CDPs can be found in the same as conventional electrorheological fluids and ionic liquids [40].

The simulation outcomes of *D* in 3D ER-CDPs at constant *N* = 500, different values of *Γ* as a function of the uniaxial electric field (*MT*) for *κ* = 2 and 3 are shown in **Figures 4** and **5**, respectively. The effect of *MT* was computed in four different

*Studies of Self Diffusion Coefficient in Electrorheological Complex Plasmas… DOI: http://dx.doi.org/10.5772/intechopen.98854*

### **Figure 3.**

*Self-diffusion coefficient as a function of external applied uniaxial electric field (0.0 <* MT *< 1.30), for constant plasma screening (*κ *= 1.4) and a number of particles (*N *= 500), considered four different ER-CDPs states (a)* Γ *= 20, (b)* Γ *= 50, (c)* Γ *= 100 and (d)* Γ *= 200.*

### **Figure 4.**

*Self-diffusion coefficient as a function of external applied uniaxial electric field (0.0 <* MT *< 1.40), for constant plasma screening (*κ *= 2) and a number of particles (*N *= 500), considered four different ER-CDPs states (a)* Γ *= 20, (b)* Γ *= 50, (c)* Γ *= 100 and (d)* Γ *= 200.*

### **Figure 5.**

*Self-diffusion coefficient as a function of external applied uniaxial electric field (0.0 <* MT *< 1.50), for constant plasma screening (*κ *= 3) and a number of particles (*N *= 500), considered four different ER-CDPs states (a)* Γ *= 30, (b)* Γ *= 75, (c)* Γ *= 250 and (d)* Γ *= 500.*

states of ER-CDPs by performed almost thirty-three simulations. The *MCT* values for *κ* = 2.0 at *Γ* = 20, *MT* > 0.45, *Γ* = 50, MCT > 0.90, *Γ* = 100, *MCT* > 1.23, and *Γ* = 200, *MCT* > 1.35. The increase in *D* under the action of external applied uniaxial electric field for *Γ* = 20, 50, 100, and 200 are 1.7, 11, 42, and 115%. The increases in the *D* for *κ* = 3.0 same *N* were observed nearly the same as said above for *κ* = 2.0. The effect of *MT* on the *D* is prominent at higher values of *Γ* parameters. It was also noted that *MT* does not significantly affect *D* lower to intermediate values of *Γ*. From **Figures 3**–**5**, we can conclude that the proposed MD simulation method worked for the limited strength of *MT*. It was noted that the system's critical values above the system go out of thermal equilibrations, and kinetic energy increased very largely up to 23 orders of magnitude, but the potential energy does not change higher order of magnitude. Below the critical strength, the potential energy of interacting Yukawa dust particles increases with *MT*.
