Non-Newtonian Dynamics with Heat Transport in Complex Systems

Abstract

Transport properties of complex system under various conditions are of practical interest in the field of science and technology. Homogenous nonequilibrium molecular dynamics (HNEMD) simulations have been employed to calculate the thermal conductivity (λ) of three-dimensional (3D) strongly coupled complex nonideal plasmas (SCCNPs) over a suitable range of plasma parameters (Γ, κ). New investigations show that the λ depending on plasma parameters and minimum value of λ exists at nearly same plasma states. In the present case, the non-Newtonian behavior is checked with different system sizes and it is found that the λ behavior is well matched with earlier numerical work. It is demonstrated that the present outcomes are more consistent than those obtained earlier known simulations. It is revealed that our outcomes can be acceptable for a low range of force field in order to find out the size of linear ranges, and it explains the nature of nonlinearity of SCCNPs. It has been shown that the measured outcomes at steady states of external field of F* (=0.005) are in acceptable agreement with previous numerical outcomes, and it showed that the deviations are within less than 12% for most of the data and depend on plasma states.

1. Introduction

The computational knowledge of thermophysical properties is very different of complex liquids as compared to the nonionic liquids. The important thermal conductivity of complex liquids is used in the heat design process as an important parameter. The estimations of thermal conductivity obtained by applying the molecular dynamics (MD) approach in liquids and crystalline solids are a difficult job due to perceived limitations of computational power [1]. Recently, in the field of science and technology, the transport properties of interacting particles in complex nonideal systems are of practical importance. A deep understanding of the interaction of complex systems is required for nano- and microstructuring of surfaces. The micron-size particles have recently been investigated in complex (dusty) plasmas, and in the physics and chemistry of plasmas, space environment, ionized gases, and material research and in the nuclear energy generation. The complex (dusty) plasmas play very important role in various technological applications, such as industrial processing of microelectronic devices, storage devices, and fuel burning, and future energy production [2]. For the explanation and understanding of these macroscopic phenomena, a comprehensive microscopic knowledge and calculation of the transport properties of complex (dusty) plasmas are required over the extensive range of plasma parameters (Г, к). Both the Coulomb coupling (Г) and Debye screening strength (к) are the dimensionless parameters, which can be used to characterize the plasma. In statistical mechanics, the microscopic dynamical origin of heat transport is a fundamental problem. Moreover, the purpose of the present work is to investigate the thermal conductivity dependences on the strength of different perturbation fields and to understand the non-Newtonian behaviors in the Yukawa liquids along with the calculations of thermal conductivity.

1.1 Dusty plasmas

Dusty plasmas are also known as nonideal complex plasmas that contain particulates of condensed matter. The dust particles may have sizes ranging from nanometers to micrometers, and typically much more massive than that of plasma ions, electrons, and neutrons. When the dust particles immersed in the plasma, they attain a high electric charge (negative charge) which makes the dusty plasmas interesting and technological important in the area of applied plasma physics. The dynamical behavior of these massive dust charge particles is much complex and occurs on considerably slower time scales, because their charge-to-mass ratio is in orders of magnitude smaller than that of the corresponding charge-to-mass ratio of either the electrons or ions. Dust particles are found in the large abundant in planetary plasmas, cosmic plasmas, plasmas in the laboratory, and plasmas near the earth. As a matter of fact, one may cogitate that except in the hottest regions of fusion plasmas, where the dust particles would not survive, most are known as dusty plasmas in the sense that some dust particles might be present. To understand this fact, one recognizes two cases in which: (1) there are just few secluded (noninteracting) dust particles, with the goal that they have nearly nothing if any impact on the plasma, and (2) there are countless number of dust particles in the plasma so that their existence really changes the properties and behavior of the plasma. In the event (1), the dust particles are charged by their interactions with the plasma yet do not change the plasma in any noticeable way. Then, again case (2) agrees to the situation in which the charge dust is a component of the plasma, subject to the collective interactions that recognize an ionized gas from a neutral gas. Case (2) is what is ordinary characterized as “dusty plasmas.” At a significantly bigger scale, it is outstanding that comets for the most part have two tails. One tail is expected to the comet’s dust particles, the other is because of ionized gas comet coma. These tails are not separate near the coma but overlap forming dusty plasma [3].

