Fluid–Solid Phase Transition in Molecular Layers Adsorbed on a

Article Cite This: J. Phys. Chem. C XXXX, XXX, XXX−XXX

pubs.acs.org/JPCC

Fluid−Solid Phase Transition in Molecular Layers Adsorbed on a Smooth Surface: A New Insight from Molecular Simulations Eugene Ustinov*

J. Phys. Chem. C Downloaded from pubs.acs.org by UNIV OF SUNDERLAND on 10/10/18. For personal use only.

Ioffe Institute, 26 Polytechnicheskaya, St. Petersburg 194021, Russian Federation ABSTRACT: Two-dimensional (2D) phase transitions in molecular layers adsorbed on various surfaces have attracted much research attention over the past decades. Special attention is paid to the solidification/melting transition using methods based on molecular simulation techniques. The main research interest is focused on parameters of the liquid−solid or fluid−solid transitions; that is, the pressure, the chemical potential, and the density change. It is implied that the liquid and solid phases coexist like vapor and liquid with a flat interface, and the transition occurs simply via gradual increase of the solid fraction at a constant pressure. However, visualization of the transition based on molecular simulation in a canonical Monte Carlo method in a relatively large simulation cell has shown that the liquid and solid phases form islands in the liquid or “lakes” in the solid of arbitrary size and shape without a distinct interface. For this reason, the aim of this study was a detailed analysis of the 2D phase transition with the emphasis on the properties of the intermediate state between the pure liquid (fluid) and the crystal. Thermodynamically, the transition phase behaves like a two-phase system: when the density changes over a certain range at a specified temperature, the tangential pressure is nearly constant. However, no energy barrier was observed at the transition from the pure liquid or fluid to the intermediate phase. The absence of the van der Waals loop provides a very efficient way to determine parameters of the 2D liquid/fluid−solid transition and the chemical potential of the crystalline phase by a thermodynamic integration or direct evaluation with the kinetic Monte Carlo simulation in a sufficiently large simulation cell. The advantage of this way is that it makes the liquid−solid transition reversible without resorting to any special procedures. origin of the “insipient triple point” phenomenon.1,13,14 In some cases, the surface corrugation causes a commensurate− incommensurate phase transition in the adsorbed molecular layer.11,13,15,16 Even if the adsorbed monolayer does not form a commensurate structure, the surface corrugation affects the contact layer properties. Thus, X-ray measurements show a continuous change of the rotational angle of the crystalline lattice of argon relative to that of graphite in the temperature range 44−50 K.17 Some authors18 argue that the experimentally observed heat capacity peak in the system Ar−graphite at 47.2 K2 has to be attributed to the rotational transition of the crystalline argon molecular layer relative to the graphite lattice rather than to the melting transition of about 50 K. Unlike experimental techniques, molecular simulation methods allow the exclusion of some secondary complicating details to evaluate the fundamental features of adsorption systems. One of the most intriguing phenomena is the melting/ freezing transition in 2D contact layers. A liquidlike and crystalline monolayer on a smooth planar surface can be easily modeled with Monte Carlo (MC) or molecular dynamics (MD) methods. However, determination of the liquid−solid phase transition parameters requires special techniques. As is known, the transition occurs when the pressure and chemical potential in both phases

1. INTRODUCTION Molecular layers of gases adsorbed on a smooth surface at sufficiently low temperatures represent clear examples of twodimensional (2D) systems. The behavior of such systems is qualitatively similar to that of their three-dimensional (3D) counterparts but has some distinctions, which have received much interest over the past few decades. The liquid−solid and the fluid−solid transitions have received particular attention. Experimentally, the phase transition is observed using calorimetric measurements of the heat capacity change in the adsorbed phase,1−3 the differential heat of adsorption,4,5 and the adsorption isotherm analysis.6−8 Some investigations of the 2D transitions are based on the neutron diffraction technique9 and X-ray scattering.10,11 A description of other scattering and resonance techniques applied to the analysis of adsorbed molecular layers can be found elsewhere.12 Experimental data have shown that, at some conditions, the 2D liquid or supercritical fluid transform into the crystalline molecular layer. A serious complication is that, in many cases, the transition is strongly affected by a periodic modulation of the surface potential due to the atomic structure of the solid. For this reason, the crystalline contact layer can be commensurate or incommensurate with the solid surface lattice. In the first case, each atom or molecule of the molecular layer is trapped in a local potential minimum (for example, in the center of the graphite hexagon) enhancing the melting temperature of the layer, sometimes even above the 2D critical temperature. This is the © XXXX American Chemical Society

Received: August 8, 2018 Revised: September 27, 2018

A

DOI: 10.1021/acs.jpcc.8b07698 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

