On the DLM/FD methods for simulating neutrally buoyant swimmer motion in non-
Newtonian shear thinning fluids

Ang Lia, Tsorng-Whay Panb,111Corresponding author: e-mail<EMAIL_ADDRESS>and Roland Glowinskib a Department of Mathematics, Lane College, Jackson,
Tennessee 38301, USA b Department of Mathematics, University of Houston,
Houston, Texas 77204, USA

Abstract In this article we discuss the generalization of a Lagrange
multiplier based fictitious domain (DLM/FD) method to simulating the motion of
neutrally buoyant particles of non-symmetric shape in non-Newtonian shear
thinning fluids. Numerical solutions of steady Poiseuille flow of non-
Newtonian shear thinning fluids are compared with the exact solutions in a
two-dimensional channel. Concerning a self-propelled swimmer formed by two
disks, the effect of shear thinning makes the swimmer moving faster and
decreases the critical Reynolds number (for the moving direction changing to
the opposite one) when decreasing the value of the power index in the Carreau-
Bird model.

keywords: Carreau-Bird model; Neutrally buoyant particle; Self-propelled
swimmer; Fictitious domain method; Operator splitting; Finite element
approximations.

##  1 Introduction

Micro-swimmers, like microorganisms, swimming in a low Reynolds number regime,
encounter stringent constraints due to the dominance of viscous over inertial
forces [1, 2]. As a result of kinetic reversibility, Purcell’s scallop theorem
[1] rules out reciprocal motion (i.e., strokes with time-reversal symmetry)
for effective locomotion in the absence of inertia in Newtonian fluids.
Another important fact is that many biological fluids, including blood,
respiratory and cervical mucus, have non-Newtonian properties, such as
viscoelasticity and shear-thinning viscosity (i.e., the viscosity decreases
non-linearly with the shear rate [3] and [4]). The physics governing
microorganism locomotion at low Reynolds numbers in Newtonian fluids is
relatively well understood, while low Reynolds propulsion in non-Newtonian
fluids remains largely unexplored [2]. Since the scallop theorem does not hold
in non-Newtonian fluids (see, e.g., [5] and [6]), it should be possible to
design and build novel micro-swimmers that can move with reciprocal motion in
non-Newtonian fluids. A study on a reciprocal sliding sphere swimmer in a
shear-thinning (non-elastic) fluid [7] suggests that propulsion is achievable
by reciprocal motion, in which backward and forward strokes occur at different
rates. Similarly, another recent study in [8] reported that a symmetric micro-
scallop, a single-hinge microswimmer, can propel in shear thickening and shear
thinning (non-elastic) fluids at low Reynolds number. Thus reciprocal swimming
mechanism can help in designing biomedical microdevices that can propel by a
simple actuation scheme in biological fluids.

In this article, we have focused on the effect of shear thinning on the
swimmer motion in a two-dimensional channel. To study such effect by direct
numerical simulation, we have adapted the Carreau-Bird viscosity model (see
[3] and [4]) for shear thinning property. For example, in [11] and [12],
particle settling in Oldroyd-B viscoelastic fluids with shear thinning
property was studied numerically by using the Carreau-Bird viscosity model.
Direct simulation of the motion of particles has been carried out by using an
arbitrary Lagrangian-Eulerian moving mesh technique with finite element
method. The power law, which is the other commonly used model, was used to
investigate the effects of both shear thickening and thinning on the motion of
scallop swimmers numerically in [8]. A finite element method with the dynamic
mesh-adaptation for moving boundaries was used to study how the scallop
swimmer can swim by reciprocal motion at low (not zero) Reynolds number. For
simulating the swimmer motion in non-Newtonian shear thinning fluids, we have
extended a Lagrange multiplier based fictitious domain (DLM/FD) method
developed for disks in [14] to simulate the motion of a neutrally buoyant non-
symmetric swimmer in non-Newtonian shear thinning fluids. This non-symmetric
swimmer of two disks with different radii is formed by connecting their mass
centers with a massless spring like the one in, e.g., [9] and [10]. To keep
the main flavor of the fictitious domain methods, i.e., the usage of fast
solvers to solve the resulting linear systems, we have splitted the diffusion
term in our proposed methodologies. Accurate computational solutions are
obtained by the proposed methods for Poiseuille flow of non-Newtonian shear
thinning fluids in a two-dimensional channel (thanks to the recently published
exact solutions for such Poiseuille flow by Griffiths in [13]). Then we have
further validated our proposed numerical methods by studying the migration of
56 disks in a two-dimensional channel, which was considered in [29]. For the
effect of shear thinning on the swimmer motion, we have obtained that, via
direct numerical simulation, the swimmer moves faster and the critical
Reynolds number (for the moving direction changing to the opposite one)
decreases when decreasing the value of the power index in the Carreau-Bird
model. The content of this article is as follows: In Section 2 we introduce a
fictitious domain formulation of the model problem associated with the
neutrally buoyant particle cases; then in Section 3 we discuss the time and
space discretizations and the related numerical methodology. We then present
and discuss the numerical results in Section 4. Concluding remarks are given
in Section 5.

##  2 A fictitious domain formulation of the model problem

Figure 1: An example of two-dimensional flow region with one rigid body.

In [14], Pan and Glowinski developed a distributed Lagrange
multiplier/fictitious domain (DLM/FD) method to simulate the motion of
neutrally buoyant disks freely moving in a Newtonian viscous incompressible
fluid in a two dimensional channel. Later they extended in [15], [16], and
[17] such DLM/FD method to simulate and investigate the motion of neutrally
buoyant balls in circular Poiseuille flows. In this article, we have extended
the above DLM/FD method to study the motion of a self-propelled swimmer in
non-Newtonian shear thinning fluids. This non-symmetric swimmer of two disks
with different radii is formed by connecting their mass centers with a
massless spring as the one in, e.g., [9] and [10]. In this section, we will
address the DLM/FD formulation for such problem. Let
$\Omega\subset{\text{I\\!R}}^{2}$ be a rectangular region. We suppose that
$\Omega$ is filled with a non-Newtonian viscous incompressible fluid of
density $\rho_{f}$ and contains a moving neutrally buoyant rigid particle $B$
centered at ${\bf G}=\\{G_{1},G_{2}\\}^{T}$ of density $\rho_{f}$ as shown in
Figure 1; the flow is modeled by the Navier-Stokes equations and the motion of
$B$ is described by the Euler-Newton’s equations. Following the DLM/FD
formulation developed in [14], we define

$\displaystyle W_{0,p}=\\{{\bf v}|{\bf v}\in(H^{1}(\Omega))^{2},\ {\bf v}={\bf
0}\ \text{\it on the top and bottom boundary of $\Omega$ and}$
$\displaystyle\hskip 55.0pt{\bf v}\ \text{\it is periodic in the $x_{1}$
direction}\\},$ $\displaystyle L_{0}^{2}=\\{q|q\in
L^{2}(\Omega),\int_{\Omega}q\,d{\bf x}=0,\\},$
$\displaystyle\Lambda_{0}(t)=\\{{\boldsymbol{\mu}}|{\boldsymbol{\mu}}\in(H^{1}(B(t)))^{2},<{\boldsymbol{\mu}},{\bf
e}_{i}>_{B(t)}=0,\
i=1,2,<{\boldsymbol{\mu}},{{\overrightarrow{Gx}}^{\perp}}>_{B(t)}=0\\}$

with ${\bf e}_{1}=\\{1,0\\}^{T}$, ${\bf e}_{2}=\\{0,1\\}^{T}$,
${{\overrightarrow{Gx}}^{\perp}}=\\{-(x_{2}-G_{2}),x_{1}-G_{1}\\}^{T}$ and
$<\cdot,\cdot>_{B(t)}$ an inner product on $\Lambda_{0}(t)$ which can be the
standard inner product on $(H^{1}(B(t)))^{2}$ (see [23], Section 5, for
further information on the choice of $<\cdot,\cdot>_{B(t)}$). Then the
fictitious domain formulation with distributed Lagrange multipliers for flow
around a freely moving neutrally buoyant particle is as follows

$\displaystyle For\ a.e.\ t>0,\ find\ {\bf u}(t)\in W_{0,p},\ p(t)\in
L_{0}^{2},\ {\bf V}_{{\bf G}}(t)\in{\text{I\\!R}}^{2},\ {\bf
G}(t)\in{\text{I\\!R}}^{2},$ $\displaystyle\omega(t)\in{\text{I\\!R}},\
{\boldsymbol{\lambda}}(t)\in\Lambda_{0}(t)\ \ \text{\it such that}$ (1)
$\displaystyle\begin{cases}\rho_{f}{\displaystyle{\int_{\Omega}}}\left[\dfrac{\partial{\bf
u}}{\partial t}+({\bf u}\cdot{\boldsymbol{\nabla}}){\bf u}\right]\cdot{\bf v}\
d{\bf x}+2{\displaystyle{\int_{\Omega}}}\eta\ {\bf D}({\bf
u})\boldsymbol{:}{\bf D}({\bf v})\ d{\bf
x}-{\displaystyle{\int_{\Omega}}}p{\boldsymbol{\nabla}}\cdot{\bf v}\,d{\bf
x}\\\ \ \ \ -<{\boldsymbol{\lambda}},{\bf
v}>_{B(t)}=\rho_{f}{\displaystyle{\int_{\Omega}}}{\bf g}\cdot{\bf v}\,d{\bf
x}+{\displaystyle{\int_{\Omega}}}{\bf F}\cdot{\bf v}\,d{\bf x},\ \forall{\bf
v}\in W_{0,p},\end{cases}$ (2)
$\displaystyle{\displaystyle{\int_{\Omega}}}q{\boldsymbol{\nabla}}\cdot{\bf
u}(t)d{\bf x}=0,\ \forall q\in L^{2}(\Omega),$ (3)
$\displaystyle<{\boldsymbol{\mu}},{\bf u}(t)>_{B(t)}=<{\boldsymbol{\mu}},{\bf
u}_{p}(t)>_{B(t)},\ \forall{\boldsymbol{\mu}}\in\Lambda_{0}(t),$ (4)
$\displaystyle\dfrac{d{\bf G}}{dt}={\bf V}_{{\bf G}},$ (5) $\displaystyle{\bf
V}_{{\bf G}}(0)={\bf V}_{{\bf G}}^{0},\ \omega(0)=\omega^{0},\ {\bf G}(0)={\bf
G}^{0}=\\{G^{0}_{1},G^{0}_{2}\\}^{T},$ (6) $\displaystyle{\bf u}({\bf
x},0)={\overline{\bf u}}_{0}({\bf x})=\begin{cases}{\bf u}_{0}({\bf x}),\
\forall{\bf x}\in{\Omega\backslash\overline{B(0)}},\\\ {\bf V}_{{\bf
G}}^{0}+\omega^{0}\\{-(x_{2}-G^{0}_{2}),x_{1}-G^{0}_{1}\\}^{T},\ \forall{\bf
x}\in\overline{B(0)},\end{cases}$

where ${\bf u}$ and $p$ denote velocity and pressure, respectively,
$\eta(\dot{\gamma})$ is the fluid viscosity, ${\boldsymbol{\lambda}}$ is a
Lagrange multiplier, ${\bf D}({\bf v})=({\boldsymbol{\nabla}}{\bf
v}+{\boldsymbol{\nabla}}{\bf v}^{T})/2$, ${\bf g}$ is gravity, ${\bf F}$ is
the pressure gradient pointing in the $x_{1}$ direction for particle moving in
a Poiseuille flow, ${\bf V}_{{\bf G}}$ is the translation velocity of the
particle $B$, and $\omega$ is the angular velocity of $B$. We suppose that the
no-slip condition holds on $\partial B$. We also use, if necessary, the
notation $\phi(t)$ for the function ${\bf x}\to\phi({\bf x},t)$. In equation
(3), ${\bf u}_{p}(t)$ is a part of actual particle motion given as follows

${\bf u}({\bf x},t)={\bf V}_{{\bf
G}}(t)+\omega(t)\\{-(x_{2}-G_{2}),x_{1}-G_{1}\\}^{T}+{\bf u}_{p}({\bf x},t).$

Thus we have ${\bf u}_{p}={\bf 0}$ if $B$ is a single piece neutrally buoyant
rigid particle like a disk suspended in fluid; but for $B$ as a swimmer formed
by two neutrally buoyant disks considered in this article, ${\bf u}_{p}(t)$ is
a given reciprocal motion of the two disks with respect to the swimmer mass
center (see Section 4). Shear thinning can be easily added to the above model
by using the Carreau-Bird viscosity model (see, [3] and [4])

