Introductory Chapter: A Novel Approach to Compute Thermal Conductivity of Complex System

1. Introduction

Thermophysical properties of complex fluids materials (dusty plasmas) have been very actively investigated both experimentally and by computer simulations. These properties describe the physical and chemical behavior of material that remains in terms of fluids formation. Among all of the properties, thermal conductivity is most important; it is an intrinsic property of fluid materials. It is a difficult property from computational point of view because it is sensitive to the internal energy of the system. Its molecular-level observations are becoming more important in dusty plasmas. Dust is abundant in nature. Dust particles coexist in plasmas and then form a complex dusty plasma. Dusty plasmas are ionized gases that contain particulates of condensed matter. Dusty plasmas are of interest as a non-Hamiltonian system of interacting particles and as a means to study generic fundamental physics of self-organization, pattern formation, phase transitions, and scaling. Their discovery has therefore opened new ways of precision investigations in many-particle physics.

2. Dusty plasmas

Dusty plasmas are ionized gases that contain particulates of condensed matter. These particles have different sizes ranging from tens of nanometers to hundreds of microns. The dynamics of these massive charged particles happens at slower timescales than the ordinary plasma ions because charge to mass ratio (e/m) is orders of magnitude smaller that the corresponding (e/m) of either the electrons or ions. They may be in the shape of rods, irregularly shaped pancakes or spheres. They are made up of dielectric, e.g., SiO₂ or Al₂O₃, or conducting materials. Even though the particles are normally solid, they might also be liquid droplets or fluffy ice crystals [1].

3. History of dusty plasmas

Despite almost a century–long history—the first observations of dust in discharges have been reported by Langmuir in 1924. After that, Lyman Spitzer along with Hannes Alfven recognized that dust in the universe was not merely a hindrance to optical observation, but that it was an essential component of the universe. One of the most exciting events in the field of dusty plasmas occurred in early 1980 during the Voyager 2 flyby of planet Saturn. In 2005, Cassini spacecraft took new and improved images of spokes with detail that would provide a better understanding of their origin. In 1992, the European spacecraft Ulysses flew by the planet Jupiter and detected the dust particles and measured their masses and impact speed. The current enormous interest in complex plasmas started in the mid-1980s, triggered by laboratory investigation of thermal conductivity of dusty plasmas [2].

4. Types of dusty plasmas

4.1. Weakly vs. strongly coupled dusty plasmas

The vital attention in the investigation of dusty plasmas is whether the particles are in the weakly or strongly coupled state. The description as weakly or strongly coupled denotes the subject of whether the particles average potential energy, due to nearest neighbor interactions, are smaller or larger than their average thermal energy. Coulomb coupling parameter, Γ is defined as the ratio of the interparticle Coulomb potential energy to the (thermal) kinetic energy of the particles:

$$\mathrm{\Gamma }=\frac{q^2}{\mathit{rT}}\mathit{exp}\left(\frac{-r}{{\mathrm{\lambda }}_D}\right)$$E1

where r = (3/4πn)^1/3 is the average interparticle spacing [3].

5. Dusty plasma physics

The brief discussion of the basic principles of dusty plasma physics is given in the following sections.

5.1. Debye shielding

The Debye length is an important physical parameter in a plasma. It provides the distance over which the influence of the electric field of an individual charged particle is felt by other charged particles (such as ions) inside the plasma. The charged particles actually rearrange themselves in order to shield all electrostatic fields within a Debye distance. In dusty plasmas, the Debye length can be defined as follows [4].

$${\mathrm{\lambda }}_D=\frac{{\mathrm{\lambda }}_{\mathit{De}}.{\mathrm{\lambda }}_{\mathit{Di}}}{\sqrt{{\mathrm{\lambda }}_{\mathit{De}}^2+{\mathrm{\lambda }}_{\mathit{Di}}^2}}$$E2

where λ_De = √KT_e /4πn_ee² and λ_Di = √KT_i /4πn_ie² are the Debye lengths associated to electrons and ions, respectively.

