Liquid-Hextic-Solid Phase Transition of a Hard-Core Lattice Gas with Third Neighbor Exclusion

TL;DR

This paper investigates the phase transitions in a two-dimensional hard-core lattice gas with third neighbor exclusion, revealing a two-step liquid-hexatic-solid transition through a novel kinetic approach and visualization techniques.

Contribution

It introduces an alternative kinetic method based on the Gibbs adsorption isotherm and random sequential adsorption with surface diffusion to analyze phase behavior.

Findings

First-order phase transition confirmed at high surface coverage.

Two-step liquid-hexatic-solid transition observed.

Visualization of local bond orientation supports phase identification.

Abstract

The determination of phase behavior and, in particular, the nature of phase transitions in two-dimensional systems is often clouded by finite size effects and by access to the appropriate thermodynamic regime. We address these issues using an alternative route to deriving the equation of state of a two-dimensional hard-core particle system, based on kinetic arguments and the Gibbs adsorption isotherm, by use of the random sequential adsorption with surface diffusion (RSAD) model. Insight into coexistence regions and phase transitions is obtained through direct visualization of the system at any fractional surface coverage via local bond orientation order. The analysis of the bond orientation correlation function for each individual configuration confirms that first-order phase transition occurs in a two-step liquid-hexatic-solid transition at high surface coverage.

Equations

$$d\Pi =kT\frac{\Theta }{A_a}d\ln C$$

$$\underset{0}{\overset{\Theta }{\int }}\frac{\Theta }{C}\frac{\partial C}{\partial \Theta }d\Theta =\frac{A_a}{kT}\Pi$$

$$K_{a}C(1-\beta (\Theta ))=K_d\Theta$$

$$\underset{0}{\overset{\Theta }{\int }}(1-\beta (\Theta ))\frac{\partial }{\partial \Theta }\left[\frac{\Theta }{1-\beta (\Theta )}\right]d\Theta =\frac{A_a}{kT}\Pi$$

$$\frac{\partial N}{\partial n}=1-\beta (\Theta )$$

$$\Psi (r_j)=\frac{1}{N_k}\sum _{k=1}^{N_k}{e}^{i6{\theta }_{jk}}$$

$$g_6(r)=\frac{⟨\sum _{k\ne j}^N\Psi {(r_j)}^{\ast }\Psi (r_k)\delta (r-\left|r_j-r_k\right|)⟩}{⟨\sum _{k\ne j}^N\delta (r-\left|r_j-r_k\right|)⟩}$$

Full text

I Introduction

The phase behavior of the two dimensional systems is a key aspect of many current research areas such as emulsion stability due to particle adsorption at the interface Pauchard et al. (2014); Liu et al. (2017), particle self-assembly into clusters Hou et al. (2019); Fortuna et al. (2010); Weber et al. (2008); Gorbunov et al. (2013), chemisorption on metal surface Bak et al. (1985); Taylor et al. (1985) and melting at an interface Bernard and Krauth (2011); Engel et al. (2013). Numerous equations of state ($EOS$) such as the Langmuir Langmuir (1918) and Volmer models have been introduced over the years to describe the adsorption behavior of such systems Rane et al. (2015, 2013). The Langmuir model is based on localized adsorption where the adsorbate molecules are smaller than the adsorption sitesFoo and Hameed (2010), but in many practical cases the adsorbate is actually larger than the adsorption site. Likewise, the Volmer model, which assumes fully delocalized adsorption (adsorbates much bigger than the adsorption sites), is equally inappropriate except, perhaps, for nanoparticles. Rather than directly prescribing the $EOS$, an alternative is the random sequential adsorption model ($RSA$), which describes the dynamics of the adsorption process by allowing objects to adsorb sequentially onto the open sites of a one or two-dimensional lattice Wang et al. (1993a). The model has been extensively used in the literature, for example by Manzi et al. Manzi et al. (2019) who discuss the adsorption of human serum albumin on the nano-structure of a black silicon surface using $RSA$. Further applications of this model include chemisorption, deposition, layered growth, vibrated granular material and the car-parking problem Lončarević et al. (2009, 2007); Talbot et al. (2000); Evans (1993); Krapivsky and Ben-Naim (1994); Wang et al. (1993b); Šćepanović et al. (2016); Tarjus and Viot (2004); Meineke and Gezelter (2001); Schaaf and Talbot (1989).