The tangential pressure was determined with the Irving− Kirkwood scheme.32 Verification of the results on modeling of the liquid/fluid− solid transition was accomplished with argon adsorption isotherms simulated at various temperatures in an elongated cell with an external potential imposed on the gas phase.27−30 The liquid/fluid−solid transition is identified by a substep in each isotherm due to the increase in the phase density, so the tangential pressure of the intermediate zone during the phase transition can be determined independently.

are the same. However, the difficulty is that the chemical potential of the crystalline phase uniformly distributed over the simulation cell cannot be directly evaluated as a function of density and temperature. This is because the pressure in the crystal at a specified temperature is not a single-valued function of its density as it explicitly depends on simulation box size, unless the box volume is extremely large. To circumvent the above restriction, some variants of thermodynamic integration have been developed that involve a number of tricks to make a phase transition reversible. This can be achieved using the idea of the Einstein crystal,19 the Wigner−Seitz cell method,20,21 and the λ integration method.22 If a reference state for the liquid−solid equilibrium is reliably established, the coexisting line can be evaluated by integration of the Gibbs−Duhem integration.23,24 In the group of “direct” methods, the liquid−solid interface is accounted explicitly.25,26 This approach had been further developed in the framework of the kinetic Monte Carlo (kMC) method by imposing a nonuniform external potential field on the system.27−30 An advantage of the technique is that the gas−liquid or gas− solid equilibria can be shifted by variation of the external field that allows one to determine the pressure−chemical potential curve for each coexisting phase. The point of intersection of the curves corresponds to the parameters of the phase transition. While the parameters of the liquid−solid transition can be evaluated with various techniques, not much is known on the structure and thermodynamic properties of the intermediate liquid−solid phase. Of particular interest is whether the transition occurs via formation of a peculiar homogeneous phase or by gradually replacing one pure phase by another: whether the van der Waals loop in the 2D liquid−solid transition exists or it is an artifact of a small simulation cell. This is an important question because, in the latter case, the liquid−solid transition parameters can be determined by a direct modeling of the intermediate phase in a sufficiently large simulation cell. These topics are addressed in this study.

3. RESULTS AND DISCUSSION 3.1. Equation of State for the 2D Liquid and Solid on Graphite. Simulations have been fulfilled for 400 particles in a canonical ensemble in a rectangular cell with the usual periodic boundary conditions. The number of MC steps was 50 × 106 both for the equilibration and sampling stages. The dimension Lz normal to the graphite surface was taken as 5σff. The ratio Ly/Lx was kept equal to 31/2/2. Figure 1 presents the tangential pressure (a) and internal energy (b) as functions of temperature along constant values of density for the crystalline monolayer. Solid lines in Figure 1 are fit with the following equations:30,31 n1

u = cT * +

n2

∑ ∑ (1 − j)bjk(T*) j ρk (1)

j=0 k=0

pT = ρT* +

n1

n2

∑ ∑ kbjk(T*) j ρk+ 1 (2)

j=0 k=1

Here, T* = RgT (kJ/mol), u is the potential energy, including the fluid−solid potential (kJ/mol), pT is the tangential pressure (mN/m), ρ is the 2D density (μmol/m2), and Rg is the gas constant. Equations 1 and 2 are in agreement with the Maxwell relation,33 so the functions u and pT are not independent:

2. METHODOLOGY AND SIMULATION DETAILS Simulations were performed for the system Ar−graphite with the kMC method27−30 in a canonical ensemble with the usual periodic boundary conditions. The liquidlike and crystalline molecular layers were modeled for 400 atoms at a number of specified densities over a wide temperature range. The Lennard−Jones (LJ) parameters for argon were taken as 0.3393 nm and 118.4 K for the collision diameter, σf f, and the potential well depth, εf f/kB (kB is the Boltzmann constant), respectively.27 The corresponding LJ parameters for the carbon atom were taken as 0.26 nm and 49.933 K.29 The p−ρ−T diagrams of the liquid and solid phases were generalized by regression equations similar to those applied to the nitrogen− graphite system.31 The idea was that, if the chemical potential is known in one point, it can be easily recalculated with the regression equation for any other temperature and density of the same phase. This is a way to prove that the chemical potential of the liquidlike phase evaluated by a direct method with the kMC scheme at any density and temperature is thermodynamically correct and consistent with the Gibbs−Duhem equation. The intermediate zone between the liquid and solid was simulated at a relatively large number of argon atoms, N = 2500. The number of cycles was 4 × 104 and 8 × 104 for equilibration and sampling, respectively. The structure of the intermediate phase was visualized in this study with contour lines that present the monolayer density distribution averaged over the sampling stage of each run.

ρ2

∂p /T ∂u + T2 T =0 ∂ρ ∂T

(3)

