Active nematic fluids on Riemannian 2-manifolds
Cuncheng Zhu***Department of Mechanical and Aerospace Engineering, University
of California San Diego, USA<EMAIL_ADDRESS>David Saintillan†††Department
of Mechanical and Aerospace Engineering, University of California San Diego,
USA<EMAIL_ADDRESS>Albert Chern‡‡‡Department of Computer Science and
Engineering, University of California San Diego, USA<EMAIL_ADDRESS>

###### Abstract

Recent advances in cell biology and experimental techniques using
reconstituted cell extracts have generated significant interest in
understanding how geometry and topology influence active fluid dynamics. In
this work, we present a comprehensive continuous theory and computational
method to explore the dynamics of active nematic fluids on arbitrary surfaces
without topological constraints. The fluid velocity and nematic order
parameter are represented as the sections of the complex line bundle of a
2-manifold. We introduce the Levi-Civita connection and surface curvature form
within the framework of complex line bundles. By adopting this geometric
approach, we introduce a gauge-invariant discretization method that preserves
the continuous local-to-global theorems in differential geometry. We establish
a nematic Laplacian on complex functions that can accommodate fractional
topological charges through the covariant derivative on the complex nematic
representation. We formulate advection of the nematic field based on a
unifying definition of the Lie derivative, resulting in a stable geometric
semi-Lagrangian discretization scheme for transport by the flow. In general,
the proposed surface-based method offers an efficient and stable means to
investigate the influence of local curvature and global topology on the 2D
hydrodynamics of active nematic systems. Moreover, the complex line
representation of the nematic field and the unifying Lie advection present a
systematic approach for generalizing our method to active $k$-atic systems.

##  1 Introduction

Active nematics, a unique class of systems defined by the interaction between
intrinsic activity and orientational order, have emerged as a captivating
paradigm to explore complex dynamical phenomena across scales. In biophysical
contexts, deciphering the dynamics of these active matter systems offers
profound insights into a myriad of processes, encompassing collective dynamics
in motile bacterial suspensions, structural organization in cytoskeletal
assemblies, the dynamics of self-propelled colloidal particles, and the
processes governing tissue morphogenesis [1, 2]. Recent experimental
discoveries and the potential for significant advances in fields like material
science and biophysics have sparked a growing interest in understanding the
intricate connections between the dynamic behavior of active nematic systems
and the geometric and topological constraints imposed by spatial confinement
[3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. However, the theoretical understanding
and technological advancements in studying these systems are currently
limited. Mathematically, there exists a fundamental connection between tensor
fields describing the nematic system with the underlying geometry and
topological characteristics of the space it lives in. This relationship can
formally described by the mathematical concept known as a fiber bundle
structure. In this paper, we introduce active nematic systems situated on
generalized Riemannian surfaces using a geometric language and the complex
line bundle formalism. Our focus is on understanding how nematic order
parameters relate to the intrinsic curvature and global topological properties
of the underlying space. We also present a straightforward and stable
computational methodology that aligns with this theoretical framework. We
expect the method to be a useful tool for practitioners studying active
nematic systems where the local curvature and global topology of the system is
important, such as in morphogenesis, and for exploring potential engineering
applications in soft robotics and morphing materials.

To characterize nematic fields, the prevailing approach in the existing
literature relies on the Landau-de Gennes $\boldsymbol{Q}$-tensor theory,
which describes the configuration of the nematic field on a surface via a
second-order symmetric tensor. However, in 2D there exists an alternative
representation through a homogeneous quadratic complex scalar function. This
approach offers distinct advantages because of the inherent constraints on its
degrees of freedom. Unlike the complex representation, the evolution of the
$\boldsymbol{Q}$-tensor needs to be confined within the realm of traceless and
symmetric matrices, which necessitates additional complications such as the
use of Lagrange multipliers during evolution [14, 15]. Furthermore, the
complex representation facilitates a natural extension to $k$-atic fields
(i.e., director fields with $k$-rotational symmetry) by simply raising the
exponent of the complex function to $k$. Recently, there has been significant
interest in modeling tissue dynamics through $k$-atic liquid crystal theory,
highlighting the practical relevance of this approach [16].

The mathematical formulation of nematohydrodynamics, or more broadly of
complex fluid flows, has been extensively explored in the literature [17, 18,
19, 20, 21]. Classic approaches primarily rely on an explicit approach,
employing coordinate-based expressions and standard vector calculus. They
often focus on the mechanics of specific physical objects in flow and root
their arguments in the principle of material frame invariance [21]. These
approaches become particularly challenging when attempting to translate
ordinary vector calculus identities to curved surfaces. However, in the realm
of differential geometry, there exists an abstract yet powerful geometric
language, exemplified by the Lie derivative, which functions independently of
coordinate systems. This language offers a natural way to comprehend various
physical phenomena and conservation laws in fluid dynamics [22]. Concepts such
as convected time derivatives in complex fluids can be effectively generalized
through the use of the Lie derivative, unifying the derivation process.
Despite its generality and utility, the adoption of this geometric language
remains relatively obscure in applied fluid mechanics and is not widely
recognized in the field of active nematics. In light of these challenges, our
objective is to bridge this knowledge gap and recast a minimal active nematic
system using the Lie derivative. This effort not only provides a unified
perspective on the various components of the active nematohydrodynamic system
but also, as we will demonstrate in the methods section, leads to a stable and
elegant numerical approach [23, 24, 25].

Similar to the Lie derivative, we also seek to establish a coordinate-free
framework for analyzing active nematic fields on curved surfaces by using the
mathematical structure called _fiber bundles_. This structure encapsulates the
concept of a family of tensor spaces parameterized by a base manifold,
providing a coordinate-free geometric language for describing the tensor
fields on the manifold as _sections_ of bundles. When each _fiber_ in the
bundle is a one-dimensional complex vector space equipped with a Riemannian
metric, which houses the $k$-atic director and velocity vector, this bundle is
considered a Hermitian complex line bundle. The theory of the complex line
bundle has been the subject of extensive research in both pure mathematics and
theoretical physics, particularly in modeling condensed matter phenomena [26].
The Levi-Civita connection, or covariant derivative, connects each fiber
within the tangent bundle. The curvature of the base manifold manifests from
the holonomy of this connection. The Levi-Civita connection on the tangent
bundle induces an associated connection onto the nematic bundle, which gives
rise to a nematic Bochner Laplacian that accommodates fractional topological
defects under the relaxation of nematic elasticity. The nuanced difference
between the nematic and vector Bochner Laplacians has often been overlooked or
inadequately explored in the existing literature. In addition to the
theoretical advancement it provides, this complex line representation also
facilitates the construction of an intrinsic complex-valued discrete Laplacian
matrix based on the local chart, instead of relying on a global
parametrization or an $\mathbb{R}^{3}$ representation in the embedded
coordinate. Overall, the formalism grounded in the complex line bundle theory
unveils a richer perspective on the representation of nematic—and more broadly
$k$-atic—fields, across arbitrary 2-dimensional manifolds.

We also show how one can formulate the Stokes equations through the Dolbeault
operators in the complex line bundle theory. In the smooth setting, the flow
of a viscous medium follows the Onsager variational principle. According to
this principle, the velocity of the viscous film minimizes a dissipation
measure known as the Rayleighian, which is a function of the strain rate
tensor. On a Riemannian manifold, the measurement of strain rate is determined
by the symmetric part of the covariant derivative, commonly referred to as the
Killing operator. With the complex structure, the covariant derivative can be
orthogonally decomposed into the holomorphic _Del_ and the anti-holomorphic
_Del bar_ operators, jointly known as the Dolbeault operators. We elucidate
that, under the incompressibility condition, the Laplacian-like operator that
determines the viscous stress on a Riemannian manifold is induced by the Del
bar operator. In addition to the Bochner Laplacian, this operator encompasses
contributions from the curvature of the manifold. This fact is acknowledged in
the literature, despite the results having been derived from the perspective
of Riemannian geometry [22, 27, 28]. While both derivations are equivalent in
this specific case, we believe that the more general formulation based on
complex line bundles can reveal new perspectives and potential for further
theoretical or computational advancements.

In the computational realm, several efforts have previously been undertaken to
solve active fluid dynamics on curved surfaces. These endeavors are primarily
founded on continuum formulations [29, 30, 31] or on the Lattice Boltzmann
method [32, 10]. In particular, we highlight the limitations of using the
streamfunction formulation for incompressible fluid mechanics on manifolds
with non-trivial topology without specialized modifications [33, 34, 35, 36,
37, 30]. However, one common challenge with velocity-pressure-based methods is
the representation of tangent vectors. In our method, we represent both the
velocity and director using an intrinsic complex-valued function alongside the
Levi-Civita connection. This approach eliminates the need for the commonly
used numerical technique of constraining the $\mathbb{R}^{3}$-valued velocity
to lie tangent to the manifold [30, 38]. Moreover, we draw parallels between
the numerical scheme for active nematics and classic problems in computational
geometry. For example, insightful analogies can be found from the resemblance
between the Stokes equations and the computation of Killing vector fields,
which aid in discerning symmetries within geometries [39, 40, 41, 42].
Similarly, there is a correspondence between the modeling of $k$-atic liquid
crystals with vector field design, which has been developed in applications
such as mesh and texture generation [43, 42, 44]. In the discretization we
present below, we integrate the finite element approximation with concepts
from discrete differential geometry (DDG). DDG adopts geometric arguments in
constructing geometric invariants, such as curvatures, allowing exact
transference of fundamental theorems—such as the Gauss-Bonnet and Poincaré-
Hopf theorems—from the smooth geometry to the discrete one.

The rest of the paper is structured as follows. We first introduce the
geometric formulation of active nematics in a continuous setting in Sec. 2.
This is followed by the development of a discretization approach that adheres
to the underlying geometric principles in Sec. 3. Results are discussed in
Sec. 4 , where we verify our numerical scheme through convergence tests and
demonstrate its computational efficiency, and where we present numerical
solutions for the dynamics of active nematics on diverse geometries,
highlighting the ability of our framework to capture the intricate dynamics
under various flow regimes and surface topologies. We conclude and discuss
potential extensions of this work in Sec. 5.

##  2 Theoretical formulation

An active fluid is a soft material whose constituent elements have the ability
to exert internal mechanical stresses. In systems consisting of microscopic
particles, like those found in many biological systems, inertia plays a
negligible role due to their small size and low velocities, resulting in very
small Reynolds numbers. In such cases, the Navier-Stokes equations governing
the hydrodynamics simplify to the Stokes equations. Here, the fluid flow is
driven by the active stresses induced by the constituent elements and is
dominated by viscous effects. To study the collective behavior of these active
fluids on large length scales, we use a continuum description of the particles
through an order parameter field. For example, a polar element can be
represented as a vector, denoted as $\boldsymbol{q}$. In active nematics,
microscopic constituents exhibit nematic symmetry, meaning there is an
equivalence relation between orientations $\boldsymbol{q}$ and
$-\boldsymbol{q}$. In 2D, this nematic order, denoted as $[\boldsymbol{q}]$,
can be equivalently represented by a special real $2\times 2$ matrix or by a
complex-valued function [26], as we explain in more detail below. The
configuration of the nematic field is governed by the nematodynamics equation,
which resembles an advection-diffusion equation acting on $[\boldsymbol{q}]$.
The diffusion part of nematodynamics stems from the alignment with neighboring
nematic elements through local steric interactions. When the nematic particles
are suspended in a fluid, they are also advected by the fluid flow and
corotate due to fluid deformations. The time scale separation between the
nematodynamics and fluid viscous dissipation allows for an effective model
based on a quasi-static coupling between nematodynamics and the Stokes
equations. When a significant amount of energy is inputted by the active
components, these systems can exhibit a wide range of phases, ranging from
well-ordered patterns to chaotic states. Because the system operates on a
curved Riemannian surface and our discretization method is devoid of global
coordinates, we first introduce the system of equations through an abstract
geometric language and then present the coordinate-based expression in a later
section.

###  2.1 Mathematical framework

In this section, we will establish the geometric framework required to
formulate the hydrodynamics of active nematics on a 2D Riemannian manifold. In
Secs. 2.1.1 and 2.1.2, we introduce the complex representation of a $k$-atic
director, specifically for the tangent vector and the nematic director. In
Sec. 2.1.3, we introduce the Lie derivative on tensor fields, which is crucial
for describing transport phenomena and deformations in fluid mechanics. The
Bochner Laplacian governs the relaxation dynamics of vector and nematic
fields. We start by introducing the Levi-Civita connection (i.e., covariant
derivative) on the tangent bundle in Sec. 2.1.4, which enables us to extend
the discussion to the nematic connection in Sec. 2.1.6 and the Bochner
Laplacian on each bundle in Section 2.1.7. In Sec. 2.1.5, through the
perspective of complex manifolds, we decompose the covariant derivative into
Dolbeault operators and establish their connections with fluid kinematics.
Local features like curvature and nematic singularities have an influence on
active nematodynamics on a manifold and are related to surface topology
through local-to-global theorems in differential geometry. In Sec. 2.1.8, we
detect and quantify charges associated with nematic topological defects by
measuring the total rotation of a director field along a closed loop using the
Levi-Civita connection. The curvature emerges as the angle holonomy of this
path integration.

####  2.1.1 Order parameter of $k$-atic director

We first introduce the general description for a $k$-atic director in the
plane, such as a vector ($k=1$) and a nematic director ($k=2$). Let $V$ be a
2D Euclidean vector space, and the multiplicative group action of complex
numbers, $\mathbb{C}^{\times}$, on $V$ is given by:

$\displaystyle\mathbb{C}^{\times}\times V\rightarrow
V,\quad(re^{\boldsymbol{\mathrm{i}}\theta},\boldsymbol{u})\mapsto\mathrm{Scale}{(\mathrm{Rotate}{(\boldsymbol{u},\theta)},r)}.$
(1)

A finite subgroup of rotation actions,
$N=\\{1,e^{\boldsymbol{\mathrm{i}}{2\pi}/{k}},\ldots,e^{\boldsymbol{\mathrm{i}}(k-1){2\pi}/{k}}\\}\subset\mathbb{C}^{\times}$,
defines an equivalence relation of $k$-rotational symmetry denoted by $\sim$.
For any $\boldsymbol{u},\boldsymbol{v}\in V$, we have:

$\displaystyle\boldsymbol{u}\sim\boldsymbol{v}\quad\mathrm{iff}\quad\exists\quad\hat{z}\in
N\quad\mathrm{s.t.}\quad\boldsymbol{v}=\hat{z}\boldsymbol{u}.$ (2)

This equivalence relation allows us to define a quotient set $V_{k}=V/N$ and
the map of taking the equivalence class $[\cdot]_{k}:V\rightarrow V_{k}$. Each
element $[\boldsymbol{q}]_{k}\in V_{k}$ is called a $k$-atic director. The
action of $\mathbb{C}^{\times}$ on $V$ naturally induces a canonical quotient
action $\mathbb{C}^{\times}_{k}=\mathbb{C}^{\times}/N$ on $V_{k}$, given by:

$\displaystyle\mathbb{C}_{k}\times V_{k}\rightarrow
V_{k},\quad(r^{k}e^{\boldsymbol{\mathrm{i}}k\theta},[\boldsymbol{u}]_{k})\mapsto[re^{i\theta}\boldsymbol{u}]_{k},$
(3)

where the rotational group is partitioned in $\theta$, with
$e^{\boldsymbol{\mathrm{i}}k\theta}\sim\hat{z}e^{\boldsymbol{\mathrm{i}}k\theta},\hat{z}\in
N$. Formally, we have the power map
$P:\mathbb{C}^{\times}\to\mathbb{C}^{\times}$,
$re^{\boldsymbol{\mathrm{i}}\theta}\mapsto
r^{k}e^{\boldsymbol{\mathrm{i}}k\theta}$, which is a homomorphism with kernel
being exactly $N$. By the first isomorphism theorem,
$\mathbb{C}^{\times}/\ker(P)\cong\mathrm{im}(P)$, the quotient group is
isomorphic to the image of the map,
$\mathbb{C}^{\times}_{k}\cong\mathrm{im}(P)=\mathbb{C}^{\times}$. If given a
unit basis vector $\boldsymbol{\mathrm{1}}\in V$ with
$|\boldsymbol{\mathrm{1}}|=1$, we can establish an isomorphism between complex
numbers $\mathbb{C}\ni\hat{q}$ and vectors $V\ni\boldsymbol{q}$ by
$\boldsymbol{q}=\hat{q}\boldsymbol{\mathrm{1}}=re^{\boldsymbol{\mathrm{i}}\theta}\boldsymbol{\mathrm{1}}\in
V$. Similarly, given a $k$-director basis $[\boldsymbol{\mathrm{1}}]_{k}\in
V_{k}$, we can establish the isomorphism between $\mathbb{C}\ni\hat{[q]}_{k}$
and its associated $k$-director
$V_{k}\ni[\boldsymbol{q}]_{k}=\hat{[q]}_{k}[\boldsymbol{\mathrm{1}}]_{k}=r^{k}e^{\boldsymbol{\mathrm{i}}k\theta}[\boldsymbol{\mathrm{1}}]_{k}=\hat{q}^{k}[\boldsymbol{\mathrm{1}}]_{k}$.
Through the chain rule, the _pushforward_ of the equivalence map, denoted as
$d[]_{k}|_{\boldsymbol{q}}$, maps an increment vector,
$\boldsymbol{q}^{\prime}\in T_{\boldsymbol{q}}V$, to an increment $k$-atic
director,
$[\boldsymbol{q}]^{\prime}_{k}=k[\boldsymbol{q}]_{k-1}\boldsymbol{q}^{\prime}\in
T_{[\boldsymbol{q}]_{k}}V_{k}$.

When focusing on a nematic field with $k=2$ and $N=\\{1,-1\\}$, an alternative
representation of the equivalence class is through a rank-$1$ symmetric tensor
$[\boldsymbol{q}]_{2}\cong\boldsymbol{q}\otimes\boldsymbol{q}\in\mathcal{Q}$,
where $\mathcal{Q}=\\{\,\boldsymbol{Q}\in
V\otimes_{\operatorname{symm}}V\mid\text{rank }\boldsymbol{Q}\leq 1\,\\}$. We
called the bijective map from the representation of complex numbers to the
matrix representation the Veronese map,

$\displaystyle\mathcal{V}:V_{2}\rightarrow\mathcal{Q},\quad[\boldsymbol{q}]_{2}\mapsto\boldsymbol{q}\otimes\boldsymbol{q}.$
(4)

For clarity, in the reminder of this paper, we will omit the subscript and
specify the bracket to denote a nematic equivalence (i.e., $[]=[]_{2}$). With
$[\boldsymbol{q}]\in V_{2}$,
$\boldsymbol{Q}=\mathcal{V}([\boldsymbol{q}])\in\mathcal{Q}$, the pushforward
of the Veronese map $d\mathcal{V}:T_{[\boldsymbol{q}]}V_{2}\rightarrow
T_{\boldsymbol{Q}}\mathcal{Q}$ maps an increment $[\boldsymbol{q}]^{\prime}$
in complex representation to the increment in matrix representation
$\boldsymbol{Q}^{\prime}$.

####  2.1.2 Vector field and nematic field on a Riemannian manifold

On a closed Riemannian 2-manifold $M$, $\partial M=\varnothing$, each point
$p\in M$ has a tangent space $T_{p}M$ that is a copy of the Euclidean plane:
$T_{p}M\cong V$. Tangent spaces at different points, $T_{p}M$, $T_{q}M$, for
$p\neq q$, are disjoint spaces, and the disjoint union of all tangent spaces
is referred to as the _tangent bundle_ , $TM=\bigcup_{p}T_{p}M$. Here the
tangent space $T_{p}M$ is also referred to as the _fiber_ of the _fiber
bundle_ $TM$ at a base point $p\in M$. A tangent vector field $\boldsymbol{q}$
is described as a _section_ of the tangent bundle,
$\boldsymbol{q}\in\Gamma(TM)$, which is an assignment of a fiber element
$\boldsymbol{q}_{p}\in T_{p}M$ at each point $p\in M$. We also denote a
section defined over a region $U\subset M$ as
$\boldsymbol{q}\in\Gamma_{U}(TM)$.

Similarly, we can describe a nematic field as a section of the nematic bundle
$L$, $[\boldsymbol{q}]\in\Gamma(L)$. Each fiber of the nematic bundle $L$ at
$p\in M$, $L_{p}$, is given by the equivalence relation,
$L_{p}\coloneqq[T_{p}M](=[T_{p}M]_{2})$ (c.f. Sec 2.1.1). In other words, the
nematic bundle $L=[TM]$ is obtained by taking fiberwise equivalence class of
the tangent bundle.

With a basis section for the tangent bundle,
$\boldsymbol{\mathrm{1}}\in\Gamma_{U}(TM),|\boldsymbol{\mathrm{1}}|=1$,
defined in a neighborhood $U\subset M$, we can construct the corresponding
nematic basis section for $L$ as $[\boldsymbol{\mathrm{1}}]\in\Gamma_{U}(L)$ .
Basis sections allow us to represent vector fields and nematic fields by
complex scalar functions, $\hat{q}\colon U\to\mathbb{C}$ and $\hat{[q]}\colon
U\to\mathbb{C}_{2}$.

Recall that a nematic vector $[\boldsymbol{q}]$ can also be represented by a
rank-1 symmetric tensor via the Veronese map of Eq. (4). We can similarly
represent a nematic field as a section of the symmetric rank-1 tensor bundle
as follows. The complex representation of the nematic field is isomorphic to a
tensor field $\boldsymbol{Q}=\boldsymbol{q}\otimes\boldsymbol{q}\in\Gamma(E)$,
where $E$ is the fiber bundle given by the set $\mathcal{Q}$ of rank $\leq 1$
in the symmetric tensor bundle $TM\otimes_{\rm symm}TM$. By populating the
fiber-wise Veronese map $\mathcal{V}:V_{2}\rightarrow\mathcal{Q}$, we get the
section-wise map $\mathcal{V}:\Gamma(L)\rightarrow\Gamma(E)$.

The Riemannian metric structure,
$\boldsymbol{g}|_{p}:T_{p}M\times_{\operatorname{symm}}T_{p}M\rightarrow\mathbb{R}$,
allows a positive-definite inner product
$\langle\cdot,\cdot\rangle=\boldsymbol{g}(\cdot,\cdot)$ between tangent vector
fields or nematic fields. For any two vector fields
$\boldsymbol{u}=\hat{u}\boldsymbol{\mathrm{1}}$,
$\boldsymbol{v}=\hat{v}\boldsymbol{\mathrm{1}}$ and two nematic fields
$[\boldsymbol{q}]=\hat{[q]}[\boldsymbol{\mathrm{1}}]\cong\boldsymbol{q}\otimes\boldsymbol{q}$,
$[\boldsymbol{p}]=\hat{[p]}[\boldsymbol{\mathrm{1}}]\cong\boldsymbol{p}\otimes\boldsymbol{p}$,
the inner products are defined using their complex representations,

$\displaystyle\langle\cdot,\cdot\rangle:\begin{cases}V\times
V\rightarrow\mathbb{R},&\langle\boldsymbol{u},\boldsymbol{v}\rangle\mapsto\Re(\hat{u}\overline{\hat{v}})\\\
V_{2}\times
V_{2}\rightarrow\mathbb{R},&\langle[\boldsymbol{q}],[\boldsymbol{p}]\rangle\mapsto\Re(\hat{[q]}\overline{\hat{[p]}})\\\
\mathcal{Q}\times\mathcal{Q}\rightarrow\mathbb{R},&\langle\boldsymbol{p}\otimes\boldsymbol{p},\boldsymbol{q}\otimes\boldsymbol{q}\rangle=\langle\boldsymbol{p},\boldsymbol{q}\rangle^{2}\mapsto\left(\Re(\hat{q}\overline{\hat{p}})\right)^{2},\end{cases}$
(5)

where $\overline{\cdot}$ denotes complex conjugation. Note that the inner
product for $\mathcal{Q}$ is the Frobenius inner product that can be extended
and denoted as $\boldsymbol{G}=\boldsymbol{g}\otimes\boldsymbol{g}:(V\otimes
V)\times(V\otimes V)\rightarrow\mathbb{R}$. A functional that maps a tangent
vector field to a scalar function, $\alpha:T_{p}M\rightarrow\mathbb{R}$, is
referred to as a _covector_ field. A natural _dual pairing_ of $\alpha$ with a
vector field $\boldsymbol{u}$ is denoted as
$\alpha\llbracket\boldsymbol{u}\rrbracket$. We denote this as
$\alpha\in\Gamma(T^{*}M)$, where $T^{*}M$ is known as the cotangent bundle,
which is the dual space to the tangent bundle. The metric provides a way to
associate a vector with its dual counterpart,
$\boldsymbol{g}|_{p}:T_{p}M\rightarrow T_{p}^{*}M$, where a vector
$\boldsymbol{u}$ is mapped to $\boldsymbol{g}(\boldsymbol{u},\cdot)$. In index
notation, this is referred to as _lowering_ an index, and it can be expressed
as $u_{i}=g_{ij}u^{j}$. Similarly, the Frobenius inner product
$\boldsymbol{G}|_{p}:(T_{p}M\otimes T_{p}M)\rightarrow(T_{p}^{*}M\otimes
T_{p}^{*}M)$ provides a way to map a $(2,0)$ tensor to a $(0,2)$ tensor, or
$P_{ij}=g_{im}g_{jn}P^{mn}$ in index notation. With metric, we can
subsequently denote the $L_{2}$ inner product of a vector or nematic field on
$M$ as
$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\cdot,\cdot\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}\coloneqq\int_{M}\langle\cdot,\cdot\rangle~{}dA:\Gamma(\circ)\times\Gamma(\circ)\rightarrow\mathbb{R}$.
The corresponding metric norm and $L_{2}$ norm are denoted as $|\cdot|$ and
$\|\cdot\|$.

