e1e-mail<EMAIL_ADDRESS>(Current affiliation is Neturen Co., Ltd.)
e2e-mail<EMAIL_ADDRESS>
11institutetext: Department of Physics and Material Science, Shimane
University, Matsue, Shimane 690-8504, Japan. 22institutetext: Graduate School
of Engineering Science, Osaka University, Toyonaka, Osaka 560-8531, Japan.

We numerically investigate the stress distribution of a paste when an external
oscillation is applied. The paste memorizes the oscillation through plastic
deformation. Due to the plastic deformation, the residual stress remains after
the oscillation, where the residual stress distribution depends on the number
of cycles in the oscillation. As this number increases, the symmetry of the
stress distribution is enhanced, which is consistent with the crack patterns
observed in the experiments using a drying paste.

# Memory effect of external oscillation on residual stress in a paste

Jumpei Moritae1,addr1  Michio Otsukie2,addr2

(Received: date / Accepted: date)

††journal: Eur. Phys. J. E

##  1 Introduction

When a paste containing powder and water is dried, cracks are formed Goehring
. Crack patterns on the surface of the paste are usually random and isotropic
in shallow containers Groisman . However, recent experiments have revealed
that the crack patterns become anisotropic on application of external fields
Nakahara05 ; Nakahara06a ; Nakahara06b ; Matsuo ; Nakayama ; Akiba ; Mal ;
Khatun ; Pauchard ; Ngo ; Lama . In particular, Nakahara and Matsuo reported
that subjecting the container to horizontal oscillation before desiccation
imprints a “memory” in the paste. This memory induces lamellar crack patterns
perpendicular to the direction of the oscillation after several days
Nakahara05 .

They also discovered that such crack patterns appear when the stress induced
by the oscillation exceeds the yield stress of the paste Nakahara05 . This
experimental result indicates that the plastic deformation caused by the
oscillation is related to the memory effect. The visualization of the plastic
deformation in the drying paste showing the memory effect supports this
conjecture Goehring .

The cracks are formed to release the tension caused by the desiccation when a
fracture criterion is satisfied Beuth ; Kitsunezaki ; Singh ; Man . Therefore,
the memory effect can be explained if plastic deformation affects one of the
factors in crack formation, namely the tension or the fracture criterion. In
this regard, recent experiments to measure stresses in drying paste with
memory effect Kitsunezaki16 ; Kitsunezaki17  deserve special attention.
According to these experiments, the fracture criterion is unchanged, but the
external oscillation induces the residual stress in the direction of the
oscillation, which increases due to the desiccation and leads to the
perpendicular crack pattern. However, the mechanism by which plastic
deformation causes residual stress is not clear in the experiments.

The stress-measuring experiments Kitsunezaki16 ; Kitsunezaki17  had been
motivated by the theoretical conjecture that anisotropy in residual stresses
caused by plastic deformation is responsible for the memory effect Otsuki ;
Ooshida08 ; Ooshida09 . However, theoretical predictions in the previous
studies seem to be incoherent. On the one hand, a quasi-linear analysis of an
elastoplastic model based on the infinitesimal strain theory was performed in
Ref. Otsuki . In this analysis, a non-uniform plastic deformation remains
after the external oscillation owing to the boundary condition at the walls of
the container. The plastic deformation causes an asymmetric stress
distribution in the direction of the oscillation, and the tension increases in
some parts of the paste. The desiccation after the oscillation enhances the
asymmetric stress, which creates cracks perpendicular to the direction of the
external oscillation. Therefore, in this analysis, the crack patterns are
asymmetric. On the other hand, in Ref. Ooshida08 ; Ooshida09 , a nonlinear
analysis based on finite strain theory is performed. This analysis shows that
a uniform plastic deformation under periodic boundary conditions increases the
tension over the entire area of the paste. This uniform increase of the
tension leads to uniform and symmetric crack patterns along the direction of
the oscillation.

However, the seemingly incoherent predictions may represent two different
limiting cases, as is suggested by a recent experiment Nakahara18 , where the
symmetric and asymmetric crack patterns occur depending on the number of
cycles in the external oscillation. For a few cycles, the asymmetric crack
pattern predicted in the quasi-linear analysis was formed. For a sufficiently
large number of cycles, the symmetric pattern predicted in the nonlinear
analysis was formed. However, the present theories cannot explain the
dependence of the crack pattern on the number of cycles. Therefore, a unified
theoretical model is required.

In this study, we formulate such a unified model, with which we numerically
examine the residual stress in a paste after external oscillation. In Sec. 2,
we review a Lagrangian description of kinematics. In Sec. 3, we derive the
time evolution equations for a paste. The memory effect of the external
oscillation is numerically investigated in Sec. 4. The setup of our simulation
is shown in Sec. 4.1. We demonstrate the relation between the plastic
deformation and the residual stress in Sec. 4.2. In Secs. 4.3 and 4.4, we
present the residual stress after the oscillation for different numbers of
cycles, respectively. The dependence of the stress distribution on the number
of cycles is examined in Sec. 4.5. We discuss our results and conclude our
paper in Sec. 5. In A, the left Cauchy–Green tensor in our setup is derived.
The stretching tensor is obtained in B. The constitutive relation of the paste
is derived in C. In D, we discuss the shear stress under steady shear in our
model.

##  2 Kinematics

The configuration of a 3-dimensional paste is represented by the current
coordinates $\bm{r}$ mapped from the reference coordinates $\bm{X}$ as

$\displaystyle(\bm{X},t)\mapsto\bm{r}(\bm{X},t)=\left[\begin{array}[]{ccc}x(X,Y,Z,t)\\\
y(X,Y,Z,t)\\\ z(X,Y,Z,t)\\\ \end{array}\right]_{\rm C},$ (4)

where $[\cdot]_{\rm C}$ denotes the representation in terms of Cartesian
components. The velocity $\bm{v}$ is given by

$\displaystyle\bm{v}(\bm{X},t)=\partial_{t}\bm{r}(\bm{X},t),$ (5)

and the acceleration is given by $\partial_{t}^{2}\bm{r}$. Here,
$\partial_{t}$ stands for the time derivative in the Lagrangian description,
and the independent variables of the paste are $\bm{X}$ and $t$. Therefore,
$\partial_{t}$ represents the derivative with respect to $t$, with $\bm{X}$
being constant.

The line element $d\bm{r}$ between two points labeled as $\bm{X}$ and
$\bm{X}+d\bm{X}$ with an infinitesimal vector $d\bm{X}$ is given by