(7)
$\dfrac{\eta-\eta_{\infty}}{\eta_{0}-\eta_{\infty}}=(1+(\lambda_{3}\dot{\gamma})^{2})^{(n-1)/2},$

where $\eta_{0}$ (resp., $\eta_{\infty}$) is the fluid viscosity at zero shear
rate (resp., infinite shear rate), $\lambda_{3}$ is the relaxation time,
$\dot{\gamma}$ is determined by the second invariant of the rate of strain
tensor $(\dot{\gamma})^{2}=({\bf D}({\bf u}):{\bf D}({\bf u}))/2$, and the
power index $n$ is in $(0,1]$ (if $n=1$, the fluid is a Newtonian fluid with
constant viscosity $\eta_{0}$).

######  Remark 1.

In (3), the rigid body motion in the region occupied by the particle is
enforced via Lagrange multipliers ${\boldsymbol{\lambda}}$. To recover the
translation velocity ${\bf V}_{{\bf G}}(t)$ and the angular velocity
$\omega(t)$, we solve the following equations

(8) $\begin{cases}<{\bf e}_{i},{\bf u}(t)-{\bf V}_{{\bf G}}(t)-\omega(t)\
{{\overrightarrow{Gx}}^{\perp}}>_{B(t)}=0,\ for\ i=1,2,\\\
<{{\overrightarrow{Gx}}^{\perp}},{\bf u}(t)-{\bf V}_{{\bf G}}(t)-\omega(t)\
{{\overrightarrow{Gx}}^{\perp}}>_{B(t)}=0.\end{cases}$

##  3 Space approximation and time discretization

Concerning the finite element based space approximation of problem (1)-(6), we
use $P_{1}$-$iso$-$P_{2}$ and $P_{1}$ finite elements for the velocity field
and pressure, respectively (see, e.g., [18] and [19](Chapter 5)). More
precisely with $h$ a space discretization step, we introduce a finite element
triangulation ${\mathcal{T}_{h}}$ of ${\overline{\Omega}}$ and then
${\mathcal{T}}_{2h}$ a triangulation twice coarser. Then we approximate
$W_{0,p}$, $L^{2}$ and $L^{2}_{0}$ by the following finite dimensional spaces

(9) $\displaystyle W_{0,h}=\\{{\bf v}_{h}|{\bf
v}_{h}\in(C^{0}(\overline{\Omega}))^{2},\ {\bf v}_{h}|_{T}\in P_{1}\times
P_{1},\ \forall T\in{\mathcal{T}}_{h},\ {\bf v}_{h}={\bf 0}\ \text{\it on
the}$ $\displaystyle\hskip 60.0pt\text{\it top and bottom boundary of $\Omega$
and}\ {\bf v}\ \text{\it is periodic in the $x_{1}$ direction}\ \\},$
$\displaystyle L^{2}_{h}=\\{q_{h}$ $\displaystyle|$ $\displaystyle q_{h}\in
C^{0}(\overline{\Omega}),\ q_{h}|_{T}\in P_{1},\ \forall
T\in{\mathcal{T}}_{2h},\ q_{h}\ \text{\it is periodic}$
$\displaystyle\text{\it in the $x_{1}$ direction}\\},$

and

(11) $L^{2}_{0,h}=\\{q_{h}|q_{h}\in L^{2}_{h},\ \int_{\Omega}q_{h}\,d{\bf
x}=0\\},$

respectively; in (9)-(11), $P_{1}$ is the space of polynomials in two
variables of degree $\leq 1$.

Figure 2: An example of set of collocation points chosen inside the disk and
at its boundary for enforcing the rigid body motion.

A finite dimensional space approximating $\Lambda_{0}(t)$ is defined in the
following: let $\\{{\bf x}_{i}\\}_{i=1}^{N}$ be a set of points covering
$\overline{B(t)}$; we define then

(12)
$\Lambda_{h}(t)=\\{{\boldsymbol{\mu}}_{h}|{\boldsymbol{\mu}}_{h}={\displaystyle\sum}_{i=1}^{N}{\boldsymbol{\mu}}_{i}\delta({\bf
x}-{\bf x}_{i}),\ {\boldsymbol{\mu}}_{i}\in{\text{I\\!R}}^{2},\ \forall
i=1,...,N\\},$

where $\delta(\cdot)$ is the Dirac measure at ${\bf x}={\bf 0}.$ Then, instead
of the inner product of $(H^{1}(B(t)))^{2}$ we shall use
$<\cdot,\cdot>_{B_{h}(t)}$ defined by

(13) $<{\boldsymbol{\mu}}_{h},{\bf
v}_{h}>_{B_{h}(t)}={\displaystyle\sum}_{i=1}^{N}{\boldsymbol{\mu}}_{i}\cdot{\bf
v}_{h}({\bf x}_{i}),\ \forall{\boldsymbol{\mu}}_{h}\in\Lambda_{h}(t),\ {\bf
v}_{h}\in W_{0,h}.$

Then we approximate $\Lambda_{0}(t)$ by

(14) $\Lambda_{0,h}(t)=\\{\mu_{h}|\mu_{h}\in\Lambda_{h}(t),\ <\mu_{h},{\bf
e}_{i}>_{B_{h}(t)}=0,\ i=1,2,\
<\mu_{h},{{\overrightarrow{Gx}}^{\perp}}>_{B_{h}(t)}=0\\}.$

A typical choice of points for defining (12) is to take the grid points of the
velocity mesh internal to the particle $B$ and whose distance to the boundary
of $B$ is greater than, e.g. $h$, and to complete with selected points from
the boundary of $B(t)$ (see Figure 2).

Using the above finite dimensional spaces leads to the following approximation
of problem (1)-(6):

$\displaystyle For\ a.e.\ t>0,\ find\ {\bf u}(t)\in W_{0,h},\ p(t)\in
L_{0,h}^{2},\ {\bf V}_{{\bf G}}(t)\in{\text{I\\!R}}^{2},\ {\bf
G}(t)\in{\text{I\\!R}}^{2},$ $\displaystyle\omega(t)\in{\text{I\\!R}},\
{\boldsymbol{\lambda}}_{h}(t)\in\Lambda_{0,h}(t)\ \ \text{\it such that}$ (15)
$\displaystyle\begin{cases}\rho_{f}{\displaystyle{\int_{\Omega}}}\left[\dfrac{\partial{\bf
u}_{h}}{\partial t}+({\bf u}_{h}\cdot{\boldsymbol{\nabla}}){\bf
u}_{h}\right]\cdot{\bf v}\,d{\bf x}+2{\displaystyle{\int_{\Omega}}}\eta\ {\bf
D}({\bf u}_{h})\boldsymbol{:}{\bf D}({\bf v})\ d{\bf x}\\\ \hskip
20.0pt-{\displaystyle{\int_{\Omega}}}p_{h}{\boldsymbol{\nabla}}\cdot{\bf
v}\,d{\bf x}-<{\boldsymbol{\lambda}}_{h},{\bf
v}>_{B_{h}(t)}={\displaystyle{\int_{\Omega}}}{\bf F}\cdot{\bf v}\,d{\bf x},\
\forall{\bf v}\in W_{0,h},\end{cases}$ (16)
$\displaystyle{\displaystyle{\int_{\Omega}}}q{\boldsymbol{\nabla}}\cdot{\bf
u}_{h}(t)d{\bf x}=0,\ \forall q\in L^{2}_{h},$ (17)
$\displaystyle<{\boldsymbol{\mu}},{\bf
u}_{h}(t)>_{B_{h}(t)}=<{\boldsymbol{\mu}},{\bf u}_{p}(t)>_{B_{h}(t)},\
\forall{\boldsymbol{\mu}}\in\Lambda_{0,h}(t),$ (18) $\displaystyle\dfrac{d{\bf
G}}{dt}={\bf V}_{{\bf G}},$ (19) $\displaystyle{\bf V}_{{\bf G}}(0)={\bf
V}_{{\bf G}}^{0},\ \omega(0)=\omega^{0},\ {\bf G}(0)={\bf
G}^{0}=\\{G^{0}_{1},G^{0}_{2}\\}^{T},$ (20) $\displaystyle{\bf u}_{h}({\bf
x},0)={\overline{\bf u}}_{0,h}({\bf x})\
(\text{with}{\boldsymbol{\nabla}}\cdot{\overline{\bf u}}_{0,h}=0).$

Applying a first order operator splitting scheme, namely the Lie scheme [20]
(see also [19] and [21] (Chapter 2)), to discretize equations (15)-(20) in
time, we obtain (after dropping some of the subscripts $h$):

(21) ${\bf u}^{0}={\overline{\bf u}}_{0,h},\ {\bf V}_{{\bf G}}^{0},\
\omega^{0},\ and\ {\bf G}^{0}\ are\ given;$

For $n\geq 0$, knowing ${\bf u}^{n}$, ${\bf V}_{{\bf G}}^{n}$, $\omega^{n}$
and ${\bf G}^{n}$, compute ${\bf u}^{n+1/6}$ and $p^{n+1/6}$ via the solution
of

(22)
$\begin{cases}\displaystyle\rho_{f}{\displaystyle{\int_{\Omega}}}\frac{{\bf
u}^{n+1/6}-{\bf u}^{n}}{\triangle t}\cdot{\bf v}\,d{\bf
x}-{\displaystyle{\int_{\Omega}}}p^{n+1/6}{\boldsymbol{\nabla}}\cdot{\bf
v}\,d{\bf x}=0,\ \forall{\bf v}\in W_{0,h},\\\
\displaystyle{\displaystyle{\int_{\Omega}}}q{\boldsymbol{\nabla}}\cdot{\bf
u}^{n+1/6}\,d{\bf x}=0,\ \forall q\in L^{2}_{h};\ {\bf u}^{n+1/6}\in W_{0,h},\
p^{n+1/6}\in L^{2}_{0,h}.\end{cases}$

Then compute ${\bf u}^{n+2/6}$ via the solution of

(23)
$\displaystyle\begin{cases}\displaystyle{\displaystyle{\int_{\Omega}}}\frac{\partial{\bf
u}}{\partial t}\cdot{\bf v}\,d{\bf x}+{\displaystyle{\int_{\Omega}}}({\bf
u}^{n+1/6}\cdot{\boldsymbol{\nabla}}){\bf u}\cdot{\bf v}\,d{\bf x}=0,\
\forall{\bf v}\in W_{0,h},\ on\ (t^{n},t^{n+1}),\\\ \displaystyle\ {\bf
u}(t^{n})={\bf u}^{n+1/6};\ \ {\bf u}(t)\in W_{0,h},\end{cases}$ (24)
$\displaystyle{\bf u}^{n+2/6}={\bf u}(t^{n+1}).$

Next, compute ${\bf u}^{n+3/6}$ via the solution of

(25)
$\begin{cases}\displaystyle\rho_{f}{\displaystyle{\int_{\Omega}}}\frac{{\bf
u}^{n+3/6}-{\bf u}^{n+2/6}}{\triangle t}\cdot{\bf v}\,d{\bf
x}+2\alpha{\displaystyle{\int_{\Omega}}}\eta(\dot{{\bf u}^{n+2/6}})\ {\bf
D}({\bf u}^{n+3/6})\boldsymbol{:}{\bf D}({\bf v})\ d{\bf x}\\\ \hskip
20.0pt={\displaystyle{\int_{\Omega}}}{\bf F}\cdot{\bf v}\,d{\bf x},\forall{\bf
v}\in W_{0,h};\ \ {\bf u}^{n+3/6}\in W_{0,h}.\end{cases}$

Now predict the position and the translation velocity of the center of mass of
the particles as follows: Take ${\bf V}_{{\bf G}}^{n+\frac{4}{6},0}={\bf
V}_{{\bf G}}^{n}$ and ${\bf G}^{n+\frac{4}{6},0}={\bf G}^{n}$; then predict
the new position of the particle via the following sub-cycling and predicting-
correcting technique:

For $k=1,\dots,N$, compute (26) $\displaystyle{\widehat{\bf V}_{{\bf
G}}}^{n+\frac{4}{6},k}={\bf V}_{{\bf G}}^{n+\frac{4}{6},k-1}+{\bf F}^{r}({\bf
G}^{n+\frac{4}{6},k-1}){\triangle t}/{2N},$ (27) $\displaystyle{\widehat{\bf
G}}^{n+\frac{4}{6},k}={\bf G}^{n+\frac{4}{6},k-1}+({\widehat{\bf V}_{{\bf
G}}}^{n+\frac{4}{6},k}+{\bf V}_{{\bf G}}^{n+\frac{4}{6},k-1}){\triangle
t}/{4N},$ (28) $\displaystyle{\bf V}_{{\bf G}}^{n+\frac{4}{6},k}={\bf V}_{{\bf
G}}^{n+\frac{4}{6},k-1}+({\bf F}^{r}({\widehat{\bf G}}^{n+\frac{4}{6},k})+{\bf
F}^{r}({\bf G}^{n+\frac{4}{6},k-1})){\triangle t}/{4N},$ (29)
$\displaystyle{\bf G}^{n+\frac{4}{6},k}={\bf G}^{n+\frac{4}{6},k-1}+({\bf
V}_{{\bf G}}^{n+\frac{4}{6},k}+{\bf V}_{{\bf
G}}^{n+\frac{4}{6},k-1}){\triangle t}/{4N},$ enddo; (30) and let ${\bf
V}_{{\bf G}}^{n+\frac{4}{6}}={\bf V}_{{\bf G}}^{n+\frac{4}{6},N},\ {\bf
G}^{n+\frac{4}{6}}={\bf G}^{n+\frac{4}{6},N}$.

