Equilibrium properties of the lattice fluid with the repulsion between the nearest neighbors on the two-level lattice with nonrectangular geometry

The equilibrium properties of the lattice fluid with the repulsion between the nearest neighbors on the two-level planar triangular lattice are investigated. The numerical results obtained from the analytical expressions are compared with the Monte Carlo simulation data. It is shown that the previously proposed diagrammatic approximation makes it possible to determine the equilibrium characteristics of the lattice fluid with the repulsion between the nearest neighbors on a two-level lattice with an accuracy comparable to the accuracy of modelling the system using the Monte Carlo method in the entire range of thermodynamic parameters. It was found that, in contrast to a similar one-level system, a lattice fluid with the repulsion between the nearest neighbors undergoes a first-order phase transition.


Introduction
The model of the lattice gas or the lattice liquid is widely used in the description of processes occurring on the surfaces and in the volumes of solids [1], as well as in the study of various electrochemical systems [2][3][4][5][6].
The main idea of this model consists in the decomposition of the original physical system into two spatial interpenetrating subsystems, one of which is rather rigidly structured and plays the role of the reference one, and in the second one the particles retain high mobility. In fact, the particles of the second subsystem move in the field created by the supporting subsystem. In this case, it is assumed that the given field does not depend on the distribution of the moving particles. Thus, it can be argued that the reference subsystem creates a potential with a spatially periodic profile for the moving particles, in which there are sufficiently deep potential energy minima that form one or other periodic lattice.
Another condition for the possibility of using the lattice fluid model to describe a real physical system is the presence of two significantly different time scales. The first of which, 1 , is determined by the average residence time of the mobile particle in the minimum of potential energy, and the second, 2 , is determined by the oscillation period of the mobile particle near this position. Hence, for the indicated time scales, we have 1 2 . The fulfillment of the above conditions, the stability of the reference system and the presence of two significantly different time scales, allows us to restrict ourselves to considering only the subsystem of moving particles, and the characteristics of the potential landscape can be considered as the input parameter of the model, assuming them to be given.
The potential energy minima will be determined as the lattice sites, each of which can be either occupied by a single particle or can be vacant. When considering a system of charged particles occupying the lattice sites, the number of possible states of each lattice site obviously increases [7] in comparison

