Thermodynamics of two-dimensional Yukawa systems across coupling regimes

TL;DR

This study systematically investigates the thermodynamics of two-dimensional Yukawa systems across various coupling regimes using molecular dynamics simulations, providing practical formulas and insights relevant to colloids, plasmas, and electrolytes.

Contribution

It offers comprehensive thermodynamic data and simple models for 2D Yukawa systems across coupling regimes, including crystalline states, using MD simulations and interpolation methods.

Findings

Accurate fits for internal energy and pressure across regimes

Application of shortest-graph interpolation for pair correlations

Finite-temperature effects modeled by peak corrections in pair correlation functions

Abstract

Thermodynamics of two-dimensional Yukawa (screened Coulomb or Debye-H\"uckel) systems is studied systematically using molecular dynamics (MD) simulations. Simulations cover very broad parameter range spanning from weakly coupled gaseous states to strongly coupled fluid and crystalline states. Important thermodynamic quantities such as internal energy and pressure are obtained and accurate physically motivated fits are proposed. This allows us to put forward simple practical expressions to describe thermodynamic properties of two-dimensional Yukawa systems. For crystals, in addition to numerical simulations, the recently developed shortest-graph interpolation method is applied to describe pair correlations and hence thermodynamic properties. It is shown that the finite-temperature effects can be accounted for by using simple correction of peaks in the pair correlation function. The corresponding correction coefficients are evaluated using MD simulation. The relevance of the obtained results in the context of colloidal systems, complex (dusty) plasmas, ions absorbed to interfaces in electrolytes is pointed out.

Tables

Table 1. Table 1: Madelung constants of the 2D Yukawa crystals (triangular lattice) for various screening parameters in the range 0.5 ≤ κ ≤ 3.0 0.5 𝜅 3.0 0.5\leq\kappa\leq 3.0

$\kappa$	$M$	$\kappa$	$M$
0.5	1.11914	1.8	0.05449
0.6	0.82503	2.0	0.03660
0.8	0.48127	2.2	0.02470
1.0	0.29709	2.4	0.01672
1.2	0.18960	2.6	0.01135
1.4	0.12357	2.8	0.00772
1.6	0.08167	3.0	0.00525

Table 2. Table 2: Reduced excess energy u ex subscript 𝑢 ex u_{\rm ex} and pressure p 𝑝 p of two-dimensional Yukawa fluids evaluated using MD simulations for various coupling ( Γ Γ \Gamma ) and screening ( κ 𝜅 \kappa ) parameters.

$\kappa =0.5$
$\Gamma$	135.420	86.7254	52.7787	32.1811	19.6073	11.9310	7.27175	4.43126	2.69848	1.64302	1.00136	0.5
$u_{\mathrm{ex}}$	152.944	98.3115	60.1901	37.0087	22.8180	14.1176	8.79838	5.51964	3.48587	2.21772	1.42021	0.76495
$p$	199.434	128.303	78.6946	48.5651	30.1485	18.8835	12.0216	7.81631	5.22964	3.63556	2.64961	1.85883
$\kappa =0.6$
$\Gamma$	140.131	89.5076	54.3171	32.9737	20.0017	12.1359	7.36665	4.47442	2.71053	1.64677	1.00106	0.5
$u_{\mathrm{ex}}$	116.984	75.1128	45.9415	28.2016	17.3768	10.7727	6.73045	4.24422	2.69421	1.72956	1.11776	0.61083
$p$	160.369	103.050	63.1652	38.9451	24.1971	15.2284	9.76528	6.42899	4.37128	3.11015	2.32663	1.69701
$\kappa =0.8$
$\Gamma$	152.277	96.5736	58.0604	34.9737	21.0334	12.6675	7.61503	4.58845	2.75830	1.66410	0.99914	0.5
$u_{\mathrm{ex}}$	74.6424	47.7340	29.0608	17.8181	10.9844	6.84185	4.30139	2.74217	1.76665	1.15293	0.75437	0.42469
$p$	112.709	72.1411	44.0441	27.1658	16.9406	10.7731	7.01845	4.73986	3.33679	2.47393	1.92983	1.49910
$\kappa =1.0$
$\Gamma$	169.071	105.975	63.1038	37.6027	22.4047	13.3361	7.94729	4.73129	2.81940	1.68034	0.99956	0.5
$u_{\mathrm{ex}}$	51.5786	32.7335	19.8556	12.1451	7.50279	4.68984	2.97702	1.91799	1.25426	0.82932	0.55059	0.31770
$p$	85.4036	54.2492	33.0215	20.3527	12.7618	8.19406	5.44279	3.76791	2.74103	2.10336	1.70075	1.38135
$\kappa =1.2$
$\Gamma$	191.126	118.398	69.6429	40.9597	24.1083	14.1893	8.34919	4.90490	2.88868	1.70019	0.99984	0.5
$u_{\mathrm{ex}}$	37.5852	23.6918	14.3026	8.72609	5.39936	3.39637	2.17547	1.41736	0.93933	0.62908	0.42281	0.24960
$p$	67.9344	42.8619	25.9838	16.0024	10.0874	6.56025	4.44041	3.15023	2.36021	1.86635	1.55301	1.30594
$\kappa =1.4$
$\Gamma$	220.172	134.441	77.9949	45.2452	26.2578	15.2219	8.83634	5.12702	2.97137	1.72440	1.00140	0.5
$u_{\mathrm{ex}}$	28.5555	17.8503	10.7244	6.53392	4.05300	2.56405	1.65932	1.09552	0.73364	0.49726	0.33718	0.20253
$p$	56.0915	35.0963	21.1892	13.0574	8.28303	5.45288	3.76392	2.73780	2.10241	1.70540	1.45171	1.25396
$\kappa =1.6$
$\Gamma$	258.433	155.296	88.6297	50.6106	28.9099	16.4928	9.41249	5.37870	3.07317	1.75217	0.99889	0.5
$u_{\mathrm{ex}}$	22.4535	13.9136	8.31218	5.05719	3.14728	2.00498	1.30903	0.87473	0.59391	0.40446	0.27520	0.16486
$p$	47.7294	29.6021	17.7849	10.9674	7.00739	4.67522	3.28559	2.44432	1.92230	1.58965	1.37647	1.15781
$\kappa =1.8$
$\Gamma$	308.935	182.395	102.261	57.3435	32.1483	18.0355	10.1029	5.67241	3.17978	1.78359	0.99997	0.5
$u_{\mathrm{ex}}$	18.1745	11.1626	6.63304	4.02868	2.51560	1.61389	1.06328	0.71747	0.49051	0.33739	0.23058	0.14359
$p$	41.6428	25.5932	15.3055	9.44338	6.07675	4.10949	2.93845	2.22906	1.78546	1.50402	1.32125	1.18748
$\kappa =2.0$
$\Gamma$	375.818	217.422	119.600	65.7745	36.1611	19.8980	10.9232	6.01199	3.30681	1.81767	1.00051	0.5
$u_{\mathrm{ex}}$	15.0964	9.17319	5.42177	3.29200	2.06139	1.33276	0.88426	0.60261	0.41513	0.28650	0.19651	0.12379
$p$	37.1333	22.5775	13.4413	8.30684	5.38337	3.68921	2.67835	2.06727	1.68347	1.43752	1.27850	1.16494
$\kappa =2.2$
$\Gamma$	463.975	262.948	141.568	76.2338	41.0173	22.0958	11.9035	6.41082	3.45056	1.85303	1.00113	0.5
$u_{\mathrm{ex}}$	12.7875	7.69994	4.52708	2.74830	1.72461	1.12217	0.75368	0.51642	0.35777	0.24734	0.17009	0.10850
$p$	33.6575	20.2710	12.0118	7.43585	4.85060	3.36425	2.48426	1.94445	1.60450	1.38520	1.24473	1.14408
$\kappa =2.4$
$\Gamma$	578.968	320.871	168.949	89.0382	46.8778	24.7092	12.9953	6.85634	3.60307	1.89919	0.99952	0.5
$u_{\mathrm{ex}}$	10.9709	6.56430	3.83850	2.33031	1.47100	0.96365	0.65089	0.44974	0.31141	0.21697	0.14862	0.09589
$p$	30.8215	18.4175	10.8648	6.74135	4.43655	3.11369	2.32748	1.84722	1.53931	1.34446	1.21673	1.12942
$\kappa =2.6$
$\Gamma$	723.656	392.384	202.051	104.080	53.5742	27.6270	14.2191	7.32182	3.76653	1.93971	1.00200	0.5
$u_{\mathrm{ex}}$	9.50055	5.63818	3.28596	1.99866	1.26783	0.83500	0.56905	0.39442	0.27600	0.19145	0.13130	0.08576
$p$	28.3633	16.8096	9.89231	6.16190	4.09049	2.90245	2.19936	1.76426	1.48858	1.30961	1.19408	1.11954
$\kappa =2.8$
$\Gamma$	893.746	474.549	239.143	120.685	60.8483	30.6642	15.4796	7.80951	3.93161	1.98042	1.00296	0.5
$u_{\mathrm{ex}}$	8.19448	4.82859	2.81518	1.71951	1.09985	0.73051	0.50093	0.35038	0.24489	0.17117	0.11700	0.07671
$p$	25.9004	15.2521	8.98792	5.63831	3.78782	2.72194	2.08856	1.69631	1.44344	1.28133	1.17497	1.10201
$\kappa =3.0$
$\Gamma$	1071.02	558.495	276.444	136.953	67.7922	33.5897	16.6383	8.22716	4.07874	2.02013	0.99949	0.5
$u_{\mathrm{ex}}$	6.93189	4.07091	2.38838	1.47193	0.95056	0.64023	0.44340	0.31146	0.21994	0.15395	0.10494	0.06958
$p$	23.1181	13.5906	8.07317	5.12679	3.49444	2.55590	1.98879	1.63334	1.40554	1.25677	1.15868	1.09682

Table 3. Table 3: Values of the anharmonic correction coefficient β 𝛽 \beta for different screening parameter κ 𝜅 \kappa .

$\kappa$	0.0	0.2	0.3	0.4	0.6	0.8	1.0	1.2	1.4	1.6	1.8	2.0	2.2	2.4	2.6	2.8	3.0
$\beta \left(\kappa \right)$	3.01	9.23	12.38	14.30	10.53	9.71	9.35	9.28	9.14	9.08	8.97	8.855	8.68	8.71	8.46	8.47	8.51

Table 4. Table 4: Reduced excess energy u ex subscript 𝑢 ex u_{\mathrm{ex}} of the 2D Yukawa crystal obtained in MD simulations for various screening parameters κ 𝜅 \kappa and reduced coupling parameters Γ m / Γ subscript Γ m Γ \Gamma_{\rm m}/\Gamma .

$\kappa$	${\Gamma }_m/\Gamma$
	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.5	1595.62	798.828	532.689	399.681	319.981	266.796	228.880	200.332	178.283
0.6	1217.36	609.282	406.628	305.117	244.469	203.938	174.914	153.188	136.267
0.8	773.025	387.104	258.328	194.074	155.484	129.733	111.343	97.5607	86.8364
1.0	529.643	265.306	177.235	133.215	106.726	89.1490	76.5169	67.1314	59.7831
1.2	382.522	191.740	128.152	96.3972	77.2970	64.6022	55.5318	48.7317	43.4438
1.4	287.408	144.232	96.4804	72.5942	58.2862	48.7586	41.9386	36.8484	32.8838
1.6	223.185	112.096	75.0671	56.5515	45.4466	38.0606	32.7681	28.8120	25.7391
1.8	178.133	89.6228	60.0889	45.3116	36.4631	30.5563	26.3521	23.1896	20.7451
2.0	145.774	73.3800	49.2712	37.2003	29.9641	25.1447	21.7011	19.1314	17.1275
2.2	121.609	61.3067	41.2021	31.1620	25.1352	21.1177	18.2517	16.1113	14.4385
2.4	102.908	51.9465	34.9672	26.4819	21.3920	17.9999	15.5706	13.7650	12.3602
2.6	87.4157	44.2324	29.8212	22.6181	18.2990	15.4212	13.3710	11.8300	10.6351
2.8	73.5771	37.3025	25.2028	19.1490	15.5271	13.1108	11.3865	10.0997	9.10597
3.0	60.2002	30.6118	20.7457	15.8118	12.8497	10.8840	9.47465	8.43053	7.65187

Table 5. Table 5: Reduced pressure (compressibility) p 𝑝 p of the 2D Yukawa crystal obtained in MD simulations for various screening parameters κ 𝜅 \kappa and reduced coupling parameters Γ m / Γ subscript Γ m Γ \Gamma_{\rm m}/\Gamma .

$\kappa$	${\Gamma }_m/\Gamma$
	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.5	2080.63	1041.70	694.789	521.370	417.454	348.100	298.669	261.442	232.679
0.6	1669.06	835.485	557.680	418.523	335.380	279.814	240.022	210.233	187.022
0.8	1168.03	585.024	390.480	293.406	235.104	196.197	168.410	147.583	131.370
1.0	878.208	440.005	294.000	221.023	177.106	147.964	127.016	111.450	99.2542
1.2	693.046	347.470	232.288	174.765	140.162	117.162	100.726	88.4011	78.8053
1.4	566.555	284.386	190.275	143.196	114.994	96.2113	82.7636	72.7234	64.8975
1.6	476.692	239.477	160.406	120.865	97.1465	81.3696	70.0608	61.6053	55.0288
1.8	410.580	206.621	138.561	104.505	84.1086	70.4915	60.7970	53.5005	47.8555
2.0	361.191	181.859	122.134	92.2267	74.2973	62.3524	53.8144	47.4405	42.4641
2.2	322.729	162.732	109.386	82.7430	66.7485	56.0825	48.4703	42.7821	38.3327
2.4	291.498	147.173	99.0847	75.0489	60.6307	51.0175	44.1300	39.0087	35.0158
2.6	263.437	133.325	89.9002	68.1935	55.1747	46.4976	40.3128	35.6615	32.0486
2.8	235.188	119.260	80.5872	61.2342	49.6540	41.9257	36.4074	32.2829	29.0897
3.0	203.533	103.516	70.1601	53.4777	43.4588	36.8063	32.0351	28.4887	25.8081

Equations

$$\phi (r)=\frac{\epsilon \lambda }{r}\exp \left(-\frac{r}{\lambda }\right),$$

$${\Gamma }_{\mathrm{m}}(\kappa )\simeq \frac{131}{1-0.388{\kappa }^2+0.138{\kappa }^3-0.0138{\kappa }^4}.$$

$$U=-T^2{\left(\frac{\partial }{\partial T}\frac{F}{T}\right)}_V,$$

$$P=-{\left(\frac{\partial F}{\partial V}\right)}_T.$$

$$\begin{array}{rl}& U=N\left(k_{\mathrm{B}}T+n\int d\mathbf{r}\phi (r)g(\mathbf{r})\right),\\ & PV=N\left(k_{\mathrm{B}}T-\frac{n}{4}\int d\mathbf{r}r{\phi }^{\prime }(r)g(\mathbf{r})\right),\end{array}$$

$$p=1+\frac{\Gamma }{2}\frac{\partial f_{\mathrm{ex}}}{\partial \Gamma }-\frac{\kappa }{2}\frac{\partial f_{\mathrm{ex}}}{\partial \kappa },f_{\mathrm{ex}}={\int }_0^{\Gamma }d{\Gamma }^{\prime }\frac{u_{\mathrm{ex}}(\kappa ,{\Gamma }^{\prime })}{{\Gamma }^{\prime }}.$$

$$g(\mathbf{r})=\frac{1}{n}\sum _{\alpha }{p}_{\alpha }(\mathbf{r}-{\mathbf{r}}_{\alpha }),$$

$$\begin{array}{rl}& p_{\alpha }(\mathbf{r})\propto \exp [-\frac{\phi (\mathbf{r}+{\mathbf{r}}_{\alpha })}{k_{\mathrm{B}}T}-b_{\alpha }({\mathbf{e}}_{\alpha }\cdot \mathbf{r})-\\ & -\frac{({\mathbf{e}}_{\alpha }\cdot \mathbf{r}{)}^2}{2a_{\Vert \alpha }^2}-\frac{{\mathbf{r}}^2-({\mathbf{e}}_{\alpha }\cdot \mathbf{r}{)}^2}{2a_{\perp \alpha }^2}].\end{array}$$

$$\begin{array}{rl}& \int d\mathbf{r}{p}_{\alpha }(\mathbf{r})=1,\int d\mathbf{r}\mathbf{r}{p}_{\alpha }(\mathbf{r})=0,\\ & \int d\mathbf{r}({\mathbf{e}}_{\alpha }\cdot \mathbf{r}{)}^2{p}_{\alpha }(\mathbf{r})={\sigma }_{\Vert \alpha }^2,\\ & \int d\mathbf{r}[{\mathbf{r}}^2-({\mathbf{e}}_{\alpha }\cdot \mathbf{r}{)}^2]p_{\alpha }(\mathbf{r})=(D-1){\sigma }_{\perp \alpha }^2,\end{array}$$

$${\sigma }_{\Vert ,\perp \alpha }^2={\tilde{\sigma }}_{\Vert ,\perp \alpha }^2\left[1+\beta (\kappa )N{\tilde{\sigma }}_1^2/V\right],$$

$$f_{\mathrm{ex}}\simeq \pi n\int \left[1-e^{-\phi (r)/k_{\mathrm{B}}T}\right]rdr.$$

$$u_{\mathrm{ex}}=u_{\mathrm{st}}+u_{\mathrm{th}},$$

$$M=-1.1061+0.5038\kappa -0.11053{\kappa }^2+0.00968{\kappa }^3+1/\kappa .$$

$$u_{\mathrm{th}}=A\ln (1+B\Gamma /{\Gamma }_{\mathrm{m}}).$$

$$f_{\mathrm{ex}}=M(\kappa )\Gamma -A{\mathrm{Li}}_2(-B\Gamma /{\Gamma }_{\mathrm{m}}),$$

$$u_{\mathrm{ex}}=A(\kappa )\ln \left[1+B(\kappa ){\Gamma }^{s(\kappa )}\right],$$

$$p_{\mathrm{ex}}={\gamma }_{\mathrm{ex}}{u}_{\mathrm{ex}}.$$

$${\gamma }_{\mathrm{ex}}(\kappa )=1+0.526\kappa +0.13{\kappa }^2-0.02{\kappa }^3.$$

Full text

PACS numbers

61.50.-f, 61.66.-f, 07.05.Tp

I Introduction

Systems of particles interacting via Yukawa (screening Coulomb, or Debye-Hückel) potential are widely found in nature. In the physics of soft matter, the Yukawa potential plays particularly important role, because it is traditionally used to describe interactions between ions in screening media (for instance, in aqueous solutions of electrolytes), charged colloidal micro- and nano-particles in various solvents and at interfaces of fluid media, as well as between charged particles in complex (dusty) plasmas. Chaudhuri et al. (2011); Ivlev et al. (2012); Löwen (1994); Morfill and Ivlev (2009) In addition, Yukawa systems represent a useful model of classical interacting particles with the softness of interaction variable in a very wide range, from extremely soft Coulomb interaction (one-component-plasma limit Baus and Hansen (1980) in the absence of screening) to the extremely hard interaction (hard spheres Mulero (2008) in the limit of very strong screening).

Various aspects of structural, dynamical, and thermodynamical properties of Yukawa systems were studied numerically, some representative examples can be found in Refs. Robbins, Kremer, and Grest, 1988; Meijer and Frenkel, 1991; Stevens and Robbins, 1993; Hamaguchi, Farouki, and Dubin, 1996, 1997; Hartmann et al., 2005; Vaulina and Koss, Adamovich; Qi et al., 2010; Heinen et al., 2011; Tolias, Ratynskaia, and de Angelis, 2014. Results of these studies found wide applications, for instance, to explain phase transitions in complex (dusty) plasmas, Nosenko et al. (2009); Hartmann et al. (2010); Klumov (2010); Khrapak et al. (2011, 2012) and in colloidal systems. Löwen (1994); Gasser (2009); Murray and Van Winkle (1987); Royall et al. (2006); Bergenholtz, Poon, and Fuchs (2003); Taffs et al. (2013) In particular, in colloidal suspensions, screened Coulomb repulsion determines various crystalline structures and their properties (see, e.g., Refs. Löwen, 1994; Gasser, 2009; Yethiraj, 2007; Kratzer and Arnold, 2015; Wang et al., 2015; Terao and Nakayama, 1999; Palberg et al., 1995). The screened Coulomb repulsion is the basic interaction for ions and microparticles in electrolytes. Levin (2002); Belloni (2000); Taffs et al. (2013); Kung, Gonzalez-Mozuelos, and de la Cruz (2010); van der Linden et al. (2013); Mohanty, Yethiraj, and Schurtenberger (2012)

Both three-dimensional (3D) and two-dimensional (2D) Yukawa systems can be of interest in the context of colloids, complex plasmas, and electrolytes. In the 2D situation the particles are normally confined to a thin layer or are located at an interface. For instance, 2D Yukawa systems of ions can arise in electrolytes at the interfaces due to the ion specific effects. In drops of aqua solutions of electrolytes, the ionic redistribution near the surface and surface trapping of anions change the surface tension (see, e.g., Refs. Wilson and Pohorille, 1991; Jungwirth, , and Tobias, 2001; Levin, dos Santos, and Diehl, 2009; Levin, 2009). In bulk aqua solutions of electrolytes, similar ion-specific effects lead to the formation of bubbles stabilized by ions, the so-called bubstons, Bunkin et al. (2012a, b, 2016); Yurchenko et al. (2016) in which the ions are adsorbed in thin layer inside the bubble surface to compensate the pressure by surface tension. Two-dimensional plasma crystals and fluids represent one of the major topics of experimental research into complex plasmas in laboratory conditions. Chu and I (1994); Thomas et al. (1994); Hayashi and Tachibana (1994); Melzer, Trottenberg, and Piel (1994); Thomas and Morfill (1996); Melzer, Homann, and Piel (1996); Schweigert et al. (1998); Morfill et al. (2004); Ivlev, Nosenko, and Röcker (2015) Thus, 2D systems of Yukawa particles occur in a rather broad range of applications and related problems are of fundamental importance.

Thermodynamic properties of 2D Yukawa systems have been of considerable continuous interest in the last couple of decades, largely in the context of complex plasmas. Totsuji et al. (2004); Hartmann et al. (2005); Vaulina and Koss (Adamovich, 2014); Khrapak et al. (2015a); Feng et al. (2016a, b) However, to the best of our knowledge, no comprehensive results across coupling regimes along with simple and reliable approximations convenient for practical use have been proposed. In the present paper, we study systematically thermodynamics of 2D Yukawa systems in a very broad parameter regime, from very weakly interacting gaseous state to strongly interacting fluid and solid states. Using MD simulations we systematically calculate the excess energy and pressure. For gases and fluids simple practical expressions for the thermodynamic properties are then proposed. For crystals, the shortest-graph interpolation method Yurchenko (2014); Khrapak et al. (2015b); Yurchenko, Kryuchkov, and Ivlev (2015, 2016) is applied to calculate the pair correlations and thermodynamic properties. The advantages of this method as well as its excellent accuracy are demonstrated. Overall, this paper describes simple and reliable tools to calculate the thermodynamic properties of 2D Yukawa systems across coupling regimes with required accuracy.

II Methods

II.1 System description

We investigate a classical system of point-like particles in the 2D geometry interacting via the pairwise repulsive Yukawa potential of the form

[TABLE]

where $\varepsilon$, and $\lambda$ are the energy and (screening) length scales of the interaction. For charged particles immersed in a plasma-like screening environment, the energy scale is $\varepsilon=Q^{2}/4\pi\epsilon_{0}\lambda$ (in SI units), where $Q$ is the charge and $\epsilon_{0}$ is the permittivity of free space. The properties of Yukawa systems are determined by the two dimensionless parameters. The first is the coupling parameter, $\Gamma=(Q^{2}/4\pi\epsilon_{0}ak_{\mathrm{B}}T)$, where $k_{\mathrm{B}}$ is the Boltzmann constant, $T$ is the temperature, $a=(\pi n)^{-1/2}$ is the 2D Wigner-Seitz radius, and $n=N/V$ is the areal density of $N$ particles occupying the 2D volume $V$. The second is the screening parameter, $\kappa=a/\lambda$. Note, that the coupling parameter is roughly the ratio of the potential energy of interaction between two neighbouring particles to their kinetic energy. The system is usually said to be in the strongly coupled state when this ratio is large, that is $\Gamma\gtrsim 1$.

When coupling increases the system forms a strongly coupled fluid phase, which can crystallize upon further increase in $\Gamma$. This fluid-solid transition can be characterized by the temperature and/or coupling parameter, $T_{\rm m}$ and $\Gamma_{\rm m}$, where the subscript “m” refers to melting. Both $T_{\rm m}$ and $\Gamma_{\rm m}$ are the functions of the screening parameter $\kappa$. The dependence $\Gamma_{\rm m}(\kappa)$ has been approximated in Ref. Hartmann et al., 2005 by the following fit:

[TABLE]

This fit describes relatively well the melting points found from the bond angular correlation analysis (see Fig. 6 of Ref. Hartmann et al., 2005) up to $\kappa=3.0$ and it should not be applied for larger $\kappa$. In the limit $\kappa=0$ the system reduces to the 2D one-component-plasma (OCP) with the Coulomb interaction. In this case $\Gamma_{\rm m}\simeq 131$ lies in the range predicted in earlier numerical simulations Gann, Chakravarty, and Chester (1979) and obtained in experiments with a classical 2D sheet of electrons Grimes and Adams (1979) (see also Ref. Khrapak and Khrapak, 2016 for a recent overview of OCP thermodynamics in 2D and 3D).

Finally, it is worth to comment on the nature of the fluid-solid phase transition in 2D Yukawa systems. Recently, it has been demonstrated that the potential softness is very important factor, which determines the melting scenario. Kapfer and Krauth (2015) For sufficiently steep repulsive interactions the hard-disk melting scenario holds: a first-order liquid-hexatic and a continuous hexatic-solid transition can be identified.  Bernard and Krauth (2011); Engel et al. (2013) For softer interactions the liquid-hexatic transition is continuous, with correlations consistent with the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) scenario. (For example, in 2D colloidal systems, hexatic phase was observed in the experiment by Zahn et al. Zahn, Lenke, and Maret (1999)) For the Yukawa potential the transition between these two scenarios occurs at about $\kappa\simeq 6$. Kapfer and Krauth (2015) Below we consider systems with $\kappa$ in the range from $0.5$ to $3.0$ (this range is particularly relevant to 2D plasma crystals and fluids in laboratory experiments Fortov et al. (2004, 2005),Rios de Anda et al. (2015)), thus belonging to the soft interaction class. In this range of $\kappa$, the hexatic phase occupies a rather narrow region on the phase diagram, Kapfer and Krauth (2015) and the study of its properties is beyond the scope of the present investigation.

