Equilibrium properties and phase diagram of two-dimensional Yukawa systems

Equilibrium properties and phase diagram of two-dimensional Yukawa systems

P. Hartmann,¹ G. J. Kalman,²Z. Donkó,¹ and K. Kutasi¹

1Research Institute for Solid State Physics and Optics of the Hungarian Academy of Sciences, H-1525 Budapest, P.O. Box 49, Hungary

2Department of Physics, Boston College, Chestnut Hill, Massachusetts 02467, USA 共Received 26 April 2005; published 22 August 2005兲

Properties of two-dimensional strongly coupled Yukawa systems are explored through molecular dynamics simulations. An effective coupling coefficient⌫^*for the liquid phase is introduced on the basis of the con- stancy of the first peak amplitude of the pair-correlation functions. Thermodynamic quantities are calculated from the pair-correlation function. The solid-liquid transition of the system is investigated through the analysis of the bond-angular order parameter. The static structure function satisfies consistency relation, attesting to the reliability of the computational method. The response is shown to be governed by the correlational part of the inverse compressibility. An analysis of the velocity autocorrelation demonstrates that this latter also exhibits a universal behavior.

DOI:10.1103/PhysRevE.72.026409 PACS number共s兲: 52.27.Gr, 05.20.⫺y, 73.21.⫺b

I. INTRODUCTION

The Yukawa 共screened Coulomb兲 potential is a widely used approximation to describe the interaction of particles in a variety of physical systems, e.g., dusty plasmas in labora- tory关1兴and space 关2兴and charged colloids关3兴. The Yukawa interaction potential 共energy兲 is ␾共r兲=共Q²/r兲exp共−␬r兲.

Many-particle systems with Yukawa interaction can fully be characterized by two dimensionless parameters:共i兲 thecou- pling parameter ⌫=␤Q²/a 共where Q is the charge of the particles,a is the Wigner–Seitz radius, and␤= 1 /k_BT is the inverse temperature兲, and 共ii兲 the screening parameter ␬^¯

=␬a.

Besides three-dimensional共3D兲systems, two-dimensional 共2D兲configurations also appear in a wide variety of physical systems. As examples of such Coulomb or Yukawa systems, respectively, the layer of electrons on the surface of liquid helium or layers of dust particles formed in low-pressure gas discharges may be mentioned关4–8兴. The dynamics of these 2D systems has extensively been investigated both theoreti- cally and by simulation techniques, e.g., 关9–12兴, but much less so insofar as the study of the equilibrium properties of such systems is concerned共see, however, a recent paper by Totsuji et al. 关10兴 and 关13–15兴兲. This paper addresses this latter issue: we report on molecular-dynamics共MD兲simula- tions through which various equilibrium properties of a 2D Yukawa system are explored. In Sec. II we describe our simulation technique. In Sec. III we display pair-correlation functions and show how most of the equilibrium properties of a 2D Yukawa system can be described in terms of an effective coupling parameter. In Secs. IV and V, respectively, we calculate the energy and pressure, and display the phase boundary between the liquid and solid phases. Sections VI and VII are devoted to the related issues of the static struc- ture factor, the static dielectric function, and the behavior of the isothermal compressibility. Finally, in Sec. VIII we present simulation results for the velocity autocorrelation function and provide a qualitative explanation of its observed features.

II. SIMULATION TECHNIQUE

The motion of the particles is traced by molecular- dynamics technique based on the PPPM 共particle-particle particle-mesh兲algorithm关16兴, which allows efficient simula- tion of systems composed of a large number of particles关17兴. We use a variant of the method modified for a Yukawa po- tential. The simulation domain is a three-dimensional cubic box with periodic boundary conditions in which the particles are constrained to move within a two-dimensional layer. In the PPPM method the interparticle force is partitioned into 共i兲a force componentF_PMthat can be calculated on a mesh 共the “mesh force”兲and共ii兲a short-range forceF_PP, which is to be applied to closely separated pairs of particles only共for more details see关16兴兲. In this way the PPPM method makes it possible to take into account periodic images of the system 共in the PM part兲, without truncating the long-range Coulomb or low-␬^¯ Yukawa potentials. 共For high␬^¯ values the PP part alone provides sufficient accuracy, in these cases the mesh part of the calculation is not used.兲In most simulations, the desired system temperature is reached by rescaling the par- ticle momenta during the initialization phase of the simula- tion. The system temperature does not show an observable drift during the subsequent measurement phase of the simu- lation that follows the initialization period. In the simulations aimed at the determination of the solid-liquid phase bound- ary the Nosé–Hoover thermostat共e.g.,关18兴兲is used to ensure a correct control of the system temperature.

The primary output data of our simulations are the pair- correlation functions 共PCFs兲 g共r兲, which contain indispens- able information about the system关19兴. The PCFs are used as input data for the calculation of the correlational energy, pressure, and compressibility of the system, and for the cal- culation of the static structure function S共k兲 as well. Addi- tional quantities, e.g., velocity autocorrelation functions 共VACF兲 are also analyzed. The solid-to-liquid transition is studied through monitoring the temperature dependence of the bond-angular order parameter共see below兲.

