STUDY OF ACTIVE NEMATICS: CONTINUUM THEORY AND

PARTICLE SIMULATIONS

By

Sheng Chen

A DISSERTATION

Submitted to

Michigan State University

in partial fulﬁllment of the requirements

for the degree of

Mechanical Engineering – Doctor of Philosophy

2021

ABSTRACT

STUDY OF ACTIVE NEMATICS: CONTINUUM THEORY AND

PARTICLE SIMULATIONS

By

Sheng Chen

The concept of ‘active matter’ refers to a system that is far away from

equilibrium. It comprises internal self-driven units that consume a local fuel

(e.g., chemical fuel) from surroundings and transforms it into mechanical

work.

It is ubiquitous not only in the macroworld such as ﬂocking of ﬁsh,

birds or animal herds, but also pervasive in the microworld. Examples include

biological swimmers such as bacteria and microalgae, synthetic colloidal

surfers actuated by chemical reactions, as well as puriﬁed biopolymers (e.g.,

microtubules (MTs)) mixed with molecular motors.

In contrast

to the

previously well-studied passive systems where the instability mainly comes

from the thermodynamic ﬂuctuations, the inherent spontaneity of the active

matter endows itself with a complex and ever-changing dynamics and

structure, the study of which is still nascent.

In this thesis, we focus primarily on a speciﬁc type of active system in the

microscale that is featured by comprising self-driven particles with elongated

shape, i.e., the rod-like particles. Such system is described as ‘active nematics’

due to its resemblance to nematic liquid crystals. The out-of-equilibrium

pattern formation of active nematics is caused by the inextricable interplay

between the short-range steric interaction and the long-range hydrodynamics.

As a result, the dynamics and structure of active nematics display hallmarks of

collective motion of particles, chaotic ﬂow structure and the concomitantly

long-ranged nematic order and motile topological defect. To gain a physical

insight

to such complex while intriguing phenomena, we develop a

coarse-grained Q-tensor continuum theory coupled with low-Reynold ﬂuid

dynamics. With the assistance of computational framework for simulating

suspensions of rigid particles in Newtonian Stokes ﬂow, we are able to

conduct large-scale particle simulations to mimic many-particle couplings.

The thesis is organized as follows: In chapter 1, we study the complex

dynamics of a two-dimensional suspension comprising non-motile active

particles conﬁned in an annulus. A coarse-grained liquid crystal model is

employed to describe the nematic structure evolution, and hydrodynamically

couples with the Stokes equation to solve for the induced active ﬂows in the

annulus.

In chapter 2, We study ﬂuid and mass transport in a dilute apolar

active suspension conﬁned in a rectangular channel. By using a Galerkin

mixed ﬁnite element method, we are able to reveal various patterns of

spontaneous coherent ﬂows that can be unidirectional, traveling-wave, and

chaotic. In chapter 3, we study the long-time rotational Brownian diffusivity

in a crowded bath of hard rods with ﬁnite aspect ratios, where the topological

constraint dominate over hydrodynamics. In chapter 4, we study the nonlinear

dynamics of an undulatory microswimmer in a quasi-2D liquid-crystal

polymer solution.

ACKNOWLEDGMENTS

In the ﬁrst place, my sincere gratitude goes to my advisor Dr. Tong Gao

for his mentoring and support throughout my PhD period and the completion

of this dissertation. His expertise in ﬂuid mechanics and mathematics lays a

sound theoretical foundation for my research at Michigan State University and

I believe will deﬁnitely have a far-reaching impact on my future study. Also

thanks for his friendly advice on my career plans and his heartwarming support

during my hard times.

Besides this, I would like to thank my thesis committee members, Dr.

Yingda Cheng, Dr. Farhad Jaberi, and Dr. André Bénard for their support and

suggestions to this dissertation, and for their great course in computational

ﬂuid mechanics. I appreciate Dr. Chang Wang for his course in perturbation

theory. The perturbation techniques I learned from Dr. Wang and my advisor

Dr. Gao accounts for a considerable proportion of my theoretical study. I also

thank Dr. Charles Petty for his multiphase ﬂuids course that is beneﬁcial to my

research and thanks Dr. Indrek Wichman for being my teaching mentor.

I am also deeply grateful to my group colleagues: Dr. Andrew Michael Hess,

and Dr. Zhaowu Lin. Their kindly assistance, expertise in different research

ﬁeld contribute to the completion of different part of this thesis.

Last but not least, I would like to thank my parents and my grandma for their

unconditional support and love. Thank Jing Ke, for being both a best friend and

a soulmate all the way through.

iv

TABLE OF CONTENTS

LIST OF FIGURES

vii

1 DYNAMICS AND STRUCTURE OF AN APOLAR ACTIVE

.

.
.

.
.

.
.

.
.

1
SUSPENSION IN AN ANNULUS
1.1 Abstract
. . . . . . . . . . . . . . . . . . . . . . . . . 1
. . . . . . .
.
.
.
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Introduction .
.
.
1.3 Mathematical model .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4.1 Dilute Suspension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4.2 Concentrated Suspension . . . . . . . . . . . . . . . . . . . . . . . . . 15
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

1.5 Conclusion .

.

.

.

.

.

.

2 TRANSPORT OF AN APOLAR ACTIVE SUSPENSION IN A

.

.
.

.
.

.
.

.
.

21
RECTANGULAR CHANNEL
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.1 Abstract
.
.
.
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.2 Introduction .
.
.
2.3 Mathematical model .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4 Analysis and simulation results . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.4.1 Spontaneous Coherent Flows . . . . . . . . . . . . . . . . . . . . . . . 27
2.4.2 Concatenation of Hydrodynamic Instabilities . . . . . . . . . . . . 28
2.4.3 Active Transport of Circular Cylinders . . . . . . . . . . . . . . . . 34
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

2.5 Conclusion .

.

.

.

.

.

.

3 SCALING LAW OF BROWNIAN ROTATION IN CONCENTRATED

.

.
.

.
.

.
.

.
.
.

3.4 Results .

HARD-ROD SUSPENSIONS
3.1 Abstract
.
.
.
.
3.2 Introduction .
.
.
3.3 Numerical method .

41
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.3.1 Complementarity formulation and equation of motion . . . . . 44
. . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.4.1 Equilibrium PDF of the Doi-Onsager Model . . . . . . . . . . . . 46
3.4.2 Measuring Dr from quadratic OACFs . . . . . . . . . . . . . . . . . 49
3.4.3 Scaling law delineated by modiﬁed tube model
. . . . . . . . . . 53
3.4.4 Approximations of the double-conﬁgurational averages . . . . 58
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

3.5 Conclusion .
3.6 Future plan .

. . . . . .

.
.

.
.

.
.

.
.

.
.

.
.

.

.

.

.

.

.

.

.

.

4 Q−TENSOR MODEL FOR UNDULATORY SWIMMING IN A

LIQUID CRYSTAL

61

v

.
.

.
.

.
.

.
.

. . . . . . . . . . . . . . . . . . . . . . . . . 61
4.1 Abstract
. . . . . . .
.
.
.
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.2 Introduction .
.
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.3 Mathematical Model .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.4 Numerical method .
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.5 Asymptotic analysis .
4.5.1 Isotropic cases .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.5.2 Nearly-aligned cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.6 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.7 Conclusion .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

.

.

.

.

.

.

.

5 PUBLICATION

BIBLIOGRAPHY

90

92

vi

LIST OF FIGURES

Figure 1.1: (a-c) Characteristic ﬂow patterns (streamline overlaid on the
colormap of vorticity Ω) and (d-f) nematic structures (nematic
director m ﬁeld overlaid on the colormap of scalar order
parameter s) when choosing R1 = 1.0 and R2 = 2.0. Insets in
(a) and (d) are the comparisons between the analytical (solid
lines) and numerical (open circles) results for normalized uθ
and Drθ components taken along the white lines when |α|
slightly goes beyond the critical value |α|cr = 1.23 for the
(g) Phase diagram of
Isotropy-Circulation instability.
(|α|, R2 − R1) when ﬁxing R1 = 1.0. The solid and dashed
lines characterize the Isotropy-Circulation instability in an
annulus and a straight channel, respectively.The dash-dot line
characterizes the Circulation-Traveling wave instability in a
straight channel for sharply-aligned rods by neglecting the
rotational diffusion.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

instabilities

Figure 1.2: (a) Growth rate Re (σ) for the long-wave (diamonds) and
ﬁnite-wavelength (lines)
function of
wavenumber k for
Isotropy-Circulation instabilities in a
periodic straight channel when choosing h = 0.5,
dT = dR = 0.025. (b) Streamlines overlaid on the colormap
of u velocity component at α = −4.0. (c) Velocity proﬁles at
α = −2.0 (umax = 0.07) and α = −4.0 (umax = 0.13),
respectively. .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

as

a

.

.

.

Figure 1.3: (a) Growth rate Re (σ) as a function of wavenumber k for
bending instabilities of the sharply-aligned Extensors when
conﬁned in a periodic straight channel when choosing
h = 0.5, dR = 0 and dT = 0.025. The maximum growth rates
at the critical wavenumbers kcr are marked by solid circles.
(b) and (c) show the vorticity ﬁeld Ω as bending instability
start to grow at α = −4.0 and α = −8.0, respectively. The
solid black scales represent the length scales π/kcr that match
the separation distance between the neighboring vortices.

. . . . 10

vii

Figure 1.4: Equilibrium solutions of nematically-aligned states at different
values of α when choosing R1 = 1.0 and R2 = 2.0, nematic
director ﬁeld overlaid on the colormap of s. Inset in panel(d):
Comparisons between the analytical (solid line) and numerical
(open symbols) results at α = −0.5,−2.0 for Drr component.
The numerical data are taken along the white line marked in
panel(c).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

.

.

.

.

.

.

Figure 1.5: Evolution of the ﬂow patterns for a concentrated Extensor
suspension when choosing α = −2.0, R1 = 1.0 and R2 = 2.0.
Insets: nematic director ﬁeld overlaid on the colormap of s.

. . 17

Figure 1.6: Flow patterns (streamline overlaid on the colormap of Ω, upper
row) and nematic structures (nematic director ﬁeld overlaid on
the colormap of s, lower row) that characterize the long-time
dynamics of a conﬁned Extensor suspension when choosing
α = −2.0, R1 = 0.75 and R2 = 2.0.

. . . . . . . . . . . . . . . . . . 17

Figure 1.7: Short-time evolution of the ﬂow (velocity vector overlaid the
colormap of on Ω, upper row) and nematic (nematic director
overlaid on the colormap of s, lower row) ﬁelds for the case
shown in ﬁgure 1.6.
In (e-h), the black arrows shows the
direction of the crack formation and +1/2 defect movement. . . 18

Figure 2.1: Typical (quasi-)steady coherent ﬂows (a,c,e) and nematic
traveling wave
orders (b,d,f) for unidirectional ﬂow (a,b),
(c,d), and chaos (e,f).
Simulations are performed when
choosing the half channel width H = 0.5, dR = dT = 0.025,
and α = −2.0,−8.0,−20.0. Comparison of simulation and
analytical solutions of the velocity proﬁle is shown in inset of
(a) where the maximum velocity is umax = 0.57. (g) Phase
diagram of (|α|, 2H). The solid line are the marginal stability
curves from the analyses of the Q-U and U-T instabilities.
The
the
geometrically-selected high-order modes. The dotted line is
our previous prediction when using a sharp-alignment
assumption ([1]).

predictions

dashed

using

lines

the

are

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

viii

the channel width is chosen as H = 0.5;

Figure 2.2: Unstable modes and growth rates Re (σ) of linear stability
analyses for Q-U (a), U-T (b), and T-C (c) instabilities.
In
in (c),
(a,b),
α = −4.0 is ﬁxed. The open diamond symbols (in panel(a) at
k = 0) represent the long-wavelength (M0) mode associated
with the unidirectional ﬂow solution. The solid, dashed, and
dotted lines in (a-c) represent the 1st (M1), 2nd (M2), and
3rd (M3) unstable modes, corresponding to the single, double,
and triple vortex arrays, respectively. The colorﬁelds for
vorticity and streamlines exhibit ﬂow patterns, and are
selected at the critical wavenumbers (circular red dots).

. . . . . 32

Figure 2.3: Active ﬂow past a ﬁxed cylinder of R = 0.1 in the channel
center: rescaled (a) u and (b) v velocity component, and the
surface distribution of (c) pressure force, (d) viscous force, and
(e) active force. In (a)-(e), the upper row are the simulation
results when choosing α = −2.0, and the bottom row shows
the normalized analytical results.

. . . . . . . . . . . . . . . . . . . . . 34

Figure 2.4: Transport of free-moving cylinders in active ﬂows. (a) FEM
results of typical time evolutions of the migration velocities
for cylinders of different sizes when choosing α = −2.0 and
H = 0.5. (b) Solid lines represent the normalized steady-state
migration velocities as functions of re-scaled radius. Symbols
are the FEM results. (c) FEM results of cylinder trajectories
(R = 0.1) in traveling wave ﬂows.
(d) FEM results of
time-averaged velocity proﬁles at different cylinder volume
fractions when choosing R = 0.1 and H = 0.5. The solid and
dashed lines correspond to α = −4.0 and −6.0, respectively. . . 38

Figure 3.1: Equilibrium orientational order parameter S as a function of
φ. The open squares represent the GCBD simulation data.
The solid lines are calculated using Ψeq in Eq. (3.12). Inset:
Instantaneous snapshots at φ = 20% and 40% for γ = 10
correspond to the I and N phase, respectively. . . . . . . . . . . . . . 47

Figure 3.2: Examples of orientational PDF for γ

=

10 at
φ = 20%, 30%, 40%. Open squares represent the GCBD data,
and the solid lines represent the results by using Ψeq in Eq.(2)
in the main text.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

.

ix

(cid:68)

(p (t) · p (0))2(cid:69)

Figure 3.3: Orientation relaxation characterized by (cid:104)p (t) · p (0)(cid:105) (a) and
(b) across φ for γ = 20. The open symbols
Schematics of conﬁgurational
represent the GCBD data.
space of p(t) are highlighted at the bottom for the different
phases.
the ﬁtting
function in Eq. (3.13). Inset: Zoom-in window highlights that
the short-time relaxation of OACF ﬁts an exponential decay
governed by a single time scale τR. . . . . . . . . . . . . . . . . . . . . 50

the solid lines represent

In panel(b),

Figure 3.4: Dr/Dr0 as functions of n(cid:96)3 across φ for various γ. The open
squares represent the data extracted from GCBD simulations
via Eq. (3.18). The solid and dashed lines represent the
scaling laws in Eq. (3.20) and (3.25), respectively. The solid
red square on the γ = 20 curve marks the location of the data
for ﬁtting β. Inset: (cid:104)|p × p(cid:48)|(cid:105)p
p(cid:48) as functions of
S,
lines)
the
quadratic-closure
[2],
respectively.
.

p(cid:48) and (cid:104)|p · p(cid:48)|(cid:105)p
(dashed
. . . . . . . . . . . . . . . . . . . 53

using Ψeq
approximation
. . . . . . . . . . . . .

and
lines)

evaluated

