# Multiscale Modeling of Sorption Kinetics

Clarissa Astuto  King Abdullah University of Science and Technology (KAUST),
4700, Thuwal, Saudi Arabia  Antonio Raudino  Department of Chemical Sciences,
University of Catania, Italy  Giovanni Russo  Department of Mathematics and
Computer Science, University of Catania, Italy

###### Abstract

In this paper we propose and validate a multiscale model for the description
of particle diffusion in presence of trapping boundaries. We start from a
drift-diffusion equation in which the drift term describes the effect of
bubble traps, and is modeled by a short range potential with an attractive
term and a repulsive core. The interaction of the particles attracted by the
bubble surface is simulated by the Lennard-Jones potential that simplifies the
capture due to the hydrophobic properties of the ions. In our model the effect
of the potential is replaced by a suitable boundary condition derived by mass
conservation and asymptotic analysis. The potential is assumed to have a range
of small size $\varepsilon$. An asymptotic expansion in the $\varepsilon$ is
considered, and the boundary conditions are obtained by retaining the lowest
order terms in the expansion.

Another aspect we investigate is saturation effect coming from high
concentrations in the proximity of the bubble surface. The validity of the
model is carefully checked with several tests in 1D, 2D and different
geometries.

##  1 Introduction

The diffusion of particles in presence of (moving or static) traps is an
interesting topic because of its applications to different fields such as
chemistry, physics and biology. An interesting case is to consider the
trapping of diffusing particles near the surface of the traps. The problem of
surfactant diffusion, that we consider in this paper, has been investigated by
several authors [1, 2, 3, 4, 5] and in each work the transport phenomena in
ionic solutions requires the calculation of the space and time evolution of
the concentration of negative (anions) and positive (cations) species
diffusing in a confined space at a time $t$, as we did in [6].

In this work we study the behaviour of one species at a time because here we
focus on the space multiscale aspects of the problem, while many authors have
focused on the time multiscale challenge of the phenomenon [4, 7]. The most
frequently used strategy to describe the dynamic properties of liquid
adsorption layers is the dynamic surface or the interfacial tension. Many
different techniques exist to cover the time interval from milliseconds to
hours and days.

Introducing the local concentration of ions $c=c(\vec{x},t)$, its time
evolution in a fluid is governed by the conservation law

$\displaystyle\frac{\partial c\left(\vec{x},t\right)}{\partial t}=-\nabla\cdot
J(\vec{x},t).$ (1)

In absence of bubbles, the flux $J$ represents the standard diffusion given by
Fick’s law $J=-D\nabla c$, where $D$ is the diffusion coefficient. In the
presence of bubbles (traps), the flux expression is augmented to include a
drift effect

$\displaystyle\displaystyle J=\ -D\left(\nabla c+\ \frac{1}{k_{B}T}c\nabla
V\right)$ (2)

where $k_{B}$ is the Boltzmann’s constant, $T$ is the absolute temperature
(assumed to be constant) and $V(x)$ is a suitable potential function that
models the attractive-repulsive behavior of the bubble with the particles. In
particular, when a particle is near the bubble, it is attracted towards its
surface, and when it is at a very short distance from the surface of the
bubble, the potential is designed to repulse the particle in order to simulate
impermeability, therefore trapping a portion of particles in the proximity of
the surface of the bubble.

The potential $V$, in presence of positive and negative diffusing particles,
may depend also on $c$ to describe the electrostatic interaction between the
two species for high concentrations. In [6] we partitioned the potential as:

$\displaystyle V=\ V_{\rm{ion-bubble}}+\ V_{\rm{ion-ion}}$ (3)

where the term $V_{\rm{ion-bubble}}$ describes the interaction between the
bubble with the ions, located at a generic position $r$ (Fig. 3), and the
(possibly oscillating) interface, while the term $V_{\rm{ion-ion}}$ accounts
for the interaction among the diffusing ions. In this work we focus on
modeling only the interaction of anions with the bubble, so we do not take
into account the term $V_{\rm{ion-ion}}$, describing the interaction with the
cations. Such term can indeed be neglected if only anions are present in the
fluid, or, more realistically, when the Debye length associated to the two
types of ions is much smaller than the size of the bubble (quasi-neutral
limit, see [6]). The more realistic case of bipolar transport, which includes
the Coulomb interaction between ions, will be subject of future investigation.

In this work we also include the study of the _saturation_ effect of the
system when high concentrations get stuck to the surface of the trap ([8, 9]).
Standard expressions for the flux satisfies a Fermi-Dirac form for the mixing
entropy: [10, 11, 12] $S_{\rm mix}=-k(c\log(c)+(1-c)\log(1-c))$ providing the
following expression for the flux

$\displaystyle\displaystyle J=\ -D\left(\nabla c+\
\frac{1}{k_{B}T}c(1-c)\nabla V\right)$ (4)

and it leads to a non-linear equation in $c$. Note that for low concentration
this expression reduces to Eq. (2). The main reason why studies take into
account saturation is that it seems to be relevant near the trap, where local
concentrations may become non negligible. Several authors have demonstrated
that it is more appropriate to use models with nonlinear saturating mobility
instead of the classical linear mobility models, both for a single species
([13, 14, 15, 16]) and also for multiple species ([13, 17]). In 4 we see the
effect of the flux defined in Eq. (4) that is different from the one acting in
dilute solutions where one can set $1-c\approx 1$.

In the linear model there is no upper bound on the concentration $c$, while in
the model which makes use of Eq. (4) the natural bound $c<1$ is always
satisfied.

Numerical results and the resulting saturation curve show how the
concentration can not go over certain values ($c\leq 1$ in our case). This
effect is known to exist when the ions are exposed to a sufficiently strong
potential (trap) but it can not be obtained in a limit of _quasi-neutrality_
[6].

Let us briefly summarize similar models dealing with the same issue that have
been widely discussed in the literature, because of numerous applications in
different problems in chemistry, physics and biology. For the sake of
simplicity, people consider a semi-infinite mono-dimensional system, i.e.:
$0<x<\infty$.

Since the potential acting on the diffusing particles is mainly located in a
narrow region near the left boundary and null otherwise (see Fig. 1), the
transport equation can be modeled by a standard diffusion equation (Eq. (5))
together with the condition that the flux of particles
$\displaystyle-J=D\left.\frac{\partial c}{\partial x}\right|_{x=0}$, entering
the left boundary, equals the time variation of surface concentration
($-J=\partial\Gamma/\partial t$) at the interface (continuity equation). So
doing, the effect of the potential is adsorbed into a boundary condition at
$x=0$:

$\displaystyle\displaystyle\frac{\partial c}{\partial t}$ $\displaystyle=$
$\displaystyle D\frac{\partial^{2}c}{\partial x^{2}}\quad{\rm for}\quad
x\in[0,+\infty[,t>0$ (5)
$\displaystyle\displaystyle\frac{\partial\Gamma}{\partial t}$ $\displaystyle=$
$\displaystyle D\left.\frac{\partial c}{\partial x}\right|_{x=0}\quad{\rm
for}\quad x=0,t>0$ (6)

and initial condition

$c(x,0)=c_{0}\quad x>0,t=0$ (7)

Eq.(5-6) must be supplemented by the further conditions:

$\displaystyle c|_{x=0}$ $\displaystyle=$ $\displaystyle f(\Gamma)$ (8)
$\displaystyle c|_{x=\infty}$ $\displaystyle=$ $\displaystyle c_{0}$ (9)

In Eq. (6), $f(\Gamma)$ is a prescribed function (see [2, 3, 4, 5, 7])
(derived by empirical or theoretical arguments) relating $c|_{x=0}$ (the
diffusant concentration near the surface satisfying the diffusion equation
(5)) to the still unknown diffusant concentration $\Gamma$ at the interface
between the two adjacent phases, Eq. (8). $\Gamma$ is an important quantity
because it is related to experimentally accessible parameters like the surface
tension. Other models in the literature (see, e.g.,[18]), are based on a
similar assumption because in addition to the bulk properties, the interface
characteristics can affect drastically the motion of substances in the liquid.
The presence of any little surfactant contamination can significantly affect
the bubble motion and shape. Sorption effect is reversible: in some part of
the surface particles may be absorbed, while in other part they may be
desorbed, according to the sign of $\partial c/\partial x$. These adsorption
and desorption fluxes at the bubble surface are governed by diffusion and lead
to the formation of a so-called diffusion boundary layer adjacent to the
surface [18].

In the simplest approximation (investigated by Sutherland a long time ago
[19]) it is assumed: $f(\Gamma)\approx K\Gamma$. Such an approximation leads
to rather cumbersome but exact analytical results (the same reference of S.S.
Dukhin as before) for the time evolution of $\Gamma$.

Our model represents an advancement over the classical approaches because:

A) Our approach allows the derivation of Sutherland’s constant $K$ from first
principles in the case of dilute solutions.

B) It can be generalized to concentrated solutions, providing a nonlinear
relation between the surface concentration and the concentration near the
trap.

C) A detailed validation of the derived model with the original drift-
diffusion model described by a potential is performed.

D) The multiscale model provides tremendous savings in computational
complexity over the model based on fully resolved potential, even if the
latter is solved adopting Adaptive Mesh Refinement, as illustrated in detail
in Appendix A.

The plan of the paper is the following: in Sec. 2 we derive the multiscale
model in 1D and 2D and then we extend the condition to 3D. In Sec. 3 we
validate the model comparing its numerical results with those obtained with
the original drift-diffusion model defined by Eqs. (1-2). Sec. 4 is devoted to
the treatment of the saturation effect, which becomes relevant for non-
negligible concentrations near the bubble surface. In the final Sec. 5 we draw
some conclusions.

##  2 Modeling multiscale problem

Modeling surfactant diffusion and trapping is very challenging because of the
multiple scales involved in space. The range of the attractive-repulsive core
of a trap is of the order of nanometers, a length that is prohibitively small
for detailed simulations.

To overcome such difficulties we propose a _multi-scale model_. The ions are
modeled by a volumetric concentration $c(x,t)$ which satisfies a simple
diffusion equation in the domain of the liquid, coupled with a superficial
one, which describes the carriers trapped on the surface. The two
concentrations are related by a mass balance equation which acts as a suitable
boundary condition for the concentration in the bulk.

The multiscale model has been tested in various geometries and dimensions $d$.
For $d=1$ the domain is $[-\varepsilon,1]$, the action of the trap is defined
with an external potential that is non negligible in a few $\varepsilon$ and
its minimum is at the origin. For $d=2$, we first consider a natural extension
of the 1D case with a trap described by a potential that is constant in
$y$-direction. Secondly, we consider a 2D region inside a square, external to
a circular disk. The arbitrary boundary of the bubble is described by a level
set function. The set of equations is then solved by a finite difference
scheme on a regular Cartesian mesh, which makes use of regular grid points
(inside the domain $\Omega$) and ghost points (neighbors of regular grid
points which are out of $\Omega$). The equations on the regular points are
obtained by finite difference discretization of the Laplacian operator, while
the equations on each ghost point are given by interpolating the numerical
solution near the ghost point and imposing the boundary condition on the
interpolant at the orthogonal projection of the point on $\partial\Omega$ (see
[20, 21] and Sec. 3 for more details).

Upon the validation in the geometries described below, it is observed that
when the range of the interaction decreases, the agreement of the multiscale
model and the detailed numerical solution of the full drift diffusion model
becomes better and better.

###  2.1 1D Multiscale model

Figure 1: Representation of the external potentials $V(x)$, on the left, and
$U(\xi)$, on the right, for $\varepsilon=0.05$. $U(\xi)$ is obtained from a
change of variable of $V(x)$, with $\xi\in[0,L+1]$. After $\xi=L+1$, the
contribution of the potential $U(\xi)$ is negligible. For this reason, on the
left the dashed line $x=\varepsilon L$ denotes the right boundary of
$\Omega_{b}^{\varepsilon}$, the domain affected by the potential. Numerically,
we notice an appropriate value for this quantity is $L=2$.

In this subsection we derive the multiscale model in 1D space dimension. The
drift-diffusion Eqs. (1)-(2) in 1D read

$\displaystyle\displaystyle\frac{\partial c}{\partial t}+\frac{\partial
J}{\partial x}=0,$ (10)

where

$\displaystyle J=-D\left(\frac{\partial c}{\partial
x}+\frac{1}{k_{B}T}cV^{\prime}\right).$ (11)

We assume that the full domain is $\Omega=[-\varepsilon,1]$ and then we
partition it into a bubble domain
$\Omega^{\varepsilon}_{b}=[-\varepsilon,L\varepsilon]$, inside which the ions
feel the effect of the potential, and the fluid domain, which is not affected
by the bubble, $\Omega^{\varepsilon}_{f}=[L\varepsilon,1],$ see Fig. 1.
Therefore, the potential has effect only in $\Omega^{\varepsilon}_{b}$, while
it is constant for $x\in\Omega^{\varepsilon}_{f}$, say $V(x)=0$.

As a prototype potential we consider the Lennard-Jones potential (LJ), which
describes attraction at long distances and repulsion at short distances due to
Van der Waals and Pauli terms, respectively. A typical shape of the potential
$V(x)$ is shown in Fig. 1 and it takes the form