More specifically, in the $RSA$ model molecules or particles are progressively added at random to an initially empty surface with the only restriction that overlap is not allowed, an assumption based physically on short-range electrostatic repulsion. As the coverage increases, the free area left for further adsorption decreases, not only because of the sites occupied by previously adsorbed molecules, but also because vacancies can be too small to allow adsorption without overlap. Without desorption or surface diffusion, adsorption kinetic rapidly slows down and coverage only asymptotically approaches the jamming limit, equivalent to random maximum packing, if the substrate is not pre-patterned. However, it has been experimentally observed Ramsden (1992) that the relaxation time scale of adsorbed particles, due to their rearrangement on the surface, can be comparable to the deposition time scale. The final configuration is comparable with a dense-packed ordered system. Furthermore, none of the above models is able to explain the ordered layering observed in adsorption of certain materials on the surface Bak et al. (1985); Taylor et al. (1985); Ramsden (1992).

A different approach is taken in the lattice gas model, which uses statistical mechanical reasoning to describe adsorbate configurations on a surface, where an adsorbate molecule may occupy one or a few adsorption sites. Applications of this model include a recent study of photo-excited Rydberg gases by Ji et al. Ji et al. (2011), where an order-disorder phase transition corresponds to the phase transition on a square lattice with first neighbor exclusion, a study of self-assembly of isophthalic acid on graphite by Lackinger et al. Lackinger et al. (2004), adsorption of selenium on Nickel surface by Bak et al. Bak et al. (1985) and chemisorption of oxygen on palladium by Zhang et al. Zhang et al. (2007). Although the lattice gas model has been studied extensively in the literature, only the single case of a triangular lattice with first neighbor exclusion was solved exactly, by Baxter Baxter (2016). For all other variants, a number of lattice gas methods have been developed over years based on various approximation methods: the matrix method of Kramer and Wannier Bellemans and Nigam (1967); Orban and Bellemans (1968); Ree and Chesnut (1966); Runnels et al. (1967); Runnels and Combs (1966); Feng et al. (2011); Rotman and Eisenberg (2009), the density (or activity) series expansion method Bellemans and Nigam (1967); Orban and Bellemans (1968); Gaunt and Fisher (1965); Gaunt (1967); Eisenberg and Baram (2005); Ushcats (2015); Ushcats et al. (2016) , the generalized Bethe method Cowley (1979); Burley (1960); Hansen-Goos and Weigt (2005), Monte Carlo simulation Feng et al. (2011); Binder and Landau (1980); Chesnut (1971); Liu and Evans (2000); Fernandes et al. (2007a); Nath and Rajesh (2014), the Rushbrooke and Scoins method Rushbrooke and Scoins (1955) and fundamental measure theory Lafuente and Cuesta (2003a). Despite all of this effort, the lattice gas model is not able to describe the adsorption isotherm of the system, and instead focuses on the equation of state and the nature of phase transition.

In order to combine the advantages of both the RSA and lattice gas models, we previously developed an alternative route to deriving the equation of state of a two-dimensional hard-core particle with first neighbor exclusion, based on kinetic arguments and the Gibbs adsorption isotherm: the $RSAD$ model Darjani et al. (2017). Here one considers a two-dimensional lattice gas in equilibrium with a three-dimensional solution of adsorbate molecules, where the equality of chemical potential throughout the system leads to:

[TABLE]

where $A_{a}$ is the interfacial area covered by a single adsorbate molecule, $\theta$ is the fractional surface coverage and $C$ is the concentration of the (three-dimensional) solution. Integrating the above equation gives:

[TABLE]

From which we see that knowledge of the adsorption isotherm, the relationship between $C(\Theta)$, bulk concentration, and fractional coverage, enables one to calculate the equation of state $\Pi(\Theta)$.