5.2. Macroscopic neutrality

Dusty plasmas are characterized as a low-temperature ionized gas whose constituents are electrons, ions, and micron-sized dust particulates. The presence of dust particles (grains) changes the plasma parameters and affects the collective processes in such plasma systems. In particular, the charged dust grains can effectively collect electrons and ions from the background plasma. Thus, in the state of equilibrium, the electron and ion densities are determined by the neutrality condition [5], which is given by

$${\mathit{en}}_{i\mathit{0}}‐{\mathit{en}}_{e\mathit{0}}+{\mathit{qn}}_{d\mathit{0}}=0$$E3

5.3. Inter-grain spacing

In a multicomponent dusty plasma, inter-grain spacing is very important to distinguish between dust in plasma and dusty plasma. Like the dust grain radius, b the inter-grain spacing, r is usually much smaller than the Debye length λ_D. For dust in plasma, r > λ_D and the dust particles are completely isolated from their neighbors. For dusty plasmas, r ≤ λ_D and the dust particles can be considered as massive point particles like multiple-charged negative (or positive) ions in a multispecies plasma where the effect of neighboring particles can be significant.

5.4. Coulomb coupling parameter

Charged dust grains can be either weakly or strong correlated depending on the strength of the Coulomb coupling parameter Eq. (1). When Γ > > 1, the dust is strongly coupled and this condition is met in several laboratory dusty plasmas, such as dust “plasma crystals”. When the dust is weakly coupled, the dispersion relation of waves is not affected by the spatial correlation of the dust grains. A dusty plasma is considered as weakly coupled if Γ < < 1 [6].

5.5. Lattice parameter

In the case of Yukawa interaction, additional screening parameter becomes necessary

$$\kappa =\frac{a_{\mathit{ws}}}{{\mathrm{\lambda }}_D}$$E4

It is how much effective and protective the shielding out behavior of a single specie against the external or internal stimuli (voltages) that affect inside the plasmas. a_ws is the Wigner-Seitz radius which can be estimated as a_ws = (4πn/3)^−1/3 [7].

6. The forces on a dust grain in a plasma

The knowledge about various forces acting on dust particles in a plasma is necessary for an understanding of their dynamics and transport.

6.1. Force of gravity, F_g

The dust particle is subject to gravity with a force that is proportional to its mass, but under microgravity condition, it is ignored.

$${\mathit{F}}_g=m\mathit{g}=\frac{4}{3}\pi b^3\rho \mathit{g}$$E5

where g is the local acceleration due to gravity and ρ is the mass density of the particle [8].

6.2. The electric force, F_e

At a location in the plasma having an electric field E, the electric force acting on dust particles of charge q is

$${\mathit{F}}_e=q\mathit{E}$$E6

The electric field is smaller in the bulk of the plasma while it is larger in sheaths next to the plasma wall boundary.

6.3. Neutral drag force, F_n

This force is produced from collisions with the background neutral gas atoms or molecules, and is proportional to the neutral pressure in the vacuum chamber. F_n is given as

$${\mathit{F}}_n={\mathit{Nm}}_n{v^2}_n\mathrm{\pi }{b}^2$$E7

where N is the neutral density, m_n is the mass of the neutral atoms (or molecules), and v_n is the average relative velocity [1].

6.4. Thermophoretic force, F_th

This force will be produced due to the temperature gradient in the neutral gas in the plasma, and it‘s direction is opposite to temperature gradient. This force is given approximately by

$$\begin{array}{l}{\mathit{F}}_{\mathit{th}}=\frac{16\sqrt{\mathrm{\pi }}}{15}\frac{b^2{k}_T}{{\nu }_{T,n}}\nabla T_n\\ \end{array}$$E8

where, v_T,n is the thermal speed of the neutral gas, k_T is the translational part of the thermal conductivity, and T_n is the temperature of the neutrals [9].

7. Transport properties

7.1. Diffusion coefficient