The coefficients bjk have been determined by the least-squares fitting for n1 = n2 = 3. The constant c is the heat capacity in units of kB at the limit of zero temperature. It turned out that c = 1.507 as a result of the fitting parameters of eqs 1 and 2, which suggests that its true value is 1.5. By taking into account the kinetic energy and three degrees of freedom, the total heat capacity in the limiting case is 3kB, which is in exact conformity with the monatomic gas theory in the classical limit.30 Note that the idea of harmonic oscillations of atoms in the crystal lattice in this case is unnecessary. Similar results obtained for the liquidlike argon monolayer are shown in Figure 2. As seen in Figures 1 and 2, the simulated values of pressure and internal energy are excellently fitted by regression eqs 1 and 2. Naturally, the coefficients are different for the liquid (fluid) and solid, but the constant c = 1.5 is taken the same in both cases. The coefficient b10 is the constant of integration, and it cannot be determined using only the p−ρ−T data, but it appears in the expression for the chemical potential: μ = T*ln(ρ) − cT *ln(T*) +

n1

n2

∑ ∑ (k + 1)bjk(T*) j ρk j=0 k=0

(4)

If the chemical potential is known at least at one point of the p−ρ−T diagram, all its other values can be determined with eq 4. B

DOI: 10.1021/acs.jpcc.8b07698 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 1. Effect of temperature on (a) the tangential pressure, pT, and (b) the specific internal energy, u, in the crystalline argon monolayer on graphite along constant values of density (numbers at curves). Surface density is incremented by 0.1 μmol/m2. Solid lines are fits with eqs 1 and 2.

Figure 2. Dependence of (a) the tangential pressure, pT, and (b) the specific internal energy, u, on temperature in the liquidlike argon monolayer on graphite along constant values of density (numbers at curves). Solid lines are fits with eqs 1 and 2. The 2D density is incremented by 0.1 μmol/m2.

Figure 3. (a) Helmholtz free energy and (b) entropy of the contact liquidlike argon layer on graphite versus surface density at different temperatures (numbers at curves).

as [u(ρ, T) − f(ρ, T)]/(kBT). The standard state corresponds to that of the ideal gas at the same temperature. The absence of systematic deviations of fitted lines from points evaluated with the molecular simulation confirms thermodynamic consistency of the results obtained with the kMC method. Note that the entropy of the liquidlike contact layer is reduced with the decrease of temperature and the increase of the 2D density that indicates the growth of ordering in the layer structure. 3.2. Intermediate Zone between the 2D Liquid and the Solid. As the density of the contact layer increases at constant temperature, the tangential pressure passes through the inflection point and then approaches a constant value, which is demonstrated in Figure 4a. The evident reason is the appearance of the crystalline phase.

This can be easily done for the liquid phase using the kMC scheme that allows for evaluation of the chemical potential directly at any density. Once the chemical potential is known, the Helmholtz free energy and the entropy change can be easily determined. The coefficients bjk in eqs 1, 2, and 4 including b10 have been determined by the least-squares minimization of the residual that includes deviations of the chemical potential, the pressure, and the internal energy together as functions of density at several values of temperature. Figure 3 presents the Helmholtz free energy and the entropy versus density of the contact layer at temperatures within the range 55−75 K. The Helmholtz free energy, f(ρ, T) was determined as the difference μ − p/ρ. The entropy change ΔS/kB was then calculated C

DOI: 10.1021/acs.jpcc.8b07698 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 4. (a) Tangential pressure and (b) the chemical potential simulated as a function of density for the 2D liquid, fluid, and the transition zone (circles) and crystal (crosses) at fixed temperatures over the range 55−75 K in the system Ar−graphite. In the density range below 11.5 μmol/m2 and for the crystalline phase the number of molecules was taken as 400. Transition zones, pure liquid, and fluid at a density larger than 11.5 μmol/m2 were simulated at N = 2500. Empty circles in panel (a) show the results of the simulation in a small simulation cell with 256 particles at 70 K.

