# Thermodynamic efficiency of interactions in self-organizing systems

Ramil Nigmatullin  Department of Physics and Astronomy, Macquarie University,
NSW 2109, Australia.  Centre for Complex Systems, Faculty of Engineering, The
University of Sydney, Sydney, NSW 2006, Australia.  Mikhail Prokopenko  Centre
for Complex Systems, Faculty of Engineering, The University of Sydney, Sydney,
NSW 2006, Australia.

(September 2020)

###### Abstract

The emergence of global order in complex systems with locally interacting
components is most striking at criticality, where small changes in control
parameters result in a sudden global re-organization. We introduce a measure
of thermodynamic efficiency of interactions in self-organizing systems, which
quantifies the change in the system’s order per unit work carried out on (or
extracted from) the system. We analytically derive the thermodynamic
efficiency of interactions for the case of quasi-static variations of control
parameters in the exactly solvable Curie-Weiss (fully connected) Ising model,
and demonstrate that this quantity diverges at the critical point of a second
order phase transition. This divergence is shown for quasi-static
perturbations in both control parameters, the external field and the coupling
strength. Our analysis formalizes an intuitive understanding of thermodynamic
efficiency across diverse self-organizing dynamics in physical, biological and
social domains.

Typically, self-organization is defined as a spontaneous formation of spatial,
temporal, spatiotemporal structures or functions in a system comprising
multiple interacting components. Importantly, a self-organizing process is
assumed to be developing in the absence of specific external controls, as
pointed out by Haken haken :

> _a system is self-organizing if it acquires a spatial, temporal or
> functional structure without specific interference from the outside. By
> ‘specific’ we mean that the structure or functioning is not impressed on the
> system, but that the system is acted upon from the outside in a non-specific
> fashion. For instance, the fluid which forms hexagons is heated from below
> in an entirely uniform fashion, and it acquires its specific structure by
> self-organization._

To explain structures that spontaneously self-organize when energy or matter
flows into a system typically describable by many variables, Haken employed
the notion of order parameters (degrees of freedom) and control parameters
synerg ; haken : slowly varying a relevant control parameter, such as
temperature of a ferromagnetic material, may induce an abrupt change, a phase
transition, in an observable order parameter, such as the net magnetization.
The emergence of global order in complex systems is most striking at
criticality, when the characteristic length and dynamical time scales of the
system diverge. A phase transition is usually accompanied by global symmetry
breaking. Crucially, in the more organized (coherent) phase of the system
dynamics, the global behavior of the system can be described by only a few
order parameters, that is, the system becomes low-dimensional as some dominant
variables “enslave” others. A canonical physical example of such coherent
dynamics, when the whole system acts in synchrony, is laser: a beam of
coherent light created out of the chaotic movement of particles haken .

In physical systems, the local interactions are usually determined by physical
laws, e.g., interactions among fluid molecules or crystal ions, while the
interactions within a biological organism may evolve over generations under
environmental selection pressures, bring survival benefits. The role of
locally interacting particles contributing to self-organizing pattern
formation in biological systems has been captured in a definition offered by
Camazine et al. camaz :

> _Self-organization is a process in which pattern at the global level of a
> system emerges solely from numerous interactions among the lower-level
> components of the system. Moreover, the rules specifying interactions among
> the system’s components are executed using only local information, without
> reference to the global pattern._

These definitions concur with many other approaches to formalize self-
organization, highlighting three important aspects bona ; polani2003 ;
Prokopenko2008 : (i) a system dynamically advances to a more organized state,
while exchanging energy, matter and/or information with the environment, but
without a specific external ordering influence; (ii) the interacting system
components have only local information, and so exchange only local
information, but exhibit long-range correlations; (iii) the increase in
organization can be observed as a more coherent global behavior.

In general, as the state of a complex system evolves, its configurational
entropy changes. The reduction (or increase) in the configurational entropy
occurs at the expense of work extracted or carried out on the system, and the
heat exported to the environment. Thus, a thermodynamic analysis of the
interactions in self-organizing systems aims to quantify the work, heat and
energy exchange between the system and the environment. One can reasonably
expect that self-organization is most thermodynamically efficient in the
vicinity of the critical points, i.e., at criticality one may expect that a
smaller amount of work extracted/done on a system can result in a larger
change of the configurational entropy. Indeed, it has been conjectured before
that a system in a self-organized low-dimensional phase with fewer available
configurations (i.e., describable by just a few order parameters and
exhibiting macroscopic stability) may be more efficient than the system in a
high-dimensional disorganized phase with more configurations.

To formalize this conjecture, Kauffman proposed a succinct principle behind
the higher efficiency of self-organized systems — the generation of
constraints during the release of energy — the constrained release channels
energy to perform some useful work, which can propagate and be used again to
create more constraints, releasing further energy and so on Kau00 . Following
a similar characterization, Carteret et al. Carteret  have shown that
available power efficiency is maximized at _critical_ Boolean networks. The
question of thermodynamic efficiency has also been proposed and studied in the
context of the cellular information-processing, from the perspective of how
close life has evolved to approach maximally efficient computation
barato2014efficiency ; Kempes2017 . Furthermore, a recent thermodynamic
analysis of a model of active matter demonstrated that the efficiency of the
collective motion diverges at the transition between disordered and coherent
collective motion Crosato2018Thermodynamics . However, the precise nature of
the divergence of the efficiency of collective motion, and its relation to the
critical exponents describing the system behavior in the vicinity of the phase
transitions remained unclear, due to the lack of analytical expressions for
the corresponding configurational probability distributions.

In this study, we offer a generic measure of thermodynamic efficiency of
interactions within self-organizing systems, aiming to clearly differentiate
between phases of system dynamics, and identify the regimes when the
efficiency is maximal. This measure is expressed by contrasting (i) the change
of organization attained within the system (i.e., change in the created order
or predictability) with (ii) the thermodynamic work involved in driving such a
change. We demonstrate that the maximal efficiency is indeed achieved _at the
critical regime_ , i.e., during the phase transition, rather than at the
macroscopically stable low-dimensional phase _per se_. The reasons for the
maximal efficiency exhibited by systems during self-organization, i.e., at a
critical regime, are articulated precisely in terms of the increased order (or
the reduction of Shannon entropy) related to the amount of the work carried
out during the transition. This measure is defined for specific
configurational changes (perturbations), rather than states or regimes — in
line with the point made by Carteret et al. Carteret  that the maximization of
power efficiency occurs at a finite displacement from equilibrium.

In developing our approach we build on information-theoretic and statistical-
mechanical methods, interpreting the process of self-organization as a
thermodynamic phenomenon, while considering the interactions within the system
as distributed information processing or distributed computation Carteret ;
barato2014efficiency ; Prokopenko-Einav ; Spinney2016-arrow ; Kempes2017 ;
Crosato2018Thermodynamics . Our aim is to develop a common understanding of
thermodynamic efficiency across multiple examples of self-organizing dynamics
in physical, biological and social domains. These phenomena include
transitions from disordered to coherent collective motion Gregoire2004 ;
Buhl2006 ; Vicsek2006 ; Szabo2007 ; mora2011biological ; bialek2012statistical
; Crosato2018Thermodynamics , phase transitions in spin systems and active
matter Montanari ; FTE ; Crosato2019 , chaos-to-order transitions in genetic
regulatory networks modeled as random Boolean networks Carteret ; Wang-Fisher
, synchronization in networks of coupled oscillators near “the edge of chaos”
Kuramoto1987StatisticalMO ; Miritello ; Kalloniatis2018 , transitions across
epidemic thresholds during contagions and cascading failures Newman1999 ;
Sander2002 ; Cupac2013 ; WangLiu2015 ; Harding2018 , critical dynamics of
urban evolution WilsonDearden ; Crosato2018Urban ; Slavko2019 , among many
others. Self-organizing criticality (SOC) SOC  is a related but distinct
phenomenon, as we are not attempting to reveal the mechanisms of self-
organization towards critical regimes, focusing instead on defining and
determining the thermodynamic efficiency of interactions in a representative
self-organizing system.

In this work, we select an abstract statistical-mechanical model (Curie-Weiss
model of interacting spins in a fully connected graph) — one of the simplest
model exhibiting a second-order phase transition — from the widely applicable
mean-field universality class. We analytically evaluate dynamics of this model
in the vicinity of a phase transition, and prove that the thermodynamic
efficiency has a power law divergence at the critical point, and compute its
critical exponent.

##  I Framework

Consider a statistical mechanical system in thermodynamic equilibrium, where
$\textbf{X}=\\{X_{1},...,X_{n}\\}$ are intensive thermodynamic quantities
which act as control parameters that can be changed externally (e.g., magnetic
field). A perturbation in the control parameter,
$\textbf{X}\rightarrow\textbf{X}+\delta\textbf{X}$, will result in a change in
thermodynamic potentials in the system including its entropy and energy. We
define the thermodynamic efficiency of interactions as

$\eta(\textbf{X};\delta\textbf{X})=\frac{1}{k_{B}}\frac{\delta S}{\delta W},$
(1)

where $\delta S$ and $\delta W$ are the change in entropy and the work
done/extracted on the system due to the perturbation $\delta\textbf{X}$.
Entropy $S$ is a configurational entropy, and thus
$\eta(\textbf{X};\delta\textbf{X})$ quantifies the reduction (increase) of
uncertainty in the state of the system that we gain per unit of work done. A
high value of $\eta$ signifies that it is energetically easy to create order
(reduce the configurational uncertainty) in the system by changing a control
parameter, whereas a low value of $\eta$ indicates that a lot of work is
needed to change the order in the system.

In practice, to evaluate $\eta(\textbf{X};\delta\textbf{X})$ we need to
specify the perturbation protocol. A change in control parameters moves the
system out of thermal equilibrium, and we need to compute the amount of work
done/extracted, $\delta W$, as the system relaxes back to its equilibrium
state. Thus, $\eta(\textbf{X};\delta\textbf{X})$ depends on how we perturb the
system, and on the master equation that describes the relaxation of the system
back to its equilibrium state. In what follows we will consider the case of a
quasi-static perturbation protocol, i.e., we assume that the perturbation is
sufficiently slow that the system effectively adjusts instantaneously to its
new equilibrium state. Helmholtz free energy, $F(\theta,\textbf{X})$, is the
most useful thermodynamic potential for analyzing the quasi-static protocols
at constant temperature. Helmholtz free energy is related to the internal
energy $U$ and entropy $S$ via equation

