# Free energy formulas for confined nematic liquid crystals based on analogies
with Kirchhoff-Routh theory in vortex dynamics

Hiroyuki Miyoshi1, Hiroki Miyazako1,2, Takaaki Nara1,2

1Graduate School of Information Science and Technology,
The University of Tokyo, Hongo, 7-3-1, Bunkyo-Ku, Tokyo, 113-8656, Japan
2Department of Information Physics and Computing,
The University of Tokyo, Hongo, 7-3-1, Bunkyo-Ku, Tokyo, 113-8656, Japan
Corresponding author<EMAIL_ADDRESS>
###### Abstract

Active nematics are influenced by alignment angle singularities called
topological defects. The localization of these defects is of major interest
for biological applications. The total distortion of alignment angles due to
defects is evaluated using Frank free energy, which is one of the criteria
used to determine the location and stability of these defects. Previous work
used the line integrals of a complex potential associated with the alignments
for the energy calculation (Miyazako and Nara, R. Soc. Open Sci., 2022), which
has a high computational cost. We propose analytical formulas for the free
energy in the presence of multiple topological defects in confined geometries.
The formulas derived here are an analogue of Kirchhoff-Routh functions in
vortex dynamics. The proposed formulas are explicit with respect to the defect
locations and conformal maps, which enables the explicit calculation of the
energy extrema. The formulas are applied to calculate the locations of defects
in so-called doublets and triplets by solving simple polynomial formulas. A
stability analysis is also conducted to detect whether defect pairs with
charges $\pm 1/2$ are stable or unstable in triplet regions. Our numerical
results are shown to match the experimental results (Ienaga et al., Soft
Matter, 2023).

_Keywords_ topological defects $\cdot$ conformal mappings $\cdot$ nematic
liquid crystals

##  1 Introduction

A spindle-shaped cell is characterized by its elongated shape. Such cells play
crucial roles in various physiological processes, including muscle
contraction, tissue organization, and cell migration. Investigations of the
dynamics of biological cells are important for elucidating fundamental
processes in development, disease progression, and tissue homeostasis. In such
investigations, spindle- and rod-shaped cells are modeled as nematic liquid
crystals, which are also known as “active nematics” [1, 2]. Active nematics
are essential not only for understanding the collective dynamics of biological
systems such as cells and tissues [3] but also for artificial biological
systems. The various phenomena observed in active nematics under diverse
scenarios have been reviewed in [4].

An important property in the field of active nematics is the existence of
topological defects [5, 2]. Topological defects are regions where the
alignment angles of nematic liquid crystals are not defined. Such
singularities influence the orientation of surrounding nematics [6] and
determine the physical behaviours of active nematics. The information of the
rotations in the surrounding alignments is referred to as the topological
charge of defects, which takes discrete values of integers or half-integers.

Topological defects play a significant role in biological systems. Maroudas-
Sacks et al. presented experimental data on topological defects in Hydra and
showed that topological defects with topological charges $-1/2$, $+1/2$, and
$+1$ regulate the process of morphogenesis [7]. Kawaguchi et al. found that
topological defects drive the aggregation of cultured murine neural progenitor
cells [8]. They demonstrated that cells concentrate around defects with a
charge of $+1/2$, and escape from defects with a charge of $-1/2$.

The behaviour of topological defects is influenced by the presence of
boundaries. Saw et al. observed topological defects with charges of $+1/2$ and
$-1/2$ in a star-shaped region to study cell death and extrusion [9]. Ienaga
et al. modified the closed domains of single circular shapes to regions
consisting of multiple discs overlapping each other (doublets and triplets) to
study the relationship between existing defects and geometries [10]. They
conjectured that the positions of defects in both doublets and triplets
linearly correlate with the distances between the centres of circles.
Furthermore, they experimentally examined the transition from a formation
where two defects with a charge of $+1/2$ exist to a state where a defect with
a charge of $-1/2$ and three defects with a charge of $+1/2$ exist as the
three discs of triplets move apart.

Based on liquid crystal theory and one-constant approximation [5],
mathematical analysis using the theory of two-dimensional (2D) harmonic
functions has been conducted to determine the location and stability of
topological defects in confined geometries. Duclos et al. investigated the
positions of topological defects in circular closed domains and confirmed both
experimentally and mathematically that two defects with a charge of $+1/2$ are
located at a distance of $5^{-1/4}R$ from the centre of the disc with radius
$R$ [11]. They modelled the Frank free energy associated with cell alignment
angles based on one-constant approximation and mathematically proved that
defects exist at the minimum of this free energy. They expressed the energy as
the surface integral of the square of the spatial gradients of orientation
excluding the vanishingly small defect core regions.

Because the alignment angles are harmonic functions in 2D, it is natural to
use a complex analysis formulation. Tarnavskyy et al. showed that topological
defects exist on the edge of a half of a lens-shaped region using conformal
mappings [12]. Chandler and Spagnolie used Schwarz-Christoffel mappings to
find topological defects on the boundary of a polygonal-shaped immersed body
[13]. The theory developed by Chandler and Spagnolie was also applied to study
two-body interactions such as body forces and torques [14]. Tang et al. showed
that Frank free energy can be classified into torque-free energy and
orientation energy and that torque determines the dynamics of topological
defects [15].

Miyazako and Nara proposed a complex potential for cell alignment angles based
on complex analysis to represent the orientational order of cells [16]. They
derived an analytical formula that describes the orientation of cells with
topological defects at arbitrary positions strictly inside arbitrary closed
domains [16]. They employed the theory of conformal mapping to compute cell
orientations in arbitrary shapes. Additionally, by utilizing the complex
Green’s theorem, they evaluated the energy of the domain as line integrals to
avoid the divergence of energy at singularities of defects. They updated
defect positions using steepest descent methods, and discussed whether the
initial configuration of defects would annihilate or stabilize in star-shaped
domains. Recently, this approach was extended to doubly connected domains [17]
using the prime functions developed by Crowdy [18].

Following Duclos et al. [11], we assume that topological defects become stable
when the Frank free energy is minimized. Energy calculations are necessary for
localizing topological defects at arbitrary positions in confined domains.
However, there are no analytical formulas for the Frank free energy in these
geometries except for the line integrals proposed by Miyazako and Nara [16],
which have a high numerical cost and require finding a suitable branch cut of
the logarithmic singularities associated with defects. This significantly
complicates the theoretical analysis of the positions of topological defects.

In this paper, we develop analytical formulas for the Frank free energy when
the radius of the defect core is small as shown in Figure 1. The explicit
formulas are derived from an expression for the complex potential,
representing the sum of each defect, excluding the contribution from the
defect itself. The proposed formulas are new in the field of liquid crystal
theory. This paper provides an analogy between the concept of Hamiltonian or
Kirchhoff-Routh functions for vortex dynamics [19, 18] and the Frank free
energy of nematic liquid crystals. The proposed formulas enable to derive the
explicit expression for the energy up to the order $\epsilon\log\epsilon$,
where $\epsilon$ is the size of defect cores. It is mathematically
demonstrated that two defects are located symmetrically at $\pm 5^{-1/4}$ in
the unit disc. We also validate the conjecture regarding the positions of
topological defects in doublets and triplets given by Ienaga et al. [10]. It
is found that our analytical results closely match the experimental results
for doublets reported by Ienaga et al., and detect the transition of the
formation of topological defects in triplets.

The rest of this paper is organized as follows. Section 2 introduces a
mathematical formulation of alignment with topological defects and the energy
of alignment angles. Section 3 presents the proposed formulas for the Frank
free energy in the unit disc and arbitrarily shaped regions. Section 4
calculates the locations of topological defects in doublets and triplets based
on the theory of quadrature domains [20] and the stability of these defects.
Section 5 gives the conclusions and provides future prospects.

Figure 1: Schematic diagram of topological defects (red circles) with charge
of $+1/2$ and resulting alignments due to defects and boundary. The alignment
angle changes by $\pi$ around the topological defects. The defects here are
assumed to have a defect core with size $\epsilon$, which is sufficiently
small compared to the outer boundary of confined domain $\Omega$.

##  2 Formulation of topological defects and Frank-free energy

###  2.1 Problem setting

Following the theory developed by Miyazako and Nara [11], this section
introduces a mathematical formulation of alignment angles in the presence of
multiple topological defects based on complex analysis. The nematic liquid
crystals are assumed to be spindle-shaped and small relative to the outer
boundary, thus allowing them to be treated as rod-like particles. They are
assumed to be aligned in the direction of a unit vector $\bm{n}_{\rm d}$ as
follows:

$\displaystyle\bm{n}_{\rm d}=(\cos\phi(u,v),\sin\phi(u,v))^{\top},$ (1)

where $\phi(u,v)$ is the angle of the field at position $(u,v)$. The direction
$\bm{n}_{\rm d}$ is called a “director”. Since $\bm{n}_{\rm d}$ and
$-\bm{n}_{\rm d}$ are indistinguishable, an arbitrary angle $\phi$ is
equivalent to $\hat{\phi}\in[-\pi/2,\pi/2)$ [16]; that is, the two angles
$\phi_{1}$ and $\phi_{2}$ are equivalent if there exists an integer
$n\in\mathbb{Z}$ such that $\phi_{1}=\phi_{2}+n\pi$.

The governing equation for the alignment angle $\phi$ is derived using the
one-constant approximation of elastic energy called Frank free energy and its
variational formula. The Frank free energy with $\phi(u,v)$ at position
$(u,v)$ is defined as the integral of a surface $\Omega$ as follows:

$\displaystyle F\equiv\frac{K}{2}\int_{\Omega}|\nabla\phi(u,v)|^{2}{\rm
d}u{\rm d}v,$ (2)

where $K$ is a Frank elastic constant and $K/2=1$ in this paper for
simplicity. The variation of the energy defined in equation (2) leads to the
energy taking a minimum when the alignment angle $\phi(u,v)$ satisfies the
following 2D Laplace’s equation:

$\displaystyle\nabla^{2}\phi(u,v)=0,\quad(u,v)\in\Omega.$ (3)

Figure 2: (a) Colour plot of ${\rm Re}[f_{\mathbb{D}}(z|\mathcal{Z})]$ with
three defects with charge of $+1/2$, and one defect with charge of $-1/2$.
Four branch cuts exist due to the logarithmic singularities in the unit disc.
(b) Angle $\phi$ on the unit circle decomposed into $\theta$ and step function
$\chi(\theta)$ associated with branch cuts minus $\pi/2$. Note that
$\theta+\chi(\theta)$ on $\partial\mathbb{D}$ is discontinous at (i)-(iv).

A complex analysis formulation of the alignment for simply connected domains
is now introduced. First, an analytical formula for the unit disc is
introduced. Let $z=x+{\rm i}y$ define a complex plane and $\mathbb{D}$ define
the unit disc. There exist topological defects at $z_{k}=x_{k}+{\rm i}y_{k}$,
$k=1,\ldots,N$, strictly inside $\mathbb{D}$, whose topological charges are
$q_{k}$. The total number of topological charges is $N$ and the charges are
integers or half-integers except $0$; that is:

$\displaystyle q_{k}\in\\{n/2|\ n\in\mathbb{Z},\ n\neq 0\\}.$ (4)

We denote the set of positions of defects as
$\mathcal{Z}\equiv\\{z_{k}|k=1,2,\ldots,N\\}$.

The alignment angle $\phi(x,y)$ created by these topological defects satisfies
the following conditions [11, 16]:

  1. 1.

Alignment angle $\phi$ satisfies 2D Laplace’s equation except at
$\mathcal{Z}$:

$\displaystyle\nabla^{2}\phi(x,y)=0,\quad(x,y)\in\mathbb{D}\backslash\mathcal{Z}.$
(5)

  2. 2.

The sum of topological charges $q_{k}$ is equal to $1$ [16]; that is:

$\displaystyle\sum_{k=1}^{N}q_{k}=1.$ (6)

  3. 3.

For simply connected domains $D_{k}\subset\mathbb{D}$ $(k=1,\ldots,N)$ that
contain only the $k$-th defect, the circulation around $D_{k}$ satisfies

$\displaystyle\oint_{\partial D_{k}}\frac{\partial\phi}{\partial s}{\rm
d}s=2\pi q_{k}.$ (7)

  4. 4.

Because the directors align along the boundary of $\mathbb{D}$, alignment
angle $\phi(x,y)$ satisfies

$\displaystyle\phi(x,y)=\nu(x,y),\quad(x,y)\in\partial\mathbb{D},$ (8)

where $\nu(x,y)$ is a tangential angle on $\partial\mathbb{D}$.

Miyazako and Nara derived an analytical formula that describes the alignment
angle that satisfies conditions (i)-(iv) in $\mathbb{D}$ [16]. The alignment
angle $\phi(x,y)$ for the unit disc is given by $\phi(x,y)={\rm
Re}[f_{\mathbb{D}}(z|\mathcal{Z})]$, where $f_{\mathbb{D}}(z|\mathcal{Z})$ is
an analytic function in domain $\mathbb{D}\backslash\mathcal{Z}$ which is
explicitly given by

$\displaystyle
f_{\mathbb{D}}(z|\mathcal{Z})=\sum_{k=1}^{N}q_{k}\mathcal{G}(z,z_{k})-\frac{\pi}{2},$
(9)

where a complex function for the cell alignment angles is defined by

$\displaystyle\mathcal{G}(z,z_{k})\equiv-{\rm
i}(\log(z-z_{k})+\log(1-\overline{z_{k}}z)),$ (10)

where $\overline{a}$ represents the complex conjugate of complex number $a$.
Topological charge $q_{k}$ satisfies condition (6). Figure 2 (a) shows an
example of $\phi$ in the presence of three topological defects with a charge
of $+1/2$ and a defect with a charge of $-1/2$.

Miyazako and Nara extended the analytical formula to simply connected
arbitrarily shaped domains using conformal mappings. Riemann’s mapping theorem
states that there exists a conformal mapping $w=g(z)$ from the unit disc
$\mathbb{D}$ in the complex $z$-plane to the region $\Omega$ in a complex
$w=u+{\rm i}v$-plane. The positions of topological defects inside $\Omega$ are
denoted as $w_{k}$, $k=1,\ldots,N$, and the set of these positions is
$\mathcal{W}\equiv\\{w_{1},w_{2},\ldots,w_{N}\\}$. The preimage of $w_{k}$ is
denoted as $z_{k}$; that is, $z_{k}=g^{-1}(w_{k})\in\mathbb{D}$,
$z_{k}\notin\partial\mathbb{D}$. Using the map $w=g(z)$, the tangential angle
on the boundary of $\Omega$ is given by

$\displaystyle\nu(\theta)={\rm Re}\left[-{\rm
i}\log\left(\frac{du}{d\theta}+{\rm i}\frac{dv}{d\theta}\right)\right]={\rm
Re}\left[-{\rm i}\log g^{\prime}(z)+\theta+\frac{\pi}{2}\right],\quad
z=e^{{\rm i}\theta},$ (11)

which gives the following complex potential of cell alignment for $\Omega$
[16, 17]:

$\displaystyle f_{\Omega}(z|\mathcal{Z})=f_{\mathbb{D}}(z|\mathcal{Z})-{\rm
i}\log g^{\prime}(z),\quad w=g(z).$ (12)

###  2.2 Line integral for energy calculations

The positions of topological defects $w_{k}$ tend to be located where the
energy defined by equation (2) satisfies the local minimum [11, 16]. Due to
the presence of topological defects at $w_{k}$, the surface integral of
$|\nabla\phi|^{2}$ is divergent and thus special care must be taken. Duclos et
al. [11] used the Fourier expansions and calculated the energy in the disc
with radius $R$ except for small regions around the defects. They showed that
the free energy of the disc is minimized when there exist two defects located
in the disc at a distance $\pm 5^{-1/4}R$ from the origin. Miyazako and Nara
modified the calculations while preventing the energy divergence of
$|\nabla\phi|^{2}$ at the defects and calculated the energy in star-shaped
domains [16]. With the use of the complex Green’s theorem, the energy (2) is
written as the following line integral:

$\displaystyle F=\int_{\Omega(\epsilon)}|\nabla\phi(u,v)|^{2}{\rm d}u{\rm
d}v=\frac{1}{2{\rm i}}\oint_{\partial\Omega(\epsilon)}\frac{\partial
f_{\Omega}}{\partial w}\overline{f_{\Omega}}dw,\quad\frac{\partial}{\partial
w}\equiv\frac{1}{2}\left(\frac{\partial}{\partial u}-{\rm
i}\frac{\partial}{\partial v}\right),$ (13)

where a small disc around the $k$-th defect with radius $\epsilon$ is defined
by

$\displaystyle D_{k}(\epsilon)\equiv\\{w_{k}+re^{{\rm i}\theta}|\ 0\leq
r\leq\epsilon,0\leq\theta<2\pi\\},$ (14)

and $\Omega(\epsilon)\equiv\Omega\backslash\cup_{k=1}^{N}D_{k}$. The path of
the line integral for the case $N=2$ is illustrated in Figure 3.

Figure 3: Line integral along $\partial(\mathbb{D}\backslash D(\epsilon))$ for
case $N=2$. Branch cuts of logarithmic singularities $z_{k}$ are chosen as
$l_{k}$ and the imaginary part of $\log(z-z_{k})$ jumps on $l_{k}$ by $2\pi$.

Using the line integral formula (13), Miyazako and Nara confirmed that for the
unit disc, the energy is minimized when defects exist at positions $z=\pm
5^{-1/4}$. They verified that when there are multiple $+1/2$ defects and
$-1/2$ defects within the unit disc, the pair of positive and negative defects
annihilate through energy minimization. They also showed that this
annihilation happens on symmetrical domains that contain multiple convex
regions such as triplets.

Although the line integral (13) is evaluated in a simple manner, the
logarithmic branch cuts of $f_{\Omega}$ exist between $\partial
D_{k}(\epsilon)$ and $\partial\Omega$. When there exist multiple defects,
special care must be taken when selecting the paths to prevent these branch
cuts from crossing.

##  3 Analytical formulas for the Frank free energy

This section proposes an analytical formula for the energy defined by (2) in
the unit disc $\mathbb{D}$ and in arbitrary domain $\Omega$. This can be done
by splitting equation (9) into a singular part that includes a logarithmic
pole at $z=z_{k}$ and a regular part around $z=z_{k}$. The concept of this
function is identical to that of the Kirchhoff-Routh functions in the field of
vortex dynamics [19, 21]. This enables the explicit evaluation of the line
integral (13), which gives an analytical formula for the energy up to the
order $\mathcal{O}(\epsilon\log\epsilon)$.

###  3.1 Analytical formula for unit disc

######  Theorem 3.1.

Assume that there exist $N$ topological defects with an integer or half-
integer charge located at $z_{k}\in\mathbb{D}$, $k=1,\ldots,N$,
$z_{k}\notin\partial\mathbb{D}$. The sum of the topological charges satisfies
(6). Then the total energy $F$ of the region $\mathbb{D}\backslash
D(\epsilon)$, $D(\epsilon)\equiv\cup_{k=1}^{N}D_{k}(\epsilon)$ is written as
follows:

$\displaystyle F=-2\pi\sum_{k=1}^{N}q_{k}^{2}\log\epsilon+2\pi
F_{d}+\mathcal{O}(\epsilon\log\epsilon),$ (15)

where the first term is independent of $z_{k}$ and $F_{d}$ only depends on the
position of the defects $z_{k}$, and is given by

$\displaystyle F_{d}=\sum_{k=1}^{N}q_{k}^{2}{\rm
Im}[\mathcal{R}(z_{k})]+\sum_{k=1}^{N}\sum_{l\neq k}^{N}q_{k}q_{l}{\rm
Im}[\mathcal{G}(z_{k},z_{l})],$ (16)

where $\mathcal{G}(.,.)$ is given in equation (10) and

$\displaystyle\mathcal{R}(z_{k})=\lim_{z\rightarrow
z_{k}}\left(\mathcal{G}(z,z_{k})+{\rm i}\log(z-z_{k})\right)=-{\rm
i}\log(1-z_{k}\overline{z_{k}})-\frac{\pi}{2}.$ (17)

The function $\mathcal{G}(z)+{\rm i}\log(z-z_{k})$ is a function analytic at
$z=z_{k}$, which is obtained by removing the contribution of the logarithmic
pole at $z=z_{k}$.

Remark 1: The first term in (15) is independent of locations $z_{k}$, which
means that $z_{k}$ only contributes the second term in (15). The minimum of
the energy is obtained by minimizing $F_{d}$.

Remark 2: When we assume that two topological defects with a charge of
$q_{1}=q_{2}=+1/2$ are located at $z=\pm r$, $0<r<1$, the energy $F_{d}$ for
this case is given by

$\displaystyle F_{d}=-\frac{1}{2}[\log(1-r^{2})+\log(1+r^{2})+\log(2r)].$ (18)

Taking derivative with respect to $r$ gives

$\displaystyle\frac{\partial F_{d}}{\partial r}=\frac{5r^{4}-1}{2r(1-r^{4})}.$
(19)

Solving the equation $\partial F_{d}/\partial r=0$, we have $r=5^{-1/4}$,
which is consistent with the results given by Duclos et al. [11].

Remark 3: Our results show that the concept of this expression is identical to
that of Hamiltonians for multiple vortices. In classical fluid dynamics,
Hamiltonians are referred to as Kirchhoff-Routh path functions, which give the
vortex trajectories (see p. 121 in [22]). The theory of Kirchhoff-Routh path
functions was extended to multiply connected domains by Lin [19] and a
mathematical formulation for multiply connected domains with prime functions
was given in a pioneering work by Crowdy [23]. The concept presented here has
been used to derive a variational formulation for vortex dynamics based on the
principle of least action and asymptotic formulas for the kinematic energy
have been obtained [24]. A Hamiltonian $\mathcal{H}$ in the presence of
multiple point vortices with circulation $\Gamma_{k}$, $k=1,\ldots,N$ is given
by [22]

