In recent years, the use of new materials for different specific applications, it is vital to understand and characterize its transport properties such as thermal transport, mass transport, and electrical transport. In the advancement of large and fast massively parallel computers, it is now possible to devise new molecular modeling methods that can reliably compute transport properties of complex materials from bottom to the top. Technological development in different areas such as electronics, automobiles and nuclear energy demands better thermal management. Conventional methods of measuring heat properties are reaching their limits. There is an immediate need to look for innovative means to achieve enhanced heat capacity and low thermal conductivity of conventional heat transfer fluid is proving to be a serious limitation. A comprehensive and systematic effort is necessary to incorporate the effective thermal conductivity of complex materials and their limitations. The real structures and geometries of multi-phase materials are so vast and vivid that one cannot use a single model to estimate the effective thermal conductivity of complex fluid materials in the whole range due to their inherent limitations. In this scenario, complex fluid materials have emerged as an attractive solution to meet these challenges.^[1]

In addition to further characteristics of complex fluid materials, the composition and dimensions of these fluid materials help in the investigation of microscopic level properties, which are difficult to measure even with state of the art experimental techniques. Such techniques are also required for studying new thermo-physical phenomena in the micro to a nanoscale device made of novel nanostructures, such as carbon electronics and nano level fluidic flow.^[2–5] Similarly, size, and dimension effects on vibrational modes for thermal conductivity in a complex crystal are extremely important especially for the variety of electronic and energy conversion technologies. The boundary scattering effects are strong for many vibrational modes because of strong anisotropy in their non-propagating nature. This non-diffusive (non-propagating) nature can be suppressed by the adequate composition of nanostructures inside the fluids materials. So, the lattice thermal conductivities get important for fluids materials, dependent on their structure, heat capacities and the equation of state.^[6] Thermal conductivities are carried out by atomic vibrations called phonons (quasi-particle traveling at the speed of sound). Heat propagation and conductivity are directly related to the time a phonon travels in a material before it collides with another phonon or defect. This characteristic time is called phonon lifetime. Shortening phonon lifetimes achieve low thermal conductivity, which is important for thermo-electronic materials. Phonon lifetime is one of the key parameters for quantifying the thermal conductivity, but assessing and measuring it is extremely challenging both experimentally and theoretically. State of the art simulation techniques is required to overcome it. Traditional computer simulations did not lead to a dramatic increase in thermal conductivity as found in experiments.^[7]

The atomic and molecular levels study of fluid dynamics is still an open challenge from the analytical and experimental point of view. The molecular dynamics (MD) simulation is the powerful tool for the investigation of transport properties of fluid dynamics at atomic or molecular scale and gives the understandings about the complex systems that cannot be directly accessible by experiments.^[8–10] MD simulation is the favorable technique over the other computational techniques for the dynamical study of materials and complex systems.^[11] The thermal properties, structures, and behavior of complex and large systems may be explored by using faster and powerful computer simulation modeling tools. Generally, two MD methods are in considerations for the estimation of thermal conductivity, namely, equilibrium MD (EMD), non-equilibrium MD (NEMD). The former method employs the Green-Kubo relations (GKRs) to calculate time correlation function of microscopic heat flux and the thermal conductivity. The method based on GKRs requires a long time to run the simulation and its computational cost limits to the systems having only a few hundred atoms. The lateral method is also a popular way to simulate the complex mechanisms, an example of it is inhomogeneous NEMD (InHNEMD) techniques, which mimic the experiment by imposing temperature differences,^[11–12] heat current,^[13–14] transient heat impulsion^[15] across the system. This technique has faster convergence than EMD method. However, both techniques have complexities in terms finite size effects, large temperature gradients, spatially inhomogeneity. The Evan’s homogenous NEMD (HNEMD) required minimum sizes and produced small statistical errors than previously discussed techniques.^[16] The HNEMD has no such complexities as discussed in above mentioned NEMD methods. The reason for preferring this scheme is that if physical walls are replaced by periodic boundary conditions (PBCs), all particles perceive similar treatment. The link between EMD and NEMD methods is due to fluctuation-dissipation theorem,^[17] which represents the relation between linear response to external perturbations and equilibrium time correlation function of the fluxes of the system. The HNEMD techniques have already been checked with emphasis on the application on heat transport problems of different fluids, thermal transport coefficients, and the autocorrelation function of heat current.^[14]