$U(\theta,\textbf{X})=\theta S(\theta,\textbf{X})+F(\theta,\textbf{X}),$

where $\theta\equiv k_{B}$T. To a first order in $\delta\textbf{X}$ the change
in internal energy, entropy and free energy induced by varying the control
parameters are $\delta U=\delta\textbf{X}\cdot\left.\nabla
U\right|_{\textbf{X}}$, $\delta S=\delta\textbf{X}\cdot\left.\nabla
S\right|_{\textbf{X}}$ and $\delta F=\delta\textbf{X}\cdot\left.\nabla
F\right|_{\textbf{X}}$. In a quasi-static process the change in free energy
can be identified with the work done on the system, $\delta F=\delta W$, and
the entropy change in the system balances the entropy exported to the
environment, $\delta S=-\delta S^{\textrm{exp}}$. Thus, for a quasi-static
protocol, the thermodynamic efficiency reduces to

$\eta(\textbf{X};\delta\textbf{X})=\frac{1}{k_{B}}\frac{\delta\textbf{X}\cdot\left.\nabla
S\right|_{\textbf{X}}}{\delta\textbf{X}\cdot\left.\nabla
F\right|_{\textbf{X}}}.$ (2)

In the case when the variation of control parameter is one-dimensional
$\textbf{X}=X$, equation (2) simplifies to

$\displaystyle\eta(X,\delta X)$
$\displaystyle=\left.\frac{1}{k_{B}}\frac{\partial S}{\partial
X}\middle/\frac{\partial F}{\partial X}\right.$ (3)
$\displaystyle=\frac{1}{k_{B}}\frac{\partial S}{\partial F}.$

Equation (3) applies for general quasistatic processes. When the system is
close to a critical point of a phase transition, expression for $\eta$ can
further be simplified using the following argument. Let $\psi$ be an extensive
quantity conjugate to $X$

$\psi=-\frac{\partial F}{\partial X}.$ (4)

Entropy is related to free energy via $S=-\partial F/\partial T$, and thus the
derivative of $S$ with respect to $X$ is

$\frac{\partial S}{\partial X}=-\frac{\partial^{2}F}{\partial X\partial
T}=-\frac{\partial^{2}F}{\partial T\partial X}=\frac{\partial\psi}{\partial
T}.$

Thus, in terms of the extensive variable conjugate to the control parameter
the efficiency given by equation (3) can be expressed as

$\eta(X,\delta X)=-\frac{1}{k_{B}}\frac{1}{\psi}\frac{\partial\psi}{\partial
T}$ (5)

If $\psi$ is an order parameter of a phase transition then near the critical
point we have $\psi=a|T-T_{c}|^{\beta}$, where $T_{c}$ is the critical
temperature, $\beta$ is the critical exponent and $a$ is non-universal
proportionality constant. Upon substitution of this expression for $\psi$ into
(5), the constant $a$ cancels and we obtain

$\eta(X,\delta X)=-\frac{1}{k_{B}}\frac{\beta}{|T-T_{c}|}.$ (6)

Equation (6) expresses the divergence of $\eta$ solely in terms of universal
exponent $\beta$. This result explains why in many thermodynamic models the
efficiency of self-organization is expected to peak near the critical point.

In many complex systems, there may not exist a readily available physical
model expressed in terms of a Hamiltonian and the expression for the order
parameter may not be evident. Nevertheless, if there is a record of samples of
the states of the system then one may still use equation (3) to compute the
efficiency of self-organization. The reason is that all of the thermodynamic
quantities in (3), expressed in terms of Gibbs probability distribution, have
a clear information-theoretic interpretation. Entropy $S$ is directly
proportional to the Shannon entropy $H$, $S=-k_{B}\sum_{x}p(x)\log
p(x)=k_{B}H$. The free energy $F$ is related to the Fisher information
$\mathcal{I}$ via equation $\mathcal{I}=\partial F^{2}/\partial^{2}X$, with
the Fisher information quantifying the sensitivity of the probability
distribution to the change in the control parameter,
$\mathcal{I}\equiv\sum_{x}(\partial\log p(x)/\partial X)^{2}p(x)$.

There are several interpretations of the Fisher information relevant to
critical dynamics and scale dependence: $\mathcal{I}$ is equivalent to the
thermodynamic metric tensor, a curvature of which diverges at phase
transitions; $\mathcal{I}$ measures the size of the fluctuations in the
collective variables around equilibrium; also, $\mathcal{I}$ is proportional
to the derivatives of the corresponding order parameters with respect to the
collective variables brody1995geometrical ; brody2003information ;
janke2004information ; crooks2007measuring ; Wang-Fisher ;
prokopenko2011relating ; Machta2013 . Substituting $\partial F/\partial
X=\int\mathcal{I}dX$ into equation (3) gives

$\eta(X)=\frac{\partial H/\partial X}{\int\mathcal{I}dX}.$ (7)

Equation (7) expresses the thermodynamic efficiency of interaction during
configurational perturbations in terms of information-theoretic quantities of
entropy and Fisher information.

Equation (7) was derived and used in Crosato2018Thermodynamics  in the
thermodynamic analysis of collective motion (e.g., swarming) exhibiting a
kinetic phase transition. Crosato et al. Crosato2018Thermodynamics  computed
the efficiency $\eta$ from the distribution $p(x)$, estimated via sampling
produced by numerical simulations of the model, consequently yielding
estimates of $H$ and $\mathcal{I}$. It was then demonstrated that $\eta$
diverges at the critical point where the swarm transitions from disordered to
coherent motion.

The notion of thermodynamic efficiency $\eta$ was also applied to the analysis
of urban transformations Crosato2018Urban , driven by quasi-static changes in
the social disposition: a control parameter characterising the attractiveness
of different areas. The thermodynamic efficiency of urban transformations was
defined as the reduction of configurational entropy resulting from the
expenditure of work. In the socioeconomic context of urban dynamics, it
expressed the ratio of the gained predictability of income flows to the amount
of work required to change the social disposition. Importantly, the efficiency
was shown to peak at a critical transition separating dispersed and
polycentric phases of urban dynamics Crosato2018Urban .

Similarly, Harding et al. Harding2018  considered thermodynamic efficiency of
quasi-static epidemic processes, defined for a value of some control parameter
(e.g., the infection transmission rate), as the ratio of the reduction in
uncertainty to the expenditure of work needed to change the parameter. On the
one hand, this could be the efficiency of an intervention process consuming
work in order to reduce the transmission rate. On the other hand, the
efficiency can be defined in terms of the pathogen emergence — a process which
increases the transmission rate, and in doing so extracts the work.
Irrespective of the interpretation, the efficiency was shown to peak at the
epidemic threshold Harding2018 .

Our contribution builds on this research, showing that according to equation
(6), the divergence of the efficiency of self-organization is generally
expected to occur at a second-order phase transition. In the following
section, we illustrate this result by explicitly computing $\eta$ in a simple
model exhibiting paramagnetic to ferromagnetic phase transition, showing that
the efficiency of self-organization peaks at the critical point when the
control parameter is either the coupling strength between the spins or the
external magnetic field.

##  II Example: Curie-Weiss Model

The energy of a wide variety statistical mechanical systems, including spin
glasses, can be written in the following form

$E(\underline{\sigma},\\{X_{i}\\})=E_{0}+\sum_{i}X_{i}E_{i}(\underline{\sigma}),$
(8)

where $\\{X_{i}\\}=\\{X_{1},X_{2},...,X_{K}\\}$ are the control parameters of
the system, $\underline{\sigma}$ denotes the microscopic state of the system
and $\theta\equiv k_{B}T$. We are working with a canonical ensemble, where the
system is in contact with a heat bath in thermal equilibrium and the average
energy is fixed. In this case, the probability of finding the system in
configuration $\underline{\sigma}$ is given by the Gibbs measure

$p(\underline{\sigma};\\{X_{i}\\})=\frac{e^{-E(\underline{\sigma},\\{X_{i}\\})/\theta}}{Z(\theta,\\{X_{i}\\})},$
(9)

where $Z=\sum_{\underline{\sigma}}e^{-E(\underline{\sigma})}/\theta$ is the
partition function. The free energy of the system is given by $F=\ln Z$. The
free energy The free energy can be used to compute any thermodynamic quantity,
in particular, the expectation of $\langle E_{1}\rangle,...,\langle
E_{K}\rangle$ are given by

$\langle E_{j}\rangle=-\frac{\partial F}{\partial X_{j}}.$ (10)

For an interacting statistical mechanical system in thermal equilibrium there
is a one to one map between the the set of control parameters
$\\{T,X_{1},...,X_{K}\\}$ and $\\{S,\langle E_{1}\rangle,...\langle
E_{K}\rangle\\}$ Wei2016 , and thus we will refer to $\psi_{i}\equiv\langle
E_{i}\rangle$ as a order parameter conjugate to the control parameter $X_{i}$.
Phase transitions are often accompanied by divergences in one or more order
parameters $\psi_{i}$ or their derivatives.

In the rest of the paper, we will focus on computing $\eta$ for a system
governed by a specific energy function of the form (8) — the Curie-Weiss (CW)
model. The CW model is a model of ferromagnetism, where each spin interacts
with all other spins via pairwise interactions and for this reason it is also
known as the fully connected Ising model. This model exhibits a second-order
phase transition at a finite critical temperature $T_{c}$. In the vicinity of
the critical point, the analytic expression to all of the thermodynamic
quantities are known, which enables the derivation of the analytic expression
for $\eta$. The phase transition from ferromagnetic to paramagnetic states in
the Curie-Weiss model belongs to the mean field universality class.

Let $N$ spins $\sigma_{i}\in\\{\pm 1\\}$ be assigned to sites
$i\in\\{1...N\\}$. A configuration of the system is given by
$\underline{\sigma}=(\sigma_{1},...,\sigma_{N})$. The energy function for the
system containing pairwise interactions between spins and in the presence of
an external magnetic field $B$ is given by