always melts at equilibrium. Therefore, in the sufficiently large system, the 2D liquid−solid transition occurs without hysteresis despite the transition being of the first order. It should also be noted that a direct modeling of the liquid− solid coexistence in the same cell with a distinct flat interface is hardly possible in the case of the 2D system. Even if the crystal and liquid are initially separated, these phases inevitably form a mixture at a sufficiently long chain of MC steps. However, the density, pressure, and chemical potential of the crystal corresponding to the phase transition at a given temperature can be determined by the point of intersection of the corresponding pT−ρ curves. Evolution of the argon contact layer structure at 70 K at different average densities is presented in Figure 5. At a density 11.85 μmol/m2 (Figure 5a), the Ar monolayer is a homogeneous dense fluid without a long-range order. Panels (b), (c), and (d) show gradual change of the intermediate phase structure with the increase of its average density. At a relatively low density (12.0 μmol/m2), it looks like a glass with small crystalline inclusions. Further increase of the density results in the formation of hexagonally packed crystalline domains separated from each other by a flexible phase which has an ordered structure resembling parallel linear chains of atoms. This is clearly seen in panels (b) and (c) of Figure 5. Such a stringlike arrangement means that atoms can move predominantly along one direction that is parallel to the symmetry axis of the nearest crystalline domain. In its turn, this means that this flexible phase is nonautonomous and does not exist individually. It should be stressed that the solid−fluid interface in the intermediate zone has a quite arbitrary shape rather than straight or circular that points to a small surface tension of the interface. At a relatively low average density, the intermediate phase resembles crystalline islands of arbitrary shape and size in the fluid as seen in Figure 5b, while at a higher density the phase is viewed as a number of liquidlike lakes in the crystal (see Figure 5d). In view of the relatively small contribution of the surface tension to the Helmholtz free energy of the intermediate phase, one can suggest that thermodynamically this phase behaves as the solid− fluid system where each fraction maintains properties of the corresponding pure phase. 3.2.1. Kinetic Monte Carlo Scheme versus the Metropolis Algorithm. The chemical potential presented in Figure 4b was simulated with the kMC method.30 This raises the question of whether the same results could be obtained with the standard Metropolis MC (MMC) scheme. To clarify this issue, the

It is worth noting that a metastable state of the 2D liquid (fluid) is not observed at all temperatures below and above the 2D critical temperature (≈63 K), which is an indication of a small energy barrier for nucleation. In part, this can be explained by a relatively large number of particles (N = 2500) and MC cycles, which was not less than 80 000. The plateau in each pressure−density curve (Figure 4a) exactly corresponds to that for the chemical potential within the same density range (see Figure 4b), which agrees with the Gibbs−Duhem equation. Figure 4a shows a striking and fundamental result that follows from the size effect. The pressure−density curve at 70 K (empty circles) was additionally simulated in a relatively small simulation box comprising only 256 LJ particles, which is 1 order of magnitude less than the number of particles taken for simulation of the intermediate zone in this study. In this case, the transition zone is highly diffused, which is always interpreted as an indication of thermodynamic instability of the system. Note that the number of MC cycles in this case was 1 × 106; that is, much higher than that taken for the larger simulation cell. The fitting curve (the dashed line in Figure 4) has a typical form of the van der Waals loop making an erroneous representation that the 2D liquid−solid transition occurs via a nucleating mechanism, which implies that there is a metastable overcompressed 2D liquid or fluid. One could think that the true transition pressure is somewhere between the liquidlike and crystal-like spinodal points; that is, in the range 20.5−23.0 mN/m. A sign that it is not the case is clearly seen in Figure 4a: the false liquidlike spinodal point does not exceed the pressure in the intermediate zone determined by simulation in the large box with 2500 LJ particles. Several runs have been performed starting from the distribution of random particles along the plane parallel to the graphite surface. However, at high densities, large crystalline domains in the monolayer rapidly appeared spontaneously. This suggests that the 2D transition from liquid (or fluid at a temperature higher than the critical temperature) to crystal is equilibrium; that is, the van der Waals loop appears as an artifact in systems that are too small. The metastable state of the crystalline phase does exist as seen in Figure 4a, but again, only when the simulation box is sufficiently small. The lines in Figure 4 corresponding to the crystalline phase were simulated for a system of 400 argon atoms on the graphite surface. It has been convincingly shown27−30 that, if the crystal coexists with the gas in an elongated cell and the interface is not curved, the crystal D

DOI: 10.1021/acs.jpcc.8b07698 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 5. Density distribution in the system Ar−graphite at 70 K. Average densities (μmol/m2) are (a) 11.85, (b) 12.0, (c) 12.15,and (d) 12.30. The contour lines have been plotted using sets of argon atom positions averaged over 2 × 108 kMC steps. The number of atoms is 2500. Panels (b), (c), and (d) refer to the intermediate zone.

Figure 6. (a) Tangential pressure and (b) the chemical potential in the system Ar−graphite simulated with the standard MMC coupled with the Widom test particle insertion method and the kMC algorithm at 60 and 70 K. Lines and crosses in panel (a) show the tangential pressure of crystal simulated with MMC.

nearly identical results for the pure liquid or fluid. Yet, within the plateau region, statistical noise of the pressure and chemical potential in the case of the standard MC scheme is markedly larger. This is a consequence of different ways of averaging. The kMC method realizes the time averaging over configurations that have different lifetimes. If, for example, one particle randomly takes energy much larger than the average energy of the system, the duration of that specific configuration is very small and therefore does not contribute much to the fluctuations of any parameter. In the conventional method based on the ensemble

