Phase diagram and universality of the Lennard-Jones gas-liquid system

University and scaling are key concepts in the study of critical systems,^1,2 including liquid-gas systems,^3–5 Ising model,^6,7 percolation model,^8–12 dimer model,¹³ etc. According to the modern theory of critical phenomena,^1,2 critical systems can be classified into different universality classes such that the systems in the same class have the same set of critical exponents.

The Lennard-Jones (LJ) particle system has been extensively studied since it is the basic model of the off-lattice systems with short-range interactions. It is widely believed that the gas-liquid phase transition of the particle system with short-range interaction belongs to the Ising universality class¹⁴ and the LJ-particle system is no exception. However, there is no definite proof that the universality class of the liquid-gas phase transition of the LJ system is same as that of the Ising universality. In order to study the universality class of the LJ system, the precise information of the phase diagram, especially the location of the critical point, is required. Aside from the critical behavior, the determination of the phase boundary itself is important to study nonequilibrium phenomena such as bubble nuclei.¹⁵ 

Recently, due to the advance in the computational power, many researchers have tried to determine the value of the critical exponents as accurately as possible, in order to confirm whether the universality class of the LJ system is really the same as that of the Ising model. Panagiotopoulos proposed the Gibbs ensemble Monte Carlo (GEMC) method to calculate the phase equilibria of the coexisting phase.¹⁶ The GEMC method involves two simulation boxes, one is in liquid phase and the other is in gas phase under a given temperature. There are three types of Monte Carlo (MC) steps; displacements of particles, changing volumes of simulation boxes, and exchanging particles between the simulation boxes. When the two systems reach their equilibrium states, we can obtain the physical quantities such as densities of gas and liquid phases, vapor pressure, and chemical potential of the coexisting state in the given temperature. While this method is simple, the system can be unstable near the criticality, e.g., the simulation box may change its phase from gas to liquid and vice versa because of the fluctuation.

Möller and Fischer proposed a different approach, NpT and test particle method.¹⁷ This method involves NpT simulation, i.e., the isobaric and isothermal condition. The chemical potential is observed as a function of temperature by Widom's test-particle method,¹⁸ and the critical point is determined by the crossing point of the chemical potential of liquid and gas branches. This method is more stable than GEMC near the critical point; since only the volume fluctuates while both density and volume fluctuate in GEMC. Additionally, Okumura and Yonezawa¹⁹ proposed NVT and test particle method which is more stable than the last two methods since this method does not involve change in volume. However, the calculation of the chemical potential with Widom's test-particle method is quite expensive, especially for larger systems. Wilding and Bruce combined the histogram reweighting (HR) techniques with the finite-size scaling (FSS) analysis.²⁰ The main idea of this method is to locate the apparent critical point from the matching condition between the ordering operator distribution and that of the Ising model with the help of the mixed-field theory. Pérez-Pellitero et al. compared the HR techniques with the direct estimation of the critical point using Binder parameter.²¹ 

While the HR techniques combined with FSS analysis are efficient to study the coexisting phase, the distribution of the order parameter of the Ising model at the critical point is required for the reference. Additionally, most of past studies determined the critical point with the help of the knowledge of the Ising universality, and therefore, it is difficult to determine whether the gas-liquid phase transition of LJ-particle system belongs to the Ising universality class.

In the present study, we propose a simple method to study the gas-liquid coexisting phase using the molecular dynamics (MD) simulations. The main purpose of our method is to determine the critical properties of the LJ gas-liquid system without a priori knowledge of the Ising universality, such as the values of the critical exponents, the value of Binder parameter at the critical point, or the distribution of the order parameter. The critical exponents of the order parameter is estimated from the liquid-gas coexisting density and that of the correlation length is estimated from the interface thickness taking into account the effect of the capillary waves. We have found that the critical exponents of the order parameter and correlation length are β = 0.3285(7) and ν = 0.63(4), respectively, which are consistent with those of the Ising universality class.^1,22–24

This paper is organized as follows. In Sec. II, we describe our model system and details of the simulation procedure. We first describe how to determine the coexisting curve in Sec. II A, then we describe how to locate the critical point using the Binder parameter with the aid of the law of rectilinear diameter in Sec. II B. In Sec. III, we present results of our simulation data and critical exponents determined from such data. Section IV contains our conclusions and discusses problems for further studies.