$E(\underline{\sigma})=-\frac{J}{N}\sum_{(ij)}\sigma_{i}\sigma_{j}-\mu
B\sum_{i=1}^{N}\sigma_{i},$ (11)

where the sum over $(ij$) runs over all of the $N(N-1)/2$. The $1/N$ scaling
in front of the spin-spin interaction term is to yield an extensive free
energy. In this model, the control parameters are $\\{J,B\\}$, which
respectively denote exchange interaction strength and externally applied
magnetic field. The probability of finding the system in configuration
$\underline{\sigma}$ is given by the Gibbs measure

$p(\underline{\sigma};T,B,J)=\frac{e^{-E(\underline{\sigma})/\theta}}{Z_{N}(\theta,B,J)},$
(12)

where $\theta\equiv k_{B}T$ and $Z_{N}$ is a partition function for the
$N$-spin system. The free energy of the $N$ spin system is given by
$F_{N}(\theta,B)=\ln Z_{N}(\theta,B)$. The thermodynamic limit is obtained by
taking $N\rightarrow\infty$. In the thermodynamic limit the free energy
density $f(\theta,B)=\lim_{N\rightarrow\infty}F_{N}(\theta,B)/N$ can have the
following analytic expression CWmodel_paper :

$f(\theta,B)=-\theta\ln 2-\theta\ln(\Phi(\theta,B))$ (13)

with

$\Phi(\theta,B)=e^{-Jy^{2}/(2\theta)}\cosh\left(\frac{Jy+B}{\theta}\right)$
(14)

Here $y$ is defined as solution to the equation

$y=\tanh\left(\frac{Jy+\mu B}{\theta}\right).$ (15)

The average magnetization per spin is the order parameter conjugate to
magnetic field and is given by $m=-(\partial f/\partial B)_{\theta}=\mu y$,
and thus the equation of state is $m=\mu\tanh[(Jm+B\mu)/\theta\mu]$. The phase
diagram can be constructed by analyzing the equation of state. The critical
point of a second-order phase transition occurs at $B=0$ and $\theta_{c}=J$.
When $B=0$ and $\theta>J$ there is only one stable solution of the equation of
state, which is $m=0$. When $B=0$ and $\theta<J$ , there are three solutions:
one unstable solution $m=0$ and two stable solution $m=\pm m^{*}$ where
$m^{*}$ is found by numerically solving the equation $m=\mu\tan(Jm/\theta)$.
Thus at $B=0$ and at the critical temperature $\theta_{c}=J$, the system
transitions from a paramagnetic disordered state where $m=0$ to a
ferromagnetic ordered state where $m=\pm m^{*}$. This transition is of second
order, since the second derivatives of $f$ with respect to both $B$ and
$\theta$ are discontinuous at $\theta_{c}$.

Having reviewed the phase change behavior of the CW model, we will now
evaluate the thermodynamic efficiency $\eta$ associated with varying the
magnetic field $B$ along a quasi-static protocol. The entropy density is
related to the free energy density via equation

$s=-\frac{\partial f(y(\theta,B),\theta,B)}{\partial\theta}.$ (16)

Using the equations (13)-(16) and (2), one can compute efficiency of self-
organization $\eta$ resulting from variation one or more control parameters
$B$, $J$ or $\theta$.

###  II.1 Varying external field, $B$

Since equation (15) does not have a closed form solution for $y(\theta,B)$, it
has to be solved numerically. Thus, for a general choice of parameters of the
Curie-Weiss model, the efficiency $\eta(\theta;\delta B)$ needs to be
evaluated numerically. The plots of derivatives of free energy and entropy
densities computed numerically by solving equation (15) are shown in Figure 1.
The thermodynamic efficiency is the ratio of these two derivatives,
$\eta=\frac{1}{k_{B}}\frac{\partial s}{\partial B}/\frac{\partial f}{\partial
B}$, which is plotted in Figure 2. As expected, the efficiency $\eta$ peaks
near the critical point of the phase transition. In the rest of this section,
we focus on behaviour of $\eta$ in the vicinity of the critical point, where
it is possible to obtain an analytic solution for all thermodynamic quantities
and study their scaling behaviour.

Near the critical point of the paramagnetic to ferromagnetic phase transition,
the average magnetization $y/\mu$ is small and thus equation of state (15) can
be approximated by a low order Taylor expansion in $y$. Keeping up to
$O(y^{3})$ the equation of state is

$K^{3}y^{3}-3y(K-1)-3h=0,$ (17)

where $K\equiv J/\theta=\theta_{c}/\theta$, $h=\mu B/\theta$. In the case of
zero magnetic field, $h=0$, the solution of (17) is

$\displaystyle y$ $\displaystyle=0,\qquad\textrm{for }t\geq 0,$ (18)
$\displaystyle=\pm\sqrt{\frac{3(K-1)}{K^{3}}}\sim\sqrt{3}\left(-t\right)^{1/2}\qquad\textrm{for
}t<0$

where $t$ is the reduced temperature $t\equiv(\theta-\theta_{c})/\theta_{c}$
and $h\equiv\mu B/\theta$. Equation (17) produces the well-known mean field
scaling law for magnetization $m\sim(-t)^{\beta}$ for $t<0$, with the critical
exponent $\beta=1/2$. Use equation (6) we arrive at $\eta(\theta,\delta
B)=-\frac{1}{2k_{B}}\frac{1}{t}$ for $t<0$.

Figure 1: Derivatives of entropy and free energy as a function of temperature
$\theta$ at zero magnetic field. The inset shows how the presence of a small
magnetic field smooths out the singularity in $\partial S/\partial B$ at the
critical point $\theta_{c}=J=1.$

In the paramagnetic case $t\geq 0$, $y=0$, both $\partial f/\partial B$ and
$\partial s/\partial B$ are zero and consequently the efficiency appears to be
undefined as it is a ratio of these derivatives. Nevertheless, in the
paramagnetic regime the derivatives of free energy and entropy can be made
finite by either adding a small external magnetic field or by considering a
finite size system. Here, we will consider the efficiency $\eta$ in the
presence of constant magnetic field $B_{0}$, which can be made arbitrary
small. In the presence of external field and when $t\gg 0$, the equation of
state (17) simplifies to $y(1-K)+h=0$, since the term $K^{3}y^{3}$ is
negligible. Thus, in this regime, $y\sim h/(1-K)=\mu B/(\theta-\theta_{c})$
and $f$ can now be evaluated using equations (13) and (14). From $f$ we
compute $\partial_{B}f$, $\partial_{B}s$ and then Taylor expand to the leading
order in $B$ to obtain

$\displaystyle\left.\frac{\partial f}{\partial B}\right|_{B=B_{0}}$
$\displaystyle=-\frac{B_{0}}{\theta_{c}t}\qquad\textrm{for }t\geq 0$ (19)
$\displaystyle\left.\frac{\partial s}{\partial B}\right|_{B=B_{0}}$
$\displaystyle=-\left.\frac{\partial f}{\partial
B\partial\theta}\right|_{B=B_{0}}=-\frac{B_{0}}{\theta_{c}^{2}t^{2}}\qquad\textrm{for
}t\geq 0.$ (20)

Now, we can evaluate scaling behavior of the thermodynamic efficiency $\eta$
around the critical point:

$\displaystyle\eta(\theta,\delta B)$
$\displaystyle=\left.\frac{1}{k_{B}}\frac{\partial s}{\partial
B}\middle/\frac{\partial f}{\partial B}\right.$ $\displaystyle=$
$\displaystyle\begin{cases}-\frac{1}{k_{B}}\frac{1}{2}t^{-1}&\qquad\textrm{for
}t<0\\\ \frac{1}{k_{B}\theta_{c}}t^{-1}&\qquad\textrm{for }t>0.\end{cases}$
(21)

A plot of $\eta$ in the vicinity of the critical point for several small
values of bias field $B_{0}$ is shown in figure 2. The curves were obtained by
numerically solving for $y$ and numerically computing the derivative of $f$
and $s$. The $|t|^{-1}$ scaling prediction agrees very well with the numerical
results. The deviations at finite $B_{0}$ and very close to the critical point
are expected as the scaling was obtained by neglecting $K^{3}y^{3}$ term in
equation (17), which is not small around $\theta=\theta_{c}$.

Figure 2: Thermodynamic efficiency $\eta(\theta,\delta B)$ as a function of
$\theta$ at several small values of $B$. The critical point is at
$\theta_{c}=1.0$ or equivalently at $t\equiv(\theta-\theta_{c})/\theta_{c}=0$.
For $t>0$, $\eta$ is undefined at $B=0$. The solid lines $-1/2t^{-1}$ for
$t<0$ and $t^{-1}$ for $t>0$ are analytic expressions for $\eta$ in the
vicinity of the critical point.

###  II.2 Varying coupling strength, $J$

We now consider computing $\eta$, when $J$ is used as a control parameter. In
this case, the relevant order parameter is
$\psi_{J}\equiv\langle\sum_{ij}\sigma_{i}\sigma_{j}\rangle$, which quantifies
the interaction energy between pairs of spins. The spins will spontaneously
align at a critical value of coupling strength $J_{c}=\theta$, and the
efficiency $\eta(J,\delta J)$ is expected to peak near the critical point.

Near the critical point there is a closed form expression for $y$, and thus we
can derive the scaling relation between $\eta$ and the reduced coupling
strength $\mathcal{J}\equiv(J-J_{c})/J_{c}$. For the ferromagnetic case,
$J>J_{c}$, inserting equation (18) into the expressions for the free energy
and entropy, taking derivatives with respect to $\mathcal{J}$ and then Taylor-
expanding in $\mathcal{J}$ to the lowest orders yields

$\displaystyle\frac{\partial f}{\partial\mathcal{J}}$
$\displaystyle=-\frac{3\mathcal{J}}{2\theta}\qquad\textrm{for }\mathcal{J}>0$
(22) $\displaystyle\frac{\partial s}{\partial\mathcal{J}}$
$\displaystyle=-\frac{3}{2\theta}\qquad\textrm{for }\mathcal{J}>0$ (23)

