Article Cite This: J. Phys. Chem. C 2018, 122, 2897−2908
pubs.acs.org/JPCC
From Simulation to Thermodynamics of Orientational Transitions in Molecular Layers: Nitrogen Contact Layer on Solids Eugene Ustinov,*,† Vitaly Gorbunov,‡ and Sergey Akimenko‡ †
Ioffe Institute, 26 Polytechnicheskaya, St. Petersburg 194021, Russian Federation Omsk State Technical University, 11 Pr Mira, Omsk 644050, Russian Federation
‡
S Supporting Information *
ABSTRACT: Orientational ordering is a well-known feature of polyatomic molecules in contact layers on solid surfaces. The majority of reported results refer to systems significantly depended on specific fluid−substrate interaction. An example is the nitrogen molecular layer commensurate with the graphite lattice at low temperatures. At the same time, it is instructive to analyze intrinsic thermodynamic properties of the molecular layer stabilized by a smooth substrate surface. In this study the nitrogen molecular layer is modeled by the Monte Carlo method with the emphasis placed on the orientational transition. The quadrupole−quadrupole interaction induces the herringbone orientational ordering which stretches the layer lattice along the glide line. The increase of temperature reduces the deformation, and its disappearance coincides with the transition from long- to short-range order. The magnitude of dispersion and electrostatic forces are always comparable even if the long-range ordering is absent, though the contribution of the latter to the internal energy is less than 17%. A new methodology of analysis of orientational transitions is developed accounting for the anisotropy of electrostatic forces inherent in the long-range orientational arrangement, which allows one to circumvent the problem of very small latent heat of the transition and provides its thermodynamic analysis. K.13,15,17,31−34 Melting of the adsorbed monolayer is believed to occur by the vacancy mechanism as a first-order phase transition. Moreover, the melting temperature increases with the surface coverage from 48 K for gas−crystal equilibrium to 82 K when the monolayer is completed.17 The situation with the orientational phase transition is quite ambiguous. Experimental measurements of the heat capacity of the adsorption layer, as well as the LEED and NMR results, are consistent with a first-order phase transition.35,36 On the other hand, theoretical approaches and computer simulations performed in the end of 1970s and over the 1980s supported both the first-order37−40 and the second-order41−43 orientational phase transition, with the temperature of the transition being significantly dependent on assumptions involved into the model. Besides, limited computational resources available at that time presented an additional challenge in reliable analysis of the orientational ordering phenomenon by a Monte Carlo method or molecular dynamics. Because of a very small entropy change and the latent heat of the orientational phase transition, information on its mechanism is still insufficient because it is mainly based on molecular simulations which are not statistically reliable or performed in small-sized systems. In the latest paper of Cai,29 advanced methods such as the histogram reweighting method, finite-size scaling, and the fourth-order cumulant method were used to ascertain the
1. INTRODUCTION A detailed understanding of driving forces and a mechanism of orientational ordering in molecular layers of nonspherical molecules possessing rotational degrees of freedom is an important task of physical chemistry.1−7 Representatives of such systems are molecular layers of diatomic molecules adsorbed on the surface of a solid. Properties of hydrogen,8−11 nitrogen,3,12−17 oxygen,18−23 and carbon monoxide3,24−26 adsorbed on the graphite surface have been studied experimentally using the low-energy electron diffraction (LEED), X-ray difraction (XRD), nuclear magnetic resonance (NMR), heat capacity, and adsorption isotherm measurements, which is comprehensively presented in the literature. A physical interpretation of the orientational phenomenon has attracted the attention of researchers for several decades. It has been well established that electrostatic forces caused, for example, by the quadrupole−quadrupole interaction play the key role in the orientational ordering at relatively low temperatures. The increase in temperature induces the orientational order− disorder phase transition at which the translational ordering of the crystalline phase remains, but the long-order orientational ordering disappears. The most studied adsorption system that clearly exhibits the orientational ordering and corresponding phase transitions is the contact layer of nitrogen on graphite.16,12,17,3,13−15 It is well-known that the nitrogen monolayer adsorbed on graphite undergoes an orientational phase transition from herringbone to orientationally disordered crystalline phase at temperature close to 28 K,3,12,13,15,16,27−30 which subsequently melts at a temperature above about 50 © 2018 American Chemical Society
Received: December 14, 2017 Revised: January 11, 2018 Published: January 24, 2018 2897
DOI: 10.1021/acs.jpcc.7b12300 J. Phys. Chem. C 2018, 122, 2897−2908
Article
The Journal of Physical Chemistry C
ular layer with a short-range orientational ordering. We rely on generalizing of the pVT data obtained by the Monte Carlo simulation in the form of a regression equation of state for those phases. For an indicator of the orientational phase transition, we take the point of loss of the tangential pressure anisotropy contributed by the quadrupole−quadrupole interaction. In this study the orientational transition in the uniform crystalline phase is modeled at a specified surface density and, therefore, constant Helmholtz free energy. The corresponding thermodynamic method and the algorithm adapted to this condition are considered in section 2.1. In the case of gas− solid/liquid coexistence on the graphite surface we relied on the kinetic Monte Carlo method50 and modeled the isothermal− isobaric orientational transition and the first-order melting transition. In this study we consider the orientational phenomenon without accounting for the quantum effects which play a significant role at low temperatures. Thus, it was shown52 that quantum fluctuations markedly reduce the orientational transition temperature. Nonetheless, it is instructive to analyze the phenomenon in the framework of classical statistical mechanics even at very low temperatures to distinguish the contribution of quantum effects depending on the temperature and molecular mass. Additionally, we consider an idealized case of rigorously twodimensional molecular layer and compare its properties including the orientational ordering with those of the contact molecular layer on the graphite surface. The latter looks like a 2D layer but is in fact a three-dimensional system, so it is of theoretical interest whether the experimentally observed molecular layer can be modeled in the framework of the twodimensional model.
orientational herringbone phase transition within the framework of anisotropic planar rotor model. This phase transition was shown to belong to the same universality class as the wellknown 3-state Potts model. This means that the orientational phase transition is a second-order transition. However, in the subsequent comprehensive work by Marx et al.,12 this conclusion was questioned on the basis of analysis of the correlation length near the phase transition and the forth-order cumulants of the order parameter. In particular, it has been established that the correlation length is finite, though unusually large. Based on these data, the orientational phase transition from the herringbone structure to the orientationally disordered crystal phase was concluded to be a fluctuation driven weak first-order transition.12,44 An obstacle to accurately determining the herringbone transition temperature and its type, in addition to the reasons mentioned above, is also related to the fact that the transition to the orientationally disordered state is smeared and occurs in the temperature range of 6−7 K.16 It is confirmed by the LEED mesurements, where a sharp decrease in the intensity of the diffraction pattern is observed between 25 and 29 K and then the intensity decreases slowly up to the melting point.45,46 Analysis of this feature of the phase transition by theoretical and computational methods supports a viewpoint that the short-range orientational order maintains above the herringbone transition.12,16,47 Thus, the only criterion by which the true temperature of the orientational phase transition can be determined is the loss of long-range orientational order.16 Partly for this reason, orientational phase transitions are usually identified by change of indirect parameters, such as one or more order parameters, correlation length, and radial distribution function, rather than purely thermodynamic criteriaequality of the chemical potential and pressure in the coexisting phases at the transition point. The reason is uncertainty in the chemical potential of crystalline phase evaluated with molecular simulations. Recently, it was shown that the chemical potential in the crystalline phase as a part of a two-phase system in an elongated cell can be accurately determined using the Monte Carlo method.15,48−50 This approach allows us to calculate the tangential pressure in the adsorbed layer as a function of the chemical potential in the coexisting phases at a specified temperature. The intersection point of these curves corresponds to the phase transition. To provide a rigorous thermodynamic description of the phases and phase transitions, the tangential pressure dependences on the chemical potential obtained with molecular simulations for the particular phases can be fitted by polynomials satisfying the Gibbs−Duhem equation. A more reliable way to determine the phase transition point is to generalize the pVT data in the entire region of existence of the corresponding phases in the form of regression equations for density, pressure, and free energy taking into account the Maxwell relation.50,51 It is worth to note that this approach was used to describe phases and phase transitions only in the simplest systems, namely, the Lennard-Jones crystal and fluid. In this paper, we extend it to systems that are more complicated. The main aim of the present paper is to reconsider the orientational phase transition from a rigorous thermodynamic viewpoint basing on molecular simulation of two crystalline nitrogen contact layers on graphite, which are not commensurate with the graphite lattice. These are the long-range orientationally ordered herringbone structure and the molec-
2. METHODOLOGY 2.1. Orientational Long-to-Short Order Transition. In this study we determined dependences of the tangential pressure and the internal energy on the surface density and temperature for three phases with the standard Monte Carlo method. The first phase is crystalline with hexagonal ordering of positions of N2 centers and the herringbone orientational ordering. The second phase is also crystalline, but without the long-range orientational ordering. The third phase is translationally and orientationally disordered. Simulations were performed for several series along specified values of the density. For each phase we approximated the array of values of the internal energy, u, and the tangential pressure, pT, as functions of the 2D density, ρ, and temperature, T, by regression equations as follows: n1
u = cT * +
n2
∑ ∑ (1 − j)bjk(T *) j ρk (1)
j=0 k=0 n1
pT = ρT * +
n2
∑ ∑ kbjk(T *) j ρk+ 1 j=0 k=1
(2)
Here T* = RgT (kJ/mol). The tangential pressure is presented in units of mN/m, which is usually used for the surface tension. The internal energy and density are expressed in μmol/m2 and kJ/mol, respectively. Equations 1 and 2 are derived from the Helmholtz free energy written as 2898
DOI: 10.1021/acs.jpcc.7b12300 J. Phys. Chem. C 2018, 122, 2897−2908
Article
The Journal of Physical Chemistry C n1
f = [ln(ρ) − 1]T * − cT * ln(T *) +
n2
μL − μS = T * ln(ρL /ρS ) − (c L − cS)T * ln(T *)
∑ ∑ bjk(T *) j ρk
n1
j=0 k=0
+ (b10 − c10)T * +
(3)
j=0 j≠1
Therefore, the pressure and internal energy satisfy the Maxwell relation ρ2
∑ (bj0 − cj0)(T *) j
n1
+
n2
∑ ∑ (k + 1)bjk(T *) j ρLk j=0 k=1 n1 n2
∂p /T ∂u =0 + T2 T ∂ρ ∂T
(4)
−
∑ ∑ (k + 1)cjk(T *) j ρSk
=0 (7)
j=0 k=1
For the same notations the chemical potential can be then written as n1
μ = T * ln(ρ) − cT * ln(T *) +
n1
pT,L − pT,S = (ρL − ρS )T * +
n2
j=1 k=1
∑ ∑ (k + 1)bjk(T *) j ρk
n1
(5)
j=0 k=0
−
n1
n1
j=0 j≠1
∑ ∑ kcjk(T *) j ρSk+ 1 = 0 (8)
Solving of the above equations gives the densities ρL and ρS of coexisting phases. Once the densities have been determined, the corresponding values of the internal energy and the tangential pressures can be easily calculated with eqs 1 and 2. It should be noted that the herringbone structure substantially distorts the symmetry of the 2D crystal. In the case of N2 molecular layer commensurate with the graphite lattice components of the tangential pressure are different relative to the herringbone glide line.15 In this study we analyze properties of incommensurate 2D layer which requires the components of the tangential pressure to be equal in all directions. To ensure this equality, the ratio Lx/Ly of the simulation box on the equilibration stage was corrected in each 106 Monte Carlo steps keeping the same area of the box LxLy. 2.2. Melting Transition. The method described in the previous section is not suitable for analysis of melting transition because the density of translationally disordered phase is markedly less than that of the 2D crystal when the two coexisting phases are at the same pressure and chemical potential. Keeping the same density prevents melting up to very high temperatures. Therefore, the chemical potential should be determined for each phase explicitly. To this end we modeled the gas−solid and gas−2D liquid in an elongated cell using the kinetic Monte Carlo method. The dense phase was kept at the center of the cell. The existence of the highly diluted gas-like phase at the cell ends allowed for the evaluation of the chemical potential of the whole gas−solid (or liquid) system with a high accuracy. Since at equilibrium the chemical potential is the same in each point of the system, the chemical potential was evaluated for the dense phase at the specified temperature and the tangential pressure calculated by the common virial route.53 Other details of this technique can be found elsewhere.15,50,54 A further scenario was as follows. For the disordered uniformly distributed phase over the simulation cell we simulated the internal energy and the tangential pressure as functions of the temperature and density with the standard Metropolis Monte Carlo method. Then we fitted those dependencies using eqs 1 and 2 to determine coefficients djk (except d10), similarly to the case described in the previous section. The constants of integration d10 and c10 were determined by a comparison of eq 5 for the chemical potential with its values obtained directly with the technique based on simulation of the gas−2D solid/liquid coexistence in the elongated cell with the kinetic Monte Carlo scheme.15,54 Once all parameters of the regression equations had been determined,
fL − fS = − (c L − cS)T * ln(T *) + (b10 − c10)T * n2
∑ (bj0 − cj0)(T *) j + ∑ ∑ (bjk − cjk)(T *) j ρk
n2
j=1 k=1
The parameter b10 is not involved in eqs 1 and 2 and therefore cannot be evaluated by simultaneous fitting u = u(T*,ρ) and pT = pT(T*,ρ). Being the constant of integration, this parameter can be determined either by a thermodynamic integration technique or using numerical experiments where the chemical potential is evaluated directly. In this study we used an alternative technique to analyze the long-to-short ordering transition. Since this is presumably a weakly first-order transition, the change of the internal energy and density is relatively small. For this reason, the transition could occur during gradual raising the temperature at a specified density at a condition very close to equilibrium. In this case the phase coexistence corresponds to equality of the specific Helmholtz free energy, f, of the two phases rather than their chemical potentials, i.e.
+
n2
∑ ∑ kbjk(T *) j ρLk+ 1
=0
j=0 k=1
(6)
Here subscripts “L” and “S” denote the long-range and shortrange order, respectively, with parameters bjk and cjk being specified for the two phases. The algorithm of evaluation of the phase transition parameters can be specified as follows. (1) Determine the coefficients cL, cS, bjk, and cjk (except b10 and c10) for the phases with long- and short-range orientational ordering, respectively, by fitting dependences of the internal energy and the tangential pressure on the temperature simulated for several series at specified densities. (2) At each density determine the temperature at which the long-range ordering disappears. This point can be controlled by merging of x- and y-components of the tangential pressure contributed by quadrupole−quadrupole interactions into a single line. (3) Using eq 6 fit the array of (T*, ρ) points simulated for the phase transition by the least-squares method to determine the difference (b10 − c10). (4) Having determined the difference (b10 − c10), the isobaric phase transition parameters at a specified temperature can be then evaluated at the condition of equality of the chemical potential and the pressure in coexisting phases, i.e. 2899
DOI: 10.1021/acs.jpcc.7b12300 J. Phys. Chem. C 2018, 122, 2897−2908
Article
The Journal of Physical Chemistry C the parameters of the melting transition including the temperature and the density drop were calculated at the tangential pressure corresponding to the 2D liquid−solid equilibrium which is always close to zero. In this study we first consider the case of nitrogen molecular layer adsorbed on graphite accounting for the surface mediation. However, we did not account for the effect of surface corrugation of graphite to elucidate intrinsic properties of the N2 monolayer. The case when the N2 contact layer is commensurate with the graphite lattice was comprehensively considered previously.15 Additionally, we have analyzed idealized 2D system without explicit accounting for the interaction of the nitrogen molecular layer with a hypothetical substrate having flat smooth surface. Thus, we assumed that the substrate just stabilizes the monolayer in the 2D form. The displacement of molecules and its rotation is allowed only along the xy-plane. This simplification allows for evaluation of main features of the 2D molecular layer composed of linear quadrupolar molecules including orientational transitions. 2.3. Model and Simulation Details. The uniform twodimensional liquid and solid were simulated with the standard Metropolis scheme in a canonical ensemble. The cutoff distance, rc, was taken as 5σ. In order to account for the tail correction, in each 106 MC steps the potential and pressure components related to each molecule were calculated with the cutoff distance of 10σ. Therefore, the tail values of the potential and pressure of the molecule were attributed to a “frozen” shell between spheres of the radius 5σ and 10σ, which refreshed in each 106 MC steps. The solid state was modeled as the crystal of 400 molecules arranged in hexagonal lattice. The density of the molecular layer was specified at a number of series with the dimensions Lx and Ly of the simulation cell being adjusted on the equilibration stage to ensure the equality of the x- and ycomponents of the tangential pressure at a given surface density. Initially the orientation of molecules was set as the herringbone arrangement with the glide line directed along the x-axis. The temperature was gradually increased at constant density up to the disappearance of the long-range orientational ordering, which was indicated by coinciding the x- and ycomponents of the electrostatic contribution to the total tangential pressure. In the proximity of the orientational phase transition the temperature was incremented by 0.1 K. The number of MC steps was 3 × 108 and 1 × 108 on the equilibration and sampling stage, respectively. Modeling of the crystalline phase with the short-range orientation order and the translationally disordered phase was performed in the simulation cell with the ratio Lx/Ly = 2/31/2. We used the two-center Lennard-Jones model for the nitrogen molecule with the distance between the atoms of 0.1098 nm. The LJ parameters for the N−N interaction were taken as 36.4 K and 0.3318 nm for the potential well depth, εNN/kB, and the collision diameter, σNN, respectively.55 The interaction potential of the molecules i and j is given by uijLJ =
∑ φij(a , b) a,b
⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σ ⎢⎜ σNN ⎟ ⎟⎥ = 4εNN ∑ ⎢⎜ (a , b) ⎟ − ⎜⎜ (NN a , b) ⎟ ⎥ r ⎝ rij ⎠ ⎥⎦ a,b ⎢ ⎣⎝ ij ⎠
pTLJ, x
pTLJ, y
(a , b)
1 = ρkBT − Sxy
∑ xij ∑ i