II.2 Computational details

To obtain the thermodynamic properties of the 2D Yukawa systems across coupling regime, extensive MD simulations have been performed. The MD simulations have been done in the $NVT$ ensemble at different temperatures using $N=64000$ particles and the Langevin thermostat. The numerical time step was chosen $\Delta t_{c}=5\times 10^{-4}\sqrt{m\lambda^{2}/\epsilon}$ for the crystalline phase and $\Delta t_{c}\sqrt{\Gamma/\Gamma_{\rm m}}$ for the fluid phase. The cutoff radius of the Yukawa potential was set equal to $15n^{-1/2}$. The simulations were run for $1.5\times 10^{6}$ time steps to equilibrate the system and obtain the equilibrium properties. In the simulation run with $\kappa=0.5$ Ewald summation was implemented.

The simulations have been performed for a number of screening parameters $\kappa$ ranging from $0.5$ to $3.0$. This corresponds to sufficiently soft interactions as discussed above. For each value of the screening parameter $\kappa$, twelve simulation runs correspond to the fluid phase and nine runs to the crystalline phase. In the fluid phase the coupling parameter ranges from $\Gamma=0.5$ to $\simeq 0.95\Gamma_{\rm m}$. In the solid phase the values corresponding to $\Gamma_{\rm m}/\Gamma=0.9,0.8,...,0.1$ are taken.