Model
The crystal plane (1, 1, 1) of the crystal with the simple cubic Bravais lattice can be considered as the basis for constructing a two-level lattice system. In this plane, the atoms of the crystal form a triangular lattice as shown in figure 1, where these atoms are represented by open squares. The interaction between the atoms of the plane forms an energy profile, the minima of which are the locations of the particles adsorbed on a given plane of the crystal, the so-called adparticles. The ensemble of these adparticles forms the system under study.
The geometric features of the original lattice make it possible to distinguish two types of sites of the lattice model. Some of them --sites -are located between two neighboring atoms of the crystal plane and are shown in figure 1 with open circles. The second type of sites --sites -are in the center of the regular triangle built on three neighboring atoms. In figure 1, -sites are shown as closed circles.
In the general case, the depths of potential wells corresponding to the -and -sites can be different where , denote the depth of the potential well corresponding to the lattice sites of the indicated type, and is the difference in the depths of the potential wells.
Since the and sites are not energetically equivalent, the probability of finding the particle in each of them can be different. For a theoretical description of this system, the original lattice can be divided into the system of two sublattices, each of which contains the sites of only one type. This allows us to speak in what follows about the and sublattices, respectively.
Thus, the constructed lattice system is a combination of two two-dimensional hexagonal sublattices. In this case, the -sublattice has a lattice constant /2, and the -sublattice has a lattice constant / √ 3, where is the distance between the atoms of the crystal plane.
It can also be noted that the first coordination number, i.e., the number of the nearest neighbor nodes is different for the sublattices introduced above. Thus, for example, the sublattice site has = 2 nearest-neighbor sites. Moreover, all these nearest neighboring sites belong to the sublattice. On the other hand, for a -site, the number of nearest sites on the sublattice is = 3.
In this paper, we consider the lattice fluid with the interaction between nearest neighbors. The energy of this interaction is assumed to be equal to 1 . It should be noted that the particles occupying the lattice sites of the same type do not interact with each other, because these nodes are not nearest neighbors.
Counting the sites on each of the sublattices is not difficult. As noted above, the -site is located between every two nearest atoms of the crystal plane. If we assume that the initial crystal plane contains atoms forming the triangular lattice with the first coordination number = 6, then the number of pairs of nearest atoms, and hence the number of -sites, is The -site lies in turn in the center of each triangular graph built on triples of the nearest atoms. The method for calculating the corresponding weight coefficients, meaning the number of different graphs of a given type that can be built on the lattice of the selected type, is described in detail in [22]. The weighting coefficient of the three-vertex graph on the triangular lattice, and hence the number of -nodes, is Thus, in the case of a lattice containing lattice sites, we obtain For further calculations, it is convenient to introduce the average first coordination number 1 which has the meaning of the average number of the nearest neighboring sites for an arbitrary lattice node. For the lattice under consideration, the coordination number 1 will not be an integer. It can be determined by the following relation . (1.5) Thus, the original lattice can be considered as the system of two nonequivalent sublattices built on the same lattice sites. For each of the sublattices, the concentration of particles 1 , and vacancies 0 can be determined, and the order parameter of the model where ( = , ) is the number of adparticles on the sublattice , and is the average concentration of adparticles ( = 1) and vacancies ( = 0) in the system.
The thermodynamic state of the constructed lattice model can be defined by specifying the set of occupation numbers taking values = 1 or = 0, if the -th lattice site is occupied by the particle or is vacant, respectively.
For a given set of occupation numbers, the potential energy of the system of particles on the lattice can be represented as: where , and denote a summation over all -, -sites and over all pairs of nearest-neighbor sites, respectively.
In what follows, an additional condition will be imposed on the considered lattice system (1.11) Compliance with condition (1.11) ensures the symmetry of the phase diagram around = 0.5.

Monte Carlo simulation algorithm
Determination of the equilibrium properties of the lattice system can be performed within the framework of the standard Metropolis algorithm [24], the application of which to the lattice fluids of various types is described in detail, for example, in [25][26][27][28].
In this case, the modelling of the equilibrium properties of the system is carried out in the grand canonical ensemble, i.e., with a variable number of particles in the system. Within the framework of this algorithm, the value of the chemical potential and the interaction energy 1 is fixed, and the concentration of particles on the lattice and their distribution are determined directly during the simulation. In this case, the initial distribution of the particles over the -site lattice can be arbitrary. After choosing an arbitrary lattice site and changing its state (adding a particle if the site is vacant, and removing it if the site is occupied), the change in the energy corresponding to the given change in the state of the site is determined as where the plus corresponds to the addition of a particle, and the minus to its removal; 1 is the number of particles occupying the nearest lattice sites; is the Kronecker symbol; index denotes the type of sublattice to which an arbitrarily chosen node belongs.
It should be noted that the maximum number of nearest neighbours for the and sublattices are different; i.e., depending on its position, the considered lattice site has a different number of the nearest neighboring nodes. To simplify the program code of the simulation procedure, the lattice system can be supplemented with ghost sites, as shown in figure 2. Ghost lattice sites are represented in this figure by closed squares. In this case, some of the introduced ghost sites coincide in position with the atoms of the crystal plane, which is reflected by the position of the closed squares inside the open ones, representing the atoms of the plane (figure 1).
The introduction of ghost sites makes it possible to transform the original lattice structure into the "square" lattice, in which each lattice site has 4 nearest neighbors. In this case, any site of the sublattice will have 2 fictitious sites among its nearest neighbors, and a site of the sublattice will have only one neighboring ghost site.
Thus, each node of the modelled system ( , ) can be in one of 3 possible states: ( , ) = 1 corresponds to a site ( , ) occupied by the particle, ( , ) = 0 -to a vacant lattice site, ( , ) = −1 -to a ghost lattice site. Obviously, during the simulation, the state of the ghost lattice sites does not change. In case of a random choice of the ghost lattice node, the selection procedure is not taken into account and a new node is selected.
If the change in the state of the lattice site leads to a decrease in the energy of the system ( 0), then the new configuration is accepted. If the energy of the system increases ( > 0), then the change in the state of the site is accepted with the probability where is the temperature; B is the Boltzmann constant.