Our simulations also provide information about the spec- tra of the longitudinal and transverse current fluctuations in

the liquid phase of the system where isotropy can be as- sumed关11兴. These spectra are obtained through the Fourier transform关20,21兴

L共k,␻兲=兩F兵␭共k,t兲其兩²,

T共k,␻兲=兩F兵␶共k,t兲其兩², 共1兲 of the microscopic quantities

␭共k,t兲=k

兺

j

vjxexp共ikxj兲,

␶共k,t兲=k

兺

_j ^{vjyexp共ikxj兲, 共2兲}

where the index jruns over all particles.

The spectra defined by共1兲serve as the basis for the analy- sis of the collective excitations of the system. We have re- ported detailed calculations on this topic in关11兴.

In the following, the simulation results are given with the length scale normalized by the 2D Wigner–Seitz radius a

= 1 /

冑

n␲ 共where n is the areal density兲, i.e.,¯r=r/a and¯k

=ka.

III. SCALING

The issue of scaling共i.e., whether only some combination of the⌫ and␬^¯ parameters rather than both of these param- eters independently or, alternatively, the ratio of the tempera- ture to the melting temperature govern the behavior of Yukawa systems兲has been addressed by several studies; the universal scaling of structural properties and transport pa- rameters has continued to receive attention for many years 关22–29兴.

In the case of three-dimensional Yukawa systems, the sug- gestion for an effective coupling parameter ⌫^*=⌫exp共−^¯␬兲 was given by Ikezi关30兴and was based on the argument that the coupling for a Yukawa potential is obtained from the Coulomb coupling parameter by multiplying the latter by exp共−␬^{¯ r¯}兲. Subsequent studies, however, have shown that dif- ferent systems characterized by the same ⌫^* do not show many similarities.关According to our 2D simulations, at cer- tain values of⌫^*=⌫exp共−^¯␬兲 共e.g., 120兲, the system may be either in the liquid or in the solid phase, depending on the value of␬^¯.兴A more recent and more meaningful definition for the effective coupling parameter was given by Vaulina,et al.关25兴, Vaulina and Vladimirov关26兴, and Fortovet al.关27兴. Their definition of⌫^*is based on the frequency of dust lat- tice waves, resulting in ⌫^*=⌫共1 +␬^¯+␬^¯2/ 2兲exp共−␬^¯兲. They have shown that within the 0艋^¯␬艋5 range, their effective coupling parameter has a nearly constant value ⌫^*=⌫m

3D

共where⌫m3D⬵175 is the coupling coefficient value at which the 3D one-component plasma melts兲along the solid-liquid phase boundary determined by Hamaguchiet al.关31兴. More- over, it was also demonstrated 关27兴 that systems with the same ⌫^* value have very similar pair-correlation functions 关25–27兴. In our earlier work 关11兴we considered defining an effective coupling parameter for a two-dimensional Yukawa system. We established a different criterion for⌫^*that relies

on associating a constant amplitude of the first peak of the PCF关g共r兲兴with a constant⌫^*value.

Our aim now is to find a properf共^¯␬兲function that allows one to partition⌫^*共⌫,␬^¯兲as

⌫^*=⌫f共^¯␬兲. 共3兲 In order to obtain the f共␬^¯兲 function we have carried out an extensive scan of the PCFs over a wide domain of the共⌫,␬^¯兲 parameter space. The pair-correlation functions of the 2D Yukawa liquid are displayed in Fig. 1共a兲for⌫= 120 and for a series of^¯␬values. It can be seen that the range of the rather pronounced order characteristic for^¯␬= 0 rapidly diminishes with increasing␬^¯. The amplitude of the first peak of the PCF can, however, be reestablished if⌫is also increased together with␬^¯. This is illustrated in Fig. 1共b兲. In fact, as Fig. 1共b兲 shows, within the range of¯rdisplayed not only the amplitude of the first peak, but theg共r¯兲functions in their entireties are nearly the same for fixed⌫^*values, even though the associ- ated␬^¯ values are quite different. A caveat, however, has to be added. As it is discussed below, the values of the static struc- ture functionS共k兲atk= 0 are sensitive to the variation of^¯␬^, even at ⌫^*= constant values. Since S共k= 0兲= 1 +n兰h共r兲dr 关where h共r兲=g共r兲− 1 and r lies in the particle plane兴, the difference must come from the tail ofg共r兲, which is, appar- ently “unscalable” and the scaling property of g共r兲 cannot extend tor→⬁.

Figure 2共a兲 shows the contours on the ⌫−^¯␬ plane that belong to constant effective coupling values⌫^*= 120, 40, and FIG. 1. Pair-correlation functions vs distance in units of the Wigner-Seitz distanceaof the 2D Yukawa liquid for共a兲⌫= 120 and different values of¯, and␬ 共b兲for a series of共⌫,¯␬兲pairs correspond- ing to a constant⌫^*= 120.