The order parameter conjugate to $\mathcal{J}$ can be defined as
$\phi\equiv-\partial f/\partial\mathcal{J}$, which according to equation (22)
is linearly proportional to $\mathcal{J}$, i.e. $\phi\sim\mathcal{J}^{\beta}$
with $\beta=1$. Using equation (6) with critical exponent $\beta=1$, we
immediately arrive at
$\eta(\theta,\delta\mathcal{J})=\frac{1}{k_{B}}\frac{1}{\mathcal{\mathcal{J}}}$
for $\mathcal{J}>0$.

In the paramagnetic case, $\mathcal{J}<0,$ the magnetization is zero in the
absence of the external magnetic field and the efficiency of interactions is
undefined since $\partial f/\partial\mathcal{J}=0$. However, in the presence
of small bias magnetic field $B_{0}$, $\eta$ can be computed, since in that
case $y\sim h/(1-K)=\mu B_{0}/(\theta-J)$. Taylor-expanding $\partial
f/\partial\mathcal{J}$ and $\partial s/\partial\mathcal{J}$ computed with this
expression for $y$ to the lowest orders in $\mathcal{J}$ gives

$\displaystyle\left.\frac{\partial f}{\partial\mathcal{J}}\right|_{B=B_{0}}$
$\displaystyle=-\frac{B_{0}^{2}\mu^{2}}{\mathcal{J}^{3}}\qquad\textrm{for
}\mathcal{J}<0$ (24) $\displaystyle\left.\frac{\partial
s}{\partial\mathcal{J}}\right|_{B=B_{0}}$
$\displaystyle=\frac{2B_{0}^{2}\mu^{2}}{\mathcal{J}^{4}}\qquad\textrm{for
}\mathcal{J}<0.$ (25)

Using equations (22)-(25) we can compute the efficiency of interactions in the
vicinity of the critical point:

$\displaystyle\eta(J,\delta J)$
$\displaystyle=\left.\frac{1}{k_{B}}\frac{\partial
s}{\partial\mathcal{J}}\middle/\frac{\partial f}{\partial\mathcal{J}}\right.$
$\displaystyle=$
$\displaystyle\begin{cases}\frac{1}{k_{B}}\mathcal{J}^{-1}&\qquad\textrm{for
}\mathcal{J}>0\\\ -\frac{1}{k_{B}}2\mathcal{J}^{-1}&\qquad\textrm{for
}\mathcal{J}<0.\end{cases}$ (26)

Figure 3: Thermodynamic efficiency $\eta(J,\delta J)$ as a function of $J$ at
several small values of $B$ at $\theta=1.0$. The critical point is at
$J_{c}=1.0$ or equivalently at $\mathcal{J}\equiv(J-J_{c})/J_{c}=0$. For
$\mathcal{J}<0$, $\eta$ is undefined at $B=0$. The solid lines
$-2\mathcal{J}^{-1}$ for $\mathcal{J}<0$ and $\mathcal{J}^{-1}$ for
$\mathcal{J}>0$ are analytic expressions for $\eta$ in the vicinity of the
criticality.

Figure 3 shows the plot of $\eta(J,\delta J)$ in the vicinity of the critical
point for several small values of bias field $B_{0}$. The dotted curves were
obtained by numerically solving for $y$ and numerically computing the
derivative of $f$ and $s$. The solid black lines indicate the
$|\mathcal{J}|^{-1}$ scaling, which agrees very well with the numerical
results.

##  III Conclusions

The increasing interest in developing a comprehensive thermodynamic framework
for studying complex system, including the process of self-organization, is
driven by several recent developments: theoretical advances in stochastic
thermodynamics Seifert2012  which enable rigorous quantitative analysis of
small and mesoscale systems; technological advances that enable measurement of
thermodynamic quantities of such systems Berut2015 ; Nisoli2018 ; SpinIce2018
; and a fusion of information-theoretic, computation-theoretic and
statistical-mechanical approaches for analyzing energy-efficiency of
information processing devices Wolpert2019 .

We have introduced a measure of thermodynamic efficiency of interactions in
self-organizing systems, which quantifies the change in the order in the
system per unit of the work done/extracted due to the changes in control
parameters. We have shown that this quantity peaks at the critical regime, by
explicitly deriving it for the exactly solvable Curie-Weiss model — a
paradigmatic model of second-order phase transitions. Quasi-static
perturbation in both control parameters, the interaction strength between
spins and the externally applied magnetic field, have been considered, and
both protocols have been shown to lead to divergence of the efficiency of
interactions at criticality.

Our work aims to support systematic thermodynamic studies of self-organization
in complex systems. In the future, it would be interesting to examine models
exhibiting other types of phase transitions in a broad class of dynamical
systems, including econo- and socio-computation Bouchaud2013 ; Cancho2003 ;
Cancho2007 ; Prokopenko2010 ; Salge2015 ; Liar2019 , and determine whether
they also maximally efficient at criticality. It will also be important to
extend the analysis to the protocols that drive the system out-of-equilibrium.
We believe that an approach to self-organization incorporating the
thermodynamic efficiency will also help in clarifying the fundamental
relationship between the structure of a complex system and its collective
behavior and function Newman2003 , as well as support efforts to
systematically control and guide the dynamics of complex systems
RevModPhys.88.035006 ; Daniels2017 .