####  2.1.3 Flow, Lie derivative, and advection

In differential geometry, the Lie derivative provides a generalized framework
for assessing changes in tensor fields encompassing functions with scalar,
vector, and covector values, along the flow defined by another vector field
[45, 46]. In the context of fluid mechanics and in particular complex fluids,
the Lie derivative generalizes the concept of directional derivative and
henceforth material derivative by capturing the dynamical behavior and
transformations of the constituent elements such as the nematic field on a
general manifold. By employing the Lie derivative, we can establish a
comprehensive definition of the advection transport phenomena. Notably, when
the Lie derivative is applied to a vector field, it recovers the well-known
upper-convected derivative. Applying the Lie derivative to the metric yields
information about the infinitesimal generation of flow velocity deformation.
This unveils the concept of viscous stress, arising from a fluid’s inherent
resistance to infinitesimal deviations from isometry. From a practical
standpoint, the Lie derivative perspective provides a unifying geometric
discretization that forgoes the use of coordinate-based expressions and
respects the underlying continuous structure.

A velocity field $\boldsymbol{u}\in\Gamma(TM)$ generates an instantaneous
flowmap $\varphi(t):M\rightarrow M$ such that the time derivative of the
flowmap is the velocity field, $\dot{\varphi}=\boldsymbol{u},\varphi(0)={\rm
id}_{M}$. The Lie derivative is defined as the rate of change of the so-called
pullback map $\varphi^{*}$ induced by this time-varying flowmap. Here, we
first motivate these concepts by applying the Lie derivative on a scalar
function. A pullback map for a scalar function pulls back a smooth function
defined on image $\varphi(M)$, $f\in C^{\infty}(\varphi(M))$, to a function
defined on the domain $M$, $\varphi^{*}f\in C^{\infty}(M)$, by precomposition,
$\varphi^{*}f=f\circ\varphi:M\rightarrow\varphi(M)\rightarrow\mathbb{R}$. As a
result, the Lie derivative on a scalar function leads to the familiar material
derivative commonly seen in continuum mechanics,

$\displaystyle\left.\frac{\partial}{\partial
t}\right|_{t=0}\varphi^{*}f\eqqcolon\frac{\partial f}{\partial
t}+\mathcal{L}_{\dot{\varphi}}f=\lim_{t\rightarrow
0}\frac{f_{\varphi(p)}-f_{p}}{t}=\dot{f}+d_{\boldsymbol{u}}f,$ (6)

where $d_{\boldsymbol{u}}f$ is the scalar directional derivative. The
directional derivative is an example of a natural dual pairing,
$d_{\boldsymbol{u}}f=df\llbracket\boldsymbol{u}\rrbracket\in\mathbb{R}$,
between the function map $df|_{p}:T_{p}M\rightarrow
T_{f(p)}\mathbb{R}\cong\mathbb{R}$, with a tangent vector $\boldsymbol{u}\in
T_{p}M$ without invoking the metric.

To define pullback for a general tensor field, we first introduce the
pushforward map or tangent map. The flowmap $\varphi$ induces a fiber-wise
pushforward map on the tangent space, $d\varphi|_{p}:T_{p}M\rightarrow
T_{\varphi(p)}M$. In differential geometry, the pushforward of a map sends a
tangent vector $\boldsymbol{q}$ in the domain to a tangent vector at the
image. In the context of continuum mechanics, the pushforward map is referred
to as the deformation gradient, where it maps a material fiber in the
Lagrangian coordinate to a fiber in the Eulerian coordinates. We therefore
conveniently denote the vector at $p$ deformed by the flowmap as
$d\varphi|_{p}(\boldsymbol{q})=\boldsymbol{F}\boldsymbol{q}$. Populating to
the fiber bundle, the pushforward map maps a tangent vector field to another,
$d\varphi:\Gamma(TM)\rightarrow\Gamma(TM)$. Because of the linearity of the
tangent map, $d\varphi(-\boldsymbol{q})=-d\varphi(\boldsymbol{q})$, we can
deduce that the deformation of the nematic field is equivalent to first
deforming its associated vector field and taking the nematic equivalence, or
$d\varphi([\boldsymbol{q}])=[d\varphi(\boldsymbol{q})]=[\boldsymbol{F}\boldsymbol{q}]$.
This is equivalent to the map on the tensorial representation of the nematic
field,
$d\varphi(\boldsymbol{q}\otimes\boldsymbol{q})=d\varphi(\boldsymbol{q})\otimes
d\varphi(\boldsymbol{q})=\boldsymbol{F}\boldsymbol{Q}\boldsymbol{F}^{\top}$.
Therefore we have the overloaded pushforward map defined over the vector field
and both representations of the nematic field,

$\displaystyle
d\varphi:\begin{cases}\Gamma(TM)\rightarrow\Gamma(TM),&d\varphi(\boldsymbol{q})=\boldsymbol{F}\boldsymbol{q}\\\
\Gamma(L)\rightarrow\Gamma(L),&d\varphi([\boldsymbol{q}])=[\boldsymbol{F}\boldsymbol{q}]\\\
\Gamma(E)\rightarrow\Gamma(E),&d\varphi(\boldsymbol{Q})=\boldsymbol{F}\boldsymbol{Q}\boldsymbol{F}^{\top}.\end{cases}$
(7)

The pullback map for a covector field pulls back a covector field on
$\varphi(M)$ to $M$, $\varphi^{*}:T^{*}_{\varphi(p)}M\rightarrow T^{*}_{p}M$,
which is the adjoint operator of the pushforward map. Therefore for covector
field $\boldsymbol{\alpha}\in\Gamma(T^{*}M)$, the pullback map is
$\varphi^{*}\boldsymbol{\alpha}\mapsto(d\varphi)^{\top}\boldsymbol{\alpha}=\boldsymbol{F}^{\top}\boldsymbol{\alpha}$.
Because the metric
$\boldsymbol{g}(\cdot,\cdot)\in\Gamma(T^{*}M\otimes_{\operatorname{symm}}T^{*}M)$
is a valence $(0,2)$ tensor, its pullback map follows the covector
transformation rule and distributes over the tensor product, or
$\varphi^{*}\boldsymbol{g}\mapsto\boldsymbol{F}^{\top}\boldsymbol{g}\boldsymbol{F}$.
To satisfy the homomorphism property of the pullback map on preserving the
dual pairing,
$\varphi^{*}(\alpha\llbracket\boldsymbol{q}\rrbracket)=\varphi^{*}\alpha\llbracket\varphi^{*}\boldsymbol{q}\rrbracket$,
a pullback map for a vector field is the inverse map of the pushforward
$\varphi^{*}|_{p}=(d\varphi|_{p})^{-1}=\boldsymbol{F}^{-1}:T_{\varphi(p)}M\rightarrow
T_{p}M$. This fiberwise pullback map pulls back a tangent vector in the
Eulerian coordinate to a tangent vector in the Lagrangian coordinate. Similar
to the pushforward map, we can populate the pullback map to sections of
bundles. With omission of the precomposition for non-scalar fields (e.g.,
${\boldsymbol{q}}\circ\varphi$ is denoted as $\boldsymbol{q}$), we can denote
pullback maps as

$\displaystyle\varphi^{*}:\begin{cases}C^{\infty}(M)\rightarrow
C^{\infty}(M),&\varphi^{*}f=f\circ\varphi\\\
\Gamma(G)\rightarrow\Gamma(G),&\varphi^{*}\boldsymbol{g}=\boldsymbol{F}^{\top}\boldsymbol{g}\boldsymbol{F}\\\
\Gamma(TM)\rightarrow\Gamma(TM),&\varphi^{*}\boldsymbol{q}=\boldsymbol{F}^{-1}\boldsymbol{q}\\\
\Gamma(L)\rightarrow\Gamma(L),&\varphi^{*}[\boldsymbol{q}]=[\boldsymbol{F}^{-1}\boldsymbol{q}]\\\
\Gamma(E)\rightarrow\Gamma(E),&\varphi^{*}\boldsymbol{Q}=\boldsymbol{F}^{-1}\boldsymbol{Q}\boldsymbol{F}^{-\top}.\end{cases}$
(8)

The Lie derivative measures the rate of change of this time-varying pullback
field at time equals zero,

$\displaystyle\left.\frac{\partial}{\partial
t}\right|_{t=0}\circ\varphi^{*}\eqqcolon\frac{\partial}{\partial
t}+\mathcal{L}_{\dot{\varphi}}:\begin{cases}\dot{f}+\mathcal{L}_{\boldsymbol{u}}f&=\lim_{t\rightarrow
0}\left(f_{\varphi(p)}-f_{p}\right)/t\\\
\mathcal{L}_{\boldsymbol{u}}\boldsymbol{g}&=\lim_{t\rightarrow
0}\left({\boldsymbol{F}^{\top}\boldsymbol{g}_{\varphi(p)}\boldsymbol{F}-\boldsymbol{g}_{p}}\right)/t\\\
\dot{\boldsymbol{q}}+\mathcal{L}_{\boldsymbol{u}}\boldsymbol{q}&=\lim_{t\rightarrow
0}\left({\boldsymbol{F}^{-1}\boldsymbol{q}_{\varphi(p)}-\boldsymbol{q}_{p}}\right)/t\\\
\dot{[\boldsymbol{q}]}+\mathcal{L}_{\boldsymbol{u}}[\boldsymbol{q}]&=\lim_{t\rightarrow
0}\left({\left[\boldsymbol{F}^{-1}\boldsymbol{q}_{\varphi(p)}\right]-[\boldsymbol{q}]_{p}}\right)/t\\\
\dot{\boldsymbol{Q}}+\mathcal{L}_{\boldsymbol{u}}\boldsymbol{Q}&=\lim_{t\rightarrow
0}\left({\boldsymbol{F}^{-1}\boldsymbol{Q}_{\varphi(p)}\boldsymbol{F}^{-\top}-\boldsymbol{Q}_{p}}\right)/t.\\\
\end{cases}$ (9)

Here, we assume that the metric is time-independent, $\dot{\boldsymbol{g}}=0$.
The Lie derivative of the metric,
$\mathcal{L}_{\boldsymbol{u}}\boldsymbol{g}$, is commonly known as the strain
rate tensor in continuum mechanics. It represents the rate of change of the
finite strain characterized by the Green-Lagrange strain tensor,
$\boldsymbol{F}^{\top}\boldsymbol{g}\boldsymbol{F}-\boldsymbol{g}$. It is
important to note that the rate of change of a pullback field generates an
increment-valued field. For instance,
$\mathcal{L}_{\boldsymbol{u}}\boldsymbol{q}$ is a $T_{\boldsymbol{q}}V$-valued
vector field. Consequently, the Lie derivative does not commute with nematic
equivalence, that is,
$\mathcal{L}_{\boldsymbol{u}}[\boldsymbol{q}]\neq[\mathcal{L}_{\boldsymbol{u}}\boldsymbol{q}]$.
Similarly, $\mathcal{L}_{\boldsymbol{u}}\boldsymbol{Q}$ is a tensor field
valued in $T_{\boldsymbol{q}}\mathcal{Q}$, and it preserves the rank-1
structure of the tensorial representation of nematics. When the pullback field
is constant, or ${\partial\varphi^{*}}/{\partial t}=0$, the section is
referred to as being _Lie advected_ by the flow,
$\dot{\varphi}=\boldsymbol{u}(t)$. In the scalar case, this recovers the
familiar material derivative. When the Lie derivative is applied to a vector
and nematic field, the deformation caused by the velocity gradient of the flow
will rotate and stretch the director fiber. This phenomenon is commonly
referred to as the upper-convected derivative or Oldroyd derivative.

