Open access

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

Written By

Aamir Shahzad, Syed Irfan Haider, Maogang He and Yan Feng

Submitted: 28 January 2018 Published: 05 September 2018

DOI: 10.5772/intechopen.75367

From the Edited Volume

Impact of Thermal Conductivity on Energy Technologies

Edited by Aamir Shahzad

Chapter metrics overview

1,238 Chapter Downloads

View Full Metrics

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., SiO2 or Al2O3, 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:


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].


where λDe = √KTe /4πnee2 and λDi = √KTi /nie2 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


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


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. aws is the Wigner-Seitz radius which can be estimated as aws = (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, Fg

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


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

6.2. The electric force, Fe

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


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, Fn

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. Fn is given as


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

6.4. Thermophoretic force, Fth

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


where, vT,n is the thermal speed of the neutral gas, kT is the translational part of the thermal conductivity, and Tn 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 JA=DΔCA, where Δ is known as gradient operator, D is known as diffusion coefficient (m2. s−1), and CA is the concentration (mole/m3). 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 Si3N4 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,


The plasma phase of Yukawa system is representing three dimensionless parameters [7], plasma coupling parameter Γ= (q2 /4πε0) (1/awsKBT), screening strength κ = awsD, and Fe(t) = (Fz) is external perturbation with its normalized value F* = Fzaws/JQZ, where JQZ is thermal heat energy along z-axis. The inverse of plasma frequency ωp = (q2/2πε0maws3)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 Lx / Ly and the dimensions of square simulation cell are Lxaws × Lyaws. 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−1. 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.


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


In this equation, rij = rirj is the position vector and Fij is the force of interaction on particle i due to j and pi represents the momentum vector of the ith particle. The energy Ei of particle i is Ei = pi / 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


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


where 0 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/ KBT.


When external force is selected parallel to the z-axis Fe(t) = δ (0, Fz), δ 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:


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 λ0 = λ/nmωpaws2 (normalized by plasma frequency) or λ* = λ/nmωEaws2 (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* = (Fzaws/JQZ). 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 14, 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 λ0 (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 λ0 (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 λ0 (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 λ0 (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 λ0 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 λ0 is well matched with earlier 2D numerical estimations [19] at intermediate Γ (= 20). It is significant to note that a constant λ0 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 λ0 is found at intermediate to higher Γ at constant F*.

Two further set of simulations is plotted to illustrate the plasma λ0 behaviors of the simulated complex dusty plasmas at a higher value of screening. For this case, Figures 3 and 4 show the normalized λ0 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.


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.


SCCDPsstrongly coupled complex dusty plasmas
EMDequilibrium molecular dynamics
HNEMDhomogeneous nonequilibrium molecular dynamics
HPMDhomogeneous perturbed molecular dynamics
MMDSmolecular mechanics dynamic simulation
DPLsdusty plasma liquids
LRTlinear response theory
GKRGreen-Kubo relation
CDPscomplex dusty plasma
PBCSperiodic boundary conditions
ΓCoulomb coupling
κscreening strength
λDDebye length
KBBoltzmann constant


  1. 1. Shukla PK, Mamun AA. Series in Plasma Physics: Introduction to Dusty Plasma Physics. Bristol: Taylor & Francis group, CRC Press; 2001; ISBN: 9780750306539 – CAT#IP39
  2. 2. Merlino RL. Dusty Plasmas and Applications in Space and Industry. Plasma Physics Applied. 2006; ISBN: 81-7895-230-0
  3. 3. Bellan PM. Fundamentals of Plasma Physics. 1st ed. UK: Cambridge University Press; 2008. ISBN-13: 9780521528009
  4. 4. Shukla PK. Dusty Plasmas: Physics, Chemistry and Technological Impacts in Plasma Processing. In: Bouchoule A. Wiley, New York. 2000; ISBN: 0471973866
  5. 5. Chen FF. Introduction to Plasma Physics and Controlled Fusion. 2nd ed. New York: Springer-Verlag. 2006; 200 p. ISBN: 9780521825689
  6. 6. Thomas H, Morfill GE, Demmel V, Goree J, Feuerbacher B, Möhlmann D. Plasma crystal: Coulomb crystallization in a dusty plasma. Physical Review Letters. 1994;73(5):652. DOI: 0031-9007/94
  7. 7. Shahzad A, He M-G. Thermal conductivity calculation of complex (dusty) plasmas. Physics of Plasmas. 2012;19(8):083707. DOI: 10.1063/1.4748526
  8. 8. Kikuchi H. Electrohydrodynamics in Dusty and Dirty Plasmas: Gravito-Electrodynamics and EHD. Dordrecht; Boston: Kluwer Academic Publishers; 2001. 207 p. ISBN: 0792368223
  9. 9. Peratt AL. Physics of the Plasma Universe, Appendix C. Dusty and Grain Plasmas. New York: Springer; 1992. ISBN: 0-387-97575-6
  10. 10. Boulos MI, Fauchais P, Pfender E. Thermal Plasmas: Fundamentals and Applications. New York: Springer Science & Business Media; 2013. DOI: 10.1007/978-1-1337-1
  11. 11. Shahzad A, He M-G. Thermal Conductivity and Non-Newtonian Behavior of Complex Plasma Liquids, A Chapter from the Book of Thermoelectrics for Power Generation—A Look at Trends in the Technology. Ri jeka, Croatia: InTech; 2016; 305 p. DOI: 10.5772/65563
  12. 12. Vladimir EF, Gregor EM. Complex and Dusty Plasmas: From Laboratory to Space. Bosa Roca, Uni ted States: Taylor & Francis Inc, CRC Press Inc.; 2010; ISBN-10: 1420083112/ISBN-13: 9781420083118
  13. 13. Faussurier G, Murillo MS. Gibbs-Bogolyubov inequality and transport properties for strongly coupled Yukawa fluids. Physical Review E. 2003;67(4):046404. DOI: 10.1103/PhysRevE.85.046405
  14. 14. Gillan MJ, Dixon M. The calculation of thermal conductivities by perturbed molecular dynamics simulation. Journal of Physics C: Solid State Physics. 1983;16(5):869. DOI: 0022-3719/83/050869
  15. 15. Mazars M. Ewald sums for Yukawa potentials in quasi-two-dimensional systems. The Journal of Chemical Physics. 2007;126(5):056101. DOI: 10.1063/1.2431371
  16. 16. Mandadapu KK, Jones RE, Papadopoulos P. A homogeneous nonequilibrium molecular dynamics method for calculating thermal conductivity with a three-body potential. The Journal of Chemical Physics. 2009;130(20):204106. DOI: 10.1063/1.3141982
  17. 17. Khrustalyov YV, Vaulina OS. Numerical simulations of thermal conductivity in dissipative two-dimensional Yukawa systems. Physical Review E. 2012;85(4):046405. DOI: 10.1103/PhysRevE.85.046405
  18. 18. Hou LJ, Piel A. Heat conduction in 2D strongly coupled dusty plasmas. Journal of Physics A: Mathematical and Theoretical. 2009;42(21):214025. DOI: 10.1088/1751-8113
  19. 19. Shahzad A, He MG. Numerical experiment of thermal conductivity in two-dimensional Yukawa liquids. Physics of Plasmas. 2015;22(12):123707. DOI: 10.1063/1.4938275

Written By

Aamir Shahzad, Syed Irfan Haider, Maogang He and Yan Feng

Submitted: 28 January 2018 Published: 05 September 2018