Thermal Conductivity of Complex Plasmas Using Novel Evan-Gillan Approach

Aamir Shahzad∗Syed Irfan HaiderMuhammad KashifMuhammad Shahzad ShifaTariq Munir and Mao-Gang He

1Molecular Modeling and Simulation Laboratory,Department of Physics,Government College University Faisalabad(GCUF),Allama Iqbal Road,Faisalabad 38040,Pakistan

2Key Laboratory of Thermo-Fluid Science and Engineering,Ministry of Education(MOE),Xi’an Jiaotong University,Xi’an 710049,China

1 Introduction

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 e ff ort 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 nonpropagating 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 fluctuationdissipation 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]Donkoet 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,at−1long-time tail exists for auto correlation function(ACF)as compared tot−3/2short-time decay for the 3D systems.Presently,the external small perturbation is applied for the 2D dusty plasma systems at the lowerintermediately and high couplings(Γ).This small perturbation can cause small deviation from the initial state and the current ACF decays faster thant−1for large simulation time step dt,and it shows a definite value of heat energy fluxJ(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(λ0)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.

2 Computational Method

2.1 Simulation Technique and Parameters

Several simulation techniques have unnecessary cost,then simplify approach could be a select e.g.molecular mechanics dynamic simulation(MMDS).It allows study molecular ensembles for thousands of atoms.The MMDS technique works as a core on a simple explanation of force between the individual atoms.Here,HNEMD approach is implemented to determine the thermal conductivity of CDPLs by applying external perturbation,which is modeled by using Yukawa potential model use for the explanation of dust particles interact with one another.Yukawa potential is used for a system of charged particles.While

2.2 HNEMD Model and Thermal Conductivity

The GRKs are the mathematical terms for transport coefficients in the form of time integral correlation functions.GKR is for hydrodynamic transport coefficient of neutral particles.It has also been used for OCCP[25]and SCCDPs.[18,27]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 it required input include timeseries for position and velocity of a dust particle.

whereArepresents the area,Tdenotes the absolute temperature,kBis Boltzmann’s constant.The relation of microscopic heat energyJQis

In this equationrij=ri−rjis the position vector andFijis the force of interaction on particleidue tojandpirepresents the momentum vector of thei-th particle.The energyEiof particleiisEi=Pi/2m+1/2Σφij,fori/=j,whereφijis the Yukawa pair potential given in Eq.(1)between particleiandj.According to,linear response theory(LRT)the perturbed equations of motion,given by Evans-Gillan in Ref.[22]define interparticle force acting on the particleiand the tensorial phase space distribution functionDi(ri,pi)describes the coupling of the system.According to non-Hamiltonian dynamics,LRT describes thisDi(ri,pi)as arbitrary phase space dynamical variable for a system moving underFe(t)[22−23]and is calculated as

whereH0is the time derivative of the total energy with respect to field dependent equation of motion and in Ref.[22]and average brackets denote the statistical average andβ=1/kBT.

The external force does a mechanical work on the system and it disturbs its equilibrium position,therefore,the Gaussian thermostat is applied in the dynamics of the system to maintain the equilibrium of the system.The dynamics of the system satisfy the condition of adiabatic incompressibility of phase space and Eq.(4)is only valid forBi(ri,pi)=0 and it is given in Ref.[22]ifDi(ri,pi)is taken as

Then,Eq.(4)withJQA=A(t)is simply related to Green-Kubo formula given in Eq.(2).It is an assumption to our system that force is sufficiently weak and the system remains homogeneous and compatible with PBCs by taking momentum derivative sum equals zero.Therefore,the response in heat energy flux is

The above Eq.(10)is the basic formula for evaluation of autocorrelation function of heat energy current by a perturbation method.The efficiency of the above formula depends on the extensive range of Yukawa plasma parameters.It is important to know the perturbation has Dirac delta function,therefore,the response of heat energy current is proportional to autocorrelation function itself rather than time integral of this function.[24]Recently,the presented authors(Shahzad and He)have reported a detail discussion on thermal conductivity calculation and Ewald-Yukawa sum for the case ofJQthat corresponds to phase space variable.[16]

3 HNEMD Results and Discussion

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λ0=λ/nmωpa2ws(normalized by plasma frequency)and orλ∗=λ/nmωEa2ws(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 perturbationP∗(=Pzaws/JQZ).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 ofP∗.In our case,the possible low value of external perturbation isP∗=0.02 at which 2D complex plasma system gives equilibrium thermal conductivity for all plasma state points.It is interesting and

When external force is selected parallel to thez-axisFe(t)=δ(0,Pz),δis Dirac delta function,the above Eq.(8)becomes significant here that this normalized steady state low value ofP∗=0.02 is very small as compared to earlier known value of external force in Ref.[23]This low value re flects more appropriate and acceptable results using presented HPMD technique than earlier used HNEMD technique.Theλ0results obtained through HPMD computer simulation are checked for the universal temperature scaling law at this reduced steady-state value ofP∗=0.02:

Equation(11)gives the simple scaling law(universal temperature law)and it is showed that the PHMD data calculated by usingλ∗=nmωEa2wsandT∗=T/Tm≡Γm/Γ(ratio of the system temperature to melting temperature),hereTmand Γmare the melting points and related detail is given in Refs.[15−19,23]Here,the unknown constants(A,B and C)are found after curve fitting to available HPMD simulation data for complex plasmas at different plasma state points(Γ,κ).

We now turn our attention to the main results obtained through the HPMD simulations.In our case,before the external perturbationP∗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 perturbationP∗.[3,23]Presently,for a possible low value of the external perturbation strength ofP∗=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 GKREMD 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λ0(Γ)are performed forN=400 particles and a sequence of four different computations are taken into account at constant perturbation ofP∗=0.02 for eachκ=1,2,and 3,respectively.Our numerical data ofλ0(Γ)are in nearly reasonable agreement with earlier numerical data of GKR-EMD,NEMD and HNEMD investigation.[4−5,23]

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]

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.

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λ0at 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λ0is well matched with earlier 2D numerical estimations[23]at intermediate Γ(=20).It is significant to note that a constantλ0is observed at intermediate to higher plasma coupling Γ at constant external perturbationP∗=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λ0is found at intermediate to higher Γ at constantP∗,con firming earlier numerical results.[4−5,23]It is interesting to note here that the existence ofλ0is present for low-intermediate to higher Γ with an increase inκand remains within a satisfactory limited statistical uncertainty,con firming previous simulation results.[23]In our simulation,the presented plasma conductivity for lower to intermediate Γ shows the existence ofλ0and it is a clear contradiction with the earlier simulation results of Donkoet al.[29]where theλ0was not found at lower Γ.

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λ0behaviors of the simulated complex dusty plasmas at higher value of screening.For this case,Fig.3 shows the normalizedλ0computed by the HPMD approach forN=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λ0at all plasma state points and it is observed that plasmaλ0found to be constant,as expected in earlier simulation results.[4]Moreover,it is noted that the measured data of plasmaλ0estimations atP∗=0.02,where plasmaλ0has equilibrium values and independent ofP∗,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 lowerNare more accurate and acceptable.

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 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λ∗(=λ/nmωEa2WS)verses various normalized temperaturesT∗(=T/Tm).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 temperatureT∗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.02302,B=−1.49422 andC=0.69373.These obtained fitting coefficients(A,BandC)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 coefficientBin 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,con firming 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κandT∗.For higher values ofκ,at intermediate values ofT∗,λ∗shows less dependence on these two variables and it is little far from functional form(temperature scaling).But at higherT∗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 ofT∗for the reduced force field strengthP∗=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.

4 Conclusions

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 temperatureT∗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.

Acknowledgments

We are very obliged to the National Advanced Computing Centre of National Centre for Physics(NCP),Pakistan and National High-Performance Computing Center(NHPCC)of Xian Jiaotong University,China for allocating computer time to test and run our MD code.

[1]Y.Feng,B.Yu,P.Xu,and M.Zou,J.Phys.D40(2007)3164.

[2]K.K.Mandadapu,R.E.Jones,and P.Papadopoulos,J.Chem.Phys.130(2009)204106.

[3]A.Shahzad and M.G.He,AIP Conf.Proc.1547(2013)173.

[4]L.J.Hou and A.Piel,J.Phys.A42(2009)214025.

[5]Y.V.Khrustalyov and O.S.Vaulina,Phys.Rev.85(2012)046405.

[6]W.Yu,D.M.France,J.L.Routbort,and S.U.Choi,Heat Trans.Eng.29(2008)432.

[7]A.J.H.McGaughey and M.Kaviany,Int.J.Heat Mass.Trans.47(2004)783.

[8]G.Ciccotti,G.Jacucci,and I.R.McDonald,J.Stat.Phys.21(1979)01.

[9]W.G.Hoover and W.T.Ashurst,Nonequilibrium Molecular Dynamics,Academic London,London(1975).

[10]D.J.Evans and G.P.Morriss,Statistical Mechanics of Non-Equilibrium Liquids,Academic London,London(1990).

[11]Z.Donko and P.Hartmann,Phys.Rev.E69(2004)016405.

[12]F.Muller-Plathe,J.Chem.Phys.106(1997)6082.

[13]J.P.Hansen and I.R.McDonald,Theory of Simple Liquids,Academic London,London(1986).

[14]D.J.Evans,Phys.Lett.A91(1982)457.

[15]R.J.Hulse,R.L.Rowley,and W.V.Wilding,Int.J.Thermo.Phys.26(2005)01.

[16]A.Shahzad and M.G.He,Contrib.Plasma.Phys.52(2012)667.

[17]R.Kubo,Rep.Prog.Phys.29(1966)255.

[18]G.Salin and M.J.Caillol,Phys.Plasmas.10(2003)1220.

[19]G.Faussurier,M.S.Murillo,and Gibbs-Bogolyubov,Phys.Rev.E67(2003)046404.

[20]Z.Donko,J.Phys.A:Math.Theor.42(2009)214029.

[21]A.Shahzad,M.G.He,S.Irfan Haider,and Y.Feng,Phys.Plasmas.24(2017)093701.

[22]A.Shahzad and M.G.He,Phys.Plasmas.19(2012)083707.

[23]A.Shahzad and M.G.He,Phys.Plasmas.22(2015)23707.

[24]M.J.Gillan and M.Dixon,J.Phys.C16(1983)869.

[25]C.Pierleon,G.Ciccotti,and B.Bernu,Euro.Phys.Lett.4(1987)1115.

[26]E.Wigner,Phys.Rev.46(2004)1002.

[27]H.Ohta and S.Hamaguchi,Phys.Plasmas7(2000)4506.

[28]T.Saigo and S.Hamaguchi,Phys.Plasmas9(2002)1210.

[29]Z.Donko,J.Goree,P.Hartmann,and B.Liu,Phys.Rev.E79(2009)026401.