Thermodynamics and phase transitions in two-dimensional Yukawa systems

JID:PLA

AID:22865 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.140; Prn:14/10/2014; 13:08] P.1 (1-5)

Physics Letters A ••• (••••) •••–•••

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

Thermodynamics and phase transitions in two-dimensional Yukawa systems O.S. Vaulina a,b , X.G. Koss a,b,∗ a b

Joint Institute for High Temperatures RAS, 125412, Izhorskaya st. 13 bd. 2, Moscow, Russia Moscow Institute of Physics and Technology, 117303, Kerchenskaya st. 1a-1, Moscow, Russia

a r t i c l e

i n f o

Article history: Received 4 June 2014 Received in revised form 18 August 2014 Accepted 1 October 2014 Available online xxxx Communicated by F. Porcelli Keywords: 2D systems Yukawa potential Phase transitions Numerical simulations Hexatic phase

a b s t r a c t The results of numerical simulations of strongly-coupled two-dimensional dissipative Yukawa systems are presented. The thermodynamic characteristics of these systems were studied, namely the internal energy, the specific heat and the entropy. For the first time, it is discovered that the considered characteristics have two singular points on the melting line; one of these points corresponds to the first-order phase transition from crystal to the hexatic phase, and another point corresponds to the second-order phase transition from the hexatic phase to the isotropic liquid. The obtained results are compared to the existing numerical and analytical data. © 2014 Elsevier B.V. All rights reserved.

Thermodynamic properties of the systems of interacting particles are of significant interest in various fields of science and technology, such as plasma physics, medical industry, physics of polymers and others [1–8]. Nowadays the special focus is on the study of the behavior of two-dimensional (2d) systems. The analysis of the physical characteristics of these systems has not only the fundamental scientific value, but is also vital in nano- and microtechnology, as well as for the development of new materials and coatings with tailored properties [5–8]. In thermodynamics of two-dimensional structures, their melting is of the greatest interest. At present, there are two main approaches describing the phase transitions in these systems, based on the analysis of formation of various topological defects. The first is KTHNY (Kosterlitz–Thouless–Halperin–Nelson–Young) theory, which predicts the two-stage transition from the crystal to the liquid state of the system via the intermediate hexatic phase [9–11]. The second is GBI (Grain-Boundary-Induced melting) theory, describing the melting of a 2d-system as the first-order transition from the crystal to the liquid without the formation of the intermediate phase [12,13]. The evidences for the KTHNY theory in 2d-systems with various interaction potentials were found both in numerical and experimental studies [14–23]. So, for example, the detailed analysis of

*

Corresponding author at: Joint Institute for High Temperatures RAS, 125412, Izhorskaya st. 13 bd. 2, Moscow, Russia. Tel.: +7 (495) 484 23 55. E-mail address: [email protected] (X.G. Koss). http://dx.doi.org/10.1016/j.physleta.2014.10.004 0375-9601/© 2014 Elsevier B.V. All rights reserved.

the dynamics of the grains interacting via the screened Coulomb potential (of Yukawa type)

φ(r ) = (e Z )2 exp(−r /λ)/r ,

(1)

where λ is the screening length, e Z is the charge of a grain, is presented in [21–23]. This model is of a special interest for the theoretical study of dusty plasma [5–8], and also is widely used for the numerical simulation of repulsion in the kinetics of interacting particles [1–8]. The numerical experiments in dissipative Yukawa systems with νfr = 0 (where νfr is the friction coefficient of particles due to their collisions with the neutrals of surrounding gas) show that the physical characteristics (e.g. the maximum/minimum ratio of the pair correlation functions, the isothermal compressibility, the heat capacities, and the diffusion constants) of these systems have two singularities [21–23]. The first corresponds to the “liquid-hexatic” phase transition and is observed when the effective coupling parameter Γ ∗ = Γh∗ ≈ 98 ± 3; the second singular point (when Γ ∗ = Γc∗ ≈ 154 ± 4) is due to the transition from the hexatic phase to the “ideal” crystal, where the diffusion coefficient of the particles D → 0. Here Γ ∗ = 1.5 (e Z )2 (1 + κ + κ 2 /2) exp(−κ )/( T rp ), where κ = rp /λ, T is the temperature of the particles measured in energy units, and rp is the mean interparticle distance. Note that Clark et al. [24] obtained similar values of Γ ∗ on the phase transitions lines for the Coulomb systems (κ = 0) in case of νfr = 0 in classical limit of the problem, using the quantum Monte Carlo method.