$\displaystyle d\bm{r}=(\partial_{i}\bm{r})dX^{i},$ (6)

where $\partial_{i}=\partial/\partial X^{i}$, and we use the Einstein notation
for the summation of repeated indices. The square of the distance between the
points in the current configuration is

$(ds)^{2}=|(\partial_{i}\bm{r})dX^{i}|^{2}=g_{ij}dX^{i}dX^{j},$ (7)

where $g_{ij}$ is a component of the Euclidean metric tensor $\bm{g}$ given by

$g_{ij}=(\partial_{i}\bm{r})\cdot(\partial_{j}\bm{r}).$ (8)

Assuming the existence of a local stress-free “natural state”, the line
element between the two points in this state is denoted as
$d\bm{r}^{\natural}$ Ooshida08 . The square of the distance between the points
in the “natural state” is given by

$(ds^{\natural})^{2}=|d\bm{r}^{\natural}|^{2}=g_{ij}^{\natural}dX^{i}dX^{j},$
(9)

where $g_{ij}^{\natural}$ is the natural metric tensor satisfying
$g_{ij}^{\natural}=g_{ji}^{\natural}$. The change in the natural metric tensor
represents the plastic deformation.

The elastic deformation is represented by the transformation from
$d\bm{r}^{\natural}$ to $d\bm{r}$:

$d\bm{r}=\mathbf{F}d\bm{r}^{\natural},$ (10)

where $\mathbf{F}$ is the deformation gradient tensor. The left Cauchy–Green
Tensor is defined as Romano ; Truesdell

$\mathbf{B}=\mathbf{F}\mathbf{F}^{T}.$ (11)

$\mathbf{B}$ is given by

$\mathbf{B}=g^{ij}_{\natural}(\partial_{i}\bm{r})\otimes(\partial_{j}\bm{r}),$
(12)

where $(g^{ij}_{\natural})$ denotes the inverse of the component matrix of the
metric tensor $(g_{ij}^{\natural})$ satisfying

$g^{\natural}_{ij}g_{\natural}^{jk}=\delta_{i}^{k}.$ (13)

See A for the derivation of Eq. (12).

The difference in velocities between two points labeled as $\bm{X}$ and
$\bm{X}+d\bm{X}$ is given by

$d\bm{v}=\left(\partial_{i}\bm{v}\right)dX^{i}.$ (14)

The time rate of the deformation is described by the velocity gradient tensor
$\mathbf{W}$, which satisfies Romano ; Truesdell

$d\bm{v}=\mathbf{W}d\bm{r}.$ (15)

The stretching tensor is defined as

$\mathbf{D}=\frac{1}{2}\left(\mathbf{W}+\mathbf{W}^{T}\right).$ (16)

As shown in B, the stretching tensor is represented as

$\mathbf{D}=-\frac{1}{2}\left(\partial_{t}g^{ij}\right)(\partial_{i}\bm{r})\otimes(\partial_{j}\bm{r}),$
(17)

where $(g^{ij})$ is the inverse of the component matrix of the metric tensor
$(g_{ij})$ satisfying

$g_{ij}g^{jk}=\delta_{i}^{k}.$ (18)

##  3 Equation of motion and constitutive relation

The equation of motion for $\bm{r}(\bm{X},t)$ is given by

$\rho\partial^{2}_{t}\bm{r}={\rm div}\bm{\sigma}+\bm{f}_{\rm e},$ (19)

where $\rho$, $\bm{\sigma}$, and $\bm{f}_{\rm e}$ are the density in the
current configuration, the Cauchy stress tensor, and the body force,
respectively. The Cauchy stress tensor consists of an elastoplastic part
$\bm{\sigma}^{\rm(EP)}$ and a viscous part $\bm{\sigma}^{\rm(V)}$ as

$\bm{\sigma}=\bm{\sigma}^{\rm(EP)}+\bm{\sigma}^{\rm(V)}.$ (20)

The elastoplastic part is given by the constitutive relation with the strain
energy per unit volume in the reference configuration $\Sigma_{\rm R}$ as

$\bm{\sigma}^{\rm(EP)}=\frac{1}{J}\frac{\partial\Sigma_{\rm
R}}{\partial\tilde{\mathbf{F}}}\ \tilde{\mathbf{F}}^{T},$ (21)

where $\tilde{\mathbf{F}}=\partial\bm{r}/\partial\bm{X}$ is the “apparent”
deformation gradient tensor characterizing the transformation from the
reference configuration to the current configuration Beatty . Here,
$J=\det\tilde{\mathbf{F}}=\sqrt{\det\bm{g}}$ is the Jacobian. The strain
energy per unit volume in the current coordinates $\Sigma$ is related to
$\Sigma_{\rm R}$ as

$\displaystyle\Sigma_{\rm R}$ $\displaystyle=$ $\displaystyle J\Sigma.$ (22)

In this study, we apply the Hadamard strain energy destrade :

$\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\left\\{\mu(I_{1}-3)+\Psi(I_{3})\right\\}/2,$ (23)
$\displaystyle\Psi(I_{3})$ $\displaystyle=$
$\displaystyle(\lambda+\mu)(I_{3}-1)-2(\lambda+2\mu)(\sqrt{I_{3}}-1),$ (24)

where $I_{1}=\rm{tr}\mathbf{B}$ and $I_{3}=\rm{det}\mathbf{B}$ are the
rotational invariants of the left Cauchy–Green strain tensor $\mathbf{B}$ with
the Lamé constants $\mu$ and $\lambda$. Note that $I_{1}$ and $I_{3}$ depend
on $\tilde{\mathbf{F}}$, as shown in C. In Ref. Ooshida08 , the incompressible
neo-Hookean model is adopted for the strain energy, but we assume Hadamard
strain energy with small compressibility (i.e. large but finite $\lambda$) to
avoid difficulties in numerical simulations. Note that the neo-Hookean strain
energy is the incompressible limit of the Hadamard strain energy. Substituting
Eq. (22) with Eqs. (23) and (24) into Eq. (21), we obtain

$\bm{\sigma}^{\rm(EP)}=\\{\Sigma+I_{3}\Psi^{\prime}(I_{3})+\mu\\}\mathbf{I}+\mu\left(\mathbf{B}-\mathbf{I}\right).$
(25)

Here, $\mathbf{I}$ is the unit tensor. See C for the derivation of Eq. (25).
For $\bm{\sigma}^{\rm(V)}$, we adopt the linear viscous stress tensor as

$\bm{\sigma}^{\rm(V)}=2\eta\mathbf{D}$ (26)

with the viscosity $\eta$ Romano .

