Thermal conductivity of strongly coupled Yukawa liquids
Z. Donko´ and P. Hartmann
Research Institute for Solid State Physics and Optics of the Hungarian Academy of Sciences, P.O. Box 49, H-1525 Budapest, Hungary 共Received 15 August 2003; published 28 January 2004兲
The thermal conductivity of strongly coupled Yukawa liquids, being relevant to dusty plasmas, is calculated from nonequilibrium molecular dynamics simulations. The calculations cover a wide range of plasma coupling (⌫) and screening () parameters and yield data which are generally in good agreement with the results of recent independent calculations. An improved analytical formula, relating the thermal conductivity to the reduced temperature and to the screening length, is proposed.
DOI: 10.1103/PhysRevE.69.016405 PACS number共s兲: 52.27.Gr, 52.25.Fi, 52.27.Lw I. INTRODUCTION
The first set of data for the transport parameters 共self- diffusion coefficient, shear and bulk viscosities, and thermal conductivity兲 of strongly coupled one-component plasmas 共OCP兲 have been derived in the 1970s from theoretical cal- culations关1–3兴and from computer simulations关4兴. The self- diffusion coefficient D was found to be a monotonically de- creasing function of the coupling strength of the plasma, ⌫
⫽(Q2/40)(1/akBT), where Q is the charge of the par- ticles, T is the temperature, a⫽(4n/3)⫺1/3 is the Wigner- Seitz 共WS兲radius, and n is the density of particles. In con- trast with the behavior of D, it was also shown that the shear viscosityexhibits a minimum at intermediate values of the coupling coefficient, ⌫⬇ 10, and that the bulk viscosity is orders of magnitude lower compared to the shear viscosity.
The first computer simulation results for the thermal conduc- tivity have been obtained only for a few⌫ values关4兴, but already suggested that has a minimum at around the same value of⌫, as the shear viscosity.
More detailed data for the shear viscosity and the thermal conductivity of the classical OCP have been obtained from recent nonequilibrium 共transient兲 molecular dynamics共MD兲 simulations关5,6兴. These simulations were based on the mea- surement of the relaxation time of the ensemble of particles after introducing a specific perturbation to the system. In the case of thermal conductivity ‘‘measurements’’ a sinusoidal temperature profile was imposed on the system and the re- laxation of this profile was monitored. For the calculations of the shear viscosity the decay of the perturbation of the ve- locity field was measured 关6兴—these simulations have con- firmed the results of the earlier calculations of , both in terms of the⌫value wherehas its minimum, and the value ofmin. Simulations aiming the determination of the thermal conductivity 关5,6兴 confirmed that has a minimum at ⌫
⬇10–15.
In the last few years increased attention has been devoted to particle systems interacting through the screened Coulomb 共Yukawa兲potential, due to the relevance of this model sys- tem to dusty plasmas. Self-diffusion properties of Yukawa systems have been explored by Ohta and Hamaguchi关7兴and Liu et al. 关8兴. The shear viscosity was studied by Sanbon- matsu and Murillo 关9兴 and Saigo and Hamaguchi关10兴. The results forpresented in Ref.关10兴in the low-screening limit confirmed the nonequilibrium MD results of Donko´ et al.
关5,6兴. Very recent comprehensive studies of Yukawa trans- port coefficients by Salin and Caillol 关11,12兴 and by Faus- surier and Murillo 关13兴 indicate current interest in the sub- ject.
The aim of this paper is to present calculations of the thermal conductivity of Yukawa systems, based on a differ- ent simulation method than those used previously. The cal- culations are carried out for a wide range of characteristic parameters of the system共wider range of⌫ than in any pre- vious work兲, and the results are compared with those derived in recent independent calculations关12,13兴.
Section II of the paper describes the simulation technique and Sec. III presents the results and their comparison with data published by other authors. The work is summarized in Sec. IV.
II. SIMULATION METHOD
The motion of the particles in the system is described by molecular dynamics simulation. N simulation particles are contained in a cubic computational box having periodic boundary conditions in all three space dimensions, and they interact through the Yukawa potential
共r兲⫽ Q2 40
exp共⫺r/D兲
r , 共1兲
where D is the Debye 共screening兲 length. The screening parameter ⫽a/D, together with the coupling parameter
⌫, fully characterizes the equilibrium system. In the→ 0 limit the interaction reduces to the 1/r Coulomb type 共OCP limit兲.
The calculation of the thermal conductivity is based on the nonequilibrium molecular dynamics 共NEMD兲 approach introduced by Mu¨ller-Plathe 关14兴. The simulation box is di- vided into several slabs along one of the axes (x axis in our case兲and two of the slabs are assigned to be the ‘‘cold slab’’
and a ‘‘hot slab,’’ as illustrated in Fig. 1. In each kth simu- lation time step the fastest particle in the cold slab and slow- est particle in the hot slab are searched for and their mo- menta are exchanged. This artificial perturbation of the system creates a heat flux between the hot and the cold slabs and gives rise to a specific spatial temperature profile throughout the system关14兴.
At the initialization of the simulation the particles are loaded at randomly chosen positions into the cubic compu- tational box. Their initial velocities are randomly directed with magnitudes sampled from a Maxwellian distribution corresponding to the specified value of ⌫. In a number of initial time steps the velocities of the particles are slightly adjusted to eliminate the possible small drift of the system, i.e., to ensure that 具vx典⫽具vy典⫽具vz典⫽0 (具•典 denotes aver- aging over the ensemble of particles兲. The perturbation giv- ing rise to the spatial temperature profile 共as explained above兲 is applied from the beginning of the simulation and the system is given sufficient time 共10 000 and 5000 time steps, respectively, in the case of N⫽1600 and 6400 par- ticles兲to relax to the stationary 共perturbed兲 state before the measurement 共data collection兲phase starts. During this ini- tial phase of the simulation the momenta of the particles are regularly rescaled in order to keep the average temperature of the system at the desired value. The number of simulation time steps in the data collection phase is 5⫻105 for N
⫽1600 and ranges between 5⫻104 and 1⫻105 for N
⫽6400. For these conditionspt 共wherep⫽
冑
Q2n/0m is the plasma frequency兲is in the order of several thousands.The above momentum-rescaling procedure is also applied during the data collection phase, however, only in every 5000th time step, to compensate for the long term accumu- lation of numerical errors. 共It is noted that the system tem- perature during the data collection phase does not rise in an observable way even without this additional thermostation.兲
spatial and temporal average of the temperature gradient 具⌬T/⌬x典between the cold and hot slabs is measured and the thermal conductivity is calculated as
⫽ E
2St
冓
⌬⌬Tx冔
, 共2兲where E is the 共cumulative兲 energy exchanged between the hot and cold slabs during the runtime of the simulation t and S⫽L2 is the cross sectional area of the simulation box. This NEMD technique has already been successfully applied to calculations of thermal conductivity of Lennard-Jones and molecular fluids 关14 –16兴.
The thermal conductivity resulting from Eq.共2兲can con- veniently be expressed in reduced units. Most common are the normalizations by the plasma frequency p or by the Einstein frequencyE关13兴, respectively
⬘
⫽ nkBpa2 *⫽ nkB
冑
3Ea2, 共3兲
whereEis the oscillation frequency of a particle around its equilibrium position when all other particles are at their equi- librium lattice sites, see, e.g., Ref. 关7兴. The Einstein fre- quency decreases with increasing , from a value E
⫽p/
冑
3 at ⫽0 关7兴. Here we follow the notation used in Ref. 关13兴, while in our earlier work 关5,6兴 and in Ref. 关12兴/nkBpa2 was denoted by*.
III. RESULTS
Our simulations are carried out with number of particles chosen to be N⫽1600 and N⫽6400. These two different settings make it possible to investigate system size effects.
The simulation box is divided into 32 slabs.
The temperature gradient establishing in the system de- pends on the frequency of perturbations. This way—by ex- changing the momenta of the two selected particles in the cold and hot slabs relatively rarely—arbitrarily small tem- perature gradients can be established, causing very small per- turbation to the system. However, at very low⌬T/⌬x noise may be dominant, thus a proper choice of the delay between momenta exchanges is a compromise between the degree of perturbation and signal to noise ratio of the quantities (E/t and⌬T/⌬x) to be measured. In Fig. 2 the T(x) profiles are plotted for different values of k 共number of time steps be- tween momenta exchanges兲, k⫽75, 150, and 300. The T(x) profiles are very nearly linear, with a slope⌬T/⌬x depend- ing on k. In the further simulations the frequency of the exchange of momenta between the selected particles in the cold and hot slabs is chosen to be between k⫽50 and 300 time steps, to result in a ⱗ⫾10% perturbation of the tem- perature profile.
The results of our calculations for the thermal conductiv- ity normalized by the plasma frequency
⬘
are shown in Fig.3 for⫽0.1, and in Fig. 4 for⫽1, 2, and 3. Figures 3 and FIG. 1. Illustration of the nonequilibrium molecular dynamics
simulation method originally proposed by Mu¨ller-Plathe 关14兴. Lighter and darker colors indicate colder and hotter parts of the simulation box. The exchange of particles’ momenta between the cold and hot slabs results in a heat flux and a temperature gradient, which is shown in the lower part of the figure.
4 also include the data published by Salin and Caillol 关12兴 and a universal functional form for the thermal conductivity proposed by Faussurier and Murillo 关13兴
*⫽0.01176T*⫹0.881
T* ⫹0.1655, 共4兲
converted to
⬘
as ⬘
⫽*(冑
3E/p), where the normal- ized temperature T*⫽⌫c/⌫, with ⌫c being the couplingstrength corresponding to melting共which, as a function of is given in Ref. 关17兴兲. This functional form was found to describe properly the scaling properties of the shear viscos- ity, as explored by Saigo and Hamaguchi 关10兴. Additionally, Fig. 3 also shows the data given in Refs. 关5,6兴for the Cou- lomb OCP, obtained with a different simulation technique.
In the limit of small共see Fig. 3兲the result of the present simulations agree well with the Coulomb OCP data 关5,6兴, FIG. 2. Spatial temperature profiles obtained at ⌫⫽20,
⫽0.1, for different number of time steps (k) between artificial changes of particle momenta in the cold and hot slabs, situated at x/L ⫽0.25 and 0.75, respectively.
FIG. 3. Reduced thermal conductivity 关normalized by the plasma frequency共3兲兴for 0⭐⭐0.1—the precise values are given in the legend. Bernu, Viellefoss, and Hansen共BVH兲 关4兴, Donko´ and Nyı´ri共DN兲 关6兴 共obtained with 8192 and 1024 particles兲, Salin and Caillol共SC兲 关12兴, Faussurier and Murillo共FM兲 关13兴, and Donko´ and Harmann共DH兲: present results共with 6400 and 1600 particles兲. Note that*⯝⬘at⬇0.
FIG. 4. Reduced thermal conductivity 关normalized by the plasma frequency共3兲兴at共a兲⫽1, 共b兲 ⫽2, and共c兲⫽3. Salin and Caillol共SC兲 关12兴 and Faussurier and Murillo 共FM兲 关13兴, and Donko´ and Harmann 共DH兲: present results 共with 6400 and 1600 particles兲.
⌫
agree in terms of the position of the minimum of
⬘
共occur- ring at⌫⬵10–15) as well as of the valuemin⬘
⬵0.4. At the lowest ⌫ values the present simulations give a slightly (⬇10%) higher ⬘
, compared to Refs. 关5,6兴. Around the minimum position, at⌫⫽10 the present value is well below the result of Bernu et al. 关4兴. At⌫ ⫽ 100 the difference is still about 20%. The data of Salin and Caillol 关12兴 lie be- tween the present data and those of Ref. 关4兴; they are 20–30 % higher compared to the present results. At⬇0 the fit 共4兲approximates well the present data, especially at interme- diate values of⌫.
The results obtained at the higher values are shown in Fig. 4. The present data indicate that the position of the mini- mum of
⬘
shifts towards higher⌫asincreases. The mini- mum value of ⬘
decreases with increasing , to min⬘
⬵0.24 at ⫽3. 共It is noted that the thermal conductivity normalized by the Einstein frequency* increases with in- creasing, to a valuemin* ⬵0.8 at⫽3.兲The data are gen- erally in a good agreement with the results of Salin and Caillol 关11,12兴 given for ⌫⭐100 range. At high and ⌫
⬎100 the thermal conductivity exhibits an increasing ten- dency with increasing ⌫, similarly to the ⬵ 0 case. The analytical formula 共4兲reproduces well the shift of the posi- tion of min
⬘
with increasing. The agreement between the calculatedmin⬘
values and those given by Eq.共4兲is not very good, in the worst case, at⫽3, the data differ by a factor of 2. The simulation results obtained with N⫽6400 and 1600 particles agree within the limits of errors at high ⌫ values, but deviate at the lower values of⌫. This deviation shows up at increasingly higher values of ⌫ when increases共from⌫⬇3 at ⫽0.1 to⌫⬇20 at⫽3).
We believe that the system size 共N兲 dependence of the results at low coupling values originates from the long tra- jectory segments where particle motion is unperturbed by the other particles. In other words, the simulation scheme is ex- pected to be valid as long as the characteristic size of the system is bigger than the length scale for interaction共signifi- cant exchange of momenta兲between the particles. This latter can be quantitatively studied by calculating the velocity au- tocorrelation function of the system Avv共already discussed in Ref. 关7兴 for Yukawa systems兲 and its decay for different⌫ and values. Figure 5 shows Avvfunctions defined as
Avv共t兲⫽具v共t兲•v共t0⫹t兲典
具v共t0兲•v共t0兲典 , 共5兲 where具典 denotes averaging over particles and different ini- tial times t0. At high⌫ values Avv exhibits a fast decay as illustrated in Fig. 5共a兲for the⫽0.1 case. The decay time of Avv, defined as Avv(t⫽)⫽1/⑀(⑀⫽2.718) increases as⌫ is decreased. The time needed for a particle to proceed a distance L共where L is the edge length of the simulation box兲 is ⌬t⫽L/v¯ , wherev¯⫽
冑
3kBT/m is the average velocity of the particles. Considering the case of N⫽1600 particles, at⌫⫽5 we findp⌬t⫽42.2, which is about an order of mag- nitude greater than the decay timeof Avv关see Fig. 5共a兲兴. At higher values of we can observe a drastic increase of, as
shown in Figs. 5共b兲and 5共c兲for⌫⫽100 and ⌫⫽1, respec- tively. At⌫⫽100p⌬t⫽189, the decay time of Avvis well below this value, even for the greatestvalue studied. How- ever,p⌬t becomes 18.9 for⌫⫽1, and this is already in the same order of magnitude as the decay time of Avv, even for the smallest value of the screening parameter, ⫽ 0.1.
Based on these arguments we find the low-⌫ limit of the applicability of the present simulation technique to be:
FIG. 5. Velocity autocorrelation functions 共a兲 for different ⌫ values at ⫽0.1; 共b兲 and 共c兲 dependence on at⌫⫽100 and⌫
⫽1, respectively.
⌫min,1600⬇3 for ⫽0.1,⌫min,1600⬇5 for ⫽1,⌫min,1600
⬇10 for⫽2, and ⌫min,1600⬇20 for⫽3. The differences between the results obtained with different system sizes in- deed show up at coupling values ⌫ⱗ⌫min. For N⫽6400 particles the lower bounds of ⌫ are expected to be⌫min,6400
⫽2⫺2/3⌫min,1600, giving ⌫min,6400⬇2 at →0 and ⌫min,6400
⬇13 at⫽3.
The form 共4兲 taken from Ref.关13兴 relates * to the re- duced temperature T*. The plot of our* data against T*, as shown in Fig. 6共a兲, indicates that the above functional form cannot give an accurate fit. While it predicts the shape of the*(T*) curve correctly and accounts for the universal behavior that the minimum of the curves at different val- ues occurs nearly at the same T* 关see Fig. 6共a兲兴, the vertical shift of the curves can only be eliminated by introducing an additional (-dependent兲 term. This way, unfortunately, the universal behavior is lost, but the resulting expression pro- vides a more reliable approximation for the thermal conduc- tivity
*⫽0.018T*⫹1.05
T* ⫹0.115⫹0.127. 共6兲
The results of the present simulation共with N⫽6400) for all
values are shown in Fig. 6共b兲, together with the above expression for*. One should keep in mind, however, that the validity of this expression has a limited⌫ range, as dis- cussed above.
IV. SUMMARY
We have presented nonequilibrium molecular dynamics calculations for the thermal conductivity of strongly coupled Yukawa liquids.
In the limit of low screening a very good agreement has been found between the result of the present simulations and our previous OCP data obtained from a different 共transient兲 molecular dynamics simulation 关5,6兴. A reasonable agree- ment is also found between our data and recent results of Salin and Caillol 关12兴 for the whole domain of the system parameters⌫ and.
Studying the decay of the velocity autocorrelation func- tion we have determined the ⌫min 共applicability兲 limit as a function of and system size. Our simulations provide reli- able measurements of in the coupling range ⌫ⲏ2 at
⫽0.1. With increasing the range of applicability of the present technique was found to be ⌫ⲏ3 at ⫽1, ⌫ⲏ7 at
⫽2, and⌫ⲏ13 at⫽3 共using N⫽6400). The (⌫,) do- main where the present method is applicable is limited by the ballistic trajectories of particles at low⌫ and high values.
At lower ⌫ values our simulation method could still work with higher number of particles, however, limitations set by computing speed did not allow us to increase the number of particles above 6400. We believe that any other 共e.g., equi- librium兲 molecular dynamics simulations may have similar limitations set by the finite number of particles, at low⌫and high values, thus in any simulation studies the system size effects should be checked carefully over the whole domain of system parameters共as done in the present work兲.
The formula共4兲proposed by Faussurier and Murillo关13兴, relating*to the reduced temperature T*and, was found to reproduce correctly the shape of the *(T*) curves at different values of . We have found, however, that an ad- ditional-dependent term has to be added to Eq.共4兲in order to give a more accurate relationship between * and T*.
ACKNOWLEDGMENTS
Support through Grant Nos. OTKA-T-34156 and MTA- NSF-OTKA-028 is gratefully acknowledged. We thank S.
Hamaguchi and T. Saigo 共Kyoto University兲 for providing their Yukawa pair correlation function data for the verifica- tion of our simulation code, and G. Ba´no´ for his comments on the manuscript.
FIG. 6.共a兲Reduced thermal conductivity*共normalized by the Einstein frequency兲as a function of the reduced temperature T*for different values;共b兲 quasiuniversal behavior found by adding a
-dependent term to Eq.共4兲, resulting in a new formula,共6兲, shown by the heavy line.
11, 1025共1975兲.
关2兴P. Vieillefosse and J.P. Hansen, Phys. Rev. A 12, 1106共1975兲. 关3兴J. Wallenborn and M. Baus, Phys. Lett. 61A, 35共1977兲; Phys.
Rev. A 18, 1737共1978兲.
关4兴B. Bernu, P. Vieillefosse, and J.P. Hansen, Phys. Lett. 63A, 301 共1977兲; B. Bernu and P. Vieillefosse, Phys. Rev. A 18, 2345 共1978兲.
关5兴Z. Donko´, B. Nyı´ri, L. Szalai, and S. Hollo´, Phys. Rev. Lett.
81, 1622共1998兲.
关6兴Z. Donko´ and B. Nyı´ri, Phys. Plasmas 7, 45共2000兲. 关7兴H. Ohta and S. Hamaguchi, Phys. Plasmas 7, 4506共2000兲. 关8兴Y.H. Liu et al., J. Phys. A 35, 9535共2002兲.
共2001兲.
关10兴T. Saigo and S. Hamaguchi, Phys. Plasmas 9, 1210共2002兲. 关11兴G. Salin and J.-M. Caillol, Phys. Rev. Lett. 88, 065002共2002兲. 关12兴G. Salin and J.-M. Caillol, Phys. Plasmas 10, 1220共2003兲. 关13兴G. Faussurier and M.S. Murillo, Phys. Rev. E 67, 046404
共2003兲.
关14兴F. Mu¨ller-Plathe, J. Chem. Phys. 106, 6082共1997兲.
关15兴D. Bedrov and G.D. Smith, J. Chem. Phys. 113, 8080共2000兲. 关16兴D. Bedrov, G.D. Smith, and T.D. Sewell, Chem. Phys. Lett.
324, 64共2000兲.
关17兴S. Hamaguchi, R.T. Farouki, and D.H.E. Dubin, Phys. Rev. E 56, 4671共1997兲.