$\displaystyle\mathcal{H}=-\frac{1}{2}\sum_{k=1}^{N}\Gamma_{k}^{2}R(z_{k})-\frac{1}{2}\sum_{k=1}^{N}\sum_{k\neq
l}^{N}\Gamma_{k}\Gamma_{l}G(z_{k},z_{l}),$ (20)

where $G(.,.)$ is a hydrodynamic Green’s function for a point vortex and
$R(.)$ is a Robin function [25] for the Green’s function $G(.,.)$ and is
defined by

$\displaystyle R(z_{k})=\lim_{z\rightarrow
z_{k}}\left(G(z,z_{k})-\frac{1}{2\pi}\log|z-z_{k}|\right).$ (21)

Note that the Green’s function for a point vortex in the unit disc is given by

$\displaystyle
G(z,z_{k})\equiv\frac{1}{2\pi}\log\left|\frac{z-z_{k}}{|z_{k}|(z-1/\overline{z_{k}})}\right|.$
(22)

The imaginary part of $\mathcal{G}(z,z_{k})$ is similar to $-2\pi G(z,z_{k})$
except for the opposite sign of the second term in (10) due to the boundary
conditions. Our result indicates that the classical approach used in vortex
dynamics can be extended to the theory of nematic liquid crystals.

###### Proof.

For simplicity, we define

$\displaystyle r_{k}(z)\equiv f_{\mathbb{D}}(z|\mathcal{Z})+{\rm
i}q_{k}\log(z-z_{k}).$ (23)

This function is analytic at $z=z_{k}$ and it is easy to see that $r_{k}(z)$
satisfies

$\displaystyle\lim_{z\rightarrow
z_{k}}r_{k}(z)=q_{k}\mathcal{R}(z_{k})+\sum_{l\neq
k}^{N}q_{l}\mathcal{G}(z_{k},z_{l})-\frac{\pi}{2}.$ (24)

The branch cuts of $\log(z-z_{k})$ are chosen such that they do not intersect
each other. The path of the line integral is shown in Figure 3. In the figure,
$l_{k}^{+}$, $l_{k}^{-}$ denote the paths of the line integral from
$\partial\mathbb{D}$ to $\partial D_{k}(\epsilon)$ used to avoid the jumps of
the branch cuts of $\log(z-z_{k})$. It is noted that for $z^{\pm}\in
l_{k}^{\pm}$, the following property is always satisfied:

$\displaystyle\log(z^{+}-z_{k})-\log(z^{-}-z_{k})=2\pi{\rm i}.$ (25)

Let $a_{k}^{+}$ and $a_{k}^{-}$ denote the intersections of the paths
$l_{k}^{+}$ and $l_{k}^{-}$ with $\partial\mathbb{D}$, and $b_{k}^{+}$ and
$b_{k}^{-}$ denote the intersections of $l_{k}^{+}$ and $l_{k}^{-}$ with
$\partial D_{k}(\epsilon)$, respectively. We choose $0\leq{\rm
arg}(a_{1}^{\pm})<{\rm arg}(a_{2}^{\pm})<\ldots<{\rm arg}(a_{N}^{\pm})<2\pi$
as shown in Figure 2. Due to the existence of branch cuts, angle $\phi$ has
jumps at these points. We define the argument on the boundary $\theta$ as
shown in Figure 2 (a). For $z=e^{{\rm i}\theta}$, $\theta\in[0,2\pi)$,

$\displaystyle{\rm Re}[f(e^{{\rm
i}\theta}|\mathcal{Z})]=\phi(\theta)=\theta+\chi(\theta)-\frac{\pi}{2},$ (26)

where $\chi(\theta)$ is a step function whose value jumps at $a_{k}^{\pm}$ as
follows:

$\displaystyle\chi(\theta)=\left\\{\begin{aligned} &0,\quad&&{\rm for}\
\theta\in[0,{\rm arg}(a^{+}_{1})],\\\ &-\sum_{k=1}^{l}2\pi q_{k}\quad&&{\rm
for}\ \theta\in[{\rm arg}(a_{l}^{-}),{\rm arg}(a^{+}_{l+1})],\quad
l=1,\ldots,N-1,\\\ &-2\pi,\quad&&{\rm for}\ \theta\in[{\rm
arg}(a^{-}_{N}),2\pi).\end{aligned}\right.$ (27)

This decomposition is illustrated in Figure 2 (b).

Using the complex Green’s theorem, the line integral of the energy (13) is
given by

$\displaystyle F$ $\displaystyle=\frac{1}{2{\rm
i}}\left[\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial
f_{\mathbb{D}}}{\partial z}\overline{f_{\mathbb{D}}}{\rm
d}z+\sum_{k=1}^{N}\left(\int_{l_{k}^{+}}\frac{\partial
f_{\mathbb{D}}}{\partial z}\overline{f_{\mathbb{D}}}{\rm
d}z-\int_{l_{k}^{-}}\frac{\partial f_{\mathbb{D}}}{\partial
z}\overline{f_{\mathbb{D}}}{\rm d}z-\oint_{\partial
D_{k}(\epsilon)}\frac{\partial f_{\mathbb{D}}}{\partial
z}\overline{f_{\mathbb{D}}}{\rm d}z\right)\right],$ (28)

where $a_{N+1}^{+}=a_{1}^{+}$. We evaluate each term in equation (28) using
the functions defined by (23).

(i) Line integral along $\partial D_{k}(\epsilon)$

The arguments of the function $\log(z-z_{k})$ at $z=b_{k}^{\pm}$ are set to
${\rm arg}(b_{k}^{\pm}-z_{k})=\beta_{k}^{\pm}$. Using the change of variables
along $\partial D_{k}(\epsilon)$ such that $z=z_{k}+\epsilon e^{{\rm
i}\theta}$, $\theta\in[\beta^{-}_{k},\beta_{k}^{+}]$ and equation (23), the
line integral of $\partial D_{k}(\epsilon)$ is evaluated as follows:

$\displaystyle-\frac{1}{2{\rm i}}\oint_{\partial
D_{k}(\epsilon)}\frac{\partial f_{\mathbb{D}}}{\partial
z}\overline{f_{\mathbb{D}}}{\rm d}z$ $\displaystyle=-\frac{1}{2{\rm
i}}\int_{\beta_{k}^{-}}^{\beta_{k}^{+}}\left[-\frac{{\rm i}q_{k}}{\epsilon
e^{{\rm i}\theta}}+\frac{\partial r_{k}}{\partial z}(z_{k}+\epsilon e^{{\rm
i}\theta})\right]\left[{\rm i}q_{k}\log(\epsilon e^{-{\rm
i}\theta})+\overline{r_{k}(z_{k}+\epsilon e^{{\rm i}\theta})}\right]{\rm
i}\epsilon e^{{\rm i}\theta}{\rm d}\theta$
$\displaystyle=-\frac{1}{2}\int_{\beta_{k}^{-}}^{\beta_{k}^{+}}\left[-{\rm
i}q_{k}\overline{r_{k}(z_{k})}+q_{k}^{2}\log\epsilon-{\rm
i}q_{k}^{2}\theta\right]{\rm d}\theta+\mathcal{O}(\epsilon\log\epsilon)$
$\displaystyle=-\pi q_{k}^{2}\log\epsilon+{\rm i}\pi
q_{k}\overline{r_{k}(z_{k})}+{\rm i}\pi^{2}q_{k}^{2}+{\rm i}\pi
q_{k}^{2}\beta_{k}^{-}+\mathcal{O}(\epsilon\log\epsilon).$ (29)

where we used that $r_{k}(z_{k}+\epsilon e^{{\rm
i}\theta})=r_{k}(z_{k})+\mathcal{O}(\epsilon)$ and
$\beta_{k}^{+}-\beta_{k}^{-}=2\pi$ from the second to third line.

(ii) Integrals along $l_{k}^{+}$ and $l_{k}^{-}$

The line integrals along $l_{k}^{\pm}$ are given by

$\displaystyle\frac{1}{2{\rm i}}\int_{l_{k}^{\pm}}\frac{\partial
f_{\mathbb{D}}}{\partial z}\overline{f_{\mathbb{D}}}{\rm d}z=\frac{1}{2{\rm
i}}\int_{a_{k}^{\pm}}^{b_{k}^{\pm}}\left(-\frac{{\rm
i}q_{k}}{z-z_{k}}+\frac{\partial r_{k}}{\partial z}\right)\left({\rm
i}q_{k}\log\overline{(z^{\pm}-z_{k})}+\overline{r_{k}(z)}\right){\rm d}z,$
(30)

where $z^{\pm}$ are points on $l_{k}^{\pm}$. The line integral along
$l_{k}^{+}$ minus the integral along $l_{k}^{-}$ becomes

$\displaystyle\frac{1}{2{\rm i}}\int_{l_{k}^{+}}\frac{\partial
f_{\mathbb{D}}}{\partial z}\overline{f_{\mathbb{D}}}{\rm d}z$
$\displaystyle-\frac{1}{2{\rm i}}\int_{l_{k}^{-}}\frac{\partial
f_{\mathbb{D}}}{\partial z}\overline{f_{\mathbb{D}}}{\rm d}z=2\pi
q_{k}\cdot\frac{1}{2{\rm i}}\int_{a_{k}^{-}}^{b_{k}^{-}}\left[-\frac{{\rm
i}q_{k}}{z-z_{k}}+\frac{\partial r_{k}}{\partial z}\right]{\rm d}z$
$\displaystyle=-{\rm i}\pi q_{k}\left[-{\rm
i}q_{k}(\log(b_{k}^{-}-z_{k})-\log(a_{k}^{-}-z_{k}))+r_{k}(b_{k}^{-})-r_{k}(a_{k}^{-})\right].$
(31)

Due to the fact that $\log(b_{k}^{-}-z_{k})=\log\epsilon+{\rm i}\beta_{k}^{-}$
and $r_{k}(b_{k}^{-})=r_{k}(z_{k})+\mathcal{O}(\epsilon)$ as
$\epsilon\rightarrow 0$, equation (31) becomes

$\displaystyle\frac{1}{2{\rm i}}\int_{l_{k}^{+}}\frac{\partial
f_{\mathbb{D}}}{\partial z}\overline{f_{\mathbb{D}}}{\rm d}z$
$\displaystyle-\frac{1}{2{\rm i}}\int_{l_{k}^{-}}\frac{\partial
f_{\mathbb{D}}}{\partial z}\overline{f_{\mathbb{D}}}{\rm d}z=-\pi
q_{k}^{2}\log\epsilon-{\rm i}\pi q_{k}^{2}\beta_{k}^{-}-{\rm i}\pi
q_{k}r_{k}(z_{k})+{\rm i}\pi
q_{k}f_{\mathbb{D}}(a_{k}^{-}|\mathcal{Z})+\mathcal{O}(\epsilon).$ (32)