13501-5
In the latter case, to accept or reject the proposed change in the state of the lattice site, a random number is generated from the interval [0; 1], and if , then the new configuration is accepted. Otherwise, the system reverts to its previous state. Repetition of the described procedure of times, where is equal to the number of particles on the lattice, forms one step of the Monte Carlo algorithm (MC step).
In this work, the modelled system contains 128 × 128 = 2 14 lattice sites. Of these, 5 × 2 11 sites are lattice sites of the two-level system, and 3 × 2 11 are the ghost lattice nodes. The complete simulation procedure consists of 30000 MC steps. Since the modelling procedure begins with a random distribution of particles on the lattice, the first 10000 MC steps are assigned to the transition of the system from an arbitrary initial state to an equilibrium state and are not taken into account in further modelling. Starting from 10001 steps of the algorithm, the number of particles in the system, the number of particles on each of the sublattices, and the number of pairs of nearest neighboring sites filled with particles are determined. On the basis of these results, after the end of the modelling procedure, the average concentration of particles in the entire system, the concentration of particles on the sublattices, the order parameter of the system, its thermodynamic factor, and the probability of detecting the two nearest lattice sites with occupied particles.
To reduce the effect of the size of the system on the obtained values of the equilibrium characteristics, we used periodic boundary conditions.

Quasi-chemical approximation
Together with the investigated lattice system, a similar reference lattice system can be considered [29][30][31]. This reference system is determined by the one-particle average potentials ( ) ( ) of the interaction of a particle ( = 1) or vacancy ( = 0) located at site of the -sublattice with some site , which is a neighbor of order ( = , ; = 1 corresponds to nearest neighbors, etc.). The potential energy of the reference system can be as where is the number of neighboring sites of order for the lattice site of the -sublattice. Summation over index implies summation over the so-called coordination shells, i.e., over complete sets of neighbors of order .
All equilibrium characteristics of the original lattice system are determined by its partition function where Sp denotes the summation over all possible particle distributions given by sets of occupation numbers { }, taking into account the normalization conditions The partition function of the initial lattice system (3.2) can be represented as where (0) and (d) are the partition function of the reference system and the diagram part of the partition function, respectively: (3.5)