The main simulation results are summarized in Tables 2-5 of the Appendix.

II.3 Thermodynamic definitions and relations

The main thermodynamic quantities which will be required below are the internal energy $U$, Helmholtz free energy $F$, and pressure $P$ of the system. The following thermodynamic definitions exist Landau and Lifshitz (2013)

[TABLE]

In addition, $U$ and $P$ can be calculated using the integral equations of state Hansen and MacDonald (2006); Frenkel and Smit (2001)

[TABLE]

where $g(\mathbf{r})$ denotes the radial distribution function, which is isotropic in gas and fluid phases and anisotropic in the crystalline phase.

We will use conventional reduced units: $u=U/Nk_{\rm B}T$, $f=F/Nk_{\rm B}T$, and $p=PV/Nk_{\rm B}T$ and divide the thermodynamic quantities into the kinetic (ideal gas) and potential (excess) components, so that $u=1+u_{\rm ex}$ (in 2D), $f=f_{\rm id}+f_{\rm ex}$, and $p=1+p_{\rm ex}$. Finally, it is useful to operate with the Yukawa system state variables $\Gamma$ and $\kappa$. In these variables the thermodynamic identities for 2D Yukawa fluids are Khrapak et al. (2015a); Semenov, Khrapak, and Thomas (2015)