(solid

by

.

.

Figure 3.5: GCBD simulations motivate a new tube model. (a) Schematic
of a virtual tube with penetrations from both the sidewall and
the two ends.
(b) Fractional contributions of the averaged
collision numbers on the sidewall (λI) and the ends (λII) for
γ = 5 and 20 as functions of φ. (c) The collision number ratio
λI/ (γλII) as functions of S for various γ.
(d) The OPCF
g(r⊥) in the nematic regime as functions of the dimensionless
transverse distance r⊥/b.

. . . . . . . . . . . . . . . . . . . . . . . . . . 56

Figure 4.1: Convergence tests of a parallel moving swimmer with the time-
dependent C.O.M. velocity U(Ux, Uy) by varying (a) domain
size (∆t = 6.25 × 10−5, h = 1/128), (b) Eulerian grid width
(Ωf = 2 × 2), and (c) Lagrangian grid spacing (Ωf = 2 × 2,
h = 1/128). These parameters are ﬁxed: P e = 10, P et = 100,
Er = 1, β = 0.005, ζ = 8, L = 1, A = 0.05, k = 2π, ω = 2π,
σb = 0.5, σs = 1500. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

Figure 4.2: Time-averaged C.O.M. speed UOB for undulatory swimming
motion in an Oldroyd-B ﬂuid. These parameters are ﬁxed:
Ωf = 8 × 8, h = 1/64, Ns = 64, ηp/ηs = 1/2, L = 4,
De = 1, and ∆t = 2.5 × 10−4. . . . . . . . . . . . . . . . . . . . . . . . 72

x

Figure 4.3: The time-dependent C.O.M. velocity U(Ux, Uy) for a parallel
moving swimmer when choosing different values of bending
stiffness σb. These parameters are ﬁxed: These parameters are
ﬁxed: Ωf = 2 × 2, ∆t = 6.25 × 10−4, h = 1/128, P e = 10,
P et = 100, Er = 1, β = 0.005, ζ = 8, L = 1, Ns = 32,
A = 0.05, k = 2π, ω = 2π, σs = 1500. . . . . . . . . . . . . . . . . . 73

Figure 4.4: I-N phase transition led anisotropic swimming.

(a) δ as a
function of ζ, with the bifurcation occurring at ζc = 4.
(b)Constant and apolar distribution of Ψ0 (φ) at ζ = 0 and 8.
(c) ULC/UN as a function of ζ for stiff (σb = 0.5) and soft
(σb = 0.005) swimmer when choosing P e = 1 and Er = 10,
for the parallel (denoted by “ (cid:107)”) and and perpendicular
(denoted by “⊥”) motions. . . . . . . . . . . . . . . . . . . . . . . . . . . 81

Figure 4.5: Speed ratio ULC/UN as a function of P e at Er = 1 and 10
obtained from numerical simulation (a-c) and asymptotic
solutions (d-f). .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

.

Figure 4.6: The instantaneous (a-d) and time-averaged (e-h) velocity ﬁeld
around a “stiff” swimmer in the body-ﬁxed frame when
choosing σb = 0.5, P e = 1, and Er = 10: (a,e) Newtonian;
(b,f) Isotropic, ζ = 2; (c,g) Nematic and parallel, ζ = 8; (d,h)
Nematic and perpendicular, ζ = 8. The vector ﬁeld and the
background color correspond to the velocity u and its
magnitude |u|.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

.

.

Figure 4.7: The time-averaged director and polymer force ﬁelds around a
stiff swimmer corresponding to Figs. 4.6(f-h). (a-c): nematic
director (cid:104)m(cid:105) superposed on the colormap of
the order
parameter S. (d-f): polymer force vector (cid:104)fp(cid:105) superimposed
on the force magnitude |(cid:104)fp(cid:105)|. . . . . . . . . . . . . . . . . . . . . . . . . 86

xi

1

DYNAMICS AND STRUCTURE OF AN APOLAR ACTIVE

SUSPENSION IN AN ANNULUS

1.1 Abstract

We study the complex dynamics of a two-dimensional suspension comprising

non-motile active particles conﬁned in an annulus. A coarse-grained liquid

crystal model is employed to describe the nematic structure evolution, and

hydrodynamically couples with the Stokes equation to solve for the induced

active ﬂows in the annulus. For dilute suspensions, coherent structures are

captured by varying particle activity and gap width, including unidirectional

circulations, traveling waves, and chaotic ﬂows. For concentrated suspensions,

the internal collective are featured by motile disclination defects and ﬂows at

ﬁnite gap widths. In particular, we observe an intriguing quasi-steady state at

certain gap widths during which +1/2-order defects oscillate around

equilibrium positions accompanying traveling-wave ﬂows

that

switch

circulating directions periodically. We perform linear stability analyses to

reveal the underlying physical mechanisms of pattern formations during a

concatenation of instabilities.

1.2 Introduction

Active matter is a class of far-away-from-equilibrium systems comprising

self-driven particles ([3, 4, 5]). When they are immersed in a ﬂuid, motile

1

microstructures exert forces upon the ambient liquid which itself acts as a

coupling medium for generation of multiscale dynamics, which presents

challenges in design, analysis, and control of novel active materials. To take

advantage of the anomalous properties (e.g.,

large-scale induced motion,

enhanced diffusion, and energy conversion), it is essential to guide active

matter to perform useful mechanical work. One way of doing this is to tune

the suspension concentration and the amount of chemical fuels ([6, 7, 8]).

Alternatively, it is possible to make use of the particle interactions, either

individually or collectively, with obstacles and geometric boundaries for more

direct manipulation.

By trapping active suspensions (such as Pusher

swimmers and Quincke rollers) within straight and curved boundaries,

experimental and numerical studies have revealed stable coherent structures

and ﬂow types ([9, 10, 11, 12, 13, 14, 15, 16, 17, 18]).

In this paper, we study the complex dynamics of a two-dimensional (2D)

apolar active suspension conﬁned in an annulus. The rod-like microparticles

are non-motile but mobile, and only advected by ﬂuid ﬂow. In the meantime,

they elongate through nearly symmetric stretching or growth to exert extensile

dipolar stresses on the solvent, and interact with each other

through

hydrodynamic coupling and steric alignment torques, which we referred to as

the “Extensor” particles ([19, 20]). Examples include bacterial cell division

([21]), microtubule (MT) bundles undergoing polarity sorting driven by

molecular motors ([7]), and tripartite Au-Pt-Au nanomotors generating surface

ﬂows due to catalytic reactions ([22, 23]). In the previous works by [19] and

[20], a macroscale liquid crystal model, i.e., “BQ”-tensor theory, has been

2

derived from a kinetic theory that describes Extensors’ ensemble dynamics.

Compared to the macro models for the motile suspensions ([24]) where the

coupled evolution equations for the suspension concentration, polarity, and

tensorial-order parameter need to be solved together, we have shown that the

apolar BQ-tensor model naturally conserves the local concentration and the

global particle numbers, and can be characterized by the evolution equation of

the second-moment tensor only. Through theoretical analyses and numerical

simulations for Extensor suspensions in both unbounded domains and

conﬁned in circular chambers, we have demonstrated that this active ﬂuid

model, while being apolar, inherits all the basic transitions and instabilities

associated with motile suspensions, and exhibits a rich set of collective

dynamics.

We examine both the dilute and concentrated Extensor suspensions conﬁned

in an annulus geometry. For the dilute cases, we identify emergent coherent

structures that resemble those observed in “polar” active suspensions of Pusher

particles ([25, 17]) by varying particle activity and annulus gap width. Linear

stability analyses are performed to reveal the underlying physical mechanisms

of a series of hydrodynamic instabilities starting from a near-isotropic state. At

a ﬁnite concentration, we show that the internal collective dynamics can be

effectively inhibited when conﬁned in thin gaps; they rise when the gap width

is large enough so that the bending instability can develop in the radial direction

to generate active nematic ﬂows. In particular, we capture an intriguing quasi-

steady state at a ﬁnite gap width featured by oscillating +1/2-order defects and

traveling-wave ﬂows that switch circulating directions periodically.

3

1.3 Mathematical model

We consider a collection of Extensors particles suspended in a Newtonian

ﬂuid. These active particles are non-motile but can elongate or stretch near

symmetrically to produce extensional ﬂows that effectively exert dipolar

stresses upon the liquid. As shown by [19, 20], the ensemble dynamics of

Extensor particles whose center-of-mass position locates at x with an
orientation p (|p| = 1) can be described by a probability distribution function
Ψ(x, p, t) through a Smoluchowski equation. The particles exert stresses ΣΣΣ on

the liquid to drive ﬂuid motion which is governed by the (dimensionless)

incompressible Stokes equation:

∇p − ∆u = ∇ · ΣΣΣ,

∇ · u = 0,

(1.1)

(1.2)

where u and p are the ﬂuid velocity and pressure, respectively. The extra stress

tensor ΣΣΣ can be written as ([26, 27]):

ΣΣΣ = αD + βS : E − 2ζβ (D · D − S : D)

(1.3)

where E = 1/2(cid:0)∇u + ∇uT(cid:1) is the rate-of-strain tensor; D =(cid:82)
S = (cid:82)

p ppΨdp and
p ppppΨdp are the second- and fourth-moment tensors respectively,
and are calculated by taking the average of p on the surface of a unit sphere.

In equation (1.3), the ﬁrst term represents a dipolar extensile stress that arises

from the local extensional ﬂows with the dimensionless coefﬁcient α < 0

quantifying particle activity ([28, 26]). The second term is a constraining

stress due to particle rigidity ([29, 30]) with an effective shape factor β > 0

characterizing particle’s aspect ratio.

The third term is due to steric

4

interactions through a Maier-Saupe potential with a strength coefﬁcient ζ

([31, 30, 27]). Note that both the constraining and steric interaction stresses

are proportional to the effective volume fraction of the suspension, and hence

cannot be neglected for concentrated suspensions ([30, 27]).

The ﬂow equation is coupled with the evolution equation of D tensor, which

is derived from the kinetic model by taking a standard moment average ([29,

19, 20]):

(cid:18)

(cid:19)

D∇ + 2E : S = 4ζ (D · D − S : D) + dT ∆D − 4dR

∂t + u · ∇D −(cid:0)∇u · D + D · ∇uT(cid:1) is an upper-convected

where D∇ = ∂D
derivative, dR and dT are the rotational and translational diffusion coefﬁcients,

D − I
2

(1.4)

respectively.

It is worthwhile to mention that in the classical Landau-de

Gennes approach for liquid crystalline ﬂuids ([32]), dR and dT are essentially

the linear terms in a Landau potential and the Frank elastic constant for a

nematic liquid crystal respectively. In nondimensionalization, we assume there
are N Extensor particles of length l and width b (b/l (cid:28) 1) distributed in a
volume V , and introduce an effective volume fraction ν = nbl2 where

n = N/V is the mean number density ([30, 19, 20]). Then we choose the
length scale lc = b/ν, the velocity scale |u0| which represents the surface
velocity due to elongation or stretching motion, and the stress scale µ|u0|/lc.
Equation (2.2) can be closed by expressing the fourth-moment tensor S in

terms of D through the so-called Bingham closure ([33, 34, 19]) which

reconstructs a distribution function ΨB(x, p, t) in terms of a traceless

5

symmetric tensor T(x, t) to yield

ΨB[T] =

(cid:82)

exp (T : pp)

p exp (T : pp) dp

.

(1.5)

relation

We

determine

tensor T numerically

by

solving

the

p ΨB[T]ppdp. Given T,

the fourth-moment

tensor S is then

D = (cid:82)
approximated by SB =(cid:82)
tensor-order parameter Q(x, t) = D/φ(x, t) − I/2 (here φ(x, t) =(cid:82)

p ΨB[T]ppppdp as ΨB is determined by assuming S
and D are co-aligned in D’s principal coordinates. We refer to the above

model as the BQ-tensor model ([19, 20]). In the following, we deﬁne the 2D

p Ψdp is
the concentration) whose maximal non- negative eigenvalue λmax of Q
satisfying 0 ≤ λmax ≤ 1/2. We call its associated unit eigenvector the nematic
director m, and 0 ≤ s = 2λmax ≤ 1 the scalar order parameter.

1.4 Results and Discussion

When simulating the conﬁned active ﬂuid motion, we solve the D dynamics

equation (2.2) together with the incompressible Navier-Stokes equation at small
Reynolds numbers (Re = 10−3) using a mixed ﬁnite element method ([35, 36,
19]). At the inner (r = R1) and outer (r = R2) boundaries, a no-slip condition
u = 0 is imposed for the velocity ﬁeld, and a no-ﬂux condition (n · ∇)D = 0
is used for the second moment tensor D. Note that in this way, we do not

enforce any particular alignment direction at the boundary but guarantees the

global particle number conservation. For all the simulation results below, we

vary α and R1,2, and ﬁx dT = dR = 0.025 (for both dilute and concentrated

suspensions), as well as ζ = 0.5, β = 0.874 (for concentrated suspensions).

6

Figure 1.1:
(a-c) Characteristic ﬂow patterns (streamline overlaid on the
colormap of vorticity Ω) and (d-f) nematic structures (nematic director m ﬁeld
overlaid on the colormap of scalar order parameter s) when choosing R1 = 1.0
and R2 = 2.0. Insets in (a) and (d) are the comparisons between the analytical
(solid lines) and numerical (open circles) results for normalized uθ and Drθ
components taken along the white lines when |α| slightly goes beyond the
critical value |α|cr = 1.23 for the Isotropy-Circulation instability. (g) Phase
diagram of (|α|, R2 − R1) when ﬁxing R1 = 1.0. The solid and dashed lines
characterize the Isotropy-Circulation instability in an annulus and a straight
channel, respectively.The dash-dot line characterizes the Circulation-Traveling
wave instability in a straight channel for sharply-aligned rods by neglecting the
rotational diffusion.

1.4.1 Dilute Suspension

We ﬁrst examine the cases of a conﬁned dilute suspension by neglecting both

the steric interaction stress (ζ = 0) and the constraining stress (β = 0), and

start simulations from an initially near-isotropic state. As shown in ﬁgure

2.1(a-f) where R1 = 1.0 and R2 = 2.0 are ﬁxed, characteristic nematic
structures and ﬂow patterns are highlighted at different regimes of |α|. As
shown in panel(a), a unidirectional circulating ﬂow ﬁrst appears when |α| goes
beyond a certain critical value αcr (in this case αcr = −1.23). While the
nematic order is still low as shown in panel(d), it clearly shows that particles

7

become aligned in the azimuthal direction due to ﬂuid shear ([26, 37]), with a
stronger alignment near rigid walls. As |α| further increases up to 10.0 in
panel(b), a circulating ﬂow pattern still dominates but the streamlines exhibit

periodic bending deformations in the radial direction to form traveling waves,

leading to counterclockwise (CCW) and clockwise (CW) vortices near the

inner and outer boundaries, respectively. Such oscillatory patterns have also

been seen in other active liquid crystal systems owing to the bending (or splay)

deformations ([38, 39, 40]). In panel(e), the ﬂow-induced alignment is further

enhanced to form low-order structures (blue color). At even higher values of
|α|, the ﬂow pattern in panel(c) becomes seemingly chaotic. The resultant
phase
nematic
([7, 41, 42, 27, 5]) in which ±1/2 disclination defects stream around, and are
constantly generated/annihilated during interactions (see panel(f)).

liquid-crystalline

typical

active

ﬁeld

exhibits

a

Note that while our model is apolar, the observed phenomena are very

similar to those reported for polar active suspensions where effects of

self-swimming motions of microparticles (e.g., bacterium) are taken into

account ([43, 25, 17]). It is essential to recognize that the nonlinear dynamics

is governed by a concatenation of hydrodynamic instabilities as α increases,

and can be characterized by a series of instabilities including Isotropy to

Circulation, Circulation to Traveling wave, and Traveling wave to Chaos. The

overview of instabilities is highlighted in a phase diagram plotted in the
(|α|, R2 − R1) plane in panel(g) where R1 is ﬁxed to be 1.0. Furthermore, the
marginal stability curve that delineates the Isotropy-Circulation instability can

be calculated analytically in the polar coordinates. To do this, we perturb the

8

Figure 1.2: (a) Growth rate Re (σ) for the long-wave (diamonds) and ﬁnite-
wavelength (lines) instabilities as a function of wavenumber k for Isotropy-
Circulation instabilities in a periodic straight channel when choosing h = 0.5,
dT = dR = 0.025.
(b) Streamlines overlaid on the colormap of u velocity
component at α = −4.0. (c) Velocity proﬁles at α = −2.0 (umax = 0.07) and
α = −4.0 (umax = 0.13), respectively.
near-isotropic solutions as D = I/2 + D(cid:48) (r) and u = u(cid:48) (r) ( (cid:28) 1) by
assuming azimuthal symmetry. Following the procedure by [9], we examine
∂t = κD(cid:48) with the growing factor κ ≥ 0.
an initially exponential growth ∂D(cid:48)
After some algebra, we ﬁnd D(cid:48)
θθ = 0, and the hydrodynamic instability
arises from the coupled linearized equations for u(cid:48)

rr = D(cid:48)

θ and D(cid:48)

rθ below:

rθ

.

4dT

(cid:17)

,

4dT

(cid:16)

=

2r

θ

rθ

r2 ∂2D(cid:48)

rθ

(1.6)

(1.7)

r + αD(cid:48)

θ

We are able to seek the solutions

rθ = A0
2r2 ,
rθ = − A0

∂r +(cid:0)γr2 − 4(cid:1) D(cid:48)
(cid:0)√
γr(cid:1) + 2A2√
(cid:0)√

∂r − u(cid:48)
∂u(cid:48)
∂r2 + r ∂D(cid:48)
where γ = − 16dR+4κ+α
D(cid:48)
θ = − A0
u(cid:48)
Nj are respectively the Bessel functions of the ﬁrst and the second kind, with
coefﬁcients Aj to be determined. By applying boundary conditions ∂D(cid:48)
∂r = 0
and u(cid:48)
θ = 0 on the two walls, we then solve the eigenvalues of γ numerically at
κ = 0 (solid line in panel(g)), as well as the velocity and D ﬁeld (insets in

γr(cid:1) + A3r, where Jj and

− A0
4γdT r2
+ 2A1√

√
∂N2(
∂r

1 + α
4γdT

√
∂J2(
∂r

γr)

γr)

and

+ A1

+ A2

γ J1

γ N1

rθ

panels(a,d)), which, again, conﬁrm our simulation results.

Next, we reveal the nature of hydrodynamic instabilities and gain physical

9

Figure 1.3: (a) Growth rate Re (σ) as a function of wavenumber k for bending
instabilities of the sharply-aligned Extensors when conﬁned in a periodic
straight channel when choosing h = 0.5, dR = 0 and dT = 0.025. The
maximum growth rates at the critical wavenumbers kcr are marked by solid
circles. (b) and (c) show the vorticity ﬁeld Ω as bending instability start to grow
at α = −4.0 and α = −8.0, respectively. The solid black scales represent the
length scales π/kcr that match the separation distance between the neighboring
vortices.

insights into formation of the observed coherent structures in ﬁgure 2.1. As

shown

below, while

analytical

approaches

are

feasible

for

the

Isotropy-Circulation and Circulation-Traveling wave transitions, the Traveling

wave-Chaos instability is highly nonlinear and hence can only be explored

numerically. Here we avoid tedious mathematical manipulations in the polar

coordinates to simplify our analysis by considering conﬁnement in a straight
periodic channel (i.e., in the limit R1,2 → ∞) where very similar collective
dynamics and hydrodynamic instabilities are observed ([17]). For the two

cases considered, we employ a homogeneous base-state solution of D0

without background ﬂow (i.e., u0 = 0, p0 = 0), and introduce perturbations

10

(D(cid:48), u(cid:48), p(cid:48)) to yield the following linearized equations:

∂t − (∇u(cid:48)) · D0 − D0 ·(cid:16)∇u(cid:48)T(cid:17)

∂D(cid:48)

+ 2E(cid:48) : S0

= dT∇2D(cid:48) − 4dRD(cid:48),
∇p(cid:48) − ∆u(cid:48) = ∇ · (αD(cid:48)) ,

∇ · u(cid:48) = 0,

(1.8)

(1.9)

(1.10)

∂y = 0 and u(cid:48) = 0 on the upper (y = h) and
with boundary conditions ∂D(cid:48)
bottom (y = −h) walls. (Note here D(cid:48)
22 = 0.) We then take Fourier
transform in the x direction of all perturbed quantities f(cid:48) such that f(cid:48) (x, y) =
˜f (k, y) exp (ikx + σt), where ˜f is the amplitude function, k is the wavenumber,

11 + D(cid:48)

and σ is the complex-valued growth rate.

Case 1: Isotropy-Circulation Instability. In this case the instability grows

from an initially isotropic state, i.e., D0 = I/2. The linearized equations in

(2.3) and (2.4) become

(cid:0)δ(2) − k2(cid:1)2

˜v + 2αk2δ ˜D11 − ikα(cid:0)δ(2) + k2(cid:1) ˜D12 = 0,
(cid:0)δ(2) − ω(cid:1) ˜D11 − δ˜v = 0,

4dT k(cid:0)δ(2) − ω(cid:1) ˜D12 + i(cid:0)δ(2) + k2(cid:1) ˜v = 0,

2dT

where δ(n) = dn

dyn represents the nth derivative of y, and ω = σ+4dR

dT

(1.11)

(1.12)

(1.13)

+ k2. We

ﬁnd that the general solution of the above equations has the vector form

( ˜D11, ˜D12, ˜v) = (B1F1 + B2F2)eky + (B3F3 + B4F4)e−ky

+ B5F5e

√

ωy + B6F6e−√

ωy + B7F7ei

√−χy + B8F8e−i

(1.14)
√−χy,

11

where χ = ω + α
4dT

, and Fj can be written as:

F2 =(cid:8)ik2y, k2y, 2idT
F4 =(cid:8)−ik2y, k2y, 2idT

(cid:2)k(k2 − ω)y + k2 + ω(cid:3)(cid:9) ,
(cid:2)k(k2 − ω)y − k2 − ω(cid:3)(cid:9) ,

F1 =(cid:8)ik, k, 2idT (k2 − ω)(cid:9) ,
F3 =(cid:8)−ik, k, 2idT (k2 − ω)(cid:9) ,
F5 =(cid:8)i(cid:0)k2 + ω(cid:1) , 2
ωk, 0(cid:9) ,
F6 =(cid:8)−i(cid:0)k2 + ω(cid:1) , 2
ωk, 0(cid:9) ,
√−χ,−χ − k2,−iαk(cid:9) ,
F7 =(cid:8)2k
√−χ,−χ − k2,−iαk(cid:9) .
F8 =(cid:8)−2k

√
√

By implementing the eight boundary conditions δ ˜D11 = δ ˜D12 = δ˜v = ˜v = 0 at
y = ±h, we obtain a linear system for Bj. The existence of non-trivial solutions
gives ω and hence σ. On the other hand, we ﬁnd that the long-wave solution has

to be solved separately due to singularity at k = 0. By making use of the fact

that the solution is unidirectional, we are able obtain the equilibrium solution:

(cid:16)(cid:113) α+4σ+16dR

(cid:17)(cid:17)

y

4dT

, which yields

(1.15)

(cid:16)

(cid:16) ˜D11, ˜D12, ˜v
(cid:17)

the growth rate

=

0, sin

y

, cos

(cid:17)

(cid:16)(cid:113) α+4σ+16dR
(cid:17)2
−(cid:16) π

4dT

2h

σ = −α
4

dT − 4dR, at k = 0.

As shown by the comparisons in ﬁgure 1.2(a), we ﬁnd the real part of the

growth rates of the long-wave modes (diamond symbols) are always larger

than those of the corresponding ﬁnite-wavelength ones (solid lines), and hence

condition: |α| > 4(cid:0) π

dominate instability. Moreover, equation (1.15) suggests a marginal stability

(cid:1)2 dT + 16dR, which recovers the theoretical predictions

2h

for unbounded suspensions in the limit h → ∞ ([27]). When making

12

connections to the annulus geometry by replacing the gap width 2h by
R2 − R1, we ﬁnd the marginal instability curve in ﬁgure 2.1(g) can be
−2 dT + 16dR but with
accurately ﬁtted as a function of the form K (R2 − R1)
a different constant K due to curved geometries (in this case K ≈ 28.7). It is
apparent that the prediction in the straight periodic channel (dashed line in

ﬁgure 2.1(g)) agrees well with that of the curved annulus, suggesting that

equation (1.15) in fact provides a reasonable estimate of the time scale of the

Isotropy-Circulation instability under a general parallel conﬁnement.

Case 2: Circulation-Traveling wave Instability. It needs to keep in mind

that the Circulation-Traveling wave instability is in fact a secondary instability

developed from a shear-induced aligned state whose analytical form, however,

is lacking. Nevertheless, we consider a reduced model where all the Extensor

particles are initially aligned in the x direction, and neglect the rotational

diffusion by choosing dR = 0 ([19]).

For such a “sharply-aligned”

the

homogeneous

conﬁguration,
simply
D0 = diag {1, 0}. After some algebra, we are able to eliminate several
unknown variables, and show that the hydrodynamic instability is determined

base-state

solution

is

by the following six-order ODE for the off-diagonal component ˜D12:

(cid:19)

δ(6) ˜D12 − k2

(cid:18) σ
(cid:18)2σ − α

dT k2 + 3

(cid:19)

dT k2 + 3

δ(4) ˜D12
δ(2) ˜D12 − k6

+ k4

(cid:19)

˜D12 = 0

(1.16)

(cid:18)σ + α
6(cid:80)

dT k2 + 1

which admits the general solution ˜D12 =

Cj exp (λjy). The eigenvalues λj

and hence the growth rate σ are then solved numerically by applying the

converted boundary conditions (cid:0)σ + dT k2(cid:1) ˜D12 − dT δ(2) ˜D12 = 0 and

j=1

13

δ ˜D12 = δ(3) ˜D12 = 0 at y = ±h, and seeking for non-trivial solutions. The real
part of the growth rates as a function of k are plotted in ﬁgure 1.3(a) at

different values of α, showing the maximum growth rates occurring at ﬁnite

critical wavenumbers kcr > 0.

In panels(b,c), we set up the corresponding

numerical simulations with the same initial conditions, and highlight the

snapshots of the vorticity ﬁeld Ω = ∂u

∂x. The simulations results show the
generation of ﬂuid jets that are perpendicular to the particle alignment

∂y − ∂v

direction due to the nematic bending instability, which leads to uniformly

spaced vortices of alternating signs. The separation distance between the

neighboring vortex pair can be accurately measured by the characteristic

length scale π/kcr.

While the above analysis is for a reduced model by assuming particle

alignment

in the x direction without background ﬂow,

it

reveals the

ﬁnite-wavelength nature of the bending instability, and provides quantitative

measurements for the intrinsic time and length scales selected by the parallel

conﬁnement. Similar to Case 1, the predicted Circulation-Traveling wave

borderline (dash-dot line in panel(g)) in the straight channel indeed agrees

with the simulation results of the annulus geometry, especially in the regime of
high |α| and small R2 − R1 where particles are strongly aligned and hence
close to the sharp alignment approximation. Beyond the traveling-wave

regime, the transition towards chaotic ﬂows is highly non-linear, and can only

be determined by numerical simulations.

14

1.4.2 Concentrated Suspension

In this

section we examine the nonlinear dynamics of concentrated

suspensions where both β and ζ are non-zero. A mean-ﬁeld torque is

introduced to govern the rotational dynamics of Extensor particles in the

kinetic model through a Maier-Saupe steric potential ([31]). As shown by

[27, 44, 19], for rod-like particles, the resultant enhanced steric interactions

spontaneously drive the systems away from an initially isotropic state to form

a nematically-aligned state when ζ > 4dR.

Compared to the dilute cases where isotropy is the only admissible

equilibrium state, we have obtained the steady-state solutions of D without
ﬂows at relatively small values of |α|, which exhibit much more complicated
nematic conﬁgurations as shown in ﬁgure 1.4. For passive particles when

choosing α = 0, ﬁgure 1.4(a) reveals that the system quickly relaxes to a

homogeneous nematic state where the equilibrium solution satisﬁes a 2D

distribution of Bingham form ([27, 19, 20])

(cid:82) 2π

ΨB(φ) =

exp [δ (ξ) cos (2φ)]

0 dφ(cid:48) exp [δ (ξ) cos (2φ(cid:48))]

(1.17)

where φ is the rod’s orientation angle, and δ is a function of coefﬁcient

ξ = 2ζ/dR. Therefore, the degree of alignment is in fact determined by the

ratio between ζ and dR. As α increases, it is seen that the nematic ﬁeld lines in

ﬁgure 1.4(b) become distorted to spiral outwards to the annulus boundaries,

reminiscent of the defect structures of active nematics when conﬁned in a
circular disk ([9, 19]). As |α| approximately goes beyond 0.5, the nematic
structures in ﬁgure 1.4(c,d) switch to become azimuthal-symmetric and

independent of the values of α, which can be solved analytically by setting the

15

Figure 1.4: Equilibrium solutions of nematically-aligned states at different
values of α when choosing R1 = 1.0 and R2 = 2.0, nematic director ﬁeld
overlaid on the colormap of s. Inset in panel(d): Comparisons between the
analytical (solid line) and numerical (open symbols) results at α = −0.5,−2.0
for Drr component. The numerical data are taken along the white line marked
in panel(c).

right-hand-side of equation (1.5) to be zero, and then seeking for the

steady-state solution of D as a function of r only. The observed variations in

the equilibrium nematic states can be illustrated by examining the system’s

transient dynamics that evolves from near-isotropy. One interesting example is
shown in ﬁgure 1.5 when choosing α = −2.0. An Isotropy-Circulation
instability ﬁrst occurs in ﬁgure 1.5(a) to reorient particles, due to the combined

effect of the circulating ﬂow and the Maier-Saupe steric interactions. However,

in this case the conﬁnement effect in the thin gap is so strong that it suppresses

the bending deformation occurring in the radial direction. As a result, the

induced circulating ﬂow gradually diminishes, and the system eventually

relaxes to an axisymmetric equilibrium state where all particles are

approximately aligned in the azimuthal direction; see insets in ﬁgure 1.5(b-d).

Increasing gap width relaxes the conﬁnement effect to facilitate the active

ﬂow generation through bending instability, which can be similarly explained

by the analysis in the dilute cases above. As shown in ﬁgure 1.6, system’s

16

Figure 1.5: Evolution of the ﬂow patterns for a concentrated Extensor
suspension when choosing α = −2.0, R1 = 1.0 and R2 = 2.0. Insets: nematic
director ﬁeld overlaid on the colormap of s.

Figure 1.6: Flow patterns (streamline overlaid on the colormap of Ω, upper
row) and nematic structures (nematic director ﬁeld overlaid on the colormap of
s, lower row) that characterize the long-time dynamics of a conﬁned Extensor
suspension when choosing α = −2.0, R1 = 0.75 and R2 = 2.0.

17

Figure 1.7: Short-time evolution of the ﬂow (velocity vector overlaid the
colormap of on Ω, upper row) and nematic (nematic director overlaid on the
colormap of s, lower row) ﬁelds for the case shown in ﬁgure 1.6.
In (e-h),
the black arrows shows the direction of the crack formation and +1/2 defect
movement.

long-time dynamics clearly exhibits an Isotropy-Circulation (panels(b,f)) and

then a Circulation-Traveling wave (panels(c,g)) instability. Compared to the

stable structures observed in the dilute cases, structure formation and

evolution in concentrated suspensions are much more complicated.

In

panel(c), it is shown that active ﬂows bend the streamlines to form traveling

waves, leading to six evenly-spaced “incipient cracks” growing from the inner

wall into the bulk regime shown in panel(g). Disclination defects of charge
±1/2 are seen to be born along these cracks, move around, and undergo
nucleation/annihilation when interacting with each other as well as with rigid

boundaries. Interestingly, the system gradually approaches a quasi steady-state

where four +1/2 defects oscillate around their equilibrium positions; see

panel(h). The ﬂow ﬁeld remains to be traveling-wave like but also periodically

18

switches directions; see panel(d). The reader is referred to the supplementary

movie S1 for the entire history of the complex dynamics.

As shown in ﬁgure 1.7, a close look at

the short-time dynamics of

oscillating +1/2 defects and the accompanying swirling ﬂow reveals the

subtle interplay among the nematic structure variation,

the active ﬂow

generation, as well as their interactions with the curved walls. The swirling

ﬂow drives the upward movement of the +1/2 disclination defect, and in the

mean time,

induces an opposite bending deformation nearby to form an

incipient crack. As shown in panels(f,g), the crack “head” keeps moving

downward to become a +1/2 defect (see panel(h)); while its “tail” gradually

merges with the up-moving +1/2 defect before it hits the wall. Corresponding

to the nematic structure change, the resultant CCW (CW) vorticity is seen to

be strengthened (weakened) when the crack and defect interact, and then

weakened (strengthened) during their annihilation. We ﬁnd such an oscillating
state is very robust when choosing R1 less than 1.5, and the gap width R2 − R1
between 1.0 and 1.5 (also see movie S2 for the similar dynamics observed in a

smaller annulus). As shown in movie S3, more complex defect dynamics and
seemingly chaotic ﬂows dominate when R2 − R1 > 2.0, resembling the active
nematic ﬂows observed in an unbounded domain ([19, 27]).

1.5 Conclusion

In this work we have studied the nonlinear dynamics and ﬂow patterns of

apolar Extensor suspensions conﬁned in an annulus geometry through

modeling and simulation. The BQ-tensor model with the Bingham closure,

19

which was coarse-grained from a mesoscale kinetic model, is used to describe

the evolution of the second-moment D tensor ([19]). The D dynamics

equation is coupled with the incompressible Navier-Stokes equation which are

solved by using a Galerkin mixed element method at low Reynolds numbers.

We have observed active ﬂow generation and spontaneous coherent structures

in both dilute and concentrated suspensions. We have gained physical insights

of the analytical structure, hydrodynamic instabilities, and characteristic time

and length scales under parallel conﬁnement by performing stability analyses

for dilute suspensions in a periodic straight channel.

In particular, we have

shown analytically the long-wave and ﬁnite-wavelength nature of

the

Isotropy-Circulation and Circulation-Traveling wave instabilities, respectively.

By using similar models and analyses, we will be able to explore more

challenging situations when active ﬂuids are conﬁned in complex geometries

at high dimensions.

20

2

TRANSPORT OF AN APOLAR ACTIVE SUSPENSION IN A

RECTANGULAR CHANNEL

2.1 Abstract

We study ﬂuid and mass transport in a dilute apolar active suspension conﬁned

in a rectangular channel. We adopt a coarse-grained active liquid crystal

model to describe collective unstable dynamics of non-motile but mobile

rod-like particles, i.e., the “Extensors”, in Stokes ﬂows. By using a Galerkin

mixed ﬁnite element method, we are able to reveal various patterns of

spontaneous coherent ﬂows that can be unidirectional, traveling-wave, and

chaotic. Comparing with our previous work by [1], we extend our analysis in

weak ﬂows to strong ﬂow regimes by solving steady unidirectional ﬂows of

ﬁnite strengths analytically, and then perform more rigorous analyses to

precisely capture the marginal stability curve of the unidirectional ﬂow to

traveling-wave instability. More interestingly, we ﬁnd that the concatenation

of instabilities may produce unique high-order unstable modes selected by

geometric conﬁnement, which can also be used to predict all the different

kinds of instabilities shown in the phase diagram.

In addition, we study

transport of a circular cylinder in active ﬂows, and have identiﬁed intriguing

trajectories by varying the cylinder size, channel width, and Extensor particle

activity. The weakening effects of active ﬂow are observed when a large

cylinder blocks the channel or when multiple cylinders are transported at a

21

high volume fraction.

2.2 Introduction

Manipulating small volumes of ﬂuids is critical in a wide range of biological

and engineering processes such as mixing, separation, targeted drug delivery,

and micropumps.

In a microﬂuidic environment that is typically within

sub-millimeter, accomplishing these functions requires us to implement

efﬁcient ﬂuid actuation methods via external controls from pressure or

electric/magnetic/optical ﬁelds, in order to overcome ﬂuid resistance. Recent

advances of active matter, a novel class of non-equilibrium systems

comprising self-driven constituents that consume a local fuel (e.g., chemical

fuel) to generate particle motions or induce ﬂows ([3, 4]), can potentially serve

as new “engines” for precisely manipulating ﬂuid in engineering applications.

In “wet” active matter systems, suspended motile microparticles effectively

exert stresses upon the ambient liquid, which itself acts as a coupling medium

for complex collective dynamics with ordering transition, ﬂuctuating density

and force/stress generation. Examples include biological swimmers such as

bacteria and microalgae ([45]),

synthetic colloidal surfers actuated by

chemical reactions ([46]), as well as puriﬁed biopolymers (e.g., microtubules

(MTs)) mixed with molecular motors (e.g. kinesin) [7, 27]. More interestingly,

in concentrated active suspensions, the enhanced particle interactions may

lead to a myriad of dynamic states such as active crystalline phases that are

featured by spontaneous coherent structures such as streaming topological

defects and coherent ﬂows ([18]).

22

Figure 2.1: Typical (quasi-)steady coherent ﬂows (a,c,e) and nematic orders
(b,d,f) for unidirectional ﬂow (a,b),
traveling wave (c,d), and chaos (e,f).
Simulations are performed when choosing the half channel width H = 0.5,
dR = dT = 0.025, and α = −2.0,−8.0,−20.0. Comparison of simulation
and analytical solutions of the velocity proﬁle is shown in inset of (a) where
the maximum velocity is umax = 0.57. (g) Phase diagram of (|α|, 2H). The
solid line are the marginal stability curves from the analyses of the Q-U and
U-T instabilities. The dashed lines are the predictions using the geometrically-
selected high-order modes. The dotted line is our previous prediction when
using a sharp-alignment assumption ([1]).

As we glean a higher level understanding of active matter, it is intriguing to

exploit their anomalous physical properties for doing useful mechanical work.

Indeed, recently several microﬂuidic experiments have shown that active

matter, when appropriately guided by boundaries, can be used to pump/mix

ﬂuids ([12, 25]), transport cargos ([47]), and power microgears ([48]). These

intriguing experimental observations have given rise to new,

as yet

unanswered questions: How will microparticle propulsion and local

interactions change the near-wall dynamics, and hence bulk coherent ﬂows?

Will active ﬂows selectively deliver cargos of particular sizes and shapes? To a

large extent,

the fundamental understandings of how collective dynamics

interact with the complex microﬂuidic environment will shed light upon novel

ﬂuid and mass transport mechanisms promised by active matter.

23

Here we study the transport mechanisms of a dilute active suspension

conﬁned by two parallel plates. We employ a reduced model for non-motile

but mobile rod-like particles that are only advected by ﬂuid ﬂow. These

particles elongate or stretch to exert extensile dipolar stresses on the solvent,

and interact with each other through hydrodynamic coupling and steric

alignment torques, which is referred to as the “Extensors”. As shown by [19],

from a mesoscale kinetic model, we are able to derive a macroscale active

liquid crystal model,

i.e., “BQ”-tensor

theory,

to describe Extensors’

ensemble dynamics by solving the evolution equation of the second moment

tensor. Using this model, we have successfully revealed unstable collective

dynamics and pattern formation of Extensor suspensions conﬁned in an

annulus geometry ([1]).

Compared with our previous work ([1]) where analyses were performed

under some ad-hoc assumptions (e.g., being isotropic or sharply-aligned) in

weak ﬂow regimes, here we make analytical efforts in strong active ﬂows by

solving the analytical solutions of unidirectional ﬂows with ﬁnite strengths.

Therefore, we are able to capture the unidirectional ﬂow to traveling wave

instability accurately. More interestingly, we show that besides the

leading-order modes that characterize the main ﬂow features, the existence of

the geometrically-selected high-order modes at

the current ﬂow regime

provides accurate predictions of the ﬂow transition at the next level. Next, we

investigate transport behaviors of circular cylinders in apolar active ﬂows, and

ﬁnd that their dynamics, when either ﬁxed or free-moving, can be analytically

solved in the weak ﬂow regime.

Furthermore, we perform numerical

24

simulations in strong ﬂow regimes to capture intriguing trajectories in

traveling wave ﬂows for cylinders of different sizes, as well as weakened

active ﬂows for transport of a large cylinder or multiple cylinders at a high

volume fraction.

2.3 Mathematical model

Here we brieﬂy describe the mathematical model for the ensemble dynamics

of a collection of Extensor particles suspended in a Newtonian ﬂuid. More

details can be found in our previous publications ([19, 20, 1]). In this model,

we assume rod-like particles are non-motile but mobile. They exert dipolar

stresses upon the ambient liquid through nearly symmetric elongation or

stretching.

The resultant active ﬂuid is governed by the following

(dimensionless) incompressible Stokes equation driven by an extra particle

stress:

∇ · u = 0, ∇p − ∆u = ∇ · (αD)

(2.1)

where u and p are the ﬂuid velocity and pressure, respectively. For dilute

suspensions,

the evolution equation of D tensor follows the evolution

(cid:18)

equation:

∂D
∂t

+ u·∇D− (∇u)· D− D· (∇u)T + 2E : S− dT ∆D + 4dR

(cid:19)
where E = 1/2(cid:0)∇u + ∇uT(cid:1) is the rate-of- strain tensor, dR and dT are

D − I
2

= 0,

(2.2)

the rotational and translational diffusion coefﬁcients, respectively. The right-

hand-side of equation (2.1) is attributed to a dipolar active stress αD whose

strength coefﬁcient α (α < 0) characterizes activity of microswimmer. Here

25

D and S are the second- and fourth-rank symmetric tensors respectively, which

are the corresponding second- and fourth-order orientation moments in terms
of the particle orientation vector p (|p| = 1), and can are deﬁned by taking
the average of p on the surface of a unit sphere via a probability distribution

function Ψ(x, p, t) such that D = (cid:82)

p ppΨdp and S = (cid:82)

p ppppΨdp ([29]).
The model is referred to as the BQ-tensor model ([19, 1]) since a 2D tensor-
order parameter is deﬁned as Q(x, t) = D− I/2 with its maximal non-negative
eigenvalue 0 ≤ λmax ≤ 1/2. We call its associated unit eigenvector the nematic
director m, and 0 ≤ s = 2λmax ≤ 1 the scalar order parameter.

In nondimensionalization, we assume that there are N Extensor particles of
length l and width b (b/l (cid:28) 1) distributed in a volume V , and hence introduce
an effective volume fraction ν = nl3 where n = N/V is the mean number
density ([49]). Then we choose the length scale lc = (nl2)−1, the velocity
scale |u0| which represents the surface velocity due to elongation or stretching
motion, and the stress scale µ|u0|/lc. Equation (2.2) can be closed by
expressing the fourth-moment tensor S in terms of D through the so-called

Bingham closure ([33, 19]).

In numerical simulation, we solve the D dynamics equation (2.2) together

with the incompressible Navier-Stokes equation at small Reynolds numbers
(10−3) using a mixed ﬁnite element method ([19, 1]). For all the simulation
results below, we ﬁx dT = dR = 0.025. To resolve the conﬁned geometries as
well as stationary and moving obstacles, a no-ﬂux condition (n · ∇)D = 0 is
used for the second moment tensor D on geometric boundaries. To solve the

velocity ﬁeld, we imposed a no-slip condition (u = 0) on the top and bottom

26

denoted

by

xc,

a

walls, i.e., y = ±H. The freely-moving cargo, whose center-of-mass position
is
that
u = U + Ω × (x − xc) is imposed on the surface, where U and ω are
respectively the translational and rotational speed, and are determined by the

rigid-body motion

follows

such

force and torque free conditions. An Arbitrary Lagrangian-Eulerian ﬁnite

element method with a moving-mesh technique is used to simulate obstacle

motions in conﬁned active ﬂows ([36, 35]).

2.4 Analysis and simulation results

2.4.1 Spontaneous Coherent Flows

Generally speaking,

the unstable dynamics of conﬁned active ﬂows are

determined by a delicate competition between the (destabilizing) active

nematic force due to aligned structures against the (stabilizing) viscous force
and topological constraints of channel walls. As |α| goes beyond a certain
critical value, the system deviates from an initially isotropic state to form

spontaneous coherent ﬂows as the pattern changes from unidirectional ﬂows to

traveling waves, and eventually turns into chaos.

As shown in Figure 2.1, we have captured three transitions, namely

quiescence to unidirectional ﬂow (Q-U), unidirectional ﬂow to traveling wave
(U-T), and traveling wave to chaos (T-C), by systematically varying |α| and
channel width. As shown in panel(a), in the weak ﬂow regimes of small |α|,
when ﬁxing the channel width, an Q-U instability ﬁrst occurs, during which

randomly-oriented microparticles become more aligned,

accompanying

development of a unidirectional ﬂow with an approximately parabolic velocity

27

proﬁle. As |α| increases, as shown in panel (b), the ﬂuid velocity ﬁeld
gradually loses symmetry as the streamlines bend in the transverse direction to

form rotating vortices, yielding wavy patterns with a clear separation of high
and lower order zones. Further increasing |α| eventually leads to chaotic ﬂows
and motile defects. The constrained collective dynamics resemble those of

Extensor particles trapped in a annulus, or Pusher particles in race-track

geometries ([25, 17, 1]).

2.4.2 Concatenation of Hydrodynamic Instabilities

Next, we go beyond the analyses of weak ﬂows ([1])

to investigate

concatenation of hydrodynamic instabilities that occur in both weak and

strong ﬂow regimes, aiming to obtain a thorough physical understanding of the

observed ﬂow transitions in the phase diagram in ﬁgure 2.1(g). We ﬁrst derive

the analytical solutions of steady-state unidirectional ﬂows with ﬁnite

strengths. To begin with, we linearize the governing equations (2.1)-(2.2) by

introducing the base-state (denoted by subscript “0”) and disturbance solutions
(denoted by superscript “,”) such that D = D0 + D(cid:48) and u = u0 + u(cid:48)
( (cid:28) 1), which lead to the following coupled equations:

∂t + u0 · ∇D(cid:48) + u(cid:48) · ∇D0 − ∇u0 · D(cid:48) − D(cid:48) · ∇(cid:0)u0(cid:1)T

∂D(cid:48)

−∇u(cid:48) · D0 − D0 · ∇u(cid:48)T

+2E0 : S(cid:48) + 2E(cid:48) : S0 = dT ∆D(cid:48) − 4dRD(cid:48),
∇ · u(cid:48) = 0, ∇p(cid:48) − ∆u(cid:48) = ∇ · (αD(cid:48)),

(2.3)

(2.4)

with no-ﬂux and no-slip boundary conditions imposed on the top and bottom
plates. We then take Fourier transform in x−direction of all disturbance

28

quantities f(cid:48) such that f(cid:48) (x, y) = ˜f (k, y) exp (ikx + σt), where ˜f is the
amplitude, k is the wavenumber, and σ is the complex-valued growth rate.

Quiescence-Unidirectional Flow Instability. Consider initially rods are

randomly oriented, corresponding to an isotropic state without ﬂows. The

base-state solutions are simply D0 = I/2 and u0 = 0 (hence E0 = 0). From

an microscopic perspective,

this state can be described by a constant

distribution function Ψ0 = 1/2π which yield the fourth-order moment

8 (δijδkl + δikδjl + δilδjk). As shown in our previous work by [1], the
ijkl = 1
S0
system ﬁrst undergoes an Q-U instability. Due to the singularity at k = 0, we

solved the long-wavelength mode separately to obtain the growth rate
Re (σ) = − α

4 −(cid:0) π
(cid:1)2 dT − 4dR (here Re denotes the real part).
2H y(cid:1) (Φ0 is a unknown coefﬁcient), as well as the other
2α sin(cid:0) π

We were also able to derive analytically both the long-wavelength

disturbance solution u(cid:48) = Φ0H cos(cid:0) π

2H y(cid:1), v(cid:48) = D(cid:48)

11 = −D(cid:48)

22 = 0 and

2H

D(cid:48)
12 = Φ0π

unstable modes which can be represented by a linear combination of eight

vector polynomials in Fourier space (see details in [1]). As shown in ﬁgure

2.2(a),

these high-order unstable modes grow more slowly and are of

ﬁnite-wavelength with the maximum growth rate occurring at k > 0.

Unidirectional Flow-Traveling Wave Instability. As |α| further increases,
the enhanced unidirectional ﬂows will realign Extensor microparticles in the

ﬂow direction. As we have previously proven ([27, 1]), extensile rod-rod

interactions with a nematically-aligned conﬁguration are hydrodynamically

unstable, and will drive a shear-induced secondary instability that grows

perpendicular to the alignment direction to bend streamlines in the transverse

29

directions and create traveling-wave patterns,

leading to counterclockwise

(CCW) and clockwise (CW) vortices near the walls.

In the strong ﬂow regime, while the exact unidirectional steady-state

solutions are unfeasible, analytical manipulation is still possible if good

approximations can be made for the rods of aligned conﬁgurations. To do this,
we expand Ψ near isotropy such that Ψ = Ψ0 + Ψ(cid:48) = A0 + A1 : (D(cid:48)) where
the two constants A0 = 1/2π and A1 = 2
2) can be determined by
p Ψdp ) and the
second-order moment D ([20]). At the leading order, the fourth-moment