The adsorption isotherm, in turn, can be obtained through kinetic arguments. At equilibrium the rates of adsorption and desorption of molecules are equal:

[TABLE]

where $K_{a}$ and $K_{d}$ are the adsorption and desorption rate constants, respectively, and $\beta(\Theta)$ is the “blocking function ”, the fraction of the surface area which is excluded from further adsorption by already adsorbed molecules. Solving for $C$ and inserting the result into the integral version of the Gibbs adsorption isotherm yields:

[TABLE]

Thus, the blocking function is the only information needed to calculate the equation of state, and we have shown previously Darjani et al. (2017) that for lattice gases the blocking function can be easily extracted from $RSAD$ model simulations.

In the $RSAD$ model, where surface diffusion is introduced in parallel with adsorption, vacancies large enough to adsorb a further particle are both created and destroyed. When diffusion is sufficiently rapid, the size distribution of vacancies no longer depends on the history of adsorption (the positions where the adsorbates first arrived on the substrate) but only on the fractional surface coverage. One of the advantages of using the $RSAD$ model is in locating the equilibrium state, which assures us that enough thermalization is present to reach the equilibrium state. Note that in this model the potential energy is effectively infinite for particle overlap, due to the repulsive interaction, which restricts the occupancy of neighbors, and is zero otherwise. The system can therefore be considered as athermal Rotman and Eisenberg (2009); Fernandes et al. (2007a); Darjani et al. (2017). Our results show that the $RSAD$ model can be used as an equilibrium model and our equation of state, the nature of our phase transition and the phase transition coverage are in excellent agreement with the only model with an exact solution in the literature Baxter (2016). From the definition of the adsorption rate, used above to define adsorption equilibrium, the blocking function can be extracted from the numerical simulations through the derivative of surface coverage with respect to the number of attempts:

[TABLE]

Here $N$ is the number of adsorbed molecules, $n$ is the number of attempts and t, defined by $nA_{a}/A$=$K_{a}C/t$, is an adimensional time defined via the adsorption rate. The latter definition will be used in practice by equating the blocking function to the rebuttal rate of adsorption attempts. Ushcats et al. Ushcats et al. (2016) used an alternative method based on the power of activity at low and high densities, which shows the importance of accounting for the holes in deriving the equation of state of lattice gases. Based on their method Ushcats et al. (2016); Ushcats (2015), hole-particle symmetry, the total interaction energy is directly related to the interaction of holes at any specific configuration.

In this paper we study the phase behavior of hard-core molecules with third neighbor exclusion on a triangular lattice and compare our results with those of Orban et al. Orban and Bellemans (1968) who studied this model previously. Hard-core molecules with extended exclusion ranges are studied extensively in the literature but mainly on a square lattice Feng et al. (2011); Rotman and Eisenberg (2009); Fernandes et al. (2007b, a); Lafuente and Cuesta (2003b); Ramola and Dhar (2012); Mandal et al. (2017); Nath and Rajesh (2014, 2016). In some experiments it is observed that lateral interactions of adsorbed particles on solid surfaces (chemisorption) follow the extended exclusion range, which is important in surface science as ordering affects the surface functionality Bak et al. (1985); Taylor et al. (1985); Zhang et al. (2007); Liu and Evans (2004). Increasing the exclusion range could also correspond to a smaller lattice site which becomes equivalent to the continuum limit when the exclusion range is significantly large Nath and Rajesh (2014); Fernandes et al. (2007a); Orban and Bellemans (1968) .

Simulation data related to the triangular lattice with third neighbor exclusion is given in the next section, where a first order liquid-hexatic-solid phase transition is obtained at high surface coverage. In the liquid state, positional and orientational correlations of particles decay exponentially, while in solid state they have a quasi-long-range orientational order. Besides these two phases there is an intermediate hexatic phase, where particles have quasi-long range orientational order and exponential positional order. We quantify the ordering using a bond orientation correlation function, ${g_{6}}(r)$, calculated from the local bond orientation order, $\Psi({r})$. In a dense system, most of the particles are surrounded by six particles, and the local bond orientation is represented by the sixfold orientation:

[TABLE]

Where $k$ is the number of the nearest neighbors of particle $j$, $\theta_{jk}$ is the angle between the line joining the centers of mass of particles $j$ and $k$ and a reference axis. The bond orientation correlation function is then defined as an average of the local bond orientation order:

[TABLE]

A power-law decay of ${g_{6}}(r)$ means that there is a quasi-long-range orientation correlation. For the system in an hexatic phase ${g_{6}}(r)\propto{r^{-\eta}}$ where $0<\eta<0.25$.

Liquid-solid transitions have been reported extensively in the literature for two-dimensional systems in lattice gas models with extended exclusion ranges. However here our simulation results reveal a liquid-hexatic-solid phase transitions for the triangular lattice cover seven sites, which have the same nature of the phase transition as the melting transition of hard-disk molecules. Although melting is studied extensively in the literature Bernard and Krauth (2011); Engel et al. (2013), there is still controversy about the nature of its phase transition that can arise from finite size effect or inefficiency of the system to equilibrate the system. In two-dimensional systems, there are three different scenarios for the melting transition of hard disks in the literature:

• •

Kosterlitz, Thouless, Halperin, Nelson and Young ($KTHNY$) scenario Young (1979); Nelson and Halperin (1979); Kosterlitz and Thouless (1973): a two-step continuous phase transition, first from liquid to hexatic phase and then from hexatic to solid phase Mak (2006); Peng et al. (2010),

• •

Two step transition: first order phase transition from liquid to hexatic phase and then continuous phase transition from hexatic to solid phase Bernard and Krauth (2011); Engel et al. (2013),

• •

First order phase transition from liquid to solid phase Alder and Wainwright (1962).

We will show in the next section how finite size effects and access to the thermodynamic regimes can bring uncertainty regarding phase behavior of the system; mainly the nature of phase transition where these issues are extensively reported in the literature about the nature of melting transition Bernard and Krauth (2011); Wierschem and Manousakis (2011).

II Simulation details

The adsorption of hard-core molecules with a third neighbor exclusion range on the triangular lattice involves the adsorption of molecules covering 7 adsorption sites in the manner as represented by a red circle in Figure. 1. As it is clear from the hexagon drawn around the adsorbate (red circle) in Figure. 1(a)-(b), this model could have two different orientations at high coverage in comparison to other adsorbates such as hard-core molecule with first neighbor exclusion Darjani et al. (2017). Here, we employ two complementary methods as in our previous work Darjani et al. (2017): an “adsorption method”, which begins from an empty lattice, and a “desorption method”, which begins with a full lattice and progressively decreases coverage. The results are expected to bracket the correct equilibrium equation of state.

In the adsorption method, molecules or particles are progressively added to an initially empty d$\times$d lattice surface where a periodic boundary condition is used to ameliorate finite size effects. The only restriction is that overlap is not allowed; an assumption based physically on short-range electrostatic repulsion. For each adsorption attempt, a random position $(x,y)$ is selected representing the center of mass of the particle. If the selected site and its neighbors are empty, adsorption is accepted. Otherwise it is rejected. Diffusion, the simultaneous movement of particles, is introduced sequentially with a predefined ratio $D$ between the number of diffusion attempt and the adsorption attempt: For $D=3$ each adsorption is followed by 3 diffusion attempts, etc. For each diffusion attempt, a previously adsorbed particle and a direction for the displacement of the particle are selected randomly; an arrow in Figure. 1 illustrates the possible direction. If moving the center of the mass of the particle to the next node along this direction does not infringe the non-overlap condition, diffusion is accepted. Otherwise it is rejected. In the RSAD model, when diffusion is fast enough, the surface layer is at internal equilibrium (even during transient adsorption) and the blocking function can be considered as a state function.

For the desorption method the lattice is initially full. For each simulation step, two particles are randomly selected and removed. Then one adsorption attempt and $D$ diffusion attempts are performed following the same procedure as for the adsorption method. The choice of the sequence ($2$ desorption events followed by $1$ adsorption) is arbitrary but answers the need at each time step to decrease coverage and add at least one particle to calculate the blocking function.