Now, compute ${\bf u}^{n+5/6}$, ${{\boldsymbol{\lambda}}}^{n+5/6}$, ${\bf
V}_{{\bf G}}^{n+5/6}$, and $\omega^{n+5/6}$ via the solution of

(31)
$\begin{cases}\displaystyle\rho_{f}{\displaystyle{\int_{\Omega}}}\frac{{\bf
u}^{n+5/6}-{\bf u}^{n+3/6}}{\triangle t}\cdot{\bf v}\,d{\bf
x}+2\beta{\displaystyle{\int_{\Omega}}}\eta(\dot{{\bf u}^{n+3/6}})\ {\bf
D}({\bf u}^{n+5/6})\boldsymbol{:}{\bf D}({\bf v})\ d{\bf x}\\\ \
=<{\boldsymbol{\lambda}},{\bf v}>_{B_{h}^{n+4/6}},\ \forall{\bf v}\in
W_{0,h},\\\ <{\boldsymbol{\mu}},{\bf
u}^{n+5/6}>_{B_{h}^{n+4/6}}=<{\boldsymbol{\mu}},{\bf
u}_{p}(t^{n+1})>_{B_{h}^{n+4/6}},\
\forall{{\boldsymbol{\mu}}}\in\Lambda_{0,h}^{n+4/6};\displaystyle\ {\bf
u}^{n+5/6}\in
W_{0,h},{{\boldsymbol{\lambda}}}^{n+5/6}\in\Lambda_{0,h}^{n+4/6},\end{cases}$

and solve for ${\bf V}_{{\bf G}}^{n+5/6}$ and $\omega^{n+5/6}$ from

(32) $\begin{cases}<{\bf e}_{i},{\bf u}^{n+5/6}-{\bf V}_{{\bf
G}}^{n+5/6}-\omega^{n+5/6}\
{\overrightarrow{G^{n+4/6}x}}^{\perp}>_{B_{h}^{n+4/6}}=0,\ for\ i=1,2,\\\
<{\overrightarrow{G^{n+4/6}x}}^{\perp},{\bf u}^{n+5/6}-{\bf V}_{{\bf
G}}^{n+5/6}-\omega^{n+5/6}\
{\overrightarrow{G^{n+4/6}x}}^{\perp}>_{B_{h}^{n+4/6}}=0,\end{cases}$

Finally, take ${\bf V}_{{\bf G}}^{n+1,0}={\bf V}_{{\bf G}}^{n+5/6}$ and ${\bf
G}^{n+1,0}={\bf G}^{n+4/6}$; then predict the final position and translation
velocity as follows:

For $k=1,\dots,N$, compute (33) $\displaystyle{\widehat{\bf V}_{{\bf
G}}}^{n+1,k}={\bf V}_{{\bf G}}^{n+1,k-1}+{\bf F}^{r}({\bf
G}^{n+1,k-1}){\triangle t}/{2N},$ (34) $\displaystyle{\widehat{\bf
G}}^{n+1,k}={\bf G}^{n+1,k-1}+({\widehat{\bf V}_{{\bf G}}}^{n+1,k}+{\bf
V}_{{\bf G}}^{n+1,k-1}){\triangle t}/{4N},$ (35) $\displaystyle{\bf V}_{{\bf
G}}^{n+1,k}={\bf V}_{{\bf G}}^{n+1,k-1}+({\bf F}^{r}({\widehat{\bf
G}}^{n+1,k})+{\bf F}^{r}({\bf G}^{n+1,k-1})){\triangle t}/{4N},$ (36)
$\displaystyle{\bf G}^{n+1,k}={\bf G}^{n+1,k-1}+({\bf V}_{{\bf
G}}^{n+1,k}+{\bf V}_{{\bf G}}^{n+1,k-1}){\triangle t}/{4N},$ enddo;

and let ${\bf V}_{{\bf G}}^{n+1}={\bf V}_{{\bf G}}^{n+1,N},\ {\bf
G}^{n+1}={\bf G}^{n+1,N}$; and set ${\bf u}^{n+1}={\bf u}^{n+5/6}$,
$\omega^{n+1}=\omega^{n+5/6}$.

In algorithm (21)-(36), we have $t^{n+s}=(n+s)\triangle t$,
$\Lambda_{0,h}^{n+s}=\Lambda_{0,h}(t^{n+s})$, $B_{h}^{n+s}$ is the region
occupied by the particle centered at ${\bf G}^{n+s}$, and ${\bf F}^{r}$ is a
short range repulsion force which prevents the particle/particle and
particle/wall penetration (see, e.g., [22, 23]). Finally, $\alpha$ and $\beta$
verify $\alpha+\beta=1$; we have chosen $\alpha=1$ and $\beta=0$ in the
numerical simulations discussed later.

###  3.1 Solutions of the subproblems (22), (23), (25) and (31)

The degenerated quasi-Stokes problem (22) is solved by a preconditioned
conjugate gradient method introduced in [24], in which the discrete elliptic
problems from the preconditioning are solved by a matrix-free fast solver from
FISHPAK by Adams et al. in [25]. The advection problem (23) for the velocity
field is solved by a wave-like equation method as in [26, 27]. To enforce the
rigid body motion inside the region occupied by the particles, we employ the
conjugate gradient method discussed in [14] for problem (31). Problem (25) is
a classical discrete elliptic problem. To solve problem (25) by the same
matrix-free fast solver, we have modified the equation as follows:

(37)
$\begin{cases}\displaystyle\rho_{f}{\displaystyle{\int_{\Omega}}}\frac{{\bf
u}^{n+3/6}-{\bf u}^{n+2/6}}{\triangle t}\cdot{\bf v}\,d{\bf
x}+\alpha\eta_{0}{\displaystyle{\int_{\Omega}}}{\boldsymbol{\nabla}}{\bf
u}^{n+3/6}{\bf:}{\boldsymbol{\nabla}}{\bf v}\,d{\bf
x}={\displaystyle{\int_{\Omega}}}{\bf F}\cdot{\bf v}\,d{\bf x}\\\ \hskip
20.0pt+\alpha\eta_{0}{\displaystyle{\int_{\Omega}}}{\boldsymbol{\nabla}}{\bf
u}^{n+2/6}{\bf:}{\boldsymbol{\nabla}}{\bf v}\,d{\bf
x}-2\alpha{\displaystyle{\int_{\Omega}}}\eta(\dot{{\bf u}^{n+2/6}})\ {\bf
D}({\bf u}^{n+2/6})\boldsymbol{:}{\bf D}({\bf v})\ d{\bf x},\\\ \forall{\bf
v}\in W_{0,h};\ \ {\bf u}^{n+3/6}\in W_{0,h}.\end{cases}$

In problem (37), the value of ${\boldsymbol{\nabla}}{\bf u}^{n+2/6}$ at each
mesh node is obtained via a super-convergent scheme discussed in [28]. In this
article, all computational results of shear thinning fluids were obtained by
algorithm (21)-(36) with the above subproblem (37) instead of problem (25).

##  4 Numerical experiments and discussion

###  4.1 Poiseuille flow of non-Newtonian shear thinning fluids

In the first validation case, we have considered the pressure driven
Poiseuille flow of a non-Newtonian shear thinning fluid. The computational
results are validated by comparing them with the exact solutions recently
reported in [13]. Since there is no particle involved in the computation, the
computational results are obtained from steps (21), (22), (23), (24), and (37)
in algorithm (21)-(36).

The computational domain is $\Omega=(0,L)\times(-H,H)=(0,8)$ $\times$
$(-1,1)$. We have used structured triangular meshes in all computations. The
mesh sizes for the velocity field are $h_{v}=1/32$, 1/64, and 1/96, while the
mesh size for the pressure is $h_{p}=2h_{v}$. The time step is chosen to be
$\triangle t=0.001$. The pressure gradient pointing in the $x_{1}$ direction
is $(F_{1},F_{2})^{T}=(-G,0)^{T}=(-1,0)^{T}$ so that the flow moves from left
to right. The flow velocity is periodic in the $x_{1}$ (horizontal) direction
and zero at the top and bottom boundary of the domain. The initial flow
velocity is ${\bf 0}$. The fluid density is 1 and the fluid viscosities are
$\eta_{0}=1$ and $\eta_{\infty}=0$. The three values of the power index $n$
considered in this case are 1/3, 1/2, and 2/3. Then as given in [13], for the
power index $n=1/3$ in (7), the exact solution $(u_{1},0)$ and the associated
relaxation time $\lambda_{3}$ of the above shear thinning Poiseuille flow
problem are given by

$\displaystyle\lambda_{3}=\sqrt[6]{2}\sqrt{3}\eta_{0}/GH=\sqrt[6]{2}\sqrt{3},$
$\displaystyle
u_{1}(Y)=\dfrac{\sqrt[3]{2}GH^{2}}{4\eta_{0}}\left[\dfrac{3}{\sqrt[3]{4}}+1-Y^{4}\right.$
$\displaystyle\hskip
60.0pt\left.+\dfrac{(3-\sqrt{1+8Y^{6}})(1+\sqrt{1+8Y^{6}})^{2/3}-2Y^{2}(3+\sqrt{1+8Y^{6}})}{4(1+\sqrt{1+8Y^{6}})^{1/3}}\right],$

where $Y=y/H$ and $-H\leq y\leq H$. For $n=1/2$, they are

$\displaystyle\lambda_{3}=\sqrt{2}\eta_{0}/(\sqrt[4]{3}GH)=\sqrt{2}/\sqrt[4]{3},$
$\displaystyle
u_{1}(Y)=\dfrac{GH^{2}(\sqrt{3+Y^{4}}-2Y^{2})(Y^{2}+\sqrt{3+Y^{4}})^{1/2}}{3\sqrt[4]{3}\eta_{0}},$

and for $n=2/3$, they are

$\displaystyle\lambda_{3}=\sqrt{3}\eta_{0}/(\sqrt[3]{2}GH)=\sqrt{3}/\sqrt[3]{2},$ $\displaystyle u_{1}(Y)=\dfrac{GH^{2}}{5\sqrt[6]{2}\eta_{0}}\left[\sqrt{2}\right.$ $\displaystyle\hskip 10.0pt\left.+\dfrac{[(1+\sqrt{1-Y^{6}})^{4/3}+Y^{4}-3Y^{2}(1+\sqrt{1-Y^{6}})^{2/3}][(1+\sqrt{1-Y^{6}})^{2/3}+Y^{2}]^{1/2}}{(1+\sqrt{1-Y^{6}})^{5/6}}\right],$ $n$ |  $h$ |  $L^{2}$-error | order |  $n$ |  $h$ |  $L^{2}$-error | order
---|---|---|---|---|---|---|---
1/3 | 1/32 | 8.510003368$\times 10^{-3}$ |  | 1/2 | 1/32 | 1.405723487$\times 10^{-4}$ |
1/3 | 1/64 | 2.475437454$\times 10^{-3}$ | 1.833 | 1/2 | 1/64 | 3.667227424$\times 10^{-5}$ | 1.938
1/3 | 1/96 | 1.157654939$\times 10^{-3}$ | 1.815 | 1/2 | 1/96 | 1.633156504$\times 10^{-5}$ | 1.959
2/3 | 1/32 | 7.847740427$\times 10^{-5}$ |  |  |  |  |
2/3 | 1/64 | 1.990780506$\times 10^{-5}$ | 2.008 |  |  |  |
2/3 | 1/96 | 8.704140007$\times 10^{-6}$ | 2.001 |  |  |  |

Table 1. $L^{2}$-errors for the steady state velocity field of Poiseuille flow
for the power indices $n=1/3$, 1/2, and 2/3.

Figure 3: Normalized horizontal velocity profiles obtained for $h=1/96$ (red
dashed line) and exact solution (blue solid line) for $n=1/3$. Figure 4:
Normalized horizontal velocity profiles obtained for $h=1/96$ (red dashed
line) and exact solution (blue solid line) for $n=1/2$. Figure 5: Normalized
horizontal velocity profiles obtained for $h=1/96$ (red dashed line) and exact
solution (blue solid line) for $n=2/3$.