In summary, the Lie derivative serves as a generalization of the concept of
advection from scalar functions to tensor fields. It allows us to extend the
idea of a directional derivative, originally defined for scalar functions, to
general tensor fields within a continuous deformation generated by a flow
represented by a vector field. This Lie derivative of the director field
remains well-defined even in the absence of the Riemannian metric. However, it
is crucial to note that, unlike directional derivatives, the Lie derivative
with respect to a _single_ vector is ill-defined due to the lack of a well-
defined deformation map. To extend the concept of single-vector directional
derivative on a manifold from scalar fields to vector or tensor fields,
additional structures such as the connection or parallel transport become
essential. The resulting derivative is known as the covariant derivative.

####  2.1.4 Levi-Civita connection on the tangent bundle

In this section, we cover the necessary background on the Levi-Civita
connection similar to other expositions on differential geometry [47]. On a
curved space, each tangent space is in general different from other tangent
spaces. The Levi-Civita connection introduces the concept of parallel
transport, allowing us to compare vectors at different points along a curved
manifold. By establishing this connection, we can define the covariant
derivative $\nabla$ that extends the concept of differentiation of a tensor
field from flat spaces to curved spaces.

Again, we first motivate the concept of covariant derivative using a scalar
function. Given a parametrized segment on $M$,
$\gamma(t):[-\epsilon,\epsilon]\rightarrow M$, with the starting point
$\gamma(0)=p$ and velocity $\gamma^{\prime}(0)=\boldsymbol{v}$, the
directional derivative of a scalar function $f:M\rightarrow\mathbb{R}$ with
the _single_ vector $\boldsymbol{v}$ is

$\displaystyle
d_{\boldsymbol{v}}f\coloneqq\left.\frac{d}{dt}f(\gamma(t))\right|_{t=0}=\lim_{t\rightarrow
0}\frac{f(\gamma(t))-f(p)}{t}.$ (10)

To extend the concept of directional derivatives, initially applicable to
scalar functions, to functions represented as tangent vectors
$\boldsymbol{u}\in\Gamma(TM)$, we must address the challenge of comparing
vectors in different tangent spaces: $\boldsymbol{u}(\gamma(t))$ in
$T_{\gamma(t)}M$ and $\boldsymbol{u}(p)$ in $T_{p}M$, where generally
$T_{\gamma(t)}M\neq T_{p}M$. To overcome this challenge, we employ the concept
of parallel transport, which is a one-parameter family of linear maps
dependent on the path $\gamma$, denoted as
$\Pi_{p\rightsquigarrow\gamma(t)}\colon T_{p}M\to T_{\gamma(t)}M$, with
$t\in[-\epsilon,\epsilon]$. For a 2-manifold $M$ equipped with a Riemannian
metric, we can intrinsically characterize parallel transport. Specifically,
when transporting along a geodesic segment $\gamma$,
$\Pi_{p\rightsquigarrow\gamma(t)}\boldsymbol{u}(p)$ maintains its length and
forms a constant angle with respect to the tangent vector
$\gamma^{\prime}(t)$. In cases where $\gamma$ is not a geodesic, the rate of
change of the angle between
$\Pi_{p\rightsquigarrow\gamma(t)}\boldsymbol{u}(p)$ and $\gamma^{\prime}(t)$
is determined by the negative geodesic curvature of $\gamma$ at $\gamma(t)$.
Using this notion of parallel transport, we define the directional covariant
derivative of $\boldsymbol{u}$ with respect to $\boldsymbol{v}$ as follows:

$\displaystyle\nabla_{\boldsymbol{v}}\boldsymbol{u}=\left.\frac{\nabla}{dt}\boldsymbol{u}(\gamma(t))\right|_{t=0}\coloneqq\lim_{t\rightarrow
0}\frac{\Pi_{\gamma(t)\rightsquigarrow
p}\boldsymbol{u}(\gamma(t))-\boldsymbol{u}(p)}{t},$ (11)

where $\Pi_{\gamma(\epsilon)\rightsquigarrow p}$ denotes the inverse of
$\Pi_{p\rightsquigarrow\gamma(\epsilon)}$. Note that the covariant derivative
is uniquely defined when $\boldsymbol{v}=\gamma^{\prime}(0)$ and is
independent of the choice of $\gamma$ as long as $\gamma(0)=p$ and
$\gamma^{\prime}(0)=\boldsymbol{v}$. The use of $\nabla$ in place of $d$
signifies the incorporation of infinitesimal parallel transport. The covariant
derivative, denoted as $\nabla\boldsymbol{u}$, is obtained by removing the
directional vector from the directional covariant derivative such that
$\nabla\boldsymbol{u}\llbracket\boldsymbol{v}\rrbracket\coloneqq\nabla_{\boldsymbol{v}}\boldsymbol{u}$
or, in index notation, as $v^{i}\nabla_{i}u^{j}$. The covariant derivative
maps a vector field to a vector-valued covector field, represented as
$\nabla:\Gamma(TM)\rightarrow\Gamma(TM\otimes T^{*}M)$.

Alternatively, when considering the surface being isometrically embedded in 3D
Euclidean space through an embedding function
$\boldsymbol{r}:M\rightarrow\mathbb{R}^{3}$, we can intuitively understand
parallel transport via the ambient 3D space. The embedding function induces a
pushforward map, $d\boldsymbol{r}|_{p}:T_{p}M\rightarrow\mathbb{R}^{3}$, which
realizes each tangent vector, $\boldsymbol{u}\in T_{p}M$, as an Euclidean
vector,
$d\boldsymbol{r}|_{p}\llbracket\boldsymbol{u}\rrbracket\in\mathbb{R}^{3}$. The
concept of parallel transport along a path $\gamma$ that connects point
$p=\gamma(0)$ to $\gamma(\epsilon)$ can be geometically described by rolling
the tangent plane at $\gamma(0)$ along the path $\gamma$ to $\gamma(\epsilon)$
against the surface without slipping. Because this rolling process is
typically well-understood in a discrete setting, it is applied to discretize
the covariant derivative in Sec. 3.3.

Recall that under a local basis section,
$\boldsymbol{\mathrm{1}}\in\Gamma_{U}(TM)$, where
$|\boldsymbol{\mathrm{1}}|=1$, a tangent vector field
$\boldsymbol{u}\in\Gamma_{U}(TM)$ can be represented as a complex field
$\hat{u}\colon U\to\mathbb{C}$, given by
$\boldsymbol{u}=\hat{u}\boldsymbol{\mathrm{1}}$. In this complex
representation, the covariant derivative can be expressed as

$\displaystyle\nabla\boldsymbol{u}=(d\hat{u})\boldsymbol{\mathrm{1}}+\hat{u}(\nabla\boldsymbol{\mathrm{1}})\eqqcolon\left((d+\boldsymbol{\mathrm{i}}\omega)\hat{u}\right)\boldsymbol{\mathrm{1}}.$
(12)

where
$d\hat{u}=d\Re{(\hat{u})}+\boldsymbol{\mathrm{i}}d{\Im{(\hat{u})}}\in\Gamma_{U}(T^{*}M;\mathbb{C})$
represents a complex-valued covector field and $\omega$ is a real-valued
covector field $\omega\in\Gamma_{U}(T^{*}M;\mathbb{R})$ defined by the
covariant derivative of the basis field
$\nabla\boldsymbol{\mathrm{1}}\eqqcolon\boldsymbol{\mathrm{i}}\omega\boldsymbol{\mathrm{1}}\in\Gamma_{U}(TM\otimes
T^{*}M)$. The components of $\omega$ are commonly referred to as the
Christoffel symbol.

####  2.1.5 Holomorphic structure with Levi-Civita connection

As being demonstrated repeatedly, a vector field, when equipped with a basis
field, can be viewed as a complex scalar field. Using the Levi-Civita
connection as the spatial derivative operator for this complex field, we can
distinguish its holomorphic and anti-holomorphic components by orthogonally
decomposing the covariant derivative into the _Del_ and _Del bar_ operators,
$\nabla=\partial^{\nabla}+\bar{\partial}^{\nabla}$,

$\displaystyle\begin{split}\bar{\partial}_{\boldsymbol{v}}^{\nabla}\boldsymbol{u}&\coloneqq\frac{\nabla_{\boldsymbol{v}}\boldsymbol{u}+\boldsymbol{\mathrm{i}}\nabla_{\boldsymbol{\mathrm{i}}\boldsymbol{v}}\boldsymbol{u}}{2}\\\
\partial_{\boldsymbol{v}}^{\nabla}\boldsymbol{u}&\coloneqq\frac{\nabla_{\boldsymbol{v}}\boldsymbol{u}-\boldsymbol{\mathrm{i}}\nabla_{\boldsymbol{\mathrm{i}}\boldsymbol{v}}\boldsymbol{u}}{2}.\end{split}$
(13)

Note that $\boldsymbol{\mathrm{i}}$ represents rotation of the vector by
$90^{\circ}$,
$\boldsymbol{\mathrm{i}}(\cdot)=\operatorname{Rotate}(\cdot,90^{\circ})$. A
vector field $\boldsymbol{u}$ is holomorphic if and only if
$\nabla\boldsymbol{u}\llbracket\boldsymbol{\mathrm{i}}\cdot\rrbracket=\boldsymbol{\mathrm{i}}\nabla\boldsymbol{u}\llbracket\cdot\rrbracket$,
or $\bar{\partial}^{\nabla}\boldsymbol{u}=0$. Its conjugate
$\partial^{\nabla}\boldsymbol{u}=0$ indicates $\boldsymbol{u}$ is
antiholomorphic.

It turns out that the holomorphic and anti-holomorphic parts play imprtant
roles when the vector field $\boldsymbol{u}$ represents the velocity field of
a surface flow. When $\boldsymbol{u}$ denotes the fluid velocity, the
holomorphic part, $\partial^{\nabla}\boldsymbol{u}$, quantifies the extent of
conformal deformation (involving infinitesimal rotation with local dilation)
of the fluid. The anti-holomorphic part,
$\bar{\partial}^{\nabla}\boldsymbol{u}$, corresponds to the symmetic and
trace-free part of $\nabla\boldsymbol{u}$. It is the strain rate tensor for
incompressible flow and measures the infinitesimal shearing that viscosity
resists. This concept will be applied in Sec. 2.2.1 when formulating the
Stokes equations on the manifold.

####  2.1.6 Levi-Civita connection on the nematic bundle

We can extend the Levi-Civita connection from the tangent bundle to the
nematic bundle. This is equivalent to defining a nematic covariant derivative,
denoted by $\nabla^{L}:\Gamma(L)\rightarrow\Gamma(L\otimes T^{*}M)$. Under a
local nematic basis $[\boldsymbol{\mathrm{1}}]\in\Gamma_{U}(L)$, the nematic
covariant derivative can be defined as:

$\displaystyle\nabla^{L}[\boldsymbol{q}]=d\hat{[q]}\boldsymbol{\mathrm{1}}+\hat{[q]}\nabla^{L}[\boldsymbol{\mathrm{1}}]\coloneqq\left((d+\boldsymbol{\mathrm{i}}2\omega)\hat{[q]}\right)[\boldsymbol{\mathrm{1}}],$
(14)

where $\omega$ is the covector field
$\boldsymbol{\mathrm{i}}\omega\boldsymbol{\mathrm{1}}\coloneqq\nabla\boldsymbol{\mathrm{1}}$
associated to the tangent bundle connection and basis. The nematic connection
$\nabla^{L}[\boldsymbol{\mathrm{1}}]\coloneqq\boldsymbol{\mathrm{i}}2\omega[\boldsymbol{\mathrm{1}}]$
reflects increased rotational speed due to nematic symmetry. This definition
is basis-independent and aligns with the algebraic construction of the Levi-
Civita connection on the tangent bundle,
$\nabla^{L}[\boldsymbol{q}]=2\boldsymbol{q}\nabla\boldsymbol{q}=\nabla(\boldsymbol{q}^{2})$.

The nematic covariant derivative $\nabla^{L}[\boldsymbol{q}]$ is consistent
with the standard covariant derivative for matrix-valued functions,
$\nabla:\Gamma(TM\otimes TM)\rightarrow\Gamma(TM\otimes TM\otimes T^{*}M)$,
when applied to
$\boldsymbol{Q}=\mathcal{V}([\boldsymbol{q}])\in\Gamma(E)\subset\Gamma(TM\otimes
TM)$. Like the Lie derivative in the advection equation (c.f. discussion on
Eqn. (9)), covariant derivatives also produce increment-valued nematic fields,
$\nabla^{L}[\boldsymbol{q}]$ and $\nabla\boldsymbol{Q}$. Via the pushforward
of the Veronase map,
$d\mathcal{V}|_{{[\boldsymbol{q}]}}:T_{{[\boldsymbol{q}]}}\Gamma(L)\rightarrow
T_{\mathcal{V}({[\boldsymbol{q}]})}\Gamma(E)$, one can confirm the equivalence
of these two derivatives,
$\nabla\boldsymbol{Q}=\nabla\boldsymbol{q}\otimes\boldsymbol{q}+\boldsymbol{q}\otimes\nabla\boldsymbol{q}=d\mathcal{V}|_{[\boldsymbol{q}]}\llbracket\nabla^{L}[\boldsymbol{q}]\rrbracket$.

####  2.1.7 Laplacian and diffusion

The concept of the covariant derivative naturally leads to a Laplacian
operator referred to as the Bochner Laplacian, which is the first variation of
the Dirichlet energy.

To begin, we introduce the Bochner Laplacian when it operates on vector
fields. In the context of liquid crystals in Euclidean space, the Dirichlet
energy of a vector field is also known as the one constant Oseen-Frank free
energy of a polar liquid crystal [48, 49]. The Dirichlet energy functional
maps a vector field to an energy measure,
$\mathcal{F}:\Gamma(TM)\rightarrow\mathbb{R}$, and can be expressed as
follows:

$\mathcal{F}(\boldsymbol{u})=\frac{1}{2}\|\nabla\boldsymbol{u}\|^{2}\overset{\partial
M=\varnothing}{=}\frac{1}{2}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\boldsymbol{u},\nabla^{*}\nabla\boldsymbol{u}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}\eqqcolon-\frac{1}{2}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\boldsymbol{u},\Delta\boldsymbol{u}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}},$
(15)

where the operator $\nabla^{*}:\Gamma(TM\otimes T^{*}M)\rightarrow\Gamma(TM)$
represents the $L_{2}$ adjoint of $\nabla$, defining the Bochner Laplacian as
$\Delta\coloneqq-\nabla^{*}\nabla:\Gamma(TM)\rightarrow\Gamma(TM)$. Note that
we use the negative semi-definite convention for the Laplacian operator.

Using the same procedure, the nematic covariant derivative is employed to
construct the nematic Dirichlet energy,
$\mathcal{F}^{L}\in\Gamma(L)\rightarrow\mathbb{R}$, which is expressed as:

$\displaystyle\mathcal{F}^{L}([\boldsymbol{q}])=\frac{1}{2}\|\nabla^{L}[\boldsymbol{q}]\|^{2}\overset{\partial
M=\varnothing}{=}\frac{1}{2}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}{[\boldsymbol{q}]},\nabla^{L*}\nabla^{L}[\boldsymbol{q}]\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}\eqqcolon-\frac{1}{2}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}{[\boldsymbol{q}]},\Delta^{L}[\boldsymbol{q}]\mathclose{\hbox{\set<EMAIL_ADDRESS>(16)

The variation of $\mathcal{F}^{L}$ leads to the definition of the nematic
Bochner Laplacian, denoted as
$\Delta^{L}\coloneqq-\nabla^{L*}\nabla^{L}:\Gamma(L)\rightarrow\Gamma(L)$.

Analogously, the standard Bochner Laplacian $\Delta:\Gamma(TM\otimes
TM)\rightarrow\Gamma(TM\otimes TM)$ defined for a general matrix field
$\boldsymbol{H}$ can be derived from the matrix Dirichlet energy,
$\mathcal{E}:\Gamma(TM\otimes TM)\rightarrow\mathbb{R}$, measured by the
Frobenius inner product,
$\mathcal{E}(\boldsymbol{H})=\frac{1}{2}\|\nabla\boldsymbol{H}\|^{2}=\frac{1}{2}\int_{M}\nabla_{i}H_{jk}\nabla^{i}H^{jk}dA$.
This Laplacian $\Delta$ can act on a nematic field
$\boldsymbol{Q}=\mathcal{V}([\boldsymbol{q}])$, but it differs from the
nematic Bochner Laplacian, $\Delta\boldsymbol{Q}\neq
d\mathcal{V}|_{[\boldsymbol{q}]}(\Delta^{L}[\boldsymbol{q}])$.

The reason for this difference is two-fold. First, nematic fields only
correspond to rank-1 matrices in the nematic subspace $E\subset TM\otimes TM$,
but the matrix Bochner Laplacian spans the entire $TM\otimes TM$. This results
in the matrix representation of nematic fields being an overparametrization
under this operator. This overparametrization can be demonstrated by
considering the matrix diffusion equation
$\dot{\boldsymbol{Q}}=\Delta\boldsymbol{Q}$. In general, a rank-1 initial
condition $\boldsymbol{Q}|_{t=0}\in\Gamma(E)$ does not remain rank-1 under
this diffusion (i.e., $\boldsymbol{Q}$ diverges from the nematic subspace
$E$). In contrast, $\dot{[\boldsymbol{q}]}=\Delta^{L}[\boldsymbol{q}]$ remains
within the nematic phase. Second, despite the equivalence of covariant
derivatives
$\nabla\boldsymbol{Q}=d\mathcal{V}|_{[\boldsymbol{q}]}\llbracket\nabla^{L}[\boldsymbol{q}]\rrbracket$
(cf. Sec. 2.1.6), their corresponding Dirichlet energies are distinct:
$\mathcal{E}(\boldsymbol{Q})\neq\mathcal{F}^{L}([\boldsymbol{q}])$. This is
because the pushforward of the Veronese map
$d\mathcal{V}_{[\boldsymbol{q}]}\colon T_{[\boldsymbol{q}]}L\to
T_{\boldsymbol{Q}}E$ does not isometrically map the intrinsic norm on
$T_{[\boldsymbol{q}]}L$ to the Frobenius norm for matrices on
$T_{\boldsymbol{Q}}E$.

As a side remark, it is possible to represent the nematic Dirichlet energy
$\mathcal{F}^{L}([\boldsymbol{q}])$ in terms of $\boldsymbol{Q}$ through an
equivalent functional $\mathcal{F}^{E}\colon\Gamma(E)\to\mathbb{R}$ such that
$\mathcal{F}^{E}(\boldsymbol{Q})=\mathcal{F}^{L}([\boldsymbol{q}])$. We can
achieve this by introducing a specific quadratic form,
$|\cdot|_{\mathcal{Q}}^{2}\coloneqq|\cdot|^{2}-2\det(\cdot)$, for
$\mathbb{R}^{2\times 2}$ matrices on $TM\otimes TM$. This quadratic form
induces a corresponding symmetric bilinear form
$\langle\cdot,\cdot\rangle_{\mathcal{Q}}=\langle\cdot,\cdot\rangle-2\det(\cdot,\cdot)$.
The equivalent nematic Dirichlet energy $\mathcal{F}^{E}$ relies on this
bilinear form and can be explicitly expressed in index notation:

$\displaystyle\mathcal{F}^{E}(\boldsymbol{Q})=\frac{1}{2}\int_{M}|\nabla\boldsymbol{Q}|_{\mathcal{Q}}^{2}~{}dA\equiv{\frac{1}{2}}\int_{M}g^{ij}(g_{km}g_{\ell
n}-\det(\boldsymbol{g})\epsilon_{km}\epsilon_{\ell
n})\nabla_{i}Q^{k\ell}\nabla_{j}Q^{mn}\,dA.$ (17)

Here, $g_{ij}$ are the components of the metric tensor $\boldsymbol{g}$ and
$\epsilon_{ij}$ represents the Levi-Civita permutation symbol.

In this paper, we adopt the more natural choices of
$\mathcal{F}^{L}([\boldsymbol{q}])$ and $\Delta^{L}$ to model the nematic
field. We use $\Delta^{E}\colon\Gamma(E)\to\Gamma(E)$ to denote the Laplacian
equivalent to $\Delta^{L}$ but acting on the matrix field representation. That
is, $\Delta^{E}$ satisfies the equivalence
$\Delta^{E}\boldsymbol{Q}=\Delta^{E}(\mathcal{V}[\boldsymbol{q}])=d\mathcal{V}\llbracket\Delta^{L}[\boldsymbol{q}]\rrbracket\in
T_{\boldsymbol{Q}}\Gamma(E)$ and is defined as,

$\displaystyle\Delta^{E}\boldsymbol{Q}\coloneqq-\frac{\delta\mathcal{F}^{E}({\boldsymbol{Q}})}{\delta{\boldsymbol{Q}|_{E}}}=-d\mathcal{V}\left\llbracket\frac{\delta\mathcal{F}({[\boldsymbol{q}]})}{\delta{[\boldsymbol{q}]}}\right\rrbracket={d\mathcal{V}\circ\Delta^{L}\circ\mathcal{V}^{-1}}(\boldsymbol{Q}).$
(18)

The subscript $|_{E}$ emphasizes the variation within the nematic subspace
$E$, rather than the entire $TM\otimes TM$.

####  2.1.8 Curvature, defects, and local-global geometry

The singularities within the nematic field, also referred to as topological
defects, are associated with points where the director field exhibits zeros
and discontinuities. These defects are commonly characterized by their index
or charge $\mathcal{Z}$. To identify these defects, we employ the covariant
derivative in closed-loop integration, a classic procedure also discussed in
other texts [47, 50]. Interestingly, poles and zeros of complex functions can
be naturally viewed as the singularies in director fields. Power functions,
$f(z)={(z-z_{0})}^{\mathcal{Z}}$, where $f:\mathbb{C}\rightarrow\mathbb{C}$,
are an archetypal demonstration of defects with a charge $\mathcal{Z}$ at
$z=z_{0}$, with $\mathcal{Z}$ being an integer or a fraction. On a manifold,
the Poincaré-Hopf theorem asserts that the total charge of defects is
determined by a topological invariant called the Euler characteristic of the
surface, thus establishing a fundamental connection between the nematic field
and the topology of the underlying surface.

The charge of a topological defect in a $k$-atic field $[\boldsymbol{q}]_{k}$
is determined by the cumulative turning angle that the director field
undergoes around the defect. Here we consider a defect at point $p\in U$,
where $U\subset M$ represents a sufficiently small disk that does not encircle
any other singularities. The boundary of the disk is denoted as
$\gamma=\partial U$. Locally in $U$ it is always possible to establish a
smooth and defect-free basis $[\boldsymbol{\mathrm{1}}]_{k}$ and subsequently
express the rotation form under the Levi-Civita connection, denoted as
$\eta\in\Gamma(T^{*}M)$, as

$\displaystyle\eta\coloneqq\frac{1}{|[\boldsymbol{q}]_{k}|^{2}}\left\langle\nabla[\boldsymbol{q}]_{k},\boldsymbol{\mathrm{i}}[\boldsymbol{q}]_{k}\right\rangle=\frac{1}{|[\boldsymbol{q}]_{k}|^{2}}\left\langle
d[\boldsymbol{q}]_{k}+\boldsymbol{\mathrm{i}}k\omega[\boldsymbol{q}]_{k},\boldsymbol{\mathrm{i}}[\boldsymbol{q}]_{k}\right\rangle=d^{\mathbb{S}^{1}}\arg(\hat{[q]}_{k})+k\omega.$
(19)

Note that the function $\arg$ necessitates a smooth basis to be well-behaved
and returns values in $\mathbb{S}^{1}$ rather than real numbers. Despite being
$\mathbb{S}^{1}$-valued by itself, the differential,
$d^{\mathbb{S}^{1}}\arg(\hat{[q]}_{k})$, can be defined as a real-valued
function. This is because, at a given point, this differential depends only on
the function’s values in a small neighborhood around the point. In such a
local context, a continuous $\mathbb{S}^{1}$-valued function can be unwrapped
to a real-valued function up to a $2\pi$-integer constant that does not affect
the differential.

The closed-loop integral of the turning rate $\eta$ defined in Eqn. (19) along
the parametrized boundary loop $\gamma(t)=\partial U(t)$, where $t\in[0,T]$
and $\gamma(0)=\gamma(T)$ is referred to as the cumulative turning angle of
the director field around a defect. The integral of the connection $\omega$ in
Eqn. (19), which is known as the holonomy angle, is equal to the negative of
the integrated Gaussian curvature: $\oint_{\gamma}\omega=-\int_{U}KdA$.
Therefore, the charge of the topological defect in a $k$-atic director field,
denoted as $\mathcal{Z}\in\mathbb{Z}/k$, admits the following relation:

$\displaystyle\mathcal{Z}_{p}([\boldsymbol{q}]_{k})\coloneqq\frac{1}{2\pi
k}\oint_{\gamma}d^{\mathbb{S}^{1}}\arg(\hat{[q]}_{k})=\frac{1}{2\pi
k}\oint_{\gamma}\eta+\frac{1}{2\pi}\int_{U}KdA.$ (20)

Note that the closed integral ensures that the orientation after integration,
$\arg(\hat{[q]}_{k})|_{t=T}$, differs only by an integer multiple of $2\pi$
from its starting value, $\arg(\hat{[q]}_{k})|_{t=0}$. As a result, the charge
is consistently an integer multiple of $1/k$ (i.e., an integer value for a
vector field and a multiple of $1/2$ for a nematic field). This quantized
property also implies that the charge of the defect remains invariant
regardless of variations in the size and shape of $\gamma$ and is solely the
property of the director field. As a special case,
$\oint_{\gamma}d^{\mathbb{S}^{1}}\arg(\hat{[q]}_{k})=0$ if $\gamma$ does not
enclose any topological defects of $\arg(\hat{[q]}_{k})$.

Owing to the Poincaré-Hopf Theorem, we can relax the local constraint on $U$.
This theorem relates the sum of charges at each point $p_{i}$ on the entire
manifold $M$ (with $\partial M=\varnothing$) to a fundamental topological
invariant known as the Euler characteristic $\chi$:

$\displaystyle\sum_{i}\mathcal{Z}_{p_{i}}([\boldsymbol{q}]_{k})=\chi(M)=2-2g,$
(21)

where $g$ represents the genus of the manifold. For instance, this theorem
establishes the absence of a nonvanishing continuous director field on a
topological 2-sphere. By replacing $U$ with manifold $M$, Eqn. (20) leads to
the Gauss-Bonnet Theorem, which states that the total Gaussian curvature is
also a topological invariant, $\int_{M}KdA=2\pi\chi(M)$.

In a nutshell, we can detect defects and calculate their charges by
establishing the cumulative turning angle through the use of the covariant
derivative. Local-to-global theorems, including the Gauss-Bonnet and Poincaré-
Hopf theorems, connect local properties like curvature of the surface and
charges of the nematic field to robust global invariants. Practically, as we
will elaborate in Sec. 3.5, these formulas translate to the discretization of
Gaussian curvature and the algorithmic computation of topological charges,
which satisfy the discrete analog of local-to-global theorems.

####  2.1.9 Expressing the Lie derivative as a covariant derivative

Under the Levi-Civita connection, the Lie derivative can be expressed in terms
of the covariant derivative. To demonstrate this, we denote the directional
derivative, a derivation associated with its corresponding vector field, as
$\boldsymbol{u}\coloneqq d_{\boldsymbol{u}}:C^{\infty}(M)\rightarrow
C^{\infty}(M)$. In index notation, it is commonly denoted as
$u^{i}\partial_{i}$. The pullback map respects the derivation structure; for
$\boldsymbol{q}\in\Gamma(TM)$, $f\in C^{\infty}(M)$, we have
$\varphi^{*}(d_{\boldsymbol{q}}f)=d_{\varphi^{*}\boldsymbol{q}}\varphi^{*}f$.
Consequently, the Lie derivative satisfies the Leibniz rule over derivation,
$\mathcal{L}_{\boldsymbol{u}}(d_{\boldsymbol{q}}f)=d_{\mathcal{L}_{\boldsymbol{u}}\boldsymbol{q}}f+d_{\boldsymbol{q}}\mathcal{L}_{\boldsymbol{u}}f$.
Since the Lie derivative of a scalar field is equivalent to the directional
derivative, we have $\mathcal{L}_{\boldsymbol{u}}\boldsymbol{q}(f)\coloneqq
d_{\mathcal{L}_{\boldsymbol{u}}\boldsymbol{q}}f=\mathcal{L}_{\boldsymbol{u}}(d_{\boldsymbol{q}}f)-d_{\boldsymbol{q}}\mathcal{L}_{\boldsymbol{u}}f=\boldsymbol{u}\boldsymbol{q}(f)-\boldsymbol{q}\boldsymbol{u}(f)\eqqcolon[\boldsymbol{u},\boldsymbol{q}](f)$,
known as the Lie bracket. The Lie bracket can in fact be expressed using the
covariant derivative of the Levi-Civita connection, referred to as the
torsion-free condition of the connection,
$T^{\nabla}\llbracket\boldsymbol{u},\boldsymbol{q}\rrbracket\coloneqq\nabla_{\boldsymbol{u}}\boldsymbol{q}-\nabla_{\boldsymbol{q}}\boldsymbol{u}-[\boldsymbol{u},\boldsymbol{q}]=0$.
Under the Levi-Civita connection, the Lie material derivative for vector
$\boldsymbol{q}$ can be expressed as

$\displaystyle\dot{\boldsymbol{q}}+\mathcal{L}_{\boldsymbol{u}}\boldsymbol{q}=\dot{\boldsymbol{q}}+\nabla_{\boldsymbol{u}}\boldsymbol{q}-\nabla_{\boldsymbol{q}}\boldsymbol{u},$
(22)

which is a familiar form of the upper-convected derivative, or Oldroyd
derivative for vectors [51]. The term $\nabla_{\boldsymbol{u}}\boldsymbol{q}$
represents the parallel transport components, while
$\nabla_{\boldsymbol{q}}\boldsymbol{u}$ accounts for the local deformation
caused by the differential flow velocity.

Through the pushforward of $k$-atic equivalence $d[]_{k}$, the Lie derivative
of the $k$-atic director field can be expressed as follows:

$\displaystyle\dot{[\boldsymbol{q}]}_{k}+\mathcal{L}_{\boldsymbol{u}}[\boldsymbol{q}]_{k}=k\boldsymbol{q}^{k-1}(\dot{\boldsymbol{q}}+\mathcal{L}_{\boldsymbol{u}}\boldsymbol{q}).$
(23)

When $k=2$, this is equivalent to the advection of the tensorial
representation of the nematic field. Applying the Leibniz rule of the Lie
derivative to the tensor product, we have
$\mathcal{L}_{\boldsymbol{u}}\boldsymbol{Q}=\mathcal{L}_{\boldsymbol{u}}\boldsymbol{q}\otimes\boldsymbol{q}+\boldsymbol{q}\otimes\mathcal{L}_{\boldsymbol{u}}\boldsymbol{q}$.
By rearranging the terms, we recover the operator for the upper-convected
derivative of the symmetric tensor $\boldsymbol{Q}$,

$\displaystyle\begin{split}\dot{\boldsymbol{Q}}+\mathcal{L}_{\boldsymbol{u}}\boldsymbol{Q}=\dot{\boldsymbol{Q}}+\nabla_{\boldsymbol{u}}\boldsymbol{q}\otimes\boldsymbol{q}+\boldsymbol{q}\otimes\nabla_{\boldsymbol{u}}\boldsymbol{q}-\nabla_{\boldsymbol{q}}\boldsymbol{u}\otimes\boldsymbol{q}-\boldsymbol{q}\otimes\nabla_{\boldsymbol{q}}\boldsymbol{u}=\dot{\boldsymbol{Q}}+\nabla_{\boldsymbol{u}}\boldsymbol{Q}-(\nabla\boldsymbol{u})\boldsymbol{Q}-\boldsymbol{Q}(\nabla\boldsymbol{u}^{\top}).\end{split}$
(24)

The dual pairing order follows the convention of matrix multiplication. In
index notation, this can be expressed as
$(\mathcal{L}_{\boldsymbol{u}}\boldsymbol{Q})^{ij}=u^{k}\nabla_{k}Q^{ij}-\nabla_{k}u^{i}Q^{kj}-Q^{ik}\nabla_{k}u^{j}$.

As demonstrated in Section 2.1.3, the strain rate tensor can be formulated as
the Lie derivative of the metric,
$\mathcal{L}_{\boldsymbol{u}}\boldsymbol{g}$. When expressing it in terms of
the covariant derivative acting on the Eulerian velocity, we can retrieve its
well-known form as the symmetric component of the velocity gradient. By
expressing
$\dot{\boldsymbol{F}}=\dot{(d\varphi)}=\varphi^{*}(\nabla\dot{\varphi})=(\nabla\boldsymbol{u})\boldsymbol{F}\in\Gamma(TM\otimes
T^{*}M)$, we can express the Lie derivative of the metric as the pullback of
the strain rate tensor,

$\displaystyle\mathcal{L}_{\boldsymbol{u}}\boldsymbol{g}=\frac{\partial}{\partial
t}(\boldsymbol{F}^{\top}\boldsymbol{g}\boldsymbol{F})=\dot{\boldsymbol{F}}^{\top}\boldsymbol{g}\boldsymbol{F}+\boldsymbol{F}^{\top}\boldsymbol{g}\dot{\boldsymbol{F}}$
$\displaystyle=\boldsymbol{F}^{\top}(\nabla\boldsymbol{u}^{\top}\boldsymbol{g}+\boldsymbol{g}\nabla\boldsymbol{u})\boldsymbol{F}=\varphi^{*}(\nabla\boldsymbol{u}^{\top}\boldsymbol{g}+\boldsymbol{g}\nabla\boldsymbol{u}).$
(25)

Since the flowmap, $\varphi:M\rightarrow M$, is instantaneously an identity
map, we can express
$\mathcal{L}_{\boldsymbol{u}}\boldsymbol{g}=\nabla\boldsymbol{u}^{\top}\boldsymbol{g}+\boldsymbol{g}\nabla\boldsymbol{u}\in\Gamma(T^{*}M\otimes_{\operatorname{symm}}T^{*}M)$,
which measures the infinitesimal deformation rate of a fluid element. This Lie
derivative defines the Killing operator,

$\displaystyle\mathcal{K}:\Gamma(TM)\rightarrow\Gamma(T^{*}M\otimes_{\operatorname{symm}}T^{*}M),\quad\mathcal{K}\colon\boldsymbol{u}\mapsto\frac{1}{2}\mathcal{L}_{\boldsymbol{u}}\boldsymbol{g}.$
(26)

A vector field $\boldsymbol{v}$ is referred to as a Killing vector field if
$\mathcal{K}\boldsymbol{v}=0$ (i.e., the flow it generates preserves the
metric).

###  2.2 Hydrodynamics of active nematic films

In this section, we will extend the mathematical framework from the previous
sections to describe the hydrodynamics of an active nematic film on a curved
surface. The system can be analyzed in two parts. First, we have the steady-
state Stokes equations, which describe how the fluid responds to the active
stress induced by the nematic configuration. Subsequently, the steady-state
fluid velocity acts as the driving force for the advection-diffusion equation
that governs the behavior of the nematic field. This equation is commonly
referred to as the nematodynamics equation. The system assumes quasi-static
coupling, meaning that the dissipation of the fluid at this length scale is
considered instantaneous.

####  2.2.1 Forced surface Stokes flow

We assume the active nematic is confined to a fluid interface whose flow is
governed by the overdamped Stokes equations. The theoretical foundation for
understanding such dissipation-driven systems can be provided by the Onsager
variational principle or the Helmholtz minimum dissipation theorem [52, 53,
54, 37, 55, 56]. The associated functional, known as the Rayleigh dissipation
function, captures the dissipation rate of the system, and Stokes flow emerges
as the stationary condition of this dissipation function. The fluid flow is
driven by the stresses exerted by the active nematic. In the present work, we
only account for the active stress and neglect passive stresses associated
with nematic elasticity. The active stress aligns with the orientation of the
nematic constituents and is represented by the stress tensor
$\boldsymbol{\sigma}=\alpha\boldsymbol{Q}\in\Gamma(TM\otimes_{\operatorname{symm}}TM)$,
which can be either contractile ($\alpha>0$) or extensile ($\alpha<0$) [57].
According to Onsager’s variational principle, the formulation of Stokes flow
involves minimizing the Rayleigh dissipation potential while satisfying the
incompressibility constraint,

$\displaystyle\min_{\boldsymbol{u}}\quad\mu\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{K}\boldsymbol{u},\mathcal{K}\boldsymbol{u}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}-\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\boldsymbol{u},\mathrm{div}^{\nabla}\boldsymbol{\sigma}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}\quad\text{subject
to}\quad\mathrm{div}~{}\boldsymbol{u}=0,$ (27)