JID:PLA

2

AID:22865 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.140; Prn:14/10/2014; 13:08] P.2 (1-5)

O.S. Vaulina, X.G. Koss / Physics Letters A ••• (••••) •••–•••

Table 1 The energy, U 0 (Γ T )−1 , for hp-lattices in Yukawa systems with various screening parameters, κ .

Fig. 1. The normalized values of the coefficients of diffusion D ∗ and viscosity ν ∗ for quasi-2d Yukawa systems [18] (solid lines). The symbols are the normalized data: (P) – for 2d-colloidal systems [25]; (!) – for 2d-disperse systems (νfr = 0, ξ → ∞, κ = 0.56) [33].

We must emphasize that one should perform the calculation of physical characteristics of 2d-systems (including spatial correlation functions) near the phase transition points and for the conditions corresponding to the hexatic phase with especial thoroughness, and the time of the simulation should be controlled for the results to be consistent [6,8,25]. This is particularly important for the systems with the small friction (νfr → 0). Thus, in [25] it was found that for the pure dispersive Yukawa systems with νfr = 0 the hexatic phase is metastable and disappears in the limit of long times. The most laboratory studies of the properties of two-dimensional systems is carried out in the colloids or in the weakly ionized gas discharge plasma, where the dissipation caused by the collisions of the charged grains with the neutral particles of the surrounding medium (νfr = 0) can substantially affect the dynamics and the conditions of formation of the ordered dust structures. We should note that, unlike the experiments with colloidal systems, demonstrating the good agreement with the statements of KTHNY theory [16–19], the experiments with dusty structures in the laboratory radio-frequency discharge did not show the compelling evidence either of KTHNY or of GBI theories [26–29]. Nevertheless, the authors of the latest works [28,29] tend towards the GBI theory describing the melting scenario of the two-dimensional dusty structures forming in RF-discharge plasma. It is usually assumed that in case of KTHNY scenario, both phase transitions in the system are of the second order [13,30]. Nevertheless, the results of some numerical and theoretical studies show that, depending on the type of pair interaction potential, both transitions could be of the first order, or one of them can be of the first, and the other – of the second order [24,25,31,32]. For the qualitative and quantitative analysis of the phase state of 2d-systems, the spatial correlation functions are usually used (pair g (r ) and orientational g 6 (r )) [4–11,23,25,33]. It should be noted that in the work [23] the behavior of g 6 (r ) function is presented, including the vicinity of the points of phase transitions. Also it is useful to study the behavior of various dynamical and structural properties of the system [21,22], for example, the diffusion and viscosity coefficients, which have the distinct singularity at Γ ∗ close to 98, which is much lower than the point of crystallization (Γc∗ ≈ 154), where D → 0, see Fig. 1 [21,34,35].

κ

U 0 (Γ T )−1

1 1.5 2 3 4

1.62532 7.59943 × 10−1 3.93543 × 10−1 1.19909 × 10−1 3.91733 × 10−2