We use the truncated Lennard-Jones potential of the form

with the well depth ɛ and atomic diameter σ.²⁵ The coefficients c₂ and c₀ are determined so that V(r_c) = V′(r_c) = 0 with the cut-off length r_c , i.e., the values of potential and force become continuously zero at the truncation point and V(r) = 0 for r > r_c. The two additional terms are introduced for good conservation of the total energy. The extra terms, which make the force continuous at the truncation point, are required to perform the MD simulations stably for long time and are unnecessary for MC simulations. In the following, we use the physical quantities reduced by σ, ɛ, and the Boltzmann constant k_B, e.g., the length is measured with the unit of σ, and so forth. We set the cut-off length as r_c = 3.0. The simulation box is a rectangular parallelepiped with sizes L_x × L_y × L_z. The periodic boundary condition is taken for all directions. We set the ratio of L_x: L_y: L_z = 1: 1: 2 in order to achieve the gas-liquid coexisting phase. The studied system sizes are L_x = 32, 64, 96, and 128. We setup the initial configuration by placing particles densely in the region z < L_z/2 and sparsely in the region z ⩾ L_z/2, and we equibrate the system with the heatbath. Nosé–Hoover heatbath is used for controlling temperature.²⁶ Time integration scheme is the second order Reversible System Propagator Algorithm (RESPA)²⁷ for isothermal simulations and the second order symplectic integrator for isoenergetic simulations. The time step is chosen to be 0.005 throughout the simulations. We check whether the system has reached the equilibrium state by observing the potential energy, since the relaxation of the potential energy is much slower than that of the temperature (see Fig. 1). The potential energy U(t) behaves as U(t) − U_eq ∝ exp (− t/τ) with the equilibrium value U_eq and the relaxation time τ. We perform thermalization at least ten times longer than the relaxation time. Then we turn off the heatbath and start observation of the density profile at and near gas-liquid phase interface and surface tension. The simulations are performed with MPI parallelized MD program.^28,29 The largest run, involving 1 583 504 particles and 5 000 000 steps, takes about 8 h using 1024 cores at SGI Altix ICE 8400 EX.

The functional form of the density profile at and near the gas-liquid interface is of the form

where ρ_l, ρ_g, and λ are the density of liquid, the density of gas, and the thickness of the interface, respectively. The typical density profile is shown in Fig. 2. By fitting Eq. (2) to the simulation results, we obtain the gas density ρ_g, the liquid density ρ_l, and the thickness λ as functions of temperature. The thickness λ is associated with the correlation length, and becomes longer as the temperature approaches to the critical point. The finite-size effect can be ignored while λ is sufficiently shorter than the system size L_x. This means that larger system size allows us to approach closer to the critical point. Typical configurations of particles below, near, and above the critical point are shown in Figs. 3(a)–3(c), respectively.

When the phase interfaces are normal to the z axis, the surface tension γ is defined by the following equation:

where π_xx, π_yy, and π_zz are the diagonal components of the stress tensor.³⁰ The integral is taken over the whole volume of the simulation box. The factor 2 of the lhs comes from the fact that there are two interfaces in the system due to the periodic condition. For the particles interacting via two body potential V(r), Eq. (3) reduces as

where r_ij is the distance between particles i and j, V′(r_ij) is the derivative of V(r_ij) with respect to r_ij, x_ij is the x component of the relative position vector between particles i and j, and so forth. The sum is taken over all pairs. Note that the kinetic term of the stress tensor is omitted since the kinetic term does not contribute to the surface tension. After the thermalization, the surface tension is sampled every 1000 steps and averaged. The observation is performed for 100 000 steps for each temperature.

The order parameter of the gas-liquid transition is the density difference Δρ ≡ ρ_l − ρ_g. The order parameter becomes zero at the critical temperature T_c, and the behavior near the criticality is

with the reduced temperature ɛ ≡ |T_c − T|/T_c and the critical exponent β. The surface tension γ and the thickness of the interface λ also show critical behavior as follows:

with the critical exponent of the correlation length ν. It is difficult to determine the critical temperature and the critical exponents simultaneously only from the above scaling relations without assuming the Ising universality. Therefore, we determine the critical point using the Binder parameter.³¹ 