tangential pressure and the chemical potential at 70 and 60 K, above and below the 2D critical temperature, were additionally calculated with the conventional Metropolis scheme coupled with the Widom test particle insertion method.34 A comparison of the results obtained with the two methods is presented in Figure 6. Simulations with the two methods were performed under equal conditions. In both cases, the number of particles and the number of cycles at the sampling stage were taken as 2500 and 8 × 104, respectively. As seen in Figure 6, both schemes provide E

DOI: 10.1021/acs.jpcc.8b07698 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 7. (a) Tangential pressure and (b) internal energy in the intermediate zone of the system Ar−graphite at constant densities (numbers in legends). Solid lines show the prediction using the idea of the liquid (fluid)−solid mixture without contribution of the interface. Gray lines in panel (a) correspond to the pressure−temperature dependence in the crystal at constant densities (numbers at curves in μmol/m2).

densities that correspond to the transition zone. The fact that those lines are extended to higher temperatures above the 2D triple point is explained by a relatively small number of argon atoms (N = 400) in the simulation box. At large numbers of particles, the metastable section of each curve substantially shortens suggesting that it fully disappears at the limit of an infinitely large cell. The internal energy of the transition zone is shown in Figure 7b. At a specified temperature, the internal energy slightly decreases with the increase of the average density due to the growth of the crystalline fraction. The solid lines are predicted by the approach. The correspondence of the simulated and predicted values of the internal energy is excellent in view of its very small change within the transition zone, which does not exceed 150 J/mol at a specified temperature. The fraction α of the crystalline phase in the intermediate zone can be determined by two ways. The first way is as follows. Once the tangential pressure at a given density ρ and a temperature is known by simulation, the densities of the pure liquid, ρL, and solid, ρS, phases can be calculated with eq 2. Then, the fraction α is

average, all configurations have equal weight. This increases the contribution of infrequent events to the averages in relatively small systems where the number of particles does not exceed several thousands. 3.3. Thermodynamic Properties of the Intermediate Zone. The main question is whether the intermediate zone can be thermodynamically viewed as a mixture of solid and liquid as individual phases. If so, the system can be described using eqs 2 and 4 for the pressure and chemical potential, respectively. At equilibrium, these variables for the solid and liquid are equal to each other. Since the coefficients bjk for the 2D liquid and cjk for the crystal are already known, the density of each phase can be determined by solving eq 2 iteratively at a given temperature and simulated pressure. By using these values and the condition of equality of partial chemical potentials, in the next step, one can determine the difference of the coefficients c10 and b10, which cannot be evaluated only with the p−ρ−T diagram. This difference, δ10 = (c10 − b10), must be the same at any temperature and average density of the intermediate phase, so it can be averaged over the array of p−ρ−T points simulated for the intermediate phase. The value of δ10 is not unique due to the interplay between δ10 and other constants. However, at a given set of coefficients bjk and cjk evaluated with the p−ρ−T diagram, the parameter δ10 appeared to be very stable. No correlation of δ10 with the temperature and density within the transition zone was found. The relative standard deviation of δ10 over 226 points proved to be 0.12 × 10−3, which is an extremely small value. The difference in the chemical potentials of the two phases is proportional to δ10T, so the condition of its equality does not need to know b10 and c10 separately. Once the coefficients bjk, cjk, and δ10 are known, the pT−T coexisting curve can be predicted using the condition of equality of the pressure and chemical potential in coexisting phases at any temperature. Figure 7 presents the simulated tangential pressure and internal energy versus temperature at constant densities within the intermediate zone. As seen in Figure 7a, 188 points of the tangential pressure in the intermediate liquid−solid phase fall onto a single black line simulated with eqs 2 and 4 using the condition of equality of partial pressures and chemical potentials and the assumption that the 2D liquid (fluid) and solid fractions in the transition zone have the same properties as those of the corresponding pure phases. The effect of the liquid−solid interface on the pressure is assumed to be insignificant. Gray lines indicate the pT−T dependence for the pure crystalline phase at the same

α = (ρ − ρL )/(ρS − ρL )

(5)

It is assumed that coefficients bjk and cjk for the liquid and solid, respectively, are known except b10 and c10. In this case, there is no need to know the chemical potential of the intermediate phase. According to the second way, the densities ρL and ρS can be calculated from two conditions: equality of chemical potentials and pressures for the pure liquid and solid. This is possible when the difference δ10 = (c10 − b10) has already been determined. Figure 8 presents the dependence of fraction α on the transition zone density along constant temperatures calculated with the two methods. Some deviations in the crystalline fraction determined by the two methods probably stem from neglecting the contribution of the liquid−solid interface to the Helmholtz free energy and the assumption that the liquid fraction in the intermediate zone has the same thermodynamic properties as the pure 2D liquid phase. It should be noted that the use of any “direct” methods based on order parameters or energy distribution to evaluate the liquid and solid fractions in the transition phase are highly unreliable, which is evident from Figure 5. 3.4. Argon Adsorption Isotherms at Low Temperatures. This section aims at a comparison of the results F