Nevertheless, the analysis of these functions does not help to find out the character of the phase transitions observed. To know the order of the phase transition we should explore the behavior of various thermodynamic functions and characteristics, such as the specific entropy, the internal energy, the specific heat, the compressibility, etc. Thermodynamic characteristics of quasi-twodimensional Yukawa systems (taking into account the displacements of particles perpendicular to the layer) close to the phase transitions were studied numerically in [22,36]. These works have confirmed the existence of two singularities near Γ ∗ = 98 and Γ ∗ = 154. However, the processed numerical data were not averaged, and were essentially discrete near the points of phase transitions, which is the reason the order of the observed phase transitions wasn’t determined. In the present work we show the results of the numerical study of thermodynamic properties of strongly coupled liquid and crystalline 2d-systems of particles which interact via the screened Coulomb (Yukawa) potential. The coupling parameter Γ ∗ was varied within the wide range, and its discretization interval was chosen to be small enough to disambiguate the behavior of the system close to the points of phase transitions. In case of isotropic pair interactions (with the known interaction energy φ ≡ φ(r )) the physical properties of coupled systems, such as the energy density U , are determined by the particle temperature T , the concentration, n, and the pair correlation function, g (r ), which can be measured experimentally or may be found from the computer simulations [4,36,37]:

U=

m 2

∞ T + (m − 1)π n

φ(r ) g (r )rm−1 dr ,

(2)

0

where m = 2, 3 is the number of dimensions in the system, and n = rp−m . For a system with a constant number of particles, a lot of thermodynamic characteristics can be obtained from the calorific equation of state U ( T , n, φ, g ) and the basic formulas of thermodynamics – for example, the heat capacity, C V = (∂ U /∂ T ) V , where V is the volume/surface area of the system, the pressure, the entropy, etc. [36]. In case of T → 0 the value of U tends to U 0 – the energy density for the crystal lattice at T = 0. For any lattice of the known type the values of U 0 may be easily computed [36]. So, for the classical triangular lattice (hp-crystal) these values are as follows:

U0 =

1  2

φ(ri j ).

(3)

i =1 j =1



√

Here r i j = a i 2 + j 2 + i j, and a = ahp ≡ rp (2/ 3 )1/2 is the step (the lattice spacing) of hp-lattice. The normalized values of energy, U 0 (Γ T )−1 , for hp-lattices in Yukawa systems with various screening parameters, κ , used in the numerical experiment, are shown in Table 1. The numerical simulation was carried out by the Langevin molecular dynamics method (using the periodic boundary conditions in two chosen directions), based on the solution of the system of N p ordinary differential equations, N p being the number

JID:PLA

AID:22865 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.140; Prn:14/10/2014; 13:08] P.3 (1-5)

O.S. Vaulina, X.G. Koss / Physics Letters A ••• (••••) •••–•••

Fig. 2. The maximum g max of the pair correlation function dependent on Γ ∗ for: (solid line) – 2d-systems (the results of the present work); (!) – quasi-2d-systems [21]. The presented error bars correspond to 4%.

of the independent particles in computational cell. The stochastic character of motion of the particles with the given kinetic temperature T was determined by the Langevin force F ran . The simulation technique is detailed in [6]. The calculations were carried out for the homogeneous 2d-Yukawa systems with the values of screening parameter κ ≡ rp /λ = 1, 1.5, 2, 3, 4 (where rp = ( N p / S )1/2 , and S is the area of the computational cell). The number of independent particles in the central computational cell N p was varied from 256 to 4096. Dependent on the number of particles, the cutoff length of potential r cut changed from 5rp to 25rp . The majority of data were obtained for N p = 1024 independent particles and rcut = 12rp , because the further increase of the number of particles did not lead to the substantial change of the numerical results. For the analysis of thermodynamic properties of the systems under study the equations of motion were solved for various values of effective parameters: the coupling parameter Γ ∗ and the scaling parameter ξ = ω∗ /νfr , where ω∗ = {2(e Z )2 (1 + κ + κ 2 /2) exp(−κ )/(rp3 π M )}1/2 , M is the mass of the particle. The value of the effective coupling parameter Γ ∗ was varied in the range from ∼5 to 275 (with the step Γ ∗ ≈ 2.5), and the value of the scaling parameter – from ξ ≈ 0.25 to ξ ≈ 4 (ξ ≈ 0.25; 0.5; 1; 2; 4), the parameters typical for the experiments in gas-discharge plasmas [5–8]. The internal energy was calculated by the expression (2) based on the computations of the pair correlation functions, g (r ). The calculations of the g (r ) functions were carried out after the system had reached the equilibrium condition, which was controlled by the comparison of the velocity distribution function with the Maxwell function for a given temperature. The duration of the calculations t sim of the g (r ) functions in equilibrium systems was ∼1000/ min{νfr , ω∗ }. For the chosen t sim the results of the calculations did not depend on the time of numerical simulation for all the values of Γ ∗ , including the range corresponding to the hexatic phase. For all cases under study the form of obtained pair correlation functions, g (r ), was determined by the value of Γ ∗ parameter for the systems with Γ ∗ > 10–15. (The disagreement between the data for all calculations was within the limits of the numerical error and did not exceed ±(1–3)%.) Similar results for 3d-systems with isotropic pair potentials and in quasi-two-dimensional structures were obtained in [5–8,21,36] for the wide ranges of values of their parameters. The results of the calculation of the maximum g max of the pair correlation function dependent on Γ ∗ are shown in Fig. 2. Deviations of the presented calculations of g max from the data, obtained earlier for the quasi-2d-systems, are most likely de-