where the first term
$\mu\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{K}\boldsymbol{u},\mathcal{K}\boldsymbol{u}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}$
quantifies the rate of viscous dissipation [27, 28, 36]. Here the divergence
operator, $\mathrm{div}\coloneqq-\mathrm{grad}^{*}:\Gamma(TM)\rightarrow
C^{\infty}$, is the negated $L_{2}$ adjoint of the gradient operator,
$\mathrm{grad}\coloneqq\boldsymbol{g}^{-1}d:C^{\infty}\rightarrow\Gamma(TM)$.

The negated adjoint of the covariant derivative is referred to as as the
covariant divergence, given by
$\text{div}^{\nabla}\coloneqq-\nabla^{*}:\Gamma(TM\otimes
TM)\rightarrow\Gamma({TM})$.

Note that the Killing operator can be factored out using the relation
$\mathcal{K}\boldsymbol{u}=\bar{\partial}^{\nabla}\boldsymbol{u}+\mathrm{div}(\boldsymbol{u})\boldsymbol{I}/2$,
where $\mathrm{div}(\boldsymbol{u})\boldsymbol{I}/2$ isolates the trace of the
strain rate tensor that characterizes fluid dilation (c.f. Sec 2.1.5).
Therefore, we can replace $\mathcal{K}$ with $\bar{\partial}^{\nabla}$ under
the incompressibility constraint and represent the constrained optimization
problem as a minimax problem:

$\displaystyle\begin{split}\min_{\boldsymbol{u}}\max_{p}~{}\mathcal{R}&=\mu\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\bar{\partial}^{\nabla}\boldsymbol{u},\bar{\partial}^{\nabla}\boldsymbol{u}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}-\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\boldsymbol{u},\mathrm{div}^{\nabla}\boldsymbol{\sigma}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}-\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}p,\text{div
}\boldsymbol{u}\mathclose{\hbox{\set<EMAIL_ADDRESS>(28)

In this expression, $p$ is the fluid pressure acting as a Lagrange multiplier
for incompressibility. On a closed surface with $\partial M=\varnothing$, the
stationary conditions with respect to velocity $\boldsymbol{u}$ and pressure
$p$ yield the incompressible, steady Stokes equations on a 2D Riemannian
manifold:

$\displaystyle\begin{split}2\mu~{}\bar{\partial}^{\nabla*}\bar{\partial}^{\nabla}\boldsymbol{u}+\mathrm{grad}~{}p-\alpha~{}\mathrm{div}^{\nabla}\boldsymbol{Q}&=0,\\\
\text{div }\boldsymbol{u}&=0.\end{split}$ (29)

The curvature contributes to the viscous term as
$-2\bar{\partial}^{\nabla*}\bar{\partial}^{\nabla}\boldsymbol{u}=(\Delta+K)\boldsymbol{u}$.
This result can be deduced from the known approach based on Riemannian
geometry that results in the relation
$-2\mathcal{K}^{*}\mathcal{K}\boldsymbol{u}=(\Delta+K+\mathrm{grad}~{}\mathrm{div})\boldsymbol{u}$
[22, 27, 28]. Additionally, in A, we provide an alternative derivation within
the framework of complex manifolds using complex differential forms.

It is worth noting that the presence of a Killing vector field indicates
symmetries in the geometry, such as the rigid body rotation on a rotationally
symmetric surface. In such case, the Stokes equations might lack full rank,
leading to non-unique solutions with its kernel represented by the Killing
vector field. In our study, we select the least-norm solution as the canonical
one to the Stokes equations. This choice is equivalent to projecting out all
modes of Killing fields in the $L_{2}$ sense. The treatment is justified by
the following. In theory, the $L_{2}$ Killing component of the velocity does
not contribute to energy dissipation, and as a result, it remains conserved in
the presence of any non-zero inertia. This is an instance of Noether’s
theorem, which states that the presence of continuous spatial symmetry,
represented by a Killing vector field, implies the conservation of momentum.
Therefore, in such ideal case, the correct solution is the least-norm solution
with the addition of the Killing mode in the initial condition. However, in
practical situations, particularly when considering an interfacial fluid
surface that separates two bulk viscous fluids, the Killing mode will
eventually dissipate energy through frictional interactions with the
substrate§§§Interactions with the bulk fluid are not considered in this
study.. In this case, the system should converges to the least-norm solution
over time. In Sec. 3.2, we will elaborate on the specific implementation of
the Killing projection.

####  2.2.2 Nematodynamics

Modeling nematodynamics involves studying the dynamics of nematic liquid
crystals, which possess both fluid-like and orientational order properties. In
this study, we are particularly interested in the sharply aligned limit of the
active nematic field. This sharply aligned phase arises when the microscopic
rotational diffusion of the nematic molecules is negligible [58]. Under these
conditions, the nematodynamics of the macroscopic order parameter becomes
identical to the dynamics of individual microscopic nematic filaments. A
commonly used approach is to utilize an advection-diffusion equation, where
advection accounts for the flow-induced alignment and rotation of the nematic
director field, while diffusion represents the smoothing of the director due
to random molecular interactions [19].

The relaxation or diffusion part of the nematodynamics equation involves the
nematic Bochner Laplacian established in Sec. 2.1.7. The constraint
$|[\boldsymbol{q}]|=1$ imposes the sharply-aligned phase. However, according
to the Poincaré-Hopf theorem, zeros of the director field are guaranteed
unless the system is on a torus with a vanishing Euler characteristic $\chi=0$
(cf. Sec. 2.1.8). As a result, the hard constraint $|[\boldsymbol{q}]|=1$ is
in general not permissible for a continuous vector field. To address this
constraint, we adopt a Ginzburg-Landau formulation and introduce a weak
constraint through a penalizing energy term
$\mathcal{G}([{\boldsymbol{q}}])=\mathcal{G}(|[\boldsymbol{q}]|)=(|[\boldsymbol{q}]|^{2}-1)^{2}/(4\varepsilon^{2})$,
where $\varepsilon$ is referred to as the coherence length scale around the
defects [26]. This combination of the Ginzburg-Landau potential and the one-
constant Frank-Oseen energy,
$\mathcal{G}([{\boldsymbol{q}}])+|\nabla^{L}[\boldsymbol{q}]|^{2}$, is
commonly referred to as the Leslie-Ericksen energy, whose variation induces a
Ginzburg-Landau modulated diffusion [59, 48, 60, 61]. Note that the
nematodynamics in the close vicinity of defects is not captured by the present
model and requires more detailed analysis [62].

The advection of nematic molecules is captured by a combination of the Lie
transportation and the “rotation-only” Lie transportation. Lie transportation,
as described in Sec. 2.1.3, models how the nematic director field is passively
transformed by the flow generated by a velocity field. The rotation-only Lie
transportation factors out the stretching part of the Lie transportation, and
thus a nematic director is only rotated while keeping its magnitude.

By decomposing the velocity gradient into symmetric
($\boldsymbol{E}=\mathcal{K}\boldsymbol{u}$, representing strain rate) and
skew-symmetric ($\boldsymbol{W}=\nabla\boldsymbol{u}-\boldsymbol{E}$,
representing vorticity) components, we can express the nematic Lie derivative
(cf. Eqn. (9) and Sec. 2.1.9) as
$\mathcal{L}_{\boldsymbol{u}}[\boldsymbol{q}]=2\boldsymbol{q}\mathcal{L}_{\boldsymbol{u}}\boldsymbol{q}$,
where
$\mathcal{L}_{\boldsymbol{u}}\boldsymbol{q}=\nabla_{\boldsymbol{u}}\boldsymbol{q}-\boldsymbol{E}\boldsymbol{q}-\boldsymbol{W}\boldsymbol{q}$.
The vorticity $\boldsymbol{W}$ and the parallel transport
$\nabla_{\boldsymbol{u}}$ contribute to infinitesimal rotational
transformations that preserve length, whereas the strain rate $\boldsymbol{E}$
affects the molecules by aligning them with the flow and causing axial
stretching or compression.

A rotational Lie derivative,
$\mathcal{L}^{\operatorname{rot}}_{\boldsymbol{u}}$, also known as the Jaumann
or corotational derivative, specifically focuses on the rotation effect
induced by the flow, (e.g.,
$\mathcal{L}^{\operatorname{rot}}_{\boldsymbol{u}}\boldsymbol{q}=\nabla_{\boldsymbol{u}}\boldsymbol{q}-\boldsymbol{W}\boldsymbol{q}$
for vector fields $\boldsymbol{q}$). Like the classical Lie derivatives (Sec.
2.1.3), the rotational Lie derivative
$\mathcal{L}^{\operatorname{rot}}_{\boldsymbol{u}}$ is defined through a
pullback operation. The only difference is that this pullback is only the
rotation component of the classical pullback operator. Recall Eqn. (8) that if
$\varphi$ is a flow map with deformation gradient $\boldsymbol{F}$, the
pullback of a vector field and respectively a nematic field is given by
$\varphi^{*}\boldsymbol{q}=\boldsymbol{F}^{-1}\boldsymbol{q}$,
$\varphi^{*}[\boldsymbol{q}]=[\boldsymbol{F}^{-1}\boldsymbol{q}]$ and
$\varphi^{*}\boldsymbol{Q}=\boldsymbol{F}^{-1}\boldsymbol{Q}\boldsymbol{F}^{-\top}$.
Now, let $\boldsymbol{R}$ denote the rotational component of $\boldsymbol{F}$
in the polar decomposition $\boldsymbol{F}=\boldsymbol{R}\boldsymbol{U}$.
Define the _rotational pullback_ $\varphi_{\operatorname{rot}}^{*}$ on vector
fields and nematic fields by
$\varphi_{\operatorname{rot}}^{*}{\boldsymbol{q}}\coloneqq\boldsymbol{R}^{-1}\boldsymbol{q}$,
$\varphi_{\operatorname{rot}}^{*}{[\boldsymbol{q}]}\coloneqq[\boldsymbol{R}^{-1}\boldsymbol{q}]$,
and
$\varphi_{\operatorname{rot}}^{*}\boldsymbol{Q}\coloneqq\boldsymbol{R}^{-1}\boldsymbol{Q}\boldsymbol{R}^{-\top}$.
Via this rotational pullback operator, define the Jaumann derivative as
$\mathcal{L}_{\boldsymbol{u}}^{\operatorname{rot}}\coloneqq{\frac{\partial}{\partial
t}}\circ\varphi_{\operatorname{rot}}^{*}$, where $\varphi$ is the
infinitesimal flow map generated by $\boldsymbol{u}$. In terms of the
gradients $\nabla\boldsymbol{u}$, $\boldsymbol{E}$, $\boldsymbol{W}$ of
$\boldsymbol{u}$,

$\mathcal{L}_{\boldsymbol{u}}^{\operatorname{rot}}\boldsymbol{q}=\nabla_{\boldsymbol{u}}\boldsymbol{q}-\boldsymbol{W}\boldsymbol{q},\quad\mathcal{L}_{\boldsymbol{u}}^{\operatorname{rot}}[\boldsymbol{q}]=2\boldsymbol{q}\mathcal{L}_{\boldsymbol{u}}^{\operatorname{rot}}\boldsymbol{q},\quad\mathcal{L}_{\boldsymbol{u}}^{\operatorname{rot}}\boldsymbol{Q}=\nabla_{\boldsymbol{u}}\boldsymbol{Q}-\boldsymbol{W}\boldsymbol{Q}+\boldsymbol{Q}\boldsymbol{W}.$
(30)

The combination of the Jaumann derivative with the traditional Lie derivative
in a weighted manner, denoted as
$\mathcal{L}^{\lambda}_{\boldsymbol{u}}\coloneqq\lambda\mathcal{L}_{\boldsymbol{u}}+(1-\lambda)\mathcal{L}_{\boldsymbol{u}}^{\operatorname{rot}}$,
$\lambda\in[0,1]$, allows for the modeling of spheroidal objects with finite
ellipticity. The parameter $\lambda\in[0,1]$ allows for a linear transition
from considering the full velocity gradient at $\lambda=1$ to only the spin
tensor at $\lambda=0$. In the literature, this non-dimensional parameter
$\lambda$ is referred to as the tumbling parameter or Bretherton’s constant,
which determines the dominance of alignment in extensional flow versus
tumbling in rotational flow [63, 64].

In addition, nematic molecules are inextensible when subjected to flow. During
advection, the magnitude of the order parameter is transported like a scalar
field without undergoing stretching. Therefore, we project out the stretching
component in the Lie derivatives $\mathcal{L}_{\boldsymbol{u}}^{\lambda}$.
Explicitly the projected Lie derivatives are given by

$\displaystyle\begin{split}&\mathcal{P}(\mathcal{L}^{\lambda}_{\boldsymbol{u}}\boldsymbol{q})\coloneqq\mathcal{L}^{\lambda}_{\boldsymbol{u}}\boldsymbol{q}-\frac{1}{|\boldsymbol{q}|^{2}}\langle\mathcal{L}^{\lambda}_{\boldsymbol{u}}\boldsymbol{q},\boldsymbol{q}\rangle\boldsymbol{q}=\nabla_{\boldsymbol{u}}\boldsymbol{q}-\boldsymbol{W}{\boldsymbol{q}}-\lambda\boldsymbol{E}{\boldsymbol{q}}+\frac{\lambda}{|\boldsymbol{q}|^{2}}\langle\boldsymbol{E}{\boldsymbol{q}},\boldsymbol{q}\rangle\boldsymbol{q},\\\
&\mathcal{P}^{L}(\mathcal{L}^{\lambda}_{\boldsymbol{u}}[\boldsymbol{q}])\coloneqq\mathcal{L}^{\lambda}_{\boldsymbol{u}}[\boldsymbol{q}]-\frac{1}{|[\boldsymbol{q}]|^{2}}\langle\mathcal{L}^{\lambda}_{\boldsymbol{u}}[\boldsymbol{q}],[\boldsymbol{q}]\rangle[\boldsymbol{q}]=2\boldsymbol{q}~{}\mathcal{P}(\mathcal{L}^{\lambda}_{\boldsymbol{u}}\boldsymbol{q}),\\\
&\mathcal{P}^{E}(\mathcal{L}^{\lambda}_{\boldsymbol{u}}\boldsymbol{Q})\coloneqq\mathcal{L}^{\lambda}_{\boldsymbol{u}}\boldsymbol{Q}-\frac{1}{|\boldsymbol{Q}|^{2}}\langle\mathcal{L}^{\lambda}_{\boldsymbol{u}}\boldsymbol{Q},\boldsymbol{Q}\rangle\boldsymbol{Q}=\nabla_{\boldsymbol{u}}\boldsymbol{Q}-(\lambda\boldsymbol{E}+\boldsymbol{W})\boldsymbol{Q}-\boldsymbol{Q}(\lambda\boldsymbol{E}-\boldsymbol{W})+\frac{2\lambda}{|\boldsymbol{Q}|}\langle\boldsymbol{E},\boldsymbol{Q}\rangle\boldsymbol{Q}.\end{split}$
(31)

Here we normalize $\boldsymbol{Q}$ based on the Frobenius matrix product
$\langle\cdot,\cdot\rangle$. It is equally valid to use the nematic norm
$\langle\cdot,\cdot\rangle_{\mathcal{Q}}$ introduced in Eqn. (17). Projections
based on both norms are equivalent because $\det(\boldsymbol{Q})=0$, meaning
$|\boldsymbol{Q}|=|\boldsymbol{Q}|_{\mathcal{Q}}$, and
$2\det(\boldsymbol{Q},\boldsymbol{Q}^{\prime})=\det(\boldsymbol{Q},\boldsymbol{Q})^{\prime}=\det(\boldsymbol{Q})^{\prime}=0$,
indicating
$\langle\boldsymbol{Q}^{\prime},\boldsymbol{Q}\rangle=\langle\boldsymbol{Q}^{\prime},\boldsymbol{Q}\rangle_{\mathcal{Q}}$.
In the literature, the generalized isometric advection equation based on Eqn.
(31) is known as the Jeffery equation [65].

By combining Jeffery advection with the Ginzburg-Landau diffusion, we obtain
the nematodynamics equations governing the evolution of the surface nematic:

$\displaystyle\begin{split}&\dot{[\boldsymbol{q}]}+\mathcal{P}^{L}(\mathcal{L}_{\boldsymbol{u}}^{\lambda}[\boldsymbol{q}])=\frac{1}{\eta}\left(\Delta^{L}[\boldsymbol{q}]-\frac{\delta\mathcal{G}}{\delta[\boldsymbol{q}]}\right),\\\
&\dot{\boldsymbol{Q}}+\mathcal{P}^{E}(\mathcal{L}_{\boldsymbol{u}}^{\lambda}\boldsymbol{Q})=\frac{1}{\eta}\left(\Delta^{E}\boldsymbol{Q}-d\mathcal{V}\left\llbracket\frac{\delta\mathcal{G}}{\delta[\boldsymbol{q}]}\right\rrbracket\right),\end{split}$
(32)

where

$\displaystyle\frac{\delta\mathcal{G}([\boldsymbol{q}])}{\delta[\boldsymbol{q}]}=\frac{1}{\varepsilon^{2}}(|[\boldsymbol{q}]|^{2}-1)[\boldsymbol{q}],$
(33)

and $\eta$ is the co-called rotational viscosity. Although not explicitly
apparent, the Gaussian curvature $K$ plays an implicit role in the Bochner
Laplacian through the Levi-Civita connection [22, 42].

####  2.2.3 System of equations

We summarize and non-dimensionalize the governing equations for the
hydrodynamics of an active nematic fluid film. Using the domain length scale,
denoted as $r$, as the characteristic length scale and the time scale
$\tau=\mu|\alpha|^{-1}$, we can scale all variables and differential
operators. This leads to the dimensionless system given by:

$(\Delta+K)\boldsymbol{u}-{\text{grad}}~{}p+\mathrm{div}^{\nabla}\boldsymbol{Q}=0,\quad\text{div}~{}\boldsymbol{u}=0,$
(34a)
$\dot{[\boldsymbol{q}]}+\mathcal{P}^{L}\mathcal{L}^{\lambda}_{\boldsymbol{u}}[\boldsymbol{q}]=\frac{1}{\mbox{\em
Pe}}\left(\Delta^{L}[\boldsymbol{q}]-\frac{1}{\epsilon^{2}}(|[\boldsymbol{q}]|^{2}-1)[\boldsymbol{q}]\right),$
(34b)

where $\lambda\in[0,1]$ is the tumbling parameter; $\mbox{\em Pe}=|\alpha|\eta
r^{2}\mu^{-1}$ is the active Peclet number, which measures the influence of
activity-driven Jeffery advection compared to nematic relaxation;
$\epsilon=\varepsilon/r$ characterizes the size of a defect core by comparing
the coherence length to the system size. Throughout the rest of the study, we
solve the non-dimensional version of the system. In the numerical results
presented below, we assume the limit $\epsilon\rightarrow 0$ and adopt a
normalization procedure to account for the Ginzburg-Landau term.

##  3 Discretization and algorithm

In this section, we will explain the construction of a discrete system for
modeling an active fluid on a 2D triangular manifold mesh. The mesh serves as
a connected and embedded representation of a continuous 2D manifold,
comprising vertices, edges, and faces.

When discretizing Eqn. (34) on a manifold, we need to discretize surface
differential operators. We begin by constructing discrete differential
operators for scalar functions using linear finite elements. Every vertex and
face is assigned with a tangent space, over which vector fields and nematic
fields are defined. Derivatives of vector and nematic fields require the Levi-
Civita connection. By leveraging the Levi-Civita connection, we obtain both
the vector and nematic Bochner Laplacian by constructing a discrete Dirichlet
energy and taking its variation. By a discrete version of the Gauss–Bonnet
Theorem, we also obtain Gaussian curvatures on vertices and faces.

Using these discrete differential operators, we solve Eqn. (34) by time
splitting (cf. Alg. 1). At each time step, the active force is calculated by
taking the covariant divergence of the stress tensor induced by the nematic
field. This active force gives rise to a fluid velocity field, which is
obtained by solving the discrete Stokes equations (Eqn. (34a)) through
Augmented-Lagrangian (AL) iteration. Next, the nematodynamics equation (Eqn.
(34b)) evolves the nematic field with the advection by the Stokes flow
velocity field, as well as the elastic diffusion. The advection and diffusion
terms are separated in a splitting procedure. We first compute the evolution
of the nematic field Lie-advected by the velocity using an explicit semi-
Lagrangian (sL) scheme. The one-step diffusion, based on the Bochner Laplacian
operator, is solved using the implicit Euler method to ensure unconditional
stability. The projection that removes the stretching component of the Lie
derivative and the constraint imposed by the Ginzburg-Landau term are
consolidated into a single normalization step. For completeness, we also
include a discrete counterpart of the theory of Sec. 2.1.8 about curvature and
defects. The resulting algorithm enables us to determine the locations and the
charges of the topological defects of the nematic field.

Algorithm 1 Hydrodynamics of active nematics

1:Inputs:

2: Surface $M$, curvature $K$, time step $\Delta t$

3:Initialize:

4: Nematic field $\boldsymbol{Q}=\mathcal{V}([\boldsymbol{q}])$

5:for $t=0$ to $T$ do

6: Stokes flow($\boldsymbol{Q}$)

7: $\boldsymbol{f}\leftarrow\mathrm{div}^{\nabla}\boldsymbol{Q}$

8:
$\boldsymbol{u},p\leftarrow\textsc{Solve}((\Delta+K)\boldsymbol{u}=\text{grad
}p-\boldsymbol{f},\text{div }\boldsymbol{u}=0)$ $\triangleright$ AL iteration
(Sec. 3.6)

9: Nematodynamics($\boldsymbol{u}$, $\boldsymbol{Q}$)

10:
$\boldsymbol{Q}\leftarrow\textsc{Solve}(\dot{\boldsymbol{Q}}+\mathcal{P}^{E}\mathcal{L}^{\lambda}_{\boldsymbol{u}}\boldsymbol{Q}=0)$
$\triangleright$ Explicit sL (Sec. 3.2)

11: $[\boldsymbol{q}]\leftarrow\mathcal{V}^{-1}(\boldsymbol{Q})$

12: $[\boldsymbol{q}]\leftarrow\textsc{Solve}(\dot{[\boldsymbol{q}]}=\mbox{\em
Pe}^{-1}\Delta^{L}[\boldsymbol{q}])$ $\triangleright$ Implicit Euler

13:
$\boldsymbol{Q}\leftarrow\mathcal{V}(\textsc{Normalize}([\boldsymbol{q}]))$

14: $t\leftarrow t+\Delta t$

15:end for

In summary, the method presented in this work offers several distinctive
features that contribute to the effectiveness and versatility of the method in
modeling and analyzing the behavior of active nematics on curved surfaces:

  * •

The calculations are carried out on a 2D triangular mesh, rather than a
volumetric grid, allowing for efficient and accurate representation of the
system.

  * •

Tensor fields and differential operators are intrinsically represented using
the discrete complex line bundle. This representation, in contrast to the
common approach of applying extra tangentiality constraints on Euclidean
tensors and operators, allows better accuracy and lower computational cost.

  * •

A generalized semi-Lagrangian method is employed to handle the Lie derivative
for both vector and tensor quantities, thereby providing a general procedure
for the advection of $k$-atic directors.

  * •

We utilize isomorphic representations of the nematic field through the
Veronese map. This allows us to adapt to the more convenient representation at
each stage of the algorithm.

###  3.1 Discrete manifold and choice of variable space

In this section, our goal is to establish notations for representing a 2D
manifold using a discrete mesh and navigating its local topology. Our focus
here lies in developing a common language for tensor calculus on a discrete
manifold, which forms the basis for discretizing the system of equations.

####  3.1.1 Representation and navigation of discrete 2D manifolds

Figure 1: Schematics illustrating key components of a 2D triangular mesh: A)
Voronoi dual mesh is constructed by circumcenters of the triangular mesh,
defining the weight of each primal element. B) Tangent space associated with
the faces of the mesh, denoted as $T_{f}M$. C) Vertex tangent space $T_{v}M$
via the geodesic polar map.