[TABLE]

II.4 The shortest-graph method

To describe the thermodynamics of 2D Yukawa crystals analytically, we employ the shortest-graph method, proposed and developed in Refs. Yurchenko, 2014; Yurchenko, Kryuchkov, and Ivlev, 2015, 2016. Following these papers, thermodynamical properties of classical crystals can be obtained very accurately from the following consideration. The anisotropic pair-correlation function $g(\mathbf{r})$ of a crystal is written in the form

[TABLE]

where the summation is over all the nodes $\alpha$, and each individual peak has the shape

[TABLE]

The normalization constant as well as the parameters $a_{\|,\perp\alpha}^{2},b_{\alpha}$ are defined by the following conditionsYurchenko, Kryuchkov, and Ivlev (2016)

[TABLE]

where $D=2$ is the spatial dimensionality and $\mathbf{e_{\alpha}}=\mathbf{r_{\alpha}}/r_{\alpha}$ is the unit vector in the direction of $\mathbf{r_{\alpha}}$, $\sigma_{\|,\perp}^{2}$ is the mean squared displacement for longitudinal and transversal directions, respectively, calculated using the finite-temperature phonon spectra, taking into account the anharmonic effects.Yurchenko, Kryuchkov, and Ivlev (2016) By using the pair correlation function $g(\mathbf{r})$ the excess energy and pressure can then be obtained. However, calculation of the finite-temperature phonon spectra is a difficult problem, which is beyond the scope of the present paper. Therefore, we propose here a simpler practical approach, which yields very accurate results and can be used for practical calculations.