Following Ref. Ooshida08 ; Ooshida09 , the plastic deformation is described by
the temporal evolution of the natural metric $g_{\natural}^{ij}$ as

$(1+\tau\partial_{t})g_{\natural}^{ij}=Kg^{ij}.$ (27)

Here, $K$ is given by

$K=\frac{3}{g^{\natural}_{ij}g^{ij}}=\frac{3}{I_{1}}$ (28)

to satisfy the incompressibility condition
$J^{\natural}=\sqrt{\det\bm{g}^{\natural}}=1$ in the natural state. It should
be noted that we assume effective incompressibility for $\bm{g}$ by adopting
Poisson’s ratio $\nu\equiv\lambda/(\lambda+\mu)/2\simeq 1/2$ in our numerical
simulations. The inverse of the relaxation time $\tau$ is given by

$\tau^{-1}(\sigma_{\rm e})=\frac{1}{\tau_{0}}{\rm
max}\left(0,1-\frac{\sigma_{\rm Y}}{\sigma_{\rm e}}\right)$ (29)

with characteristic time $\tau_{0}$ and tensile yield strength $\sigma_{\rm
Y}$. The equivalent tensile stress $\sigma_{\rm e}$ is given by

$\displaystyle\sigma_{\rm e}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\sigma_{xx}^{\rm(EP)}-\sigma_{yy}^{\rm(EP)})^{2}$
(30)
$\displaystyle+\frac{1}{2}(\sigma_{yy}^{\rm(EP)}-\sigma_{zz}^{\rm(EP)})^{2}$
$\displaystyle+\frac{1}{2}(\sigma_{zz}^{\rm(EP)}-\sigma_{xx}^{\rm(EP)})^{2},$
$\displaystyle+3\left\\{\left(\sigma_{xz}^{\rm(EP)}\right)^{2}+\left(\sigma_{yz}^{\rm(EP)}\right)^{2}+\left(\sigma_{zx}^{\rm(EP)}\right)^{2}\right\\}.$

Equations (27) and (29) indicate that the plastic deformation associated with
the change of the natural metric $g_{\natural}^{ij}$ occurs when the von Mises
yield criterion $\sigma_{\rm e}=\sigma_{\rm Y}$ is satisfied Jones . See D for
the flow curve of our model under steady shear.

##  4 Memory of external oscillation

In this section, we numerically investigate the memory effect of the external
oscillation. In Sec. 4.1, we explain our setup. In Sec. 4.2, we demonstrate
how plastic deformation affects the residual stress. The residual stress for
different numbers of cycles is shown in Secs. 4.3 and 4.4. In Sec. 5, we
discuss the dependence of the symmetry of the stress distribution on the
number of cycles.

###  4.1 Setup

Let us consider a paste of thickness $H$ and width $2L$ in a container, as
shown in Fig. 1. The center of the container is at $X=0$, and the bottom is at
$Z=0$. We assume plane strain deformation so that the current configuration is
given by

$\displaystyle\bm{r}(X,Z,t)=\left[\begin{array}[]{ccc}X+u(X,Z,t)\\\ Y\\\
Z+w(X,Z,t)\\\ \end{array}\right]_{\rm C}$ (34)

with displacements $u$ and $w$ in the $X$ and $Z$ directions, respectively.
The body force is given by

$\displaystyle\bm{f}_{\rm e}=-\rho\left[\begin{array}[]{c}A(t)\\\ 0\\\ G\\\
\end{array}\right]_{\rm C},$ (38)

where $G$ is the gravitational acceleration and $A(t)$ is the horizontal
acceleration due to the external oscillation of the container.

Figure 1: Schematic illustration of a paste in a container.

Assuming that the incompressibility of the natural state
($\det\mathbf{g}^{\natural}=1$), $\mathbf{g}^{\natural}$ can be expressed by
two parameters, and we set

$\displaystyle\mathbf{g}^{\natural}=\left[\begin{array}[]{ccc}e^{-\alpha}&0&\beta\\\
0&1&0\\\ \beta&0&e^{\alpha}(1+\beta^{2})\\\ \end{array}\right]$ (42)

and

$\displaystyle\mathbf{g}_{\natural}=\left[\begin{array}[]{ccc}e^{\alpha}(1+\beta^{2})&0&-\beta\\\
0&1&0\\\ -\beta&0&e^{-\alpha}\\\ \end{array}\right],$ (46)

where $\beta$ represents the plastic shear strain, whereas the plastic normal
strain is characterized by $\alpha$. $e^{-\alpha}$ represents a normal
component of the natural metric tensor in the direction of the external
oscillation.

The initial configuration is given by $u=w=\alpha=\beta=0$ at $t=-T_{\rm I}$
with $g_{ij}^{\natural}=\delta_{ij}$. For $-T_{\rm I}\leq t<0$, we set
$A(t)=0$ and relax the configuration under gravitational acceleration. We
apply the external oscillation as

$\displaystyle A(t)=-A_{x}\cos\frac{2\pi t}{T}$ (47)

for $0\leq t<N_{\rm c}T$ with the maximum acceleration $A_{x}$, the period
$T$, and the number of cycles as $N_{\rm c}$. After $N_{\rm c}$ cycles of the
oscillation, we set $A(t)=0$ and relax the system until $t=T_{\rm w}$. The
profiles of the stress and plastic deformation shown below are obtained at
$t=T_{\rm w}$.

We assume the following no-slip boundary conditions at the bottom of the
container:

$\displaystyle u|_{Z=0}=w|_{Z=0}=0,$ (48)

while the stress applied to the free surface at the top of the paste is given
by

$\displaystyle\sigma_{xz}|_{Z=H}=0,\ \ \sigma_{zz}|_{Z=H}=-p_{0},$ (49)

where $p_{0}$ is the atmospheric pressure. For the lateral walls, we assume
that the paste does not leave the wall as

$\displaystyle u|_{X=\pm L}=0,$ (50)

but it freely slides on the wall as

$\displaystyle\sigma_{zx}|_{X=\pm L}=0.$ (51)