The numerical results obtained with $h=1/96$ and the actual profiles are shown
in Figures 3, 4, and 5 for the steady state horizontal profile of Poiseuille
flow. The $L^{2}$-errors of steady numerical solutions in Table 1 show that
the expected error order for a $P_{1}$ finite element approximation has been
obtained.

(a)

(b) 1

Figure 6: Migration of 56 neutrally buoyant identical particles in the
pressure driven flow of a Newtonian fluid without shear thinning ($n=1.0$).
(a) Velocity profile of the fluid without particles (top curve) and velocities
of particles (bottom curve); (b) particle positions in the channel at time
$t=100$.

Figure 7: Velocity profile of the fluid without particles and velocities of
the 56 particles: $n=0.7$ at $t=42.46$ (top), $n=0.5$ at $t=31.78$ (middle),
and $n=0.4$ at $t=27.41$ (bottom).

Figure 8: Positions of the 56 neutrally buoyant particles in the pressure
driven flow of a shear thinning fluid: $n=0.7$ at $t=42.46$ (top), $n=0.5$ at
$t=31.78$ (middle), and $n=0.4$ at $t=27.41$ (bottom).

###  4.2 On the migration of 56 disks in the Poiseuille flow of a shear
thinning fluid

The second validation case we consider concerns the motion and cross stream
migration of 56 neutrally buoyant identical circular cylinders in the pressure
driven Poiseuille flow of a on-Newtonian shear thinning fluid, which was a
test problem studied numerically in [29]. The proposed numerical scheme has
been validated by comparing our results with those reported in [29]. The
computational domain is $\Omega=(0,21)$ $\times$ $(0,10)$. The initial flow
velocity and particle velocities are ${\bf 0}$. The pressure drop is
$F_{1}=-2$ so that the maximum horizontal speed is 25 when there is no
particle for the Poiseuille flow of a Newtonian fluid with the viscosity
$\eta_{0}=1$. For non-Newtonian shear thinning fluids modeled by (7), the
viscosity at infinity shear rate is $\eta_{\infty}=0.1\eta_{0}$ and the
relaxation time is $\lambda_{3}=1$. The values of the power index $n$ are 0.4,
0.5, and 0.7. The fluid and particle densities are 1. The disk diameter $d$ is
1 and ${\bf u}_{p}(t)$ is ${\bf 0}$ for each disk. The mesh size for the
velocity field is $h_{v}=1/16$, while the mesh size for the pressure is
$h_{p}=2h_{v}$. The time step is $\triangle t=0.001$.

For the case of a Newtonian fluid without shear thinning ($n=1$), the
migration of the disks is shown in Figure 6. There are no disks at the center
line and the disks tend to accumulate a distance of 0.6 from the center line
of the channel, which is known as the Segre-Silberberg effect (see [30] and
[31]). In general, lubrication forces move the particles away from the channel
walls while the curvature of the velocity profile moves the particles away
from the centerline of the channel. For a small neutrally buoyant ball in a
plane Poiseuille flow via the matched asymptotic expansion methods, it was
found that the of the fluid without particles and the velocities of the
particles at $t=100$ are shown in Figure 6. At $t=100$, the maximum particle
velocity is 14.0864, while the maximum fluid velocity without particles is 25
so the slip velocity is 10.9136 near the centerline. In [29] the maximum
particle velocity is 15 near the centerline and the maximum slip velocity is
10. Both results are in good agreement.

The velocity profile with shear thinning is quite different from the one at
$n=1$ as shown in Figure 7. The effects of shear thinning can be increased by
decreasing the power index $n$. The maximum particle velocities are 23.0664,
29.9960, and 30.1694, while the maximum fluid velocities are 47.3956, 84.5356,
and 111.9036 for $n=0.7$, 0.5, and 0.4, respectively. Hence we obtain the slip
velocities in Table 2 by following the approach given in [29]. Their slip
velocities are 23, 54, and 81.4 for $n=0.7$, 0.5, and 0.4, respectively,
implying our numerical results are in a good agreement with those in [29]. Due
to the Segre-Silberberg effect, the disks separate into two groups (see
Figures 7 and 8), the top group and the bottom one according to the particle
positions shown in Figure 8. Following these particle separations, there are
three particle free zones, the one in the middle (called middle zone), the one
next to the top wall (called the top boundary zone), and the one next to the
bottom wall (the bottom boundary zone), which can be identified easily from
the particle velocity plots in Figure 7. The sizes of these particle free
zones are called the gaps and shown in Table 2, in which the average boundary
gap is the mean of the two boundary gaps. Due to the shear thinning effect,
the disk accumulation is enhanced so that the middle gap of particle free zone
increases and the average boundary gap decreases, when decreasing the value of
the power index $n$. Our results are consistent with those obtained in [29]
qualitatively.

$n$ | Slip velocity | Average boundary gap | Middle gap
---|---|---|---
1 | 10.9136 | 0.95167 | 1.60403
0.7 | 24.3292 | 0.85906 | 1.67752
0.5 | 54.5396 | 0.83388 | 1.80585
0.4 | 81.7342 | 0.83313 | 2.00206

Table 2. The slip velocity, average of gap sizes next to two walls, and gap
size in the center of the channel for the power indices $n=1/3$, 1/2, 2/3, and
1.

###  4.3 On the motion of a self-propelled swimmer

Even Purcell’s scallop theorem [1] rules out reciprocal motion (i.e., strokes
with time-reversal symmetry) for effective locomotion in the absence of
inertia ($Re=0$) in Newtonian fluids, the motion of a simple reciprocal model
swimmer (an asymmetric dumbbell) in a Newtonian fluid at intermediate Reynolds
numbers has been studied extensively in [9] and [10]. In this section, we have
considered the case of a neutrally buoyant self-propelled swimmer freely
moving in a non-Newtonian shear thinning fluid. Such swimmer is formed by two
disks of different sizes connecting by a massless spring at their mass centers
(see Figure 9). Hence it is an asymmetric long body, but symmetric with
respect to the line segment connecting two disk mass centers. The reciprocal
motion of two disks with respect to the swimmer mass center is prescribed (see
Figure 9 and equation (38)). Due to the asymmetric disk motion with respect to
the swimmer mass center, swimmer can move in either directions parallel to the
string connecting the two disk mass centers. It is interesting to find out how
the shear thinning impacts the swimmer motion. Following the units used in
[9], we assume in this section all dimensional quantities are in the physical
MKS units.

Figure 9: An example of the relative position of the two disks: The reciprocal
motion with respect to the swimmer mass center at different time over one
period is $d(t)=d_{0}+A\sin(t)$. Disk centers of mass (blue *) are indicated.

Figure 10: Average horizontal velocity $V_{1}$ of a swimmer versus $Re$ for
the power index $n=1$ (top) and 0.8 (bottom).

Figure 11: Velocity field next to the swimmer of two same size disks and
swimmer position at $t=8$ (top) and 8.05 (bottom) for $Re=15$ and $n=1$.
Figure 12: Histories of the horizontal $x_{1}$-position of the swimmer mass
center for $Re=5$ and 15 and $n=1$ where the swimmer is formed by two
different size disks (see Figure 13).

The computational domain is $\Omega=(-4,4)$ $\times$ $(-4,4)$. The flow
velocity is ${\bf 0}$ initially and zero at the boundary of domain $\Omega$.
The densities of fluid and swimmer are $10^{3}$, which is the water density.
The initial position of the swimmer mass center is at (0,0). The distance
between the two disk centers is $d(t)=d_{0}+A\sin(\omega t)$ as in [9] where
the equilibrium distance between the two disk centers is $d_{0}=10r_{0}$, the
amplitude is $A=r_{0}$, and $r_{0}=0.15$. The angular frequency is
$\omega=20\pi$ (so the frequency is 10). The massless spring connecting the
two disks is located in the direction parallel to the $x_{1}$ direction. The
radius of the larger disk $D_{M}$ is $R_{M}=2r_{0}=0.3$ and that of the
smaller disk $D_{m}$ is $R_{m}=r_{0}=0.15$. The amplitude $A$ is
$A=A_{M}+A_{m}$ where the component for the larger disk is
$A_{M}=\dfrac{AR_{m}^{2}}{R_{M}^{2}+R_{m}^{2}}$ and the one for the smaller
disk is $A_{m}=\dfrac{AR_{M}^{2}}{R_{M}^{2}+R_{m}^{2}}$. The Reynolds number
is defined by $Re=A_{m}R_{m}\omega/\nu$ (as in [9]) where $\nu$ is the kinetic
viscosity of the fluid at zero shear rate. We vary the value of $Re$ to obtain
the fluid viscosity at zero shear rate. We have used structured triangular
meshes in all simulations. The mesh size for the velocity field is
$h_{v}=1/128$, while the mesh size for pressure is $h_{p}=2h_{v}$. The time
step is $\triangle t=0.001$. Concerning the shear thinning, the ratio of
viscosity at infinity shear rate and the one at zero shear rate is 0.1 and the
relaxation time is $\lambda_{3}=1$. The values of the power index $n$ are 0.5,
0.6, 0.7, 0.8, and 0.9.

Figure 13: Swimmer position at $t=0$ (top), velocity field next to the swimmer
and swimmer position at $t=20$ for $Re=5$ (middle) and 15 (bottom) and $n=1$.

In this section, we have mainly studied the effect of the power index $n$ and
of the Reynolds number on the moving direction and speed of the swimmer. In
order to focus on the swimmer moving direction, the swimmer is only allowed to
freely move in the $x_{1}$-direction. To accommodate this restricted motion
and the Dirichlet boundary condition, we have to modify the spaces defined
previously as follows

$\displaystyle W_{0}$ $\displaystyle=$ $\displaystyle\\{{\bf v}|{\bf
v}\in(H^{1}(\Omega))^{2},\ {\bf v}={\bf 0}\ \text{\it on the boundary of
$\Omega$.}\\},$ $\displaystyle\Lambda_{0}(t)$ $\displaystyle=$
$\displaystyle\\{{\boldsymbol{\mu}}|{\boldsymbol{\mu}}\in(H^{1}(B(t)))^{2},<{\boldsymbol{\mu}},{\bf
e}_{1}>_{B(t)}=0\\},$ $\displaystyle W_{0,h}$ $\displaystyle=$
$\displaystyle\\{{\bf v}_{h}|{\bf v}_{h}\in(C^{0}(\overline{\Omega}))^{2},\
{\bf v}_{h}|_{T}\in P_{1}\times P_{1},\forall T\in{\mathcal{T}}_{h},$
$\displaystyle{\bf v}_{h}={\bf 0}\ \text{\it on the boundary}\ \text{\it of
$\Omega$}\\},$ $\displaystyle L^{2}_{0,h}$ $\displaystyle=$
$\displaystyle\\{q_{h}|q_{h}\in C^{0}(\overline{\Omega}),\ q_{h}|_{T}\in
P_{1},\ \forall T\in{\mathcal{T}}_{2h},\ \int_{\Omega}q_{h}\,d{\bf x}=0\\},$
$\displaystyle\Lambda_{0,h}(t)$ $\displaystyle=$
$\displaystyle\\{\mu_{h}|\mu_{h}\in\Lambda_{h}(t),\ <\mu_{h},{\bf
e}_{1}>_{B_{h}(t)}=0\\}.$

Following the above given restricted motion, the reciprocal motion with
respect to the swimmer mass center is

(38) ${\bf u}_{p}(t)=\begin{cases}(A_{M}\omega\cos(\omega t),0)^{T},\
\forall{\bf x}\in D_{M},\\\ (-A_{m}\omega\cos(\omega t),0)^{T},\ \forall{\bf
x}\in D_{m}.\\\ \end{cases}$

Figure 14: Average horizontal speed $V_{1}$ of a swimmer versus $Re$ for the
power index $n=0.5$, 0.6, 0.7, 0.8, 0.9, and 1.

Figure 15: Velocity field next to the swimmer and swimmer position at $t=20$
for $Re=15$ and $n=0.7$ (top), $n=0.8$ (middle), and $n=0.9$ (bottom).

Applying algorithm (21)-(36) (with the subproblem (37) instead of problem (25)
for shear thinning fluids) and above discrete spaces, we have studied
numerically the motion of two disk swimmer in the $x_{1}$ direction in
Newtonian and non-Newtonian fluids. We first show the average velocities of
swimmer for different values of $Re$ and two power index values, $n=1$ and
0.8, obtained by the velocity mesh sizes $h=1/128$, 1/256, and 1/512 and time
step $\triangle t=0.001$. Each average value was computed from the last twenty
periods of horizontal velocity for $0\leq t\leq 20$. All those average
velocities for each Reynolds number are in a good agreement as shown in Figure
10.