1.2 Dusty plasma parameters

There are two basic dimensionless parameters which are used for the analysis of transport coefficients in 2D and 3D dusty plasma systems, and which are responsible for mass transfer and phase state in nonideal dissipative systems [4, 5, 6]. These dimensionless parameters are known as effective Coulomb coupling (Г^*) parameter and scaling parameter (ξ). These parameters (Г^*, ξ) are responsible for transport and structural processes for nonideal systems [5]. Screening parameter (к) is the third one which is the important for the classification of dusty plasmas. Here, we discuss only two parameters, Coulomb coupling (Г) and the Debye screening strength (к).

1.2.1 Coulomb coupling parameter (Г)

The Coulomb coupling (Г) parameter is the ratio between the interparticle potential energy (P.E) to kinetic energy (K.E), and mathematically, it is written as, $\mathrm{Г}=\frac{P.E}{K.E}$ . This Coulomb coupling (Г) parameter is used for the classification of strongly coupled complex (dusty) plasmas and weakly coupled complex (dusty) plasmas.

1.2.2 Screening parameter (к)

Another important parameter of dusty plasmas is the screening strength (к), which is the ratio of interparticle distance to the Debye length, and mathematically, it is written as $к=\frac{a}{{\lambda }_D}$ , here “a” is the Wigner-Seitz (WS) radius and “λ_D” is the Debye length. The screening parameter (к) is also used for the classification of dusty plasmas.

1.3 Types of dusty plasmas

On the basis of Coulomb coupling (Г) parameter, the complex (dusty) plasmas are classified into two classes: one is called weakly coupled (ideal) complex (dusty) plasmas and the other is called strongly coupled (nonideal) complex (dusty) plasmas. These are succinctly discussed below.

1.3.1 Weakly coupled complex (dusty) plasmas

For weakly coupled complex (dusty) plasmas (WCCDP_S), the Coulomb coupling parameter (Г) is less than the unity (Г < 1), and it also called the weakly coupled ideal plasmas because Columbic collisions are negligible. Weakly coupled plasmas, like a gas, have no structure because Coulombic interactions are negligible between the particles and particle motion is like molecular motion in gases and particles have nearly random positions with respect to nearest neighbors [7]. For WCCDP_S, the K.E (thermal energy) is much larger than the Coulomb interparticle P.E of nearest particles.

1.3.2 Strongly coupled complex (dusty) plasmas

When the Coulomb coupling parameter (Г) is greater than the unity (Г > 1), then the plasma is known as strongly coupled (nonideal) complex plasmas. For weak-to-intermediate Coulomb coupling (Г) values, the SCCDP_S can have structure of liquids, and structure of solids for higher values of Г. Furthermore, SCCDP_S are the collection of free microparticles that interact with each other with a strong Coulomb repulsion force and have structure at microscopic scale for the arrangement of particles. The particle motion in SCCDP_S resembles that in the liquids or solids, and particles remain in relative fixed positions with respect to neighboring particles because of the strong Coulomb interactions present between the charged particles [7]. In SCCDP_S, the Coulomb interaction P.E is much larger than the K.E of nearby particles.

1.4 Formation and growth of dust particles in a plasma

An innovative feature of plasmas is that comprise chemically energetic species which grow the dust particles. It is true for plasma appliances that are using in the plasma processing semiconductor engineering, in which combinations of gases such as oxygen (O₂), argon (Ar), and silane (SiH₄) are castoff in the assembly and figments of microelectronic chips. Dust particles could also produce through sputtering, arcing, electron bombardment, etc., where the atoms or molecules from the walls of chamber or electrodes come out and immersed in the plasmatic system through a different mechanism, called plasma-material interaction [8]. Plasma processing devices are employed for the production of silicon wafers characteristically used as parallel plate electrode, in which 13.56 MHz radio frequency (rf) power is connected to the lower electrode to create the plasma. Etching process involves a reactive species such as silane along with a buffer gas like Ar. Plasma-aided gas phase chemical reactions produced silicon hydride (SiH₂) by the reaction {e⁻ + SiH₄ → (SiH₄)* → SiH₂ + 2H}. The vibrationally excited state is produced by the collisions of SiH₄ with electrons which then dissociates into SiH₂ [9]. The particle that grows in plasmas passes through certain phases like nucleation, coagulation, and surface growth.