HARTMANNet al. PHYSICAL REVIEW E72, 026409共2005兲

10. It can be seen that these lines have approximately the same shape; thus, they can be scaled to a single universal line, as shown in Fig. 2共b兲, which displays the dependence of the ratio ⌫/⌫^* on ␬^¯ for the chosen values of ⌫^*. At high values of⌫^*, the ratio⌫/⌫^*depends only on␬^¯, the partition- ing given in 共3兲 is indeed possible, and f共^¯␬兲 can be fitted with the aid of the following:

f共␬^¯兲= 1 +f₂␬^¯2+f₃␬^¯3+f₄^¯␬⁴, with

f₂= − 0.388; f₃= 0.138; f₄= − 0.0138. 共4兲 It can be seen in Fig. 2共b兲that the results can accurately be represented by Eq.共4兲except for the lowest⌫^*= 10 value.

IV. THERMODYNAMIC QUANTITIES

The energy共per particle兲of the systemE consists of the thermal partE₀, the positive Hartree partE_H, and the nega- tive correlational partE_c:

␤E₀= 1

␤E_H=␤ⁿ

2

冕

^{␾共r兲dr=␤Q2}␬^␲n=

⌫

␬

¯

␤^Ec=␤ⁿ

2

冕

^{h共r兲␾共r兲dr=⌫}

冕

0

⬁h共r¯兲e^{−␬¯ r¯}dr¯. 共5兲

Similarly, the pressure consists of the ideal, the positive Hartree and the negative correlational contributions

␤^P0=n

␤P_H=␤ⁿ²

2

冕

^{␾共r兲dr=␲Q}␬²ⁿ²⁼␤^⌫

␬

¯n

␤^Pc= −␤ⁿ²

4

冕

^r⳵␾⳵^{共r兲r h共r兲dr=} n⌫

2

冕

0

⬁

¯r

冋

^␬¯+¯1r

册

^{e−␬¯ r¯h共r¯兲dr¯.}

共6兲 The correlation energy per particle关calculated according to共5兲兴 and the correlational part of the pressure关calculated according to 共6兲兴 can now be obtained from the PCF. The former is plotted in Fig. 3共a兲as a function of⌫, for different values of ␬^¯. Plotting the same data as a function of ^¯␬ ^for constant values of⌫^*shows thatE_cdepends only slightly on

⌫^*in the strong coupling domain, but varies strongly with^¯␬, as expected, since onlyh共r兲—but not␾共r兲—is scaled within the integral. The data can be approximated as

␤E_c=⌫关b共␬^¯兲+c共␬^¯兲⌫^*−2/3兴, with

b共␬^¯兲=b₀+b₁␬^¯+b₂^¯␬²+b₃^¯␬³+b₄^¯␬⁴, and FIG. 2. 共a兲 Constant effective coupling共⌫^*兲 lines on the⌫−␬¯

plane. The dashed line represents Ikezi’s definition⌫^*=⌫exp共−␬¯兲 for⌫^*= 10.共b兲Dependence of the ratio⌫/⌫^*on␬¯. The symbols are data taken from共a兲, while the solid line is a fit according to共3兲and 共4兲.

FIG. 3. 共a兲 Correlation energy per particle of the 2D Yukawa liquid as a function of⌫, for selected values of¯␬. Lines: present results, symbols: data of Totsujiet al.关10兴.共b兲Correlation energy per particle as a function of␬¯, for constant values of the effective coupling parameter⌫^*.

c共␬^¯兲=c₀+c₁␬^¯+c2␬^¯2+c3␬^¯3+c4␬^¯4, 共7兲 where b₀= −1.103; b₁= 0.505; b₂= −0.107; b₃= 0.00686; b₄

= 0.0005 andc₀= 0.384;c₁= −0.036;c₂= −0.052;c₃= 0.0176;

c₄= 0.00165.

Our data shown in Fig. 3共a兲are in an excellent agreement with the energy values recently calculated by Totsuji et al.

关10兴 and at ␬= 0 with the energy values given for the 2D one-component plasma共OCP兲 关5,10,32兴.

The correlational part of the pressure is plotted in Fig.

4共a兲as a function of⌫, for different values of^¯␬. Similarly to the energy, the data are in an excellent agreement with those of Totsujiet al.关10兴. In the 5艋⌫艋120 and 0.5艋^¯␬艋3 in- terval the correlational part of the pressure共Pc兲can be fitted using the form

␤P_c=n⌫共b₀

⬘

+b₁

⬘

^¯␬兲, 共8兲 whereb₀

⬘

= −0.5638 andb₁

⬘

^{= 0.09367.}

It is interesting to note that while in the Coulomb case the pressure and the energy are linked through the virial theorem as P_c=¹₂nE_c, no such relationship exists for ^¯␬⫽0 关see 共6兲 and Fig. 4共b兲兴.

V. PHASE DIAGRAM

The issue of the number of distinct phases in a 2D Cou- lomb or Yukawa systems has been, for some time, a matter of intense controversy关33兴. The Kosterlitz–Thouless–Halperin–