## References

  * [1] H. Haken.  Information and Self-Organization: A Macroscopic Approach to Complex Systems.  Springer-Verlag, 1988.
  * [2] H. Haken.  Synergetics, an Introduction: Nonequilibrium Phase Transitions and Self-Organization in Physics, Chemistry, and Biology.  Springer-Verlag, 3rd rev. enl. ed., New York, 1983.
  * [3] S. Camazine, J.-L. Deneubourg, N. R. Franks, J. Sneyd, G. Theraulaz, and E. Bonabeau.  Self-Organization in Biological Systems.  Princeton University Press, 2001.
  * [4] E. Bonabeau, G. Theraulaz, J.-L. Deneubourg, and S. Camazine.  Self-organisation in social insects.  Trends in Ecology and Evolution, 12(5):188–193, 1997.
  * [5] D. Polani.  Measuring self-organization via observers.  In W. Banzhaf, T. Christaller, P. Dittrich, J. T. Kim, and J. Ziegler, editors, Advances in Artificial Life – Proceedings of the 7th European Conference on Artificial Life (ECAL), Dortmund, pages 667–675. Springer-Verlag, Heidelberg, 2003.
  * [6] M. Prokopenko.  Design vs. Self-organization, pages 3–21.  Springer London, London, 2nd edition, 2008.
  * [7] S. A. Kauffman.  Investigations.  Oxford University Press, Oxford, 2000.
  * [8] H. Carteret, K. Rose, and S. Kauffman.  Maximum power efficiency and criticality in random boolean networks.  Physical Review Letters, 101:218702, 12 2008.
  * [9] A.C. Barato, D. Hartich, and U. Seifert.  Efficiency of cellular information processing.  New Journal of Physics, 16(10):103024, 2014.
  * [10] C. Kempes, D. Wolpert, Z. Cohen, and Juan Pérez-Mercader.  The thermodynamic efficiency of computations made in cells across the range of life.  Philosophical Transactions of The Royal Society A Mathematical Physical and Engineering Sciences, 375, 06 2017.
  * [11] E. Crosato, R. E. Spinney, R. Nigmatullin, J. T. Lizier, and M. Prokopenko.  Thermodynamics and computation during collective motion near criticality.  Physical Review E, 97:012120, 2018.
  * [12] M. Prokopenko and I. Einav.  Information thermodynamics of near-equilibrium computation.  Physical Review E, 91:062143, Jun 2015.
  * [13] R. E. Spinney, J. T. Lizier, and M. Prokopenko.  Transfer entropy in physical systems and the arrow of time.  Physical Review E, 94:022135, 2016.
  * [14] G. Grégoire and H. Chaté.  Onset of collective and cohesive motion.  Physical Review Letters, 92(2):025702, 2004.
  * [15] J. Buhl, D. Sumpter, I. Couzin, J.J. Hale, E. Despland, E. Miller, and S.J. Simpson.  From disorder to order in marching locusts.  Science, 312:1402–6, 07 2006.
  * [16] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Sochet.  Novel type of phase transition in a system of self-driven particles.  Physical Review Letters, 75:1226, 12 2006.
  * [17] B. Szabo, G. Szollosi, B. Gönci, Z. Jurányi, D. Selmeczi, and T. Vicsek.  Phase transition in the collective migration of tissue cells: Experiment and model.  Physical review. E, Statistical, nonlinear, and soft matter physics, 74:061908, 01 2007.
  * [18] T. Mora and W. Bialek.  Are biological systems poised at criticality?  Journal of Statistical Physics, 144(2):268–302, 2011.
  * [19] W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A.M. Walczak.  Statistical mechanics for natural flocks of birds.  Proceedings of the National Academy of Sciences, 109(13):4786–4791, 2012.
  * [20] M. Mezard and A. Montanari.  Information, Physics, and Computation.  Oxford University Press, Inc. (New York), 2009.
  * [21] M. Prokopenko, L. Barnett, M. Harré, J. Lizier, O. Obst, and X. R. Wang.  Fisher transfer entropy: Quantifying the gain in transient sensitivity.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science, 471:20150610, 12 2015.
  * [22] E. Crosato, M. Prokopenko, and R. E. Spinney.  Irreversibility and emergent structure in active matter.  Physical Review E, 100:042613, 2019.
  * [23] X.R. Wang, J.T. Lizier, and M. Prokopenko.  Fisher information at the edge of chaos in random Boolean networks.  Artificial Life, 17(4):315–329, 2011.
  * [24] Yoshiki K. and Ikuko N.  Statistical macrodynamics of large dynamical systems. case of a phase transition in oscillator communities.  Journal of Statistical Physics, 49:569–605, 1987.
  * [25] G. Miritello, A. Pluchino, and A. Rapisarda.  Central limit behavior in the Kuramoto model at the ‘edge of chaos’.  Physica A: Statistical Mechanics and its Applications, 388:4818–4826, 12 2009.
  * [26] A. Kalloniatis, M. Zuparic, and M. Prokopenko.  Fisher information and criticality in the Kuramoto model of nonidentical oscillators.  Physical Review E, 98:022302, 08 2018.
  * [27] M. Newman and D. Watts.  Scaling and percolation in the small-world network model.  Physical Review E, 60:7332–42, 01 1999.
  * [28] L. Sander, C. Warren, I. Sokolov, C. Simon, and J. Koopman.  Percolation on heterogeneous networks as a model for epidemics.  Mathematical Biosciences, 180:293–305, 11 2002.
  * [29] V. Cupac, J. Lizier, and M. Prokopenko.  Comparing dynamics of cascading failures between network-centric and power flow models.  International Journal of Electrical Power and Energy Systems, 49:369–379, 2013.
  * [30] W. Wang, Q.-H. Liu, L.-F. Zhong, M. Tang, H. Gao, and H. Stanley.  Predicting the epidemic threshold of the susceptible-infected-recovered model.  Scientific Reports, 6:24676, 12 2015.
  * [31] N. Harding, R. Nigmatullin, and M. Prokopenko.  Thermodynamic efficiency of contagions: A statistical mechanical analysis of the SIS epidemic model.  Interface Focus, 8:20180036, 12 2018.
  * [32] A. Wilson and J. Dearden.  Phase transitions and path dependence in urban evolution.  Journal of Geographical Systems, 13:1–16, 2011.
  * [33] E. Crosato, R. Nigmatullin, and M. Prokopenko.  On critical dynamics and thermodynamic efficiency of urban transformations.  Royal Society Open Science, 5(10):180863, 2018.
  * [34] B. Slavko, K. Glavatskiy, and M. Prokopenko.  Dynamic resettlement as a mechanism of phase transitions in urban configurations.  Physical Review E, 99:042143, 04 2019.
  * [35] P. Bak, C. Tang, and K. Wiesenfeld.  Self-organized criticality: An explanation of the 1/f noise.  Physical Review Letters, 59(4):381–384, 1987.
  * [36] D.C. Brody and N. Rivier.  Geometrical aspects of statistical mechanics.  Physical Review E, 51:1006–1011, 1995.
  * [37] D.C. Brody and A. Ritz.  Information geometry of finite Ising models.  Journal of Geometry and Physics, 47(2):207–220, 2003.
  * [38] W. Janke, D.A. Johnston, and R. Kenna.  Information geometry and phase transitions.  Physica A: Statistical Mechanics and its Applications, 336(1–2):181–186, 2004.
  * [39] G.E. Crooks.  Measuring thermodynamic length.  Physical Review Letters, 99:100602, 2007.
  * [40] M. Prokopenko, J.T. Lizier, O. Obst, and X.R. Wang.  Relating Fisher information to order parameters.  Physical Review E, 84:041116, 2011.
  * [41] B. B. Machta, R. Chachra, M. K. Transtrum, and J. P. Sethna.  Parameter space compression underlies emergent theories and predictive models.  Science, 342(6158):604–607, 2013.
  * [42] B.-B. Wei.  New J. Phys., 18:113035, 2016.
  * [43] M. Kochmanski, T. Paszkiewicz, and S. Wolski.  Curie-Weiss magnet — a simple model of phase transition.  European Journal of Physics, 34:1555, 2013.
  * [44] U. Seifert.  Stochastic thermodynamics, fluctuation theorems, and molecular machines.  Reports on Progress in Physics, 75:126001, 2012.
  * [45] Antoine Bérut, Artyom Petrosyan, and Sergio Ciliberto.  Information and thermodynamics: Experimental verification of landauer’s erasure principle.  Journal of Statistical Mechanics Theory and Experiment, 2015:P06015, 2015.
  * [46] Cristiano Nisoli.  Write it as you like it.  Nature Nanotechnology, 13:5–6, 2018.
  * [47] Yuyang L., Francesco C., Mohammed S., Joseph S., Daniel G., J. D. Watts, A. M. Albrecht, Andreas S., Karin D., Cristiano N., and Peter S.  Classical topological order in the kinetics of artificial spin ice.  Nature Physics, 14(7):723–727, 2018.
  * [48] D. H. Wolpert.  The stochastic thermodynamics of computation.  Journal of Physics A, 52:193001, 2019.
  * [49] Jean-Philippe Bouchaud.  Crises and collective socio-economic phenomena: Simple models and challenges.  J Stat Phys, 2013.
  * [50] R. Ferrer-i Cancho and R. Sole.  Least effort and the origins of scaling in human language.  Proceedings of the National Academy of Sciences of the United States of America, 100:788–91, 03 2003.
  * [51] R. Ferrer-i Cancho and A. Diaz-Guilera.  The global minima of the communicative energy of natural communication systems.  Journal of Statistical Mechanics: Theory and Experiment, 2007:P06009, 06 2007.
  * [52] M. Prokopenko, N. Ay, O. Obst, and D. Polani.  Phase transitions in least-effort communications.  Journal of Statistical Mechanics: Theory and Experiment, 2010:P11025, 11 2010.
  * [53] C. Salge, N. Ay, D. Polani, and M. Prokopenko.  Zipf’s law: Balancing signal usage cost and communication efficiency.  PLoS ONE, 10:e0139475, 10 2015.
  * [54] M. Prokopenko, M. Harré, J. Lizier, F. Boschetti, P. Peppas, and S. Kauffman.  Self-referential basis of undecidable dynamics: from the Liar paradox and the halting problem to the edge of chaos.  Physics of Life Reviews, 2019.
  * [55] M. E. J. Newman.  The structure and function of complex networks.  SIAM Review, 45:167–256, 2003.
  * [56] Y.-Y. Liu and A.-L. Barabási.  Control principles of complex systems.  Review of Modern Physics, 88:035006, Sep 2016.
  * [57] Bryan C. Daniels, David C. Krakauer, and Jessica C. Flack.  Control of finite critical behaviour in a small-scale social system.  Nature Communications, 8(1):14301, 2017.


# Effects of charge and size on the coadsorption of counterionic colloids in
Gibbs monolayers

J. M. Gómez-Verdú  $|$ , B. Martínez-Haya  $|$ , A. Cuetos  $|$ Corresponding
author<EMAIL_ADDRESS>
(Received June 14, 2023, in final form July 30, 2023)

###### Abstract

This study uses a coarse-grained Monte Carlo algorithm to model and simulate
the coadsorption of a binary mixture of counterionic colloids in Gibbs
monolayers. These monolayers form at a idealized air-water interface, with one
non-soluble species confined at the interface and the second one partially
soluble in the aqueous phase. The investigation focuses on the effect of
colloidal size and charge on the thermodynamics and microstructure of the
monolayer. We find that the composition of the monolayer evolves non-trivially
with surface coverage, depending on the balance of steric and electrostatic
forces. When the electrostatic interactions are weak, the soluble species is
expelled from the monolayer upon compression, yielding a phase behaviour
particularly sensitive to the relative size of the soluble and non-soluble
colloids. By contrast, strong electrostatic interactions favour the
stabilization of the soluble particles in the monolayer and the formation of
quasi-equimolar fluids, with only a weak dependence on particle size. The
combination of these phenomena results in the formation of a number of two-
dimensional mesoscopic arrangements in the monolayer, ranging from diluted
gas-phase behaviour to domains of aggregates and percolates, and to incipient
crystalline structures.

Key words: Gibbs monolayers, Monte Carlo, equations of state, aggregation

###### Abstract

У цьому дослідженні застосовується огрублений алгоритм Монте-Карло для
моделювання коадсорбції бінарної суміші протиіонних колоїдів у моношарах
Гіббса. Моношари утворюються на ідеалізованій межі розділу ‘‘повітря-вода’’,
при цьому нерозчинна речовина обмежена межею розділу, а друга речовина
частково розчинна у водній фазі. Дослідження закцентовано на вивченні впливу
розміру колоїдної частинки та її заряду на термодинаміку та мікроструктуру
моношару. З’ясовано, що склад моношару міняється нетривіальним чином зі зміною
покриття поверхні в залежності від балансу стеричних та електростатичних сил.
Коли електростатичні взаємодії слабкі, розчинна речовина викидається з
моношару під час його стиснення, формуючи фазову поведінку, яка особливо
чутлива до відношення розмірів розчинних і нерозчинних колоїдів. Навпаки,
сильні електростатичні взаємодії сприяють стабілізації розчинних частинок у
моношарі та утворенню квазіеквімолярних плинів з незначною залежністю від
розміру частинок. Поєднання цих явищ призводить до утворення низки двовимірних
мезоскопічних структур у моношарі в діапазоні від розрідженої газової фази до
областей, що складаються з агрегатів та перколятів, і далі аж до зародків
кристалічних структур.

Ключов слова: моношари Гіббса, метод Монте-Карло, рівняння стану, агрегація

##  1 Introduction

Computer simulation has consolidated as a fundamental tool for the
rationalization and technological exploitation of the properties of molecular
liquids, soft condensed matter, and materials in general. An important step in
using computer simulation tools to investigate the properties of any system is
to obtain a suitable model for the problem at hand. There is no single recipe
for this task, and the level of complexity of the model to be used will depend
on the properties and time and length scales to be studied. Very detailed
microscopic models may provide quantitative information about some properties,
but they usually make it difficult to explore the general physical
characteristics that guide the behaviour of the system. This latter
information is typically more accessible with simplified models, which may
however compromise the acquisition of reliable quantitative and atomistic
information. The choice and design of appropriate models for a specific study
is an art in itself, which Professor Jaroslav Ilnytskyi has excelled
throughout his career, in particular within the framework of complex molecular
fluids [1, 2].