$\displaystyle
V(x)=E\left(\left(\frac{x+\varepsilon}{\varepsilon}\right)^{-12}-2\left(\frac{x+\varepsilon}{\varepsilon}\right)^{-6}\right).$

where $\varepsilon$ denotes the range of the potential and $E$ represents the
depth of the well (see Fig. 1, left panel). It is convenient to adopt a non
dimensional form of the potential, expressed as a function of the rescaled
variable $\xi=1+x/\varepsilon\in[0,L+1]$:

$U(\xi)=\phi\left(\xi^{-12}-2\xi^{-6}\right)$ (12)

where $\displaystyle\phi={E}/{k_{B}T}$ represents the ratio between the depth,
$E$, of the potential well of the bubble and $k_{B}T$. Typical values of
$\phi$ should range between 10 and 20 [22]. $L$ is a scaled distance beyond
which the potential is negligible.

Let us assume that the trap is located at $x=0$ and its ”thickness” is of
order of $\varepsilon$. The original model requires to solve Eq. (10) in
$\Omega^{\varepsilon}=\Omega^{\varepsilon}_{b}\cup\Omega^{\varepsilon}_{f}$
with initial condition $c=c_{0}(x)$ when $t=0$ and with boundary conditions

$J=0\quad{\rm at}\>x=-\varepsilon,{\rm and}\>x=1$ (13)

to simulate wall effects. If the potential diverges at $x=-\varepsilon$ such a
condition is not necessary, since the potential wall prevents the particles
from approaching the boundary. This condition can be used in case we simplify
the description of the potential by an attractive well and a vertical
boundary.

The mathematical model proposed so far presents a multiscale nature, since the
influence of the bubble on the particles is simulated in
$\Omega^{\varepsilon}_{b}$ that is very small compared to the rest of the
domain, leading to some extra computational effort to ensure that the correct
behaviour is captured. In order to overcome this difficulty, we propose a
multiscale model for $\varepsilon\ll 1$ in which the domain of the problem is
reduced to $\Omega^{\varepsilon}_{f}$ and the effect of $V(x)$ on
$\Omega^{\varepsilon}_{b}$ is approximated by a suitable boundary condition at
the origin obtained as follows. Let us assume that the potential
$V_{\varepsilon}(x)$ depends on $\varepsilon$ but maintains the same
functional form, i.e. there exist two functions:
${U}\colon[0,L+1]\rightarrow\mathbb{R}$ (e.g. the LJ function defined in Eq.
(12)) and $\mathcal{U}\colon\mathbb{R}\rightarrow\mathbb{R}$ that do not
depend on $\varepsilon$ and such that

