Liquid-Hextic-Solid Phase Transition of a Hard-Core Lattice Gas with Third Neighbor Exclusion
Shaghayegh Darjani, Joel Koplik, Sanjoy Banerjee, Vincent Pauchard
LLiquid-Hexatic-Solid Phase Transition of a Hard-Core Lattice Gas with ThirdNeighbor Exclusion
Shaghayegh Darjani , † , Joel Koplik , ∗ Sanjoy Banerjee , and Vincent Pauchard Energy Institute and Department of Chemical Engineering,City College of New York, New York, NY 10031, USA. Benjamin Levich Institute and Department of Physics, City College of New York, NY 10031, USA. Benjamin Levich Institute and Department of Chemical Engineering,City College of the City University of New York, New York, New York 10031, USA (Dated: August 16, 2019)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 ther-modynamic regime. We address these issues using an alternative route to deriving the equation ofstate of a two-dimensional hard-core particle system, based on kinetic arguments and the Gibbsadsorption isotherm, by use of the random sequential adsorption with surface diffusion (
RSAD )model. Insight into coexistence regions and phase transitions is obtained through direct visualiza-tion of the system at any fractional surface coverage via local bond orientation order. The analysis ofthe bond orientation correlation function for each individual configuration confirms that first-orderphase transition occurs in a two-step liquid-hexatic-solid transition at high surface coverage.
I. INTRODUCTION
The phase behavior of the two dimensional systemsis a key aspect of many current research areas such asemulsion stability due to particle adsorption at the in-terface [1, 2], particle self-assembly into clusters [3–6],chemisorption on metal surface [7, 8] and melting at aninterface [9, 10]. Numerous equations of state (
EOS )such as the Langmuir [11] and Volmer models have beenintroduced over the years to describe the adsorption be-havior of such systems [12, 13]. The Langmuir modelis based on localized adsorption where the adsorbatemolecules are smaller than the adsorption sites[14], but inmany practical cases the adsorbate is actually larger thanthe adsorption site. Likewise, the Volmer model, whichassumes fully delocalized adsorption (adsorbates muchbigger than the adsorption sites), is equally inappropri-ate except, perhaps, for nanoparticles. Rather than di-rectly prescribing the
EOS , an alternative is the randomsequential adsorption model (
RSA ), which describes thedynamics of the adsorption process by allowing objectsto adsorb sequentially onto the open sites of a one or two-dimensional lattice [15]. The model has been extensivelyused in the literature, for example by Manzi et al. [16]who discuss the adsorption of human serum albumin onthe nano-structure of a black silicon surface using
RSA .Further applications of this model include chemisorption,deposition, layered growth, vibrated granular materialand the car-parking problem [17–26].More specifically, in the
RSA model molecules or par-ticles are progressively added at random to an initiallyempty surface with the only restriction that overlap isnot allowed, an assumption based physically on short-range electrostatic repulsion. As the coverage increases, ∗ [email protected] the free area left for further adsorption decreases, notonly because of the sites occupied by previously adsorbedmolecules, but also because vacancies can be too smallto allow adsorption without overlap. Without desorp-tion or surface diffusion, adsorption kinetic rapidly slowsdown and coverage only asymptotically approaches thejamming limit, equivalent to random maximum pack-ing, if the substrate is not pre-patterned. However, ithas been experimentally observed [27] that the relaxationtime scale of adsorbed particles, due to their rearrange-ment on the surface, can be comparable to the depositiontime scale. The final configuration is comparable with adense-packed ordered system. Furthermore, none of theabove models is able to explain the ordered layering ob-served in adsorption of certain materials on the surface[7, 8, 27].A different approach is taken in the lattice gas model,which uses statistical mechanical reasoning to describeadsorbate configurations on a surface, where an adsor-bate molecule may occupy one or a few adsorption sites.Applications of this model include a recent study ofphoto-excited Rydberg gases by Ji et al. [28], wherean order-disorder phase transition corresponds to thephase transition on a square lattice with first neighborexclusion, a study of self-assembly of isophthalic acid ongraphite by Lackinger et al. [29], adsorption of seleniumon Nickel surface by Bak et al. [7] and chemisorption ofoxygen on palladium by Zhang et al. [30]. Althoughthe lattice gas model has been studied extensively inthe literature, only the single case of a triangular lat-tice with first neighbor exclusion was solved exactly, byBaxter [31]. For all other variants, a number of latticegas methods have been developed over years based onvarious approximation methods: the matrix method ofKramer and Wannier [32–38], the density (or activity)series expansion method [32, 33, 39–43] , the generalizedBethe method [44–46], Monte Carlo simulation [37, 47– a r X i v : . [ c ond - m a t . s o f t ] A ug RSAD model [54]. Here one consid-ers a two-dimensional lattice gas in equilibrium with athree-dimensional solution of adsorbate molecules, wherethe equality of chemical potential throughout the systemleads to: d Π = kT Θ A a d ln C (1)where A a is the interfacial area covered by a single adsor-bate molecule, θ is the fractional surface coverage and C is the concentration of the (three-dimensional) solution.Integrating the above equation gives: Θ (cid:90) Θ C ∂C∂ Θ d Θ = A a kT Π (2)From which we see that knowledge of the adsorptionisotherm, the relationship between C (Θ), bulk concen-tration, and fractional coverage, enables one to calculatethe equation of state Π(Θ).The adsorption isotherm, in turn, can be obtainedthrough kinetic arguments. At equilibrium the rates ofadsorption and desorption of molecules are equal: K a C (1 − β (Θ)) = K d Θ (3)where K a and K d are the adsorption and desorption rateconstants, respectively, and β (Θ) is the “blocking func-tion”, the fraction of the surface area which is excludedfrom further adsorption by already adsorbed molecules.Solving for C and inserting the result into the integralversion of the Gibbs adsorption isotherm yields: Θ (cid:90) (1 − β (Θ)) ∂∂ Θ (cid:20) Θ1 − β (Θ) (cid:21) d Θ = A a kT Π (4)Thus, the blocking function is the only informationneeded to calculate the equation of state, and we haveshown previously [54] that for lattice gases the blockingfunction can be easily extracted from
RSAD model sim-ulations.In the
RSAD model, where surface diffusion is intro-duced in parallel with adsorption, vacancies large enoughto adsorb a further particle are both created and de-stroyed. When diffusion is sufficiently rapid, the sizedistribution of vacancies no longer depends on the his- tory of adsorption (the positions where the adsorbatesfirst arrived on the substrate) but only on the fractionalsurface coverage. One of the advantages of using the
RSAD model is in locating the equilibrium state, whichassures us that enough thermalization is present to reachthe equilibrium state. Note that in this model the po-tential energy is effectively infinite for particle overlap,due to the repulsive interaction, which restricts the oc-cupancy of neighbors, and is zero otherwise. The systemcan therefore be considered as athermal [38, 50, 54]. Ourresults show that the
RSAD model can be used as anequilibrium model and our equation of state, the natureof our phase transition and the phase transition coverageare in excellent agreement with the only model with anexact solution in the literature [31]. From the definitionof the adsorption rate, used above to define adsorptionequilibrium, the blocking function can be extracted fromthe numerical simulations through the derivative of sur-face coverage with respect to the number of attempts: ∂N∂n = 1 − β (Θ) (5)Here N is the number of adsorbed molecules, n is thenumber 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 equatingthe blocking function to the rebuttal rate of adsorptionattempts. Ushcats et al. [43] used an alternative methodbased on the power of activity at low and high densities,which shows the importance of accounting for the holes inderiving the equation of state of lattice gases. Based ontheir method [42, 43], hole-particle symmetry, the totalinteraction energy is directly related to the interaction ofholes at any specific configuration.In this paper we study the phase behavior of hard-coremolecules with third neighbor exclusion on a triangularlattice and compare our results with those of Orban etal. [33] who studied this model previously. Hard-coremolecules with extended exclusion ranges are studied ex-tensively in the literature but mainly on a square lat-tice [37, 38, 50, 51, 55–59]. In some experiments it isobserved that lateral interactions of adsorbed particleson solid surfaces (chemisorption) follow the extended ex-clusion range, which is important in surface science asordering affects the surface functionality [7, 8, 30, 60].Increasing the exclusion range could also correspond toa smaller lattice site which becomes equivalent to thecontinuum limit when the exclusion range is significantlylarge [33, 50, 51] .Simulation data related to the triangular lattice withthird neighbor exclusion is given in the next section,where a first order liquid-hexatic-solid phase transitionis obtained at high surface coverage. In the liquid state,positional and orientational correlations of particles de-cay exponentially, while in solid state they have a quasi-long-range orientational order. Besides these two phasesthere is an intermediate hexatic phase, where particleshave quasi-long range orientational order and exponen-tial positional order. We quantify the ordering using abond orientation correlation function, g ( r ), calculatedfrom the local bond orientation order, Ψ( r ). In a densesystem, most of the particles are surrounded by six par-ticles, and the local bond orientation is represented bythe sixfold orientation:Ψ( r j ) = 1 N k N k (cid:88) k =1 e i θ jk (6)Where k is the number of the nearest neighbors of particle j , θ jk is the angle between the line joining the centersof mass of particles j and k and a reference axis. Thebond orientation correlation function is then defined asan average of the local bond orientation order: g ( r ) = (cid:42) N (cid:80) k (cid:54) = j Ψ( r j ) ∗ Ψ( r k ) δ ( r − | r j − r k | ) (cid:43)(cid:42) N (cid:80) k (cid:54) = j δ ( r − | r j − r k | ) (cid:43) (7)A power-law decay of g ( r ) means that there is a quasi-long-range orientation correlation. For the system in anhexatic phase g ( r ) ∝ r − η where 0 < η < . • Kosterlitz, Thouless, Halperin, Nelson and Young(
KT HN Y ) scenario [61–63]: a two-step continuousphase transition, first from liquid to hexatic phaseand then from hexatic to solid phase [64, 65], • Two step transition: first order phase transitionfrom liquid to hexatic phase and then continuousphase transition from hexatic to solid phase [9, 10], • First order phase transition from liquid to solidphase [66].We will show in the next section how finite size effectsand access to the thermodynamic regimes can bring un-certainty regarding phase behavior of the system; mainlythe nature of phase transition where these issues are ex-tensively reported in the literature about the nature ofmelting transition [9, 67].
II. SIMULATION DETAILS
The adsorption of hard-core molecules with a thirdneighbor exclusion range on the triangular lattice in-volves the adsorption of molecules covering 7 adsorp-tion sites in the manner as represented by a red circlein Figure. 1. As it is clear from the hexagon drawnaround the adsorbate (red circle) in Figure. 1(a)-(b), thismodel could have two different orientations at high cov-erage in comparison to other adsorbates such as hard-core molecule with first neighbor exclusion [54]. Here,we employ two complementary methods as in our pre-vious work [54]: an “adsorption method”, which beginsfrom an empty lattice, and a “desorption method”, whichbegins with a full lattice and progressively decreases cov-erage. The results are expected to bracket the correctequilibrium equation of state. (a) (b)
FIG. 1: Triangular lattice with first, second and third neighborexclusion where each adsorbate covers 7 sites. The center of theadsorbate is represented by a circle, arrows indicate possible dis-placements of particles, and stars represent the sites where thecenter of other particles are not allowed to adsorb. (a) and (b)show two different configuration at the maximum close packing.
In the adsorption method, molecules or particles areprogressively added to an initially empty d × d latticesurface where a periodic boundary condition is used toameliorate finite size effects. The only restriction is thatoverlap is not allowed; an assumption based physically onshort-range electrostatic repulsion. For each adsorptionattempt, a random position ( x, y ) is selected representingthe center of mass of the particle. If the selected site andits neighbors are empty, adsorption is accepted. Other-wise it is rejected. Diffusion, the simultaneous movementof particles, is introduced sequentially with a predefinedratio D between the number of diffusion attempt andthe adsorption attempt: For D = 3 each adsorption isfollowed by 3 diffusion attempts, etc. For each diffu-sion attempt, a previously adsorbed particle and a di-rection for the displacement of the particle are selectedrandomly; an arrow in Figure. 1 illustrates the possibledirection. If moving the center of the mass of the particleto the next node along this direction does not infringe thenon-overlap condition, diffusion is accepted. Otherwiseit is rejected. In the RSAD model, when diffusion is fastenough, the surface layer is at internal equilibrium (evenduring transient adsorption) and the blocking functioncan be considered as a state function.For the desorption method the lattice is initially full.For each simulation step, two particles are randomly se-lected and removed. Then one adsorption attempt and D diffusion attempts are performed following the same pro-cedure as for the adsorption method. The choice of thesequence (2 desorption events followed by 1 adsorption)is arbitrary but answers the need at each time step to de-crease coverage and add at least one particle to calculatethe blocking function.For both adsorption and desorption methods, theblocking function is extracted from the success rate ofadsorption attempts. 500 runs are performed, and anensemble average is used to reduce the noise arising fromthe numerical calculation of the derivative of the cov-erage. The blocking function is fitted with a polyno-mial function before it is used to generate the adsorptionisotherm. The latter is inserted into the Gibbs adsorptionisotherm equation to obtain the equation of state.In this work we also used a relaxation method in or-der to track the structure of particles by time and proofour hypothesis about expecting the correct equilibriumequation of state starting from either the adsorption ordesorption method. In this method, either the adsorptionor desorption method can be used to reach a specific cov-erage. In the second step, one adsorption attempt and D diffusion attempts are performed following the same pro-cedure as for the adsorption method. In order to keep thecoverage constant, if the adsorption attempt is successful,one particle is randomly selected and removed. 1500 runsare performed to extract the success rate of adsorptionattempts or, in other words, the blocking function. III. RESULTS
Accessing the thermodynamic regime is an initial steptoward studying the phase behavior of the system [9, 10]where some algorithms are insufficient to equilibrate thesystem in order to study the phase transition of the sys-tem in a reliable manner [64]. The effect of surface diffu-sion for a lattice size d = 105 is studied in Figure. 2. Ini-tially, when the system is diluted, all of the curves regard-less of their methods or their surface diffusion overlap inthe low surface coverage as presented in Figure. 2(a). Thefluctuation in the phase transition region is much largerthan the pure liquid and solid phase, where the differenceis maximized in the middle of phase transition as illus-trated in Figure2(b) [66]. As presented in Figure2 (b),at high surface coverage probability of success of adsorb-ing a new particle initially decreases due to the cagingeffect. However by increasing the ordering of particles,this caging effect will be diminished in order to maximizethe available surface for accepting the new incoming par-ticles. The system reaches the equilibrium state whentwo curves overlap in the whole range of fractional sur-face coverage, so accessing the thermodynamic regimewill be apparent. (a) (b) FIG. 2: (a) Effect of surface diffusion on adsorption rate versussurface coverage for d = 105. Ads and Des refer to adsorptionand desorption methods, respectively. (b) The inset expands thehigh-coverage region where sensitivity to surface diffusion ap-pears. Figure. 3 represents the relaxation of blocking functionof hard-core molecules with third neighbor exclusion ona triangular lattice for D = 0 . d = 196 and θ = 0 . FIG. 3: Relaxation dynamics of hard core molecule with thirdneighbor exclusion on triangular lattice for surface coverage of θ = 0 .
915 and d = 196. D = 0 .
01 in both first and second steps.Adsorption and desorption method refer to the initial configura-tion to reach the surface coverage of 0 . Finite size effects are another important key issue infinding the equation of state in lattice simulations. Fig-ure. 4 presents the variation in adsorption rate for tri-angular lattices covering 7 sites of sizes 105 to 266. Forsurface diffusion of D = 1, the same blocking functionis obtained from both adsorption and desorption meth- (a) (b) FIG. 4: (a) Effect of lattice size on adsorption rate versus sur-face coverage for D = 1. Ads and Des refer to adsorption anddesorption methods, respectively. (b) the inset expands the high-coverage region where sensitivity to lattice size appears. (a) (b) FIG. 5: (a) Comparison between our EOS where d = 196 and D = 80 with Orban and Bellman (Matrix method and low andhigh density method) [33] for hard core molecules on a triangularlattice covers 7 sites. (b) The inset shows a magnified view oferror bar between adsorption and desorption method for ( d =196, D = 80) and ( d = 105, D = 80) in the phase transition region. ods at low surface coverage as presented in Figure. 4(a).Initially, in the vicinity of the phase transition, by in-creasing the lattice size dimension, the adsorption rateincreases for both the adsorption and desorption methodas illustrated in Figure. 4(b). Exactly before all of thecurves overlap around full coverage, the adsorption ratedecreases by increasing the lattice dimension for bothmethods. The results in Figure. 4(b) show that the ad-sorption rate is significantly different for d = 105 in com-parison to d = 196 and d = 266. As a result d = 105is 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 theimportance of using a large enough system. One of theadvantages of using the RSAD method is we know howbig our system should be to ensure that the results areat the same time accurate and computationally less ex-pensive.Although an exact solution of this model does not existin the literature, meaning there is no certain agreement (a) (b)
FIG. 6: (a) Analysis of phase transition region where d = 196 and D = 80, (b) d = 105 and D = 80. upon the equation of state and its phase transition, wecan compare our results with the analytic calculation ofOrban and Bellmans [33] for d = 196 and D = 80 inFigure. 5(a). This paper used two different methods, thematrix method, based on the sequence of exact solutionsfor lattices of infinite length and increasing finite width,and the series expansion method, based on knowing thefinal structure at close packing and constructing the den-sity and activity series to find the surface pressure. A firstorder transition was found for both the matrix and theseries expansion methods, where the transition occurs atsurface pressure of 4 . − . ± . ± .
75, where no phase transition wasfound.As illustrated in Figure. 5(a), at low surface coveragethere is no difference between the reported equations ofstate. However, in the lower part of the phase transitionzone, our equation of states shows lower surface pres-sure than Orban and Bellman’s series expansion method(low and high density) but higher surface pressure thantheir matrix method. On the contrary, in the upperpart of the phase transition region, the matrix method g ( r ) ∝ g ( r ) ∝ r − g ( r ) ∝ g ( r ) ∝ g ( r ) ∝ g ( r ) ∝ g ( r ) ∝ g ( r ) ∝ r − g ( r ) ∝ g ( r ) ∝ r − (a) (b) (c-1) (c-2) (d-1) (d-2) (e-1) (e-2) (f-2) (f-1) FIG. 7: Bond orientation correlation function of the hard core molecule with third neighbor exclusion based on relaxation method where d = 196 and D = 0 .
01 in both step of relaxation at different surface coverage (a) θ = 0 .
75, (b) θ = 0 .
85, (c-1) and (c-2) θ = 0 . θ = 0 . θ = 0 . θ = 0 .
98. Before relaxation refer to the configuration obtained fromadsorption method at D = 0 .
01 in the first step of relaxation method , and after relaxation refer to the equilibrium configuration. shows higher surface pressure than what we found, butthe series expansion method overlaps with both of ourmethods. The thermodynamically-stable surface pres-sure loop is observed in our simulation results as is re-ported extensively in the literature due to finite size ef-fects [9, 10, 68, 69]. Surface pressure can create a ther-modynamically stable loop at equilibrium for a finite sizesystem, but this loop will disappear at infinite size andthe coexistence zone will be flat. Creation of the flatpressure is visible in Figure. 5(b) for d = 196 around asurface pressure of 4 . ± .
05. Moreover, based on thelower part of the phase transition in Figure. 4(b), theadsorption rate tends to increase by increasing the lat-tice dimension, which means less surface pressure will beexpected for a larger lattice size. Conversely, the ad-sorption rate tends to decrease for both the adsorptionand desorption methods in the upper part of the phasetransition, which means higher surface pressure will beexpected for a larger lattice size. This behavior indicatesthe tendency of the system towards flatness at an infi-nite size. From the equality shown in the hatched area of Figure. 5(b), the horizontal surface pressure is obtainedfrom the Maxwell construction, where overlapping of thisconstruction line with the flat region of the equation ofstate confirms the tendency of the system to be flat atinfinite size. Figure. 5(b) shows that by increasing thelattice dimension, surface pressure decreases in the lowerpart of the phase transition and increases in the upperpart of phase transition. Bernard et al [9, 10] reportedthe same trend for the equation of state of the meltingtransition of hard disks by increasing the number of parti-cles where they reported a two step transition; first orderliquid-hexatic and continuous hexatic-solid transition fortheir system.The phase transition zone is studied in Figure. 6 basedon the derivative of surface pressure with respect to sur-face coverage in the adsorption method. One can seefrom Figure. 6(b) that for an insufficiently large system,the first order liquid-solid phase transition is obtainedfor d = 105 and D = 80; however, a different phasetransition is obtained for d = 196 and D = 80 in Fig-ure. 6(a), which indicates the importance of using a suf- Before relaxation After relaxation Θ = . Θ = . Θ = . Before relaxation After relaxation Θ = . Θ = . Θ = . ( a ) ( b ) ( c ) ( d ) ( e ) (f ) FIG. 8: Local bond orientation order based on relaxation method at low surface diffusion D = 0 .
01 and different surface coverage: (a) θ = 0 .
75, (b) θ = 0 .
85, (c) θ = 0 . θ = 0 . θ = 0 . θ = 0 .
98. Before relaxation refer to the configuration obtainedfrom adsorption method at D = 0 .
01 in the first step of relaxation method , and after relaxation refer to the equilibrium configuration. ficiently large system. Phase transition peaks obtainedfrom Figure. 6(a) for d = 196, D = 80 are analyzedthrough a bond orientation correlation function, g6(r),in Figure 7 based on the relaxation method for each in-dividual configuration. For additional insight into thephase transition region, each of these configurations isvisualized via the local bond orientation order functionΨ( r ), based on the relaxation method, in Figure. 8. Thefirst peak in the phase transition curve corresponds tosurface coverage of 0 . g ( r ) before and after relaxation atsurface coverage of 0 .
75 in Figure. 7(a) suggests that thesystem is in a pure liquid regime below a surface cover-age of 0 . . d = 105 and d = 196 through fractional surfacecoverages up to 0 . . . .
869 to see the behavior of the system below andabove the transition value. For surface coverage of 0 . g ( r ) decays exponentially at the same rate before and af-ter relaxation in Figure. 7(b); more peaks were observedafter relaxation, which is an indication of the creation ofa more ordered phase in the system. Creation of a moreordered phase is also confirmed by the local bond orien-tation order seen in Figure. 8(b), indicated by more redspread through the surface. For surface coverage of 0 . g ( r ) still decays exponentially, but at different rates be-fore and after relaxation as illustrated in Figure. 7 (c-1)and (c-2), respectively. The local bond orientation orderin Figure. 8(c) shows that ordered particles prefer to stickto each other. The third peak in Figure. 6 correspondsto surface coverage of 0 . g ( r ) decays witha power law with θ = 0 .
25 as illustrated in Figure. 7(d-2) , which is an indication of a hexatic phase. The localbond orientation order shown in Figure. 8(d) reveals that t=0 t=21 t=211 t=2113 t=2958 (a) (b)
FIG. 9: Local bond orientation order for surface coverage of 0.915 with D=0.01 and d=196, (a) for adsorption method over time, (b) fordesorption method over time. at this coverage all the particles tend to stick togetherand they are between the mobile particles, as was alsoreported in reference [3]. The tendency of ordered parti-cles to stick together causes the probability of success toincrease in Figure. 2(b) due to maximizing the availablefree surface for accepting the new incoming particles. Byincreasing the surface coverage from 0.915 to 0.963, η de-creases from 0 .
25 to 0 .
08 after relaxation as illustratedin Figure. 7(e-2), which indicates that the fourth peakin Figure. 6 corresponds to a first order transition fromhexatic phase to solid phase in equilibrium state.To further illustrate the relaxation dynamics towardan equilibrium state, the local bond orientation orderparameter of particles at surface coverage of 0 .
915 wastracked over time for two different initial configurationsobtained from adsorption and desorption methods atvery low surface diffusion ( D = 0 . d = 196). Fig-ure. 9(a) indicates that the density should be increasedvery slowly when we start from an empty lattice, or elsethe system will be locked into a glassy configuration [66].A glassy state is reported in reference [3] during rapidcompression of a system composed of spherical particles.For the desorption method, Figure. 9(b), the system isinitially ordered, so the density should be decreased veryslowly or else the system will be metastable [66]. As timepasses in the equilibrium state, phase separation occurs,where ordered particles tends to cluster and create morespace for further adsorption. IV. CONCLUSION
In this paper we studied the phase behavior of a hard-core molecules with third neighbor exclusions in a trian-gular lattice with
RSAD simulations, in order to derive the equation of state of a two-dimensional hard-core par-ticle based on kinetic arguments and the Gibbs adsorp-tion isotherm. We compared our results with Orban andBellman [33], who used a matrix method and a seriesexpansion methods in low and high densities, and foundonly partial agreement. Our results show that the systemis in a pure liquid regime below surface coverage of 0 . .
877 and 0 . . g ( r ) decays algebraically afterrelaxation with a power-law exponent η = 0 .
25, which isan indication of a hexatic phase. Our simulation resultsreveal that as the surface coverage increases, after relax-ation, the surface particles tend to form tightly-packedclusters in like-oriented domains, while the remainingmobile particles have more random orientations. By in-creasing the surface coverage above 0 . η decreasesfrom 0 .
25 toward zero after relaxation, where the systemundergoes a first order transition from a hexatic phase toa solid phase.One of the advantages of using the
RSAD model isbeing able to locate the equilibrium state, which assuresus that adequate thermalization and finite size are be-ing used to reach the equilibrium state. Moreover, subtledetails of the clustering structure, through direct visu-alization of the system using the relaxation method atany fractional surface coverage, provides insight regard-ing coexistence regions and phase transitions.
ACKNOWLEDGMENTS
This research was supported by the National ScienceFoundation under grant
N o. [1] Vincent Pauchard, Jayant P Rane, and Sanjoy Banerjee,“Asphaltene-laden interfaces form soft glassy layers incontraction experiments: A mechanism for coalescenceblocking,” Langmuir , 12795–12803 (2014).[2] Fang Liu, Shaghayegh Darjani, Nelya Akhmetkhanova,Charles Maldarelli, Sanjoy Banerjee, and VincentPauchard, “Mixture effect on the dilatation rheology ofasphaltenes-laden interfaces,” Langmuir , 1927–1942(2017).[3] Zhanglin Hou, Kun Zhao, Yiwu Zong, and Thomas GMason, “Phase behavior of two-dimensional browniansystems of corner-rounded hexagons,” Physical ReviewMaterials , 015601 (2019).[4] Sara Fortuna, David L Cheung, and Alessandro Troisi,“Hexagonal lattice model of the patterns formed byhydrogen-bonded molecules on the surface,” The Jour-nal of Physical Chemistry B , 1849–1858 (2010).[5] UK Weber, VM Burlakov, LMA Perdigao, RHJ Fawcett,PH Beton, NR Champness, JH Jefferson, GAD Briggs,and DG Pettifor, “Role of interaction anisotropy in theformation and stability of molecular templates,” Physicalreview letters , 156101 (2008).[6] VA Gorbunov, SS Akimenko, AV Myshlyavtsev,VF Fefelov, and MD Myshlyavtseva, “Adsorption oftriangular-shaped molecules with directional nearest-neighbor interactions on a triangular lattice,” Adsorption , 571–580 (2013).[7] Per Bak, P Kleban, WN Unertl, J Ochab, G Akinci,NC Bartelt, and TL Einstein, “Phase diagram of se-lenium adsorbed on the ni (100) surface: A physical re-alization of the ashkin-teller model,” Physical review let-ters , 1539 (1985).[8] David E Taylor, Ellen D Williams, Robert L Park,NC Bartelt, and TL Einstein, “Two-dimensional order-ing of chlorine on ag (100),” Physical Review B , 4653(1985).[9] Etienne P Bernard and Werner Krauth, “Two-step melt-ing in two dimensions: First-order liquid-hexatic transi-tion,” Physical review letters , 155704 (2011).[10] Michael Engel, Joshua A Anderson, Sharon C Glotzer,Masaharu Isobe, Etienne P Bernard, and WernerKrauth, “Hard-disk equation of state: First-order liquid-hexatic transition in two dimensions with three simula-tion methods,” Physical Review E , 042134 (2013).[11] Irving Langmuir, “The adsorption of gases on plane sur-faces of glass, mica and platinum.” Journal of the Amer-ican Chemical society , 1361–1403 (1918).[12] Jayant P Rane, Sharli Zarkar, Vincent Pauchard,Oliver C Mullins, Dane Christie, A Ballard Andrews, An-drew E Pomerantz, and Sanjoy Banerjee, “Applicabilityof the langmuir equation of state for asphaltene adsorp-tion at the oil–water interface: Coal-derived, petroleum,and synthetic asphaltenes,” Energy & Fuels , 3584–3590 (2015).[13] Jayant P Rane, Vincent Pauchard, Alexander Couzis,and Sanjoy Banerjee, “Interfacial rheology of asphaltenesat oil–water interfaces and interpretation of the equationof state,” Langmuir , 4750–4759 (2013). [14] Keng Yuen Foo and Bassim H Hameed, “Insights intothe modeling of adsorption isotherm systems,” Chemicalengineering journal , 2–10 (2010).[15] Jian-Sheng Wang, Peter Nielaba, and Vladimir Priv-man, “Collective effects in random sequential adsorptionof diffusing hard squares,” Modern Physics Letters B ,189–196 (1993).[16] Berardo M Manzi, Marco Werner, Elena P Ivanova, Rus-sell J Crawford, and Vladimir A Baulin, “simulationsof protein adsorption on nanostructured surfaces,” Sci-entific reports (2019).[17] I Lonˇcarevi´c, Lj Budinski-Petkovi´c, SB Vrhovac, andA Beli´c, “Adsorption, desorption, and diffusion of k-merson a one-dimensional lattice,” Physical Review E ,021115 (2009).[18] I Lonˇcarevi´c, Lj Budinski-Petkovi´c, and SB Vrhovac,“Reversible random sequential adsorption of mixtureson a triangular lattice,” Physical Review E , 031104(2007).[19] J Talbot, G Tarjus, PR Van Tassel, and P Viot, “Fromcar parking to protein adsorption: an overview of se-quential adsorption processes,” Colloids and Surfaces A:Physicochemical and Engineering Aspects , 287–324(2000).[20] James W Evans, “Random and cooperative sequentialadsorption,” Reviews of modern physics , 1281 (1993).[21] PL Krapivsky and E Ben-Naim, “Collective properties ofadsorption–desorption processes,” The Journal of chem-ical physics , 6778–6782 (1994).[22] Jian-Sheng Wang, Peter Nielaba, and Vladimir Privman,“Locally frozen defects in random sequential adsorptionwith diffusional relaxation,” Physica A: Statistical Me-chanics and its Applications , 527–538 (1993).[23] JR ˇS´cepanovi´c, D Stojiljkovi´c, ZM Jakˇsi´c, Lj Budinski-Petkovi´c, and SB Vrhovac, “Response properties in theadsorption–desorption model on a triangular lattice,”Physica A: Statistical Mechanics and its Applications , 213–226 (2016).[24] G Tarjus and P Viot, “Statistical mechanical descriptionof the parking-lot model for vibrated granular materials,”Physical Review E , 011307 (2004).[25] Matthew A Meineke and J Daniel Gezelter, “Randomsequential adsorption model for the differential coverageof gold (111) surfaces by two related silicon phthalocya-nines,” The Journal of Physical Chemistry B , 6515–6519 (2001).[26] P Schaaf and J Talbot, “Surface exclusion effects in ad-sorption processes,” The Journal of chemical physics ,4401–4409 (1989).[27] Jeremy J Ramsden, “Observation of anomalous diffusionof proteins near surfaces,” The Journal of Physical Chem-istry , 3388–3391 (1992).[28] S Ji, C Ates, and I Lesanovsky, “Two-dimensional ryd-berg gases and the quantum hard-squares model,” Phys-ical review letters , 060406 (2011).[29] Markus Lackinger, Stefan Griessl, Thomas Markert, Fer-dinand Jamitzky, and Wolfgang M Heckl, “Self-assemblyof benzene- dicarboxylic acid isomers at the liquid solid interface: steric aspects of hydrogen bonding,” The Jour-nal of Physical Chemistry B , 13652–13655 (2004).[30] Yongsheng Zhang, Volker Blum, and Karsten Reuter,“Accuracy of first-principles lateral interactions: Oxygenat pd (100),” Physical Review B , 235406 (2007).[31] Rodney J Baxter, Exactly solved models in statistical me-chanics (Elsevier, 2016).[32] Andr´e Bellemans and RK Nigam, “Phase transitions intwo-dimensional lattice gases of hard-square molecules,”The Journal of Chemical Physics , 2922–2935 (1967).[33] J Orban and A Bellemans, “Phase transitions in two-dimensional lattice gases of hard-core molecules. the tri-angular lattice,” The Journal of Chemical Physics ,363–370 (1968).[34] Francis H Ree and Dwayne A Chesnut, “Phase transitionof a hard-core lattice gas. the square lattice with nearest-neighbor exclusion,” The Journal of Chemical Physics , 3983–4003 (1966).[35] LK Runnels, LL Combs, and James P Salvant, “Exactfinite method of lattice statistics. ii. honeycomb-latticegas of hard molecules,” The Journal of Chemical Physics , 4015–4020 (1967).[36] LK Runnels and LL Combs, “Exact finite method of lat-tice statistics. i. square and triangular lattice gases ofhard molecules,” The Journal of Chemical Physics ,2482–2492 (1966).[37] Xiaomei Feng, Henk WJ Bl¨ote, and Bernard Nienhuis,“Lattice gas with nearest-and next-to-nearest-neighborexclusion,” Physical Review E , 061153 (2011).[38] Z Rotman and E Eisenberg, “Critical exponents fromcluster coefficients,” Physical Review E , 031126(2009).[39] D. S. Gaunt and M. E. Fisher, “Hard-sphere lattice gases.i. plane-square lattice,” J. Chem. Phys. , 2840–2863(1965).[40] David S Gaunt, “Hard-sphere lattice gases. ii. plane-triangular and three-dimensional lattices,” The Journalof Chemical Physics , 3237–3259 (1967).[41] E Eisenberg and A Baram, “A first-order phase transitionand a super-cooled fluid in a two-dimensional lattice gasmodel,” EPL (Europhysics Letters) , 900 (2005).[42] MV Ushcats, “High-density equation of state for a latticegas,” Physical Review E , 052144 (2015).[43] MV Ushcats, LA Bulavin, VM Sysoev, and SJ Ushcats,“Virial and high-density expansions for the lee-yang lat-tice gas,” Physical Review E , 012143 (2016).[44] ER Cowley, “A bethe model for phase transitions in thehard sphere lattice gas,” The Journal of Chemical Physics , 458–461 (1979).[45] DM Burley, “A lattice model of a classical hard spheregas,” Proceedings of the Physical Society , 262 (1960).[46] Hendrik Hansen-Goos and Martin Weigt, “A hard-spheremodel on generalized bethe lattices: Statics,” Journalof Statistical Mechanics: Theory and Experiment ,P04006 (2005).[47] K Binder and DP Landau, “Phase diagrams and crit-ical behavior in ising square lattices with nearest-andnext-nearest-neighbor interactions,” Physical Review B , 1941 (1980).[48] Dwayne A Chesnut, “A monte carlo study of the thermo-dynamic properties of “hard hexagons” on the triangularlattice,” Journal of Computational Physics , 409–434(1971). [49] Da-Jiang Liu and James W Evans, “Ordering and per-colation transitions for hard squares: Equilibrium versusnonequilibrium models for adsorbed layers with c (2 ×
2) superlattice ordering,” Physical Review B , 2134(2000).[50] Heitor C Marques Fernandes, Jeferson J Arenzon, andYan Levin, “Monte carlo simulations of two-dimensionalhard core lattice gases,” The Journal of chemical physics , 114508 (2007).[51] Trisha Nath and R Rajesh, “Multiple phase transitionsin extended hard-core lattice gas models in two dimen-sions,” Physical Review E , 012120 (2014).[52] George Stanley Rushbrooke and HI Scoins, “On the isingproblem and mayer’s cluster sums,” Proceedings of theRoyal Society of London. Series A. Mathematical andPhysical Sciences , 74–90 (1955).[53] Luis Lafuente and Jos´e A Cuesta, “Density functionaltheory for nearest-neighbor exclusion lattice gases in twoand three dimensions,” Physical Review E , 066120(2003).[54] Shaghayegh Darjani, Joel Koplik, and VincentPauchard, “Extracting the equation of state of latticegases from random sequential adsorption simulations bymeans of the gibbs adsorption isotherm,” Physical Re-view E , 052803 (2017).[55] Heitor C Marques Fernandes, Yan Levin, and Jefer-son J Arenzon, “Equation of state for hard-square latticegases,” Physical Review E , 052101 (2007).[56] Luis Lafuente and Jos´e A Cuesta, “Phase behavior ofhard-core lattice gases: A fundamental measure ap-proach,” The Journal of chemical physics , 10832–10843 (2003).[57] Kabir Ramola and Deepak Dhar, “A high density per-turbation expansion for the hard square lattice gas,” in AIP Conference Proceedings , Vol. 1447 (AIP, 2012) pp.135–136.[58] Dipanjan Mandal, Trisha Nath, and R Rajesh, “Esti-mating the critical parameters of the hard square latticegas model,” Journal of Statistical Mechanics: Theory andExperiment , 043201 (2017).[59] Trisha Nath and R Rajesh, “The high density phase ofthe k-nn hard core lattice gas model,” Journal of Statis-tical Mechanics: Theory and Experiment , 073203(2016).[60] Da-Jiang Liu and JW Evans, “Lattice-gas modeling ofthe formation and ordering of oxygen adlayers on pd (10 0),” Surface science , 13–26 (2004).[61] AP Young, “Melting and the vector coulomb gas in twodimensions,” Physical Review B , 1855 (1979).[62] David R Nelson and BI Halperin, “Dislocation-mediatedmelting in two dimensions,” Physical Review B , 2457(1979).[63] John Michael Kosterlitz and David James Thouless,“Ordering, metastability and phase transitions in two-dimensional systems,” Journal of Physics C: Solid StatePhysics , 1181 (1973).[64] CH Mak, “Large-scale simulations of the two-dimensionalmelting of hard disks,” Physical Review E , 065104(2006).[65] Yi Peng, Ziren Wang, Ahmed M Alsayed, Arjun G Yodh,and Yilong Han, “Melting of colloidal crystal films,”Physical review letters , 205703 (2010).[66] BJ Alder and TE Wainwright, “Phase transition in elas-tic disks,” Physical Review , 359 (1962). [67] Keola Wierschem and Efstratios Manousakis, “Simula-tion of melting of two-dimensional lennard-jones solids,”Physical Review B , 214108 (2011).[68] Joseph E Mayer and Wm W Wood, “Interfacial tensioneffects in finite, periodic, two-dimensional systems,” TheJournal of Chemical Physics , 4268–4274 (1965). [69] Manuel Schrader, Peter Virnau, and Kurt Binder, “Sim-ulation of vapor-liquid coexistence in finite volumes: Amethod to compute the surface free energy of droplets,”Physical Review E79