# Towards Explaining Adversarial Examples Phenomenon in Artificial Neural
Networks ††thanks: Preprint submitted to ICPR 2020.

1st Ramin Barati  Computer engineering Department
Amirkabir university of Technology
Tehran, Iran
<EMAIL_ADDRESS>2nd Reza Safabakhsh  Computer engineering Department
Amirkabir university of Technology
Tehran, Iran
<EMAIL_ADDRESS>3rd Mohammad Rahmati  Computer engineering Department
Amirkabir university of Technology
Tehran, Iran
<EMAIL_ADDRESS>
###### Abstract

In this paper, we study the adversarial examples existence and adversarial
training from the standpoint of convergence and provide evidence that
pointwise convergence in ANNs can explain these observations. The main
contribution of our proposal is that it relates the objective of the evasion
attacks and adversarial training with concepts already defined in learning
theory. Also, we extend and unify some of the other proposals in the
literature and provide alternative explanations on the observations made in
those proposals. Through different experiments, we demonstrate that the
framework is valuable in the study of the phenomenon and is applicable to
real-world problems.

###### Index Terms:

adversarial attack, robustness, artificial neural network, classifier,
learning theory, supervised learning, adversarial training

##  I Introduction

With the recent advancements in processing power and algorithms, machine
learning has become one of the more popular tools in many industries. It is
arguable that most of this attention is due to the combined effect of deep
learning and big training sets. The effectiveness of applying deep learning to
common machine learning tasks like image classification, object recognition,
natural language processing and speech recognition has fortified the position
of multi-layer perceptron (MLP) among the standard tools of solving machine
learning problems.

Considering these achievements, it would be natural to use artificial neural
networks (ANNs) in applications like autonomous driving where safety and
security is of importance. However, it was shown by Szegedy et al. that deep
networks are sensitive to adversarial examples[1]. Adversarial examples are
points close to a natural example in which the classifier output changes
drastically. In some cases, the difference between the adversarial example and
the natural example is not perceivable to humans.

In the near future, it is expected that many instances of MLPs are embedded in
the physical world. However, research has shown that adversarial samples can
be misused in physical world as well[2]. For example, the input signal can be
changed in a way that a speech recognition system order a wrong action[3] or
attack a self-driving car autonomous driving control[4]. As a result, many
concerns has been raised in regards to the use of ANNs in safety-critical
applications[5].

In practice, ANNs are considered black-box models. Although these models work
well in many tasks, they have proved to be difficult to explain and interpret.
Providing an explanation for the existence of adversarial examples is crucial
in minimizing their effect. But, a combination of difficulties in explaining
the inner workings of ANNs and visualization of high dimensional spaces has
made the problem hard to approach. Explaining the existence of adversarial
examples is still an open question and we refer the reader to [5] for a more
comprehensive study of research done on other aspects of this phenomenon.

In this paper, we describe a new perspective on the reasons behind the
existence of adversarial examples and discuss the consequences of our
approach. This paper is organized in five sections. After introduction, in
Section II, we provide a literature study of the relevant proposed
explanations of this phenomenon. Next, in Section III, we will describe the
proposed perspective and will continue by providing evidence in support of the
proposal in Section IV. We will finish the paper with a discussion of the
consequences and possible next steps in Section V.

##  II Related work

There have been multiple proposals for explaining the existence of adversarial
examples phenomenon. The first explanation was proposed in [1]. Szegedy et al.
viewed the adversarial examples as pockets in the input space which have a low
probability of being observed and correctly classified. This view has been
further described with an analogy between the real numbers as the natural
samples and the rational numbers as the adversarial examples. In other words,
There is a rational number in the neighborhood of every real number even
though the probability of observing a rational number in an interval of the
real line is zero. Nevertheless, as mentioned in [6], the argument does not
justify why would a classifier show such a behavior and it is only considered
as a side effect of nonlinearity of ANNs. This view is further undermined when
it was shown in [7] that adversarial examples are not isolated points and they
form dense regions in the input space.

Next, linearity hypothesis was proposed in [8]. In contrast with the low
probability pockets perspective, the linearity perspective relates the
adversarial phenomenon to a side effect of linear classifiers in high
dimensions. According to the linearity perspective, nonlinear classifiers like
ANNs are showing the phenomenon because they are trained with algorithms that
prefer linear models. Using this theory, Goodfellow et al. developed the Fast
Gradient Search Method (FGSM) which was able to easily produce adversarial
examples that would also transfer to other networks.

This attack is constructed as follows. Consider the inner product of the
weights $w$ and an adversarial example $\tilde{x}=x+\eta$.

$\langle w,\tilde{x}\rangle=\langle w,x\rangle+\langle w,\eta\rangle$ (1)

According to (1), if we choose $\eta=\mathrm{sign}(w)$, the change in the
result of the inner product would be maximized, constrained by $l^{\infty}$
norm. According to [8], since the sensory precision is fairly limited, the
input domain is discrete. But, on a computer, we have to represent these
quantities using floating-point precision numbers. If the input dimension is
large enough, the accumulated error of floating-point arithmetic in the output
could be considerable even if $\eta$ was actually small. Linearity means that
the trained classifiers show a similar behavior and would be fooled if we
choose $\eta$ to be the sign of the gradient of the loss function with respect
to $x$.

Linearity perspective explains adversarial transfer as the side effect of
converging to the optimal linear classifier. According to [6], this theory has
two predictions. First, that all linear classifiers show adversarial
phenomenon and second, that the effect would get worse as we increase the
dimensions of the input space. Nevertheless, Tanay and Griffin showed that
both of these predictions are wrong[6]. They constructed a linear classifier
which did not show the adversarial phenomenon and they showed that the effect
of the adversaries does not increase as the dimensions of the problem
increases.

The proposal of Ilyas et al. is also of interest in our discussion[9]. They
explain the existence of adversarial examples by attributing them to features
that are predictive, but nonrobust. A useful but nonrobust feature is a
feature that is highly predictive of the true label on the empirical
distribution of samples and labels, but if we add adversarial perturbations to
samples, it would not be as useful anymore. That is, for a distribution
$\mathcal{D}$, adversarial perturbations $\Delta(x)$ and a feature $f$ with a
positive correlation with the label $y$,

$\mathbb{E}_{\mathcal{D}}[y.f(x)]>\mathbb{E}_{\mathcal{D}}[y.f(x+\Delta(x))].$
(2)

The authors show that these features consistently exist in standard datasets
and tie the phenomenon they observed to a misalignment between the human-
specified notion of robustness and the inherent geometry of data. To support
their proposal, authors describe a process in which robust and nonrobust
features could be extracted from the training samples and showed that removal
of nonrobust features from the training samples would exchange accuracy with a
decrease in the adversarial phenomenon effects.

There are other proposal and studies as well. For example Shamir et al.
proposed that the phenomenon is a natural consequence of geometry of
$\mathbb{R}^{n}$ with Hamming metric[10], or Tanay and Griffin proposed that
the adversarial phenomenon can be reproduced by changing the angle of the
decision boundary manifold and the manifold that the natural samples reside
on[6]. These studies strongly suggest that the adversarial phenomenon is
somehow influenced by the training set and the way classifiers are defined and
trained. To the best of our knowledge, there is no approach to the problem in
which all these different observations and ideas unite. In the following
sections, we describe a new perspective on the cause of the phenomenon that
can potentially clarify the reasons behind the existence of adversarial
phenomenon and at the same time justify the observations made in the
literature.

##  III Proposed notion

In this section, we describe a setting in which the adversarial samples
phenomenon occur. Our proposal is based on the proposition that convergence of
ANNs to the true classifier is pointwise, not uniform. In the following, we
discuss the proposition and describe its role in producing the adversarial
examples phenomenon.

###  III-A Pointwise convergence

Consider the training set $S=\\{(0,-),(1,+)\\}$. This training set consists of
two samples from $[0,1]$ interval. Suppose that we want to find the maximum
margin classifier of $S$. To do so, we need to choose a set of features first.
Here, we use the Bernstein basis polynomials as features. The $n+1$ Bernstein
basis polynomials of degree $n$ are defined as

$b_{i,n}(x)={\binom{n}{i}}x^{i}\left(1-x\right)^{n-i},\quad i=0,\ldots,n.$ (3)

Then, we choose support vector classifiers $F_{n}$ as our hypothesis space,

$F_{n}(x;\alpha)=\alpha_{0}\sum_{i=0}^{n}b_{i,n}(0)b_{i,n}(x)+\alpha_{1}\sum_{i=0}^{n}b_{i,n}(1)b_{i,n}(x).$
(4)

It could be seen from (3) that $b_{i,n}(0)=\delta_{i,0}$ and
$b_{i,n}(1)=\delta_{i,n}$ where $\delta$ is the Kronecker delta function

$\delta_{i,j}={\begin{cases}0&{\text{if }}i\neq j,\\\ 1&{\text{if
}}i=j.\end{cases}}$

As a result, the maximum margin classifier $f_{n}(x)$ of degree $n$ would be
defined as

$f_{n}(x)=x^{n}-(1-x)^{n}.$ (5)

The above equation defines a sequence of maximum margin classifiers. For
example, $f_{1}(x)=2x-1$ is the linear maximum margin classifier of $S$. It
could be checked that $f_{n}$ is robust and does not show any adversarial
regions.

To illustrate the phenomenon, we pose the same classification task on $S$, but
use the first $n+1$ shifted Chebyshev basis polynomials as features instead.
Unlike Bernstein basis polynomials, Chebyshev basis polynomials do not
guarantee uniform convergence and reveal the error caused by pointwise
convergence. Shifted Chebyshev polynomials of the first kind are defined as

$T_{i}^{*}(x)=T_{i}(2x-1).$

and Chebyshev basis polynomials of the first kind are defined through the
recurrence relation

$\displaystyle T_{0}(x)$ $\displaystyle=1,$ $\displaystyle T_{1}(x)$
$\displaystyle=x,$ $\displaystyle T_{i+1}(x)$
$\displaystyle=2xT_{i}(x)-T_{i-1}(x).$

Similar to the case of Bernstein basis polynomials, the support vector
classifiers $G_{n}$ are defined by

$G_{n}(x;\alpha)=\alpha_{0}\sum_{i=0}^{n}T_{i}^{*}(0)T_{i}^{*}(x)+\alpha_{1}\sum_{i=0}^{n}T_{i}^{*}(1)T_{i}^{*}(x).$
(6)

It is known that $T_{i}^{*}(1)=1$ and $T_{i}^{*}(0)=(-1)^{i}$. As a result,
the maximum margin classifier $g_{n}(x)$ of degree $n$ would be defined as

$g_{n}(x)=\frac{2}{n+1}\sum_{i=0}^{[\frac{n-1}{2}]}T_{2i+1}^{*}(x).$ (7)

Figure 1 compares a few members of $f_{n}$ and $g_{n}$. It could be seen that
points exist that maximize or minimize the estimation error
$g_{n}(x)-f_{n}(x)$. The result is oscillations in the output of $g_{n}$. Due
to the sign of $g_{n}$ alternating between every root, predictions of $g_{n}$
agrees with $f_{n}$ at most half of the time. For a large $n$, the points that
do not agree become dense in the points that do. Nonetheless, for a small $n$,
they produce measurable regions in the input space.

Figure 1: A comparison of the Bernstein and Chebyshev maximum margin
classifiers of $S$. It could be seen that convergence of Chebyshev polynomials
is not uniform.

To put the phenomenon in learning theory terms, the hypothesis space
$\mathcal{P}$ of polynomials is nonuniform learnable. Hence, the hypothesis
space $\mathcal{P}_{n}$ of polynomials of degree $n$ has the uniform
convergence property. Bernstein polynomials of degree $n$ are testament to
this fact since they guarantee uniform convergence. A truncated Chebyshev
polynomial on the other hand only guarantees pointwise convergence. As a
result, even though as a hypothesis class it enjoys the uniform convergence
property, a Chebyshev polynomial falls under the notion of consistency as a
learning rule . We refer the reader to [11] for an overview of nonuniform
learnability and consistency.