$\displaystyle\displaystyle\mathcal{U}(\xi)=\left\\{\begin{array}[]{l}+\infty\quad\xi\leq
0\\\ U(\xi)\quad\xi\in[0,L+1]\\\ 0\quad\quad\xi>L+1\end{array}\right.$ (17)

It is clear that the solution $c$ and the flux $J$ depend on $\varepsilon$ as
well, and then the original problem can be written as

$\displaystyle\displaystyle\frac{\partial c_{\varepsilon}}{\partial
t}+\frac{\partial J_{\varepsilon}}{\partial x}$ $\displaystyle=$
$\displaystyle 0\quad\rm{in}\,\Omega^{\varepsilon}$ (18)
$\displaystyle\displaystyle J_{\varepsilon}$ $\displaystyle=$
$\displaystyle-D\left(\frac{\partial c_{\varepsilon}}{\partial
x}+\frac{1}{k_{B}T}\,c_{\varepsilon}\,V^{\prime}_{\varepsilon}\right)$ (19)

For $x\in\Omega^{\varepsilon}_{b}$ we can use the scaled variable
$\displaystyle\xi=1+\frac{x}{\varepsilon}\in[0,L+1]$ and then the flux in
$\Omega^{\varepsilon}_{b}$ can be written as

$J_{\varepsilon}=-D\frac{1}{\varepsilon}\left(\frac{\partial
c_{\varepsilon}}{\partial\xi}+c_{\varepsilon}\,U^{\prime}\right).$ (20)

Rewriting (18) in the rescaled variable and using Eq. (20) we obtain

$\frac{\partial c_{\varepsilon}}{\partial
t}=-\frac{1}{\varepsilon}\frac{\partial
J_{\varepsilon}}{\partial\xi}=\frac{D}{\varepsilon^{2}}\frac{\partial}{\partial\xi}\left(\frac{\partial
c_{\varepsilon}}{\partial\xi}+c_{\varepsilon}U^{\prime}(\xi)\right).$ (21)

The range of the scaled variable $\xi$ does not depend on $\varepsilon$ and
then we can assume that $c_{\varepsilon}(\xi,t)$ has the following expansion
in $\Omega^{\varepsilon}_{b}$:

${c_{\varepsilon}(\xi,t)=c^{(0)}(\xi,t)+\varepsilon^{2}c^{(1)}(\xi,t)+\varepsilon^{4}c^{(2)}(\xi,t)+\cdots.}$
(22)

Inserting the expansion (22) into Eq. (20) we obtain the following expansion
for the flux:

$\varepsilon
J_{\varepsilon}=J^{(0)}+\varepsilon^{2}J^{(1)}+\varepsilon^{4}J^{(2)}+\cdots,$
(23)

with

$J^{(k)}=-D\left(\frac{\partial
c^{(k)}}{\partial\xi}+c^{(k)}U^{\prime}(\xi)\right),\quad k\geq 0.$

Using expansion (23) in (21) we obtain, to the various orders in
$\varepsilon$:

$\displaystyle O(\varepsilon^{-1}):$
$\displaystyle\phantom{\quad\frac{\partial c^{(k)}}{\partial
t}+}\frac{\partial J^{(0)}}{\partial\xi}$ $\displaystyle=0,$ (24)
$\displaystyle O(\varepsilon^{k}):$ $\displaystyle\quad\frac{\partial
c^{(k)}}{\partial t}+\frac{\partial J^{(k)}}{\partial\xi}$
$\displaystyle=0,\quad k\geq 0.$ (25)

Eq. (24) states that the lowest order flux $J^{(0)}$ is constant.
Condition $\lim_{x\to-\varepsilon}J(x,t)=0$, which is a consequence of the
singularity of the potential, gives:

$\frac{\partial c^{(0)}}{\partial\xi}+U^{\prime}(\xi)c^{(0)}=0.$ (26)

Therefore, it is possible to integrate the equation

$\frac{1}{c^{(0)}}\frac{\partial c^{(0)}}{\partial\xi}=-U^{\prime}(\xi)$ (27)

whose solution is

$c^{(0)}(\xi,t)=c^{(0)}(L+1,t)\exp\left(-U(\xi)+U(L+1)\right)=c^{(0)}(L+1,t)\exp\left(-U(\xi)\right)$
(28)

since $U(L+1)=0$.

Now we assume that $V_{x}(\varepsilon L)=0$, which is compatible with having a
smooth potential with compact support. Notice that, strictly speaking, LJ
potential is not with compact support, however it is negligible small,
together with its derivative, at distances which are sufficiently larger than
its range $\varepsilon$ (see Table 1 and Fig. 2.)

Integrating (10) in the range of the potential, using the approximation
$c(x,t)=c^{(0)}(\xi,t),\,x\in\varepsilon[-1,L]$, the zero boundary condition
(13) at $x=-\varepsilon$ and that $V_{\varepsilon}(\varepsilon
L)=V^{\prime}_{\varepsilon}(\varepsilon L)=0$, we obtain

$\frac{d}{dt}\int_{-\varepsilon}^{\varepsilon
L}c(x,t)\,dx+J_{\varepsilon}(\varepsilon L)=0,$

for which it follows

$\displaystyle\frac{\partial c(\varepsilon L,t)}{\partial
t}\;\varepsilon\int_{0}^{L+1}\exp\left(-U(\xi)\right)d\xi-D\frac{\partial
c(\varepsilon L,t)}{\partial x}$ $\displaystyle=$ $\displaystyle 0$

that represents a boundary condition at $x=\varepsilon L$. Notice that, out of
the range of the potential, the flux $J$ is given by just the diffusion term,
thus Eq. (25) will not be used since for sufficiently small values of
$\varepsilon$ it is $c(x,t)\approx c^{(0)}(x,t)$, and Eq. (24) is used to
provide an effective boundary condition for the diffusion equation.

Summarizing, the multiscale model can be stated as a limit model for
$\varepsilon\to 0$ as follows:

$\displaystyle\frac{\partial c}{\partial t}$ $\displaystyle=$ $\displaystyle
D\frac{\partial^{2}c}{\partial x^{2}}\quad{\rm{in}}\,x\in[0,1]$ (29)
$\displaystyle\frac{\partial c}{\partial x}$ $\displaystyle=$ $\displaystyle
0\quad{\rm{at}}\,x=1$ (30) $\displaystyle M\frac{\partial c}{\partial t}$
$\displaystyle=$ $\displaystyle D\frac{\partial c}{\partial
x}\quad{\rm{at}}\,x=0$ (31)

where

$M=\varepsilon\int_{0}^{L+1}\exp\left(-U(\xi)\right)d\xi.$ (32)

As $\varepsilon\to 0$ one could let $L\to\infty$ in such a way that
$L\varepsilon\to 0$, so in the expression for $M$ the integral can be written
from $0$ to $\infty$. In practice we use a finite value for $L$. Except where
otherwise stated, we used $L=2$. In Fig. 2 we see that this is a convenient
value for the parameter $L$, since the quantity $\mathcal{I}_{L}(\phi)$
approximately does not change increasing $L$. Also in Table 1 we see an
agreement with the choice of the parameter $L$. We show that the action of the
potential $U(\xi)$ is negligible in the domain $[L+1,\infty]$, with $L=2$.

We observe that, if the potential does not depend on $\varepsilon$,
$M\rightarrow 0$ as $\varepsilon\rightarrow 0$ and then the condition (31)
reduces to a zero Neumann boundary condition. Therefore, in validating the
model, we shall assume that $\varepsilon\to 0$ with $M$ fixed and we explain
how to perform such a limit.

There are several ways to choose the values of the parameters $M$,
$\varepsilon$, and of the function $U(\xi)$ (that describes the shape of the
bubble interaction with the ions) in Eq. (32). As we said before, since $M\to
0$ such as $\varepsilon\to 0$, we may consider a potential
$U=U^{\varepsilon}(\xi)$ which scales in such a way that $M$ is finite as
$\varepsilon$ approaches $\tilde{\varepsilon}$, where $\tilde{\varepsilon}$ is
a typical length of the attractive-repulsive core of the bubble and it is
prohibitively small for a detailed calculation using Eq. (2). In practice, we
can assume that $M$ is non negligible when $\varepsilon\ll l$, where $l$ is a
typical length scale of the physical setup. Here we show how to choose the
parameters such that $M$ is constant when $\varepsilon\to\tilde{\varepsilon}$.
The role of $\tilde{\varepsilon}$ will be evident in Sec. 4 where we describe
the effect of saturation.

Once we define $U(\xi)$ in Eq. (12), the resulting expression for $M$ from Eq.
(32) is

$\displaystyle M$ $\displaystyle=$
$\displaystyle\varepsilon\mathcal{I}_{L}(\phi)$ (33)

where

$\mathcal{I}_{L}(\phi)=\int_{0}^{L+1}\exp\left(-\phi\left(\xi^{-12}-2\xi^{-6}\right)\right)d\xi.$
(34)

and it is represented in Fig. 2. Notice that, since the integrand in (34) is
exponentially large in $\phi$ over a finite interval, then we can choose
$(\varepsilon,\phi)$ such that $M$ is constant by solving the non linear
equation for $\phi$, $\mathcal{I}(\phi)=M/\varepsilon$.

$L$ | 2 | 4 | 6 | 8
---|---|---|---|---
ratio | $8\times 10^{-9}$ | $10^{-9}$ | $2\times 10^{-10}$ | $8\times 10^{-11}$
Table 1: In this table we define ratio :=
${\int_{L+1}^{\infty}U(\xi)d\xi}/{\int_{0}^{\infty}U(\xi)d\xi}$, for different
values of $L$, and $\phi=6$. Here we show that we obtain a good approximation
of the integral of the potential $U(\xi)$ with a few units of $L$. Figure 2:
Plot of the profile of $\mathcal{I}(\phi)$, defined in Eq. (34), with M = 3
and for different values of $L\in\\{2,4,6,8,10\\}$. This quantity increases
with $\phi$ and we see that, for sufficiently large values of $\phi$, say
$\phi>6$, the values of $\mathcal{I}(\phi)$ are almost independent on $L$.

###  2.2  Extension to more dimensions

Figure 3: Scheme of the boundary condition in 2D where $\hat{n}$ is the
outgoing normal vector to the domain $\Omega$, represented by the white region
outside the light gray bubble, and $\varepsilon$ is the thickness parameter of
the bubble that may be a function of space. The contour lines represent the
iso-potential lines. The white region near the bubble surface represents the
local minima of the potential well, and represents the boundary $\Gamma_{B}$
of the bubble. The potential changes in the direction $\hat{n}$, while it is
constant along the tangential direction $\hat{\tau}$.

In this section we provide a unified treatment of 2D and 3D case, making use
of the projection operator on the line (2D case) or the surface (3D case),
which delimits the bubble ${\mathcal{B}}$.

We assume that the bubble $\mathcal{B}$ is implicitly defined by a level set
function $\Phi(\vec{x})$, $\vec{x}\in\mathbb{R}^{d}$, $d=2,3$, and assume
$\Phi(\vec{x})<0$ inside the bubble and $\Phi(\vec{x})>0$ in $\Omega$. For
example $\Phi$ may be the signed distance from $\Sigma_{B}$ (see [23]), where
$\Sigma_{B}$ denotes the boundary of the bubble $\mathcal{B}$. The unit normal
$\hat{n}$ to $\Sigma_{B}$ is given by
$\hat{n}(\vec{x})=\nabla\Phi(\vec{x})/|\nabla\Phi(\vec{x})|$, $\forall\vec{x}$
on $\Sigma_{B}$. Note that the $\hat{n}$ is naturally defined everywhere in
$\Omega\cup\mathcal{B}$. We denote by $n_{i},i=1,\ldots,d$ the Cartesian
components of $\hat{n}$, and we shall denote by
$h_{ij}=\delta_{ij}-n_{i}n_{j}$ the projection operator on the tangent plane
on $\Sigma_{B}$. Finally, we denote by $\tilde{\partial}_{i}\equiv
h_{ij}\partial_{j}$ the gradient operator on the tangent plane to
$\Sigma_{B}$. We adopt standard Einstein’s convention of summation over
repeated indices. Furthermore, we assume that the surface $\Sigma_{B}$ is
smooth, and its curvature $\kappa$ times the range of the potential
$\tilde{\varepsilon}$ is much smaller than one.

For points near the boundary, we decompose the flux vector
$J_{i},i=1,\ldots,d$ into a normal and tangential component, i.e.

$\displaystyle
J_{i}=\underbrace{n_{i}J_{n}}_{J^{n}_{i}}+\underbrace{h_{ij}J_{j}}_{J^{\tau}_{i}}$
(35)

where $J_{n}=n_{i}J_{i}$ and $J^{\tau}_{i}=h_{ij}J_{j}$.

The divergence of the flux is given by:

$\nabla\cdot\vec{J}=\partial_{i}J_{i}$ (36)

If we substitute Eq.(35) in Eq.(36) we obtain

$\nabla\cdot\vec{J}=\partial_{i}\left(n_{i}J_{n}\right)+\partial_{i}\left(h_{ij}J_{j}\right).$
(37)

Now the first term in the last equation can be written as

$\displaystyle\partial_{i}\left(n_{i}J_{n}\right)$
$\displaystyle=n_{i}\partial_{i}J_{n}+J_{n}\underbrace{\partial_{i}n_{i}}_{\chi_{ii}}$
(38)

where $\chi_{ij}=\partial_{i}n_{j}$ denotes the second fundamental form of the
surface $\Sigma_{B}$, and its trace $\chi_{ii}$ gives the curvature of the
surface $\Sigma_{B}$:

$\displaystyle\chi_{ii}=\kappa=\left\\{\begin{array}[]{rl}\displaystyle
1/R&{\rm in}\,2D\\\ \kappa_{1}+\kappa_{2}&{\rm in}\,3D,\end{array}\right.$
(41)

where $R$ is the local radius of curvature, $\kappa_{1}$ and $\kappa_{2}$
denote the Gauss principal curvatures.

Since the gradient of the potential is parallel to the normal, then one has

$J^{\tau}_{i}=-D\tilde{\partial}_{i}c=-Dh_{ij}\partial_{j}c$ (42)

At this point we substitute the quantities (38-42) in Eq.(36) and the
expression for the divergence of the flux becomes:

$\nabla\cdot\vec{J}=n_{i}\partial_{i}J_{n}+J_{n}\chi_{ii}-D\Delta_{\perp}c$
(43)

where
$\Delta_{\perp}c=\tilde{\partial}_{i}\tilde{\partial}_{i}c=h_{ij}\partial_{i}\partial_{j}c-\chi_{ii}n_{j}\partial_{j}c$
denotes the surface Laplacian of the concentration (see Appendix B).

The evolution equation of for the concentration near the bubble surface is
then given by

$\frac{\partial c}{\partial t}=-\frac{\partial J_{n}}{\partial
n}-J_{n}\chi_{ii}+D\Delta_{\perp}c.$ (44)

Following the procedure adopted in 1D, let $\vec{x}$ be a point in
$\Sigma_{B}$ and integrate along the normal direction:

$\displaystyle\frac{\partial}{\partial t}\int_{-\varepsilon}^{\varepsilon
L}c(\vec{x}+r\hat{n})\,dr$
$\displaystyle=-J_{n}(\vec{x}+L\varepsilon\hat{n})-\int_{-\varepsilon}^{\varepsilon
L}\chi_{ii}J_{n}(\vec{x}+r\hat{n})\,dr$
$\displaystyle+D\Delta_{\perp}\int_{-\varepsilon}^{\varepsilon
L}c(\vec{x}+r\hat{n})\,dr$

where we dropped the term $J_{n}(\vec{x}-\varepsilon\hat{n})$ which vanishes
because of the repulsive core.

Notice that the term $\int_{-\varepsilon}^{\varepsilon
L}\chi_{ii}J_{n}(\vec{x}+r\hat{n})\,dr=(1+L)\varepsilon\langle\chi_{ii}J_{n}\rangle$
can be neglected since we assume that the range of the potential is much
smaller than the radius of curvature (2D) or to the inverse of the mean
curvature (3D).

Now we consider a change of variable from $\vec{x}$ to $\xi$ in the
integration intervals, and obtain

$\displaystyle\varepsilon\frac{\partial}{\partial
t}\int_{0}^{L+1}c(\vec{x}+\varepsilon\xi\hat{n})\,dr$
$\displaystyle=-J_{n}(\vec{x}+L\varepsilon\hat{n})+\varepsilon
D\Delta_{\perp}\int_{0}^{L+1}c(\vec{x}+\varepsilon\xi\hat{n})\,d\xi$

Let us consider the normal components of the flux $J_{n}$:

$J_{n}=-D\left(\frac{\partial c}{\partial{n}}+\frac{1}{k_{B}T}\frac{\partial
V}{\partial{n}}c\right),$ (45)

since $\displaystyle V(x)/k_{B}T=U\left(\xi\right)$, it means that we have the
following relation between the two expressions of the potentials

$\displaystyle\frac{\partial(V/k_{B}T)}{\partial{n}}=\frac{1}{\varepsilon}U^{\prime}(\xi)$

We notice, as in the 1D case, that there is also a dependence on $\xi$ of the
solution inside the layer, such that
$c=c(\vec{x}+\varepsilon\xi\hat{n},t)=\tilde{c}(\xi,t)$, thus we can put
$1/\varepsilon$ as common factor:

$J_{n}=-D\frac{1}{\varepsilon}\left(\frac{\partial
c}{\partial\xi}+\frac{\partial U}{\partial\xi}c\right)$ (46)

where we omitted the $\tilde{\cdot}$. Following the same argument we used in
1D to derive Eq. (22), we perform a formal expansion in $\varepsilon$ of the
solution. To the lowest order in $\varepsilon$ for the normal component of the
flux we have $\displaystyle{J_{n}^{(0)}}=0$, which implies

$\frac{\partial c^{(0)}}{\partial\xi}+\frac{\partial
U}{\partial\xi}c^{(0)}=0.$ (47)

Therefore we obtain the expression for the concentration to order zero in
$\varepsilon$:

$c^{(0)}(\vec{x}+\varepsilon\xi\hat{n})=c^{(0)}(\vec{x}+\varepsilon
L\hat{n})\exp\left(-U(\xi)\right).$ (48)

Using this expression, we can compute the line density (2D) or the surface
density (3D) of entrapped ions:

$\mathcal{C}(\vec{x})=\varepsilon\int_{0}^{L+1}c(\vec{x}+\varepsilon\xi\hat{n})d\xi\simeq
Mc(\vec{x})$ (49)

with $\displaystyle M=\varepsilon\int_{0}^{L+1}\exp\left(-U(\xi)\right)d\xi$
and where we assumed $U(L+1)=U^{\prime}(L+1)=0$.

At the end, in 2D (3D), the system and the new boundary condition for the
concentration becomes:

$\displaystyle\displaystyle\frac{\partial c}{\partial t}$ $\displaystyle=$
$\displaystyle D\Delta c\quad\rm{in}\>\Omega$ (50)
$\displaystyle\displaystyle\frac{\partial c}{\partial n}$ $\displaystyle=$
$\displaystyle 0\quad\rm{on}\,\partial\Omega\backslash\Sigma_{B}$
$\displaystyle\displaystyle M\frac{\partial c}{\partial t}$ $\displaystyle=$
$\displaystyle MD\Delta_{\perp}c+D\frac{\partial c}{\partial
n}\quad\rm{on}\,\Sigma_{B}$ (51)

because $J_{n}=-D\partial c/\partial n$ just out of the range of the
potential.

Notice that in 2D the Laplace-Beltrami operator reduces to the second
derivative with respect to the arclength of the boundary:111Here by
$\partial/\partial\tau$ we denote the derivative on $\Gamma_{B}$, i.e. the
derivative along the arclength that parametrizes the curve, likewise
$\partial^{2}/\partial\tau^{2}$ denotes the second derivative along
$\Gamma_{B}$, not the second derivative along the tangent direction. See Fig.
3 and Appendix B.

$\Delta_{\perp}c=\frac{\partial^{2}c}{\partial\tau^{2}}.$ (52)

##  3 Validation of the model and results

In this section we validate the model in 1D and 2D, considering different
geometries. In all numerical tests, when dealing with the full model, we
choose a space discretization that fully resolves the $\varepsilon$-scale.

We start considering the effect of the LJ potential in 1D and the agreement
with the multiscale model. To validate the model we compare the mass entrapped
in the hole for the full model, with the concentration calculated in the
boundary and compare the profiles of the fully resolved solution and of the
solution of the multiscale model in the bulk.

Secondly, we consider a 2D system in a squared domain with a potential that is
constant in the $y$-direction and concentrated on the left edge of the domain.

Finally we consider a more realistic domain in 2D, with a circular hole
describing the attractive bubble surface at the center of a squared domain.

###  3.1 Validation 1D

The domain in 1D is
$\Omega^{\varepsilon}=[-\varepsilon,1]=\Omega^{\varepsilon}_{b}\cup\Omega^{\varepsilon}_{f},\,\Omega^{\varepsilon}_{b}=[-\varepsilon,L\varepsilon],\,\Omega^{\varepsilon}_{f}=[L\varepsilon,1]$.
In all our calculations we used $L=2$, as we explain in Section 2.1. The
computational domain $\Omega^{\varepsilon}_{h}$ is then discretized by a
uniform Cartesian mesh with spatial step $h$. The concentration $c_{i}$ is
defined at the center of the cells with
$x_{i}=-\varepsilon+(i-1/2)h\in\\{1,\dots,N_{x}\\},hN_{x}=1+\varepsilon$. We
choose a cell centered discretization in order to guarantee the exact
conservation of the total volume of the solute $v=\sum_{i}c_{i}h$ deriving by
the zero boundary condition for the flux. The scheme is second order accurate
and it is stable even in presence of a drift term, provided the so called
_mesh Péclet number_ is smaller than 2 [24], defined as follows:

${\rm pec}\equiv\max_{x}\left|{\partial_{x}U}\right|h<{2}.$ (53)

For the multiscale model the domain is $\Omega^{0}=[0,1]$, and again we
discretize the computational domain $\Omega^{0}_{h}$ as follows:
$x_{i}=(i-1/2)h\in\\{0,\dots,N\\},hN=1$. Internal equations are discretized by
finite differences by the classical three points approximation of the second
derivative, while the Eq. (31) is discretized to second order as

$M\frac{\partial}{\partial t}\frac{c_{0}+c_{1}}{2}=D\frac{c_{1}-c_{0}}{h}.$

In this case we do not have restriction defined in Eq. (53) on the space step
because $U=0$ in $\Omega^{0}$.

If we denote with $c^{\varepsilon}$ and $c^{0}$ the solutions of the full and
multiscale models, respectively, the problems (10)-(11) and (29)-(32) can be
summarized as

$\frac{\partial c^{\varepsilon}}{\partial
t}=Q^{\varepsilon}(c^{\varepsilon}),$ (54)

where we distinguish between the cases $\varepsilon>0$ and $\varepsilon=0$:

$\displaystyle Q^{\varepsilon}(c^{\varepsilon})=\left\\{\begin{matrix}D\Delta
c^{\varepsilon}&{\rm in}\,\Omega^{\varepsilon}_{f}\\\ \displaystyle
D\left(\Delta
c^{\varepsilon}+\frac{1}{k_{B}T}\nabla\cdot(c^{\varepsilon}V^{\prime})\right)&\rm{in}\,\Omega^{\varepsilon}_{b}\end{matrix}\right.\qquad{\rm
if}\>\varepsilon>0$ (55) $\displaystyle
Q^{0}(c^{0})=\left\\{\begin{matrix}D\Delta c^{0}&{\rm in}\,\Omega^{0},\\\
\displaystyle DM^{-1}\frac{\partial c^{0}}{\partial x}\quad&{\rm
on}\,x=0\end{matrix}\right.\qquad{\rm if}\>\varepsilon=0$ (56)

Eq. (54) is discretized in time by using the Crank-Nicolson method, which is
second order accurate:

$\displaystyle\frac{c^{\varepsilon,n+1}-c^{\varepsilon,n}}{\Delta t}$
$\displaystyle=$
$\displaystyle\frac{1}{2}Q^{\varepsilon}(c^{\varepsilon,n})+\frac{1}{2}Q^{\varepsilon}(c^{\varepsilon,n+1})$
(57)

where $\Delta t$ is the time step, and we solve the system for
$c^{\varepsilon,n+1}$. Eq. (54) is supplemented by the following initial
condition:

$c(t=0,x)=c_{\rm
in}(x)=\frac{v_{0}}{\sqrt{2\sigma^{2}\pi}}\exp\left(-\left(x-x_{m}\right)^{2}/2\sigma^{2}\right),\qquad\sigma\in\mathbb{R},\quad
x_{m}\in(0,1),$ (58)

and $v_{0}$ denotes the total volume.

The method that we consider takes into account two parameters, the range of
the potential $\varepsilon$ and the discretization step $h=\Delta x$ (we chose
$\Delta t\propto\Delta x$). In this section first we see the effect of the LJ
potential, described in Eq. (12), on the concentration (see Fig. 4, left). As
$\varepsilon$ decreases, the agreement between the two models gets better and
better, as we see in Fig. 4 (right).

In order to perform a quantitative comparison between the two models, we
calculate the difference in the amount of solute near the surface, which we
denote by $e_{S}$, and the norm of the concentration difference in the bulk,
denoted by $e_{B}$ (scaled by the total amount of the solute $v_{0}$), as
functions of time:

$\displaystyle\displaystyle{\rm
e_{S}}=\frac{\left|\int_{-\varepsilon}^{\varepsilon L}c^{\varepsilon}(x,t)dx-
Mc^{0}(x=0,t)\right|}{v_{0}},\quad{\rm e_{B}}=\frac{\int_{\varepsilon
L}^{1}\left|c^{\varepsilon}(x,t)-c^{0}(x,t)\right|dx}{v_{0}}$ (59)
$\displaystyle{\rm
error}=\frac{|\int_{-\varepsilon}^{L\varepsilon}c^{\varepsilon}(x)dx-
Mc^{0}(x=0,t_{\rm final})|}{|Mc^{0}(x=0,t_{\rm final})|}$ (60)

The results are reported in Fig. 6. As we can see, the discrepancy between the
two approaches is approximately proportional to the range of the potential
$\varepsilon$.

Finally in 1D we show that the error between models does not depend on the
space discretization. In Fig. 7 we see the error is almost independent of the
space step $\Delta x$, especially when $\Delta x\ll\varepsilon$.

Figure 4: Here we see two aspects of the same plot. On the left we see the
effect of the attractive potential that becomes stronger decreasing the value
of $\varepsilon$. We consider a zoom-in in $x$-direction to focus on the
length of the entrapped mass. In the plot on the right we have a zoom-in in
$y$-direction with an additional case. The dashed red line represents the
solution of the multiscale model, $c^{0}$. In this panels we see the agreement
of the two models improves with $\varepsilon\to 0$. Parameters of the tests:
$\Delta x=10^{-4},\,t_{\rm final}=0.05,\,M=3$ and for the initial condition
defined in Eq. (58) we choose $\sigma=0.1$, $v_{0}=10^{-6}$, $x_{m}=0.5$.
Parameters of the tests: $\Delta x=10^{-4},\,t_{\rm final}=0.05,\,M=3$ and for
the initial condition defined in Eq. (58) we choose $\sigma=0.1$,
$v_{0}=10^{-6}$, $x_{m}=0.5$.

Figure 5:  Here we see the effect of the potential with two different values
of $M$, on the left $M=5$ and on the right $M=10$. We see that for a higer
value of $M$, (but same values of $\varepsilon$) the peak of the solution is
higher. Parameters of the tests: $\Delta x=10^{-4},\,t_{\rm final}=0.05,$ and
for the initial condition defined in Eq. (58) we choose $\sigma=0.1$,
$v_{0}=10^{-6}$, $x_{m}=0.5$.

Figure 6: Errors defined in Eq. (59) for the entrapped mass near the surface
(left panel) and in the bulk (right panel). We see the errors decrease with
$\varepsilon$ as expected. Parameters of the simulations: $\Delta x=1.82\times
10^{-4},\,t_{final}=0.05$ and for the initial condition defined in Eq. (58):
$\sigma=0.2$, $v_{0}=10^{-6}$, $x_{m}=0.5$.

Figure 7: In this figure we show how the error, defined in Eq. (60), does not
depend on the spatial discretization, once the $\varepsilon$-scale is
resolved. It also shows that error increase linearly with $\varepsilon$. This
is why we have the same behaviour for all values of $\Delta x$ (left panel)
and the errors are constants as functions of $\Delta x$ (right panel).

###  3.2 Validation in 2D

In this section we define a very simple domain in 2D. It is a direct extension
of the 1D case, with the potential defined as a function $V:=V(x,y)$ but it is
constant in $y$-direction. In this case the domain is defined as
$\Omega^{\varepsilon}=[-\varepsilon,1]\times[0,1]$ and the potential is
negligible for $x>L\varepsilon$ with $L=2$ as we do in 1D.

The computational domain $\Omega^{\varepsilon}_{h}$ is then discretized by a
uniform Cartesian mesh with spatial steps $\Delta x$ and $\Delta y$ in $x$ and
$y$ directions, respectively. The concentration $c_{ij}$ is defined at the
center of the cells with $x_{i}=-\varepsilon+(i-1/2)\Delta
x\in\\{1,\dots,N_{x}\\},\Delta xN_{x}=1+\varepsilon$ and $y_{j}=(j-1/2)\Delta
y\in\\{1,\dots,N_{y}\\},\Delta yN_{y}=1$. We choose again a cell centered
discretization, and we make sure that the stability condition imposed by the
_mesh Péclet number_ always holds.

For the multiscale model the domain is $\Omega^{0}=[0,1]\times[0,1]$, and
again we discretize the computational domain $\Omega^{0}_{h}$ as follows:
$x_{i}=(i-1/2)\Delta x\in\\{0,\dots,N\\},\Delta xN=1$ and $y_{i}=(i-1/2)\Delta
x\in\\{1,\dots,N\\},\Delta yN=1$, and $h=\Delta x=\Delta y$.

About the discretization in time, in both cases, full and multiscale model, we
choose the Alternate Direction Implicit (ADI) method [24], which is basically
as accurate as Cranck-Nicholson, but it is much mode efficient, since it
requires to solve only tri-diagonal systems at each time step. In Fig. 8
panels (a,b,c) we see the effect of the potential $V$ as $\varepsilon\to 0$
and in panel (d) we see a one dimensional plot as a summary of the three
previous cases ( black, blue and green lines) and the solution of the
multiscale model (dashed red line).

Figure 8: We see the effect of the attractive potential that becomes stronger
and stronger decreasing the value of $\varepsilon$ in panels (a,b,c). In panel
(d) we have a 1D plot, that is a section of the previous one for $y=0.5$ with
an additional line. The dashed red line represents the solution of the
multiscale model, $c^{0}$. The agreement of the two models improves with
$\varepsilon\to 0$. Grid resolution of the simulations: $\Delta x=\Delta
t=3.33\times 10^{-4}$, $t_{final}=0.05s$ and the initial condition is defined
in Eq. (58) with $x_{m}=y_{m}=0.5,\sigma=0.1,v_{0}=10^{-6}$.

###  3.3 Validation 2D in presence of a circular bubble

Here we consider a more realistic 2D domain, in which the bubble is described
as a circular disk centered at the origin, for the multiscale model, while for
the full problem the potential is defined as a function
$V(x,y):=\mathcal{U}\left({(r-R)}/{\varepsilon}\right)$ where $\mathcal{U}$ is
a function $\mathcal{U}\colon\mathbb{R}\to\mathbb{R}$ such that

$\displaystyle\mathcal{U}(\xi)=\left\\{\begin{array}[]{l}+\infty\quad\xi<-1\\\
U(\xi)\quad\xi\in[-1,L]\\\ 0\quad\quad\xi>L\end{array}\right.$

with

$\displaystyle\xi=\frac{r-R}{\varepsilon}=\frac{\sqrt{x^{2}+y^{2}}-R}{\varepsilon}$

where $R$ is the radius of the bubble and $r$ is the distance from the origin,
and $U(\xi)$ is given by Eq. (12). The equation for $c=c^{0}$ in the domain
$\Omega^{0}$ is the following

$\frac{\partial c}{\displaystyle\partial t}=\nabla\cdot\left(D\nabla
c\right)\text{ in }\Omega^{0}.$ (61)

The fluid is contained in a box
$\mathcal{S}=(-a,a)^{2}\subseteq\mathbb{R}^{2}$, with $a>0$. Wall boundary
conditions are prescribed on the boundary
$\Gamma_{\mathcal{S}}=\partial\mathcal{S}$ (zero Neumann conditions for the
solution):

$\frac{\partial c}{\partial{n}}=0\quad\text{ on }\Gamma_{\mathcal{S}}.$ (62)

Figure 9: Representation of the domain on the left and classification of
inside grid points (green points), ghost points (red squares) and inactive
points (blue circles) on the right.

About the discretization in space for the full problem, it is analogue to the
previous case, described in 3.2 and the expression for the initial condition
is the following:

$c(t=0,x,y)=c_{\rm
in}(x,y)=\frac{2v_{0}}{\sqrt{2\sigma^{2}\pi}}\exp\left(\left(-\left(x-x_{m}\right)^{2}+\left(y-y_{m}\right)^{2}\right)/2\sigma^{2}\right),\quad\sigma\in\mathbb{R}$
(63)

In the presence of a bubble, the fluid domain is represented by
$\Omega^{0}=\mathcal{S}\backslash\mathcal{B}$, where $\mathcal{B}$ is the
region occupied by the bubble and represented by a circle centred in the
origin and with radius $R$ such that $0<R<a$. We choose
$(x_{m},y_{m})\in\Omega^{0}$.

The computational domain $\mathcal{S}$ is discretized through a uniform
Cartesian mesh with spatial step $h=2a/N=\Delta x=\Delta y$, where $N^{2}$ is
the number of cells. We use a cell-centered discretization, therefore the set
of grid points is
$\mathcal{S}_{h}=\\{(x_{i},y_{j})=(-a-h/2+ih,-a-h/2+jh),(i,j)\in\\{1,\cdots,N\\}^{2}\\}$.
Within the set of grid points we define the set of internal points
$\Omega^{0}_{h}=\mathcal{S}_{h}\cap\Omega^{0}$, the set of bubble points
$\mathcal{B}_{h}=\mathcal{S}_{h}\cap\mathcal{B}$ and the set of ghost points
$\mathcal{G}_{h}$:

$(x_{i},y_{j})\in\mathcal{G}_{h}\iff(x_{i},y_{j})\in\mathcal{B}_{h}\text{ and
}\\{(x_{i}\pm h,y_{j}),(x_{i},y_{j}\pm h)\\}\cap\Omega^{0}_{h}\neq\emptyset.$
(64)

The remaining grid points that are neither inside nor ghost points are called
inactive points. See Fig. 9 (right panel) for a classification of inside,
ghost, and inactive points.

For an accurate description of the discretization of the problem we refer to
[25], where we present a ghost-point method to solve drift-diffusion equations
in domains with circular holes, in 2D and 3D, with boundary conditions defined
in Eq. (32), even in presence of a moving bubble. A suitable geometric
multigrid approach is adopted to solve the problem, based on a proper
relaxation of the boundary conditions on the ghost points, with a relaxation
parameter chosen in order to achieve the convergence of the iterative method
[20]. Here we briefly describe the space and time discretization.

The problem (57) is discretized in space, leading to a linear system

$\left(I_{h}-\frac{k}{2}Q_{h}\right)c_{h}^{n+1}=\left(I_{h}+\frac{k}{2}Q_{h}\right)c_{h}^{n}$
(65)

to be solved at each time step, where $I_{h}$ and $Q_{h}$ are the
$(N_{I}+N_{G})\times(N_{I}+N_{G})$ matrices representing the discretization of
the identity operator $I$ and the operator $Q$ in Eqs. (55-56).

While for the full model we use an ADI discretization in time, for the
multiscale model we use a Crank-Nicolson method. The reason of this choice
comes from the shape of $\Gamma_{B}$, that makes the discretization in time
depending on the coordinates, losing the benefits of the ADI scheme.

Figure 10: Representation of the five-point stencil for the discretization of
internal points $P_{ij}=(x_{i},y_{j})$ (left panel) and the reduced stencil
when $P_{ij}$ is close to the wall $\Gamma_{\mathcal{S}}$ (right panel). In
the latter case the stencil is composed by four points.

If $P_{ij}=(x_{i},y_{j})\in\Omega^{0}_{h}$ is an internal grid point (as in
Fig. 10, left panel), then the equation of the linear system is obtained from
the discretization of the internal Eq. (61) and the standard central
difference on a five-point stencil is used to discretize the Laplace operator
on $(x_{i},y_{j})$. If $G=(x_{i},y_{j})\in\mathcal{G}_{h}$ is a ghost point,
then we discretize the boundary condition (47), following a ghost-point
approach similar to the one proposed in [20].

Figure 11: In these plots we see the effect of the attractive potential in a
more appropriate 2D case. Here the bubble has a reasonable circular shape. As
before, in the first three panels (a,b,c) we show the solution
$c^{\varepsilon}$ for three different values of $\varepsilon$. In panel (d) we
have a 1D plot, that is a section of the previous ones for $y=0.0$ with two
additional lines. The dashed black line represents the solution of the
multiscale model, $c^{0}$, with the same space resolution of the other three
cases. The dashed red line represents $c^{0}$ with a larger space step. The
agreement of the two models improves with $\varepsilon\to 0$. Grid resolution:
$\Delta x=0.001$ and $\Delta x=0.0025$ for the dashed red line,
$t_{final}=0.05$ and $M=3$. The initial condition is defined in Eq. (63) with
$x_{m}=1$, $y_{m}=0,\sigma=0.1$ and $v_{0}=10^{-6}$.

##  4 Modeling the saturation effect

In this section we describe how to extend the model to include nonlinear
effects due to the saturation when the concentration is not negligibly small.
Such effect may be relevant near the trap, because that is where the
surfactant concentrates.

###  4.1 Saturation: 1D model

We derive a boundary-problem from the drift diffusion non linear equation
using the same asymptotic analysis as in the case of a single species
multiscale model. Starting with a one-dimensional model we impose that the
coefficient of $O(\varepsilon^{-1})$ in the expansion of the flux is zero:

$\frac{\partial c(x)}{\partial x}+\frac{1}{k_{B}T}V^{\prime}(x)c(x)(1-c(x))=0$
(66)

This equation can be integrated giving

$c(x)=\frac{c_{B}\tilde{f}(x)}{1-c_{B}(1-\tilde{f}(x))}$ (67)

where $c_{B}=c(\varepsilon L)$ is the concentration at $x=\varepsilon L$, i.e.
at the right boundary of the support of the potential and

$\tilde{f}(x)=\exp\left(\frac{1}{k_{B}T}(V(\varepsilon L)-V(x))\right)$ (68)

denotes the classical Boltzmann factor. The expression for the entrapped ion
surface concentration $C_{B}$ is:

$C_{B}(c_{B})=\int_{-\varepsilon}^{\varepsilon L}c(x)\,dx=M(c_{B})c_{B}$ (69)

where now the coefficient $M$ in front of $c_{B}$ depends on the concentration
$c_{B}$ just out of the potential well:

$\displaystyle M(c_{B})$
$\displaystyle=\varepsilon\int_{0}^{L+1}\frac{f(\xi,\phi)}{1-c_{B}(1-f(\xi))}\,d\xi,\quad
f(\xi)\equiv\tilde{f}(x)$ (70) $\displaystyle f(\xi,\phi)$
$\displaystyle=\exp\left(-U(\xi;\phi)\right)$ (71)

where we assume $U(L+1)=0$, and

$\mathcal{I}_{L}(\phi,c_{B})\equiv\frac{M(c_{B})}{\varepsilon}=\int_{0}^{L+1}\frac{f(\xi,\phi)}{1-c_{B}+c_{B}f(\xi,\phi)}\,d\xi$
(72)

Notice that now we cannot choose $M$ arbitrarily, since the integrand in (72)
is now bounded by $1/c_{B}$. In particular

$\mathcal{I}_{L}(\phi,c_{B})<\frac{L+1}{c_{B}}.$

This means that when taking into account the effect of saturation, two
potentials with the same value of $M(0)$ will have different behavior, since a
potential with a wider and less deep well can accommodate more particles, and
therefore the effect of saturation will be less evident than in the case of a
potential with the deeper and narrower well.

The quantity $\mathcal{I}(\phi,c_{B})$ approximately represents the average
amplification of the concentration in the potential well. In Fig. 12 (left) we
plot the amplification $\mathcal{I}(\phi,c_{B})$ for a wide range of
concentrations.

Figure 12: On the left we plot $\mathcal{I}(\phi,c_{B})$ for different values
of $\phi$. If $M$ is fixed, when $\phi$ increases, $\varepsilon$ decreases.
The result is that $\mathcal{I}(\phi,c_{B})$ is constant when $c_{B}\to 0$ and
it increases as $\varepsilon$ decreases. On the right we plot
$M(c_{B})/\varepsilon$ for $\phi=10$ (thick black line), together with the
Taylor expansion (73), truncated at various order $k$. Also in Table 2 we show
the oscillating behaviour of these quantities. The value $L=4$ has been used
here.

For low densities, this factor increases a lot with the parameter
$\phi=E/(k_{B}T)$, while when the density increases, the amplification is
reduced, and becomes the width of the well as $c_{B}$ approaches 1.

For moderate but non negligible values of $c_{B}$ a truncated power series
expansion of $M(c_{B})$ could be used:

$M(c_{B})=\varepsilon\sum_{k=0}^{\infty}(c_{B})^{k}M_{k},$ (73)

with

$M_{k}=\int_{0}^{L+1}\,f(\xi)(1-f(\xi))^{k}\,d\xi.$ (74)

However such an expression will not be particularly useful because the
coefficients $M_{k}$ increase very rapidly with $k$. In Table 2 we see a few
values of $M_{k}$ for $k=0,\ldots,4$. In Fig. (12) (right panel) we show the
expression of $M(c_{B})/\varepsilon$ for $\phi=10$, together with its Taylor
expansion for several values of the order $k$ of the Taylor expansion. It is
evident that such an expansion is not useful for concentrations bigger that
about $3\times 10^{-5}$.

$k$ | 0 | 1 | 2 | 3 | 4
---|---|---|---|---|---
$M_{k}$ | 2.9065e+04 | -7.0352e+09 | 1.9534e+15 | -5.7704e+20 | 1.7634e+26
Table 2: The first few coefficients $M_{k}$ (Eq. (74)) of the power series
(73). $\varepsilon=10^{-4}$.

A more useful approximation is given by a rational function where numerator
and denominator are both polynomials. In Fig. 13 (left) we see the comparison
between the exact value of $\mathcal{I}(\phi,c_{B})$ and a rational
approximation. Of course better approximation is obtained using rational
functions of higher degree.

Figure 13: On the left we plot $M(c_{B})/\varepsilon$ for $\phi=10$ (blue
line), together with the rational approximation $R_{3,4}(c)=P_{3}(c)/Q_{4}(c)$
(dashed red line). In order to have a better approximation we can use rational
fractions with higher degree. On the right we show the values of $c_{\rm
max}(c_{B})$ for $\phi\in\\{6,8,10,12,14\\}$. If $M$ is fixed and we increase
$\phi$, the value of $\varepsilon$ goes to 0 and the maximum value of the
concentration, $c_{\rm max}(\phi,c_{B})$, goes very close to 1.

Once the trapped charge $C_{B}$ is expressed in terms of $c_{B}$ near the
boundary, the boundary condition for the multiscale model with non negligible
saturation effect becomes:

$\displaystyle\frac{dC_{B}}{dt}+\left.D\frac{\partial c}{\partial
n}\right|_{x=0}=0$ (75)

where $C_{B}=M(c_{B})c_{B}$.

###  4.2 Results

In this section we show various tests to validate the saturation effect
described in Eq. (70).

In Fig. 12 (left panel) we show how the quantity $\mathcal{I}_{L}(\phi,c_{B})$
depends on $c_{B}$ for different values of $\phi$. We can easily see that it
does not depend on $c_{B}$ when $c_{B}$ is sufficiently low and also it
becomes stronger and stronger when we increase the value of $\phi$.

The product $c_{\max}(\phi,c_{B})\equiv c_{B}\mathcal{I}_{L}(\phi,c_{B})$
provides the concentration near the bubble. In Fig. 13 (right panel) we show
how such concentration reaches 1 very easily when $\phi$ increases. We plot
the quantity $c_{B}\mathcal{I}_{L}(\phi,c_{B})$, and, as expected, we see the
saturation effect becomes more effective when $\phi$ increases. In Fig. 14 we
see the fraction of entrapped mass for the Eqs.(75) versus time, for different
values of $\varepsilon$. In this figure we see that the concentration
saturates earlier when $\varepsilon$ decreases.

Figure 14: Left: Here we show the fraction of entrapped mass versus time for
different values of $\varepsilon$. We see the solution saturates at earlier
times when $\varepsilon$ decreases. Grid resolution chosen: $\Delta x=10^{-4}$
and for the initial condition defined in Eq. (58) we choose
$v_{0}=10^{-6},\sigma=0.2$ and $x_{m}=0.5$.

Figure 15: Profiles of the solution at the circular boundary, as function of
the angle $\theta$. We see the saturation effect is almost negligible for very
low concentrations, but it becomes more evident increasing the initial volume,
$v_{0}$. In these plots we show the results of the solution with (black line)
and without (dashed line) saturation. Grid resolution: $\Delta x=0.002$,
$t_{final}=0.15$ and $M(c_{B}=0)=3$. The values of the initial volume are
described in each plot.

##  5 Conclusions

In this paper we present a multiscale model for the description of the
interaction between a reversible trap and a surfactant. The starting point is
a drift-diffusion model for a single carrier, in which particles are driven by
the gradient of the potential describing the effect of the bubble. The
multiscale model is based on the assumption that the range of the potential is
much smaller than typical macroscopic length scales, such as, for example, the
radius of a spherical trap. It is shown that near the trap the particles are
in local equilibrium, and therefore the concentration satisfies a Boltzmann-
type distribution. This allows a drastic simplification of the system, which
is then described by just a diffusion equation, with a peculiar evolutionary
time dependent solution derived by the continuity equation of the trapped
charge and the relation of local equilibrium, and which does not require the
resolution of the small scales of the potential.

The multiscale model solved on a uniform grid provides a tremendous reduction
in the computational complexity compared to the full model, even if the latter
is solved on using Adaptive Mesh Refinement, as is clearly illustrated in
Appendix A.

Saturation effect is also taken into account. The model is carefully validated
numerically against a detailed numerical solution of a fully resolved model
with a potential of width $\varepsilon$. It is shown that the agreement
between the multiscale model and the original drift-diffusion model becomes
better and better as $\varepsilon$ is decreased. As basic potential we adopt
Lennard-Jones type potential, which is uniquely determined once one specifies
the depth of the potential well, that we call $E$, and its width
$\varepsilon$. We observe that for very dilute solutions the effect of the
potential is summarized by a single quantity that is a function of $E/k_{B}T$
and $\varepsilon$, because the saturation effect is negligible (for a fixed
$\varepsilon$ and suitably small concentration). However we show that the
effect of saturation may become relevant even for moderate concentration,
since near the trap the concentration increases by orders of magnitude. This
means that one has to provide realistic values of $\varepsilon$ (that we call
$\tilde{\varepsilon}$ to distinguish it from the mathematical parameter used
in the asymptotic expansion of the solution) and $E$, and the model is able to
capture the behavior of the diffusant without resolving the small scale
$\tilde{\varepsilon}$. This multiscale model is applied to study the evolution
of the surfactant in presence of a time dependent bubble [25]. A physically
realistic model should be based on the transport equation for both cations and
ions, coupled by a self-consistent Poisson equation for the Coulomb potential
[6]. The current multiscale model could be used to describe carrier evolution
under the condition of the so-called quasi-neutral limit [6, 26], in which the
Coulomb interaction is so strong that the two carrier charge densities are
practically overlapped, and move as a single carrier with effective mobility
and diffusion coefficient. A multiscale two-carrier model which takes into
account Coulomb interaction far from the quasi-neutral limit is current under
investigation.

##  Appendix A Comparison with Adaptive Mesh Refinement

In the derivation of the multiscale model we perform an expansion in the
parameter $\varepsilon$ which in our paper represents the non-dimensional
thickness of the potential well near the bubble, and is given by the ratio of
the actual thickness of the potential well and a typical size of the device.
We can imagine that the device is of the order of a few millimeters, while a
reasonable thickness of the potential well is of the order of nanometers, this
giving a non-dimensional thickness of about $\tilde{\varepsilon}\approx
10^{-6}$. Following the argument in [27], in this appendix we give a rough
estimate of the number of degrees of freedom which is necessary to fully
resolve all the scales of the problem, comparing two different approaches: one
based on the original drift-diffusion model in which a small grid is used only
near the bubble to resolve the small scale of the thin potential of the
bubble, while a coarse grid is adopted in the bulk, where the evolution is
governed by a diffusion equation, and the other one which is based on the
multiscale model, in which the effect of the potential is modeled by the
evolutionary boundary condition for the trapped charge, and a uniform coarse
grid is adopted in the domain. We consider the problem in 1D, 2D and 3D. Let
us assume that the bulk is a domain discretized by a grid of (non-dimensional)
mesh size $H$, while the region influenced by the potential has to be
discretized by a grid of meshsize $h$. Let us denote by $N=1/H$ the number of
grid points per direction needed to fully resolve the bulk, and
$n=\tilde{\varepsilon}/h$ the number of grid points required to resolve the
thickness of the potential. For simplicity we assume the domain is a unit cube
in dimension $d$, i.e. $\Omega=[-0.5,0.5]^{d}$, and the bubble is a $d$-sphere
of radius $R<0.5$.

  1. 1.

In one space dimension one needs $n$ grid points to discretize the region
influenced by the potential, and $N$ grid points to discretize the bulk.

  2. 2.

In two space dimensions we assume we have small regions of size $h^{2}$ inside
the potential and large regions of size $H^{2}$ in the bulk. The number of
regions inside the potential is then given by the ratio between the area
influenced by the potential and the area of the small cell: $2\pi
R\tilde{\varepsilon}/h^{2}$. On the other hand, the number of elements needed
to model the boundary conditions in the multiscale model is given by $2\pi
R/H=2\pi RN$. The number of elements in the bulk is approximately $N^{2}$ in
both cases.

  3. 3.

In three space dimensions the number of elements in the region influenced by
the potential is given by the volume of such region, $4\pi
R^{2}\tilde{\varepsilon}$, divided by the volume $h^{3}$ of each small
element.

We summarize in Table 3 the number of degrees of freedom which are needed to
fully resolve the problem in the two approaches (AMR and Multiscale). Here we
neglected the effect of factors $2\pi R$ and $4\pi R^{2}$ which are $O(1)$.

d | AMR #dof | Multiscale #dof
---|---|---
1 | $n+N$ | $110$ | $1+N$ | $101$
2 | $n^{2}/\tilde{\varepsilon}+N^{2}$ | $10^{8}$ | $N+N^{2}$ | $10^{4}$
3 | $n^{3}\tilde{\varepsilon}^{2}+N^{2}$ | $10^{15}$ | $N^{2}+N^{3}$ | $10^{6}$
Table 3: Order of magnitude of the number of degrees of freedom needed to
resolve the problem in dimension $d$, using either full model discretized by
adaptive mesh refinement (AMR) or discretized by the new multiscale model. The
number in columns 3 and 5 refer to the specific assumption $n=10$, $N=100$ and
$\tilde{\varepsilon}=10^{-6}$.

From these considerations it follows that the multiscale model allows an
improvement of efficiency by orders of magnitude over a detailed direct
computation performed by resolving all the small scales of the problem.

##  Appendix B Laplace-Beltrami operator

Here we prove that

$\partial_{i}\tilde{\partial}_{i}c=\tilde{\partial}_{i}\tilde{\partial}_{i}c=:\Delta_{\perp}c$
(76)

where $\Delta_{\perp}$ is the Laplace-Beltrami operator. We also provide a
representation of $\Delta_{\perp}$ in Cartesian coordinates.

First, observe that

$\partial_{i}=\tilde{\partial}_{i}+n_{i}n_{j}\partial_{j},$

therefore

$\partial_{i}\tilde{\partial}_{i}c=\tilde{\partial}_{i}\tilde{\partial}_{i}c+n_{i}n_{j}\partial_{j}\tilde{\partial}_{i}c.$

Now consider the term

$n_{i}n_{j}\partial_{j}\tilde{\partial}_{i}c=n_{i}\frac{\partial}{\partial
n}\tilde{\partial}_{i}c.$

Since $n_{i}\tilde{\partial}_{i}c\equiv 0$, it is

$0=\frac{\partial}{\partial n}(n_{i}\tilde{\partial}_{i}c)=\frac{\partial
n_{i}}{\partial n}\tilde{\partial}_{i}c+n_{i}\frac{\partial}{\partial
n}\tilde{\partial}_{i}c.$

Now $\partial n_{i}/\partial n=n_{j}\partial_{j}n_{i}=0$, therefore

$n_{i}\frac{\partial}{\partial n}\tilde{\partial}_{i}c=0$

from which (76) follows.

Now we give a representation for $\Delta_{\perp}$. We start observing that
$\tilde{\partial}_{i}c=h_{ij}\partial_{j}c$ is orthogonal to $\Sigma_{B}$,
since $n_{i}h_{ij}=0$. It is

$\displaystyle\Delta_{\perp}c=\tilde{\partial}_{i}\tilde{\partial}_{i}c$
$\displaystyle=\tilde{\partial}_{i}(h_{ij}\partial_{j}c)=\tilde{\partial}_{i}(\partial_{i}c-n_{i}n_{j}\partial_{j}c)$
(77)
$\displaystyle=h_{ij}\partial^{2}_{ij}c-\tilde{\partial}_{i}n_{i}\frac{\partial
c}{\partial
n}-n_{i}(\tilde{\partial}_{i}n_{j})\partial_{j}c-n_{i}n_{j}\tilde{\partial}_{i}\partial_{j}c$
(78)

The last two terms are zero, because $n_{i}\chi_{ij}=0$ and
$n_{i}\tilde{\partial}_{i}\partial_{j}c=0$.

Therefore it follows:

$\Delta_{\perp}c=\tilde{\partial}_{i}\tilde{\partial}_{i}c=h_{ij}\partial^{2}_{ij}c-\chi_{ii}n_{j}\partial_{j}c$

We recall that the trace of the second fundamental form of a 2D surface in 3D
is the sum of the two principal curvatures, i.e. twice the mean curvature of
the surface.

### Laplace-Beltrami operator in two dimensions

In two space dimensions (see Fig.3) expressions considerably simplify. Here we
prove Eq. (52), i.e.

$\frac{\partial^{2}c}{\partial\tau^{2}}\equiv\tau_{j}\partial_{j}(\tau_{i}\partial_{i}c)=h_{ij}\partial^{2}_{ij}c-\chi_{ii}\frac{\partial
c}{\partial n}\equiv\Delta_{\perp}c.$ (79)

Indeed it is

$\frac{\partial^{2}c}{\partial\tau^{2}}=\tau_{j}\partial_{j}(\tau_{i}\partial_{i}c)=\tau_{i}\tau_{j}\partial^{2}_{ij}c+\tau_{j}(\partial_{j}\tau_{i})\partial_{i}c=\tau_{i}\tau_{j}\partial^{2}_{ij}c-\frac{1}{R}n_{i}\partial_{i}c$
(80)

the latter equality is a consequence of the relation

$\frac{\partial\hat{\tau}}{\partial\tau}=-\frac{1}{R}\hat{n}$

where $R$ denotes the local radius of curvature. Finally we observe that in 2D
it is

$h_{ij}\equiv\delta_{ij}-n_{i}n_{j}=\tau_{i}\tau_{j}$ (81)

as it can be easily shown by inspection checking that

$n_{i}n_{j}+\tau_{i}\tau_{j}=\delta_{ij}$

Using Eqs. (80) and (81) in (79) completes the proof of relation (52).

## References

  * [1] H.C. Berg and E.M. Purcell.  Physics of chemoreception.  Biophysical Journal, 1977.
  * [2] Haim Diamant, Gil Ariel, and David Andelman.  Kinetics of surfactant adsorption: the free energy approach.  Colloids and Surfaces A: Physicochemical and Engineering Aspects, 183-185:259–276, 2001.
  * [3] R. Miller, E.V. Aksenenko, and V.B. Fainerman.  Dynamic interfacial tension of surfactant solutions.  Advances in Colloid and Interface Science, 247:115–129, 2017.  Dominique Langevin Festschrift: Four Decades Opening Gates in Colloid and Interface Science.
  * [4] Reinhard Miller, Paul Joos, and Valentin B. Fainerman.  Dynamic surface and interfacial tensions of surfactant and polymer solutions.  Advances in Colloid and Interface Science, 49:249–302, 1994.
  * [5] F.W Wiegel.  Diffusion and the physics of chemoreception.  Physics Reports, 95(5):283–319, 1983.
  * [6] Antonio Raudino, Antonio Grassi, Giuseppe Lombardo, Giovanni Russo, Clarissa Astuto, and Mario Corti.  Anomalous sorption kinetics of self-interacting particles by a spherical trap.  Communications in Computational Physics, 31(3):707–738, 2022.
  * [7] C.J. Beverung, C.J. Radke, and H.W. Blanch.  Protein adsorption at the oil/water interface: characterization of adsorption kinetics by dynamic interfacial tension measurements.  Biophysical Chemistry, 81(1):59–80, 1999.
  * [8] A. F. H. Ward and L. Tordai.  Time-dependence of boundary tensions of solutions i. the role of diffusion in time-effects.  jcp, 1946.
  * [9] J. Liu and U. Messow.  Diffusion-controlled adsorption kinetics at the air/solution interface.  Colloid and Polymer Science, 2000.
  * [10] Pierre-Henri Chavanis.  Generalized stochastic fokker-planck equations.  Entropy, 17(5):3205–3252, 2015.
  * [11] Nicolas Martzel and Claude Aslangul.  Mean-field treatment of the many-body fokker-planck equation.  Journal of Physics A: Mathematical and General, 34(50):11225–11240, 2001.
  * [12] D. Marin, M. A. Ribeiro, H. V. Ribeiro, and E. K. Lenzi.  A nonlinear fokker-planck equation approach for interacting systems: Anomalous diffusion and tsallis statistics.  Physics Letters A, 382(29):1903–1907, 2018.
  * [13] Matthew J. Simpson, Kerry A. Landman, and Barry D. Hughes.  Multi-species simple exclusion processes.  Physica A: Statistical Mechanics and its Applications, 388(4):399–406, 2009.
  * [14] M J Simpson, B D Hughes, and K A Landman.  Diffusing populations: Ghosts or folks?  Australasian Journal of Engineering Education, 15(2):59–68, 2009\.
  * [15] Giambattista Giacomin and Joel L. Lebowitz.  Phase segregation dynamics in particle systems with long range interactions. i. macroscopic limits.  Journal of Statistical Physics, 87(1-2):37–61, Apr 1997.
  * [16] Kevin Painter and Thomas Hillen.  Volume-filling and quorum-sensing in models for chemosensitive movement.  Canadian Applied Mathematics Quarterly, 10(4):501–544, 2002.
  * [17] Kevin J. Painter.  Continuous models for cell migration in tissues and applications to cell sorting via differential chemotaxis.  Bulletin of Mathematical Biology, 71(5):1117–1147, 2009.
  * [18] S. S. Dukhin, G. Kretzschmar, and R. Miller.  Dynamics of adsorption at liquid interfaces: theory, experiment, application.  Elsevier, Amsterdam, 1995.
  * [19] K. L. Sutherland.  The kinetics of adsorption at liquid surfaces.  Australian Journal of Scientific Research. Series A., 5:683–696, 1952.
  * [20] Armando Coco and Giovanni Russo.  Finite-difference ghost-point multigrid methods on cartesian grids for elliptic problems in arbitrary domains.  Journal of Computational Physics, 241:464 – 501, 2013.
  * [21] Armando Coco and Giovanni Russo.  Second order finite-difference ghost-point multigrid methods for elliptic problems with discontinuous coefficients on an arbitrary interface.  Journal of Computational Physics, 361:299–330, 2018.
  * [22] Douglas H Everett and John C Powl.  Adsorption in slit-like and cylindrical micropores in the henry’s law region. a model for the microporosity of carbons.  Journal of the Chemical Society, Faraday Transactions 1: Physical Chemistry in Condensed Phases, 72:619–636, 1976.
  * [23] Giovanni Russo and Peter Smereka.  A remark on computing distance functions.  Journal of computational physics, 163(1):51–67, 2000.
  * [24] P. Wesseling.  Principles of computational fluid dynamics.  Springer, 2001.
  * [25] Clarissa Astuto, Armando Coco, and Giovanni Russo.  A finite-difference ghost-point multigrid method for multi-scale modelling of sorption kinetics of a surfactant past an oscillating bubble.  Journal of Computational Physics, page 111880, 2023.
  * [26] Ansgar Jungel and Yue-Jun Peng.  A hierarchy of hydrodynamic models for plasmas. quasi-neutral limits in the drift-diffusion equations.  Asymptotic Analysis, 28, 2000.
  * [27] Matteo Semplice, Armando Coco, and Giovanni Russo.  Adaptive mesh refinement for hyperbolic systems based on third-order compact weno reconstruction.  Journal of Scientific Computing, 66(2):692–724, 2016.


WGMwgm QEqe EPep PMSpms BECbec DEde

xxx [1] [cor1]Corresponding author

# The role of intra- and inter-group Matthew effect in the social dilemma of
public goods games

Chaoqian Wang<EMAIL_ADDRESS>Department of Computational and Data
Sciences, George Mason University, Fairfax, VA 22030, USA

###### Abstract

The Matthew effect describes the phenomenon where the rich tend to get richer.
Such a success-driven mechanism has been studied in spatial public goods games
in an inter-group way, where each individual’s social power is enhanced across
all groups. For instance, factors like knowledge can exert an advantage across
various social contexts. In contrast, certain factors, especially local
material goods, only enhance advantages within their current group. Building
on this, we further explore the intra-group Matthew effect where the
enhancement of social power is calculated separately in each group. Our
findings indicate that the intra-group Matthew effect sustains cooperation
more at high productivity, while the inter-group Matthew effect promotes
cooperation at low productivity. Moreover, the mixture of the intra- and
inter-group Matthew effect harms cooperation. This study provides insights
into addressing social dilemmas by adjusting wealth accumulation across
diverse social groups.

###### keywords:

Matthew effect Public goods game Evolutionary game theory Social dilemma

##  1 Introduction

Matthew 25:29 (New Revised Standard Version) states, “For to all those who
have, more will be given, and they will have an abundance; but from those who
have nothing, even what they have will be taken away.” In the real world, it
is evident that wealthy individuals have greater opportunities to become
richer, while the poor have limited prospects of improving their financial
situations. This phenomenon is known as the Matthew effect [1, 2]. The role of
the Matthew effect has been studied in the social dilemma of evolutionary
public goods games [3].

Social dilemmas refer to the conflict between individual and collective
interests [4]. While cooperation (also known as prosocial behavior) benefits
everyone in a group, individual defection can benefit oneself at the expense
of others, resulting in a reduction in average benefits for all group members,
leading to the tragedy of the commons [5]. However, cooperation can thrive in
human society [6], and several explanations have been proposed for this
phenomenon [7], including network reciprocity, which has been studied by
researchers from evolutionary biology to statistical physics [8] and from
computer simulations [9, 10, 11] to theoretical analyses [12, 13, 14].

It is a natural idea to consider the role of the Matthew effect in multiplayer
games instead of two-player games because the former involves the allocation
of public goods and the enhancement of wealth. The emergence of cooperation in
multiplayer games has become increasingly popular over the past few decades
[15, 16, 17], and group interactions [18] play a significant role in higher-
order population structures [19, 20, 21, 22]. The fundamental model type is
the public goods game, in which cooperators contribute to the group, producing
more than the cost, but defectors can share the fruits equally [23]. In other
words, in a classic public goods game group, everyone shares an equal amount
of fruits. However, introducing the Matthew effect mentioned earlier, those
receiving more can share more. It was found that this can have an intriguing
effect on the cooperation level: the Matthew effect in distributing fruits
promotes cooperation but also preserves defection [3, 24]. Similar to this
success-driven mechanism in distribution, other success-driven mechanisms in
different aspects of evolutionary dynamics have been studied, such as success-
driven migration [25, 26, 27, 28], success-driven participation [29], success-
driven group formation [30], and success-driven multigames [31]. Opposite
assumptions from the Matthew effect have also been investigated, such as the
situation where rich people contribute more [32].

A widely accepted notion is that inter-group reciprocity, where an individual
participates in multiple groups and engages in games in each, plays a crucial
role in spatial public goods games [33, 34, 35]. When coupling groups,
apparently similar things can differ [36]. Previous research has investigated
various inter-group mechanisms in public goods games, such as group
competition [37] and inter-group selection [38, 39]. Some studies have also
examined scenarios where individuals use different strategies with different
neighbors [40, 41].

Similarly, as an additional mechanism, the Matthew effect can also be
considered separately in different groups. On the one hand, an individual can
accumulate social power in different groups independently. On the other hand,
an individual can also possess a unified social power through all its groups,
as previous research has shown [3]. More generally, the intra- and inter-group
social power can be blended, and the impact of their proportions on the role
they play can be explored. The idea of differentiating the intra- and inter-
group Matthew effects captures a detail of the real world: people are involved
in multiple social groups and interact with members in each. For example, a
student attends different classes at a university, and their advantages accrue
separately or jointly in various classes, depending on the relevance of the
knowledge. Similarly, a multinational corporation operates in different
countries, and its advantages in each country may accumulate separately or
jointly, depending on the trade policies.

In this study, we establish a model to describe the intra- and inter-group
Matthew effects in the spatial public goods game and investigate their roles
separately and collectively. It should be noted that we use synchronous
updates [9, 10, 42, 43] in this work. Although most previous studies adopt
asynchronous updates, the synchronous approach is more intuitive in describing
specific mechanisms [44, 45]. The results do not differ significantly in
quality [46].

##  2 Model

The spatial public goods game is conducted on an $L\times L$ square lattice,
where each node is occupied by an agent. The agent set is denoted by
$\mathcal{N}$, and its size is the population $N$, which yields
$|\mathcal{N}|\equiv N=L^{2}$. Each agent $i\in\mathcal{N}$ is the center of a
group $\Omega_{i}$, which contains the centering agent $i$ and its four
nearest neighbors $j\in\Omega_{i}\setminus\\{i\\}$. In this work, we set the
group size to $G=5$. Each agent $i$ is involved in $G$ groups, centered on
itself and its neighbors, $\Omega_{j}$ for $j\in\Omega_{i}$.

The payoff received by agent $i$ from the group centered on $j$ at time $t$ is
denoted by $\pi_{i,t}^{j}$, calculated as follows:

$\pi_{i,t}^{j}=M_{i,j,t}r\sum_{k\in\Omega_{j}}s_{k,t}c-s_{i,t}c,$ (1)

where $s_{i,t}$ denotes the strategy of agent $i$ at time $t$. If $s_{i,t}=1$,
agent $i$ cooperates by contributing $c$ ($c>0$) to the group $j$, while if
$s_{i,t}=0$, agent $i$ defects and makes no contribution. The contributions to
group $j$ by players $k\in\Omega_{j}$ are enlarged by the productivity factor
$r$ ($r>1$) and allocated according to the Matthew effect $M_{i,j,t}$, which
we will define later.

The average payoff of agent $i$ received from its $G$ groups at time $t$ is
denoted by $\pi_{i,t}$, yielding

$\pi_{i,t}=\frac{1}{G}\sum_{j\in\Omega_{i}}\pi_{i,t}^{j}.$ (2)

At each time step $t$, the agents’ strategies $s_{i,t}$ are updated
synchronously. For each agent $i$, one of its neighbors
$j\in\Omega_{i}\setminus\\{i\\}$ is selected randomly, and agent $i$ adopts
$j$’s strategy $s_{j,t}$ with probability

$W(s_{i,t+1}\leftarrow
s_{j,t})=\frac{1}{1+\mathrm{e}^{-(\pi_{j,t}-\pi_{i,t})/\kappa}},$ (3)

where $\kappa$ is the noise parameter, set to $\kappa=0.1$ in this work. If
the update is not successful, agent $i$ keeps its current strategy,
$s_{i,t+1}\leftarrow s_{i,t}$.

To define the Matthew effect, we introduce the concept of social power, which
is determined by an agent’s payoff, both within and between groups. First, an
agent’s intra-group social power varies across groups. Specifically, agent
$i$’s intra-group power in group $j$ at time $t+1$ is based on its payoff from
group $j$ at time $t$, as given by Eq. (4):

$P_{i,j,t+1}^{(\text{intra})}=\pi_{i,t}^{j}.$ (4)

Second, an agent’s inter-group social power is the average of its intra-group
power across all of its $G$ groups. In other words, it is simply the agent’s
payoff at time $t$, as expressed by Eq. (5):

$P_{i,t+1}^{(\text{inter})}=\frac{1}{G}\sum_{j\in\Omega_{i}}P_{i,j,t+1}^{(\text{intra})}=\pi_{i,t}.$
(5)

The social power of agent $i$ is a weighted combination of its intra- and
inter-group power, as given by Eq. (6):

$P_{i,j,t+1}=aP_{i,j,t+1}^{(\text{intra})}+(1-a)P_{i,t+1}^{(\text{inter})},$
(6)

where $0\leq a\leq 1$ determines the relative weight of intra-group power
versus inter-group power.

The Matthew effect, which determines the allocation of public goods to agent
$i$ in group $j$, is proportional to the agent’s social power. Specifically,
the higher an agent’s social power, the more public goods are allocated to it.
This is captured by Eq. (7):

$M_{i,j,t+1}=\frac{\mathrm{e}^{wP_{i,j,t+1}}}{\sum_{k\in\Omega_{j}}\mathrm{e}^{wP_{k,j,t+1}}},$
(7)

where $w\geq 0$ determines the strength of the Matthew effect. A larger value
of $w$ corresponds to a stronger Matthew effect, in which agents with higher
payoffs receive a larger share of the public good. When $w=0$, Eq. (7) reduces
to $M_{i,j,t}\equiv 1/G$ for all $i\in\mathcal{N}$, $j\in\Omega_{i}$, and
$t\in\mathbf{N}^{*}$, which corresponds to the classic public goods game.

In our model, the Matthew effect is self-cumulative, as described by Eqs.
(4)–(5). Specifically, the social power at time $t$ affects the allocation of
public goods at time $t+1$, which in turn affects the social power at time
$t+1$ and the allocation at time $t+2$, and so on. This automatic iteration
process obviates the need for storing data from every previous step.

##  3 Results and discussion

The fixed secondary parameters of this study are population size $N=100\times
100=10^{4}$, group size $G=5$, and selection noise $\kappa=0.1$. The main
parameters are productivity $r$, the strength of the Matthew effect $w$, and
the proportion of the intra-group Matthew effects $a$, whose values are
flexible.

At the start of the simulation, each agent $i\in\mathcal{N}$ is assigned a
random strategy of cooperation ($s_{i,0}=1$) or defection ($s_{i,0}=0$), with
the proportion of cooperators being approximately $1/2$. The social power
$P_{i,j,0}$ for each agent in their respective groups $j\in\Omega_{i}$ is
initialized to 0, resulting in $M_{i,j,0}=1/G$ and an equal distribution of
public goods at the first time step.

The system is then allowed to evolve from $t=0$ to $t=10^{5}$, and the
proportion of cooperators $\sum_{i\in\mathcal{N}}s_{i,t}/N$ is calculated at
each time step $t$. The average proportion of cooperators in the last $2\times
10^{4}$ time steps is used to define the cooperation level ($\rho_{C}$) under
the assigned parameter values, $\rho_{C}=\sum_{t=8\times
10^{4}+1}^{10^{5}}\sum_{i\in\mathcal{N}}s_{i,t}/N/(2\times 10^{4})$.

###  3.1 The role of the inter-group Matthew effect

We begin by examining a special case where only the inter-group Matthew effect
is present. When $a=0$, the level of cooperation, $\rho_{C}$, is plotted as a
function of productivity $r$ and the strength of the Matthew effect $w$ in
Fig. 1. To facilitate a comprehensive presentation of the results, we have
rescaled the $w$-axis logarithmically for $w>1$, while the scale in $0<w<1$ is
linear. When $w=0$, the results coincide with those of the classic synchronous
spatial public goods game. However, when $w>0$, we can observe the impact of
the inter-group Matthew effect’s strength by increasing $w$.

Figure 1: Cooperation level $\rho_{C}$ as a function of productivity $r$ and
the strength of inter-group Matthew effect $w$, with fixed intra-group
proportion $a=0$. The $w$-axis scale is linear for $0\leq w\leq 1$ and
logarithmic for $10^{0}\leq w\leq 10^{3}$. The results are qualitatively
similar to the ones in [3]: a strong inter-group Matthew effect promotes the
emergence of cooperation when productivity is low, but also preserves
defection when productivity is high. Meanwhile, a moderate inter-group Matthew
effect disfavors cooperation most.

As illustrated in Fig. 1, the inter-group Matthew effect promotes cooperation
at low productivity ($r\lesssim 4$), with cooperation emerging at $w\gtrsim
1$. However, at intermediate or high productivity ($r\gtrsim 4$), the inter-
group Matthew effect inhibits cooperation. Interestingly, the inhibition of
cooperation is most significant when the inter-group Matthew effect is
moderate. As the inter-group Matthew effect becomes stronger, the cooperation
level rebounds, albeit remaining lower than in the absence of the Matthew
effect. Notably, for intermediate productivity ($4\lesssim r\lesssim 6$),
cooperation disappears entirely when the effect strength is moderate. Our
findings are consistent with previous work [3] that investigated the pure
inter-group Matthew effect and concluded that it promotes cooperation at low
productivity but sustains defection at high productivity, albeit for a
different implementation of our approach in modeling.

###  3.2 The role of the intra-group Matthew effect

Next, we analyze the impact of the intra-group Matthew effect by setting $a=1$
and examine the relationship between the intra-group effect’s strength and the
cooperation level. Fig. 2 illustrates the cooperation level as a function of
productivity and the strength of the intra-group Matthew effect, with the
$w$-axis scaled similarly to Fig. 1. When the intra-group Matthew effect is
present ($w>0$), we can observe its role by increasing $w$.

Figure 2: Cooperation level $\rho_{C}$ as a function of productivity $r$ and
the strength of intra-group Matthew effect $w$, with fixed intra-group
proportion $a=1$. The $w$-axis scale is linear for $0\leq w\leq 1$ and
logarithmic for $10^{0}\leq w\leq 10^{3}$. A strong intra-group Matthew effect
impedes cooperation, and, compared to this, a moderate intra-group Matthew
effect can sustain cooperation more.

In contrast to the inter-group Matthew effect, a certain level of intra-group
Matthew effect ($w\gtrsim 0.2$) can only promote cooperation at high
productivity ($r\gtrsim 5$). At moderate strength ($0.2\lesssim w\lesssim 1$),
neither the intra-group nor the inter-group Matthew effect can promote
cooperation at low productivity, but the intra-group Matthew effect can
sustain higher levels of cooperation at high productivity. At higher strength
($w\gtrsim 10^{0}$), the intra-group Matthew effect cannot promote cooperation
at low productivity as effectively as the inter-group effect, but it can
maintain higher levels of cooperation at high productivity. In conclusion,
while the intra-group Matthew effect cannot promote cooperation at low
productivity, it can retain more cooperation at high productivity compared to
the inter-group Matthew effect. This advantage is most pronounced at a
moderate strength of the Matthew effect.

To explain how the intra-group Matthew effect maintains cooperation better
than the inter-group effect at high productivity, we can imagine the interface
between cooperation and defection clusters and consider the role of the intra-
group Matthew effect. In a public goods game within a single group, defectors
receive higher payoffs than cooperators because of the extra cost $c$ incurred
by cooperators. Consequently, when allocation in the next moment within a
group depends solely on the current payoff, individuals who adopt the
defective strategy are always allocated a larger proportion $M_{i,j,t+1}$ of
the public good than cooperative individuals. Under the high productivity of
the classic model, cooperative clusters expand rapidly. Defectors who can
allocate more public goods in the next moment mostly become cooperators in the
next moment. As a result, they are allocated the same proportion of public
goods in the group facing the defection region and more public goods in the
group facing the cooperation region due to the effect of defection in the
previous moment. This transformation of defectors into cooperators encourages
adjacent defectors to become more willing to transform into cooperators, which
preserves the rapid expansion of cooperative clusters.

###  3.3 The combination of the intra- and inter-group Matthew effect

Finally, we investigate the combined effects of the intra- and inter-group
Matthew effects, where the intra-group social power is represented as a
proportion of $a$ and the inter-group social power as a proportion of $1-a$.
We use Fig. 3 to display the cooperation level as a function of productivity
$r$ and the proportion of the intra-group Matthew effect $a$ at different
strengths of the Matthew effect $w$.

Figure 3: Cooperation level $\rho_{C}$ as a function of productivity $r$ and
the intra-group proportions $a$ at different strength of intra-group Matthew
effect $w$. The values of $w$ are presented in the subtitles. Only the inter-
group Matthew effect plays a role when $a=0$. Only the intra-group Matthew
effect plays a role when $a=1$. The intra-group and inter-group Matthew
effects are mixed when $0<a<1$. Generally speaking, cooperation is disfavored
when the intra- and inter-group Matthew effects are both present.

When the strength of the Matthew effect is moderate [Fig. 3(b), (c), and (d)],
we find that the intra-group Matthew effect has a dominant range where the
cooperation level is high, particularly in the region where $a$ is large. The
cooperation level varies continuously within this region. Similarly, the
inter-group Matthew effect has its dominant range of smaller $a$. The
cooperation level varies sharply on the adjoining interface of the dominant
range of both effects, which is the most prominent in Fig. 3(b).

In contrast, when the strength of the Matthew effect is large, we observe that
the intra-group Matthew effect only manifests when it exists alone ($a=1$). In
other cases where $a<1$, the change in cooperation level is dominated by the
inter-group Matthew effect, which varies continuously at low cooperation
levels. Overall, in each panel of Fig. 3, the cooperation level is high at the
ends of the $a$-axis and low in the middle. These results suggest that a
combination of intra- and inter-group Matthew effects is detrimental to
cooperation.

##  4 Conclusion

In conclusion, this study expands upon the concept introduced in [3] and
examines the impact of the intra- and inter-group Matthew effect separately
and in combination. Our findings indicate that the inter-group Matthew effect
inhibits cooperation at a moderate level but can sustain it to some degree
when it is strong. In comparison, the intra-group Matthew effect encourages
cooperation, particularly when productivity is high, due to the rapid spread
of cooperative strategies at cluster boundaries. However, the combination of
both intra- and inter-group Matthew effects has a negative impact on
cooperation.

The numerical simulations in this study have the potential to provide insights
for policy analysis on wealth distribution. As previously noted, individuals
participate in different games within various social relationships, and the
Matthew effect is a factor that cannot be ignored in the real world. To
address this issue, we suggest adjusting the role of the wealth enhancement
effect based on the current level of social productivity. At low productivity,
the Matthew effect should be transferable between different social
relationships, and success in one group should facilitate success in all
groups. At high productivity, the Matthew effect should not be correlated
across different social relationships, and success in one group should only
increase the likelihood of success in the same group. In any case, we
recommend minimizing the mixing of these two approaches.

## Declaration of competing interest

None.

## Code availability

The Matlab code for the numerical simulations is available upon request from
the corresponding author.

## References

  * Merton [1968] R. K. Merton,  The Matthew effect in science: The reward and communication systems of science are considered.,  Science 159 (1968) 56–63.
  * Perc [2014] M. Perc,  The Matthew effect in empirical data,  Journal of The Royal Society Interface 11 (2014) 20140378.
  * Perc [2011] M. Perc,  Success-driven distribution of public goods promotes cooperation but preserves defection,  Physical Review E 84 (2011) 037102.
  * Sigmund [2010] K. Sigmund, The calculus of selfishness, Princeton University Press, 2010.
  * Hardin [1968] G. Hardin,  The tragedy of the commons,  Science 162 (1968) 1243–1248.
  * Nowak [2006a] M. A. Nowak, Evolutionary dynamics: exploring the equations of life, Harvard University Press, 2006a.
  * Nowak [2006b] M. A. Nowak,  Five rules for the evolution of cooperation,  Science 314 (2006b) 1560–1563.
  * Perc et al. [2017] M. Perc, J. J. Jordan, D. G. Rand, Z. Wang, S. Boccaletti, A. Szolnoki,  Statistical physics of human cooperation,  Physics Reports 687 (2017) 1–51.
  * Nowak and May [1992] M. A. Nowak, R. M. May,  Evolutionary games and spatial chaos,  Nature 359 (1992) 826–829.
  * Nowak and May [1993] M. A. Nowak, R. M. May,  The spatial dilemmas of evolution,  International Journal of Bifurcation and Chaos 3 (1993) 35–78.
  * Szabó and Tőke [1998] G. Szabó, C. Tőke,  Evolutionary prisoner’s dilemma game on a square lattice,  Physical Review E 58 (1998) 69.
  * Ohtsuki et al. [2006] H. Ohtsuki, C. Hauert, E. Lieberman, M. A. Nowak,  A simple rule for the evolution of cooperation on graphs and social networks,  Nature 441 (2006) 502–505.
  * Lieberman et al. [2005] E. Lieberman, C. Hauert, M. A. Nowak,  Evolutionary dynamics on graphs,  Nature 433 (2005) 312–316.
  * Allen et al. [2017] B. Allen, G. Lippner, Y.-T. Chen, B. Fotouhi, N. Momeni, S.-T. Yau, M. A. Nowak,  Evolutionary dynamics on any population structure,  Nature 544 (2017) 227–230.
  * Chen et al. [2017] W. Chen, C. Gracia-Lázaro, Z. Li, L. Wang, Y. Moreno,  Evolutionary dynamics of $n$-person Hawk-Dove games,  Scientific Reports 7 (2017) 1–10.
  * Luo et al. [2021] Q. Luo, L. Liu, X. Chen,  Evolutionary dynamics of cooperation in the $n$-person stag hunt game,  Physica D: Nonlinear Phenomena 424 (2021) 132943.
  * Wang and Szolnoki [2022] C. Wang, A. Szolnoki,  Involution game with spatio-temporal heterogeneity of social resources,  Applied Mathematics and Computation 430 (2022) 127307.
  * Perc et al. [2013] M. Perc, J. Gómez-Gardenes, A. Szolnoki, L. M. Floría, Y. Moreno,  Evolutionary dynamics of group interactions on structured populations: a review,  Journal of the Royal Society Interface 10 (2013) 20120997.
  * Burgio et al. [2020] G. Burgio, J. T. Matamalas, S. Gómez, A. Arenas,  Evolution of cooperation in the presence of higher-order interactions: From networks to hypergraphs,  Entropy 22 (2020) 744.
  * Alvarez-Rodriguez et al. [2021] U. Alvarez-Rodriguez, F. Battiston, G. F. de Arruda, Y. Moreno, M. Perc, V. Latora,  Evolutionary dynamics of higher-order interactions in social networks,  Nature Human Behaviour 5 (2021) 586–595.
  * Pan et al. [2023] J. Pan, L. Zhang, W. Han, C. Huang,  Heterogeneous investment promotes cooperation in spatial public goods game on hypergraphs,  Physica A: Statistical Mechanics and its Applications 609 (2023) 128400\.
  * Duh et al. [2023] M. Duh, M. Gosak, M. Perc,  Unexpected paths to cooperation on tied hyperbolic networks,  Europhysics Letters 142 (2023) 62002.
  * Szabó and Hauert [2002] G. Szabó, C. Hauert,  Phase transitions and volunteering in spatial public goods games,  Physical Review Letters 89 (2002) 118101.
  * Hu et al. [2023] K. Hu, L. Shi, Y. Tao, M. Perc,  Cumulative advantage is a double-edge sword for cooperation,  Europhysics Letters 142 (2023) 21001.
  * Helbing and Yu [2009] D. Helbing, W. Yu,  The outbreak of cooperation among success-driven individuals under noisy conditions,  Proceedings of the National Academy of Sciences 106 (2009) 3680–3685.
  * Jiang et al. [2010] L.-L. Jiang, W.-X. Wang, Y.-C. Lai, B.-H. Wang,  Role of adaptive migration in promoting cooperation in spatial games,  Physical Review E 81 (2010) 036108.
  * Liu et al. [2012] Y. Liu, X. Chen, L. Zhang, F. Tao, L. Wang,  Does migration cost influence cooperation among success-driven individuals?,  Chaos, Solitons & Fractals 45 (2012) 1301–1308.
  * Buesser et al. [2013] P. Buesser, M. Tomassini, A. Antonioni,  Opportunistic migration in spatial evolutionary games,  Physical Review E 88 (2013) 042806.
  * Chen and Szolnoki [2016] X. Chen, A. Szolnoki,  Individual wealth-based selection supports cooperation in spatial public goods games,  Scientific Reports 6 (2016) 1–8.
  * Szolnoki and Chen [2016] A. Szolnoki, X. Chen,  Cooperation driven by success-driven group formation,  Physical Review E 94 (2016) 042311.
  * Szolnoki and Perc [2014] A. Szolnoki, M. Perc,  Coevolutionary success-driven multigames,  Europhysics Letters 108 (2014) 28004.
  * Tu [2018] J. Tu,  Contribution inequality in the spatial public goods game: Should the rich contribute more?,  Physica A: Statistical Mechanics and its Applications 496 (2018) 9–14.
  * Nowak et al. [2010] M. A. Nowak, C. E. Tarnita, T. Antal,  Evolutionary dynamics in structured populations,  Philosophical Transactions of the Royal Society B: Biological Sciences 365 (2010) 19–30.
  * Su et al. [2018] Q. Su, L. Wang, H. E. Stanley,  Understanding spatial public goods games on three-layer networks,  New Journal of Physics 20 (2018) 103030.
  * Su et al. [2019] Q. Su, A. Li, L. Wang, H. Eugene Stanley,  Spatial reciprocity in the evolution of cooperation,  Proceedings of the Royal Society B 286 (2019) 20190041.
  * Wang and Szolnoki [2022] C. Wang, A. Szolnoki,  A reversed form of public goods game: equivalence and difference,  New Journal of Physics 24 (2022) 123030.
  * Cheikbossian [2012] G. Cheikbossian,  The collective action problem: Within-group cooperation and between-group competition in a repeated rent-seeking game,  Games and Economic Behavior 74 (2012) 68–82.
  * Eckel et al. [2016] C. C. Eckel, E. Fatas, S. Godoy, R. K. Wilson,  Group-level selection increases cooperation in the public goods game,  PloS One 11 (2016) e0157840.
  * Wang et al. [2021] J. Wang, W. Xu, W. Chen, F. Yu, J. He,  Inter-group selection of strategy promotes cooperation in public goods game,  Physica A: Statistical Mechanics and its Applications 583 (2021) 126292\.
  * Su et al. [2016] Q. Su, A. Li, L. Zhou, L. Wang,  Interactive diversity promotes the evolution of cooperation in structured populations,  New Journal of Physics 18 (2016) 103007.
  * Su et al. [2017] Q. Su, A. Li, L. Wang,  Evolutionary dynamics under interactive diversity,  New Journal of Physics 19 (2017) 103023.
  * Wang et al. [2021] C. Wang, Q. Pan, X. Ju, M. He,  Public goods game with the interdependence of different cooperative strategies,  Chaos, Solitons & Fractals 146 (2021) 110871.
  * Wang and Huang [2022] C. Wang, C. Huang,  Between local and global strategy updating in public goods game,  Physica A: Statistical Mechanics and its Applications 606 (2022) 128097\.
  * Liu et al. [2010] Y. Liu, Z. Li, X. Chen, L. Wang,  Memory-based prisoner’s dilemma on square lattices,  Physica A: Statistical Mechanics and its Applications 389 (2010) 2390–2396.
  * Dong et al. [2019] Y. Dong, H. Xu, S. Fan,  Memory-based stag hunt game on regular lattices,  Physica A: Statistical Mechanics and its Applications 519 (2019) 247–255.
  * Shu et al. [2018] F. Shu, X. Liu, K. Fang, H. Chen,  Memory-based snowdrift game on a square lattice,  Physica A: Statistical Mechanics and its Applications 496 (2018) 15–26.