In this study, we use the unit mass, length, and time as $m=\rho H^{3}$,
$l=H$, and $\tau=\sqrt{H/G}$, respectively. The parameter values are set as
$L/H=10$, $A_{x}/G=6.0\times 10^{-2}$, $T/\tau=31.25$,
$\nu=\lambda/(\lambda+\mu)/2=0.4999$, $\tau_{0}/\tau=0.94$, $\mu/(\rho
HG)=7.3\times 10^{-2}$, $\sigma_{\rm Y}/(\rho HG)=4.7\times 10^{-2}$,
$\eta/(\rho\sqrt{GH^{3}})=0.14$, $T_{I}/\tau=100$, and $T_{w}/\tau=300$ based
on experiments Nakahara05 , except for $\mu$ and $\eta$. It should be noted
that $\tau_{0}$ is estimated from the flow curve under steady shear, which is
shown in D. The atmospheric pressure $p_{0}$ is set to $0$ because the
numerical results shown below do not depend on the value of $p_{0}$. We have
checked that our numerical results do not depend on the value of Poisson’s
ratio $\nu\equiv\lambda/(\lambda+\mu)/2$ for $\nu\geq 0.4999$. We adopt the
finite-difference method with the time interval as $\Delta
t/\sqrt{H/G}=5.0\times 10^{-5}$ and the spatial mesh size as $\Delta
x/(2L)=\Delta z/(2L)=0.025$.

###  4.2 Effect of plastic deformation on residual stress

In this section, we discuss the normal component of the deviatoric stress

$\displaystyle s_{xx}=\sigma_{xx}-(\sigma_{xx}+\sigma_{yy}+\sigma_{zz})/3$
(52)

in the direction of the oscillation because its increase is essential for the
formation of cracks. In the linear approximation of $\partial_{i}u$, $\alpha$,
and $\beta$ while ignoring $\partial_{i}w$, $s_{xx}$ is given by

$\displaystyle s_{xx}=s_{xx}^{\rm(L)}+s_{xx}^{\rm(N)}$ (53)

with

$\displaystyle s_{xx}^{\rm(L)}=2\mu u_{X}$ (54)

and

$\displaystyle s_{xx}^{\rm(N)}=2\mu\alpha.$ (55)

Here, we abbreviate $\partial u/\partial X$ as $u_{X}$. Finite $u_{X}$ due to
the plastic deformation leads to the change of $s_{xx}$ through Eqs. (53) and
(54), which corresponds to the mechanism of the memory effect assumed in the
quasi-linear analysis Otsuki . The plastic deformation characterized by
$\alpha$ induces the increase of $s_{xx}$ through Eqs. (53) and (55), which is
consistent with the scenario proposed in the nonlinear analysis Ooshida08 ;
Ooshida09 .

###  4.3 Residual stress for $N_{c}=1/2$

In Fig. 2, we plot $s_{xx}$ as a function of $(X,Z)$ after the oscillation
with $N_{c}=1/2$ (i.e. half a cycle). Here, $(X,Z)$ indicates the initial
position. During the oscillation, the equivalent tensile stress $\sigma_{\rm
e}$ exceeds the tensile yield strength $\sigma_{\rm Y}$, and plastic
deformation occurs at the bottom of the paste. The residual stress remains as
a memory effect due to the plastic deformation. The residual stress
distribution is almost symmetric with respect to the inversion of $X$ (i.e.,
$X\to-X$) and positive near the bottom, while it is asymmetric near the
surface of the paste. The asymmetric pattern near the surface is consistent
with the stress distribution predicted in the quasi-linear analysis Otsuki ,
which induces the asymmetric crack patterns observed in experiments with a few
cycles of oscillation Nakahara18 .

Figure 2: $s_{xx}$ as a function of $(X,Z)$ after the external oscillation
with $N_{c}=1/2$.

Figure 3 displays $\alpha$ as a function of $(X,Z)$ after the external
oscillation with $N_{c}=1/2$. $\alpha$ is positive near the center of the
bottom because the equivalent tensile stress $\sigma_{\rm e}$ at the center of
the bottom maximizes and exceeds the yield stress $\sigma_{Y}$ under the
horizontal force owing to the external oscillation. The symmetric $\alpha$ in
Eq. (55) explains the nearly symmetric distribution of $s_{xx}$ near the
bottom as shown in Fig. 2.

Figure 3: $\alpha$ as a function of $(X,Z)$ after the external oscillation
with $N_{c}=1/2$.

In Fig. 4 (a), we show $\beta$ as a function of $(X,Z)$ after the external
oscillation for $N_{c}=1/2$. For $N_{c}=1/2$, the external force is applied in
the negative $X$-axis direction before the external oscillation is stopped,
which leads to the negative shear strain $\partial_{Z}u<0$ getting associated
with the plastic shear strain $\beta<0$ near the center of the bottom, as
shown in Fig. 4 (a). The plastic shear strain $\beta$ causes negative
displacement $u<0$ in the $X$-axis direction except for the boundaries, which
results in the asymmetric $u_{X}$ as shown in 4 (b). This asymmetric $u_{X}$
in Eq. (54) explains the asymmetric residual stress near the surface in Fig.
2.

Figure 4:  (a) $\beta$ as a function of $(X,Z)$ after the external oscillation
with $N_{c}=1/2$. (b) $u_{X}$ as a function of $(X,Z)$ after the external
oscillation with $N_{c}=1/2$.

###  4.4 Residual stress for $N_{c}=2$

In Fig. 5, we plot $s_{xx}$ as a function of $(X,Z)$ after the external
oscillation with $N_{c}=2$ (i.e. two cycles). $s_{xx}$ near the bottom for
$N_{c}=2$ is higher than that for $N_{c}=1/2$. The distribution of $s_{xx}$
becomes almost symmetric even near the surface, which is consistent with the
symmetric crack patterns shown in experiments with sufficient cycles of the
external oscillation Nakahara18 .

Figure 5: $s_{xx}$ as a function of $(X,Z)$ after the external oscillation
with $N_{c}=2$.

Figure 6 displays $\alpha$ as a function of $(X,Z)$ after the oscillation with
$N_{c}=2$. As shown in D, $\alpha$ monotonically increases under shear.
Therefore, $\alpha$ for $N_{c}=2$ becomes larger than that for $N_{c}=1/2$,
which leads to an increase in $s_{xx}$ near the bottom, as shown in Fig. 5.

Figure 6: $\alpha$ as a function of $(X,Z)$ after the external oscillation
with $N_{c}=2$.

In Fig. 7 (a), we plot $\beta$ for $N_{c}=2$ as a function of $(X,Z)$. For
$N_{c}=2$, the plastic deformation is accumulated during the oscillation and
the distribution of $\beta$ becomes nearly antisymmetric, which leads to the
symmetric $u_{X}$ as shown in 4 (b). This symmetric $u_{X}$ explains the
symmetric stress distribution near the surface, as shown in Fig. 5.

Figure 7:  (a) $\beta$ as a function of $(X,Z)$ after the external oscillation
with $N_{c}=2$. (b) $u_{X}$ as a function of $(X,Z)$ after the external
oscillation with $N_{c}=2$.