DOI: 10.1021/acs.jpcc.8b07698 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 8. Fraction α of the crystalline phase in the intermediate zone as a function of its density at constant temperatures. Circles denote values determined by the tangential pressure and virial eq 2 for the pressure of pure liquid and solid phases. Solid lines are predicted dependences corresponding to the condition of equality of the chemical potential and tangential pressure for pure liquid and solid phases.

Figure 9. Selected argon adsorption isotherms on graphite simulated with a canonical kMC in an elongated cell at different temperatures. The dashed line is the experimental Ar−graphite adsorption isotherm.35

substep that corresponds to the 2D liquid/fluid−solid transition. The bulk pressure has been determined from the chemical potential assuming that the bulk phase is close to the ideal gas due to extremely low pressures. The chemical potential and the tangential pressure at the substeps have been determined quite reliably, which enables the comparison of the parameters of the liquid/fluid−solid transition with those evaluated in the previous section. This can be done taking into account that the chemical potential for the pure liquid or fluid phase is already known, which means that the coefficient b10 is also known. On the other hand, the difference b10 − c10 has been evaluated independently from properties of the intermediate zone, which eventually determines the coefficient c10. This completes formulation of the model based on virial expansion of the chemical potential in the form of eq 4 that enables the description of both liquid/fluid and crystalline phases. The only question left is whether the coefficient b10 determined for the pure liquidlike monolayer can be used for the intermediate zone. Figure 10 presents the temperature dependence of the tangential pressure and chemical potential in the intermediate transition zone determined (i) by a direct simulation, (ii) analytically with eqs 2 and 4, and (iii) from the adsorption isotherms simulated in the elongated cell with the external potential. Figure 10 demonstrates excellent correspondence between the tangential pressure and chemical potential determined by three different ways. The important result is that the chemical potential in the intermediate zone evaluated by a direct kMC simulation coincides with that for the substeps on adsorption isotherms. The intermediate zone is not only dense but also includes the crystalline phase. Nevertheless, the kMC method provides reliable simulation of the chemical potential of the dense phases. Among other things, the kMC method does not need any external procedures or auxiliary methods to determine the chemical potential, making it promising in modeling phase transitions in more complex systems.

obtained by simulation of the argon molecular layer uniformly distributed over the xy plane of the simulation cell with those obtained by simulation in an elongated cell with a positive external potential imposed on its ends.27−30 The external potential splits the system into gas in the vicinity of the cell ends and a dense fluid, liquid, or solid in the central part of the cell. The total chemical potential of the system that includes the external potential is the same over the cell, but the intrinsic chemical potential of the gas phase is lower than that of the dense phase. Variation of the external potential shifts the gas− solid or gas−fluid equilibria providing a possibility to model the adsorption isotherm over a wide range of chemical potential. The advantage of this technique is that the presence of the gas phase allows for evaluation of the chemical potential even with the standard MMC scheme coupled with the Widom particle insertion method. Having determined the chemical potential of the gas phase, we thereby determine one for the dense phase, because the chemical potential of the coexisting phases is the same. An additional advantage of this technique is that, in the case of the crystalline phase in the center of the cell, the tangential pressure components can be easily equilibrated by tuning the cell width on the equilibration stage of the run. This ensures a thermodynamically correct crystalline density corresponding to the chemical potential of the system at a given temperature. Several selected argon adsorption isotherms simulated with a canonical kMC scheme using the above technique are presented in Figure 9. There is an important feature of the low-temperature adsorption isotherms shown in Figure 8 that should be emphasized. The amount adsorbed increases with the chemical potential even when the second and further molecular layers are nearly absent and the contact layer is crystallized. Thus, at 40 K, the 2D monolayer density grows by 10% with the bulk pressure from the point of the gas−solid transition to the onset of the second layer. This is a requirement of the Gibbs−Duhem equation: the increase of the chemical potential results in the increase of the tangential pressure, but the latter is only possible with the density growtha. Starting from 50 K, that is, slightly above the 2D triple-point temperature, each simulated adsorption isotherm exhibits a

4. CONCLUSIONS In the present study, the crystallization phenomenon in the molecular layer of argon adsorbed on graphite was investigated by means of the kMC method. Special attention was paid to the structure and thermodynamic properties of the intermediate G

DOI: 10.1021/acs.jpcc.8b07698 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 10. (a) Tangential pressure and (b) chemical potential on the line of liquid/fluid−solid coexistence. Crosses designate direct simulation with the kMC scheme in the series presented in Figure 4. The solid lines indicate predictions with eqs 2 and 4 at equality of the chemical potential and pressure of the coexisting phases. White circles correspond to (a) the tangential pressure and (b) the chemical potential at the substeps on the argon adsorption isotherms presented in Figure 9.

