# Optimal s-boxes against alternative operations

Marco Calderini  Department of Mathematics - University of Trento
<EMAIL_ADDRESS>,  Roberto Civino  DISIM - University of L’Aquila
<EMAIL_ADDRESS>and  Riccardo Invernizzi  COSIC - KU Leuven
<EMAIL_ADDRESS>
###### Abstract.

Civino et al. have characterised diffusion layers that expose an SPN to
vulnerability from differential cryptanalysis when employing alternative
operations coming from groups isomorphic to the translation group on the
message space. In this study, we present a classification of diffusion layers
that exhibit linearity in _parallel_ alternative operations for ciphers with
4-bit s-boxes, enabling the possibility of an alternative differential attack
simultaneously targeting all the s-boxes within the block. Furthermore, we
investigate the differential behaviour with respect to alternative operations
for all classes of optimal 4-bit s-boxes, as defined by Leander and Poschmann
(2007). Our examination reveals that certain classes contain weak permutations
w.r.t. alternative differential attacks, and we leverage these vulnerabilities
to execute a series of experiments.

###### Key words and phrases:

Differential cryptanalysis; alternative operations; 4-bit s-boxes

###### 2010 Mathematics Subject Classification:

20B35, 94A60, 68P25

M. Calderini and R. Civino are members of INdAM-GNSAGA (Italy). R. Civino is
partially funded by the Centre of excellence ExEMERGE at University of
L’Aquila. R. Invernizzi is supported by the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement ISOCRYPT - No. 101020788) and by CyberSecurity Research
Flanders with reference number VR20192203.

##  1\. Introduction and preliminaries

Differential cryptanalysis, originally introduced by Biham and Shamir in the
late 1980s [BS91] and subsequently generalised [Wag99, Knu95, BCJW02, BBS05],
has become one of the cornerstones for evaluating the robustness of various
symmetric primitives. The fundamental premise of differential cryptanalysis is
that analysing the differences (differentials) between pairs of plaintexts and
the corresponding ciphertexts can unveil undesired biases. While differentials
can be calculated with respect to any difference operator, regardless of which
operation is responsible for performing the sum with the round key during
encryption, it is usual for the two operations to coincide. For this reason,
classical differential cryptanalysis of a cipher in which the key is xor-ed to
the state is typically performed by studying the distribution of xor-
differentials, whose propagation is traditionally prevented by the combined
action of the linear diffusion layer and the s-box layer. In particular,
s-boxes are pivotal for ensuring the security of almost all contemporary block
ciphers, serving as the primary non-linear component within the cipher,
particularly in the case of SPNs. Equally relevant, the efficiency of a cipher
is significantly influenced by the size of the s-boxes. In practical
scenarios, s-boxes typically have a size of 4 or 8 bits, with 4 being the most
popular choice for ciphers designed to operate on power-constrained devices
[BAK98, BKL+07, SIH+11, BBI+15]. It is clear that the selection of appropriate
s-boxes is critical to fortify the cipher against various types of attacks. In
this sense, Leander and Poschmann have classified 4-bit s-boxes which are
optimal w.r.t. standard criteria that guarantee poor propagation of xor-
differentials [LP07].

A recent line of research is focused on the study of alternative difference
operators for the differential cryptanalysis of xor-based ciphers [CBS19,
CCS21, Teş22, CCI24]. These new operators are designed to induce a novel
operation with respect to which differentials are computed. Within this
approach, a large class of possible alternative operations has been studied,
all of which have in common that they are induced by a group of translations
isomorphic to the group of translations acting on the message space by means
of the xor addition with the key. In the context of an SPN, where the
encrypted message is generated by iterating through a sequence of s-box
layers, (xor)-linear diffusion, and xor-based key addition layers, altering
the differential operator yields a dual impact. On one hand, it is highly
probable that differentials traverse the s-box layer more effectively, given
that its non-linearity is maximised with respect to xor. On the other hand,
differentials do not deterministically propagate through the diffusion layer,
as observed in classical scenarios. This pivotal limitation effectively
restricts the success of the attack only to cases where the target layer is
linear not only concerning xor but also with respect to the operation under
consideration for computing differentials.

A first successful attempt based on the study of the alternative differential
properties of a xor-based toy cipher of the SPN family has shown that it is
possible to highlight a bias in the distribution of the differences calculated
compared to an alternative operation which is instead not detectable by means
of the standard xor-differential-based approach [CBS19]. The target cipher
featured five 3-bit s-boxes and the operation used to perform the attack acted
as the xor on the last four s-boxes, while on the first one matched with one
of the alternative sums defined by Calderini et al. [CCS21], coming from
another translation groups. The advantage of employing an alternative
operation in this case was only derived from the benefit induced by a single
s-box. In a more recent experimental approach [CCI24], we showed that better
results in a similar context can be obtained using an _alternative parallel
operation_ , in which every s-box can be targeted. In this case, the diffusion
layer of the cipher was determined through an algorithm, ensuring that it
adheres to the constraint of linearity with respect to both xor and the target
operation.