(iii) Line integral along the unit circle

First we note that $\partial f_{\mathbb{D}}/\partial z$ is analytic everywhere
in $\mathbb{D}\backslash\mathcal{Z}$. Using equation (26), the line integral
of the unit circle becomes

$\displaystyle\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial
f_{\mathbb{D}}}{\partial z}\overline{f_{\mathbb{D}}}{\rm d}z=\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial
f_{\mathbb{D}}}{\partial z}(2\theta+2\chi(\theta)-\pi-f_{\mathbb{D}}){\rm
d}z.$ (33)

Let the argument of $a_{k}^{\pm}$ define $\alpha_{k}^{\pm}$. Note that
$\chi(\alpha_{k}^{+})-\chi(\alpha_{k}^{-})=-2\pi q_{k}$. The first two terms
in (33) are

$\displaystyle\frac{1}{{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}(\theta+\chi(\theta))\frac{\partial
f_{\mathbb{D}}}{\partial z}{\rm d}z=-{\rm
i}\sum_{k=1}^{N}[(\theta+\chi(\theta))f_{\mathbb{D}}(e^{{\rm
i}\theta}|\mathcal{Z})]_{\theta=\alpha_{k}^{-}}^{\theta=\alpha_{k+1}^{+}}+{\rm
i}\sum_{k=1}^{N}\int_{\alpha_{k}^{-}}^{\alpha_{k+1}^{+}}f_{\mathbb{D}}(e^{{\rm
i}\theta}|\mathcal{Z}){\rm d}\theta$ (34) $\displaystyle=-{\rm
i}\sum_{k=1}^{N}(2\pi
q_{k}f_{\mathbb{D}}(a_{k}^{-}|\mathcal{Z})+4\pi^{2}q_{k}^{2}+2\pi
q_{k}(\chi(\alpha_{k}^{-})+\alpha_{k}^{+}))+{\rm
i}\sum_{k=1}^{N}(2\pi^{2}q_{k}^{2}+2\pi
q_{k}(\chi(\alpha_{k}^{-})+\alpha_{k}^{+}))-{\rm i}\pi^{2}$
$\displaystyle=-{\rm i}\sum_{k=1}^{N}(2\pi
q_{k}f_{\mathbb{D}}(a_{k}^{-}|\mathcal{Z})+2\pi^{2}q_{k}^{2})-{\rm i}\pi^{2},$
(35)

where the second term in (34) is calculated as

$\displaystyle{\rm
i}\sum_{k=1}^{N}\int_{\alpha_{k}^{-}}^{\alpha_{k+1}^{+}}f_{\mathbb{D}}(e^{{\rm
i}\theta}|\mathcal{Z}){\rm d}\theta$ $\displaystyle={\rm
i}\sum_{k=1}^{N}\int_{\alpha_{k}^{-}}^{\alpha_{k+1}^{+}}{\rm
Re}[f_{\mathbb{D}}(e^{{\rm i}\theta}|\mathcal{Z})]{\rm d}\theta=-{\rm
i}\sum_{k=1}^{N}\int_{\alpha_{k}^{-}}^{\alpha_{k+1}^{+}}\left(\theta+\chi(\theta)-\frac{\pi}{2}\right){\rm
d}\theta.$ (36)

For the first equality in (36), we used the fact that when $|t|<1$,

$\displaystyle\int_{0}^{2\pi}\log|e^{{\rm
i}\theta}-t|d\theta=\int_{0}^{2\pi}\log|1-\overline{t}e^{{\rm
i}\theta}|d\theta=0,$ (37)

which can be proven by the residue theorem for $|t|<1$:

$\displaystyle\oint_{\partial\mathbb{D}}\frac{\log(1-tz)}{z}{\rm d}z=0.$ (38)

The third term in (33) is

$\displaystyle-\frac{\pi}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial
f_{\mathbb{D}}}{\partial z}{\rm d}z=-\frac{\pi}{2{\rm
i}}\sum_{k=1}^{N}\left[f_{\mathbb{D}}(z|\mathcal{Z})\right]_{z=a_{k}^{-}}^{z=a_{k+1}^{+}}={\rm
i}\pi^{2}\sum_{k=1}^{N}q_{k}={\rm i}\pi^{2}.$ (39)

where we used the fact that $f_{\mathbb{D}}$ shifts $2\pi q_{k}$ at
$z=a_{k}^{\pm}$. The fourth term in (33) is

$\displaystyle-\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial
f_{\mathbb{D}}}{\partial z}f_{\mathbb{D}}{\rm d}z$
$\displaystyle=-\frac{1}{4{\rm
i}}\sum_{k=1}^{N}\left[f^{2}_{\mathbb{D}}(z|\mathcal{Z})\right]_{z=a_{k}^{-}}^{z=a_{k+1}^{+}}=-\frac{1}{4{\rm
i}}\sum_{k=1}^{N}2\pi q_{k}(2f_{\mathbb{D}}(a_{k}^{-}|\mathcal{Z})+2\pi
q_{k}).$

Hence, the line integral around $\partial\mathbb{D}$ in (33) is given by

$\displaystyle\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial
f_{\mathbb{D}}}{\partial z}\overline{f_{\mathbb{D}}}{\rm d}z$
$\displaystyle=-{\rm
i}\pi\sum_{k=1}^{N}(q_{k}f_{\mathbb{D}}(a_{k}^{-}|\mathcal{Z})+\pi
q_{k}^{2}).$ (40)

Combining (29), (32), and (40), we have

$\displaystyle F=-2\pi\sum_{k=1}^{N}q_{k}^{2}\log\epsilon+2\pi
F_{d}+\mathcal{O}(\epsilon\log\epsilon),\quad F_{d}=\sum_{k=1}^{N}q_{k}{\rm
Im}[r_{k}(z_{k})].$ (41)

Using the formula (24), the final expression (15) is derived. ∎

###  3.2 Analytical formula for arbitrarily shaped domain

The theory of energy calculation on the unit disc shown in the previous
section is utilized to derive an explicit formula for energy in arbitrarily
shaped regions.

######  Theorem 3.2.

We consider an arbitrarily shaped simply connected domain $\Omega$ in the
$w$-plane mapped from the unit disc $\mathbb{D}$ by a conformal mapping
represented as $w=g(z)$, where the boundary of the domain is sufficiently
smooth and $g^{\prime}(z)\neq 0$ for $z\in\mathbb{D}$. Assume that there exist
$N$ topological defects with an integer or half-integer charge located at
$w_{k}\in\Omega$, $k=1,\ldots,N$, $w_{k}\notin\partial\Omega$, and define the
set of topological defects as $\mathcal{W}\equiv\\{w_{k}|k=1,\ldots,N\\}$. We
also define the set of topological defects
$\mathcal{Z}\equiv\\{z_{k}|k=1,\ldots,N\\}$, $z_{k}=g^{-1}(w_{k})$, mapped
from $\mathcal{W}$ and the small disc with radius $\epsilon$ which encloses
only $w_{k}$ as $D_{k}(\epsilon)\equiv\\{w_{k}+re^{{\rm i}\theta}|\ 0\leq
r\leq\epsilon,0\leq\theta\leq 2\pi\\}$. The sum of topological charges
satisfies (6). The total energy $F$ of the region $\Omega\backslash
D(\epsilon)$, $D(\epsilon)\equiv\cup_{k=1}^{N}D_{k}(\epsilon)$, is written as
follows:

$\displaystyle F=-2\pi\sum_{k=1}^{N}q_{k}^{2}\log\epsilon+2\pi
F_{d}+F_{g}+\mathcal{O}(\epsilon\log\epsilon),$ (42)

where

$\displaystyle
F_{d}\equiv\sum_{k=1}^{N}\left(q_{k}^{2}-2q_{k}\right)\log|g^{\prime}(z_{k})|+\sum_{k=1}^{N}q_{k}^{2}{\rm
Im}[\mathcal{R}(z_{k})]+\sum_{k=1}^{N}\sum_{l\neq k}^{N}q_{k}q_{l}{\rm
Im}[\mathcal{G}(z_{k},z_{l})],$ (43) $\displaystyle
F_{g}\equiv\int_{\mathbb{D}}\left|\frac{g^{\prime\prime}(z)}{g^{\prime}(z)}\right|^{2}{\rm
d}S+4\pi\log|g^{\prime}(0)|.$ (44)

The energy $F_{d}$ depends on $z_{k}$ and the map $g(z)$, whereas $F_{g}$ is
independent of $z_{k}$. The energy $F_{g}$ can be seen as the total energy of
the map $g(z)$ in $\mathbb{D}$.

Figure 4: Line integral along $\partial(\Omega\backslash D(\epsilon))$ for
case $N=2$. Branch cuts of logarithmic singularities $z_{k}$ are chosen as
$l_{k}$ and argument of $\log(z-z_{k})$ jumps on $l_{k}$ by $2\pi$. The path
of integrals in the $w$-plane is transformed into the line integral in the
$z$-plane under a conformal map $w=g(z)$.

###### Proof.

Since the energy defined by equation (2) is conformally invariant, the Frank
free energy in the domain $\Omega$ is expressed by line integrals:

$\displaystyle F$ $\displaystyle=\frac{1}{2{\rm
i}}\oint_{\partial(\Omega\backslash D(\epsilon))}\frac{\partial
f_{\Omega}}{\partial
w}(g^{-1}(w)|\mathcal{Z})\overline{f_{\Omega}(g^{-1}(w)|\mathcal{Z})}{\rm
d}w=\frac{1}{2{\rm
i}}\oint_{\partial(\mathbb{D}\backslash\tilde{D}(\epsilon))}\frac{\partial
f_{\Omega}}{\partial z}(z|\mathcal{Z})\overline{f_{\Omega}(z|\mathcal{Z})}{\rm
d}z$ $\displaystyle=\frac{1}{2{\rm
i}}\left[\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial
f_{\Omega}}{\partial z}\overline{f_{\Omega}}{\rm
d}z+\sum_{k=1}^{N}\left(\int_{l_{k}^{+}}\frac{\partial f_{\Omega}}{\partial
z}\overline{f_{\Omega}}{\rm d}z-\int_{l_{k}^{-}}\frac{\partial
f_{\Omega}}{\partial z}\overline{f_{\Omega}}{\rm
d}z-\oint_{\partial\tilde{D}_{k}(\epsilon)}\frac{\partial f_{\Omega}}{\partial
z}\overline{f_{\Omega}}{\rm d}z\right)\right],$ (45)

where the points $a_{k}^{\pm}$, $b_{k}^{\pm}$, and the paths $l_{k}^{\pm}$ are
chosen in the same manner as that in the unit disc case and
$\tilde{D}_{k}(\epsilon)\equiv g^{-1}(D_{k}(\epsilon))$. Figure 4 shows the
paths of the line integrals. The point on $\partial D_{k}(\epsilon)$ is
parameterized as $w=w_{k}+\epsilon e^{{\rm i}\theta}$, and
$\theta\in[0,2\pi]$, which is transformed into