Let us consider a swimmer of two same size disks, i.e.,
$R_{m}=R_{M}=2r_{0}=0.3$, suspended in a Newtonian fluid ($n=1$). Due to the
symmetry of reciprocal motion of the two same size disks (see Figure 11), the
swimmer mass center does not move as expected, at least for $Re=1$, 3, 5, 7,
9, 11, 13, and 15, the cases we have considered. But for the case of two
different size disks in a Newtonian fluid with $R_{M}=2r_{0}$ and
$R_{m}=r_{0}$ for $r_{0}=0.15$, the swimmer can move in the positive or
negative $x_{1}$-direction up to the Reynolds number. As shown in Figure 12,
the swimmer mass center oscillates with period 0.1 but moves in the positive
(resp., negative) $x_{1}$-direction for $Re=15$ (resp., $Re=5$). There is a
critical Reynolds number so that the swimmer moves to the left (resp., right)
for $Re<Re_{c}$ (resp., $Re>Re_{c}$) if the larger disk is located on the
right side (see the plots in Figures 13 and 14). The critical Reynolds number
is about $Re_{c}=12.5$ for the swimmer shown in Figure 13.

But when swimming in a non-Newtonian shear thinning fluid, the average moving
speed ($V_{1}$) per period in the positive-$x_{1}$ direction increases when
decreasing the value of the power index $n$ as shown in Figures 14 for the
power index values, $n=0.5$, 0.6, 0.7, 0.8, 0.9, and 1. Interestingly the
average swimming velocity is also slow down when the swimmer moves in the
minus-$x_{1}$ direction when decreasing the value of the power index $n$. Thus
the swimmer position is further to the right when comparing its position at
$t=20$ as in Figure 15. Our computational results show that the shear thinning
does have a strong effect on the motion of a swimmer made of two different
size disks and the critical Reynolds number decreases for a non-Newtonian
shear thinning fluid as the power index value decreases.

##  5 Conclusions

In this article, a distributed Lagrange multiplier/fictitious domain method
has been developed for simulating a non-symmetric (two-disk) particle moving
freely in non-Newtonian shear thinning fluids. The Carreau-Bird model is
adapted for the fluid shear thinning property. Our numerical methods have been
validated by comparing our computed solutions with the recently published
exact solutions for Poiseuille flows of non-Newtonian fluids in two dimensions
(in ref. [13]). Accurate numerical solutions have been obtained and the
$L^{2}$-errors of the velocity field show the expected 2nd order for $P_{1}$
finite elements. For many neutrally buoyant particle cross stream migration in
a Poiseuille flow of non-Newtonian shear thinning fluids, the particles
migrate away from the wall and center of the channel. The size of the boundary
gap decreases and that of the middle gap increases when reducing the power
index value. These results are consistent qualitatively with those reported in
[29]. Concerning the two-disk swimmer, the effect of shear thinning does let
the swimmer move faster when moving in the positive $x_{1}$-direction, i.e.,
from the smaller to larger disk, and the critical Reynolds number (for
changing the moving direction from the left to the right) decreases for
smaller value of the power index.

Acknowledgments

We acknowledge the helpful comments and suggestions of Howard H. Hu
(University of Pennsylvania).

## References

  * [1] Purcell EM. Life at low Reynolds number. Am J Phys. 1977;45(1):3-11.
  * [2] Lauga E, Powers TR. The hydrodynamics of swimming microorganisms. Rep Prog Phys. 2009;72(9):096601.
  * [3] Carreau PJ. Rheological Equations from Molecular Network Theories. Transactions of the Society of Rheology. 1972;16(1):99-127.
  * [4] Bird RB, Armstrong RC, Hassager O. Dynamics of polymeric liquids, Volume 1 Fluid Mechanics. 2 ed. New York, NY:John Wiley and Sons; 1987.
  * [5] Lauga E. Life at high Deborah number. Europhysics Letters. 2009;86(6):64001.
  * [6] Normand T, Lauga E. Flapping motion and force generation in a viscoelastic fluid. Phys Rev E. 78 (2008) 061907.
  * [7] Montenegro-Johnson TD, Smith DJ, Loghin D. Physics of rheologically enhanced propulsion: different strokes in generalized Stokes. Phys Fluids. 2013;25(8):081903.
  * [8] Qiu T, Lee TC, Mark AG, Morozov KI, Münster R, Mierka O,et al. Swimming by reciprocal motion at low Reynolds number. Nature communications. 2014;5(1):1-8.
  * [9] Dombrowski T, Jones SK, Katsikis G, Bhalla APS, Griffith BE, Klotsa D. Kinematics of a simple reciprocal model swimmer at intermediate Reynolds numbers. Phys Rev Fluids. 2020;5(6):063103.
  * [10] Dombrowski T, Klotsa D. Transition in swimming direction in a model self-propelled inertial swimmer. Phys Rev Fluids. 2019;4(2):021101.
  * [11] Huang PY, Feng J, Hu HH, Joseph DD. Direct simulation of the motion of solid particles in Couette and Poiseuille flows of viscoelastic fluids. J Fluid Mech. 1997;343:73-94.
  * [12] Huang PY, Hu HH, Joseph DD. Direct simulation of the motion of elliptic particles in Oldroyd-B fluids. J Fluid Mech. 1998;362:297-325.
  * [13] Griffiths PT. Non-Newtonian channel flow—exact solutions. IMA J Appl Math. 2020;85(2):263-279.
  * [14] Pan T-W, Glowinski R. Direct simulation of the motion of neutrally buoyant circular cylinders in plane Poiseuille flow. J Comp Phys. 2002;181(1):260-279.
  * [15] Pan T-W, Glowinski R. Direct simulation of the motion of neutrally buoyant balls in a three-dimensional Poiseuille flow. C R Mecanique, Acad Sci Paris. 2005;333(12):884-895.
  * [16] Yang BH, Wang J, Joseph JJ, Hu HH, Pan T-W, Glowinski R. Migration of a sphere in tube flow. J Fluid Mech. 2005;540:109-131.
  * [17] Pan T-W, Li A, Glowinski R. Numerical study of equilibrium radial positions of neutrally buoyant balls in circular Poiseuille flows. Phys Fluids. 2021;33(3):033301.
  * [18] Bristeau MO, Glowinski R, Periaux J. Numerical methods for the Navier-Stokes equations. Applications to the simulation of compressible and incompressible viscous flow. Computer Physics Reports. 1987;6(1-6):73-187.
  * [19] Glowinski R. Finite element methods for incompressible viscous flow. In: Ciarlet PG, Lions JL eds. Handbook of Numerical Analysis, vol. IX. Amsterdam: North-Holland; 2003:3-1176.
  * [20] Chorin AJ, Hughes TJR, McCracken MF, Marsden JE. Product formulas and numerical algorithms. Comm Pure Appl Math. 1978;31(2):205-256.
  * [21] Glowinski R, Osher SJ, Yin W eds. Operator-Splitting in Communications and Imaging, Sciences, and Engineering. Gewerbestrasse, Switzerland:Springer;2016.
  * [22] Glowinski R, Pan T-W, Hesla T, Joseph DD. A distributed Lagrange multiplier/fictitious domain method for particulate flows. Int J Multiphase Flow. 1999;25(5):755-794.
  * [23] Glowinski R, Pan T-W, Hesla T, Joseph DD, Periaux J. A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: Application to particulate flow. J Comput Phys. 2001;169(2):363-426.
  * [24] Glowinski R, Pan T-W, Periaux J. Distributed Lagrange multiplier methods for incompressible flow around moving rigid bodies. Comput Methods Appl Mech Engrg. 1998;151(1-2):181-194.
  * [25] J. Adams, P. Swarztrauber, and R. Sweet, FISHPAK: A package of FORTRAN subprograms for the solution of separable elliptic partial differential equations (The National Center for Atmospheric Research, Boulder, CO, 1980).
  * [26] Dean EJ, Glowinski R. A wave equation approach to the numerical solution of the Navier-Stokes equations for incompressible viscous flow. C R Acad Sci Paris, Série 1-Mathematics. 1997;325(7):783-791.
  * [27] Pan T-W, Glowinski R. A projection/wave-like equation method for the numerical simulation of incompressible viscous fluid flow modeled by the Navier-Stokes equations. Computational Fluid Dynamics Journal. 2000;9:28-42.
  * [28] Whiteman JR, Goodsell G. A survey of gradient superconvergence for finite element approximation to second order elliptic problems on triangular tetrahedral meshes. In: Whiteman JR, ed. The Mathematics of Finite Elements and Applications VII. London:Academic Press; 1991:55-74.
  * [29] Huang PY, Joseph DD. Effects of shear thinning on migration of neutrally buoyant particles in pressure driven flow of Newtonian and viscoelastic fluids. J. Non-Newtonian Fluid Mech. 2000;90(2-3):159-185.
  * [30] Segré G., Silberberg A. Radial particle displacements in Poiseuille flow of suspensions. Nature. 1961;189(4760):209-210.
  * [31] Segré G, Silberberg A. Behavior of macroscopic rigid spheres in Poiseuille flow. Part I, Determination of local concentration by statistical analysis of particle passages through crossed light beams. J Fluid Mech. 1962;14(1):115-157.


Microscopy and Microanalysis

Paper

[$\ast$]Corresponding author<EMAIL_ADDRESS>
# Calibrating coordinate system alignment in a scanning transmission electron
microscope using a digital twin

Dieter Weber0000-0001-6635-9567  David Landers0000-0002-3295-9473  Chen
Huang0000-0002-7864-8427  Emanuela Liberti0000-0001-5306-7936  Emiliya
Poghosyan0000-0002-0961-731X  Matthew Bryan0000-0001-9134-384X  Alexander
Clausen0000-0002-9555-7455  Daniel G. Stroppa0000-0002-7711-1839  Angus I.
Kirkland0000-0001-8068-1990  Elisabeth Müller0000-0001-8239-7109  Andrew
Stewart0000-0002-3081-5644  Rafal E. Dunin-Borkowski0000-0001-8082-0647  Ernst
Ruska-Centre for Microscopy and Spectroscopy with Electrons, Forschungszentrum
Jülich, 52425 Jülich, Germany  Department of Chemistry, University College
London, 20 Gordon Street, London, WC1H 0AJ, UK  LETI, CEA, 17 Av. des Martyrs,
38054 Grenoble, France  Rosalind Franklin Institute, Harwell Science and
Innovation Campus, Didcot, OX11 0QX, UK  Electron Microscopy Facility at PSI,
Paul Scherrer Institute, 5232 Villigen, Switzerland  DECTRIS AG, Taefernweg 1,
5405 Baden-Daettwil, Switzerland  Department of Materials, University of
Oxford, Oxford, OX1 3PH, UK

(2023)

###### Abstract

In four-dimensional scanning transmission electron microscopy (4D STEM) a
focused beam is scanned over a specimen and a diffraction pattern is recorded
at each position using a pixelated detector. During the experiment, it must be
ensured that the scan coordinate system of the beam is correctly calibrated
relative to the detector coordinate system. Various simplified and approximate
models are used implicitly and explicitly for understanding and analyzing the
recorded data, requiring translation between the physical reality of the
instrument and the abstractions used in data interpretation. Here, we
introduce a calibration method where interactive live data processing in
combination with a digital twin is used to match a set of models and their
parameters with the action of a real-world instrument.

###### keywords:

scanning transmission electron microscopy, 4D STEM, metadata

##  1 Introduction

For several 4D scanning transmission electron microscopy (STEM) data analysis
methods such as differential phase contrast (DPC) (Hamilton & Sheppard, 1984;
Lazić et al., 2016) or center of mass (CoM) (Lazić & Bosch, 2017), strain
(Usuda et al., 2005; Ozdol et al., 2015) and orientation (Viladot et al.,
2013; Rollett & Barmak, 2014) mapping, and ptychography (Hoppe, 1969; Yang et
al., 2017) it is essential that the relative orientation between the scan
coordinate system and the detector coordinate system is known so that
directions on the detector can be put into correspondence with directions in
the scan coordinate system (Ning et al., 2022). Figure 1 shows the impact of
an incorrect rotation calibration on integrated CoM reconstruction for
illustration. This calibration is surprisingly difficult and error-prone in
practice, and it was often not required in STEM previously, since conventional
detectors such as those used for bright field and annular dark field STEM are
rotationally symmetric. Consequently, there has been substantial recent work
in both theory and practical implementation to calibrate a microscope’s
coordinate system as well as validate the calibration of a given 4D STEM
dataset.