The diffusion coefficient is the proportionality constant between molar fluxes which is ${\mathrm{J}}_A=-D\mathrm{\Delta }{\mathrm{C}}_A$, where $\mathrm{\Delta }$ is known as gradient operator, D is known as diffusion coefficient (m². s⁻¹), and C_A is the concentration (mole/m³). In the complex (dusty) plasma, the mass, momentum, and energy are transported through dust particles.

7.2. Shear viscosity

Shear viscosity is the measure of force between different layers of fluid. It is the dynamical property of a material such as liquid, solids, gas and dusty plasmas. For liquid, it is familiar thickness, for example, honey and water have different viscosity. The ideal fluid has no resistance between layers in shear stress. The superfluid has zero viscosity at very low temperature. All other liquids have positive viscosity and are said to be viscous.

7.3. Thermal conductivity

The most vital property of the dusty plasma liquids (DPLs) is thermal conductivity and it is due to the internal energy of the molecules. The energy transference of the DPLs and its dependence on the applied external field can be checked by the thermal conductivity. The complete analysis of it is important and necessary for the designing and manufacturing of numerous heat flow devices. The applications in thermo-electronic devices e.g., in semiconductor systems, the phonon thermal conductivity has got a special attention. Due to the atoms’ oscillations, phonons of different wavelengths and frequencies are created in solids, and these phonons would disappear when the oscillations stop. Especially, nanotechnology (nanomaterials) requires the accurate calculation of heat transport features [10, 11]. Molecular simulation is recognized as substantial for micro- and nano-scale heat transport phenomena. Furthermore, the new advances require the complete explanations of phase change and the heat, mass transport in micrometer to nanometer scale regimes.

8. Applications

8.1. Dust is a good thing

In the present era, scientists do not consider the dust as an undesirable pollutant; interestingly its positive impacts lead in manufacturing, designing of new devices and direct new developments in material science. In the plasma-chemical mechanism, fine dust particles are also important and having useful properties related to their size and composition. A few of the applications are given below:

1. The efficiency and lifetime of silicon solar cell was increased by the incorporation of amorphous hydrogenated silicon particles (a-Si:H) with the nanocrystalline silicon particles which grow in silane plasmas.

2. Thin films of TiN in an amorphous Si₃N₄ matrix prepared by PECVD (plasma enhanced chemical vapor deposition) have enormously high hardness and elastic modulus. A thin film coating is applied to materials to improve surface properties.

3. Diamond whiskers made up by etching in radio frequency (rf) plasmas improve electron field emission. The reactive ion etching process used in rf plasma devices effectively sharpen the micro tips of diamond [12].

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

At first, it was presumed that the semiconductor surfaces were contaminated during handling of the wafers. To lighten the problem, all fabrication steps were done in clean rooms. Yet even with the best state-of-the art clean rooms, semiconductor wafers showed evidence of contamination. It drove out that the surfaces were being soiled by dust particles generated within the processing plasmas. Complex plasma-enhanced chemical reactions take place within these discharges that produce and grow dust particles. Several experimental devices are used to measure the presence and growth of dust particles such as transmission electron microscope (TEM), scanning electron microscope (SEM), laser light scattering etc. Also, theoretical and computational works are directed to investigate the dust formation, growth, charging, and transport.

9. Computational method

9.1. Simulation technique and parameters

Molecular mechanics dynamic simulation (MMDS) is a simplified approach as compared to other techniques. It allows study of molecular ensembles for thousands of atoms. The MMDS technique works as a core on a simple explanation of force between the individual atoms. Here, HPMD approach is implemented to determine the thermal conductivity of CDPs by applying external perturbation which is modeled by using Yukawa potential model use for the explanation of dust particles interacting with one another. Yukawa potential is used for a system of charged particles. While Green-Kubo relation applies to neutral particles,

$$\varphi \left(\left|\mathbf{r}\right|\right)=\frac{q^2}{4{\mathrm{\pi \epsilon }}_0}\frac{e^{-\left|\mathbf{r}\right|/{\mathrm{\lambda }}_D}}{\left|\mathbf{r}\right|},$$E9