###  4.5 Dependence of asymmetry on $N_{c}$

The asymmetric part of the deviatoric stress is given by

$\displaystyle
s_{xx}^{\rm(a)}(X,Z)=\left\\{s_{xx}(X,Z)-s_{xx}(-X,Z)\right\\}/2.$ (56)

Here, we define a parameter

$A=\frac{\int_{-L}^{L}dX\int_{0}^{H}dZ|s_{xx}^{\rm(a)}|^{2}}{\int_{-L}^{L}dX\int_{0}^{H}dZ|s_{xx}|^{2}}$
(57)

characterizing the asymmetry of the stress distribution. In Fig. 8, we show
$A$ after the external oscillation as a function of $N_{c}$. $A$ decreases
with increasing $N_{c}$, which is consistent with the symmetry change of the
crack patterns in the experiments Nakahara18 . Figure 8 indicates that the
nonlinear effect becomes dominant for $N_{c}\geq 1$. A similar dependence is
obtained for different values of $\mu$. The sign of the residual stress in the
quasi-linear analysis is reversed when the direction of the external force is
reversed Otsuki , while the stress field is independent of the direction of
the external force in the nonlinear analysis Ooshida08 ; Ooshida09 . Hence, we
consider that the quasi-linear effect is gradually canceled when the external
force is applied in both directions for $N_{c}\geq 1$, which leads to the fast
decrease of $A$ in Fig. 8.

Figure 8: $A$ as a function of $N_{c}$ after the external oscillation.

##  5 Summary and discussions

In this study, we numerically investigated the residual stress of a paste
after an external oscillation based on an elastoplastic model. The residual
stress remains as a memory of the oscillation, which leads to crack patterns
perpendicular to its direction Nakahara05 ; Nakahara06a . The residual stress
distribution is asymmetric when the number of cycles in the oscillation is
small. The symmetry of the residual stress is enhanced by increasing the
number of cycles, which is consistent with the results of the experiments in
Ref. Nakahara18 .

The plastic deformation induced by the oscillation remains until the formation
of cracks Goehring . The desiccation after the oscillation enhances the
residual stress caused by the plastic deformation as shown in Ref.
Kitsunezaki16 ; Kitsunezaki17 . Thus, we expect that the dependence of the
stress symmetry on the number of cycles can be detected in an experiment using
the method in Ref. Kitsunezaki16 ; Kitsunezaki17 . Such an experiment will
verify the validity of our theory.

Different crack patterns appear depending on the types of powders and ways of
applying external forces Nakahara06b ; Matsuo ; Nakayama ; Akiba , which is
not explained by the present theory. In our analysis, the process of crack
formation and the displacement in the $y$ direction are neglected because we
have restricted our attention to the residual stress leading to the
perpendicular crack when an external force is applied in one direction. An
extension of the models of drying crack patterns Kitsunezaki ; Ito14a ; Ito14b
; Halasz  incorporating the effect of the plastic deformation discussed in
this study will help understand these crack patterns.

Recent experiments using micro-focus X-ray computerized tomography have
revealed that the configuration of microscopic particles in a paste becomes
anisotropic after the external oscillation Kitsunezaki17a . This anisotropy is
considered as a microscopic memory effect. The relationship between the
microscopic memory effect and the macroscopic residual stress discussed in
this study is unclear. An extension of microscopic theories such as the mode-
coupling theory Miyazaki ; Fuchs ; Ballauff ; Mohan ; Fritschi ; Mohan17 ;
Moghimi  and the pair distribution function theory Otsuki06 , predicting the
macroscopic constitutive relation, may clarify this relationship.

###### Acknowledgements.

The authors thank A. Nakahara, S. Kitsunezaki, R. Tarumi, H. Hayakawa, and T.
Ooshida for fruitful discussions. This work was supported by JSPS KAKENHI
Grant Numbers JP16H04025, JP19K03670, and JP21H01006.

## Author contribution statement

J.M. carried out the simulations. Both authors developed the theory and wrote
the manuscript.

##  Appendix A Left Cauchy–Green Tensor

In this Appendix, we derive Eq. (12) using the method used in Ref. Ooshida08 .
We first represent the natural metric tensor $g_{ij}^{\natural}$ as a dot
product of orthogonal vectors. $\bm{g}^{\natural}$, satisfying Eq. (9), is a
positive definite symmetric matrix with positive eigenvalues
$\lambda_{1}^{2}$, $\lambda_{2}^{2}$, and $\lambda_{3}^{2}$, whose
corresponding eigenvectors are $\bm{q}_{1}$, $\bm{q}_{2}$, and $\bm{q}_{3}$,
respectively. A matrix $\mathbf{Q}=[\bm{q}_{1}\ \bm{q}_{2}\ \bm{q}_{3}]$ is
given by

$\displaystyle\bm{g}^{\natural}\mathbf{Q}=\mathbf{Q}\left[\begin{array}[]{ccc}\lambda_{1}^{2}&0&0\\\
0&\lambda_{2}^{2}&0\\\ 0&0&\lambda_{3}^{2}\\\ \end{array}\right]$ (61)

and

$\displaystyle\mathbf{Q}^{-1}=\mathbf{Q}^{T},$ (62)

which indicates

$\displaystyle\bm{g}^{\natural}=\mathbf{Q}\left[\begin{array}[]{ccc}\lambda_{1}^{2}&0&0\\\
0&\lambda_{2}^{2}&0\\\ 0&0&\lambda_{3}^{2}\\\
\end{array}\right]\mathbf{Q}^{T}.$ (66)

It should be noted that $\mathbf{Q}$ is an orthogonal matrix. Here, we define
a symmetric matrix

$\displaystyle\mathbf{P}=\mathbf{Q}\left[\begin{array}[]{ccc}\lambda_{1}&0&0\\\
0&\lambda_{2}&0\\\ 0&0&\lambda_{3}\\\ \end{array}\right]\mathbf{Q}^{T}.$ (70)

The natural metric tensor $\bm{g}^{\natural}$ satisfies

$\displaystyle\bm{g}^{\natural}=\mathbf{P}\mathbf{P}.$ (71)

Here, we introduce independent vectors $\bm{p}_{1}$, $\bm{p}_{2}$, and
$\bm{p}_{2}$ satisfying

$\displaystyle\mathbf{P}=[\bm{p}_{1}\ \bm{p}_{2}\ \bm{p}_{3}].$ (72)

The component of $\bm{g}^{\natural}$ is given by