Due to the anharmonicity of phonon spectra at finite temperatures, the second-order term becomes more significant in the temperature expansion of the mean-squared displacements $\sigma^{2}$. To account for this effect, we propose the anharmonic correction of the mean-squared displacements

[TABLE]

where the tildes denote the mean-squared displacement calculated using zero-temperature phonon spectra (see Ref.Yurchenko, Kryuchkov, and Ivlev, 2015), $\widetilde{\sigma}_{1}^{2}$ is the total mean-squared displacement for the nearest neighbours, and we have introduced the anharmonic correction coefficient $\beta(\kappa)$, which does not depend on the temperature and should be found using MD simulations for different screening parameters. The correction given by Eq.(9) conserves the ratio $\sigma_{\|}^{2}/\sigma_{\perp}^{2}$ between the mean-squared displacements in the longitudinal and transversal directions. A posteriori comparison with MD results proves that this assumption allows to obtain excellent accuracy.

III Results

III.1 Weakly-coupled fluids

A simple and physically transparent approach to the thermodynamics of weakly coupled Yukawa systems for small deviations from the ideal gas behavior is to calculate the second virial coefficient. This has recently been shown to work well in 3D Yukawa systems. Khrapak (2016) In the 2D geometry the excess free energy is expressed in this approximation as