For both adsorption and desorption methods, the blocking function is extracted from the success rate of adsorption attempts. $500$ runs are performed, and an ensemble average is used to reduce the noise arising from the numerical calculation of the derivative of the coverage. The blocking function is fitted with a polynomial function before it is used to generate the adsorption isotherm. The latter is inserted into the Gibbs adsorption isotherm equation to obtain the equation of state.

In this work we also used a relaxation method in order to track the structure of particles by time and proof our hypothesis about expecting the correct equilibrium equation of state starting from either the adsorption or desorption method. In this method, either the adsorption or desorption method can be used to reach a specific coverage. In the second step, one adsorption attempt and $D$ diffusion attempts are performed following the same procedure as for the adsorption method. In order to keep the coverage constant, if the adsorption attempt is successful, one particle is randomly selected and removed. $1500$ runs are performed to extract the success rate of adsorption attempts or, in other words, the blocking function.

III Results

Accessing the thermodynamic regime is an initial step toward studying the phase behavior of the system Bernard and Krauth (2011); Engel et al. (2013) where some algorithms are insufficient to equilibrate the system in order to study the phase transition of the system in a reliable manner Mak (2006). The effect of surface diffusion for a lattice size $d=105$ is studied in Figure. 2. Initially, when the system is diluted, all of the curves regardless of their methods or their surface diffusion overlap in the low surface coverage as presented in Figure. 2(a). The fluctuation in the phase transition region is much larger than the pure liquid and solid phase, where the difference is maximized in the middle of phase transition as illustrated in Figure2(b) Alder and Wainwright (1962). As presented in Figure2 (b), at high surface coverage probability of success of adsorbing a new particle initially decreases due to the caging effect. However by increasing the ordering of particles, this caging effect will be diminished in order to maximize the available surface for accepting the new incoming particles. The system reaches the equilibrium state when two curves overlap in the whole range of fractional surface coverage, so accessing the thermodynamic regime will be apparent.

Figure. 3 represents the relaxation of blocking function of hard-core molecules with third neighbor exclusion on a triangular lattice for $D=0.01$, $d=196$ and $\theta=0.915$. This figure confirms that the equilibrium state is a function of a blocking function and does not depend on initial configuration. The blocking function of the desorption method increases in time until it reaches the equilibrium value, whereas the adsorption method shows the opposite trend. Simulation data shows that the desorption method reaches the equilibrium value faster than the adsorption method.

Finite size effects are another important key issue in finding the equation of state in lattice simulations. Figure. 4 presents the variation in adsorption rate for triangular lattices covering $7$ sites of sizes $105$ to $266$. For surface diffusion of $D=1$, the same blocking function is obtained from both adsorption and desorption methods at low surface coverage as presented in Figure. 4(a). Initially, in the vicinity of the phase transition, by increasing the lattice size dimension, the adsorption rate increases for both the adsorption and desorption method as illustrated in Figure. 4(b). Exactly before all of the curves overlap around full coverage, the adsorption rate decreases by increasing the lattice dimension for both methods. The results in Figure. 4(b) show that the adsorption rate is significantly different for $d=105$ in comparison to $d=196$ and $d=266$. As a result $d=105$ is not a reliable lattice for finding the equation of state. Later on we will discuss this case and compare it with $d=196$ in Figure. 5(b) and Figure. 6 only to show the importance of using a large enough system. One of the advantages of using the $RSAD$ method is we know how big our system should be to ensure that the results are at the same time accurate and computationally less expensive.

Although an exact solution of this model does not exist in the literature, meaning there is no certain agreement upon the equation of state and its phase transition, we can compare our results with the analytic calculation of Orban and Bellmans Orban and Bellemans (1968) for $d=196$ and $D=80$ in Figure. 5(a). This paper used two different methods, the matrix method, based on the sequence of exact solutions for lattices of infinite length and increasing finite width, and the series expansion method, based on knowing the final structure at close packing and constructing the density and activity series to find the surface pressure. A first order transition was found for both the matrix and the series expansion methods, where the transition occurs at surface pressure of $4.5-5$ for surface coverage of ${\rm{0.81}}\pm{\rm{3}}$ and ${\rm{0.98}}\pm{\rm{1}}$ for fluid and solid regimes respectively. This system is also studied by Chestnut Orban and Bellemans (1968) using a Monte Carlo technique, but the study was limited to surface coverages lower than $0.75$, where no phase transition was found.