$\displaystyle g_{ij}^{\natural}=\bm{p}_{i}\cdot\bm{p}_{j}.$ (73)

From Eqs. (9) and (73), we obtain

$\displaystyle d\bm{r}^{\natural}=\bm{p}_{i}dX^{i}.$ (74)

Because $\bm{p}_{1}$, $\bm{p}_{2}$, and $\bm{p}_{3}$ are independent, there
exist dual vectors $\bm{p}_{*}^{1}$, $\bm{p}_{*}^{2}$, and $\bm{p}_{*}^{3}$
satisfying $\bm{p}_{i}\cdot\bm{p}_{*}^{j}=\delta_{i}^{j}$. It should be noted
that the inverse of the component matrix of the metric tensor
$(g_{ij}^{\natural})$ satisfying Eq. (13) is represented by

$\displaystyle g^{ij}_{\natural}=\bm{p}_{*}^{i}\cdot\bm{p}_{*}^{j}.$ (75)

Taking a dot product of Eq. (74) and $\bm{p}_{*}^{i}$, we find

$\displaystyle dX^{i}=\bm{p}_{*}^{i}\cdot d\bm{r}^{\natural}.$ (76)

Substituting this equation into Eq. (6), we obtain

$\displaystyle d\bm{r}=(\partial_{i}\bm{r})\bm{p}_{*}^{i}\cdot
d\bm{r}^{\natural}=\left\\{(\partial_{i}\bm{r})\otimes\bm{p}_{*}^{i}\right\\}d\bm{r}^{\natural}$
(77)

From Eqs. (10) and (77), we find

$\displaystyle\mathbf{F}=(\partial_{i}\bm{r})\otimes\bm{p}_{*}^{i}.$ (78)

Substituting this equation into Eq. (11) with Eq. (75), we obtain Eq. (12).

##  Appendix B Stretching Tensor

In this Appendix, we derive Eq. (17). First, we introduce

$\displaystyle\nabla X^{i}=\frac{\partial X^{i}}{\partial\bm{r}},$ (79)

which satisfies

$\displaystyle(\partial_{i}\bm{r})\cdot\nabla X^{j}=\delta_{i}^{j}$ (80)

and

$\displaystyle(\partial_{i}\bm{r})\otimes\nabla X^{i}=\mathbf{I}.$ (81)

A component of the inverse matrix of $\bm{g}$ satisfying Eq. (18) is given by
the dot product of $\left\\{\nabla X^{i}\right\\}$:

$\displaystyle g^{ij}=\nabla X^{i}\cdot\nabla X^{j}.$ (82)

Using Eq. (81) with Eq. (82), we obtain

$\displaystyle\nabla X^{i}$ $\displaystyle=$ $\displaystyle\mathbf{I}\nabla
X^{i}$ (83) $\displaystyle=$
$\displaystyle\left\\{(\partial_{j}\bm{r})\otimes\nabla X^{j}\right\\}\nabla
X^{i},$ $\displaystyle=$ $\displaystyle\left\\{\nabla X^{i}\cdot\nabla
X^{j}\right\\}(\partial_{j}\bm{r}),$ $\displaystyle=$ $\displaystyle
g^{ij}\partial_{j}\bm{r}.$

Substituting Eq. (83) into Eq. (81),

$\displaystyle
g^{ij}(\partial_{i}\bm{r})\otimes(\partial_{j}\bm{r})=\mathbf{I}.$ (84)

Taking the dot product of Eq. (6) with $\nabla X^{i}$, we obtain

$\displaystyle dX^{i}=\nabla X^{i}\cdot d\bm{r}.$ (85)

Substituting this equation into Eq. (14), we obtain

$\displaystyle d\bm{v}=(\partial_{i}\bm{v})(\nabla X^{i}\cdot
d\bm{r})=\left\\{(\partial_{i}\bm{v})\otimes\nabla X^{i}\right\\}d\bm{r}.$
(86)

Comparing this equation with Eq. (15), we find

$\displaystyle\bm{W}=(\partial_{i}\bm{v})\otimes\nabla X^{i}$ (87)

Substituting Eq. (87) into Eq. (16) with Eq. (87), we obtain

$\displaystyle\bm{D}$ $\displaystyle=$
$\displaystyle\frac{1}{2}g^{ij}\left\\{(\partial_{i}\bm{v})\otimes(\partial_{j}\bm{r})+(\partial_{i}\bm{r})\otimes(\partial_{j}\bm{v})\right\\}$
(88) $\displaystyle=$
$\displaystyle\frac{1}{2}g^{ij}\partial_{t}\left\\{(\partial_{i}\bm{r})\otimes(\partial_{j}\bm{r})\right\\}$

Differentiating Eq. (84) by $t$, we find

$g^{ij}\partial_{t}\left\\{(\partial_{i}\bm{r})\otimes(\partial_{j}\bm{r})\right\\}=-\left(\partial_{t}g^{ij}\right)\left\\{(\partial_{i}\bm{r})\otimes(\partial_{j}\bm{r})\right\\}.$
(89)

Substituting this equation into Eq. (88), we obtain Eq. (17).

##  Appendix C Elasto-plastic part of stress tensor

In this appendix, we derive Eq. (25). Substituting Eq. (22) into Eq. (21), we
obtain

$\displaystyle\bm{\sigma}^{\rm(EP)}=\frac{\Sigma}{J}\frac{\partial
J}{\partial\tilde{\mathbf{F}}}\
\tilde{\mathbf{F}}^{T}+\frac{\partial\Sigma}{\partial\tilde{\mathbf{F}}}\
\tilde{\mathbf{F}}^{T}.$ (90)

Using a formula

$\displaystyle\frac{\partial\det{\mathbf{A}}}{\partial\mathbf{A}}\
{\mathbf{A}}^{T}=\left(\det{\mathbf{A}}\right)\mathbf{I}$ (91)

for a tensor $\mathbf{A}$, $J=\det\tilde{\mathbf{F}}$ satisfies

$\displaystyle\frac{\Sigma}{J}\frac{\partial J}{\partial\tilde{\mathbf{F}}}\
\tilde{\mathbf{F}}^{T}=\Sigma\mathbf{I}.$ (92)

The strain energy per unit volume in the current coordinates $\Sigma$ is a
function of the rotational invariants $I_{1}$ and $I_{3}$ of the left
Cauchy–Green strain tensor $\mathbf{B}$. Hence, the Cartesian component of the
second term in Eq. (90) satisfies