[TABLE]

The excess energy and pressure can be readily obtained from the excess free energy. We compare the values $u_{\rm ex}$ at a fixed coupling parameter $\Gamma=0.5$ obtained from Eq. (10) and computed using MD simulations in Fig. 1. The agreement is satisfactory: in the range of $\kappa$ investigated the deviations are within several percent. The agreement naturally improves with increasing $\kappa$, because at a fixed $\Gamma$ the actual interaction strength weakens as $\kappa$ increases.

III.2 Strongly-coupled fluids

The excess energy and pressure of the 2D Yukawa fluids have been determined using MD simulations in a wide range of coupling and screening parameters. The results are summarized in the Table 2 of the Appendix. Here we describe simple analytical approximations, which can be used to evaluate the energy and pressure for practical purposes.

In the strongly coupled fluid regime it is helpful to divide the thermodynamic quantities, such as energy and pressure, into static and thermal contributions. The static contribution corresponds to the value of internal energy when the particles are frozen in some regular configuration and the thermal corrections arise due to the deviations of the particles from these fixed positions (due to thermal motion). Of course, such a division is only meaningful when the regular structure is specified. For crystals, the obvious choice is a corresponding lattice sum (Madelung energy). For fluids this choice is also meaningful and we use it here (Note, that in 3D Yukawa system a slightly different definition of the static fluid energy is traditionally employed. Khrapak (2016); Khrapak et al. (2014))

The excess internal energy is thus a sum of the static and thermal contributions,

[TABLE]

where $u_{\rm st}=M\Gamma$ and $M$ is the Madelung constant. The values of the Madelung constant for 2D Yukawa systems in the regime of relatively weak screening, $0.5\leq\kappa\leq 3.0$, are tabulated in Table 1. The dependence $M(\kappa)$ can be fitted using a functional form similar to that proposed by Totsuji *et al.*Totsuji et al. (2004)

[TABLE]

The last term in (12) accounts for the absence of neutralizing background in our case (but present in Ref. Totsuji et al., 2004), the energy of this background being simply $-\Gamma/\kappa$. The fit is chosen in such a way that when $\kappa\rightarrow 0$ and the neutralizing background is introduced, the Madelung constant is reduced to the well known value of the triangular lattice sum of the 2D one-component-plasma (OCP) with Coulomb interactions, $M_{\rm OCP}\simeq-1.1061$. This fit is accurate to within a tiny fraction of percent for $\kappa\lesssim 1.0$ and to within $\sim 1\%$ when screening becomes stronger ($\kappa\sim 3$).

The thermal part of the excess energy is expected to exhibit a quasi-universal scaling with respect to the reduced coupling parameter $\Gamma/\Gamma_{\rm m}$. This is a general property of classical particle systems with sufficiently soft interactions, which was first pointed out by Rosenfeld and Tarazona (RT scaling) for 3D systems. Rosenfeld and Tarazona (1998); Rosenfeld (2000) In the context of 3D Yukawa systems, the RT scaling has been proven to be very useful in Refs. Khrapak et al., 2015b; Khrapak, 2016; Khrapak and Thomas, 2015a, b The emergence of RT scaling analogue for 2D systems has been discussed in the context of OCP with Coulomb and logarithmic interactions, Yukawa systems near the OCP limit, and inverse-power-law interactions. Khrapak and Khrapak (2016); Khrapak et al. (2015a) The dependence of $u_{\rm th}$ on $\Gamma/\Gamma_{\rm m}$ in the strongly coupled regime is displayed in Fig. 2. The quasi-universality is well pronounced, although there is clearly some systematic tendency of decreasing the value of $u_{\rm th}$ with $\kappa$ at the same value of $\Gamma/\Gamma_{\rm m}$. This tendency is expected when the potential steepness increases (see e.g. Fig. 4 from Ref. Khrapak et al., 2015a). Overall, the data points corresponding to the dependence $u_{\rm th}(\Gamma/\Gamma_{\rm m})$ are confined to a relatively narrow range. The important point is that towards the side of soft interactions (sufficiently small $\kappa$ in our case), the static component of the internal energy is dominant over the thermal one. For example, at $\kappa=1$ the thermal component contributes only to about $2\%$ of the total excess energy near the fluid-solid phase transition. Therefore, even moderately accurate fits for $u_{\rm th}$ allow to obtain high accuracy with respect to the total excess energy $u_{\rm ex}$.