A closed, oriented triangular mesh $M$ of arbitrary genus $g$ consists of the
sets of vertices, edges, and faces, denoted as $\\{\mathfrak{V,E,F}\\}$. Let
$|\mathfrak{V}|$, $|\mathfrak{E}|$, and $|\mathfrak{F}|$ denote the numbers of
vertices, edges, and faces, respectively. To indicate the local topology of
the mesh, we use the symbols $\prec$ and $\succ$, which signify the relation
“incident to”. The orientation of the symbol indicates the dimensionality
relationship. For example, for incident $v\in\mathfrak{V}$,
$e\in\mathfrak{E}$, and $f\in\mathfrak{F}$, we write $v\prec e\prec f$. When
using the summation sign, the iterating variable is placed first. For example,
$\sum_{f\succ v}$ means summing over faces $f$ that are adjacent to a given
vertex $v$. The directed edge vector associated with the edge
$e\in\mathfrak{E}$ in face $f\in\mathfrak{F}$, where $e\prec f$, is denoted as
$\boldsymbol{l}_{e\prec f}$. This vector is positively-oriented (i.e.,
counterclockwise) in face $f$. In particular, this edge has opposite
orientations in the two faces that share it. Similarly we can also denote the
edge vector associated with a vertex as $\boldsymbol{l}_{v\prec f}$, which has
the positive orientation of pointing outward from the vertex. Lastly, the
angle of the interior corner of a triangle $f$ at a vertex $v\prec f$ is
denoted by $\angle_{v\prec f}$, $\angle_{f\succ v}$, $\angle_{e\prec f}$, or
$\angle_{f\succ e}$, where $e$ is the edge opposite to vertex $v$ across
triangle $f$.

Using the circumcenters of the mesh faces, we can construct a _Voronoi dual
mesh_. In this dual mesh, each edge, known as the dual edge (illustrated by
Fig. 1A), orthogonally bisects the corresponding primal edge in the primal
mesh. Note that the dual edge is intrinsically straight, but it may be kinked
when the mesh is embedded due to the dihedral angle formed by neighboring
faces. We denote the directed dual edge as $\boldsymbol{l}^{*}_{e\prec
f}=w_{e}\boldsymbol{\mathrm{i}}\boldsymbol{l}_{e\prec f}$, where
$w_{e}={|\boldsymbol{l}_{e}|}/{|\boldsymbol{l}^{*}_{e}|}$ is the length ratio
between the dual edge and the primal edge referred to as the _cotangent
weight_ ¶¶¶since $w_{e}={\frac{1}{2}}(\cot(\angle_{e\prec
f_{1}})+\cot(\angle_{e\prec f_{2}}))$, where $f_{1}\neq f_{2}\succ e$ are the
two incident faces of $e$., and $\boldsymbol{\mathrm{i}}$ representes
$90^{\circ}$ counter-clockwise rotation in face $f$. Analogous to the
continuous notation that $C^{\infty}(M)$ denotes the space of continuous
functions on $M$, we use
$C_{\mathfrak{V}}(M)=\\{h\colon\mathfrak{V}\to\mathbb{R}\\}$ to denote the
space of functions defined over vertices. The function evaluation at vertex
$v\in\mathfrak{V}$ is denoted as $h_{v}$.

The notion of a tangent space to a mesh includes the tangent space $T_{v}M$ at
a vertex $v$ and the face tangent space $T_{f}M$ at a face $f$. A vector field
defined on faces is a _face section_ of the tangent bundle,
$\Gamma_{\mathfrak{F}}(TM)$, and similarly, a vector field defined on vertices
is a _vertex section_ denoted as $\Gamma_{\mathfrak{V}}(TM)$. A vector basis
section is represented by $\boldsymbol{\mathrm{1}}\in\Gamma(TM)$, where the
basis vector at a vertex $v$ is denoted as $\boldsymbol{\mathrm{1}}_{v}$, and
at a face $f$ as $\boldsymbol{\mathrm{1}}_{f}$. Given a basis section, a
vector field can be represented by a complex-valued function. For example, a
face-based velocity field, $\boldsymbol{u}\in\Gamma_{\mathfrak{F}}(TM)$, can
be represented as $\hat{u}:\mathfrak{F}\rightarrow\mathbb{C}$. The same
representation holds for a nematic field with its basis section denoted as
$[\boldsymbol{\mathrm{1}}]\in\Gamma(L)$. Specifically, on the face tangent
space $T_{f}M$, we start by selecting an arbitrary positively-oriented
normalized edge vector (e.g., the first edge $e_{1}$ in the data structure) as
the basis, denoted as $\boldsymbol{\mathrm{1}}_{f}=\boldsymbol{l}_{e_{1}\prec
f}/|\boldsymbol{l}_{e_{1}\prec f}|\in T_{f}M$ and
$[\boldsymbol{\mathrm{1}}]_{f}=[\boldsymbol{\mathrm{1}}_{f}]\in L_{f}$ (cf.
Fig. 1B). Then, we can associate all vectors and nematic director in this
tangent space with a complex number. In particular, we represent the phase
angle of the oriented edge $\boldsymbol{l}_{e\prec f}=\hat{l}_{e\prec
f}\boldsymbol{\mathrm{1}}_{f}$ as $\arg(\boldsymbol{l}_{e\prec
f})\coloneqq\arg(\hat{l}_{e\prec f})$. This phase angle is used to compute the
discrete Levi-Civita connection in Sec. 3.1.2.

On a vertex tangent space $T_{v}M$, which is cone-shaped, we can construct an
intrinsic representation referred to as the _geodesic polar map_ (cf. Fig.
1C). To achieve this, we first need to flatten the cone. However, due to
inherent curvature, a cone generally cannot be flattened perfectly due to
angle defects or excesses. In other words, the sum of all corner angles of the
faces incident to the vertex does not equal $2\pi$ (i.e., $\sum_{f\succ
v}\angle_{f\succ v}\neq 2\pi$). To perform an intrinsic projection, we scale
the corner angles by the ratio $r_{v}=(2\pi)/\sum_{f\succ v}\angle_{f\succ
v}$. By applying this scaling factor to all interior angles at vertex $v$, we
obtain a _geodesic angle_ for each corner, denoted as $\tilde{\angle}_{f\succ
v}=r_{v}\angle_{f\succ v}$. Similar to the face tangent space, we choose an
arbitrary positively-oriented edge vector $\boldsymbol{l}_{e_{1}\prec v}$ at
the vertex as the vector basis and its corresponding nematic basis,
represented by $\boldsymbol{\mathrm{1}}_{v}\coloneqq\boldsymbol{l}_{e_{1}\prec
v}/|\boldsymbol{l}_{e_{1}\prec v}|\in T_{v}M$ and
$[\boldsymbol{\mathrm{1}}]_{v}=[\boldsymbol{\mathrm{1}}_{v}]\in L_{v}$. We
accumulate the geodesic corner angles and obtain the phase angle
$\arg(\boldsymbol{l}_{e\succ v})=\arg(\hat{l}_{e\succ v})$ for each edge $e$
in this vertex polar plane. Consequently, the edge vector in $T_{v}M$ can be
represented by a complex number $\boldsymbol{l}_{e\succ v}=\hat{l}_{e\succ
v}\boldsymbol{\mathrm{1}}_{v}$.

####  3.1.2 Levi-Civita connection

We equip the discrete manifold with a Levi-Civita connection (cf. Sec. 2.1.4)
by locally constructing the discrete parallel transport. This includes the
parallel transport from the tangent space of one vertex to its neighboring
vertex, from a face to its neighboring face, as well as from a vertex to its
neighboring face.

Figure 2: Schematics depicting intrinsic parallel transport on a 2D triangular
mesh: A) face-face transportation by flattening the dihedral angle within a
two-triangle stencil. B) vertex-vertex transportation preserves angles along
geodesic paths within the geodesic polar map.

The face-face and vertex-vertex transportations can be defined intrinsically
(i.e., invariant under isometric deformation). Therefore, parallel transport
of a vector $\boldsymbol{u}$ from face $f_{i}$ to face $f_{j}$ along dual edge
$e^{*}$, $e\prec f_{i},f_{j}$, can be conceptualized as if the dihedral angle
at $e$ were flattened (cf. Fig. 2A). In this context, parallel transport
simplifies to a basis change across face $f_{i}$ to face $f_{j}$,

$\displaystyle\Pi_{i\rightsquigarrow j}:T_{f_{i}}M\rightarrow
T_{f_{j}}M,\quad\Pi_{i\rightsquigarrow
j}(\boldsymbol{u}_{i})=\Pi_{i\rightsquigarrow
j}(\hat{u}_{i}\boldsymbol{\mathrm{1}}_{i})=e^{-\boldsymbol{\mathrm{i}}\Omega_{e^{*}}}\hat{u}_{i}\boldsymbol{\mathrm{1}}_{j},$
(35)

where the angle-valued connection along the dual edge,
$\Omega_{e^{*}}=\arg(\boldsymbol{l}_{e\prec
f_{j}})-\arg(\boldsymbol{l}_{e\prec f_{i}})-\pi$, accounts for the basis
change,
$\boldsymbol{\mathrm{1}}_{j}=e^{\boldsymbol{\mathrm{i}}\Omega_{e^{*}}}\boldsymbol{\mathrm{1}}_{i}$
(cf. Fig. 2A and the continuous theory discussed around Eqn. (12)).

Analogous to the continuous theory that parallel transport can be achieved by
preserving the angle along a geodesic path (cf. Sec. 2.1.4), the vertex-vertex
transportation is intrinsically defined by maintaining the angle with respect
to the connecting edge $e$ between neighboring vertices $v_{i},v_{j}\prec e$,
illustrated by Fig. 2B. As seen in the continuous theory (cf. Sec. 2.1.6), the
nematic field requires amplifying the Levi-Civita connection on tangent bundle
$\Omega$ to $2\Omega$. Subsequently, we can express the parallel transport of
a nematic director $[\boldsymbol{q}]_{i}\in L_{v_{i}}$ to $L_{v_{j}}$ along
$\boldsymbol{l}_{e}$ as

$\displaystyle\Pi_{i\rightsquigarrow j}:L_{v_{i}}\rightarrow
L_{v_{j}},\quad\Pi_{i\rightsquigarrow
j}([\boldsymbol{q}]_{i})=\Pi_{i\rightsquigarrow
j}(\hat{[q]}_{i}[\boldsymbol{\mathrm{1}}]_{i})=e^{-\boldsymbol{\mathrm{i}}2\Omega_{e}}\hat{[q]}_{i}[\boldsymbol{\mathrm{1}}]_{j}.$
(36)

Here, the primal edge connection, $\Omega_{e}=\arg(\boldsymbol{l}_{e\succ
v_{j}})-\arg(\boldsymbol{l}_{e\succ v_{i}})-\pi$, quantifies the basis change
across $T_{v_{i}}M$ and $T_{v_{j}}M$ (cf. Fig. 2B).

The vertex-face transportation is achieved using the $\mathbb{R}^{3}$
embedding. Given the embedding of the mesh and the realization of the tangent
vector in $\mathbb{R}^{3}$, we can utilize a 3D dihedral rotation to map from
one vertex tangent space to a neighboring face tangent space and vice versa
(cf. the continuous theory in Sec. 2.1.4). To illustrate, consider a vertex
unit normal, $\boldsymbol{n}_{v}\in T_{v}^{\perp}M$,$|\boldsymbol{n}_{v}|=1$,
and a face unit normal, $\boldsymbol{n}_{f}\in T_{f}^{\perp}M$,
$|\boldsymbol{n}_{v}|=1$. The dihedral rotation from $\boldsymbol{n}_{v}$ to
$\boldsymbol{n}_{f}$ is defined as the rotation $\boldsymbol{R}_{vf}$ of with
angle $\cos^{-1}\langle\boldsymbol{n}_{v},\boldsymbol{n}_{f}\rangle$ about the
axis
$\boldsymbol{l}={\boldsymbol{n}_{v}\times\boldsymbol{n}_{f}}/{|\boldsymbol{n}_{v}\times\boldsymbol{n}_{f}|}$.
Define $\Pi_{v\rightsquigarrow
f}\coloneqq\boldsymbol{R}_{vf}:T_{v}M\rightarrow T_{f}M$. Note that the
definition of $\boldsymbol{n}_{v}$ is not unique. In practice, we use an
angle-weighted vertex normal $\boldsymbol{n}_{v}\coloneqq(\sum_{f\succ
v}\tilde{\angle}_{f\succ v}\boldsymbol{n}_{f})/(2\pi)$.