We present a computer simulation study inspired by the work of prof. Ilnytsky,
of the thermodynamics and structural properties of a multicomponent fluid of
charged particles dispersed in an interfacial monolayer. More specifically, we
investigate generic cases of Gibbs monolayers representing dispersed
amphiphilic particles at an air-water interface [3, 4]. These systems began to
be studied in modern times already in the 18th century, from the celebrated
observations of films of olive oil on water carried out by Benjamin Franklin
[5] and later revisited by Lord Rayleigh [6]. Subsequently, Agnes Pockels [6,
7] designed the famous trough to measure the surface tension of monolayers,
which was eventually improved by Irving Langmuir and his assistant Katherine
Blodgett, implementing the efficient transfer of the monolayer to a solid
substrate [8]. A monolayer is said to be of Langmuir-type if all the molecules
are completely insoluble, while it is said to be a Gibbs monolayer if some of
its constituents are partially soluble in the bulk of at least one of the
phases that share the interface. In a general situation, the species belonging
to the insoluble fraction will remain confined at the interface, while those
belonging to the soluble fraction will be distributed between the fluid phases
and the interface [9]. Gibbs monolayers find a broad range of applications,
often related to sensing devices, where interfacial changes are induced when
species in solution selectively bind to non-soluble colloidal receptors hosted
by the monolayer [10, 11]. A nice example of this type of application is the
interfacial binding of biomacromolecules, such as DNA fragments, by cationic
surfactants. For such systems, it has been shown that the compression of the
monolayer affects the concentration of the soluble biomolecules in the
monolayer, a thermodynamic effect that is not yet well understood [11, 12].
Gibbs monolayers have also a potential to be utilized in the synthesis of
thin-film topologies, facilitated by counterionic template frameworks that are
incorporated to the interface from solution. An example of this is the
fabrication of Langmuir-Blodgett layered structures based on guanidinium
surfactants employing soluble carboxylate and phosphate moieties [10].

Not surprisingly, Gibbs monolayers have attracted much attention and have
motivated important developments from the field of molecular and colloidal
simulation. In this context, different levels of approximation have been
proposed, ranging from models with a highly detailed atomistic description of
the molecules and particles involved, including water [13, 14, 15], to coarse-
grained models where entire parts of the molecules are modelled as beads,
using various effective force fields [16, 17]. In many studies, the monolayer
is examined as a three-dimensional system, simultaneously modelling and
simulating the interface and the two coexisting phases. In this way, it is
possible to obtain valuable insights into the microscopic properties of the
system, especially into the conformations that favour colloidal interactions
[13, 17, 15]. However, due to the high computational cost, these models are
usually restricted to small systems and short timescales, which often
precludes the access to thermodynamic and mesoscopic structural properties.

In an attempt to complement the information obtained with more detailed
modelling studies, we have recently proposed a simple coarse-grained model for
the study of Gibbs monolayers [18]. Within our approach, only the interface
itself is simulated as a two-dimensional system with interacting particles
represented as flat disks. Amphiphilic and soluble particles can diffuse
within the interface, but only the latter ones are capable of moving in and
out of the interface. The interchange of soluble particles between the
interface and the bulk phase is simulated with Grand Canonical Monte Carlo
moves. Another aspect of our model is that the interaction between charged
particles is modelled as a dipolar interaction, an assumption that resembles
the different screening of colloidal charge induced by the two immiscible
phases that conform the interface, e.g., the air-water interfaces [19, 20].
With these ingredients, we recently studied the phase behaviour and the
structural properties of monolayers of mutually attractive soluble and non-
soluble colloids with the same size and absolute dipolar charge (but opposite
sign) [18]. For that system, an interesting sequence of disorder fluid,
cluster fluid, and crystalline phases was found, confirming the potential
richness of the structural landscape of Gibbs monolayers. In this article we
extend the investigation by introducing asymmetry in colloidal size in the two
complementary cases where either the soluble or the non-soluble component of
the mixture is larger in size. In this way, we intend to isolate the size
effects on interfacial behaviour from those of other colloidal properties.

##  2 Model and methods

The colloidal model and the Monte Carlo simulation methodology are similar to
those employed in our previous work [18]. The Gibbs monolayer at the air-water
interface is modelled using a simplified coarse-grained approach, which
represents the system as a two-dimensional fluid, comprising a binary mixture
of disk particles with different diameters, interacting through both steric
and electrostatic interactions. One of the colloidal species (hereafter
denoted as $S$) is assumed to be soluble, while its counterionic counterpart
(denoted as $M$) is non-soluble (or amphiphilic) and it is therefore confined
at the interface. Though the bulk of the liquid phase is not explicitly
considered, it is taken into account on the assumption that it is in
thermodynamic equilibrium with the interface. Thus, it is assumed that the
chemical potential of the soluble $S$ particles, denoted as $\mu_{S}$, is the
same in the bulk and at the interface. This chemical potential is dependent on
the concentration of the particles in solution, $C_{S}$, through the
expression $\mu_{S}=\mu^{\text{ref}}_{S}+k_{\text{B}}T\ln\gamma C_{S}$, where
$\mu^{\text{ref}}_{S}$ is the chemical potential in a reference state and
$\gamma$ is the activity coefficient ($\gamma=1$ for an ideal solution). Here,
$k_{\text{B}}$ is the Boltzmann constant, and $T$ is the temperature. The
chemical potential modulates the incorporation of $S$ particles from the
solution to the monolayer and this is emulated by means of Monte Carlo
simulations in the $\mu_{s}N_{M}\Pi T$ ensemble, where $N_{M}$ denotes the
fixed number of non-soluble particles in the monolayer and $\Pi$ is the
surface pressure. The Monte Carlo algorithm samples the equilibrium
configurations with a combination of random displacements of $S$ or $M$
particles, with the insertion or removal of $S$ particles, and changes in the
surface area of the monolayers. The acceptance or rejection of these move
attempts will depend on the interaction between the particles, the chemical
potential $\mu_{S}$ and the surface pressure for volume changes, following
general acceptance rules based on the Metropolis algorithm [21].

The simulations were run on systems of $N_{M}=$ 4000 non-soluble $M$ colloids,
while the number of soluble $S$ colloids is allowed to change throughout the
simulation without constraints. Each Monte Carlo cycle involved $N_{M}$
attempts for random displacements of a particle ($S$ or $M$), for changes in
the surface area, and for the insertion or deletion of a particle $S$ in or
from the interface. Using the values optimized in our previous work, the
displacement of a random particle was attempted with frequency $w_{d}=0.765$,
the changes of surface with frequency $w_{a}=0.0005$, while the
insertion/deletion moves are both attempted with frequency $w_{i}=0.230$. A
typical simulation at a given surface pressure and chemical potential starts
from a random initial configuration of low surface density. To equilibrate the
system, about $1-2\cdot 10^{6}$ Monte Carlo cycles were applied. Additional
$5\cdot 10^{5}$ Monte Carlo cycles were then run to obtain the averages of the
observables of interest.

The choice of the interaction potential between the particles is a key aspect
of the model. The interaction between charged particles at an air-water
interface is a far from trivial problem that has attracted the interest of
many researchers [19, 22, 23, 24, 20]. A long-range $r^{-3}$ dipole-dipole
interaction has been proposed to describe the interaction between charged
interfacial particles, as a consequence of different shielding of the
Coulombic forces in the air and water sides of the interface. Here, we
consider such dipolar interaction, on top of a short-range steric repulsion
between the particles. Other types of interaction, such as capillary
interactions, have been pointed out as relevant [20], but are not considered
in our model. With all these ingredients, the interfacial colloids interact
through a pair potential $U(r)=U_{LJ}(r)+U_{D}(r)$, expressed as the sum of
short-range steric repulsions and long-range effective dipole interactions
[18].

On the one hand, the short-range interaction is represented by a truncated and
shifted Lennard-Jones potentials of the form

