## 1. Introduction

The exact numerical investigation of transport properties of complex liquids is a fundamental research task in the field of thermophysics, as various transport data are closely related with setup and the confirmation of equations of state. A reliable knowledge of transport data is also important for optimization of technological processes and apparatus design in various engineering and science fields (incl. thermoelectric devices) and, in particular, when provision of precise data for parameters of heat, mass, and momentum transport is required [1–3]. In thermophysical properties of fluids, chemical properties remain unaffected, but physical properties of material are changed by variable temperature, composition, and pressure. These properties of simple and complex liquids explain the phase transition [4]. These fluids can be examined experimentally, theoretically, and by simulation techniques. Thermophysical properties (thermodynamics and transport coefficients) include thermal conductivity, thermal expansion, thermal radiative properties, thermal diffusivity, enthalpy internal energy, Joule-Thomson coefficients, and heat capacity, as well as, thermal diffusion coefficients, mass coefficients, viscosity, speed of sound, and interfacial and surface tension. Thermophysical properties of gases and liquids, such as hydrogen, H_{2}; oxygen, O_{2}; nitrogen, N_{2}; and water, H_{2}O are different from ideal gas at high pressure and low temperature. Specific models are required for the calculation of these properties in the widest range of pressure and temperature. Different fluids, such as gases and liquids are used as a power generation source in different power plants. For example, heavy water, steam, air, and different gases are used for power generation in nuclear power plants, gas turbine plants, and internal combustion plants. Also, for cooling in refrigerators and fast nuclear reactors, ammonia and sodium in liquid phase are used as a cooling agent.

### 1.2. Dusty plasma

Nowadays dusty plasma refers to as complex plasma in analogy to the condensed matter field of “complex liquids” in soft matter (colloidal suspensions, polymers, surfactants, etc.). The dust particles combine physics of nonideal plasmas and condensed matter, and this field has played an important role in both newly system designs and advance development micro- and nanotechnology. This complex plasma system has four components, i.e., ions, electrons, neutral atoms, and dust particles with high charges, as compared to other species, which are responsible for the extraordinary plasma properties. The study of complex (dusty) plasmas reveals rich variety of interesting phenomena and extends knowledge on fundamental aspects of plasma physics at the microscopic level. Among these, the freezing (gaseous-liquid-solid) phase transition is of particular interest. Complex plasma is called strongly coupled plasma, in which thermal energy (kinetic energy) of nearest neighbors is much smaller than their Coulomb interaction potential energy, whereas plasma is called weakly coupled when Coulomb interparticle potential energy of nearest particles is much smaller than their kinetic energy [5–7].

Plasma is the fourth state of matter, and usually, it is said, that there are three states of matter, but another state was also found to exist, named as plasma. Irving Langmuir (American physicist) defined plasma as “it is a quasi-neutral gas of charged and neutral particles, which exhibits collective behavior,” and he got the Nobel Prize in 1927 because firstly he was using the term plasma [8]. In this definition, quasi-neutral means that plasma is electrically neutral and has approximately equal ion and electron density (*n*_{i} ≈ *n*_{e}). The term collective behavior shows, that due to Coulomb potential or electric field, plasma’s particles collide with each other. Simply, plasma is an ionized gas. In a gas, sufficient energy is given to eject free electrons from atoms or molecules, and, as a result, ions and electrons, both species, coexist. There are several ways to convert a gas into plasma, but all include pumping the gas with energy. For instance, plasma can be created due to a spark in a gas. Usually, neutral plasmas are relatively hot, such as the solar corona (1000000 K), a candle flame (1000 K), or the ionosphere around our planet (300 K). In this universe, the main source of plasma is sun in which the large number of electrons of hydrogen and helium molecules is removed. So, the sun is a great big ball of plasma like other stars.

#### 1.2.1. Types of plasma

There are different types of plasma that are described by many characteristics, such as temperature, degree of ionization, and density.

### 1.2.1.1. Cold plasma

