PHYSICAL REVIEW E
108
, 024610 (2023)
Editors’ Suggestion
Wildebeest herds on rolling hills: Flocking on arbitrary curved surfaces
Christina L. Hueschen
,
1
,
*
Alexander R. Dunn
,
1
and Rob Phillips
2
,
3
,
†
1
Department of Chemical Engineering, Stanford University, Palo Alto, California 94305, USA
2
Department of Physics, California Institute of Technology, Pasadena, California 91125, USA
3
Division of Biology and Biological Engineering, California Institute of Technology, Pasadena, California 91125, USA
(Received 11 June 2022; accepted 10 July 2023; published 24 August 2023)
The collective behavior of active agents, whether herds of wildebeest or microscopic actin filaments propelled
by molecular motors, is an exciting frontier in biological and soft matter physics. Almost three decades ago,
Toner and Tu developed a continuum theory of the collective action of flocks, or herds, that helped launch the
modern field of active matter. One challenge faced when applying continuum active matter theories to living
phenomena is the complex geometric structure of biological environments. Both macroscopic and microscopic
herds move on asymmetric curved surfaces, like undulating grass plains or the surface layers of cells or embryos,
which can render problems analytically intractable. In this paper, we present a formulation of the Toner-Tu
flocking theory that uses the finite element method to solve the governing equations on arbitrary curved surfaces.
First, we test the developed formalism and its numerical implementation in channel flow with scattering obstacles
and on cylindrical and spherical surfaces, comparing our results to analytical solutions. We then progress to
surfaces with arbitrary curvature, moving beyond previously accessible problems to explore herding behavior on
a variety of landscapes. This approach allows the investigation of transients and dynamic solutions not revealed
by analytic methods. It also enables versatile incorporation of new geometries and boundary conditions and
efficient sweeps of parameter space. Looking forward, the paper presented here lays the groundwork for a
dialogue between Toner-Tu theory and data on collective motion in biologically relevant geometries, from drone
footage of migrating animal herds to movies of microscopic cytoskeletal flows within cells.
DOI:
10.1103/PhysRevE.108.024610
I. INTRODUCTION
The beautiful collective motions of flocking or herding
animals have mesmerized the human mind for millennia [
1
]
and inspired the modern field of study of active matter [
2
–
4
].
In recent decades, active matter theories have been used to
describe collective motion in living systems across nearly a
billionfold difference in scales [
5
], from starling flocks wheel-
ingoverRomeat10m
/
s[
6
], to 50
μ
m
/
s flows of cytoplasm
in cm-size internodal cells of the algae
Chara
[
7
], to 0
.
1
μ
m
/
s
flows of actin filaments and myosin molecules in 50
μ
mdevel-
oping
C. elegans
worm embryos [
8
]. In these macroscopic and
microscopic contexts, beautiful self-organization phenomena
often take place on surfaces. The collective motion of flock-
ing sheep or kilometer-wide migrating wildebeest herds are
influenced by the hills, valleys, or canyons of the landscape
on which they travel. Similarly, the flows of actin and myosin
molecules mentioned above occur in thin layers at the surface
of the embryo or cell, and are thus constrained by its surface
*
Corresponding author: chueschen@gmail.com
†
Corresponding author: phillips@pboc.caltech.edu
Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International
license. Further
distribution of this work must maintain attribution to the author(s)
and the published article’s title, journal citation, and DOI.
topology and shape. While modern techniques such as dy-
namical drone imaging [
9
] and high-resolution fluorescence
microscopy [
10
,
11
] enable experimental measurements of
these animal, cell, or molecule kinematics, a dialogue between
measurement and theory requires predictions that incorporate
the complex shapes of the real-world surfaces on which they
move.
In this paper, we present a general curved-surface for-
mulation and numerical implementation of a minimalistic
continuum theory of flocking or herding active matter with
the hope that it will prove useful to others interested in ex-
ploring continuum theory predictions on complex geometries.
For active matter systems such as bacterial swarms [
12
],
active colloidal fluids [
13
], self-propelled rods [
14
], or pu-
rified cytoskeletal networks [
15
], studying the contribution
of engineered confinement geometries to emergent patterns
has proven a fruitful path [
16
–
20
]. Understanding the con-
tribution of
biological
geometries to pattern formation is an
exciting direction of growth [
21
–
26
]. Indeed, our own interest
in solving active matter theories on arbitrary curved surfaces
was initially inspired by a need to predict emergent actin
polarity patterns at the surface of single-celled
Toxoplasma
gondii
parasites, whose gliding motility is driven and directed
by this surface actin layer [
27
]. These cells have a beautiful
but complex shape which lacks the symmetries that license
traditional analytic approaches.
We focus on a classic continuum active matter model orig-
inally developed by Toner and Tu [
3
,
4
,
28
], inspired by the
2470-0045/2023/108(2)/024610(28)
024610-1
Published by the American Physical Society
HUESCHEN, DUNN AND PHILLIPS
PHYSICAL REVIEW E
108
, 024610 (2023)
work of Vicsek
et al.
[
2
], to describe the collective behavior
of flocking or herding animals. The Toner-Tu theory helped
launch the modern field of active matter [
29
] and can describe
collections of dry, polar, and self-propelled agents at any
length scale. We develop a general curved surface framework
for the Toner-Tu theory and implement it in the finite ele-
ment setting, enabling its convenient use on complex surfaces.
While flocking theory—and, in particular, the herding of
wildebeests—serves as our example in this paper, the general
curved surface formulation and finite element method (FEM)
approach presented here may prove useful for any continuum
active matter theory. There is a rich and beautiful literature as-
sociated with formulating field equations on arbitrary curved
surfaces [
30
–
36
]. It has been said of these formulations, “The
complexity of the equations may explain why they are so
often written but never solved for arbitrary surfaces.” [
34
].
Because of our interest in collective motions on biological sur-
faces, we could not afford to avoid solving these equations on
such surfaces. We used a finite element approach, which
can be flexibly adjusted to an arbitrary choice of geometry
and permits an exploration of active, self-organized solutions
predicted by the Toner-Tu theory. With our general surface
formulation and finite element implementation in hand, we
test our framework and approach on cylindrical and spherical
geometries for which analytic solutions exist as well. Satisfied
by the agreement between our numerical and analytic results,
we explore flocking phenomena on a broad collection of
surfaces.
The remainder of the paper is organized as follows. In
Sec.
II
, we show how the flocking theory of Toner and Tu
can be recast in a general surface form to describe flocking
motions on arbitrary curved surfaces. In Sec.
III
, we con-
struct a finite element surface formulation that respects the
low symmetry of realistic curved geometries while permitting
numerical analysis of the dynamics. In Sec.
IV
,weturnto
the great wildebeest herds for inspiration and perform simple
parameter estimates for Toner-Tu wildebeest herding. We then
explore the dimensionless ratios that appear when recasting
the Toner-Tu equations in dimensionless form. In Sec.
V
, with
the full surface formulation and its numerical implementation
in hand, we explore scaling relationships and changes of dy-
namical state that arise for herds in channels with scattering
obstacles, and we use cylindrical and spherical surfaces to
compare our FEM results to corresponding analytic solutions.
Our approach also allows the exploration of transients and
dynamic solutions not revealed by analytic methods. Finally,
in Sec.
VI
, we playfully use the curved-space formalism and
its finite element implementation for case studies of wilde-
beest herding on landscapes of rolling hills, and we note that
this paper serves as a theoretical foundation to solve flocking
theories on real-world curved surfaces of biological interest.
II. FLOCKING THEORY FOR ARBITRARY
CURVED SURFACES
A. A minimal Toner-Tu theory in the plane
The first step in the development of continuum theories of
active matter involves selecting the relevant field variables.
Following the classic work of Toner and Tu [
3
,
4
] and inspired
by observations of herds of land animals like wildebeests,
we consider a two-dimensional density field
ρ
(
r
,
t
) and a
corresponding two-dimensional velocity field
v
(
r
,
t
), which
captures both the orientation and speed of a polar agent.
Throughout this paper, wildebeest herds will serve as inspi-
ration and example, although in a sense the term wildebeest
is our shorthand for a self-propelled polar agent at any length
scale. We also note that the Toner-Tu model used here de-
scribes
dry
systems, in which momentum is not conserved
and hydrodynamic coupling between the agents is negligible
relative to frictional drag. With variables
ρ
and
v
in hand,
our next step is to write the partial differential equations that
describe the spatiotemporal evolution of those field variables.
The first governing equation is the continuity equation given
by
∂ρ
∂
t
+
∂
(
ρ
v
i
)
∂
x
i
=
0
,
(1)
which allows modulations in the density field but enforces
mass conservation, forbidding wildebeest birth or death in
the midst of herding phenomena. Note that we are using the
Einstein summation convention, which tells us to sum over all
repeated indices. For example,
a
·
b
=
3
i
=
1
a
i
b
i
=
a
i
b
i
.The
second governing equation is a minimal representation of the
dynamics of the
v
field offered by Toner and Tu [
3
] and is
given by
∂
v
i
∂
t
=
[
α
(
ρ
−
ρ
c
)
−
β
v
j
v
j
]
v
i
−
σ
∂ρ
∂
x
i
+
D
∇
2
v
i
−
λ
v
j
∂
v
i
∂
x
j
,
(2)
where
ρ
c
is the critical density above which the herd moves
coherently, the ratio of
α
(
ρ
−
ρ
c
) and
β
sets the wildebeest
mean-field speed
v
pref
=
√
α
(
ρ
−
ρ
c
)
/β
,theterm
σ∂ρ/∂
x
i
provides an effective pressure,
D
tunes wildebeest alignment
and speed matching with neighbors, and
λ
tunes velocity self-
advection. We can conceptualize each term as an update rule
that computes incremental changes in the velocity vector at
each point in space and at every step in time. In Fig.
1
and in
the remainder of this section, we seek to provide an intuitive
interpretation of each term and its contribution to updating the
velocity field.
The first term on the right-hand side of Eq. (
2
), the pre-
ferred speed term, pushes the velocity magnitude toward the
characteristic speed of a wildebeest,
v
pref
[
37
]. The second
is based on an equation of state originally used by Toner
and Tu to relate density and pressure [
4
]. This pressure term
punishes gradients in density, adjusting velocity to flatten the
density field. The third, the neighbor coupling term, provides
a smoothing or diffusion of velocity (orientation and speed)
that reflects coordination between nearby wildebeests. Inter-
estingly, if wildebeests adopt the rule that a given wildebeest
averages the difference between its own velocity and that of
its neighbors, mathematically, the result is a Laplacian term
like that seen in the minimal model [
37
]. The final term
has analogy to the gradient component of the material time
derivative in the Navier-Stokes equations. In essence, the ve-
locity field advects itself; wildebeests move along a direction
dictated by their orientation and velocity, and they bring that
orientation and velocity with them. In the case of pure velocity
self-advection,
λ
=
1, but this is not necessarily the case for
024610-2
WILDEBEEST HERDS ON ROLLING HILLS: FLOCKING ...
PHYSICAL REVIEW E
108
, 024610 (2023)
(a)(b)
(c)(d)
FIG. 1. Building intuition for the physics of the velocity field in the Toner-Tu theory. The different terms in the Toner-Tu analysis are
represented graphically for a one-dimensional herd. We can think of each term as providing an update to the current velocity. In each case, an
example of a simple velocity profile (blue) or density profile [graphic in (b)] is shown alongside the corresponding velocity update (green).
(b) The preferred speed term increases the velocity of wildebeests that are moving too slow and decreases the velocity of wildebeests that are
moving too fast. (b) The pressure term punishes gradients in density, adjusting velocity to flatten these gradients. Given mass conservation in
a finite system, the pressure term alone would lead to a steady state of uniform density
ρ
=
ρ
0
. (c) The neighbor coupling term captures the
velocity adjustment made by a wildebeest to better match its neighbors, smoothing out differences in velocity. A given wildebeest (middle
black dot) adjusts its velocity by an amount represented by the green arrow: an average of the difference between its velocity and its two
neighbors’ velocities. This averaging of two differences is mathematically analogous to taking a local second derivative or local Laplacian,
∇
2
v
. (d) Finally, the advection term ensures that the filament velocity field is swept along according to its own velocities, just as fluid velocity
is self-advected in the Navier-Stokes equations. Mathematically, the velocity adjustment needed to ensure velocity advection is a function of
the spatial gradient of velocity (
−
∇
v
, how mismatched in velocity nearby wildebeests are) and the velocity itself (
v
, how fast that mismatch is
carried along by the herd).
active flocks [
37
,
38
]. In the case of wildebeest herds, we can
conceptualize this effect by considering that
λ
=
1
+
ξ
, where
ξ
reflects a behavioral response to gradients in velocity; for
example, wildebeests may resist running quickly into a steep
gradient of decreasing velocity, and may slow down. For a full
pedagogical derivation of the Toner-Tu model, we recommend
a series of lectures by Toner [
37
] in which he uses symmetry
arguments to infer what terms should be kept in a complete
continuum description of herding and also provides intuitive
arguments about what these terms mean.
B. Formulating the theory for arbitrary curved surfaces
Inspired by a desire to make contact with phenomena of
the natural world, like wildebeests navigating an undulat-
ing landscape, sheep flocks crossing hilly pastures, and the
surfaces flows of flocking actin we study in Ref. [
27
], we
sought to solve the minimal Toner-Tu theory presented above
on complex and asymmetric surface geometries. We consider
here flocking on nondeformable surfaces, but we refer the
interested reader to earlier work on surface hydrodynamics
024610-3
HUESCHEN, DUNN AND PHILLIPS
PHYSICAL REVIEW E
108
, 024610 (2023)
FIG. 2. Illustration of the projection operator essential to the general surface formulation and finite element treatment of the Toner-Tu
equations on curved surfaces. The surface of interest is represented by a collection of nodes (vertices of triangles) and a field of surface normal
vectors,
n
(black arrows). An arbitrary vector
v
is projected onto the tangent plane (shown in gold), permitting a decomposition of the form
v
=
v
+
v
⊥
.
in the completely general case in which the surface itself can
evolve over time [
32
–
34
,
36
]. Predicting flocking or herding
behavior on arbitrary surfaces requires us to reformulate the
theory in a more general way that accounts for curvature.
Instead of basing our formulation on a parameterized surface
and the intrinsic differential geometry tools that this approach
licenses, as done in Ref. [
39
], from the outset we have in mind
arbitrary surfaces that can be represented using finite element
meshes and described by a field of local normal vectors,
n
.
Our choice and use of this extrinsic differential geometry
and FEM approach was aided by the work of many, includ-
ing Refs. [
33
–
35
,
40
–
45
] and Chap. 3 of the Supplemental
Material for Ref. [
46
]. The finite element setting permits us
the convenience of carrying out the mathematics in the full
three-dimensional setting of
R
3
, while using our knowledge
of the normal vectors everywhere on the surface of interest
to project our governing equations onto the surface. While
velocity is described by the three-dimensional vector
v
, both
v
and the scalar
ρ
are defined only on the surface. At every point
on the surface, derivatives evaluated in the usual
R
3
way are
projected onto the tangent plane using the local normal. For
an insightful description of this extrinsic differential geometry
approach to handling curved surfaces and its mathematical
equivalence to the intrinsic differential geometry strategy, we
recommend Chaps. 22 and 23 of Needham [
47
]. Central to
the extrinsic geometry approach is the projection operator,
defined as
projection onto tangent plane
=
P
=
I
−
n
⊗
n
,
(3)
where
I
is the identity matrix and
n
⊗
n
is the outer product of
the surface normal vector as shown in Fig.
2
. To make sense
of this expression mathematically, we recall that the outer
product
n
⊗
n
is defined through its action on a vector
v
as
(
n
⊗
n
)
v
=
n
(
n
·
v
)
.
(4)
We can write
v
=
v
+
v
⊥
, where
v
is the component of
v
in the tangent plane of the surface and
v
⊥
is normal to the
surface. The action of
I
−
n
⊗
n
on a vector
v
is given by
(
I
−
n
⊗
n
)
v
=
v
−
n
(
n
·
v
)
=
v
−
v
⊥
=
v
,
(5)
as illustrated in Fig.
2
. We can write the projection operator
in component form as
P
ij
=
δ
ij
−
n
i
n
j
,
(6)
recalling that the normal vector to the surface is given by
n
=
(
n
1
,
n
2
,
n
3
), or in full component form as
P
=
⎡
⎢
⎣
1
−
n
2
1
−
n
1
n
2
−
n
1
n
3
−
n
1
n
2
1
−
n
2
2
−
n
2
n
3
−
n
1
n
3
−
n
2
n
3
1
−
n
2
3
⎤
⎥
⎦
.
(7)
Using the projection operator allows us to perform calcula-
tions in the ordinary three-dimensional space within which
the surface of interest is embedded, but then to pick off only
the pieces of the resulting vectors that live within the surface.
In the remainder of this section, we examine each term
in the minimal Toner-Tu equations and translate it into its
projected form. In Appendix Sec. 1, we define the projected
surface form of calculus operators as an additional reference
for the reader. First, to modify the preferred speed term to
its curved-surface implementation, we note that it is now
v
,
the in-plane velocity, that has a privileged magnitude. This
magnitude is imposed through the condition
preferred speed term
=
[
α
(
ρ
−
ρ
c
)
−
β
v
j
v
j
]
v
i
.
(8)
If
|
v
|
is either larger or smaller than the privileged value
v
pref
=
√
α
(
ρ
−
ρ
c
)
/β
, this term will adjust the velocity to-
wards that steady-state magnitude.
Next, we consider the pressure term, whose physical origin
comes from a model of pressure in powers of density of the
form introduced by Toner and Tu [
3
],
P
(
ρ
)
=
∞
n
=
1
σ
n
(
ρ
−
ρ
0
)
n
,
(9)
where
ρ
0
is the mean density. In the minimal Toner-Tu
theory we adopt here, only the first-order term in that expan-
sion is kept, with the notational simplification that
σ
1
=
σ
,
024610-4
WILDEBEEST HERDS ON ROLLING HILLS: FLOCKING ...
PHYSICAL REVIEW E
108
, 024610 (2023)
resulting in
pressure term
=−
∂
P
∂
x
i
=−
σ
∂ρ
∂
x
i
.
(10)
We note that while
ρ
0
is not explicitly present in the gradient
of the first-order pressure term, for the case considered here
of a finite surface on which total density is conserved, this
pressure term effectively maintains a density range centered
around the mean density
ρ
0
established by our choice of
initial condition. The curved-space version of the pressure
term requires projecting the full 3D gradient of the density
onto the tangent plane, using the projection operator defined
in Eq. (
3
). We follow Jankuhn
et al.
[
34
] in introducing the
notation
∇
for the projected surface gradient operator, where
the subscript
indicates that the gradient is evaluated on the
surface of interest. Mathematically, this amounts to computing
∇
ρ
=
(
I
−
n
⊗
n
)
∇
ρ,
(11)
where
∇
is the ordinary, three-dimensional gradient in Carte-
sian coordinates. This can be rewritten in component form
as
(
∇
ρ
)
i
=
[
(
I
−
n
⊗
n
)
∇
ρ
]
i
=
(
δ
ij
−
n
i
n
j
)
∂ρ
∂
x
j
=
∂ρ
∂
x
i
−
n
i
n
j
∂ρ
∂
x
j
.
(12)
The next term in the Toner-Tu equations that we consider
in its curved-space format is the advection term in Eq. (
2
),
namely,
advection term
=−
λ
v
j
∂
v
i
∂
x
j
.
(13)
In direct notation, the curved-space version of this term has
the form
advection term
=−
λ
(
v
·
∇
)
v
,
(14)
which involves the curved-space version of the velocity gra-
dient tensor. In Appendix Sec. 1, we describe how to compute
this tensor, which stated simply is
∇
v
=
P
(
∇
v
)
P
,
(15)
where
P
is the projection operator defined in Eq. (
3
). In indi-
cial notation, this leads to the result
(
∇
v
)
ij
=
[
P
(
∇
v
)
P
]
ij
=
P
ik
∂
v
k
∂
x
l
P
lj
.
(16)
Invoking the definition of the projection operator from Eq. (
6
),
this expression simplifies to the result
(
∇
v
)
ij
=
[
P
(
∇
v
)
P
]
ij
=
∂
v
i
∂
x
j
−
n
l
n
j
∂
v
i
∂
x
l
−
n
i
n
k
∂
v
k
∂
x
j
−
n
l
n
j
∂
v
k
∂
x
l
.
(17)
Thus, the curved-space advection term can be written as
−
λ
[
P
(
∇
v
)
P
]
ij
v
j
=−
λ
∂
v
i
∂
x
j
−
n
l
n
j
∂
v
i
∂
x
l
−
n
i
n
k
∂
v
k
∂
x
j
−
n
l
n
j
∂
v
k
∂
x
l
v
j
.
(18)
With these examples of tangent-plane calculus established,
we now turn to the most tricky of the terms in the Toner-
Tu equations, the surface-projected version of the neighbor
coupling term,
D
∇
2
v
. First, we recall that the Laplacian of
a vector field in normal 3D Cartesian space is defined as
∇
2
v
=
∇
·
(
∇
v
). Further, we note that
∇
v
is itself a tensor. As
shown by Jankuhn
et al.
[
34
], the surface-projected version of
the Laplacian term (see their Eq. (3.16) for the surface Navier-
Stokes equations for the tangential velocity on a stationary
surface), is therefore given by
∇
2
v
flat Toner-Tu
=
div(
∇
v
)
in the plane
−→
P
div
(
G
(
v
))
curved surface
.
(19)
To begin to unpack this expression, we note that
P
is the
projection operator defined in Eq. (
3
) and that the tensor
G
is the surface velocity gradient,
G
(
v
)
=
∇
v
,
(20)
already presented in Eq. (
17
) and repeated here in component
form as
G
ij
=
∂
v
i
∂
x
j
−
n
l
n
j
∂
v
i
∂
x
l
−
n
i
n
k
∂
v
k
∂
x
j
−
n
l
n
j
∂
v
k
∂
x
l
.
(21)
We note that Jankuhn
et al.
[
34
] use a symmetrized version
of the velocity gradient since they are deriving the surface
versions of the Navier-Stokes equations and are thus taking
the divergence of a stress. For the Toner-Tu case of interest
here, the psychological viscosity that comes from neighboring
wildebeests in the herd comparing their velocities is equiva-
lent to the divergence of the velocity gradient itself. We then
invoke the definition of the surface-projected divergence of a
tensor
G
presented in Appendix
A1
, giving rise to a vector of
the form
div
(
G
)
l
=
∂
G
lj
∂
x
j
−
n
j
n
k
∂
G
lj
∂
x
k
.
(22)
We also introduce the shorthand notation for this divergence
as
∂
G
lj
∂
x
j
=
∂
G
lj
∂
x
j
−
n
j
n
k
∂
G
lj
∂
x
k
(23)
to simplify some of the complex expressions to follow. Fi-
nally, we can write the
i
th component of the surface version
of the Toner-Tu term
D
∇
2
v
in indicial notation as
D
[
P
div
(
G
(
v
))]
i
=
DP
il
∂
G
lj
∂
x
j
−
n
j
n
k
∂
G
lj
∂
x
k
.
(24)
We now assemble the results of Eqs. (
8
), (
12
), (
18
), and
(
24
) to construct a complete curved-surface formulation of the
024610-5
HUESCHEN, DUNN AND PHILLIPS
PHYSICAL REVIEW E
108
, 024610 (2023)
minimal Toner-Tu equations. Specifically, we have
∂
v
i
∂
t
=
[
α
(
ρ
−
ρ
c
)
−
β
v
j
v
j
]
v
i
−
σ
∂ρ
∂
x
i
−
n
i
n
j
∂ρ
∂
x
j
+
DP
il
∂
G
lj
∂
x
j
−
n
j
n
k
∂
G
lj
∂
x
k
−
λ
[
P
(
∇
v
)
P
]
ij
v
j
,
(25)
which can be streamlined to the alternative form
∂
v
i
∂
t
=
[
α
(
ρ
−
ρ
c
)
−
β
v
j
v
j
]
v
i
−
σ
∂ρ
∂
x
i
+
DP
il
∂
G
lj
∂
x
j
−
λ
v
j
∂
v
i
∂
x
j
−
n
i
n
k
∂
v
k
∂
x
j
.
(26)
Similarly, the curved-space formulation of the governing
equation for density can be written as
∂ρ
∂
t
=−
∂
(
ρ
v
i
)
∂
x
i
+
n
i
n
j
∂
(
ρ
v
i
)
∂
x
j
=−
∂ρ
v
i
∂
x
i
.
(27)
We now have a complete formulation of our minimal Toner-Tu
equations for the general surface context, requiring only a
description of the surface in the language of normal vectors.
We next turn to the implementation of this general surface for-
mulation in a fashion consistent with finite element treatments
on arbitrary surfaces.
III. FORMULATING THE SURFACE TONER-TU
EQUATIONS FOR NUMERICAL IMPLEMENTATION
Numerically solving active matter equations on complex
surfaces relevant to the living world presents a practical
challenge. It precludes the differential geometric formalism
used to describe parameterized surfaces, which would lead
to equations featuring covariant derivatives [e.g., Eq. (
81
)in
Sec.
VB
]. In the finite element setting, the surface of interest
is represented by a collection of nodes and corresponding
surface normals, as shown in Fig.
2
. Using these surface nor-
mals, we perform surface projections of the full three-space
derivatives following Jankuhn
et al.
[
34
]. In the previous sec-
tion, we presented a formal statement in Eq. (
25
) for handling
the minimal Toner-Tu model on an arbitrary surface, using
these surface normals. We now recast those curved-surface
equations once more, in a fashion consonant with a FEM
solver. To formulate the equations in an expression convenient
for the FEM, we aim to rewrite the Toner-Tu equations in the
form
∂
v
∂
t
+
∇
·
J
(
v
)
=
f
(
v
)
.
(28)
Here, the whole formulation comes down to the definitions
of
J
(
v
)
and
f
(
v
)
. The fluxlike quantity
J
(
v
)
is a 3
×
3matrix
defined such that
∇
·
J
(
v
)
=
∂
∂
x
1
,
∂
∂
x
2
,
∂
∂
x
3
⎡
⎢
⎣
J
11
J
12
J
13
J
21
J
22
J
23
J
31
J
32
J
33
⎤
⎥
⎦
,
(29)
recalling that the divergence of a second-rank tensor is a
vector. In indicial notation, this can be written as
(
∇
·
J
(
v
)
)
i
=
∂
J
(
v
)
ji
∂
x
j
.
(30)
We need to define the components of the tensor
J
(
v
)
such that
they yield the correct Toner-Tu terms. As we will see below,
J
(
v
)
is not symmetric despite its superficial resemblance to a
stress tensor. To set notation and make sure that the strategy is
clear, we begin by demonstrating how to implement the flat-
space version of the Toner-Tu equations in a form consistent
with Eq. (
2
). If we define the force term
f
(
v
)
as
f
(
v
)
i
=
[
α
(
ρ
−
ρ
c
)
−
β
v
j
v
j
]
v
i
−
σ
∂ρ
∂
x
i
−
λ
v
j
∂
v
i
∂
x
j
,
(31)
then the remaining terms are captured if we define the tensor
J
(
v
)
as
J
(
v
)
=
⎡
⎢
⎢
⎣
−
D
∂
v
1
∂
x
1
−
D
∂
v
2
∂
x
1
−
D
∂
v
3
∂
x
1
−
D
∂
v
1
∂
x
2
−
D
∂
v
2
∂
x
2
−
D
∂
v
3
∂
x
2
−
D
∂
v
1
∂
x
3
−
D
∂
v
2
∂
x
3
−
D
∂
v
3
∂
x
3
⎤
⎥
⎥
⎦
.
(32)
Considering the one-component of velocity, we see that
(
∇
·
J
(
v
)
)
1
=
∂
J
j
1
∂
x
j
=
∂
J
11
∂
x
1
+
∂
J
21
∂
x
2
+
∂
J
31
∂
x
3
.
(33)
Plugging in the components of
J
(
v
)
gives the result
(
∇
·
J
(
v
)
)
1
=
∂
∂
x
1
−
D
∂
v
1
∂
x
1
+
∂
∂
x
2
−
D
∂
v
1
∂
x
2
+
∂
∂
x
3
−
D
∂
v
1
∂
x
3
.
(34)
Evaluating these derivatives leads to the result
(
∇
·
J
(
v
)
)
1
=−
D
∂
2
v
1
∂
x
2
1
−
D
∂
2
v
1
∂
x
2
2
−
D
∂
2
v
1
∂
x
2
3
,
(35)
as expected from the original Toner-Tu equations. Combin-
ing the contributions from Eqs. (
31
) and (
32
)intheform
∂
v
/∂
t
+
∇
·
J
(
v
)
=
f
(
v
)
, we recover Eq. (
2
) precisely, as we
set out to do. By defining the quantities
J
(
v
)
and
f
(
v
)
,wehave
successfully reframed the flat-space Toner-Tu equations in a
format that will be conveniently implemented in the finite
element setting. We now need to tackle the more demanding
formulation for an arbitrary curved surface.
Abstractly, our finite element version of the curved-space
Toner-Tu equations is written as
∂
v
∂
t
=−
div
J
(
v
)
+
f
(
v
)
,
(36)
where div
is the surface projected version of the divergence
introduced in Eq. (
22
). We now introduce the definition
J
(
v
)
ji
=−
DP
il
G
lj
,
(37)
where we recall
P
from Eq. (
3
) and
G
from Eq. (
20
). Given
this definition, we can attempt to write our Toner-Tu equations
as
∂
v
i
∂
t
=−
D
∂
(
−
P
il
G
lj
)
∂
x
j
+
f
(
v
)
i
,
(38)
024610-6
WILDEBEEST HERDS ON ROLLING HILLS: FLOCKING ...
PHYSICAL REVIEW E
108
, 024610 (2023)
where we use the shorthand notation introduced in Eq. (
23
).
This can be rewritten as
∂
v
i
∂
t
=
DP
il
∂
G
lj
∂
x
j
P
div
(
∇
v
)
+
D
∂
P
il
∂
x
j
G
lj
unwanted term
+
f
(
v
)
i
.
(39)
Unfortunately, as written, these equations contain an extra
term that is not present in the Toner-Tu formulation. To re-
move this unwanted extra term, we must introduce a new term
in force
f
(
v
)
that subtracts the unwanted term:
f
fict
i
=−
D
∂
P
il
∂
x
j
G
lj
.
(40)
We view this as a mathematical trick that allows us to use a
formalism convenient for FEM analysis and have not sought
to find a physical interpretation of this force in the way that
ideas such as the Coriolis force arises in mechanics. In this
case, our
real
Toner-Tu equations can be written as
∂
v
i
∂
t
=
DP
il
∂
G
lj
∂
x
j
P
div
(
∇
v
)
+
D
∂
P
il
∂
x
j
G
lj
unwanted term
+
f
(
v
)
i
+
f
fict
i
.
(41)
We now have precisely the equations we want, cast in the form
we will use in the finite element setting.
As we saw above, the force term
f
(
v
)
+
f
fict
needs to account for those terms that are not present in
∇
·
J
(
v
)
and to subtract
terms that are present in
∇
·
J
(
v
)
but unwanted. To that end, the one-component of
f
(
v
)
+
f
fict
takes the form
f
(
v
)
1
+
f
fict
1
=
(
α
(
ρ
−
ρ
c
)
−
β
((
v
1
)
2
+
(
v
2
)
2
+
(
v
3
)
2
))
v
1
−
σ
∂ρ
∂
x
1
−
λ
v
1
∂
v
1
∂
x
1
+
v
2
∂
v
1
∂
x
2
+
v
3
∂
v
1
∂
x
3
−
v
1
n
1
n
1
∂
v
1
∂
x
1
−
v
1
n
1
n
2
∂
v
2
∂
x
1
−
v
1
n
1
n
3
∂
v
3
∂
x
1
−
v
2
n
1
n
1
∂
v
1
∂
x
2
−
v
2
n
1
n
2
∂
v
2
∂
x
2
−
v
2
n
1
n
3
∂
v
3
∂
x
2
−
v
3
n
1
n
1
∂
v
1
∂
x
3
−
v
3
n
1
n
2
∂
v
2
∂
x
3
−
v
3
n
1
n
3
∂
v
3
∂
x
3
−
D
∂
P
11
∂
x
1
G
11
+
∂
P
11
∂
x
2
G
12
+
∂
P
11
∂
x
3
G
13
−
D
∂
P
12
∂
x
1
G
21
+
∂
P
12
∂
x
2
G
22
+
∂
P
12
∂
x
3
G
23
−
D
∂
P
13
∂
x
1
G
31
+
∂
P
13
∂
x
2
G
32
+
∂
P
13
∂
x
3
G
33
.
(42)
Note that the
α
and
β
terms capture the preferred speed contribution, the
σ
term captures the pressure, the terms multiplied by
λ
capture the advection contribution to the Toner-Tu equations, and the final terms involving
P
and
G
subtract off the fictitious
force. We can repeat a similar analysis for the two- and three-components of the force as
f
(
v
)
2
+
f
fict
2
=
(
α
(
ρ
−
ρ
c
)
−
β
((
v
1
)
2
+
(
v
2
)
2
+
(
v
3
)
2
))
v
2
−
σ
∂ρ
∂
x
2
−
λ
v
1
∂
v
2
∂
x
1
+
v
2
∂
v
2
∂
x
2
+
v
3
∂
v
2
∂
x
3
−
v
1
n
2
n
1
∂
v
1
∂
x
1
−
v
1
n
2
n
2
∂
v
2
∂
x
1
−
v
1
n
2
n
3
∂
v
3
∂
x
1
−
v
2
n
2
n
1
∂
v
1
∂
x
2
−
v
2
n
2
n
2
∂
v
2
∂
x
2
−
v
2
n
2
n
3
∂
v
3
∂
x
2
−
v
3
n
2
n
1
∂
v
1
∂
x
3
−
v
3
n
2
n
2
∂
v
2
∂
x
3
−
v
3
n
2
n
3
∂
v
3
∂
x
3
−
D
∂
P
21
∂
x
1
G
11
+
∂
P
21
∂
x
2
G
12
+
∂
P
21
∂
x
3
G
13
−
D
∂
P
22
∂
x
1
G
21
+
∂
P
22
∂
x
2
G
22
+
∂
P
22
∂
x
3
G
23
−
D
∂
P
23
∂
x
1
G
31
+
∂
P
23
∂
x
2
G
32
+
∂
P
23
∂
x
3
G
33
(43)
and
f
(
v
)
3
+
f
fict
3
=
(
α
(
ρ
−
ρ
c
)
−
β
((
v
1
)
2
+
(
v
2
)
2
+
(
v
3
)
2
))
v
3
−
σ
∂ρ
∂
x
3
−
λ
v
1
∂
v
3
∂
x
1
+
v
2
∂
v
3
∂
x
2
+
v
3
∂
v
3
∂
x
3
−
v
1
n
3
n
1
∂
v
1
∂
x
1
−
v
1
n
3
n
2
∂
v
2
∂
x
1
−
v
1
n
3
n
3
∂
v
3
∂
x
1
−
v
2
n
3
n
1
∂
v
1
∂
x
2
−
v
2
n
3
n
2
∂
v
2
∂
x
2
−
v
2
n
3
n
3
∂
v
3
∂
x
2
−
v
3
n
3
n
1
∂
v
1
∂
x
3
−
v
3
n
3
n
2
∂
v
2
∂
x
3
−
v
3
n
3
n
3
∂
v
3
∂
x
3
−
D
∂
P
31
∂
x
1
G
11
+
∂
P
31
∂
x
2
G
12
+
∂
P
31
∂
x
3
G
13
−
D
∂
P
32
∂
x
1
G
21
+
∂
P
32
∂
x
2
G
22
+
∂
P
32
∂
x
3
G
23
−
D
∂
P
33
∂
x
1
G
31
+
∂
P
33
∂
x
2
G
32
+
∂
P
33
∂
x
3
G
33
.
(44)
By setting
∂
v
/∂
t
equal to the sum of
−
div
J
(
v
)
and the
f
(
v
)
+
f
fict
terms, we have fully reproduced the curved-space version of
the Toner-Tu equations.
024610-7