Kinetic Equations of Physicochemical Processes with Allowance for Multi-Particle Effects in the Lattice Gas Model

A way of deriving kinetic equations of physicochemical processes in dense phases is developed on the basis of the discrete–continuous description of the spatial distribution of components in the lattice gas model (LGM), with allowance for multi-particle effects. The emergence of multi-particle effects is associated with the simultaneous influence of all neighbors on the rate of the elementary stage with the participation of a given particle. They include multi-particle potentials of interaction, including quantum–chemical energy calculations, the effect the configurations of neighboring molecules have on the internal motion of the central particle, and the effects of the indirect correlation of interacting particles that occurs for any potential of pair interaction, assuming the internal motions of particles do not depend on the local configurations of neighbors. Multi-particle effects take models beyond the quasi-chemical approximation, which reflects direct correlations of interacting particles through pair distribution functions, and require the use of correlation functions for a larger number of particles in describing their kinetics. The rates of elementary one- and two-node stages are calculated within the theory of absolute rates of reactions in non-ideal reaction systems. Ways of calculating approximate rates of the elementary stages of mono- and bimolecular processes are discussed, along with the possibilities of generalizing the derived equations.


INTRODUCTION
The development of a new associative model of the liquid state [1] has revealed a number of questions in the statistical theory of the condensed state that failed earlier to receive the attention they deserve [1]. Intermolecular pair interaction potentials of components are renormalized into multi-particle potentials under the effect of anharmonic vibrations in bound associates. This results in effective multi-particle potentials that depend on the temperature and concentrations of components. The emergence of effective potentials dependent on temperature and density has long been known from the theory of metal melts [2] and the theory of fluids [3]. We are referring here to effective pair potentials obtained by averaging over the spatial distribution of interacting particles without considering the effects of their vibrations. Directly allowing for vibrations produces effective multi-particle potentials. The effect intermolecular interaction has on the frequencies of interparticle vibration changes the average distance between neighboring particles, along with the depth of the potential well and its curvature near the minimum. A change in energy alters the rates of elementary reactions in liquid phases and their concentration dependences. The associative fluid model was formulated on the basis of the discrete-continuous description of the spatial distribution of components in the LGM [4][5][6]. This theory corresponds to the best developed approach to many of the physicochemical processes in condensed phases that occur in three aggregate states, including large homogeneous [7][8][9][10][11][12][13][14][15] and numerous heterogeneous systems [13][14][15][16][17][18][19][20][21][22][23][24]. Apart from situations with spatial phase separation (the formation of heterogeneous systems), these can contain systems with different spatial matrices formed by one component, and between which other components can locate. These systems contain, e.g., solid porous and non-porous adsorbents (membranes) that absorb vapor or fluid molecules; solid absorbents that absorb atoms, ions, and molecules from other phases into their interstitial sites; and the dissolution of diverse mixtures in complex dispersion phases.
The modern kinetic theory of condensed phases is based on that of absolute reaction rates of non-ideal systems [6,25]. Strictly speaking, multi-particle effects arise automatically in kinetics due to the simultaneous effect all neighboring molecules have on the rate of the elementary stage with the participation of the relevant central particle. The central particle can