Nelson–Young 共KTHNY兲 theory of the 1970s 共see, e.g., 关34兴兲suggested that in the course of the transition from the liquid phase to the crystalline phase with quasi-long-range

positional order and a long-range orientational order, an in- termediate “hexatic” phase exists without a long-range posi- tional order but dominated by a quasi-long-range orienta- tional order.

The transitions between the different phases can be char- acterized by the behavior of the bond-angular order param- eterG_⌰共e.g.,关35兴兲. For a system with hexagonal symmetry, the bond-angular order parameter关35,36兴has the form

G_⍜= 1

N

冏 ^兺

^{k=1N 16m=1}

^兺

^{6 exp共i6⌰k,m兲}

冏

^{2, 共9兲}

where the indiceskruns over all particles of the system, m runs over the six nearest neighbors neighbors of thekth par- ticle, and ⌰k,m is the angle between a predefined 共e.g., x兲 direction and the vector connecting thekth andmth particles.

For the determination of the solid-liquid phase diagram we follow Schweigertet al.关36兴who investigated systems char- acterized by different interaction potentials and inferred that the solid-to-liquid transition can be identified by a drop of the bond-angular order parameter below a value of G_⌰

⬵0.45;共see also关13,14兴兲.

The melting “experiment” of the 2D Yukawa layer is il- lustrated in Fig. 5. At the beginning of the simulation, a particle configuration with random particle positions is ini- tialized with an initial temperature well above the expected value of melting temperature. The temperature of the system is slowly decreased afterward, and ample time is given to the system to reach an ordered configuration 共as shown in the inset of Fig. 5兲. Subsequently, the temperature is slowly in- creased and the bond-angular order parameterG_⌰ is calcu- lated according to共9兲in each time step. As a consequence of the increasing temperature, first we observe a slow decay of G_⌰共from an initial value close to 1.0, indicating nearly per- fect hexagonal order兲 and later on, when the temperature reaches a certain value,G_⌰suddenly drops to⬇0, indicating FIG. 4.共a兲Correlational part of the pressure共␤P_c/n兲as a func-

tion of⌫ for␬¯= 0.5, 1, 2, and 3. The dashed line shows the theo- retical behavior of the pure Coulomb OCP.共b兲The ratio P_c/共nE_c兲 as a function of⌫.

FIG. 5. Illustration of the “melting experiment:” time depen- dence of the bond-angular order parameterG_⌰and system tempera- tureT, obtained at¯␬= 2. The sudden decay ofG_⌰below the 0.45 value 关36兴—marking the solid → liquid transition—occurs at ⌫m

= 384. The inset shows a snapshot of particle positions recorded right before the temperature starts to increase.

HARTMANNet al. PHYSICAL REVIEW E72, 026409共2005兲

an abrupt loss of the long-range orientational order in the system. We identify this event, taking place at⌫=⌫m, as the solid to liquid transition. The temperature control of the sys- tem is realized by the Nosé–Hoover algorithm 共see, e.g., 关18兴兲. To obtain the⌫m−␬^¯ phase diagram of the system, the simulations are carried out for a series of screening param- eter values. All results presented here have been obtained usingN= 1600 particles.

The⌫m−^¯␬phase boundary obtained from simulations il- lustrated above is plotted in Fig. 6. At␬^¯= 0 the simulations closely reproduce the well-known value ⌫mCoulomb⬵137 for the 2D OCP. The figure also shows the⌫values calculated from共4兲, assuming⌫^*= 131. We find an excellent agreement between the two sets of data. This agreement shows that the first peak amplitude of the PCF is nearly constant along the melting line of 2D Yukawa systems, regardless of the value of^¯␬, as already pointed out before. Studies of Vaulinaet al.

关25兴, and Vaulina and Vladimirov关26兴, and Fortovet al.关27兴 have actually reached a similar conclusion for 3D Yukawa systems.

An early attempt to theoretically determine the solid- liquid phase boundary for a 2D Yukawa system is due to Peeters and Wu 关15兴. Their calculation is based on the KTHNY theory of dislocation-mediated melting关34,35兴. The structure of the phase diagram obtained in this work is not corroborated by our MD results: Peeters and Wu关15兴find a significantly共⬇40%兲lower melting⌫than is shown in Fig.

6. The possible presence of a hexatic phase共as predicted by the KTHNY theory兲was not investigated in关15兴.

VI. STATIC STRUCTURE FUNCTION

The static structure function is derived from the PCFs as

S共k兲= 1 + 2␲ⁿ

冕

0

⬁

h共r兲rJ0共kr兲dr. 共10兲

Representative plots of S共k兲 obtained from the PCF through共10兲are displayed in Figs. 7共a兲and 7共b兲 for the⌫^*