In laboratory, in the positive column of a glow discharge tube, there exists plasma in which the same number of ions and electrons is present. When gas pressure is low, collision between electrons and gas molecules is not frequent. So, nonthermal equilibrium between energy of electrons and gas molecules does not exist. So, energy of electrons is very high as compared to gas molecules and the motion of gas molecules can be ignored. We have *T*_{e} ≫ *T*_{i} ≫ *T*_{g}, where *T*_{e}, *T*_{i}, and *T*_{g} represent temperature of electron, ion, and gas molecules and such type of plasma is known as cold plasma. In cold plasma, the magnetic force can be ignored and only the electric force is considered to act on the particles. In cold plasma technique, cold gases are used to disinfect surfaces of packaging or food products. Vegetative microorganisms and spores on packaging materials can be inactivate at temperature below 40°C. This process can have a clear advantage compared to heat treatment for temperature-sensitive products. Also it can reduce the amount of water used for disinfection of packaging materials. Because cold plasma is in the form of gas, so, the irregularly shaped packages, such as bottles can be treated easily as compared to UV or pulsed light where shadowing occurs.

### 1.2.1.2. Hot plasma

When gas pressure is high in the discharge tube, then electrons collide with gas molecules very frequently and thermal equilibrium exists between electrons and gas molecules. We have *T*_{e} ≅ *T*_{i}. Such type of plasma is known as hot plasma and it is also known as thermal plasma. Hot plasma is one which approaches to a local thermodynamic equilibrium (LTE). Atmospheric arcs, sparks, and flames are used for the production of such type of plasmas.

### 1.2.1.3. Ultracold plasma

If plasma occurs at temperature as low as 1 K, then such type of plasma is known as ultracold plasma, and it can be formed by photoionizing laser-cooled atoms and pulsed lasers. In ultracold plasmas, the particles are strongly interacting because their thermal energy is less than Coulomb energy between neighboring particles [9].

### 1.2.1.4. Ideal plasma

There are mainly two types of plasma according to plasmas’ ideality and properties study, nonideal plasmas (weakly coupled and strongly coupled plasmas) and ideal plasmas (very weakly coupled plasmas). Whenever the kinetic energy of plasma is much larger than the potential energy and plasma has a low temperature and high density, then such type of plasma is known as ideal plasma. Ideal plasma is one in which Coulomb collisions are negligible. If the average distance among the interacting particles is large, then the interaction potential can be ignored due to this large-distance ideal plasma that does not have any arrangement of particles [2].

### 1.2.1.5. Nonideal (complex) plasma

Nonideal plasmas are often found in nature, as well as, in technological services. They can be shown as electron plasma in solid and liquid metals and electrolytes, the superdense plasma of the matter of white dwarfs, the sun and the interiors (deep layers) of the giant planets of the solar system, and astrophysical objects, whose structure and evolution are defined by plasma characteristics [2]. Further examples of nonideal plasmas are brown dwarfs, laser-generated plasmas, capillary discharges, plasma-opening switches, high-power electrical fuses, exploding wires, etc. On the bases of Coulomb coupling, nonideal plasma can be divided into two families: weakly coupled plasma (WCP, Г < 1) and strongly coupled plasma (SCP, Г ≥ 1).

Nonideal complex plasmas are found in daily life and can be found in processing industries to manufacture many products that we deal in our everyday life directly or indirectly at moderate temperature, such as plastic bags, automobile bumpers, airplane turbine blades, artificial joints, and, most importantly, in semiconductor circuits. Moreover, nonideal plasmas (terrestrial plasmas) are not hard to find. They occur in gas-discharge lighting, such as neon lighting used for commercial purposes and fluorescent lamps, for instance, compact fluorescence light sources, which have a higher performance than the traditional incandescent light sources, a variety of laboratory experiments, and a growing array of industrial processes. Modern display methods contain plasma screens, in which small plasma discharges are used to stimulate a phosphor layer, which then emits light [2, 7].

#### 1.2.2. Complex (dusty) plasma