The plasma phase of Yukawa system is representing three dimensionless parameters [7], plasma coupling parameter Γ= (q² /4πε₀) (1/a_wsK_BT), screening strength κ = a_ws /λ_D, and F_e(t) = (F_z) is external perturbation with its normalized value F^* = F_za_ws/J_QZ, where J_QZ is thermal heat energy along z-axis. The inverse of plasma frequency ω_p = (q²/2πε₀ma_ws³)^1/2 characterizes timescale. Simulations are performed for N = 400–14,400 particles in canonical ensemble with PBCs and minimum image convention of Yukawa particles. In our case, most of the simulations are performed with N = 400 particles. The particles are placed in a unit cell with edge length L_x / L_y and the dimensions of square simulation cell are L_xa_ws × L_ya_ws. The equations of motions for N-Yukawa dust particles are integrated through the predictor–corrector algorithm with simulation time step of Δt = 0.001ω_p⁻¹. In our case, the conductivity calculations are reported for a wide range of plasma coupling (1 ≤ Γ ≤ 100) and screening parameters (1 ≤ κ ≤ 4) of 2D Yukawa system at constant normalized external perturbation F^*.

9.2. HNEMD model and thermal conductivity

The Green-Kubo relations (GKRs) are the mathematical terms for transport coefficients in the form of time integral correlation functions. GKR is for hydrodynamic transport coefficient of neutral particles. This formula gives linear response expression for thermal conductivity. It enables our calculations using a time-series record of motion of individual dust particles. For thermal transport coefficient, it is a time integral of the correlation function of the microscopic flux of heat energy and where the required input includes time series for position and velocity of a dust particle.

$$\mathrm{\lambda }=\frac{1}{2K_{\mathrm{B}}{\mathit{AT}}^2}\underset{0}{\overset{\infty }{\int }}⟨{\mathbf{J}}_Q\left(t\right).{\mathbf{J}}_Q\left(0\right)⟩\mathit{dt},$$E10

where A represents the area, T denotes the absolute temperature, K_B is Boltzmann’s constant. The relation of microscopic heat energy J_Q is

$${\mathbf{J}}_{Q}A={\displaystyle \sum _{i=1}^N}\left[{\mathrm{E}}_i-\frac{1}{2}\sum _{i=\mathrm{j}}\left(\mathit{r}i-\mathit{r}j\right).\mathit{F}\mathit{ij}\right].\frac{{\mathit{p}}_i}{m},$$E11

In this equation, r_ij = r_i– r_j is the position vector and 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 particle i is E_i = p_i / 2 m + ½ Σɸ_ij, for i ≠ j, where ɸ_ij is the Yukawa pair potential given in Eq. (9) between particle i and j. According to linear response theory (LRT), the perturbed equations of motion, given by Evans-Gillan [7] is

$$\begin{array}{l}{\dot{\mathit{r}}}_i=\frac{{\mathit{p}}_i}{m}\\ {\dot{\mathit{p}}}_i={\mathit{F}}_i+{\mathbf{D}}_i\left({\mathit{r}}_i,{\mathit{p}}_i\right).{\mathit{F}}_e\left(t\right)\\ \alpha =-\frac{{\displaystyle {\displaystyle \sum _{i=1}^N}}\left[{\mathit{F}}_i+{\mathbf{D}}_i\left({\mathit{r}}_i,{\mathit{p}}_i\right).{\mathit{F}}_e\left(t\right)\right].{\mathit{p}}_i}{{\displaystyle {\displaystyle \sum _{i=1}^N}}{p}_i^2/m_i}\end{array}$$E12

The tensorial phase space distribution function D_i (r_i, p_i) describes the coupling of the system. In the generalization of LRT to a system moving according to non-Hamiltonian dynamics, the response of D_i (r_i, p_i) [13] is

$$⟨A\left(t\right)⟩-⟨A\left(0\right)⟩=\beta \underset{0}{\overset{t}{\int }}⟨A\left(t\right){\dot{H}}_0\left(t^{\prime }\right)⟩\mathit{dt},$$E13