We propose that the phenomenon occurs when a universally consistent learning
rule on a nonuniform learnable hypothesis class does not guarantee uniform
convergence. According to the definition of consistency, sample complexity of
a consistent learning rule depends on the generating distribution of the data
as well as the hypothesis. As a result, there could be distributions of data
that maximize or minimize the sample complexity of the learning rule for any
hypothesis. To make it precise, the probability $\delta$ of the classifier
making an error approaches zero, but the maximum possible error $\epsilon$ of
the classifier does not vanish as fast as $\delta$ as we draw more samples
from a sub-optimal training distribution.

To showcase the proposal, we keep the degree of the Chebyshev polynomial
constant and instead increase the number of training points. The training
points get labeled according to the nearest neighbour classifier of $S$. We
construct two training sets, $E_{m}$ and $C_{m}$. $E_{m}$ is constructed using
an equispaced grid of $m$ points, $C_{m}$ on the other hand use a Chebyshev
grid as samples.

The Chebyshev maximum margin classifier of degree 30 of $E_{m}$ and $C_{m}$ is
depicted in Figure 2 for a few choices of $m$. It could be seen in the figure
that in case of $E_{m}$, uniform convergence does not occur even for
relatively large $m$. In contrast, classifiers trained on $C_{m}$ converge in
the expected number of samples. As a matter of fact, Chebyshev polynomials
guarantee uniform convergence on a Chebyshev grid, and hence, do not suffer
from the phenomenon.

Figure 2: The trained Chebyshev maximum margin classifiers of degree 30 on
$C_{m}$ (top row) and $E_{m}$ (bottom row). We expect that the predictor
converge for $m>30$ on the equispaced grid. In practice though, it does not
converge even for relatively large $m$.

###  III-B Optimal points of a hypothesis class

In the previous section, we demonstrated that Chebyshev grid optimizes the
sample complexity of Chebyshev polynomials. These points seem to be reflecting
the position of the critical points of the estimation error
$g_{n}(x)-f_{n}(x)$. Thus, we define the optimal training points of a
hypothesis class in a way that reflects this observation.

Optimal training points: Let $\mathcal{X}$ be a domain set, let $\mathcal{H}$
be a hypothesis class and let $A$ be a universally consistent learning rule
with respect to $\mathcal{H}$. For every $X\subset\mathcal{X}$ and every
$h\in\mathcal{H}$, let

$\hat{h}_{X}=A(\\{(x,h(x))\,|\,x\in X\\}).$

The optimal training points of $A$ is a solution to the following problem plus
the boundary $\partial\mathcal{X}$ of $\mathcal{X}$, X∫_H∑_x∈X∥∇(^h_X(x)-
h(x))∥ dh X is feasible

The intuition behind the definition of the optimal points is that by
minimizing the empirical risk on the critical points of the error function, we
are effectively minimizing the maximum norm of the true error.

Computing the objective of (III-B) is not tractable for any practical purpose.
However, steps could be taken to calculate an approximation. First, we can
utilize the loss function $L$ as a surrogate for the objective of (III-B). We
know that the loss function is proportional to the magnitude of the error
function,

$L(h(x),h(x))\propto|\hat{h}_{X}(x)-h(x)|.$

Hence, points that maximize $L(h(x),h(x))$ would be critical points of
$\hat{h}_{X}(x)-h(x)$ and consequently would minimize
$\|\nabla(\hat{h}_{X}(x)-h(x))\|$. Of course, this is only a heuristic and
critical points of the error function could exist that are local minimas of
$|\hat{h}_{X}(x)-h(x)|$. Second, we approximate the integral in (III-B) by
restricting the hypothesis class. We call the solution to this approximation
the hard training points of the restricted class.

Hard training points of $H$: The hard training points of $A$ with respect to a
distribution $H$ on $\mathcal{H}$ is a solution to the following problem,
XE_H[∑_x∈XL(h(x),h(x))] X is feasible

Looking at (III-B), we can see that it is similar to the program that results
in adversarial examples. Indeed, if we reduce $H$ to a specific $h$, we would
find an approximation to (III-B) that closely resembles the objective of
evasion attacks.

Adversarial training points of $h$: The adversarial training points of $A$
with respect to $h\in\mathcal{H}$ is a solution to the following problem,
X∑_x∈XL(h(x),h(x)) X is feasible

There are two differences between the proposed definition of the adversarial
training points of $h$ and the standard definition used in adversarial
training. First, our definition does not use the true labels in computing the
loss. While the standard definition may lead to stronger adversaries, our
proposal makes the independence of the phenomenon from the true predictor
explicit. Second, Our definition does not require the points to be close to
natural samples.

##  IV Experiments

In section III, we introduced a framework in which adversarial examples are
produced. Here, we back the propositions by providing empirical evidence in
their support. In our first experiment, we show that a Chebyshev grid does
optimize (III-B). Then, we demonstrate how (III-B) could be used to
adversarially train a MLP.

###  IV-A Hard training points of Chebyshev polynomials

To exhibit that the definition of hard training points is consistent with the
observations made in section III about Chebyshev polynomials, we now derive
and visualize the objective of (III-B) for a special family of polynomials in
a maximum margin classification setting.

It is known that the hinge loss function is the loss of the maximum margin
classifier. Thus, we are looking for hard training points of a shifted
Chebyshev polynomial of degree $n$ with respect to the hinge loss function.
XE_H[∑_x∈Xmax(0,1-h(x)^2)] 0 ≤x_1¡…¡x_m ≤1

The $\max$ operator in the objective of (IV-A) would make the objective hard
to manipulate. As a result, we restrict ourselves to a family of polynomials
with range in the $[-1,1]$ interval. Shifted Chebyshev basis polynomials are
examples of such polynomials. Using shifted Chebyshev basis polynomials as
$H$, the problem will change to a form that is easier to analyse and compute.
X∑_x∈X1n+1∑_i=0^n 1-T_i^*(x)^2 0 ≤x_1¡…¡x_m ≤1

From (IV-A) it could be deduced that the local minimas of
$\sum_{i=0}^{n}T_{i}^{*}(x)^{2}$ solve (IV-A). Next, it could be argued that
the optimal solution is close to the roots of $T_{n}^{*}(x)$, which are
exactly the points of a Chebyshev grid.

Figure 3 compares the position of Chebyshev nodes with
$\frac{1}{n+1}\sum_{i=0}^{n}1-T_{i}^{*}(x)^{2}$. It could be checked in the
figure that Chebyshev nodes are very close to the optimal points. The two
nodes on the two end of $[0,1]$ are also critical points of the error function
because they belong to the boundary of the interval. In our experiments,
removal of the two boundary points has a considerable negative effect on the
trained classifier. This shows that we should consider the boundary
$\partial\mathcal{X}$ of the domain set as well.

Figure 3: The expected value of loss on $H$ in (IV-A) for $n=10$. The points
mark the value of the expectation at each node of a Chebyshev grid.

###  IV-B Optimal training points in MLPs

We showed that the proposed notion can explain the phenomenon in case of
polynomials. Here, we analyze MLPs through our framework. To this end, we
sample a random MLP classifier and visualize the adversarial and hard training
points objectives with respect to that network.

The domain set of the network is $[-1,1]^{2}$. The dimensions of the network
are $2\times 100\times 100\times 100\times 5$. The weights of the network are
sampled independently from a standard normal distribution. The activation
function of the hidden layers is $\tanh$ and the output of the network is
activated through softmax. We have chosen the cross entropy loss as our
surrogate for the error function.

The objective of (III-B) is straightforward to compute. To compute (III-B), we
need to choose a distribution over the hypothesis class of MLPs. Following the
case of Chebyshev polynomials, we use each neuron in the feature layer of the
MLP as the hypothesis set. In other words, we create a copy of the network,
but replace the last layer with an identity transform and compute (III-B) for
the new network. For this particular example, the size of the resulting
network would be $2\times 100\times 100\times 100\times 100$.

Figure 4 shows the results. By comparing the output of the original network
with the adversarial objective of the network, we can see that the adversarial
objective is sensitive to the position of the decision boundary. But, we can
also identify some regions of the objective that get lighter or darker without
respecting the position of the decision boundary. According to our proposal,
if we sample the training points proportional to the adversarial objective,
the sample complexity of finding this exact hypothesis would be minimized.

Contrasting the adversarial objective and the hard objective in Figure 4
reveals that a hypothesis agnostic set of points would not be as helpful as
for the case of polynomials. In other words, while we could reach a definitive
set of optimal points for polynomials, the optimal points of MLPs are not
unique even up to the layers of a single network. As a result, even though the
proposed notion can explain the existence and abundance of adversarial
examples in MLPs, it cannot further explain their transfer between different
architectures of MLPs. Specially, we could marginalize out the hypothesis
space from the sample complexity in polynomials, but the same could not be
said about MLPs.

Figure 4: A comparison between the output, the adversarial objective and hard
objective of a random network. The adversarial objective is representative of
the decision boundary. On the other hand, the hard objective suggests a very
different distribution of optimal points.

###  IV-C Adversarial training of a MLP

In this section, we show that sampling according to the hard objective
described in the previous experiment does indeed optimize the sample
complexity of the robust hypothesis. To do so, we adversarially train a MLP on
a subset of MNIST using the hard training points of the feature layer. Next,
we train a MLP with the same architecture but this time use the whole of MNIST
as the training set with no adversarial training. Then, we compare the
performance of the networks in both adversarial and nonadversarial settings.
To show that the results are consistent across multiple real-world datasets,
we conduct the same experiments for Fashion-MNIST[12] and report the results.

Here we use $784\times 512\times 512\times 10$ networks with ReLU activation.
The robust network is trained using a third of the full datasets. We have
trained each network for 10 epochs. In each step of the adversarial training,
we find the hard training points close to the batch with regards to the
weights of the feature layer in that step. To be more precise, we use the
current batch as the initial points and take a few gradient ascend steps in
the direction of (III-B). Then, we train the network on this new batch of
samples as well.

Figure 5 shows the accuracy of all networks on the test set as a function of
the magnitude of the adversarial perturbation. A $L^{2}$ normalized one-step
gradient based attack is used for this experiment.

(a) MNIST

(b) Fashion-MNIST

Figure 5: The accuracy of MLPs compared to the magnitude of the adversarial
perturbation. The data shows that the adversarial examples of the
adversarially trained network does transfer to the reference network, but not
vice versa.

It could be seen that the adversarially trained network gets a comparable
performance with the reference network on the natural samples set even though
it is trained on a third of the samples. The adversarially trained network is
also more robust compared to the reference network if we compare their
performance in the direct and transfer attack settings. The results show that
adversarial training using (III-B) would make the trained MLP more robust
against adversarial attacks by attacking random hypothesises that use the same
feature layer instead of attacking the hypothesis itself.

##  V Conclusion and future work

In this paper, we introduced a framework to explain the adversarial examples
phenomenon. We defined the adversarial examples as the critical points of the
error function. We showed that this principle can fully explain adversarial
examples existence, training and transfer for pointwise converging
polynomials.

Unfortunately, this approach does not enjoy the same success in case of MLPs.
We could demonstrate that MLPs do follow the principle in case of existence
and training, but the same principle proved to be insufficient in explaining
transferable adversarial examples. Nevertheless, the results suggest that the
phenomenon would disappear if the learning process guarantee uniform
convergence of MLPs.

With regards to other proposals in literature, our definition is more aligned
with the low probability pockets perspective. We constructed the first
instance of a classifier with adversarial examples that are dense in the input
domain. While the low probability pockets perspective relates the phenomenon
to highly nonlinear nature of ANNs, we have shown that the phenomenon does
occur for low degree polynomials as well, which would produce adversarial
regions in turn.

With respect to the linearity hypothesis, our analysis does not confirm the
view that the phenomenon is a result of representing the input by floating-
point numbers. We also introduced an alternative explanation for the success
of gradient based attacks that maximize the training loss function. While our
proposal does not confirm that adversarial examples transfer is due to
converging to the optimal linear classifier, we share the idea that the
transfer is rooted in properties of the optimal classifier.

The proposed notion is similar to the nonrobust features perspective in that
we also relate the phenomenon to a form of optimal training points. Our
proposal can accommodate a notion of pointwise converging features similar to
nonrobust features as well. On the other hand, we do not support the idea that
adversarial examples transfer in MLPs is caused by patterns in the input that
are imperceptible to humans. In case of MLPs, we did not observe any
hypothesis independent relation between the distribution of the optimal points
and the domain set in our experiments.