CHEMICAL KINETICS AND CATALYSIS
be single when the stage proceeds at one site (monomolecular stages) or two particles when it proceeds in two neighboring sites (bimolecular stages). Due to the simultaneous effect of all neighbors of the central reagents, the rate of the elementary stage depends on the state of the entire multi-particle configuration formed by the central particle and its neighbors. The interaction between them form a reaction cluster and the local activation energy of the stage. The rate depends on the approximation used when considering interparticle interaction. The theory most widely accepted today uses the so-called quasi-chemical approximation (QCA) with allowance for direct correlations between interacting particles [6,19]. This approximation provides a self-consistent description of equilibrium and kinetics when considering the pair interaction potential between components, so long as the internal motions of molecules in LGM cells are independent of the configurations of neighbors.
The aim of this work is to extend this theory to considering multi-particle effects for a discrete description of the spatial distribution of components in the LGM. It should be noted that apart from the above associative fluid model, multi-particle effects also emerge when different types of multi-particle interaction potentials are used, including quantum-chemical calculations of a system's energy (with optimization of the geometry of complexes [26,27]), and in the familiar means of the liquid state theory (integral equations) that describe indirect correlations. Integral equations underlie the current theory of the liquid state described by correlation functions (CFs) with a wide spectrum of interaction potentials [3,4,28,29]. The same CFs are used in discrete interpretations of the LGM [30][31][32]. The need for more detailed consideration of contributions from different molecular configurations, rather than simply allowing for pair contributions from interacting particles, is common to the above situations.
The procedure for deriving kinetic equations for the above mono-and bimolecular elementary stages in condensed phases was developed in [6,25,[33][34][35][36][37]. In this work, we go beyond the pair description of correlation effects, which alters the structure of kinetic equations and ways of calculating approximate rates of elementary stages. Multi-particle effects take models beyond the use of the QCA by employing pair CFs and requiring them for a larger number of particles in describing kinetics (below, we use the term correlator in addition to the abbreviation CF). It should also be noted that initial logic of using CFs to describe the full set of local particle configurations remains the same [6,19,[33][34][35][36][37]. The change in LGM parameters when considering multi-particle effects uses the same set of configurations, so the way in which relationships between CFs are constructed is not altered.
Below, we formulate a model that describes interactions with nearest neighbors in the cluster approach (CA) [6,19,33,[37][38][39]. The CA introduces correlators of the first and second types. The first type of correlators has a simple relationship with the local energy of the system, and all thermodynamic functions are expressed through correlators of the second type. Kinetic equations are deduced for correlators of the second type, and the rates of elementary stages are expressed via correlators of the first type. The structure of kinetic equations is derived below without considering the multi-particle parameters of the LGM. Specific features of deriving new kinetic equations are discussed. General expressions for the rates of the elementary steps of mono-and bimolecular stages are then considered, along with proof of the self-consistency of the derived expressions for forward and backward one-and two-node stages with an equilibrium description of the reaction systems. We then discuss ways of approximating expressions for the rates of elementary stages under the condition of local equilibrium, along with the relationship between calculations of the derived rates of elementary stages and thermodynamic functions of the system beyond the QCA.

MODEL AND CORRELATION FUNCTIONS
In the LGM, system volume V is divided into cells about the size of a molecule ( ). Each cell is numbered with index f, where 1 ≤ f ≤ М; and М = V/ is the number of cells in the system. The occupancy of each site is characterized by value . If site f is occupied by a molecule of type i, = 1. If the site contains another type of molecule, = 0. The number of components in the system is s, and 1 ≤ i ≤ s. The values satisfy conditions (1) and (2) , where is the Kronecker delta. These two conditions mean that each cell (or site) is occupied by a molecule of only one kind.
Let the lattice structure have z nearest neighbors. For simplicity in considering the interaction between molecules, we limit ourselves to the distance of the nearest neighbors. We must consider the multi-particle parameters of the LGM, which include the occurrence of multi-particle interaction potentials and the effect different configurations of neighboring molecules have on the internal motions of molecules.
The total energy of the system (Hamiltonian Н) is written [6] as the sum over all sites of the system: , where is the Hamiltonian for cluster K 1 consisting of one central site f and its z nearest neighbors g 1 …g z , with summands of two types of contributions from the internal motion of molecule i RUSSIAN  (1) where symbol for the parameter of multi-particle interaction retains the meaning of the bond between central particle i in cell f and j k in cell g k , with 1 ≤ k ≤ z. However, the energy of the pair bond is modified by the effect of the other neighbors from g k = 1 to , excluding given site g k . The sum in (1) is taken over all neighbors g 1 …g z containing particles j 1 …j z , respectively. Here, we give abbreviated relation . In other words, the states of occupancy for neighboring particles j 1 …j z change on a cluster with fixed central molecule i in all neighboring sites g 1 …g z . The total number of K 1 cluster configurations is . We use the grand canonical ensemble to describe equilibrium. Chemical potential μ i of molecule i is included on the value: where is the potential of a lattice system at site f, and is the statistical sum of molecule i in a dense phase. For the chemical potential in the gas phase, , where is the reference point of the former: = 0 for a free particle in the gas phase or = ε i in the ad-or absorbed state, and β = (k B T) −1 . Here, P i is the partial pressure of molecule I. The value reflects the ratio of the statistical sum of molecule i in cell f of the lattice system, and in the gas phase: = / is the local constant of retention (for problems of adsorption and absorption, this is the Henry constant at an additional non-zero potential of an ad-or absorbent) of particle i in cell f; and is the statistical sum of a particle in gas. A similar relation of internal motions is given for the motion of component i in a dense phase: In the CA used to derive equations for the equilibrium distribution of components, we use two types of correlators [6,38] to avoid finding the statistical weights of individual configurations in the statistical sum of system Q. Correlators of the first type distinguish local energy characteristics of configurations, while those of the second type are needed to calculate the thermodynamic characteristics of the system. The first correlator of the second type is thus related to local partial occupancies θ I of system sites and/or chemical potential (isotherm) μ i . The second correlators represent probabilities θ ij of the presence of particles i and j in neighboring lattice sites, through which we can determine, e.g., heats of sorption or mixing, heat capacities, and activity coefficients. Introducing two types of correlators allows us to reduce deriving a set of equations for correlators of the second type to an operation using the ratios of their probability (instead of determining the statistical weights of configurations).
Correlators of the first type are defined through cluster Hamiltonians (1) as (3) Exponent (3) contains a cluster local internal energy for a particular configuration ij 1 …j z of cluster K 1 . The ratio of correlators of the first type with different central sites allows us to exclude unknown constant and the statistical sum of system Q, which gives (4) Formula (4) determines the pair relationship between correlators of the first type with identical states of all neighbors j 1 …j z through local energies of a cluster with different occupancy states of central sites А and В.
Correlators of the second type are obtained from correlators of the first type by averaging over the occupancy states of some of the nearest sites. They are defined as (5) The relation for correlators of the first and second types is written as (6) On the right, summation is done over the occupancy states of …j z sites that are missing from subscripts of the correlator on the left. To derive equations for the equilibrium distribution of components we must associate correlators of the second type derived for different central particles [6,38]. The need to move to relationships between correlators of the second type with different central sites is due to all approximations being constructed through correlators of the second type. By substituting Eq. (4) in the right side of (6) for central component i = A expressed through central component B, we have the relation (7) in which correlators of the second type with center А are expressed through correlators of the first type with center В.
Subsequent derivation of equations for the equilibrium distribution of components requires substituting correlator of the first type in (7) for correlators of the second type with the same central particle В ( ). After finding relation , they are substituted into the right side of the set of Eqs. (7), producing the sought expression for correlators of the second type with different central components: On the left, there is a correlator of the second type with center А; on the right, there is a correlator of the second type with center В.
To find relation , we must move from the coordinate description of the distribution of neighbors in different configurations (5) and (6) to their energy contributions from particles of different kinds. If we use the description of configurations of j 1 …j z neighbors through numbers of bonds n ij between center i and a neighbor of the j type, we can rewrite formula (6) for particular configurations by grouping correlators of the first type with fixed numbers n ij (configurations with at least n ij neighbors of the j type in this case correspond to correlators of the second type) [6,38]: (8) where the upper limits of the sums are defined as , .
Here, are the numbers of transitions of correlators of the first type needed to obtain a correlator of the second type with dimensionality n. Symbol σ({n}) denotes the arrangement the number of particles {n} = (the braces signify the full set of types), as does symbol σ({n + k}) corresponding to neighboring particles of different kinds around central particle В. Symbol σ({n}) is of a vector nature in a multicomponent system. ...
,..., s n n Inverting the set of Eq. (8) gives expressions for correlators of the first type through correlators of the second type: (9) Derived relation (9) is substituted into the right side of the set of Eqs. (7), producing the sought expression for correlators of the second type with different central components, which describes their equilibrium distribution. This set of equations reflects multi-particle effects and covers all configurations of neighbors. It is therefore cumbersome.
As the simplest example, let us illustrate the use of symbols σ in a binary system (particles А and В) when the arrangement of neighbors is clearly associated with irreducible configurations of neighbors [6,25]. If we omit the subscript of site numbers, and instead denote through n the number of components А among z neighbors around central particle i = A and B, instead of expressions (6) and (8), we have (10) where sum σ(n + k) is taken over all possible configurations at a given number of particles А (n + k) that occur as they are added to configuration σ(n): , is the number of combinations. The right side of Eq. (10) contains correlators with brackets in subscripts, which correspond to correlators of the first type with a fixed number of particles А and their arrangement.
Using the CA [6,39] when considering multi-particle LGM parameters for a binary mixture, we obtain the following expression for the equilibrium distribution of components (only was considered earlier): (11) where in, e.g., the QCA, , where is the conventional probability of one particle А being near another particle А; is the concentration of component А; and is the probability of two particles А being in neighboring sites.
Multi-particle effects generally require the use of numerical calculations in which deriving the coefficients of is a not complicated procedure.
Correlators of the second type with different dimensionalities are related to one another by normalization conditions that are satisfied for any characteristic time scales τ: (12) Correlators (12) characterize the probability of a cluster consisting of m particles, in which a particle of type i n is at a site with number f n at time τ; 1 ≤ n ≤ m, where т is the dimensionality of the correlator.

KINETIC EQUATIONS FOR CORRELATORS
The evolution of the system is described by kinetic equations for correlators of the second type with different dimensionality. The approach to deriving kinetic equations in the LGM was developed in [6,37,40]. Let the considered process consist of many stages, and we denote by α the number of the stage of the elementary process.
Let us denote through the total set (or the full list) of γ f i angle values of all sites of lattice М, 1 ≤ i n ≤ s, which unambiguously define the total configuration of particle arrangements on the lattice at time τ. Kinetic equations for correlators of the system upon moving from initial state {I} to final state {II} are written as (13) where particle types in the first summand on the right side correspond to the occupancy states of lattice sites in final state {II}. Here, W α ({I} → {II}) is the probability of elementary process α occurring, with the system moving from initial state {I} to final state {II} by time τ as a result. In formula (13), the sum is taken over different kinds of direct processes (index α) and all inverse processes {II} in which the occupancy state of each site of the system changes. Probabilities is the total energy of a lattice system in state {I}. Equation (13) and probabilities of transitions W α are traditionally written in the Mark approximation, where it is assumed that relaxation of the internal degrees of freedom of all particles proceeds more rapidly than changes in the occupancy states of different sites of a lattice system. This condition is modified in the associative fluid model so that the internal states of molecules with different numbers Equations (13) are the start for obtaining kinetic equations of processes that occur in condensed phases. In the right side of the kinetic equation for the local concentration of molecules (where are first order correlators (12) and m = 1) has summands for the rates of elementary one-site processes i ↔ b (where h ∈ z f ), and for the rates of elementary two-site processes i + j α ↔ b + d α (h ∈ z) in the nearest f and g sites [6,29,32]: (14) These summands contain higher order CFs than the first, and Eq. (14) cannot be closed without losing the correlation effects.
The next equation of the system for pair correlators ( ) = ) is written as (15) where the second summand in describes stage i + m ↔ b + c in neighboring f and h sites. Equations (15) contain summands (α) and (α) corresponding to one-(i ↔ b) and two-site (i + m ↔ b + c) reactions of particle i with neighboring particle j. However, particle j does not participate in elementary process α itself; rather, it changes the activation barrier for reacting particle i in the one-site stage and reacting particles i and m in the two-site stage.
As an example, we present the explicit form of kinetic equations for the first two correlators of a binary system in which there is the first stage of nondissociative adsorption-desorption (α = 1), and jumps of particle А to neighboring vacancies V (α = 2) The first brackets correspond to the one-site stage of adsorption-desorption on the central f site (sum- describes the desorption stage), and the second brackets (symbol (2) corresponds to jumps, j = A, V; h ∈ z f ) describes the stage of the А particle jump from the central f site to site h ( ) and vice versa. The value is found from normalization condition (12).
The equation of a system for pair CF is written as (17) Reacting particle А or V is identified in the superscript in parentheses. The first two brackets correspond to the one-site stage of adsorption-desorption in central f and g sites. The second particle of the pair (the superscript without parentheses) affects the course of the process, but the occupancy state of the site with this particle does not change. Third summand corresponds to jumps of particle А from site f to neighboring vacant sites h, ; symbol means there is no site g with neighboring particle А among the nearest neighbors of site f.
. The fourth summand has a similar form: particle А at site g exchanges h with its neighbors, , around site g, except for site f: .
Pair CFs have normalization relations for the occupancy states of both first site f and second site g: , . Note that to describe the evolution of pair CFs it is not enough to know one pair function if (e.g., there is the density gradient along the line/axis passing through fg). We also need the second equation for pair CFs in, e.g., the form (18) The first two summands have the same meaning as in Eq. (17): they are elementary stages at sites f and g, and they have a fixed state of particle А at site f and vacancy V at site g. The third summand describes the particle exchange at the central pair of sites fg (superscripts have no parentheses, since both particles participate in the stage). The last two summands have the structure of expressions and : (1) To consider indirect correlations along with multi-particle effects, we must to move to the equations of system (13) with large correlators. In Eqs. (14) and (15), the dimensionality of correlators on the right is larger than that of correlators on the left [6,25,37,40]. This dimensionality equals the number of all neighboring molecules that simultaneously affect the rate of this stage. The minimum dimensionality of the correlator in (14) is the dimensionality of reaction cluster K 1 = 1+ z for onesite stages (for two-site stages, the dimension of the reaction cluster is around K 2 = 2z). We recall that in Bogolyubov's kinetic hierarchy [41], the dimensionality of correlators in the right side of kinetic equations grows in steps of 1. This differentiates the structure of LGM kinetic equations from the traditional form of kinetic equations of Bogolyubov's hierarchies [41][42][43][44], and the problem of the self-consistent calculation of high dimensionality correlators up to the K 1 dimension arises in the kinetic theory. Aspects of the procedure for deriving kinetic equations for highest correlators is can be discussed using the example of a triple correlator in a binary system.

Triple correlator
is the first correlator that describes indirect correlations. As above, we let the kinetic equation for its evolution also be formed in two types of stages: adsorption-desorption (α = 1) and jumps to the neighboring sites (α = 2) Here, first summand reflects the contributions from one-site adsorption and desorption process at sites fgh (three parentheses describe the participation of particles А in different sites fgh in this triple correlator) (20) Second summand in (19)   , and the parentheses in the superscripts indicate at which sites fξ jumps take place. The potential of other particles А at sites gh affect the rate of the elementary jump stage between sites fξ, but they remain in their sites. Similarly, and .
The triple correlator has four conditions of normalization (12), , which associate them with four pair CFs, , This means that in addition to the equation for triple correlator (19), there must be other kinetic equations containing one and two vacancies in the superscripts: (21) and The structure of these summands is written as when moving from (17) to (18). In addition to jumps from neighboring sites ξ, there are summands with direсt jumps between the central sites for each of three types of correlators at sites fgh. For the and , summands are written as (22) Here, one-site stages occur in three sites fgh of the correlator. The character of the arrangement of sites fgh is important for contributions from the summands with jumps. For instance, the three neighboring sites in square lattice (z = 4, d = 2) can be either on one line or form an angle, so sites f and h cannot be nearest sites. We therefore have one jump between sites g and h for z = 4: while for planar hexagonal lattice (z = 6, d = 2) three sites can be nearest ones. This results in the second summand in : In both cases, the three last summands in the expressions describe the exchange between particles А (2) and neighboring vacancies (i.e., they provide an expansion of the dimensionality of correlators, along with the effect of the lateral interaction potential): The same procedure is used in the successive increase in correlator dimensionality n at sites (f 1 f 2 … f n ). Below, we present the structure of kinetic equations in the general form only for particles А and the vacancies among them (occupancy states of correlators of dimensionality n at sites f 1 f 2 … f n are denoted as The presented structure of set of kinetic equations (25) assumes the additional use of n(s − 1) normalizations (12) for correlators with dimensionality n and different numbers of vacancies in the superscripts. The system contains contributions from one-site stages and particle jumps to vacant sites , which grow along with n. These jumps occur at both sites of this correlator with dimensionality n (summand ) and sites of the correlator with dimensionality (n + 1) (summand ). This set of equations for adsorption of a non-dissociating molecule can be easily rewritten for dissociative adsorption-desorption when the elementary stage proceeds at two neighboring sites. All summands from stage α = 1 at one site in this case disappear. In their place emerge summands from stage α = 1 at two sites. The principle of deriving the kinetic equations is retained for multicomponent mixtures and the growing number of elementary stages α, but the system becomes much more cumbersome.
To close the equations for correlators we need the explicit form of expressions for the rates of elementary one-and two-site stages that occur in reaction clusters with sizes K 1 = 1+ z and K 2 = 2z. All rates of the elementary stages in dense phases and are calculated using the theory of the absolute rates of reactions in non-ideal reaction systems [6,25,[33][34][35][36].

RATES OF ONE-SITE STAGES
In the transitional state, an activated complex (AC) interacts with neighboring particles. The multi-particle parameter of interaction between the AC and neighbors ( ) formed by particle А  ij j j ij fg fg g g K with neighboring particles j 1 , … j z at sites g 1 ,…,g z differing from parameter of interaction between initial reagent i from which the AC is formed and neighboring particle j. When there is no AC of one-site reaction А → В, the total energy of the system is described by Eq. (1). The formation of the AC changes the local energy of the system by a value equal to a change in the energy of the cluster Hamiltonian (26) where The probability of AC reaction А → В occurring at site f of reaction cluster K 1 is that of a complex event: particle А occurring at site f upon a change in the system's energy by the value of . It is expressed as (27) The exponential dependence in (27) is due to the assumption that AC has chemical potential and is in equilibrium with a reaction subsystem consisting of initial reagent A and the fixed arrangement of all neighboring particles around A, interaction with which determines the activation barrier height. We are talking here about the local equilibrium (reagent-AC), though the distribution of particles over lattice sites and the equilibrium distribution can differ as much as possible. As shown above (see the text between formulas (1) and (2)), the chemical potential is associated with internal motions of AC in through AC statistical sum . By distinguishing the translational degree of freedom along the reaction coordinate from the AC statistical sum and expressing the rate of the motion along it using characteristics of the activation barrier , where all values have their usual meaning [6,33,[45][46][47] ( is the transmission coefficient; is the average rate of the AC's motion along the vertex of the activation barrier with length ; is the concentration of AC at the given stage), we obtain an expression for the rate of an elementary process at one site of the system with correlators of the first type (3): (28) where index k at neighboring site g k is associated with the same rearrangement of correlators into a sum over  (8) and (9). The rate constant of elementary stages at a fixed distribution of neighbors K 1 , which determines the local activation energy through the energies of the AC and the initial reagent (h is the Planck constant) was introduced in (28): (29) In the right side of (29) there is a difference between the chemical potentials of the AC and the initial reagent to reflect within the LGM the occurrence of a particle from the gas phase not contained in the LGM components. If the AC is formed by components of the lattice system, and Eq. (29) have a common meaning: . If a gas molecule participates in the formation of the AC (e.g., vacancy V and gas molecule A participate at the stage of adsorption), , where is the chemical potential of molecules А in the gas phase, , and its partial pressure appears in formula (29) [6,19,25]: .
The multi-particle nature of statistical sums and in (29) does not allow us to transform Eq. (28) into the traditional rate constant of the elementary stage in the form , where is the rate constant of the one-site stage in an ideal system and is the concentration component of the rate that considers a change in the activation energy because of the non-ideality of the reaction system (as for the pair potential [6,25,33]), though when there is no interaction. With a weak dependence, however, the ratios of statistical sums in (29) from local configurations in the reaction cluster (as when considering indirect correlations at the potential of pair interaction), the structure of (29) allows us to move to the form in (28).
For a binary system of expressions (28) and (29) and using (n, σ) to determine the configurations, the rates of the stages are written through correlators of the first type, as above: (30) Parameters of interaction on a homogeneous lattice do not depend on the site numbers and can be character- ized by the number of particles A in the coordination sphere and their configurations σ(n). In the right side of kinetic equations for all correlators n > 1, summands of one-site stages of the (15) type for n = 2 or for any n > 2 (25) appear in addition to summands . The meaning of these summands is that neighboring particle j at site g is in the immediate environment of reacting particle i, and it changes the activation barrier height for the reaction of particle i. The particle itself does not participate in the reaction, and the occupancy state of site g does not change. Summands are in expression for one-site stages that can occur in all n sites of a correlator of this dimensionality. All neighboring sites that are nearest to the reagent and which retain their occupancy states can contain particles affecting the activation barrier height. In these situations, Eq. (28) for the rates retains the correlator dimensionality of n, but a contribution from energy appears in (28), depending on the approximation used to consider correlations between interacting particles: where the prime mark before the sum over all states of occupancy of sites g 1 … g z means that averaging is not performed over one of the neighboring g sites. Particle j makes a constant contribution to the rate constant through . For example [6], we introduce Eq. (28a) for the pair potential in the QCA when n = 2, s = 2, have the form = , = exp(βδ )/ , where = , , , where is the parameter of pair interaction between the AC and initial neighbors. With LGM multi-particle parameters, Eq. (28a) depends on the approximation that is used and correlator dimensionality n.
The derived expressions for the elementary onesite reaction rates must have the property of a self-consistent description of the equilibrium state in the kinetics path to equilibrium. Ideal reaction systems always have this property, since the theory is of a oneparticle nature. Problems arise when considering interparticle interactions responsible for the non-ideal behavior of real systems.
The essence of self-consistency is that by equalizing the rates of forward and backward reactions, we obtain the same equation for the equilibrium distribution of reaction mixture components, which immediately follows from the derivation of the equation in the equilibrium theory. From a formal viewpoint, this  means that by equalizing the rates of opposite elementary stages containing AC parameters (the AC's internal statistical sum and energies of interaction between molecules of the system and the AC), we obtain an expression that depends only on the molecular parameters of system components, rather than these AC parameters.
To prove this property, let us consider one-site stage А ↔ В and write Eq. (4) for correlators of the first type with all indices and central components А and В in the form (31) Note that AC properties at the top of the barrier do not depend on the direction of the reaction (i.e., the statistical sum and the AC parameters of potential interaction formed from reagents А and В are the same: and . Let us divide Eq. (31) on the left and right, respectively, by the same value of the exponent with indices and . On each side we then have an expression for the local rate of a one-site stage for a fixed configuration of neighbors K 1 (32) After the identical averaging of Eq. (32) on the right and on the left over all configurations of neighbors {j} on cluster K 1 = , it can be rewritten in the form of (28) with regard to (29). Otherwise we obtain the equality of rates of stage А ↔ В in opposite directions .
RATES OF TWO-SITE STAGES Let us consider two-site stages of the form А + С → В + D at two central sites f and h [6,25,33,35]. The AC of the two-site stage occupies two sites f and h. Neighbors of reaction cluster K 2 are all neighbors g located at the shortest distance from one of the sites of the central pair of sites f and h with reagents А + С. Dimensionality K 2 is 2(z − 1) + 2 = 2z, if neighboring (z − 1) sites for central sites f and h do not intersect. If they do intersect, dimensionality K 2 is less than 2z. (For a 2D lattice where z = 6, dimensionality K 2 is 2(z − 1).) To simplify denotation, below we consider dimensionality K 2 to be 2z for all lattice structures and refer to the total number of sites of the nearest neighbors g as 2z (1 ≤ g ≤ 2z), though it can be smaller (e.g., 2(z − 1)).
In analogy with above denotaion for onesite clusters K 1 , we accept for the states of all neighbors of cluster K 2 and introduce a system of correlators of the first type, constructed at two central sites f and h: (33) The meaning of all values introduced is similar to that of the above parameters for one-site clusters K 1 with initial Hamiltonian (1). For simplicity, we retain the form of the denotation of LGM multi-particle parameters that describe internal motions and interparticle interactions . Neighboring sites g, 1 ≤ g ≤ 2z, are located asymmetrically around f and h, and two distances are required between neighbor g and two central sites f and h [6,25,33]. One index k, 1 ≤ k ≤ 2z, is used to simplify the denotations in (33). Here, , where the last summand means that energy of bond fh is considered twice in both and . The value of corresponds to bond fh, so it depends on the model of multi-particle interactions. The simplest version of symmetrization is , which is an approximation of the quantum chemical calculations of bonds in [26,27].
The ratio of correlators of the first type and different pairs of central sites allows us to exclude the unknown constant and statistical sum of system Q, yielding an expression with a structure similar to that of Eq. (4) for clusters K 1 : (34) Correlators of the second type are defined as (35)  is the concentration of the AC of the considered two-site stage, we obtain an expression for the elementary rate of the process through correlators of the first type per system site: (39) where the value that appears when the local contribution from the AC is considered is the energy of reagents А and C before their transformation into the AC (the fh bonding energy is internal for the AC and is not written in explicit form. As above, the symbol of index k at neighboring site g k in the exponent is due to the rearrangement of correlators into the sum over different occupancy states of all neighbors. Here, the rate constant of an elementary stage is also introduced at the fixed distribution of neighbors K 2 defining the local activation energy through the energies of the AC and initial reagents (where h is the Planck constant): (40) As in (29), the difference of chemical potentials of AC and the initial reagent is introduced into the right side of (40) to reflect within the LGM the presence of a particle from the gas phase that is not one of its components. If the AC is formed from lattice system components, and Eq. (40) have a common meaning: . If a gas molecule participates in the formation of the AC (e.g., two neighboring vacancies V and gas molecule А 2 participate in the stage of dissociative adsorption, , where is the chemical potential of molecules А 2 in the gas phase), then , , and partial pressure appears in formula (40) [6,25,27]: .
In the right sides of kinetic equations for all correlators n > 1, summands appear for two-site stages of type (15) when n = 2 and for when n > 2 (25), in addition to . As above, particle j in the immediate environment of reacting particle i does not participate in the reaction, but its presence changes the activation barrier height for the reaction of parti-  where the prime mark before the sum over all occupancy states of sites g 1 … g 2z in (39a) means that averaging does not include one of neighboring sites g, 1 ≤ g ≤ 2z. Particle j makes a constant contribution to the rate constant (40) through . For example [6], let us consider Eq. (39a) for the pair potential in QCA at n = 2, s = 2. Function has the form , = , where all values are defined above. With LGM multi-particle parameters, Eq. (39a) depends on the approximation that is used and correlator dimensionality n.
For correlators with two central sites we can prove the self-consistency of the description of forward and backward stage rates and local equilibrium. Let us consider two-site stage А + С ↔ В + D, and for central components (А + C) and (В + D) we write an expression for correlators of the first type (34) in the equivalent form (41) Note that AC properties at the top of the barrier do not depend on the reaction direction (i.e., statistical sums and potential interaction parameters of the AC formed from reagents А and В are the same: and . We divide equality (41) on the left and right, respectively, by the same value of the exponent with indices and . On each side we then have an expression for the local one-site stage rate at a fixed configuration of neighbors K 2 : Equation (42) corresponds to reaction cluster K 2 as a whole, and we have neither the isolation of contributions from clusters K 1 nor the values of and . By doing the same averaging over all configurations of neighbors {j} on cluster K 2 = on the right and left sides of Eq. (42), we can rewrite (42) in the form of (39). Otherwise, we obtain the equality of rates of stage А + С ↔ В + D in opposite directions: . In this section, we do not present an analog of Eq. (9) for correlators of the second type with two central sites for multi-particle potentials because it is cumbersome. In [35], we gave an example of applying this equation and deriving functions of indirect correlations for dissociative adsorption for the potential of pair interaction. We also proved the compatibility of equilibrium systems of equations for correlators of the second type constructed at one and two central sites.

CLOSURE OF KINETIC EQUATIONS AND RATES OF ELEMENTARY STAGES
The derived structure of the system of kinetic equations for correlators from unary (n = 1) to size n contains contributions (28) from reaction correlators of dimensionalities K 1 and K 2 on the right for one-site stages and (39) with correlators K 2 and larger due to two-site stages. To close a system of kinetic equations we must express and , where is the set of independent correlators of the second type that describe molecular distributions for the selected level of the accuracy of approximation. This set of correlators is identical for contributions (28) and (39).
The dimensionality of reaction correlators K 1 and K 2 always exceeds that of unknown coupled correlators n on the left, and these equations are derived for their evolution. Kinetic equations for unary (concentrations) and pair DFs are usually considered. For the potential of pair interaction, the QCA corresponds to adequate closure of the reaction correlators [6,25]. This way of deriving kinetic equations permits generalization from a rigid lattice structure to a soft (deformed) lattice structure [1,48]. This allows correct calculations of pressure with regard to its experimentally measured periods of relaxation with respect to that of the chemical potential. The concept of F mod allows us to consider coupled correlators of any dimensionalities from n = 2, 3, … to K 1 and further. The use of F mod combines the introduced model representations of multi-particle internal motions in the LGM [1] in both rigid and soft lattices.
The ideas of the cluster approach do not depend on the type of the potential or the statistical sum of internal motions of system components. They operate with a full set of system configurations and can be applied to its non-equilibrium states. The self-consistent equations for the equilibrium and non-equilibrium spatial molecular distributions derived in the LGM with regard to the potential of multi-particle interaction and the dependence of internal molecular motions on the configurations of neighboring molecules (and their composition and structure in bound associates [1]) greatly expand the possibilities of calculating the kinetics of physicochemical processes. These equations can be applied not only to vapor-liquid phases (as in [1]) but to all heterogeneous threeaggregate systems as well, including adsorption systems with contributions from the anharmonism of molecular vibrations [20][21][22][23][24][57][58][59][60].

CONFLICT OF INTEREST
The author declares he has no conflicts of interest.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.