####  3.1.3 Space of system variables and operator mappings

In this section, we specify the choice of space where the discretization
variable resides, and the mapping spaces imposed by the differential
operators. In Secs. 3.2-3.4, we will recapitulate the specific definitions of
the operators.

As we mentioned in the previous section, it is convenient to define tangent
spaces over both the vertices and faces of the mesh. In our approach, we
discretize the Stokes equations on the faces of the mesh and the
nematodynamics equation on the vertices. Specifically, we work with a face-
based fluid velocity denoted as $\boldsymbol{u}\in\Gamma_{\mathfrak{F}}(TM)$
and a vertex-based nematic director represented as
$[\boldsymbol{q}]\in\Gamma_{\mathfrak{V}}(L)$, alongside the corresponding
matrix representation
$\boldsymbol{Q}=\boldsymbol{q}\otimes\boldsymbol{q}\in\Gamma_{\mathfrak{V}}(E)$.

Subsequently, we develop the discrete differential operators associated with
these choices. In the context of the Stokes equations, we construct the
Bochner Laplacian, denoted as $\Delta:\Gamma(TM)\rightarrow\Gamma(TM)$,
utilizing the face-face Levi-Civita connection. The gradient operator,
$\mathrm{grad}:C_{\mathfrak{V}}(M)\rightarrow\Gamma_{\mathfrak{F}}(TM)$, maps
scalar vertex data to vector data on faces. The divergence operator, the
adjoint of the gradient, maps vectors from faces to scalar measurements at
vertices, $\mathrm{div}:\Gamma_{\mathfrak{F}}(TM)\rightarrow
C_{\mathfrak{V}}(M)$. Similarly, in the context of the nematodynamics
equation, we establish the nematic Bochner Laplacian as
$\Delta^{L}[\boldsymbol{q}]:\Gamma_{\mathfrak{V}}(L)\rightarrow\Gamma_{\mathfrak{V}}(L)$,
using the vertex-vertex Levi-Civita connection. We also introduce a vertex-
based Lie derivative, denoted as
$\mathcal{L}_{\boldsymbol{u}}:\Gamma_{\mathfrak{V}}(E)\rightarrow\Gamma_{\mathfrak{V}}(E)$,
to advect the nematic director.

###  3.2 Semi-Lagrangian Lie advection and nematodynamics

The sL method is an explicit space-time integrator for advection equations
similar to an upwinding scheme [66]. Here we introduce the sL method
generalized to any Lie advection equation for tensor fields as follows. Given
a velocity field, we construct a backward flowmap by tracing the velocity
field upwind. This flowmap is used to pullback the tensor field from the
previous time step as the update rule. The principle can be applied to tensors
of arbitrary $(r,s)$ type (e.g., the valence-$k$ tensor that represents the
$k$-atic field [16]), serving as a discrete counterpart to the continuous
theory discussed in Sec. 2.1.3. The resultant scheme, summerized in Alg. 2,
captures the geometric structure of nematic advection and enhance numerical
stability.

Concretely, given a time-independent velocity field
$\boldsymbol{u}\in\Gamma_{\mathfrak{V}}(TM)$, we construct a finite-time
backward flowmap $\Psi_{\Delta t}$ by integrating the ODE,
$\partial_{\tau}\Psi_{\tau}=-\boldsymbol{u}(\Psi_{\tau})$, $\tau\in[0,\Delta
t]$, with the initial condition $\Psi_{0}=\mathrm{id}_{M}$. In practice, we
represent $\Psi_{\tau}$ and $\boldsymbol{u}$ as $\mathbb{R}^{3}$-valued
functions using surface embeddings and implement the integration using RK4. To
evaluate $\boldsymbol{u}\circ\Psi_{\tau}$, we use a closest-point projection
followed by a linear interpolation in a triangle. Recall in the continuous
setting that Lie advection
$\partial_{t}{[\boldsymbol{q}]}+\mathcal{L}_{\boldsymbol{u}}[\boldsymbol{q}]=0$
is characterized by a constant pullback field
$\partial_{t}\varphi^{*}[\boldsymbol{q}]=0$ (cf. Eqn. (9)). With the backward
flowmap $\Psi_{\Delta t}$ and its gradient $d\Psi_{\Delta t}$, a discrete Lie
advection is analogously represented by updating the nematic field
$[\boldsymbol{q}]$ using its pullback value $\Psi_{\Delta
t}^{*}[\boldsymbol{q}]$. The procedure will be made explicit in following
subsections.

####  3.2.1 Deformation gradient of backward flowmap

By evaluating the finite-time flowmap over vertices (i.e., $\Psi_{\Delta t}$
such that $M_{\Delta t}=\Psi_{\Delta t}M$), we can construct the deformation
gradient $d\Psi_{\Delta t}|_{f}$ (realized in $\mathbb{R}^{3}$) on each
triangle face $f$ as follows:

$\displaystyle d\Psi_{\Delta
t}:\Gamma_{\mathfrak{F}}(TM)\rightarrow\Gamma_{\mathfrak{F}}(T_{\Psi_{\Delta
t}}M_{\Delta t}),\quad d\Psi_{\Delta
t}|_{f}\coloneqq\begin{bmatrix}\underset{|}{\overset{|}{\boldsymbol{l}_{1}}}&\underset{|}{\overset{|}{\boldsymbol{l}_{2}}}&\underset{|}{\overset{|}{\boldsymbol{n}}}\end{bmatrix}_{\Psi_{\Delta
t}(f)}\begin{bmatrix}\underset{|}{\overset{|}{\boldsymbol{l}_{1}}}&\underset{|}{\overset{|}{\boldsymbol{l}_{2}}}&\underset{|}{\overset{|}{\boldsymbol{n}}}\end{bmatrix}_{f}^{-1}.$
(37)

These column vectors, illustrated by Fig. 3, represent the embedded edge
vectors $\boldsymbol{l}\in\mathbb{R}^{3}$ and the unit face normal
$\boldsymbol{n}\in\mathbb{R}^{3}$ before and after the backward flowmap
$\Psi_{\Delta t}$. A vertex-based deformation gradient
$\boldsymbol{F}:\Gamma_{\mathfrak{V}}(TM)\rightarrow\Gamma_{\mathfrak{V}}(T_{\Psi_{\Delta
t}}M_{\Delta t})$ is obtained by averaging deformation gradient on incident
faces, $\boldsymbol{F}|_{v}=\frac{1}{|f\succ v|}\sum_{f\succ v}d\Psi_{\Delta
t}|_{f}$.

Figure 3: The sL Lie advection is performed on a triangular mesh using a
backward finite time flowmap $\Psi$. This flowmap engenders a piecewise
constant deformation gradient $d\Psi|_{f\succ v}$ defined over the triangle
faces. By averaging this face-based deformation gradient around the one-ring
neighborhood of each vertex, we derive a vertex-based deformation gradient
$d\Psi|_{v}$. The nematic Lie advection is modeled by updating
$\boldsymbol{Q}_{v}$ with
$d\Psi|_{v}^{-1}\boldsymbol{Q}_{\Psi(v)}d\Psi|_{v}^{-\top}$.

####  3.2.2 Pullback of nematic field

Recall that the nematic field can be equivalently expressed using the matrix
representation
$\boldsymbol{Q}=\mathcal{V}([\boldsymbol{q}])\in\Gamma_{\mathfrak{V}}(TM\otimes
TM)$. With $\boldsymbol{Q}$ realized as an $\mathbb{R}^{3\times 3}$-valued
function, its Lie advection,
$\partial_{t}{\boldsymbol{Q}}+\mathcal{L}_{\boldsymbol{u}}\boldsymbol{Q}=0$,
with $t\in[0,\Delta t]$, $\boldsymbol{Q}|_{t=0}=\boldsymbol{Q}_{0}$, can be
modeled by

$\displaystyle\boldsymbol{Q}_{\Delta t}=\Psi_{\Delta
t}^{*}\boldsymbol{Q}_{0}=\boldsymbol{F}^{-1}(\boldsymbol{Q}_{0}\circ\Psi_{\Delta
t})\boldsymbol{F}^{-\top},$ (38)

where $\boldsymbol{Q}_{\Delta t}$ denotes the updated value of
$\boldsymbol{Q}$ for the subsequent time step. Because the backward flowmap
$\Psi_{\Delta t}$ can land anywhere on $M$, evaluating
$\boldsymbol{Q}_{0}\circ\Psi_{\Delta t}$ requires interpolation, such as
trilinear interpolation utilizing the face-barycentric representation of
$\Psi_{\Delta t}$. However, this interpolation may not preserve the rank-1
structure of $\boldsymbol{Q}$. To address this, we perform SVD and and extract
the leading component of the interpolated value, ensuring that
$\boldsymbol{Q}_{0}\circ\Psi_{\Delta t}=\boldsymbol{q}_{\Psi_{\Delta
t}}\otimes\boldsymbol{q}_{\Psi_{\Delta t}}$, and thereby
$\boldsymbol{Q}_{\Delta t}$ remains rank-1. Lastly, by projecting
$\boldsymbol{Q}_{\Delta t}$ onto the surface using tangent bases, we recover
the intrinsic representation of $\boldsymbol{Q}_{\Delta t}$ and its
corresponding complex representation $[\boldsymbol{q}]_{\Delta
t}=\mathcal{V}^{-1}(\boldsymbol{Q}_{\Delta t})$.

To model the Jeffery equation
$\dot{\boldsymbol{Q}}+\mathcal{P}^{E}(\mathcal{L}^{\lambda}_{\boldsymbol{u}}\boldsymbol{Q})=0$
(cf. Eqn. (31)), we can include the additional rotational Lie advection
$\mathcal{L}_{\boldsymbol{u}}^{\operatorname{rot}}$ and normalization as
follows:

$\displaystyle\boldsymbol{Q}_{\Delta
t}=\frac{\lambda\boldsymbol{F}^{-1}(\boldsymbol{Q}_{0}\circ\Psi_{\Delta
t})\boldsymbol{F}^{-\top}+(1-\lambda)\boldsymbol{R}^{-1}(\boldsymbol{Q}_{0}\circ\Psi_{\Delta
t})\boldsymbol{R}^{-\top}}{|\lambda\boldsymbol{F}^{-1}(\boldsymbol{Q}_{0}\circ\Psi_{\Delta
t})\boldsymbol{F}^{-\top}+(1-\lambda)\boldsymbol{R}^{-1}(\boldsymbol{Q}_{0}\circ\Psi_{\Delta
t})\boldsymbol{R}^{-\top}|},$ (39)

where $\boldsymbol{R}$ isolates rotation in $\boldsymbol{F}$ through polar
decomposition $\boldsymbol{F}=\boldsymbol{R}\boldsymbol{U}$.

As a side remark, in the special case where the transported field $h$ is a
scalar function, its advection
$\partial_{t}h+\mathcal{L}_{\boldsymbol{u}}h=\partial_{t}h+\nabla_{\boldsymbol{u}}h=0$
simplies to the traditional sL method, $h_{\Delta t}=\Psi_{\Delta
t}^{*}h_{0}=h_{0}\circ\Psi_{\Delta t}$.

Algorithm 2 sL method for nematic Lie advection over time interval $[0,\Delta
t]$

1:Inputs:

2: Surface mesh $M$, nematic field
$\boldsymbol{Q}_{0}=\mathcal{V}([\boldsymbol{q}]_{0})$, time step $\Delta t$,
flow velocity $\boldsymbol{u}$

3:for each vertex $v\in\mathfrak{V}$, do

4: $\Psi_{\Delta
t}\leftarrow\textsc{RK4}({d\Psi_{\tau}}/{d\tau}=-\boldsymbol{u}(\Psi_{\tau}))$
$\triangleright$ Backward flowmap

5: $\boldsymbol{Q}_{\Psi_{\Delta
t}}\leftarrow\textsc{Interp}(\boldsymbol{Q}_{0}\circ\Psi_{\Delta t})$

6: $\boldsymbol{Q}_{\Psi_{\Delta
t}}\leftarrow\textsc{SVD}(\boldsymbol{Q}_{\Psi_{\Delta t}})$

7:end for

8:for each face $f\in\mathfrak{F}$ do $\triangleright$ Deformation gradient

9: $d\Psi_{\Delta
t}\leftarrow[\boldsymbol{l}_{1},\boldsymbol{l}_{2},\boldsymbol{n}]_{\Psi_{\Delta
t}}[\boldsymbol{l}_{1},\boldsymbol{l}_{2},\boldsymbol{n}]^{-1}$

10:end for

11:for each vertex $v\in\mathfrak{V}$, do$\triangleright$ Pullback

12: $\boldsymbol{F}\leftarrow\textsc{Avg}(d\Psi_{\Delta t}|_{f\succ v})$

13:
$\boldsymbol{Q}_{v}\leftarrow\boldsymbol{F}^{-1}\boldsymbol{Q}_{\Psi_{\Delta
t}}\boldsymbol{F}^{-\top}$

14:end for

####  3.2.3 Discussion of sL advection scheme

Although the traditional scalar sL method guarantees unconditional stability,
the extension to tensor Lie advection in general does not [24]. The stability
depends on the conditioning of the pullback map. However, by incorporating a
diffusion term into the sL advection through time-splitting (cf. Alg. 1) and
scaling the system with respect to advection, the time step for advection
remains fixed regardless of the Péclet number. Because the diffusion equation,
$\dot{[\boldsymbol{q}]}=\mbox{\em Pe}^{-1}\Delta^{L}[\boldsymbol{q}]$, is
solved using an implicit method (cf. Alg. 1), the algorithm is numerically
stable with high Péclet number. This also highlights some well-known issues
associated with time-spliting and the sL method. First, time-spliting
introduces errors, which can be reduced using techniques such as the Strang
method [66]. Second, sL method based on linear interpolation has significant
numerical diffusion. Improvements can be made by implementing higher-order
interpolation and techniques such as the McCormack scheme [66]. However, these
challenges are common to the sL method and are not unique to the
generalization to tensor-valued functions. Despite these challenges, the
current minimal implementation is sufficient for our purpose of studying
active nematics on curved surfaces.

###  3.3 Gradient and (covaraint) divergence

We used the linear finite element hat function $\Phi_{v}$ to discretize the
gradient for a scalar function $p\in C_{\mathfrak{V}}(M)$ as follows:

$\displaystyle\mathrm{grad}:C_{\mathfrak{V}}(M)\rightarrow\Gamma_{\mathfrak{F}}(TM),\quad
A_{f}(\mathrm{grad}~{}p)_{f}=A_{f}\sum_{v\prec
f}p_{v}(\mathrm{grad}_{f}\Phi_{v})=\frac{1}{2}\sum_{v\prec
f}p_{v}\boldsymbol{\mathrm{i}}\boldsymbol{l}_{v\prec f},$ (40)

where $\boldsymbol{\mathrm{i}}\boldsymbol{l}_{v\prec f}$ denotes the
$90^{\circ}$ rotation of edge vector $\boldsymbol{l}_{v\prec f}$ and $A_{f}$
is the face area (cf. Sec. 3.1.1).

The divergence operator is given by the negated adjoint operation of the
gradient. The adjoint of Eqn. (40) applied to a face-based vector field
$\boldsymbol{u}\in\Gamma_{\mathfrak{F}}(TM)$ is a measure on each vertex
$(\mathrm{grad}^{*}\boldsymbol{u})_{v}=-{1\over 2}\sum_{f\succ
v}\langle\boldsymbol{\mathrm{i}}\boldsymbol{l}_{v\prec
f},\boldsymbol{u}_{f}\rangle$, representing the total divergence at vertex
$v$. Using a reference vertex area measure $A_{v}$ we define the divergence
operator as

$\displaystyle\mathrm{div}:\Gamma_{\mathfrak{F}}(TM)\rightarrow
C_{\mathfrak{V}}(M),\quad
A_{v}(\mathrm{div}~{}\boldsymbol{u})_{v}=-\frac{1}{2}\sum_{f\succ
v}\langle\boldsymbol{\mathrm{i}}\boldsymbol{l}_{v\prec
f},\boldsymbol{u}_{f}\rangle.$ (41)

Note that the discrete equation involving the divergence operator is often
written as in the weak form, where divergence is treated as a measure. In
particular these discrete equations are independent of the specific choice of
the vertex area measure $A_{v}$. In cases where pointwise evaluation is
needed, such as in visualization, we use $A_{v}=\sum_{f\succ v}A_{f}/3$.

The divergence operator, represented by Eqn. (41), can also be understood as a
finite-volume discretization applied to the dual cells of the mesh [66]. We
can analogously obtain a vertex-based divergence
$\mathrm{div}:\Gamma_{\mathfrak{V}}(TM)\rightarrow C_{\mathfrak{F}}(M)$ by
applying finite volume on the primal mesh.

By incorporating an additional vertex-face parallel transport,
$\Pi_{v\rightsquigarrow f}:T_{v}M\rightarrow T_{f}M$, (cf. Sec. 3.1.2), we can
extend such construction to discretize covariant divergence:

$\displaystyle\mathrm{div}^{\nabla}:\Gamma_{\mathfrak{V}}(TM\otimes
TM)\rightarrow\Gamma_{\mathfrak{F}}(TM),\quad
A_{f}(\mathrm{div}^{\nabla}~{}\boldsymbol{Q})_{f}=-\frac{1}{2}\sum_{v\prec
f}\Pi_{v\rightsquigarrow
f}\left[\boldsymbol{Q}_{v}\cdot\Pi^{-1}_{v\rightsquigarrow
f}(\boldsymbol{\mathrm{i}}\boldsymbol{l}_{v\prec f})\right].$ (42)