Our framework is heavily inspired by analogous concepts in approximation
theory. For now, we have analysed both the training points and the hypothesis
class by introducing a definition for optimal training points and revealing
the role of convergence properties of the hypothesis class in the phenomenon.
If approximation theory is any indicator of the path forward, our next step
should be to analyse the classifier independently from the training points. In
other words, we guess that there is something special about generalizing
decision boundaries that is causing the transfer of the adversarial examples
in MLPs. In future work, we would focus on finding extensions to the framework
that would help explain the transfer of adversarial examples in MLPs.

## References

  * [1] C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. Goodfellow, and R. Fergus, “Intriguing properties of neural networks,” _arXiv preprint arXiv:1312.6199_ , 2013.
  * [2] A. Kurakin, I. J. Goodfellow, and S. Bengio, “Adversarial examples in the physical world,” _CoRR_ , vol. abs/1607.02533, 2016. [Online]. Available: http://arxiv.org/abs/1607.02533
  * [3] L. Schönherr, K. Kohls, S. Zeiler, T. Holz, and D. Kolossa, “Adversarial attacks against automatic speech recognition systems via psychoacoustic hiding,” in _Network and Distributed System Security Symposium (NDSS)_ , 2019\.
  * [4] A. Boloor, X. He, C. D. Gill, Y. Vorobeychik, and X. Zhang, “Simple physical adversarial examples against end-to-end autonomous driving models,” _CoRR_ , vol. abs/1903.05157, 2019. [Online]. Available: http://arxiv.org/abs/1903.05157
  * [5] X. Yuan, P. He, Q. Zhu, and X. Li, “Adversarial examples: Attacks and defenses for deep learning,” _IEEE Transactions on Neural Networks and Learning Systems_ , p. 1–20, 2019. [Online]. Available: http://dx.doi.org/10.1109/TNNLS.2018.2886017
  * [6] T. Tanay and L. Griffin, “A boundary tilting persepective on the phenomenon of adversarial examples,” _arXiv preprint arXiv:1608.07690_ , 2016.
  * [7] P. Tabacof and E. Valle, “Exploring the space of adversarial images,” _2016 International Joint Conference on Neural Networks (IJCNN)_ , Jul 2016\. [Online]. Available: http://dx.doi.org/10.1109/IJCNN.2016.7727230
  * [8] I. J. Goodfellow, J. Shlens, and C. Szegedy, “Explaining and harnessing adversarial examples,” _arXiv preprint arXiv:1412.6572_ , 2014.
  * [9] A. Ilyas, S. Santurkar, D. Tsipras, L. Engstrom, B. Tran, and A. Madry, “Adversarial examples are not bugs, they are features,” _arXiv preprint arXiv:1905.02175_ , 2019.
  * [10] A. Shamir, I. Safran, E. Ronen, and O. Dunkelman, “A simple explanation for the existence of adversarial examples with small hamming distance,” _CoRR_ , vol. abs/1901.10861, 2019. [Online]. Available: http://arxiv.org/abs/1901.10861
  * [11] S. Shalev-Shwartz and S. Ben-David, _Understanding machine learning: From theory to algorithms_. Cambridge university press, 2014.
  * [12] H. Xiao, K. Rasul, and R. Vollgraf. (2017) Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms.


# Aggregation and disaggregation processes in clusters of particles: simple
numerical and theoretical insights of the competition in 2D geometries

Louis-Vincent Bouthier  Corresponding author
<EMAIL_ADDRESS>Romain Castellani  Groupe CFL, CEMEF, Mines
Paris, PSL Research University, 1 Rue Claude Daunesse, 06904 Sophia Antipolis,
France  Sébastien Manneville  ENSL, CNRS, Laboratoire de physique, F-69342
Lyon, France  Institut Universitaire de France (IUF)  Arnaud Poulesquen  CEA,
DES, ISEC, DE2D, SEAD, LCBC, Université de Montpellier, Marcoule, France  Rudy
Valette  Elie Hachem  Groupe CFL, CEMEF, Mines Paris, PSL Research University,
1 Rue Claude Daunesse, 06904 Sophia Antipolis, France

###### Abstract

Aggregation and disaggregation of clusters of attractive particles under flow
are studied from numerical and theoretical points of view. Two-dimensional
molecular dynamics simulations of both Couette and Poiseuille flows highlight
the growth of the average steady-state cluster size as a power law of the
adhesion number, a dimensionless number that quantifies the ratio of
attractive forces to shear stress. Such a power-law scaling results from the
competition between aggregation and disaggregation processes, as already
reported in the literature. Here, we rationalize this behavior through a model
based on an energy function, which minimization yields the power-law exponent
in terms of the cluster fractal dimension, in good agreement with the present
simulations and with previous works.

Colloids, Clusters, Aggregation, Disaggregation, 2D Poiseuille flow, 2D
Couette flow

##  I Introduction

Colloidal suspensions composed of attractive particles, that aggregate into
clusters when dispersed into a liquid medium, have raised great interest for a
number of decades due to their wide range of applications, from paints and
coatings to food products or construction materials [1, 2, 3, 4, 5, 6, 7, 8].
The large variety of attractive forces driving cluster aggregation, including
van der Waals interactions, depletion or capillary forces, combined with
stabilizing repulsive forces, such as electrostatic or steric interactions,
lead to a complex phase diagram involving equilibrium liquid or phase-
separated regimes as well as out-of-equilibrium gel or glassy phases [9, 2, 3,
10]. Beyond their phase behaviour at rest, understanding how colloidal
suspensions respond to deformation and flow is pivotal for their processing
and their practical use. Thus, a huge amount of work has been devoted to the
rheology of colloidal suspensions and their structure under external
mechanical solicitations [2, 3, 1, 11, 12]. In particular, one key question
concerns the evolution of particle clusters under shear.

Colloidal clusters generically show a ramified structure characterized by a
fractal dimension $D$ linking the number of particles $n(l)$ in the cluster to
its size $l$ through $n(l)\propto l^{D}$ [13, 14, 15]. Other structural
parameters may be identified, such as the chemical dimension related to the
number of bonds in the shortest path in a cluster [16], the connectivity or
the coordination number [17]. Since the seminal work of Weitz et al. [18, 19]
on diffusion-limited aggregation, a large number of experimental studies have
characterized the microscopic mechanisms underpinning cluster formation at
rest, either from direct visualisation [20, 21], or from scattering techniques
[22]. In cases where colloidal aggregation leads to the formation of a space-
spanning cluster network, i.e. to a colloidal gel, these measurements have
provided key information on both the maximum cluster size $\ell$ and their
fractal dimension $D$ [23]. When shear is applied to a flocculated suspension,
aggregates break up and form a suspension of isolated clusters. The steady-
state cluster size distribution is a key observable and the dependence of
$\ell$ and $D$ on the shear stress applied to the suspension has been probed
by combining light, x-ray, and neutron scattering and rheometry [24, 25, 26,
27]. Typical results are that (i) $\ell$ decreases as a power-law of the shear
stress $\sigma$, $\ell\propto\sigma^{-m}$, with an exponent $m$ ranging from
0.2 to 0.5 [28, 29], and (ii) $D$ increases from 1.8–2.2 depending on the
specific aggregation mechanism to 2.4–2.7 under shear, a phenomenon known as
“shear-induced cluster densification” [24, 30].

Furthermore, thanks to the growth of computational capabilities over the past
decade, large-scale numerical simulations have been developed to model
suspensions, both at rest and under shear, with up to several millions of
colloidal particles [31, 32], sometimes taking into account the solvent
through hydrodynamic interactions [33, 34, 35, 36, 37]. The model suspensions
rely on various idealized two-body interaction potentials [38], which in some
cases, include bending stiffness [39, 40, 17], or solid friction [41, 42, 43].
One of the main advantages of numerical simulations is that they yield a
complete picture of the cluster microstructure. Another modelling approach
referred to as “population balance model” (PBM), relies on kinetic equations
for “classes” of clusters of given sizes, which mimic the competition between
aggregation and disaggregation processes through the use of various kernels
[44, 45, 46, 47]. Although some ingredients of the kernels can be
theoretically predicted thanks to Smoluchowksi equations, in practice, PBM
often require a significant amount of empirical tuning of the kernel
parameters [15, 48].

Strong theoretical predictions have also been derived from analytical models
that relate microscopic parameters, such as the interparticle bond stiffness,
to macroscopic quantities, such as the elastic modulus $G^{\prime}$ at rest,
the yield strain $\gamma_{y}$ or the dynamic viscosity $\eta$ under shear [49,
12, 50]. In particular, one popular model by Wessel & Ball [11] predicts that
$\ell\propto\sigma^{-1/3}$ from a force balance under the assumption that the
clusters behave hydrodynamically like compact spheres. While this model
provides a correct estimate of the exponent $m$ for various experiments [26],
it does not account for the range of exponents reported in the literature [51,
52, 27], most probably due to the stringent assumption of hydrodynamically
compact clusters. More refined models have considered soft, highly deformable
clusters leading to $m=1/2$ rather than $m=1/3$ for rigid aggregates [53].
Yet, to the best of our knowledge, there is no general theoretical framework
that may describe the behaviour of colloidal clusters under an external
stress.

The goal of the present paper is to propose very general theoretical arguments
to predict the steady-state size $\ell$ of colloidal aggregates submitted to
an external solicitation. Starting from dimensional analysis, we show in
Section II that the competition between cluster aggregation and disaggregation
can be captured through the minimization of an energy function. In Section
III, we then provide evidence for the existence of an energy minimum in simple
numerical simulations. We proceed to detail an analytical model in Section IV,
which yields a power-law for $\ell$ as a function of the adhesion number, a
dimensionless number that quantifies the ratio of attractive forces to shear
stress. Finally, this model is compared to the coagulation-fragmentation
approach and to previous experimental and numerical findings in Section V.
Conclusions and open questions are drawn in Section VI.

##  II General framework

Before describing the general approach based on a grand-canonical free energy,
we start with simple dimensional arguments to justify that the competition
between cluster aggregation led by attractive forces and disaggregation driven
by shear forces depends on a single dimensionless group, namely the adhesion
number $\mathrm{Ad}$, once the attraction range is fixed. As discussed in
Refs. [41, 54], attractive forces may be estimated by $Ua/\delta^{2}$, with
$U$ the depth of the interaction potential between two particles, $a$ the
particle radius, and $\delta$ the range of interaction, which can be taken as
the center-to-center distance between two particles at equilibrium or as the
width of the potential well. Disaggregating forces, on the other hand, may be
estimated by $\sigma a^{2}$, with $\sigma$ the external stress exerted on the
clusters. One chooses here, as a distinction from the literature, to consider
the stress rather than the shear rate because (i) whenever the shear rate is
involved in similar definitions, it appears as multiplied by the solvent
viscosity, therefore as a shear stress, and (ii) one is convinced that stress
drives the disaggregation, like in plasticity or fracture [55, 56, 57, 58],
rather than a kinematic quantity. Following Ref. [54], the adhesion number is
defined as the ratio of attractive forces to disaggregating forces:

$\mathrm{Ad}=\frac{U}{\sigma a\delta^{2}}\,.$ (1)

The condition $\mathrm{Ad}\ll 1$ implies that hydrodynamic forces dominate,
while $\mathrm{Ad}\gg 1$ indicates that attractive forces are predominant. It
is important to note that the choice of this adhesion number contains some
degree of arbitrariness. Indeed, other similar dimensionless groups may be
built more generally by replacing $a\delta^{2}$ by
$a^{\alpha}\delta^{3-\alpha}$, with $\alpha\in\left[0,3\right]$, in Eq. (1).
The precise choice depends on whether one considers energies ($\alpha=0$),
forces ($\alpha=1$), or stiffnesses ($\alpha=2$), i.e., energies per unit
surface, or energies per unit volume ($\alpha=3$). Finally, according to the
Vaschy-Buckingham theorem [59, 60, 61], any characteristic length that depends
on $U$, $\sigma$, $a$, and $\delta$, such as the maximum cluster size $\ell$,
may be expressed as $\ell/a=\mathcal{F}\left(\mathrm{Ad},\delta/a\right)$.
Therefore, a combination of both $\mathrm{Ad}$ and $\delta/a$ is expected.
Also note that the adhesion number simply corresponds to the inverse of the
“Mason number,” $\mathrm{Mn}=\sigma a^{2}\delta/U$, a dimensionless group
popular in the rheology community that quantifies the ratio of shearing forces
to attractive forces [36, 62, 63].