Three fits are shown in Fig. 2. The upper (lower) curve corresponds to the data portion for $\kappa=0.5$ ($\kappa=3.0$). The intermediate curve has been obtained using the entire massive of the data points (corresponding to the parameter regime shown). It can be considered as representative for strongly coupled 2D Yukawa fluids in the vicinity of the freezing transition. The functional form of the fit is the same as used previously Khrapak et al. (2015a)

[TABLE]

The use of the coefficients $A=0.257$ and $B=195.4$ determined here would somewhat improve previous approximations.

The excess free energy can be routinely calculated using the model for the excess energy formulated above and the second of Eqs. (5). The resulting expression is rather simple,

[TABLE]

where ${\rm Li}_{2}(z)=\int_{z}^{0}dt\ln(1-t)/t$ is dilogarithm. Note that in deriving Eq. (14), the thermodynamic integration over the coupling parameter from 0 to $\Gamma$ has been performed, while Eq. (13) is strictly speaking not applicable at $\Gamma\ll 1$. The correct procedure would be to start thermodynamic integration from some small, but finite value $\Gamma_{0}$, and then add the constant $f_{\rm ex}(\Gamma_{0})$ evaluated using Eq. (10). However, since the actual contribution from the weakly coupling regime is small, Eq. (14) remains rather accurate at strong coupling and we use it here.

The calculation of pressure from the excess free energy is straightforward, but rather cumbersome in the considered case. This is because the differentiation with respect to $\kappa$ is involved, and the two fits for $M(\kappa)$ and $\Gamma_{\rm m}(\kappa)$ are present. For this reason, the explicit expression for $p$ is not displayed. We verified that near freezing (at $\Gamma/\Gamma_{\rm m}\simeq 0.95$) the derived expression yields the pressures which deviate from the exact MD results by $\sim 0.001\%$ at $\kappa=0.5$, $~{}\sim 0.1\%$ at $\kappa=1.0$, and $\sim 1\%$ at $\kappa=2.0-2.8$. The accuracy drops at the highest value $\kappa=3.0$. This is not surprising, since the fits for $M(\kappa)$ and $\Gamma_{\rm m}(\kappa)$ are only applicable for $\kappa\lesssim 3.0$ and, therefore, derivatives from these fits at $\kappa=3.0$ can produce significant errors.

We also found out that if better accuracy is required, the data for the excess thermal energy can be fitted by the following slightly modified expression

[TABLE]

where $A$ and $B$ are now assumed $\kappa$-dependent and a $\kappa$-dependent exponent $s$ is introduced. Based on all the data points obtained in MD simulations the following relations are identified: $A(\kappa)=0.35708+0.09397\kappa$, $B(\kappa)=1.65491\exp(-0.76911\kappa)$, $s(\kappa)=0.68838-0.05183\kappa$. Some representative examples are shown in Fig. 3. The fit of Eq. (15) is clearly more accurate and can be used in the regime of weaker coupling, compared to the simple form (13). However, it is also less practical in evaluating thermodynamic parameters other than the excess internal energy.

III.3 Relation between excess pressure and energy

It is sometimes advantageous to operate with an equation of state written in the form of relation between the pressure and internal energy of the system. For soft purely repulsive potentials a simplest formulation of this kind can be written as

[TABLE]

Here the parameter $\gamma_{\rm ex}$ generally depends both on the temperature and density, that is both on $\Gamma$ and $\kappa$ for Yukawa systems. Note that the parameter $\gamma_{\rm ex}$ introduced in this way is not directly related to the conventional definitions of either the density scaling exponent or Grüneisen parameter. Hummel et al. (2015) Nevertheless, it may be helpful in characterizing the softness of the repulsive potential. We remind that for inverse-power-law (IPL) repulsive potentials of the form $\varphi(r)\propto r^{-\alpha}$ the relation between the excess pressure and energy is particularly simple, $p_{\rm ex}=\tfrac{\alpha}{2}u_{\rm ex}$ in 2D. Thus, an “effective IPL exponent” may be associated with the quantity $2\gamma_{\rm ex}$.

Having approximations for both $p_{\rm ex}$ and $u_{\rm ex}$ for 2D Yukawa fluids we can easily estimate the value of $\gamma_{\rm ex}$. The corresponding plot of $\gamma_{\rm ex}$ as a function of Yukawa systems state variables $\kappa$ and $\Gamma/\Gamma_{\rm m}$ is shown in Fig. 4. To produce this plot, Eq. (13) for the thermal component of the excess energy has been used. Figure 4 shows that in the strongly coupled regime $\gamma_{\rm ex}$ is very weakly dependent on the coupling strength (temperature), but exhibits considerable dependence on $\kappa$ (density). Using the exact MD results for $p_{\rm ex}/u_{\rm ex}$ in the vicinity of the fluid-solid phase transition ($\Gamma/\Gamma_{\rm m}\simeq 0.95$) we have obtained a representative dependence $\gamma_{\rm ex}(\kappa)$ in the strongly coupled regime:

[TABLE]

Importantly, $\gamma_{\rm ex}\rightarrow 1$ as $\kappa\rightarrow 0$. This seems counter-intuitive at first, because one would naturally expect $\gamma_{\rm ex}=\tfrac{1}{2}$ in the OCP Coulomb interaction limit in 2D. The difference is attributed to the presence of the neutralizing background in the OCP model. In the limit of very soft interaction, the energy and pressure are dominated by their static contributions. As $\kappa\rightarrow 0$, the dominant contribution is the Madelung energy, so that $f_{\rm ex}\sim u_{\rm ex}\sim M\Gamma\sim\Gamma/\kappa$ (without background). This implies $p_{\rm ex}=\tfrac{\Gamma}{2}(\partial f_{\rm ex}/\partial\Gamma)-\tfrac{\kappa}{2}(\partial f_{\rm ex}/\partial\kappa)\sim\Gamma/\kappa\sim u_{\rm ex}$. In the presence of neutralizing background the term $\Gamma/\kappa$ disappears and we have $f_{\rm ex}\sim u_{\rm ex}\sim M_{\rm OCP}\Gamma$. This yields $p_{\rm ex}\sim\tfrac{1}{2}M_{\rm OCP}\Gamma\sim\tfrac{1}{2}u_{\rm ex}$. This consideration demonstrates that the Yukawa systems in the limit $\kappa\rightarrow 0$ are not fully equivalent to the Coulomb systems with the neutralizing background.