13501-6
Angular brackets with index 0 denote the averaging over the states of the reference system, the distribution functions of which, due to the one-particle structure of its potential energy, can be written as the product of the average concentrations of particles and vacancies.
The partition function of the reference lattice system can be calculated analytically On the other hand, the diagram part of the partition function can be represented as where ( ) are the Mayer-like functions [32] ( where sites and are neighbors of order ; is the interaction energy of particles occupying the -th order neighbors; , = , . Equation (3.7) can be rewritten as This series allows for a clear geometric interpretation. Each of its members is associated with the diagram (graph) consisting of the lines (link) corresponding to the functions ( ) connecting the lattice nodes and . The type of each link is determined by the index , i.e., the degree of proximity of nodes and . When averaging each diagram, the summation is performed over the lattice sites, leading to the appearance of the weight coefficient (weight) of the diagram, equal to the number of ways of placing the diagram on the lattice. The weight of the diagram is determined by the geometric properties of both the diagram itself and the lattice on which the diagram is placed. Equation (3.4) makes it possible to represent the free energy of the initial lattice per one lattice site as the sum of two terms where (0) and (d) are the reference and diagram parts of the free energy, respectively. We have where is the -th coordination number for the sublattice , i.e., the number of neighboring sites of order on the sublattice ( 1 = , 1 = ). To calculate the diagram part of the free energy, the corresponding partition function can be represented as

13501-7
This makes it possible to determine the diagram part of the free energy (d) by selecting from expansion (3.9) diagrams with weight coefficients linear in . Renormalization using the mean potentials of the Mayer functions improves the convergence of the expansion of the free energy of the system. The mean potentials themselves are determined from the principle of minimum susceptibility of the free energy to their variations [29] ( ) ( ) = 0. (3.14) In turn, the order parameter (1.7) can be found from the condition for the extremality of the free energy [26,31] ( , ) , ( ) = 0, (3.15) which is equivalent to the condition of equality of chemical potentials on the sublattices. Holding various parts of the expansion of the diagrammatic part of the free energy, one can construct several approximate methods for calculating the free energy of a lattice system. The quasi-chemical approximation (QChA) corresponds to taking into account the contribution of only two-vertex graphs in the diagrammatic part of the free energy.
The use of the principle of minimum susceptibility leads to the following system of equations for determining the mean potentials of the nearest neighbors [22] If the order parameter of the system is given, then the system of equations (3.16) has an analytical solution, the substitution of which in relation (3.10) leads to the following expression for the free energy of the lattice fluid in the quasi-chemical approximation, where Calculations show that, within the framework of the QChA, the range of the average potentials of the reference system coincides with the range of the interaction potential in the initial lattice system. Therefore, in the case of the lattice system with the interaction of the nearest neighbors, the average potentials for the second and more distant neighbors will be equal to zero In [22,29] it was shown that when the average potentials are determined in accordance with relation (3.16), the diagram part of the free energy (d) vanishes. Thus, in the framework of QChA, the free energy of the lattice model under study turns out to be equal to the free energy of the reference system (0) .

Diagram approximation
The most obvious way to refine the results of QChA is to consider additional terms in the expansion of the diagram part of the free energy. The implicitly proposed strategy can be realized in the framework of the so-called diagram approximation (DA) [30].
The main idea of this approximation is that the average potentials of the reference system are taken to be equal to their values obtained within the framework of the quasi-chemical approximation considered above. In this case, for the diagrammatic part of the free energy, an additional assumption is put forward that the contribution to the free energy of the entire diagram series is proportional to the contribution of the simplest nonzero graph in the quasi-chemical approximation.
In accordance with the properties of QChA, such a graph will be the simplest ring graph, i.e., a closed irreducible connected diagram containing the minimum possible number of links and vertices. In the case under consideration, such a graph contains 12 vertices, 6 of which belong to the -sublattice and 6 belong to the -sublattice and has a weight coefficient equal to 1/5. The graphically indicated object is shown in figure 1. Considering relations (3.16), the expression for the diagrammatic part of the free energy can be written as [22]  .

(4.2)
The information about the free energy of the lattice fluid makes it possible to determine any of its equilibrium properties. In particular, the chemical potential of the system , the thermodynamic factor , and the probability of two nearest sites to be occupied by particles (1; 1) are defined as (1; 1) = 2 . (4.5)

Results and discussion
The fulfillment of the condition (1.11) leads to the fact that in the case of the system with repulsion of the nearest neighbors, the number of shallower sites exceeds the number of deeper sites.

13501-9
At the same time, the critical parameter of the system turns out to be equal to the critical parameter of the similar system with attraction between particles [33]   Comparison of the results of analytical calculations and modelling showed that the parameter included in relation (4.2), which ensures the equality of the critical parameters of the system found during its modelling and calculated within the framework of the diagram approximation, is equal to 1.10308. It can be noted that its value differs very insignificantly from the analogous value obtained for a system with attraction between nearest neighbors [33]. The isotherms of the chemical potential calculated for temperatures above and below the critical temperature are show in figure 3. Figure 4 shows the isotherms of the chemical potential at a temperature of 0.8 determined within the framework of the quasi-chemical and diagram approximation.
It can be seen from the comparison (figure 3) that the results of the diagram approximation are in good agreement with the results of MC simulation in the entire range of the considered temperatures and concentrations. As expected, in the two-phase region, the approximate approaches give a Van der

13501-11
Waals loop (see figure 4), while the simulation of the system demonstrates the constancy of the chemical potential. This allows us to conclude that a first-order phase transition is observed in the case of a two-level lattice.
Maxwell's construction allows you to determine the points of phase transitions. The phase diagram of the system is shown in figure 5. The phase diagram turns out to be much narrower compared to the system with attractive interaction [33]. In general, it can be noted that the thermodynamic characteristics do not qualitatively reveal the new properties in comparison with a similar system with the attraction of the nearest neighbors.
By contrast, the structural properties -the order parameter (figure 6), the thermodynamic factor (figure 7) and the correlation function (figure 8) -show the new features qualitatively.
For example, the order parameter is much larger than for the fluid with attraction between particles [33], and decreases monotonously with increasing temperature. In this case, in contrast to the one-level system [22], the order parameter turns out to be nonzero even at temperatures above the critical one.
One can also see a sharp increase in the thermodynamic factor at particle concentrations near 0.4 and 0.6. This behavior of the thermodynamic factor is explained by the fact that at a sufficiently low temperature and at concentrations less than 0.4, deeper -sites are predominantly filled. At = 0.4, the -sublattice is almost filled, and the transition of the particle to the -sublattice is energetically unfavorable both because the -sites are shallower and because of the repulsion between the particles.
Thus, at a concentration of = 0.4, the mobility of particles sharply decreases. A similar situation is observed at a concentration of 0.6 with the only difference that, due to the interaction between particles at this concentration, the -sites are predominantly filled, and the sublattice is almost completely empty.
The features in the dependence of the correlation functions on the concentration shown in figure 8 can be explained in a similar way. The sharp increase in the correlation function in the concentration range of 0.4 and 0.6 for temperatures close to the critical one is due to an equally sharp decrease in the concentration of particles on one of the sublattices.

Conclusion
Comparison of the results of the developed approximate methods and simulation data shows that these results are in full qualitative and quantitative agreement in the entire range of temperature and concentration variation. This makes it possible to use, for example, the diagram approximation to study the thermodynamic and structural characteristics of the lattice fluid on a two-level lattice, since DA does not require powerful computers as compared with MC method.
The proposed quasi-chemical approximation refinement method is not the only possible one. Another strategy can be implemented within the framework of the so-called self-consistent diagram approximation. The main idea of the latter is to consider a reference system in which the range of the average potentials exceeds the range of the interaction potentials [29]. However, this approach leads to a significant complication of the obtained relations, which complicates their subsequent analysis.
It can also be noted that in the case of the system with repulsive interaction, the transition from a one-level to a two-level system changes the kind of phase transition in it.
To explain this fact, the isotherm of the chemical potential at low temperature (see curve 1 in figure 3) must be considered together with the corresponding dependence of the order parameter versus concentration. The number of particles on the lattice increases monotonously with an increase in the chemical potential at low temperatures and low concentrations. In this case, as follows from figure 6 (see curve 1), deep -sites are predominantly filled. Obviously, such distribution of particles on the lattice is possible up to the concentration of 0.4. When this limit is reached, all -sites will be occupied and additional particles can be placed on shallow -sites. In this case, each particle on -sublattice will receive 2 nearest neighbors. Therefore, due to the interaction between particles, the concentration of particles increases with a jump to 0.6 when the chemical potential reaches a critical value. This corresponds to the situation when all the -sites are occupied, and the -sites are vacant. There is a plateau on the isotherm of the chemical potential which corresponds to the phase transitions of the first order.
Thus, the phase transition itself corresponds to the transition between two different ordered states of the system. In the first of them, the deep -sublattice is predominantly occupied, and in the second, the shallow -sublattice is predominantly occupied. This conclusion is fully confirmed by the analysis of snap-shots of the simulation procedure. An increase in temperature reduces the ordering of the system. However, this ordering is essential at a sufficiently high temperature 2 (see curve 5 in figure 6) and does not disappear even at extremely high temperatures (10 3 ).