Dust is present everywhere in the universe and mostly it is present in solid form. It is also present in gaseous form, which is often ionized, and thus the dust coexists with plasma and forms “dusty plasma.” In dusty plasma, dust particles are immersed in plasma, in which ions, electrons, and neutrals are present. These dust particles are charged and then affected by electric or magnetic fields and can cause different changes in the properties of plasma. The presence of dust component gives rise to new plasma phenomena and allows study of fundamental aspects of plasma physics at the microscopic level. Dust particles are charged due to the interaction between dust particles and the surrounding plasmas. Due to this interaction, grains are charged very rapidly. The charge on grains depends on the flow of ions and electrons. These charged grains enhance plasma environments, for example, setting up space charges. Also, to determine the charge on dust grain, it is assumed, that a spherically symmetric isolated dust grain is injected in plasma and only the effect of ion and electron is considered. Moreover, there are many other charging processes, such as secondary emission, electron emission, thermionic emission, field emission, radioactivity, and impact ionization. Complex plasma is condensed plasma characterized by strong interaction between existing molecules and atoms; it is also called strongly coupled complex plasma. Dusty plasma is complex plasma which includes many components: ions, electrons, neutral particle, and dust particles. Last 20–25 years, strongly coupled plasmas were mainly studied theoretically, due to lack of suitable laboratory tools and equipment. However, experimental strongly coupled plasma studies became more common with the discovery of ways to find dusty plasma [10], laser-cooled ion plasmas in a penning trap [11], and ultracold neutral plasmas [12]. Plasma systems can be treated theoretically in a straightforward way in the extreme limits of both weakly coupled and strongly coupled plasmas [2].

The main goals of this chapter are to study thermal conductivity (*λ*) along with lattice correlation (long-range crystalline order) at the corresponding plasma states and to extend the set of plasma states (Γ, *κ*) by using the same method as introduced by Shahzad and He [3] with different system sizes (*N*). The effect of external force field strengths on Ψ(*t*) and corresponding *λ* values of complex Yukawa liquids of dust particles is another interesting task that is calculated under near-equilibrium condition.

## 2. HNEMD model and simulation approach

In this section, we will introduce theoretical background needed in this work. We start by introducing the model system, which is used in our HNEMD simulations. We consider a cubic box of edge length *L* and have *N* number of particles or millions of atoms. In MD technique, the range of particle number is chosen as *N* = 500–1000 [13]. Periodic boundary conditions (PBCs) are used for selection of the size or dimension of the box. Practically, PBCs avoid the surface size effects. The particles present in the box interact with each other with known interaction potential. This potential may be Yukawa, Coulomb, and Lennard-Jones potential depends on the type of the system considered.

The particles interact through screened Coulomb potential, which depends on the physical parameters and the background plasmas. Average interparticle interaction is frequently considered to be isotropic and basically repulsive and approximated by Yukawa interaction potential [1–3]. Yukawa model has been employed in many physical and chemical systems (for instance, biological and pharmaceutical sciences, colloidal and ionic systems, space and environment sciences, physics and chemistry of polymers and materials, etc.) [1–8]. In the present case, the interaction of potential energy of particles is in Yukawa form:

where charge on dust particle is *Q*_{d}, magnitude of interparticle division is *r*, and Debye length *λ*_{D} accounts for the screening of interaction by other plasma species. Due to the long-range interaction between the particles, Yukawa potential energy cannot be solved directly. Ewald sum method is used with periodic boundary conditions to calculate Yukawa potential energy, force, and heat energy current. Improved nonequilibrium molecular dynamics (NEMD) method proposed by Evans has been employed to estimate thermal conductivity of strongly coupled complex (dusty) plasmas. During the simulation of Yukawa systems, sufficient number of particles *N* is to be selected to study the size effect of system. It comes to know, that there is no effect of system size on thermal conductivity or on any other properties under limited statistical uncertainties. Negative divergence of Yukawa potential *F* = (−∇*φ*) and Newton’s equation of motion integrated by predictor-corrector algorithm are used for calculation of force exerted by the particles on each other or in minimum image convection [14].

In HNEMD technique, in order to measure thermal conductivity and nonlinearity of NICDPs, the system will be perturbed by applying the external field along *z*-axis. In three-dimensional systems, we use standard Green-Kubo relations (GKR) for calculation of thermal conductivity coefficient of uncharged particles [15]:

Here, **J**_{Q} is vector of heat flux, *V* is the volume, *T* is the temperature of system, and *k*_{B} is Boltzmann constant. At microscopic level, heat flux vector has value:

In this equation, r_{ij} = r_{i} − r_{j} is the position vector, F_{ij} is the force of interaction on particle *i* due to *j*, and p_{i} represents the momentum vector of the *i*th particle. The energy *E*_{i} of *i*th particle is given by:

where ∅_{ij} is Yukawa pair potential between particles *i* and *j*. According to non-Hamiltonian dynamics, generalization of linear response theory, proposed by Evans, for system moving with equations of motion and recently detailed understanding of Ewald-Yukawa sums [1–3] allows to present thermal conductivity as:

where J_{Qz} is *z*-component of heat energy flux vector for strongly coupled complex (dusty) plasma liquids. Thermal conductivity coefficient for charged particle of plasma according to GKR is given by Eq. (5), and further detail is provided in Ref. [3]. Plasma states of Yukawa systems can be illustrated fully by three reduced parameters: plasma coupling parameter Γ = (*Q*^{2}/4πε_{0})(*a*_{ws}*k*_{B}*T*), screening parameter *κ* ≡ *a*_{ws}/*λ*_{D}, and reduced external force *F*^{ext} for HNEMD model, where *a*_{ws} is Wigner-Seitz radius, *Q* is charge on dust particle, and *ε*_{0} is permittivity of free space [1–3]. Gaussian thermostat is used to control temperature of systems [14]. Simulation time step is d*t* = 0.001/*ω*_{p}, where *ω*_{p} = (*nQ*^{2}/ε_{0}*m*)^{1/2} is dust plasma frequency with *m* is dust particle’s mass and *n* is number density. Reported simulations are performed between 3.0 × 10^{5}/*ω*_{p} and 1.5 × 10^{5}/*ω*_{p} time units in the series of data recording of thermal conductivity (*λ*) [16, 17].

## 3. Computer simulation outcomes

### 3.1. Particle lattice correlation

The structural information of Yukawa system is given by lattice correlation. For the calculation of lattice correlation, density of given material at point *r* can be calculated as:

where *ρ* is density of system, *N* is number of particles, *δ* represents distribution of particles, and *r*_{j} is the position of particle corresponding to particle at position *r*. Eq. (6) gives information about the system being in ordered state and then it may be in solid or crystal form depending on *ρ*. The lattice correlation equation according to Fourier transform is:

System arrangement (ordered or disordered) is calculated by simulation based on Eq. (7). When value of lattice correlation approaches to |Ψ| ≈ 1, then the system will be in ordered state, and if value becomes |Ψ| ≈ 0, then the system will be in liquid or gas (nonideal gas) state. In Eq. (7), *k* is lattice correlation vector for ordered state, and its value is different for different lattice structures. Its value for face-centered cubic (FCC) is *k* = 2π/(1,−1,1)*l*, for body-centered cubic (BCC) is *k* = 2π/(1,0,1)*l*, and for simple cubic (SC) is *k* = 2π/(1,0,0)*l*; here, *l* is edge length [14].

Lattice correlation was examined in 3D NICDPs in the limit of appropriate constant near equilibrium external force field strength *F*^{ext} = 0.005. **Figure 1** illustrates lattice correlation in NICDPs versus simulation time at normalized *F*^{ext} = 0.005. In this case, additional parameter includes the heat energy flux J_{Q} and the external force field strength *F*^{ext}(*t*) = (0,0,*F*_{z}) is selected along *z*-axis, in the limit of *t* → ∞, and its normalized value *F*^{ext} = (*F*_{z})(*a*_{ws}/J_{Q}) [17].

### 3.2. Normalized thermal conductivity

We now turn attention to the key results obtained through HNEMD simulations. Obtained computer-simulated data confirm, that thermal conductivity of Yukawa system can be calculated with satisfactory statistics by an extended HNEMD approach. **Figures 2**–**5** display the main results calculated from HNEMD method for various plasma states for Yukawa liquids at *κ* = 1.4 and *κ* = 2 and *κ* = 4, respectively. HNEMD simulation is used to compute the thermal conductivity normalized by plasma frequency (*ω*_{p}) as *λ*_{0} = *λ*/*nk*_{B}*ω*_{p}*a*_{ws}, or by Einstein frequency (*ω*_{E}) as *λ*^{*} = *λ*/√3*nk*_{B}*ω*_{E}*a*_{ws} of YDPLs, at the normalized external field strength. These normalizations of transport properties, including *λ*_{0}, were widely used in earlier studies of one-component complex plasma (OCCP) [18] and NICDPs [1–3, 19–21]. HNEMD method is employed to investigate *λ*_{0} of 3D NICDPs at reduced external force field *F*^{ext} = 0.005 over suitable domain of plasma parameters of coupling (1 ≤ Γ ≤ 300) and screening (1 ≤ *κ* ≤ 4).