1.5 Dust particle in the plasma

Dust particles are found everywhere in the entire universe with different shapes and sizes and mostly found in the solid form but also found in the liquid and gaseous ionized form. When the dust particles coexist with the plasma, then “dusty (complex) plasma” is formed. Dust particles acquire an electric negative charge, when these dust particles are immersed in the plasma and then affected by electric and magnetic fields and plasma properties are changed. Moreover, these dust particles attain electric negative charge (typically depends on the flow of ions and electrons) very fast due to the interactions between the dust particles and the nearby plasmas.

In recent years, dusty plasmas have opened up an entirely new field of research of science and technology by investigating transport properties (thermal conductivity, shear viscosity, and diffusion coefficient) of dusty plasmas in the laboratory.

1.5.1 Charge on dust particle

A lot of mechanisms are adopted to produce charge on a dust particle. If all mechanisms are considered at once, then the measurement of the equilibrium charge condition on a grain becomes very difficult. For the electron temperature T_e and ions temperature T_i, the flux of ions and electrons has an individual thermal velocity v_te. The thermal velocity of electron is higher than that of the heavier ions which have a minute thermal velocity v_ti. These motions develop the charge Q on the grain and make its surface potential (ϕ_s) negative. The charge on a grain fluctuates continuously due to the collective current of electron and ion, i.e., dQ/dt = I_e + I_i and equilibrium of charge occurs under condition I_e + I_i = 0 or ϕ_s = −2.49kT/e for an electron-ion plasma. The charge itself is associated to the surface potential by Q = Cϕ_s, in which C tells about the capacitance of a grain in a plasma [10]. The dust particle density differs from the density of electrons and ions because in normal plasma, neutral n₀ exists. If the primary electrons are very energetic, then they release subordinate electrons, which make the positive potential surface. The absorption of plasma ions also makes the positive potential surface. The transition of primary electrons and ions influenced by the surface potential of the grain depends on the relative velocity between the plasma species and the grain. Electrons are repelled and the grain current will be decreased if surface potential is negative. Electrons show attraction and the grain current increases for positive potential surface.

1.5.2 Size of dust particle

In dusty plasma, the dust particles can have any shape and can be made of either dielectric or conducting materials. The size of dust particles is much larger than the size of electrons and ions of plasmas which is in microns or tens of nanometers. So, the dust particles can easily be seen without any microscope. The typical size range of dust particles is from 100 nm up to about 100 μm. For experimental studies, dust particles that are distributed into plasmas are generally plastic or glass particles (commonly used particle is melamine formaldehyde). They are spherical in shape with a very narrow distribution of diameters. For example, the diameter of a normally used particle may be 3.50 ± 0.05 μm and a mass $\sim 3\times {10}^{-11}$  kg. Such particles are named monodisperse. Fine powders, such as aluminum silicate (kaolin) with a broad size distribution ranging from the submicron to tens of microns, are used in some experiments. Such powders contain particles having laminar shapes with jagged edges.

1.6 Forces acting on dust particles in a plasma

When the dust particles are immersed in the plasma, then various forces are acting on the dust particles which are significant because of their dynamics and transport characteristics. These various forces determine that where the particles are trapped or not, and these are sensitive to the position of dust particles within the plasmas.

1.6.1 Force of gravity (F_g)

The dust particles which are under a gravitational force are proportional to the mass of dust particle, and if dust particles are under microgravity condition, then it must be ignored.

$$F_g={\mathit{mg}}_d=\frac{4}{3}\pi a^3{\rho }_{d}g$$ E1

where “g” describes the gravitational acceleration and “ρ_d” represents the mass density of dust particles. Mostly, its value for solid materials is ρ_d = 2000 kg/m³. This force can be neglected for submicron particles, but for the micron sized or larger particles, this can be considered as the dominant force that typically confines the time during which the particle resides in the plasmas.

1.6.2 The electric force (F_e)

For plasma having an electric field (E), the electric force acting on dust species is written as:

$$F_e=\mathit{QE}$$ E2

This force is much larger in the bulk of the plasmas, while it is much smaller in the sheaths next to the plasma edge. For example, for characteristic radio frequency (rf) parallel plate glow discharge plasmas, particles reside beneath the negative electrode where the downward electric field offers an upward electric force that stabilizes the weight of the particle with force.

1.6.3 Neutral drag force (F_n)

This force is generated due to impacts of dust particles with the neutral gas species (atoms and molecules) and it is proportional to the neutral pressure in the vacuum chamber. Mathematically, it is written as:

$$F_n={\mathit{Nm}}_n{v^2}_{\mathit{dn}}\pi a^2$$ E3

where N defines the density of neutral species, m_n denotes the mass of neutral species, and v_dn represents the average relative velocity between the neutral elements and dust species. The resulting damping force also acts on the dust particles if the dust particles drift with drag force in the opposite direction to its motion.

1.6.4 Thermophoretic force (F_th)

The thermophoretic force occurs from the effect of temperature gradient in the neutral gas in the plasmas, and it is in the opposite direction to the temperature gradient. This force occurs due to transfer of momentum by the gas molecules from the hotter region to the colder portion of the gas. It can be written as:

$$F_{\mathit{th}}=\frac{16\sqrt{\pi }}{15}\frac{a^2{\kappa }_T}{{\nu }_{T,n}}\nabla T_n$$ E4

where v_T,n describes the thermal speed of the neutral gas of plasma, κ_T defines the translational effects in the λ, and T_n tells about the temperature of the neutrals. It can be occurred in the discharge by heating one of the electrodes. The force of gravity acting on the dust particles in plasma is balanced by the temperature gradient [11].

1.7 Application of dusty plasmas in industry

Plasma-based material processing technologies are extensively employed in the designing and commercialization of very large-scale integrated circuits (VLSI). Usually, chemically reactive plasmas are useful to sputter, etch, or otherwise alter the surface characteristics especially for silicon. Surface characteristics are approximately done at length scale of 0.2 μm wide and 4 μm deep in silicon films by such kind of mechanisms. The presence of dust is of critical concern to the microelectronic industry since particle contamination of semiconductor materials was estimated to account for more than 50% of device failures. Dust contamination diminishes the yield and recital characteristics of fabricated devices. Simply, the dust particles fall into the surface topographies of semiconductors either interfering with the etching process, preventing the adhesion of thin films or contaminating the final products. The occurrence of even the smallest dust particles became a crucial problem as the microelectronics industry moved to smaller and smaller structures.

1.7.1 Dust is a good thing

In those days, it was investigated that dust particles in plasma can have very interesting and useful properties, e.g., very small sizes, uniform size distribution, and chemical activity. There are many applications of plasma-produced particles. For example, large and active surface in catalysis is profitable. They are also essential in ceramic industry for sintering, in the modern technology of composite materials, and in fabrication of hard coatings [12] and solar cells [13]. Also, by injecting particles in plasma can furnish unique objects, like coated or layered grains with desired surface structure, color, and fluorescent properties. These particles are used as toners in copying machines [14] or in some optical devices [15].

1.7.2 Dust in plasma processing devices (dust is a bad thing)

During the last decade, in microelectronics industry, dust particles become the major cause of contamination and reduce the yield and performance of fabricated devices. In early 1990s, more than 50% devices were failed due to the particle contamination. The adhesion of thin films was reduced due to the submicron particles deposited on the surface, also causes dislocations. In semiconductor technology, the elimination of even smallest dust particles has become an urgent issue to develop smaller structures and thin films. Firstly, it was thought that the cause of most of the contamination is that the processed surfaces were handling in the clean rooms but soon it was seemed that plasma is the major source of dust particles and causes the loss of costly wafers. Now, the dust contamination is well controlled.

2. HNEMD algorithm and computational technique

We start, as usual, the Green-Kubo relations (GKR_S) for the hydrodynamic transport coefficients of uncharged particles [16]. It is well-known form and has been shown the standard GKR_S of fluids to the YDPL_S [17, 18, 19, 20, 21, 22]. The typical GKR_S used for the estimation of thermal conductivity of interacting dust particles for YDPL_S:

$$\mathrm{\lambda }=\frac{1}{3k_{\mathrm{B}}{\mathit{VT}}^2}\underset{0}{\overset{\infty }{\int }}⟨{\mathbf{J}}_Q\left(t\right).{\mathbf{J}}_Q\left(0\right)⟩\mathit{dt}$$ E5

where in Eq. (5), k_B is the Boltzmann’s constant, V is the system volume, T is the absolute temperature, and J_Q is the heat flux vector. The expression for the microscopic heat flux vector J_Q can be given by:

$${\mathbf{J}}_{Q}V=\sum _{i=1}^N{E}_i\frac{{\mathbf{p}}_i}{m}-\frac{1}{2}\sum _{i\ne j}\left({\mathbf{r}}_i-{\mathbf{r}}_j\right)\left(\frac{{\mathbf{p}}_i}{m}.{\mathbf{F}}_{\mathit{ij}}\right)$$ E6

In the above expression, F_ij is the total interparticle force on particle i due to j, r_ij = r_i − r_j are the position vectors, and P_i is the momentum vector of the ith particle. E_i is the energy of particle i and is given by the expression as:

$$E_i=\frac{{{\mathbf{p}}^2}_i}{2m}+\frac{1}{2}\sum _{i\ne j}{\varphi }_{\mathit{ij}}$$ E7

where ${\varphi }_{\mathit{ij}}$ is the Yukawa pair potential between particle i and j. Evans [23] has developed the non-Hamiltonian linear response theory (LRT) used for a moving system representing the equation of motion:

$${\stackrel{•}{\mathbf{r}}}_i=\frac{{\mathbf{p}}_i}{m}$$ E8

$${\stackrel{•}{\mathbf{p}}}_i=\sum _{J=1}^N{\mathbf{F}}_i+{\mathbf{D}}_i\left({\mathbf{r}}_i,{\mathbf{p}}_i\right).{\mathbf{F}}_e\left(t\right)-\alpha {\mathbf{p}}_i$$ E9

where ${\mathbf{F}}_i=\left(-\partial {\varphi }_{\mathit{ij}}/\partial r_i\right)$ is the total interparticle force acting on particle i and ${\mathbf{D}}_i={\mathbf{D}}_i\left({\mathbf{r}}_i,{\mathbf{p}}_i\right)$ is the tensor phase variable that describes the coupling of system to the fictitious external force field ${\mathbf{F}}_e\left(t\right)$ . A mechanical work is performed through the external applied force field ${\mathbf{F}}_e\left(t\right)$ , and thus, the equilibrium cannot be maintained. In the above Eq. (9), $\alpha$ is the Gaussian thermostat multiplier that maintains the system temperature and it is given as [16, 23]:

$$\alpha =\frac{{\sum }_{i=1}^N\left[{\mathbf{F}}_i+{\mathbf{D}}_i\left({\mathbf{r}}_i,{\mathbf{p}}_i\right).{\mathbf{F}}_e\left(t\right)\right].{\mathbf{p}}_i}{{\sum }_{i=1}^N{p}_i^2/m_i}$$ E10

The external force field parallel to the z-axis is of the form ${\mathbf{F}}_e\left(t\right)=\left(0,0,F_Z\right)$ ; therefore, the thermal conductivity is calculated as:

$$\begin{array}{ll}\mathrm{\lambda }& =\frac{V}{3k_B{T}^2}\underset{0}{\overset{\infty }{\int }}⟨{\mathbf{J}}_{Qz}\left(t\right){\mathbf{J}}_{Qz}\left(0\right)⟩\mathit{dt}\\ & =\underset{F_z\to 0}{lim}\underset{t\to \infty }{lim}\frac{-⟨{\mathbf{J}}_{Q_z}\left(t\right)⟩}{{\mathit{TF}}_z}\end{array}$$ E11

where ${\mathbf{J}}_{Q_Z}\left(t\right)$ is the z-component of the heat flux vector and the external force field ${\mathbf{F}}_e\left(t\right)=\left(F_Z\right)$ .

3. Numerical results and discussion