= 40 and 10 values, respectively. As the S共k兲 functions are derived from the PCFs we expect theS共k兲functions belong- ing to the same ⌫^* to show similarity. On the other hand, S共k= 0兲 is expected to increase with increasing ␬^¯: for the Coulomb caseS共k= 0兲= 0, whereas Yukawa systems are char- acterized by S共k= 0兲⬎0. As this effect is also emphasized with decreasing⌫, the increase ofS共k= 0兲is clearly observed in the case of⌫^*= 10, while it is less pronounced at⌫^*= 40.

The value ofS共k兲fork→0 is governed by the compress- ibility sum rule for ␹^¯, the “screened” 共total兲 static density response function, which can be written in the form

␹

¯共k→0兲= − ␤ⁿ

1 +L_c+Mk¯², 共11兲 whereL_cis the correlational part of the inverse compressibil- ity andM represents the correlation-induced contribution of O共k²兲 in the ␹^¯共k→0兲 expansion. From the fluctuation- dissipation relationship

S共k兲= − ^¯␹共k兲/␤ⁿ

1 −␾共k兲^¯␹共k兲 共12兲 关␾共k兲= 2␲^Q2/

冑

_␬²+k²is the Fourier transform of the Yukawa potential energy兴, one finds

S共k¯兲=S₀+S₂¯k², where FIG. 6. ⌫m as a function of¯␬ as obtained from the “melting

experiments”共symbols兲and the⌫^*= 131 line.

FIG. 7. Static structure functions of the 2D Yukawa liquid for 共a兲⌫^*= 40 and共b兲⌫^*= 10.

S₀=

冉

^2⌫¯␬ ^{+ 1 +Lc}

冊

^{−1 and}

S₂= 2⌫

冉

²␬^{¯⌫+ 1 +Lc}

冊

²

冉

^21¯␬^3−e

冊

^{. 共13兲}

The constante appearing in 共13兲is linked to the M coeffi- cient in共11兲and, as such, cannot be predicted without adopt- ing some model for the dynamics of the strongly coupled systems. If we espouse the extended version of the quasilo- calized charge approximation共QLCA兲that allows the analy- sis of the static response 关37兴, e can be identified as the coefficient of the O共k⁴兲 term in the dynamical matrixD共k兲 关11兴, i.e.,

D共k兲 k² →␻p

2a²共−d+ek¯²兲. 共14兲

Our calculations show thateis small in comparison to the other terms in共13兲, and its precise value does not substan- tially affect the results of calculations ensuing from 共13兲.

Nevertheless, a comparison of theS₀andS₂coefficients共see Fig. 8兲shows that inserting the QLCA value ofeinS₂well reproduces the relationship betweenS₀andS₂as required by 共13兲.

The correlational part of the inverse compressibilityL_cis obtained from the pressure through the relation L_c

=␤共⳵^Pc/⳵n兲. Based on the fitting formula共8兲L_cbecomes

L_c=␤⳵^Pc

⳵^{n =}

冉

^32b0

^⬘

^+b1

^⬘

^¯␬

冊

^⌫=共− 0.8458 + 0.09367␬^¯兲⌫. 共15兲 If, on the other hand, the static structure function is known,L_ccan be determined directly fromS₀=S共k= 0兲 as

L_c= 1 S₀−2⌫

␬

¯ − 1. 共16兲

We note that the difficulty in this latter method is the strong dependence of S共k兲 for small k values on the fluctuations appearing in the simulation due to the finite number of par- ticles. This makes long runs and extensive averaging neces- sary to achieve data with an acceptable level of scattering.

The comparison of the results obtained from the two in- dependent methods provides an important check on the con- sistency and accuracy of the computation. The outcomes of the two calculations are compared in Fig. 9. A strong coin- cidence of the two sets of results, especially for larger ␬^¯ values, verifies the consistency of the computational proce- dure.

Figure 9 shows that the absolute value of the共negative兲L_c decreases with increasing^¯␬values. We note that, in any case, the total compressibility includes the Hartree termL_H, and is always positive, as required by the usual equilibrium stability criterion.

VII. DIELECTRIC FUNCTION

The static dielectric function␧共k兲is an important quantity governing the screening of an impurity and the tendency of the system to develop charge-density waves. It is obtained from S共k兲 via the fluctuation-dissipation relation 共11兲 by identifying

␧共k兲= 1 −␾共k兲␹^¯共k兲, 共17兲 thus providing

␧共k兲= 1 1 − 2⌫ 1

冑

¯k²+^¯␬^2S共k兲

. 共18兲

Figure 10 shows the behavior of the static␧共k兲for various combinations of ⌫ and ␬^¯ values. For high ⌫ values ␧共k兲 exhibits the familiar “inverted U” behavior characteristic for FIG. 8. Ratio of the zeroth- and second-order terms of the static

structure functionS共k→0兲 expansion. Plotted are the values of the S₁/S₀²⌫ expression as functions of⌫ for¯␬= 0.5, 1, 2, and 3. The lines show the predictions of the QLCA theory关the dashed lines show calculations ignoring termein共13兲兴.

FIG. 9. Correlational inverse compressibilityL_Cas a function of