In this paper, we establish a general result that, in the context of an SPN
with 4-bit s-boxes, characterises all xor-linear maps that are concurrently
linear with respect to a parallel alternative operation (Sec. 2). This finding
enables the execution of a differential attack wherein each s-box affected by
a non-trivial differential contributes to the final differential probability
with increased efficacy compared to the conventional xor differentials.
Additionally, differentials propagate deterministically through the linear
layer in this scenario. Moreover, we examine all possible alternative
operations on 4 bits and investigate the differential properties of optimal
4-bit s-boxes, following the classification outlined by Leander and Poschmann
(a comparable methodology, albeit in the context of modular addition, was
recently employed by Zajac and Jókay [ZJ20]. Our analysis demonstrates that
each class comprises potentially weak permutations (Sec. 3). When coupled with
a diffusion layer as described earlier, these permutations have the potential
to render the cipher susceptible to differential attacks with alternative
operations. To substantiate our findings, we conclude the paper by presenting
experimental results on a family of toy SPNs (Sec. 4).

###  1.1. Notation

Let $V$ be an $n$-dimensional vector space over
$\operatorname{\mathbb{F}}_{2}$ which represents the message space. We write
$V=V_{1}\oplus V_{2}\oplus\cdots\oplus V_{b}$, where each $V_{j}$ is
isomorphic to a vector space $B$ such that $\dim(B)=s$ on which every s-box
acts. Thefore we have $n=sb$. We denote by $\\{e_{i}\\}_{i=1}^{n}$ the
canonical basis of $V$. If $G$ is any finite group acting on $V$, for each
$g\in G$ and $v\in V$ we denote the action of $g$ on $v$ as $vg$, i.e. we use
postfix notation for every function evaluation. We denote by
$\operatorname{Sym}(V)$ the symmetric group acting on $V$, i.e. the group of
all permutation on the message space, by $\operatorname{GL}(V,+)$ the group of
linear transformations, and by $\operatorname{AGL}(V,+)$ the group of affine
permutations. The identity matrix of size $l$ is denoted by $\mathbb{1}_{l}$
and the zero matrix of size $l\times h$ is denoted by $\mathbb{0}_{l,h}$, or
simply $\mathbb{0}_{l}$ if $l=h$. We finally denote by $T_{+}$ the group of
translations on $V$, i.e. $T_{+}\vcentcolon=\\{\sigma_{a}\mid a\in V,x\mapsto
x+a\\}<\operatorname{Sym}(V)$. We remind that the translation $\sigma_{k}$
acts on a vector $x$ in the same way the key-addition layer of an SPN acts
xor-ing the round key $k$ to the message $x$, i.e. $x\sigma_{k}=x+k$.

###  1.2. Preliminaries on alternative operations

An alternative operation on $V$ can be defined given any 2-elementary abelian
regular subgroup $T<\operatorname{AGL}(V,+)$, that we can write as
$T=\\{\tau_{a}\mid a\in V\\}$, where $\tau_{a}$ is the unique element in $T$
which maps $0$ into $a$. Consequently, for all $a,b\in$V, we can define
$a\circ b\vcentcolon=a\tau_{b}$, resulting in $(V,\circ)$ forming an additive
group. The operation $\circ$ induces a vector space structure on $V$, with the
corresponding group of translation being $T_{\circ}=T$. Additionally, for each
$a\in V$, there exists $M_{a}\in\operatorname{GL}(V,+)$ such that
$\tau_{a}=M_{a}\sigma_{a}$, meaning that for every $x\in V$,

$x\circ a=x\tau_{a}=xM_{a}+a.$

It is also assumed throughout that $T_{+}<\operatorname{AGL}(V,\circ)$, where
$\operatorname{AGL}(V,\circ)$ is the normaliser in $\operatorname{Sym}(V)$ of
$T_{\circ}$ (i.e., the group of affine permutations w.r.t. $\circ$). This
crucial technical assumption renders the key-addition layer an affine operator
concerning the new operation, enabling the prediction of how the key addition
affects the differentials with a reasonable probability. Further details on
this aspect, which may not be directly relevant to the scope of the current
paper, can be found in Civino et al. [CBS19]. In this context, we define the
_weak keys subspace_ as

$W_{\circ}\vcentcolon=\\{a\mid a\in V,\sigma_{a}=\tau_{a}\\}=\\{k\mid k\in V,\
\forall x\in V\ x\circ k=x+k\\}.$

$W_{\circ}$ is a vector subspace of both $(V,+)$ and $(V,\circ)$. It is known
[CDVS06, CCS21] that $W_{\circ}$ is non empty and that

$2-(n\bmod 2)\leq\dim(W_{\circ})\leq n-2.$ (1)

Moreover, up to conjugation we can always assume $W_{\circ}$ to be the span of
the last $d$ canonical vectors of $V$ [CCS21]. This allows to represent the
new sum in a canonical way [CCS21]: for each $a\in V$ there exists a matrix
$E_{a}\in\operatorname{\mathbb{F}}_{2}^{(n-d)\times d}$ such that

$M_{a}=\begin{pmatrix}\mathds{1}_{n-d}&E_{a}\\\
\mathbb{0}_{d,n-d}&\mathds{1}_{d}\end{pmatrix}.$ (2)

Fixing such an operation as above is therefore equivalent to defining the
matrices

$M_{e_{i}}=\begin{pmatrix}\mathds{1}_{n-d}&E_{e_{i}}\\\
\mathbb{0}_{d,n-d}&\mathds{1}_{d}\end{pmatrix}=\left(\begin{array}[]{@{}c|c@{}}\mathds{1}_{n-d}&\begin{matrix}\textbf{b}_{i,1}\\\
\vdots\\\ \textbf{b}_{i,n-d}\end{matrix}\\\
\hline\cr\mathbb{0}_{d,n-d}&\mathds{1}_{d}\end{array}\right)$

for $1\leq i\leq n$, where
$\textbf{b}_{i,j}\in\operatorname{\mathbb{F}}_{2}^{d}$. The assumptions on
$T_{\circ}$ and on $W_{\circ}$ imply that $E_{e_{i}}=0$ for $n-d+1\leq i\leq
n$, $\textbf{b}_{i,i}=\textbf{0}$ and $\textbf{b}_{i,j}=\textbf{b}_{j,i}$. In
conclusion, the following result characterises the criteria that the vectors
$\textbf{b}_{i,j}$ must adhere to in order to define an alternative operation
as previously described.

######  Theorem 1.1 ([CBS19]).

Let $T_{\circ}<\operatorname{AGL}(V,+)$ be 2-elementary, abelian, and regular,
and let $d\leq n-2$. The operation $\circ$ induced by $T_{\circ}$ is such that
$d=\dim(W_{\circ})$, $T_{+}<\operatorname{AGL}(V,\circ)$, and
$W_{\circ}=\operatorname{Span}\\{e_{n-d+1},\dots,e_{n}\\}$ if and only if the
matrix
$\Theta_{\circ}\in(\operatorname{\mathbb{F}}_{2^{d}})^{(n-d)\times(n-d)}$
defined as

$\Theta_{\circ}\vcentcolon=\begin{pmatrix}\mathbf{b}_{1,1}&\mathbf{b}_{1,2}&\cdots&\mathbf{b}_{n-d,1}\\\
\mathbf{b}_{2,1}&\mathbf{b}_{2,2}&\cdots&\mathbf{b}_{n-d,2}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{b}_{n-d,1}&\mathbf{b}_{n-d,2}&\cdots&\mathbf{b}_{n-d,d}\end{pmatrix}$

is zero-diagonal, symmetric and no $\operatorname{\mathbb{F}}_{2}$-linear
combination of its columns is the null vector. The matrix $\Theta_{\circ}$ is
also called the _defining matrix_ for $\circ$.

In the subsequent discussion, the term _alternative operation_ refers to an
additive law $\circ$ on $V$ as defined above.

##  2\. Parallel operations and their automorphism groups

Let $\circ$ be an alternative operation on the block-sized space $V$. As
outlined in the introduction, if $\lambda\in\operatorname{GL}(V,+)$ represents
a (xor)-linear diffusion layer, and $\Delta\in V$ is an input difference
traversing $\lambda$, predicting the output difference with respect to
$\circ$, i.e.,

$x\lambda\circ(x\circ\Delta)\lambda,$

becomes inherently challenging without additional assumptions on $\lambda$
that ensure a sufficiently high predictive probability. For this reason, the
examination of the following object becomes crucial: in cryptographic terms,
it contains potential diffusion layers that allow differentials, whether
computed with respect to xor or $\circ$, to propagate with a probability of 1.

######  Definition 2.1.

Let $\circ$ be an alternative operation on $V$. Let us define

$H_{\circ}:=\\{f\in\operatorname{GL}(V,+)\mid\forall a,b\in V:\,\,\,(a\circ
b)f=af\circ bf\\}$

to be the subgroup of $\operatorname{GL}(V,+)$ of permutations that are linear
w.r.t. the operation $\circ$. More precisely, denoting by
$\operatorname{AGL}(V,\circ)$ the normaliser in $\operatorname{Sym}(V)$ of
$T_{\circ}$ and by $\operatorname{GL}(V,\circ)$ the stabiliser of $0$ in
$\operatorname{AGL}(V,\circ)$, we have $H_{\circ}$ =
$\operatorname{GL}(V,+)\cap\operatorname{GL}(V,\circ)$.

The structure of the group $H_{\circ}$ in its most general case has not been
understood yet. This work addresses this challenge in a specific scenario,
guided by assumptions that are deemed reasonable within the context of
differential cryptanalysis.

#### Assumption 1: $\circ$ is a parallel operation.

While the operation $\circ$ could, in theory, be defined on the entire message
space $V$, studying the differential properties of the s-box layer, considered
as a function with $2^{n}$ inputs, is impractical for standard-size ciphers.
For this reason, we focus on operations applied in a _parallel_ way to each
s-box-sized block, i.e., $\circ=(\circ_{1},\circ_{2},\dots,\circ_{b})$, where
for each $1\leq j\leq b$, $\circ_{j}$ is an operation on $V_{j}$. In this
scenario, every operation is acting independently on the s-box space $B$,
regardless of the others. This motivates the following definition.

######  Definition 2.2.

Let $\circ$ be an alternative operation on $V$. We say that $\circ$ is
_parallel_ if for each $1\leq j\leq b$ there exists an alternative operation
$\circ_{j}$ on $V_{j}$ such that for each $x,y\in V$ we have

$x\circ y=\begin{pmatrix}x_{1}\\\ \vdots\\\
x_{b}\end{pmatrix}\circ\begin{pmatrix}y_{1}\\\ \vdots\\\
y_{b}\end{pmatrix}=\begin{pmatrix}x_{1}\circ_{1}y_{1}\\\ \vdots\\\
x_{b}\circ_{b}y_{b}\end{pmatrix},$

where $x=(x_{1},x_{2},\dots,x_{b}),y=(y_{1},y_{2},\dots,y_{b})$ and each
component belongs to the s-box-sized space, i.e., $x_{j},y_{j}\in V_{j}\cong
B$ for $1\leq j\leq b$.

In the notation of Sec. 1.2, up to a block matrix conjugation, we can assume
that every element $x\in V$ is associated to a translation
$\tau_{x}=M_{x}\sigma_{x}$, with

$M_{x}=\begin{pmatrix}M_{x_{1}}^{\circ_{1}}&\cdots&\mathbb{0}\\\
\vdots&\ddots&\vdots\\\ \mathbb{0}&\cdots&M_{x_{b}}^{\circ_{b}}\end{pmatrix}$

where $M_{x_{i}}^{\circ_{i}}$ is the matrix associated to the translation
$\tau_{x_{i}}$ with respect to the sum $\circ_{i}$, as defined in Eq. (2).
Notice that it can be assumed, without loss of generality, that all the
operations $\circ_{j}$ coincide.

#### Assumption 2: $\dim(W_{\circ_{j}})=s-2$.

According to Eq. (1), every operation $\circ_{j}$ defined at the s-box level
must satisfy the bound $\dim(W_{\circ_{j}})\leq s-2$, being $s=\dim(B)$. The
situation where the (upper) bound is reached holds particular interest for
several reasons, as elaborated further in Civino et al. [CBS19]. Notably,

  * •

if the s-box size $s$ is four, the case where $\dim(W_{\circ_{j}})=2$ is the
sole possibility;

  * •

the considered case stands today as the only one for which the structure of
$H_{\circ_{j}}$ is well understood.

For the reader’s convenience, we present the classification result for
$H_{\circ_{j}}$ obtained by Civino et al. in the considered case.
Additionally, it is worth recalling that, according to Theorem 1.1, any
$\circ_{j}$ for which $\dim(W_{\circ_{j}})=s-2$ is determined by a single non-
null vector $\mathbf{b}\in(\operatorname{\mathbb{F}}_{2})^{s-2}$.

######  Theorem 2.3 ([CBS19]).

Let $\circ_{j}$ be an alternative operation such that
$d=\dim(W_{\circ_{j}})=s-2$ defined by a vector
$\textbf{b}\in(\operatorname{\mathbb{F}}_{2})^{s-2}$, and let
$\lambda\in(\operatorname{\mathbb{F}}_{2})^{s\times s}$. The following are
equivalent:

  * •

$\lambda\in H_{\circ_{j}}$;

  * •

there exist $A\in\operatorname{GL}((\operatorname{\mathbb{F}}_{2})^{2},+)$,
$D\in\operatorname{GL}((\operatorname{\mathbb{F}}_{2})^{d},+)$, and
$B\in(\operatorname{\mathbb{F}}_{2})^{2\times d}$ such that

$\lambda=\begin{pmatrix}A&B\\\ \mathbb{0}_{d,2}&D\end{pmatrix}$

and $\textbf{b}D=\textbf{b}$.

We are now prepared to present the first novel contribution of this work,
wherein we characterise the group $H_{\circ}$ for a parallel operation
$\circ=(\circ_{1},\circ_{2},\dots,\circ_{b})$ with components at the s-box
level satisfying $\dim(W_{\circ_{j}})=s-2$. For the sake of simplicity and
without losing generaly, we assume that the $b$ operations at the s-box level
coincide.

######  Theorem 2.4.

Let $\circ=(\circ_{1},\circ_{2},\dots,\circ_{b})$ be a parallel alternative
operation on $V$ such that for each $1\leq j\leq b$ $\circ_{j}$ is an
alternative operation on $V_{j}$. Let us assume that every $\circ_{j}$ is such
that $\dim(W_{\circ_{j}})=s-2$ and it is defined by a vector
$\textbf{b}\in(\operatorname{\mathbb{F}}_{2})^{s-2}$. Let
$\lambda\in(\operatorname{\mathbb{F}}_{2})^{n\times n}$. Then, $\lambda\in
H_{\circ}$ if and only if it can be represented in the block form

$\lambda=\left(\begin{array}[]{@{}c|c|c@{}}\begin{matrix}A_{11}&B_{11}\\\
C_{11}&D_{11}\end{matrix}&\cdots&\begin{matrix}A_{1b}&B_{1b}\\\
C_{1b}&D_{1b}\end{matrix}\\\ \hline\cr\vdots&\ddots&\vdots\\\
\hline\cr\begin{matrix}A_{b1}&B_{b1}\\\
C_{b1}&D_{b1}\end{matrix}&\cdots&\begin{matrix}A_{bb}&B_{bb}\\\
C_{bb}&D_{bb}\end{matrix}\\\ \end{array}\right),$

where

  1. (1)

$A_{ij}\in(\operatorname{\mathbb{F}}_{2})^{2\times 2}$ such that for each row
and each column of blocks there exists one and only one non-zero $A_{ij}$;
moreover, all the non-zero $A_{ij}$ are invertible;

  2. (2)

$B_{ij}\in(\operatorname{\mathbb{F}}_{2})^{2\times(s-2)}$;

  3. (3)

$C_{ij}=\mathbb{0}_{(s-2)\times 2}$;

  4. (4)

$D_{ij}\in(\operatorname{\mathbb{F}}_{2})^{(s-2)\times(s-2)}$ such that if
$A_{ij}$ is zero, then $\textbf{b}D_{ij}=\textbf{0}$, and if $A_{ij}$ is
invertible, then $\textbf{b}D_{ij}=\textbf{b}$. Moreover, the matrix $D$
defined by

$D\vcentcolon=\begin{pmatrix}D_{11}&\cdots&D_{1b}\\\ \vdots&\ddots&\vdots\\\
D_{b1}&\cdots&D_{bb}\end{pmatrix}$

is invertible.

###### Proof.

The proof involves standard linear algebra techniques, but its extensive and
laborious nature necessitates omission due to page limitations. ∎

##  3\. Differential properties of optimal s-boxes

In this section we delve into the examination of the differential properties
exhibited by all possible 4-bit permutations, with respect to all possible
alternative operations defined as in Sec. 1.2. In particular, we set $s=4$ and
therefore consider $B=\mathbb{F}_{2}^{4}$. We begin by acknowledging that,
despite the compact size of the space, the count of alternative operations on
$B$ is considerable:

######  Proposition 3.1 ([CCS21]).

There exist 105 different elementary abelian regular subgroups groups
$T_{\circ}$ in $AGL(\mathbb{F}_{2}^{4},+)$. Furthermore, each of them
satisfies $T_{+}<AGL(\mathbb{F}_{2}^{4},\circ)$ and $\dim W_{\circ}=s-2=2$.

We recall that given a permutation $f\in\operatorname{Sym}(B)$ we can define

$\delta_{f}(a,b)=\\#\\{x\in B\mid xf+(x+a)f=b\\}.$

The _differential uniformity_ of $f$ is defined as $\delta_{f}:=\max_{a\neq
0}\delta_{f}(a,b)$ and it represent the primary metric to consider when
assessing the resistance of an s-box to differential cryptanalysis [Nyb93].

Several cryptographic properties, including differential uniformity, are
preserved under affine equivalence for vectorial Boolean functions. Two
functions, denoted as $f$ and $g$, are considered _affine equivalent_ if there
exist two affine permutations, $\alpha$ and $\beta$, in
$\operatorname{AGL}(V,+)$ such that $g=\beta f\alpha$.

Leander and Poschmann [LP07] provided a comprehensive classification (up to
affine equivalence) of permutations over $B=\mathbb{F}_{2}^{4}$. They
identified 16 classes with _optimal_ cryptographic properties. All 16 classes
exhibit a classical differential uniformity equal to 4, which represents the
best possible value for s-boxes in $\operatorname{Sym}(B)$. The
representatives of the 16 classes are listed in Table 1, where each vector is
interpreted as a binary number, most significant bit first.

Table 1. Optimal 4-bit permutations according to Leander and Poschmann | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | 3${}_{\text{x}}$ | 4${}_{\text{x}}$ | 5${}_{\text{x}}$ | 6${}_{\text{x}}$ | 7${}_{\text{x}}$ | 8${}_{\text{x}}$ | 9${}_{\text{x}}$ | A${}_{\text{x}}$ | B${}_{\text{x}}$ | C${}_{\text{x}}$ | D${}_{\text{x}}$ | E${}_{\text{x}}$ | F${}_{\text{x}}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$G_{0}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | B${}_{\text{x}}$ | C${}_{\text{x}}$ | 9${}_{\text{x}}$ | 3${}_{\text{x}}$ | E${}_{\text{x}}$ | A${}_{\text{x}}$ | 5${}_{\text{x}}$
$G_{1}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | B${}_{\text{x}}$ | E${}_{\text{x}}$ | 3${}_{\text{x}}$ | 5${}_{\text{x}}$ | 9${}_{\text{x}}$ | A${}_{\text{x}}$ | 12${}_{\text{x}}$
$G_{2}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | B${}_{\text{x}}$ | E${}_{\text{x}}$ | 3${}_{\text{x}}$ | A${}_{\text{x}}$ | C${}_{\text{x}}$ | 5${}_{\text{x}}$ | 9${}_{\text{x}}$
$G_{3}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | C${}_{\text{x}}$ | 5${}_{\text{x}}$ | 3${}_{\text{x}}$ | A${}_{\text{x}}$ | E${}_{\text{x}}$ | B${}_{\text{x}}$ | 9${}_{\text{x}}$
$G_{4}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | C${}_{\text{x}}$ | 9${}_{\text{x}}$ | B${}_{\text{x}}$ | A${}_{\text{x}}$ | E${}_{\text{x}}$ | 5${}_{\text{x}}$ | 3${}_{\text{x}}$
$G_{5}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | C${}_{\text{x}}$ | B${}_{\text{x}}$ | 9${}_{\text{x}}$ | A${}_{\text{x}}$ | E${}_{\text{x}}$ | 3${}_{\text{x}}$ | 5${}_{\text{x}}$
$G_{6}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | C${}_{\text{x}}$ | B${}_{\text{x}}$ | 9${}_{\text{x}}$ | A${}_{\text{x}}$ | E${}_{\text{x}}$ | 5${}_{\text{x}}$ | 3${}_{\text{x}}$
$G_{7}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | C${}_{\text{x}}$ | E${}_{\text{x}}$ | B${}_{\text{x}}$ | A${}_{\text{x}}$ | 9${}_{\text{x}}$ | 3${}_{\text{x}}$ | 5${}_{\text{x}}$
$G_{8}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | E${}_{\text{x}}$ | 9${}_{\text{x}}$ | 5${}_{\text{x}}$ | A${}_{\text{x}}$ | B${}_{\text{x}}$ | 3${}_{\text{x}}$ | 12${}_{\text{x}}$
$G_{9}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | E${}_{\text{x}}$ | B${}_{\text{x}}$ | 3${}_{\text{x}}$ | 5${}_{\text{x}}$ | 9${}_{\text{x}}$ | A${}_{\text{x}}$ | 12${}_{\text{x}}$
$G_{10}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | E${}_{\text{x}}$ | B${}_{\text{x}}$ | 5${}_{\text{x}}$ | A${}_{\text{x}}$ | 9${}_{\text{x}}$ | 3${}_{\text{x}}$ | 12${}_{\text{x}}$
$G_{11}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | E${}_{\text{x}}$ | B${}_{\text{x}}$ | A${}_{\text{x}}$ | 5${}_{\text{x}}$ | 9${}_{\text{x}}$ | C${}_{\text{x}}$ | 3${}_{\text{x}}$
$G_{12}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | E${}_{\text{x}}$ | B${}_{\text{x}}$ | A${}_{\text{x}}$ | 9${}_{\text{x}}$ | 3${}_{\text{x}}$ | C${}_{\text{x}}$ | 5${}_{\text{x}}$
$G_{13}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | E${}_{\text{x}}$ | C${}_{\text{x}}$ | 9${}_{\text{x}}$ | 5${}_{\text{x}}$ | B${}_{\text{x}}$ | A${}_{\text{x}}$ | 3${}_{\text{x}}$
$G_{14}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | E${}_{\text{x}}$ | C${}_{\text{x}}$ | B${}_{\text{x}}$ | 3${}_{\text{x}}$ | 9${}_{\text{x}}$ | 5${}_{\text{x}}$ | 10${}_{\text{x}}$
$G_{15}$ | 0${}_{\text{x}}$ | 1${}_{\text{x}}$ | 2${}_{\text{x}}$ | D${}_{\text{x}}$ | 4${}_{\text{x}}$ | 7${}_{\text{x}}$ | F${}_{\text{x}}$ | 6${}_{\text{x}}$ | 8${}_{\text{x}}$ | E${}_{\text{x}}$ | C${}_{\text{x}}$ | B${}_{\text{x}}$ | 9${}_{\text{x}}$ | 3${}_{\text{x}}$ | A${}_{\text{x}}$ | 5${}_{\text{x}}$

###  3.1. Dealing with affine equivalence

Our goal is to analyse the differential uniformity of each optimal s-box
class, with respect to every alternative operation $\circ$ on $B$. The
definitions given above can be generalised in the obvious way setting
$\delta^{\circ}_{f}(a,b)=\\#\\{x\in B\mid xf\circ(x\circ a)f=b\\}$ and calling
_$\circ$ -differential uniformity_ of $f$ the value
$\delta^{\circ}_{f}:=\max_{a\neq 0}\delta^{\circ}_{f}(a,b)$.

It is noteworthy that, unlike in the case of classic differential uniformity,
the value of $\delta^{\circ}_{f}$ is not invariant under affine equivalence.
However, verifying the $\circ$-differential uniformity of $g_{2}G_{i}g_{1}$
for any optimal class and every pair $g_{1},g_{2}\in\operatorname{AGL}(V,+)$
would be impractical. Therefore, a reduction in the number of permutations to
be checked is necessary, and for this purpose, we make the following
observations. First, similar to the classical case, the $\circ$-differential
uniformity is preserved under affine transformations w.r.t. $\circ$.

######  Proposition 3.2.

Given $f\in\operatorname{Sym}(B)$ and
$g_{1},g_{2}\in\operatorname{AGL}(B,\circ)$ we have

$\delta_{g_{1}fg_{2}}^{\circ}(a,b)=\delta_{f}^{\circ}(g_{2}(a),g_{1}^{-1}(b)).$

Moreover, Proposition 3.1 establishes that for any $\circ$ derived from a
translation group in $\operatorname{AGL}(B,+)$, the $+$-translations are
affine with respect to $\circ$. This initial observation allows us to narrow
down the analysis to $g_{2}G_{i}g_{1}$ with
$g_{1},g_{2}\in\operatorname{GL}(B,+)$, which still remains impractical.
Furthermore, considering that
$H_{\circ}=\operatorname{GL}(B,+)\cap\operatorname{GL}(B,\circ)$, Proposition
3.2 establishes that left and right multiplication by elements in $H_{\circ}$
preserves both $\circ$ and $+$-differential uniformity. It is noteworthy that
during this process, the rows of the matrix containing all the
$\delta^{\circ}_{f}(a,b)$ ($\mathrm{DDT}^{\circ}$) are merely shuffled,
thereby preserving the highest element of each row. Therefore, the following
conclusion can be easily obtained.

######  Proposition 3.3.

Let $g_{1},g_{2}\in\operatorname{GL}(B,+)$ and $f\in\operatorname{Sym}(B)$.
For any $g_{1}^{\prime}\in g_{1}H_{\circ}$ and $g_{2}^{\prime}\in
H_{\circ}g_{2}$ we have

$\delta^{\circ}_{g_{2}fg_{1}}=\delta^{\circ}_{g_{2}^{\prime}fg_{1}^{\prime}}.$

###### Proof.

Take $h_{1},h_{2}\in H_{\circ}$ such that $g^{\prime}_{1}=g_{1}h_{1}$ and
$g_{2}^{\prime}=h_{2}g_{2}$. Then,

$xg^{\prime}_{2}fg_{1}^{\prime}\circ(x\circ
a)g^{\prime}_{2}fg_{1}^{\prime}=xh_{2}g_{2}fg_{1}\circ(xh_{2}\circ
ah_{2}g_{2}fg_{1})h_{1},$

implying that
$\delta^{\circ}_{g^{\prime}_{2}fg_{1}^{\prime}}(a,b)=\delta^{\circ}_{g_{2}fg_{1}}(ah_{2},bh_{1}^{-1})$.
So,
$\delta^{\circ}_{g^{\prime}_{2}fg_{1}^{\prime}}=\delta^{\circ}_{g_{2}fg_{1}}$.∎∎

The final proposition allows us to focus solely on $g_{1}$ and $g_{2}$ within
the left and right cosets of $H_{\circ}$. These reductions facilitate the
analysis of the potential $\circ$-differential uniformities attainable across
all classes of optimal permutations for the 105 conceivable alternative sums
defined over $B$. For each of the 105 alternative operations, we
systematically explored each of the 16 classes, following the described
procedure, and we recorded the $\circ$-differential uniformity for every
candidate. To streamline the presentation, we calculated the average across
the 105 operations and presented the consolidated results in Tab. 2.

Table 2. Avg. number of functions with given $\circ$-differential uniformity Class $\delta^{\circ}$ | 2 | 4 | 6 | 8 | 10 | 12 | 14 | 16
---|---|---|---|---|---|---|---|---
$G_{0}$ | 0 | 914 | 7842 | 3463 | 420 | 19 | 0 | 14
$G_{1}$ | 0 | 1019 | 10352 | 4226 | 560 | 0 | 0 | 18
$G_{2}$ | 0 | 1003 | 8604 | 3805 | 462 | 21 | 0 | 16
$G_{3}$ | 0 | 16733 | 117740 | 27639 | 1779 | 0 | 0 | 0
$G_{4}$ | 0 | 1101 | 9295 | 2715 | 179 | 0 | 0 | 0
$G_{5}$ | 0 | 2479 | 24135 | 5402 | 639 | 0 | 0 | 0
$G_{6}$ | 0 | 1632 | 10842 | 3071 | 218 | 0 | 0 | 0
$G_{7}$ | 0 | 1257 | 10679 | 2994 | 186 | 28 | 0 | 0
$G_{8}$ | 0 | 1691 | 12821 | 6113 | 583 | 93 | 0 | 24
$G_{9}$ | 0 | 1228 | 7734 | 2693 | 154 | 39 | 0 | 0
$G_{10}$ | 0 | 1228 | 8063 | 2763 | 166 | 41 | 0 | 0
$G_{11}$ | 0 | 1637 | 9940 | 2941 | 214 | 0 | 0 | 0
$G_{12}$ | 0 | 2541 | 16832 | 5308 | 352 | 0 | 0 | 0
$G_{13}$ | 0 | 1124 | 9520 | 2416 | 217 | 15 | 0 | 0
$G_{14}$ | 0 | 1207 | 7641 | 2584 | 160 | 51 | 0 | 0
$G_{15}$ | 0 | 1227 | 7776 | 2630 | 163 | 52 | 0 | 0

In our examination, we observe that if, for a given operation $\circ$, certain
elements within an affine equivalence class yield a $\circ$-differential
uniformity $\delta$, then this value $\delta$ is achieved by some element in
the entire class for all alternative operations. Our analysis reveals that
certain optimal functions may exhibit the highest differential uniformity (16)
for alternative operations, specifically the classes $G_{0}$ (containing,
e.g., the s-box S1 of Serpent [BAK98]), $G_{1}$ (containing, e.g., the s-box
of Present [BKL+07]), $G_{2}$, and $G_{8}$. Conversely, the classes $G_{3}$,
$G_{4}$, $G_{5}$, $G_{6}$, $G_{11}$, and $G_{12}$ demonstrate more favorable
behavior concerning alternative operations.

##  4\. Experiments on a 16-bit block cipher with 4-bit s-boxes

In this concluding section, we aim to apply the results obtained above to a
family of (toy) ciphers. These ciphers may exhibit security under classical
differential cryptanalysis but reveal vulnerabilities to the alternative
differential approach.

In our experiments, we set $V=\operatorname{\mathbb{F}}_{2}^{16}$, $n=4$, and
$s=4$, defining $\circ$ as the parallel sum by applying the alternative
operation defined by the vector $\mathbf{b}=(0,1)$ to each $4$-bit block.
Moreover, all our ciphers will feature the 4-bit permutation
$\gamma:\operatorname{\mathbb{F}}_{2}^{4}\rightarrow\operatorname{\mathbb{F}}_{2}^{4}$
defined by the sequence (0${}_{\text{x}}$, E${}_{\text{x}}$, B${}_{\text{x}}$,
1${}_{\text{x}}$, 7${}_{\text{x}}$, C${}_{\text{x}}$, 9${}_{\text{x}}$,
6${}_{\text{x}}$, D${}_{\text{x}}$, 3${}_{\text{x}}$, 4${}_{\text{x}}$,
F${}_{\text{x}}$, 2${}_{\text{x}}$, 8${}_{\text{x}}$, A${}_{\text{x}}$,
5${}_{\text{x}}$) as its s-box. Precisely, four copies of $\gamma$ will act on
the 16-bit block. Notice that the s-box $\gamma\in$ belongs to $G_{0}$ and has
$\delta_{\gamma}=4$ and $\delta^{\circ}_{\gamma}=16$.

In all the experiments described below, we consider the SPN whose $i$-th round
is obtained by the composition of the parallel application of the s-box
$\gamma$ on every 4-bit block, a ‘diffusion layer’ $\lambda$ sampled random
from $H_{\circ}$, and the xor with the $i$-th random round key. We study the
difference propagation in the cipher in a long-key scenario, i.e., the key-
schedule selects a random long key $k\in\mathbb{F}_{2}^{16r}$ where $r$ is the
number of rounds. To avoid potential bias from a specific key choice, we
conduct our experiments by averaging over $2^{15}$ random long-key
generations. This approach gives us a reliable estimate of the expected
differential probability for the best differentials in this cipher.

In 150 distinct executions, spanning a range of rounds from 3 to 10, we
calculated the discrepancy between the most effective $\circ$-trail and
$+$-trail. To manage computational resources, our focus was narrowed down to
input differences with a Hamming weight of 1.

The results are depicted in Fig. 1, where each dot represents an individual
simulation. The $x$ axis corresponds to the negative logarithm of the
probability of the best $\circ$ differential, while the $y$ axis represents
the difference between that value and the negative logarithm of the
probability of the best $+$ differential. Darker dots indicate a higher number
of rounds, as explained in the legend. Notably, about half of the dots lie
above zero, suggesting that the best $\circ$ differential consistently
outperforms the best $+$ differential until they become indistinguishable.
Interestingly, this convergence often occurs when the $\circ$ probability is
already very close to 16, providing potential candidates for our distinguisher
attack.

Figure 1. Comparison of $+$ and $\circ$ trails for random mixing layers

## References

  * [BAK98] Eli Biham, Ross Anderson, and Lars Knudsen, _Serpent: A new block cipher proposal_ , International workshop on fast software encryption, Springer, 1998, pp. 222–238.
  * [BBI+15] Subhadeep Banik, Andrey Bogdanov, Takanori Isobe, Kyoji Shibutani, Harunaga Hiwatari, Toru Akishita, and Francesco Regazzoni, _Midori: A block cipher for low energy_ , Advances in Cryptology–ASIACRYPT 2015: 21st International Conference on the Theory and Application of Cryptology and Information Security, Auckland, New Zealand, November 29–December 3, 2015, Proceedings, Part II 21, Springer, 2015, pp. 411–436.
  * [BBS05] Eli Biham, Alex Biryukov, and Adi Shamir, _Cryptanalysis of Skipjack reduced to 31 rounds using impossible differentials_ , Journal of Cryptology 18 (2005), 291–311.
  * [BCJW02] Nikita Borisov, Monica Chew, Rob Johnson, and David Wagner, _Multiplicative differentials_ , Fast Software Encryption: 9th International Workshop, FSE 2002 Leuven, Belgium, February 4–6, 2002 Revised Papers 9, Springer, 2002, pp. 17–33.
  * [BKL+07] Andrey Bogdanov, Lars R Knudsen, Gregor Leander, Christof Paar, Axel Poschmann, Matthew JB Robshaw, Yannick Seurin, and Charlotte Vikkelsoe, _Present: An ultra-lightweight block cipher_ , Cryptographic Hardware and Embedded Systems-CHES 2007: 9th International Workshop, Vienna, Austria, September 10-13, 2007. Proceedings 9, Springer, 2007, pp. 450–466.
  * [BS91] Eli Biham and Adi Shamir, _Differential cryptanalysis of DES-like cryptosystems_ , Journal of CRYPTOLOGY 4 (1991), 3–72.
  * [CBS19] Roberto Civino, Céline Blondeau, and Massimiliano Sala, _Differential attacks: using alternative operations_ , Designs, Codes and Cryptography 87 (2019), 225–247.
  * [CCI24] Marco Calderini, Roberto Civino, and Riccardo Invernizzi, _Differential experiments using parallel alternative operations_ , Journal of Mathematical Cryptology 18 (2024), no. 1, 20230030.
  * [CCS21] Marco Calderini, Roberto Civino, and Massimiliano Sala, _On properties of translation groups in the affine general linear group with applications to cryptography_ , Journal of Algebra 569 (2021), 658–680.
  * [CDVS06] A. Caranti, F. Dalla Volta, and M. Sala, _Abelian regular subgroups of the affine group and radical rings_ , Publ. Math. Debrecen 69 (2006), no. 3, 297–308.
  * [Knu95] Lars R Knudsen, _Truncated and higher order differentials_ , Fast Software Encryption: Second International Workshop Leuven, Belgium, December 14–16, 1994 Proceedings 2, Springer, 1995, pp. 196–211.
  * [LP07] Gregor Leander and Axel Poschmann, _On the classification of 4 bit s-boxes_ , Arithmetic of Finite Fields: First International Workshop, WAIFI 2007, Madrid, Spain, June 21-22, 2007. Proceedings 1, Springer, 2007, pp. 159–176.
  * [Nyb93] Kaisa Nyberg, _Differentially uniform mappings for cryptography_ , Workshop on the Theory and Application of of Cryptographic Techniques, Springer, 1993, pp. 55–64.
  * [SIH+11] Kyoji Shibutani, Takanori Isobe, Harunaga Hiwatari, Atsushi Mitsuda, Toru Akishita, and Taizo Shirai, _Piccolo: an ultra-lightweight blockcipher_ , Cryptographic Hardware and Embedded Systems–CHES 2011: 13th International Workshop, Nara, Japan, September 28–October 1, 2011. Proceedings 13, Springer, 2011, pp. 342–357.
  * [Teş22] George Teşeleanu, _The security of quasigroups based substitution permutation networks_ , International Conference on Information Technology and Communications Security, Springer, 2022, pp. 306–319.
  * [Wag99] David Wagner, _The boomerang attack_ , International Workshop on Fast Software Encryption, Springer, 1999, pp. 156–170.
  * [ZJ20] Pavol Zajac and Matúš Jókay, _Cryptographic properties of small bijective s-boxes with respect to modular addition_ , Cryptography and Communications 12 (2020), 947–963.


# Particle-scale modeling of the drying characteristics of colloidal
suspensions

Rei Tatsumi<EMAIL_ADDRESS>Environmental Science Center, The
University of Tokyo, Tokyo 113-8656, Japan  Osamu Koike  Yukio Yamaguchi
Product Innovation Association, Tokyo 113-8656, Japan  Yoshiko Tsuji1,
Department of Chemical System Engineering, The University of Tokyo, Bunkyo,
Tokyo 113-8656, Japan

###### Abstract

During drying of colloidal suspensions, colloidal particles can form
concentrated particle layers beneath the receding free surface. The drying
rate can gradually decrease with the growth of the particle layers. We
construct a model to investigate how such drying characteristics is affected
by interactions between particles. In this model, the formation of the
particle layers is described by Langevin dynamics simulations, and the drying
rate is evaluated from the permeation resistance of the particle layers. We
show that the decrease in the drying rate is suppressed when the particles
form aggregates by attractive interactions. The present model would enable us
to predict and control the drying characteristics through the character of
colloidal particles.

##  I Introduction

Colloidal suspensions are coated on substrates and upon drying produce various
functional materials, whose quality is determined by the microstructure
composed of the colloidal particles. During drying of colloidal suspensions,
the receding free surface (the liquid–air interface) induces the particles to
form structures, thereby causing a decrease in the drying rate. dryc1 ; dryc2
; dryc3  The prediction and control of such drying characteristics is
therefore required to improve material quality as well as to reduce drying
time.

Figure 1 schematically illustrates the structure formation corresponding to
the drying characteristics of colloidal suspensions. During the constant rate
period, the drying rate is controlled by vapor diffusion from the free surface
where the liquid evaporates. The receding free surface sweeps the particles to
form concentrated particle layers, and then the falling rate period begins.
The condition that this structure formation occurs can be quantified by the
particle drying Péclet number as
$\mathrm{Pe}=\tau_{\mathrm{D}}/\tau_{\mathrm{E}}>1$, where $\tau_{\mathrm{D}}$
and $\tau_{\mathrm{E}}$ are the time scales of particle diffusion by the
Brownian motion and the recession of the free surface, respectively. clyr1 ;
clyr2 ; sdry1  The continuous evaporation requires liquid transport through
the pores of the particle layer to the free surface. With the growth of the
particle layer, the permeation resistance of the particle layer increases. The
drying rate is thus controlled by the liquid transport through the particle
layer and decreases gradually. Such decrease in the drying rate has been
quantitatively observed in experiments of the unidirectional drying of
colloidal suspensions confined in a thin rectangular cell. In these
experiments, the drying rate can be evaluated by observing the decrease in the
volume of the suspension in the drying cell. udry1 ; udry2 ; udry3 ; udry4 ;
udry5 ; udry6  When the particle layer recedes with the free surface and grows
to reach the substrate with further drying, the free surface goes into the
particle layer. The drying rate is then controlled by vapor diffusion through
the particle layer and will further decrease, as observed for granular
materials. In this stage, liquid flow can be driven by a capillary pressure
due to menisci formed among particles. This liquid transport would suppress
the decrease in the drying rate. pors1 ; pors2

We focus on the decrease in the drying rate due to the growth of concentrated
particle layers. The character of colloidal particles would affect the
structure of the particle layers. For example, in the drying of droplets of
colloidal suspensions, the morphology of the dried grains changes depending on
interactions between the particles. sdry2 ; sdry3  As for the drying of binary
colloidal mixtures containing particles of different sizes, segregation of
smaller particles to the top surface can occur depending on the mixing ratio,
the particle size ratio, and drying rates. seg1 ; seg2 ; seg3  Such structural
differences would be reflected in the drying characteristics.

In this study, we construct a model to investigate how the drying
characteristics of colloidal suspensions is affected by interactions between
the colloidal particles. We first derive an analytical expression of the
relationship between the drying rate and the permeation resistance of
concentrated particle layers. We then perform Langevin dynamics simulations
that describe the Brownian motion of colloidal particles by stochastic
differential equations. The receding free surface is also considered in the
simulations to induce the formation of concentrated particle layers. During
the simulations, the drying rate is evaluated from the permeation resistance
of the particle layers by use of the derived analytical expression. We assume
the DLVO potential between charged colloidal particles and tune the potential
from repulsive to attractive by varying ionic strength. We show that the
decrease in the drying rate is suppressed when the potential energy barrier is
so low that the particles form aggregates. We discuss how the interaction
between particles affect the drying characteristics in terms of the structural
differences of the particle layers.

Figure 1:  Sketch of the drying characteristics of colloidal suspensions.

##  II Model

###  II.1 Drying rates

Figure 2:  Drying rate as a function of permeation resistance at different
humidities. The profiles are given by Eq. (4) for $\varphi<1$ and Eq. (11) for
$\varphi\rightarrow 1$.  Figure 3:  Characteristic thickness of concentrated
particle layers as a function of particle diameter. The values scaled by $1\
\mathrm{mm}$ and the particle diameter $d$ are exhibited in the left and right
axes, respectively. The characteristic thickness
$L_{0}=R_{0}/\rho_{\mathrm{c}}$ is calculated by Eqs. (7) and (12) with
setting the parameters: $\eta=1\times 10^{-3}\ \mathrm{Pa\ s}$, $T=300\
\mathrm{K}$, $V_{\mathrm{m}}=3\times 10^{-29}\ \mathrm{m^{3}}$,
$p_{0}=3.5\times 10^{3}\ \mathrm{Pa}$, $K=8\times 10^{-2}\ \mathrm{m/s}$, and
$\Phi_{\mathrm{c}}=0.64$. From Eq. (6), this parameter setting corresponds to
$U_{0}=1\times 10^{-6}\ \mathrm{m/s}$ for $\varphi=0.5$.

The drying rate $U$ represents the decrease rate of liquid due to evaporation
and equals the mass transfer rate of the vapor from the free surface as

$\displaystyle U=\frac{KV_{\mathrm{m}}}{k_{\mathrm{B}}T}(p-p_{\infty}),$ (1)

where $K$ is the mass transfer coefficient of the vapor, $V_{\mathrm{m}}$ is
the molecular volume of the liquid, $T$ is the thermodynamic temperature,
$k_{\mathrm{B}}$ is the Boltzmann constant. Equation (1) describes the
diffusion of the vapor from the free surface to the ambient air stream. The
stagnant region in vicinity to the free surface is considered as the boundary
layer where the vapor diffuses. Based on Fick’s law, the mass transfer rate is
proportional to the difference in the vapor pressure of saturation $p$ and
that of the ambient air $p_{\infty}$.

We then consider the drying of colloidal suspensions where concentrated
particle layers form as illustrated in Fig. 1. Evaporation continues by liquid
transport through the particle layer to the free surface. This transport flow
is originated by capillarity due to the curvature of the menisci formed among
particles. The permeation resistance of the particle layer induces a pressure
drop $\Delta P$, which is balanced by the capillary pressure. Due to the
pressure drop, the vapor pressure at the free surface decreases from $p_{0}$
as

$\displaystyle p=p_{0}\exp\left(-\frac{V_{\mathrm{m}}\Delta
P}{k_{\mathrm{B}}T}\right),$ (2)

where $V_{\mathrm{m}}\Delta P$ gives the chemical potential drop of the
liquid. Equation (2) is known as the Kelvin equation. Since the liquid flow
rate equals the drying rate due to mass balance, the pressure drop is related
to the drying rate by Darcy’s law as

$\displaystyle\Delta P=\eta RU,$ (3)

where $R$ is the permeation resistance of the particle layer, and $\eta$ is
the liquid viscosity. Combining Eqs. (1) – (3) yields an equation that relates
the drying rate to the permeation resistance as

$\displaystyle\tilde{U}=\frac{e^{-(1-\varphi)\tilde{R}\tilde{U}}-\varphi}{1-\varphi},$
(4)

where the relative humidity is given by

$\displaystyle\varphi=\frac{p_{\infty}}{p_{0}}.$ (5)

In Eq. (4), the drying rate and the permeation resistance are given as
dimensionless quantities. The drying rate $U$ is scaled by that without
particle layers $U_{0}$ as

$\displaystyle\tilde{U}=\frac{U}{U_{0}},\ \ \ \ \
U_{0}=\frac{KV_{\mathrm{m}}}{k_{\mathrm{B}}T}(p_{0}-p_{\infty}).$ (6)

The permeation resistance $R$ is scaled by a characteristic resistance $R_{0}$
defined as

$\displaystyle\tilde{R}=\frac{R}{R_{0}},\ \ \ \ \
R_{0}=\frac{(k_{\mathrm{B}}T)^{2}}{\eta
KV_{\mathrm{m}}^{2}p_{0}}=\frac{1}{\eta
U_{0}}\frac{k_{\mathrm{B}}T}{V_{\mathrm{m}}}(1-\varphi).$ (7)

As compared to Eq. (3), the definition of $R_{0}$ indicates that the flow rate
through the particle layer with a permeation resistance of $R_{0}$ equals
$U_{0}$ when a pressure of $(k_{\mathrm{B}}T/V_{\mathrm{m}})(1-\varphi)$ is
applied. In the previous study, Eq. (4) was derived with ignoring humidity as
$\varphi=0$ and was approximately solved with assuming $V_{\mathrm{m}}\Delta
P/(k_{\mathrm{B}}T)\ll 1$, udry5  although this assumption only holds for
$\tilde{R}\ll 1$.

The solution of Eq. (4) is expressed by the Lambert $W$ function that
satisfies $x=W(xe^{x})$ as

$\displaystyle\tilde{U}=\frac{1}{1-\varphi}\left[\frac{1}{\tilde{R}}W(\tilde{R}e^{\varphi\tilde{R}})-\varphi\right].$
(8)

As described in Fig. 2, the drying rate decreases with an increase in the
permeation resistance. The asymptotic behavior is different depending on the
humidity. The asymptotic form for $\varphi>0$ decreases proportional to
$\tilde{R}^{-1}$ as

$\displaystyle\tilde{U}\sim-\frac{\ln\varphi}{1-\varphi}\frac{1}{\tilde{R}}\ \
\ \ (\tilde{R}\rightarrow\infty),$ (9)

while that for $\varphi=0$ decreases more slowly as

$\displaystyle\tilde{U}\sim\frac{\ln\tilde{R}}{\tilde{R}}\ \ \ \
(\tilde{R}\rightarrow\infty).$ (10)

Combining Eqs. (3) and (9) shows that the pressure drop $\Delta P$ approaches
to $-(k_{\mathrm{B}}T/V_{\mathrm{m}})\ln\varphi$ as
$\tilde{R}\rightarrow\infty$. If this limiting pressure drop is larger than
the maximum capillary pressure achieved by the maximum meniscus curvature, the
free surface will invade into the particle layer during drying. cap

Equation (8) becomes a simple expression in the limit of $\varphi\rightarrow
1$ as

$\displaystyle\lim_{\varphi\rightarrow 1}\tilde{U}=\frac{1}{1+\tilde{R}}.$
(11)

This expression indicates that the drying rate equals the rate of the liquid
flow through the series resistance of $R$ and $R_{0}$ when a pressure of
$(k_{\mathrm{B}}T/V_{\mathrm{m}})(1-\varphi)$ is applied. From this
expression, the characteristic resistance $R_{0}$ can be interpreted as the
resistance of vapor transport in the boundary layer. As noted in B, Eq. (11)
has the same form as the Ruth’s filtration equation that describes the
temporal variation of the flow rate in dead-end filtration. The Ruth’s
filtration equation consider the series permeation resistance of the filter
media and the cake deposited onto the filter. filt

When we assume that the particle layer is spatially uniform, the permeation
resistance is proportional to the thickness of the layer $L$ as
$R=\rho_{\mathrm{c}}L$ with a constant resistivity $\rho_{\mathrm{c}}$. We
define a characteristic thickness as $L_{0}=R_{0}/\rho_{\mathrm{c}}$. The
drying rate will decrease significantly when the particle layer grows as
$L>L_{0}$. The permeation resistivity can be evaluated by the Kozeny–Carman
equation given by

$\displaystyle\rho_{\mathrm{c}}=\frac{180}{d^{2}}\frac{\Phi^{2}_{\mathrm{c}}}{(1-\Phi_{\mathrm{c}})^{3}},$
(12)

where $d$ is the particle diameter, and $\Phi_{\mathrm{c}}$ is the volume
fraction of the particles in the particle layer. With assuming the uniformity
of the particle layer, the present model gives an analytical expression of the
drying characteristics as derived in B. As shown in Fig. 3, the characteristic
thickness $L_{0}$ increases proportional to $d^{2}$. We thus expect that the
decrease in the drying rate will not be obvious for large particles. In fact,
in the previous experiments of unidirectional drying in confined systems, the
decrease in the drying rate was observed only for the particles of $d\lesssim
40\ \mathrm{nm}$. udry1 ; udry2 ; udry3 ; udry4 ; udry5 ; udry6  The condition
that the drying rate will obviously decrease during drying can be estimated as
$(\Phi_{0}/\Phi_{\mathrm{c}})H_{0}>L_{0}$, where $\Phi_{0}$ and $H_{0}$ are
the initial particle volume fraction and thickness of the coating film of the
colloidal suspension, respectively. In this estimation,
$(\Phi_{0}/\Phi_{\mathrm{c}})H_{0}$ equals the particle layer thickness when
all the particles in the suspension form the particle layer. When the
particles of $d=10\ \mathrm{nm}$ form a particle layer with
$\Phi_{\mathrm{c}}=0.64$ (random close packing rcp ), the characteristic
thickness is $L_{0}=4.3\ \mathrm{mm}$. The decrease in the drying rate will be
obvious when the initial film thickness satisfies $H_{0}>27.5\ \mathrm{mm}$
for $\Phi_{0}=0.1$.

###  II.2 Formation of concentrated particle layers

We consider a suspension containing spherical colloidal particles of a
diameter $d$ and a mass $m$. The suspension is coated on a plane substrate
where we set the $z$-coordinate along the vertical direction. The coating film
is bounded by two parallel planes: the top free surface at $z=H$ and the
bottom substrate at $z=0$. We describe the Brownian motion of the particles by
the Langevin equations and solve them numerically. For the $i$-th particle,
the time evolution of the velocity $\bm{v}_{i}$ and the position $\bm{r}_{i}$
is expressed by snap1 ; seg2 ; snap2 ; snap3

$\displaystyle m\frac{\mathrm{d}\bm{v}_{i}}{\mathrm{d}t}$
$\displaystyle=-\zeta\bm{v}_{i}+\bm{F}_{i}^{\mathrm{R}}+\sum_{j}(\bm{F}_{ij}^{\mathrm{cnt}}+\bm{F}_{ij}^{\mathrm{DLVO}})+\bm{F}_{i}^{\mathrm{cpl}},$
$\displaystyle\frac{\mathrm{d}\bm{r}_{i}}{\mathrm{d}t}$
$\displaystyle=\bm{v}_{i}.$ (13)

As the influence of the ambient liquid, the hydrodynamic drag
$-\zeta\bm{v}_{i}$ and the random force $\bm{F}_{i}^{\mathrm{R}}$ due to
thermal fluctuations are considered, while hydrodynamic interactions among the
particles are neglected. The hydrodynamic drag is expressed by Stokes’ law as
$\zeta=3\pi\eta d$. The three components of the random force are given as the
stochastic variables obeying independent Gaussian distributions satisfying

$\displaystyle\langle\bm{F}_{i}^{\mathrm{R}}(t)\rangle=\bm{0},\ \ \ \ \
\langle\bm{F}_{i}^{\mathrm{R}}(t)\bm{F}_{i}^{\mathrm{R}}(0)\rangle=2k_{\mathrm{B}}T\zeta\delta(t)\bm{I}.$
(14)

The forces $-\zeta\bm{v}_{i}$ and $\bm{F}_{i}^{\mathrm{R}}$ cause the Brownian
motion of the particle, whose diffusion coefficient is given by the
Stokes–Einstein relation $D=k_{\mathrm{B}}T/\zeta$.

The contact force $\bm{F}_{ij}^{\mathrm{cnt}}$ describes collisions between
the particles and is given by the Voigt model with the Herzian contact theory:
Herz

$\displaystyle\bm{F}_{ij}^{\mathrm{cnt}}=\left[\frac{1}{3}E^{\ast}d^{1/2}\delta_{ij}^{3/2}-\Gamma(\bm{v}_{i}-\bm{v}_{j})\cdot\hat{\bm{n}}_{ij}\right]\hat{\bm{n}}_{ij}.$
(15)

The summation on the index $j$ includes the force from the substrate as well
as the other particles. The contact force acts in the direction normal to the
tangent plane at the contact point. This direction is expressed by the unit
vector $\hat{\bm{n}}_{ij}=\bm{r}_{ij}/|\bm{r}_{ij}|$ with
$\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$. The normal relative displacement of the
contact point is given by $\delta_{ij}=\max(0,d-|\bm{r}_{ij}|)$. The
longitudinal elastic modulus $E^{\ast}$ is given by the Young’s modulus $E$
and the Poisson ratio $\nu$ as $E^{\ast}=E/(1-\nu^{2})$. We assume that the
damping coefficient is related to the elastic modulus as
$\Gamma=\lambda(mE^{\ast}d^{1/2}\delta_{ij}^{1/2})^{1/2}$ to describe
partially inelastic collisions. cnt  The contact force from the substrate is
expressed by Eq. (15) with substitution of $E^{\ast}\rightarrow
2^{1/2}E^{\ast}$ and $\delta_{ij}=\max(0,d/2-z_{i})$, where
$z_{i}=\hat{\bm{z}}\cdot\bm{r}_{i}$ is the $z$-coordinate of the $i$-th
particle. When an adhesive force $F_{\mathrm{a}}$ acts between the particles,
there is an equilibrium displacement given by
$\delta_{0}=[3F_{\mathrm{a}}/(E^{\ast}d^{1/2})]^{2/3}$. Considering small
oscillation around $\delta_{0}$, the elastic force is approximated by a linear
restoring force with a spring constant of $E^{\ast}(d\delta_{0})^{1/2}/2$.
Critical damping is realized when the factor of the damping coefficient
$\Gamma$ satisfies $\lambda=1$, which we set in the present study.

The DLVO force $\bm{F}_{ij}^{\mathrm{DLVO}}$ describes the force acting
between charged particles in liquid. This force is given by the DLVO potential
$V$ as dlvo

$\displaystyle\bm{F}_{ij}^{\mathrm{DLVO}}=\left.-\frac{\mathrm{d}V(h)}{\mathrm{d}h}\right|_{h=h_{ij}}\hat{\bm{n}}_{ij},$
(16) $\displaystyle
V(h)=-\frac{Ad}{24h}+\pi\varepsilon_{\mathrm{r}}\varepsilon_{0}d\psi_{\mathrm{eff}}^{2}e^{-\kappa
h},$ (17)

where $h$ is the surface separation between the particles. The DLVO potential
is the combination of the van der Waals attraction (first term) and the
electric double layer repulsion (second term). The magnitude of the van der
Waals attraction is provided by the Hamaker constant $A$. The electric double
layer repulsion is characterized by the effective surface potential
$\psi_{\mathrm{eff}}$ and the Debye parameter $\kappa$, and they are given by

$\displaystyle\psi_{\mathrm{eff}}=\frac{4k_{\mathrm{B}}T}{Z\mathrm{e}}\tanh\left(\frac{Z\mathrm{e}\psi_{0}}{4k_{\mathrm{B}}T}\right),\
\ \
\kappa=\left(\frac{2N_{\mathrm{A}}\mathrm{e}^{2}I}{\varepsilon_{\mathrm{r}}\varepsilon_{0}k_{\mathrm{B}}T}\right)^{1/2},$
(18)

where $\mathrm{e}$ is the elementary charge, $N_{\mathrm{A}}$ is the Avogadro
constant, $\varepsilon_{0}$ is the electric constant,
$\varepsilon_{\mathrm{r}}$ is the relative permittivity, $Z$ is the ion
valence, and $\psi_{0}$ is the surface electric potential. The expression of
$\psi_{\mathrm{eff}}$ in Eq. (18) is that for valence symmetric electrolytes.
In this case, the ionic strength $I$ equals the electrolyte concentration
multiplied by a factor of $Z^{2}$. The inverse of the Debye parameter
$\kappa^{-1}$ serves as a measure of the thickness of the electric double
layer. Note that the expression of Eq. (17) is the approximate form for
$\kappa^{-1}\ll h\ll d$, and we use it for simplicity. We introduce a cut off
distance $h_{\mathrm{c}}$ within which the DLVO force becomes constant, and
hence the adhesive force between the particles is given by
$F_{\mathrm{a}}=(\mathrm{d}V/\mathrm{d}h)_{h=h_{\mathrm{c}}}$. The cut off
distance can be interpreted as the atomic length scale where the continuum
picture, which is assumed in Eq. (17), breaks down. dlvo  We calculate the
DLVO force with setting $h_{ij}=\max(h_{\mathrm{c}},|\bm{r}_{ij}|-d)$.

The capillary force $\bm{F}_{i}^{\mathrm{cpl}}$ acts to restore the particles
to the position satisfying an equilibrium contact angle with the free surface.
When the particles prefer to be perfectly wet, the capillary force is
expressed as

$\displaystyle\bm{F}_{i}^{\mathrm{cpl}}=-4\pi\gamma\delta_{i}^{\mathrm{fs}}\left(1-\frac{\delta_{i}^{\mathrm{fs}}}{d}\right)\hat{\bm{z}},$
(19)

where $\gamma$ is the surface tension of the liquid. The protruding length of
the particle from the free surface is given by
$\delta_{i}^{\mathrm{fs}}=\max[0,d/2-(H-z_{i})]$. To describe the evaporation
of the liquid, the free surface moves in the $-z$ direction with a rate of
$U$. The particles are swept by the free surface through the capillary force,
thereby forming concentrated particle layers. Although the ion concentration
and the temperature can change during drying in real systems, we assume that
they are uniformly constant and that the DLVO potential does not change.

According to Eq. (1), the drying rate $U$ changes with reflecting the
formation of concentrated particle layers. We define the particle layer as the
cluster of the particles moving with the free surface. The clustering
particles are detected as the particles contacting each other, i.e.
$\delta_{ij}<0$. When one of the particles that compose the cluster satisfies
$\delta_{i}^{\mathrm{fs}}<0$, we regard the cluster as the particle layer. The
permeation resistance of the particle layer is calculated by integrating the
resistivity of the cross-section of the computational domain vertical to the
$z$-axis:

$\displaystyle R=\int_{0}^{H}\rho(z)\mathrm{d}z,\ \ \ \
\rho(z)=\frac{80}{[D_{\mathrm{H}}(z)]^{2}}\frac{S}{S_{\mathrm{f}}(z)},$ (20)

where $S$ is the cross-sectional area of the computational domain,
$S_{\mathrm{f}}$ is that of the liquid flow path, and $D_{\mathrm{H}}$ is the
hydraulic diameter. Since the cross-section of the flow path through the
particle layer is quite complex shape, the resistivity is evaluated by
regarding the flow path as a cylindrical channel with a diameter of
$D_{\mathrm{H}}$. resist  The evaluation of the resistivity in Eq. (20) is an
analogy to the resistivity of a cylindrical channel whose diameter is $D$,
i.e. $32/D^{2}$, which is derived from the Hagen–Poiseuille equation. The
factor $S/S_{\mathrm{f}}$ appears to consider the variation of the flow rate
inversely proportional to the cross-sectional area of the flow path. The
factor 80 reflects the channel tortuosity and is selected such that the
Kozeny–Carman equation Eq. (12) can be derived from Eq. (20) as indicated
below. The hydraulic diameter $D_{\mathrm{H}}$ is given by the cross-sectional
area of the flow path $S_{\mathrm{f}}$ and the wetted perimeter of the cross-
section $L_{\mathrm{f}}$ as resist

$\displaystyle D_{\mathrm{H}}=\frac{4S_{\mathrm{f}}}{L_{\mathrm{f}}},$ (21)
$\displaystyle S_{\mathrm{f}}(z)=S-\frac{\pi}{4}d^{2}\sum_{i\in
P}\left[1-\frac{4(z-z_{i})^{2}}{d^{2}}\right]\Theta\left(1-\frac{2|z-z_{i}|}{d}\right),$
(22) $\displaystyle L_{\mathrm{f}}(z)=\pi d\sum_{i\in
P}\Theta\left(1-\frac{2|z-z_{i}|}{d}\right),$ (23)

where the Heaviside step function $\Theta$ appears to indicate the region
where the particles occupy: $|z_{i}-z|\leq d/2$. The index in the summation
$i\in P$ indicates the particles that compose the particle layer. Assuming the
uniform distribution of the particles with a volume fraction of
$\Phi_{\mathrm{c}}$, we obtain $S_{\mathrm{f}}=(1-\Phi_{\mathrm{c}})S$ and
$L_{\mathrm{f}}=6\Phi_{\mathrm{c}}S/d$, which lead to the Kozeny–Carman
equation given by Eq. (12).

##  III Results and discussion

Table 1:  Time scales in the motion of particles.

---
Time scale | Ratio to $\tau_{\mathrm{hyd}}$
Migration (Characteristic length: $d$)
Evaporation |  $\tau_{\mathrm{E}}=d/U_{0}$ | $\ 25$
Diffusion |  $\tau_{\mathrm{D}}=d^{2}/D$ | $\ 1.0\times 10^{4}$
Velocity change
Hydrodynamic drag |  $\tau_{\mathrm{hyd}}=m/\zeta$ | $\ 1$
Contact force |  $\tau_{\mathrm{cnt}}=(m/E^{\ast})^{1/2}(d\delta_{0})^{-1/4}$ | $\ 1.0\times 10^{-1}$
Capillary force |  $\tau_{\mathrm{cpl}}=(m/\gamma)^{1/2}$ | $\ 3.1\times 10^{-1}$

Table 2:  Dimensionless Debye parameter $\kappa d$ and potential barrier
$\Delta E$ of DLVO potentials for different ionic strengths.

---
$I/(\mathrm{mol/L})$ | $\ \ \kappa d\ \ \ $ | $\Delta E/(k_{\mathrm{B}}T)$
0.01 | 6.5 | $15.8\ \ \ $
0.03 | 11.3 | $13.7\ \ \ $
0.1 | 20.5 | $10.7\ \ \ $
0.3 | 35.6 | $7.5\ \ \ $
1 | 64.9 | $3.5\ \ \ $

Figure 4:  DLVO potentials for different ionic strengths, as functions of
surface separation between particles.  Figure 5:  Potential barrier (left
axis) and the average contact number of particles for the initial
configurations (right axis), as functions of ionic strength. Corresponding
snapshots of the initial configurations are shown in upper row, where the
contact number of the particles is indicated by their color.

Numerical simulations are performed with varying the DLVO potential. We
investigate the effects of the interaction between the particles on the drying
characteristics. To describe a part of the coating film far from edges,
periodic boundary conditions are applied in the directions parallel to the
substrate with a side length of $18d$. In the initial coating film of the
thickness of $50d$, the volume fraction of the particles is set as
$\Phi_{0}=0.1$.

The particle motion described by Eq. (13) contains multiple time scales as
listed on Table 1. The time integration of Eq. (13) is performed by use of the
velocity Verlet algorithm. The time increment is set as
$0.2\tau_{\mathrm{cnt}}$ to be smaller than the shortest time scale in Table
1, i.e. $\tau_{\mathrm{cnt}}$. The drying Péclet number is
$\mathrm{Pe}=\tau_{\mathrm{D}}/\tau_{\mathrm{E}}=400$, which corresponds to
$U_{0}=0.44\ \mathrm{m/s}$ for an aqueous suspension with $d=20\ \mathrm{nm}$,
$\eta=1\times 10^{-3}\ \mathrm{Pa\ s}$, and $T=300\ \mathrm{K}$. This drying
rate is larger by a factor of about $10^{6}$ than actual situations. The
drying Péclet number represents also the ratio of the particle diameter to the
diffusion length of the particles concentrated by the receding free surface,
clyr1 ; seg3  i.e. $\mathrm{Pe}=d/(D/U_{0})$. The current setting of
$\mathrm{Pe}$ corresponds to reducing the length scale in the drying direction
to perform simulations in smaller computational domains. We correspondingly
set the characteristic thickness given in Sec. II.1 as $L_{0}/d=2$, which is
about $10^{-6}$ times smaller than actual situations as demonstrated in Fig.
3. The humidity is set as $\varphi=0.5$.

The dimensionless parameters of the DLVO potential are set as
$A/(24k_{\mathrm{B}}T)=0.1$ and
$\pi\varepsilon_{\mathrm{r}}\varepsilon_{0}d\Psi_{\mathrm{eff}}^{2}/(k_{\mathrm{B}}T)=23$.
This setting corresponds to considering the silica particles in water with
$d=20\ \mathrm{nm}$, $T=300\ \mathrm{K}$, $\psi_{0}=-50\ \mathrm{mV}$, and
$A=1\times 10^{-20}\ \mathrm{J}$. hamak  As shown in Table 2 and Fig. 4, we
change the other dimensionless parameter $\kappa d$ by changing the ionic
strength according to Eq. (18). We define the potential barrier $\Delta E$ as
the difference between the maximum and the local minimum (potential well) of
the DLVO potential. The potential barrier decreases with an increase in the
ionic strength. The cut off length $h_{\mathrm{c}}$ is set to satisfy
$\delta_{0}/d=0.02$. In the present simulations, we consider the DLVO force
only between the particles with the separation of $h/d\leq 1$, where the DLVO
potential sufficiently decays since all the settings in Table 2 satisfies
$(\kappa d)^{-1}<1$. The initial configurations exhibited in Fig. 5 are
obtained as follows: the particles are firstly distributed randomly so that
the surface separations between the particles become more than that of the
potential maximum, and then the Brownian motion of the particles are solved
during a time of $5\tau_{\mathrm{D}}$. The particles aggregate at $I\geq 0.3\
\mathrm{mol/L}$ where $\Delta E/(k_{\mathrm{B}}T)\lesssim 10$ is satisfied.
The average contact number of the particles rises sharply accompanied by the
aggregation.

Figure 6:  Snapshots of the particle configurations during drying. The ionic
strengths are (a) $0.1\ \mathrm{mol/L}$ and (b) $1\ \mathrm{mol/L}$. The
particles composing the concentrated layers are indicated by the red color.
Figure 7:  Time variation of (a) coating film thickness and (b) drying rate
for different ionic strengths. The dotted line represents the time change of
the film thickness containing no particles.

As shown in Fig. 6, the receding free surface sweeps the particles to form the
concentrated particle layers. Figure 7 thus indicates that the drying rate
decreases with the growth of the particle layers. If there are no particles,
the thickness of the liquid film decreases linearly with a constant drying
rate. The drying rate decreases most slowly at the highest ionic strength of
$I=1\ \mathrm{mol/L}$ where the particles aggregate. Comparing the cases of
$I\leq 0.1\ \mathrm{mol/L}$ where the particles are stably dispersed, the
drying rate is initially larger at the higher ionic strength of $I=0.1\
\mathrm{mol/L}$, while it finally becomes larger at the lower ionic strength
of $I=0.01\ \mathrm{mol/L}$. This difference reflects the permeation
resistance as discussed below.

Figure 8:  Distributions of particle volume fraction at the times of
$t/\tau_{\mathrm{E}}=0,5,10,30,50$ for different ionic strengths. The broken
lines represent the position of the free surface.

Figure 8 demonstrates the growth of the particle layers as the increase in the
particle volume fraction beneath the free surface, where the clear oscillatory
profiles can be found except at $I=1\ \mathrm{mol/L}$. This oscillatory
profile indicates the stratification of the particles. The DLVO potential
barrier hinders the particles to contact each other, thereby allowing the
concentrated particles to rearrange to stratify. When the potential barrier is
too low so that the particles aggregate at $I=1\ \mathrm{mol/L}$, the
particles cannot easily rearrange and form loose particle layers, thus the
oscillatory profiles do not appear.

Figure 9:  Time variation of (a) the permeation resistance and (b) the
resistivity of concentrated particle layers for different ionic strengths. The
resistivity is scaled by that evaluated by the Kozeny–Carman equation Eq. (12)
with $\Phi_{\mathrm{c}}=0.64$.

The permeation resistance increases with the growth of the concentrated
particle layers as shown in Fig. 9(a). The increase of the permeation
resistance results in the decrease of the drying rate through Eq. (4). The
larger permeation resistance is resulted from the larger particle volume
fraction of the particle layers. In the cases of $I\leq 0.1\ \mathrm{mol/L}$,
the permeation resistance increases sharply at $t/\tau_{\mathrm{E}}\approx 5$,
which seems to correspond to the time when the first layer of the particle
layer is formed. In fact, this time is roughly evaluated as
$t/\tau_{\mathrm{E}}=(\Phi_{\mathrm{c}}-\Phi_{0})/\Phi_{0}$, which equals
$5.4$ with assuming $\Phi_{\mathrm{c}}=0.64$. This evaluation is derived from
the particle mass balance of $(\Phi_{\mathrm{c}}-\Phi_{0})d=\Phi_{0}U_{0}t$.
Comparing the cases of $I\leq 0.1\ \mathrm{mol/L}$, the permeation resistance
is initially larger for the smaller ionic strength of $I=0.01\
\mathrm{mol/L}$. In Fig. 8, the initial distributions at
$t/\tau_{\mathrm{E}}=0$ show that the stronger repulsive interaction between
the particles results in the larger particle volume fraction beneath the free
surface, and thus the permeation resistance becomes larger. Such increased
accumulation of particles near a wall by repulsive interactions between
particles was also shown in the previous study. wall  However, the permeation
resistance finally becomes larger at the larger ionic strength of $I=0.1\
\mathrm{mol/L}$, corresponding to the larger particle volume fraction of the
particle layer as demonstrated by the profiles at $t/\tau_{\mathrm{E}}\geq 10$
in Fig. 8. This larger particle volume fraction can be attributed to the
smaller potential barrier and potential range, which would allow the particles
to be densely packed.

In Fig. 9(b), the average resistivity is evaluated as
$\bar{\rho}=R/[2(H-\langle z\rangle_{P})]$, where $\langle z\rangle_{P}$ is
the mean $z$-coordinate of the particles that compose the particle layers, and
$2(H-\langle z\rangle_{P})$ gives an estimation of the thickness of the
particle layers. The resistivity also increases with time but converges in the
cases of the stably dispersed particles at $I\leq 0.1\ \mathrm{mol/L}$.
Although the particle volume fraction of the particle layer is not uniform
especially near the free surface as shown in Fig. 8, the convergence of the
resistivity suggests that the particle layer gradually become uniform with
distance from the free surface. However, such convergence is not observed at
$I=1\ \mathrm{mol/L}$, where the particle volume fraction of the particle
layer continuously increases with time. When we can observe a convergence of
the resistivity in the Langevin dynamics simulations, combining this evaluated
resistivity and the continuum model given in B would enable us to evaluate the
drying characteristics of coating films whose thickness is larger than the
simulation scale.

Figure 10:  Drying rate and the particle volume fraction of concentrated
particle layers at $t/\tau_{\mathrm{E}}=50$, as functions of ionic strength.

Figure 10 summarizes the effect of the interaction between the particles on
the drying characteristics. The particle volume fraction of the particle layer
$\bar{\Phi}$ is evaluated as the mean particle volume fraction in the region
within the depth of $3d$ from the free surface, i.e. $[H-3d,H]$. As long as
the particles stably dispersed at $I\leq 0.1\ \mathrm{mol/L}$, decreasing the
potential barrier by increasing the ionic strength results in a little
decrease in the drying rate with a little increase in the particle volume
fraction of the particle layer. On the other hand, when the particles
aggregate at the higher ionic strength of $I>0.1\ \mathrm{mol/L}$, we find
obvious change that the particle layer becomes loose and the drying rate
increases. In the previous experimental study of drying of colloidal
suspension droplets, it was also indicated that the morphology of the dried
grains can be changed by forming aggregates, due to the decrease in the
permeation resistance of the particle layers. sdry3

##  IV Conclusions

We constructed a model that can calculate the drying characteristics by use of
Langevin dynamics simulations. This model enables us to evaluate the drying
characteristics of colloidal suspensions from the character of colloidal
particles such as interactions. The present calculation of this model shows
that the formation of aggregates by attractive interactions suppresses the
decrease in the drying rate. This suppression is due to the loose structure of
the concentrated particle layers. The aggregation is thus preferable to reduce
the drying time, while the final dried structure would be loose and
disordered. snap3  To improve the density and order of structures with
reducing the drying time by aggregation, further investigations considering
the effects of morphology and strength of the aggregates would be required.
Note that the present finding will be true for filtration. The aggregated
particles are expected to form loose cake, thereby suppressing the decrease in
the filtering flux. Both filtration and drying are operations of solid–liquid
separation where liquid is removed through an interface.

## Appendix

###  A Formulation of drying rates in terms of chemical potential

Figure A1:  Sketch of the chemical potential profile of the liquid. The
liquid-to-vapor phase transition occurs at the free surface. The chemical
potential drops in two steps: permeation through the concentrated particle
layer ($\Delta\mu_{\mathrm{p}}$) and diffusion through the boundary layer
($\Delta\mu_{\mathrm{b}}$).

The evaporation can also be considered to be driven by the difference in
chemical potential, instead of the difference in vapor pressure as indicated
in Eq. (1). As shown in Fig. A1, when a concentrated particle layer is formed
beneath the free surface, the chemical potential of the liquid drops in two
steps. The total chemical potential drop is $-k_{\mathrm{B}}T\ln\varphi$ that
equals the chemical potential difference between the vapor of the pressure
$p_{0}$ and $p_{\infty}$:

$\displaystyle\Delta\tilde{\mu}_{\mathrm{b}}+\Delta\tilde{\mu}_{\mathrm{p}}=\Delta\tilde{\mu}=-\ln\varphi,$
(A1)

where the chemical potential is scaled by $k_{\mathrm{B}}T$ and dimensionless
quantities are denoted by the tilde. The chemical potential drops through the
boundary layer and the particle layer are given as

$\displaystyle\Delta\tilde{\mu}_{\mathrm{b}}=\ln\frac{\tilde{p}}{\varphi},$
(A2) $\displaystyle\Delta\tilde{\mu}_{\mathrm{p}}=\frac{V_{\mathrm{m}}\Delta
P}{k_{\mathrm{B}}T},$ (A3)

respectively, where we denote $\tilde{p}=p/p_{0}$. The drying rate is
expressed to be proportional to the chemical potential drops in each step as

$\displaystyle\tilde{U}=\Lambda_{\mathrm{b}}(\tilde{p})\Delta\tilde{\mu}_{\mathrm{b}}=\Lambda_{\mathrm{p}}\Delta\tilde{\mu}_{\mathrm{p}}.$
(A4)

Comparing Eq. (A4) with Eqs. (1) and (3), the transport coefficients are
expressed as

$\displaystyle\Lambda_{\mathrm{b}}(\tilde{p})=\frac{\tilde{p}-\varphi}{(1-\varphi)\ln(\tilde{p}/\varphi)},$
(A5) $\displaystyle\Lambda_{\mathrm{p}}=\frac{1}{(1-\varphi)\tilde{R}}.$ (A6)

Note that the former depends on the vapor pressure at the free surface $p$.
Combining Eqs. (A1) and (A4), we finally obtain an expression of the
relationship between the drying rate and the total chemical potential drop as

$\displaystyle\tilde{U}=\left[\frac{1}{\Lambda_{\mathrm{b}}(\tilde{p})}+\frac{1}{\Lambda_{\mathrm{p}}}\right]^{-1}\Delta\tilde{\mu}.$
(A7)

In this expression, the coefficient of $\Delta\tilde{\mu}$ represents the
overall transport coefficient. In the limit of $\varphi\rightarrow 1$,
$\Lambda_{\mathrm{b}}\sim(1-\varphi)^{-1}$ and
$\Delta\tilde{\mu}\sim(1-\varphi)$ result in Eq. (11) that indicates a flow
rate through serial resistance.

###  B Drying characteristics associated with the growth of uniform particle
layers

Figure A2:  Time variation of drying rate at different humidities. The
profiles are given by Eq. (A12) for $0<\varphi<1$; Eq. (A16) for $\varphi=0$;
and Eq. (A18) for $\varphi\rightarrow 1$.  Figure A3:  Time variation of
permeation resistance at different humidities. This describes the growth of
the concentrated particle layer. The profiles are given by combining Eqs. (4)
and (A12) for $0<\varphi<1$; combining Eqs. (4) and (A16) for $\varphi=0$; and
Eq. (A19) for $\varphi\rightarrow 1$

Assuming that the density of the particle layer is spatially uniform and does
not change in time, the growth of the particle layer can be described by a
one-dimensional differential equation as

$\displaystyle\frac{\mathrm{d}L}{\mathrm{d}t}=\alpha U,\ \ \ \ \
\alpha=\frac{\Phi_{0}}{\Phi_{\mathrm{c}}-\Phi_{0}}.$ (A8)

This equation is derived from the particle mass balance:
$(\Phi_{\mathrm{c}}-\Phi_{0})\mathrm{d}L=\Phi_{0}U\mathrm{d}t$. With the
relation $R=\rho_{\mathrm{c}}L$, a dimensionless form of Eq. (A8) is obtained
as

$\displaystyle\frac{\mathrm{d}\tilde{R}}{\mathrm{d}\tilde{t}}=\tilde{U},$ (A9)

where the time is scaled as

$\displaystyle\tilde{t}=\frac{t}{t_{0}},\ \ \ \ \ t_{0}=\frac{L_{0}}{\alpha
U_{0}}.$ (A10)

An explicit expression of $\tilde{R}$ as a function of $\tilde{U}$ is obtained
from Eq. (4). Substituting this explicit expression to Eq. (A9) yields the
following differential equation:

$\displaystyle\frac{\mathrm{d}\tilde{U}}{\mathrm{d}\tilde{t}}=\tilde{U}^{3}\left\\{\frac{\ln[(1-\varphi)\tilde{U}+\varphi]}{1-\varphi}-\frac{\tilde{U}}{(1-\varphi)\tilde{U}+\varphi}\right\\}^{-1}.$
(A11)

The solution of Eq. (A11) is obtained as

$\displaystyle\tilde{t}=$
$\displaystyle-\frac{\ln[(1-\varphi)\tilde{U}+\varphi]}{2(1-\varphi)\tilde{U}^{2}}+\frac{1}{2\varphi\tilde{U}}$
$\displaystyle-\frac{1-\varphi}{2\varphi^{2}}\ln\left[1+\frac{\varphi}{(1-\varphi)\tilde{U}}\right]-\frac{1}{2\varphi}\left[1+\frac{1-\varphi}{\varphi}\ln(1-\varphi)\right].$
(A12)

As shown in Figs. A2 and A3, the drying rate decreases with time, and thus the
increasing rate of the permeation resistance is reduced. Unlike the simulation
results of $I\leq 0.1\ \mathrm{mol/L}$ in Fig. 9(a), sharp increases at an
early stage do not appear. This is because the model in this section assumes
the growth of spatially uniform particle layers even when their thickness is
less than the particle diameter. The influence of the humidity does not appear
in the initial behavior of $\tilde{R}\sim\tilde{t}$ as $\tilde{t}\rightarrow
0$, but it appears with time. The long-time asymptotic behavior for
$\varphi>0$ is

$\displaystyle\tilde{U}\sim\left[\frac{-\ln\varphi}{2(1-\varphi)}\right]^{1/2}\tilde{t}^{-1/2}\
\ \ \ (\tilde{t}\rightarrow\infty),$ (A13)
$\displaystyle\tilde{R}\sim\left(\frac{-2\ln\varphi}{1-\varphi}\right)^{1/2}\tilde{t}^{1/2}\
\ \ \ (\tilde{t}\rightarrow\infty),$ (A14)

while a different asymptotic behavior is derived for $\varphi=0$ as follows:

$\displaystyle\frac{\mathrm{d}\tilde{U}}{\mathrm{d}\tilde{t}}=\frac{\tilde{U}^{3}}{\ln\tilde{U}-1},$
(A15)
$\displaystyle\tilde{t}=\frac{1}{4}\left(\frac{1-2\ln\tilde{U}}{\tilde{U}^{2}}-1\right)\sim-\frac{1}{2}\frac{\ln\tilde{U}}{\tilde{U}^{2}}\
\ \ \ (\tilde{t}\rightarrow\infty).$ (A16)

This expression indicates that the drying rate decreases more slowly than in
the cases of $\varphi>0$.

In the limit of $\varphi\rightarrow 1$, Eq. (A11) becomes a simple form as

$\displaystyle\frac{\mathrm{d}\tilde{U}}{\mathrm{d}\tilde{t}}=-\tilde{U}^{3}.$
(A17)

The solution of Eq. (A17) is obtained as

$\displaystyle\tilde{U}=(1+2\tilde{t})^{-1/2},$ (A18)
$\displaystyle\tilde{R}=(1+2\tilde{t})^{1/2}-1.$ (A19)

This result corresponds to the Ruth’s filtration equation that describes the
temporal variation of the flow rate in dead-end filtration. filt  The Ruth’s
filtration equation is derived from Darcy’s law with considering the series
permeation resistance of the filter media and the cake deposited onto the
filter. In the limit of $\varphi\rightarrow 1$, the present model also becomes
an expression indicating series resistance as Eq. (11).

In the previous experiments of unidirectional drying in confined systems,
udry1 ; udry2 ; udry3 ; udry4 ; udry5  the change in the time dependence of
the particle layer thickness from $t$ to $t^{1/2}$ was observed. To explain
this behavior, equations that have the same form as Eq. (A19) were suggested,
but the derivation of them was empirical or approximative. udry1 ; udry5  With
considering the influence of humidity, which was ignored in the previous
studies, udry1 ; udry5  the present model yields the analytical expressions of
the drying characteristics and shows that the previously suggested expressions
hold only for the limit of $\varphi\rightarrow 1$. Figures A2 and A3 shows
that the humidity affects the drying characteristics especially for long-time
behavior. The present model would be helpful to explain the drying
characteristics observed in the experiments.

## Acknowledgements

This work was carried out under the project of Products Innovation
Association, “Structure of NAno Particles (SNAP) study group,” and was
supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI
Grant No. 19K15335 and 21K04747.

## References

## References

  * (1) S. Erkselius, L. Wadsö, O. J. Karlsson, J. Colloid Interface Sci., 2008, 317, 83.
  * (2) Y. Komoda, R. Kimura, K. Niga and H. Suzuki, Drying Technol., 2011, 29, 1037.
  * (3) S. Lim, K. H. Ahn and M. Yamamura, Langmuir, 2013, 29, 8233.
  * (4) T. Okuzono, K. Ozawa and M. Doi, Phys. Rev. Lett., 2006, 97, 136103.
  * (5) A. Osman, L. Goehring, A. Patti, H. Stitt and N. Shokri, Ind. Eng. Chem. Res., 2017, 56, 10506.
  * (6) M. Léang, D. Lairez, F. Cousin, F. Giorgiutti-Dauphiné, L. Pauchard and L.-T. Lee, Langmuir, 2019, 35, 2692.
  * (7) E. R. Dufresne, E. I. Corwin, N. A. Greenblatt, J. Ashmore, D. Y. Wang, A. D. Dinsmore, J. X. Cheng, X. Xie, J. W. Hutchinson and D. A. Weitz, Phys. Rev. Lett., 2003, 91, 224501.
  * (8) E. R. Dufresne, D. J. Stark, N. A. Greenblatt, J. X. Cheng, J. W. Hutchinson, L. Mahadevan and D. A. Weitz, Langmuir, 2006, 22, 7144.
  * (9) A. Sarkar and M. S. Tirumkudulu, Soft Matter, 2011, 7, 8816.
  * (10) S. Inasawa and Y. Yamaguchi, Soft Matter, 2012, 8, 2416.
  * (11) P. Lidon and J.-B. Salmon, Soft Matter, 2014, 10, 4151.
  * (12) K. Abe and S. Inasawa, Phys. Chem. Chem. Phys., 2018, 20, 8935.
  * (13) P. Lehmann, S. Assouline and D. Or, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 056309.
  * (14) J. Thiery, S. Rodts, D. A. Weitz and P. Coussot, Phys. Rev. Fluids, 2017, 2, 074201.
  * (15) N. Tsapis, E. R. Dufresne, S. S. Sinha, C. S. Riera, J. W. Hutchinson, L. Mahadevan and D. A. Weitz, Phys. Rev. Lett., 2005, 94, 018302.
  * (16) É. Lintingre, G. Ducouret, F. Lequeux, L. Olanier, T. Périé and L. Talini, Soft Matter, 2015, 11, 3660.
  * (17) M. Schulz and J. L. Keddie, Soft Matter, 2018, 14, 6181.
  * (18) R. Tatsumi, T. Iwao, O. Koike, Y. Yamaguchi and Y. Tsuji, Appl. Phys. Lett., 2018, 112, 053702.
  * (19) R. Tatsumi, O. Koike, Y. Yamaguchi and Y. Tsuji, J. Chem. Phys., 2020, 153, 164902.
  * (20) O. Vincent, A. Szenicer and A. D. Stroock, Soft Matter, 2016, 12, 6656.
  * (21) B. F. Ruth, Ind. Eng. Chem., 1935, 27, 708, 806.
  * (22) R. D. Kamien and A. J. Liu, Phys. Rev. Lett., 2007, 99, 155501.
  * (23) M. Fujita and Y. Yamaguchi, J. Chem. Eng. Jpn., 2006, 39, 83.
  * (24) O. Koike, S. Ohta, M. Fujita and Y. Yamaguchi, Jpn. J. Appl. Phys., 2008, 47, 8124.
  * (25) S. Usune, T. Takahashi, M. Kubo, E. Shoji, T. Tsukada, O. Koike, R. Tatsumi, M. Fujita and T. Adschiri, J. Chem. Eng. Jpn., 2019, 52, 680.
  * (26) S. P. Timoshenko and J. N. Goodier, Theory of Elasticity, McGraw-Hill, New York, 1970.
  * (27) Y. Tsuji, T. Tanaka and T. Ishida, Powder Technol., 1992, 71, 239.
  * (28) J. N. Israelachvili, Intermolecular and Surface Forces, 3rd ed., Academic Press, Burlington, MA, 2011.
  * (29) R. B. Bird, W. E. Stewart, and E. N. Lightfoot, Transport Phenomena, revised 2nd ed., Wiley, New York, 2006.
  * (30) H.-J. Butt, B. Cappella and M. Kappl, Surf. Sci. Rep., 2005, 59,1.
  * (31) S. Karanikas, J. Dzubiella, A. Moncho-Jordá and A. A. Louis, J. Chem. Phys., 2008, 128, 204704.