where Ḣ₀ is the time derivative of the total energy with respect to field-dependent equation of motion [7] and average brackets denote the statistical average and β = 1/ K_BT.

$$⟨{\mathbf{J}}_Q\left(t\right)⟩=\mathit{\beta A}\underset{0}{\overset{t}{\int }}\mathit{dt}⟨{\mathbf{J}}_{Q_Z}\left(t^{\prime }\right).{\mathbf{J}}_{Q_Z}\left(0\right)⟩.{\mathit{F}}_e\left(t\right),$$E14

When external force is selected parallel to the z-axis F_e(t) = δ (0, F_z), δ is Dirac delta function. Due to this function, the response of heat energy current is proportional to autocorrelation function itself rather than time integral of this function [14]. The reduced thermal conductivity has the following form:

$$\mathrm{\lambda }=\frac{1}{{\mathit{ATF}}_Z}\underset{0}{\overset{\infty }{\int }}⟨{\mathbf{J}}_{Q_Z}\left(t^{\prime }\right)\mathbf{.}{\mathbf{J}}_{Q_Z}\left(0\right)⟩\mathit{dt},$$E15

Eq. (15) is the basic formula for evaluation of autocorrelation function of heat energy current by a perturbation method. Here it is important to discuss some other factors that are in association with the thermal conductivity i.e., Ewald sum. It is used to measure force, Yukawa potential energy, heat energy current (GKR). In this scheme, original interaction potential is divided into two parts: the long-range part that converges quickly in reciprocal space, a short-range interaction that converges quickly in the real-space part of Ewald-Yukawa potential [15].

10. HNEMD results and discussion

In this section, the thermal conductivity calculations are obtained through homogenous perturbed MD (HPMD) simulations, using Eq. (15), for 2D complex dusty plasma systems. The thermal conductivity is compared here with appropriate frequency normalization in the limit of a suitable equilibrium low value of normalized external perturbation, for an absolute range of plasma coupling (Γ ≥ 1) and screening strength (κ ≥ 1). For 2D case, the thermal conductivity of complex dusty plasmas may be represented as λ₀ = λ/nmω_pa_ws² (normalized by plasma frequency) or λ^* = λ/nmω_Ea_ws² (normalized by plasma frequency). This improved HPMD approach to 2D strongly coupled plasmas enables it possible to compute all the possible range of plasma states (Γ, κ) at a constant value of normalized perturbation F^* = (F_za_ws/J_QZ). In our case, the possible low value of external perturbation is F^* = 0.02 at which 2D complex plasma system gives equilibrium thermal conductivity for all plasma state points. Before the external perturbation F^* is switched on, the system is equilibrated using the Gaussian thermostat which generates the canonical ensemble given in Eq. (12). In practice, it is necessary for the MD system to be thermostated for the removal of additional heat that is generated due to work done by the external perturbation F* [16]. The results obtained through present HPMD approach are shown in Figures 1–4, where we have traced the plasma thermal conductivity through a computation of usual Yukawa particles in 2D within the strongly coupled regime for different screening parameters of κ = 1, 2, 3, and 4, respectively.

Figure 1.

Comparison of results obtained from Yukawa thermal conductivity λ₀ (normalized by ω_p) as a function of plasma coupling Γ (1 ≤ Γ ≤ 50) (system temperature) for SCCDPs at 𝜅 = 1.4.

Figure 2.

Comparison of results obtained from Yukawa thermal conductivity λ₀ (normalized by ω_p) as a function of plasma coupling Γ (1 ≤ Γ ≤ 50) (system temperature) for SCCDPs at 𝜅 = 2.

Figure 3.

Comparison of results obtained from Yukawa thermal conductivity λ₀ (normalized by ω_p) as a function of plasma coupling Γ (1 ≤ Γ ≤ 50) (system temperature) for SCCDPs at 𝜅 = 3.

Figure 4.