In order to compute the fluctuation of density from NVT or NVE simulations, we adopt the block distribution analysis.³² The main idea of this method is to divide the system into small cells, and each cell exhibits approximately the grandcanonical ensemble, while cells are not completely independent from each other. Suppose there are a cubic system with linear dimension L containing N particles. The simulation box with linear dimension L is divided in to cells of size L/M_b, where M_b is an integer. Each cell contains different number of particles, and therefore, has different density ρ_i, i being the label of the cells. The moments of density required for the calculation of Binder parameter are defined to be

Using the above definition, we can define the Binder parameters as follows:

where ⟨⋅ ⋅ ⋅⟩ denotes ensemble average. With different dividing number M_b, we can obtain the Binder parameters for different volume sizes from a single run for a sufficiently large system.

The crossing point of the Binder parameter gives the critical point. Since we perform NVT ensemble, the density, and the temperature should be given simultaneously. Therefore, we assume the law of rectilinear diameters³³ that the average density shows linear dependence on temperature as (ρ_g + ρ_l)/2 = aT + b. While some fluids exhibit strong deviations from the above linear relation, simple pure fluids such as oxygen and xeon follow this law within high accuracy.³⁴ The studied system is a cube with linear dimension L = 128. First the system is thermalized for 200 000 steps with heatbath, and the Binder parameters are observed every 1000 steps. The observation is performed for 100 000 steps. The Binder parameters are calculated for M_b = 16, 24, and 32.

While the critical exponent can be also estimated from the finite-size analysis of Binder parameter, it is better to estimate the critical exponent β from the liquid-gas coexisting density. Because the coexisting density can be determined more accurately than Binder parameter, since the coexisting density is the first order moment of the density while Binder parameter contains the higher moments. Therefore, we use Binder parameter only for locating the critical point and we estimate the critical exponent from the liquid-gas coexisting density.

We first calculate the phase diagram, and determine the coefficients of the linear relation. Then we calculate the Binder parameter on the line in order to determine the critical temperature. Using the obtained critical temperature, we determine the critical exponents.

The gas-liquid phase boundary obtained from the simulations is shown in Fig. 4. The obtained gas and liquid coexisting densities do not exhibit finite-size effect, i.e., the values of different system sizes are within the statistical errors. We also include the experimental data of neon which is denoted by “Exp.”³⁵ The obtained results for the largest systems L_x = 128 are listed in Table. I. We confirm that the law of rectilinear diameters (ρ_g + ρ_l)/2 = aT + b is well satisfied. We determine the coefficients to be a = −0.195(1) and b = 0.533(1). The obtained line is shown as the dashed line in Fig. 4. Note that we extend the line into the supercritical phase while the law of rectilinear diameters does not make sense outside the coexisting phase, since we want to calculate Binder parameter in the region over the critical point in order to have the intersection point.

The Binder parameters are calculated on this line and they are shown in Fig. 5. The Binder parameter does not intersect exactly at one point as reported before,²¹ but we can determine the critical temperature with reasonable accuracy as T_c = 1.100(5). Using the critical point, we determine the value of the critical exponent β from the slope of ln Δρ/ln ɛ, as shown in Fig. 6. The obtained value β = 0.3285(7) is consistent with the reported value β = 0.3269(6) of the spontaneous magnetization¹ of the three-dimensional Ising model on the simple cubic lattice obtained by Talapov and Blöte²² or β = 0.325(5) by Ito et al.²³ of the Ising universality class. The critical behavior of the surface tension is shown in Fig. 7. The critical exponent of the correlation length is determined to be ν = 0.640(4) which is slightly larger than that of the Ising universality class ν = 0.63002(10).²⁴ 

We also study the critical behavior of the interface thickness. We find that the interface thickness exhibits strong finite-size effect, as shown in Fig. 8(a), while the finite-size effect on the surface tension is almost negligible. These finite-size dependencies come from the interface broadening due to the capillary waves.^36,37 According to the capillary wave theory, the system size dependence of the interface thickness λ is expressed to be