In order to describe the competition between aggregation and disaggregation,
we consider a grand-canonical ensemble with a population of clusters of mass
$k\in\mathbb{N}^{*}$ associated to a number of clusters $n_{k}$. Each cluster
of mass $k$ has an associated energy $\mathcal{E}\left(k\right)$, which we
seek to determine. The number of primary particles is not fixed and is
related, for each cluster mass $k$, to a chemical potential $\alpha_{k}/\beta$
with $\beta=1/(k_{B}T)$, $k_{B}$ the Boltzmann constant, and $T$ the
temperature. Note that the use of a canonical ensemble may seem more
appropriate for a problem with a fixed number of particles. However, the
calculation of the canonical partition function leads to the use of the
complete Bell polynomials, where each variable is
$\mathrm{e}^{-\beta\mathcal{E}\left(k\right)}$. Inverting the relation is not
straightforward, so that computing the distribution of clusters of mass $k$ is
very cumbersome within a canonical framework. Here, thanks to the grand-
canonical formulation, $\mathcal{E}\left(k\right)$, as well as the energy
$n_{k}\mathcal{E}\left(k\right)$ associated to all clusters of mass $k$, may
be computed rather easily. The grand canonical partition function reads:

$\displaystyle\Xi$
$\displaystyle=\sum_{\left(n_{k}\right)_{k\in\mathbb{N}^{*}}\in\mathbb{N}^{\mathbb{N}^{*}}}\exp\left(-\sum_{k\in\mathbb{N}^{*}}\left(\beta
n_{k}\mathcal{E}\left(k\right)-n_{k}k\alpha_{k}\right)\right)$ (2)
$\displaystyle=\prod_{k\in\mathbb{N}^{*}}\left(1-\mathrm{e}^{-\beta\mathcal{E}\left(k\right)+k\alpha_{k}}\right)^{-1}.$
(3)

The convergence of the series in Eq. (2) is guaranteed if
$\alpha_{k}k<\beta\mathcal{E}\left(k\right)$ for all $k\in\mathbb{N}^{*}$.
Moreover, it follows from Eq. (2) that the distribution of each level
population $n_{k}$ is a geometric distribution with a parameter
$\mathrm{e}^{-\beta\mathcal{E}\left(k\right)+k\alpha_{k}}$. Since the average
number of particles $\left\langle N\right\rangle$ is related to the average
population $\left(\left\langle n_{k}\right\rangle\right)_{k\in\mathbb{N}^{*}}$
of clusters of mass $k$ through $\left\langle
N\right\rangle=\sum_{k\in\mathbb{N}^{*}}k\left\langle n_{k}\right\rangle$, one
may compute the average number of clusters of mass $k$ and the standard
deviation through:

$\displaystyle\left\langle n_{k}\right\rangle$
$\displaystyle=-\frac{1}{k}\frac{\partial\ln\Xi}{\partial\alpha_{k}}=\left(\mathrm{e}^{\beta\mathcal{E}\left(k\right)-\alpha_{k}k}-1\right)^{-1},$
(4) $\displaystyle\Delta n_{k}$
$\displaystyle=\frac{1}{k}\sqrt{\frac{\partial^{2}\ln\Xi}{\partial\alpha_{k}^{2}}}=\sqrt{\left\langle
n_{k}\right\rangle\left(1+\left\langle n_{k}\right\rangle\right)}$ (5)
$\displaystyle=\frac{1}{2}\,\mathrm{csch}\left(\frac{\beta\mathcal{E}\left(k\right)-\alpha_{k}k}{2}\right),$
(6)

which boils down to the Bose-Einstein statistics.

In practice, in an experiment or a simulation, a initial sample of
$\left(n_{k}\right)_{k\in\mathbb{N}^{*}}$ is chosen, which evolves in time
according to external conditions. It is clear from the above expressions that
the $\left(n_{k}\right)_{k\in\mathbb{N}^{*}}$ should end up being centered
around the averages $\left(\left\langle
n_{k}\right\rangle\right)_{k\in\mathbb{N}^{*}}$ with a dispersion
$\left(\Delta n_{k}\right)_{k\in\mathbb{N}^{*}}$. Therefore, a higher number
of clusters of mass $k$ also brings a higher dispersion around this value.
Since numerical simulations yield a probability distribution function
$f\left(k\right)$ for the number of cluster of mass $k$ , which is assumed to
be equal to the average distribution $\left\langle n_{k}\right\rangle$, i.e.,
$f\left(k\right)=\left\langle n_{k}\right\rangle$, one can find the average
potential energy of the clusters and the standard deviation respectively
through:

$\displaystyle\beta\mathcal{E}\left(k\right)-\alpha_{k}k$
$\displaystyle=\ln\left(1+\frac{1}{f\left(k\right)}\right)$ (7)
$\displaystyle\Delta\left(\beta\mathcal{E}\left(k\right)-\alpha_{k}k\right)$
$\displaystyle=\left(f\left(k\right)\left(1+f\left(k\right)\right)\right)^{-\frac{1}{2}}$
(8)
$\displaystyle=2\,\sinh\left(\frac{\beta\mathcal{E}\left(k\right)-\alpha_{k}k}{2}\right)$
(9) $\displaystyle>\beta\mathcal{E}\left(k\right)-\alpha_{k}k\,.$ (10)

##  III Numerical simulation

###  III.1 Numerical scheme and analysis

To get some insight of the competition between aggregation and disaggregation
processes, we turn to simple molecular dynamics simulations of two-dimensional
Couette and Poiseuille flows based on the LAMMPS library [64]. We use reduced
units, where the unit size is given by the particle radius. The simulation box
is of size $\mathcal{L}\times h$ with periodic boundary conditions along the
$x$ direction and solid boundaries along the $y$ direction located at $y=0$
and $y=h$. The length of the box $\mathcal{L}$ is fixed to $\mathcal{L}=620$,
and the width $h$ is either 103 or 206, much larger than the particle size.
Each boundary is constituted of one layer of particles of unit size which
positions are fixed and that interact with bulk particles through a repulsive
Yukawa potential $w_{s}(r)=100\mathrm{e}^{-r}/r$, with $r$ the distance
between two particles and a cut-off distance $r=5$ to save computation time.
The interaction potential between two particles in the bulk is a classical
12-6 Lennard-Jones potential $w(r)=4\left(r^{-12}-r^{-6}\right)$ with a cut-
off distance $r=40$, again to save computation time.

First, the system is initialised with a particle surface fraction of
$\phi=0.12$. Using a canonical formulation, the temperature is fixed at
$T=0.01$. The time step is taken as $\Delta t=0.01$, which is small enough to
account realistically for temporal variations while keeping the computation
time reasonably low. Particles are distributed over a square lattice, and
their initial velocities are chosen according to a Maxwell distribution for
the given temperature. The system is then let to evolve for a duration $5,000$
time units in order to create the initial cluster distribution. More
precisely, the equation of motion for each particle $i\in\left\llbracket
1,N\right\rrbracket$, with $N$ the total number of particle, is

$\displaystyle\frac{\mathrm{d}^{2}\bm{r}_{i}}{\mathrm{d}t^{2}}$
$\displaystyle=-\frac{\partial\mathcal{W}}{\partial\bm{r}_{i}}\left(\left(\bm{r}_{j}\right)_{j\in\left\llbracket
1,N\right\rrbracket}\right)+\bm{B}_{i}$ (11)
$\displaystyle\mathcal{W}\left(\left(\bm{r}_{j}\right)_{j\in\left\llbracket
1,N\right\rrbracket}\right)$
$\displaystyle=\sum_{k=1}^{N}\sum_{j=k+1}^{N}w\left(\left|\bm{r}_{k}-\bm{r}_{j}\right|\right)$
(12) $\displaystyle w\left(r\right)$
$\displaystyle=4\left(r^{-12}-r^{-6}\right)\,,$ (13)

with $\bm{r}_{i}$ the position of particle $i$ and $\bm{B}_{i}$ a Brownian
white noise for particle $i$.

Figure 1: Simulated particle distributions for $V=10$, $h=206$, and $C=1$
after the preparation step prior to shearing (a) and at after shearing in the
Couette geometry (b). Rendering using _Ovito_ [65]. (c) Cluster detection and
various estimates of the cluster size as defined in the text for the top right
particle distribution.

In a second step, a drag force
$\bm{F}_{i}=C\left(\bm{u}-\frac{\mathrm{d}\bm{r}_{i}}{\mathrm{d}t}\right)$ is
applied on each particle $i\in\left\llbracket 1,N\right\rrbracket$, with $C=1$
the drag coefficient 111Additional simulations were carried out with various
values of $C$ and the main results remained similar up to a rescaling of the
time scale., $\bm{u}=Vy/h\bm{e}_{x}$ for the 2D Couette flow or
$\bm{u}=4V\left(1-y/h\right)y/h\bm{e}_{x}$ for the 2D Poiseuille flow, where
$V$ is the maximum flow velocity. Here, a micro-canonical formulation is used
and the time step is set to $\Delta t=0.001$. The simulation is run for a
duration of $1,000$ time units. The equation of motion for each particle $i$
thus reads:

$\displaystyle\frac{\mathrm{d}^{2}\bm{r}_{i}}{\mathrm{d}t^{2}}$
$\displaystyle=-\frac{\partial\mathcal{W}}{\partial\bm{r}_{i}}\left(\left(\bm{r}_{j}\right)_{j\in\left\llbracket
1,N\right\rrbracket}\right)+F\left(\bm{r}_{i}\right)\bm{e}_{x}$ (14)
$\displaystyle F\left(\bm{r}_{i}\right)$
$\displaystyle=C\left(-\frac{\mathrm{d}\bm{r}_{i}}{\mathrm{d}t}\cdot\bm{e}_{x}+\right.$
(15)
$\displaystyle\left.V\frac{\bm{r}_{i}\cdot\bm{e}_{y}}{h}\begin{cases}1&\text{2D
Couette flow}\\\ 4\left(1-\frac{\bm{r}_{i}\cdot\bm{e}_{y}}{h}\right)&\text{2D
Poiseuille flow}\end{cases}\right)\,.$ (16)

The parameters investigated in the present work are
$\left(V,h\right)\in\left\\{0.1,0.3,1,3,10,30,100\right\\}\times\left\\{103,206\right\\}$.
The numerical scheme used here is a Velocity Verlet algorithm [67]. The
computations output are (i) the position of the particles at each time step,
and (ii) the clusters to which the particles belong based on a connected-
component algorithm [68] with a distance threshold of 1.4, consistently with
the literature [69, 39, 70]. Moreover, we checked that for thresholds ranging
from 1.1 to 2.0, the distribution of neighbours does not change, so that the
results are not sensitive to the specific choice of threshold. The size of
cluster number $I$ is quantified according to the following estimates:

  * •

the radius of gyration $R_{g}^{I}$ given by

$\displaystyle R_{g}^{I}$
$\displaystyle=\sqrt{\frac{1}{\left|P\left(I\right)\right|}\sum_{k\in
P\left(I\right)}\left\lVert\bm{r}_{k}-\overline{\bm{r}}\right\rVert^{2}},$
(17) $\displaystyle\overline{\bm{r}}$
$\displaystyle=\frac{1}{\left|P\left(I\right)\right|}\sum_{k\in
P\left(I\right)}\bm{r}_{k},$ (18)

with $\bm{r}_{k}$ the position of the particles in the cluster,
$P\left(I\right)$ the set of particles in cluster $I$ and
$\left|P\left(I\right)\right|$ the number of particles in cluster $I$.

  * •

the Feret radii in the $x$ and $y$ directions given by $\left(\max_{k\in
P\left(I\right)}\bm{r}_{k}\cdot\bm{e}\right)/2-\left(\min_{k\in
P\left(I\right)}\bm{r}_{k}\cdot\bm{e}\right)/2$ with $\bm{e}$ the unit vector
in the $x$ and $y$ directions respectively.

  * •

the half maximum chord length given by $\max_{\left(k,l\right)\in
P\left(I\right)^{2}}\left\lVert\bm{r}_{k}-\bm{r}_{l}\right\rVert/2$.