The area of complex dusty plasmas has rapidly grown over the past two decades to become the major area of study in the field of plasma physics. The results of strongly coupled complex dusty plasmas (SCCDPs) have been calculated for three-dimensional (3D) thermal conductivity (λ) by Salin and Caillol,^[18] Faussurier and Murillo,^[19] Donko et al.^[20] Donko and Hartmann^[11] and presented authors of Shahzad and He,^[21–22] as well as two-dimensional (2D) estimations, have been explored by Hou and Piel,^[4] Khrustalyov and Vaulina^[5] and presented authors of Shahzad and He.^[23] Moreover, for the 2D systems, a t⁻¹ long-time tail exists for auto correlation function (ACF) as compared to t^−3/2 short-time decay for the 3D systems. Presently, the external small perturbation is applied for the 2D dusty plasma systems at the lower-intermediately and high couplings (Γ). This small perturbation can cause small deviation from the initial state and the current ACF decays faster than t⁻¹ for large simulation time step dt, and it shows a definite value of heat energy flux J(t).^[10,20] The current ACF oscillation, retardation of perturbation, and damping of external field generate smooth decay in the vicinity of 2D dusty plasmas. Consequently, integral converges when attempting to compute the transport coefficients by using Green-Kubo relation. The main objective of current work is to investigate the preliminary normalized (λ₀) of 2D complex fluid material in strongly coupled complex regime through an improved Evan-Gillan HNEMD algorithm at constant external perturbation. The Gillan and Dixon^[24] have also used this modified approach for LJ liquids to measure the autocorrelation function of microscopic heat current and thermal coefficients with weak external perturbations. This modified algorithm has already been used for transport coefficients of one component plasma (OCCP),^[25] ionic liquids and for the investigation of simple fluids,^[13–14] rheological issues of Yukawa liquids,^[11] and semiconductor systems.^[2] Therefore, this method is considered as a best computational tool in the limit of zero applied external perturbation. The extensive HNEMD simulations are performed to study the performance of algorithms and to compare the results obtained from EMD and NEMD simulations for a wide range of plasma parameters (Γ, κ) than those used for formerly used for Yukawa liquids.

In this section, the thermal conductivity calculations are obtained through homogenous perturbed MD (HPMD) simulations, using Eq. (10), 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 (Γ ≥ 10) and screening strength (κ ≥ 1). For 2D case, the thermal conductivity of complex dusty plasmas may be represented as (normalized by plasma frequency) and or (normalized by plasma frequency). These types of normalizations have been used usually for the earlier studies of OCCP^[25] and CDPLs^[18–19] for estimating thermal conductivity. Especially for the 3D strongly coupled system the Einstein frequency decreases by increasing screening parameter.^[22,28] This improved HPMD approach to 2D strongly coupled plasmas enables it possible to compute all the possible range of plasma states (Γ, κ) at constant value of normalized perturbation P* (= P_za_ws/J_QZ). For results reported here, we have checked and varied the following parameters including system size (N), normalized perturbation (P*), thermostat (α), simulations total run time, simulation step size (dt), and Debye screening (κ), Coulomb coupling (system temperature ≡ 1/Γ) for the investigation of plasma thermal conductivity. Different sequences of HPMD simulations are performed for various suitable low values of normalized external perturbation in order to find appropriate value of P*. In our case, the possible low value of external perturbation is P* = 0.02 at which 2D complex plasma system gives equilibrium thermal conductivity for all plasma state points. It is interesting and significant here that this normalized steady state low value of P* = 0.02 is very small as compared to earlier known value of external force in Ref. [23] This low value reflects more appropriate and acceptable results using presented HPMD technique than earlier used HNEMD technique. The λ₀ results obtained through HPMD computer simulation are checked for the universal temperature scaling law at this reduced steady-state value of P* = 0.02:

We now turn our attention to the main results obtained through the HPMD simulations. In our case, before the external perturbation P* is switched on, the system is equilibrated using the Gaussian thermostat, which generates the canonical ensemble given in Refs. [22–23]. 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 P*.^[3,23] Presently, for a possible low value of the external perturbation strength of P* = 0.02 (steady state value) is to be chosen for the estimation of equilibrium thermal conductivity at all plasma states of Γ (≡ 10, 100) and κ (≡ 1, 3). The results obtained through present HPMD approach are shown in Figs. 1–3, 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, and 3 respectively. These figures show our key results along with the earlier numerical estimations taken from 2D GKR-EMD of dissipative Yukawa systems of Khrustalyov and Vaulina^[5] as well as the 2D NEMD results of Hou and Piel^[4] and the previous measurements taken from the 2D homogenous NEMD computations of Shahzad and He^[23] nearly at the same data points. Our HPMD data are in practically good agreement with the previous numerical computations based on different methods that yield better measurements for plasma thermal conductivity. It is observed from figures that the measured thermal conductivity has lower values as compared to earlier estimated values at nearly same plasma state points. The presented simulation results λ₀(Γ) are performed for N = 400 particles and a sequence of four different computations are taken into account at constant perturbation of P* = 0.02 for each κ = 1, 2, and 3, respectively. Our numerical data of λ₀ (Γ) are in nearly reasonable agreement with earlier numerical data of GKR-EMD, NEMD and HNEMD investigation.^[4–5,23]

Fig. 1

Figure Option ViewDownloadNew Window

Fig. 1 Comparison of results obtained from Yukawa thermal conductivity λ0 (normalized by ωp) as a function of plasma coupling Γ (system temperature) for SCCDPs at κ = 1. Our 2D HPMD simulation results: present data (for N = 400 particles) and simulation results for the 2D HNEMD obtained by the Shahzad and He,[23] 2D NEMD (Brownian dynamics) results of Hou and Piel,[4] GKR-EMD of Khrustalyov and Vaulina at scaling factors of ζ = 1, 0.25 and ∞.[5]

Figures 1 and 2 show the thermal conductivity, normalized by the plasma frequency (ω_p), as a function of Coulomb coupling (system temperature = 1/Γ) for the cases of κ = 1 and 2, respectively. For both cases, our simulations covering the appropriate range of Coulomb coupling parameter i.e. from the nearly liquid state to strongly coupled states, depending on different κ values. The presented simulation data are generally in fair agreement at nearly same plasma parameters and figures show overall the same trends as in the earlier numerical methods of 2D Yukawa liquids.^[4–5,23] It is observed that our investigation of λ₀ at low value of Γ (= 10) is definitely higher than that of NEMD of Hou and Piel^[4] and GKR-EMD estimations of Khrustalyov and Vaulina^[5] but slightly higher that HNMED (N = 1024) simulations Shahzad and He.^[23] It is noted that our result for 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, and 100), the present results lie closer to earlier 2D NEMD simulations^[4] and HNMED (N = 4096) computations^[23] but slightly less than 2D dissipative Yukawa GKR-EMD numerical results.^[5] For both cases, it can be seen that the presented λ₀ is well matched with earlier 2D numerical estimations^[23] at intermediate Γ(= 20). It is significant to note that a constant λ₀ is observed at intermediate to higher plasma coupling Γ at constant external perturbation P* = 0.02, however, it is observed that a very slightly decreasing behavior is observed at higher Γ, contrary to earlier simulations of Shahzad and He.^[23] But it is examined that a constant λ₀ is found at intermediate to higher Γ at constant P*, confirming earlier numerical results.^[4–5,23] It is interesting to note here that the existence of λ₀ is present for low-intermediate to higher Γ with an increase in κ and remains within a satisfactory limited statistical uncertainty, confirming previous simulation results.^[23] In our simulation, the presented plasma conductivity for lower to intermediate Γ shows the existence of λ₀ and it is a clear contradiction with the earlier simulation results of Donko et al.^[29] where the λ₀ was not found at lower Γ.

Fig. 2

	Figure Option ViewDownloadNew Window
	Fig. 2 Comparison of results obtained from Yukawa thermal conductivity λ0 (normalized by ωp) as a function of plasma coupling Γ (system temperature) for SCCDPs at κ = 2. For details, see the caption of Fig. 1.

Fig. 3

	Figure Option ViewDownloadNew Window
	Fig. 3 Comparison of results obtained from Yukawa thermal conductivity λ0 (normalized by ωp) as a function of plasma coupling Γ (system temperature) for SCCDPs at κ = 3. For details, see the caption of Fig. 1.