3

Fig. 3. The values of δ U vs. Γ ∗ for 2d-Yukawa systems with various κ : (×) – 1; (!) – 1.5; (P) – 2; (E) – 3; (1) – 4. The line is the approximation [36]. The presented error bars correspond to 3%. (Here we have omitted some calculated results for the figure to be more demonstrative, see the text.)

Fig. 4. The averaged values of δ U (bold solid line) vs. Γ ∗ for 2d-Yukawa systems with various κ = 1–4.

termined by the spatial partitioning of the analyzed area of the structure, which is more fine in the present work. The calculations of thermodynamic characteristics (U , C V , S) vs. Γ ∗ , using the obtained g (r ) functions, are presented in Figs. 3–6 for various screening parameters, κ . The illustrations of normalized functions U (Γ ∗ ) are presented in Fig. 3. It can be easily seen that the obtained normalized values of δ U = (U − U 0 − mT /2)/ T are completely determined by the value of Γ ∗ parameter in accordance with the semi-empirical approach, suggested in [36] for the coupled liquid-like systems with Γ ∗ lying in the range Γc∗ > Γ ∗ > 10. The comparison of the averaged (by all κ and ξ ) functions δ U (Γ ∗ ) with the averaged data obtained for the quasi-2d-structures [36], is shown in Fig. 4. It can be easily noticed that, despite the good agreement of our data with the approximation [36], the numerical curves have the singularities close to Γ ∗ = 98 (point of inflexion) and close to Γ ∗ = 154 (jump), which aren’t reflected in the behavior of the analytical curve, see Figs. 3, 4. Remember, that in the systems with a constant volume the criterion of the first-order phase transition is the jump of the specific entropy of the system, and the criterion of the second order phase

JID:PLA

AID:22865 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.140; Prn:14/10/2014; 13:08] P.4 (1-5)

O.S. Vaulina, X.G. Koss / Physics Letters A ••• (••••) •••–•••

4

Fig. 5. The (C V − 1) values vs. Γ ∗ for 2d-Yukawa systems (solid line), obtained from the averaged δ U values. Dotted curves are the lines of phase transitions.

Fig. 6. The excess entropy  S vs. Γ ∗ for 2d-Yukawa systems, obtained from the averaged δ U values.

transition is the jump of its derivatives with respect to the temperIn this case the entropy can be presented ature and the pressure.   in the form S = dT C V / T + const ≡ U / T − U /( T Γ ∗ )dΓ ∗ + const, and its absolute variation  S = S (Γ ∗ )–S (0) – in the form

Γ S = δU −

∗

δU

dΓ ∗

Γ∗

(4)

.

0



We have

∂S ∂T





= V

∂( S ) ∂T

 ≡ V

CV T