In order to infer statistical estimations, each size distribution is further
weighted by the number of particles in each cluster. Such weighting is needed
because the number of clusters is not constant. Therefore, because the total
number of particles is constant, weighting by the mass of each cluster allows
one to recover the number of particles when integrating over the whole
distribution. Finally, thanks to the reduced units, $U/a\delta^{2}=1$ and the
shear stress is $\sigma=CV/h=V/h$ here due to $C=1$ so that the adhesion
number simply reads $\mathrm{Ad}=h/V$.

###  III.2 Simulation results

Figure 1 shows typical particle distributions computed after the preparation
step prior to shearing (a), and after application of shear in the Couette
geometry (b). It appears clearly that the system starts from a space-spanning
ramified structure and evolves toward dense, isolated clusters under shear.
Moreover, as shown in Fig. 1(c), the different estimates for the cluster size
yield consistent values. In the following, for the sake of simplicity, we
shall focus only on the weighted average of the radius of gyration to estimate
the cluster size $\ell$.

Figure 2: Energy functional $\beta\mathcal{E}\left(k\right)-\alpha_{k}k$
computed from from Eq. (7) as a function of the number $k$ of particles in a
cluster for (left) Couette flow and (right) Poiseuille flow. Colors correspond
to the adhesion number $\mathrm{Ad}$ as indicated in the legend. Results
obtained from simulations performed with $h=103$ or $h=206$.

Following the general framework introduced in Sect. II, we compute the
distribution of the cluster mass and the dimensionless energy
$\beta\mathcal{E}\left(k\right)-\alpha_{k}k$ based on Eq. (7). Figure 2 shows
this energy functional plotted against the mass $k$ of particles within a
cluster for adhesion numbers ranging from 1 to about $2,000$. First, looking
at the small- and large-mass limits, it appears that $\lim_{k\to
0^{+}}\beta\mathcal{E}\left(k\right)-\alpha_{k}k=\lim_{k\to+\infty}\beta\mathcal{E}\left(k\right)-\alpha_{k}k=+\infty$,
which indicates that extreme masses are not accessible to the system. Second,
in all cases, there exists a global minimum of the energy functional that
shifts towards larger values of $k$ as the adhesion number $\mathrm{Ad}$ is
increased. As expected intuitively, this suggests that the average cluster
mass increases with $\mathrm{Ad}$. Moreover, the steeper slope
$\partial_{k}\left(\beta\mathcal{E}\left(k\right)-\alpha_{k}k\right)$ of the
energy functional on the right side of the global minimum than on the left
side indicates that the system reaches the energy minimum more rapidly when
starting from large masses than from small masses. This confirms the intuition
that disaggregation processes are much faster than aggregation processes.
Third, while the value of the global minimum energy does not show any clear
trend with $\mathrm{Ad}$, there may exist several local minima in
$\beta\mathcal{E}\left(k\right)-\alpha_{k}k$. This means that several
metastable states may occur. These states may disappear when increasing the
simulation duration or the system size. Still, this shows that a rather
polydisperse population of clusters may be found, at least transiently. This
may also be related to the dispersion in the energy minimum illustrated
previously.

Figure 3: (a) Mass $k_{\min}$ of the clusters corresponding to the energy minimum and (b) cluster size $\ell$ defined as the weighted average of the radius of gyration of the clusters, as a function of the adhesion number $\mathrm{Ad}$ for Couette flow ($\bullet$ ) and Poiseuille flow ($\square$) and for two different width $h=103$ (blue) and $h=206$ (orange). Solid and dotted lines correspond to the best power-law fits respectively for the Couette and the Poiseuille flow, and computed over $\mathrm{Ad}\in\left[2,200\right]$. (c) Cluster size $\ell$ as a function of the equilibrium mass $k_{\min}$. The black solid line corresponds to $\ell\propto k_{\min}^{1/2}$. Same symbols and colors as in (a) and (b).  Geometry | $h$ | Variable | Exponent | Prefactor
---|---|---|---|---
Couette | 103 | $\ell/a$ | $0.346\pm 0.009$ | $6.1\pm 0.2$
Couette | 103 | $k_{\min}$ | $0.4\pm 0.4$ | $400\pm 700$
Couette | 206 | $\ell/a$ | $0.57\pm 0.03$ | $3.5\pm 0.4$
Couette | 206 | $k_{\min}$ | $0.9\pm 0.1$ | $60\pm 30$
Poiseuille | 103 | $\ell/a$ | $0.271\pm 0.007$ | $6.1\pm 0.2$
Poiseuille | 103 | $k_{\min}$ | $0.5\pm 0.1$ | $170\pm 70$
Poiseuille | 206 | $\ell/a$ | $0.48\pm 0.03$ | $3.9\pm 0.5$
Poiseuille | 206 | $k_{\min}$ | $1.0\pm 0.1$ | $40\pm 20$
Table 1: Exponents and prefactors of the power-law fits in Fig. 3 for the
different geometries, widths, and sizes $\ell$ and $k_{\min}$.

Figure 3(a) and Fig. 3(b) respectively show the “equilibrium mass” $k_{\min}$
of the clusters corresponding to the energy minimum, i.e.,
$k_{\min}=\arg\min_{k\in\mathbb{R}_{+}}\left(\beta\mathcal{E}\left(k\right)-\alpha_{k}k\right)$,
and the cluster size $\ell/a$ taken as the weighted average of their radius of
gyration plotted against the adhesion number $\mathrm{Ad}$. For both
observables, a power-law regime is identified over almost two decades in
adhesion numbers. The exponents inferred from power-law fits for
$\mathrm{Ad}\in\left[2,200\right]$ are reported in Table 1. For
$\mathrm{Ad}\gtrsim 100$, a saturation is observed in $k_{\min}$ together with
large variations, for both the Couette and Poiseuille flows. This behavior is
most probably linked to finite-size effects as the average cluster size
becomes comparable or larger than the system width $h$. Consistently with Fig.
3(a), the cluster size $\ell$ increases as a power-law of $\mathrm{Ad}$ [Fig.
3(b)]. There, although the data for $\ell$ do not show such a strong
saturation as for $k_{\min}$, significant deviations from power-law behavior
are still observed for $\mathrm{Ad}\gtrsim 100$.

Moreover, the exponents for the dependence of $\ell$ with $\mathrm{Ad}$ seem
to depend significantly on the geometry, with values 0.35 and 0.57 for the
Couette flow and 0.27 and 0.48 for the Poiseuille flow, respectively for
$h=103$ and $h=206$ (see Table 1). Similarly, the corresponding exponents for
$k_{\min}$, namely 0.4 and 0.9 for Couette flow, and 0.5 and 1.0 for
Poiseuille flow, also differ for the two values of $h$. At this stage, the
reason for such an influence of $h$ in the simple two-dimensional simulations
remains unclear, and future work should focus on more realistic interaction
potentials and three-dimensional geometries.

Still, it is interesting to note that $\ell/a$ and $k_{\min}$ are not expected
to have the same power-law behavior with the adhesion number because of their
intrinsic relationship. More specifically, considering the framework of
fractal clusters, the number of particles in a cluster $k$ is related to the
geometrical size $l$ through $k=\left(l/a\right)^{D}$ with $D$ the fractal
dimension. In the two-dimensional case, $D$ should fall into the range
$\left[1,2\right]$. Figure 3(c) shows that $\ell/a\propto k_{\min}^{1/2}$,
which is compatible with a fractal dimension $D\simeq 2$, i.e., with almost
dense aggregates in two dimensions geometries. The fact that very compact
aggregates are obtained in the simulations is confirmed visually by looking at
Fig. 1(b). Another confirmation comes from the ratio of the exponent for
$k_{\min}$ and that for $\ell/a$ in Table 1, which is also found to be close
to 2 (except for the Couette flow with $h=103$ but the scatter of the
$k_{\min}$ data in this latter case makes it difficult to conclude).

##  IV Proposition of model

###  IV.1 Description of the disaggregation and reaggregation processes

Figure 4: Sketch of the disaggregation and reaggregation processes under a
uniform external stress $\sigma$. The numbers represent different sizes of
aggregates with different levels of energy. The aggregates are disaggregating
and reaggregating according to the stress solicitation. Transition from state
1 to states 2 and 4 corresponds to a “fragile” rupture, while transition from
state 1 to state 3 is an “erosion” process.

The precise internal structure of the aggregates is really complex, and
depends upon many different parameters including the volume fraction, the
nature of the interparticle forces and of the solvent, temperature, and
chemical environment. Hence, we shall consider the aggregates as a continuum
without further internal details. We consider a suspension of particles of
size $a$ that interact through an attractive potential so that they gather
into aggregates as sketched in Fig. 4. Following a statistical approach as in
Refs. [48, 71, 72, 73, 74], we assume that there exists a distribution
$f\left(t,l,a,U,\delta,\sigma\right)$ giving the population of aggregates of
size $l$ at time $t$ made of particles of size $a$ interacting through a pair
potential $U$ over a distance $\delta$ under a stress solicitation $\sigma$.
Such a distribution may be computed through a coagulation-fragmentation
equation as in Refs. [48, 71, 74]. Yet, this approach requires to identify
coagulation kernels and fragmentation coefficients, which involves much effort
for theoretical and/or numerical validation. Without ignoring the power of
such tools, we aim at a simpler approach through the present statistical
approach.

We start by imposing the conservation of the total number of particles $N$,
which reads:

$N=\int_{0}^{+\infty}n\left(\frac{l}{a}\right)f\left(t,l,a,U,\delta,\sigma\right)\,\mathrm{d}l=\mathrm{constant}\,,$
(19)

where $n\left(l/a\right)$ denotes the number of particles per cluster of
relative size $l/a$. The average size of the aggregates is thus given by:

$\overline{l}\left(t,a,U,\delta,\sigma\right)=\frac{a}{N}\int_{0}^{+\infty}n\left(\frac{l}{a}\right)\frac{l}{a}f\left(t,l,a,U,\delta,\sigma\right)\,\mathrm{d}l\,.$
(20)

From a statistical point of view, and assuming that the aggregates are
submitted to a uniform external stress $\sigma$ far from any boundary,
aggregates should evolve from one state to another as sketched in Fig. 4. If
the stress is sufficiently large to break some initial aggregate (state 1), a
disaggregation occurs either due to “fragile” rupture into two pieces of
similar sizes (bottom part of Fig. 4, state 4) or due to “erosion,” where
small pieces detach from the initial aggregate (top part of Fig. 4, states 2
and 3). If the larger pieces can still be broken down (states 2 and 4), the
process continues. However, if the pieces become too small (state 3), the
interparticle attraction dominates and reaggregation occurs. In the process,
the aggregates progressively decrease their global energy, until a minimum is
reached at long times $t\to+\infty$. Therefore, the steady state (state 5)
eventually corresponds to the optimum of all possible sizes, which results
from a dynamical equilibrium between disaggregation and reaggregation
processes. More formally, when a steady state is reached, the average
aggregate size $\ell$ is given by

$\ell\left(a,U,\delta,\sigma\right)=\lim_{t\to+\infty}\overline{l}\left(t,a,U,\delta,\sigma\right).$
(21)

In practice, since $t$ remains finite, we note that the longer the final time,
the smaller the spread of the distribution around the steady-state size. An
additional comment is that, keeping the previous set of variables for the
distribution $f$, an equivalent form using the Vaschy-Buckingham theorem [59,
60, 61] is found by replacing $f\left(t,l,a,U,\delta,\sigma\right)\mathrm{d}l$
with
$\tilde{f}\left(l/a,\mathrm{Ad},\delta/a\right)\mathrm{d}\left(l/a\right)$,
where $t$ can be discarded based on unit independence. Time becomes relevant,
however, when the viscosity $\eta$ of the suspending liquid or any other time-
related quantity is considered.

###  IV.2 Analytical formulation of the model