where C is a system-size independent constant for a fixed temperature. Note that the temperature T and system-size L used here are dimensionless values. The system-size dependence of the interface thickness for T = 0.9 is shown in the inset of Fig. 8(a). From this finite-size dependence, we estimate the surface tension as a function of temperature. The critical behavior of the surface tension obtained from the interface thickness is shown in Fig. 8(b). The critical exponent of the correlation length is determined to be ν = 0.63(4) which is consistent with the value of the Ising universality class.

We determine the phase boundary between gas and liquid phases by achieving gas-liquid coexisting state with MD. This method is simple and provides the phase boundary between liquid and gas phases with high accuracy. Since we have gas-liquid interface in the system, we can check the finite-size effect on the obtained results; if the linear length of the system L_z is much longer than the interface thickness λ, then the coexisting densities of liquid and gas are expected to be free from the finite-size effect. In our experience, the finite-size effect on the coexisting density is vanished when L_z > 10λ. We determine the critical temperature from the block density analysis of the Binder parameter, which allows us to perform finite-size analysis only from one simulation run with large system size. Additionally, this method does not require the distribution of order parameter of the Ising model at the criticality, and therefore, we do not have to consider the mixed-field theory.

The critical exponents β and ν are determined from the behaviors of the order parameter and the surface tension. While the critical exponents of the order parameter β is consistent with the Ising universality class, the exponent of correlation length ν obtained from the surface tension is found to be slightly larger than that of the Ising universality. While this deviation can be due to the capillary waves, it is difficult to remove the effect since the surface tension does not exhibit apparent finite-size effect. Potoff and Panagiotopoulos studied the critical behavior of the surface tension using HR technique and FSS analyses, and reported ν = 0.71(4) which is much larger than the Ising universality.³⁸ However, they did not consider the effect from the capillary waves. Considering the effect of the capillary waves, we obtained the critical exponent of the correlation length ν = 0.63(4) which is consistent with the Ising universality. Note that the statistical error became larger than that of the value obtained directly from the stress tensor using Eq. (4). This is because the surface tension analysis with considering the capillary waves requires two steps; first measuring the surface thicknesses for various system sizes and then extracting the surface tension from its size dependence. While the obtained values of exponents are consistent with that of the Ising universality class, it is still difficult to declare that the universality class of the LJ system belongs to the Ising universality since the accuracy is insufficient especially for the exponent of the correlation length. Some techniques are required which can estimate the surface tension accurately. Also the distances from the criticality is insufficient enough to discuss the universality class precisely. In order to get closer to the critical point, larger and longer simulations are necessary. It should be one of the further issues.

Some of the past studies reported the effect from the correction to scaling.^20,21 They first obtained the apparent critical temperatures for several system size, and they took the infinite size limit considering the correction. However, we did not find any sign of the correction to scaling in Fig. 6. It is difficult to determine the value of the correction to scaling exponent if it exists. While the obtained values of critical exponents are consistent with the Ising universality without considering the correction, it would become more important in order to investigate more accurately.

Very recently, Ma and Hu^39–43 used MD to study relaxation and aggregation of polymer chains, in which neighboring monomers of a polymer chain are connected by rigid bonds^39,41 or springs with various spring constants.^40,42 They found that when the bending-angle depending potential and the torsion-angle depending potential are zero or very small, polymer chains tend to aggregate.^41–43 Such results are useful for understanding the mechanism of protein aggregation^41,42 related to neurodegenerative diseases, including Alzheimer's disease, Huntington's disease, and Parkinson's disease. One may argue that polymer chains can aggregate only when they are inside the phase boundary. It is of interest to extend the method of the present paper to calculate the phase diagram to calculate phase diagrams of polymer chains of various chain lengths.⁴⁴ Such phase diagrams will be useful for understanding under what temperature and density conditions, polymer chains would aggregate.

The computation was carried out by using facilities of the Supercomputer Center, Institute for Solid State Physics, University of Tokyo, and Information Technology Center, Nagoya University. We would like to thank K. Binder, N. Kawashima, S. Todo, and Y. Tomita for helpful discussions. This work was partially supported by Grants-in-Aid for Scientific Research (Contract No. 23740287), by KAUST GRP (KUK-I1-005-04), by Grants NSC 100-2112-M-001-003-MY2, and by NCTS (North).