Figure 1: Impact of the rotation calibration on an integrated CoM (Lazić &
Bosch, 2017) reconstruction of a 4D STEM dataset of $\text{SmB}_{6}$, recorded
with a DECTRIS ARINA on a Thermo Fisher Spectra 200: a) The calibration of the
scan rotation is off by 90°. b) Correct calibration of the scan rotation.

Three methods are established to perform this alignment. First, the deflection
distribution around atom columns in atomic-resolution STEM can be analyzed.
Since the deflection is, in a thin specimen approximation, proportional to the
projected local field, a beam of negative electrons should be attracted to the
positively charged atom columns (Shibata et al., 2012). Furthermore, an
electrostatic field should be free of curl. The coordinate system alignment
that minimizes curl and ensures a negative divergence at atom positions is
assumed to be the correct one for this method (Savitzky et al., 2021).

Second, the self-consistency of the data in a ptychography reconstruction can
be analyzed, such as described in (Ning et al., 2022).

In practice, however, real-world specimens are usually not thin enough to
assume a simple direct dependence between displacement and field without
multiple scattering. In particular, slight misalignments between the atom
column axis and the beam axis as well as optical aberrations might lead to
apparent distortions (Bürger et al., 2020). Furthermore, atom columns are only
reliably resolved if the instrument aberrations are very low and a high
spatial resolution is achieved. Not every scanning transmission electron
microscope is capable of this, and it is not always desirable to tune and
operate the microscope in this mode. Additionally, methods based on analyzing
field distribution around atom columns are only applicable to thin crystalline
specimens.

The third and more conventional approach relies on over- or underfocusing the
beam so that the crossover above resp. below the specimen acts as a point
source. In this configuration, a shadow image of the specimen is projected
onto the detector. If this shadow image is compared with a STEM image, the
transformation between STEM image coordinates and shadow image, i.e. detector
coordinates, can be determined. Instead of a separate STEM image, a virtual
detector image can be generated from a 4D STEM dataset using a very small
virtual detector, for example a single detector pixel, so that it has a high
depth of field. In that case both detector shadow image and STEM image are
taken from the same 4D STEM dataset (Hu et al., 2023). Underfocus can be used
instead of overfocus, which leads to an inverted detector image compared to
overfocus.

This method relies on finding a field of view on the specimen without mirror
or rotational symmetry and adjusting the microscope so that it is clearly
recognizable, including exposing with sufficient dose. However, the matching
procedure can be error-prone if performed manually without the help of a
dedicated software stack. As an example, the direction of rotation is easily
reversed if the procedure is performed manually, an additional handedness
change could be missed, and different software or hardware uses different
coordinate system conventions. If conventional STEM and 4D STEM use a
different acquisition software and/or scan generator on the same microscope,
which is often the case, the coordinate system of a conventional STEM image
might differ from the coordinate system of a 4D STEM dataset.

Furthermore, different software might interpret and display the data
differently. As an example, some software displays the first pixel of a
detector frame in the lower left corner with the y axis pointing upwards,
while other image processing software shows the first pixel in the upper left
with the y axis pointing downwards. That may introduce inadvertent errors if
data is inspected visually with one software to calibrate the alignment, but
then processed numerically with another.

To eliminate such error sources, the calibration should be performed with the
same acquisition system and software as the actual acquisition to ensure
consistency. Furthermore, the calibration method and analysis methods should
be validated together to interpret data and parameters in a consistent way.

Here we present a method that is derived from the defocus method. It uses
automated data processing and a digital twin of the microscope to superimpose
all shadow images in an over- or underfocused 4D STEM dataset. If the
transformation by the digital twin matches the actual transformation by the
microscope, a sharp image is obtained. In case of a mismatch between
microscope and twin, the images are not superimposed correctly and the result
is blurred. This blurring is independent of symmetries in the specimen, since
the translation of the projected images on the detector caused by scanning the
beam breaks any symmetry.

##  2 Methods

The shadow imaging in the overfocused state is modeled with a linear ray
tracing approximation. All rays are assumed to originate from a focus point
above the specimen, pass through the specimen, and then hit the detector
plane. Underfocus, where the focus point is below the specimen, is represented
with a negative overfocus value that leads to the correct result without
separate adaptation of this model. Scanning the beam shifts the position of
the focus point laterally relative to the specimen, but doesn’t shift the beam
position on the detector or lead to a beam tilt in the specimen plane. The
specimen is assumed to be thin and oriented orthogonal to the optical axis,
i.e. not tilted or warped. All coordinate systems are right-handed by default.

Figure 2: Abstract model of the projection with overfocused STEM with the
parameters and coordinate axes labeled. Figure 3: Visualization of the ray
paths in the digital twin. The figure is generated using TEMGYM’s
visualization routine from the model used in the calculations with the actual
dataset (Figures 5-9). The overfocus and the scan step are exaggerated by a
factor of $10^{4}$ to make them visible. Note how this model is closer to the
physical reality of a scanning transmission electron microscope than the
abstract model in Figure 2. TEMGYM Basic translates the parameters of the
abstract model into parameters for simulated optical elements such as
deflectors and lenses to yield results equivalent to the abstract model.

With these assumptions in place, the electron-optical setup of this imaging
mode is described with the following parameters:

  * •

Overfocus Virtual vertical distance of the beam focus over the specimen plane
under the assumption of straight propagation between focus point and specimen.
Underfocus is specified as a negative overfocus value.

  * •

Scan Pixel Size Lateral distance between scan points. The scan is assumed to
be rectangular with uniform step size in both directions.

  * •

Camera Length Virtual vertical distance between specimen and detector under
the assumption of straight ray propagation between specimen and detector. In
real microscopes the electrons pass through a set of projection lenses that
allow adjusting this virtual distance.

  * •

Detector Pixel Size Size resp. pitch of detector pixels. The detector is
assumed to be rectangular with equidistant square pixels without gaps.

  * •

Y and X Center Pixel position on the detector that is defined as the beam
center, i.e. where the ray from the focus point straight down orthogonal
through the specimen plane hits the detector.

  * •

Scan Rotation Rotation between the scan coordinate system in scan pixel units
and the detector pixel coordinate system in pixel units. Physically, this is
composed of three components: Actual rotation of the scan, rotation imposed by
the projection system, and rotation of the detector. In the simple model used
here these three rotations are cumulative and can be described with a single
parameter.

  * •

Flip Y Invert the detector row index (y) axis. This encodes a handedness
change between detector and scan. Any combination of axis inversion and
rotation for both detector and scan can be described with a combination of
Scan Rotation and Flip Y.

Figure 2 shows a schematic of this model of the experiment with labels for the
parameters.

The parameters are chosen to correspond with the usual parameters that are
available from a microscope’s user interface. However, it is not guaranteed
that they are also interpreted in the same way. As an example, for Thermo
Fisher microscopes the parameter “Scan Rotation” describes a clockwise
rotation of the STEM image on the screen, while our implementation interprets
it as a clockwise rotation of the scan coordinate system around the beam (Z)
axis, following the convention (Clausen et al., 2019) in the LiberTEM software
(Clausen et al., 2023b). That means the parameter “Scan Rotation” rotates in
opposite direction between the Thermo Fisher software and LiberTEM.

“Flip Y” was chosen to describe a handedness difference since it makes it easy
to account for an inverted Y axis of a detector.

The parameters Overfocus $o$, Scan Pixel Size $s$, Camera Length $l$ and
Detector Pixel Size $d$ are not independent, but their ratios define a single
scaling factor $M=\frac{sl}{do}$ between scan pixel coordinates and detector
pixel coordinates. Since Scan Pixel Size, Camera Length and Detector Pixel
Size are typically calibrated by other means or are fixed physical distances,
Overfocus is a natural choice to adjust the scaling factor.

TEMGYM Basic (Landers et al., 2023c) was used to construct a digital twin for
a STEM instrument (Figure 3) that can accept these parameters and interprets
them in a consistent way to match the abstraction shown in Figure 2. This
digital twin can calculate both the detector pixel position and the specimen
pixel position of rays as a function of diffraction angle and scan position.
With this information, the shadow image projected onto the detector can be
transformed into images in scan pixel coordinates. This model assumes that
changing the focus doesn’t affect rotation and handedness.

Figure 4: Illustration of the superposition of projected shadow images: a)
Projection at first scan position. b) Projection at second scan position. c)
Projected image from a). d) Projected image from b). e) Incorrect
superposition of the projected images: The features don’t align. This will
lead to blurring if many images from different scan positions are superimposed
incorrectly. f) Correct superposition: The features align and the
superposition of many images from different scan positions creates a sharp
image of the object.

If the digital twin describes the action of the microscope correctly,
superimposing the images from an overfocused 4D STEM experiment after
transforming into scan pixel coordinates will form a sharp image of the
specimen (Figure 4). This superimposed image will also correspond to a virtual
STEM detector image formed from a single detector pixel, for example the
central ray like described in (Hu et al., 2023). However, if the
transformation of the digital twin differs from the actual microscope, the
images will not be superimposed correctly, leading to a blurred result.

By adjusting the parameters of the digital twin while continuously performing
the superposition operation, a user or algorithm can optimize the model
parameters until the result is as sharp as possible. A possible metric for
sharpness resp. blurriness was described by Crete et al. (2007). The process
of optimizing sharpness is similar to conventional focusing and two-fold
astigmatism adjustment, meaning it is intuitive for microscope users.

The data processing routine is implemented as a LiberTEM (Clausen et al.,
2023b) user-defined function (UDF) called “OverfocusUDF” that performs the
superposition and extracts the additional control images. This allows for fast
offline and live processing together with live visualization and dynamic
parameter updates.

For a smooth user interaction the superposition should be performed very
quickly. TEMGYM Basic does trace rays at high speed, but it is still too slow
to determine the target scan coordinates for each data point in a 4D STEM
dataset in real time. For that reason a linear transformation from detector
pixel position and scan position to specimen pixel position is derived from a
set of sample rays using a least square optimization. This transformation is
applied using an optimized kernel written in Numba.

In many cases the initial parameters will be far off their ideal values,
resulting in a featureless superposition image. Since multiple parameters have
to match at least approximately in order to form an image with recognizable
features, two helper tranformations were implemented:

By selecting a single detector pixel and plotting its value as a function of
scan position, a virtual detector image with very large depth of focus can be
obtained, as described in (Hu et al., 2023) (Figures 5 \- 8, label
“OverfocusUDF: point”). As long as the dose and detector resolution is
sufficient, this will be a recognizable image of the specimen in scan
coordinates. Second, a single shadow image on the detector can be transformed
to scan coordinates without superimposing it with shadow images from other
scan positions (Figures 5 \- 8 label “OverfocusUDF: selected”).

If the parameters are adjusted with sufficient accuracy (Figures 8 and 9),
these two images will look the same. By visually comparing them, a user can
roughly adjust the overfocus value by comparing the scale, and the rotation
and handedness by observing a unique feature on the specimen (Figures 5 – 7).
This coarse adjustment will usually be accurate enough to observe features in
the superimposed image (Figure 7), allowing for further incremental
optimization of all parameters following a protocol similar to optimizing
focus and astigmatism.

The model currently doesn’t account for descan error. However, the presence of
descan error is checked by summing up detector images without transformation.
Without descan error, the beam should have the same position on the detector
independent of the scan position. That means that a sum of all detector images
should yield a sharp image of the aperture that forms the primary beam.

Additionally, a loss function suitable for numerical optimization can be
derived from the calculation of the blurriness of the superimposed image as a
function of Scan Rotation and Overfocus. A utility function is included with
the “OverfocusUDF” that creates such a loss function by wrapping other
parameters for executing it, such as LiberTEM Context, dataset and the
remaining model parameters, in a closure (Sussman & Steele Jr., 1998). The
returned function accepts only the Scan Rotation and Overfocus as a single
vector of two values, and maps them from a range of [(-10, 10), (-10, 10)] to
a sensible interval of $\pm 10\deg$ resp. $\pm 25\%$ around their starting
value for conditioning. This ensures that the parameters are in a convenient
range for numerical optimizers and have roughly the same impact on blurriness.
That way, common optimizer implementations work with their default parameters.
See section “Data and code availability” for a link to the source code and
examples.

The test dataset presented in Figures 5 \- 9 was obtained using a gold cross-
grating replica sample. Data acquisition was conducted on a JEOL GRANDARM2
equipped with a MerlinEM 4R (quad-chip Medipix) system. The physical pixels on
the detector chip are 55 µm wide. The scanning coils of the microscope are
synchronized with the detector through TTL signals via a BNC cable. The scan
array consists of 64 by 64 steps, with each step measuring 12.5 nm. The
frames, each comprising 512 by 512 pixels, were each exposed for 1ms,
culminating in a total acquisition time of approximately 4 seconds.