Let us consider a suspension of particles of diameter $a$ gathered in
aggregates. The bond between each particle involves an energy $U$ and an
interparticle distance of separation $\delta$. The system is submitted to a
uniform stress $\sigma$. One aims at estimating the steady-state size of
stable aggregates $\ell$ as a function of $a$, $U$, $\delta$, and $\sigma$.
Following Eggersdorfer  _et al._ [42], for a dense aggregate of size $l$, the
applied mechanical energy per unit area is $\sigma l$. However, aggregates are
not completely dense and one should more generally account for their fractal
nature. As already introduced above in Sect. III.2, the number of particles in
an aggregate of size $l$ is proportional to $\left(l/a\right)^{D}$ with $D$
the fractal dimension of the clusters. Hence, on the one hand, the mechanical
energy per unit of effective surface of the aggregate $M\left(l\right)$ is
proportional to $\sigma a\left(l/a\right)^{D+1}$. On the other hand, if one
isolates an aggregate of size $l$, the energy per unit surface liberated due
to broken bonds $E$ is $E\left(l\right)=U\delta^{-2}$, which is independent of
the size $l$ [54]. Indeed, considering an intermolecular potential
$w\left(d\right)$, with $d$ the intermolecular distance, one can compute the
interparticle potential per unit area through:

$W\left(d\right)=2\pi\varrho^{2}\intop_{d}^{+\infty}\intop_{y}^{+\infty}\intop_{0}^{+\infty}rw\left(\sqrt{r^{2}+z^{2}}\right)\,\mathrm{d}r\mathrm{d}z\mathrm{d}y\,,$
(22)

where $\varrho$ is the molecular density inside a particle. Defining $\delta$
as the distance that satisfies $W^{\prime}\left(\delta\right)=0$, or
equivalently
$\delta=\arg\min_{d\in\mathbb{R}_{+}^{*}}W\left(d\right),$222Taking some usual
examples leads to
$\delta\approx\arg\min_{d\in\mathbb{R}_{+}^{*}}w\left(d\right)$, and the
energy $U=\delta^{2}W\left(\mathfrak{\delta}\right)$,
$E\left(l\right)=U\delta^{-2}$ indeed corresponds to the energy per unit
surface liberated due to broken bonds. This energy is independent of the size
of the aggregate because one may assume that, on the boundary of the
aggregate, the number of particles per unit surface does not depend on $l$ and
is only related to the structure, which is assumed to be fixed in steady
state. Finally, the steady-state size $\ell$ corresponds to the size for which
the mechanical energy balances that due to broken bonds, i.e.,
$M\left(\ell\right)=E\left(\ell\right)$, which leads to:

$\frac{\ell}{a}=\left(\frac{U}{\sigma
a\delta^{2}}\right)^{\frac{1}{1+D}}=\mathrm{Ad}^{\frac{1}{1+D}}.$ (23)

The same result can be obtained by minimizing the total energy
$\mathcal{E}\left(l\right)=l^{2}\left(M\left(l\right)-E\left(l\right)\right)$.
Equation (23) also agrees with Eq. (21) through the Vaschy-Buckingham theorem
[59, 60, 61]. Indeed, Eq. (21) may be rewritten as:

$\frac{\ell}{a}=G\left(\mathrm{Ad},\frac{\delta}{a}\right).$ (24)

Adding that $U/\delta^{2}$ provides all the information about the interaction
potential, one has $\partial_{\delta/a}G=0$ and thus

$\frac{\ell}{a}=G\left(\mathrm{Ad}\right)\,,$ (25)

consistently with Eq. (23).

##  V Discussion

In this section, we discuss the theoretical approach and numerical results in
light of the literature. We start by comparing the model with the classical
coagulation-fragmentation approach. Then, the power-law scaling predicted for
the cluster size, $\ell\propto\mathrm{Ad}^{1/(1+D)}$, is confronted to the
present simulations and to previous experimental results.

###  V.1 Comparison with a coagulation-fragmentation model

The model proposed above in Sect. IV may be compared to the coagulation-
fragmentation approach introduced in the literature more than thirty years ago
[76] and subsequently enriched over the years, e.g., through the population
balance models [48, 71, 47, 46]. In such an approach, the probability density
function $f$ of clusters of size $x$ at time $t$ obeys the following dynamical
equation:

$\frac{\partial f}{\partial t}\left(t,x,\bullet\right)=\\\
\frac{1}{2}\int_{0}^{x}K\left(y,x-y,\bullet\right)f\left(t,y,\bullet\right)f\left(t,x-y,\bullet\right)\,\mathrm{d}y\\\
-\frac{1}{2}\int_{0}^{x}F\left(y,x-y,\bullet\right)f\left(t,x,\bullet\right)\,\mathrm{d}y\\\
-\int_{0}^{+\infty}K\left(x,y,\bullet\right)f\left(t,x,\bullet\right)f\left(t,y,\bullet\right)\,\mathrm{d}y\\\
+\int_{0}^{+\infty}F\left(x,y,\bullet\right)f\left(t,x+y,\bullet\right)\,\mathrm{d}y\,,$
(26)

with $\bullet=\left(a,U,\delta,\sigma\right)$, $K$ the aggregation kernel, and
$F$ the fragmentation kernel. Considering the long-term behaviour, it can be
shown that $f$ takes the general form [48, 71, 76]:

$f\left(t,x,\bullet\right)=\frac{1}{s^{2}\left(t,\bullet\right)}\,\varphi\left(\frac{x}{s\left(t,\bullet\right)},\bullet\right),$
(27)

where $s$ and $\varphi$ are two functions that depend only on time and size
respectively, and such that
$\int_{\mathbb{R}_{+}}x\varphi\left(x,\bullet\right)\,\mathrm{d}x=N$, the
fixed total number of particles. The function $s$ corresponds to the average
mass of the clusters according to the distribution $f$. It can be related to a
size through the cluster fractal dimension $D$ by $s\propto l^{D}$. Then,
integrating the previous expression over the size as in Eq. 20, one computes
the steady-state average size as:

$\overline{l}\left(t,\bullet\right)=as^{\frac{1}{D}}\left(t,\bullet\right)\,,$
(28)

with

$\displaystyle\lim_{t\to+\infty}s\left(t,\bullet\right)$
$\displaystyle=\left(N\frac{A\left(\bullet\right)}{B\left(\bullet\right)}\right)^{\frac{1}{\alpha+2-\lambda}}$
(29) $\displaystyle A\left(\bullet\right)$
$\displaystyle=\iint_{\mathbb{R}_{+}^{2}}xy\varphi\left(x,\bullet\right)\varphi\left(y,\bullet\right)K\left(x,y,\bullet\right)\,\mathrm{d}x\mathrm{d}y$
(30) $\displaystyle B\left(\bullet\right)$
$\displaystyle=\iint_{\mathbb{R}_{+}^{2}}xy\varphi\left(x+y,\bullet\right)F\left(x,y,\bullet\right)\,\mathrm{d}x\mathrm{d}y\,,$
(31)

where $\lambda$ and $\alpha$ are the respective homogeneity coefficients of
$K$ and $F$ 333The homegeneity coefficients are defined as $K\left(\xi x,\xi
y\right)=\xi^{\lambda}K\left(x,y\right)$ and $F\left(\xi x,\xi
y\right)=\xi^{\alpha}F\left(x,y\right)$ for all
$\left(x,y,\xi\right)\in\mathbb{R}_{+}^{3}$., assuming $\alpha+2>\lambda$. It
follows from Eq. 21 that

$\frac{\ell}{a}=\left(N\frac{A\left(\bullet\right)}{B\left(\bullet\right)}\right)^{\frac{1}{D(\alpha+2-\lambda)}}.$
(32)

Finally, identifying with Eq. (23), one gets

$\frac{U}{\sigma
a\delta^{2}}=\left(N\frac{A\left(\bullet\right)}{B\left(\bullet\right)}\right)^{\frac{1+D}{D(\alpha+2-\lambda)}}\,.$
(33)

Most of the physical quantities appear as multiplicative factors, in the sense
that, for example, the coagulation kernel $K$ is usually built with a factor
$U/\eta$ with $\eta$ the viscosity of the fluid and no other dependence on
physical quantities besides the variables $x$ and $y$. Therefore, the ratio
$NA\left(\bullet\right)/B\left(\bullet\right)$ can be considered as the
product of a dimensionless number built in a similar manner as the adhesion
number and another factor that depends only on the shape of respective kernels
without involving any additional physical parameter. Therefore, in order to
keep Eq. (33) true in general, one should impose that the exponent is 1, which
leads to

$1=D\left(\alpha+1-\lambda\right)\,.$ (34)

This is an important result which, to our knowledge, has not been reported in
the literature before. Indeed, Eq. (34) allows one to relate the first
mechanical approach to population balance models, where the homogeneity
coefficients may seem disconnected from physical and measurable quantities.
This result also emphasizes the fact that the adhesion number is a correct
measure of the relative importance of aggregation versus disaggregation.

We note that aggregation kernels are pretty well covered, either in terms of
theoretical solutions [48, 71, 78, 79] with simple sums or products, or in
terms of a physical construction, e.g., based on collisions, thermal
fluctuations, and diffusion [80, 81]. Fragmentation kernels, however, are more
poorly controlled, in the sense that some theoretical solutions impose some
strong conditions on these kernels without much physical justification [48,
71, 47, 46], so that the expressions of the fragmentation kernels remain
mostly empirical or semi-empirical [82, 83]. Therefore, although the use of
coagulation-fragmentation equations is well established, relating some of the
main terms to physical phenomena, such as the interaction potential and the
flow stresses, is an important step yet to be fully achieved.

###  V.2 Comparison with simulations and experiments

Relating the model proposed in Sect. IV to the simulations of Sect. III, we
expect the exponents in Table 1 to be linked to the fractal dimension
respectively by $1/\left(1+D\right)$ for the steady-state size $\ell$, and
$D/\left(1+D\right)$ for the equilibrium mass $k_{\min}$. The broad
variability of the exponents does not allow to properly extract a fractal
dimension from the simulations. Nevertheless, as discussed above in Sect.
III.2, the shear-induced clusters are almost dense, so that we may assume a
fractal dimension $D$ close to 2. Such a compactness most probably results
from the specific interaction potential, i.e. a 12-6 Lennard-Jones potential,
which is a central-force potential and is likely to lead to clusters with
droplet-like shapes. We note that $D=2$ would be consistent with the exponent
$\beta=\left(1+D\right)^{-1}=1/3$ expected for $\ell$ as a function of
$\mathrm{Ad}$, at least in the smaller geometry, since Table 1 reports
exponents of 0.35 and 0.27 for $h=103$.

Moreover, the results may be compared to the recent three-dimensional
simulations under simple shear flow by Ruan et al. [43], which also report
rather dense shear-induced clusters at steady state. Focusing both on the
cluster size and on the average number of particles in a cluster as in the
present work, the authors extract the cluster fractal dimension $D$ and show
that $\ell/a$ scales as a power-law of the shear stress. The exponents $\beta$
for $\ell$ reported in Ref. [43] are replotted as a function of $D$ as blue
circles in Fig. 5 and show good agreement with the prediction
$\beta=\left(1+D\right)^{-1}$. Note that this is also consistent with the
empirical expression proposed in Ref. [43], namely $\beta=S/(S+\langle
D\rangle)$, where $\langle D\rangle$ denotes the average cluster fractal
dimension, and the fitting parameter $S$ is found to be close to 1.

On the experimental side, a number of works have reported results assessing
the size of aggregates under the application of a mechanical stress. However,
only a handful of papers explicitly state the values of the various parameters
$a$, $U$, $\sigma$, $\delta$, and $D$. Table 2 and Fig. 5 gather a selection
of such previous works. First, Table 2 shows that Eq. (23) predicts a typical
cluster size $\ell_{\mathrm{th}}$ which is of the same order as the
experimental size $\ell_{\mathrm{exp}}$. All the parameters fall into commonly
known ranges and justify that the previous approach may be a good proxy to
evaluate the most stable cluster size. Second, Fig. 5 probes the sensitivity
of $\ell/a$ with the adhesion number $\mathrm{Ad}$ by plotting the power-law
exponent $\beta$ in $\ell/a=\mathrm{Ad}^{\beta}$ as a function of $D$. Except
for a few points that lie far from the average estimation, most exponents are
gathered along the prediction of Eq. (23), consistently with the numerical
results of Ref. [43], which provides strong support for the approach described
in Sect. IV.