One further set of simulations is plotted to illustrate the plasma λ₀ behaviors of the simulated complex dusty plasmas at higher value of screening. For this case, Fig. 3 shows the normalized λ₀ computed by the HPMD approach for N = 400 at κ = 3 and a sequence of different simulations is performed. This figure shows that our results are satisfactory agreement with various simulation data sets and the uncertainties inherent to the different earlier approaches are comparable. It is depicted from this figure that the present results lie close to the earlier 2D NEMD results of Hou and Piel^[4] at intermediate to higher Γ(= 20, 100). At lower value of Γ, our simulation result is slightly higher than earlier HNMED simulation result, however, it is definitely higher than earlier numerical results of NEMD, GKR-EMD. Moreover, it is examined that the presence of normalized plasma λ₀ at all plasma state points and it is observed that plasma λ₀ found to be constant, as expected in earlier simulation results.^[4] Moreover, it is noted that the measured data of plasma λ₀ estimations at P* = 0.02, where plasma λ₀ has equilibrium values and independent of P*, are within limited statistical uncertainties. The overall HPMD simulation data obtained with lower system size (N = 400) are revealed to be well matched within statistical limits of errors at lower, intermediate and higher Γ states, however, some data points are deviate at the lower Γ states. This deviation grows up suddenly at lower Γ states and intermediate screening κ = 2. This deviation of our numerical result from previous data point is still acceptable, for all cases. It is demonstrated from all figures with comparisons of earlier results that the presented results through HPMD approach with lower N are more accurate and acceptable.

We determine the universal behavior of complex dusty plasma in which normalized conductivity λ* follows a temperature scaling law. Figure 4 shows the variation of normalized plasma verses various normalized temperatures T* (= T/T_m). It is observed that these measured results are in good agreement to the former reported results for 2D Yukawa liquids.^[11] In our case, T* is plotted along horizontal axis and λ* is plotted along vertical axis as shown in figure. This figure displays the variation of normalized (by Einstein ω_E) λ* for different normalized temperature T* at κ = 1, 2, and 3. The bold line, shown in Fig. 4, is obtained by fitting curve of the functional form (temperature scaling law) given in Eq. (11) reported in Refs. [19, 23] with the coefficients: A = 0.023 02, B = −1.494 22 and C = 0.693 73. These obtained fitting coefficients (A, B and C) for the dimensionless plasma thermal conductivity given in Eq. (11), λ* = AT* + B/T* + C, are measured from presented HPMD simulation data display in Figs. 1–3. It is observed that there is dispersion of obtained data of normalized λ* shown in Fig. 4. The scattering of these data from bold line suggested one possible reason that this may be happen due to high negative value of coefficient B in the functional fit of Eq. (11) in comparison to the previous EMD, HNMED measurements. It is noted that the bold fitting line is nearly exact fitting, confirming earlier numerical results. Moreover, Einstein frequency (ω_E) is much more important than plasma frequency (ω_p) because the distribution of data along solid line explains more accurately the physical significance of dusty plasma thermal conductivity. It is observed from Fig. 4 that the λ* of dusty plasma is close to functional form for low value of κ and T*. For higher values of κ, at intermediate values of T*, λ* shows less dependence on these two variables and it is little far from functional form (temperature scaling). But at higher T* the conductivity of dusty plasmas is close to functional form, for κ = 3 in 2D case. It is concluded that present results show the less growing behavior of dusty plasmas λ* with the increase of normalized temperature and screening. The plotted functional form demonstrates the correct universal behavior at three different values of screening on the extensive range of T* for the reduced force field strength P* = 0.02. The result measured by plasmas λ*(T*), employing the HPMD technique, give empirical fitting and the plot shows better fit compared to the prior results.

Fig. 4

	Figure Option ViewDownloadNew Window
	Fig. 4 Variation of normalized plasma conductivity (λ*) by Einstein frequency (ωE) with normalized temperature (T*) for complex dusty plasmas system at different κ = 1, 2, and 3. The bold line is computed by employing simple functional form of: λ* = AT* + B/T* + C, representing the universal temperature law for the 2D complex dusty plasmas.[23]

We have estimated thermal conductivity of the 2D strongly coupled complex Yukawa liquid using improved Evan-Gillan HPMD approach for suitable range of plasma parameters of screening lengths κ (=1, 3) and Coulomb couplings Γ (=10, 100). Nonequilibrium molecular dynamics method uses the thermal response of heat energy current to calculate the preliminary results of plasma thermal conductivity. Our presented method is better than earlier HNEMD and NEMD methods because the very small value of external perturbation (P* = 0.02) is only imposed on several individual particles each time step. We have shown that normalized plasma λ* as the function of normalized temperature T* follows simple temperature scaling law. 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 (P*) can be varied to examine how effectively this improved HPMD algorithm calculates the thermal conductivities of Yukawa and other Coulomb systems.