$\displaystyle g^{-1}(w_{k}+\epsilon e^{{\rm
i}\theta})=z_{k}+\frac{1}{g^{\prime}(z_{k})}\epsilon e^{{\rm
i}\theta}+\mathcal{O}(\epsilon^{2})$ (46)

under the inverse of the conformal map of $g(z)$. This implies that the disc
at $w=w_{k}=g(z_{k})$ with radius $\epsilon$ is mapped to the disc at
$z=z_{k}$ with radius $\tilde{\epsilon}\equiv\epsilon/|g^{\prime}(z_{k})|$ up
to $\mathcal{O}(\epsilon^{2})$.

Now we use the same technique as that used for the unit disc. Before proving
this, it is convenient to define $h(z)\equiv-{\rm i}\log g^{\prime}(z)$ and a
function that excludes only the singularity at $z=z_{k}$ as

$\displaystyle\tilde{r}_{k}(z)\equiv f_{\Omega}(z|\mathcal{Z})+{\rm
i}\log(z-z_{k})=r_{k}(z)+h(z).$ (47)

(i) Line integral along $\partial\tilde{D}_{k}(\epsilon)$

The integral along $\partial\tilde{D}_{k}(\epsilon)$ follows the calculations
for the unit disc except for two simple modifications, namely, the inner
radius $\epsilon$ on the $w$-plane is $\tilde{\epsilon}$ on the $z$-plane and
the function $r_{k}(z)$ analytic at $z=z_{k}$ is $\tilde{r}_{k}(z)$ in this
case. Following equation (29), the line integral around
$\tilde{D}_{k}(\epsilon)$ is

$\displaystyle-\frac{1}{2{\rm
i}}\oint_{\partial\tilde{D}_{k}(\epsilon)}\frac{\partial f_{\Omega}}{\partial
z}\overline{f_{\Omega}}dz=-\pi q_{k}^{2}\log\tilde{\epsilon}+{\rm i}\pi
q_{k}\overline{\tilde{r}_{k}(z_{k})}+{\rm i}\pi^{2}q_{k}^{2}+{\rm i}\pi
q_{k}^{2}\beta_{k}^{-}+\mathcal{O}(\epsilon\log\epsilon)$ $\displaystyle=-\pi
q_{k}^{2}\log\epsilon+\pi q_{k}^{2}\log|g^{\prime}(z_{k})|+{\rm i}\pi
q_{k}\overline{r_{k}(z_{k})}-\pi q_{k}\log\overline{g^{\prime}(z_{k})}+{\rm
i}\pi^{2}q_{k}^{2}+{\rm i}\pi
q_{k}^{2}\beta_{k}^{-}+\mathcal{O}(\epsilon\log\epsilon).$

(ii) Integrals along branch cuts $l_{k}^{+}$ and $l_{k}^{-}$

This integral is calculated using the same technique as that used for the unit
disc:

$\displaystyle\frac{1}{2{\rm i}}\int_{l_{k}^{+}}\frac{\partial
f_{\Omega}}{\partial z}\overline{f_{\Omega}}{\rm d}z-\frac{1}{2{\rm
i}}\int_{l_{k}^{-}}\frac{\partial f_{\Omega}}{\partial
z}\overline{f_{\Omega}}{\rm d}z$ $\displaystyle=-\pi q_{k}^{2}\log\epsilon+\pi
q_{k}^{2}\log|g^{\prime}(z_{k})|-{\rm i}\pi q_{k}^{2}\beta_{k}^{-}-{\rm i}\pi
q_{k}r_{k}(z_{k})-\pi q_{k}\log g^{\prime}(z_{k})+{\rm i}\pi
q_{k}f_{\Omega}(a_{k}^{-}|\mathcal{Z})+\mathcal{O}(\epsilon).$

(iii) Line integral along the unit circle

The line integral along the unit circle is given by

$\displaystyle\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial
f_{\Omega}}{\partial z}\overline{f_{\Omega}}{\rm d}z$
$\displaystyle=\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\left(\frac{\partial
f_{\mathbb{D}}}{\partial z}+\frac{\partial h}{\partial
z}\right)\overline{(f_{\mathbb{D}}+h)}{\rm d}z$ $\displaystyle=\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\left(\frac{\partial
f_{\mathbb{D}}}{\partial z}\overline{f_{\mathbb{D}}}+\frac{\partial
f_{\mathbb{D}}}{\partial z}\overline{h}+\frac{\partial h}{\partial
z}\overline{f_{\mathbb{D}}}+\frac{\partial h}{\partial
z}\overline{h}\right){\rm d}z.$ (48)

Using the same argument as that for the unit disc, the first term in (48) is
identical to (40); that is,

$\displaystyle\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial
f_{\mathbb{D}}}{\partial z}\overline{f_{\mathbb{D}}}{\rm d}z=-{\rm
i}\pi\sum_{k=1}^{N}(q_{k}f_{\mathbb{D}}(a_{k}^{-}|\mathcal{Z})+\pi
q_{k}^{2}).$ (49)

The last term in (48) is

$\displaystyle\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial h}{\partial
z}\overline{h}{\rm
d}z=\int_{\mathbb{D}}\left|\frac{g^{\prime\prime}}{g^{\prime}}\right|^{2}{\rm
d}x{\rm d}y,$ (50)

where the complex Green’s theorem and the fact that $h(z)$ is analytic in
$\mathbb{D}$ are used. The second term in (48) is

$\displaystyle\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial
f_{\mathbb{D}}}{\partial z}\overline{h}{\rm
d}z=-\frac{1}{2}\oint_{\partial\mathbb{D}}\sum_{k=1}^{N}q_{k}\left(\frac{1}{z-z_{k}}-\frac{\overline{z_{k}}}{1-\overline{z_{k}}z}\right)\overline{h}{\rm
d}z={\rm i}\pi\sum_{k=1}^{N}q_{k}\overline{h(z_{k})}-2{\rm
i}\pi\overline{h(0)},$ (51)

where we used $\overline{z}=1/z$, ${\rm d}\overline{z}=-1/z^{2}{\rm d}z$ on
$\partial\mathbb{D}$, and the following residue theorem:

$\displaystyle\frac{1}{2}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\overline{h(z)}}{z-z_{k}}{\rm
d}z=\frac{1}{2}\oint_{\partial\mathbb{D}}\frac{\overline{h(z)}}{z-z_{k}}{\rm
d}z=\overline{\frac{1}{2}\oint_{\partial\mathbb{D}}\frac{h(z)}{1/z-\overline{z_{k}}}\left(-\frac{1}{z^{2}}\right){\rm
d}z}={\rm i}\pi\overline{h(0)},$ (52)
$\displaystyle-\frac{\overline{z_{k}}}{2}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\overline{h(z)}}{1-\overline{z_{k}}z}{\rm
d}z=-\frac{\overline{z_{k}}}{2}\overline{\oint_{\partial\mathbb{D}}\frac{h(z)}{1-z_{k}/z}\left(-\frac{1}{z^{2}}\right){\rm
d}z}=-{\rm i}\pi\overline{h(z_{k})}+{\rm i}\pi\overline{h(0)}.$ (53)

The third term is calculated based on the property that
$\overline{f_{\mathbb{D}}}=2\theta+2\chi(\theta)-\pi-f_{\mathbb{D}}$ on the
unit circle:

$\displaystyle\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial h}{\partial
z}\overline{f_{\mathbb{D}}}{\rm d}z=\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial h}{\partial
z}(2\theta+2\chi(\theta)-\pi-f_{\mathbb{D}}){\rm d}z.$ (54)

The first two terms in (54) are explicitly calculated as follows:

$\displaystyle\frac{1}{{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial h}{\partial
z}(\theta+\chi(\theta)){\rm d}z=-2{\rm
i}\pi\sum_{k=1}^{N}q_{k}h(a_{k}^{-})+{\rm i}\int_{0}^{2\pi}h(e^{{\rm
i}\theta}){\rm d}\theta=-2{\rm i}\pi\sum_{k=1}^{N}q_{k}h(a_{k}^{-})+2{\rm
i}\pi h(0).$

The third term in (54) is zero because $h(z)$ is analytic in $\mathbb{D}$. The
fourth term in (54) is

$\displaystyle-\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial h}{\partial
z}f_{\mathbb{D}}{\rm d}z$ $\displaystyle=-\frac{1}{2{\rm
i}}\sum_{k=1}^{N}[h(z)f_{\mathbb{D}}]_{a_{k}^{-}}^{a_{k+1}^{+}}+\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial
f_{\mathbb{D}}}{\partial z}h{\rm d}z$ $\displaystyle={\rm
i}\pi\sum_{k=1}^{N}q_{k}h(a_{k}^{-})-{\rm i}\pi\sum_{k=1}^{N}q_{k}h(z_{k}),$
(55)

where we used the residue theorem. Hence, the line integral along the unit
circle is given by

$\displaystyle\frac{1}{2{\rm
i}}\sum_{k=1}^{N}\int_{a_{k}^{-}}^{a_{k+1}^{+}}\frac{\partial
f_{\Omega}}{\partial z}\overline{f_{\Omega}}{\rm d}z$ $\displaystyle=-{\rm
i}\pi\sum_{k=1}^{N}q_{k}h(a_{k}^{-})-2{\rm i}\pi\sum_{k=1}^{N}q_{k}{\rm
Im}[h(z_{k})]+4\pi\log|g^{\prime}(0)|$ $\displaystyle-{\rm
i}\pi\sum_{k=1}^{N}q_{k}f_{\mathbb{D}}(a_{k}^{-}|\mathcal{Z})-{\rm
i}\pi^{2}\sum_{k=1}^{N}q_{k}^{2}+\int_{\mathbb{D}}\left|\frac{g^{\prime\prime}}{g^{\prime}}\right|^{2}{\rm
d}x{\rm d}y.$ (56)

Combining (i), (ii), and (iii), the energy for $f_{\Omega}(z|\mathcal{Z})$
given by (42) is derived.

∎

##  4 Topological defects localization on quadrature domains

Section 3 presented the explicit formula for the Frank free energy of
arbitrarily shaped domains. The complex potential $f_{\Omega}$ was split into
analytic and logarithmic parts and the line integrals along boundaries were
evaluated up to $\epsilon\log\epsilon$. The formulas derived here are
explicit, written as logarithmic singularities plus the contribution from
conformal maps. This section applies the formulas (42) to validate a
conjecture given by Ienaga et al. [10], namely, that the positions of defects
in both doublets and triplets linearly correlate with the distances between
the centres of the overlapping discs.

The simplest extension of localizing topological defects with boundaries, as
done by Duclos et al. [11], is to estimate defect positions in regions with
multiple discs that intersect each other.  When two or three discs with the
same radius overlap at a distance $\Delta$, they are called “doublets” and
“triplets”, respectively. These geometries exhibit ferromagnetic and
antiferromagnetic order in bacterial vortex lattices, depending on the
distance between the centres of the discs relative to the radius of the discs
[26]. Ienaga et al. created regions with two or three circular areas and
conducted cell alignment experiments within these regions to investigate the
number and positioning of defects [10]. They hypothesized a proportional
relationship between the distance $\Delta$ between the centres of two
overlapping discs and the distance $D_{++}$ between two $+1/2$ defects within
those doublets. They conjectured the following equation for doublets through
fitting:

$\displaystyle D_{++}=2\left(D_{++,0}+\frac{\Delta}{2}\right),\quad
D_{++,0}=0.58.$ (57)

For triplets, it was experimentally confirmed that when $\Delta$ is large, a
topological defect with a charge of $-1/2$ exists at the centre, and three
topological defects with a charge of $+1/2$ are located on each disc. They
also hypothesized that the distance $D_{+-}$ between the $-1/2$ defect and the
$+1/2$ defect, and the distance between disc centres $\Delta$, have a
proportional relationship with a slope of $1/\sqrt{3}$.

In this paper, we analytically verify these hypotheses using the proposed
formula and the mapping of quadrature domains [20]. The conformal mapping for
a region with $M$-overlapping unit discs can be approximated by the following
rational function:

$\displaystyle
w=g(z)=\frac{a^{2M}-1}{a^{M}-a^{M-1}+1}\cdot\frac{z}{a^{M}-z^{M}},$ (58)

where the real parameter $a$ should be greater than $1$ and $M$ is a positive
integer. This region is known as the quadrature domains. Such domains have
been extensively studied in fields such as Hele-Shaw flows [18, 20].

Figure 5: Examples of quadrature domains with various values of $M$ and $a$.
Black dots correspond to $w=g(\omega_{M}^{m}/a)$,
$\omega_{M}\equiv\exp(2\pi{\rm i}/M)$, $m=0,1,\ldots,M-1$, which are
considered as the centres of the discs. Red circles are unit discs with
centres at the black dots. (i) $a=1.05$, $M=2$, (ii) $a=1.4$, $M=3$, and (iii)
$a=1.5$, $M=5$.  Figure 6: Example of quadrature domain with $M=3$ and two
stationary topological defects. Two topological defects with a charge of
$+1/2$ exist inside the unit disc. Because of the potential $f_{\Omega}$,
topological defects with a charge of $+1/2$ are also located at
$1/\overline{z_{n}}$, $n=1,2$. For $M=3$, $g^{\prime}(z)$ has three simple
zeros and three double poles outside the unit disc, which indicates that the
correction term $h(z)=-{\rm i}\log g^{\prime}(z)$ produces topological defects
with a charge of $+1$ and $-2$. The problem of finding the stationary
positions of defects is equivalent to finding an equilibrium of ten
topological defects in the infinite plane.

Quadrature domains are applied to calculate topological defects for the first
time in the field of active nematics. The advantage of using the rational
function (58) is that it enables the boundary of doublets or triplets to be
represented by changing a single parameter $a$. Figure 5 shows the results of
mappings for various values of $M$ and $a$. Through this mapping, the origin
of the unit disc in the $z$-plane is mapped to the origin of the $w$-plane,
corresponding to the centre of the region in the $w$-plane. The first and
second derivatives of this function are given as follows:

$\displaystyle
g^{\prime}(z)=\frac{a^{2M}-1}{a^{M}-a^{M-1}+1}\cdot\frac{(a^{M}+(M-1)z^{M})}{(a^{M}-z^{M})^{2}},$
(59) $\displaystyle
g^{\prime\prime}(z)=\frac{a^{2M}-1}{a^{M}-a^{M-1}+1}\cdot\frac{Mz^{M-1}[(M-1)z^{M}+(M+1)a^{M}]}{(a^{M}-z^{M})^{3}}.$
(60)

The mapping $g(z)$ may exhibit self-intersections for certain value of $a$.
The condition for no self-intersections is that the turning number of $g(z)$
around the unit disc is $1$, that is,

$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}\frac{{\rm d}}{{\rm d}\theta}{\rm
Im}\left[\log g^{\prime}(e^{{\rm i}\theta})\right]{\rm d}\theta=1,$ (61)

which results in no zeros of $g^{\prime}(z)$ being in the unit disc, i.e.,
$a>(M-1)^{1/M}$.

The centres of the discs are well approximated by $g(\omega^{m}_{M}/a)$,
$m=0,1,\ldots,M-1$, $\omega_{M}=\exp(2\pi{\rm i}/M)$. The coefficients of the
mapping in (58) are chosen such that the radius $R=g(1)-g(1/a)$ equals $1$.

The advantage of the map is that $g(z)$ is rational, allowing the Frank free
energy in the presence of topological defects to be expressed as logarithm of
certain rational functions. This enables the calculation of stationary defect
locations and stability using standard calculus. We note that

$\displaystyle-{\rm i}\log g^{\prime}(z)=-{\rm i}\log(M-1)-{\rm
i}\sum_{m=1}^{M}\log(z+a\omega_{M}^{m}/(M-1)^{1/M})+2{\rm
i}\sum_{m=1}^{M}\log(z-a\omega_{M}^{m}),$ (62)

where $\omega_{M}\equiv\exp(2\pi{\rm i}/M)$. This indicates that the problem
of finding two stationary defects considered here is to find the stability of
four topological defects with a charge of $+1/2$, $M$ topological defects with
a charge of $+1$ located at $z=-a\omega_{M}^{m}/(M-1)^{1/M}$,
$m=0,\ldots,M-1$, and topological defects with a charge of $-2$ located at
$z=a\omega^{m}_{M}$, $m=0,\ldots,M-1$ in the infinite plane. The configuration
for $M=3$ is illustrated in Figure 6.

###  4.1 $M=2$: Doublet domain

We first investigate the defect positions and their alignment angles in
doublet regions. The condition for no self-intersections is $a>1$ and

$\displaystyle g(z)=\frac{a^{4}-1}{a^{2}-a+1}\cdot\frac{z}{a^{2}-z^{2}}.$ (63)

Since the centres of the discs of doublet regions are represented as $g(\pm
1/a)$, the distance between the centres of the two discs is given by

$\displaystyle\Delta=2g(1/a)=\frac{2a}{a^{2}-a+1}.$ (64)

Using the symmetry in the region where the two discs intersect, we assume that
there exist topological defects with charges $q_{1}=q_{2}=+1/2$ at two
symmetric points $\pm g(s)$ along the $u$-axis of the $w=u+{\rm i}v$-plane,
where $0<s<1$. In other words, the positions of topological defects in the
predomain are $z=\pm s$. Using equation (42), the free energy that depends
solely on the positions of the defects is explicitly given by

$\displaystyle F_{d}$
$\displaystyle=-\frac{1}{2}\log[2s(1-s^{4})]-\frac{3}{4}\log|g^{\prime}(s)|-\frac{3}{4}\log|g^{\prime}(-s)|.$
(65)

To find the extrema of the free energy, consider the derivative with respect
to $s$:

$\displaystyle\frac{\partial F_{d}}{\partial
s}=-\frac{3g^{\prime\prime}(s)}{4g^{\prime}(s)}+\frac{3g^{\prime\prime}(-s)}{4g^{\prime}(-s)}+\frac{5s^{4}-1}{2s(1-s^{4})}=\frac{3s(s^{2}+3a^{2})}{s^{4}-a^{4}}+\frac{5s^{4}-1}{2s(1-s^{4})}=0.$
(66)

This gives an $8$th-degree polynomial equation with only even powers of $s$:

$\displaystyle
p_{2}(t)=t^{4}+18a^{2}t^{3}+5(a^{4}-1)t^{2}-18a^{2}t-a^{4}=0,\quad t=s^{2}.$
(67)

Using the fact that $p_{2}(0)=-a^{4}<0$ and $p_{2}(1)=4(a^{4}-1)>0$,
$p_{2}(t)$ has at least one real solution in the interval $0<t<1$, and

$\displaystyle p_{2}^{\prime}(t)=4t^{3}+54a^{2}t^{2}+10(a^{4}-1)t-18a^{2},$
(68) $\displaystyle
p_{2}^{\prime\prime}(t)=12t^{2}+108a^{2}t+10(a^{4}-1)>0,\quad 0<t<1,$ (69)

which indicates that $p_{2}(t)$ has only one real solution in the interval
$0<t<1$. We denote this solution as $t_{0}$ and define $s_{0}^{2}=t_{0}$,
$s_{0}>0$. We solve the polynomial equation (67) numerically between $0<s<1$.
We refer to this method as the polynomial method.

Figure 7: Colour map of $\phi$ created by $+1/2$ defect pairs (red circles) in
quadrature domain with $M=2$. The defect locations are determined through
energy minimization under the assumption that two $+1/2$ defects symmetrically
located at $z=\pm s$ on the real axis of the predomain. Black dots correspond
to $w=g(1/a)$, which are assumed to be the centres of the discs of doublets
based on numerical experiments.

Figure 7 shows the locations of topological defects and colour maps of $\phi$
in $M=2$ with (i) $a=1.05$, (ii) $=1.5$, and $a=2.1$. The defect locations are
determined by using the polynomial method. Figure 8 shows numerical results
for the distance between two defects $D_{++}$ for various values of $\Delta$.
The solid black line represents the results obtained using the polynomial
method. The results conjectured by [10] (red dashed line) are in good
agreement with our numerical results. However, it can be observed that the
formula expected by Ienaga et al. is slightly shifted towards smaller values.

Figure 8: (Left) Relationship between distance between two discs $\Delta$ and
distance between two defects $D_{++}$. Black line is the numerical results
obtained using the polynomial method and the red dotted line corresponds to
the conjecture given in [10]. The polynomial is well approximated by the blue
dotted line given by asymptotic expansion near $\Delta=2$. (Right)
Relationship between $\Delta$ and $D_{++}$, $D_{+-}$ for triplets. We checked
that the slope of the black line at $\Delta=1.9$ is $0.7232$, which is close
to the value of $1/\sqrt{3}$ conjectured in [10]. Our numerical results are
within the error bars in [10].

###  4.2 Asymptotic expansions of $D_{++}$ and $\Delta$ for $M=2$

We demonstrate that the distance between two defects is proportional to the
distance between the centres of the two circles in the limit $a\rightarrow 1$
(i.e., when the two discs are almost touching each other). For $a=1+k$ with a
small parameter $k$, the distance $\Delta$ is expanded with respect to $k$ as
follows:

$\displaystyle\Delta=\frac{2a}{a^{2}-a+1}=2(1-k^{2}+k^{3})+\mathcal{O}(k^{4}).$
(70)

The expansion of the root of the polynomial (67) is given by

$\displaystyle
s_{0}=1-\frac{1}{5}k+\frac{3}{25}k^{2}-\frac{99}{625}k^{3}+\frac{32}{125}k^{4}+\mathcal{O}(k^{5}),$
(71)

which gives