comparing with the zeroth (i.e., concentration φ = (cid:82)

π (pp − I

tensor

S0

ijkl

=

is

then

calculated

as

p pipjpkplΨ0dp

=

− 1
24 (δijδkl + δikδjl + δilδjk) +

S

(cid:82)

1

6 (δijDkl + δikDjl + δilDjk + δjkDil + δjlDik + δklDij).

After substituting S0 into equation (2.2) and eliminating x−dependency, we

can simplify the governing equations at steady state as

δ(2)D0

11 = 4

δ(2)D0

12 =

(cid:19)

(cid:18)
(cid:18)16dR + 3α

11 − 1
D0
2

dR
dT

4dT

(cid:1)2 ,

(cid:0)D0
(cid:19)

12

+

α
dT

D0
11

− α
D0
12,
dT
δ(1)u0 = −αD0
12,

(2.5)

(2.6)

(2.7)

where δ(n) = dn

dyn represents the nth derivative of y. As suggested by the
analytical disturbance solutions, we ﬁnd it is reasonable to assume the nematic

solutions

22 = B0 + B1 cos(cid:0) π

ﬁeld
11 = 1 − D0
D0
unknown coefﬁcients Bi. Then we can solve the velocity components as
u0 = − 2HαB2

2H y) and v0 = 0. By using the expressions of D0 and u0

12 = −B2 sin(cid:0) π

2H y(cid:1), with the

H y(cid:1) and D0

trigonometric

cos( π

form

to

be

of

π

30

and applying boundary conditions at the two walls, we are able to obtain Bi

(cid:113)− 8dRγ1γ2

α2γ3

B2 =

analytically after some algebra as B0 = 1

, B1 = − dRγ2
H 2 dT + α + 16dR, and
H 2 dT + 6dR (hence B1 (cid:28) B0). As shown in the inset of ﬁgure 2.1(a) for
the velocity proﬁle, this solution agrees well with the simulation results. In

2 + γ1γ2
αγ3
4H 2 dT + dR, γ2 = π2

, where γ1 = π2

γ3 = π2

, and

αγ3

addition, it can be easily veriﬁed that setting γ2 = 0 recovers the isotropic

solution, and yields the marginal stability curve of the Q-U instability in the

phase diagram.

Substituting the base-state solutions into (2.3) and (2.4) yields the linearized

equations in Fourier space as:

π

(cid:0)dT δ(2) − ω + 2iHkαC2
+(cid:0)ikC2 sin( π
−(cid:0) 1
(cid:0)dT δ(2) − ω + 2iHkαC2
+(cid:0) π
−(cid:0) i

2H C2 cos( π

2δ(1) + i

π

H sin( π

2H sin( π

2H y)(cid:1) ˜D11 + αC2
H y)(cid:1) ˜v
2H y)δ(2)(cid:1) ˜v = 0,
2H y)(cid:1) ˜D12 − αC2
4)δ(2)(cid:1) ˜v = 0,

2H sin( π
H y) + C0 − 1

cos( π

2H y) + πC1

k C2 sin( π

cos( π

2H y) + ik(C1 cos( π
H y) + C0 − 3

(cid:0)δ(4) − 2k2δ(2) + k4(cid:1) ˜v = −2k2αδ(1) ˜D11 + (αkiδ(2) + αk3i) ˜D12, (2.10)

k (C1 cos( π

2H y) ˜D12

2H y) ˜D11

4)(cid:1)

(2.8)

(2.9)

where ω = σ + 4dR + dT k2. By sweeping values of α and H, we numerically

solve the eigenvalue problem to obtain the growth rate σ for non-trivial
solutions. The Re (σ) − k curves in ﬁgure 2.2(b) show that the U-T instability
is dominated by unstable modes of a ﬁnite wavelength whose maximum

growth rate occurs at a critical wavenumber kcr > 0. Also as shown in ﬁgure

2.1(b), the predicted characteristic length π/kcr (solid blue line) agrees with

31

Figure 2.2: Unstable modes and growth rates Re (σ) of linear stability analyses
for Q-U (a), U-T (b), and T-C (c) instabilities. In (a,b), the channel width is
chosen as H = 0.5; in (c), α = −4.0 is ﬁxed. The open diamond symbols (in
panel(a) at k = 0) represent the long-wavelength (M0) mode associated with
the unidirectional ﬂow solution. The solid, dashed, and dotted lines in (a-c)
represent the 1st (M1), 2nd (M2), and 3rd (M3) unstable modes, corresponding
to the single, double, and triple vortex arrays, respectively. The colorﬁelds for
vorticity and streamlines exhibit ﬂow patterns, and are selected at the critical
wavenumbers (circular red dots).

the separation distance between the neighboring CW and CCW vortices from

simulation results. Comparing to our previous results of “sharply-aligned”

rods by assuming perfect alignment and neglecting the rotational diffusion

([1]),

the newly predicted U-T transition agrees much better with the

numerical results.

Geometrically-Selected Unstable Modes. A thorough understanding of the

phase diagram in ﬁgure 2.1 requires uncovering the transition mechanisms of

the T-C instability. Nevertheless, directly performing stability analysis is

unfeasible due to the lack of analytical traveling-wave solutions at the base

state. Besides the leading-order modes, we ﬁnd that

the selection and

development of high-order modes in different ﬂow transitions not only reﬂect

the effects of geometrical conﬁnement but also are “concatenated” in the sense

the emergent high-order modes captured in each instability can be used to

32

predict the next-level ﬂow transition.

We start from the lowest

level Q-U transition with a leading-order

unidirectional solution (denoted by M0 mode) as shown in ﬁgure 2.1(a).
Stability analyses show that when H (|α|) is ﬁxed,
associated with ﬁnite-wavelength instability as predicted in ﬁgure 2.2(a) only
occur when |α| (H) goes beyond a certain critical value. Interestingly, we also
ﬁnd these ﬁnite-wavelength modes coexist with M0, and are all appear to be

the unstable modes

wavy. The corresponding vorticity maps show one vortex array in the channel

center (denoted as M1), reminiscent of the traveling wave solutions. When

comparing our results with the phase diagram, we ﬁnd that examining the

high-order modes of the disturbance solution yields a critical curve (black

dashed line) of (αcr, kcr) that agree well with the previous determined

marginal stability curve of U-T transition.

To further verify our conjecture, we explore the T-C transition by examining

the high-order modes captured in the U-T analysis in ﬁgure 2.2 (b,c). As an

example shown in Figure 2.2(b), we ﬁx H = 0.5 and examine unstable modes
by systematically increasing |α|. As 3.5 < |α| < 10, the stability analyses
yield only one branch of solution, corresponding to the leading-order traveling
wave solutions (black solid line, M1 mode). When |α| goes beyond 10, the
eigenvalue search picks another branch of solutions (marked by dashed line)

that represents the high-order modes featured by double vortex arrays (denoted
by M2). Likewise, as shown in panel(c), when ﬁxing α = −4 and increasing
H from 0.45 (i.e., U-T transition) to approximately 1.5, more complex patterns

are captured. Especially in a wide channel, double (M2) and even triple vortex

33

Figure 2.3: Active ﬂow past a ﬁxed cylinder of R = 0.1 in the channel center:
rescaled (a) u and (b) v velocity component, and the surface distribution of (c)
pressure force, (d) viscous force, and (e) active force. In (a)-(e), the upper row
are the simulation results when choosing α = −2.0, and the bottom row shows
the normalized analytical results.

arrays (denoted by M3) exhibit, which suggests that the chaotic motions are

in fact resultant from nonlinear interactions between multiple unstable modes
with different length and time scales. The transition occurs at Hcr ≈ 0.9. When
sweeping values of α and H, the obtained critical curve that characterizes the

mode switch from M1 to M2 and M3 indeed agree with the borderline for the

T-C transition on the phase diagram.

2.4.3 Active Transport of Circular Cylinders

Next, we study how macro obstacles transport in the spontaneous coherent

ﬂows in the rectangular channel, which will be determined by the resultant

hydrodynamic and nematic forces. To facilitate analytical manipulation and

analysis, we consider the case of circular cylinder with possibly the simplest

geometry.

Fixed Cylinder. We ﬁrst consider the case of a circular cylinder ﬁxed in the

channel center. As shown by [19], it is convenient to introduce a 2D stream

as
function ϕ(x, y)
u(cid:48) = (u(cid:48), v(cid:48)) = (∂ϕ/∂y,−∂ϕ/∂x). Substituting ϕ(x, y) into the governing

disturbance

velocity

ﬁeld

to

represent

the

34

equations and assuming isotropy, we can derive a reduced linearized equation

∂t + (cid:0) α

∂∆2ϕ

4 + 4dR

for ϕ as:

bi-harmonic operator using a scalar function such that ∆2ϕ = ψ, leading to a

steady-state equation of an elliptic form:

(cid:1) ∆2ϕ − dT ∆3ϕ = 0. We then replace the
(cid:18)α + 16dR

(cid:19)(cid:19)

(2.11)

4dT

ψ = 0,

(cid:18)

∆ −

which will be solved in the half domain x ≥ 0, −H ≤ y ≤ H, and on the
cylinder surface x2 + y2 = R2 with R the cylinder radius.

We ﬁrst seek the general solution of ψ, and then solve an inhomogeneous

Poisson equation to obtain the scalar stream function ϕ. After some algebraic

manipulations, we have found that the solution can be written as:



(cid:80)N

n=1



+Re

ϕ(ˆx,ˆy)

H = 2C0

π sin(cid:0) π
2 ˆy(cid:1)
 cosh ξn sinh ξn ˆy
 sinh ηn cosh ηn ˆy

−ˆy sinh ξn cosh ξn ˆy

−ˆy cosh ηn sinh ηn ˆy

 eiξn ˆx
 eiηn ˆx

Cn

+Dn


 ,

(2.12)

where ˆx = x/H and ˆy = y/H, and Cn and Dn are unknown coefﬁcients. In

the above, the ﬁrst part is attributed to the long-wavelength weak ﬂow solution

of the Q-U instability without cylinder. The second part satisﬁes the no-slip

condition on the top and bottom walls by using an expansion of

Papkovich-Fadle (PF) eigenfunctions ([50]). The corresponding eigenvalues

ξn and ηn
(n = 1, 2, 3...) are the complex roots of the equations:
ξ − sinh ξ cosh ξ = 0 and η + sinh η cosh η = 0, with positive real and
imaginary parts. The roots are ordered as 0 < Imξ1 < Imξ2... and

35

0 < Imη1 < Imη2... where “Im” denotes the imaginary part. Applying the

no-slip condition on the cylinder surface enable us to solve for Cn and Dn
(n > 0) in terms of C0. The solution of ϕ (and hence u(cid:48)) at x < 0 can be
obtained by using the symmetry condition ϕ (ˆx, ˆy) = ϕ (−ˆx, ˆy). Next, we
solve for the disturbance nematic ﬁeld D(cid:48) through an elliptic equation:

dT ∆D(cid:48) − 4dRD(cid:48) = −1
4

which yields the solution for D

,

11 = −D(cid:48)
D(cid:48)

22

−ˆy cosh ηn sinh ηn ˆy

(cid:16)∇u(cid:48) + (∇u(cid:48))T(cid:17)

 ,
 eiηn ˆx
 eiξn ˆx

 .
 eiξn ˆx
 eiηn ˆx

 sinh ηn cosh ηn ˆy
 cosh ξn sinh ξn ˆy
2αH sin(cid:0) π
2 ˆy(cid:1)
 cosh ξn sinh ξn ˆy
 sinh ηn cosh ηn ˆy

−ˆy cosh ηn sinh ηn ˆy

−ˆy sinh ξn cosh ξn ˆy

−ˆy sinh ξn cosh ξn ˆy

D(cid:48)
12 = πC0

(2.13)

(2.14)

(2.15)







= Re

N(cid:80)

n=1

En

+Fn

+Re

N(cid:80)

n=1

Gn

+Hn

with the unknown coefﬁcients (En, Fn, Gn, Hn) being solved by applying the
no-ﬂux condition on the cylinder surface. In the end, D(cid:48) at x < 0 is obtained
by observing that D(cid:48)

11 (−ˆx, ˆy) = −D(cid:48)

11 (ˆx, ˆy) and D(cid:48)

12 (ˆx, ˆy) = D(cid:48)

12 (−ˆx, ˆy).

In ﬁgure 2.3 we exhibit the steady-state ﬂow (from left to right) ﬁelds and

surface force distributions, and compare the results obtained by using FEM

and the weak-ﬂow analysis.

In general, it is observed that the simulation

36

results are in good agreement with the analytical predictions. Not surprisingly,

the ﬂow ﬁeld in panels(a,b) resembles the classical case of an incoming ﬂow

past a cylinder. To highlight the near surface dynamics, we break the traction

exerted on the cylinder surface into the contributions from pressure ﬁeld
(panel(c), fp = −pn), viscous force (panel(d), fv =
active (nematic) force (panel(e),

fa = α (D − I/2) · n),

(cid:16)∇u + (∇u)T(cid:17) · n), and

respectively.

Performing surface integration of traction yields a net rightward force exerted

on the cylinder, which is largely attributed to the pressure and viscous forces

due to their fore-aft asymmetry. The net nematic force is approximately zero

due to the near-symmetric structure as estimated by the weak-ﬂow analysis,

although the structure will become more asymmetric in a ﬁnite unidirectional

ﬂow which has to be resolved by numerical simulations.

Free-Moving Cylinder. When released in the active ﬂow, the cylinder will

be dragged into motion. As shown in ﬁgure 2.4(a), in the weak ﬂow regime
where |α| is small, it is seen that the cylinder initially accelerates to reach a
maximum, and then relaxes to a steady-state migration speed. Such an initial

pulsative behavior is due to the hyperbolic nature of equation (2.2), somewhat

reminiscent of viscoelastic ﬂuids of Maxwell type where similar start-up

behaviors are observed due to the propagation of stress waves ([51]).

Furthermore, linearity of the weak-ﬂow governing equations (2.11) and (2.13)

permits analytical manipulations to obtain the steady-state cylinder migration
velocity U in x−direction by decomposing the total force into two parts as
Fx = U F (1)
respectively represents the net force

x . Here F (1)
x

x + C0F (2)

and F (2)
x

exerted on a cylinder moving in a quiescent Stokes ﬂow at a constant speed

37

Figure 2.4: Transport of free-moving cylinders in active ﬂows. (a) FEM results
of typical time evolutions of the migration velocities for cylinders of different
sizes when choosing α = −2.0 and H = 0.5. (b) Solid lines represent the
normalized steady-state migration velocities as functions of re-scaled radius.
Symbols are the FEM results. (c) FEM results of cylinder trajectories (R = 0.1)
in traveling wave ﬂows. (d) FEM results of time-averaged velocity proﬁles at
different cylinder volume fractions when choosing R = 0.1 and H = 0.5. The
solid and dashed lines correspond to α = −4.0 and −6.0, respectively.

U = 1, and on a ﬁxed cylinder by setting the unknown coefﬁcient C0 = 1.

After calculating F (1)
x

by following a similar solution procedure via PF

eigenfunction expansion and substituting the previously computed results for

the ﬁxed cylinder when choosing the critical values αcr = −(cid:0) π

(cid:1)2 dT − 16dR

H

x /F (1)
x

for a free-moving cylinder (i.e., Fx = 0). In ﬁgure 2.4(b),

associated with the I-A instability, the migration velocity is then obtained by
U = −C0F (2)
we
rescale
steady-state migration velocities by the maximum
u−component of the unidirectional ﬂow (umax), and plot them in terms of the
is re-scaled by the channel width 2H. Not
dimensionless radius that

the

surprisingly, we have found that the all the analytical results computed at

different channel widths collapse into one master curve, and, again, are in

agreement with numerical simulation results.

In the strong ﬂow regimes, the nonlinear dynamics of the free-moving

cylinder need to be solved numerically. For cylinders whose size is much

38

smaller than the channel width, we expect that they will approximately follow

the coherent ﬂow patterns. As R increases,

the cylinder motion may

signiﬁcantly impact topological structures of active ﬂows, and hence may lead

to rich dynamical behaviors. Indeed, as shown in ﬁgure 2.4(c) (also see movie

S1) for three typical trajectories in traveling wave ﬂows, a small cylinder

oscillates periodically between the two parallel plates by following the wavy

ﬂow patterns. On the other hand, a ﬁnite-size cylinder, whose diameter is close

to the half-wavelength of the traveling waves, seems to hardly catch the ﬂow

direction variations.

Instead, the cylinder “sticks” to one plate such that it

periodically hits and bounces off the boundary.

Increasing the cylinder size

may show some blockage effects to impede active ﬂow development, and lead

to a much slower migration velocity as more viscous damping force exerted on

the cylinder. The blockage effect is further examined by studying transport

behaviors of multiple cylinders in active ﬂows. In ﬁgure 2.4 and movie S2, we

show examples of the time-averaged velocity proﬁles of the resultant active at

various cylinder volume fractions φ computed at R = 0.1 and H = 0.5, where

signiﬁcant ﬂow damping occurs when φ goes beyond 10% whilst the mixture

still shows active transport behavior at a even higher volume fraction (e.g.,

φ = 30%).

2.5 Conclusion

In this work we have combined theoretical analysis and numerical simulations

in studying ﬂuid and obstacle transport in a dilute 2D apolar active ﬂuid

conﬁned in a rectangular channel by using the BQ-tensor model ([19]).

39

Compared to our previous work by [1], we have not only revealed the coherent

structures in different parameter regimes but also, more importantly, precisely

predicted the ﬂow transitions that cover the entire phase diagram.

To

accomplish this goal, we ﬁrst derived analytical solutions of steady

unidirectional ﬂows with ﬁnite strengths, which permit us to perform more

rigorous stability analyses for the U-T instability without ad-hoc assumptions

that were used by [1]. Moreover, by examining the existence of concatenated

high-order wave modes from stability analyses at different levels, we were

able to delineate the borderlines for both the U-T and T-C transitions. When a

circular cylinder is placed in the apolar active ﬂuid, we have been able to

analytically solve for the ﬂow ﬁelds and particle dynamics in the weak-ﬂow

regime for cylinders that is either ﬁxed or free-moving.

In strong coherent

ﬂows of traveling waves, we have observed intriguing trajectories for cylinders

of different sizes. Suppression of active ﬂow development has been observed

when cylinder diameter is comparable of channel size, as well as when

multiple cylinders are transported in active ﬂows at a high volume fraction.

We expect this study will shed light upon the new physics and novel transport

phenomena of geometrically conﬁned active ﬂuids, which potentially suggest

novel mechanisms of ﬂuid manipulation in microﬂuidic devices.

40

3

SCALING LAW OF BROWNIAN ROTATION IN CONCENTRATED

HARD-ROD SUSPENSIONS

3.1 Abstract

Self-diffusion in rod suspensions are subject to topological constraints because

of steric interactions. This topological effect may be anisotropic when rods are

nematically-aligned with their neighbors, and remains a challenge to existing

theories. Via a classical Doi-Onsager kinetic model with the Maier-Saupe

steric potential, we characterize the long-time rotational Brownian diffusivity

Dr for dense suspensions of hard rods of ﬁnite aspect ratios, based on

quadratic orientation auto-correlation functions. Furthermore, we show that

the computed non-monotonic scalings of Dr as a function of volume fraction

can be accurately predicted by a new tube model, which extends Doi and

Edwards’ semi-dilute theory from the isotropic phase to the nematic phase.

3.2 Introduction

Understanding diffusive transport behaviors of concentrated suspensions of

rodlike microparticles (e.g., polymers, macromolecules, and colloids) in dense

systems is of crucial importance in biology, physics, chemical engineering,

and material sciences where the many-body interactions lead to intriguing

pattern formation and material properties [52, 53, 54]. For example, the

short-range steric interactions may dominate over hydrodynamic (or other

41

long-range interactions) at high volume fractions to drive isotropic-anisotropic

order transition [32, 29, 55].

Various models have been developed to reveal the physical mechanisms of

steric interactions in rod suspensions [56, 29, 55], where each rod may

frequently collide with multiple neighbors. This challenges the theories of

cluster-expansion type [57, 58], because of the fast-growing complexity of

such expansion beyond the pairwise-level, i.e., beyond B2 in virial expansion

of free energy. Alternatively, one may construct mean-ﬁeld type models by

embedding one rod (a probe) into the environment of sterically interacting

neighbors. The classical Doi and Edwards (DE hereafter) tube model [29]

views such topologically-conﬁned motion a certain “reptation” dynamics as

the probe slithers among entangled rod meshwork [59, 29]. Because of the

anisotropic shape of rods, the motion of this probe is conﬁned within an

imaginary tube formed by its neighbors, and typically exhibits restricted

rotational motion (characterized by the diffusion coefﬁcient Dr) within the

tube. When estimating the small rotation angle of the probe before it
disengages from the tube, DE’s model predicts Dr/Dr0 ∼ (n(cid:96)3)−2 in the
semi-dilute isotropic regime, where n is the rod number density, (cid:96) is the rod

length, and Dr0 is the rotational diffusivity of an isolated rod with no

neighbors. It has been veriﬁed in several experiments and particle simulation

studies [59, 60, 61, 62, 63, 64]. Nevertheless, we emphasize that most of the

prior studies focus on the dilute or semi-dilute regime with Dr being extracted

from the relaxation time of orientation auto-correlation functions (OACFs)

which typically follows a simple exponential decay. However, the OACFs in

42

the nematic regime may exhibit much more complex behaviors due to phase

transitions, which prevent them from being ﬁtted into analytically trackable

forms that permit clear deﬁnitions of the relaxation time scale(s). There is a

general lack of accurate theoretical models for measuring Dr in the ordered

phases (nematic or smectic) with strong steric conﬁnement effect.

In this Letter, we characterize the long-time Brownian rotational diffusivity

Dr for dense suspensions of rods with ﬁnite aspect ratios via a Doi-Onsager

model where a Maier-Saupe (MS) potential

is employed to capture the

(I)sotropic-(N)ematic phase transition. When the hydrodynamic interactions

are neglected, we show that the measured non-monotonic Dr can be explained

by a modiﬁed DE’s tube model, covering the entire range of volume fraction

from dilute systems in the isotropic phase to dense networks in the nematic

regime. By using a Geometric Constrained Brownian Dynamics (GCBD)

algorithm [65], we build our theory based on extensive numerical simulation

data of N Brownian spherocylinders, each with length (cid:96), diameter b, and

orientation p, in a cubic domain of length L. This method veriﬁed the classic

Brownian spherocylinders without hydrodynamic interactions investigated by

(cid:16) b3

6 + (cid:96)b2

4

(cid:17)

[66] (see more details in Refs [65]). Our simulations cover a wide range of

aspect ratio γ = (cid:96)/b and extend to volume fraction φ = πn

over

50% (i.e., close to the Nematic-Smectic phase transition [66]), where

n = N/L3 is the particle number density. As shown by the open squares in

Fig. 3.1, we perform ensemble average and calculate the (global) orientational

order parameter S = 3
2N

2 and ˆn the average orientation of
spherocylinders (i.e., the director [32, 29]), which captures ﬁrst-order I-N

(cid:80)N
i=1 (pi · ˆn)2 − 1

43

transition as bifurcations of S(φ) curves.

3.3 Numerical method

The simulation method based on geometric constraints in this work was

developed by [65]. The method is applicable for general smooth-shaped rigid

bodies but we focus on spherocylinders in this work. We summarize it here

brieﬂy.

3.3.1 Complementarity formulation and equation of motion
In general, the velocity U = (U1, Ω1, U2, Ω2, . . . ) of all spherocylinders in
the system can be decomposed as the collision velocity U c and non-collisional
velocity U nc:

U = U nc + U c.

(3.1)

In this work the non-collisional motion refers to the translational and rotational

Brownian motion.

The constraints are simply that the minimal separation distance Φα(C)
between each pair α of close rigid objects remain non-negative for all
conﬁgurations C. In total there are nc = O(N ) close particle-particle pairs that
may potentially collide. For each pair α, we denote the collision force

magnitude between this pair of particles as γα. The pair (Φα, γα) must satisfy

one of two conditions:

No contact: Φα > 0,

γα = 0.

In contact: Φα = 0,

γα > 0.

(3.2)

(3.3)

44

Mathematically, this yields a complementarity problem, written vertically over

the set of pairs as:

0 ≤ Φ ⊥ γ ≥ 0,

(3.4)

where Φ = (Φ0, Φ1, ...) ∈ Rnc and γ = (γ0, γ1, ...) ∈ Rnc denote the
collections of minimal distances and contact force magnitudes for each α.

For N rigid particles, let Dα ∈ R6N be a sparse column vector mapping the
magnitude γα to the collision force (and torque) vector applied to each particle.
We deﬁne a sparse matrix D ∈ R6N×nc as a collection of all Dα:
D = [D0 D1 . . . Dnc] ∈ R6N×nc, F c = Dγ.

(3.5)

For each pair α, Dα deﬁned in this way has 12 non-zero entries for

non-spherical shapes, corresponding to 3 translational and 3 rotational degrees

of freedom for each particle. For the special cases of two spheres, Dα has 6

non-zero entries, as normal collision forces induce no torques on spheres. The

entries of each Dα are explicitly given in [67]. For rigid particles, there is an
important relation between the transpose of D and Φ:

DT = (∇CΦ) .

(3.6)

Euler angles.

Combining (3.4) and (3.5), we reach a differential variational inequality

(DVI) for N particles and all constraints.

˙C = U nc + MDγ,
0 ≤ Φ(C) ⊥ γ ≥ 0.

(3.7a)

(3.7b)

45

Here U nc, M, D, and Φ depend only on the geometric conﬁguration C. This
DVI is solvable and integrable over time once a relation between the
conﬁguration C and the collision force γ is supplied, that is, a timestepping
scheme. Higher order schemes such as the Runge-Kutta scheme can be used,

but for simplicity of presentation we derive a ﬁrst order Euler integration

scheme here.

With given timestep ∆t,

(3.7)

can be

linearized as

a Linear

Complementarity Problem and then converted to a Constrained Quadratic

Programming and be efﬁciently solved with gradient descent methods.

Detailed derivation, discussion, and veriﬁcation of our algorithm can be found

in [65]. In this study, we simulate spherocylinders of length (cid:96) = 1µm and
aspect ratio 5 ≤ γ ≤ 50 for volume fractions 1% ≤ φ ≤ 50%, in cubic
computation domain with the length varying between the maximum value

L = 10(cid:96) (for γ = 5) and the minimum value L = 4(cid:96) (for γ = 50). The time
step is chosen as ∆t = 1 ∼ 5 × 10−5s.

3.4 Results

3.4.1 Equilibrium PDF of the Doi-Onsager Model

A classical Doi-Onsager kinetic model [29] accurately describes the above

mentioned critical behaviors via a time evolution equation

(cid:20)

(cid:21)

= DrRRR ·

∂Ψ
∂t

RRRΨ +

Ψ
kBT

RRRUs (p)

,

(3.8)

where Ψ is a probability distribution function (PDF) for rod’s orientation (unit)
vector p, and RRR = p × (∂/∂p) is a rotational operator.
characterizes the long-time (or “collective”) Brownian rotational diffusion for

In the above, Dr

46

Figure 3.1: Equilibrium orientational order parameter S as a function of φ.
The open squares represent the GCBD simulation data. The solid lines are
calculated using Ψeq in Eq. (3.12). Inset: Instantaneous snapshots at φ = 20%
and 40% for γ = 10 correspond to the I and N phase, respectively.

an equilibrium rod system [29, 64]. A mean-ﬁeld torque is imposed on the rod’s
rotational dynamics via the MS steric potential [31] Us (p) = −ζ0pp : Q due
to the crowdedness nature of densely packed rods, where Q = (cid:104)pp − I/3(cid:105) is
the nematic-order-parameter tensor with (cid:104)(cid:105) a conﬁgurational average of p, and
ζ0 is a strength coefﬁcient. At equilibrium (i.e., ∂Ψ
∂t = 0), we rewrite the MS
s = −ζ0pp : Qeq with Qeq = S (ˆnˆn − I/3), and then balance

potential as U (eq)

the rotational diffusion and the steric torque to obtain the equilibrium PDF.

The rod system can be described by a probability distribution function Ψ (p)

as a function of the orientation vector p, evolved through the Smoluchowski

equation (1) in the main text. The equilibrium solution satisﬁes

(cid:79)p ln Ψeq =

2ζ0
kBT

(I − pp) · Q · p

(3.9)

In an isotropic state, Ψ = 1

4π with the coefﬁcient ζ0 = 0.

In the nematic

regime, we assume the rods are distributed axisymmetrically around the

47

φS00.10.20.30.40.500.20.40.60.81γ = 7γ = 5γ = 10γ = 15γ = 20φ = 20%(γ = 10)nNIφ = 40%(γ = 10)∧γ = 50Figure 3.2: Examples of orientational PDF for γ = 10 at φ = 20%, 30%, 40%.
Open squares represent the GCBD data, and the solid lines represent the results
by using Ψeq in Eq.(2) in the main text.

director ˆn, yielding [30]

(cid:19)

ˆn ˆn − I
3

(cid:18)
Ψ(cid:0)3 cos2 θ − 1(cid:1) sin θdθdφ.

,

(3.10)

(3.11)

Q = S

(cid:90) 2π

(cid:90) π

0

0

S =

1
2

where the order parameter S can be expressed as

By substituting (3.10) and (3.11) into (3.9), we can solve Ψ (θ) in spherical

coordinates as

Ψeq (θ) =

2π(cid:82) π

with the coefﬁcient χ =

Sζ0
2kBT

exp (χ cos (2θ))

0 exp (χ cos (2θ(cid:48))) sin θ(cid:48)dθ(cid:48) ,
. For a given φ and γ, χ can be measured

(3.12)

through the following steps: (1) take log on both sides of (3.12), yielding

48

θ/πΨeq(θ)00.20.40.60.810123I-N40%20%30%φΨeqGCBDIsotropicIsotropicIsotropicNematiclog Ψeq = χ cos (2θ) + const; (2) construct a discrete orientation probability

distribution from particle simulation by dividing θ uniformly into small
intervals (cid:52)θ = π/400; (3) take log on the discrete orientation probability
distribution and convert it to be a function of cos (2θ) and use the linear ﬁtting

to get χ. After χ is measured, a continuous PDF can be constructed based on

(3.12). As examples shown in Fig. 3.2 for γ = 10 cases, the reconstructed PDF

Ψeq accurately represent the numerical data.

3.4.2 Measuring Dr from quadratic OACFs

We examine the orientation relaxation process of a free probe in an

equilibrium particle bath. As examples of γ = 20 shown in Fig. 3.3(a), direct
long-time (cid:104)p (t) · p (0)(cid:105) often show over-complicated
evolution behaviors as the system deviates from the isotropic cases (e.g.,

measurements of

φ = 10%) whose OACFs decay exponentially to zero [29]. During the I-N

transition (e.g., φ = 17% and 20%), however, the functions may exhibit

extraordinary slow decays after the initial transient [68]. As an alternative

(cid:68)

(p (t) · p (0))2(cid:69)

approach shown in Fig. 3.3(b), when measuring the quadratic form

(or equivalently (cid:104)P2 (p (t) · p (0))(cid:105), P2

represents the

Legendre polynomial of order 2 [69, 29, 63]), we observe that for all cases, the

OACFs approach constant values after the initial decay, which can be

accurately ﬁtted into a universal form

(cid:68)

(p (t) · p (0))2(cid:69)

= (1 − A∞) exp

(cid:19)

(cid:18)

− t
τR

+ A∞

(3.13)

with a (single) correlation time τR and a steady-state value A∞. In the isotropic

regime (φ < 15%), A∞ = 1/3 which has been proven analytically [29]; while

49

(cid:68)

(p (t) · p (0))2(cid:69)

Figure 3.3: Orientation relaxation characterized by (cid:104)p (t) · p (0)(cid:105) (a) and
(b) across φ for γ = 20. The open symbols represent the
GCBD data. Schematics of conﬁgurational space of p(t) are highlighted at the
bottom for the different phases. In panel(b), the solid lines represent the ﬁtting
function in Eq. (3.13). Inset: Zoom-in window highlights that the short-time
relaxation of OACF ﬁts an exponential decay governed by a single time scale
τR.
in the deep nematic regime (φ > 20%), A∞ → 1 and τR → 0, suggesting
strong orientation correlations.

These behaviors can be qualitatively understood by examining the long-time

orientation distribution of p(t) with respect to p(0). As sketched in the bottom

of Fig. 3.3(a), the conﬁgurational space of p(t) exhibits spherical symmetry
with the formed angle α = arccos (p (t) · p (0)) ∈ [0, π]. In the deep nematic

50

<p(t)⋅p(0)>00.511.5200.20.40.60.81φ = 20%I-Nφ = 40%φ = 17%φ = 15%φ = 10%IN(a)I-NI-Np(0)p(0)p(t)p(t)II-Np(0)p(t)Np(t)<(p(t)⋅p(0))2>00.511.520.40.60.81φ = 20%I-NINφ = 40%φ = 17%φ = 15%φ = 10%IN(b)I-NI-Np(0)p(0)p(0)p(0)p(t)p(t)1/300.10.20.40.60.81time (seconds)regime, it approximately reduces to a small and narrow cone-shaped domain

with an acute angle α when the probe rotation is strongly conﬁned around p(0).

The cone shape can be described by its maximum opening angle estimated
as αmax ∼ π/2 − arcsin(S). During the I-N transition where steric effects
are relatively weak, p(t) can simultaneously fall into two large cones that are

symmetric about α = π/2 but with unequal probabilities as marked by different

colors for the cones. A biased-distribution of p(t) towards α = 0 will cause a

slower relaxation process (e.g., φ = 17%) towards a ﬁnite equilibrium value as
t → ∞. On the other hand, as shown in panel(b), a quadratic OACF maintains
. Hence the
the fore-aft symmetry, i.e.,

(−p (t) · p (0))2(cid:69)

(p (t) · p (0))2(cid:69)

(cid:68)

(cid:68)

=

conﬁgurational domain of p(t) automatically switches from a half sphere in the

isotropic regime to a single cone domain in the nematic regime, leading to all

positive correlations at equilibrium (i.e., A∞ > 0).

Next, we construct a mathematical model to connect the quadratic OACF

(cid:68)

(p (t) · p (0))2(cid:69)

with the Brownian diffusivity Dr deﬁned in Eq. (3.8). We follow DE’s theory

=(cid:82) (pp : p(cid:48)p(cid:48)) G (p, p(cid:48); t) Ψeq (p(cid:48)) dpdp(cid:48)

[69, 29] to deﬁne
where a Green’s function G (p, p(cid:48); t) represents a conditional (or joint)
probability that the probe is in the direction p at time t, given that it was in the

direction p(0) at t = 0. It can be derived from the Smoluchowski equation that
the time evolution of G (p, p(cid:48); t) follows ∂G
, with
the initial condition G (p, p(cid:48); t = 0) = δ (p − p(cid:48)). Then the time derivative of

kBTRRRUs

∂t = DrRRR ·(cid:16)RRRG + G

(cid:17)

51

(cid:68)

(p (t) · p (0))2(cid:69)
(cid:68)
(p (t) · p (0))2(cid:69)
(cid:90) (cid:18)

∂
∂t

becomes

=

Dr

RRR2 (pp) : p(cid:48)p(cid:48) − 1
kBT

(cid:19)

GΨeq (p(cid:48)) dpdp(cid:48)

(RRR (pp) : p(cid:48)p(cid:48)) · RRRUs

(3.14)
with the ﬁrst term RRR2 (pp) : p(cid:48)p(cid:48) = −6pp : p(cid:48)p(cid:48) + 2. Substitute expression of
U (eq)

into the second term, it derives

s

(cid:0)( ˆn · p(cid:48))( ˆn · p)(p · p(cid:48)) − ( ˆn · p)2 (pp : p(cid:48)p(cid:48))(cid:1)

(RRR (pp) : p(cid:48)p(cid:48)) · RRRUs =

− 1
kBT

4ζ0S
kBT

(3.15)

Then we write the triple product term as ( ˆn· p(cid:48))( ˆn· p)(p· p(cid:48)) = |( ˆn· p(cid:48))( ˆn·
√

p)||p · p(cid:48)| ≈ ( ˆn · p)2|p · p(cid:48)|. A Taylor expansion of the function f (x) =
to the ﬁrst order near x = 1 yields |p · p(cid:48)| ≈ 1
2. The term
( ˆn· p)2 can approximated through the steady state value of (cid:104)pp : ˆn ˆn(cid:105), which is
3. Substitute these into Eq.(3.15)

(cid:1) : ˆn ˆn and it gives ( ˆn · p)2 ≈ 2

(cid:0)Qeq + I

2(pp : p(cid:48)p(cid:48)) + 1

3S + 1
and Eq.(3.16), it gives the evolution over time

x

3

(cid:68)

(p (t) · p (0))2(cid:69) − (q +

(cid:20)
(cid:1) captures the intrinsic I-N transition. After

(1 + q)

(3.16)

(cid:21)

1
3

)

= −6Dr

9kBT S(cid:0)S + 1

2

(p (t) · p (0))2(cid:69)

(cid:68)

∂
∂t

with q (S) = 2ζ0

(cid:68)

(p (t) · p (0))2(cid:69)

integral, we have the time-dependent solution

=

2

3 (1 + q)

exp [−6 (1 + q) Drt] +

3q + 1
3 (1 + q)

.

(3.17)

Therefore, comparing Eq. (3.13) and (3.17) yields A∞ = 3q+1

3(1+q) and

Dr =

1 − A∞
4τR

52

(3.18)

Figure 3.4: Dr/Dr0 as functions of n(cid:96)3 across φ for various γ. The open squares
represent the data extracted from GCBD simulations via Eq. (3.18). The solid
and dashed lines represent the scaling laws in Eq. (3.20) and (3.25), respectively.
The solid red square on the γ = 20 curve marks the location of the data for
ﬁtting β.
p(cid:48) as functions of S, evaluated by
using Ψeq (solid lines) and the quadratic-closure approximation (dashed lines)
[2], respectively.

p(cid:48) and (cid:104)|p · p(cid:48)|(cid:105)p

Inset: (cid:104)|p × p(cid:48)|(cid:105)p

which enable us to extract Dr from OACF measurements. Note that it recovers

the classical results A∞ = 1/3 in the isotropic regime when q = 0 [29, 63].

For all cases, we present the ratio Dr/Dr0 (open squares) as functions of a

dimensionless number density n(cid:96)3 in Fig. 3.4. Although Dr/Dr0 decays

monotonically across the volume fraction for short rods (e.g., γ = 5), reverse

growths are seen for rods of ﬁnite and large aspect ratios (γ > 10) during the

I-N transition. In the deep nematic regime, Dr/Dr0 becomes to be descending

again, suggesting that

the probe’s self-diffusion motion has been fully

constrained.

3.4.3 Scaling law delineated by modiﬁed tube model

To understand the complex scaling laws of Dr, as sketched in Fig. 3.5(a), we

envision the rotation of a probe to be conﬁned by surrounding neighbors that

53

n    3Dr / Dr010-11001011021030.20.40.60.81tube dilationγ = 5βγ = 10γ = 50γ = 20-2S00.20.40.60.8100.20.40.60.81〈|p⋅p’|〉〈|p×p’|〉are nematically aligned. Inspired by DE’s classical theory [29], we enclose the

probe in a virtual tube. Unlike DE’s original semi-dilute model which only

considers one effective side-by-side interaction between an inﬁnitesimally-thin

probe and its nearest neighbor, we simultaneously track multiple interactions

(contacts) attributed from two types of collisions, side-side collisions (denoted

by “type I”) and end-side/end-end collisions (denoted by “type II”). Intuitively,

type II collisions, which were neglected in DE’s model, may signiﬁcantly

change the rotational motion because the collision forces have longer moment

arms when it occurs at the probe’s end.

The conceptual

image in Fig. 3.5(a) of a restricted collision-driven

rotational motion is supported by a series of GCBD simulation data that reveal

intriguing geometric features. As shown in Fig. 3.5(b), we examine the

fractional contributions of the two types of collisions relative to the total

number of collisions, i.e., λI,II = NI,II/ (NI + NII), where NI,II represent

the time-average collision numbers. For short rods (e.g., γ = 5), the type II

collisions dominate in both the isotropic and nematic regimes; while for

relatively long rods (e.g., γ = 20), the type I collisions are more frequent in

the isotropic regime, but

they are signiﬁcantly reduced during the I-N

transition, and become much less compared to the type II collisions in the

nematic regime. Not surprisingly, all data for a re-scaled collision number

ratio λI/ (γλII) approximately collapse onto a master curve since the collision

numbers are purely determined by the surface areas of the sidewall and the

ends. In panel(d), we deﬁne an orientational pair-correlation function (OPCF)

in

the

rod’s

transverse

direction

[70,

71,

72]

54

(cid:80)

i(cid:54)=j

(cid:68)
(cid:80)

(cid:16)
r⊥−
√

|rij|2−(rij·pi)2(cid:17)(cid:69)
√
|rij|2−(rij·pi)2(cid:17)(cid:69)

(cid:68)

(cid:16)

P2(pi·pj)δ
r⊥−

i(cid:54)=j

δ

g (r⊥) =

where rij

connects

the

center-of-mass positions of the ith and jth rods, and r⊥ measured the

projected distance between the rods in perpendicular direction of pi. Again,
all data for a re-scaled OPCF g(r⊥)−gmin
gmax−gmin
and minimum (when r⊥ > 4b) values of g, especially for rods of high aspect

by using the maximum (when r⊥ = b)

ratios, approximately fall onto a master curve. These shape-independent

features provides guidance of how to construct analytical models that

accurately reﬂect the relative contributions of two types of steric interactions.

As illustrated in Fig. 3.5(a), we assume a probe oriented in the p-direction

has a thickness b and length (cid:96), sitting among nematically-aligned neighbors

whose orientation distribution can be described by Ψeq. We estimate the

number of penetrations into a virtual tube that has the same rod length (cid:96). As

suggested by Fig. 3.5(d), the tube radius aT approximately varies between b

and 3b where the OPCFs exhibit signiﬁcant correlations. Following DE’s

approach [29], we consider the number of rods that penetrate a small surface

a

be

unit

normal

vector n can

as
area ∆S with
∆N = n(cid:96)∆S(cid:104)|p(cid:48) · n|(cid:105) where the conﬁgurational average is performed in a
representative volume along the p(cid:48)−direction. For the sidewall intersections
(i.e., n⊥p) arising from the type I collisions, we reuse DE’s calculation [29] to
sum ∆N over half of the sidewall area SI = πaT (cid:96) by taking into account the

calculated

rod,

vicinal

from a
p(cid:48) = 2naT (cid:96)2(cid:104)|p × p(cid:48)|(cid:105)p

double-penetration
NI = SIn(cid:96)(cid:104)|p · n|(cid:105)p
p(cid:48) (see the derivation details in Ref
[29]). Here (cid:104)(cid:105)p
p(cid:48) represents conﬁgurational averages over both p and p(cid:48).
Likewise, the number of rods intersecting from the two ends (i.e., n//p and the

number

obtain

and

the

55

Figure 3.5: GCBD simulations motivate a new tube model. (a) Schematic of
a virtual tube with penetrations from both the sidewall and the two ends. (b)
Fractional contributions of the averaged collision numbers on the sidewall (λI)
and the ends (λII) for γ = 5 and 20 as functions of φ. (c) The collision number
ratio λI/ (γλII) as functions of S for various γ. (d) The OPCF g(r⊥) in the
nematic regime as functions of the dimensionless transverse distance r⊥/b.

total area SII = 2πa2
NII = SIIn(cid:96)(cid:104)|p · n|(cid:105)p
p(cid:48) = 2nπa2
rods only penetrate the tube once.

T ) can be calculated in a straightforward manner as
p(cid:48) by assuming that the surrounding

T (cid:96)(cid:104)|p · p(cid:48)|(cid:105)p

To avoid double-counting the neighboring rods that penetrate from the side

and one end simultaneously, we calculate the total penetration number as

Np = ωINI + ωIINII
= 2ωInaT (cid:96)2(cid:104)|p × p(cid:48)|(cid:105)p

p(cid:48) + 2ωIInπa2

T (cid:96)(cid:104)|p · p(cid:48)|(cid:105)p
p(cid:48).

(3.19)

Since the calculations of NI,II have incorporated the shape effect, according

56

(a)φλI, II00.10.20.30.40.50.20.40.60.8(b)λI, γ = 5λII, γ = 5λI, γ = 20λII, γ = 20SλI / (γλII)00.20.40.60.8100.050.10.15(c)γ = 5γ = 10γ = 20γ = 50r⊥/ b(g - gmin)/(gmax - gmin)1234500.20.40.60.81(d)γ = 50, φ = 40%γ = 5, φ = 40%γ = 10, φ = 30%γ = 10, φ = 40%γ = 20, φ = 30%γ = 20, φ = 40%γ = 50, φ = 30%to Fig. 3.5(c), we assume the weights ωI,II (S) are shape-independent, and

will be selected such that NI and NII coexist in the isotropic regime but NI

vanishes in the “sharply-aligned” limit where all rods are perfectly aligned

(i.e., S = 1). We assume the probe rotates a small angle δφ as it moves for a
distance of rod length (cid:96) during the “reptation” time scale τ0 ∝ D−1
r0 , which
leads to Dr = δφ2/τ0 = βδφ2Dr0 with β a shape-dependent constant.

Considering the fact that the probe encounters these Np penetrations with

equal

Dr/Dr0 = β(cid:0)n(cid:96)3(cid:1)−2(cid:104)

probability,

we

estimate

ωI (cid:104)|p × p(cid:48)|(cid:105)p

p(cid:48) + πωII

δφ

(cid:0) aT

(cid:96)

=

(cid:1)(cid:104)|p · p(cid:48)|(cid:105)p

p(cid:48)

aT

(cid:105)−2

Np(cid:96),

yielding

. The β value

for each γ case can be numerically ﬁtted by matching Dr/Dr0 with the GCBD

data on a point at the onset of the I-N transition (as one example, see the solid

red square marked on the γ = 20 curve in Fig. 3.4) [60, 61, 73, 29, 74].

After β is ﬁxed, we ﬁnd that choosing ωI = (cid:104)|p × p(cid:48)|(cid:105)p

p(cid:48), ωII = (cid:104)|p · p(cid:48)|(cid:105)p
p(cid:48),
and the tube radius aT ≈ 2b leads to excellent agreements between the GCBD
data. The analytical form can be written as

= β(cid:0)n(cid:96)3(cid:1)−2

(cid:20)(cid:16)(cid:104)|p × p(cid:48)|(cid:105)p

p(cid:48)

(cid:17)2

(cid:16)(cid:104)|p · p(cid:48)|(cid:105)p

p(cid:48)

(cid:17)2(cid:21)−2

+

2π
γ

Dr
Dr0

(3.20)

where the two double-conﬁgurational-averages are evaluated by using Ψeq.

The inset of Fig. 3.4 reveals the intrinsic competing mechanism for the the

non-monotonic critical behaviors of Dr. Mathematically, we see that the

enhanced rotational diffusions observed during the I-N transition are caused
by the divergent behaviors of Eq. (3.20) as (cid:104)|p × p(cid:48)|(cid:105)p
p(cid:48) approaches zero at
large S which suggests weakened side-side steric interactions as the probe

reptates. Such “tube dilation” effect is regulated by the enhanced end-steric
effects in the nematic regime as (cid:104)|p · p(cid:48)|(cid:105)p

p(cid:48) dominates, especially for relatively

57

short rods with small γ.

3.4.4 Approximations of the double-conﬁgurational averages

We perform analytical calculation of
(cid:104)|p × p(cid:48)|(cid:105)p

p(cid:48) and (cid:104)|p · p(cid:48)|(cid:105)p

p(cid:48). We ﬁrst deﬁne m = pp − I/3, and rewrite

the two conﬁgurational averages

|p · p(cid:48)| =(cid:112)pp : p(cid:48)p(cid:48) =
|p × p(cid:48)| =(cid:112)1 − pp : p(cid:48)p(cid:48) =

(cid:114)1
(cid:114)2

3

3

√

(cid:114)

1 + 3m : m(cid:48),

1 − 3
2

m : m(cid:48)

(3.21)

(cid:114)1

(cid:18)

3

For the ﬁrst expression in (3.21), we perform Taylor expansions in terms of
m : m(cid:48), leading to

|p · p(cid:48)| =

1 +

3
2

m : m(cid:48) − 9
8

(m : m(cid:48))2 +

27
16

(m : m(cid:48))3 + ...

(3.22)

To proceed, we employ a quadratic closure to express the high-order-moments
in terms of the lower-order moments, leading to (cid:104)
Hence we have

mm...mm(cid:105)p =

QQ...QQ.

N

(cid:122)

(cid:125)(cid:124)

N

(cid:122)

(cid:125)(cid:124)

(cid:123)

(cid:123)

(cid:19)

(cid:114)1

(cid:18)

3

≈

1 +

3
2

Q : Q(cid:48) − 9
8

(cid:104)|p · p(cid:48)|(cid:105)p
p(cid:48)
(Q : Q(cid:48))2 +
√

(cid:114)1

1 + 2S2.

=

3

(cid:19)

(Q : Q(cid:48))3 + ...

27
16

Similarly, the second term in (3.21) can be calculated as

(cid:104)|p × p(cid:48)|(cid:105)p

p(cid:48) ≈

1 − S2.

Then Eq. (3.20) can be rewritten as

≈ β(cid:0)n(cid:96)3(cid:1)−2

(cid:20)(cid:18)

Dr
Dr0

(cid:19)

− 1

(cid:21)−2

.

S2

(3.23)

(3.24)

(3.25)

(cid:114)2
(cid:19)

3

(cid:112)
(cid:18)2π

γ

1 +

+

π
γ

58

As shown by the dashed lines in Fig. 3.4, the above equation provides accurate
yet convenient form without using Ψeq. For thin rods (γ → ∞), Eq. (3.25)
captures the DE’s semi-dilute theory in the isotropic regime when S = 0. It also
recovers DE’s correction in the nematic regime when S > 0 and γ → ∞ by
considering the tube-dilation effect only [75, 29, 76], which, however, becomes
singular in the limit S → 1. We can further predict that such dilation effects
only occur for rods with γ > 2π as the resultant negative sign of the S2 term

directly causes reverse growing. This can be veriﬁed by the γ = 7 case where

Dr/Dr0 approximately follows the(cid:0)n(cid:96)3(cid:1)−2 law in the nematic regime.

3.5 Conclusion

In summary, we have built a Doi-Onsager model for characterizing the

long-time Brownian diffusivity Dr in the nematic regime for concentrated rod

suspensions. The non-monotonic critical behaviors of Dr are accurately

captured via a new tube model. We have shown that when appropriately

incorporating steric hindrances around a probe rod, DE’s tube concept can

elegantly yet accurately describe topologically-constrained rod rotation in a

concentrated meshwork without scrutinizing detailed rod-rod interactions.

While hydrodynamic interactions are neglected, we expect our scaling law

remains accurate for dense rod suspensions where the short-range steric

shielding dominates. We anticipate this model will provide a new angle in

characterizing Brownian rotations in a wide variety of computational and

experimental studies of soft condensed matter.

59

3.6 Future plan

There is still a lack of a thorough understanding of the long-time translational

Brownian diffusivity Dt for dense suspensions of hard rods of ﬁnite aspect

ratios. It is meaningful to explore the scaling of Dt as a function of volume

fraction φ. Besides, it may be a good try to expand to the active system with

hydrodynamics effect, to see what the steric interactions bring to the existing

Taylor dispersion theory, which describes an effect in ﬂuid mechanics where a

shear ﬂow can increase the effective diffusivity of particles.

60

4

Q−TENSOR MODEL FOR UNDULATORY SWIMMING IN A

LIQUID CRYSTAL

4.1 Abstract

Microorganisms take full advantage of anisotropic drag forces when

swimming in low-Reynolds-number ﬂows where the inertia effect

is

negligible. Besides adopting various swimming gaits,

they may exploit

anisotropic properties of the coupling liquid medium for efﬁcient swimming

or locomotion.

In this paper, we build a computation model to study the

nonlinear dynamics of an undulatory microswimmer

in a quasi-2D

liquid-crystal polymer solution. We describe the polymer ﬂuid phase using a
hydrodynamic Q−tensor model for rigid microrods, and treat the swimmer as
a ﬁnite-length ﬁber whose undulatory motion is actuated by imposing

propagating traveling waves on the body curvature.

The ﬂuid-structure

interactions are resolved via an immersed boundary method. Compared to the

dynamics in Newtonian ﬂuids, we observe rich non-Newtonian behaviors of

both enhanced and retarded swimming motions when calculating the

swimmer’s mean speed. We reveal the propulsion mechanism by analyzing the

near-body ﬂow ﬁelds and polymeric force distributions,

together with

asymptotic analysis for an idealized model built upon Taylor’s swimming

sheet.

61

4.2 Introduction

Many organisms live in microﬂuidic environments, either biological or

synthetic, where the ﬂuid inertia is negligible.

In the so-called Stokes (or

creeping) ﬂows, Purcell’s Scallop theorem explains

that performing

time-reversible motions cannot generate directional swimming or locomotion

due to kinematic reversibility ([77]).

Instead, microswimmers often adopt

special propulsion strategies (swimming gaits). Especially for those with

slender shapes or having thin appendages (e.g., C. elegans, E. coli, B. subtilis),

it is well-understood that they can perform undulation or ﬂagellar beating to

generate directional motions by breaking symmetry ([78]). On the other hand,

people have conducted extensive research to reveal rich dynamics and new

propulsion mechanisms of microorganisms that utilize the complex ﬂuids’

intrinsic anisotropy arising from their non-Newtonian physical and rheological

properties. For example, while the extra particle stresses from shear-thinning

viscoelastic (e.g., Oldroyd B) ﬂuids generally impose additional hindrance

effects, it has been found that undulatory swimmers such as C. elegans may

leverage their ﬁnite body-length and polymeric stress relaxation to achieve a

higher swimming speed than that in the Newtonian ﬂuids ([79, 80, 81]).

Furthermore, there has been a growing interest in exploring non-equilibrium

physics of biosynthetic active materials “powered” by many-swimmer

interactions in complex ﬂuids that are capable of exploiting collective

dynamics for useful mechanical work ([3, 5]). Of particular interest is the

study of motile bacteria in lyotropic nematic liquid crystals (LC) wherein the

extra stress generation is directly determined by the LC molecule orientation

62

(or

the so-called orientation order),

leading to intriguing optical and

mechanical properties ([82, 83]).

So far, there have been several models proposed to study the dynamics of

microswimmers in LCs. [82, 84] constructed a hydrodynamic Ericksen-Leslie

(EL) model to solve for the orientation variation of the director ﬁeld in the LC

with motile bacteria (B. subtilis) with prescribed translational motions, and

illustrated how the induced shear determines the ordering patterns (or

birefringent cloud) around the swimmer. [85] studied the motion of Taylor’s

swimming sheet with imposed small-amplitude travelling waves in LC via a

similar EL model. When enforcing the local LC orientation to be tangential to

the undulatory body (i.e., tangential anchoring), they derived the asymptotic

solutions, which, for the ﬁrst time, mathematically reveal how the undulatory

gaits are controled by LC’s material parameters ([86, 87]). Later [88] used the

Immersed Boundary (IB) method to simulate a ﬁnite-length undulatory

swimmer in LCs conﬁned between two parallel walls, when subjected to

large-amplitude actuation and the tangential anchoring condition, and only

found retarded motions compared to the swimming speed in the Newtonian

ﬂuid. However, as they pointed out, the EL model cannot be robustly extended

to more general cases, and may break down when applying other types (e.g.,

homeotropic) of anchoring conditions.

Motivated by the prior

studies and to further

seek quantitative

understanding of undulatory swimming mechanisms in LCs, we present a new
computational framework which is built upon a minimal Q−tensor model
coarse-grained from the Doi-Onsager kinetic model for rigid microrods ([29]).

63

The swimmer is treated as a ﬁnite-length elastic ﬁber whose undulation is

activated by imposing traveling waves on its body curvature. An IB method is

employed to solve the ﬂuid/structure interactions (FSIs) between the swimmer

and LCs.

Interestingly, our simulations reveal both enhanced and retarded

swimming motions in the nematic regime, corresponding to the scenarios

when the swimming direction is parallel and perpendicular to the LC’s

alignment directions, respectively. These observations are different from [88]

who only see retarded motions. Our numerical results are qualitatively

conﬁrmed by performing the asymptotic analysis of Taylor’s swimming sheet

model, and can be further understood by revealing the unique near-body

polymer force distributions.

4.3 Mathematical Model

In the kinetic model by [29], a probability distribution function (PDF)

Ψ(x, p, t) is employed to describe the ensemble dynamic of microrods in

terms of the rod’s center-of-mass (C.O.M.) position x and orientation p
(|p| = 1). When simulating the macro scale hydrodynamic interactions
between rods and ambient ﬂows, it is often preferred to use coarse-grained

PDEs that are derived via proper moment closure to eliminate p due to their

low computation costs.

Consider rodlike LC molecules with fore-aft

symmetry and uniform length (cid:96). In the ﬂuid domain Ωf, a two-dimensional

hydrodynamic LC model can be derived for the second-moment

tensor

D =(cid:82)

p ppΨdp as ([29]):
∇
D +2E : S = 4ζdr(D · D − D : S) − 4dr

(cid:19)

(cid:18)

D − I
2

+ dt∆D,

(4.1)

64

where

and S = (cid:82)

∇
D = ∂D

∂t + u · ∇D − (D · ∇u + ∇uT · D) is the so-called

(cid:0)∇u + ∇uT(cid:1) is the strain-rate tensor,

upper-convected time derivative, E = 1
2

p ppppΨdp is the fourth-moment tensor. The ﬁrst term on the
right-hand-side (R.H.S.) represents the steric alignment effects arising from
the Maier-Saupe (MS) potential UM S (p) = −ζpp : D, with ζ the
dimensionless strength coefﬁcient ([31]). The diffusivity coefﬁcients dt and dr

the so-called Q-tensor,

originate from the translational and rotational Brownian motion, respectively.
The above equation is a minimal “Q−tensor” model since the trace-free
is deﬁned by
(normalized) tensor-order parameter,
Q(x, t) = c−1D − I/2, with c(x, t) = 1 the concentration ﬁeld (i.e.,
zeroth-moment of Ψ) that remains to be homogeneous constant.
In 2D, the
tensor Q has a maximal nonnegative eigenvalue ρ satisfying 0 ≤ ρ ≤ 1
2.
Assuming that ρ is isolated, then we call its associated unit eigenvector the
director m, and 0 ≤ S = 2ρ ≤ 1 the orientation order parameter. The
incompressible ﬂuid velocity u is incorporated via a mean-ﬁeld fashion, and is
governed by the forced Stokes equation: ∇p − µf ∆u = ∇ · σp + fe, where µf
is the ﬂuid viscosity, and fe(x, t) represents the force exerted by the motile

swimmer. The extra stress due to the LC polymer solution is deﬁned as

([75, 89, 76])

(cid:18)

(cid:19)

D − I
2

− 2νkBT ζ (D · D − S : D) +

νkBT β0
2dr(ν(cid:96)3)2 E : S (4.2)
σp = 2νkBT
where ν ∼ (cid:96)−3 is an effective volume fraction, kBT is the thermal energy
scale, β0

(cid:0)ν(cid:96)3(cid:1)−2 represents a crowdedness factor with β0 being an empirical

coefﬁcient. In this study, we will choose a small empirical crowdedness factor

(cid:0)ν(cid:96)3(cid:1)−2 ∼ 10−3 − 10−4 as suggested by [89] and [76].

β0

65

In the Lagrangian frame (ΩL) for the undulatory swimmer of length Ls, its

kinematics can be described by the parametric form X (s, t), with s the local
arc length s ∈ [0, Ls]. The actual shape X (s, t) is determined by penalizing
deviations from actuation imposed on the body curvature κ0 (s, t):

κ0 (s, t) = −k2A0 sin (ks − ωt) ,

(4.3)

starting from an initial conﬁguration given by X(s, t = 0) = (s, A0(sin(ks))).

Equation (4.3) describe the right-propagating traveling waves with amplitude

A0, wavenumber k, and angular frequency ω. Imposing actuation in Eq. (4.3)

produces elastic forces Fe (X) along the body, and effectively drive periodic

shape changes (or swimming gaits). Following [90], the Lagrangian body

force can be derived by performing the variational derivative upon the elastic
energy E, i.e., Fe (X) = − δE[X(s,t)]
δX . Here the total elastic energy E [X]
include the contributions from both stretching (denoted by subscript s) and

bending (denoted by subscript b) deformation

(cid:12)(cid:12)(cid:12)(cid:12) − 1
(cid:19)2

ds +

σ0
b
2

(cid:90)

ΩL

(cid:18)∂2X
∂s2 · n − κ0

(cid:19)2

ds (4.4)

E [X (s, t)] =

σ0
s
2

(cid:90)

ΩL

∂s

(cid:18)(cid:12)(cid:12)(cid:12)(cid:12)∂X
(cid:90)

where n denotes the local normal direction. Using the IB method, the Eulerian

forcing term fe can be calculated by

fe (x, t) =

Fe (s, t) δ (x − X (s, t)) ds,

(4.5)

where δ denotes the Dirac delta function that maps between the Eulerian (Ωf)

ΩL

and Lagrangian (ΩL) domain. Likewise,

the Lagrangian velocity can be

interpolated from the nearby Eulerian grids via

∂X
∂t

= U (s, t) =

u (x, t) δ (X (s, t) − x) dx.

(4.6)

(cid:90)

ΩE

66

For nondimensionalization, we choose the wavelength λ = 2πk−1 the
length scale, wave speed ω/k the velocity scale, 2νkBT the polymer stress

scale, which lead to the dimensionless equations in the LC phase after proper

re-scaling:

∇ · u = 0, ∇p − ∆u = Er∇ · τp + fe,

− ζ (D · D − S : D) + βE : S,

(cid:18)

(cid:19)

(4.7)

(4.8)

(cid:18)

(cid:19)

D − I
2

τp =
∇
D +2E : S =

+

∆D.

ζ
P e

(D · D − D : S) − 1
P e

Here we deﬁne two P´eclet numbers, P e = ω
8πdr

D − I
2
(∼ O(1) − O(10)) and
dtk2 (∼ O(10) − O(100)), which represent the ratio of the time scales
P et = 2πω
for rod’s rotation and transport over that of undulation (i.e., ω−1), respectively.
Especially P e characterizes the time evolution of LC’s orientation ﬁeld

1
P et

(4.9)

(so-called nematic elasticity).

In comparison, we only consider minimal or

negligible translational diffusion effect, and choose P et to be at least one order

of magnitude higher than P e. The LC-ﬂuid coupling is characterized by the

Ericksen number Er = 4πνkBT
µf ω

([91, 85]). Higher (lower) Er values

And β = β0ω

correspond to the stronger (weaker) ﬂow modiﬁcation on the LC conﬁguration.
8πdr(ν(cid:96)3)2 is an effective crowdness coefﬁcient. Equations (4.3)-(4.4)
and the two mapping functions (4.5) and (4.6) retain the same mathematical
forms, with the dimensionless parameters σs = σ0
s 2πk
ω2

b k3
σb = σ0
2πω2.

,

We ﬁx the wavenumber k = 2π, and set the dimensionless swimmer length

L be the same as the wavelength, i.e., L = Lsk/2π = 1, and the amplitude

A = A0k = 0.05. We choose a large stretching stiffness σs = O(103) to

enforce inextensibility, but vary the bending stiffness σb between 0.005 and 0.5.

67

We then use a pseudo-spectral method to solve the above equations.

4.4 Numerical method

In the solid phase, the elastic body force Fe(Fx, Fy) is calculated explicitly at

each time step. Following the IB approaches in Refs [92, 93], we discretize the

variational derivative of the bending energy functional, and derive the following



∆s3

− 1

∆s

∆s

Fy,i = σs
∆s

(cid:17)(cid:16) Xi+1−Xi
(cid:17)(cid:16) Yi+1−Yi

|Xi+1−Xi|

|Xi+1−Xi|

Fx,i = σs
∆s
− σb

(cid:17)(cid:16) Xi−Xi−1
(cid:17)(cid:16) Yi−Yi−1

|Xi−Xi−1|

component form for the ith node at the position Xi (Xi, Yi):

(cid:17) −(cid:16)|Xi−Xi−1|
(cid:17) −(cid:16)|Xi−Xi−1|

(cid:104)(cid:16)|Xi+1−Xi|
(cid:0)(κ − κ0)i−1 (Yi−1 − Yi−2) − (κ − κ0)i (Yi+1 − Yi−1) + (κ − κ0)i+1 (Yi+2 − Yi+1)(cid:1)
(cid:104)(cid:16)|Xi+1−Xi|
(cid:0)(κ − κ0)i−1 (Xi−1 − Xi−2) − (κ − κ0)i (Xi+1 − Xi−1) + (κ − κ0)i+1 (Xi+2 − Xi+1)(cid:1)
(4.10)
where the curvature κ is deﬁned as κi = ˆez · (Xi+1 − Xi) × (Xi − Xi−1)/∆s3,
and ∆s = L/Ns. When mapping between the Lagrangian and Eulerian (x)

(cid:17)(cid:105)
(cid:17)(cid:105)

|Xi−Xi−1|

+ σb
∆s3

∆s

∆s

− 1

− 1

− 1

meshes, we use a four-point smoothed Dirac delta function φ(r) [92, 94]:

(cid:18)x − X

(cid:19)

(cid:18)y − Y

(cid:19)

h

φ

h

δ(x − X) =

1
h2 φ

(4.11)

(4.12)

where delta function φ(r) is constructed as



0,

1
8

1
8

φ(r) =

(cid:16)
(cid:16)

5 − 2|r| −(cid:112)−7 + 12|r| − 4r2
(cid:17)
3 − 2|r| +(cid:112)1 + 4|r| − 4r2

,

(cid:17)

,

|r| ≥ 2,
2 ≥ |r| ≥ 1,
1 ≥ |r| ≥ 0.

In the LC phase, we ﬁrst solve the constitutive equations for the velocity u

and the second-moment tensor D. Here we use the method by Vaithianathan

68

Figure 4.1: Convergence tests of a parallel moving swimmer with the time-
dependent C.O.M. velocity U(Ux, Uy) by varying (a) domain size (∆t = 6.25×
10−5, h = 1/128), (b) Eulerian grid width (Ωf = 2 × 2), and (c) Lagrangian
grid spacing (Ωf = 2 × 2, h = 1/128). These parameters are ﬁxed: P e = 10,
P et = 100, Er = 1, β = 0.005, ζ = 8, L = 1, A = 0.05, k = 2π, ω = 2π,
σb = 0.5, σs = 1500.

and Collins [95] to enforce the symmetric positive deﬁniteness of D whose
diagonal components are bounded, i.e., tr(D) = D11 + D22 = 1 where 0 ≤
D11, D22 ≤ 1. A Cholesky decomposition



 l11

l21

0 l22

 =

 l2

11

l11l21

21 + l2
l2
22

l11l21

 (4.13)

l21

l22

D = L · LT =

 l11 0
components of L again via

is used, where L is a lower triangular matrix. Then we transform the

l11 = 1

π arctan η + 1
2,

l22 = 1

π arctan γ + 1
2,

l21 = 2

π arctan χ.

(4.14)

69

timevelocity component234-0.6-0.4-0.200.20.40.6Ux, Ωf=2×2Uy, Ωf=2×2Ux, Ωf=4×4Uy, Ωf=4×4(a)time234-0.6-0.4-0.200.20.40.6Ux, h=1/64, Δt=2.5×10-4Uy, h=1/64, Δt=2.5×10-4Ux, h=1/128, Δt=6.25×10-5Uy, h=1/128, Δt=6.25×10-5(b)time234-0.6-0.4-0.200.20.40.6Ux, Ns=32, Δt=2.5×10-4Uy, Ns=32, Δt=2.5×10-4Ux, Ns=128, Δt=2×10-5Uy, Ns=128, Δt=2×10-5(c)Substituting Eqs. (4.13-4.14) into the evolution equation of D, we get

∂y = π(cid:0)1 + η2(cid:1) l11

∂η

∂ux
∂x + l21

∂ux

∂y − 1
∆D11 − 1

2P e

2l11

+ 1
P et

1
2l11

l11 − 1

2l11

(cid:16)

K11

(cid:17)

(cid:17)



(cid:17)





∂η
∂t + ux

∂η
∂x + uy

= π
2

∂χ
∂t + ux

∂χ
∂y

∂χ
∂x + uy

(cid:0)1 + χ2(cid:1) l11


∂γ
∂x + uy

∂γ
∂y

l22

∂γ
∂t + ux

= π(cid:0)1 + γ2(cid:1)

(cid:16)

(cid:16) l2

∂uy

∂x + l2

− 1

P et

l21
2l2
11

∂ux

22
l11

∂y − l21
∆D11 + 1
P et

∂ux

2P e

∂x − 1
∆D12 + l21
2l2
11

1
l11

l21 + l21
2l2
11

K11 − 1

l11

K12

− l22 + 1

2l22

21
2l22l2
11

∂uy

∂y − l22l21

l11

2P e

∂ux
∂y + 1
∆D11 − 1

P et

l2
21
2l22l2
11

+ 1
P et
− l2

21
2l22l2
11

K11 + l21
l22l11

K12 − 1

2l22

K22

l21

l22l11

∆D12 + 1
P et

1
2l22

∆D22

(4.15)
P e(D · D − D : S). After solving the D ﬁeld, we
where K = 2E : S − ζ
evaluate the polymer stress τp, and then directly compute the 2D velocity ﬁeld

u = (ux, uy) in the Fourier space as:

(cid:16)

(cid:17) · ˜ftotal

I − ˆξ ˆξ

(4.16)

˜u =

1
ξ2

where ftotal = Er∇·τp +fe, ˆξ is the wavevector with |ξ| = ξ ( ˆξ = ξ/ξ). All the
equations in the LC phase is solved using a pseudo-spectral method, with the
time integration being evaluated using the 4−th order Runge-Kutta method.For
all simulations performed in this study, we choose the computation parameters
as: Ωf = 2 × 2, h = 1
128, ∆t = 10−5 − 10−4, P e = 0.1 − 10, P et = 10 − 100,
Er = 1 − 10, ζ = 0 − 10, β = 0.0005 − 0.005, Ns = 32, L = 1, A = 0.05,
k = 2π, ω = 2π, σb = 0.005 − 0.5, σs = 3000σb.

70

As the examples of time-dependent velocity components shown in Fig. 4.1(a-

c), we have veriﬁed the numerical results by performing convergence tests for

a parallel moving swimmer. As shown in Fig. 4.2, we have performed another

benchmark study for the undulatory swimming motions in an Oldroyd-B (OB)

ﬂuid where the dimensionless Deborah (De) number, playing the similar role as

P e in the LC cases, is deﬁned as the wave frequency by the OB ﬂuid relaxation

time. We measured the mean C.O.M. swimming speed UOB of the swimmer

that has a long body length L = 4, and compared the speed ratio with the

numerical data by Salazar et al. [96] and the asymptotic results for Taylor’s

swimming sheet by Lauga [97]

UOB
UN

=

1 + ηs

ηs+ηp

De2

1 + De2

,

(4.17)

where ηs and ηp respectively represent the solvent and polymer contribution to

the viscosity. The Newtonian speed UN can be derived as

UN =

1
2

(

ω
k

)(Ak)2 + O (Ak)4

(4.18)

where k = 2π is the wavenumber. We ﬁnd our results agree well with the

previous studies.

In the third test, we have examined the impact of varying bending stiffness

σb when choosing high values. As shown in Fig. 4.3 for the time-dependent

velocity, only small differences are seen when σb goes beyond 0.5. In these

scenarios, the swimmer can quickly respond to the imposed target curvature,

and well follow the traveling-wave actuation. Hence, in the main text, we

choose σb = 0.5 for the cases of a stiff swimmer.

71

Figure 4.2: Time-averaged C.O.M. speed UOB for undulatory swimming
motion in an Oldroyd-B ﬂuid. These parameters are ﬁxed: Ωf = 8 × 8,
h = 1/64, Ns = 64, ηp/ηs = 1/2, L = 4, De = 1, and ∆t = 2.5 × 10−4.

4.5 Asymptotic analysis

In the moving frame of the swimmer, we consider the vertical displacement of

an inﬁnitely-long wavy sheet with the described traveling-wave motion
y(x, t) = A sin(kx − ωt). When choosing 1/k the length scale, 1/ω the time
scale, and ω/k the velocity scale, the dimensionless form can be written as
y(x, t) = ε sin(x − t). Here we assume a small amplitude ε = Ak (cid:28) 1.
Following Ref[97], we adopt a streamfunction ϕ(x, y, t), and hence the
velocity components can be computed as ux = ∂ϕ/∂y, uy = −∂ϕ/∂x with
the incompressibility condition being satisﬁed. The boundary conditions for

ϕ(x, y, t) arise from conditions at inﬁnity and on the undulatory sheet with a
steady speed −ULCˆex. Then the far-ﬁeld condition at y = ∞ reads

∇ϕ|(x,∞) = ULCˆey.

72

(4.19)

Amplitude AUOB00.010.020.030.040.050.0600.010.020.030.040.050.06Lauga  (analytical)Salazar et al. (numerical)present (numerical)Figure 4.3: The time-dependent C.O.M. velocity U(Ux, Uy) for a parallel
moving swimmer when choosing different values of bending stiffness σb. These
parameters are ﬁxed: These parameters are ﬁxed: Ωf = 2×2, ∆t = 6.25×10−4,
h = 1/128, P e = 10, P et = 100, Er = 1, β = 0.005, ζ = 8, L = 1, Ns = 32,
A = 0.05, k = 2π, ω = 2π, σs = 1500.

On the swimming sheet, the no-slip velocity condition is imposed as

∇ϕ|(x,ε sin(x−t)) = ε cos(x − t)ˆex.

Recalling the forced Stokes equation

∇p = ∆u + Er∇ · τp,

where polymer stress, when neglecting β (β ≈ 0), is given as

(cid:18)

(cid:19)

And the D evolution equation in the LC phase

τp =

D − I
2

− ζ(D · D − D : S).

 ∂D
∂t + u · ∇D −(cid:0)D · ∇u + ∇uT · D(cid:1) + 2E : S

P e(D · D − D : S)

= − 1

P e

(cid:0)D − I

(cid:1) + ζ

2

73

(4.20)

(4.21)

(4.22)

(4.23)

timevelocity component23456-0.8-0.6-0.4-0.200.20.4Ux, σb=0.4Uy, σb=0.4Ux, σb=0.5Uy, σb=0.5Ux, σb=0.6Uy, σb=0.6When applying the curl on the both sides of Eq. (4.21), we have

∇ × (∇ · τp) =

∇4ϕ

1
Er

Next, we expand all the variables with ε to the second order, i.e.,

ϕ = εϕ(1) + ε2ϕ(2) + O(cid:0)ε3(cid:1) ,

τ = τ (0) + ετ (1) + ε2τ (2) + O(ε3),

D = D(0) + εD(1) + ε2D(2) + O(ε3),

LC + O(cid:0)ε3(cid:1) .

ULC = εU (1)

LC + ε2U (2)

After some manipulations, we can derive the following governing equations:

(4.24)

(4.25)

(4.26)

(4.27)

(4.28)

(4.29)

(4.30)

∂D(0)

0th−order:

D(0) · ∇u(0) + ∇u(0)T · D(0)(cid:17)
∂t + u(0) · ∇D(0) −(cid:16)
(cid:1) + ζ
(cid:0)D(0) − I
(cid:1) − ζ(D(0) · D(0) − D(0) : S(0)).
p =(cid:0)D(0) − I

P e(D(0) · D(0) − D(0) : S(0)),

− 1

τ (0)

P e

2

2

+ 2E(0) : S(0) =

∂D(1)

1st−order:
∂t + u(0) · ∇D(1) + D(1) · ∇D(0)

D(0) · ∇u(1) + D(1) · ∇u(0) + ∇u(0)T · D(1) + ∇u(1)T · D(0)(cid:17)
−(cid:16)
+2(cid:0)E(0) : S(1) + E(1) : S(0)(cid:1)

(cid:0)D(0) · D(1) + D(1) · D(0) − D(0) : S(1) − D(1) : S(0)(cid:1) ,

= − 1

P eD(1) + ζ

P e

(4.31)
p = D(1) − ζ(D(0) · D(1) + D(1) · D(0) − D(0) : S(1) − D(1) : S(0)). (4.32)
τ (1)

74

+ 2(cid:0)E(0) : S(2) + E(1) : S(1) + E(2) : S(0)(cid:1)

P e

∂D(2)

P eD(2) + ζ

P e

= − 1
− ζ

2nd−order:
∂t + u(0) · ∇D(2) + u(1) · ∇D(1) + u(2) · ∇D(0)

(cid:0)D(0) · D(2) + D(1) · D(1) + D(2) · D(0)(cid:1)

D(0) · ∇u(2) + D(1) · ∇u(1) + D(2) · ∇u(0) + ∇u(0)T · D(2)(cid:17)
−(cid:16)
−(cid:16)∇u(1)T · D(1) + ∇u(2)T · D(0)(cid:17)
(cid:0)D(0) : S(2) + D(1) : S(1) + D(2) : S(0)(cid:1) ,
−ζ(cid:0)D(0) · D(2) + D(1) · D(1) + D(2) · D(0) − D(0) : S(2)(cid:1)
−ζ(cid:0)D(1) : S(1) + D(2) : S(0)(cid:1) .
∇ϕ(1)(cid:12)(cid:12)(cid:12)(x,∞)
∇ϕ(1)(cid:12)(cid:12)(cid:12)(x,0)

= cos(x − t)ˆex

τ (2)
p = D(2)

= U (1)

LCˆey

At the 0th−order, we can solve for homogeneous solutions. Now the boundary
conditions at the 1st−order become:

(4.33)

(4.34)

(4.35)

(4.36)

Note that in the above, instead of being satisﬁed exactly along the wavy body,
the no-slip boundary condition is projected onto the x−axis, i.e., at y = 0. And
at the 2nd−order, they take the form

∇ϕ(2)(cid:12)(cid:12)(cid:12)(x,∞)

∇ϕ(2)(cid:12)(cid:12)(cid:12)(x,0)

= U (2)

LCˆey

(cid:18)∂ϕ(1)

(cid:19)

|(x,0).

(4.37)

(4.38)

∂y

= − sin(x − t)∇

Isotropic cases

4.5.1
In the isotropic regime where 0 ≤ ζ < ζc, we close the 4th−moment tensor S
by expanding the PDF near the isotropy (i.e., Ψ → 1
2π in the limit of random

75

conﬁguration) via

which leads to

(cid:18)

pp − I
2

(cid:19)

: D,

Ψ =

1
2π

+

2
π

Sijkl = − 1
+ 1

24 (δijδkl + δikδjl + δilδjk)

6 (δijDkl + δikDjl + δilDjk + δjkDil + δjlDik + δklDij) .

p = 0, and D(0) = diag(cid:0) 1

2, 1

2

At the 0th−order, we obtain u(0) = 0, τ (0)
1st−order, we substituting D(0) into Eq. (4.32) and obtain

p = (1 − ζ/4)D(1).
τ (1)

When substitute Eq. (4.41) into Eq. (4.31), we can get

1

1 − ζ/4

∂τ (1)
p
∂t

= − 1
P e

τ (1)
p +

1
2

E(1).

(4.39)

(4.40)

(cid:1). At the

(4.41)

(4.42)

Then we ﬁrst take the divergence and then apply the curl on the both sides of

4

Er(4 − ζ)

∂
∂t

+ (

1

ErP e

+

1
4

)

∇4ϕ(1) = 0

(4.43)

Given the boundary conditions in Eqs. (4.35-4.36), it is straightforward to solve
for the 1st−order solutions as:

ϕ(1)(x, y, t) = (1 + y)e−y sin(x − t),

U (1)

LC = 0.

(4.44)

(4.45)

Next, at the 2nd−order, we ﬁrst evaluate Eq.(4.34) to obtain

p = (1 − ζ
τ (2)
4

)D(2) − ζ(D(1) · D(1) − D(1) : S(1))

(4.46)

The Eq. (4.33) can be rewritten as

1−ζ/4

1

−(cid:16)

p

∂τ (2)
∂t + 1

P eτ (2)

D(1) · ∇u(1) + ∇u(1)T · D(1)(cid:17) − ζ

2E(2) = −2E(1) : S(1) − u(1) · ∇D(1)
∂(D(1)·D(1)−D(1):S(1))

p − 1

,

1−ζ/4

∂t

76

Eq. (4.42) to obtain(cid:20)

(cid:21)

which allows us to derive the equation for ϕ(2) as

[

4

Er(4−ζ)

∂

−(cid:16)

∂t + ( 1

ErP e + 1

D(1) · ∇u(1) + ∇u(1)T · D(1)(cid:17) − ζ

4)]∇4ϕ(2) = ∇ × {∇ · [−2E(1) : S(1) − u(1) · ∇D(1)
]} (4.47)

∂(D(1)·D(1)−D(1):S(1))

1−ζ/4

∂t

Note that to solve for mean speed ULC, it doesn’t require that we obtain the time-

dependent solution of Eq. (4.47). Instead, we can perform the time averaging
(“(cid:104)(cid:105)”) as the following:

(

1

ErP e

+

1
4

)∇4(cid:104)ϕ(2)(cid:105) =

4P e2

(4 − ζ)2 + 16P e2 (8y2 − 24y + 8)e−2y

(4.48)

with boundary conditions in Eqs. (4.37-4.38) becoming

∂(cid:104)ϕ(2)(cid:105)

(∞) = U (2)
LC,

∂y
∂(cid:104)ϕ(2)(cid:105)

We can obtain

∂y

(0) = 1
2.

(cid:104)ϕ(2)(cid:105)(x, y) = a0(−y2 +

d
dy

)e−2y + (

1
2

1
2

− 1
2

a0)

(4.49)

(4.50)

(4.51)

where a0 = ErP e
4+ErP e
at the 2nd−order as

16P e2

(4−ζ)2+16P e2. So we can evaluate the mean swimming speed

ULC
UN

=

U (2)
LC
U (2)
N

= 1 −

16ErP e3

(4 + ErP e) [(4 − ζ)2 + 16P e2]

(4.52)

where UN = U (2)

N ε2 = 1/2ε2 is the swimming velocity in a Newtonian ﬂuid.
Note that in Doi’s theory, the LC polymer contribution to the zero-shear-rate

viscosity can be effectively deﬁned as [98]

µp
µf

= α(S)ErP e

77

(4.53)

where S =(cid:112)(2D : D − 1)/2 is the order parameter, α (S) is a concentration

coefﬁcient, and µp represents the polymer contribution to the viscosity. Hence,

when ζ = 0, we estimate S = 0 and α(S) = 1, leading to

(cid:16) µf

(cid:17)

µf +µp/4
1 + P e2

P e2

.

(4.54)

1 +

ULC
UN

=

4.5.2 Nearly-aligned cases

In the nematic regime, we adopt a classical quadratic closure to approximate S

[29] as

S = DD,

(4.55)

(cid:16)

ζ > ζc.

corresponding to the strong alignment conﬁguration with large ζ values when
the 0th−order, we solve the steady-state equation for
At
(cid:20)
11 , 1 − D(0)
(cid:16)
D(0)

D(0) = diag

(cid:26)

(cid:19)

(cid:18)

(cid:17)

via

11

2 −

D(0)
11

2

D(0)
11

+

1 − D(0)

ζ

= 0,

(4.56)

−

11 − 1
D(0)
2

which yields two equilibrium solution. One takes the form

D(0)

11 =

>

1
2

.

(4.57)

This solution represents that the LC’s principal alignment direction is aligned
with the x−axis. When ﬁxing the swimmer’s moving direction to be
horizontal, the above solution corresponds to the parallel swimming cases.

Another solution is

D(0)

11 =

(4.58)
meaning the principal alignment direction of LC is along the y−axis. Hence,
this solution corresponds to the swimming direction is perpendicular to the

2ζ

<

,

1
2

(cid:17)2(cid:21)
(cid:27)
ζ +(cid:112)ζ 2 − 2ζ

D(0)
11

11

2ζ

ζ −(cid:112)ζ 2 − 2ζ

LC’s alignment direction.

78

 1 0

0 −1

 .

At the 1st−order, it is straightforward to show that

Substituting Eq. (4.59) into Eq. (4.31), we have

11

p = (ζ − 2)D(1)
τ (1)
(cid:17)
(cid:16)
(cid:16)
(cid:17)
11 − 1
D(0)
11 − 1
D(0)

∂τ (1)
∂t + 4
p,11
∂t − 4
∂τ (1)
p,22

 1

ζ−2

D(0)

11 E(1)

11 + 1

P eτ (1)

p,11 = 0

P eτ (1)
Applying the same trick in Eq. (4.24), we can further derive

11 + 1

11 E(1)

D(0)

1
ζ−2

p,22 = 0

(

1
ζ − 2

1
Er

∂
∂t

+

1

ErP e

)∇4ϕ(1) +

4
ζ

∂4ϕ(1)
∂x2∂y2 = 0

To facilitate further analytical manipulations, we consider the strong alignment
cases when ζ is large (ζ (cid:29) 1), which allows us to neglect the second term, and
obtain the same 1st − order solution as Eq. (4.44).

At the 2nd−order, using the equilibrium solutions of D(0), we can show that

p = [(ζ − 2)D(2)
τ (2)

11 + ζ(3D(1)

11

+ D(1)
12

2

)(2D(0)

11 − 1)]

 1 0

0 −1



2

1 0

 0 1

2(cid:17)(cid:16)
2(cid:17)(cid:16)

+2ζD(1)

11 D(1)

12 (2D(0)

11 − 1)

which lead to

 1

ζ−2

1
ζ−2

∂
∂t

∂
∂t

(cid:104)
(cid:104)

p,11 − ζ
τ (2)
τ (2)
p,22 + ζ

(cid:16)
(cid:16)

2

2

3D(1)
11

+ D(1)
12

3D(1)
11

+ D(1)
12

(4.59)

(4.60)

(4.61)

,

(4.62)

(cid:17)(cid:105)
(cid:17)(cid:105)
11 − 1
11 − 1

2D(0)

2D(0)

+ G11 + 1

P eτ (2)

p,11 = 0

+ G22 + 1

P eτ (2)

p,22 = 0

(4.63)

79

where

G = u(1) · ∇D(1)

D(0) · ∇u(2) + D(1) · ∇u(1) + ∇u(1)T · D(1) + ∇u(2)T · D(0)(cid:17)
−(cid:16)
+2(cid:0)E(1) : S(1) + E(2) : S(0)(cid:1) .
(cid:34)
(cid:16) 1

(cid:17)(cid:105)
ErP e∇4ϕ(2) + ∂2(G22−G11)

(cid:17)
11 −1)

∇4ϕ(2) + 2

In the end, we can derive the linear equation for ϕ(2) as

(cid:17)(cid:104)

+D(1)
12
∂x∂y

∂2(cid:16)

(cid:16)

(cid:35)

(2D(0)

2(cid:17)

3D(1)
11

1
ζ−2

∂
∂t

1
Er

∂x∂y

(cid:16)

ζ

(cid:17)(cid:16) ∂2
∂x2 − ∂2

∂y2

2ζD(1)

11 D(1)

12

2D(0)

+ 1
11 − 1

=

ζ−2

∂

∂t + 1

P e

2

(4.64)

Similar to the isotropic case, here we take the time averaging of Eq. (4.64) to

obtain

1

(cid:16)
ErP e∇4(cid:104)ϕ(2)(cid:105) =
2(ζ−2)
2D(0)
P e2+(ζ−2)2

(cid:17)(cid:104)
11 − 1
(cid:16)

(cid:16)

4y2 +

8D(0)

(cid:17)
11 − 12
(cid:17)

y −(cid:16)
(cid:32)

e−2y +

When applying the boundary conditions (4.37-4.38), we obtain

d
dy

(cid:104)ϕ(2)(cid:105)(x, y) = −a1
2

y2 + 2D(0)

11 y + D(0)

11

1
2

+

a1D(0)
11

2

(cid:17)(cid:105)

e−2y.

8D(0)

11 − 6

(cid:33)

(4.65)

(4.66)

2(ζ−2)

where a1 = ErP e
speed ratio at the 2nd−order as

(ζ−2)2+P e2 (2D(0)

11 − 1). So we can eventually solve for the

ULC
UN

= 1 +

4ErP e(ζ − 2)
(ζ − 2)2 + P e2 (D(0)

11 − 1
2

)D(0)
11

(4.67)

4.6 Results and Discussion

Without the swimmer, the pure LC phase is described by a homogeneous

distribution of D0 =(cid:82)

p ppΨ0dp, i.e., obtained from taking the conﬁgurational

80

Figure 4.4: I-N phase transition led anisotropic swimming. (a) δ as a function of
ζ, with the bifurcation occurring at ζc = 4. (b)Constant and apolar distribution
of Ψ0 (φ) at ζ = 0 and 8. (c) ULC/UN as a function of ζ for stiff (σb = 0.5)
and soft (σb = 0.005) swimmer when choosing P e = 1 and Er = 10, for the
parallel (denoted by “(cid:107)”) and and perpendicular (denoted by “⊥”) motions.

average of the equilibrium solution of the original Doi’s kinetic model:

(cid:82) 2π

Ψ0 =

exp [δ cos (2φ)]

0 dφ(cid:48) exp [δ cos (2φ(cid:48))]

,

δ =

ζ
4

(cid:82) 2π
(cid:82) 2π
0 dφ(cid:48) cos (2φ(cid:48)) exp [δ cos (2φ(cid:48))]

0 dφ(cid:48) exp [δ cos (2φ(cid:48))]

.

(4.68)

As shown in Fig. 1(a), δ (ζ) admits a bifurcation at ζc = 4, and has two

solution branches for the (I)sotropic and (N)ematic states. In Fig. 4.4(b), the

symmetric form of Ψ0 (φ) suggests that the rod’s primary (or principal)
orientation direction is aligned with the x−axis when φ = 0 and φ in the
nematic regime (ζ ≥ ζc), compared to the constant solution in the isotropic
regime (0 ≤ ζ < ζc). At each time step in simulation, we reconstruct Ψ0 from
p ppppΨ0dp following a
quasi-equilibrium ansatz that D and S co-align in principal directions in the

the D ﬁeld, and then approximate S (x, t) ≈ (cid:82)

eigenspace. This method is also referred to as the Bingham closure ([33, 76]),

which has proven to be accurate in PDF reconstruction, and has been widely

used in simulating both equilibrium and non-equilibrium LCs ([76, 99, 100]).

When there are microswimmers moving in a LC, curiously, one may ask

81

ζULC/UN02468100.60.91.21.51.8//, stiff⊥, stiff//, soft⊥, soft(c)ζδ0246801234nematicζcisotropic(a)φΨ000.20.40.60.81  0                       π                    2π ζ = 8ζ = 0(b)how do they mechanically respond to this anisotropic polymer ﬂuid? To

answer this question, we consider undulatory swimming motions either
parallel with or perpendicular to LC alignment direction (i.e., x−direction),
for both the “stiff” (σb = 0.5) and “soft” (σb = 0.005) swimmers. When we

impose traveling waves that are moving rightward, the swimmer produces

leftward undulatory motions. We measure the mean C.O.M. velocity ULC by

performing time averaging over ﬁve undulation periods after the swimming

speed reaches quasi-steady states, and compare with that for the Newtonian

speed UN. As shown in Fig. 4.4(c), the speed ratio ULC/UN as a function of ζ

clearly show non-Newtonian behaviors, especially in the nematic regime

where similar bifurcations occur, corresponding the motions that are parallel
(denoted by symbol “//”) and perpendicular (denoted by symbol “⊥”) to the
x−axis. More speciﬁcally, we observe enhanced swimming speeds when the
swimmer moves parallel in the nematic regime; while retarded motions are

seen in the isotropic regime, and in the most nematic regime during

perpendicular motions. Also, while qualitatively similar behaviors are

consistently observed for the stiff and soft cases in the parameter space, it is

seen that

the stiff swimmer can generally achieve faster mean speeds,

exhibiting more signiﬁcant speed gain and loss.

Figures 4.5(a-c) show the variations of ULC/UN with respect to P e when

ﬁxing Er. When P e is small, which corresponds to a slow undulation (or

equivalently, quick rotational diffusion of LC molecules), ULC/UN is close to
1, suggesting Newtonian-like behaviors. This is because when P e → 0, the
Brownian diffusion in the LC phase dominates over convection, which not only

82

Figure 4.5: Speed ratio ULC/UN as a function of P e at Er = 1 and 10 obtained
from numerical simulation (a-c) and asymptotic solutions (d-f).

drives the LC structure towards equilibrium states (i.e., the R.H.S. of Eq. (4.9)

is approximately 0.) but also produces smaller and smaller τp. On the other

hand, at a high P e where convection dominates during fast acutation, the LC’s

conﬁguration is modulated by ﬂow via constraining stress, i.e., E : S, due

to rod’s rigidity, which leads to more complicated behaviors. In the isotropic

regime, ULC/UN monotonically decrease. Nevertheless, in the nematic regime,

ULC/UN varies non-monotonically, with the maximum (minimum) occurring
at small P e (∼ O(1)) during the parallel (perpendicular) swimming motions.
To understand the numerical observations, we build an idealized linear

model of a one-dimensional,

inﬁnitely-long sheet prescribed with small

amplitude traveling waves, i.e., the so-called Taylor’s swimming sheet ([101]).

Its vertical displacement can be described in the moving coordinate as
y(x, t) = ε sin(x − t) with amplitude ε (cid:28) 1. Then we simplify the above
governing equations by employing a stream function ϕ to replace the

83

ULC/UN02468100.80.91Er=1, stiffEr=10, stiffEr=1, softEr=10, soft(a)ζ= 2024681011.21.41.6Er=1, stiffEr=10, stiffEr=1, softEr=10, soft(b)ζ= 8, //02468100.70.80.91Er=1, stiffEr=10, stiffEr=1, softEr=10, soft(c)ζ= 8, ⊥Pe024681012345678910Er=1Er=10(e)ζ= 8, //Pe02468100.40.60.81Er=1Er=10(f)ζ= 8, ⊥PeULC/UN024681000.20.40.60.81Er=1Er=10(d)ζ= 2incompressible ﬂuid velocity via u = ∇ × (ϕˆez) (ˆez is the unit vector
denoting the out-of-place direction). By following the solution procedure that

is almost identical to [102] and [85, 88], we impose the no-slip condition on

the wavy sheet, and perform asymptotic analyses by expanding all

the

variables in terms of ε in both the isotropic and nematic regime. Moreover, we

neglect the crowdness effect (i.e., β = 0) in τp and the translational Brownian
t → 0), and adopt different closure methods for S in the
diffusion (i.e., P e−1
some
isotropic and nematic regime to facilitate derivation.

After

manipulations (see more derivation details in the supplemental material), we

ﬁnd directional motions only occur in the second order, i.e., ULC = U (2)

LCε2,

and the speed ratio can be calculated as:

(cid:18)

(cid:19)

,

16ErP e3

= 1 −

ULC
(4 + ErP e) [(4 − ζ)2 + 16P e2]
UN
4ErP e(ζ − 2)
ULC
D(0)
(ζ − 2)2 + P e2
= 1 +
11 ,
UN
√
(cid:16)
11 = ζ±
ζ 2−2ζ
2ζ
11 , 1 − D(0)
D(0)

11 − 1
D(0)
2

(cid:17)

11

where D(0)

D(0) = diag

(Isotropic)

(Nematic)

(4.69)

(4.70)

is the diagonal component of the equilibrium solution

, and the plus (minus) sign corresponds to parallel

(perpendicular) swimming motion with respect to the alignment direction.

As shown in Figs. 4.5(d-f), the above asymptotic solutions predict the

qualitatively similar behaviors of ULC/UN as the numerical results, and reveal

different scaling behaviors. In the isotropic regime, the monotonic decay as a

function of P e at a given Er is reminiscent of the behaviors in the viscoelastic

(e.g., Oldroyd-B) ﬂuids ([102, 103]). This can be witnessed by introducing an

effective polymer viscosity µp ([29, 89]) as µp
µf

= α (S)ErP e, where ErP e

deﬁnes the effective polymer contribution to the zero-shear-rate viscosity, and

84

Figure 4.6: The instantaneous (a-d) and time-averaged (e-h) velocity ﬁeld
around a “stiff” swimmer in the body-ﬁxed frame when choosing σb = 0.5,
P e = 1, and Er = 10: (a,e) Newtonian; (b,f) Isotropic, ζ = 2; (c,g) Nematic
and parallel, ζ = 8; (d,h) Nematic and perpendicular, ζ = 8. The vector ﬁeld
and the background color correspond to the velocity u and its magnitude |u|.

the estimated coefﬁcient α characterizes the dependence on the order
parameter S. When ζ → 0, we estimate S ≈ 0 and α(S) ≈ 1, and hence can
rewrite Eq. (4.69) as

(cid:16) µf

(cid:17)

1 +

ULC
UN

=

µf +µp/4
1 + P e2

P e2

,

(ζ → 0)

(4.71)

resembling the asymptotic solution by [102] for Taylor’s swimming sheet in

viscoelastic ﬂuids. In the nematic regime, from Eq. (4.70), apparently the speed

1

gain and loss are seen to be directly determined by the LC’s horizontal (D(0)
11 >
→
2) and vertical (D(0)
1 in the limit of P e → ∞, suggesting diminishing non-Newtonian behaviors
for very fast undulations in sharply-aligned LCs.

2) alignment direction. Moreover, in both cases, ULC
UN

11 < 1

Next, we take a closer look at the ﬂow patterns to uncover the nonlinear

effect due to the ﬁnite length of the swimmer. As shown in Figs. 4.6(a-d), we

85

Figure 4.7: The time-averaged director and polymer force ﬁelds around a
(a-c): nematic director (cid:104)m(cid:105)
stiff swimmer corresponding to Figs. 4.6(f-h).
superposed on the colormap of the order parameter S. (d-f): polymer force
vector (cid:104)fp(cid:105) superimposed on the force magnitude |(cid:104)fp(cid:105)|.
compare the simulation results for instantaneous ﬂow ﬁelds at the end of each

period, after the time-periodic swimming motions become quasi-steady. While

the isotropic (ζ = 2) case in panel(b) resembles the Newtonian one in panel(a),

the parallel (perpendicular) swimming case in the nematic regime in panel(c)

(panel(d)) has slightly wider (narrower) circulation zones near the two ends.

Indeed, when averaging over three to ﬁve periods, we consistently ﬁnd the

ﬂow patterns shown in Figs. 4.6(e-h) distinctive features. Compared to the

Newtonian case in panel(e), in the isotropic regime, the LC ﬂow in panel(f),

while similar, appears to be slightly weakened. More interestingly, we ﬁnd

that the mean ﬂow in the nematic regime can be extensile (panel(g)) and

contractile (panel(h)), corresponding to the parallel and perpendicular motions,

respectively.

Figures 4.7(a-c) reveal the characteristics of the nematic ﬁeld. Without

86

imposing particular anchoring conditions for D along the wavy body, we ﬁnd
the time-averaged (denoted by “(cid:104)(cid:105)”) director distribution (cid:104)m(cid:105) approximately
follow the global equilibrium ordering, i.e., random (panel(a)) and aligned

(panels(b,c)). But due to the ﬁnite length of the swimmer, the resultant LC

structures lack left-right symmetry. The disturbances on the nematic ﬁeld are

concentrated in the areas around the body’s rear side without propagating

further away from where strong steric interactions (i.e., characterized by large

values of ζ) enforce local alignment. Also, when comparing the two scenarios

in the nematic regime in panel(b) and (c), it appears that parallel motions can

impact larger areas. This is because the vertical oscillations during undulation
are a lot faster (∼ 5 − 10 times) than the resultant horizontal migration, and
can lead to stronger momentum exchange along the y−direction. Hence, it is
the LC molecules
easier for the shear-induced force tends to re-orient

perpendicularly, which leaves larger low-order zones during parallel motions.
Further examination of the time-averaged polymer force, (cid:104)fp(cid:105) = (cid:104)∇ · τp(cid:105), in
Fig. 4.7 reveal more precise propulsion mechanisms. Since the periodic elastic
body force have a zero mean, i.e., (cid:104)fe(cid:105) = 0, from Eq. (4.7), the near-body ﬂuid
ﬂows are driven by (cid:104)fp(cid:105) at a ﬁnite Er. As typically shown in panel(d),
isotropic LCs can only produce small (cid:104)fp(cid:105) without showing clear correlations
with the swimming direction. In comparison, nematic LCs apparently create
much stronger and left-right asymmetrical (cid:104)fp(cid:105) in panel(e) and (f), with the
largest values at the “tail”. And they are clearly seen to be extensile and

contractile, respectively, suggesting to be the driving forces of the ﬂow

patterns in Fig. 4.6(g) and (h). More intriguingly, unlike the isotropic cases,

87

the near-body distributions of (cid:104)fp(cid:105) in the nematic regime highly correlate with
the swimming direction. In panel(e), the projected force vectors at y = 0 all

approximately point leftwards to “push” the swimmer forward; while in panel
(f), the projected force along the x−direction, while weaker, appears to be
overall pointing rightward, and “pull” the swimmer backward. In fact, when

considering more general scenarios when there are disturbances (e.g., small

oscillations) occurring in nematic LCs, our

results suggest anisotropic

propagation of mechanical signals as the induced force transmissions are in

favor of nematic alignment directions.

4.7 Conclusion

This work has built a Q−tensor model to study the nonlinear dynamics of
undulatory swimming motions in LCs in the limit of vanishing Reynolds

numbers. Combining direct numerical simulations and asymptotic analysis,

we have revealed both enhanced and retarded swimming behaviors for a

ﬂexible swimmer moving in LC solutions which may exhibit different

orientation orders. Compared to those top-down macro models of EL or

Landau-de Gennes (LG) type ([32]), ours is derived bottom-up, and hence has

accurate microscopic origins but fewer computation parameters. Moreover, it

avoids to use any orientational derivatives, and hence does not require

applying ad-hoc anchoring conditions to enforce LC molecule alignment on

the boundaries or interfaces, which can signiﬁcantly reduce complexity in

constructing numerical algorithm and facilitate asymptotic analytical analysis.

Our results suggest that besides adopting wavy body undulations to break

88

symmetry, it is possible for a ﬁnite-length ﬂexible swimmer to make use of the

unique near-body ﬂuid dynamics and polymer force generation to adjust the

migration speed. We expect to apply the same methodology to generally study

microorganism’s motion, both individually and collectively,

in anisotropic

complex-ﬂuids environments.

89

5

PUBLICATION

[1] Lin, ZW., Chen, S., and Gao, T., 2021. Q-tensor model for undulatory
swimming in a liquid crystal. submitted to Journal of Fluid Mechanics.

[2] Chen, S., Yan, W. and Gao, T., 2020. Scaling law of Brownian rotation in
dense hard-rod suspensions. Physical Review E, 102(1), p.012608.

[3] Chen, S., Gao, P. and Gao, T., 2018. Dynamics and structure of an apolar
active suspension in an annulus. Journal of Fluid Mechanics, 835, p.393.

[4] Chen, S. and Gao, T,. Transport of an Apolar Active Suspension in a
Rectangular Channel. ﬁnished

90

BIBLIOGRAPHY

91

BIBLIOGRAPHY

[1]

[2]

[3]

S. Chen, P. Gao, and T. Gao. Dynamics and structure of an apolar active
suspension in annulus. J. Fluid Mech., 835:393–405, 2018.

See supplementary information.

S. Ramaswamy. The mechanics and statistics of active matter. Ann. Rev.
Cond. Matt. Phys., 1:323–345, 2010.

[4] M. Marchetti, J. Joanny, S. Ramaswamy, T. Liverpool, J. Prost, M. Rao,
and R. Simha. Hydrodynamics of soft active matter. Rev. Mod. Phys.,
85:1143–1189, 2013.

[5] M. Shelley. The dynamics of microtubule/motor-protein assemblies in

biology and physics. Annu. Rev. Fluid Mech., 48:487–506, 2016.

[6] A. Sokolov, I. Aranson, J. Kessler, and R. Goldstein. Concentration
dependence of the collective dynamics of swimming bacteria. Phys. Rev.
Lett., 98:158102, 2007.

[7]

T. Sanchez, D. Chen, S. DeCamp, M. Heymann, and Z. Dogic.
Spontaneous motion in hierarchically assembled active matter. Nature,
491:431–434, 2012.

[8] G. Henkin, S. J. DeCamp, D. T. N. Chen, T. Sanchez, and Z. Dogic.
Tunable dynamics of microtubule-based active isotropic gels. Phil. Trans.
R. Soc. A, 372(2029), 2014.

[9]

F. Woodhouse and R. Goldstein. Spontaneous circulation of conﬁned
active suspensions. Phys. Rev. Lett., 109:168105, 2012.

[10] M. Ravnik and J. Yeomans. Conﬁned active nematic ﬂow in cylindrical

capillaries. Phys. Rev. Lett., 110:026001, 2013.

[11] H. Wioland, F. Woodhouse, J. Dunkel, J. Kessler, and R. Goldstein.
Conﬁnement stabilizes a bacterial suspension into a spiral vortex. Phys.
Rev. Lett., 110:268102, 2013.

[12] A. Bricard, J. Caussin, N. Desreumaux, O. Dauchot, and D. Bartolo.
Emergence of macroscopic directed motion in populations of motile
colloids. Nature, 503:95–98, 2013.

92

[13] B. Ezhilan and D. Saintillan. Transport of a dilute active suspension in

pressure-driven channel ﬂow. J. Fluid Mech., 777:482–522, 2015.

[14] A.C.H. Tsang and E. Kanso.

Density shock waves in conﬁned

microswimmers. Phys. Rev. Lett., 116:048101, 2016.

[15] Pau Guillamat, Jordi Ignés-Mullol, and Francesc Sagués. Control of
active liquid crystals with a magnetic ﬁeld. Proceedings of the National
Academy of Sciences, 113(20):5498–5502, 2016.

[16] Pau Guillamat, Jordi Ignés-Mullol, Suraj Shankar, M Cristina Marchetti,
and Francesc Sagués. Probing the shear viscosity of an active nematic
ﬁlm. Physical review E, 94(6):060602, 2016.

[17] M. Theillard, R. Alonso-Matilla, and D. Saintillan. Geometric control of

active collective motion. Soft Matter, 13:363–375, 2017.

[18] K. Wu, J. Hishamunda, D. Chen, S. DeCamp, Y. Chang, A. Fernández-
Nieves, S. Fraden, and Z. Dogic. Transition from turbulent to coherent
ﬂows in conﬁned three-dimensional active ﬂuids. Science, 355(6331),
2017.

[19] T. Gao, M. Betterton, A. Jhang, and M. Shelley. Analytical structure,
dynamics, and coarse-graining of a kinetic model of an active ﬂuid. Phys.
Rev. Fluids, 2:093302, 2017.

[20] T. Gao and Z. Li. Self-driven droplet powered by active nematics. Phys.

Rev. Lett., 119:108002, 2017.

[21] D. W. Adams and J. Errington. Bacterial cell division: assembly,
maintenance and disassembly of the z ring. Nat. Rev. Micro., 7:642–653,
2009.

[22] M. Davies Wykes, J. Palacci, T. Adachi, L. Ristroph, X. Zhong, M. Ward,
J. Zhang, and M. Shelley. Dynamic self-assembly of microscale rotors
and swimmers. Soft Matter, page accepted, 2016.

[23] E. Jewell, W. Wang, and T. Malloukl. Catalytically driven assembly
of trisegmented metallic nanorods and polystyrene tracer particles. Soft
Matter, 12:2501–2504, 2016.

[24] D. Saintillan and M. Shelley. Active suspensions and their nonlinear

models. C. R. Phys., 14:497–517, 2013.

[25] H. Wioland, E. Lushi, and R. Goldstein. Directed collective motion of

bacteria under channel conﬁnement. New J. Phys., 18, 2016.

93

[26] D. Saintillan and M. Shelley. Instabilities, pattern formation, and mixing

in active suspensions. Phys. Fluids, 20:123304, 2008.

[27] T. Gao, R. Blackwell, M. Glaser, M. Betterton, and M. Shelley.
Multiscale polar theory of microtubule and motor-protein assemblies.
Phys. Rev. Lett., 114:048101, 2015.

[28] C. Hohenegger and M. Shelley. Dynamics of complex bio-ﬂuids. Oxford

University Press, 2011.

[29] M. Doi and S. Edwards. The theory of polymer dynamics, volume 73.

Oxford University Press, USA, 1988.

[30] B. Ezhilan, M. Shelley, and D. Saintillan.

Instabilities and nonlinear
dynamics of concentrated active suspensions. Phys. Fluids, 25:070607,
2013.

[31] W. Maier and A. Saupe.

Eine einfache molekulare theorie des
nematischen kristallinﬂüssigen zustandes. Zeit. Nat. Teil A, 13:564,
1958.

[32] P. G. DeGennes. The physics of liquid crystals. Oxford University Press,

1974.

[33] Christopher Bingham. An antipodally symmetric distribution on the

sphere. The Annals of Statistics, pages 1201–1225, 1974.

[34] C. Chaubal and L. Leal. A closure approximation for liquid-crystalline
J. Rheol.,

polymer models based on parametric density estimation.
42:177–201, 1998.

[35] T. Gao and H. H. Hu. Deformation of elastic particles in viscous shear

ﬂow. J. Comput. Phys., 228:2132–2151, 2009.

[36] H. H. Hu, M. Y. Zhu, and N. Patankar. Direct numerical simulations of
ﬂuid solid systems using the arbitrary lagrangian eulerian technique. J.
Comput. Phys., 169:427–462, 2001.

[37] Sumesh P Thampi, Amin Doostmohammadi, Ramin Golestanian, and
Intrinsic free energy in active nematics. EPL

Julia M Yeomans.
(Europhysics Letters), 112(2):28004, 2015.

[38] R. Voituriez, J. Joanny, and J. Prost. Spontaneous ﬂow transition in active

polar gels. Europhys. Lett., 69:404–10, 2005.

94

[39] L. Giomi, L. Mahadevan, B. Chakraborty, and M. Hagan. Excitable

patterns in active nematics. Phys. Rev. Lett., 106:218101, 2011.

[40] L Giomi, L Mahadevan, B Chakraborty, and MF Hagan. Banding,
excitability and chaos in active nematic suspensions. Nonlinearity,
25(8):2245, 2012.

[41] L. Giomi, M. Bowick, X. Ma, and M. Marchetti. Defect annihilation and

proliferation in active nematics. Phys. Rev. Lett., 110:228101, 2013.

[42] S. Thampi, R. Golestanian, and J. Yeomans. Instabilities and topological

defects in active nematics. Europhys. Lett., 105:18001, 2014.

[43] E. Lushi, H. Wioland, and R. Goldstein. Fluid ﬂows generated by
swimming bacteria drive self-organization in conﬁned ﬂuid suspensions.
Proceedings of the National Academy of Sciences, 111:9733–9738, 2014.

[44] T. Gao, R. Blackwell, M. Glaser, M. Betterton, and M. Shelley.
Multiscale modeling and simulation of microtubule/motor protein
assemblies. Phys. Rev. E, 92:062709, 2015.

[45] D. Koch and G. Subramanian. Collective hydrodynamics of swimming
living ﬂuids. Annu. Rev. Fluid Mech., 43:637–659,

microorganisms:
2011.

[46] D. Takagi, A. Braunschweig, J. Zhang, and M. Shelley. Dispersion of
self-propelled rods undergoing ﬂuctuation-driven ﬂips. Phys. Rev. Lett.,
110:038301, 2013.

[47] A. Kaiser, A. Peshkov, A. Sokolov, B. ten Hagen, H. Löwen, and
I. Aranson. Transport powered by bacterial turbulence. Phys. Rev. Lett.,
112:158101, 2014.

[48] A. Sokolov, M. Apodaca, B. Grzybowski, and I. Aranson. Swimming
bacteria power microscopic gears. Proc. Natl. Acad. Sci., 107:969–974,
2010.

[49] D. Saintillan and M. Shelley. Instabilities and pattern formation in active
particle suspensions: Kinetic theory and continuum simulations. Phys.
Rev. Lett., 100:178103, 2008.

[50]

J. Jeong and C. Jang. Slow motion of a circular cylinder in a plane
poiseuille ﬂow in a microchannel. Phys. Fluids, 154:123104, 2014.

95

[51] A. Duarte, A. Miranda, and P. Oliveira. Numerical and analytical
modeling of unsteady viscoelastic ﬂows: The start-up and pulsating test
case problems. J. Non-Newtonian Fluid Mech., 154:153–169, 2008.

[52] S. Woltman, G. Jay, and G. Crawford. Liquid-crystal materials ﬁnd a

new order in biomedical applications. Nat. Mater., 6:929–938, 2007.

[53] T. White and D. Broer. Smectic phase in suspensions of gapped dna

duplexes. Nat. Mater., 14:1087–1098, 2015.

[54] M. Salamonczyk, J. Zhang, G. Portale, C. Zhu, E. Kentzinger, J. Gleeson,
A. Jakli, C. De Michele, J. Dhont, S. Sprunt, and E. Stiakakis. Smectic
phase in suspensions of gapped dna duplexes. Nat. Commun., 7:13358,
2016.

[55] G. Vroege and H. Lekkerkerker. Phase transitions in lyotropic colloidal

and polymer liquid crystals. Rep. Prog. Phys., 55:1241–1309, 1992.

[56] M. Fixman. Dynamics of semidilute polymer rods: An alternative to

cages. Phys. Rev. Lett., 55:2429–2432, 1985.

[57] A. Perera and G. N. Patey. The solution of the hypernetted-chain and
J.

Percus–Yevick approximations for ﬂuids of hard spherocylinders.
Chem. Phys., 89(9):5861–5868, 1988.

[58] Hartmut Graf and Hartmut Löwen. Density functional theory for hard
spherocylinders: Phase transitions in the bulk and in the presence of
external ﬁelds. J. Physics: Condens. Matter, 11(6):1435, 1999.

[59] H. Hervet, L. Léger, and F. Rondelez. Self-diffusion in polymer solutions:

A test for scaling and reptation. Phys. Rev. Lett., 42:1681–1684, 1979.

[60]

J. F. Maguire, J. P. McTague, and F. Rondelez. Rotational diffusion of
sterically interacting rodlike macromolecules. Phys. Rev. Lett., 45:1891–
1894, 1980.

[61] K. Zero and R. Pecora.

Rotational and translational diffusion in
semidilute solutions of rigid-rod macromolecules. Macromolecules,
15:87–93, 1982.

[62] P. Cobb and J. Butler. Simulations of concentrated suspensions of rigid
ﬁbers: Relationship between short-time diffusivities and the long-time
rotational diffusion. J. Chem. Phys., 123:054908, 2005.

96

[63] Y. Tao, W. den Otter, J. Padding, J. Dhont, and W. Briels. Brownian
dynamics simulations of the self- and collective rotational diffusion
coefﬁcients of rigid long thin rods. J. Chem. Phys., 122:244903, 2005.

[64] Y. Tao, W. den Otter, J. Dhont, and W. Briels.

Isotropic-nematic
spinodals of rigid long thin rodlike colloids by event-driven brownian
dynamics simulations. J. Chem. Phys., 124:134906, 2006.

[65] W. Yan, H. Zhang, and M. Shelley. Computing collision stress in
assemblies of active spherocylinders: Applications of a fast and generic
geometric method. J. Chem. Phys., 150:064109, 2019.

[66] P.G. Bolhuis, A. Stroobants, D. Frenkel, and H.N.W. Lekkerkerker.
Numerical study of the phase behavior of rodlike colloids with attractive
interactions. Journal of chemical physics, 107(5):1551–1564, 1997.

[67] Mihai Anitescu, James F Cremer, and Florian A Potra. Formulating
Journal of Structural

three-dimensional contact dynamics problems.
Mechanics, 24(4):405–437, 1996.

[68] A. Perera, S. Ravichandran, M. Moreau, and B. Bagchi. Single particle
and collective orientational relaxation in an anisotropic liquid near the
isotropic–nematic transition. J. Chem. Phys., 106:1280–1283, 1997.

[69] B. Berne and R. Pecora. Dynamic light scattering. Wiley, New York,

1976.

[70] H. Löwen. Brownian dynamics of hard spherocylinders. Phys. Rev. E,

14:1232–1242, 1994.

[71] D. Saintillan and M. Shelley. Orientational order and instabilities in

suspensions of self-locomoting rods. Phys. Rev. Lett., 99, 2007.

[72] Stavros D Peroukidis, Ken Lichtner, and Sabine HL Klapp. Tunable
structures of mixtures of magnetic particles in liquid-crystalline matrices.
Soft Matter, 11(30):5999–6008, 2015.

[73] D. Frenkel and J. Maguire. Molecular dynamics study of the dynamical
properties of an assembly of inﬁnitely thin hard rods. Mol. Phys., 49:503–
541, 1983.

[74]

J. Rallison.
interacting particles. J. Fluid Mech., 186:471–500, 1988.

Brownian diffusion in concentrated suspensions of

97

[75] Masao Doi. Molecular dynamics and rheological properties of
concentrated solutions of rodlike polymers in isotropic and liquid
crystalline phases. Journal of Polymer Science: Polymer Physics Edition,
19(2):229–243, 1981.

[76]

J Feng, CV Chaubal, and LG Leal. Closure approximations for the doi
theory: Which to use in simulating complex ﬂows of liquid-crystalline
polymers? Journal of Rheology, 42(5):1095–1119, 1998.

[77] E. M. Purcell. Life at low reynolds number. Am. J. Phys., 45:3–11, 1977.

[78] E. Lauga and T. Powers.

The hydrodynamics of swimming

microorganisms. Reports on Progress in Physics, 72(9):096601, 2009.

[79]

J. Teran, L. Fauci, and M. Shelley. Viscoelastic ﬂuid response can
increase the speed and efﬁciency of a free swimmer. Phys. Rev. Lett.,
104:038101, 2010.

[80] B. Thomases and R. Guy. Mechanisms of elastic enhancement and
hindrance for ﬁnite-length undulatory swimmers in viscoelastic ﬂuids.
Phys. Rev. Lett., 113:098102, 2014.

[81] D. Salazar, A. Roma, and H. Ceniceros. Numerical study of an
inextensible, ﬁnite swimmer in stokesian viscoelastic ﬂow. Phys. Fluids,
28:063101, 2016.

[82] S. Zhou, A. Sokolov, O. Lavrentovich, and I. Aranson. Living liquid

crystals. Proc. Nat. Acad. Sci., 111:1265–1270, 2014.

[83] O. Lavrentovich. Active colloids in liquid crystals. Current Opinion in

Colloid & Interface Science, 21:97 – 109, 2016.

[84] S. Zhou, O. Tovkach, D. Golovaty, A. Sokolov, I. Aranson, and
O. Lavrentovich. Dynamic states of swimming bacteria in a nematic
liquid crystal cell with homeotropic alignment. New Journal of Physics,
19:055006, 2017.

[85] M. Krieger, S. Spagnolie, and T. Powers. Microscale locomotion in a

nematic liquid crystal. Soft Matter, 11:9115–9125, 2015.

[86] M. Krieger, M. Dias, and T. Powers. Minimal model for transient

swimming in a liquid crystal. Eur. Phys. J. E, 38:94, 2015.

[87]

J. Shi and T. Powers. Swimming in an anisotropic ﬂuid: How speed
depends on alignment angle. Phys. Rev. Fluids, 2:123102, 2017.

98

[88] M. Krieger, S. Spagnolie, and T. Powers. Swimming with small and
large amplitude waves in a conﬁned liquid crystal. J. Non-Newtonian
Fluid Mech., 273:104185, 2019.
J. Feng and L. Leal. Simulating complex ﬂows of liquid-crystalline
polymers using the doi theory. J. Rheol., 41:1317–1335, 1997.

[89]

[90] L. Fauci and C. Peskin. A computational model of aquatic animal

locomotion. Journal of Computational Physics, 77:85–108, 1988.

[91] R. G. Larson. The Structure and Rheology of Complex Fluids. Oxford

University Press, New York, 1999.

[92] Pilhwa Lee and Charles W Wolgemuth. An immersed boundary method
for two-phase ﬂuids and gels and the swimming of caenorhabditis
elegans through viscoelastic ﬂuids. Physics of Fluids, 28(1):011901,
2016.

[93] Lisa J Fauci and Charles S Peskin. A computational model of aquatic
animal locomotion. Journal of Computational Physics, 77(1):85–108,
1988.

[94] Charles S Peskin. The immersed boundary method. Acta numerica,

11:479–517, 2002.

[95] T Vaithianathan and Lance R Collins.

Numerical approach to
simulating turbulent ﬂow of a viscoelastic polymer solution. Journal
of Computational Physics, 187(1):1–21, 2003.

[96] Daniel Salazar, Alexandre M Roma, and Hector D Ceniceros. Numerical
study of an inextensible, ﬁnite swimmer in stokesian viscoelastic ﬂow.
Physics of Fluids, 28(6):063101, 2016.

[97] Eric Lauga.

Propulsion in a viscoelastic ﬂuid. Physics of Fluids,

[98]

19(8):083104, 2007.
J Feng and LG Leal. Simulating complex ﬂows of liquid-crystalline
polymers using the doi theory. Journal of Rheology, 41(6):1317–1335,
1997.

[99] T. Gao and Z. Li. Self-driven droplet powered by active nematics. Phys.

Rev. Lett., 119:108002, 2017.

[100] T. Gao, M. Betterton, A. Jhang, and M. Shelley. Analytical structure,
dynamics, and coarse-graining of a kinetic model of an active ﬂuid. Phys.
Rev. Fluids, 2:093302, 2017.

99

[101] G.I. Taylor. Analysis of the swimming of microscopic organisms. Proc.

R. Soc. Lond. Ser. A, 209:447, 1951.

[102] Eric Lauga.

Propulsion in a viscoelastic ﬂuid.

Phys. Fluids,

19(8):083104, 2007.

[103] X. Shen and P. Arratia. Undulatory swimming in viscoelastic ﬂuids. Phys.

Rev. Lett., 106:208101, 2011.

100