III.4 Crystals

In a series of MD simulations for 2D Yukawa crystals, in addition to evaluate the excess energy and pressure (which are summarized in Tables 4 and 5 of the Appendix), the mean squared displacements were calculated to find the anharmonic correction coefficient $\beta$. The resulting dependence $\beta(\kappa)$ is shown in Figure 5 (the corresponding values are also tabulated in Table 3 of the Appendix for completeness). The inset in Fig. 5 presents the radial (isotropic) pair correlation function, $g(r)\propto\int{d\varphi\;g(\mathbf{r})}$, and demonstrates excellent representation of the short- and long-distance correlations. The obtained anharmonic correction coefficient $\beta(\kappa)$ allows to calculate analytically pair correlation function and then the excess energy, pressure and other thermodynamic parameters by the thermodynamic integration with the help of the expressions given in Sec. II.3.

It is worth to point out the following observation: In the limit $\kappa\rightarrow 0$, the Yukawa interaction tends to the unscreened Coulomb interaction $\varphi\propto r^{-1}$. According to our previous MD simulations, Yurchenko, Kryuchkov, and Ivlev (2015) the finite-temperature phononic spectra differ weakly from zero-temperature ones for IPL potentials, $\varphi\propto r^{-\alpha}$. Therefore, in the OCP limit ($\kappa=0$ and $\alpha=1$) we should obtain the smallest values of $\beta(\kappa)$. This is indeed observed in Fig. 5.

In figure 6 we plot the reduced pressure versus the reduced excess energy of 2D Yukawa fluids and solids. Symbols are the MD results, the solid and dashed curves correspond to the shortest-graph method [with found anharmonic correction coefficient $\beta(\kappa)$] for the crystalline phase and the proposed fit by Eq.(15) for the fluid phase, respectively. Excellent agreement is observed.

III.5 Accuracy

The relative difference between the excess energies calculated using the shortest-graph method and those evaluated using direct MD simulations in the solid phase amounts to $\simeq 5\times 10^{-5}$, which is comparable to the values reported earlier. Yurchenko, Kryuchkov, and Ivlev (2016) The accurate fit of Eq. (15) yields the relative error in the excess energy smaller than $5\times 10^{-4}$ and $2\times 10^{-3}$ for 72% and 95% of the examined fluids data points, respectively. Maximal relative deviation, $5\times 10^{-3}$, is observed near the melting line at large values of the screening parameter $\kappa$. A simpler fit of Eq. (13) is applicable when the relative deviations within $\lesssim 1\%$ are acceptable.

In addition, we can compare our results with those recently reported in Refs. Feng et al., 2016a, b, where fits for the pressure of 2D Yukawa fluids in the $(\kappa,\Gamma)$ parameter space have been proposed. The case $\kappa=0.5$ received special attention and a simple two-term fit has been proposed based on the results of a MD simulation, Feng et al. (2016a) $p=1.53\Gamma+1.33$. We plot our MD results along with the fit of Eq. (13) and the fit from Ref. Feng et al., 2016a in Fig. 7. One can see that the fit from Ref. Feng et al., 2016a overestimates the pressure systematically at high values of $\Gamma$. At the strongest coupling in the fluid phase studied in this work, $\Gamma=135.42$, the present MD simulation yields $p=199.434$, fit by Eq. (13) yields $p=199.432$, while the fit from Ref. Feng et al., 2016a yields $p=208.523$. On the other hand, the previous model for 2D Yukawa systems in the OCP (weakly screening) limit discussed in Refs. Khrapak et al., 2015a; Semenov, Khrapak, and Thomas, 2015 yields $p=199.445$, providing confidence in the accuracy of the present results. The reasons for deviations in Ref. Feng et al., 2016a have to be identified.

IV Conclusion

We studied thermodynamics of 2D classical Yukawa systems across coupling regimes, from the weakly non-ideal gas to the strongly coupled fluid and crystalline phases. Careful analysis of the extensive MD simulation results allowed us to propose simple and physically suitable expressions for the internal energy and pressure of the studied systems. These expressions can be used to estimate other thermodynamic properties using the conventional thermodynamic relations.

For weakly non-ideal gases, virial expansion with the two first terms retained was shown to provide reasonably good estimate of the excess energy at sufficiently weak coupling, $\Gamma\lesssim 1$.

For the strongly-coupled fluid phase, we made use of the 2D analogue of the 3D Rosenfeld-Tarazona quasi-universal scaling of the thermal component of the excess energy. This quasi-universal scaling was shown to be particularly useful on approaching the fluid-solid phase transition. Deviations from the quasi-universal behaviour have been discussed and quantified.

To calculate thermodynamic properties of the crystalline phase, we employed the shortest-graph method for pair correlation functions. To account for the effects of finite-temperature phonon spectra, we proposed a simple way to correct the mean squared displacements of correlation peaks for different nodes. The coefficient of anharmonic correction was evaluated in MD simulations and then used in analytical estimates. The efficiency and accuracy of the approach was documented.

The results of this paper can be useful for thermodynamic calculations related to various phenomena in 2D and quasi-2D Yukawa fluids and solids in a broad range of parameters. In particular, this includes colloidal systems, complex (dusty) plasmas, and aqua solutions of electrolytes for ionic fluids.

Acknowledgements.

The numerical simulations are supported by Russian Science Foundation, Project No. 14-29-00277. Post-processing is supported by Russian Foundation for Basic Researches (Projects Nos. 15-38-21009 and 16-38-00952). The theoretical study by the shortest-graph method is supported by Russian Science Foundation, Project No. 14-43-00053. The present position of SAK at Aix Marseille University is supported by the A*MIDEX project (Nr. ANR-11-IDEX-0001-02) funded by the French Government “Investissements d’Avenir” program managed by the French National Research Agency (ANR).

Appendix A MD results

In the Appendix, we summarize main results from MD simulations performed in this study. Table 2 reports the reduced excess energies and pressures at different state points in the fluid phase. Table 3 summarizes the values of the anharmonic correction coefficient $\beta$ evaluated using MD simulations of the crystalline phase. Finally, Tables 4 and 5 report the excess energies and pressures in the crystalline phase.