$\displaystyle D_{++}$
$\displaystyle=2g(s_{0})=\frac{10}{3}-\frac{158}{75}k^{2}+\frac{2576}{1125}k^{3}+\mathcal{O}(k^{4})$
$\displaystyle=\frac{92}{75}+\frac{79}{75}\Delta+\frac{206}{1125}k^{3}+\mathcal{O}(k^{4}).$
(72)

It is important to note that $92/75=1.227\sim 2D_{++,0}=1.16$, which well
agrees with the formula conjectured by Ienaga et al. It is also observed that
$D_{++}$ and $\Delta$ are proportional for sufficiently small values of
parameter $k$ and that the approximation accuracy is good because the
coefficient of $k^{3}$ is small. The proportional equation, which includes
only the first two terms in (72) here, is referred to as asymptotic expansion.

The blue dotted line in Figure 8 shows the result of asymptotic expansion.
This line well matches the numerical results for wide ranges of $\Delta$.

###  4.3 $M=3$: Triplet domain

Ienaga et al. also investigated the positions of defects and the distances
between the centres of three overlapping discs (i.e., triplets) [10]. They
experimentally demonstrated that, in this case as well, there exists a
proportionality between the positions of defects and the circular regions.

For $M=3$,

$\displaystyle
g(z)=\frac{a^{6}-1}{a^{3}-a^{2}+1}\cdot\frac{z}{a^{3}-z^{3}},\quad
g^{\prime}(z)=\frac{a^{6}-1}{a^{3}-a^{2}+1}\cdot\frac{a^{3}+2z^{3}}{(a^{3}-z^{3})^{2}},$
(73)

where $a>2^{1/3}$ due to the condition for no self-intersections. Figure 9
plots different shapes created by the mapping for various $a$. Ienaga et al.
examined the intersections of discs and the number of defects. They observed
that when the distances between the discs are sufficiently large, a defect
with a charge of $-1/2$ forms at the center, and $+1/2$ defects form at
positions $w=\hat{r}\exp(2\pi m{\rm i}/3)$, $\hat{r}\in\mathbb{R}$, $m=0,1,2$.
Furthermore, when the three discs intersect, they experimentally found that
the three $+1/2$ defects no longer exist; instead, one defect annihilates with
the central $-1/2$ defect, resulting in two $+1/2$ defects forming along the
$x$-axis or at positions symmetric to the $x$-axis.

The positions of these defects are algebraically determined by assuming
symmetry for the domain. To investigate whether the resulting defects are
stable, the Hessian of the energy at these positions is examined and their
stability is evaluated based on whether the Hessian is positive definite. For
the calculation of the Hessian, the following quantities are precomputed:

$\displaystyle\frac{g^{\prime\prime}}{g^{\prime}}=\frac{6z^{2}}{a^{3}+2z^{3}}+\frac{6z^{2}}{a^{3}-z^{3}},\quad\left(\frac{g^{\prime\prime}}{g^{\prime}}\right)^{\prime}=\frac{12z(a^{3}-z^{3})}{(a^{3}+2z^{3})^{2}}+\frac{6z(2a^{3}+z^{3})}{(a^{3}-z^{3})^{2}}.$
(74)

The distance between discs is given by

$\displaystyle\Delta=\sqrt{3}g(1/a)=\frac{\sqrt{3}a^{2}}{a^{3}-a^{2}+1}.$ (75)

Figure 9: Defect positions in triplets. Black dots in (a) correspond to the
centres of the circles. Red circles are topological defects with a charge of
$+1/2$ and blue squares are defects with a charge of $-1/2$. (a) Case where
two defects exist on the real axis, (b) case where defects are located at
positions symmetric about the $x$-axis, and (c) case where defects with a
charge of $+1/2$ are located at each disc and the one defect with a charge of
$-1/2$ is located at the origin.

#### Two defects (i): two defects located on real axis of predomain

Assume that there exist two defects on the centre line of the quadrature
domain. The coordinates of these defects in the predomain are defined by
$z_{1}=r_{1},\ z_{2}=r_{2}$, $-1<r_{1}<r_{2}<1$. The energy $F_{d}$ is
explicitly given by

$\displaystyle F_{d}$
$\displaystyle=-\frac{1}{4}\left[\log(1-r_{1}^{2})+\log(1-r_{2}^{2})+2\log(r_{2}-r_{1})+2\log(1-r_{1}r_{2})\right]-\frac{3}{4}[\log
g^{\prime}(r_{1})+\log g^{\prime}(r_{2})].$

The derivative of $F_{d}$ with respect to $r_{1}$ and $r_{2}$ gives