The Supplementary Material (Weber et al., 2023) contains a full screen capture
video with explanation of the very first live adjustment using this method.
The test specimen was a semiconductor heterostructure in the bump-bonding
region of a PCB, where a random area was selected for the preparation of a TEM
lamella using a focused ion beam (FIB) instrument. A prototype of the DECTRIS
ARINA detector mounted on a probe corrected JEOL JEM-ARM200F transmission
electron microscope was used to record the live data. The microscope was
operated at 200kV acceleration voltage in STEM mode with a probe of 28.2 mrad
semi-convergence angle and 116 pA beam current. The ARINA detector was
operated at 50 kHz frame rate with (2x2) binning which resulted in 96x96
detector pixels (detector pixel size: 100 um). This was performed using a
prototype implementation of the method described here. Most notably, it didn’t
use TEMGYM and accelerated raytracing yet, but a preliminary manual ray
tracing implementation of the model described above. The detector data was not
recorded due to the ephemeral nature of live processing.

##  3 Results

A reliable manual alignment could be found in less than 15 min on the very
first attempt with this approach using a prototype implementation. See the
Supplementary Material at (Weber et al., 2023) for a screen capture video of
this first attempt with explanations, and a detailed description. Note how
live processing allowed adjustment of microscope parameters and specimen
movement interactively. Observing specimen movement allowed rough calibration
of scale, rotation and handedness even on a detector with a low number of
pixels such as the DECTRIS ARINA, which can only show a small region of
interest at coarse spatial resolution. We also note that the high frame rate
of the ARINA detector gave a sufficient scan repetition rate to observe the
effect of changes in near real time. The rotation and handedness parameters
obtained with the digital twin were validated to correspond to the parameters
of LiberTEM’s center of mass (CoM) analysis (Clausen et al., 2023a) and the
single side band (SSB) implementation in the Ptychography 4.0 project (Weber
et al., 2021). This enabled successful ptychographic reconstruction using this
SSB implementation in unrelated experiments following the adjustment in the
screen capture video. The results for Overfocus, Scan Pixel Size, Camera
Length and Detector Pixel Size have not been validated experimentally yet.

The Figures 5-8 show the manual adjustment procedure with the test dataset
described in the previous section using the most recent software version. The
dataset is available at (Weber et al., 2023).

Figure 5: Starting condition with scale and alignment miscalibrated. Control
interfaces for the calibration parameters are shown at the bottom. The field
“Blur metric” shows the blurriness of the “OverfocusUDF: shifted_sum” image as
calculated by the algorithm by Crete et al. (2007). The plot “OverfocusUDF:
point” shows the trace of a single detector pixel as a function of scan
position. “OverfocusUDF: shifted_sum” is the main adjustment plot that shows
all detector images superimposed in specimen coordinates according to the
transformation calculated by TEMGYM that is derived from the input parameters
below the plots. If the transformation by TEMGYM matches the actual
transformation by the microscope, this plot is a sharp image of the specimen.
Note that it is blurred here and doesn’t show any specimen features, meaning
the parameters are not correct yet. “SumUDF: intensity” shows a sum of all
untransformed diffraction patterns. The plot “OverfocusUDF: selected” shows
the diffraction pattern at the center of the scan grid rotated and scaled to
scan coordinates according to the current settings. Structures of the specimen
are only recognizable in “OverfocusUDF: point”, while “OverfocusUDF: selected”
magnifies the frame so much with the current settings that only a few pixels
are in the field of view. The sharp image of the beam-forming aperture shows
that the pivot point is adjusted correctly since the aperture is at the same
position for each scan point. “OverfocusUDF: sum” shows the sum of partially
transformed diffraction patterns where only “Scan Rotation” and “Flip Y” are
applied, in contrast to “SumUDF: intensity” that shows the sum of the
untransformed diffraction patterns where the position of the beam on the
detector can be judged. Figure 6: Compared to Figure 5, the “Overfocus”
parameter was adjusted so that “OverfocusUDF: point” and “OverfocusUDF:
selected” have approximately the same scale. Figure 7: Compared to Figure 6,
the “Scan Rotation” and “Flip Y” parameters were adjusted so that
“OverfocusUDF: point” and “OverfocusUDF: selected” show approximately the same
image. Blurred specimen features are starting to appear in “OverfocusUDF:
shifted_sum”. Figure 8: Compared to Figure 7, the “Overfocus” and “Scan
Rotation” parameters are fine-tuned so that “OverfocusUDF: shifted_sum”
appears sharp. Note how a small change in “Overfocus” and “Scan Rotation”
compared to Figure 7 causes a large difference in “OverfocusUDF: shifted_sum”,
demonstrating the sensitivity of this method. Figure 9: Compared to Figure 8,
the “Overfocus” and “Scan Rotation” parameters are numerically optimized by
minimizing the blur metric using global optimization within an interval around
the values in Figure 8.

The processing takes 1.5 s for a single pass over the test file used in
Figures 5-9, a MIB 4D STEM dataset of size 64 $\times$ 64 $\times$ 512
$\times$ 512 with 16 bits per pixel, corresponding to a data rate of about 1.4
GB/s or 2730 frames per second, on an AMD EPYC 7F72 with 24 physical cores and
256 GB RAM using LiberTEM 0.13 and the default Dask executor. Processing takes
1.38 s, corresponding to an uncompressed data rate of about 3.5 GB/s or 47,000
frames per second, for a file of size 256 $\times$ 256 $\times$ 192 $\times$
192 with 16 bits per pixel in the compressed HDF5 format generated by the
DECTRIS ARINA. This is fast enough to follow typical 4D STEM data rates live.
Live processing allows adjustment of both microscope parameters, such as
specimen region, magnification and focus, and model parameters. However, many
repeated scans may degrade a specimen through radiation damage or
contamination.

A loss function using the blurriness metric (Crete et al., 2007) as
implemented in the Scikit-Image (van der Walt et al., 2014) function
skimage.measure.blur_effect was created using the helper function described in
the previous section. It was minimized numerically using the “simplicial
homology global optimization” (SHGO) algorithm (Endres et al., 2018) with
“constrained optimization by linear approximation” (COBYLA, Powell (1994)) for
local search, as implemented in the SciPy (Virtanen et al., 2020) function
scipy.optimize.shgo. The result is shown in Figure 9. The fast processing of
LiberTEM and the OverfocusUDF allows numerical optimization within about 12 s
from Figure 8 to Figure 9 on the system described above. The optimization
converged reliably to near optimal values from various starting values,
provided the optimum was within the optimization bounds around the starting
value as described in the previous section. The optimization was implemented
to update the plots live, which allows observing the convergence behavior and
checking if the optimization result is sensible.

The following workflow turned out to be practical:

  1. 1.

Select a specimen region, scan pixel size and other microscope parameters that
show good contrast and recognizable features in the “OverfocusUDF: point”
plot. The method only works at low or medium magnification, not at atomic
resolution, since it doesn’t take diffraction into account.

  2. 2.

Set the Scan Pixel Size, Camera Length and Detector Pixel Size to their
correct values, as far as known.

  3. 3.

Adjust the focus of the microscope so that specimen features are recognizable
on the detector. Ideally, the features on the detector have roughly the same
scale as in the “OverfocusUDF: point” plot.

  4. 4.

Set the Overfocus value according to the overfocus given by the microscope.
Specimen features should be recognizable in both the “OverfocusUDF: point” and
“OverfocusUDF: selected” plots now.

  5. 5.

If possible, adjust Overfocus so that features in “OverfocusUDF: point” and
“OverfocusUDF: selected” have roughly the same scale. This works well on
grids, for example.

  6. 6.

Adjust Scan Rotation and Flip Y so that movements of the specimen or low-
symmetry specimen features are consistent between “OverfocusUDF: point” and
“OverfocusUDF: selected”. A blurred image of the specimen should appear in
“OverfocusUDF: shifted_sum” now.

  7. 7.

Confirm that the scale is approximately the same between “OverfocusUDF: point”
and “OverfocusUDF: selected”.

  8. 8.

Fine-tune Scan Rotation and Overfocus until “OverfocusUDF: shifted_sum” is
sharp. Confirm that Flip Y is set correctly if the image remains blurred.

  9. 9.

Optionally, refine with numerical optimization. This works once “OverfocusUDF:
shifted_sum” appears only moderately blurred, meaning the starting values are
already close enough.

##  4 Discussion

Live processing was helpful for finding a suitable specimen region and
microscope parameters as well as coarse adjustment of Scan Rotation and Flip
Y, since observing the effect of specimen movement interactively helps to
orient the user. For fine adjustments it is advantageous to record a dataset
after coarse adjustment of microscope and alignment parameters that show
recognizable features in the superposition image. The user can then continue
to adjust on such an offline dataset. LiberTEM allows using offline and live
data interchangeably, which helps the implementation of such a workflow. This
avoids excessive exposure of the specimen. Furthermore, recording such an
overfocused calibration dataset before real measurements can allow later
validation of the parameters for data interpretation.

Currently, the implementation uses TEMGYM Basic for the digital twin, i.e. a
matrix optics approximation. However, the method is not limited to this: The
ray tracing engine is exchangeable as long as it accepts the same parameters,
meaning a more advanced model such as TEMGYM Advanced (Landers et al., 2023a,
b) could be used. This software traces electron trajectories through
electromagnetic fields and enables the inclusion of 3rd order and higher
aberration effects from the objective lens. If the internal interpolation for
quick mapping of detector frames is adapted to allow non-linear relations
between detector pixel, scan position and specimen pixel, this approach could
also be extended to take aberrations into account. Furthermore, the digital
twin could be extended to match each optical element of the electron
microscope exactly, as already demonstrated in Landers et al. (2023c).

The method provides a clear validation for parameters that leaves little room
for error or interpretation: If the resulting image is sharp, the parameters
of microscope and digital twin are consistent. It can be performed as part of
a 4D STEM acquisition workflow without switching detector or scan engine.
Furthermore, it is sufficiently fast to perform many checks during a session,
for example after a change of microscope parameters. In particular if a good
initial guess is available, fine adjustment is quick and intuitive, and can
even be performed automatically. Since the image used for fine calibration is
composed of many detector frames, it has a better signal to noise ratio than
individual detector frames or the plot of the central pixel, allowing quick
fine calibration at a higher scan rate and low beam current if a good starting
value is known.

The current user interface in a Jupyter notebook works as a proof of concept.
In the future, a stand-alone tool could be developed or the method could be
integrated into the microscope control software. In particular, the actual
microscope parameters could be read out digitally and applied to the digital
twin. If microscope calibration values instead of model parameters are
optimized, it would allow the microscope to operate directly in units that
correspond to the digital twin.

The approach is also extensible to other parameters of the microscope that
change the detector image in a systematic way and allow superposition of
transformed detector images to form a sharp image of an object in the beam,
such as a specimen or aperture. As a disadvantage, it only works with 4D STEM
detectors, not conventional or segmented STEM detectors.

It should be investigated if changing the focus strongly in real microscopes
can affect the observed rotation or other parameters.

In this paper we already show how parameters can be transformed between three
different models or digital twins that describe the projection performed by
the microscope: TEMGYM Basic uses ray tracing through a simplified STEM model
based on a linear approximation for an objective lens and double beam
deflectors (Figure 3). The parameters exposed towards the user (controls in
Figure 5-9) describe a simpler model where the beam is not deflected by double
deflectors or focused by a lens, but where the specimen coordinate system is
simply shifted and a focus position is specified (Figure 2). Internally, the
“OverfocusUDF” uses a simple ray transfer matrix that directly performs a
mapping from detector coordinates and scan position to specimen pixel
coordinates for optimal performance.

Since each model is linear, and provides a good approximation of the actual
instrument’s behaviour under the chosen conditions, these are equivalent
descriptions of the same transformation that are specified in different ways
to match the context in which they are used. It also allows matching them
exactly using computational methods. The user interface aids the user in
finding approximate parameters that match the actual projection by the
instrument. Numerical optimization can then be used to fine-tune this result
further. In effect, this determines the instrument parameters in a reliable
and precise way.