This section provides the outcomes of thermal conductivity of 3D complex dusty plasmas by using HNEMD simulations over suitable plasma couplings Γ (≡1, 200) and screening strengths κ (≡1.4, 4) at constant external force strength of F_ext (≡0.005). It is noted that we have already reported our similar results with higher system sizes [13] and with different varying force fields [13]. In this present work, we have reported our HNEMD outcomes for different low to intermediate system sizes at fixed force field.

Figures 1, 2,3, 4 show our main outcomes of plasma thermal conductivity (λ) by employing HNEMD approach. Here, the thermal conductivity is normalized by plasma frequency (ω_p) as λ₀ = λ/nk_Bω_pa_ws, or by the Einstein frequency (ω_E) as λ^* = λ/√3nk_Bω_Ea_ws of SCCNPs, at the normalized external field strength F^* = (F_z)(a_ws/J_Q), where a_ws is radius of Wigner-Seitz (WS) radius with n being the equilibrium particle number density, k_B is Boltzmann constant and J_Q. It should be noted that these normalizations have been employed for classical Coulomb one-component plasmas (COCPs) [24] and SCCNPs [7].

Figure 1.

Comparison of normalized plasma thermal conductivity (λ₀), computed by various numerical approaches for plasma coupling states Γ ≡ (1, 200). Results investigated by Shahzad and He for homogenous nonequilibrium MD (HNEMD) [13], Salin and Caillol for equilibrium MD (EMD) [21], Donko and Hartmann for inhomogenous NEMD [17]. Shahzad and He for homogenous perturbed MD (HPMD) [18] and Faussurier and Murillo for variance procedure (VP) [25], (a) for N = 256, (b) N = 1372 and at κ = 1.4.

Figure 2.

Comparison of normalized plasma thermal conductivity (λ₀), computed by various numerical approaches for plasma coupling states Γ ≡ (1, 300). Results investigated by Shahzad and He for homogenous nonequilibrium MD (HNEMD) [13], Salin and Caillol for equilibrium MD (EMD) [21], Donko and Hartmann for inhomogenous NEMD [17]. Shahzad and He for homogenous perturbed MD (HPMD) [18] and Faussurier and Murillo for variance procedure (VP) [25], (a) for N = 500, (b) N = 864 and at κ = 3.

Figure 3.

Comparison of normalized plasma thermal conductivity (λ₀), computed by various numerical approaches for plasma coupling states Γ ≡ (1, 300). Results investigated by Shahzad and He for homogenous nonequilibrium MD (HNEMD) [13], Salin and Caillol for equilibrium MD (EMD) [21], Shahzad and He for homogenous perturbed MD (HPMD) [18] and Faussurier and Murillo for variance procedure (VP) [25], (a) for N = 500, (b) N = 864 and at κ = 4.

Figure 4.

Trend of normalized potential energy (P.E.) with four values of screening parameters (κ = 1, 2, 3, and 4) for different coupling states Γ ≡ (1, 5, 10, 20, 50, 100, 200, and 300) and at N = 256.

Diverse series of the plasma λ₀ subsequent to a decreasing series of external force field F^* are computed to establish the linear system of the SCCNPs under the influence of the normalized force field strength. The current HNEMD outcomes allow investigation for the complete range of plasma parameters (Γ, κ) with variation of external force field F^*. In this work, a feasible suitable value of the external force field strength F^* (=0.005) for the computation of the steady state values of the plasma normalized thermal conductivity is to be selected, for small varying practical system size. This feasible suitable external force field provides the steady state plasma thermal conductivity estimations, which are satisfactory over the whole range of the plasma state points (Γ, κ).