The Euclidean dot product, denoted as $\cdot$, denotes a contraction operation
on one slot of the $(2,0)$ tensor (i.e.,
$\boldsymbol{Q}\cdot\boldsymbol{q}\equiv Q^{ij}q_{j}$). In the case of
$\boldsymbol{Q}=\boldsymbol{q}\otimes\boldsymbol{q}$ being rank-1, the
covariant divergence can be explicitly expressed,
$A_{f}(\mathrm{div}^{\nabla}~{}\boldsymbol{Q})_{f}=-\frac{1}{2}\sum_{v\prec
f}\langle\Pi_{v\rightsquigarrow
f}(\boldsymbol{q}_{v}),\boldsymbol{\mathrm{i}}\boldsymbol{l}_{v\prec
f}\rangle\Pi_{v\rightsquigarrow f}(\boldsymbol{q}_{v})$. In practice, these
discrete operators, including those in Sec. 3.4, are organized and computed
using sparse matrix multiplication.

###  3.4 Bochner Laplacian

Here, we construct the discretization of Bochner Laplacian using the covariant
derivative and covariant divergence. We will first discuss the covariant
derivative and Bochner Laplacian for face-based vector fields and then proceed
with a similar construction for vertex-based nematic fields.

Using the face-face parallel transport $\Pi_{i\rightsquigarrow
j},i,j\in\mathfrak{F}$ across edge $e\prec i,j$ (cf. Eqn. (35)), discrete
directional covariant derivative $\nabla_{\boldsymbol{l}^{*}}$ maps a face-
based vector field
$\boldsymbol{u}=\hat{u}\boldsymbol{\mathrm{1}}\in\Gamma_{\mathfrak{F}}(TM)$ to
an edge-based vector field,

$\displaystyle\nabla_{\boldsymbol{l}^{*}}:\Gamma_{\mathfrak{F}}(TM)\rightarrow\Gamma_{\mathfrak{E}}(TM),\quad|\boldsymbol{l}_{e}^{*}|(\nabla_{\boldsymbol{l}^{*}}\boldsymbol{u})_{e}=(\boldsymbol{u}_{i}-\Pi_{i\rightsquigarrow
j}(\boldsymbol{u}_{j}))_{i,j\succ
e}=\left((\hat{u}_{j}-e^{-\boldsymbol{\mathrm{i}}\Omega_{e^{*}}}\hat{u}_{i})\boldsymbol{\mathrm{1}}_{j}\right)_{i,j\succ
e}.$ (43)

By taking the $L_{2}$ norm with the edge area as the diamond area formed by
$\boldsymbol{l}$ and $\boldsymbol{l}^{*}$,
$A_{e}=|\boldsymbol{l}||\boldsymbol{l}^{*}|$, we can construct the discrete
Dirichlet energy as:

$\displaystyle
2\mathcal{F}(\boldsymbol{u})=\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\nabla_{\boldsymbol{l}^{*}}\boldsymbol{u},\nabla_{\boldsymbol{l}^{*}}\boldsymbol{u}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}=\sum_{e\in\mathfrak{E}}A_{e}|\nabla_{\boldsymbol{l}^{*}}\boldsymbol{u}|_{e}^{2}=\sum_{e\in\mathfrak{E}}\frac{A_{e}}{|\boldsymbol{l}_{e}^{*}|^{2}}|\boldsymbol{u}_{i}-\Pi_{i\rightsquigarrow
j}(\boldsymbol{u}_{j})|_{i,j\succ
e}^{2}=\sum_{e\in\mathfrak{E}}\frac{1}{w_{e}}|\hat{u}_{j}-e^{-\boldsymbol{\mathrm{i}}\Omega_{e^{*}}}\hat{u}_{i}|_{i,j\succ
e}^{2}.$ (44)

As mentioned in Sec. 3.1.1, the edge length ratio
$w_{e}=|\boldsymbol{l}^{*}|/|\boldsymbol{l}|=(\cot\angle_{e\prec
i}+\cot\angle_{e\prec j})/2$ is commonly referred to as the cotangent weight,
where $\angle_{e\prec i}$ refers to the corner angle opposite to the edge $e$
at face $i$.

The discrete Bochner Laplacian can be obtained by taking the negated variation
of the discrete Dirichlet energy,

$\displaystyle\Delta:\Gamma_{\mathfrak{F}}(TM)\rightarrow\Gamma_{\mathfrak{F}}(TM),\quad
A_{f}(\Delta\boldsymbol{u})_{f}=-\left(\frac{\delta\mathcal{E}}{\delta\boldsymbol{u}}\right)_{f}=-\boldsymbol{\mathrm{1}}_{f}\sum_{e\prec
f}\frac{1}{w_{e}}(\hat{u}_{f}-e^{-\boldsymbol{\mathrm{i}}\Omega_{e^{*}}}\hat{u}_{i})_{i,f\succ
e}.$ (45)

Here, face $i$ is in the one-ring neighbor of face $f$, and $i$ and $f$ share
the edge $e$.

Similarly, using the vertex-vertex transportation $\Pi_{i\rightsquigarrow
j},i,j\in\mathfrak{V}$ across edge $e\succ i,j$ (cf. Eqn. (36)), we can
construct the covariant derivative for a vertex-based nematic field as
follows:

$\displaystyle\nabla^{L}_{\boldsymbol{l}}:\Gamma_{\mathfrak{V}}(TM)\rightarrow\Gamma_{\mathfrak{E}}(TM),\quad|\boldsymbol{l}_{e}|(\nabla_{\boldsymbol{l}}[\boldsymbol{q}])_{e}=([\boldsymbol{q}]_{i}-\Pi_{i\rightsquigarrow
j}([\boldsymbol{q}]_{j}))_{i,j\prec
e}=\left((\hat{q}_{j}-e^{-\boldsymbol{\mathrm{i}}2\Omega_{e}}\hat{q}_{i})[\boldsymbol{\mathrm{1}}]_{j}\right)_{i,j\prec
e}.$ (46)

The nematic Bochner Laplacian acting on the vertex-based nematic field
$[\boldsymbol{q}]=\hat{[q]}[\boldsymbol{\mathrm{1}}]\in\Gamma_{\mathfrak{V}}(L)$
is therefore:

$\displaystyle\Delta^{L}:\Gamma_{\mathfrak{V}}(L)\rightarrow\Gamma_{\mathfrak{V}}(L),\quad
A_{v}(\Delta[\boldsymbol{q}])_{v}=-[\boldsymbol{\mathrm{1}}]_{v}\sum_{e\prec
v}w_{e}(\hat{[q]}_{v}-e^{-\boldsymbol{\mathrm{i}}2\Omega_{e}}\hat{[q]}_{i})_{i,f\prec
e},$ (47)

where vertex $i$ is in the one-ring neighbor of vertex $v$, connected by edge
$e$.

Strictly speaking, there is an ambiguity in the definition of the directional
covariant derivative (Eqns. (43) and (46)) because it is equally valid to
transport a vector (nematic director) from face (vertex) $j$ to face (vertex)
$i$ before applying finite difference. However, the discretization of the
Dirichlet energy, and thus the Bochner Laplacian (Eqns. (45) and (47)), relies
only on the norm of the covariant derivative, which remains invariant,
$|\Pi_{j\rightsquigarrow
i}(\boldsymbol{u}_{j})-\boldsymbol{u}_{i}|=|\boldsymbol{u}_{j}-\Pi_{i\rightsquigarrow
j}(\boldsymbol{u}_{i})|$, $|\Pi_{j\rightsquigarrow
i}([\boldsymbol{q}]_{j})-[\boldsymbol{q}]_{i}|=|[\boldsymbol{q}]_{j}-\Pi_{i\rightsquigarrow
j}([\boldsymbol{q}]_{i})|$.

In summary, the intrinsic definition of the discrete Bochner Laplacian emerges
naturally from continuous theory of complex line bundles. Contrasting the
previous covariant divergence (Eqn. (42)) that relies on an extrinsic parallel
transport achieved through dihedral rotations between the tangent spaces of
vertices and faces, the Bochner Laplacian is constructed purely intrinsically
by vertex-vertex or face-face parallel transport. Compared to common
approaches that rely on the embedding structure in $\mathbb{R}^{3}$ and a
projection operator from $\mathbb{R}^{3}$ to the manifold, the discrete
complex line bundle approach has notably improved efficiency. This approach
reduces the degree of freedom per element from 3 to 2, and the complex nature
furthermore takes advantage of the optimized complex arithmetic handling in
established numerical linear algebra libraries. The corresponding complex-
valued Laplace matrix is Hermitian and negative-semidefinite.

###  3.5 Curvature, defects, and local-global geometry

Here we introduce a discretization of the face Gaussian curvature and an
algorithm to detect topological defects and compute their charges. This serves
as the discrete counterpart to the continuous theory in Sec. 2.1.8. The
construction of curvature and charge leads to discrete analogs of local-to-
global theorems including the Poincaré-Hopf and Gauss-Bonnet theorems.

Given a vertex-based $k$-atic field
$[\boldsymbol{q}]_{k}:\mathfrak{V}\rightarrow V_{k}$, we can utilize vertex-
vertex parallel transport $\Pi_{i\rightarrow j}$ (cf. Eqn. (36)) with
$i,j\in\mathfrak{V}$ along the incident edge $e\succ i,j$ to measure the
finite rotation $|\boldsymbol{l}_{e}|\eta_{e}$ (cf. Eqn. (19)) of
$[\boldsymbol{q}]_{k}$ along $e$ as follows:

$\displaystyle|\boldsymbol{l}_{e}|\eta_{e}=\arg\left(\frac{[\boldsymbol{q}]_{k,j}}{\Pi_{i\rightsquigarrow
j}[\boldsymbol{q}]_{k,i}}\right)=\arg\left(\frac{\hat{[\boldsymbol{q}]}_{k,j}}{e^{-\boldsymbol{\mathrm{i}}k\Omega_{e}}\hat{[\boldsymbol{q}]}_{k,i}}\right).$
(48)

Assuming the $k$-atic field $[\boldsymbol{q}]_{k}$ is continuous between
vertices $i$ and $j$, we choose the principal branch of the $\arg$ function so
that $|\boldsymbol{l}_{e}|\eta_{e}$ lies within $(-\pi,\pi]$.

Recall that in the continuous setting, Gaussian curvature quantifies holonomy,
$\int_{U}KdA=-\int_{\partial U}\omega$ (cf. Sec. 2.1.8). Analogously, discrete
Gaussian curvature $K\in C_{\mathfrak{F}}(M)$ can be represented by closed-
loop summation of the primal edge connection $\Omega_{e}$ as

$\displaystyle A_{f}K_{f}\coloneqq-\sum_{e\prec f}\Omega_{e}\mod
2\pi\in(-\pi,\pi].$ (49)

Selecting branch $(-\pi,\pi]$ for $A_{f}K_{f}$ is justified by the small
integration domain of face $f$. The discrete curvature is also the angle
defect at face $f$, given by $A_{f}K_{f}=\sum_{v\prec
f}\tilde{\angle}_{v}-\pi$, where $\tilde{\angle}_{v}$ is the geodesic corner
angle (cf. Sec. 3.1). The discrete Gauss-Bonnet theorem states,
$\sum_{f}A_{f}K_{f}=2\pi\chi(M)$, where
$\chi(M)=|\mathfrak{V}|-|\mathfrak{E}|+|\mathfrak{F}|=2-2g$ represents the
Euler characteristic of the polygonal mesh $M$.

Analogous to Eqn. (20), given a $k$-atic field $[\boldsymbol{q}]_{k}$, its
topological charge is a fraction-valued function
$\mathcal{Z}([\boldsymbol{q}]_{k}):\mathfrak{F}\rightarrow\mathbb{Z}/k$ that
can be computed for each face $f$ as follows:

$\displaystyle\mathcal{Z}_{f}([\boldsymbol{q}]_{k})$
$\displaystyle=\frac{1}{2\pi k}\sum_{e\prec
f}|\boldsymbol{l}_{e}|\eta_{e}+\frac{1}{2\pi}A_{f}K_{f}.$ (50)

The charge $\mathcal{Z}_{f}$ at face $f$ is zero when the adjacent $k-$atic
directors $[\boldsymbol{q}]_{k,v\prec f}$ are continuous, and nonzero in the
presence of defects. The charge function satisfies the discrete Poincaré-Hopf
theorem exactly, $\sum_{f}\mathcal{Z}_{f}([\boldsymbol{q}]_{k})=\chi(M)$.

###  3.6 Augmented Lagrangian method for Stokes equations

We utilize a classic AL iteration to solve the incompressible Stokes flow
[67]. Recall that the Stokes equations can be viewed as the stationary
condition for a divergence-constrained variational problem (cf. Sec. 2.2.1).
The AL approach augmented the original Rayleighian dissipation functional with
a quadratic penalty term on divergence. The augmented Rayleighian is given by:

$\displaystyle\begin{split}\min_{\boldsymbol{u}}\max_{p}~{}\mathcal{R}_{\mathrm{AL}}&=\mathcal{R}(\boldsymbol{u},p)+\frac{k}{2}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\text{div
}\boldsymbol{u},\text{div
}\boldsymbol{u}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}},\end{split}$
(51)

where the stationary condition is
$(\Delta+K+k~{}\mathrm{grad}\circ\mathrm{div})\boldsymbol{u}-\text{grad
}p+\boldsymbol{f}=0$. Here, the penalty coefficient $k$ enforces the
incompressibility constraint during the iterative solution process. Adjusting
the penalty parameter allows effective control over the enforcement of the
constraint. It is worth noting that physically this penalty term is also
included in the viscosity term, which provides resistance to expansion and
compression using the same mechanism through bulk viscosity. Shown in Alg. 3,
the solve process involves iteratively solving the modified Stokes equations.
At each iteration, the velocity and pressure fields are updated until a
convergence criterion is met. Practically, increasing the value of $k$ makes
the fluid very stiff in response to compressibility, which leads to
convergence in a small number of iterations. However, a high value of $k$ also
increases the spectral condition number of the matrix
$A=-(\Delta+K+k~{}\text{grad}\circ\text{div})$, which increases the time per
iteration [67]. An alternative method is to use the Uzawa method combined with
pressure projection [67]. We implement the Poisson solve using the Python
linear algebra library scipy.sparse, which utilizes a sparse matrix
representation and solves the system using the conjugate gradient method.
Compared to a direct solver, the conjugate gradient is suitable for medium to
large system sizes, as it effectively handles the sparse and symmetric
positive-semidefinite nature of the matrices involved.

Algorithm 3 AL subroutine for incompressible Stokes equations

1:Inputs:

2: Penalty coefficient $k$, tolerance $\varepsilon$

3:Initialize:

4: Initial guess $p$

5:while $\|\mathrm{div}~{}\boldsymbol{u}\|>\varepsilon$ do

6:
$\boldsymbol{u}\leftarrow\textsc{Solve}(-(\Delta+K+k~{}\mathrm{grad}\circ\mathrm{div})\boldsymbol{u}=\boldsymbol{f}-\text{grad
}p)$ $\triangleright$ Poisson solve (A x = b)

7: $p\leftarrow p+k~{}\mathrm{div}~{}\boldsymbol{u}$

8:end while

9:

10:if exist $N$ distinct $\boldsymbol{\omega}_{i}$ then$\triangleright$ Remove
the Killing component

11: for $i=1$ to $N$ do

12:
$\boldsymbol{u}\leftarrow\boldsymbol{u}-\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\boldsymbol{u},\boldsymbol{\omega}_{i}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}\boldsymbol{\omega}_{i}$

13: end for

14:end if

As discussed in Sec. 2.2.1, when dealing with the Stokes equations in
symmetric geometries that possess a Killing vector field, the solution is not
unique. Despite this, iterative methods such as the conjugate gradient method
can still converge to one feasible solution. To ensure that we obtain
continuous solution during time evolution, we consistently select the solution
with the minimal $L_{2}$ norm by projecting out components that lie in the
space of the Killing vector fields. To achieve this, we precompute the set of
Killing bases of the geometry, denoted as $\boldsymbol{\omega}_{i}$, before we
begin the evolution of the active nematic system. The bases are determined by
solving the nonvanishing vector field that satisfies
$(\Delta+K+\mathrm{grad}\circ\mathrm{div})\boldsymbol{\omega}_{i}=0$. In
practice, we solve this equation computationally through the smallest
eigenvalue problem and extract eigenfunctions with near-zero eigenvalues. Here
each Killing basis is normalized such that
$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\boldsymbol{\omega}_{i},\boldsymbol{\omega}_{j}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}=\delta_{ij}$.
As outlined in Alg. 3, we project out components spanned by these bases after
the AL iteration. This protocol is carried out in the Sec. 4 for examples on
analytical shapes. Note that the number of distinct Killing vector fields of a
given geometry is often apparent from its symmetry prior to performing the
actual computation. For example, a sphere exhibits three distinct Killing
vector fields, and a torus displays one, corresponding respectively to their
modes of rotational symmetries.

###  3.7 Implementation

Along with the theoretical and algorithmic development, we also provide an
implementation using SideFx Houdini, a freely-accessible software under its
educational license. Our implementation is accessible through a public
repository on Github: https://github.com/CunchengZhu/Riemannian-active-
nematics-2024.git. SideFx Houdini is primarily known as a 3D animation and
special effects software. However, it can be effectively utilized for custom
algorithm implementation by leveraging its half-edge data structure of
triangular meshes, as well as its VEX and Python programming capabilities. VEX
is a C-like shader language that operates in a vectorized, multithreaded
manner. The Python module enables access to the numerical linear algebra
libraries Scipy. Additionally, the software provides visualization and basic
geometric processing features such as mesh parametrization and interpolation
to facilitate rapid prototyping.

##  4 Results and discussion

In the preceding sections, we established the continuous formulation, its
corresponding discretization, and the numerical methods employed. In this
section, we demonstrate the robustness of the method and its potential in
addressing various biophysically relevant scenarios. First, we conduct a
convergence analysis spanning the various subroutines of our methodology. This
analysis not only attests to the method’s accuracy in approximating the smooth
solution but also illustrates the performance of the algorithms. In addition,
we explore the efficiency and scalability of the methods by presenting the
computational-time scaling characteristics of the major components within the
framework. We have chosen to benchmark these methods on simple geometries,
where a mix of numerical experiments and analytical comparisons are possible.
However, the application extends beyond idealized geometries. To demonstrate
the generality of these methods on arbitrary geometries and topologies, we
apply them to solve the full nematodynamics equations in biophysically
relevant examples. In this section, through the error analysis, scalability