As illustrated in Figure. 5(a), at low surface coverage there is no difference between the reported equations of state. However, in the lower part of the phase transition zone, our equation of states shows lower surface pressure than Orban and Bellman’s series expansion method (low and high density) but higher surface pressure than their matrix method. On the contrary, in the upper part of the phase transition region, the matrix method shows higher surface pressure than what we found, but the series expansion method overlaps with both of our methods. The thermodynamically-stable surface pressure loop is observed in our simulation results as is reported extensively in the literature due to finite size effects Bernard and Krauth (2011); Engel et al. (2013); Mayer and Wood (1965); Schrader et al. (2009). Surface pressure can create a thermodynamically stable loop at equilibrium for a finite size system, but this loop will disappear at infinite size and the coexistence zone will be flat. Creation of the flat pressure is visible in Figure. 5(b) for $d=196$ around a surface pressure of ${\rm{4.5}}\pm{\rm{0.05}}$. Moreover, based on the lower part of the phase transition in Figure. 4(b), the adsorption rate tends to increase by increasing the lattice dimension, which means less surface pressure will be expected for a larger lattice size. Conversely, the adsorption rate tends to decrease for both the adsorption and desorption methods in the upper part of the phase transition, which means higher surface pressure will be expected for a larger lattice size. This behavior indicates the tendency of the system towards flatness at an infinite size. From the equality shown in the hatched area of Figure. 5(b), the horizontal surface pressure is obtained from the Maxwell construction, where overlapping of this construction line with the flat region of the equation of state confirms the tendency of the system to be flat at infinite size. Figure. 5(b) shows that by increasing the lattice dimension, surface pressure decreases in the lower part of the phase transition and increases in the upper part of phase transition. Bernard et al Bernard and Krauth (2011); Engel et al. (2013) reported the same trend for the equation of state of the melting transition of hard disks by increasing the number of particles where they reported a two step transition; first order liquid-hexatic and continuous hexatic-solid transition for their system.

The phase transition zone is studied in Figure. 6 based on the derivative of surface pressure with respect to surface coverage in the adsorption method. One can see from Figure. 6(b) that for an insufficiently large system, the first order liquid-solid phase transition is obtained for $d=105$ and $D=80$; however, a different phase transition is obtained for $d=196$ and $D=80$ in Figure. 6(a), which indicates the importance of using a sufficiently large system. Phase transition peaks obtained from Figure. 6(a) for $d=196$, $D=80$ are analyzed through a bond orientation correlation function, g6(r), in Figure 7 based on the relaxation method for each individual configuration. For additional insight into the phase transition region, each of these configurations is visualized via the local bond orientation order function $\Psi({r})$, based on the relaxation method, in Figure. 8. The first peak in the phase transition curve corresponds to surface coverage of $0.826$. Since we find the same exponential decay of ${g_{6}}(r)$ before and after relaxation at surface coverage of $0.75$ in Figure. 7(a) suggests that the system is in a pure liquid regime below a surface coverage of $0.826$. Furthermore, the homogeneous distribution of local bond orientation order over the surface in Figure. 8(a) confirms the presence of a liquid phase below surface coverage of $0.826$. In Figure. 6 we see that lattice size does not have any impact on the liquid phase since the same behavior from liquid regime is obtained for both $d=105$ and $d=196$ through fractional surface coverages up to $0.826$.