Figure 5: Power-law exponent $\beta$ in $\ell/a=\mathrm{Ad}^{\beta}$ as a function of the fractal dimension $D$ of the aggregates extracted from Refs. [43, 27, 84, 85, 86, 51, 26, 52, 87, 88, 89, 90, 11, 28, 29, 42, 41] ($\bullet$) and compared to the model prediction $\beta=\left(1+D\right)^{-1}$ (solid line).  Source | [91, 92, 93, 36] | [21, 94, 95, 96]
---|---|---
$\ell_{\mathrm{exp}}$ (nm) | 500 | 60
$a$ (nm) | 150 | 6.5
$U/k_{B}T$ | 20 | 60
$\sigma$ (MPa) | 0.1 | 2
$\delta\,\left(\mathring{\mathrm{A}}\right)$ | 7 | 2
$D$ | 2.6 | 1.88
$\ell_{\mathrm{th}}$ (nm) | 300 | 55
Table 2: Experimental parameters and measurements of the aggregate size from
the literature. $\ell_{\mathrm{exp}}$ is the direct measurement of the
aggregate size, while $\ell_{\mathrm{th}}$ is the aggregate size estimated
using Eq. (23) based on the parameters defined in the text and reported in the
references of the first column.

Finally, Fig. 6 provides a sensitivity study of the value of $\ell$ upon the
different parameters of Eq. (23) based on Table 2. It is clear that the most
critical parameters are $D$ and $\delta$ as expected from Eq. (23). Therefore,
particular attention must be considered to assess accurately these parameters.
Yet, when one of these two parameters are unknown, the model can be used to
estimate $D$ or $\delta$ with good accuracy based on measurements of $\ell$.

Figure 6: Sensitivity study of $\ell$ predicted by Eq. (23) upon relative
variations of the different control parameters around the values listed in the
first line of Table 2 except for the fractal dimension, which is varied from 1
to 3.

##  VI Conclusion

We have investigated simple two-dimensional numerical simulations and proposed
a theoretical description of aggregation and disaggregation processes in
colloidal dispersions under shear. We have shown that the classically reported
power-law increase of the average cluster size with the adhesion number can be
rationalized through a dynamical balance of energy between aggregation and
disaggregation. In particular, the algebraic dependence of the power-law
exponent with the cluster fractal dimension predicted by the approach appears
to be in good agreement with simulations and previous experiments. Still, one
should remain careful about a direct inference of one simple general formula
as provided here, since many additional parameters may be involved in the size
selection, including time through aging phenomena, spatial confinement, and
individual cluster dynamics.

###### Acknowledgements.

The authors wish to thank Emanuela Del Gado, Thibaut Divoux, and Thomas Gibaud
for fruitful discussions.