phase appearing at the fluid−solid transition. Several new results have been established, which can be summarized as follows. • The tangential pressure of the 2D liquid or fluid grows with the 2D density at a specified temperature up to a maximal value followed by a plateau. The van der Waals loop is absent and appears as an artifact when simulation is accomplished in a relatively small simulation cell. • The intermediate zone corresponding to the pressure− density plateau is viewed as crystalline islands in liquid or liquid “lakes” in the crystalline phase at a higher average density without a clear interface of a distinct shape. This is an indication of a very small liquid−solid surface tension, which explains the absence of the metastable overcompressed liquid due to a small nucleation barrier. • The structure of the liquid phase in the intermediate zone has a linear arrangement suggesting that this phase is not autonomous. Nevertheless, the intermediate phase as a whole behaves as a two-phase liquid−solid system, with the thermodynamic properties of the liquid and solid being close to those of individual phases. This feature allows one to predict the pressure−temperature coexisting line. • The kMC method enables reliable evaluation of the chemical potential of the 2D fluid, the liquid, and the intermediate zone at any density. As a whole, this work has provided new insight into the microscopic origin of the 2D liquid−solid transition owing to simulations in a relatively large simulation box with the kMC method. The outcome of the approach is a relatively simple scheme of determination of the phase transition parameters and the thermodynamic properties of the 2D liquid, the fluid, and the solid.

■

ensemble when the adsorbed phase is crystallized. If the simulation box is unchanged, then the monolayer density and tangential pressure remain the same at any given chemical potential. This is because the lattice constant and the box size are strongly inter-related. Therefore, any simulated adsorption isotherm that includes sections with a constant amount adsorbed is thermodynamically incorrect.

■

REFERENCES

(1) Migone, A. D.; Chan, M. H. W.; Niskanen, K. J.; Griffiths, R. B. Incipient triple point for N2 adsorbed on graphite. J. Phys. C: Solid State Phys. 1983, 16, L1115−L1120. (2) Migone, A. D.; Li, Z. R.; Chan, M. H. W. Melting transition of submonolayer Ar adsorbed on graphite. Phys. Rev. Lett. 1984, 53, 810− 813. (3) Chan, M. H. W.; Migone, A. D.; Miner, K. D.; Li, Z. R. Thermodynamic study of phase transitions of monolayer N2 on graphite. Phys. Rev. B: Condens. Matter Mater. Phys. 1984, 30, 2681− 2694. (4) Rouquerol, J.; Partyka, S.; Rouquerol, F. Calorimetric evidence for a bidimensional phase change in the monolayer of nitrogen or argon adsorbed on graphite at 77 K. J. Chem. Soc., Faraday Trans. 1 1977, 73, 306−314. (5) Grillet, Y.; Rouquerol, F.; Rouquerol, J. Two-dimensional freezing of nitrogen or argon on differently graphitized carbons. J. Colloid Interface Sci. 1979, 70, 239−244. (6) Hoja, R.; Marx, R. High-resolution adsorption isotherm measurement of adsorbed krypton near the commensurate-incommensurate transition. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 34, 7823−7827. (7) Zhang, Q. M.; Larese, J. Z. Melting of monolayer argon adsorbed on a graphite substrate. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43, 938−946. (8) Choi, B. I.; Nham, H. S.; Youn, H. S. Phase diagram of a physisorbed krypton monolayer on graphite using an ellipsometric technique. J. Korean Phys. Soc. 2008, 53, 3262−3266. (9) Taub, H.; Carneiro, K.; Kjems, J. K.; Passell, L.; McTague, J. P. Neutron scattering study of 36Ar monolayer films adsorbed on graphite. Phys. Rev. B 1977, 16, 4551−4568. (10) McTague, J. P.; Als-Nielsen, J.; Bohr, J.; Nielsen, M. Synchrotron x-ray study of melting in submonolayer Ar and other rare-gas films on graphite. Phys. Rev. B: Condens. Matter Mater. Phys. 1982, 25, 7765− 7772. (11) Nielsen, M.; Als-Nielsen, J.; Bohr, J.; McTague, J. P.; Moncton, D. E.; Stephens, P. W. Melting studies of argon on ZYX graphite. Phys. Rev. B: Condens. Matter Mater. Phys. 1987, 35, 1419−1425.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]ffe.ru. ORCID

Eugene Ustinov: 0000-0003-4668-1326 Notes

The author declares no competing financial interest.

■

ADDITIONAL NOTE Special caution is required in simulation of low-temperature adsorption isotherms on a smooth surface in a grand canonical a

H

DOI: 10.1021/acs.jpcc.8b07698 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

(34) Widom, B. Some topics in the theory of fluids. J. Chem. Phys. 1963, 39, 2808−2812. (35) Gardner, L.; Kruk, M.; Jaroniec, M. Reference data for argon adsorption on graphitized and nongraphitized carbon blacks. J. Phys. Chem. B 2001, 105, 12516−12523.