The transition between positive and negative slopes occurs at surface coverage of $0.864$, located between the first and second peaks, so we analyze surface coverages of $0.85$ and $0.869$ to see the behavior of the system below and above the transition value. For surface coverage of $0.85$, ${g_{6}}(r)$ decays exponentially at the same rate before and after relaxation in Figure. 7(b); more peaks were observed after relaxation, which is an indication of the creation of a more ordered phase in the system. Creation of a more ordered phase is also confirmed by the local bond orientation order seen in Figure. 8(b), indicated by more red spread through the surface. For surface coverage of $0.869$, ${g_{6}}(r)$ still decays exponentially, but at different rates before and after relaxation as illustrated in Figure. 7 (c-1) and (c-2), respectively. The local bond orientation order in Figure. 8(c) shows that ordered particles prefer to stick to each other. The third peak in Figure. 6 corresponds to surface coverage of $0.915$. At this coverage the system decays exponentially before relaxation as presented in Figure. 7(d-1) and Figure. 8(d) shows that the cluster of ordered phase spreads through the surface and indicates a glassy state. After relaxation, ${g_{6}}(r)$ decays with a power law with $\theta=0.25$ as illustrated in Figure. 7(d-2) , which is an indication of a hexatic phase. The local bond orientation order shown in Figure. 8(d) reveals that at this coverage all the particles tend to stick together and they are between the mobile particles, as was also reported in reference Hou et al. (2019). The tendency of ordered particles to stick together causes the probability of success to increase in Figure. 2(b) due to maximizing the available free surface for accepting the new incoming particles. By increasing the surface coverage from 0.915 to 0.963, $\eta$ decreases from $0.25$ to $0.08$ after relaxation as illustrated in Figure. 7(e-2), which indicates that the fourth peak in Figure. 6 corresponds to a first order transition from hexatic phase to solid phase in equilibrium state.

To further illustrate the relaxation dynamics toward an equilibrium state, the local bond orientation order parameter of particles at surface coverage of $0.915$ was tracked over time for two different initial configurations obtained from adsorption and desorption methods at very low surface diffusion ($D=0.01$, $d=196$). Figure. 9(a) indicates that the density should be increased very slowly when we start from an empty lattice, or else the system will be locked into a glassy configuration Alder and Wainwright (1962). A glassy state is reported in reference Hou et al. (2019) during rapid compression of a system composed of spherical particles. For the desorption method, Figure. 9(b), the system is initially ordered, so the density should be decreased very slowly or else the system will be metastable Alder and Wainwright (1962). As time passes in the equilibrium state, phase separation occurs, where ordered particles tends to cluster and create more space for further adsorption.

IV Conclusion

In this paper we studied the phase behavior of a hard-core molecules with third neighbor exclusions in a triangular lattice with $RSAD$ simulations, in order to derive the equation of state of a two-dimensional hard-core particle based on kinetic arguments and the Gibbs adsorption isotherm. We compared our results with Orban and Bellman Orban and Bellemans (1968), who used a matrix method and a series expansion methods in low and high densities, and found only partial agreement. Our results show that the system is in a pure liquid regime below surface coverage of $0.826$. Increasing the surface coverage will create a more ordered phase where a first order liquid-hexatic phase transition occurs between surface coverage of $0.877$ and $0.915$. At surface coverage of $0.915$, ${g_{6}}(r)$ decays algebraically after relaxation with a power-law exponent $\eta=0.25$, which is an indication of a hexatic phase. Our simulation results reveal that as the surface coverage increases, after relaxation, the surface particles tend to form tightly-packed clusters in like-oriented domains, while the remaining mobile particles have more random orientations. By increasing the surface coverage above $0.915$, $\eta$ decreases from $0.25$ toward zero after relaxation, where the system undergoes a first order transition from a hexatic phase to a solid phase.

One of the advantages of using the $RSAD$ model is being able to locate the equilibrium state, which assures us that adequate thermalization and finite size are being used to reach the equilibrium state. Moreover, subtle details of the clustering structure, through direct visualization of the system using the relaxation method at any fractional surface coverage, provides insight regarding coexistence regions and phase transitions.

Acknowledgments

This research was supported by the National Science Foundation under grant $No.1743794$, PIRE: Investigation of Multi-Scale, Multi-Phase Phenomena in Complex Fluids for the Energy Industries.