## References

  * Macosko [1994] C. W. Macosko,  _Rheology : Principles, Measurements and Applications_ , edited by C. W. Macosko (John Wiley and Sons, 1994).
  * Mewis and Wagner [2012] J. Mewis and N. J. Wagner,  _Colloidal Suspension Rheology_, edited by J. Mewis and N. J. Wagner (Cambridge Univesity Press, 2012).
  * Wagner and Mewis [2021] N. J. Wagner and J. Mewis,  _Theory and Applications of Colloidal Suspension Rheology_, edited by N. J. Wagner and J. Mewis (Cambridge University Press, 2021).
  * Hengl  _et al._ [2014] N. Hengl, Y. Jin, F. Pignon, S. Baup, R. Mollard, N. Gondrexon, A. Magnin, L. Michot, and E. Paineau, A new way to apply ultrasound in cross-flow ultrafiltration: Application to colloidal suspensions, Ultrasonics Sonochemistry 21, 1018 (2014).
  * Ioannidou  _et al._ [2016] K. Ioannidou, K. J. Krakowiak, M. Bauchy, C. G. Hoover, E. Masoero, S. Yip, F.-J. Ulm, P. Levitz, R. J.-M. Pellenq, and E. Del Gado, Mesoscale texture of cement hydrates, Proceedings of the National Academy of Sciences of the United States of America 113, 2029 (2016).
  * Awad  _et al._ [2012] T. S. Awad, H. A. Moharram, O. E. Shaltout, and M. M. Youssef, Applications of ultrasound in analysis, processing and quality control of food: A review, Food Research International 48, 410 (2012).
  * Knorr  _et al._ [2004] D. Knorr, M. Zenker, V. Heinz, and D.-U. Lee, Applications and potential of ultrasonics in food processing, Trends in Food Science & Technology 15, 261 (2004).
  * Chandrapal  _et al._ [2012] J. Chandrapal, C. Oliver, S. Kentish, and M. Ashokkumar, Ultrasonics in food processing, Ultrasonics Sonochemistry 19, 975 (2012).
  * Bonn  _et al._ [2017] D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, Yield stress materials in soft condensed matter, Reviews of Modern Physics 89, 1 (2017), 035005.
  * Vassileva  _et al._ [2005] N. D. Vassileva, D. van den Ende, F. Mugele, and J. Mellema, Capillary Forces between spherical Particles Floating at a Liquid-Liquid Interface, Langmuir 21, 11190 (2005).
  * Wessel and Ball [1992] R. Wessel and R. C. Ball, Fractal aggregates and gels in shear flow, Physical Review A 46, R3008 (1992).
  * Shih  _et al._ [1990] W.-H. Shih, W. Y. Shih, S.-I. Kim, J. Liu, and I. A. Aksay, Scaling behavior of the elastic properties of colloidal gels, Physical Review A 42, 4772 (1990).
  * Sorensen [2001] C. M. Sorensen, Light Scattering by Fractal Aggregates: A Review, Aerosol Science Technology 35, 648 (2001).
  * Lin  _et al._ [1990] M. Y. Lin, R. Klein, H. M. Lindsay, D. A. Weitz, R. C. Ball, and P. Meakin, The Structure of Fractal Colloidal Aggregates of Finite Extent, Journal of Colloid and Interface Science 137, 263 (1990).
  * Lazzari  _et al._ [2016] S. Lazzari, L. Nicoud, B. Jaquet, M. Lattuada, and M. Morbidelli, Fractal-like structures in colloid science, Advances in Colloid and Interface Science 235, 1 (2016).
  * Grassberger [1985] P. Grassberger, On the spreading of two-dimensional percolation, Journal of Physics A: Mathematical and General 18, L215 (1985).
  * Bantawa  _et al._ [2021] M. Bantawa, W. A. Fontaine-Seiler, P. D. Olmsted, and E. Del Gado, Microscopic interactions and emerging elasticity in model soft particulate gels, Journal of Physics: Condensed Matter 33, 1 (2021).
  * Weitz and Oliveria [1984] D. A. Weitz and M. Oliveria, Fractal Structures Formed by Kinetic Aggregation of Aqueous Gold Colloids, Physical Review Letters 52, 1433 (1984).
  * Weitz  _et al._ [1985] D. A. Weitz, J. S. Huang, M. Y. Lin, and J. Sung, Limits of the Fractal Dimension for Irreversible Kinetic Aggregation of Gold Colloids, Physical Review Letters 54, 1416 (1985).
  * Hoekstra  _et al._ [2003] H. Hoekstra, J. Vermant, J. Mewis, and G. Fuller, Flow-induced anisotropy and reversible aggregation in two-dimensional suspensions, Langmuir 19, 9134 (2003).
  * Nguyen  _et al._ [2011] V. S. Nguyen, D. Rouxel, R. Hadji, B. Vincent, and Y. Fort, Effect of ultrasonication and dispersion stability on the cluster size of alumina nanoscale particles in aqueous solutions, Ultrasonics Sonochemistry 18, 382 (2011).
  * Richards  _et al._ [2017] J. J. Richards, J. B. Hipp, J. K. Riley, N. J. Wagner, and P. D. Butler, Clustering and Percolation in Suspensions of Carbon Black, Langmuir 33, 12260 (2017).
  * Masschaele  _et al._ [2009] K. Masschaele, J. Fransaer, and J. Vermant, Direct visualization of yielding in model two-dimensional colloidal gels subjected to shear flow, Journal of Rheology 53, 1437 (2009).
  * Hoekstra  _et al._ [2005] H. Hoekstra, J. Mewis, T. Narayanan, and J. Vermant, Multi length scale analysis of the microstructure in sticky sphere dispersions during shear flow, Langmuir 21, 11017 (2005).
  * Zaccone  _et al._ [2009] A. Zaccone, M. Soos, M. Lattuada, H. Wu, M. U. BÄBLER, and M. Morbidelli, Breakup of dense colloidal aggregates under hydrodynamic stresses, Physical Review E 79, 061401 (2009).
  * Torres  _et al._ [1991a] F. E. Torres, W. B. Russel, and W. R. Schowalter, Floc Structure and Growth Kinetics for Rapid Shear Coagulationof Polystyrene Colloids, Journal of Colloid and Interface Science 142, 554 (1991a).
  * Hunter and Frayne [1980] R. J. Hunter and J. Frayne, Flow Behavior of Coagulated Colloidal Sols, Journal of Colloid and Interface Science 76, 107 (1980).
  * Higashitani and Iimura [1998] K. Higashitani and K. Iimura, Two-Dimensional Simulation of the Breakup Process of Aggregates in Shear and Elongational Flows, Journal of Colloid and Interface Science 204, 320 (1998).
  * Higashitani  _et al._ [2001] K. Higashitani, K. Iimura, and H. Sanda, Simulation of deformation and breakupof large aggregates in Flows of viscous Fluids, Chemical Engineering Science 56, 2927 (2001).
  * Massaro  _et al._ [2020] R. Massaro, G. Colombo, P. Van Puyvelde, and J. Vermant, Viscoelastic cluster densification in sheared colloidal gels, Soft Matter 16, 2437 (2020).
  * Brady and Bossis [1985] J. F. Brady and G. Bossis, The rheology of concentrated suspensions of spheres in simple shear flow by numerical simulation, Journal of Fluid Mechanics 155, 105 (1985).
  * Morris [2009] J. F. Morris, A review of microstructure in concentrated suspensions and its implications for rheology and bulk flow, Rheologica Acta 48, 909 (2009).
  * Varga  _et al._ [2015] Z. Varga, G. Wang, and J. Swan, The hydrodynamics of colloidal gelation, Soft Matter 11, 9009 (2015).
  * Varga and Swan [2015] Z. Varga and J. W. Swan, Linear viscoelasticity of attractive colloidal dispersions, Journal of Rheology 59, 1271 (2015), https://doi.org/10.1122/1.4928951 .
  * Varga and Swan [2018] Z. Varga and J. W. Swan, Large scale anisotropies in sheared colloidal gels, Journal of Rheology 62, 405 (2018), https://doi.org/10.1122/1.5003364 .
  * Varga  _et al._ [2019] Z. Varga, V. Grenard, S. Pecorario, N. Taberlet, V. Dolique, S. Manneville, T. Divoux, G. H. McKinley, and J. W. Swan, Hydrodynamics control shear-induced pattern formation in attractive suspensions, Proceedings of the National Academy of Sciences 116, 12193 (2019).
  * Lorenzo and Marco [2022] T. Lorenzo and L. Marco, Brownian Dynamics simulations of shear-induced aggregation of charged colloidal particles in the presence of hydrodynamic interactions, Journal of Colloid and Interface Science 10.1016/j.jcis.2022.05.047 (2022).
  * Jungblut  _et al._ [2019] S. Jungblut, J.-O. Joswig, and A. Eychmüller, Diffusion- and reaction-limited cluster aggregation revisited, Physical Chemistry Chemical Physics 21, 5723 (2019).
  * Colombo and Del Gado [2014a] J. Colombo and E. Del Gado, Stress localization, stiffening, and yielding in a model colloidal gel, Journal of Rheology 58, 1089 (2014a).
  * Bouzid  _et al._ [2018] M. Bouzid, B. Keshavarz, M. Geri, T. Divoux, E. Del Gado, and G. H. McKinley, Computing the linear viscoelastic properties of soft gels using an optimally windowed chirp protocol, Journal of Rheology 62, 1037 (2018).
  * Kimbonguila Manounou and Rémond [2014] A. Kimbonguila Manounou and S. Rémond, Discrete element modeling of the microstructure of fine particle agglomerates in sheared dilute suspension, Physica A 412, 66 (2014).
  * Eggersdorfer  _et al._ [2010] M. L. Eggersdorfer, D. Kadau, H. J. Herrmann, and S. E. Pratsinis, Fragmentation and restructuring of soft-agglomerates under shear, Journal of Colloid and Interface Science 342, 261 (2010).
  * Ruan  _et al._ [2020] X. Ruan, S. Chen, and S. Li, Structural evolution and breakage of dense agglomerates in shear flowand Taylor-Green vortex, Chemical Engineering Science 211, 115261 (2020).
  * Soos  _et al._ [2006] M. Soos, J. Sefcik, and M. Morbidelli, Investigation of aggregation, breakage and restructuring kinetics of colloidal dispersions in turbulent flows by population balance modeling and static light scattering, Chemical Engineering Science 61, 2349 (2006).
  * Soos  _et al._ [2007] M. Soos, L. Wang, R. O. Fox, J. Sefcik, and M. Morbidelli, Population balance modeling of aggregation and breakage in turbulent taylor–couette flow, Journal of Colloid and Interface Science 307, 433 (2007).
  * Puisto  _et al._ [2012] A. Puisto, X. Illa, M. Mohtaschemi, and M. J. Alava, Modeling the viscosity and aggregation of suspensions of highly anisotropic nanoparticles, The European Physical Journal E 35, 1 (2012).
  * Lattuada  _et al._ [2016] M. Lattuada, A. Zaccone, H. Wu, and M. Morbidelli, Population-balance description of shear-induced clustering, gelation and suspension viscosity in sheared DLVO colloids, Soft Matter 12, 5313 (2016).
  * Banasiak  _et al._ [2020a] J. Banasiak, W. Lamb, and P. Laurençot,  _Analytic Methods for Coagulation-Fragmentation Models, Volume I_ , edited by J. Banasiak, W. Lamb, and P. Laurençot, Monographs and Research Notes in Mathematics, Vol. 1 (CRC Press, CRC PressTaylor & Francis Group6000 Broken Sound Parkway NW, Suite 300Boca Raton, FL 33487-2742, 2020).
  * Kantor and Witten [1984] Y. Kantor and T. A. Witten, Mechanical stability of tenuous objects, Journal de Physique Lettres 45, 675 (1984).
  * Mellema  _et al._ [2002] M. Mellema, J. H. J. van Opheusde, and T. van Vliet, Categorization of rheological scaling models for particlegels applied to casein gels, Journal of Rheology 46, 11 (2002).
  * Brakalov [1987] L. B. Brakalov, Connection between the orthokinetic coagulation capture efficiency of aggregates and their maximum size, Chemical Engineering Science 42, 2373 (1987).
  * Torres  _et al._ [1991b] F. E. Torres, W. B. Russel, and W. R. Schowalter, Simulation of Coagulation in Viscous Flows, Journal of Colloid and Interface Science 145, 51 (1991b).
  * Snabre and Mills [1996] P. Snabre and P. Mills, I. Rheology of Weakly Flocculated Suspensions of Rigid Particles, Journal de Physique III 6, 1811 (1996).
  * Marshall and Li [2014] J. S. Marshall and S. Li,  _Adhesive Particle Flow: A Discrete-Element Approach_, 1st ed., edited by J. S. Marshall and S. Li (Cambridge University Press, 32 Avenue of the Americas, NY 10013-2473, USA, 2014).
  * von Mises [1913] R. von Mises, Mechanik der festen Körper im plastisch-deformablen Zustand, Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen 1, 582 (1913).
  * Griffith [1921] A. A. Griffith, VI. The phenomena of rupture and flow in solids, Philosophical Transactions of the Royal Society A 221, 582 (1921).
  * Irwin [1957] G. R. Irwin, Analysis of Stresses and Strains Near the End of a Crack Traversing a Plate, Journal of Applied Mechanics 24, 361 (1957).
  * Creton and Ciccotti [2016] C. Creton and M. Ciccotti, Fracture and adhesion of soft materials: a review, Reports on Progress in Physics 79, 046601 (2016).
  * Buckingham [1914] E. Buckingham, On Physically Similar Systems; Illustrations of the Use of Dimensional Equations, Physical Review 4, 345 (1914).
  * Buckingham [1915a] E. Buckingham, The Principle of Similitude, Nature 96, 396 (1915a).
  * Buckingham [1915b] E. Buckingham, Model experiments and the forms of empirical equations, in  _Transactions of the American Society of Mechanical Engineers_ , Vol. 37, American Society of Mechanical Engineers (American Society of Mechanical Engineers, 1915) pp. 263–296.
  * Jamali  _et al._ [2020] S. Jamali, R. C. Armstrong, and G. H. McKinley, Time-rate-transformation framework for targeted assembly of short-range attractive colloidal suspensions, Materials Today Advances 5, 100026 (2020).
  * Nabizadeh and Jamali [2021] M. Nabizadeh and S. Jamali, Life and death of colloidal bonds control the rate-dependent rheology of gels, Nature Communications 12, 1 (2021).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Computer Physics Communications 271, 108171 (2022).
  * Stukowski [2010] A. Stukowski, Visualization and analysis of atomistic simulation data with OVITO-the Open Visualization Tool, Modelling and Simulation in Materials Science and Engineering 18, 015012 (2010).
  * Note [1] Additional simulations were carried out with various values of $C$ and the main results remained similar up to a rescaling of the time scale.
  * Swope  _et al._ [1982] W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson, A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters, The Journal of Chemical Physics 76, 637 (1982), https://doi.org/10.1063/1.442716 .
  * Pearce [2005] D. J. Pearce,  _An improved algorithm for finding the strongly connected components of a directed graph_, techreport (Victoria University, 2005).
  * Colombo  _et al._ [2013] J. Colombo, A. Widmer-Cooper, and E. Del Gado, Microscopic Picture of Cooperative Processes in Restructuring Gel Networks, Physical Review Letters 110, 198301 (2013).
  * Colombo and Del Gado [2014b] J. Colombo and E. Del Gado, Self-assembly and cooperative dynamics of a model colloidal gel network, Soft Matter 10, 4003 (2014b).
  * Banasiak  _et al._ [2020b] J. Banasiak, W. Lamb, and P. Laurençot,  _Analytic Methods for Coagulation-Fragmentation Models, Volume II_ , edited by J. Banasiak, W. Lamb, and P. Laurençot, Monographs and Research Notes in Mathematics, Vol. 2 (CRC Press, CRC PressTaylor & Francis Group6000 Broken Sound Parkway NW, Suite 300Boca Raton, FL 33487-2742, 2020).
  * Golse [2005] F. Golse, Handbook of Differential Equations: Evolutionary Equations (Elsevier-North-Holland, 2005) Chap. The Boltzmann Equation and Its Hydrodynamic Limits, pp. 159–301.
  * Alexeev [2004] B. V. Alexeev,  _Generalized Boltzmann Physical Kinetics_ , edited by B. V. Alexeev (Elsevier, 2004).
  * Stadnichuk  _et al._ [2015] V. Stadnichuk, A. Bodrova, and N. Brilliantov, Smoluchowski aggregation fragmentation equations: Fast numerical method to find steady-state solutions, International Journal of Modern Physics B 29, 1 (2015).
  * Note [2] Taking some usual examples leads to $\delta\approx\mathop{arg}\nolimits\mathop{min}\displaylimits_{d\in\mathbb{R}_{+}^{*}}w\left(d\right)$.
  * Sorensen  _et al._ [1987] C. M. Sorensen, H. X. Zhang, and T. W. Taylor, Cluster-size evolution in a coagulation-fragmentation system, Physical Review Letters 59, 363 (1987).
  * Note [3] The homegeneity coefficients are defined as $K\left(\xi x,\xi y\right)=\xi^{\lambda}K\left(x,y\right)$ and $F\left(\xi x,\xi y\right)=\xi^{\alpha}F\left(x,y\right)$ for all $\left(x,y,\xi\right)\in\mathbb{R}_{+}^{3}$.
  * Wattis [2006] J. A. D. Wattis, An introduction to mathematical models of coagulation–fragmentation processes: A discrete deterministic mean-field approach, Physica D 222, 1 (2006).
  * Spicer and Pratsinis [1996] P. T. Spicer and S. E. Pratsinis, Coagulation and fragmentation: Universal steady-state particle-size distribution, AIChE Journal 42, 1612 (1996).
  * Kryven  _et al._ [2014] I. Kryven, S. Lazzari, and G. Storti, Population Balance Modeling of Aggregation and Coalescence in Colloidal Systems, Macromolecular Theory and Simulations 23, 170 (2014), https://onlinelibrary.wiley.com/doi/pdf/10.1002/mats.201300140 .
  * Barthelmes  _et al._ [2003] G. Barthelmes, S. E. Pratsinis, and H. Buggisch, Particle size distributions and viscosity of suspensions undergoing shear-induced coagulation and fragmentation, Chemical Engineering Science 58, 2893 (2003).
  * Delichatsios and Probstein [1976] M. A. Delichatsios and R. F. Probstein, The Effect of Coalescence on the Average Drop Size in Liquid-Liquid Dispersions, Industrial & Engineering Chemistry Fundamentals 15, 134 (1976), https://doi.org/10.1021/i160058a010 .
  * Kusters [1991] K. A. Kusters,  _The influence of turbulence on aggregation of small particles in agitated vessels_ , phdthesis, Chemical Engineering and Chemistry (1991).
  * Sonntag and Russel [1986] R. C. Sonntag and W. B. Russel, Structure and Breakup of Flocs Subjected to Fluid Stresses: I. Shear Experiments, Journal of Colloid and Interface Science 115, 390 (1986).
  * Sonntag and Russel [1987a] R. C. Sonntag and W. B. Russel, Structure and Breakup of Flocs Subjected to Fluid Stresses: II. Theory, Journal of Colloid and Interface Science 115, 378 (1987a).
  * Sonntag and Russel [1987b] R. C. Sonntag and W. B. Russel, Structure and Breakup of Flocs Subjected to Fluid Stresses: III. Converging Flow, Journal of Colloid and Interface Science 113, 399 (1987b).
  * Harshe  _et al._ [2011] Y. M. Harshe, M. Lattuada, and M. Soos, Experimental and Modeling Study of Breakage and Restructuring of Open and Dense Colloidal Aggregates, Langmuir 27, 5739 (2011).
  * Potanin [1991] A. A. Potanin, On the Mechanism of Aggregation in the Shear Flow of Suspensions, Journal of Colloid and Interface Science 145, 140 (1991).
  * Potanin [1992] A. A. Potanin, On the model of colloid aggregates and aggregating colloids, Journal of Chemical Physics 96, 9191 (1992).
  * Potanin and Russel [1996] A. A. Potanin and W. B. Russel, Fractal model of consolidation of weakly aggregated colloidal dispersions, Physical Review E 53, 3702 (1996).
  * Gibaud  _et al._ [2020a] T. Gibaud, N. Dagès, P. Lidon, G. Jung, L. C. Ahouré, M. Sztucki, A. Poulesquen, N. Hengl, F. Pignon, and S. Manneville, Rheoacoustic Gels: Tuning Mechanical and Flow Properties of Colloidal Gels with Ultrasonic Vibrations, Physical Review X 10, 1 (2020a), 011028, 1905.07282 .
  * Gibaud  _et al._ [2020b] T. Gibaud, T. Divoux, and S. Manneville, Encyclopedia of Complexity and Systems Science (Springer, 2020) Chap. Nonlinear mechanics of colloidal gels: creep, fatigue and shear-induced yielding, pp. 1–24.
  * Dagés  _et al._ [2021] N. Dagés, P. Lidon, G. Jung, F. Pignon, S. Manneville, and T. Gibaud, Mechanics and structure of carbon black gels under high-power ultrasound, Journal of Rheology 65, 477 (2021).
  * Klimchitskaya  _et al._ [2000] G. L. Klimchitskaya, U. Mohideen, and V. M. Mostepanenko, Casimir and Van der Waals force between two plates or a sphere (lens) above a plate made of real metals, Physical Review A 61, 062107 (2000).
  * Visser [1972] J. Visser, On Hamaker constants: A comparison between Hamaker constants and Lifshitz-Van der Waals constants, Advances in Colloid and Inferface Science 3, 331 (1972).
  * Waite  _et al._ [2001] T. D. Waite, J. K. Cleaver, and J. K. Beattie, Aggregation Kinetics and Fractal Structure of $\gamma$-Alumina Assemblages, Journal of Colloid and Interface Science 241, 333 (2001).