$U_{LJ}=\left\\{\begin{array}[]{cc}4\epsilon\left[\left(\sigma_{ij}/r\right)^{12}-\left(\sigma_{ij}/r\right)^{6}+1/4\right]&~{}~{}r\leqslant\sqrt[6]{2}\,\sigma_{ij},\\\
0&~{}~{}r>\sqrt[6]{2}\,\sigma_{ij}.\end{array}\right.$ (2.1)

Here, $\sigma_{ij}=0.5(\sigma_{i}+\sigma_{j})$ is the contact distance between
the interacting particles, with diameters $\sigma_{i}$ and $\sigma_{j}$,
respectively. $r$ is the distance between particles. Note that since the
potential is truncated at the minimum of the Lennard-Jones potential well, it
is a purely repulsive contribution to the net interaction.

On the other hand, the dipolar contribution to the interaction potential reads
as follows

$U_{D}=\left\\{\begin{array}[]{cc}\epsilon\Delta_{i}\Delta_{j}\left[(\frac{\sigma}{r})^{3}-(\frac{\sigma}{r_{s}})^{3}\right]&~{}~{}r\leqslant
r_{s},\\\ 0&~{}~{}r>r_{s},\end{array}\right.$ (2.2)

where $\sigma$ and $\epsilon$ are the unit of length and energy, respectively.
$r_{s}$ is the cutoff distance of the dipole interaction, which as in our
previous article is set as $r_{s}=30\sigma$. $\Delta_{i}$ and $\Delta_{j}$ are
the dimensionless dipolar charge of particles $i$ and $j$. The values assigned
to these dipole moments in this study are $3$, $6$ or $9$ in reduced units, as
commented below.

##  3 Results

In our previous work [18], we focused on the reference situation of $M$ and
$S$ colloidal particles with equal dipolar charge modulus and diameter. Here,
we consider situations in which the colloidal particles $M$ and $S$ are
asymmetric in size. Two complementary cases are addressed in which either the
non-soluble $M$ particles or, alternatively, the soluble $S$ particles are
larger in size. We specifically considered the case with $\sigma_{M}=\sigma$
and $\sigma_{S}=3\sigma$, as well as the reverse situation where
$\sigma_{M}=3\sigma$ and $\sigma_{S}=\sigma$. The temperature and the chemical
potential of the soluble particles were set to
$T^{*}=k_{\text{B}}T/\epsilon=1$ and $\mu^{*}=\mu/\epsilon=-8$, respectively,
while surface pressure was varied in the range
$\Pi^{*}=\Pi\sigma^{2}/\epsilon=0.001-3$. Our previous study showed that
within this range of parameters the Gibbs monolayers display a rich phase
diagram.

###  3.1 Thermodynamic study of the monolayers

The equations of state of the systems under study are characterized in figure
1, in terms of the evolution of the surface pressure
$\Pi^{*}=\Pi\sigma^{2}/\epsilon$, the molar fraction of soluble colloids
$X_{s}$ and the internal energy per particle in reduced units,
$u^{*}=U/N\epsilon$, as a function of the inverse coverage,
$\theta^{-1}=A/A_{p}$. Note that $\theta^{-1}$ represents the interfacial area
per particle; $A$ is the total area of the interface, while $A_{p}=N_{M}\cdot
a_{M}+N_{S}\cdot a_{S}$ is the effective area physically occupied by the
particles at the interface, with $a_{i}=\piup\sigma^{2}_{i}/4$ the area of a
particle of species $i$.

Figure 1: (Colour online) Monte Carlo isotherms for systems with colloids of
diameters $\sigma_{S}=3,\sigma_{M}=1$ (black lines and circles) and
$\sigma_{S}=1,\sigma_{M}=3$ (red lines and circles) and dipole moments
$\Delta=3$ (panels a, d and g), $6$ (panels b, e and h) and $9$ (panels c, f
and i). The top (a to c), middle (d to f) and bottom (g to i) panels plot the
surface pressure $\Pi^{*}=\Pi\sigma^{2}/\epsilon$, the mole fraction of $S$
colloids $X_{s}$ and the energy per particle $u^{*}=U/N\epsilon$,
respectively, as a function of the inverse of the surface coverage,
$\theta^{-1}$. The horizontal dash-dotted lines indicate the equimolarity of
the oppositely charged colloids ($X_{s}=0.5$).

The main finding of the present study is that for weak dipolar interactions
($\Delta=3$), the monolayer tends to be equimolar in the two colloidal
components in the high dilution limit ($X_{s}\approx X_{m}\approx 0.5$, as
$\theta^{-1}\rightarrow\infty$), whereas it becomes progressively enriched in
the non-soluble component $M$ upon compression. The soluble component $S$
tends to be expelled from the monolayer, so that at high coverage the
monolayer becomes essentially a monocomponent system of $M$ colloids
($X_{m}\rightarrow$1, $X_{s}\rightarrow 0$, as $\theta^{-1}\rightarrow 0$, see
figure 1d). A consequence of this behaviour is that the system displays a
marked asymmetry with respect to the size ratio of the $M$ and $S$ colloids.
Figures 1a and 1g show that when the smaller colloid is the insoluble
component ($\sigma_{M}=1$), the monolayer displays a steep increase in
pressure and internal energy with increasing coverage, at $\theta^{-1}<10$. A
similar increase occurs at much higher coverages, onset at $\theta^{-1}<2$,
when the insoluble colloid is the larger one ($\sigma_{M}=3$). In general
terms, at any given surface coverage away from the dilution limit, much higher
pressures and internal energies are associated with the $\sigma_{M}=1$ fluid
in comparison to the $\sigma_{M}=3$ fluid.

The exclusion of the $S$ particles at high coverage just described for
$\Delta=3$, can be rationalized in terms of an increase in entropy associated
with the larger area available to the intefacial particles. At a higher
dipolar charge, such entropic effects are balanced by the attractive
interaction between $S$ and $M$ particles, which favours the presence of $S$
colloids in the monolayer. Hence, the thermodynamic behavior of the monolayer
changes qualitatively as the dipolar interactions become significant. This
also brings the system to a much more symmetric dependence on colloid size
ratio, meaning that it becomes progressively less relevant whether the larger
colloids belong to the soluble component or to the insoluble one. For the
present study, the isotherms with $\sigma_{M}=3$ ($\sigma_{S}=1$) and
$\sigma_{M}=1$ ($\sigma_{S}=3$) are similar for $\Delta=6$ and become almost
overlapping for $\Delta=9$ (see figure 1). Interestingly, the strong dipole
interactions largely favour the compression of the monolayer, which proceeds
with small changes in pressure and with an internal energy $u^{*}$, that
remains negative throughout the full isotherms and becomes more negative upon
compression. Note that this is in strong contrast with the increasingly net
positive $u^{*}$ values attained for $\Delta=3$. Negative values of $u^{*}$
are associated with dominant attractive interactions, indicating that an
appreciable fraction of $S$ colloids remains at the interface at all surface
coverages.

For $\Delta=6$ at intermediate dilutions ($\theta^{-1}>10$), $u^{*}$ is weakly
dependent on the surface coverage [figure 1(h)]. In this range, entropic
effects dominate, and the $S$ particles are still sterically ejected from the
system upon compression, yielding a decrease in $X_{s}$ qualitatively similar
to that found for $\Delta=3$ [figure 1(e)]. Interestingly however, the
$S$-particles re-enter the monolayer at sufficently high surface coverages,
correlating with a steep decrease in the internal energy. In the limit of high
coverage, the molar fraction reaches the values similar to those obtained in
the high dilution regime, $X_{s}\approx 0.4$. For the strongly interacting
$\Delta=9$ monolayer, the molar fraction of $S$ colloids is larger than $0.5$
at high dilution, indicating that the incorporation of an excess of soluble
colloids to the monolayer is favored under these conditions. The composition
of the monolayer approaches then equimolarity $X_{s}\approx 0.5$ at high
packing fractions.

###  3.2 Microscopic structure of the monolayers

Figure 2: (Colour online) Representative snapshots of microscopic
configurations of the Gibbs monolayers explored in this study. Red and blue
disks represent the nonsoluble ($M$) and the soluble ($S$) colloids,
respectively. Snapshots with labels A, B and C correspond to situations where
$\sigma_{M}=1$ and $\sigma_{S}=3$, while labels D, E and F indicate the
reciprocal cases, with $\sigma_{M}=3$ and $\sigma_{S}=1$. The number in the
labels indicates the dimensionless dipolar charge in each case, $\Delta=3$,
$6$ or $9$. The composition of the monolayer in each state can be inferred
from the representation of $X_{s}$ in the bottom panel, which shows an extract
of the data included in panels (c), (e) and (f) of figure 1.

The thermodynamic properties described so far are intimately related to the
microscopic and mesoscopic structures adopted by the particles within the
monolayer. The topology and relative stability of different interfacial
structural arrangements are strongly conditioned by the incorporation of
soluble colloids $S$ from the bulk phase. Figure 2 provides an overview of
different structures observed in our work. The ejection of S colloids from the
monolayer at $\Delta=3$ can be appreciated in the sequences of configurations
from 3A to 3C ($\sigma_{M}=1$, $\sigma_{S}=3$) and from 3D to 3F
($\sigma_{M}=3$, $\sigma_{S}=1$). In the high density configurations 3C and
3F, the monolayer consists of a monocomponent fluid of $M$ colloids. A trend
towards the hexatic crystalline order characteristic of the two-dimensional
discotic particle fluid is observed at high pressures and surface coverages
(configuration 3F), although a consolidated crystalline phase is not found
within the thermodynamic range included in our study for $\Delta=3$. When the
dipolar charge is increased, the particle aggregation and cluster formation
tend to occur even in the high dilution limit. This phenomenon is particularly
pronounced for $\Delta=9$, (configurations 9A and 9D) where linear aggregates
of several particles are observed, whereas for $\Delta=6$ (panels 6A and 6D)
mainly $M-S$ dimers are found. As the concentration increases, the formation
of clusters is more evident and percolated arrangements are formed. In the
more dense monolayers of the strongly interacting colloids, the percolates
colapse to an incipient square crystalline arrangements.

Figure 3: (Colour online) Radial functions associated with the $M-M$
($g_{MM}$) and $S-S$ ($g_{SS}$) self-correlations between alike-charged
colloids (black solid and blue dotted dashed curves, respectively) and with
the $M-S$ cross-correlations ($g_{MS}$) between oppositely charged colloids
(red dashed curves) for the states 6C (top left), 6F (top right), 9C (bottom
left) and 9F (bottom right) described in figure 2. Note that in 6C and 9C
cases $\sigma_{M}=1$ and $\sigma_{S}=3$, while for 6F ad 9F $\sigma_{M}=3$ and
$\sigma_{S}=1$.

Despite the abundance of structural defects, we obtain the evidence for square
cristalline symmetry in the high coverage limit (e.g. configuration 9F in
figure 2). The analysis of the corresponding radial distribution functions
confirms these findings. Figure 3 depicts an illustrative set of radial $M-M$,
$S-S$ and $M-S$ correlation functions ($g_{MM}$, $g_{SS}$ and $g_{MS}$
respectively). Configuration 9F is representative of the formation of a
crystalline arrangement; the position of the peaks in the correlation
functions indicates that the structure corresponds to a square lattice, in
which the smaller colloids are located in the intersticial sites generated by
the larger colloids. The distribution functions for configurations 6C, 6F and
9C display a strong correlation at short and medium distances, resembling the
domains of crystalline structure, which however does not survive at long
range. These configurations deserve further investigation to ensure that they
correspond to true equilibrium states and are not trapped in a metastable
configuration. Note that our simulations were performed by compressing the
monolayer from the high dilution limit. Metastability could in particular
explain the sharp drop in $X_{S}$ observed for the $\Delta=6$ and $9$ cases at
high coverage, and the differences in the equation of state observed in this
range while comparing the monolayers qith $\sigma_{M}=3$ vs. $\sigma_{M}=1$
(figure 1).

Figure 4: (Colour online) Fraction of colloids in clusters $n_{C}$ as a
function of surface coverages $\theta^{-1}$, for $\Delta=3$ (black lines and
circles), $6$ (red lines and squares) and $9$ (blue lines and triangles). Open
symbols indicate the cases with $\sigma_{M}=1$ and $\sigma_{S}=3$ while solid
symbols are the reciprocal case where $\sigma_{M}=3$ and $\sigma_{S}=1$.

Colloidal aggregation in the strongly interacting monolayers is intriguing.
Clustering is actually a common phenomenon in the fluids of charged colloids
at air-water interfaces, as known from experiments [25, 26]. The transition
from a regime dominated by small aggregates to either mesoscopic aggregation
or long-range percolated arrangements is subtle and challenging to predict.
Here, we assess the extent of aggregation by analyzing the interconnection
between neighbouring colloids as a function of surface coverage. To this end,
we considered that two particles belong to the same cluster if the distance
between them is smaller than $\sigma_{ij}+0.5\sigma$, where $\sigma_{ij}$ is
the contact distance between the particles. Within this framework, two types
of ensemble-averaged magnitudes were employed to characterize the clustering
of particles. Firstly, the distribution, $F(i)=i\cdot C(i)/N_{T}$, represents
the fraction of particles belonging to a cluster of size $i$, where $C(i)$ is
the average number of clusters of size $i$ and $N_{T}$ is the average number
of particles at the interface. Secondly, the clustering paremeter, $n_{C}$,
represents the fraction of colloids embedded in clusters of any size,
including dimers and larger clusters ($n_{C}=\sum_{i>1}F(i)$).

Figure 4 depicts the dependence of $n_{C}$ on $\theta^{-1}$ for the fluids
under study. Consistently, $n_{C}\rightarrow 1$ at high coverages, and it
decreases steadily as the system approaches the high dilution limit. Cluster
dissociation as the area accessible to the colloids increases, mainly driven
by entropy, is clearly apparent for $\Delta=3$, with $n_{C}$ reaching the
values close to 0 at high dilution. For stronger dipolar interactions,
clustering is still noticeable in the expanded monolayers. In the case of
$\Delta=6$, even for $\theta^{-1}>$ 100, around 10% of the colloidal particles
are incorporated to a cluster. For $\Delta=9$, the particle aggregation at
high dilution still reaches remarkable values for $n_{C}$ above 80%. The
domains of clustering and percolation are characterized in greater detail in
figure 5, which shows the cluster distributions $F(i)$ for relevant
configurations of $\Delta=9$ monolayer. We find that colloids are distributed
in clusters with a broad distribution of sizes. It can be noted that even at
high dilution (configurations 9A and 9D), the fraction of isolated colloid
monomers remains below 10%, while dimers amount to almost 20%, with the
majority of the remaining 70% corresponding to clusters with sizes up to
$i\approx 20$. In the percolation domain (configurations 9B and 9E) not more
than 0.5% of the colloids stay as monomers and the cluster distribution
broadens significantly, extending to sizes $i\approx 50-100$. Notably, a
coincident aggregation behaviour is found for the monolayers with
($\sigma_{M}=1$, $\sigma_{S}=3$) and with ($\sigma_{M}=3$, $\sigma_{S}=1$)
when cluster distributions are compared at similar surface coverages (9A vs.
9D and 9A vs. 9D, in figure 5). These results confirm the idea that for
intermediate and high values of $\Delta$, the presence of aggregates is one of
the main characteristics of Gibbs monolayers. Similar results were obtained
for colloidal particles of equal size [18]. One question we leave open for
future work is the stability of the clusters, whether they have long lifetimes
or, on the contrary, they are continuously forming and fracturing. Answering
this question demands simulation techniques other than the Monte Carlo
algorithm used here, such as molecular dynamics simulations.

Figure 5: (Colour online) Fraction of colloids in clusters of a given size $i$
for some states of the monolayer with $\Delta=9$ (see labels in figure 4).
Solid lines correspond to states with $\sigma_{M}=1$, while dashed lines
correspond to states with $\sigma_{M}=3$.

##  4 Conclusions and final remarks

This work presents a computer simulation study of the structural and
thermodynamic behaviour of Gibbs monolayers in the case where the soluble and
non-soluble colloids have different diameters. As such, it is an extension of
a previous work devoted to the case of colloids with equal size [18]. This
study also corroborates the diversity of phases and structures that are
possible in Gibbs monolayers depending on colloidal parameters such as size
and dipole charge, as well as thermodynamic conditions. Most of these phases
and structures consist of disordered fluids with the particles forming
aggregates.

It is interesting to see how the two complementary situations that we have
studied, namely cases in which the non-soluble colloids $M$ are smaller
($\sigma_{M}=1$ and $\sigma_{S}=3$), or larger ($\sigma_{M}=3$ and
$\sigma_{S}=1$), display the most marked differences with each other in the
case of weak colloid-colloid interactions (represented here by a dipolar
moment $\Delta=3$). This finding can be traced back to the fact that the
soluble colloids $S$ are expelled from the monolayer with increasingly surface
coverage, eventually leading to an interfacial a monocomponent fluid of $M$
particles. Stronger colloid interactions ($\Delta$= $6$, $9$) promote the
stabilization of $S$ colloids at the interface, leading to comparable
fractions of soluble and non-soluble colloids at the monolayer. Under these
conditions, most of these phases and structures consist of disordered fluids
with the particles forming aggregates. In very dilute systems we find gas
phases of free particles and dimers. In denser systems, progressively larger
clusters are formed, eventually leading to percolated arrangements. In the
limit of high coverages, we find evidence of square crystalline structures,
but also of fluids with short-range order which, due to the presence of
defects, the correlation between particles is lost at long distances.

In short, in this work we have confirmed the potential complexity of Gibbs
monolayers as a function of the system parameters, and the capacity of the
outlined simulation methodology to explore it. We believe that the results
obtained may be useful for a better understanding of the design and tailoring
of this type of two-dimensional fluids for specific technological
applications.

## Acknowledgments

The authors acknowledge support from Consejería de Transformación Económica,
Industria, Conocimiento y Universidades de la Junta de Andalucía/FEDER
(projects P20-00816 and P20-01258), and from the Spanish Ministerio de Ciencia
e Innovación and FEDER (Projects no. PID2021-126121NB-I00,
PID2019-110430GB-C22 and TED2021-130683B-C21). We are thankful to C3UPO for
the HPC facilities provided.

## References

  * [1] Ilnytskyi J., Wilson M. R., Comput. Phys. Commun., 2001, 134, No. 1, 23–32, doi:10.1016/S0010-4655(00)00187-9.
  * [2] Hughes Z. E., Stimson L. M., Slim H., Lintuvuori J. S., Ilnytskyi J. M., Wilson M. R., Comput. Phys. Commun., 2008, 178, No. 10, 724–731, doi:10.1016/j.cpc.2008.01.047.
  * [3] Myers D., Surfaces, Interfaces, and Colloids, John Wiley & Sons, Ltd, 1999, doi:10.1002/0471234990.fmatter_indsub.
  * [4] Binks B. P., Langmuir, 2017, 33, No. 28, 6947–6963, doi:10.1021/acs.langmuir.7b00860.
  * [5] Tanford C., Ben Franklin Stilled the Waves : An Informal History of Pouring Oil on Water With Reflections on the Ups and Downs of Scientific Life in General, Duke University Press, Durham, NC, 1989.
  * [6] Rayleigh, Nature, 1891, 43, No. 1115, 437–439, doi:10.1038/043437c0.
  * [7] Poskels A., Nature, 1892, 46, 418–419, doi:10.1038/046418e0.
  * [8] Langmuir I., Trans. Faraday Soc., 1920, 15, 62–74, doi:10.1039/TF9201500062.
  * [9] Melzer V., Vollhardt D., Brezesinski G., Möhwald H., J. Phys. Chem. B, 1998, 102, No. 3, 591–597, doi:10.1021/jp972860g.
  * [10] Ariga K., Yamauchi Y., Mori T., Hill J. P., Adv. Mater., 2013, 25, No. 45, 6477–6512, doi:10.1002/adma.201302283.
  * [11] Chen Q., Kang X., Li R., Du X., Shang Y., Liu H., Hu Y., Langmuir, 2012, 28, No. 7, 3429–3438, doi:10.1021/la204089u.
  * [12] Du X., Wang Y., Ding Y., Guo R., Langmuir, 2007, 23, No. 15, 8142–8149, doi:10.1021/la700955f.
  * [13] Giner Casares J. J., Camacho L., Martín-Romero M. T., López Cascales J. J., ChemPhysChem, 2010, 11, No. 10, 2241–2247, doi:10.1002/cphc.201000048.
  * [14] Wang D., de Jong D. H., Rühling A., Lesch V., Shimizu K., Wulff S., Heuer A., Glorius F., Galla H.-J., Langmuir, 2016, 32, No. 48, 12579–12592, doi:10.1021/acs.langmuir.6b02496.
  * [15] Bai Y., Wen W., Gao Y., Cui W., Sun Y., Yan P., J. Mol. Liq., 2022, 368, 120804, doi:10.1016/j.molliq.2022.120804.
  * [16] Düchs D., Schmid F., J. Phys.: Condens. Matter, 2001, 13, No. 21, 4853, doi:10.1088/0953-8984/13/21/313.
  * [17] González-Castro C. A., Ramírez-Santiago G., Phys. Rev. E, 2015, 91, 032409, doi:10.1103/PhysRevE.91.032409.
  * [18] Cuetos A., Morillo N., Martínez-Haya B., Langmuir, 2020, 36, No. 11, 2877–2885, doi:10.1021/acs.langmuir.9b03886.
  * [19] Hurd A. J., J. Phys. A: Math. Gen., 1985, 18, No. 16, L1055, doi:10.1088/0305-4470/18/16/011.
  * [20] Bleibel J., Domínguez A., Oettel M., Eur. Phys. J. Spec. Top., 2013, 222, 3071–3087, doi:10.1140/epjst/e2013-02076-9.
  * [21] Frenkel D., Smit B., Understanding Molecular Simulation, Academic Press, San Diego, 2002, doi:10.1016/B978-012267351-1/50008-0.
  * [22] Majee A., Schmetzer T., Bier M., Phys. Rev. E, 2018, 97, 042611, doi:10.1103/PhysRevE.97.042611.
  * [23] Chen W., Tan S., Zhou Y., Ng T.-K., Ford W. T., Tong P., Phys. Rev. E, 2009, 79, 041403, doi:10.1103/PhysRevE.79.041403.
  * [24] Lian Z., J. Chem. Phys., 2016, 145, No. 1, doi:10.1063/1.4954913.
  * [25] Suzuki D., Horigome K., J. Phys. Chem. B, 2013, 117, No. 30, 9073–9082, doi:10.1021/jp4035166.
  * [26] Nallamilli T., Ragothaman S., Basavaraj M. G., J. Colloid Interface Sci., 2017, 486, 325–336, doi:10.1016/j.jcis.2016.10.009.

Ukrainian 0

Вплив заряду та розміру на коадсорбцію протиіонних колоїдів у моношарах Гіббса
[Ж. М. Гомес-Верду, Б. Мартінес-Айя, А. Куетос]Ж. М. Гомес-Верду, Б. Мартінес-
Айя, А. Куетос

Факультет фізичних, хімічних та природніх систем, Університет Пабло де
Олавіде, 41013 Севілья, Іспанія