to notice two singularities on the curve C V (Γ ∗ ) near Γ ∗ = 98 (the jump) and Γ ∗ = 154 (the fast growth). The results of the calculations of the excess entropy  S are presented in Fig. 6. The curve  S (Γ ∗ ) was obtained with the help of Eq. (4) and the extrapolation of the numerical data in the range Γ ∗ < 5. Like δ U (Γ ∗ ) curve, it has two singular points denoting the phase transitions, as it was predicted by Eq. (4). Note that the finiteness of the peak on C V curve near the solidhexatic transition and the difference between the observed jump of the C V value near the fluid-hexatic transition from the classic Heaviside step function, as well as the smoothed features of the curve  S (Γ ∗ ), are obviously caused by the mathematical peculiarities of integration/differentiation of discontinuous data. The rougher is the discretization of analyzed data, the wider is the observed peak, the smaller is its amplitude, and the singularities obtained by the integration become less distinctive. Let us take a look at the results of the recent experiments with 2d-colloidal systems of micrometer-sized superparamagnetic particles [38] and the numerical data obtained by Monte Carlo simulations for the disperse systems (νfr = 0, ξ → ∞) with the pair interaction of particles [30]. They demonstrate the substantial growth of the specific heat near the solid-hexatic transformation; however, the authors didn’t count this fact as the sign of the firstorder phase transition. Note also that the authors of the mentioned works didn’t identify the second singularity of the specific heat corresponding to the transition from the hexatic to the isotropic liquid, though they have presented the convincing proofs of the existence of the intermediate hexatic phase. Finally, here we have presented new results of numerical investigation of thermodynamic functions and characteristics, such as the internal energy, the entropy and the specific heat, for strongly coupled liquid and crystal systems of the particles interacting via the Yukawa potential. The simulations were carried out for the two-dimensional structures within the wide range of parameters corresponding to the conditions of experiments in the laboratory dusty plasma (for the systems with the effective coupling parameter Γ ∗ > 5). We have found that the mentioned thermodynamic characteristics have two singularities within the melting of the considered systems near Γ ∗ = 98 and Γ ∗ = 154. The behavior of the thermodynamic characteristics close to these points shows that the phase transitions observed in their proximity are of the second and of the first order, respectively. The reported results can be considered as one more evidence of the two-stage melting scenario of two-dimensional Yukawa systems. The comparison of the obtained data to the existing experimental, numerical and analytical data is presented. The results of the present work may be easily adapted for the systems with the wide range of isotropic pair potentials that are interesting in plasma physics as well as in the medicine, biology, physics of polymers and other colloidal systems. Acknowledgements

.

(5)

Therefore, the jump on the δ U (Γ ∗ ) curve leads to the jump of the entropy, which is the sign of the first-order phase transition. Note that the specific heat C V of the system should tend to infinity in the ideal case. When the curve δ U (Γ ∗ ) has the point of inflexion, the dependence of the entropy S on Γ ∗ also has the inflexion point, which leads to the jump of the derivatives of these functions on temperature (that is, the jump of C V ) and is the sign of the second-order phase transition. The specific heat, C V = (∂ U /∂ T ) V , obtained by the differentiation of the averaged data on δ U (Γ ∗ ), is shown in Fig. 5. It is easy

This work was partially supported by the Russian Foundation for Basic Research (grants No. 13-08-00263-a, No. 14-08-31633mol-a), by the Ministry of Education and Science of the Russian Federation (Presidential scholarship for X.G. Koss) and by the Program of the Presidium of RAS. References [1] H.Z. Cummins, E.R. Pike (Eds.), Photon Correlation and Light Beating Spectroscopy, Plenum, New York, 1974. [2] B. Pullman (Ed.), Intermolecular Interactions: From Diatomics to Biopolymers, Wiley Interscience, Chichester, 1978.

JID:PLA

AID:22865 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.140; Prn:14/10/2014; 13:08] P.5 (1-5)

O.S. Vaulina, X.G. Koss / Physics Letters A ••• (••••) •••–•••