⌫for different␬¯values. Lines show data based on the equation-of- state calculation关through共15兲兴and symbols show points calculated using the structure functionS共k兲 关through共16兲兴. The dashed lines show the theoretical behavior of the pure Coulombic OCP.

HARTMANNet al. PHYSICAL REVIEW E72, 026409共2005兲

strongly coupled Coulomb systems共see, e.g., 关37兴兲, with a negative pole around¯k⬇2⌫. In contrast to Coulomb sys- tems, however,␧共k= 0兲is finite and its value is governed by the value ofL_c.

With increasing ⌫^* values, the maximum of ␧共k兲 in the inverted-U region approaches zero 共see Fig. 10兲. This indi- cates the tendency of the system to develop charge-density waves, i.e., an oscillatory screening response to an impurity 关38兴.

The position of the singularity ¯k

pole of ␧共k¯兲 for low ⌫ values and the negative maximum of␧共k¯兲as a function of⌫^* are displayed in Figs. 10共b兲and 10共c兲. From共11兲and共17兲it is not difficult to show that the pole at¯k=¯k^* disappears and

␧共k¯兲becomes a monotonic positive definite function of¯kfor

⌫values where 1 +L_c艌0.

VIII. VELOCITY AUTOCORRELATION FUNCTION The velocity autocorrelation functionZ共t兲is obtained ac- cording to

Z共t兲=具v共t兲v共0兲典

具兩v共0兲兩²典 , 共19兲

where the average is taken over theNparticles and different initial times.

Representative velocity autocorrelation functions共VACF兲 共19兲obtained at⌫^*= 120 for series of^¯␬values are displayed in Figs. 11共a兲–11共d兲. TheZ共t兲functions共plotted against␻pt兲 exhibit marked oscillations. Here ␻p=

冑

2␲^{Q2n/ma is the} nominal 2D plasma frequency. Such marked oscillations with a period nearly independent of ⌫ in the two-dimensional strongly coupled classical electron layer—similar to those of the three-dimensional one-component plasma关39兴—have al- ready been found by Kaliaet al.关40兴and Hansenet al.关41兴. In the 3D case the oscillations of the VACFs were indeed expected on the basis of the possible coupling between the single particle motion and long-wavelength plasmons whose frequencies are almost independent ofk. In 2D, in contrast, it is well known that the plasmon frequency 共␻兲 depends on the wave numberk and in the limit ofk→0 the mode fre- quency␻→0. The␻共k兲 dispersion curve, however, flattens at higher wave numbers, and this latter behavior has been identified to lead to an appreciable coupling between the single particle motion and collective excitations.

Our results show that the dominant frequency in the Cou- lomb case amounts ␻⬵0.9 ␻p 关Fig. 11共a兲兴, and this value decreases with increasing^¯␬关Figs. 11共b兲–11共d兲兴. Figure 11共e兲 shows theZ共t兲functions replotted against␻Etinstead of␻pt, where␻Eis the Einstein frequency, defined as the oscillation frequency of a test particle in the frozen environment of all other particles situated at lattice sites.␻Eas a function of^¯␬is FIG. 10. 共a兲Dielectric function␧共k兲for⌫= 5,¯␬= 0.5, 1, 2, and

3.共b兲The position of the singularity¯k

poleof␧共¯k兲for low⌫values.

共c兲The negative maximum of␧共¯k兲as a function of⌫^*.

FIG. 11. 共a–d兲 Velocity autocorrelation functions obtained at

⌫^*= 120 and a series of ¯␬ values 共␻p is the nominal 2D plasma frequency兲.共e兲The same data replotted against time normalized by the共␬¯-dependent兲Einstein frequency␻E.

obtained from a formula based on the QLCA method关42兴, using MD-generated pair-correlation functions. Figure 11共e兲 shows that the oscillation frequency共when given in units of

␻E兲 is nearly independent of ^¯␬ and in the vicinity of ␻

⬵1.3␻E, indicating a universality of the dynamical behavior of the systems investigated.

The Fourier transforms¯Z共␻兲of the VACFs shown in Figs.

11共a兲–11共c兲are portrayed in Fig. 12共a兲. The dominant peaks in the spectra, shifting toward lower frequencies with in- creasing^¯␬, correspond to the high-frequency oscillations of theZ共t兲functions共easily observed visually兲. As discussed in previous studies共see, e.g.,关23,43兴兲, these peaks are related to longitudinal current fluctuations, while the broad features at low frequencies are connected to the共slower兲decay ofZ共t兲.

The latter characterize the transverse current fluctuations and are related to diffusion properties of the system. The spectra of the longitudinal and transverse current fluctuations, L共k¯,␻兲andT共k¯,␻兲, respectively, are displayed in the form of color maps in Figs. 12共b兲 and 12共c兲, for the ⌫= 360, ␬^{¯= 2} case. The spectra of the longitudinal current fluctuations show that at small¯k the mode frequency increases linearly with¯k, then within a relatively wide range of¯k the mode frequency in near to␻^/␻p⬇0.45, in correspondence with the peak of¯Z共␻兲 shown in Fig. 12共a兲. The T共k¯,␻兲 spectra 关see Fig. 12共c兲兴 for any¯k are broader, compared to the L共k¯,␻兲 spectra; the fluctuations in the transverse currents are distrib- uted over a rather broad frequency domain, again in agree- ment with the behavior of the corresponding¯Z共␻兲function.