Comparison of results obtained from Yukawa thermal conductivity λ₀ (normalized by ω_p) as a function of plasma coupling Γ (1 ≤ Γ ≤ 50) (system temperature) for SCCDPs at 𝜅 = 4.

Figures 1 and 2 show the thermal conductivity for the cases of κ = 1.4 and 2, respectively. For both cases, our simulations cover the appropriate range of Coulomb coupling parameter i.e., from the nearly liquid state to strongly coupled states. It is observed that our investigation of λ₀ at low value of Γ (= 10) is definitely higher than that of GKR-EMD estimations of Khrustalyov and Vaulina [17] and NEMD of Hou and Piel [18] but for κ = 1.4 results are slightly higher than HNMED (N = 1024) simulations of Shahzad and He [19]. It is noted that our result for the low value of Γ shows that particle–particle interactions are very weak and particles have maximum kinetic energy and the effectiveness of screening parameter is large. At intermediate to higher Γ (= 20, 50), the present results lie closer to earlier 2D NEMD simulations [18] and HNMED (N = 4096) computations [19] but slightly less than 2D dissipative Yukawa GKR-EMD numerical results [17]. For both cases, it can be seen that the presented λ₀ is well matched with earlier 2D numerical estimations [19] at intermediate Γ (= 20). It is significant to note that a constant λ₀ is observed at intermediate to higher plasma coupling Γ at constant external perturbation F^* = 0.02; however, it is observed that a very slightly decreasing behavior is observed at higher Γ, contrary to earlier simulations of Shahzad and He [19]. But it is examined that a constant λ₀ is found at intermediate to higher Γ at constant F*.

Two further set of simulations is plotted to illustrate the plasma λ₀ behaviors of the simulated complex dusty plasmas at a higher value of screening. For this case, Figures 3 and 4 show the normalized λ₀ computed by the HPMD approach for N = 400 at κ = 3 and 4 and a sequence of different simulations are performed. It is characterized by these figures that the present results lie close to the earlier 2D NEMD results of Hou and Piel [18] at intermediate to higher Γ (= 20, 50). For the κ = 3 at a lower value of Γ, our simulation result is slightly higher than earlier HNMED simulation result, however, for both cases at a lower value of Γ, it is definitely higher than earlier numerical results of NEMD, GKR-EMD.

11. Conclusions

The improved Evan-Gillan HPMD method is used to investigate the thermal conductivity of the 2D strongly coupled complex Yukawa liquid for a suitable range of plasma parameters of screening lengths κ (=1, 4) and Coulomb couplings Γ (=1, 100). Nonequilibrium molecular dynamics method uses the thermal response of heat energy current to calculate the preliminary results of plasma thermal conductivity. The presented method is better than earlier HNEMD and NEMD methods because the very small value of external perturbation (F^* = 0.02) is only imposed on several individual particles each time step. It is concluded that the present approach for evaluating the thermal conductivity from homogenous PMD method yields consistent results and this method is quite accurate and much faster than the previous EMD and NEMD methods. For future work, the system size (N) and external perturbation strength (F^*) can be varied to examine how effectively this improved HPMD algorithm calculates the thermal conductivities of Yukawa and other Coulomb systems. It is suggested that the presented HPMD technique based on Ewald summation described here can be used to explore the ionic and dipolar materials.

Acknowledgments

The authors thank the National Advanced Computing Centre of National Centre for Physics (NCP), Pakistan and National High-Performance Computing Center (NHPCC) of Xian Jiaotong University, P.R. China for allocating computer time to test and run our MD code.

Abbreviations

SCCDPs	strongly coupled complex dusty plasmas
EMD	equilibrium molecular dynamics
HNEMD	homogeneous nonequilibrium molecular dynamics
HPMD	homogeneous perturbed molecular dynamics
MMDS	molecular mechanics dynamic simulation
DPLs	dusty plasma liquids
LRT	linear response theory
GKR	Green-Kubo relation
CDPs	complex dusty plasma
PBCS	periodic boundary conditions
Γ	Coulomb coupling
κ	screening strength
λD	Debye length
KB	Boltzmann constant