(12) Vilches, O. E. Phase transitions in monomolecular layer films phisisorbed on crystalline surfaces. Annu. Rev. Phys. Chem. 1980, 31, 463−490. (13) Butler, D. M.; Litzinger, J. A.; Stewart, G. A.; Griffiths, R. B. Heat capacity of krypton physisorbed on graphite. Phys. Rev. Lett. 1979, 42, 1289−1292. (14) Hansen, F. Y. Effects on the structure of monolayer and submonolayer fluid nitrogen films by the corrugation in the holding potential of nitrogen molecules. J. Chem. Phys. 2001, 115, 1529−1537. (15) Stephens, P. W.; Heiney, P. A.; Birgeneau, R. J.; Horn, P. M.; Moncton, D. E.; Brown, G. S. High-resolution x-ray-scattering study of the commensurate − incommensurate transition of monolayer Kr on graphite. Phys. Rev. B: Condens. Matter Mater. Phys. 1984, 29, 3512− 3532. (16) Specht, E. D.; Mak, A.; Peters, C.; Sutton, M.; Birgeneau, R. J.; D’Amico, K. L.; Moncton, D. E.; Nagler, S. E.; Horn, P. M. Phase diagram and phase transitions of krypton on graphite in the extended monolayer regime. Z. Phys. B: Condens. Matter 1987, 69, 347−377. (17) D’Amico, K. L.; Bohr, J.; Moncton, D. E.; Gibbs, D. Melting and orientational epitaxy in argon and xenon monolayers on graphite. Phys. Rev. B: Condens. Matter Mater. Phys. 1990, 41, 4368−4376. (18) Flenner, E.; Etters, R. D. Properties of argon adlayers deposited on graphite from Monte Carlo calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 125419. (19) Frenkel, D.; Ladd, A. J. C. New Monte Carlo method to compute the free energy of arbitrary solids. Application to the fcc and hcp phases of hard spheres. J. Chem. Phys. 1984, 81, 3188−3193. (20) Hoover, W. G.; Ree, F. H. Use of computer experiments to locate the melting transition and calculate the entropy in the solid phase. J. Chem. Phys. 1967, 47, 4873−4878. (21) Hansen, J. P.; Verlet, L. Phase transitions of the Lennard-Jones system. Phys. Rev. 1969, 184, 151−161. (22) Grochola, G. Constrained fluid λ-integration: Constructing a reversible thermodynamic path between the solid and liquid state. J. Chem. Phys. 2004, 120, 2122−2126. (23) Kofke, D. A. Direct evaluation of phase coexistence by molecular simulation via integration along the saturation line. J. Chem. Phys. 1993, 98, 4149−4162. (24) Agrawal, R.; Kofke, D. A. Thermodynamic and structural properties of model systems at solid-fluid coexistence. Mol. Phys. 1995, 85, 43−59. (25) Tepper, H. L.; Briels, W. J. Crystallization and melting in the Lennrd-Jones system: Equilibration, relaxation, and long-time dynamics of the moving interface. J. Chem. Phys. 2001, 115, 9434−9443. (26) Pedersen, U. R.; Hummel, F.; Kresse, G.; Kahl, G.; Dellago, Ch. Computing Gibbs free energy differences by interface pinning. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 094101. (27) Ustinov, E. A. Phase equilibrium in argon films stabilized by homogeneous surfaces and thermodynamics of two-stage melting transitions. J. Chem. Phys. 2014, 140, 074706. (28) Ustinov, E. A. Improved modeling of two-dimensional transitions in dense phases on crystalline surfaces. Krypton − graphite system. J. Chem. Phys. 2015, 142, 074701. (29) Ustinov, E. A. Effect of crystallization and surface potential on the nitrogen adsorption isotherm on graphite: A refined Monte Carlo simulation. Carbon 2016, 100, 52−63. (30) Ustinov, E. A. Efficient chemical potential evaluation with kinetic Monte Carlo method and non-uniform external potential: LennardJones fluid, liquid, and solid. J. Chem. Phys. 2017, 147, 014105. (31) Ustinov, E.; Gorbunov, V.; Akimenko, S. From simulation to thermodynamics of orientational transitions in molecular layers: Nitrogen contact layer on solids. J. Phys. Chem. C 2018, 122 (5), 2897−2908. (32) Irving, J. H.; Kirkwood, J. G. The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics. J. Chem. Phys. 1950, 18, 817−829. (33) van der Hoef, M. A. Free energy of the Lennard-Jones solid. J. Chem. Phys. 2000, 113, 8142−8148. I

DOI: 10.1021/acs.jpcc.8b07698 J. Phys. Chem. C XXXX, XXX, XXX−XXX