Different sequences of *λ*_{0} corresponding to decreasing sequence of *F*^{ext} are calculated to determine the linear regime of YDPLs under the action of reduced force field strength. The earlier data of thermal conductivities of YDPLs has been limited to the small values of plasma parameters at different *F*^{ext}. The present HNEMD simulation enables study over the whole domain (Γ, *κ*) of plasma-state parameters with variation of *F*^{ext}. In this case, possible low value of the force field strength *F*^{ext} = 0.005 for determination of the near-equilibrium values of Yukawa thermal conductivity is to be chosen, for small reasonable system size. This possible reasonable external force field gives the near-equilibrium thermal conductivity measurements, which are acceptable for the whole domain of plasma parameters (Γ, *κ*).

**Figures 2**–**5** show, that measured thermal conductivity is in satisfactory agreement with earlier HPMD simulations by Shahzad and He [1], inhomogeneous NEMD computations by Donko and Hartmann [20], and EMD measurements by Salin and Caillol [19]. The present results are also higher than Salin and Caillol [19] results at lower Γ = 1 and 2, for *κ* = 1.4 (*N* = 500 and 864). The minimum value of *λ*_{0} is *λ*_{min} ≈ 0.4533 at Γ = 20 and *κ* = 1.4. Deviation of the present data from the earlier calculated results based on different techniques of EMD, HPMD, and NEMD is also calculated. It is observed, that the results of *λ*_{0} are within the range of ~7–50% for EMD, ~10–16% for NEMD, and ~10–40% for HPMD. Moreover, **Figure 3(a)** shows, that obtained thermal conductivity at *F*^{ext} = 0.005 for *N* = 864 is in good agreement with HPMD simulation of Shahzad and He [1], but it is noted, that our results are slightly greater than EMD of Salin and Caillol [19], NEMD of Donko and Hartmann [20], and VP of Faussurier and Murillo [21] at lower values of Γ. Deviation of the present data from EMD, HPMD, and inhomogeneous NEMD is within the range of ~4–14%, ~3–16%, and ~4–30%. The minimum value of *λ*_{0} = 0.4331 at Γ = 20 and *κ* = 1.4.

This comparison shows, that our data remain within the limited statistical uncertainty range. Panels (b) of **Figures 2**–**5** compare the present simulation results of thermal conductivity, normalized by reference data, calculated here from HNEMD approach for different sets of external force field strengths with reference set of data and earlier known simulation data of HPMD, EMD, inhomogeneous NEMD, and VP techniques [1, 19–21]. A series of different sequences of HNEMD simulations are performed with *N* = 500 and 864 at various normalized *F*^{ext} values for varying screening parameters. In our case, results measured by Shahzad and He [3] with *N* = 13500 at *F*^{ext} = 0.005 are taken as reference set of data *λ*_{REF} and are shown panels (b) of the respective figures. It is important, that computationally noted linear regime is traced between 0.001 ≤ *F*^{ext} and ≤ 0.009, which depends on the plasma state points (Γ, *κ*). It observed, that differences between the simulation data calculated by different authors are large at some state points and differences with their own present data sets are much smaller. Plasma’s thermal conductivity is generally overpredicted within ~3–17% (~8–15%), ~5–30% (~5–38%), and ~5–40% (~3–17%) relative to the data of EMD by Salin and Caillol [19], inhomogeneous NEMD by Donko and Hartman [20], and HPMD by Shahzad and He [1], respectively, for *N* = 500 (864). It is concluded, that most of data points of presented results fall under ±15 range around the reference data.