$\displaystyle\left(\frac{\partial\Sigma}{\partial\tilde{\mathbf{F}}}\
\tilde{\mathbf{F}}^{T}\right)_{ij}=\left\\{\frac{\partial\Sigma}{\partial
I_{1}}\frac{\partial I_{1}}{\partial B_{lm}}+\frac{\partial\Sigma}{\partial
I_{3}}\frac{\partial I_{3}}{\partial B_{lm}}\right\\}\frac{\partial
B_{lm}}{\partial\tilde{F}_{ik}}\ \tilde{F}_{jk}.$ (93)

Here, we obtain

$\displaystyle\frac{\partial I_{1}}{\partial\mathbf{B}}=\mathbf{I}$ (94)

and

$\displaystyle\frac{\partial I_{3}}{\partial\mathbf{B}}\
=I_{3}{\mathbf{B}}^{-1}$ (95)

from Eq. (91) and $\mathbf{B}^{T}=\mathbf{B}$. Using Eqs. (11), (75), and
(78), we find

$\displaystyle\mathbf{B}=\tilde{\mathbf{F}}{\bm{g}}_{\natural}\tilde{\mathbf{F}}^{T}.$
(96)

From Eq. (96), we obtain

$\displaystyle\frac{\partial B_{lm}}{\partial\tilde{F}_{ik}}\
\tilde{F}_{jk}=\delta_{il}B_{jm}+\delta_{im}B_{lj}.$ (97)

Substituting Eqs. (94), (95), and (97) into Eq. (93), we derive

$\displaystyle\frac{\partial\Sigma}{\partial\tilde{\mathbf{F}}}\tilde{\mathbf{F}}^{T}=2\left\\{\frac{\partial\Sigma}{\partial
I_{1}}\mathbf{B}+I_{3}\frac{\partial\Sigma}{\partial
I_{3}}\mathbf{I}\right\\}.$ (98)

Substituting Eqs. (92) and (98) into Eq. (90), we obtain

$\displaystyle\bm{\sigma}^{\rm(EP)}=\Sigma\mathbf{I}+2\left\\{\frac{\partial\Sigma}{\partial
I_{1}}\mathbf{B}+I_{3}\frac{\partial\Sigma}{\partial
I_{3}}\mathbf{I}\right\\}.$ (99)

Substituting Eq. (23) into this equation, we derive Eq. (25).

##  Appendix D Stress under uniform steady shear

In this appendix, we discuss the rheological properties described by Eqs. (25)
and (27) with Eq. (12) under simple steady shear:

$\displaystyle\bm{r}(X,Y,Z,t)=\left[\begin{array}[]{ccc}X+\gamma(t)Z\\\ Y\\\
Z\\\ \end{array}\right]_{\rm C},$ (103)

where the shear strain is given by

$\displaystyle\gamma(t)=\dot{\gamma}t$ (104)

with the shear rate $\dot{\gamma}$. The natural metric tensor
$\mathbf{g}^{\natural}$ is represented by Eq. (42) with $\alpha=\beta=0$ at
$t=0$.

Figure 9:  (a) $\sigma_{xz}^{\rm(EP)}$ against the shear strain $\gamma$ for
$\sigma_{\rm Y}/\mu=1.0$ with different $\dot{\gamma}\tau_{0}$. (b)
$\sigma^{\rm(EP)}_{xz}$ against the shear strain rate $\dot{\gamma}$ for
$\sigma_{\rm Y}/\mu=1.0$ with different $\dot{\gamma}\tau_{0}$ in the steady
state.

Substituting Eqs. (103) and (46) into Eqs. (25) and (27) with Eq. (12), we
obtain the time evolution equations for the shear stress
$\sigma_{xz}^{\rm(EP)}$ as

$\displaystyle\partial_{t}\sigma_{xz}^{\rm(EP)}$ $\displaystyle=$
$\displaystyle\mu e^{-\alpha}\dot{\gamma}-\tau^{-1}(\sigma_{\rm
e})\sigma_{xz}^{\rm(EP)}$ (105) $\displaystyle\partial_{t}\alpha$
$\displaystyle=$ $\displaystyle\tau^{-1}(\sigma_{\rm
e})\left(1-\frac{3}{I_{1}}e^{\alpha}\right)$ (106)
$\displaystyle\partial_{t}\beta$ $\displaystyle=$
$\displaystyle\tau^{-1}(\sigma_{\rm
e})\left(\frac{3}{I_{1}}\gamma-\beta\right)$ (107)

with

$\displaystyle
I_{1}=3+2(\cosh\alpha-1)+e^{\alpha}\left(\sigma_{xz}^{\rm(EP)}/\mu\right)^{2}.$
(108)

Here, $\sigma_{\rm e}$ is given by

$\displaystyle\sigma_{\rm
e}=\sqrt{3\left(\sigma_{xz}^{\rm(EP)}\right)^{2}+\Xi^{2}}$ (109)

with

$\displaystyle\Xi^{2}$ $\displaystyle=$
$\displaystyle\frac{\mu^{2}e^{2\alpha}}{2}\left\\{1+\left(\frac{\sigma_{xz}^{\rm(EP)}}{\mu}\right)^{2}-e^{-\alpha}\right\\}^{2}$
(110) $\displaystyle+\frac{\mu^{2}}{2}(e^{-\alpha}-1)^{2}$
$\displaystyle+\frac{\mu^{2}}{2}\left\\{e^{\alpha}+e^{\alpha}\left(\frac{\sigma_{xz}^{\rm(EP)}}{\mu}\right)^{2}-e^{-\alpha}\right\\}^{2}.$

It should be noted that

$\displaystyle\partial_{t}\alpha$ $\displaystyle\simeq$
$\displaystyle\tau^{-1}(\sigma_{\rm
e})\left\\{-\alpha+\frac{1}{3}\left(\frac{\sigma_{xz}^{\rm(EP)}}{\mu}\right)^{2}\right\\}$
(111)

for $\alpha\ll 1$ and $\sigma_{xz}^{\rm(EP)}\ll\mu$, which indicates that
$\alpha$ monotonically increases under steady shear.

In Fig. 9(a), we plot the shear stress $\sigma^{\rm(EP)}_{xz}$ against the
shear strain $\gamma$ for $\sigma_{\rm Y}/\mu=1.0$ with different
$\dot{\gamma}\tau_{0}$. For sufficiently small $\gamma$,
$\sigma^{\rm(EP)}_{xz}$ is almost proportional to $\gamma$. As $\gamma$
increases, $\sigma^{\rm(EP)}_{xz}$ exhibits a peak and converges to a steady
state. $\sigma^{\rm(EP)}_{xz}$ and $\gamma$ at the peak increase with
increasing $\dot{\gamma}$. This behavior is qualitatively similar to the
stress–strain curve in experiments Aken  and numerical simulations Varnik .