The observed features ofZ¯共␻兲indeed support the model of coupling between single particle motion and collective exci- tations in the 2D system.

IX. SUMMARY

Two-dimensional strongly coupled Yukawa systems have been investigated by molecular-dynamics simulations. We have shown that an effective coupling coefficient⌫^* can be defined, based on the constancy of the first peak amplitude of pair-correlation functions 共PCFs兲. It was demonstrated that systems characterized by different共⌫,^¯␬兲 pairs belonging to the same⌫^*have very nearly the same PCFs. We have cal- culated thermodynamic quantities and have found good agreement with recent theoretical calculations关10兴.

The results of “melting experiments” on systems charac- terized by a wide range of the screening parameter␬^{¯ showed} that the solid-liquid phase boundary can be given by ⌫^*

=⌫m

OCP to a good approximation.

The satisfaction of theoretically established relationships between the static structure function S共k兲 and the correla- tional inverse compressibilityLchas been verified, indicating the high degree of reliability of the computational procedure.

We have analyzed the ensuing behavior of the static dielec- tric function and the predicted screening properties of the system.

The qualitative behavior of the velocity autocorrelation function共VACF兲can be well understood in terms of the dy- namical charge and current fluctuation spectra. The VACF has also been found to exhibit a universal behavior.

ACKNOWLEDGMENTS

This work has been supported by the Hungarian Fund for Scientific Research Grants No. OTKA-T-34156 and OTKA- T-48389, the European Commission Grant No. MERG-CT- 2004-502887, MTA-OTKA-90/46140, NSF Grant No.

PHYS-0206695, DOE Grant No. DE-FG02-03ER5471, and NSF Grant No. PHYS-0514619.

FIG. 12. 共Color online兲 共a兲 Fourier transforms ¯Z共␻兲 of the VACFs obtained at⌫^*= 120, for selected values of¯␬.共b兲longitudi- nalL共¯k,␻兲and共c兲transverseT共¯k,␻兲current fluctuations obtained at⌫= 360,␬¯= 2.共The color coding of the amplitude is logarithmic;

it only intends to illustrate qualitative features.兲

HARTMANNet al. PHYSICAL REVIEW E72, 026409共2005兲

关1兴A. Melzer, V. A. Schweigert, I. V. Schweigert, A. Homann, S.

Peters, and A. Piel, Phys. Rev. E 54, R46共1996兲. J. B. Pieper, J. Goree, and R. A. Quinn, J. Vac. Sci. Technol. 14, 519 共1996兲; M. Zuzic, A. V. Ivlev, J. Goree, G. E. Morfill, H. M.

Thomas, H. Rothermel, U. Konopka, R. Sutterlin, and D. D.

Goldbeck, Phys. Rev. Lett. 85, 4064共2000兲; S. Robertson, A.

A. S. Gulbis, J. Collwell, and M. Horányi, Phys. Plasmas 10, 3874共2003兲; Z. Sternovsky, M. Lampe, S. Robertson, IEEE Trans. Plasma Sci. 32, 632共2004兲.

关2兴E. Grun, H. A. Zook, M. Gabuhlet al., Nature共London兲 362, 428共1993兲; H. Kruger, M. Horányi, and E. Grun, Geophys.

Res. Lett. 30, art. 1058共2003兲; C. E. Krauss, M. Horányi, and S. Robertson, New J. Phys. 5, art. 70共2003兲.

关3兴H. Löwen, J. Phys. A36, 5827共2003兲; A. P. Hynninen and M.

Dijkstra, J. Phys.: Condens. Matter 15, S3557共2003兲; S. Auer and D. Frenkel,ibid. 14, 7667共2002兲.

关4兴C. C. Grimes and G. Adams, Phys. Rev. Lett. 36, 145共1976兲. 关5兴R. C. Gann, S. Chakravarty, and G. V. Chester, Phys. Rev. B

20, 326共1979兲.

关6兴R. K. Kalia, P. Vashishta, and S. W. de Leeuw, Phys. Rev. B 23, 4794共1981兲.

关7兴H. Totsuji, C. Totsuji, and K. Tsuruta, Phys. Rev. E 64, 066402共2001兲.

关8兴S. Nunomura, D. Samsonov, and J. Goree, Phys. Rev. Lett.84, 5141共2000兲; V. Nosenko, S. Nunomura, and J. Goree,ibid.

88, 215002 共2002兲; S. Nunomura, J. Goree, S. Hu, X. Wang, and A. Bhattacharjee, Phys. Rev. E 65, 066402 共2002兲; H.

Thomas, G. E. Morfill, V. Demmel, J. Goree, B. Feuerbacher, and D. Mohlmann, Phys. Rev. Lett. 73, 652 共1994兲; S. K.

