1. Introduction
Understanding the statistical mechanics of structured particles with arbitrary size and shape in external fields remains a major theoretical challenge, largely due to the complex entropic contributions arising from particle configurations at finite densities. Even for simplified models, such as linear particles with hard-core interactions on regular lattices, the problem is analytically intractable. This difficulty stems from spatial correlations among allowed particle configurations, which complicate the calculation of thermodynamic potentials. These correlations underlie various emergent collective behaviors, including nematic ordering in systems of linear
k-mers [
1], and entropy-driven competition in multicomponent mixtures. Exact solutions have been found only in a few special cases, such as dimers on square lattices [
2] and hexagons on regular lattices [
3].
In continuum systems, this problem has been extensively studied. In three-dimensional colloidal suspensions, Onsager famously demonstrated that elongated molecules undergo a phase transition from an isotropic to a nematic phase [
1]. In two dimensions, although continuous rotational symmetry cannot be spontaneously broken, a Kosterlitz–Thouless transition occurs, characterized by a power-law decay in orientational correlations [
4,
5].
In contrast, the case of hard-core particles on lattices is less well understood. Early work by Flory [
6] and Huggins [
7] initiated the study of rigid rods, or
k-mers, modeled as linear arrangements of
k identical units occupying contiguous lattice sites. These rods interact solely through hard-core exclusion, meaning that no site may be occupied by more than one unit.
The Flory–Huggins (
) theory, developed independently by Flory [
6] and Huggins [
7], generalizes the theory of binary liquid mixtures or dilute polymer solutions on lattices. In the lattice–gas framework, the adsorption of
k-mers on homogeneous surfaces is formally analogous to polymer–solvent binary solutions.
Extensive efforts have been made to assess
theory against experimental results, the theory being completely satisfactory in a qualitative, or semi-quantitative way. There is no doubt that this simple theory contains the essential features which distinguish high polymer solutions from ordinary solutions of small molecules. Modified forms of the
approximation have been also proposed. A comprehensive discussion on this subject is included in the book by Des Cloizeaux and Jannink [
8].
The
statistics, given for the packing of molecules of arbitrary shape but isotropic distribution, provides a natural foundation onto which the effect of the orientation of the ad-molecules can be added. Following this line of thought, DiMarzio [
9] developed an approximate method of counting the number of ways,
, to pack together linear polymer molecules of arbitrary shape and of arbitrary orientations. Accordingly,
was evaluated as a function of the number of molecules in each permitted direction. These permitted directions can be continuous so that
is derived as a function of the continuous function
which gives the density of rods lying in the solid angle
, or the permitted directions can be discrete so that
is the the number of ways to pack molecules onto a lattice. Based on the detailed knowledge of the orientations of the molecules, the various types (nematic, smetic, and cholestic) of liquid crystals were argued for and the reasons for their existence were ascertained. In the case of allowing only those orientations for which the molecules fit exactly onto the lattice is that for the case of an isotropic distribution the value of
reduces to the earlier result by Guggenheim [
10], now known as the Guggenheim–DiMarzio (
) approximation.
In the 2000s, two novel approaches were proposed for describing multisite adsorption. The first, developed by Ramirez-Pastor et al. [
11], introduced the Extension Ansatz (
) model for linear adsorbates on homogeneous surfaces, based on exact one-dimensional thermodynamic expressions and their generalization to higher dimensions. The second, the Fractional Statistical Theory of Adsorption (
) [
12,
13], incorporates the internal configuration of the adsorbed molecule as a model parameter.
generalizes Haldane’s fractional exclusion statistics [
14,
15], originally developed for quantum systems, to describe classical polyatomic adsorption at gas–solid interfaces.
Comparisons with simulation data [
11] have shown that the
approximation agrees well at low surface coverage, while the
model performs better at high coverage. These insights led to the development of the Semiempirical (
) Model for Polyatomic Adsorption [
11,
16], a hybrid model combining exact 1D results with
approximations, weighted appropriately.
More recently, the Multiple Exclusion (
) statistics framework was introduced to describe classical systems in which particles access spatially correlated states [
17,
18].
statistics accounts for situations in which multiple particles simultaneously exclude access to a common state, an intrinsic feature of non-monomeric particles on a lattice. The uncorrelated limit of
statistics recovers both the Haldane-Wu and
formalisms. This approach was further extended in Ref. [
19] to mixtures of particles with arbitrary shapes and sizes, allowing for analytical expressions of thermodynamic quantities in terms of coverage and species densities.
Despite the number of studies dealing with the adsorption of polyatomics on discrete lattices, there are many aspects which are still outstanding. While the problem can be easily and precisely defined, exact solutions for adsorbed correlated particles, such as
k-mers, have historically proved elusive, with results being limited to one-dimensional substrates [
20] and a few shapes in dimensions greater than one. The classic example of such a model is the lattice-gas of dimers (
) [
2,
21,
22,
23,
24,
25,
26,
27,
28]. A review on the entropy of fully packed dimers on planar lattices may be found in Ref. [
29].
The inherent complexity of the
k-mer problem is further increased when attempting to obtain approximate solutions for the thermodynamic functions of systems that, in addition to allowing multiple site occupancy, also involve lateral interactions among adsorbed molecules and/or surface heterogeneity. In this context, simple solvable models of adsorption on homogeneous surfaces serve as valuable foundations for developing alternative approaches to more complex cases involving interacting adsorbates [
30,
31] and heterogeneous surfaces [
32,
33,
34,
35,
36,
37].
In this work, we present a comprehensive overview of foundational and recent theoretical developments in the modeling of structured particle adsorption on regular lattices (commonly referred to as multisite occupancy adsorption). We focus on how particle geometry and size affect the configurational entropy of the adsorbed layer, an aspect that has rarely been systematically treated in thermodynamic models. Understanding entropic effects in polyatomic systems is particularly relevant for applications such as alkane and hydrocarbon adsorption, which are key to petrochemical separation technologies.
The paper is structured as follows:
Section 2 examines the thermodynamics of one-dimensional lattice gases composed of interacting and non-interacting linear particles, covering both single-species and mixture adsorption, including monolayer and multilayer regimes.
Section 3 presents theoretical approximations for non-interacting polyatomic species in two dimensions, including the
and
models, the
extension of one-dimensional results, the
framework based on fractional statistics, the Occupation Balance (
) approximation, and the
model. Adsorption of single and multicomponent species is discussed in both monolayer and multilayer contexts.
Section 4 explores two-dimensional lattice gases of interacting structured species via mean-field and quasi-chemical approaches. Intermolecular interactions give rise to possible phase transitions.
Section 5 introduces the
statistics framework for classical lattice gases of arbitrarily shaped particles, generalizing the formalism of multiple exclusion statistics presented in Ref. [
18].
Section 6 extends
statistics to multicomponent systems, analytically describing the exclusion spectra in terms of lattice coverage and species densities. This appears as a suitable framework to address complex lattice gases mixtures where spatial state correlations are significant to understand their phase behavior.
Section 7 discusses applications of the main theoretical models developed in this review, comparing model predictions with Monte Carlo simulations and experimental data.
Section 8 focuses on computational methods. The statistics of polyatomics is also a very demanding problem from a computational point of view. Whereas for monomer particles (
) the thermal equilibrium is quickly reached using standard adsorption–desorption MC algorithms, the relaxation time for large particles increases very quickly as the density increases. Consequently, MC simulations are very time consuming at high density and produce artefacts related to non-accurate equilibrium states. In order to cope with these difficulties, efficient MC simulations based on cluster moves were developed in the literature. The use of these techniques has made it possible to investigate the behavior of the system at high densities. In
Section 8, the main computational algorithms of interest for the study of adsorption problems involving multiple-site occupancy are presented. Most of these algorithms have been used throughout the present work. Finally,
Section 9 presents our conclusions and future perspectives.
5. Latest Developments, Part I: Multiple Exclusion Statistics for Spatially Correlated Single Species
Recently, it has been proposed that the complex problem of interacting particles with arbitrary size and shape can be addressed using concepts from fractional statistics, extending them to classical systems of particles with spatially correlated states [
12]. This approach reframes the counting of allowed configurations in terms of exclusion principles that go beyond Pauli-type constraints, providing a new route to approximate thermodynamic functions in structured lattice gases.
We briefly revisit the foundational arguments and structure of the formalism for a single particle species in a homogeneous external field. The generalization to mixtures of species will be presented in the following section of this review (
Section 6).
5.1. Multiple Exclusion Statistics Formalism
We consider the equilibrium of identical particles distributed over a set of states within a volume
V. These states are, in general, spatially correlated—meaning that the allowed equilibrium positions (i.e., the state spectrum) are distributed over space with correlation lengths smaller than the size of a particle
8. As a result, a particle occupying a particular state not only excludes that state but also prevents occupation of additional states due to the spatial extent of the particle relative to the distribution of states.
In systems where each particle excludes a constant number
g of states, this process has a classical analogy to the quantum exclusion principle, as introduced in Ref. [
14]. However, due to the spatial arrangement of states (or the lattice topology, in the case of lattice systems), multiple particles may simultaneously exclude the same state. This occurs, for example, in regular lattices with structured particles spanning more than one site. We refer to this statistical phenomenon as multiple exclusion (
) of states, which significantly influences the entropy and thermodynamic behavior of the system.
Since
is inherently configuration-dependent, we construct an approximate expression for the partition function using a state-counting ansatz that captures the configuration-dependent
while allowing for analytical and numerical thermodynamic analysis. In the special case of configuration-independent, constant statistical exclusion,
statistics reduces to the well-known fractional statistics introduced by Wu [
15]. Our objective is to develop a thermodynamic framework to describe lattice gases of structured particles through the lens of state exclusion, and to further extend it to mixtures of hard-core particles of arbitrary shape and size on a lattice.
Consider N identical particles with access to G single-particle states within volume V, at temperature T, and single-particle energy . The canonical partition function is given by , where is the Hamiltonian of the configuration, , is the total number of configurations, and the internal partition function.
Assuming particles exclude states among the
G available, which are not independent (i.e., spatially correlated), the number of configurations is approximated by
, where
is the number of states available to the
particle once the previous
have been added to the system. This expression is exact when the states are independent and the exclusion is constant, as in Ref. [
15]. For correlated states, it is approximate, since
generally depends on configuration.
A statistical ansatz has been formulated to evaluate and, in the thermodynamic limit, the density of accessible states per particle , where is the occupation number.
Applying Stirling’s approximation
, the Helmholtz free energy
leads to the intensive functions per state
and
, with
The chemical potential
takes the form
or, alternatively, as a ratio between occupied and unoccupied states:
where
,
, and
is the fraction of unoccupied states. For simplicity, we assume
.
If the system exchanges particles with a reservoir at and T, the time evolution of the occupation number is , where and are the probabilities of a state being empty or occupied, and , are the transition rates.
At equilibrium,
, which implies
Since
, Equation (
287) yields the expression
which is approximate, as long as
is not exact.
We now define the exclusion spectrum function [
17]
which gives the average number of excluded states per particle at occupation
n9.
Additionally, the exclusion density function is defined through the excluded fraction , such that . Then, measures the number of states excluded by inserting a particle at occupation n.
The normalization defines , where g is the number of states excluded per particle at saturation. Thus, corresponds to the exclusion caused by a single isolated particle, and .
From Equations (
287) and (
290), we obtain:
where
. These spectral functions provide a route to obtain detailed statistical information about the exclusion spectrum from thermodynamic observables such as the adsorption isotherm
, as will be further explored in
Section 5.7. All quantities are expressed in terms of the occupation number
n, which facilitates interpretation in this framework. They can be converted to lattice coverage
by the transformation
.
5.2. States Counting Ansatz: Density of States
The exact configuration counting for a general problem of this nature remains an open challenge in statistical mechanics and constitutes a demanding analytical task. Here, we introduce a self-consistent thermodynamic approximation, based on a new
statistics inspired by fractional statistics ideas [
14,
15], to describe classical systems of a single species with spatially correlated states—such as structured particles on regular lattices. We develop a general approximation for the density of states
that accounts for statistical correlations among states. Applications to the phase behavior of the classical problem of
k-mers on square lattices show that this
formulation predicts simulation results with remarkable accuracy. Additional results for squares and rectangles on square lattices have also been discussed [
18]. The extension to mixtures and more detailed treatment of
k-mer transitions have been recently addressed [
19] which we summarize in the next section.
All thermodynamic and exclusion functions introduced in
Section 5.1 are determined by the density of states. We now present the approximation for
.
The state counting scheme for a single species is summarized as follows. Let G be the total number of states available to a single particle within volume V. As we successively add identical particles from 1 to N, each occupies one state and excludes others. Importantly, the number of states excluded per particle changes with N, due to spatial correlations among the states—resulting in multiple exclusion () as previously discussed.
The following recursion relations define the number of available states for the
particle:
,
, ...,
, where
represents the state occupied by the particle plus
, the number of states excluded uniquely by the
particle [
17]
10.
In the thermodynamic limit
,
, and
, the exclusion
depends only on
n. Along with the recursion, we introduce a counting ansatz to evaluate
, assuming
, where
is a system-dependent parameter called the exclusion correlation parameter
11. Here,
represents the fraction of states still accessible to the
particle, and
resembles a mean-field-like approximation over single-particle states.
Thus, the recursion becomes:
In the limit with , retaining the first term of the sum yields , which describes the fraction of states (relative to total G) accessible to a particle at occupation n.
More generally,
takes the form
, where constants
,
are determined by boundary conditions:
and
at saturation. This gives
and
, with
being the state density at maximum occupation
. Hence,
For systems where the entropy per state vanishes at saturation (i.e., ), such as symmetric k-mers on a 1D lattice, one has . However, in most structured-particle systems on lattices, , and therefore .
For independent particles with uncorrelated states,
and
, and Equation (
294) reduces to
, corresponding to Haldane-Wu fractional statistics for constant state exclusion
g per particle. While Haldane’s generalization of the Pauli principle was originally introduced for quantum particles with
[
14], structured classical particles with excluded-volume interactions can behave as super-fermions with
.
To better illustrate how the
framework generalizes known statistics, from Equations (
289) and (
287), the
distribution function
can be written as
where
and
satisfies
In the limiting cases: for independent states (
,
) and
,
, so
, and Equation (
295) becomes the Fermi-Dirac distribution. For
,
,
, recovering the Bose-Einstein statistics. For constant exclusion
, with
, Equation (
296) reduces to Wu’s equation
, and Equation (
295) yields the Haldane-Wu fractional
g-statistics [
14,
15].
In general, for particles with spatially correlated states (as studied here), one finds
,
,
, and Equation (
296) must be solved for
n at a given chemical potential. In this work, we alternatively compute
as a function of
n directly using Equation (
286).
5.3. Density of States Parameters
In this section, we examine the statistical and physical interpretation of the parameters involved in the density of states function , and discuss how these parameters are practically determined through the relationship between lattice topology and particle structure (size and shape) via the thermodynamic limits of the exclusion spectrum.
In
Section 5.2, the constant
was introduced so that the density of states satisfies the boundary condition
. Physically,
represents the ratio between the number of states available to a particle at saturation and the number of states available to an isolated particle. Generally, for a phase of structured particles on a lattice at saturation, a finite number of configurations per particle is expected, so
. The saturation entropy per state
is related to
via Equations (
285) and (
294), yielding
The parameter
(or alternatively, the saturation entropy
) is the only free parameter of the
statistics needed to describe a wide class of complex lattice gases. However, in the analysis of the
k-mers problem on the square lattice (developed in subsequent sections),
is not treated as a free parameter. Instead, it is fixed for each value of
k using Equation (
297) to match the Monte Carlo values of
reported in [
90,
91].
As an example, for k-mers on a square lattice, assuming that the entropy vanishes at saturation () implies . Within this minimal approximation, the formalism predicts an isotropic-nematic (I-N) transition at intermediate coverage for . However, both the I-N and a high-density nematic-isotropic (N-I) transition arise only for , even for small positive values of . This behavior is latter discussed in the following section.
Next, we determine the exclusion correlation parameter from the lattice and particle characteristics. For a given particle-lattice system, the total number of distinguishable single-particle states G and the number of particles at saturation (i.e., the maximum number of particles that fit without overlap) are first computed. Then, gives the number of states excluded per particle at full coverage. For instance, for rod-like k-mers on a 1D lattice with M sites, , , and . On a 2D square lattice with sites, (two orientations per site), and , leading to .
Since each particle occupies one state, the occupation number (i.e., fraction of occupied states) relates to lattice coverage .
The exclusion correlation parameter
can be determined from the configurational boundary condition at infinite dilution. In this limit,
represents the number of states excluded by an isolated particle, denoted
. From the definitions of
,
,
, and
[Equations (
289)–(
294)], we find that
, which gives
This fundamental
equation relates model parameters to the number of states excluded by an isolated particle. Since
is known for a given species on a lattice, Equation (
298) can be solved to determine
. The solution has an analytical form:
where
is the principal branch of the Lambert function
12 and
. Details of the Lambert function are given in the next section.
Moreover, , so the limits of at and provide full characterization of state exclusion.
In 1D systems of ideal
k-mers, we have
,
, and
, which yield
for any
k, reproducing the exact results of Ref. [
20].
In contrast, for straight rigid
k-mers on square lattices with
and
, it follows that
, so
. The number of excluded states is
for
, while
,
for monomers (
). Here,
accounts for exclusion across the particle, and
for exclusion along it. The solution for
becomes:
with
.
Statistically,
originates from the
state counting ansatz in
Section 5.2. The exponential decay term
in
becomes
in terms of lattice coverage
. The ratio
thus defines the typical coverage at which the
term in the density of states decays. For instance, in the
k-mer model discussed next, the isotropic-nematic transition occurs around the point where
, indicating that most of the single particle states are excluded and so are most of the isotropic-phase configurations. Hence, the ratio
has a clear statistical and physical meaning.
5.4. Lambert Function
The Lambert function
, introduced in 1758 [
92], is defined by the equation
, with notable values
,
, and
as
.
The solution of Equation (
298),
yields:
with
, valid for
,
, and
. The function
denotes the principal (positive) branch. In particular,
for
. For
k-mers in 1D,
,
, which results in
.
For equations of the form
, with
, a substitution
transforms it into
, with solution
, yielding
. Specifically, for
k-mers, setting
,
, and
, we obtain:
valid for
, where
is the principal branch.
5.5. Entropy of k-Mers on Square Lattices: Orientational Phase Transitions
This section examines the phase behavior of straight rigid
k-mers adsorbed on the square lattice. This system was initially studied in Ref. [
96], where Monte Carlo simulations provided strong numerical evidence that nematic order emerges at intermediate densities for
, beyond a critical density
. Additionally, using high-density expansions, Ghosh and Dhar [
96] offered a qualitative description of a second phase transition, from a nematic to a non-nematic state, occurring at a critical density
for large
k.
Building on the seminal work of Ghosh and Dhar [
96], numerous studies have explored the phase transitions in systems of long, straight, rigid rods on two-dimensional lattices with discrete allowed orientations [
97,
98,
99,
100,
101,
102,
103,
104,
105,
106]. These investigations showed that for
, no phase transition occurs. However, for
, increasing the density leads to three distinct phases: a low-density disordered (isotropic) phase, an intermediate-density nematic phase, and a high-density disordered phase with no orientational order. The threshold value
depends on the lattice geometry:
for square [
96,
97] and triangular [
98] lattices, and
for honeycomb lattices [
99]. The intermediate-density nematic phase, characterized by a large domain of parallel
k-mers, is separated from the low-density isotropic state by a continuous phase transition at a finite critical density
. This first transition, commonly referred to as the isotropic–nematic (I–N) phase transition, belongs to the two-dimensional Ising universality class for square lattices [
97], and to the three-state Potts universality class for triangular [
97] and honeycomb [
99] lattices. In all three lattice types, the critical density associated with the I–N transition,
, follows a power-law scaling of the form
[
98]. The existence of this first transition has also been rigorously proven [
100].
While the second transition has been less understood, Ref. [
101] suggested it is continuous for
, but more recent findings [
106] argue it is first-order, showing MC evidence of phase coexistence for
.
Let us now return to the case of square lattices, which is the subject of this section. For
, the first transition occurs at
,
, and the second transition at
,
[
101,
104].
We focus on the entropy as a function of density, as derived from
statistics for both isotropic and nematic phases as
k increases. Two levels of approximation are considered: (i) a first-order approximation, where
is constant and the entropy of the isotropic phase at saturation is matched to Monte Carlo (MC) values, i.e.,
and
is fixed using Equation (
297). Additionally, the effect of assuming vanishing entropy at saturation (
,
) on the transitions and critical coverages is discussed; (ii) a second-order approximation, where
is a slowly decaying linear function of the density, as discussed in footnote
11. In both cases, the phase transitions and corresponding critical densities are predicted, with the first-order approximation yielding reasonable agreement with MC data and the second-order approximation showing very accurate results.
The prediction of transitions and their critical points is obtained by analytically evaluating the entropy per site as a function of
n for a fully aligned (nematic) phase,
, and for an isotropic phase,
, at the same occupation
n. The factor 2 accounts for the two possible orientations (states) per site in the isotropic phase, in contrast to one in the aligned nematic phase. The entropy per site as a function of coverage
is presented in
Figure 16 for
, and 8, based on Equations (
285), (
294), (
297), and (
298).
For
k-mers in a fully aligned nematic phase, the parameters are
,
,
,
, and
according to Equation (
298). For the isotropic phase, we have
,
, and
;
and
are computed using Equations (
297) and (
298). The entropy values
are taken from [
90,
91,
107] for
to 10, and for
,
and the corresponding values of
were reported in [
19]. The critical coverages
are determined by solving the equation
.
This equation yields two solutions
and
only for
; for
there is no solution other than the trivial equality at zero coverage. This behavior is shown in
Figure 16.
In particular, for
,
Figure 16(b) shows that
for
, indicating that the isotropic phase is the only stable phase. The same holds for
. Notably, for
, the smallest entropy difference occurs near
.
For , there are two critical coverages: and in the first-order approximation. This dual-transition pattern is found for all , consistent with the results discussed, where a mixture of cross-excluding of differently oriented species is considered. Although the nature of the transitions is not addressed here, there it finds that the I→N transition is continuous, indicating that this formalism does not matches a typical mean-field approach. For the sake of reference, in Bethe lattices with coordination q, a first-order transition occurs for depending on q [?] for this nematic transition.
In the second-order approximation,
is assumed to decay slowly with density as
, with
and
. For
, the critical coverages for
are
and
, which agree remarkably well with MC values
and
[
101].
Figure 17 shows the variation of
with
k. As
,
and
.
This model predicts a sequence of phases: isotropic at low/intermediate density, nematic at higher density, and isotropic again at high coverage, in agreement with MC results for . Even assuming vanishing saturation entropy, the formalism still predicts an I-N transition for , though no N-I transition occurs.
Figure 18 compares
results with MC data for
and 8, showing good agreement in both phases, especially at intermediate coverage. Discrepancies at high density explain differences in
, since it equals the derivative of entropy with respect to
(in
units).
For
, Equation (
294) becomes exact only in 1D. In 2D systems,
due to allowed local configurations at high coverage [
90,
91]. Thus,
underestimates entropy at high density. This distinction is critical for understanding the N-I transition for
, as seen in the inset of
Figure 18.
While already captures high coverage entropy via , an empirical correction improves quantitative agreement. We define , where . The term matches the MC saturation entropy, ensures , and the exponential captures high coverage behavior.
Appropriate values of
and
reproduce MC results in both isotropic and nematic regimes. Since
, the correction is significant only very close to saturation. Saturation entropy values are taken as
for
, and
for
[
90,
91,
107].
This empirical correction also accurately reproduces the entropy of more complex particles such as trimers (straight or bent) and triangles on triangular lattices [
110].
Concerning with the coverage dependence of the chemical potential from
statistics, following, we reproduce the results
for two very different k-mer’s size,
and
as illustrative examples in
Figure 19 [
18].
The variables are expressed in terms of coverage
, via the substitution
. For further details on the simulations, see Refs. [
17,
101,
111,
112].
In
Figure 19, MC data for
k-mers are compared with two analytical predictions. The dashed line corresponds to the limiting case where entropy vanishes at saturation, i.e.,
,
. When the correct value
is used, the
approximation (through Equation (
294)) provides better agreement with MC at high coverage, as seen in the solid line.
The discrepancy at high coverage for
between theory and MC is attributable to the behavior of the entropy
near
in Equation (
285). A qualitatively similar result is obtained using the empirical form
introduced in the previous section.
For
, MC simulations show that
k-mers undergo a continuous I-N transition at intermediate
[
96], consistent with predictions from the
model based on entropy comparison between nematic and isotropic phases. Despite this, Equation (
286) still gives a fair approximation for
if a value of
smaller than that from Equation (
300) is used.
This is understood as follows: nematic ordering forms compact bundles of neighboring k-mers. For a fixed particle number N, such alignment leads to more multiply excluded states and fewer total excluded states per particle, effectively reducing the parameter . A simple illustration of this can be made by comparing the number of states excluded when two k-mers are perpendicular vs. parallel and aligned.
For instance, in
Figure 19, the case
with
yields a good fit to MC data. This value of
is significantly smaller than the isotropic value
obtained from Equation (
300) reflecting bundled-like configurations of the lattice gas.
The MC simulation results presented in this section were obtained using the efficient algorithm introduced by Kundu et al. [
101,
111,
112] and described in
Section 8.5. Simulations were performed on
square lattices with periodic boundary conditions. The ratio
was set to 120, a value for which finite-size effects were found to be negligible. Equilibrium was typically reached after
MCSs, and observables were computed by averaging over
configurations.
Illustrative results for
from
statistics for squares and rectangles on square lattices were presented and compared with fast-relaxation Monte Carlo simulations in Ref. [
18]. A comprehensive study of the various phases formed by rectangles on square lattices, and other hard-core lattice gases, can be found in Refs. [
111,
113,
114,
115].
5.6. Exclusion Spectrum Functions
A singular outcome of this formalism is the thermodynamic characterization of the configuration space or state exclusion spectrum through the exclusion per particle frequency function and the cumulative exclusion per particle function , referred to collectively as exclusion spectrum functions.
From Equations (
291) and (
292), both
and
can be expressed in terms of the density dependence of the chemical potential, providing a thermodynamic description of equilibrium particle configurations.
Illustrative results for
and
for
k-mers, squares, and rectangles in the isotropic phase are shown in
Figure 20 and
Figure 21. Analytical results are shown as lines, and MC simulation data are indicated with symbols:
rectangles (triangles),
k-mers (circles),
squares (diamonds), and
k-mers (squares). The simulations were performed following the same scheme and using the same parameters as those used in
Section 5.5.
For
with
, both
and
decrease rapidly with increasing coverage:
and
. In contrast, in a 1D lattice with
, the decay is slower:
and
(not shown in
Figure 20 and
Figure 21 for clarity).
In the other cases shown: for rectangles, , , , ; for squares, , , , ; for k-mers, , , , . A good agreement is observed between analytical predictions and MC data for both and , particularly for . At saturation, , since each particle excludes on average g states, and , as all single-particle states are either occupied or excluded.
As shown in
Figure 20 and
Figure 21,
statistics provides an accurate description of the exclusion spectrum functions across the full range of
. In particular, the results for
k-mers show excellent agreement for both small and large
k. Moreover, particles with higher
g tend to exhibit larger values of
, regardless of shape. In contrast,
captures more detailed configuration-specific features, as evidenced in
Figure 21. For instance, while isolated
k-mers exclude more states than
squares, there exists an intermediate coverage range
where they exclude fewer states per particle. This indicates local alignment among
k-mers at high
, reducing exclusion relative to a disordered configuration.
These results confirm that statistics captures the thermodynamic signature of configurational exclusion with remarkable accuracy across all densities. A more in-depth analysis of the exclusion spectrum and its behavior near phase transitions is provided discussed in following sections.
5.7. Adsorption of Polyatomics: Relation Between Exclusion Functions, Thermodynamic Observables, and Adsorption Field Topology
In this section, we explore potential applications of the exclusion spectrum functions defined in
Section 5.1 in connection with experimental thermodynamic measurements. This relation provides insight into how adsorbed particles occupy and exclude states based on the spatial distribution of local minima in the adsorption field—here generically referred to as the adsorption field topology.
The average exclusion spectrum function
connects a configurational property, related to the spatial correlation of states and influenced by particle geometry, to the density dependence of a thermodynamic observable such as the chemical potential. From the relation
or equivalently , it becomes apparent that these exclusion functions can, in principle, be inferred from experiments via the dependence of on n.
A more refined exclusion description is provided by the frequency function
. Introducing
, we have
and
. Hence, the analytical or experimental form of
encapsulates the configurational exclusion information. Thus, for a given shape and size of adsorbate molecule the number of states at very low coverage could be determined and the spatial arrangement of the adsorption potential minima could be inferred, so called here, the adsorption potential topology [
18]. Additionally, the complete configuration changes on density are embodied in
and
through
. A more detailed experimental analysis of adsorption isotherms on well-defined particle–substrate systems is needed to assess the feasibility and value of this configurational framework, though that lies beyond the scope of the present review.
6. Latest Developments, Part II: Multiple Exclusion Statistics Formulation for Mixtures
In this section the statistical thermodynamics framework is extended to describe mixtures of particles with arbitrary size and shape, each having a spectrum of topologically correlated states and subject to statistical exclusion. A generalized distribution is obtained from a configuration space ansatz recently proposed for single species, accounting for the multiple exclusion phenomenon, where correlated states can be simultaneously excluded by more than one particle. Statistical exclusion on correlated state spectra is characterized by parameters , which are self-consistently determined. Self- and cross-exclusion spectral functions and are introduced to describe the density-dependent exclusion behavior. In the limit of uncorrelated states, the formalism recovers Haldane’s statistics and Wu’s distribution for single species and for mixtures of mutually excluding species with constant exclusion.
The formalism is latter applied to k-mers on the square lattice, modeled as a mixture of two orthogonally oriented, self- and cross-excluding pseudo-species. This approach offers a general and consistent framework for entropy-complex lattice gases. It reproduces k-mer phase transitions and provides access to configurational information through the exclusion spectrum functions. We summarize here basis of this approach for mixtures and the relevant analytical predictions for rigid k-mers on the square lattice are discussed in the section devoted to applications.
6.1. State Counting Approximation and Density of States for Mixtures with Multiple State Exclusion
The general self-consistent formulation for the thermodynamics of mixtures consisting of an arbitrary number of species in a volume
V is summarized in the following subsections. Each species are assumed to exclude accessible states to itself and to the others, a phenomenon intensified by spatial correlations among the states—referred to as Multiple Exclusion Statistics (
) [
17]. This leads to a particularly challenging statistical problem, especially in lattice models involving linear or arbitrarily shaped particles.
Figure 22 illustrates
through a ternary mixture of monomers, dimers, and tetramers. The single-species theory of Ref. [
18] was extended to mixtures and state density and exclusion distribution functions were derived [
19] in order to formulate an analytical thermodynamic model applicable to
k-mers on lattices with correlated spectra as well to other particle shapes.
We define the self- and cross-exclusion parameters and based on the number of states excluded per particle at saturation. Let be the total states available to an isolated particle of species i and its number in V. The occupation number is , with its maximum, and .
Cross-exclusion is rather more subtle-quantified by , where are the non-excluded states of i when j saturates the system. Expressed in terms of fractions: , with and .
The canonical partition function is , where is the energy per particle(eventually, due to an external field such as the interaction with the lattice)and henceforth. The configurational term captures how particles distribute over their respective sets of accessible states.
We define
as the number of states available to a particle of species
i given occupation vector
(and analogously for
). Generalizing the form introduced in [
18] for single species, the configuration count is:
Using Stirling’s approximation and defining
, the Helmholtz free energy becomes:
In the thermodynamic limit, we define the free energy per state of species
i as
. This leads to:
where
.
The entropy per state (in units of
) follows from Equation (
306):
In general, depends not only on the occupation vector but also on the microstates, which are inaccessible analytically. Thus, we postulate a functional form based on average occupations—sufficient for thermodynamic descriptions. The next section formalizes the derivation of based on a counting ansatz and pairwise exclusion analysis.
It is worth noting that the expression for
is exact only if particles can occupy states that are completely independent of one another. This is the case in the Haldane-Wu
g-statistics framework discussed in Ref. [
18], where each particle excludes
g states regardless of
N or the specific configuration. For a single species, this leads to
with constant
g, while for a mixture, the generalized form becomes
with constants
. In general, however,
in Equation (
304) is only approximate, since the actual number of available states for species
i should depend not only on
but also on the specific microscopic configuration of the ensemble; that is, it should be configuration-dependent.
Yet, since exact configuration counting is intractable, we develop a general approximation for based on a state counting ansatz. In the thermodynamic limit, this approach yields the density of states as a function only of the average occupation numbers at equilibrium. This can be interpreted as the effective density of states corresponding to typical equilibrium configurations which contribute most significantly to the system’s entropy.
We now show the derivation of the functional form of . Let denote the set of available states for a single particle of species i, with cardinality . We define the tuple , whose components correspond to the total number of states for each species.
To quantify how the presence of other species modifies the state space of species i, we denote by the number of available states for a particle of species i when the system contains particles in total, satisfying when for all j. We first isolate the effect of species j on species i by defining under the condition , with only species j present.
The recursion relation introduced in
Section 5 for single species can be extended to such pair interactions. The function
is then defined recursively as follows:
in general,
where
denotes the number of states of species
i excluded by the
-th particle of species
j added to the system.
Following the analogy with the single-species case, we posit that . Here, the term 1 accounts for the exclusion of at least one state, while represents the additional number of excluded states due to spatial correlations between the state spectra of species i and j.
The correlation term
is defined by the state counting ansatz [
17] as
. Substituting this into the recursive definition yields:
and so forth. The general expression becomes:
valid for
.
For convenience, we introduce in Equation (
309) the rescaled exclusion parameters
, and similarly
, where
and
.
By retaining only the leading term of the summation in Equation (
310) and taking the thermodynamic limit
, with
,
, and
, we obtain:
which yields
, in agreement with the approximation proposed in Ref. [
17]. Here,
represents the fraction of states of species
i that remain available in the presence of a concentration
of species
j, under the condition that all other species are absent (
for
). Statistically, this function encodes the depletion of the state spectrum of species
i due to the presence of species
j, and can be interpreted as a pairwise statistical interaction function.
The function
must satisfy the boundary conditions
and
. To ensure this, we define constants
and
such that
. Imposing the boundary conditions yields
and
. Therefore, the final expression for the pair function becomes:
Because of the spatial self- and cross-correlations, state multiple exclusions occur for a given microscopical configuration of the statistical ensemble (as shown in
Figure 22), then, the total fraction of states excluded to a given species
i by the others species is not merely
, being
the fraction of states of
i excluded by
j at
. Accordingly, the fraction of states for a particle of species
i when all the species are coexisting, namely
, corresponds to the ratio between the number of states in the intersection sets
and the total number of states
. Ultimately, the total fraction of states for a single particle of species
i at occupation
of all the species is given by
, which can be approximated by
.
Figure 23.
Symbolic representation of the species i’s states set (whole framed area) whose elements are the states accessible to species i when for , being its cardinality. represent the sets of states of a particle of species i excluded by the particles of species , respectively (shown generically and by the areas filled by vertical and horizontal lines), being their cardinalities, respectively. The states occupied by species i are represented by the set (oblique lines area). represents the set of states for particles of species i not excluded by particles of species j. The intersection is the set of states for a particle of species i non-excluded by any of the species (dark gray area), with cardinality and fraction .
Figure 23.
Symbolic representation of the species i’s states set (whole framed area) whose elements are the states accessible to species i when for , being its cardinality. represent the sets of states of a particle of species i excluded by the particles of species , respectively (shown generically and by the areas filled by vertical and horizontal lines), being their cardinalities, respectively. The states occupied by species i are represented by the set (oblique lines area). represents the set of states for particles of species i not excluded by particles of species j. The intersection is the set of states for a particle of species i non-excluded by any of the species (dark gray area), with cardinality and fraction .
It is worth noticing that the fraction
can also be interpreted as the probability for a state of species
i being non-excluded by particles of species
j. Accordingly,
represents the probability that a state of species
i is simultaneously non-excluded by all species. Then, assuming in a first approximation that the pairs cross-exclusion events are independent,
can be written as
where the functions
are given by Equation (
311).
In order to finally obtain
as defined in Equation (
306) from
in Equation (
312), and by denoting
, it is required that
satisfies
and
13, where
is the maximum occupation of species
i when other species are at densities
.
Introducing normalization constants
and
so that
, and using
, the boundary conditions are fulfilled by setting
and
. Then, the density of states reads
for
, with
given by Equation (
311).
The Helmholtz free energy for the generalized mixture on spatially correlated states can now be computed from Equations (
311), (
313), and (
306), providing the full thermodynamic behavior.
As will be discussed in the applications, the quantities
can be predetermined in model systems such as
k-mer mixtures on regular lattices from the system symmetry
14.
Finally, note that the set cardinalities in this derivation do not depend on specific microstates, but represent effective values for configurations at equilibrium minimizing the Helmholtz free energy at given T and V (i.e., given ).
The parameters
are consistently determined from particle and lattice properties (size, shape, connectivity), and from thermodynamic boundary conditions by generalizing the analysis developed introduced for single-component systems [
18] .
6.2. Mixtures Statistical Thermodynamics
From the Helmholtz free energy in Equation (
306) we derive the density dependence of the chemical potential
,
consequently, from Equations (
306) and (
314),
where, for the sake of shortness, the explicit dependence of
on
is implicit in Equation (
315) and
. We write Equation (
315) more conveniently as
with
for
which straightforwardly give the chemical potentials
as a function of the species state occupation numbers
, or specie’s density. The Equation (
316) also represents a system of
s-coupled equations whose solutions are the species occupation numbers
for given chemical potentials
.
By defining
, Equation (
316) can be rewritten as
or
which are the coupled equations within the
statistics from which the equilibrium distributions
can be determined.
It is worth noticing that Equation (
318) reduces to Wu’s distribution for fractional exclusion statistics [
15] when the species’s states are spatially uncorrelated. An alternative picture of what the Wu’s limiting case of spatially uncorrelated species’s states means here is that the numbers of self-excluded and cross-excluded states per particle are constants in Wu’s formalism [
15] and density dependent in
statistics.
6.3. State Exclusion Spectrum Functions: Determination of Exclusion Correlation Parameters
This section is devoted to determining the parameters from the thermodynamic limits of the exclusion spectrum functions, which quantify the average cumulative number of self-excluded and cross-excluded states per particle as functions of the occupation numbers . For model systems where both the size/shape of the particles and the spatial distribution of accessible states (e.g., the lattice geometry) are known, the values of can be determined within the statistics framework, enabling a complete thermodynamic description.
Assume that species
in volume
V can exchange particles with a reservoir at temperature
T and chemical potentials
. The time evolution of the mean occupation number
follows:
where
and
denote the average fractions of empty and occupied states for species
i, and
,
are the respective transition rates. At equilibrium,
and the detailed balance condition gives
. Since
, then from Equation (
316), we have
Substituting from Equation (
316), we obtain:
The generalized exclusion spectrum function is defined as
, which measures the total average fraction of states excluded to a particle of species
i at given occupations
(including both self- and cross-exclusion). The cumulative number of excluded states of species
i per particle of species
j, i.e., the spectrum function
, is defined as:
Additionally, the rate of excluded states per particle due to self- or cross-exclusion is given by the partial derivatives:
Let
,
, and
, which yield:
From Equations (
320) and (
322), it can be shown:
Solving for the self-exclusion parameters
,
where
refers to Lambert function (see
Section 5.4).
For the cross-exclusion case
For uncorrelated states (
),
and
, as in Wu’s formalism [
15].
The exclusion spectrum functions
are related to density and chemical potential through
Although this formulation links measurable thermodynamic quantities to exclusion spectrum functions, its experimental application, especially for mixtures, requires further work. In dicussing applications further on in
Section 7.8, we illustrate how
and
can be computed for
k-mers on a square lattice, their behaviour through the
k-mers transitions and usefulness display and characterize the order of transitions.
6.4. The k-Mers Problem as a Mixture Model: Basic Definitions
Preliminarly to applications of
statistics, which we will discuss in
Section 7.8, we introduce here the problem of
k-mers on square lattice of
M sites rationalized as a mixture of two differently oriented species.
As discussed in
Section 5.5, the adsorption of straight rigid
k-mers on square lattices for
exhibits two distinct phase transitions: (1) a continuous, entropy-driven isotropic-to-nematic (I-N) transition occurring at intermediate surface coverage [
96,
101,
104], and (2) a nematic-to-isotropic (N-I) transition taking place at densities approaching lattice saturation [
106].
In Ref. [
18]
k-mers on the square lattice has been modeled as a binary mixture of species aligned along the horizontal (
) and vertical (
) lattice directions, denoted as species 1 and 2, respectively. Both species occupy
k consecutive lattice sites along their respective directions. According to our definitions:
,
,
, and
. The exclusion parameters are
, and the saturation values satisfy
.
From Equation (
313), the saturation occupations under coexistence are:
The self-exclusion at infinite dilution is , leading to . For cross-exclusion, each k-mer excludes states orthogonal to its direction, of which are shared. Thus, .
Using Equation (
326), the cross-exclusion correlations are
where
is the Lambert function. Solving this yields
for
,
for
,
for
, and
for
.
Based upon this elementary definition of the mixture parameter, the thermodynamic and exclusion functions in the
statistics, a much comprehensive treatment of the problem is given in
Section 7.8, leading to the entropy surface, equilibrium paths, density branches, order parameters, transition critical points and state exclusion spectrum are obtained for various values of
k.
7. Applications
In the present section, we analyze the scope and limitations of the theoretical models developed in previous sections by comparing them with Monte Carlo simulations and experimental data available in the literature.
7.1. Two-Dimensional Adsorption: Comparison Between Theory and Monte Carlo Simulations
In this section, adsorption isotherms are calculated for the theoretical models introduced in
Section 3 (
,
,
,
, and
) and compared both among themselves and against Monte Carlo simulations performed within the grand canonical ensemble framework [see
Section 8]. These comparisons are conducted for honeycomb, square, and triangular lattices.
Monte Carlo simulations were carried out on honeycomb, square, and triangular lattices of size , with , and 150, respectively, using periodic boundary conditions. This lattice size ensures that finite-size effects are negligible. In addition, MCSs.
We begin by discussing some fundamental features of the adsorption isotherms.
Figure 24 presents a comparison between the exact adsorption isotherm for monomers and the simulated isotherms for dimers on honeycomb, square, and triangular lattices. As observed, the particle–vacancy symmetry, which holds for monoatomic adsorbates, breaks down when
. Furthermore, while the dimer adsorption isotherms appear similar across the different lattice types, the curves shift to lower values of
as the connectivity
increases. In other words, for a given value of
, the equilibrium surface coverage rises with increasing
. This behavior can be explained using the following relation:
which is valid for linear
k-mers at low concentrations [see Equation (
149)]. As the chemical potential increases, this effect diminishes, and consequently, the slope of the isotherms decreases with increasing
.
We now consider the case of linear adsorbates larger than dimers. For honeycomb lattices,
k-mers adsorb as described in
Section 3.1.1. When a site is selected, there are six possible equilibrium orientations for a single
k-mer (
) at extremely low coverage, resulting in a total number of
k-uples equal to
(i.e.,
), as is also the case for triangular lattices.
Based on these conditions, extensive simulations were conducted for linear adsorbates with
k ranging from 2 to 10. As illustrative examples,
Figure 25(a),
Figure 26(a), and
Figure 27(a) compare simulation isotherms with theoretical predictions for 6-mers on honeycomb, square, and triangular lattices, respectively. In all cases, the theoretical models agree well with simulations at low coverages, but deviate significantly as the surface coverage increases.
The discrepancies between simulation and theory can be quantified by the percentage reduced coverage, defined as [
73]
where
(
) represents the coverage obtained by using MC simulation (analytical approach). Each pair of values (
) is obtained at fixed
.
Figure 25b,
Figure 26b, and
Figure 27b show how
varies with surface coverage for the different lattice types. The performance of each theoretical model is as follows: the
model (dashed line) shows very good agreement with simulation results, with minimal discrepancies. Both the
(dash-dot-dot line) and
(dash-dot line) models tend to underestimate the coverage across the entire range. The
model (dotted line) performs poorly at intermediate coverages but improves at high coverages. With respect to lattice connectivity, the accuracy of
and
improves as
decreases, while the opposite trend is observed for
and
. The behavior of the
and
models supports the formulation of the semi-empirical (
) isotherm (solid line) given in Equation (
155). This trend is also illustrated in
Figure 28,
Figure 29, and
Figure 30, which show the percentage reduced coverage for the
model as a function of concentration for various values of
and
k.
To better interpret the data shown in
Figure 28,
Figure 29 and
Figure 30, we consider two summary metrics: (1) the average absolute difference between simulation and theoretical coverage,
; and (2) the maximum value of the percentage reduced coverage,
. These metrics are presented in
Figure 31. Several conclusions can be drawn: (i) theoretical models generally perform better on square lattices; (ii) both
and
remain nearly constant for
k values between 2 and 8; and (iii) both quantities increase for
. Finally, since the values of
remain below
, the
model can be considered a reliable approximation for describing multisite occupancy adsorption, at least for the
k values analyzed here.
7.2. Two-Dimensional Adsorption of Binary Mixtures: Comparison Between Theory and Monte Carlo Simulations
In this section, we analyze the main characteristics of the theoretical approximations developed in
Section 3.3 in comparison with MC simulation results and experimental data.
We consider a gas mixture composed of rigid rods of lengths
k and
l, with each component present in equal molar fraction in the gas phase. Adsorption occurs on a regular, homogeneous lattice with connectivity
(square lattice) and 6 (triangular lattice). The parameters used in the HPTMC simulations (see
Section 8.3) were:
, and
MCSs. For simplicity, we assume the standard chemical potentials are zero and that the equilibrium constants
for all species.
The partial adsorption isotherms for
,
, and two values of
are shown in
Figure 32(a) and
Figure 32(b), respectively. Symbols denote MC simulation data, while the lines represent various theoretical approaches as indicated.
A characteristic feature observed in binary mixtures of polyatomic species is evident in both figures: at higher pressures, the smaller species are displaced by the larger ones. This phenomenon, known as adsorption preference reversal (APR), arises from entropic competition between adsorbed species. A detailed investigation of the APR effect, focusing on the impact of molecular size difference, is presented in
Section 7.4.
In
Figure 32(a), the most significant discrepancies between simulation and theoretical predictions occur in the partial adsorption isotherm of the larger species. As clearly shown in
Figure 32(b), the magnitude of these deviations depends on the level of approximation used in evaluating
[
11].
To quantitatively assess the agreement between simulation and analytical results, we use the reduced coverage error defined in previous section [
73],
where
represents the value of the total coverage obtained by using the MC simulation (analytical approach). Each pair (
,
) corresponds to the same value of pressure
P.
The results obtained from Equation (
332) are shown in
Figure 32c,d. In
Figure 32c, better agreement is observed at surface coverage above
, where only dimers are present, a case where all theoretical models perform well [
16]. Conversely, for larger adsorbates, classical approaches fail to accurately describe the adsorption behavior across the entire range of coverage, as illustrated in
Figure 32d. In contrast, the
approximation yields satisfactory results, with errors remaining below
in both cases.
The same analysis presented in
Figure 32 is repeated for triangular lattices, with the corresponding results shown in
Figure 33. In most cases, the error increases with lattice connectivity. However, the
approximation consistently yields an error below
, indicating that it is a highly reliable method for modeling binary mixture adsorption with multisite occupancy, at least for the molecular sizes considered in this study.
To complete our analysis, we examine partial isotherms for mixtures with varying sizes of
l, keeping
fixed. For each pair (
k,
l), the discrepancies between theoretical predictions and simulation results are assessed using the average error across the full coverage range, defined as
where
R is the total number of isotherm points (or the number of replicas, as described in
Section 8.3).
Figure 34 illustrates the dependence of
on the size of the
l-mers and the lattice geometry for the various theoretical approaches evaluated in this work. As shown,
increases monotonically with
l, indicating that the divergence between MC simulations and analytical models becomes more significant for larger adsorbates. In contrast, the error associated with the
approximation remains nearly constant, around
for square lattices and
for triangular lattices. This excellent agreement across all values of
l highlights the robustness of the
method and underscores the importance of accurately computing the correction function [Equation (
187)] for understanding the adsorption behavior of rigid rod mixtures.
Finally, to evaluate the applicability of the proposed model, we analyzed experimental data extracted from Ref. [
116]. Specifically, adsorption isotherms of hydrocarbon mixtures—methane–ethane and ethane–propylene—on activated carbon (AC-40) at
C were examined using the
adsorption model introduced in this work. Since the experimental data were reported as adsorbed amount (in
) versus pressure (in mmHg), the theoretical isotherms were reformulated in terms of pressure
P and adsorbed amount
g to enable direct comparison and fitting.
Assuming equilibrium between the adsorbed phase and an ideal gas-phase mixture, the chemical potentials were related to the system pressure and molar fractions. Additionally, the coverage was defined as , where represents the maximum adsorption capacity of the surface.
Following a common approach in the literature, a “bead segment” chain model was employed in which each CH
n group (bead) occupies one adsorption site on the surface. Accordingly, each hydrocarbon species C
m was modeled as a rigid rod of length
[
117].
Within this framework, the experimental isotherms for methane–ethane and ethane–propylene mixtures at
C and varying molar fractions were fitted using a single value of
and temperature-dependent equilibrium constants
as adjustable parameters. The results of the fitting procedure are shown in
Figure 35, and the corresponding parameter values are summarized in
Table 1. A very good agreement is observed between experimental data (symbols) and theoretical predictions (solid lines).
While a more extensive analysis of experimental adsorption isotherms is still needed, these results suggest that the theory provides a promising and accurate framework for describing the adsorption thermodynamics of interacting polyatomic species.
7.3. Two-Dimensional Adsorption of Interacting k-Mers: Comparison Between Theory and Monte Carlo Simulations
In this section, we analyze the main features of the thermodynamic functions derived from the models presented in
Section 4.1 (Bragg-Williams Approximation,
) and
Section 4.2 (Quasi-Chemical Approximation,
), in comparison with Monte Carlo simulation results for a lattice-gas of interacting dimers on honeycomb, square, and triangular lattices.
15
As in
Section 7.1, simulations were conducted on honeycomb, square, and triangular lattices of size
, with
, 144, and 150, respectively, using periodic boundary conditions. Moreover, the lattice size
L was carefully selected to avoid perturbation of the adlayer structure.
Representative adsorption isotherms obtained from Monte Carlo simulations in the grand canonical ensemble (symbols), along with their comparison to the
(solid lines) and
(dashed lines), are shown in
Figure 36,
Figure 37, and
Figure 38 for honeycomb, square, and triangular lattices, respectively.
For attractive interactions [
Figure 36a,
Figure 37a, and
Figure 38a], as temperature decreases, a first-order phase transition occurs, evidenced by the discontinuity in the simulated isotherms and the appearance of characteristic loops in the theoretical curves. This behavior, experimentally observed in many systems [
80], corresponds to a low-coverage lattice-gas phase coexisting with a higher-coverage “lattice-fluid” phase. The lattice-fluid can be considered as a diluted version of the registered
phase where all lattice sites are occupied except for some vacancies.
This two-dimensional gas-to-liquid condensation closely resembles that of a monomeric lattice gas with attractive interactions. However, for k-mers, the particle-vacancy symmetry (valid for monomers) is broken, resulting in adsorption isotherms that are asymmetric with respect to .
In the case of repulsive interactions [
Figure 36b,
Figure 37b, and
Figure 38b], the isotherms exhibit more complex features due to the formation of ordered structures in the adsorbed layer. These ordered arrangements are indicative of subcritical behavior, where continuous phase transitions occur from disordered to ordered phases [
118,
119]. At high temperatures, isotherms remain featureless, but at low temperatures, they display distinct steps corresponding to the emergence of ordered phases. The specific form of these steps depends strongly on the lattice connectivity. As the chemical potential
increases and the surface coverage
spans from 0 to 1, two ordered phases are typically observed: (1) a low-coverage ordered phase (LCOP), characterized by site occupancies of
,
, and
for honeycomb, square, and triangular lattices, respectively; and (2) a high-coverage ordered phase (HCOP), with
site occupancy across all three geometries. Snapshots of the LCOP [part a)] and HCOP [part b)] configurations for each lattice type are presented in
Figure 39,
Figure 40, and
Figure 41. For a detailed discussion of these phases, refer to Refs. [
118,
119].
Under attractive interactions, both theoretical models yield qualitatively similar results, and the isotherms from
and
are nearly indistinguishable. However, it is known that isotherms derived from fundamentally different approximations can appear deceptively similar [
67]. To better assess the accuracy of each model, we use the absolute error in the chemical potential,
, which is defined as
where
(
) represents the chemical potential obtained by using MC simulation (analytical approach). Each pair of values (
,
) is obtained at fixed
.
As an example,
Figure 42(a) presents
for three representative attractive interaction strengths: squares for
, triangles for
, and circles for
. Solid and open symbols correspond to
and
results, respectively. In all cases,
outperforms
.
The corresponding analysis for repulsive interactions is shown in
Figure 42b, which includes
(squares),
(triangles), and
(circles). Again, solid and open symbols represent
and
, respectively. Here, differences between the two models are both quantitative and qualitative. While
fails to predict any ordered structures,
captures the formation of a pronounced plateau at low temperature. This critical coverage,
, appearing between the LCOP and HCOP, depends on both lattice geometry and adsorbate size. The adsorbate configuration at
can be interpreted as a mixture of LCOP and HCOP phases [see part (c) in
Figure 39,
Figure 40, and
Figure 41].
The curves in
Figure 42 correspond to a honeycomb lattice. However, the behavior of
for square and triangular lattices is very similar (data are not shown here for sake of simplicity).
To quantify the overall deviation between theory and simulation across the full coverage range, we define the integral error
as:
Figure 43 shows
for all lattice geometries and a broad range of
values. Several key conclusions can be drawn:
In all cases,
provides a significantly better fit to the simulation data than
. This is particularly true for repulsive interactions, where
shows large discrepancies, while
remains the simplest yet effective model for describing multisite occupancy adsorption.
The value of
increases with lattice connectivity. This may be attributed to a loss in accuracy of
as
increases [
73].
There exists a broad range of interaction strengths (
) for which
matches the simulation data extremely well. Notably, most surface science experiments fall within this range of interaction energies.
Therefore, not only represents a clear improvement over the in modeling k-mer adsorption but also provides a solid theoretical framework and compact expressions for the interpretation of thermodynamic adsorption data of polyatomic species—such as alkanes, alkenes, and other hydrocarbons—on regular surfaces.
7.4. Application of to the Adsorption of and in Zeolites 13X and 5A: Determination of the Adsorption Configuration
One interesting application of the theoretical framework presented in
Section 3.2.2 on a lattice gas model involves the interpretation of experimental adsorption isotherms for propane [
120] and oxygen [
121,
122] in 13X and 5A zeolites, as well as in simulation-based systems. In our approach, we employed Equation (
125) under two main assumptions:
since
constant, if a single molecule has
m distinct ways to adsorb per lattice site at zero coverage, then the presence of an adsorbed
k-mer, occupying
sites, effectively excludes
states from being accessible to other molecules [thus,
]; and
the energetic contribution from adsorbate-adsorbate interactions is accounted for using a mean-field approximation, as described in
Section 4.1. This analysis highlights the physical interpretation of the parameters
g and
a, linking them to the spatial configuration of adsorbed molecules and the geometric structure of the surface.
Because the experimental data are presented as adsorbed volume
v, against pressure
p, we rewrite Equation (
125) in the more convenient form:
where
(
is the volume corresponding to monolayer completion);
;
is the equilibrium constant
; and
is the mean-field term. In addition,
can be associated with the isosteric heat of adsorption
.
Figure 44 presents the adsorption isotherms of propane (
) in
zeolite. Solid lines represent theoretical predictions using the
model, while symbols show experimental data from Ref. [
120]. Following conventional modeling, alkane chains are treated as "bead segments," where each methyl group corresponds to one adsorption site. Accordingly, propane is modeled as a trimer with
. Given that the propane molecule (
) is relatively large compared to the cavity diameter (
), it likely adsorbs along a preferred orientation. Otherwise, accommodating 5–6 molecules per cavity would be unfeasible. We thus assign
(with
and
, mimicking a 1D configuration). The best fit to the experimental data over the full temperature and pressure range was obtained by simultaneously optimizing
,
and
w (see
Table 2).
Consistent with experimental findings, the resulting is slightly under 6 molecules per cavity. A value of suggests that some molecules may partially span across the cavity windows. Regarding the lateral interaction parameter w, its ratio to the known liquid-phase interaction energy for propane is about , indicating that each molecule interacts, on average, with 2.5 neighbors at full coverage. This supports the approximation of a quasi-one-dimensional system.
Figure 45 also shows strong agreement between theory and experimental data for a system with non-monotonic adsorption behavior. In this case, the derivative of the adsorbed volume with respect to pressure was fitted, again yielding excellent correlation.
To quantify model accuracy, we define the deviation
D as the average relative discrepancy (in percent) between theoretical (
) and experimental (
) data:
where
i runs over the total set of data.
For propane on 13X,
D was found to be
(
Table 2), which is within the bounds of experimental error, underscoring the robustness of the
approach. In contrast, a previous model from Ref. [
120] required eight parameters to describe similar systems. Here, the complexity of polyatomic adsorption is captured through a single, physically meaningful parameter
g, reflecting the spatial configuration of the adsorbate.
We now turn to the oxygen adsorption isotherms in 5A zeolite, shown in
Figure 46. Symbols denote experimental measurements, while lines indicate theoretical curves from Equation (
336). Experimental data were taken from two different sources in the literature, from Miller et al. [
121] and Danner et al. [
122]. In the first set of data (empty symbols) [
121], the amount adsorbed was measured in units of the number of molecules per cavity. In the other case (full symbols) [
122], the amount adsorbed was reported in units of
per
of adsorbent.
Figure 46.
Adsorption isotherms for
adsorbed in
zeolite. Empty and full symbols correspond to data from Refs. [
121,
122], respectively. Lines correspond to the adsorption isotherm function of Equation (
336).
Figure 46.
Adsorption isotherms for
adsorbed in
zeolite. Empty and full symbols correspond to data from Refs. [
121,
122], respectively. Lines correspond to the adsorption isotherm function of Equation (
336).
The fitting process involved two steps:
based on previous numerical simulations [
123], we fix
(
and
). Under these considerations, analytical isotherms in
Figure 46 were obtained by multiple fitting the set of parameters
,
and
w as in
Figure 44. The value obtained for
w is in excellent agreement with the simulational calculation of
w in Ref. [
123]. With respect to
, it was not possible from the work of Razmus et al. [
124] to estimate
in order to compare with the one from Equation (
336). However,
was independently validated through a second stage of fitting;
the values of
g,
and
w arising from
were fixed. Then,
, set as to fit the experimental isotherm measured by Danner et al. [
122], agrees with the monolayer volume reported in Ref. [
122]. The deviation between experimental and fitting curves was
(see
Table 2). Based in the consistency of this analysis,
appears to adsorb flat with two possible orientations on a two-dimensional layer defined by the cavity’s inner surface.
The isosteric heat of adsorption
for both systems was derived from the slope of the plot of
versus
, as shown in
Figure 47. The results, which are presented in
Table 2, align well with reported experimental values for propane in 13X [
120] and oxygen in 5A [
124].
Finally, in order to illustrate the applicability and versatility of
to describe systems more complex than the one in the experiments analyzed here, we show in Figure (
Figure 48) the fit (solid lines) to numerical isotherms (symbols) of dimers, flexible trimers and flexible tetramers adsorbed flat on a square lattice. Solid lines represent the best fitting to computer experiments in the crudest approximation [
constant, Equation (
125)] to the general isotherm of Equation (
124). The values obtained for
g in all cases are very consistent to the ideal value
. It is worth mentioning that, at this elementary degree of approximation
is already more accurate than the classical Flory’s theory [
6] of adsorbed chains as well as the multisite-adsorption approaches of Refs. [
20,
33].
For the case of dimers, the inset in Figure (
48) shows a comparison between the configurational entropy per site,
s, versus coverage from simulation and the corresponding ones obtained from
[being
and
] for different values of
g. As it was calculated for the isotherm, the best fit corresponds to
(solid line). The three curves in dashed (dotted) lines correspond to increments (decrements) of
,
and
with respect to
. As it can be clearly visualized, small variations of
g provide notable differences in the entropy curves. Thus, the statistical exclusion parameter
g results highly sensitive to the adparticle’s spatial configuration. The more compact the configuration of the segments attached to the surface sites the smaller
g. For instance,
g may vary from
(
) for straight trimers (tetramers), to
(
) for flexible trimers (tetramers).
In summary, the study presented in this section demonstrates how the concept of statistical exclusion (
) serves as a powerful tool to quantify the entropy associated with polyatomic adsorption, effectively capturing both the configuration and interaction of adsorbates. Furthermore, Equation (
124) lays the groundwork for analyzing configurational changes as a function of coverage [through configurational spectroscopy
] from thermodynamic data. Overall,
offers a compact and insightful framework for interpreting adsorption phenomena, from simple gases to complex hydrocarbons on structured surfaces.
7.5. Adsorption of Methane-Ethane Mixtures in Zeolites: Reversal Adsorption Phenomena
The adsorption of molecular mixtures presents a significantly greater challenge than that of pure components, both experimentally and theoretically. While the number of adsorbed molecules in a pure gas system can be accurately determined by measuring the weight gain of the adsorbent sample, studying mixtures requires additional experimental techniques to determine the composition of the adsorbed phase. This complexity is one of the primary reasons for the scarcity of experimental data on the adsorption of polyatomic mixtures.
Nonetheless, several simulation studies have investigated the behavior of hydrocarbon mixtures [
125,
126,
127,
128,
129,
130,
131,
132]. A particularly striking phenomenon has been observed in methane–ethane mixtures [
125,
126]: at low pressures, the adsorbed phase is dominated by ethane, but as pressure increases, methane gradually displaces ethane. This inversion in adsorption preference has also been reported in mixtures of linear hydrocarbons adsorbed in silicalite [
127,
128,
129,
130,
131,
132], as well as in carbon nanotube bundles [
133] and metal-organic frameworks [
134]. In all these cases, selectivity shifts from favoring the larger molecule at low pressure to favoring the smaller one at higher pressure—a behavior known as Adsorption Preference Reversal (APR) [
51].
From a theoretical standpoint, understanding the origin of APR is not straightforward. Within the framework of single-site adsorption (where each molecule occupies a single lattice site), studies by Ayache et al. [
51] and Dunne et al. [
52] have shown that competition between species under repulsive lateral interactions can cause one component to replace the other in the adsorbed phase. In their approaches, Ayache et al. [
51] used a mean-field approximation, while Dunne et al. [
52] performed exact calculations. However, hydrocarbon molecules adsorbed on solid surfaces should be regarded under the light of a multisite-adsorption model, in order to properly account for the effects of configurational entropy (
k-mer size and flexibility) on the thermodynamics of the adlayer. Entropic effects arising from differences in molecular structure and density are expected to drive a range of complex behaviors and phase transitions in such systems.
Understanding how molecular size and structure affect the thermodynamic properties of the adsorbed layer is of fundamental importance in statistical physics. Furthermore, developing accurate theoretical descriptions of these systems has practical relevance in designing and optimizing separation processes in petrochemical applications. In this context, we present the first exact model for adsorption of molecular mixtures in zeolites that accounts for multisite occupancy. Using a rigorous statistical thermodynamic approach, we analyze mixtures of s-mers and k-mers (molecules occupying s and k lattice sites, respectively) on one-dimensional substrates. Our results reveal that the APR phenomenon emerges naturally from differences in molecular size (i.e., number of occupied sites) between the two species.
Following the line of Ref. [
51], we start calculating the parameter
, which is obtained from the equilibrium equations [Equations (
33) and (
34) in
Section 2.1.2]:
In order to understand the basic phenomenology, we consider in the first place a monomer-monomer mixture (
), with equimolar amounts of each kind of molecules in the gas phase (
). Under these conditions, Equation (
338) can be written as:
In this case,
const. Accordingly, the partial adsorption isotherms corresponding to the species
s and
k do not intersect and the APR phenomenon does not occur. This situation is reflected in
Figure 49, where an equimolar monomer-monomer mixture has been studied for
and
. The inset shows a similar analysis, where the values of the adsorption energies per particle were taken from Ref. [
51]:
(
J and
T=250 K) and
(
J and
T=250 K). In addition, the values of
and
were obtained by using Equation (
168) (
Section 2.1.2) with
(
) equal to the molecular mass of methane (ethane). Thus,
uma,
uma (where 1 uma
kg) [
135] and, consequently,
and
.
The results presented in
Figure 49 align with earlier studies on monomer–monomer mixtures [
51], where the occurrence of APR required introducing a complex set of lateral interactions.
We now turn to the more general case of an equimolar mixture composed of
s-mers and
k-mers. As indicated by Equation (
338), when both species have the same size (
), the individual adsorption isotherms do not intersect, and the resulting curves resemble those observed in
Figure 49.
However, when
, the system exhibits a broader range of behaviors that depend on both the value of the interaction parameter
A and the size relationship between the two species. Without loss of generality, we consider the case where
. In evaluating the role of
A, we rely on two key physical assumptions regarding the adsorption of linear alkanes: the heats of adsorption
are attractive and
increase (in absolute value) linearly with the length of the hydrocarbon chain. It is important to emphasize that these are not merely theoretical assumptions—they reflect well-established experimental findings on the energetics of alkane adsorption in zeolites, supported by extensive literature evidence [
117,
136,
137,
138].
Under these considerations, the argument of the exponential in Equation (
338) is negative,
A varies between 0 and 1 and there exists a value of coverage,
(
), at which the partial isotherms coincide (
) and, consequently, the APR phenomenon occurs. The value of
can be obtained from Equation (
338) by simple algebra:
The crossing point separates two adsorption regimes. Thus, for , is larger than . This tendency is reverted for , where is larger than . Then, the existence of the point is directly related to the displacement of the species k by the species s, and, consequently, to the presence of APR phenomenon.
This behavior is clearly illustrated in
Figure 50, which examines a monomer–dimer mixture (
and
) with equal concentrations of each species in the gas phase. The parameter values for
,
,
and
are the same as those given in the inset of
Figure 49. This monomer–dimer case is particularly insightful because:
it represents the simplest model of binary adsorption involving species of different sizes;
it captures the essential characteristics of multisite-occupancy adsorption; and
it serves as a useful model for interpreting experimental adsorption of methane–ethane mixtures in silicalite. Following a widely used approach, we adopt the bead-segment model, where each methyl group is represented by a single unit (or “bead”) equivalent in size to one adsorption site. Within this framework, methane and ethane are modeled as
and
, respectively.
Given that ethane molecules possess a stronger affinity for the surface (higher adsorption energy), they adsorb preferentially at low pressures. However, as pressure increases, adsorption of the smaller methane molecules becomes more favorable, eventually displacing the previously adsorbed ethane. As a result, the partial isotherms intersect at a critical coverage , indicating the occurrence of APR phenomenon. In this specific case, .
The exact results presented here offer a significant contribution to understanding the APR mechanism, demonstrating that it naturally arises from size disparity (or more precisely, the number of occupied sites) between different adsorbed species. It’s worth highlighting that in earlier theoretical treatments, such as the monomer–monomer model presented in Ref. [
51], reproducing APR required a highly parameterized model involving six variables (adsorption energies for methane and ethane, pairwise interaction energies among methane–methane, methane–ethane, and ethane–ethane, plus temperature). This complexity led to poorly constrained parameters and only qualitative conclusions. In fact, extremely small changes in parameter values (much smaller than typical experimental uncertainties) produced vastly different predictions. For example, a mere 5 percent change in ethane’s adsorption energy could determine whether APR appears or not in the model of Ref. [
51].
Supporting this point, Equation (
340) demonstrates that APR can emerge solely due to size asymmetry (
), even in the absence of any energetic contributions or lateral interactions (
and
). This reveals that introducing complex lateral interactions in previous models may be seen as an artificial or effective way to account for purely geometric or steric effects through energetic means.
In summary, this study leads to two key conclusions: the entropic effects arising from the non-spherical nature of polyatomic adsorbates play a crucial role in surface phenomena and cannot be neglected in comparison with monoatomic adsorption, and failing to incorporate the polyatomic nature of adsorbates into thermodynamic models can lead to a serious misinterpretation of experimental data, particularly regarding phenomena like APR.
7.6. Alkanes Adsorbed in Carbon Nanotubes Bundles: Surface Area Characterization
The multilayer adsorption theory for polyatomic molecules, discussed in
Section 2.1.3, is utilized here to evaluate the surface area of a HiPco single-walled carbon nanotube (SWCNT) sample. This analysis is based on adsorption isotherms obtained using a series of four alkanes—methane through butane.
The SWCNT sample analyzed in this study was procured from Carbon Nanotechnology Inc. (CNI). No post-processing was applied beyond evacuating the sample under vacuum (better than
Torr) for 72 hours before each measurement cycle. Although the purification method used by the manufacturer may have opened some nanotube ends, these openings were effectively sealed by chemical groups introduced during that same process [
139]. To unblock these caps, heating to at least 650 K under vacuum is necessary [
140,
141], a procedure not undertaken in this work. As a result, internal adsorption was excluded, and the sample effectively behaved as a non-porous substrate. Adsorption primarily occurred on the external surfaces of nanotube bundles and, to a lesser extent, on large-diameter, defect-related interstitial sites [
142,
143].
The isotherms were recorded using a custom-built volumetric adsorption setup [
144]. Low temperatures were achieved using a helium closed-cycle refrigeration system, and temperature regulation was maintained with two controllers. Pressure measurements were conducted using three capacitance manometers, with upper limits of 1, 10, and 1000 Torr, respectively. These gauges were housed in the room-temperature gas-handling system. Data acquisition and gas dosing were managed by a LabView-based program developed in-house. All gases employed were of ultra-high purity, sourced from Matheson Gas.
To determine the surface area, we calculated the product of the molecular area of the adsorbate [
145,
146,
147] and the monolayer capacity of the substrate. The monolayer capacity was derived from experimental adsorption isotherms using two approaches: the BET method [
54,
60] and the point B method [
148,
149]. The specific surface area was then calculated by dividing the total surface area by the sample’s mass.
The point B method involves plotting the adsorption isotherm as adsorbed amount versus pressure on a linear scale. Initially, coverage increases steeply with pressure, eventually reaching a pronounced inflection point. Beyond this point, coverage increases more gradually and linearly, until multilayer adsorption begins. Point B is identified as the lowest pressure at which a linear extrapolation aligns with this intermediate linear region between the completion of the first layer and the start of the second [
148,
149]. The coverage at point B corresponds to monolayer completion, as illustrated in
Figure 51. Generally, this method yields a slightly higher monolayer capacity than the BET equation.
Figure 52 shows the isotherms obtained for the HiPco nanotube sample using methane, ethane, propane, and butane. Different temperatures were used for each alkane due to their differing adsorption energies. Methane requires low temperatures where other alkanes exhibit vapor pressures too low for reliable measurement in our setup [
150,
151]. Conversely, the higher temperatures suitable for measuring butane adsorption produce pressures for methane that are too high to resolve the monolayer region effectively [
151].
To ensure comparability across different gases, isotherm temperatures were scaled by each adsorbate’s bulk critical temperature (), aligning the scaled temperatures across the data sets.
Table 3 compiles the molecular areas of the alkanes used. They were obtained from published neutron scattering results for alkane films adsorbed on graphite [
145,
146,
147].
Figure 53 presents plots of the linearized BET equation for the four different adsorbates used. The linearized BET equation is
In this context, P refers to the pressure at a specific point on the adsorption isotherm, is the saturated vapor pressure of the adsorbate at the isotherm temperature, n denotes the number of molecules adsorbed at pressure P, represents the monolayer capacity of the substrate, and C is a constant that reflects the interaction strength between the adsorbate and the substrate.
The BET equation is typically applied over a range of relative pressures from approximately
0.05 to
0.3 [
54,
60,
149]. This is the pressure range used for plotting the isotherm data in
Figure 52, which yielded the results shown in
Figure 53.
Figure 54 displays the specific surface area values for the HiPco nanotube sample, as determined using both the point B method and the BET equation for different adsorbates. Notably, both methods reveal the same trend: as the molecular length of the adsorbate increases, the specific surface area calculated for the same sample decreases. Within the range of molecule sizes examined, this decrease is approximately
for each method.
In an earlier study, we measured the specific surface area of a SWCNT sample using several spherical adsorbates (neon, argon, methane, and xenon) and found that the value was essentially independent of the size of the adsorbate [
152]. That analysis employed the point B method to determine monolayer capacity. These results differ markedly from the findings reported here for the linear alkane series (methane through butane). The key distinction lies in molecular geometry: unlike the spherical adsorbates used previously, the alkanes are linear, which enhances their tendency for multisite occupancy during adsorption.
We now proceed to analyze the experimental data using the theoretical framework for multilayer adsorption of polyatomic molecules, as developed in
Section 2.1.3 and
Section 3.4. The model outlined in
Section 2.1.3 will be referred to hereafter as the "modified BET model," or simply the MBET model. As previously mentioned, in the one-dimensional (1-D) case, the MBET equation reduces to the standard BET equation when considering monomer adsorption (
). The explicit expression derived for the adsorption of dimers (
) in 1-D within the MBET framework is given by:
The variables appearing in Equation (
342) carry the same meaning as those defined in Equation (
341). Equation (
342) is applicable for determining the monolayer capacity and specific surface area in the case of ethane adsorption, where ethane molecules can be effectively treated as dimers. Unlike the standard BET equation, the MBET expression for dimers is nonlinear. When applied to ethane, the MBET model yields a larger surface area than the value obtained using the BET equation on the same dataset. Notably, the MBET model for dimers involves only two adjustable parameters—the same number as in the BET formulation.
As previously mentioned, for systems with in two dimensions (2-D), no exact analytical expressions exist for adsorption isotherms, i.e., for the fractional coverage n as a function of the relative pressure . To address this challenge, we adopted two different strategies to extract monolayer capacity from experimental data:
- (1)
-
In this first approach (procedure A), we used the one-dimensional MBET equations for all four adsorbates to account for the linear geometry of the molecules. For methane, we employed the standard BET equation, which corresponds to the MBET expression for monomers. For ethane, we applied the exact 1-D MBET formula for dimers [Equation (
342)]. For propane and butane, we utilized the same dimer equation but adjusted its parameters to fit the experimental data in the low-pressure, low-coverage regime—specifically, the same region typically used in BET analysis.
Using Procedure A, the calculated specific surface areas for ethane, propane, and butane were consistently higher than those obtained using the BET method. In addition to producing improved results, this method remains relatively straightforward to implement.
However, we also observed a consistent trend: as the length of the alkane chain increases, the derived specific surface area decreases. This behavior is illustrated in
Figure 6, which includes data obtained from both this and other approaches.
Although not perfect, the results from Procedure A represent a clear enhancement over the standard BET approach.
- 2)
-
The second strategy (procedure B) involved fitting the isotherm data for all four adsorbates to the approximate MBET expression developed for the two-dimensional case [Equation (
209) in
Section 3.4]
In this formulation, is a constant that reflects both the interaction strength between the adsorbate and the substrate, as well as the connectivity of the adsorption lattice. The parameter k denotes the number of units in the k-mer molecule.
The experimental adsorption isotherms are fitted to this model using the appropriate
k value for each adsorbate—1 for methane, 2 for ethane, 3 for propane, and 4 for butane—within the same low-pressure range typically employed for BET analysis.
Figure 55 shows the fit obtained for butane, which demonstrates excellent agreement with the data; similarly accurate fits were achieved for the other three alkanes. The fits yield values for the monolayer capacity,
, which are subsequently used to compute the specific surface area of the sample.
Figure 56 presents the key findings of this study. It shows the specific surface area values of the substrate calculated using the four methods discussed: the point B method, the BET equation, the 1-D MBET model (Procedure A), and the 2-D MBET model (Procedure B). These values are plotted against the number of carbon atoms in the alkane adsorbates used in the isotherm measurements. Among the four approaches, only the 2-D MBET model (Procedure B) produces surface area values that remain nearly constant across the range of adsorbates. In contrast, the other three methods show a decreasing trend in surface area with increasing alkane chain length.
This result has practical significance, as it offers a consistent approach for determining surface area using linear adsorbates. Specifically, the 2-D MBET model allows for surface area measurements with linear molecules that closely mirror the consistency typically achieved when using spherical adsorbates with the BET or point B methods.
From a more fundamental standpoint, our findings underscore the significance of the additional entropy contribution that emerges in monolayer films of linear molecules due to multisite occupancy.
To summarize, in this section we presented adsorption isotherm measurements for a series of four alkanes (methane, ethane, propane, and butane). These data were then used to estimate the specific surface area of a single substrate using four different methodologies: the traditional BET equation, the point B method, and two versions (1-D and 2-D) of a more recent model for the adsorption of linear molecules, known as the MBET method. Our key conclusion is that the 2-D MBET approach (procedure B) yields consistent surface area values across different linear adsorbates, in contrast to the other methods evaluated (point B, BET, or the 1-D MBET model), which do not show such consistency.
This result holds both practical and theoretical significance. Practically, it offers a reliable method for determining specific surface area using longer linear molecules, producing values that align well with those obtained from spherical adsorbates using the BET equation. This, in turn, enhances the comparability of results obtained for the same substrate across different types of adsorbates.
7.7. Crystal Growth from Aqueous Solution in the Presence of Structured Impurities
The process of crystal growth plays a vital role in both biological systems and industrial technologies. It begins with the formation of stable nuclei, followed by the gradual addition of atoms or molecules onto the crystal’s surface. These growth units—comprising atoms or molecules identical to those in the crystal—migrate from the bulk solution and attach to specific surface sites, eventually leading to the formation of macroscopic structures with well-defined surfaces.
It is well established that the adsorption of foreign atomic or molecular species, known as impurities, can significantly influence the rate of crystal growth [
153,
154,
155]. This phenomenon has been widely utilized to control crystal morphology and to enhance the properties of crystalline materials, powders, and granulated substances. Understanding the mechanisms by which impurities exert their effects is therefore of considerable importance.
Numerous studies have explored the influence of impurities on both the growth rate and morphology of single crystals [
156,
157,
158,
159,
160,
161]. Theoretical models addressing this influence typically assume that impurities (ions, atoms, or molecules) adsorb at specific surface features such as kinks, steps, and terraces during crystal growth. Early theoretical efforts in this area were made by Bliznakov [
162,
163,
164,
165], who focused on impurity adsorption at step edges, and by Cabrera and Vermilyea [
166], who considered adsorption on flat surface terraces. Experimental data on crystal growth rates have supported the predictions of these models [
160].
In a key contribution, Davey and Mullin (DM) [
167] developed a theoretical framework to describe how impurities affect the growth rate of crystals from aqueous solutions. Their model posits that impurity adsorption is limited to a thin layer adjacent to the crystal surface, and assumes no lateral or vertical interaction between adsorbed impurity molecules. Within this context, the velocity of step movement,
V in the presence of impurities is given by:
where
is the step velocity in the absence of impurities, and
denotes the fraction of surface sites occupied by impurities. Consequently,
reflects the proportion of available sites for the adsorption of pure growth units. As
approaches 1, the surface becomes fully covered with impurity species, and crystal growth effectively halts.
The ability of an impurity to inhibit crystal growth is also believed to depend on factors such as its size, shape, and orientation—a phenomenon referred to as stereochemical effect. To incorporate these aspects, Kubota and Mullin (KM) [
168] proposed an improved kinetic model. Their approach considers impurity adsorption along step edges and introduces a parameter,
, representing the efficiency of impurity blockage. The step velocity in this model is given by:
This equation shows that the overall impact of impurities on growth rate is governed by two factors: the surface coverage by impurities (
) and the effectiveness of each adsorbed molecule (
) in hindering growth.
The relationship between impurity concentration (
) and surface coverage (
) is typically described using adsorption isotherms. The Langmuir isotherm [
169] (or its extended versions [
170]) is often used due to its mathematical simplicity. Then,
where
is the Langmuir constant. Introducing Equation (
346) in Equations (
344) and (
345), the relative step velocity can be written in terms of the impurity concentration:
and
The models in Equations (
347) and (
348) assume that each impurity molecule occupies a single site on the surface and that the spatial arrangement of those sites is irrelevant. As a result, they cannot differentiate between varying impurity structures or lattice geometries. This limitation has led to the development of more advanced theoretical approaches that account for the equilibrium adsorption of polyatomic species with diverse molecular configurations.
In this context, we aim to present a new generalization of the DM model [
167] for crystal growth rates that incorporates both the size and shape of impurities, as well as the surface geometry. As in the KM model, the theoretical equations are based on the pinning mechanism of Cabrera and Vermilyea for the inhibition of step advancement [
166]. Building upon the model discussed in
Section 3.2.1, we employ the
adsorption isotherm [Equation (
120)] to describe the surface coverage of impurities
.
Before starting the comparison, and in order to analyze adsorption from liquid solutions, it is convenient to write the theoretical adsorption isotherm given in
Section 3.2.1 in a more appropriate form. By using the equilibrium condition
, Equation (
120) adopts the form
where
is the chemical potential of an ideal solution. Hence,
where
is the standard chemical potential, which is defined as
where
h is the Planck constant (
) and
m is the mass of the
k-mer. Taking into account that the mass of a
k-mer is directly proportional to the mass of the monomer unit:
with
the mass of the monomer unit. Then, Equation (
350) can be rewritten as:
and
where
is the impurity concentration,
and
. This expression relates the chemical potential in solution with the impurity sizes and the concentration in the bulk. Introducing Equation (
353) in Equation (), the following expression is obtained,
where
.
Equation (
118) is valid for
(
for
).
Now, replacing in Equation (
344)
by
, the relative step velocity
for impurity molecules of different shape and size can be obtained. For this purpose, it is necessary first to rewrite Equation (
354) in the form
(
as a function of
). This procedure was performed for linear impurities [
for
, and
for
] with
k ranging between
and
adsorbed on square lattices (
)
16. The resulting expressions are as follows:
with
. In addition,
where
where
.
The model described in Equations (
355–
358) represents an advancement over the previous KM model [
168], in which the effects of impurity size and shape were incorporated through a single proportionality constant, or effectiveness factor,
(
). In contrast to the KM approach, the new expressions for
allow us to (1) explicitly capture the entropic contributions of non-spherical impurities, and (2) investigate the influence of the spatial arrangement of adsorbed species on crystal growth, using parameters with clear physical significance. These parameters can be derived from thermodynamic measurements and are directly related to the spatial configuration of the adsorbed impurity molecules. In the following paragraphs, our theoretical predictions will be compared with Monte Carlo (MC) simulations and applied to interpret experimental data on the growth rates of the 100 crystal faces of KBr as a function of impurity concentration.
We will begin by discussing the behaviour of the relative step velocity
as a function of the concentration
for different values of the impurity size. Thus,
Figure 57 shows
for two cases:
, solid line, Equation (
355); and
, dashed line, Equation (
358). In both cases,
.
As shown in
Figure 57, the curve corresponding to
exhibits a steeper decline at low impurity concentrations and remains consistently below the
curve throughout the low- and medium-concentration regimes. As anticipated, at high impurity concentrations
, the curves converge and the step velocity approaches zero—though this regime is omitted from the figure for clarity. These results demonstrate that the model defined by Equations (
355–
358) effectively captures the essential behavior of the system. Specifically, in the low to intermediate concentration range, the ability of an impurity to suppress crystal growth is strongly influenced by its size and shape, with larger or more complex impurities exhibiting greater effectiveness. However, as
increases and the surface becomes nearly saturated with impurities, the role of size and shape becomes less significant, and the growth suppression effect levels off.
To further assess the applicability of Equations (
355–
358), we now examine an experimental case study, as illustrated in
Figure 58. This analysis provides an opportunity to evaluate both the practical relevance and the predictive power of the theoretical framework proposed in this work.
The experimental data are taken from the study by Bliznakov and Nikolaeva [
171], in which the impact of aliphatic carboxylic acid impurities on the relative growth rate of the 100 faces of KBr crystals was investigated. Their findings revealed two key observations: (1) the relative growth rate of the crystal face decreases with increasing impurity concentration, and (2) the extent of this decrease depends on the size of the impurity, as determined by the number of carbon atoms in the carboxylic acid.
Subsequently, Kubota and Mullin [
168] modeled these experimental results using their kinetic scheme [Equation (
345)]. The KM model provided a good fit to the data, using the parameters
and
listed in
Table 4.
Figure 58 presents this comparison: dashed lines represent the theoretical predictions of the KM model, while the experimental data are indicated by symbols (squares for HCOOH, triangles for CH
3COOH, diamonds for C
2H
5COOH, and circles for C
3H
7COOH). Kubota and Mullin concluded that the relatively constant values of
suggest that adsorption occurs primarily through the carboxyl group. Regarding the effectiveness factor
, they proposed that its increase with the number of carbon atoms could be attributed to the growing size or bulkiness of the impurity molecules.
We now turn to applying our multisite adsorption model to the experimental data reported in Ref. [
171], with the aim of offering an alternative perspective on the effects of impurity size. Leveraging the inclusion of the impurity size parameter
k in Equations (
355–
358), we assign specific values of
k based on the number of carbon atoms in the carboxylic acid molecules: HCOOH corresponds to
, CH
3COOH to
, C
2H
5COOH to
, and C
3H
7COOH to
.
It is important to emphasize that this assignment of
to 4 is a practical method to account for differences in impurity size and does not imply any specific assumptions about the adsorption mechanism.
17 Furthermore, this approach enables a fully analytical treatment of the problem via the explicit expressions provided in Equations (
355)–(
358). An alternative strategy would be to treat the size parameter
k as a fitting variable. However, the absence of a general analytical expression for
makes such an analysis significantly more challenging.
On the other hand, the values of were obtained using a standard least-squares fitting procedure for each experimental isotherm. Four cases were considered:
Case HCOOH: Equation (
355) was used as fitting function, with
and
as adjustable parameter.
Case CH
3COOH: Equation (
356) was used as fitting function, with
and
as adjustable parameter.
Case C
2H
5COOH: Equation (
357) was used as fitting function, with
and
as adjustable parameter.
Case C
3H
7COOH: Equation (
358) was used as fitting function, with
and
as adjustable parameter.
The results are shown in
Figure 58 (solid lines) and the fitting values of
are collected in
Table 4. Clearly, Equations (
355–
358) provide a very good agreement with the experimental data.
The values of
obtained (see the fifth column in
Table 4) are notably consistent across the different impurities, supporting the conclusion that adsorption occurs primarily through the carboxylic functional groups. Additionally, the analysis based on Equations (
355–
358) strengthens the interpretation proposed by Kubota and Mullin [
168]. Specifically, although the carboxylic acids are unlikely to lie flat on the surface, variations in impurity size (i.e., steric effects) have a significant impact on the behavior of the relative step velocity,
. As anticipated in Ref. [
168], the effectiveness of growth inhibition increases with the molecular size of the impurity.
In summary, this section has introduced a novel theoretical framework for analyzing crystal growth from aqueous solutions in the presence of structured impurities. The model is built upon physically grounded equations, involving parameters with clear physical interpretations. These parameters are experimentally accessible via thermodynamic measurements and directly reflect the spatial configuration of the impurity molecules in the adsorbed state.
Theoretical predictions have been validated against Monte Carlo simulations and successfully applied to describe experimental data on the relative growth rates of the 100 faces of KBr crystals in the presence of various aliphatic carboxylic acids. The findings confirm and further support previous conclusions regarding the critical role of impurity size in modulating crystal growth behavior [
168].
7.8. Application to k-Mers Phase Transitions
In this section, we apply the multiple exclusion (
) statistics to the problem of linear
k-mers adsorbed on a square lattice with
M sites. In this system, two phase transitions have been observed for
[
96]: (1) an entropy-driven isotropic-to-nematic (I-N) continuous phase transition at intermediate coverage [
96,
101,
104], and (2) a nematic-to-isotropic (N-I) discontinuous transition at densities close to lattice completion [
106].
It is important to highlight that for large k-mers ( in the case of square lattices), spontaneous orientational order emerges from purely entropic effects, without the need for explicit attractive interactions between the rods. This behavior exemplifies how entropy alone can induce phase transitions in systems constrained by hard-core exclusion.
The statistics provide a natural and powerful framework to model the thermodynamics of these transitions, as they inherently incorporate the excluded volume effects and the reduced configurational entropy associated with rod alignment at higher densities.
7.8.1. Basic Definitions
As discussed in
Section 6, the
k-mers problem on the square lattice is modeled as a mixture of two species, each aligned along one of the two lattice directions: horizontal (
) and vertical (
), referred to as
-mers and
-mers, respectively. We assign
-mers
and
-mers
. Both species occupy
k consecutive sites along their respective axes.
According to the definitions established earlier, the following identifications hold: , , , , , , , , , and . Furthermore, the saturation densities are and , with all exclusion coefficients equal: . It is important to note that, since and , the states available to each species are restricted to those lying along their characteristic direction.
Regarding the saturation occupation numbers
and
satisfying
(according to Equation (
313)), a fully covered lattice must fulfill
18 Given a vertical occupation number
for
-mers, the maximum horizontal occupation number is
and analogously,
For
k-mers on the square lattice, since
, these simplify to
The self-exclusion per particle for each species at infinite dilution corresponds to the number of states excluded by an isolated
k-mer along its axis:
identical to the 1D case. From Equation (
324), the solutions yield
.
The number of cross-excluded states between species 1 and 2 is
The first term,
, accounts for the initial total exclusion of perpendicular states by an isolated
k-mer, while the second term,
, corrects for the non-independent states between directions. Two isolated
k-mers of species 1 and 2 exclude together
states out of the
total cross-states. By symmetry,
, and thus
This follows formally from the general relation in Equation (
323), where
. Therefore,
.
From the solutions to Equations (
324) and (
326) with
, the exclusion correlation parameters are
and
where
denotes the Lambert function, taking its main branch for
or the lower branch otherwise.
The numerical solutions for are: for (), for (), for (), and for ().
7.8.2. Entropy Surface, Equilibrium Path, and Order Parameter
The Helmholtz free energy per site, , can be fully represented in the space . The relation between the mean occupation number n, which denotes the state occupation number of the whole lattice, and the occupation numbers of each species on their respective sublattices is given by and with , so that . The maximum value that n can reach is , corresponding to a line in the plane such that , named the saturation line.
For each value of n, the function is minimized numerically, and the points where the Helmholtz free energy attains local minima are obtained.
Since the model is athermal, only excluded-volume interactions are present and
, leading to a vanishing internal energy. Therefore, the minima of
correspond to maxima of the entropy per site in units of
. An order parameter
is defined based on the difference between the species occupation numbers as
Here, indicates isotropic k-mer distribution, while reflects nematic ordering.
For , all solutions satisfy for any n, implying and absence of phase separation: both species distribute equally.
In contrast, for , phase separation occurs beyond a critical density . A nematic phase emerges, where one direction becomes favored. The particular case of deserves a detailed discussion.
Although the model predicts phase separation at very high coverage for
(see
Figure 59), this behavior results from the Helmholtz free energy vanishing at saturation, i.e.,
and
, according to Equations (
306) and (
307), owing to the restrictive boundary conditions
imposed at saturation. It is worth noticing, that also for
the order parameter do not show any transition to an isotropic high-coverage regime as expected from MC simulations. However, if the entropy reaches a finite value
at full coverage instead of vanishing, this anomalous phase separation for
would not occur and, as shown latter, the high-coverage nematic-to-isotropic transition does occur for
.
By adopting either a more general form for the density of states (as introduced in
Section 7.8.3) or by adding an
ad hoc high-density correction to the entropy surface, this spurious behavior is eliminated. Consequently, the system does not exhibit phase transitions for
, consistent with prior single-species analyses of nematic transitions discussed earlier in
Section 5. For
, the entropy at high coverage is exceedingly close to, but still higher than, that of the nematic regime. Thus, a highly accurate approximation is required to capture the MC simulation result that nematic ordering only occurs for
. The subtle dependence of entropy at high coverage critically determines the absence of nematic transition for
and its appearance for
[
101,
172].
For
, there exists a critical state occupation
such that for
, phase separation (
) occurs and
k-mers preferentially align along a lattice direction, forming a nematic phase. Even under the restrictive condition of vanishing entropy at saturation, the model predicts, for
,
,
, a critical state occupation
, corresponding to
, and for
,
, a critical value
,
. These values are very close to the known MC simulation results
,
[
101], and also agree with earlier estimates reviewed in
Section 5 [
18].
Remarkably, the complex many-body correlations are captured by a single cross-exclusion parameter
, yielding surprisingly accurate predictions. The input boundary values that best match MC results,
(
), are close to, but slightly lower than, the analytical solution value
(
) obtained from Equation (
326). This highlights the already significant accuracy of the first-order approximation assuming constant
with vanishing entropy constraints, although it also suggests that
may vary slowly with
n, as proposed in the second-order approximation of
statistics discussed in the previous section.
Thus, the critical density predicted by the model is highly sensitive to the cross-exclusion parameter
. In this elaborated approach compared to one presented in
Section 5, the first-order approximation is kept, setting
constant, since many essential conclusions can already be drawn from it. Overall, the present formalism reveals that the emergence of nematic order for
is mainly driven by cross-excluded volume effects encoded in
, with critical values finely tuned by its magnitude.
To better match the observed behavior at high coverage, an empirical entropy contribution was introduced to the theoretical entropy derived previously. This correction accounts for the finite entropy observed in simulations even at full lattice coverage, as the lattice exhibit local rearrangements of k-mers patches.
The proposed empirical correction is given by:
where
Here, is the mean occupation, the maximum occupation, and the saturation entropy at full coverage. The parameters , , and control the behavior near saturation. Specifically, follows with fitted from MC simulations: for to respectively, and for .
The correction satisfies: , for isotropic full coverage, and for fully aligned k-mers at saturation (either , or vice versa), thus recovering the one-dimensional limit.
Consequently, the corrected entropy surface
becomes:
with
, again, being the entropy surface vanishing at the saturation line
. This correction becomes significant at high occupations (
), especially for small
values.
Values , –, and – were adjusted to reproduce the chemical potentials and transition coverages observed in simulations, as discussed later.
Figure 60 shows the resulting order parameter
as a function of coverage
for
, 8, 12, and 14, when the high-coverage empirical entropy correction
is included. The high coverage entropy contribution not only predicts the isotropic-nematic transition at intermediate densities but also captures the nematic-to-isotropic transition at high densities, as observed in MC simulations.
The critical densities for the isotropic-nematic transition obtained from the model are very close to simulation results. For instance, for
, the model predicts
(
), and for
,
(
). These values are in excellent agreement with MC simulations reported in Refs. [
101,
104].
Furthermore, the model predicts a sharp drop in the order parameter for at (corresponding to ), associated with the nematic-to-isotropic transition. Predictions for and are also provided, yielding and , respectively. Although simulation data for these higher k-values are scarce, the model suggests that the critical coverages should be slightly lower, as already observed for smaller k.
Although the curves in
Figure 60 appear to show continuous transitions at high density, a more detailed analysis (discussed later in
Section 7.8.4) reveals that the nematic-to-isotropic transition at high coverage is indeed of first order, as indicated by a discontinuous jump in the chemical potential.
Additionally, the full Helmholtz free energy surface
can be derived analytically.
Figure 61 displays this surface for
. The equilibrium states follow the black solid curve shown over the surface. For
, the equilibrium corresponds to isotropic configurations with
. As
n increases beyond
, two symmetric branches emerge corresponding to nematic ordering along either the horizontal or vertical direction. These branches are symmetric with respect to the isotropic bisector
and represent coexistence between a high-density aligned phase and a low-density dilute phase.
Specifically, for , two minima of exist at points and with , symmetrically located with respect to the line . If , the high-density phase is oriented along the vertical direction; otherwise, it is aligned horizontally. This analytical description of the density branches in the nematic phase constitutes an important result of the present formalism.
The Helmholtz free energy difference between the isotropic and nematic phases is very small, indicating that the nematic phase is only weakly more stable. It is also highly sensitive to perturbations, such as weak additional interactions, which could further stabilize or destabilize the nematic phase depending on their nature.
The discussed results and predictions remain valid even when the empirical entropy correction is not explicitly included. As shown later, the general formulation of statistics allows one to obtain equivalent behavior by simply relaxing the constraint of zero entropy at saturation, setting nonzero values for the density of states at full coverage.
Thus, the empirical entropy correction serves primarily as a practical and illustrative tool for matching MC data and providing physical insight, while the general formalism ensures the robustness of the theoretical predictions.
7.8.3. Generalized Density of States Function in Multiple Exclusion Statistics
Beyond the analysis of the high-coverage phase transitions conducted in the previous section by introducing an empirical entropy correction in the
plane through Equations (
370)–(
371), an analogous behavior naturally arises from the
statistics by adopting a less restrictive boundary condition for the density of states functions in Equation (
313). Instead of enforcing
, corresponding to a fully oriented phase at saturation (as in the previous approximation), we now consider
.
Indeed, the strict condition is valid only when the lattice is saturated entirely by -mers or -mers, with occupations or , respectively. However, for other saturation states along the line , where both types of k-mers coexist, the density of states at saturation is expected to remain finite.
To demonstrate that similar results regarding the order parameter and phase transitions emerge from this generalization, we adopt the following expressions for the density of states functions:
where
and
denote the finite density of states per particle for species 1 and 2 at the saturation line, corresponding respectively to
and
.
It should be notice that
increase smoothly along the saturation line from zero for a fully oriented phase,
or
, up to
for an isotropic phase
. Moreover,
and
are made to vary according to a functional form analogous to Equation (
370) along the saturation line, being the first term
.
where the equations are meant to be valid for
along the saturation line. Thus, given
,
is fixed or, conversely, given
, then,
.
The values
for an isotropic saturated phase takes values:
for
to
, respectively, as introduced in [
18] of this series to match the MC values for the entropy at full coverage [
90,
91,
107].
Thus, the entropy functions from Equation (
372) and the one from
statistics, Equation (
307), for the general density of states forms, Equations (
373) and (
374) evaluated along the saturation line cases differ each other by less than
of
, except at the fully aligned or fully isotropic limits where they take the identical desired values. Ultimately, the two free energy surfaces involved in the calculations compared in this section take very approximately the same values all along the saturation line in the plane
.
As we will discuss in the next section, adopting the generalized density of states functions from Equation (
373) leads to an entropy surface that captures both qualitatively and quantitatively the phase transitions observed in MC simulations. It provides an analytical description of the nematic branches, characterizes the order of the transitions, and matches the first-order nature of the high-density transition recently confirmed for
(and very possibly for
) in Ref. [
106].
7.8.4. Nematic-Phase Density Branches and Phase Transitions
Regarding adsorption isotherms, this analytical treatment provides both the low- and high-density branches, which to our knowledge had not been obtained previously from a purely analytical approach for the k-mers model on the square lattice. These results are compared with corresponding branches obtained from fast-relaxation MC simulations.
As in previous
Section 5.5, the MC simulations were conducted in the grand canonical ensemble using the algorithm developed by Kundu et al. [
101,
111,
112], which efficiently handles the sampling slowdown at high densities. For a fixed value of
, the number of
k-mers on the lattice evolves via non-local insertion and deletion of
k-mers. Briefly, each MC step consists of removing all horizontal
k-mers while keeping vertical ones fixed. The probability of forming an empty horizontal segment of a given length is calculated exactly and stored, and empty segments are subsequently filled with
k-mers according to these probabilities. An analogous procedure is applied in the vertical direction. The insertion-deletion moves satisfy detailed balance and the algorithm is ergodic. Simulations were performed on
square lattices with
and periodic boundary conditions, showing negligible finite size effects. Equilibrium was typically reached after
MC steps. The MC results for the density branches, displayed in
Figure 62 and
Figure 63, were obtained by averaging the coverage, or state occupation number, over
MC steps, regardless of the orientation of
k-mers in each configuration.
The analytical densities for the high- and low-density phases,
and
, were obtained as functions of the chemical potential of the gas-phase,
, by solving for
and
the pair of coupled equations
[Equations (
314) and (
315)], assuming equilibrium between the two
k-mer species and the
k-mer gas at reservoir. The mean occupation and lattice coverage are given by
and
, where
and
represent the coverage along the horizontal and vertical directions, respectively.
Analytical results are compared to MC simulations for
and
in
Figure 62 and
Figure 63. Two theoretical approximations are shown: (1) the results obtained by adding the empirical high-density entropy correction , Equation (
370), to the basic entropy surface corresponding to density of states functions vanishing at full coverage [Equation (
313)] resulting in entropy surface of Equation (
372), displayed as green lines; and (2) the results from the entropy entropy surface corresponding to generalized density of states functions [Equation (
373)], shown as blue lines. As seen, both approaches yield very similar results. It worth noticing that, we many times refer to entropy surfaces instead of Helmholtz Free Energy surface given that
in this problem because
.
The chemical potential versus density curves (adsorption isotherms) are very well reproduced up to saturation. At low coverage, an isotropic phase prevails; as density increases, the nematic branches emerge, and finally a transition to a disordered isotropic phase occurs. The nematic branches and coexistence region are qualitatively well captured.
The isotropic–nematic transition is clear for
(
Figure 63) and
(
Figure 62), with critical points approximately at
,
(
) and
,
(
). The transition is continuous, as evidenced by the continuous chemical potential dependence, matching MC results.
No transition occurs for
(
Figure 64). In that case, theoretical and MC results show that
for all densities, confirming the absence of ordering. The inset shows
as an example, where the branching is much more pronounced.
At high coverage, a nematic-to-isotropic transition is predicted. For
, the branches collapse around
, corresponding to a coverage jump
, indicating a first-order transition. The critical chemical potential estimated by Maxwell construction is
, close to the prediction
from Ref. [
106]. While the predicted coverage jump (
) is smaller than the value
reported by Shah et al. [
106], the trend is consistent.
For , the behavior is more sensitive: a small discontinuous jump appears, depending critically on . Analytical predictions give () for , and () for .
Overall, while a more detailed finite-size scaling analysis is needed from MC simulations, these results suggest a first-order nematic-to-isotropic transition for , consistent with MC results and previous theoretical work, and also fore although less conclusively.
7.8.5. State Exclusion Spectrum Functions of k-Mers: Coverage Dependence
The self-exclusion and cross-exclusion spectrum functions, denoted and , respectively, formally related to the density dependence of the chemical potential, offer a novel characterization of statistical exclusions arising from spatial state correlations. These functions provide insight into the statistical behavior of the system during the isotropic–nematic phase transition and throughout the phase coexistence. Furthermore, they open the possibility of experimentally probing particle spatial ordering through purely thermodynamic measurements, which is not easily inferred from standard adsorption isotherms, as we can already observe from the model k-mers isotherms.
We now focus on the statistical characterization of phase behavior in the
k-mers system through the self-exclusion and cross-exclusion per particle frequency functions,
,
,
, and
, as well as the cumulative average state exclusion per particle functions,
,
,
, and
, introduced in
Section 6.3 through Equations (
321) and (
322).
For , all exclusion functions decrease monotonically with coverage from their infinite dilution limits, and , down to zero. Moreover, and at all densities, indicating an isotropic phase where the average number of self-excluded and cross-excludes states by a - or -oriented k-mer are identical for each species. As a result, the average number of states excluded per particle decreases with increasing coverage, although the decrease is not linear due to the statistical effects.
For , the exclusion functions , , , and also decrease from up to the critical coverage , maintaining and . At , they split into two distinct branches, remaining continuous but with a discontinuous first derivative. If we assume, for simplicity, that the high-density nematic phase aligns along the direction (species 2), then represents the self-exclusion along per -oriented particle, while corresponds to the same along for -oriented particles. Similarly, refers to cross-exclusion in by a -oriented particle, and is the reciprocal.
Figure 65.
Self- and cross-exclusion frequency functions: red (
)
(lower branch) and
(upper branch); blue (
)
(lower) and
(upper); black (
)
and
; green (
)
and
. Analytical results from Equation (
322) with density of states from Equation (
373).
Figure 65.
Self- and cross-exclusion frequency functions: red (
)
(lower branch) and
(upper branch); blue (
)
(lower) and
(upper); black (
)
and
; green (
)
and
. Analytical results from Equation (
322) with density of states from Equation (
373).
After the transition at , the behavior is as follows: the decreasing branch of indicates stronger alignment and denser packing along ; the decreasing reflects compact transversal packing; the branch of increases sharply beyond , indicating more dispersed -oriented particles; and reaches a maximum and then decreases slightly as density increases.
Since and for , the -oriented phase is always more aligned and compact than the one.
While all
vanish at full coverage, the cumulative exclusion functions
, defined in
Section 6.3, are highly sensitive to changes around the transitions. They split at
and show a discontinuous jump at the high-density transition, consistent with its first-order character.
Figure 66.
Average exclusion per particle spectrum functions for and . Red and black: ; blue and green: . Analytical results from statistics formalism.
Figure 66.
Average exclusion per particle spectrum functions for and . Red and black: ; blue and green: . Analytical results from statistics formalism.
Thus, the functions provide sensitive and valuable analytical tools to trace phase transitions and characterize lattice configurations. Importantly, and can be easily computed in MC simulations.
Regarding the cross-exclusion parameters
, they quantify the strength of statistical interaction between
- and
-oriented
k-mers. The leading term
shows that
controls the exponential decay of available states for
particles as
particles fill the lattice. Changing variables to
and
in the isotropic phase, we have
thus
where for . For , , thus no nematic transition is expected.
Qualitatively, in terms of state exclusion statistics, it is the strong depletion of available states for a given orientation, driven by the density of particles in the perpendicular direction, that induces reordering into a nematic phase coexisting with a less dense transverse phase.
Summarizing, the exclusion statistics formulation introduced in
Section 6 was applied to address the challenging problem of
k-mers on the square lattice due to existence of at least two phase transitions for large
. The system was modeled statistically as a mixture of two self-excluding and cross-excluding differently oriented species, or, in equivalent words, as identical particles occupying two different set of states observing self and cross-state exclusion between them because of their topological correlations. The analytical
statistics predicts the a continuum isotropic-nematic transition at intermediate coverage only for
and, particularly significant, a first-order nematic-isotropic transition at high coverage only for
, and less clearly for
as well although not as conclusively. The theoretical model reproduces qualitatively and quantitatively well the density dependence of the chemical potential, the lattice coverage and chemical potential at the transitions, it predicts the density branches at the nematic regime as it gives the phase-coexistence line. The formalism provides an accurate approximation for the whole Helmholtz free energy and entropy surfaces in terms of the state occupation numbers,
and
, of horizontal and vertical
k-mers, from which the equilibrium paths are obtained by minimization at fixed average lattice occupation
n. The low and high density branches of the adsorption isotherms are accurately reproduced compared to MC simulations. A unified framework was developed by introducing proper thermodynamic averaged state exclusion functions and by expressing them in terms of the
k-mers chemical potential and density in order to determine the fundamental statistical correlation parameters from their boundary properties.
8. Monte Carlo Simulation Method Applied to the Problem of Adsorption with Multisite Occupancy
In the following sections, we present different computational algorithms of interest for the study of adsorption problems involving multiple-site occupancy. Most of these algorithms have been used throughout the present work.
8.1. Metropolis MC Algorithms for Adsorption of Interacting k-Mers
8.1.1. Grand Canonical Ensemble
In order to simulate the adsorption/desorption process of
k-mers in the grand canonical ensemble, we use a generalized MC algorithm based on the Metropolis scheme of transition probabilities [
173]. The method starts with a system of
M sites (and connectivity
) at temperature
T and pressure
P (or chemical potential
). The simulation process is a repetition of the following elementary steps (MCS):
- i)
set the value of the chemical potential and the temperature T;
- ii)
choose randomly a linear k-uple of nearest-neighbor sites.
- iii)
if the
k sites selected in
are empty, an attempt is made to deposit a rod with probability
; if the
k sites selected in
are occupied by units belonging to the same
k-mer, an attempt is made to desorb this
k-mer with probability
and otherwise, the attempt is rejected. Here,
and
represent the probabilities of transition from a state with
N particles to a new state with
or
particles, respectively. Following the Metropolis scheme [
174], these probabilities are given by
, where
is the difference between the Hamiltonians of the final and initial states.
- iv)
repeat steps M times.
The first MCS of each run are discarded to allow for equilibrium and the next m MCS are used to compute averages. In the low-temperature regime, where ordered phases are expected to develop, displacement (diffusional relaxation) of adparticles to nearest-neighbor positions, by either jumps along the k-mer axis or reptation by rotation around the k-mer end, must be allowed in order to reach equilibrium in a reasonable time.
8.1.2. Canonical Ensemble
In the canonical ensemble, the thermodynamic equilibrium is reached by following Kawasaki’s dynamics [
173], generalized to deal with polyatomic molecules. The algorithm to carry out an elementary Monte Carlo step (MCS), is the following:
Given a square lattice of M equivalent adsorption sites:
- i)
set the value of the temperature T;
- ii)
set the value of the coverage, , by adsorbing linear molecules onto the lattice, each molecule occupying k adsorption sites;
- iii)
a
k-mer and a linear
k-uple of empty sites are randomly selected, and their positions are established. Then, an attempt is made to interchange its occupancy state with probability given by the Metropolis rule [
174]:
where
is the difference between the Hamiltonians of the final and initial states;
- iv)
a k-mer is randomly selected. Then, a displacement to nearest-neighbor positions is attempted (following the Metropolis scheme), by either jumps along the k-mer axis or reptation by rotation around a unity of the k-mer. This procedure (diffusional relaxation) must be allowed in order to reach equilibrium in a reasonable time; and
- v)
repeat steps iii)-iv) M times.
As in the case of grand canonical MC simulations presented in
Section 8.1.1, the equilibrium state was reached after discarding
MCS, and the averages were taken over the next
m MCS.
8.2. Parallel Tempering MC Algorithm for Adsorption of Interacting k-Mers
To simulate
k-mers of increasing length, the standard grand canonical Monte Carlo (MC) algorithm (
Section 8.1.1) typically suffers from slow dynamics, particularly at medium to high coverage. This limitation necessitates the use of alternative algorithms. One such method is the hyper-parallel tempering Monte Carlo (HPTMC) [
175] algorithm. The HPTMC method consists in generating a compound system of
R non-interacting replicas of the system under study. Each replica is associated with a gas pressure
, taken from a set of properly selected pressures
[
176]. To determine the number of sampled pressures we used an acceptance ratio of
for the swapping move for each pair of replicas. Once the values of the gas pressure or the chemical potential are established, the simulation process consist in two major subroutines:
replica-update and
replica-exchange.
Replica-update.
The adsorption-desorption procedure is as follows: (i) One out of R replicas is randomly selected; (ii) a linear k-uple of nearest-neighbor sites is selected. Then, if the k sites are empty, an attempt is made to deposit a rod with probability ; if the k sites are occupied by units belonging to the same k-mer, an attempt is made to desorb this k-mer with probability and otherwise, the attempt is rejected. As in previous section, .
Replica-exchange.
Exchange of two configurations
and
, corresponding to the
i-th and
j-th replicas, respectively, is tried and accepted with probability,
, where
in a nonthermal grand canonical ensemble is given by,
8.3. Parallel Tempering MC Algorithm for Adsorption of Binary Mixtures of Interacting Species of Polyatomics
The above algorithms, in particular HPTMC, can be applied to the case of mixtures of two species. In order to simulate the adsorption of binary mixtures of
k-mers and
l-mers we consider a substrate consisting in a regular lattice with connectivity
c. The HPTMC method consists in generating a compound system of
R non-interacting replicas of the system under study. Each replica is associated with a gas pressure
, taken from a set of properly selected pressures
[
176]. To determine the number of sampled pressures we used an acceptance ratio of
for the swapping move for each pair of replicas. Once the values of the gas mixture pressure and molar fractions
are set, the chemical potential of each species is obtained an ideal gas mixture, i.e,
, where
is the standard chemical potential at temperature
T and
.
Under these considerations, the simulation process consist in two major subroutines: replica-update and replica-exchange.
Replica-update.
The adsorption-desorption procedure is as follows: (i) One out of R replicas is randomly selected; (ii) the species x is selected with equal probability from the two species, k and l; (iii) a linear x-uple of nearest-neighbor sites is selected. Then, if the x sites are empty, an attempt is made to deposit a rod with probability ; if the x sites are occupied by units belonging to the same x-mer, an attempt is made to desorb this x-mer with probability and otherwise, the attempt is rejected.
Replica-exchange.
Exchange of two configurations
and
, corresponding to the
i-th and
j-th replicas, respectively, is tried and accepted with probability,
, where
in a nonthermal grand canonical ensemble is given by,
The complete simulation procedure is the following: (1) replica-update, (2) replica-exchange, and (3) repeat from step (a) times. This is the elementary step in the simulation process or Monte Carlo step (MCs). Typically, the equilibrium state is reached after discarding the first MCs. Then, the next m MCs are used to compute averages.
For each value of pressure
, the corresponding surface coverages are determinate by simple averages,
where
represents the state of the replica
j-th at the Monte Carlo time
t.
8.4. Improving the Update Algorithm Through the Use of Lists of Full and Empty k-Uples
By using lists of full and empty k-uples it is possible to improve the performance of previous algorithms. The process of updating (replica-update) involves the random selection of a linear k-uple in the system. This k-uple may turn out to be occupied (by exactly one k-mer or rod), it may be completely empty, or it may have a partial occupation. In the last case the current update process must be rejected. It is possible to perform this random selection in a rejection-free manner using both full and empty k-uples lists. The cost of managing and using the lists is negligible compared to the advantage of avoiding rejection caused by the selection of partially occupied k-uples.
In the case of HPTMC, the procedure is as follows: (i) one out of R replicas is randomly selected; (ii) an element from the compound list of full and empty k-uplas, corresponding to the replica selected in (i), is randomly selected. This element is precisely an empty k-upla or a k-upla occupied by a k-mer. If the k-uple is empty, an attempt is made to deposit a rod with probability ; if the k-uple is occupied, an attempt is made to desorb this k-mer with probability . Finally, if any of the changes have been accepted, the lists of k-uples are updated accordingly.
8.5. Non-Local Update Kundu’s Algorithm for Adsorption of Non-Interacting Large k-Mers (Only Excluded Volume Interaction)
Simulations of
k-mers lattice gases were carried out in the Grand Canonical Ensemble through an efficient algorithm introduced by Kundu et al. [
101,
111,
112] to overcome the sampling slowdown at very high density due to the jamming effects. The temperature, chemical potential
and system’s size are held fixed and the number of particles on the lattice is allowed to fluctuate through non-local changes, i.e, insertion and deletion of many
k-mers at a time (in contrast to the standard Metropolis rule used in previous algorithms).
For a given configuration of
k-mers on the square lattice, one MCS is fulfilled by removing all the horizontal
k-mers and keeping the vertical ones. This process is shown in
Figure 67a,b. Thus, if the system is observed in the horizontal direction, it consists of segments of empty sites of different lengths, separated from each other by vertical
k-mers (see for example row
j in
Figure 67b). The probabilities of having horizontal segments of unoccupied sites can be exactly calculated and stored for all the possible segment sizes (this is essentially a 1D problem). All the horizontal segments are then reoccupied by
k-mers and empty sites with the probabilities accordingly determined (see
Figure 68). Finally, all these steps are repeated in the vertical direction.
Figure 67 shows the process of removal of horizontal (red)
k-mers, and the identification of three segments in a row.
Figure 68 shows the process of occupation of a
segment. The algorithm can be easily generalized to other geometries and dimensions. The detailed discussion of the algorithm can be found in the original work Refs. [
101,
111,
112]. The algorithm has proved to be ergodic, it satisfies the Detailed Balance Principle and equilibrium is reached after typically
MCSs.
8.6. Thermodynamic Integration Method in Canonical Ensemble: Artificial Hamiltonian Method
The advantages of using Monte Carlo simulation to calculate thermal averages of thermodynamic observables are well known [
177]. The estimation of certain quantities such as total energy, energy fluctuations, correlation functions, etc., is rather straightforward from averaging over a large enough number of instantaneous configurations (states) of a thermodynamic system. However, free energy and entropy, in general, cannot be directly computed. In order to calculate free energy and entropy, various methods have been developed. Namely, thermodynamic integration method (TIM) [
74,
75,
177,
178,
179,
180], Ma’s method of coincidence counting of states along the trajectory in phase space [
181], “stochastic models" method of Alexandrowicz [
182], “local states" method of Meirovitch [
183], “multistage sampling" and “umbrella sampling" of Valleau et al. [
184,
185,
186,
187], method of Salsburg [
188], method of Yip et al. [
189] (which is an optimized combination of coupling parameter and adiabatic switching formalisms), etc.
Among the methods mentioned in the last paragraph, the TIM is one of the most widely used and practically applicable. In the following, we briefly describe this method.
Given a lattice-gas of
N interacting particles on a regular lattice with
M sites at temperature
T, from the basic relationship
it follows
where
U is the mean total energy of the system.
is readily calculated if (reference state) is known, given that the integral in the second term can be accurately estimated by MC simulation. In practice, the calculation of S in a reference state can be rigorously accomplished by analytical methods only in a very few cases. Although the entropy of some particular states is trivially known (for example ), this is often computationally inconvenient since it would require the simulation of a thermodynamic open system to get the entropy of a state at finite density. Alternatively, integration can be carried out through a thermodynamic path of a closed (mechanically isolated) system along a constant density path, if a proper reference state is defined for which can be directly computed.
In the case monomers (
), the determination of the entropy in the reference state is trivial. In fact, for a monoatomic lattice-gas
The last equation holds for any finite value of the lateral interactions between the ad-particles.
Since cannot be exactly calculated for k-mer adsorption () by analytical means, in the following we present a general numerical methodology to obtain the entropy of generalized lattice-gas in a reference state.
If an artificial lattice-gas is defined from the system of interest (henceforth referred to as the original system) such that it fulfills the condition,
Then, the integral in Equation (
382) can be separated into two terms. Thus,
where
and
U are the mean total energy of the artificial and original system, respectively (both integrals can be evaluated by MC in the canonical ensemble). The general definition of the artificial reference system follows.
Let us assume the original system to be a discrete system of N particles on M sites with Hamiltonian ; where finite , is the potential energy in the configuration among the set of accessible configurations . The original system can only have access to those configurations within ; the total amount of configurations in is (in a lattice gas of N monomers with single-site occupancy of M sites, ).
The Hamiltonian for the artificial system, , follows from:
Definition 1: , is defined as finite , where and have analogous meanings to those given above for U and , respectively. The equalities ensure that the set of accessible configurations for the original system and the artificial system are equal (although , the energy of the configurations in the artificial system may be, in general, different from the ones in the original system).
Definition 2: The potential energy of the accessible configurations (
) for the artificial system take the following values
Definition 2 means that a given configuration (the th) is selected arbitrarily from and defined as the non-degenerate ground state of the artificial lattice-gas; hence . In practice, the configuration can be easily defined.
An example for adsorbed dimers follows in order to make this point clear. Let us consider adsorbed dimers on an homogeneous square lattice with
and interaction between NN dimer’s heads as shown in
Figure 69 (original system). For this system there is no rigorous expression of
for
in the thermodynamic limit (
,
constant).
To build up an artificial system fulfilling the definitions 1 and 2, we follow the steps:
The number of particles, size and geometry of the lattice is kept as in the original system.
The interaction energy between NN units is set to zero.
An adsorption energy is introduced for the lattice sites (representing, for each site, the interaction between the lattice and the unit of the dimer adsorbed on it in the artificial system), in such a way that two types of sites are defined, strong and weak, with energies
and
, respectively, being
. For
N adsorbed dimers we choose
strong sites conveniently on the lattice. For instance, in
Figure 70a a possible distribution of strong and weak sites is depicted, where circles and squares are sites of energy
and
, respectively.
It is assumed that dimers in a particular direction are energetically favored. This is formally handled by introducing a virtual external field such that the interaction energy between the dimers and the field is
if the
dimer is vertically aligned and
otherwise (this choice is obviously arbitrary). Care must be taken if periodic boundary conditions are applied to ensure that the state of minimum energy is unique. Then, the Hamiltonian of the artificial system for this example is given by
where
if the site is strong and
if the site is weak.
Thus, the ground state of the artificial system is the one shown in
Figure 70(b), which is non-degenerate, giving
.
The calculation of
through Equation (
386) is straightforward and computationally simple, since the temperature dependence of
and
is evaluated at constant coverage for various values of
T following the standard procedure of MC simulation in the canonical ensemble (see
Section 8.1.2) (based on the Metropolis scheme [
174]). Then,
and
are obtained as simple averages, spline-fitted and numerically integrated. It should be mentioned that
is calculated by using the hamiltonian of Equation (
388) and
is calculated by using the original hamiltonian. Two typical curves of
versus
u are shown in
Figure 71, for attractive and repulsive dimers on a square lattice.
The strategy described above is applicable to a wide class of lattice-gas systems.
9. Conclusions and Future Perspectives
This review has systematically explored the thermodynamic behavior and statistical mechanics of adsorbed rigid rods (or k-mers) on regular lattices, focusing on equilibrium properties as a function of density, temperature, and particle size. Through a sequence of progressively refined approaches—ranging from exact results in one dimension to numerical simulations and mean-field approximations in higher dimensions, and culminating in the analytical development of Multiple Exclusion (ME) statistics—the article provides a broad and coherent framework for understanding collective adsorption phenomena governed by excluded volume effects.
In
Section 2, the review focused on one-dimensional (1D) systems, which offer the unique opportunity of exact solvability. The classic lattice versions of hard rod systems are revisited, emphasizing their exact entropy, equation of state, and chemical potential. These models serve not only as rigorous benchmarks but also as reference systems to test approximate theories. Thus single-species and mixtures on the lattice can be rigorously treated and phenomena such as reversal-adsorption in binary mixtures can be well characterized and understood in the light of the model system. The treatment of
k-mers in the multilayer regime bears an elucidating relation with analogous experimental realizations such as adsorption of linear molecules (alkanes, alkenes) on quasi-regular surfaces of carbon nanotubes leading to demonstrate that molecular the size/shape entropic contribution to free energy cannot be oversimplified in interpreting adsorption isotherms to determine the specific surface of the adsorbent.
The inclusion of interacting k-mers through lateral attractions or repulsions further enriches the 1D landscape, allowing for a better understanding of the features of the coverage dependence of entropy. Notably, the exact solutions for lattice k-mers lay the groundwork for understanding the role of orientational and configurational constraints in low-dimensional systems.
The 1D treatment sets the basis for extending the understanding to higher dimensions by introducing coarse-grained approaches that effectively capture key thermodynamic quantities. Here, the notion of excluded volume per particle is generalized, leading to approximate expressions for entropy and pressure that remain remarkably accurate across a broad range of densities. A number of approximations were reviewed, ranging from the configurational dimensional ansatz to extend the 1D analytical forms to higher dimensions through the EA approximation, the conceptualization of the challenging problem of the structured particles lattice gas as one isomorphic to fractional exclusion statistics in quantum systems with independent states through FSTA, up to the introduction of the most elaborate theoretical approach, named Multiple Exclusion Statistics, embodying the complexity emerging from spatial correlations between particle states into a set of simple self-exclusion and cross-exclusion correlation parameters that can be consistently determined within the theory from limiting configurational conditions of the system. .
The article then progresses to address more intricate systems in higher dimensions, where exact solutions are generally not accessible. Here, mean-field approximations and phenomenological entropy-based models are employed to describe the emergence of phase transitions, particularly the isotropic-to-nematic transition observed for elongated particles. These transitions are governed not only by steric constraints but also by emergent orientational ordering, which is absent in 1D. The entropy per site as a function of coverage reveals non-trivial behavior, including entropy-driven phase separations and coexistence regions, which can be tracked even in approximate analytical formulations.
Section 4 presents a series of heuristic and mean-field-like methods aimed at incorporating orientational degrees of freedom and spatial correlations in 2D and 3D lattices. The simple Bragg-Williams-type approximations provide qualitative insights into isotropic-nematic transitions and the role of
k-mer length. Improved approaches, such as those based on lattice free volume estimates or effective excluded area models like the Quasi-Chemical approach, offer better agreement with simulation data, particularly. However, these approximations often neglect inter-particle correlations and fail to describe critical behavior accurately.
The review takes a major conceptual leap in
Section 5 and
Section 6, where the statistical description of the system is reframed through the lens of generalized exclusion statistics, particularly in the form of Multiple Exclusion (ME) statistics. This approach introduces the idea of state-counting based on effective exclusion rules that go beyond simple geometric packing. By systematically encoding how a particle’s presence affects the availability of nearby states, the ME formalism captures both entropic interactions and emergent correlations. This represents a profound shift from traditional lattice gas models toward a more abstract and generalizable statistical mechanics framework.
One of the most significant achievements of this formalism is the derivation of spectral exclusion functions, which encode the exclusion parameters between particle states. These functions allow for the construction of thermodynamic quantities—such as the chemical potential and entropy from the statistical exclusion rules, circumventing the need for explicit partition functions. Furthermore, the ability to reproduce not only exact 1D results but also to approximate higher-dimensional behavior with improved accuracy underscores the power of the ME approach.
A key theoretical development discussed is the formalization of ME statistics to include multi-component and geometrically complex particle mixtures. The introduction of state self-exclusion and state cross-exclusion functions and , as generalizations of chemical potential derivatives with respect to partial densities, allows for a thermodynamic characterization of complex exclusion scenarios. These spectral functions encode how the presence of particles of one species affects the availability of configurations for others, thus extending the exclusion principle into a spatially resolved thermodynamic descriptor.
It is also showcased the application of these ideas to realistic systems (
Section 7), such as adsorption on heterogeneous surfaces, adsorption of alkanes on carbon nanotubes as a prototype experimental realization to show the simplest multilayer adsorption model generalizing the pioneer work of BET to determine adsorption energy and specific area of adsorbent when probe molecules are not ideally spherical. Furthermore, application to crystal growth from aqueous solutions in presence of structured impurities were exhibited as applications were the impurity size and shape must be properly accounted for in the thermodynamic potentials in order to understand the growth rate dependence on impurity concentration.
The application examples concluded with an in-depth analysis of the long-standing k-mer problem on the square lattice, approached through the advanced framework of Multiple Exclusion Statistics theory. This problem, which still presents several open questions, such as the lack of a unified theoretical explanation for the emergence of a nematic transition at the critical length , and the nature of the high-density transition back to an isotropic phase near saturation, was revisited from the alternative and more comprehensive perspective provided by the ME formalism. By modeling the system as a mixture of two species of particles with distinct orientations, the free energy and entropy surfaces were analytically approximated across the full range of relative occupations along lattice directions. This approach enabled a detailed characterization of density branches, critical points, transition orders, and exclusion spectrum functions, offering new insights into the entropy–density relationship. The analysis reveals that is the minimal rod length required to undergo a nematic transition, while , although close, does not satisfy the necessary conditions and remains in an isotropic phase even at high coverage.
Furthermore, the review explores the novel outcomes of ME statistics such as the potential of extracting exclusion spectral functions from experimental data, inverting thermodynamic observables to infer structural information. This opens the possibility of using adsorption isotherms and fluctuations spectra not just as phenomenological descriptors, but as quantitative probes of configurational entropy and microstate topology. Remarkably, the ME statistics theoretical framework provides a bridge between observable thermodynamic quantities (e.g., chemical potential versus density curves) and latent spatial correlations or ordering tendencies that are otherwise hidden. This fact opens new experimental possibilities: exclusion spectra may eventually be inferred from measurements of adsorption isotherms or compressibility, providing insight into microscopic ordering without relying on direct imaging.
Ultimately, the ME formalism proves to be a unifying language capable of handling complex interactions and geometry-induced constraints. It not only serves as a powerful lens for interpreting equilibrium properties of complex particle systems, but also sets the stage for future theoretical and experimental explorations into the geometry and thermodynamics of constrained configuration spaces.
Several open questions and research avenues arise from this unified perspective. First, the extension of ME statistics to non-lattice systems, including continuous 2D and 3D domains with quenched disorder or curved geometries, remains largely unexplored.
Its generality and compatibility with empirical data make it not only a theoretical construct but also a bridge toward experimentally accessible quantities. From a computational standpoint, there is room for developing inverse statistical mechanical methods that reconstruct exclusion spectral functions directly from experimental data or Monte Carlo simulations. This would establish a concrete protocol for extracting statistical fingerprints from real systems, with applications in surface science, porous media, and biological adsorption.
Secondly, dynamical aspects such as adsorption/desorption kinetics, diffusion on fluctuating energy landscapes, or driven systems under external fields, are fertile grounds for applying and testing the ME framework beyond equilibrium.
Looking ahead, the ME framework is poised to play a central role in the statistical mechanics of systems with constrained configurations, such as crowding in biological environments, adsorption in porous media, or active matter with limited motility space.
Finally, there is potential to link the ME statistics approach with field-theoretic and renormalization group methods, particularly in systems where critical behavior or universality classes may be modified by spatial exclusion. The interplay between exclusion-driven entropy and geometric frustration also presents an intriguing challenge, especially in contexts where topology and boundary effects play a dominant role.
A whole
Section 8 was devoted to review the Monte Carlo techniques from the early Metropolis state sampling to the new highly efficient non-local configuration update to overcome the sampling slowdown to reach equilibrium in lattice gases with state exclusion, particularly at high density. Furthermore, the artificial-Hamiltonian technique was revisited as powerful tool to calculate entropy at arbitrary density from elementary thermodynamic relations taking advantage of efficient numerical sampling of equilibrium configurations.
Monte Carlo (MC) simulations serve as an essential tool to validate approximate analytical methods for single components and mixtures as well as exploring the influence of the lattice geometry on the thermodynamic potentials. MC studies have uncovered rich phase behavior, including continuous and discontinuous isotropic-nematic transitions, as well as layering and jamming effects near full coverage. These results have helped to clarify the limitations of earlier mean-field approximations and motivated the development of more refined theories. Importantly, the numerical data have also inspired empirical functional forms for entropy and chemical potential, which have been instrumental in bridging the gap between simulations and analytical descriptions. The multiple exclusion problem opens new questions concerning potential simulation techniques in systems with strong state correlations and geometrical constrains exploiting the topological properties of the configuration space and graph theory.
In summary, this article accomplishes several goals: It consolidates a wide range of exact, approximate, and numerical results across spatial dimensions and k-mer lengths, it clarifies the domain of validity and limitations of various theoretical approaches, it reviews the recently presented Multiple Exclusion Statistics (ME) as a promising analytical tool that reconciles geometric exclusion under topological constrains and correlations with thermodynamic consistency and it sets new perspectives for the application and testing of this newest analytical tool to a broader scenario of systems where spatial correlations of particle states dominate.
The synthesis of exact models, approximate theories, numerical validation, and generalized statistics presented in this review thus offers a comprehensive and forward-looking perspective on k-mer adsorption and excluded volume systems, with broader implications across condensed matter, materials science, and statistical physics.
Figure 1.
Surface coverage
versus relative chemical potential
. Comparison between adsorption isotherms for
k-mers in one dimension from Flory-Huggins’s approximation (dashed line) and exact solution (solid line) from Equation (
15) for dimers, 4-mers, 10-mers and 100-mers (curves from top to bottom). MC simulation for
and
in one dimension are shown in full circles.
Figure 1.
Surface coverage
versus relative chemical potential
. Comparison between adsorption isotherms for
k-mers in one dimension from Flory-Huggins’s approximation (dashed line) and exact solution (solid line) from Equation (
15) for dimers, 4-mers, 10-mers and 100-mers (curves from top to bottom). MC simulation for
and
in one dimension are shown in full circles.
Figure 2.
Molar entropy (in units of ) versus surface coverage. Comparison between exact results (solid line) and Flory-Huggins’s approach (dashed line) for monomers (), dimers () and polymers ( and ). Curves from top to bottom (the case is common for both results).
Figure 2.
Molar entropy (in units of ) versus surface coverage. Comparison between exact results (solid line) and Flory-Huggins’s approach (dashed line) for monomers (), dimers () and polymers ( and ). Curves from top to bottom (the case is common for both results).
Figure 3.
Spreading pressure in units of versus surface coverage. Comparison between exact results (solid line) and Flory-Huggins’s approach (dashed line) for monomers (), dimers () and polymers ( and ). Curves from top to bottom (the case is common for both results).
Figure 3.
Spreading pressure in units of versus surface coverage. Comparison between exact results (solid line) and Flory-Huggins’s approach (dashed line) for monomers (), dimers () and polymers ( and ). Curves from top to bottom (the case is common for both results).
Figure 4.
Total and partial adsorption isotherms for a monomer-dimer binary mixture adsorbed on a one-dimensional lattice with and . Symbols and solid lines represent MC and theoretical data, respectively.
Figure 4.
Total and partial adsorption isotherms for a monomer-dimer binary mixture adsorbed on a one-dimensional lattice with and . Symbols and solid lines represent MC and theoretical data, respectively.
Figure 5.
Cartoon representing the lattice gas model of k-mer adsorption in the multilayer regime (the case of trimer adsorption is depicted). Adsorption sites are represented by crosses on the adsorbent’s surface. The adsorbate molecule is represented by black beds connected by solid bonds.
Figure 5.
Cartoon representing the lattice gas model of k-mer adsorption in the multilayer regime (the case of trimer adsorption is depicted). Adsorption sites are represented by crosses on the adsorbent’s surface. The adsorbate molecule is represented by black beds connected by solid bonds.
Figure 6.
Adsorption isotherms for different values of the parameter
c. The solid line represents the isotherm for dimers obtained in the present work [Equation (
60)], and the dashed line does for the BET’s isotherm (monomers). The pairs of curves from bottom to top correspond to
and 100, in the range
of
in
Figure 6a) and
of
in
Figure 6b).
Figure 6.
Adsorption isotherms for different values of the parameter
c. The solid line represents the isotherm for dimers obtained in the present work [Equation (
60)], and the dashed line does for the BET’s isotherm (monomers). The pairs of curves from bottom to top correspond to
and 100, in the range
of
in
Figure 6a) and
of
in
Figure 6b).
Figure 7.
versus
for
(in arbitrary units), and
(curves from bottom to top). Solid and dashed curves represent the same models as in
Figure 6.
Figure 7.
versus
for
(in arbitrary units), and
(curves from bottom to top). Solid and dashed curves represent the same models as in
Figure 6.
Figure 8.
Adsorption isotherms for interacting (repulsive as well attractive) dimers and 10-mers; (a)
,
; (b)
,
; (c)
,
; (d)
,
. The full circles represent results from Monte Carlo Simulation in the lattice-gas model. Comparison between experimental isotherm of CH
4 from Ref. [
64] (open triangles,
K; open diamonds,
K) and theoretical isotherm [from Equation (
73)] of dimers in the lattice-gas approximation (full triangles,
K,
Kcal/mol; full diamonds
K,
Kcal/mol).
Figure 8.
Adsorption isotherms for interacting (repulsive as well attractive) dimers and 10-mers; (a)
,
; (b)
,
; (c)
,
; (d)
,
. The full circles represent results from Monte Carlo Simulation in the lattice-gas model. Comparison between experimental isotherm of CH
4 from Ref. [
64] (open triangles,
K; open diamonds,
K) and theoretical isotherm [from Equation (
73)] of dimers in the lattice-gas approximation (full triangles,
K,
Kcal/mol; full diamonds
K,
Kcal/mol).
Figure 9.
Entropy per site, , versus lattice coverage, , for interacting (repulsive as well attractive) dimers and 10-mers; (a) , ; (b) , , (c) , ; (d) , ; (e) , ; (f) , . Symbols and lines represent MC simulation data and theoretical results, respectively.
Figure 9.
Entropy per site, , versus lattice coverage, , for interacting (repulsive as well attractive) dimers and 10-mers; (a) , ; (b) , , (c) , ; (d) , ; (e) , ; (f) , . Symbols and lines represent MC simulation data and theoretical results, respectively.
Figure 10.
Specific heat in units of (at constants volume) versus at . (a) monomers (); (b) dimers (); (c) trimers (); (d) 10-mers (). The cases (a) and (b) correspond to the left hand side axis and the cases (c) and (d) correspond to the right hand side axis.
Figure 10.
Specific heat in units of (at constants volume) versus at . (a) monomers (); (b) dimers (); (c) trimers (); (d) 10-mers (). The cases (a) and (b) correspond to the left hand side axis and the cases (c) and (d) correspond to the right hand side axis.
Figure 11.
Schematic representation of a lattice-gas of dimers (, blue circles) and trimers (, red circles) adsorbed on a one-dimensional lattice. The figure shows different types of pairs of sites: a) 00, b) , c) , d) , e) and f) .
Figure 11.
Schematic representation of a lattice-gas of dimers (, blue circles) and trimers (, red circles) adsorbed on a one-dimensional lattice. The figure shows different types of pairs of sites: a) 00, b) , c) , d) , e) and f) .
Figure 12.
Rules for the mapping , from the original lattice of polyatomics to an effective lattice of monomers .
Figure 12.
Rules for the mapping , from the original lattice of polyatomics to an effective lattice of monomers .
Figure 13.
Partial adsorption isotherms for the 1D case, , and different values of the lateral interaction as indicated. Symbols correspond to MC data whereas lines represent the theoretical results.
Figure 13.
Partial adsorption isotherms for the 1D case, , and different values of the lateral interaction as indicated. Symbols correspond to MC data whereas lines represent the theoretical results.
Figure 14.
Linear tetramers adsorbed on (a) square, (b) triangular and (c) honeycomb lattices. Full and empty circles represent tetramer units and empty sites, respectively.
Figure 14.
Linear tetramers adsorbed on (a) square, (b) triangular and (c) honeycomb lattices. Full and empty circles represent tetramer units and empty sites, respectively.
Figure 15.
Schematic representation of a lattice-gas of dimers (, blue circles) and trimers (, red circles) adsorbed on a square lattice (). The figure shows different types of pairs of sites: a) , b) , c) , d) 00, e) and f) .
Figure 15.
Schematic representation of a lattice-gas of dimers (, blue circles) and trimers (, red circles) adsorbed on a square lattice (). The figure shows different types of pairs of sites: a) , b) , c) , d) 00, e) and f) .
Figure 16.
Entropy per site vs. lattice coverage
for
. Dashed lines: isotropic phase (I). Solid lines: fully aligned nematic phase (N). First-order
approximation through Equations (
285), (
294), (
300).
Figure 16.
Entropy per site vs. lattice coverage
for
. Dashed lines: isotropic phase (I). Solid lines: fully aligned nematic phase (N). First-order
approximation through Equations (
285), (
294), (
300).
Figure 17.
Critical coverages (lower branch) and (upper branch) as functions of k. Open squares: first-order approximation. Solid squares: second-order approximation.
Figure 17.
Critical coverages (lower branch) and (upper branch) as functions of k. Open squares: first-order approximation. Solid squares: second-order approximation.
Figure 18.
Entropy per site vs. lattice coverage
for
and 8. Symbols: MC data [
109]. Dashed:
model (isotropic phase) with
(
),
(
). Solid:
result for nematic phase (
). Dotted: empirical correction
. Inset: high coverage behavior for
.
Figure 18.
Entropy per site vs. lattice coverage
for
and 8. Symbols: MC data [
109]. Dashed:
model (isotropic phase) with
(
),
(
). Solid:
result for nematic phase (
). Dotted: empirical correction
. Inset: high coverage behavior for
.
Figure 19.
Lattice coverage versus for k-mers on a square lattice. : , , , (dashed line), (solid line); : , , , (dotted line). Symbols denote MC simulation results: squares for , circles for .
Figure 19.
Lattice coverage versus for k-mers on a square lattice. : , , , (dashed line), (solid line); : , , , (dotted line). Symbols denote MC simulation results: squares for , circles for .
Figure 20.
Exclusion spectrum function
vs. lattice coverage
. Lines:
theory from Equation (
291); dot-dashed (
), dashed (
), solid (
rectangles), dotted (
squares). Symbols: MC data as indicated.
Figure 20.
Exclusion spectrum function
vs. lattice coverage
. Lines:
theory from Equation (
291); dot-dashed (
), dashed (
), solid (
rectangles), dotted (
squares). Symbols: MC data as indicated.
Figure 21.
Same as
Figure 20, but for the exclusion per particle frequency function
.
Figure 21.
Same as
Figure 20, but for the exclusion per particle frequency function
.
Figure 22.
A schematic representation of a ternary mixture: a tetramer (blue), a dimer (red), and a monomer (green). The dotted lines indicate examples of multiply excluded states within the spectrum of (a) tetramers, (b) dimers due to tetramers and monomers, and (c) dimers via self-exclusion. This emphasizes the multiplicity of state exclusions even at finite densities.
Figure 22.
A schematic representation of a ternary mixture: a tetramer (blue), a dimer (red), and a monomer (green). The dotted lines indicate examples of multiply excluded states within the spectrum of (a) tetramers, (b) dimers due to tetramers and monomers, and (c) dimers via self-exclusion. This emphasizes the multiplicity of state exclusions even at finite densities.
Figure 24.
Comparison between the exact adsorption isotherm of monomers and the simulation adsorption isotherms of dimers on honeycomb, square and triangular lattices.
Figure 24.
Comparison between the exact adsorption isotherm of monomers and the simulation adsorption isotherms of dimers on honeycomb, square and triangular lattices.
Figure 25.
Adsorption isotherms of 6-mers on a honeycomb lattice. Symbols represent MC results, and lines correspond to different approaches (see inset). Percentage reduced coverage, , versus surface coverage. The symbols are as in part .
Figure 25.
Adsorption isotherms of 6-mers on a honeycomb lattice. Symbols represent MC results, and lines correspond to different approaches (see inset). Percentage reduced coverage, , versus surface coverage. The symbols are as in part .
Figure 26.
As
Figure 5 for a square lattice.
Figure 26.
As
Figure 5 for a square lattice.
Figure 27.
As
Figure 5 for a triangular lattice.
Figure 27.
As
Figure 5 for a triangular lattice.
Figure 28.
Percentage reduced coverage versus concentration for k-mers adsorbed on a honeycomb lattice and approximation. Symbols are indicated in the inset.
Figure 28.
Percentage reduced coverage versus concentration for k-mers adsorbed on a honeycomb lattice and approximation. Symbols are indicated in the inset.
Figure 29.
As
Figure 28 for a square lattice.
Figure 29.
As
Figure 28 for a square lattice.
Figure 30.
As
Figure 28 for a triangular lattice.
Figure 30.
As
Figure 28 for a triangular lattice.
Figure 31.
a) Average percentage reduced coverage , as a function of k for different connectivities. b) As in part a) for the maximum percentage reduced coverage . Hexagons, squares and triangles correspond to data obtained for honeycomb, square an triangular lattices, respectively.
Figure 31.
a) Average percentage reduced coverage , as a function of k for different connectivities. b) As in part a) for the maximum percentage reduced coverage . Hexagons, squares and triangles correspond to data obtained for honeycomb, square an triangular lattices, respectively.
Figure 32.
Adsorption isotherms for rigid molecules on a square lattice: (a) and ; (b) and . (c)[(d)] Reduced coverage error vs. the total surface coverage for the data in part (a)[(b)]. Symbols represent MC results, and lines correspond to different theoretical approaches as indicated.
Figure 32.
Adsorption isotherms for rigid molecules on a square lattice: (a) and ; (b) and . (c)[(d)] Reduced coverage error vs. the total surface coverage for the data in part (a)[(b)]. Symbols represent MC results, and lines correspond to different theoretical approaches as indicated.
Figure 33.
Same as
Figure 32 for triangular lattices.
Figure 33.
Same as
Figure 32 for triangular lattices.
Figure 34.
Average error for fixed k and different values of l: (a) square lattice, and (b) triangular lattice, . The meaning of the lines is indicated in the figure.
Figure 34.
Average error for fixed k and different values of l: (a) square lattice, and (b) triangular lattice, . The meaning of the lines is indicated in the figure.
Figure 35.
Comparison between experimental and theoretical adsorption isotherms for a binary mixture: (a) CH4/C2H6 and (b) C2H6/C3H6 adsorbed on a commercial activated carbon. Symbols represent experimental data from Ref. [
116] and lines correspond to results from Equation (
193). The parameters used in the fitting procedure are listed in
Table 1.
Figure 35.
Comparison between experimental and theoretical adsorption isotherms for a binary mixture: (a) CH4/C2H6 and (b) C2H6/C3H6 adsorbed on a commercial activated carbon. Symbols represent experimental data from Ref. [
116] and lines correspond to results from Equation (
193). The parameters used in the fitting procedure are listed in
Table 1.
Figure 36.
Adsorption isotherms for homonuclear dimers adsorbed on a honeycomb lattice with nearest-neighbor interactions. Symbols, solid lines and dashed lines represent results from Monte Carlo simulations, and , respectively. The dotted lines are included in the figure as a guide for the eyes. a) Attractive case: full circles, ; open circles, ; open squares, ; open up triangles, ; open diamonds, and open down triangles, . b) Repulsive case: full circles, ; full hexagons, ; full squares, ; full up triangles, ; full diamonds, ; full down triangles, and full stars, . Inset: Adsorption isotherms from for and different values of k as indicated.
Figure 36.
Adsorption isotherms for homonuclear dimers adsorbed on a honeycomb lattice with nearest-neighbor interactions. Symbols, solid lines and dashed lines represent results from Monte Carlo simulations, and , respectively. The dotted lines are included in the figure as a guide for the eyes. a) Attractive case: full circles, ; open circles, ; open squares, ; open up triangles, ; open diamonds, and open down triangles, . b) Repulsive case: full circles, ; full hexagons, ; full squares, ; full up triangles, ; full diamonds, ; full down triangles, and full stars, . Inset: Adsorption isotherms from for and different values of k as indicated.
Figure 37.
Adsorption isotherms for homonuclear dimers adsorbed on a square lattice with nearest-neighbor interactions. Symbols, solid lines and dashed lines represent results from Monte Carlo simulations, and , respectively. The dotted lines are included in the figure as a guide for the eyes. a) Attractive case: full circles, ; open squares, ; open up triangles, ; open diamonds, and open down triangles, . b) Repulsive case: full circles, ; full squares, ; full up triangles, ; full diamonds, and full down triangles, . Inset: Adsorption isotherms from for and different values of k as indicated.
Figure 37.
Adsorption isotherms for homonuclear dimers adsorbed on a square lattice with nearest-neighbor interactions. Symbols, solid lines and dashed lines represent results from Monte Carlo simulations, and , respectively. The dotted lines are included in the figure as a guide for the eyes. a) Attractive case: full circles, ; open squares, ; open up triangles, ; open diamonds, and open down triangles, . b) Repulsive case: full circles, ; full squares, ; full up triangles, ; full diamonds, and full down triangles, . Inset: Adsorption isotherms from for and different values of k as indicated.
Figure 38.
Adsorption isotherms for homonuclear dimers adsorbed on a triangular lattice with nearest-neighbor interactions. Symbols, solid lines and dashed lines represent results from Monte Carlo simulations, and , respectively. The dotted lines are included in the figure as a guide for the eyes. a) Attractive case: full circles, ; open circles, ; open squares, ; open diamonds, and open up triangles, . b) Repulsive case: full circles, ; full squares, ; full up triangles, ; full diamonds, and full down triangles, . Inset: Adsorption isotherms from for and different values of k as indicated.
Figure 38.
Adsorption isotherms for homonuclear dimers adsorbed on a triangular lattice with nearest-neighbor interactions. Symbols, solid lines and dashed lines represent results from Monte Carlo simulations, and , respectively. The dotted lines are included in the figure as a guide for the eyes. a) Attractive case: full circles, ; open circles, ; open squares, ; open diamonds, and open up triangles, . b) Repulsive case: full circles, ; full squares, ; full up triangles, ; full diamonds, and full down triangles, . Inset: Adsorption isotherms from for and different values of k as indicated.
Figure 39.
Snapshots of the ordered phases corresponding to repulsive dimers adsorbed on a honeycomb lattice. a) Low-coverage ordered structure (LCOP); b) high-coverage ordered structure (HCOP) and c) LCOP-HCOP mixture according to the predictions of QCA.
Figure 39.
Snapshots of the ordered phases corresponding to repulsive dimers adsorbed on a honeycomb lattice. a) Low-coverage ordered structure (LCOP); b) high-coverage ordered structure (HCOP) and c) LCOP-HCOP mixture according to the predictions of QCA.
Figure 41.
As
Figure 39 for triangular lattices.
Figure 41.
As
Figure 39 for triangular lattices.
Figure 42.
Absolute error , versus surface coverage for adsorption isotherms of dimers. The symbology is as follows: a) Squares, ; triangles, and circles, . b) Squares, ; triangles, and circles, . Full and open symbols correspond to comparisons with and , respectively.
Figure 42.
Absolute error , versus surface coverage for adsorption isotherms of dimers. The symbology is as follows: a) Squares, ; triangles, and circles, . b) Squares, ; triangles, and circles, . Full and open symbols correspond to comparisons with and , respectively.
Figure 43.
Integral error , versus lateral interaction (in units) for different geometries as indicated.
Figure 43.
Integral error , versus lateral interaction (in units) for different geometries as indicated.
Figure 44.
Adsorption isotherms for
adsorbed in
zeolite fitted by
. Symbols correspond to data from Ref. [
120] and lines represent theoretical results from Equation (
336).
Figure 44.
Adsorption isotherms for
adsorbed in
zeolite fitted by
. Symbols correspond to data from Ref. [
120] and lines represent theoretical results from Equation (
336).
Figure 45.
Derivative of the adsorbed amount versus pressure (
) for the same set of data plotted in
Figure 44.
Figure 45.
Derivative of the adsorbed amount versus pressure (
) for the same set of data plotted in
Figure 44.
Figure 47.
Temperature dependence of equilibrium constant,
, arising from fitting. Squares, red circles and black circle correspond to fitting from Refs. [
120,
121,
122], respectively.
reported in
Table 2 is the absolute value of the slope of the solid line.
Figure 47.
Temperature dependence of equilibrium constant,
, arising from fitting. Squares, red circles and black circle correspond to fitting from Refs. [
120,
121,
122], respectively.
reported in
Table 2 is the absolute value of the slope of the solid line.
Figure 48.
Comparison between Monte Carlo simulations of dimers, trimers and tetramers adsorbed on square lattices and theoretical isotherm from
[Equation (
125)].
g values from best fitting are indicated in the figure for each case. Inset: Configurational entropy versus coverage for dimers adsorbed on square lattices. Symbols represent simulation data and lines provide theoretical results for different values of
g as indicated in the text.
Figure 48.
Comparison between Monte Carlo simulations of dimers, trimers and tetramers adsorbed on square lattices and theoretical isotherm from
[Equation (
125)].
g values from best fitting are indicated in the figure for each case. Inset: Configurational entropy versus coverage for dimers adsorbed on square lattices. Symbols represent simulation data and lines provide theoretical results for different values of
g as indicated in the text.
Figure 49.
Partial adsorption isotherms for an equimolar monomer()-monomer() mixture on an one-dimensional lattice. Parameter values (figure): and . Parameter values (inset): , , and .
Figure 49.
Partial adsorption isotherms for an equimolar monomer()-monomer() mixture on an one-dimensional lattice. Parameter values (figure): and . Parameter values (inset): , , and .
Figure 50.
Partial adsorption isotherms for an equimolar methane()-ethane() mixture on an one-dimensional lattice. Parameter values: , , and .
Figure 50.
Partial adsorption isotherms for an equimolar methane()-ethane() mixture on an one-dimensional lattice. Parameter values: , , and .
Figure 51.
The adsorption isotherm for methane on single-walled carbon nanotubes is shown, with the coverage (expressed in cm3 Torr/g) plotted on the Y-axis against the relative pressure on the X-axis. Initially, there is a steep rise in coverage at low pressures, followed by a region where the coverage increases linearly with pressure. The point at which the isotherm first diverges from this linear behavior (indicated by the arrow aligned with the X-axis in the figure) is identified as point B. This point is interpreted as corresponding to the completion of the monolayer.
Figure 51.
The adsorption isotherm for methane on single-walled carbon nanotubes is shown, with the coverage (expressed in cm3 Torr/g) plotted on the Y-axis against the relative pressure on the X-axis. Initially, there is a steep rise in coverage at low pressures, followed by a region where the coverage increases linearly with pressure. The point at which the isotherm first diverges from this linear behavior (indicated by the arrow aligned with the X-axis in the figure) is identified as point B. This point is interpreted as corresponding to the completion of the monolayer.
Figure 52.
Adsorption isotherms for methane (77.3K), ethane (165 K), propane (190 K), and butane (220 K) adsorption on single-walled carbon nanotubes. The coverage in cm3 Torr/g (Y axis) is presented as a function of relative pressure (X axis).
Figure 52.
Adsorption isotherms for methane (77.3K), ethane (165 K), propane (190 K), and butane (220 K) adsorption on single-walled carbon nanotubes. The coverage in cm3 Torr/g (Y axis) is presented as a function of relative pressure (X axis).
Figure 53.
BET analysis for the adsorption isotherms for methane, ethane, propane, and butane on single-walled carbon nanotubes shown in
Figure 52.
Figure 53.
BET analysis for the adsorption isotherms for methane, ethane, propane, and butane on single-walled carbon nanotubes shown in
Figure 52.
Figure 54.
Specific surface area of single-walled carbon nanotubes computed using the BET and the point B methods. The specific surface area in m2/g (Y axis) is presented as a function of number of carbon atoms in the adsorbate (X axis).
Figure 54.
Specific surface area of single-walled carbon nanotubes computed using the BET and the point B methods. The specific surface area in m2/g (Y axis) is presented as a function of number of carbon atoms in the adsorbate (X axis).
Figure 55.
Fit of the low-pressure region of the butane isotherm data to Equation (
343) in the text (obtained in the MBET approach). The value of the monolayer capacity is extracted from this fit.
Figure 55.
Fit of the low-pressure region of the butane isotherm data to Equation (
343) in the text (obtained in the MBET approach). The value of the monolayer capacity is extracted from this fit.
Figure 56.
Specific surface area values for single-walled carbon nanotubes were calculated using four different methods: BET, point B, procedure A, and procedure B. The results, expressed in m2/g (Y-axis), are plotted as a function of the number of carbon atoms in the adsorbate molecules (X-axis). To maintain clarity in the figure, error bars are shown only for the values obtained using procedure B.
Figure 56.
Specific surface area values for single-walled carbon nanotubes were calculated using four different methods: BET, point B, procedure A, and procedure B. The results, expressed in m2/g (Y-axis), are plotted as a function of the number of carbon atoms in the adsorbate molecules (X-axis). To maintain clarity in the figure, error bars are shown only for the values obtained using procedure B.
Figure 57.
Relative step velocity
as a function of the concentration
for monomers [
, red line, Equation (
355)] and tetramers [
, blue line, Equation (
358)]. In the two cases,
.
Figure 57.
Relative step velocity
as a function of the concentration
for monomers [
, red line, Equation (
355)] and tetramers [
, blue line, Equation (
358)]. In the two cases,
.
Figure 58.
Relative growth rates of the
faces of KBr crystals as a function of impurity concentration. Aliphatic carboxylic acids, HCOOH, CH
3COOH, C
2H
5COOH and C
3H
7COOH were used as impurities. Symbols correspond to experimental data [
171], solid lines represent results from Equations (
355–
358) and dashed lines correspond to KM model [
168]. The parameters used in the theoretical models are listed in
Table 4.
Figure 58.
Relative growth rates of the
faces of KBr crystals as a function of impurity concentration. Aliphatic carboxylic acids, HCOOH, CH
3COOH, C
2H
5COOH and C
3H
7COOH were used as impurities. Symbols correspond to experimental data [
171], solid lines represent results from Equations (
355–
358) and dashed lines correspond to KM model [
168]. The parameters used in the theoretical models are listed in
Table 4.
Figure 59.
Order parameter for , , (dotted); , (dashed); , (solid); , , (dash-dotted). Critical densities for the nematic transition are very approximately, , , , , respectively. All cases correspond to densities of states functions , [ consequently, also the entropy surface ] restricted to vanish at the saturation line in the plane .
Figure 59.
Order parameter for , , (dotted); , (dashed); , (solid); , , (dash-dotted). Critical densities for the nematic transition are very approximately, , , , , respectively. All cases correspond to densities of states functions , [ consequently, also the entropy surface ] restricted to vanish at the saturation line in the plane .
Figure 60.
Order parameter
versus
for
(solid),
;
(dashed),
;
(short-dashed),
;
(dotted),
. Parameters of the empirical high density entropy term
form Equations (??)-(??),
,
,
. The superimposed full symbols for
and
correspond to the analytical results from the density of states function Equation (
373) of the
statistics for:
,
and
,
.
Figure 60.
Order parameter
versus
for
(solid),
;
(dashed),
;
(short-dashed),
;
(dotted),
. Parameters of the empirical high density entropy term
form Equations (??)-(??),
,
,
. The superimposed full symbols for
and
correspond to the analytical results from the density of states function Equation (
373) of the
statistics for:
,
and
,
.
Figure 61.
Full Helmholtz free energy surface for as a function of the state occupation numbers and along the two lattice directions. Black dots represent the equilibrium states. Branches represents the high and low density phases of the nematic transition.
Figure 61.
Full Helmholtz free energy surface for as a function of the state occupation numbers and along the two lattice directions. Black dots represent the equilibrium states. Branches represents the high and low density phases of the nematic transition.
Figure 62.
Lattice coverage
versus
for
. Red and blue lines are theoretical results from the
formalism with the generalized density of states Equation (
373). The high-density (upper blue) and low-density (lower blue) branches, along with the total coverage (red), are shown. The inset highlights the first-order transition at high coverage. Green lines (overlapping) correspond to the empirical entropy correction. Symbols are MC simulation data. Parameters:
,
,
,
,
,
,
,
,
.
Figure 62.
Lattice coverage
versus
for
. Red and blue lines are theoretical results from the
formalism with the generalized density of states Equation (
373). The high-density (upper blue) and low-density (lower blue) branches, along with the total coverage (red), are shown. The inset highlights the first-order transition at high coverage. Green lines (overlapping) correspond to the empirical entropy correction. Symbols are MC simulation data. Parameters:
,
,
,
,
,
,
,
,
.
Figure 63.
Same as
Figure 62 but for
. Parameters:
,
,
,
,
,
,
,
,
.
Figure 63.
Same as
Figure 62 but for
. Parameters:
,
,
,
,
,
,
,
,
.
Figure 64.
Coverage
vs.
for
(main figure) and
(inset). Lines are theoretical
statistics results with generalized density of states Equation (
373). MC data are shown as symbols.
Figure 64.
Coverage
vs.
for
(main figure) and
(inset). Lines are theoretical
statistics results with generalized density of states Equation (
373). MC data are shown as symbols.
Figure 67.
Example of a system of tetramers on a square lattice. (a) Current state. (b) Horizontal tetramers removed. (c) Detail of row j with three segments (only segments with are showed).
Figure 67.
Example of a system of tetramers on a square lattice. (a) Current state. (b) Horizontal tetramers removed. (c) Detail of row j with three segments (only segments with are showed).
Figure 68.
Refilling process of a segment by tetramers. Each branch is made with precalculated probabilities of the 1D case that depend on T and . Red represents a deposited k-mer and grey represents a deposited empty site.
Figure 68.
Refilling process of a segment by tetramers. Each branch is made with precalculated probabilities of the 1D case that depend on T and . Red represents a deposited k-mer and grey represents a deposited empty site.
Figure 69.
Dimers () adsorbed on a square lattice of sites. Typical units are labeled. This represents the system of interest (original).
Figure 69.
Dimers () adsorbed on a square lattice of sites. Typical units are labeled. This represents the system of interest (original).
Figure 70.
a) Square lattice of
sites representing the lattice of the artificial system; strong and weak sites are symbolized by circles and squares respectively. b) Configuration of
dimers in the lowest energy state (ground state) according to the artificial Hamiltonian of Equation (
388).
Figure 70.
a) Square lattice of
sites representing the lattice of the artificial system; strong and weak sites are symbolized by circles and squares respectively. b) Configuration of
dimers in the lowest energy state (ground state) according to the artificial Hamiltonian of Equation (
388).
Figure 71.
Mean total energy per site (in units of the interaction energy w) for dimers on a square lattice with nearest-neighbor interaction energy w at fixed coverage . a) Open circles and top x-axis correspond to attractive dimers. b) Full circles and bottom x-axis correspond to repulsive dimers. Simulations were carried out in the canonical ensemble and symbols represent averages over typically MC configurations, after equilibration steps.
Figure 71.
Mean total energy per site (in units of the interaction energy w) for dimers on a square lattice with nearest-neighbor interaction energy w at fixed coverage . a) Open circles and top x-axis correspond to attractive dimers. b) Full circles and bottom x-axis correspond to repulsive dimers. Simulations were carried out in the canonical ensemble and symbols represent averages over typically MC configurations, after equilibration steps.
Table 1.
Parameters used in the fitting of
Figure 35.
Table 1.
Parameters used in the fitting of
Figure 35.
| |
/
|
/
|
|
(mol. / g) |
0.41 |
0.41 |
|
(mmHg−1) |
3.5 |
0.48 |
|
(mmHg−1) |
1.6 |
0.6 |
Table 2.
Table of fitting parameters of data in
Figure 44 and
Figure 46.
,
,
and
are expressed in kcal/mol (absolute values given).
is expressed in molecules/cavity [(a)] and
g. of adsorbent [(b)] for data from Refs. [
121] and [
122], respectively. (c) and (d) represent experimental values from Refs. [
121] and [
120], respectively; (e), simulation data from Ref. [
123] and (f),
interaction energy in the liquid phase [
120].
Table 2.
Table of fitting parameters of data in
Figure 44 and
Figure 46.
,
,
and
are expressed in kcal/mol (absolute values given).
is expressed in molecules/cavity [(a)] and
g. of adsorbent [(b)] for data from Refs. [
121] and [
122], respectively. (c) and (d) represent experimental values from Refs. [
121] and [
120], respectively; (e), simulation data from Ref. [
123] and (f),
interaction energy in the liquid phase [
120].
| System |
k |
m |
g |
|
|
|
|
|
D(%) |
|
2 |
2 |
4 |
|
3.10 |
|
0.72 |
|
5.60 |
|
3 |
1 |
3 |
5.75 |
6.94 |
|
1.27 |
|
2.08 |
Table 3.
Single-walled nanotubes samples and gases used in the isotherm measurements.
Table 3.
Single-walled nanotubes samples and gases used in the isotherm measurements.
| sample |
type |
weight (g) |
gas |
area/mol on graphite () |
isotherm temp (K) |
| SWNTs |
HiPco |
0.1727 |
methane |
15.4 [145] |
77 |
| SWNTs |
HiPco |
0.325 |
ethane |
21 [146] |
165 |
| SWNTs |
HiPco |
0.325 |
propane |
28.8 [147] |
190 |
| SWNTs |
HiPco |
0.325 |
butane |
32.7 [146] |
220 |
Table 4.
The table shows the parameters used in the theoretical curves of
Figure 58.
Table 4.
The table shows the parameters used in the theoretical curves of
Figure 58.
| |
Kubota and Mullin [168] |
This work |
| Impurity |
|
|
k |
|
|
0.502 |
12.8 |
1 [Equation (355)] |
=3.33 |
|
0.532 |
39.6 |
2 [Equation (356)] |
=2.20 |
|
0.697 |
35.0 |
3 [Equation (357)] |
=3.39 |
|
0.850 |
35.6 |
4 [Equation (358)] |
=5.38 |