Figure 9(b) displays the shear stress $\sigma_{xz}^{\rm(EP)}$ against the
shear strain rate $\dot{\gamma}$ for $\sigma_{\rm Y}/\mu=1.0$ in the steady
state. $\sigma^{\rm(EP)}_{xz}$ monotonically increases and obeys
$\sigma_{xz}^{\rm(EP)}=\frac{\sigma_{\rm
Y}}{\sqrt{3}}+\mu\dot{\gamma}\tau_{0}$, which is consistent with the behavior
of a Bingham plastic fluid Bingham .

## References

  * (1) L. Goehring, A. Nakahara, T. Dutta, S. Kitsunezaki, S. Tarafdar, Deciccation cracks and their patterns: formation and modeling in science and nature (Wiley-VCH, Weinheim, 2015).
  * (2) A. Groisman and E. Kaplan, Europhys. Lett. 25, 415 (1994).
  * (3) A. Nakahara and Y. Matsuo, J. Phys. Soc. Jpn. 74, 1362 (2005).
  * (4) A. Nakahara and Y. Matsuo, J. Stat. Mech. 2006, P07016 (2006).
  * (5) A. Nakahara and Y. Matsuo, Phys. Rev. E 74, 045102(R) (2006).
  * (6) Y. Matsuo and A. Nakahara, J. Phys. Soc. Jpn. 81, 024801 (2012).
  * (7) H. Nakayama, Y. Matsuo, T. Ooshida, and A. Nakahara, Eur. Phys. J. E 36, 1 (2013).
  * (8) Y. Akiba and H. Shima, J. Phys. Soc. Jpn. 88, 024001 (2019).
  * (9) D. Mal, S. Sinha, T. R. Middya, and S. Tarafdar, Phys. A 384, 182 (2007).
  * (10) T. Khatun, T. Dutta, and S. Tarafdar, Langmuir 29, 15535 (2013).
  * (11) L. Pauchard, F. Elias, P. Boltenhagen, A. Cebers, J. C. Bacri, Phys. Rev. E 77, 021402 (2008).
  * (12) A. T. Ngo, J. Richardi, and M. P. Pileni, J. Phys. Chem. B 112, 14409 (2008).
  * (13) H. Lama, V. R. Dugyala, M. G. Basavaraj, and D. K. Satapathy, Phys. Rev. E 94, 012618 (2016).
  * (14) J. L. Jr. Beuth, Int. J. Solids Struct. 29, 1657 (1992).
  * (15) S. Kitsunezaki, Phys. Rev. E 60, 6449 (1999).
  * (16) K. B. Singh and M. S. Tirumkudulu, Phys. Rev. Lett. 98, 218302 (2007).
  * (17) W. Man and W. B. Russel, Phys. Rev. Lett. 100, 198302 (2008).
  * (18) S. Kitsunezaki, A. Nakahara, and Y. Matsuo, EPL 114, 64002 (2016).
  * (19) S. Kitsunezaki, Y. Matsuo, and A. Nakahara, J. Chem. Eng. Jpn. 50, 775 (2017).
  * (20) M. Otsuki, Phys. Rev. E 72, 046115 (2005).
  * (21) T. Ooshida, Phys. Rev. E 77, 061501 (2008).
  * (22) T. Ooshida, J. Phys. Soc. Jpn. 78, 104801 (2009).
  * (23) A. Nakahara, T. Hiraoka, R. Hayashi, Y. Matsuo, and S. Kitsunezaki, Phil. Trans. R. Soc. A 377, 20170395 (2018).
  * (24) A. Romano and A. Marasco, Continuum Mechanics using Mathematica: Fundamentals, Methods, and Applications (Brikhauser Basel, Basel, 2014).
  * (25) C. Truesdell and W. Noll, The Non-Linear Field Theories of Mechanics (Springer-Verlag Berlin Heidelberg, Berlin, 2004).
  * (26) M. F. Beatty, Appl. Mech. Rev. 40, 1699 (1987).
  * (27) M. Destrade and G. Saccomandi, Phys. Rev. E 72, 016620 (2005).
  * (28) R. M. Jones, Deformation Theory of Plasticity (Bull Ridge Corporation, Blacksburg, 2009).
  * (29) S.-i Ito and S. Yukawa, Phys. Rev. E 90, 042909 (2014).
  * (30) S.-i Ito and S. Yukawa, J. Phys. Soc. Jap. 83, 124005 (2014).
  * (31) Z. Halász, A. Nakahara, S. Kitsunezaki, and F. Kun, Phys. Rev. E 96, 033006 (2017).
  * (32) S. Kitsunezaki, A. Sasaki, A. Nishimoto, T. Mizuguchi, Y. Matsuo and A. Nakahara, Eur. Phys. J. E 40, 88 (2017).
  * (33) K. Miyazaki and D. R. Reichman, Phys. Rev. E 66, 050501(R) (2002).
  * (34) M. Fuchs and M. E. Cates, Phys. Rev. Lett. 89, 248304 (2002).
  * (35) M. Ballauff, J. M. Brader, S. U. Egelhaaf, M. Fuchs, J. Horbach, N. Koumakis, M. Krüger, M. Laurati, K. J. Mutch, G. Petekidis, M. Siebenbürger, Th. Voigtmann, and J. Zausch, Phys. Rev. Lett. 110, 215701 (2013).
  * (36) L. Mohan, R. T. Bonnecaze, and M. Cloitre, Phys. Rev. Lett. 111, 268301 (2013).
  * (37) S. Fritschi, M. Fuchs, and Th. Voigtmann, Soft Matter 10, 4822 (2014).
  * (38) L. Mohan, M. Cloitre, and R. T. Bonnecaze, J. Rheol. 59, 63 (2015).
  * (39) E. Moghimi, A. R. Jacob, and G. Petekidis, Soft Matter 13, 7824 (2017).
  * (40) M. Otsuki and S. Sasa, J. Sta. Mech.: Theory Exp. 2006, L10004 (2006).
  * (41) B. van Aken, P. de Hey, and J. Sietsma, Mater. Sci. Eng., A 278, 247 (2000).
  * (42) F. Varnik, L. Bocquet, and J.-L. Barrat, J. Chem. Phys. 120, 2788 (2004).
  * (43) E. C. Bingham, Fluidity and Plasticity (McGraw-Hill, New York, 1922).