Figures 1, 2, 3 display that the computed plasma thermal conductivity is in acceptable conformity with former HNEMD investigations by Shahzad and He [13], EMD calculations of Salin and Caillol [21], inhomogenous NEMD estimations of Donkó and Hartmann [17], homogenous perturbed molecular dynamics simulations (HPMD) measurements of Shahzad and He, and theoretical predictions of Faussurier and Murillo for variance procedure (VP) [18, 25]. It can be seen from Figure 1 that our results are slightly lower as compared to earlier known numerical results based on different numerical techniques, at lower Γ. However, the present results are well matched with earlier results for intermediate to higher Γ at two different system sizes N = (256 and 1372) and it is clearly shown that our results are very close EMD and HNEMD results. It is observed from Figure 2 that HNEMD results are between EMD (at lower N) results and HNEMD (at higher N) computations at low value of Γ but our outcomes are satisfactorily matched with earlier results at intermediate and higher Γ. It can be seen from Figure 2 that our HNEMD data for low intermediate system size it is increasing behavior at higher Γ, confirming earlier HPMD results [18]. It is observed that the deviation of data from earlier known measured data based on different methods of EMD, HPMD, and InNEMD is also computed and the outcomes of plasma λ₀ are within range of ∼3–22% for EMD, ∼7–20% for NEMD, and ∼10–35% for HPMD. It is noted that some of data points are far away from present data that are not mentioned here but most of data points are within limited statistical range, as expected. At higher screening κ = 4, it is examined from Figure 3 that the present results are definitely lower as compared to earlier EMD computations of Salin and Caillol and HNEMD estimations at higher N of Shahzad and He. Moreover, it can be noted that the present outcomes are slightly lower at intermediate Γ and well matched at higher Γ, confirming earlier results.

It is suggested from these figures that measured outcomes are satisfactory well matched with previous outcomes at intermediate to high Γ; however, some results diverge at the lower Γ points but all within statistical unlimited uncertainty range. Figures 1, 2, 3 show that the presented HNEMD method may precisely calculate the plasma thermal conductivity of strongly coupled complex plasmas. We have shown that the present method has good performance and its accuracy is very close to earlier EMD and InHNEMD methods. It is concluded that our outcomes depend on the plasma parameters of Coulomb coupling and Debye screening strength, confirming earlier simulations. Moreover, it is shown that the position of minimum value of thermal conductivity shifts toward higher Γ with an increase of screening κ, as expected. Presently, we have demonstrated our results for a wide range of plasma parameters, ranging from nonideal gaseous state to strongly coupled range. It is noted that the extended HNMED method is excellent for lower system sizes with constant external force field strength, where signal-to-noise ratio is acceptable for equilibrium plasma thermal conductivity.

4. Summary

Plasma thermal conductivity of SCCNPs system was computed over a suitable domain of plasma couplings (1 ≤ Γ ≤ 300) and screening strength (1.5 ≤ κ ≤ 4) by employing constant external force field strength through HNEMD approach. It is shown that our HNEMD outcomes are in reasonable agreement with the earlier outcomes measured from EMD, HNEMD, InHNEMD, HPMD, and VP approaches for SCCNPs. New computations show that the minimum values of thermal conductivity exist at same values of plasma coupling Γ and it shifts toward higher Γ with an increase of screening κ, as expected in earlier numerical approaches. It has been revealed that the plasma thermal conductivity depends on plasma parameters (Γ, κ) in 3D complex dusty systems that illustrate earlier results of SCCNPs. In this study, the HNEMD method is a mostly dominant numerical approach, which occupies fast computations of plasma thermal conductivity, for small to intermediate system sizes. This chapter provides the understanding and investigation of nonlinear regime of the SCCNPs for a suitable low value of external force field strength. In future, thermal conductivity of complex plasma can be calculated by applying external magnetic field or an electric field strength and it can be applied to other systems (Coulomb, polymer, or ionic).

Acknowledgments

The authors thank Z. Donkó (Hungarian Academy of Sciences) for providing his thermal conductivity data of Yukawa Liquids for the comparisons of our simulation results, and useful discussions. We are grateful to the National Advanced Computing Center of National Center of Physics (NCP), Pakistan, for allocating computer time to test and run our MD code.

Abbreviations

SCCNPs	strongly coupled complex nonideal plasmas
Γ	Coulomb coupling
κ	Debye screening length
F*	external force field strength
HNEMD	homogenous nonequilibrium molecular dynamics
NEMD	nonequilibrium molecular dynamics
MD	molecular dynamics
InHNEMD	inhomogenous nonequilibrium molecular dynamics
SCP	strongly coupled plasma
EMD	equilibrium molecular dynamics
λ	thermal conductivity
λ0	normalized thermal conductivity
PBCs	periodic boundary conditions
VP	variance procedure
HPMD	homogenous perturbed MD
N	number of particles