For the cases of *κ* = 2 and *κ* = 4, the minimum value of *λ*_{0} is *λ*_{min} ≈ 0.3413 (for *κ* = 2) at Γ = 50 and *λ*_{min} ≈ 0.2218 (for *κ* = 4) at Γ = 200. The comparison of the present data with the earlier investigated results obtained from different techniques of EMD, HPMD, and NEMD is also shown in respective figure panels (a) for *κ* = 2 and *κ* = 4. It is observed, that the results of *λ*_{0} are within the range ~1–15% (~2–13%, for *κ* = 4) for EMD, ~2–20% (~3–13%, for *κ* = 4) for HPMD, ~15–22% for inhomogeneous NEMD, and ~1–28% for HPMD. It is observed, that for *κ* = 4, thermal conductivity agrees relatively well with the simulation data in general within ~5–35% for EMD and ~2–23% for HPMD.

It is concluded from figures, that obtained results agree well with earlier results at intermediate and high Γ values; however, some data points deviate at the lower Γ values. **Figures 2**–**5** depict, that extended HNEMD approach can accurately predict thermal conductivity of Yukawa system (dusty plasma). We have used the present developed homogeneous NEMD method, which has an excellent performance; its accuracy is comparable to that of EMD and inhomogeneous NEMD techniques. The first conclusion from above **Figure 5** is that thermal conductivity depends on plasma parameters Γ and *κ*, confirming the earlier numerical results. Furthermore, it is examined, that the position of minimum value of *λ*_{min} shifts toward higher Γ with increase in *κ*, as expected. The minimum value of *λ*_{0} decreases with increasing *κ*, as *λ*_{0} = 0.4533 at *κ* = 1 to *λ*_{0} = 0.2218 at *κ* = 4 and *λ*_{0} = 0.4331 at *κ* = 1 to *λ*_{0} = 0.2099 at *κ* = 4 for *N* = 500 and 864, respectively. Our numerical data are considerably more comprehensive covering the full range of coupling strengths from nearly nonideal gaseous state to strongly coupled liquid (SCL) as shown in **Figure 1**. The present approach has shown excellent results of *λ*_{0} at steady-state value of *F*^{ext} at lower, intermediate, and higher coupling values, and it also gives more comparable performance of normalized thermal conductivity at lower *N*. This approach yields the practical accuracy, compared to EMD, for relatively small total simulation times due to the act of finite nonzero external force field; the signal-to-noise ratio of thermal response is high. It is significant from our simulation data, that the external force field strength of *F*^{ext} increases with increase in *κ*.

## 4. Summary

Thermal conductivity of NICDP system was investigated for wide range of Coulomb coupling parameter (1 ≤ Γ ≤ 300) and screening parameter (1 ≤ *κ* ≤ 4), applying external force field, by using HNEMD method. This HNEMD simulation method reveals, that our present results are in good agreement with the earlier results obtained by equilibrium MD and homogeneous and inhomogeneous NEMD simulations for NICDPs. It is confirmed, that lattice correlation is not affected by system size, while lattice correlation decreases with increment of *κ* and at high temperature (1/Γ). It is confirmed from presented HNEMD simulation results, in which normalized thermal conductivities follow simple universal (temperature) scaling law. Observations show, that the minimum value of thermal conductivity shifts toward higher Coulomb couplings with increasing screening strength, confirming earlier numerical results. It has been shown, that thermal conductivity depends on both temperature (plasma coupling) and density (screening) in 3D Yukawa systems, which shows previous data for NICDPs. Presented HNEMD approach was particularly a powerful numerical technique, which involves fast calculation of thermal conductivity, on small and intermediate system sizes, in contexts, applicable to the study of matter at microscopic level for 3D NICDPs. The second major contribution of this chapter is that it provides, for the first time, understanding and determination of non-Newtonian behavior of Yukawa liquid. It is important to note, that thermal conductivity is in increasing behavior at higher values of external force field. These indications show, that increasing behavior of field dependence of thermal conductivity decreases with increasing *κ* and decreasing Γ. Described simulation technique can be helpful for estimation of thermal conductivity regularities in novel liquid, organic (polymer), and multiphase solid-state thermoelectric materials.