Zhdanov, D. Samsonov, and G. E. Morfill, Phys. Rev. E 66, 026411共2002兲.

关9兴M. S. Murillo and D. O. Gericke, J. Phys. A 36, 6273共2003兲. 关10兴H. Totsuji, M. S. Liman, C. Totsuji, and K. Tsuruta, Phys. Rev.

E 70, 016405共2004兲.

关11兴G. J. Kalman, P. Hartmann, Z. Donkó, and M. Rosenberg, Phys. Rev. Lett. 92, 065001共2004兲.

关12兴Z. Donkó, P. Hartmann, and G. J. Kalman, Phys. Rev. E 69, 065401共R兲 共2004兲.

关13兴H. Löwen, J. Phys. C 4, 10105共1992兲.

关14兴G. E. Morfill, H. M. Thomas, U. Konopka, and M. Zuzic, Phys. Plasmas 6, 1769共1999兲.

关15兴F. M. Peeters and X. Wu, Phys. Rev. A 35, 3109共1987兲. 关16兴R. W. Hockney and J. W. Eastwood,Computer Simulation Us-

ing Particles共McGraw-Hill, New York, 1981兲.

关17兴N. David and S. M. Hooker, Phys. Rev. E 68, 056401共2003兲. 关18兴D. Frenkel and B. Smit,Understanding Molecular Dynamics

Simulations共Academic, New York, 2001兲.

关19兴P. Hopkins, A. J. Archer, and R. Evans, Phys. Rev. E 71, 027401共2005兲.

关20兴J.-P. Hansen, I. R. McDonald, and E. L. Pollock, Phys. Rev. A 11, 1025共1975兲.

关21兴S. Hamaguchi, Plasmas Ions 2, 57共1999兲.

关22兴Y. Rosenfeld, J. Phys.: Condens. Matter 13, L39共2001兲. 关23兴H. Ohta and S. Hamaguchi, Phys. Plasmas 7, 4506共2000兲. 关24兴T. Saigo and S. Hamaguchi, Phys. Plasmas 9, 1210共2002兲. 关25兴O. Vaulina, S. Khrapak, and G. Morfill, Phys. Rev. E 66,

016404共2002兲.

关26兴O. Vaulina and S. V. Vladimirov, Phys. Plasmas 9, 835共2002兲. 关27兴V. E. Fortov, O. S. Vaulina, O. F. Petrov, V. I. Molotkov, A. M.

Lipaev,V. M. Torchinsky, H. M. Thomas, G. E. Morfill, S. A.

Khrapak, Yu. P. Semenov, A. I. Ivanov, S. K. Krikalev, A. Yu.

Kalery, S. V. Zaletin, and Yu. P.. Gidzenko, Phys. Rev. Lett.

90, 245005共2003兲.

关28兴G. Faussurier and M. S. Murillo, Phys. Rev. E 67, 046404 共2003兲.

关29兴G. Faussurier, Phys. Rev. E 69, 066402共2004兲. 关30兴H. Ikezi, Phys. Fluids 29, 1764共1986兲.

关31兴S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin, Phys. Rev.

E 56, 4671共1997兲.

关32兴K. I. Golden, G. Kalman, and P. Wyns, Phys. Rev. A 46, 3463 共1992兲.

关33兴S. T. Chui, Phys. Rev. Lett. 48, 933共1982兲; G. Aeppli and R.

Bruinsma,ibid. 53, 2133共1984兲; K. J. Strandburg, Rev. Mod.

Phys.60, 161共1988兲; Ch. Simon, I. Rosenman, F. Batallan, C.

Lartigue, and J. F. Legrand, Phys. Rev. B 45, 2694共1992兲; V.

N. Ryzhov and E. E. Tareyeva,ibid. 51, 8789共1995兲. 关34兴J. M. Kosterlitz and D. J. Thouless, J. Phys. C6, 1181共1973兲. 关35兴B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. 41, 121

共1978兲.

关36兴I. V. Schweigert, V. A. Schweigert, and F. M. Peeters, Phys.

Rev. Lett. 82, 5293共1999兲.

关37兴K. I. Golden, H. Mahassen, and G. J. Kalman, Phys. Rev. E 70, 026406共2004兲.

关38兴S. Kyrkos and G. J. Kalman, J. Phys. A 36, 6235共2003兲. 关39兴J. P. Hansen, I. R. McDonald, and E. L. Pollock, Phys. Rev. A

11, 1025共1975兲.

关40兴R. K. Kalia, P. Vashishta, S. W. de Leeuw, and A. Rahman, J.

Phys. C 14, L991共1981兲.

关41兴J. P. Hansen, D. Levesque, and J. J. Weis, Phys. Rev. Lett. 43, 979共1979兲.

关42兴K. I. Golden and G. J. Kalman, Phys. Plasmas 7, 14共2000兲. 关43兴P. Schmidt, G. Zwicknagel, P.-G. Reinhard, and C. Toepffer,

Phys. Rev. E 56, 7310共1997兲.