[3] A.A. Ovchinnikov, S.F. Timashev, A.A. Belyy, Kinetics of Diffusion Controlled Chemical Processes, Nova Science Publishers, Commack, NY, 1989. [4] N.H. March, M.P. Tosi, Introduction to Liquid State Physics, World Scientific, 1995. [5] S.V. Vladimirov, K. Ostrikov, A.A. Samarian, Physics and Applications of Complex Plasmas, Imperial College, London, 2005. [6] V.E. Fortov, G.E. Morfill (Eds.), Complex and Dusty Plasmas, CRC Press, 2010. [7] A. Ivlev, H. Lowen, G. Morfill, C.P. Royall, Complex Plasmas and Colloidal Dispersions: Particle-Resolved Studies of Classical Liquids and Solids, World Scientific, Singapore, 2012. [8] B.A. Klumov, Phys. Usp. 53 (2010) 1053. [9] D.R. Nelson, B.I. Halperin, Phys. Rev. B 19 (1979) 2457. [10] J.M. Kosterlitz, D.J. Thouless, J. Phys. C 6 (1973) 1181. [11] A.P. Young, Phys. Rev. B 19 (1979) 1855. [12] S.T. Chui, Phys. Rev. B 28 (1983) 178. [13] K. Strandburg, Rev. Mod. Phys. 60 (1988) 161. [14] A. Jaster, Phys. Rev. E 59 (1999) 2594. [15] D.C. Glattli, et al., Phys. Rev. Lett. 60 (1988) 420. [16] C.A. Murray, R.A. Wenk, Phys. Rev. Lett. 62 (1989) 1643. [17] A.H. Marcus, S.A. Rice, Phys. Rev. Lett. 77 (1996) 2577. [18] R.E. Kusner, et al., Phys. Rev. Lett. 73 (1994) 3113. [19] K.J. Naidoo, J. Schnitker, J. Chem. Phys. 100 (1994) 3114;

[20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38]

5

W.K. Qi, et al., J. Chem. Phys. 133 (2010) 234508; A.H. Marcus, S.A. Rice, Phys. Rev. E 55 (1997) 637. K. Zahn, G. Maret, Phys. Rev. Lett. 85 (2000) 3656. O.S. Vaulina, I.E. Drangevski, Phys. Scr. 73 (2006) 577. O.S. Vaulina, et al., Phys. Rev. Lett. 97 (2006) 195001. O.S. Vaulina, E.V. Vasilieva, Phys. Lett. A 378 (2014) 719–722. B.K. Clark, M. Casula, D.M. Ceperley, Phys. Rev. Lett. 103 (2009) 055701. A. Derzsi, A.Zs. Kovacs, Z. Donko, P. Hartmann, Phys. Plasmas 21 (2014) 023706. R.A. Quinn, C. Cui, J. Goree, J.B. Pieper, Phys. Rev. E 53 (1996) 2049. A. Melzer, et al., Phys. Rev. E 53 (1996) 2757. C.A. Knapek, D. Samsonov, Zhdanov, et al., Phys. Rev. Lett. 98 (2007) 015004. V. Nosenko, S.K. Zhdanov, A.V. Ivlev, et al., Phys. Rev. Lett. 103 (2009) 015001. M. Mazars, arXiv:1301.1571 [cond-mat.stat-mech], 2013. P. Bladon, D. Frenkel, Phys. Rev. Lett. 74 (1995) 2519. T. Chou, D.R. Nelson, Phys. Rev. E 53 (1996) 2560. P. Hartmann, G.J. Kalman, Z. Donko, J. Phys. A, Math. Gen. 39 (2006) 4485–4491. H. Lowen, J. Phys. Condens. Matter 4 (1992) 10105. B. Liu, J. Goree, Phys. Rev. Lett. 94 (2005) 185002. O.S. Vaulina, X.G. Koss, Yu.V. Khrustalyov, O.F. Petrov, V.E. Fortov, Phys. Rev. E 82 (2010) 056411. N.K. Ailawadi, Phys. Rep. 57 (1980) 241. S. Deutschlander, A.M. Puertas, G. Maret, P. Keim, arXiv:1312.5511v1, Dec 2013.