$\displaystyle\left\\{\begin{aligned} \frac{\partial F_{d}}{\partial
r_{1}}&=\frac{r_{1}}{2(1-r_{1}^{2})}-\frac{1}{2(r_{2}-r_{1})}+\frac{r_{2}}{2(1-r_{1}r_{2})}-\frac{3}{4}\left[\frac{6r_{1}^{2}}{a^{3}+2r_{1}^{3}}-\frac{6r_{1}^{2}}{r_{1}^{3}-a^{3}}\right]=0,\\\
\frac{\partial F_{d}}{\partial
r_{2}}&=\frac{r_{2}}{2(1-r_{2}^{2})}-\frac{1}{2(r_{1}-r_{2})}+\frac{r_{1}}{2(1-r_{1}r_{2})}-\frac{3}{4}\left[\frac{6r_{2}^{2}}{a^{3}+2r_{2}^{3}}-\frac{6r_{2}^{2}}{r_{2}^{3}-a^{3}}\right]=0.\end{aligned}\right.$
(76)

The solution (76) for $r_{1}$ and $r_{2}$ are obtained by the simplex method
implemented in the standard matlab solvers [27]. The Hessian of the
coordinates $z_{1}=r_{1}$ and $z_{2}=r_{2}$ is calculated as follows:

$\displaystyle H=\begin{pmatrix}J&\mathcal{O}\\\
\mathcal{O}&J\end{pmatrix}\begin{pmatrix}H_{11}&H_{12}\\\
H_{21}&H_{22}\end{pmatrix}\begin{pmatrix}J&\mathcal{O}\\\
\mathcal{O}&J\end{pmatrix}^{\top},$ (77)

where each component is given by

$\displaystyle\mathcal{J}\equiv\begin{pmatrix}1&1\\\ {\rm i}&-{\rm
i}\end{pmatrix},\quad
H_{kl}\equiv\begin{pmatrix}\dfrac{\partial^{2}F_{d}}{\partial z_{k}\partial
z_{l}}&\dfrac{\partial^{2}F_{d}}{\partial
z_{k}\partial\overline{z_{l}}}\vspace{3pt}\\\
\dfrac{\partial^{2}F_{d}}{\partial\overline{z_{k}}\partial
z_{l}}&\dfrac{\partial^{2}F_{d}}{\partial\overline{z_{k}}\partial\overline{z_{l}}}\end{pmatrix}.$
(78)

When the Hessian $H$ is positive definite, the coordinates $z=r_{1}$ and
$z=r_{2}$ are stable.

The black solid line in Figure 10 shows the smallest eigenvalue of the Hessian
$H$. Our numerical simulations show that the defects become unstable when
$\Delta$ becomes larger than approximately $0.809$.

Figure 10: Smallest eigenvalue of Hessian matrix with respect to $\Delta$.
When $\Delta$ is less than $0.809$, two defects on the centre line are stable.
When $\Delta>1.481$, four defects are stable. The formation of two defect
pairs symmetrically positioned along the real axis is shown to be always
stable with respect to $\Delta$.

#### Two defects (ii): two defect pairs symmetrically positioned along real
axis

Assuming symmetry about the real axis, two defect positions in the predomain
are set to $z_{1}=re^{{\rm i}\theta}$, $z_{2}=re^{-{\rm i}\theta}$. The
conditions of $r$ and $\theta$ for the two defect locations are as follows:

$\displaystyle\left\\{\begin{aligned} &\frac{\partial F_{d}}{\partial
z_{1}}=-\frac{1}{4}\left(-\frac{re^{-{\rm
i}\theta}}{1-r^{2}}+\frac{1}{r(e^{{\rm i}\theta}-e^{-{\rm
i}\theta})}+\frac{-re^{-{\rm i}\theta}}{1-r^{2}e^{2{\rm
i}\theta}}\right)-\frac{3}{8}\frac{g^{\prime\prime}(re^{{\rm
i}\theta})}{g^{\prime}(re^{{\rm i}\theta})}=0,\\\ &\frac{\partial
F_{d}}{\partial z_{2}}=-\frac{1}{4}\left(-\frac{re^{{\rm
i}\theta}}{1-r^{2}}+\frac{1}{r(e^{-{\rm i}\theta}-e^{{\rm
i}\theta})}+\frac{-re^{{\rm i}\theta}}{1-r^{2}e^{-2{\rm
i}\theta}}\right)-\frac{3}{8}\frac{g^{\prime\prime}(re^{-{\rm
i}\theta})}{g^{\prime}(re^{-{\rm i}\theta})}=0.\end{aligned}\right.$ (79)

By solving (79) for $r$ and $\theta$ numerically subject to $0<r<1$ and
$0<\theta<\pi$, we can find a unique solution. Using this solution, the
Hessian is calculated in the same manner.

The red dashed line in Figure 10 shows the smallest eigenvalue of the Hessian
matrix $H$ in this case. The smallest eigenvalue for two defects pairs
symmetrically positioned along the real axis is always positive, which means
that this formation is always stable in this triplet model.

#### Four defects

Assume that $q_{0}=-1/2$ exists in the centre of triplets and three
topological defects with charge $q_{k}=+1/2$, $k=1,2,3$, exist at symmetric
positions. The coordinates of these defects are defined by $z_{0}=0$,
$z_{1}=s$, $z_{2}=s\omega$, $z_{3}=s\omega^{2}$, $\omega=e^{2\pi{\rm i}/3}$,
$s>0$. The energy that depends only on $z_{k}$ is given by

$\displaystyle
F_{d}=-\frac{3}{4}\left[3\log\frac{a^{3}+2s^{3}}{(a^{3}-s^{3})^{2}}+\log(1-s^{6})\right]-\frac{9}{4}\log\left(\frac{a^{6}-1}{a^{3}-a^{2}+1}\right)+\frac{5}{4}\log|g^{\prime}(0)|-\frac{3}{4}\log
3.$

The extrema of $F_{d}$ with respect to $s$ is obtained by

$\displaystyle\frac{\partial F_{d}}{\partial
s}=-\frac{3}{4}\left[\frac{18s^{2}}{a^{3}+2s^{3}}+\frac{18s^{2}}{a^{3}-s^{3}}-\frac{6s^{5}}{1-s^{6}}\right]=-\frac{9}{2}s^{2}\left[\frac{3(2a^{3}+s^{3})}{(a^{3}+2s^{3})(a^{3}-s^{3})}-\frac{s^{3}}{1-s^{6}}\right]=0,$

which yields the following $9$-th-order polynomial formula with respect to
$s$:

$\displaystyle p_{3}(t)=t^{3}+7a^{3}t^{2}-(3-a^{6})t-6a^{3}=0,\quad t=s^{3}.$
(80)

Note that $p_{3}(0)=-6a^{3}<0$, $p_{3}(1)=a^{6}+a^{3}-2>0$, and

$\displaystyle p_{3}^{\prime}(t)=3t^{2}+14a^{3}t-3+a^{6},\quad
p_{3}^{\prime\prime}(t)=6t+14a^{3}>0,$ (81)

which means that $p_{3}(t)=0$ has a root in the range $0<t<1$.

The blue dotted line in Figure 10 shows the smallest eigenvalue of the Hessian
matrices for the three cases. When $\Delta>1.481$, the smallest eigenvalue is
positive, which means the formation of three topological defects becomes
stable when $\Delta$ becomes larger than approximately $1.481$.

##  5 Discussion and conclusion

This paper provides analytical formulas for Frank free energy in the presence
of topological defects based on analogies to vortex dynamics. The energy
representation was obtained using the approach of Hamiltonians for vortex
dynamics. The analytical formulas derived in this paper were used to obtain
the positions of defects in doublet and triplet regions. The proposed formulas
are explicit and give a more accurate condition for $\Delta$ when the
transition of formation of topological defects occurs.

Our analytical results closely match the experimental results reported by
Ienaga et al., particularly in the case of doublets [10]. However, the defect
distances $D_{++}$ obtained using the actual experiments conducted by Ienaga
et al. were found to be shorter than those obtained from our analytical
methods. This suggests that the influence of cells must be considered, as we
assumed cells to behave as nematic liquid crystals.

Although the present paper focuses on simply connected domains, it is
important to note that Miyazako and Sakajo derived an analytical formula for
cell alignments in doubly connected domains [17]. The derivation of the
explicit formulas for energy in doubly connected domains is an ongoing work.
Using the analytical energy formulas and the theory of quadrature domains in
doubly connected domains described in [18], it might become possible to study
ferromagnetic and antiferromagnetic order in bacterial vortex lattices in more
complicated geometries, such as three merging circles with an inner holes.

Data Accessibility: This article has no additional data.

Authors’ Contributions: All authors contributed to the content of this paper.

Competing Interests: We declare no competing interests.

Funding: The first author is supported by JSPS KAKENHI Grant No. JP24KJ0041.
The second author is partially supported by JSPS KAKENHI, Grant No. JP23H00086
and JP23H04406.

Acknowledgements: The first author thanks Prof. Darren G. Crowdy at Imperial
College London for fruitful discussions.

##  Appendix A Hessian

We derive the Hessian with complex coordinates $z_{k}$ and $\overline{z_{k}}$.
Because the derivatives with respect to $z_{k}$ and $\overline{z_{k}}$ are
given by

$\displaystyle\frac{\partial}{\partial
z_{k}}=\frac{1}{2}\left(\frac{\partial}{\partial x_{k}}-{\rm
i}\frac{\partial}{\partial y_{k}}\right),\
\frac{\partial}{\partial\overline{z_{k}}}=\frac{1}{2}\left(\frac{\partial}{\partial
x_{k}}+{\rm i}\frac{\partial}{\partial y_{k}}\right),$ (82)

we have

$\displaystyle\frac{\partial^{2}}{\partial x_{k}\partial
x_{l}}=\left(\frac{\partial}{\partial
z_{k}}+\frac{\partial}{\partial\overline{z_{k}}}\right)\left(\frac{\partial}{\partial
z_{l}}+\frac{\partial}{\partial\overline{z_{l}}}\right)=\frac{\partial^{2}}{\partial
z_{k}\partial z_{l}}+\frac{\partial^{2}}{\partial\overline{z_{k}}\partial
z_{l}}+\frac{\partial^{2}}{\partial
z_{k}\partial\overline{z_{l}}}+\frac{\partial^{2}}{\partial\overline{z_{k}}\partial\overline{z_{l}}},$
(83) $\displaystyle\frac{\partial^{2}}{\partial x_{k}\partial y_{l}}={\rm
i}\left(\frac{\partial}{\partial
z_{k}}+\frac{\partial}{\partial\overline{z_{k}}}\right)\left(\frac{\partial}{\partial
z_{l}}-\frac{\partial}{\partial\overline{z_{l}}}\right)={\rm
i}\left(\frac{\partial^{2}}{\partial z_{k}\partial
z_{l}}+\frac{\partial^{2}}{\partial\overline{z_{k}}\partial
z_{l}}-\frac{\partial^{2}}{\partial
z_{k}\partial\overline{z_{l}}}-\frac{\partial^{2}}{\partial\overline{z_{k}}\partial\overline{z_{l}}}\right),$
(84) $\displaystyle\frac{\partial^{2}}{\partial y_{k}\partial
y_{l}}=-\left(\frac{\partial}{\partial
z_{k}}-\frac{\partial}{\partial\overline{z_{k}}}\right)\left(\frac{\partial}{\partial
z_{l}}-\frac{\partial}{\partial\overline{z_{l}}}\right)=-\frac{\partial^{2}}{\partial
z_{k}\partial z_{l}}+\frac{\partial^{2}}{\partial\overline{z_{k}}\partial
z_{l}}+\frac{\partial^{2}}{\partial
z_{k}\partial\overline{z_{l}}}-\frac{\partial^{2}}{\partial\overline{z_{k}}\partial\overline{z_{l}}}.$
(85)

Using this transformation, the Hessian with respect to $x_{k}$ and $y_{l}$,
$k,l=1,\ldots,N$, is given by

$\displaystyle\begin{pmatrix}\dfrac{\partial^{2}}{\partial x_{k}\partial
x_{l}}&\dfrac{\partial^{2}}{\partial x_{k}\partial y_{l}}\\\
\dfrac{\partial^{2}}{\partial y_{k}\partial
x_{l}}&\dfrac{\partial^{2}}{\partial y_{k}\partial
y_{l}}\end{pmatrix}=\mathcal{J}\begin{pmatrix}\dfrac{\partial^{2}}{\partial
z_{k}\partial z_{l}}&\dfrac{\partial^{2}}{\partial
z_{k}\partial\overline{z_{l}}}\\\
\dfrac{\partial^{2}}{\partial\overline{z_{k}}\partial
z_{l}}&\dfrac{\partial^{2}}{\partial\overline{z_{k}}\partial\overline{z_{l}}}\end{pmatrix}\mathcal{J}^{\top},\quad\mathcal{J}\equiv\begin{pmatrix}1&1\\\
{\rm i}&-{\rm i}\end{pmatrix}.$ (86)

## References

  * [1] J. Zhao, U. Gulan, T. Horie, N. Ohmura, J. Han, C. Yang, J. Kong, S. Wang, and B. B. Xu.  Advances in biological liquid crystals.  Small, 15(18):1900019, 2019.
  * [2] A. Doostmohammadi and B. Ladoux.  Physics of liquid crystals in cell biology.  Trends in cell biology, 32(2):140–150, 2022.
  * [3] S. J. Woltman, G. D. Jay, and G. P. Crawford.  Liquid-crystal materials find a new order in biomedical applications.  Nat. Mater., 6(12):929–938, December 2007.
  * [4] M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha.  Hydrodynamics of soft active matter.  RMP, 85(3):1143, 2013.
  * [5] P.-G. De Gennes and J. Prost.  The physics of liquid crystals.  Number 83. OUP, 1993.
  * [6] I. S. Aranson.  Topological defects in active liquid crystals.  Phys.-Uspekhi, 62(9):892, 2019.
  * [7] Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun, and K. Keren.  Topological defects in the nematic order of actin fibres as organization centres of hydra morphogenesis.  Nat. Phys., 17(2):251–259, 2021.
  * [8] K. Kawaguchi, R. Kageyama, and M. Sano.  Topological defects control collective dynamics in neural progenitor cell cultures.  Nature, 545(7654):327–331, April 2017.
  * [9] T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans, and B. Ladoux.  Topological defects in epithelia govern cell death and extrusion.  Nature, 544(7649):212–216, 2017.
  * [10] R. Ienaga, K. Beppu, and Y. T. Maeda.  Geometric confinement guides topological defect pairings and emergent flow in nematic cell populations.  Soft Matter, 19(26):5016–5028, July 2023.
  * [11] G. Duclos, C. Erlenkämper, J.-F. Joanny, and P. Silberzan.  Topological defects in confined populations of spindle-shaped cells.  Nat. Phys., 13(1):58–62, September 2016.
  * [12] O. S. Tarnavskyy, M. F. Ledney, and A. I. Lesiuk.  Generalised technique for calculation of plane director profiles in bounded nematic liquid crystals.  Liq. Cryst., 45(5):641–648, 2018.
  * [13] T. G. J. Chandler and S. E. Spagnolie.  A nematic liquid crystal with an immersed body: equilibrium, stress and paradox.  J. Fluid Mech., 967:A19, 2023.
  * [14] T. G. J. Chandler and S. E. Spagnolie.  Exact and approximate solutions for elastic interactions in a nematic liquid crystal.  arXiv preprint arXiv:2311.17708, 2023.
  * [15] X. Tang and J. V. Selinger.  Orientation of topological defects in 2D nematic liquid crystals.  Soft Matter, 13(32):5481–5490, August 2017.
  * [16] H. Miyazako and T. Nara.  Explicit calculation method for cell alignment in non-circular geometries.  R. Soc. Open Sci., 9(1):211663, January 2022.
  * [17] H. Miyazako and T. Sakajo.  Defect pairs in nematic cell alignment on doubly connected domains.  Proc. Roy. Soc. A (accepted), 2024.
  * [18] D. G. Crowdy.  Solving problems in multiply connected domains.  SIAM, 2020.
  * [19] C. C. Lin.  On the motion of vortices in two dimensions: I. existence of the kirchhoff-routh function.  PNAS, 27(12):570–575, 1941.
  * [20] B. Gustafsson and H. S. Shapiro.  What is a quadrature domain? Springer, 2005.
  * [21] P. G. Saffman.  Vortex dynamics.  CUP, 1995.
  * [22] P. K. Newton and MF Platzer.  $n$-vortex problem: Analytical techniques.  Appl. Mech. Rev., 55(1):B15–B16, 2002.
  * [23] D. G. Crowdy and J. Marshall.  Analytical formulae for the kirchhoff–routh path function in multiply connected domains.  Proc. Roy. Soc. A, 461(2060):2477–2501, 2005.
  * [24] N. M. Khalifa and H. E. Taha.  Vortex dynamics: A variational approach using the principle of least action.  Phys. Rev. Fluids, 9(3):034701, 2024.
  * [25] M. Flucher.  Vortex motion in two dimensional hydrodynamics.  Variational Problems with Concentration, pages 131–149, 1999.
  * [26] H. Wioland, F. G. Woodhouse, J. Dunkel, and R. E. Goldstein.  Ferromagnetic and antiferromagnetic order in bacterial vortex lattices.  Nat. Phys., 12(4):341–345, 2016.
  * [27] J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright.  Convergence properties of the nelder–mead simplex method in low dimensions.  SIOPT, 9(1):112–147, 1998.