Currently, the commonly used implicit models and parameters for metadata
schemas and data analysis only approximate the action of the instrument in
specific optical modes such as transmission electron microscopy (TEM), STEM,
electron energy loss spectroscopy (EELS), etc. The model in Figure 2 can
describe strongly defocused STEM, but can’t model the projection for in-focus
STEM since it doesn’t consider diffraction, for example. If a more
comprehensive or even full physical model is used to describe the projection
by the microscope, the approach described here can allow transformation of
physical parameters, such as objective lens current, into parameters for any
simplified abstract model, such as the STEM model here. It also allows
transformation of parameters between different models, for example to
determine the parameters on a different microscope to perform an equivalent
projection. A set of comprehensive model, parameters and calibration can
describe the action of the instrument in every detail and in every state in a
reproducible way, provided the model is sophisticated enough. At the same
time, corresponding parameters for any simplified model can be derived using,
for example, a numerical solver. That allows to use the same metadata schema
and model framework for any analysis performed with any TEM. That means a
comprehensive digital twin constitutes a potential foundation for universal
metadata for electron microscopy.

##  5 Code and data availability

The source code for the adjustment code as well as the Jupyter notebook to run
it on offline data can be found at https://github.com/LiberTEM/Microscope-
Calibration/tree/37066766ab446153ec2ebb3370cb879905cc410b.

The TEMGYM Basic version used here is
https://github.com/TemGym/TemGym/tree/bc28049225489b6ca75bc3377a4c64f53d77eb31.

The dataset used in the figures for offline calibration, a screen capture
video and additional explanation, as well as a runnable Apptainer image with
the complete software stack for reproducibility can be found at
https://doi.org/10.5281/zenodo.10418769.

##  6 Acronyms

STEM

    scanning transmission electron microscopy
DPC

    differential phase contrast
CoM

    center of mass
SSB

    single side band
UDF

    LiberTEM (Clausen et al., 2023b) user-defined function

##  7 Financial support

Swiss National Science Foundation is kindly acknowledged for co-funding the
electron microscope (R’Equip Project 206021_177020) used herein.

We gratefully acknowledge support by the AIDAS project.

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under grant agreement No 823717 – ESTEEM3.

We gratefully acknowledge support by DECTRIS for providing access to a
prototype of the ARINA detector at PSI.

##  8 Competing interests

The authors declare no competing interests.

## References

  * Bürger et al. (2020) Bürger, J., Riedl, T. & Lindner, J. (2020). Influence of lens aberrations, specimen thickness and tilt on differential phase contrast STEM images, Ultramicroscopy 219, 113118.
  * Clausen et al. (2019) Clausen, A., Weber, D., Bryan, M., Ruzaeva, K., Migunov, V., Baburajan, A., Bahuleyan, A., Caron, J., Chandra, R., Dey, S., Halder, S., Katz, D.S., Levin, B.D., Nord, M., Ophus, C., Peter, S., Puskás, L., Schyndel van, J., Shin, J., Sunku, S., Ånes, H.W., Müller-Caspary, K. & Dunin-Borkowski, R.E. (2019). LiberTEM concepts: Coordinate system, URL https://libertem.github.io/LiberTEM/concepts.html#coordinate-system.
  * Clausen et al. (2023a) Clausen, A., Weber, D., Bryan, M., Ruzaeva, K., Migunov, V., Baburajan, A., Bahuleyan, A., Caron, J., Chandra, R., Dey, S., Halder, S., Katz, D.S., Levin, B.D., Nord, M., Ophus, C., Peter, S., Puskás, L., Schyndel van, J., Shin, J., Sunku, S., Ånes, H.W., Müller-Caspary, K. & Dunin-Borkowski, R.E. (2023a). LiberTEM API reference: ComUDF.with_params, URL https://libertem.github.io/LiberTEM/reference/udf.html#libertem.udf.com.CoMUDF.with_params.
  * Clausen et al. (2023b) Clausen, A., Weber, D., Bryan, M., Ruzaeva, K., Migunov, V., Baburajan, A., Bahuleyan, A., Caron, J., Chandra, R., Dey, S., Halder, S., Katz, D.S., Levin, B.D., Nord, M., Ophus, C., Peter, S., Puskás, L., Schyndel van, J., Shin, J., Sunku, S., Ånes, H.W., Müller-Caspary, K. & Dunin-Borkowski, R.E. (2023b). LiberTEM/LiberTEM: 0.13.1.
  * Crete et al. (2007) Crete, F., Dolmiere, T., Ladret, P. & Nicolas, M. (2007). The blur effect: perception and estimation with a new no-reference perceptual blur metric, B.E. Rogowitz, T.N. Pappas & S.J. Daly (eds.), Human Vision and Electronic Imaging XII, SPIE.
  * Endres et al. (2018) Endres, S.C., Sandrock, C. & Focke, W.W. (2018). A simplicial homology algorithm for Lipschitz optimisation, Journal of Global Optimization 72, 181–217.
  * Hamilton & Sheppard (1984) Hamilton, D.K. & Sheppard, C.J.R. (1984). Differential phase contrast in scanning optical microscopy, Journal of Microscopy 133, 27–39.
  * Hoppe (1969) Hoppe, W. (1969). Beugung im inhomogenen Primärstrahlwellenfeld. I. Prinzip einer Phasenmessung von Elektronenbeungungsinterferenzen, Acta Crystallographica Section A 25, 495–501.
  * Hu et al. (2023) Hu, Y., Gao, S., Wu, X., Pei, X., Huang, F., Mao, W., Zhang, W., Horne, A., Gu, Z. & Wang, P. (2023). 4D-Explorer: A visual software for 4D-STEM data processing and image reconstruction, arXiv .
  * Landers et al. (2023a) Landers, D., Clancy, I., Dunin-Borkowski, R.E., Weber, D. & Stewart, A. (2023a). TEMGYM advanced: Software for electron lens aberrations and parallelised electron ray tracing, Ultramicroscopy 250, 113738.
  * Landers et al. (2023b) Landers, D., Clancy, I., Dunin-Borkowski, R.E., Weber, D. & Stewart, A.A. (2023b). TEMGYM advanced – NanoMi lens characterisation, Micron 169, 103450.
  * Landers et al. (2023c) Landers, D., Clancy, I., Weber, D., Dunin-Borkowski, R.E. & Stewart, A. (2023c). TEMGYM basic: transmission electron microscopy simulation software for teaching and training of microscope operation, Journal of Applied Crystallography 56, 1267–1276.
  * Lazić & Bosch (2017) Lazić, I. & Bosch, E.G. (2017). Analytical review of direct STEM imaging techniques for thin samples, Advances in Imaging and Electron Physics, 75–184, Elsevier.
  * Lazić et al. (2016) Lazić, I., Bosch, E.G. & Lazar, S. (2016). Phase contrast STEM for thin samples: Integrated differential phase contrast, Ultramicroscopy 160, 265–280.
  * Ning et al. (2022) Ning, S., Xu, W., Ma, Y., Loh, L., Pennycook, T.J., Zhou, W., Zhang, F., Bosman, M., Pennycook, S.J., He, Q. & Loh, N.D. (2022). Accurate and robust calibration of the uniform affine transformation between scan-camera coordinates for atom-resolved in-focus 4D-STEM datasets, Microscopy and Microanalysis 28, 622–632.
  * Ozdol et al. (2015) Ozdol, V.B., Gammer, C., Jin, X.G., Ercius, P., Ophus, C., Ciston, J. & Minor, A.M. (2015). Strain mapping at nanometer resolution using advanced nano-beam electron diffraction, Applied Physics Letters 106.
  * Powell (1994) Powell, M.J.D. (1994). A Direct Search Optimization Method That Models the Objective and Constraint Functions by Linear Interpolation, 51–67, Springer Netherlands.
  * Rollett & Barmak (2014) Rollett, A.D. & Barmak, K. (2014). Orientation Mapping, 1113–1141, Elsevier.
  * Savitzky et al. (2021) Savitzky, B.H., Zeltmann, S.E., Hughes, L.A., Brown, H.G., Zhao, S., Pelz, P.M., Pekin, T.C., Barnard, E.S., Donohue, J., DaCosta, L.R., Kennedy, E., Xie, Y., Janish, M.T., Schneider, M.M., Herring, P., Gopal, C., Anapolsky, A., Dhall, R., Bustillo, K.C., Ercius, P., Scott, M.C., Ciston, J., Minor, A.M. & Ophus, C. (2021). py4DSTEM: A software package for four-dimensional scanning transmission electron microscopy data analysis, Microscopy and Microanalysis 27, 712–743.
  * Shibata et al. (2012) Shibata, N., Findlay, S.D., Kohno, Y., Sawada, H., Kondo, Y. & Ikuhara, Y. (2012). Differential phase-contrast microscopy at atomic resolution, Nature Physics 8, 611–615.
  * Sussman & Steele Jr. (1998) Sussman, G.J. & Steele Jr., G.L. (1998). Scheme: A interpreter for extended lambda calculus, Higher Order Symbolic Computation 11, 405–439.
  * Usuda et al. (2005) Usuda, K., Numata, T., Irisawa, T., Hirashita, N. & Takagi, S. (2005). Strain characterization in SOI and strained-Si on SGOI MOSFET channel using nano-beam electron diffraction (NBD), Materials Science and Engineering: B 124–125, 143–147.
  * van der Walt et al. (2014) van der Walt, S., Schönberger, J.L., Nunez-Iglesias, J., Boulogne, F., Warner, J.D., Yager, N., Gouillart, E. & Yu, T. (2014). scikit-image: image processing in Python, PeerJ 2, e453.
  * Viladot et al. (2013) Viladot, D., Véron, M., Gemmi, M., Peiró, F., Portillo, J., Estradé, S., Mendoza, J., Llorca-Isern, N. & Nicolopoulos, S. (2013). Orientation and phase mapping in the transmission electron microscope using precession-assisted diffraction spot recognition: state-of-the-art results, Journal of Microscopy 252, 23–34.
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S.J., Brett, M., Wilson, J., Millman, K.J., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., Carey, C.J., Polat, İ., Feng, Y., Moore, E.W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, C.R., Archibald, A.M., Ribeiro, A.H., Pedregosa, F., van Mulbregt, P., Vijaykumar, A., Bardelli, A.P., Rothberg, A., Hilboll, A., Kloeckner, A., Scopatz, A., Lee, A., Rokem, A., Woods, C.N., Fulton, C., Masson, C., Häggström, C., Fitzgerald, C., Nicholson, D.A., Hagen, D.R., Pasechnik, D.V., Olivetti, E., Martin, E., Wieser, E., Silva, F., Lenders, F., Wilhelm, F., Young, G., Price, G.A., Ingold, G.L., Allen, G.E., Lee, G.R., Audren, H., Probst, I., Dietrich, J.P., Silterra, J., Webber, J.T., Slavič, J., Nothman, J., Buchner, J., Kulick, J., Schönberger, J.L., de Miranda Cardoso, J.V., Reimer, J., Harrington, J., Rodríguez, J.L.C., Nunez-Iglesias, J., Kuczynski, J., Tritz, K., Thoma, M., Newville, M., Kümmerer, M., Bolingbroke, M., Tartre, M., Pak, M., Smith, N.J., Nowaczyk, N., Shebanov, N., Pavlyk, O., Brodtkorb, P.A., Lee, P., McGibbon, R.T., Feldbauer, R., Lewis, S., Tygier, S., Sievert, S., Vigna, S., Peterson, S., More, S., Pudlik, T., Oshima, T., Pingel, T.J., Robitaille, T.P., Spura, T., Jones, T.R., Cera, T., Leslie, T., Zito, T., Krauss, T., Upadhyay, U., Halchenko, Y.O. & Vázquez-Baeza, Y. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python, Nature Methods 17, 261–272.
  * Weber et al. (2023) Weber, D., Landers, D., Huang, C., Liberti, E., Poghosyan, E., Bryan, M., Clausen, A., Stroppa, D.G., Kirkland, A.I., Müller, E., Stewart, A. & Dunin-Borkowski, R.E. (2023). Supplementary material for: Calibrating coordinate system alignment in a scanning transmission electron microscope using a digital twin., URL https://doi.org/10.5281/zenodo.10418769.
  * Weber et al. (2021) Weber, D., Lesnichaia, A., Strauch, A., Clausen, A., Bangun, A., Melnyk, O., Meissner, H., Ehrig, S., Wendt, R., Sukumaran, M., Wollgarten, M., zu Castell, W. & Müller-Caspary, K. (2021). Ptychography 4.0: 0.1.0.
  * Yang et al. (2017) Yang, H., MacLaren, I., Jones, L., Martinez, G.T., Simson, M., Huth, M., Ryll, H., Soltau, H., Sagawa, R., Kondo, Y., Ophus, C., Ercius, P., Jin, L., Kovács, A. & Nellist, P.D. (2017). Electron ptychographic phase imaging of light elements in crystalline materials using Wigner distribution deconvolution, Ultramicroscopy 180, 173–179.

  *[STEM]: scanning transmission electron microscopy
  *[DPC]: differential phase contrast
  *[CoM]: center of mass
  *[UDF]: LiberTEM (, ) user-defined function
  *[SSB]: single side band