e1e-mail<EMAIL_ADDRESS>e2e-mail:
<EMAIL_ADDRESS>
11institutetext: Institut Curie, Université PSL, Sorbonne Université, CNRS
UMR168, Laboratoire Physico Chimie Curie, 75005 Paris, France 22institutetext:
Collège de France, 75005 Paris, France

Topological defects are at the root of the large-scale organization of liquid
crystals. In two-dimensional active nematics, two classes of topological
defects of charges $\pm 1/2$ are known to play a major role due to active
stresses. Despite this importance, few analytical results have been obtained
on the flow-field and active-stress patterns around active topological
defects. Using the generic hydrodynamic theory of active systems, we
investigate the flow and stress patterns around these topological defects in
unbounded, two-dimensional active nematics. Under generic assumptions, we
derive analytically the spontaneous velocity and stall force of self-advected
defects in the presence of both shear and rotational viscosities. Applying our
formalism to the dynamics of monolayers of elongated cells at confluence, we
show that the non-conservation of cell number generically increases the self-
advection velocity and could provide an explanation for their observed role in
cellular extrusion and multilayering. We finally investigate numerically the
influence of the Ericksen stress. Our work paves the way to a generic study of
the role of topological defects in active nematics, and in particular in
monolayers of elongated cells.

# Spontaneous flow created by active topological defects

Louis Brézinaddr1,addr2  Thomas Rislere2,addr1  Jean-François
Joannye1,addr1,addr2

(February 1, 2022)

Topological defects are singularities in the topology of the order-parameter
field in a phase of broken continuous symmetry MerminTopological1979 ;
kleman1983 . They were first observed in nematic liquid crystals as early as
1904 by Lehman lehmann1904  and sparked theoretical interest later by Frank
and Kleman for their role in the mechanical and optical properties of liquid
crystals kleman1983 ; kleman1989 . In two-dimensional nematic liquids,
topological defects are point-like singularities of the director-field
orientation. They are classified by the defect strength or topological charge
$S$, which is equal to the number of full rotations of the liquid-crystal
director on a closed path of total angle $2\pi$ around the defect core
OswaldNematic2005 ; HarthTopological2020 . Due to the nematic symmetry, this
topological charge is either an integer or a half-integer. The regions that
surround a topological defect present strong gradients of the director-field
orientation, with high energetic costs. Because of their relatively lower
energy, the defects that are mostly observed in passive nematic liquid
crystals have topological charge $S=\pm 1/2$, which is the lowest possible
charge in absolute value.

Largely driven by biophysical applications, much attention has been attracted
more recently to active matter JulicherActive2007 ; marchetti2013 . In active
nematics, numerical simulations giomi2013 ; thampi2013  as well as analytical
analyses pismen2013  have revealed a very rich dynamics, dominated by the
creation, motion, and annihilation of defects. Topological defects in living
systems have been first observed in 1968 by Elsdale in monolayers of
fibroblast cells elsdale1968  and, since then, the topological charge of the
defects has helped to identify the nematic or polar nature of cellular tissues
kemkemer2000 ; duclos2017a ; saw2018 ; blanch-mercader2020 ; blanch-
mercader2020a . Recent works looked at the influence of topological defects on
the onset of three-dimensional morphogenesis maroudas-sacks2021 ;
metselaar2019 . Some material properties can also be extracted from the
observation of the orientation of the director field around topological
defects blanch-mercader2020 ; blanch-mercader2020a .

The properties of topological defects in active nematic liquid crystals are
very different from the properties of their equilibrium counterparts. The
active stress in an active nematic is proportional to the nematic
orientational tensor. As regions surrounding a topological defect have
pronounced nematic orientation gradients, they correspond to regions of strong
active force density, which induce local flows around the defect. The
structure of this flow depends on the topological charge of the defect, but
defects in general are associated with vortices of the flow field giomi2013 ;
thampi2013 . If the director field around the defect breaks a left-right
symmetry along a particular direction, and in the absence of any external
force, the flow created by active stress gradients drags the defect along that
direction with a finite velocity, limited by viscous dissipation. Defects of
topological charge $S=-1/2$ have a three-fold symmetry. Consequently, they
have no preferred direction of motion and passively diffuse. On the contrary,
defects of topological charge $S=+1/2$ break a ‘head-tail’ symmetry. In active
nematics, they generically move along their head-tail axis, with a direction
dictated by their contractile (motion toward the tail) or extensile (motion
toward the head) nature.

Experimental observations on nematic cell monolayers show that topological
defects are preferential sites of cell extrusion saw2017 ; kawaguchi2017  and
multilayering sarkar2021 . Similar observations have been made in bacterial
monolayers copenhagen2021 . In a recent work by T. Sarkar and colleagues
sarkar2021  on the spontaneous organization of cell monolayers on a solid
substrate, the authors measured the flow patterns around topological defects
of charges $\pm 1/2$. Some of the $+1/2$ topological defects were found to be
pinned at specific locations on the substrate, while others had a spontaneous
motion as expected in active nematics.

In the current paper, we study theoretically the flows generated by active,
two-dimensional topological defects of charges $S=\pm 1/2$. One key ingredient
of the hydrodynamics of nematic liquid crystals is the backflow, the
hydrodynamic motion induced by the rotation of the director field degennes-
prost . Analytical results on topological defects have been obtained by
neglecting this backflow giomi2014 ; Ronning2022 , while numerical studies
have considered the coupling between the flow and the director fields in
passive toth2002  as well as active giomi2013 ; giomi2014  nematic defects. We
provide here analytical solutions that account explicitly for these effects.
We compute the force necessary to stall an active topological defect of charge
$S=+1/2$ and investigate the link between the two-dimensional flow patterns
and the trigger of multilayering. We limit our study to a single, isolated
defect in an infinite plane. Thereby, we rely on analytical solutions of the
active hydrodynamic equations that depend on the bulk properties of the
cellular tissue only. Performing a perturbation analysis around the passive
defect orientation allows us to obtain analytical solutions for weakly active
systems. We also estimate the size of the spatial domain over which this
hypothesis is valid.

The paper is organized as follows: section 1 introduces the hydrodynamic
equations of active nematics by giving the constitutive equations and
computing the active force density generated by a topological defect of charge
$S=\pm 1/2$. We then study in section 2 the dynamics of the defect, assuming
that the nematic orientation is unperturbed compared to that of a passive
defect of the same charge. We focus on $+1/2$ defects, but computations for
$-1/2$ defects are detailed in E. We compute the flow field and the self-
advection velocity of the defect in section 2.1. In section 2.2, we compute
the force necessary to pin a topological defect to a specific location on the
substrate. In section 2.3, we study the effect of inhomogeneous cell divisions
and cell deaths or extrusions on the self-advection velocity of the defect, in
the case where the rate of cell production is linearly coupled to the
differential pressure field. In section 3, we address the role of the
rotational viscosity, which induces feedback of the flow on the nematic
orientation. We first ignore the elastic Ericksen stress and set the flow-
alignment parameter to zero in section 3.1. This allows for the computation of
the flow and the self-advection velocity to linear order in activity
analytically. Within the same hypothesis, the crucial role of the backflow on
the stall force is detailed in section 3.2. The first-order correction to the
director orientation as well as the spatial domain of validity of our
calculations are discussed in section 3.3. Finally, we investigate numerically
the role of the Ericksen stress on the self-advection velocity in section 3.4.

##  1 Hydrodynamic equations

We consider an active nematic deep in the nematic phase where the order
parameter has a constant modulus. To describe its dynamics, we make use of the
hydrodynamic theory of active nematic gels kruse2005 ; joanny2007 ;
marchetti2013 . This very general framework is based on symmetries and
conservation laws joanny2009  and has been shown to adequately describe
monolayers of elongated cells that self-organize into a nematic phase
duclos2014 ; duclos2018b ; saw2018 ; alert2020 . More specifically, we aim at
describing the dynamics of nematic topological defects as observed
experimentally in monolayers of myoblast cells sarkar2021 . Because we are
interested in monolayers of cells at confluence, we consider the active
nematic as a one-constituent dense phase, with a fixed cell density.

Within the active hydrodynamic framework, the monolayer is described by a
velocity field $\mathbf{v}$ and a director field $\mathbf{p}$. The two-
dimensional Frank-Oseen free energy $F$ associated to gradients of the
director-field orientation reads, in the one-constant approximation where the
splay and bend constants are equal degennes-prost :

$F=\int{\mathop{}\mathopen{}\mathrm{d}x\>\mathop{}\mathopen{}\mathrm{d}y\left[\frac{K}{2}\textrm{Tr}\left[\left(\nabla\otimes\mathbf{p}\right)\cdot\left(\nabla\otimes\mathbf{p}\right)^{\rm
T}\right]-\frac{1}{2}h_{\parallel}^{0}\,\mathbf{p}^{2}\right]}\,,$ (1)

where $\otimes$ denotes the tensorial product, and $h_{\parallel}^{0}$ is a
Lagrange multiplier to ensure that the director is a unit vector.

The evolution of the velocity and director fields is described by two coupled
vectorial equations. The first equation is the force-balance equation, which
reads, at steady state:

$\nabla\cdot\bm{\sigma}=\mathbf{0}\,.$ (2)

Inertia has been neglected since the Reynolds number for typical cellular
systems is much smaller than one.

The total stress tensor is further split into a passive $\bm{\sigma}^{\rm p}$
and an active $\bm{\sigma}^{\rm a}$ contribution. The passive contribution is
given by the hydrodynamic theory of passive nematics degennes-prost ;
martin1972 ; forster1971  as:

$\displaystyle\bm{\sigma}^{\rm p}=$ $\displaystyle
2\eta\,\tilde{\mathbf{u}}+\frac{\nu}{2}\left(\mathbf{h}\otimes\mathbf{p}+\mathbf{p}\otimes\mathbf{h}-(\mathbf{p}\cdot\mathbf{h})\,\mathbf{1}\right)$
$\displaystyle+\frac{1}{2}\left(\mathbf{h}\otimes\mathbf{p}-\mathbf{p}\otimes\mathbf{h}\right)+\tilde{\bm{\sigma}}^{\rm
E}-P\,\mathbf{1}\,.$ (3)

Here, $\eta$ is the shear viscosity; $\tilde{\mathbf{u}}$ is the traceless
part of the symmetric part of the velocity-gradient tensor
$\mathbf{u}=(\nabla\otimes\mathbf{v}+(\nabla\otimes\mathbf{v})^{\rm T})/2$;
$\nu$ is the so-called flow-alignment parameter; $\mathbf{h}=-\delta
F/\delta\mathbf{p}$ is the orientational field; $\tilde{\bm{\sigma}}^{\rm E}$
is the traceless part of the Ericksen stress tensor $\bm{\sigma}^{\rm E}$,
which generalizes thermodynamic pressure for anisotropic systems degennes-
prost ; chaikin1995 ; and $P$ is the pressure, which is a Lagrange multiplier
that ensures the incompressibility condition $\nabla\cdot\mathbf{v}=0$, with
$\mathbf{1}$ the identity tensor. The elastic force density associated with
the Ericksen stress is obtained from the Gibbs-Duhem relation julicher2018 ;
joanny2007 ; degennes-prost

$\nabla\cdot\bm{\sigma}^{\rm E}=-(\nabla\otimes\mathbf{p})\cdot\mathbf{h}\,.$
(4)

The active contribution is linked to the existence of a local cell alignment
along the direction of the director field $\mathbf{p}$. It reads

$\bm{\sigma}^{\rm
a}=-\zeta\Delta\mu\left(\mathbf{p}\otimes\mathbf{p}-\frac{\mathbf{1}}{2}\right)\,,$
(5)

where $\zeta\Delta\mu$ is a scalar quantity that measures the activity of the
system. It is the product of a difference in chemical potential $\Delta\mu$,
which is positive, and an Onsager coefficient $\zeta$, which can have either a
positive or a negative value. A positive value of $\zeta$ corresponds to an
extensile active stress, in which the cells push along their long axis. The
stress $\zeta\Delta\mu$ is negative in a contractile system where the cells
pull along their long axis.

The second equation describes the evolution of the director field
$\mathbf{p}$:

$\frac{\mathop{}\mathopen{}\mathrm{D}\mathbf{p}}{\mathop{}\mathopen{}\mathrm{D}t}=\frac{1}{\gamma}\mathbf{h}-\nu\,\mathbf{u}\cdot\mathbf{p}\,.$
(6)

Here,
$\mathop{}\mathopen{}\mathrm{D}\mathbf{p}/\mathop{}\mathopen{}\mathrm{D}t=\partial_{t}\mathbf{p}+(\mathbf{v}\cdot\nabla)\mathbf{p}+\bm{\omega}\cdot\mathbf{p}$
is the co-moving co-rotational derivative of the director, where
$\bm{\omega}=(\nabla\otimes\mathbf{v}-(\nabla\otimes\mathbf{v})^{\rm T})/2$ is
the vorticity tensor, and $\gamma$ is the rotational viscosity. The
orientational field $\mathbf{h}=-\delta F/\delta\mathbf{p}$ is associated to
changes of the free energy with respect to the director $\mathbf{p}$. It is
parallel to $\mathbf{p}$ in a non-flowing steady state. It is convenient to
introduce the components of the orientational field, parallel and
perpendicular to $\mathbf{p}$, $h_{\parallel}$ and $h_{\perp}$, respectively.
The component $h_{\perp}$ controls the orientation of the director
$\mathbf{p}$, whereas the component $h_{\parallel}$ controls the modulus of
$\mathbf{p}$ kruse2005 .

In the following, we define the orientation angle $\varphi$ of the director
$\mathbf{p}$ as
$\mathbf{p}=\cos{\varphi}\,\mathbf{e}_{x}+\sin{\varphi}\,\mathbf{e}_{y}$ in
the fixed reference frame $\\{\mathbf{e}_{x},\mathbf{e}_{y}\\}$. For passive
nematics, the equilibrium configurations correspond to minima of the Frank-
Oseen free energy given in eq. (1), which leads to $\Delta\varphi=0$. This
equation has the solution $\varphi(r,\theta)=S\,\theta$ for a topological
defect of charge $S$, using polar coordinates $(r,\theta)$ centered at the
defect singularity. Defects with the lowest absolute topological charge $S=\pm
1/2$ are the most stable. They are shown in fig. 1.

Figure 1: Passive orientation of the director field (solid lines) and
resulting active force density (arrows) in the contractile case
$\zeta\Delta\mu<0$, for $+1/2$ (a) and $-1/2$ (b) topological defects. The
origin of the polar angle $\theta$ is set by $\mathbf{e}_{x}$.

In free monolayers of nematic tissue, they are the only observed defects.

For a director field corresponding to a passive defect, using the polar
coordinates $(r,\theta)$ in the Cartesian basis of vectors
$\\{\mathbf{e}_{x},\mathbf{e}_{y}\\}$, the active force density reads:

$\displaystyle\mathbf{f}^{a}=\nabla\cdot\bm{\sigma}^{\rm
a}=-\frac{\zeta\Delta\mu}{2r}\begin{cases}\mathbf{e}_{x}\ ,&S=+1/2\\\
-\cos{2\theta}\ \mathbf{e}_{x}+\sin 2\theta\ \mathbf{e}_{y}\
,&S=-1/2\end{cases}$ (7)

As represented in fig. 1, the active force density for a $+1/2$ topological
defect is along the axis of symmetry. This drives a self-advection of the
defect ‘from head to tail’ in the contractile case and ‘from tail to head’ in
the extensile case. For $-1/2$ defects, because of the three-fold symmetry of
the director field, the net total force on a disk of radius $R$ vanishes, such
that the default, even tough active, can only diffuse.

##  2 Topological defects in the limit of vanishing rotational viscosity

Throughout this section, we consider the case of vanishing rotational
viscosity. In this case, the molecular field vanishes to satisfy eq. (6), so
that $\delta F/\delta\mathbf{p}=\mathbf{0}$ and the director orientation
satisfies the equilibrium condition for passive nematics $\Delta\varphi=0$,
illustrated in fig. 1. The force density associated to the Ericksen stress as
given by eq. (4) vanishes. The Ericksen stress can then be ignored and the
passive stress as given by eq. (1) reduces to
${\bm{\sigma}}^{p}=2\eta\,\tilde{\mathbf{u}}-P\,\mathbf{1}$.

###  2.1 Flow field and self-advection velocity of a $+1/2$ defect

Under this assumption, we now compute the velocity field of an isolated $+1/2$
topological defect in an infinite domain. We write the equations of motion in
the reference frame of the defect. Force balance reads

$\eta\Delta\mathbf{v}-\mathbf{\nabla}P-\xi(\mathbf{v}+\mathbf{v_{0}})-\frac{\zeta\Delta\mu}{2r}\mathbf{e}_{x}=0\
,$ (8)

where $\xi$ is the friction coefficient of the cell layer with the underlying
substrate and $\mathbf{v}_{0}$ the self-advection velocity of the defect with
respect to the substrate. This bulk equation is supplemented by two boundary
conditions on the velocity field: it vanishes at the core of the defect in the
reference frame of the defect $\mathbf{v}(0,\theta)=\mathbf{0}$, and it
vanishes at infinity in the reference frame of the substrate
$\mathbf{v}(\infty,\theta)=-\mathbf{v}_{0}$.

Considering the incompressibility condition $\nabla\cdot\mathbf{v}=0$, eq. (8)
is solved by introducing the stream function $\psi(r,\theta)$, defined by
$\mathbf{v}=\mathbf{\nabla}\times(\psi\ \mathbf{e}_{z})$, where
$\mathbf{e}_{z}=\mathbf{e}_{x}\times\mathbf{e}_{y}$. We further
adimensionalize the equations using the length scale $L=\sqrt{\eta/\xi}$ and
the time scale $\tau=2\eta/(|\zeta|\Delta\mu)$. In the following, all
dimensionless quantities are denoted with a tilde. Taking the curl of eq. (8),
the dimensionless stream function $\tilde{\psi}=(\tau/L^{2})\psi$ satisfies

$\Delta\left[\Delta\tilde{\psi}-\tilde{\psi}\right]=-s\frac{\sin\theta}{\tilde{r}^{2}}\
,$ (9)

where $\tilde{r}=r/L$, $s=\operatorname{sign}(\zeta)$ and $\Delta$ denotes the
Laplace operator with respect to the reduced variable $\tilde{r}$.

A $+1/2$ topological defect is symmetric with respect to the $x$-axis.
Therefore, the radial component of the velocity field is an even function of
the polar angle $\theta$, and the stream function an odd function of that
variable. Given that eq. (9) is linear in $\tilde{\psi}$ with a forcing term
in $\sin{\theta}$, only the Fourier mode $n=1$ contributes non trivially to
the solution, and the stream function can be written as
$\tilde{\psi}(\tilde{r},\theta)=\tilde{\psi}(\tilde{r})\sin\theta$. The
dimensionless velocity field has then polar components
$\tilde{v}_{r}(\tilde{r},\theta)=\tilde{v}_{r}(\tilde{r})\cos\theta$ and
$\tilde{v}_{\theta}(\tilde{r},\theta)=\tilde{v}_{\theta}(\tilde{r})\sin\theta$,
which can be derived directly from the generic solution of eq. (9) for
$\tilde{\psi}$, as detailed in B.1. Imposing that the solution is non-
divergent at infinity and vanishing in $\mathbf{r}=\mathbf{0}$, the solution
for $\tilde{v}_{r}$ and $\tilde{v}_{\theta}$ read:

$\displaystyle\tilde{v}_{r}(\tilde{r})=$
$\displaystyle-\frac{s}{\tilde{r}}\left\\{I_{1}(\tilde{r})\int_{\tilde{r}}^{+\infty}{K_{1}(u)(-\frac{\pi}{4}u^{2}+u)\>\mathop{}\mathopen{}\mathrm{d}u}\right.$
$\displaystyle\qquad+\left.K_{1}(\tilde{r})\int_{0}^{\tilde{r}}{I_{1}(u)(-\frac{\pi}{4}u^{2}+u)\>\mathop{}\mathopen{}\mathrm{d}u}\right\\}$
(10) $\displaystyle\tilde{v}_{\theta}(\tilde{r})=$ $\displaystyle
s\left\\{\left(I_{0}(\tilde{r})-\frac{I_{1}(\tilde{r})}{\tilde{r}}\right)\int_{\tilde{r}}^{+\infty}{K_{1}(u)(-\frac{\pi}{4}u^{2}+u)\>\mathop{}\mathopen{}\mathrm{d}u}\right.$
$\displaystyle\left.-\left(K_{0}(\tilde{r})+\frac{K_{1}(\tilde{r})}{\tilde{r}}\right)\int_{0}^{\tilde{r}}{I_{1}(u)(-\frac{\pi}{4}u^{2}+u)\>\mathop{}\mathopen{}\mathrm{d}u}\right\\}\
,$ (11)

where $I_{1}(u)$ and $K_{1}(u)$ are the modified Bessel functions of the first
and second kind, respectively abramowitz1964  (see A). The self-advection
velocity can then be computed. Taking the limits of eq. (2.1) at infinity
leads to

$\mathbf{v}_{0}=-\frac{\pi}{4}\times\frac{\zeta\Delta\mu}{2\sqrt{\xi\eta}}\
\mathbf{e}_{x}\ ,$ (12)

in physical units. This result is similar to that of ref. Ronning2022 . A
representation of this velocity field as well as of its amplitude as a
function of $\tilde{r}$ is shown in fig. 2a,c.

Figure 2: a,b Flow structure around active $+1/2$ (a) and $-1/2$ (b)
topological defects. Thin, solid black lines represent the orientation of the
director field. Thick, arrowed blue lines are the directed streamlines of the
velocity field in the reference frame of the substrate for a contractile
nematic. Vortices around $+1/2$ defects have a rotational direction that
depends on the sign of $\zeta$. The three-fold symmetry for $-1/2$ defects
does not create self-advection contrary to $+1/2$ defects. c,d Dimensionless
azimuthal (orange) and radial (blue) velocities as a function of the
dimensionless distance to the core $\tilde{r}$ in the reference frame of the
defect. For $+1/2$ defects (c), this velocity vanishes at the core and tends
to a finite limit at infinity. For $-1/2$ defects (d), this velocity vanishes
both at the core and at infinity.

Corresponding expressions for $-1/2$ defects are given in E.1 and plotted in
fig. 2b,d. The solution (2.1) depends only on intrinsic properties of the
system and can therefore be directly compared with experimental measurements
to determine intrinsic parameters in experimental systems.

Because of hydrodynamic screening, the velocity field decays as $1/r$ at
infinity in the substrate reference frame, as does the active-force density.
The total active force diverges linearly with system size but is balanced by
the total friction force. These divergences exist here because we consider a
single defect in an unbounded space. For a pair of defects of opposite charges
$+1/2$ and $-1/2$, for example, the active-force density and the velocity
field decrease as $1/r^{2}$. The pressure field is obtained from the force-
balance eq. (8):

$P(r,\theta)=-\frac{\zeta\Delta\mu}{2}\cos\theta+P_{0}\ ,$ (13)

where $P_{0}$ is a reference pressure. The active stress creates a
differential pressure between the head and the tail of the defect. In the
contractile case ($\zeta\Delta\mu<0$), the pressure is higher at the tail of
the defect, and it is the opposite in the extensile case.

###  2.2 Stall force of a $+1/2$ defect

Sarkar _et al._ observed two classes of $+1/2$ topological defects in
confluent monolayers of C2C12 mouse myoblasts sarkar2021 : motile topological
defects with a flow characteristic of self-advected defects, as presented in
sec. 2.1, and stalled, non-motile defects. These defects where not immobilized
by jamming, as it can occur in other cellular systems garcia2015 , as
significant flows were observed around them. In addition, a one-to-one
correspondence was observed between non-motile defects and the initiation of a
second layer of cells due to extrusion of cells from the first layer, a
process referred to as multilayering. Immotile $+1/2$ defects in an active
monolayer can only exist in the presence of an external force that balances
the total active force, other than the friction drag that structurally
vanishes for a vanishing cell-velocity field.

We determine in this section the local external force applied on the defect
core necessary to pin an active, $+1/2$ defect. We introduce a core region of
the defect of finite size $a$ of the order of the correlation length of the
nematic order degennes-prost , which we consider as a region where there is no
nematic order. We therefore model the monolayer as divided into two regions:
perfect nematic order with $\mathbf{p}^{2}=1$ as in section 2.1 for $r>a$, and
isotropic with $\mathbf{p}^{2}=0$ for $r<a$. This assumption allows for an
analytical treatment and is justified when looking at an isolated defect over
large distances compared to the characteristic length $L$: since the total
active force diverges with system size, the active force that would result
from the weak nematic order inside the finite core region is negligible
compared to the active force from the bulk region in a large system. Indeed,
recent work show that, for isolated defects, the self-advection velocity given
by eq. (12) is valid even when considering the effect of the weak nematic
order inside the core Ronning2022 .

Let us call $\mathbf{f}=-f\mathbf{e}_{x}$ the external force applied to the
core region necessary to stall the defect. We hypothesize further that, in the
core region $r<a$, the external force-density has the same mathematical form
as the active force in eq. (8) to allow for a similar treatment. Force balance
then reads:

$\displaystyle\eta\Delta\mathbf{v}-\mathbf{\nabla}P-\xi\mathbf{v}-\frac{\zeta\Delta\mu}{2r}\mathbf{e}_{x}$
$\displaystyle=0\qquad r>a$ (14)
$\displaystyle\eta\Delta\mathbf{v}-\mathbf{\nabla}P-\xi\mathbf{v}-\frac{f}{2\pi
ar}\mathbf{e}_{x}$ $\displaystyle=0\qquad r<a\ .$ (15)

At the boundary $r=a$ between these two regions, the global flow and stress
fields are continuous. This boundary condition leads to an analytical
treatment under the hypothesis that the radius $a$ of the core region is much
smaller than the hydrodynamic screening length $L=\sqrt{\eta/\xi}$, as
detailed in C.1. To leading order in $a/L$, we obtain

$f=\sqrt{\frac{\eta}{\xi}}\frac{\pi^{2}\zeta\Delta\mu}{\log\left(a\sqrt{\xi/\eta}\right)}\
.$ (16)

A simple way to understand the result of eq. (16) is to consider the defect
velocity in an active nematic submitted to a force localized at the core of
the defect. As the hydrodynamic equations are linear, the velocity is the sum
of two contributions: an active contribution $\mathbf{v}_{0}$ proportional to
the active stress $\zeta\Delta\mu$ and a passive contribution proportional to
the applied force $\mathbf{f}$:

$\mathbf{v}_{\rm defect}=\mathbf{v}_{0}+\chi\mathbf{f}\ ,$ (17)

where $\chi$ is the two-dimensional mobility of the defect. To linear order,
this mobility is that of a defect in a passive nematic. A naive guess is given
by the Saffman-Delbrück result for the mobility of a disk of radius $a$ in a
two-dimensional fluid of shear viscosity $\eta$: $1/\chi\sim\eta/\log(L/a)$,
where the large- and short-scale cutoffs are the screening length $L$ and the
defect core size $a$, respectively. Using the active self-advection velocity
given by eq. (12) leads to eq. (16), up to a numerical prefactor.

###  2.3 Effect of cell division/extrusion on the self-advection of a $+1/2$
defect

We now look at the role of the non-conservation of the cell number on the
active flow created by a topological defect. This non-conservation stems from
cell-division and cell-death processes, as well as from cellular extrusion
from the monolayer. We introduce an effective cell-proliferation rate $k$,
which accounts for these three processes. With incompressible cells, the
continuity equation then reads:

$\nabla\cdot\mathbf{v}=k\ .$ (18)

There is evidence that the effective cell-proliferation rate is influenced by
the cellular mechanical environment, and particularly by tissue pressure
basan2009 ; montel2011 ; delarue2013 ; delarue2014 . At a specific pressure
called the homeostatic pressure $P_{\rm h}$, cell divisions and cell deaths
and extrusions balance on average, such that $k(P_{\rm h})=0$. Far from the
defect, the tissue is in its homeostatic state $P(r=\infty)=P_{\rm h}$. In the
vicinity of the homeostatic state, the pressure-dependent division rate reads
to linear order:

$k(P)=-\frac{1}{\kappa}(P-P_{\rm h})\ ,$ (19)

where $\kappa$ is a phenomenological coefficient. It is positive to ensure the
monolayer stability and can be interpreted as an effective, long-term, bulk
viscosity.

To solve the force-balance equation (8) in this context, we write the velocity
field as a sum of a curl-free part and a divergence-free part, also known as
the Helmholtz decomposition: $\mathbf{v}=\mathbf{\nabla}\times(\psi\
\mathbf{e}_{z})+\mathbf{\nabla}\phi$. Taking the curl and the divergence of
eq. (8) leads to:

$\displaystyle\Delta\left[\Delta\tilde{\psi}-\tilde{\psi}\right]$
$\displaystyle=-s\,\frac{\sin\theta}{\tilde{r}^{2}}$ (20)
$\displaystyle\Delta\left[\Delta\tilde{\phi}-\frac{\eta}{\eta+\kappa}\tilde{\phi}\right]$
$\displaystyle=-s\,\frac{\eta}{\eta+\kappa}\frac{\cos{\theta}}{\tilde{r}^{2}}\
,$ (21)

using the same adimensional variables and functions as in Section 2.1. These
equations are solved in D in a similar way as in sec. 2.1 to obtain the
velocity field. The result for the self-advection velocity of the defect
reads:

$\mathbf{v}_{0}=-\frac{\pi}{4}\left(1+\sqrt{\frac{\eta}{\eta+\kappa}}\right)\times\frac{\zeta\Delta\mu}{2\sqrt{\xi\eta}}\
\mathbf{e}_{x}\ .$ (22)

The pressure-dependent division rate promotes motion of the defect in the same
direction as what is dictated by the active force. We recover the
incompressible limit of eq. (12) in the case of an infinite $\kappa$, that is
when the cell number is conserved. In the other limit $\kappa=0$,
corresponding to an infinite response of the effective proliferation rate to
pressure variations, the velocity of the defect is multiplied by two. In a
contractile nematic, active forces associated with a $+1/2$ defect create a
low-pressure environment at the head of the defect, which promotes cell
division, and the opposite situation at the tail, as illustrated in fig. 3a.

Figure 3: a,b. Colormap of the dimensionless divergence of the velocity—or
equivalently of the net division rate $k$—around a $+1/2$ (a) and a $-1/2$ (b)
defect. The director-field orientation is represented by solid, thin lines. In
the contractile case represented here ($\zeta<0$), there is an asymmetric
divergence profile for the $+1/2$ defect, with a positive divergence (orange)
at the head and a negative one (blue) at the tail. For the $-1/2$ defect,
there is a three-fold symmetry alternating between a positive (orange) and a
negative (blue) divergence around the defect. c,d. Radial dependence of the
dimensionless divergence for a $+1/2$ (c,
$\nabla\cdot\tilde{\mathbf{v}}=D(\tilde{r})\cos\theta$) and a $-1/2$ (d,
$\nabla\cdot\tilde{\mathbf{v}}^{-}=D^{-}(\tilde{r})\cos 3\theta$) defect. The
amplitude of this quantity is maximum at the core of the defect and decays
with the distance to the core.

At large distances, the divergence of the velocity field vanishes, as
illustrated in fig. 3c. This pattern of division and death/extrusion further
induces motion form head to tail, similarly to the active force. For extensile
systems, the situation is reversed. For a $-1/2$ defect, the pattern of
effective cell divisions follows a three-fold symmetry. Corresponding
expressions are given in E.2 and illustrated in fig. 3b,d. Interestingly, such
a divergence profile has been observed in collective self-organizations of
bacteria copenhagen2021 .

##  3 Coupling between the velocity field and orientation gradients of the
nematic director

All calculations in section 2 are performed in the limit of vanishing
rotational viscosity $\gamma$. In this limit, the orientation of the director
is obtained by minimizing the free energy eq. (1). Consequently, the velocity
field was computed with the fixed director orientation $\varphi=\theta/2$,
corresponding to the orientation field of a passive defect at equilibrium. We
now relax this hypothesis to look at the coupling between the flow field and
the orientation field of the director. We therefore expect that activity
influences the orientation of the director, further modifying the flow field
with respect to the results of section 2. Since no global analytical solution
can be found to this new system of equations, we compute the velocity field in
perturbation in activity. We limit our investigation to first-order terms in
activity. Since the flow field is already of order one in activity with the
passive defect orientation $\varphi=\theta/2$, the first-order modifications
of the orientation contribute only to a second-order term in this velocity
field. The first-order contribution of activity to the velocity field can
therefore be computed using the passive defect orientation $\varphi=\theta/2$.

As a first approximation, we consider the case where the flow-alignment
parameter $\nu$ vanishes, and we further neglect the Ericksen stress tensor
$\bm{\sigma}^{\rm E}$. The passive contribution to the stress tensor then
reads

$\bm{\sigma}^{\rm
p}=2\eta\,\tilde{\mathbf{u}}+\frac{1}{2}\left(\mathbf{h}\otimes\mathbf{p}-\mathbf{p}\otimes\mathbf{h}\right)-P\,\mathbf{1}\
,$ (23)

and the molecular field $\mathbf{h}$, given by eq. (6), reads at steady state:

$\mathbf{h}=\gamma\left[(\mathbf{v}\cdot\nabla)\mathbf{p}+\bm{\omega}\cdot\mathbf{p}\right]\,.$
(24)

###  3.1 Flow field and self-advection velocity

The active and passive contributions to the stress tensor as given by eqs. (5)
and (23) lead to the force-balance condition:

$\eta\Delta\mathbf{v}+\frac{1}{2}\mathbf{\nabla}\times(h_{\perp}\mathbf{e}_{z})-\mathbf{\nabla}P-\xi(\mathbf{v}+\mathbf{v_{0}})-\frac{\zeta\Delta\mu}{2r}\mathbf{e}_{x}=0\
.$ (25)

The perpendicular component of the molecular field $h_{\perp}$ reads, using
the equilibrium orientational field $\varphi=\theta/2$:

$h_{\perp}=\frac{\gamma}{2}\left(\frac{v_{\theta}}{r}-\omega_{\rm s}\right)\
,$ (26)

where $\omega_{\rm s}=(\nabla\times\mathbf{v})\cdot\mathbf{e}_{z}$ is the
vorticity. With this expression, the curl of the force balance equation (25)
leads to:

$\Delta\left[\Delta\bar{\psi}-\frac{\lambda}{\bar{r}}\partial_{\bar{r}}\bar{\psi}-\bar{\psi}\right]=-s\frac{\sin\theta}{\bar{r}^{2}}\,,$
(27)

where we have adimensionalized the equation using the same procedure as in
sec. 2.1, except with an effective viscosity $\bar{\eta}=\eta+\gamma/4$, and
where $\lambda=\gamma/(4\bar{\eta})$. More precisely, we have defined a new
lengthscale $\bar{L}=\sqrt{\bar{\eta}/\xi}$ and a new timescale
$\bar{\tau}=2\bar{\eta}/(|\zeta|\Delta\mu)$, defining the dimensionless
variable $\bar{r}=r/\bar{L}$ and the dimensionless field
$\bar{\psi}=(\bar{\tau}/\bar{L}^{2})\psi$. The term inside the brackets of the
Laplace operator in eq. (27) is a transformed version of the modified Bessel
equation bowman1958 . Similarly to what was done in sec. 2.1 and as detailed
in B.2, we determine the velocity field. The resulting self-advection velocity
reads:

$\mathbf{v_{0}}=-c(\lambda)\times\frac{\zeta\Delta\mu}{2\sqrt{\xi\eta}}\mathbf{e}_{x}\
,$ (28)

where

$c(\lambda)=\frac{\Gamma\left[1-\beta/2\right]\
\Gamma\left[1+1/(2\beta)\right]}{2\sqrt{\lambda}\
\Gamma\left[3/2-\beta/2\right]\ \Gamma\left[3/2+1/(2\beta)\right]}$ (29)

with $\Gamma$ the Euler gamma function abramowitz1964  and
$\beta=\sqrt{1+\lambda^{2}/4}+\lambda/2$.

We plot in fig. 4 the dimensionless amplitude $\tilde{v}_{0}$ of the self-
advection velocity as a function of $\lambda$, both not considering and
considering the Ericksen stress (see section 3.4).

Figure 4: Self-advection velocity of a $+1/2$ defect in units of
$\zeta\Delta\mu/(2\sqrt{\xi\eta})$ and as a function of $\lambda$. The dotted
horizontal line shows the value $\pi/4$, corresponding to eq. (12) when
$\lambda=0$. Two cases are presented: without the Ericksen contribution to the
stress (blue curve), as determined by eq. (28), and with the Ericksen
contribution (orange curve), determined numerically (see section 3.4).

The rotational viscosity overall decreases this amplitude, which is expected
since it corresponds to an additive source of dissipation.

###  3.2 Stall force

Proceeding similarly as in section 2.2, the in-core force balance is unchanged
as compared with eq. (15), since there is no nematic order within the core
region, and the out-core force balance condition is given by eq. (25).
Imposing continuity of the flow and stress fields at the boundary $r=a$
between the in-core and out-core regions as detailed in C.2 gives, to leading
order in $a/\bar{L}$:

$f=-\sqrt{\frac{\bar{\eta}}{\xi}}\ \zeta\Delta\mu\
g(\lambda)\left(a\sqrt{\frac{\xi}{\bar{\eta}}}\right)^{\beta-1}\ ,$ (30)

where

$\displaystyle g(\lambda)=$ $\displaystyle 2^{\beta-\lambda}\
\Gamma\left[1-\beta/2\right]\ \Gamma\left[1+1/(2\beta)\right]$
$\displaystyle\times
6\pi(2-\lambda)\frac{\left[\lambda(4-\lambda(1+\beta))+24(\beta-\lambda/2-1)\right]}{(48-20\lambda^{2}+\lambda^{4})\Gamma(1+\beta-\lambda/2)}\,.$
(31)

The force vanishes as a power law of the core size $a$ with an exponent
$\beta-1$, which varies between 0 and $(\sqrt{5}-1)/2$ as $\lambda$ varies
between 0 and 1. Compared to the previous case $\gamma=0$, the stall force
here is comparatively smaller in the limit of a small core size compared to
the hydrodynamic length $\bar{L}$, with a power-law dependence on the core
size $a$ rather than a logarithmic one.

###  3.3 First-order active correction to the director orientation

We can now investigate the effect of the flow on the orientation of the
director, to first order in activity. With a non-vanishing rotational
viscosity $\gamma$, the perpendicular component of the molecular field is
given by eq. (26). Given the definition of the molecular field and the
expression of the free energy eq. (1), we obtain, in terms of the stream
function $\psi$,

$K\Delta\varphi=\frac{\gamma}{2}\left[\Delta\psi-\frac{\partial_{r}\psi}{r}\right]\
.$ (32)

For reasons similar to those presented in sec. 2.1, the stream function
remains of the form $\psi(r,\theta)=\psi(r)\sin\theta$, and we obtain a
similar form for the term proportional to the rotational viscosity $\gamma$ in
the expression of $\varphi(r,\theta)$:

$\displaystyle\varphi(r,\theta)=\frac{\theta}{2}+\frac{\gamma}{2K}\left\\{rv_{r}(r)\right.$
$\displaystyle\quad+\left.\frac{1}{2}\left(r\int_{0}^{r}{\frac{v_{\theta}(u)}{u}\>\mathop{}\\!\mathrm{d}u}-\frac{1}{r}\int_{0}^{r}{uv_{\theta}(u)\>\mathop{}\\!\mathrm{d}u}\right)\right\\}\sin\theta\
.$ (33)

Let us now investigate the behavior of the correction
$\delta\varphi(r,\theta)$ to the passive orientation
$\varphi^{0}(r,\theta)=\theta/2$ at large distances $r\gg\bar{L}$.
Asymptotically, we get:

$\delta\varphi(r,\theta)\underset{r\gg\bar{L}}{\approx}\frac{\gamma}{4K}r\ln{\left(\frac{r}{r_{0}}\right)}v_{0}\,,$
(34)

where $r_{0}$ is an arbitrary length that has been introduced based on
dimensional analysis. Since this expression diverges with $r$, the solution
given by eq. (3.3) is only valid up to a characteristic lengthscale
$L_{\text{defect}}\propto K/(\gamma v_{0})$. Using the scaling of eq. (28), we
define

$L_{\text{defect}}=\frac{\sqrt{\xi(\eta+\gamma/4)}}{\gamma/4}\
\frac{K}{|\zeta\Delta\mu|}=\frac{L_{a}^{2}}{\lambda\bar{L}}\,,$ (35)

where $L_{\rm a}=\sqrt{K/|\zeta\Delta\mu|}$ is the active length,
corresponding to the scale above which uniform bands of active nematics become
unstable voituriez2005 ; edwards2009  and active turbulence is expected
giomi2015 ; doostmohammadi2018 . This result gives an insight into the
interplay between active and hydrodynamic stresses. Note that our description
is only valid if $L_{\rm a}\gg\bar{L}$. In this limit, eq. (34) is valid at
distances $r$ from the core region such that $\bar{L}\ll r\ll L_{\rm a}$. If
this is not the case, the homogeneous nematic system does not have a stable
steady state and the passive orientation $\varphi^{0}(r,\theta)=\theta/2$ is
spontaneously destroyed over length scales of order $\bar{L}$.

###  3.4 Contribution of the Ericksen stress tensor to the self-advection
velocity

Until now, we have only considered the coupling between the velocity field and
the nematic order through the torques associated with the antisymmetric
component of the stress tensor. In the expression of the stress tensor, we
have however neglected the Ericksen stress tensor $\bm{\sigma}^{\rm E}$ as
well as the flow-alignment coupling term proportional to the parameter $\nu$.
In perturbation theory around a homogeneous steady state, neglecting the
Ericksen stress is in general justified by the fact that it is second order in
the gradients of the director orientation, and therefore nonlinear in activity
degennes-prost . Here, however, the Ericksen stress is relevant because the
ground state upon which activity sets in has a non-homogeneous director
orientation $\varphi=\theta/2$. The force density resulting from the Ericksen
stress is given by the Gibbs-Duhem relation (4), which here reads

$\nabla\cdot\bm{\sigma}^{\rm E}=-\frac{h_{\perp}}{2r}\ \mathbf{e}_{\theta}\,.$
(36)

Accounting for the Ericksen stress while keeping $\nu=0$, the curl of the
force-balance condition leads to

$\Delta\left[\Delta\bar{\psi}-\frac{\lambda}{\bar{r}}\partial_{\bar{r}}\bar{\psi}-\bar{\psi}\right]+\frac{\lambda}{\bar{r}}\partial_{\bar{r}}\left[\Delta\bar{\psi}-\frac{\partial_{\bar{r}}\bar{\psi}}{r}\right]=-s\frac{\sin\theta}{\bar{r}^{2}}\,,$
(37)

with the same dimensionless units as in eq. (27). Contrary to the previous
equations for the stream function in this paper, this equation has no direct
analytical solution. Using asymptotic expansions close to the core and at
large distances, we can however determine the velocity of the defect using a
shooting method. We plot in fig. 4 the resulting self-advection velocity of
the defect as a function of $\lambda$ (orange curve), in comparison to the
analytical solution obtained in the absence of Ericksen stress (blue curve):
there is a decrease in amplitude due to the Ericksen stress, which is expected
given that it represents an additional source of rigidity.

##  4 Discussion

In this manuscript, we have studied the flows generated by $\pm 1/2$
topological defects of an active nematic system in two dimensions, which
interacts with its underlying substrate via viscous drag. The active nematic
system lies deep in the nematic phase, with an order parameter of constant,
maximal amplitude. Our study is inspired by monolayers of elongated cells at
confluence. In the case of a vanishing rotational viscosity, we derived
analytical expressions for the flow field. In the case of a $+1/2$ topological
defect, we determined its self-advection velocity and stall force. We showed
that the non-conservation of cell number leads to an increase in the self-
advection velocity. We then investigated the effect of the rotational
viscosity, which couples the velocity field to the nematic order. Taking the
flow-alignment parameter to zero and in the absence of the Ericksen
contribution to the stress, we derived analytical expressions of the self-
advection velocity and stall force, making use of the passive, equilibrium
configuration of the director around the defect. We finally studied the
validity domain of the small departure from the equilibrium nematic-order
configuration, as well as the influence of the Ericksen stress.

Our results show that the self-advection velocity is linear in activity for a
$+1/2$ defect, as it should, as well as its corresponding stall force. The
self-advection velocity of a $+1/2$ defect decreases in the presence of a non-
vanishing rotational viscosity, which adds an extra source of dissipation.
Accordingly, the stall force is also smaller, with a qualitatively different
dependence on the size of the core region of the defect, where the nematic
order vanishes. On the contrary, the non-conservation of the cell number
increases the self-advection velocity, when coupled linearly to the departure
from a homeostatic isotropic stress. Finally, accounting for the Ericksen
stress lowers even further that velocity, with a similar qualitative
interpretation.

Several assumptions used in our study can be questioned. A first restriction
consisted in studying an isolated, topological defect in an infinite domain,
as defects in an actual, free monolayers come into pairs to ensure a vanishing
total topological charge. This assumption allowed for an analytical treatment
that does not depend on finite-domain boundary conditions. Importantly,
section 3.3 provided a characteristic system size, within which our
calculations are valid, and beyond which active effects drive the nematic far
from its passive, equilibrium configuration. At a qualitative level, the
presence of other defects introduces another length scale in the problem,
namely the characteristic distance $d$ between defects. If this distance is
such that $d\gg\bar{L}$, our results remain valid as the defects interact only
very weakly. The interactions between defects becomes relevant in the limit
where $d\lesssim\bar{L}$. In this case, one should consider the nematic layer
as a gas of interacting defects shankar2018 . Within our approach, computing
the flow created by a pair of $+1/2$ and $-1/2$ defects is the next logical
step. A second restriction was to consider only partially the coupling between
flow and orientation. In particular, for a non-vanishing rotational viscosity,
we have limited our study to the case of a vanishing flow-alignment parameter.
Beyond this approximation, no analytical solutions to the flow equations could
be found. Finally, considering proliferation and extrusion is of prominent
relevance for living systems, and we provided a first computation of this
effect on the self-advection velocity of $+1/2$ defects.

The initial goal of this work was to explain the observation of Sarkar _et
al._ in ref sarkar2021  of motionless defects that are preferential sites for
multilayer formation. We computed the self-advection velocity of motile
defects and the force necessary to stall them. However, the origin of the
stalling and the mechanism by which another layer of cells forms at the defect
is still elusive. Accounting for pressure-dependent proliferation shows that,
on average, more extrusion is to be expected compared to the rest of the
monolayer at the tail of $+1/2$, contractile defects. Recent studies by Vafa &
Mahadevan vafa2021  and by Hoffman _et al._ hoffmann2021  show that, when
considering a deformable surface, there is an out-of-plane force at
topological defects that could be responsible for extrusion and multilayering.

###### Acknowledgements.

We thank P. Silberzan and T. Sarkar from Institut Curie for sharing their
experimental results with us and useful discussions. LB received a PhD
fellowship from the doctoral school Physique en Ile-de-France (EDPIF) and
support by the Collège de France foundation and Institut Curie. This work
received support from the grants ANR-11-LABX-0038, ANR-10-IDEX-0001-02.

## References

  * (1) N. D. Mermin, “The topological theory of defects in ordered media,” Reviews of Modern Physics, vol. 51, pp. 591–648, July 1979.
  * (2) M. Kleman, Points, Lines, and Walls: In Liquid Crystals, Magnetic Systems, and Various Ordered Media.  J Wiley, 1983.
  * (3) O. Lehmann, Flüssigkristalle.  Engelmann, 1904.
  * (4) M. Kleman, “Defects in liquid crystals,” Reports on Progress in Physics, vol. 52, pp. 555–654, May 1989.
  * (5) P. Oswald and P. Pieranski, Nematic and Cholesteric Liquid Crystals: Concepts and Physical Properties Illustrated by Experiments.  Boca Raton: CRC Press, Feb. 2005.
  * (6) K. Harth and R. Stannarius, “Topological Point Defects of Liquid Crystals in Quasi-Two-Dimensional Geometries,” Frontiers in Physics, vol. 8, p. 112, 2020.
  * (7) F. Jülicher, K. Kruse, J. Prost, and J. F. Joanny, “Active behavior of the Cytoskeleton,” Physics Reports, vol. 449, pp. 3–28, Sept. 2007.
  * (8) M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, “Hydrodynamics of soft active matter,” Reviews of Modern Physics, vol. 85, pp. 1143–1189, July 2013.
  * (9) L. Giomi, M. J. Bowick, X. Ma, and M. C. Marchetti, “Defect annihilation and proliferation in active nematics,” Physical Review Letters, vol. 110, p. 228101, May 2013.
  * (10) S. P. Thampi, R. Golestanian, and J. M. Yeomans, “Velocity Correlations in an Active Nematic,” Physical Review Letters, vol. 111, p. 118101, Sept. 2013.
  * (11) L. M. Pismen, “Dynamics of defects in an active nematic layer,” Physical Review E, vol. 88, p. 050502, Nov. 2013.
  * (12) T. Elsdale, “Parallel orientation of fibroblasts in vitro,” Experimental Cell Research, vol. 51, pp. 439–450, July 1968.
  * (13) R. Kemkemer, D. Kling, D. Kaufmann, and H. Gruler, “Elastic properties of nematoid arrangements formed by amoeboid cells,” The European Physical Journal E, vol. 1, no. 2, p. 215, 2000.
  * (14) G. Duclos, C. Erlenkämper, J.-F. Joanny, and P. Silberzan, “Topological defects in confined populations of spindle-shaped cells,” Nature Physics, vol. 13, pp. 58–62, Jan. 2017.
  * (15) T. B. Saw, W. Xi, B. Ladoux, and C. T. Lim, “Biological Tissues as Active Nematic Liquid Crystals,” Advanced Materials, vol. 30, p. 1802579, Nov. 2018.
  * (16) C. Blanch-Mercader, P. Guillamat, A. Roux, and K. Kruse, “Integer topological defects of cell monolayers – mechanics and flows,” Physical Review E, vol. 103, p. 012405, June 2021.
  * (17) C. Blanch-Mercader, P. Guillamat, A. Roux, and K. Kruse, “Quantifying material properties of cell monolayers by analyzing integer topological defects,” Physical Rview Letters, vol. 103, p. 028101, June 2021.
  * (18) 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,” Nature Physics, vol. 17, pp. 251–259, Feb. 2021.
  * (19) L. Metselaar, J. M. Yeomans, and A. Doostmohammadi, “Topology and Morphology of Self-Deforming Active Shells,” Physical Review Letters, vol. 123, p. 208001, Nov. 2019.
  * (20) 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, vol. 544, pp. 212–216, Apr. 2017.
  * (21) K. Kawaguchi, R. Kageyama, and M. Sano, “Topological defects control collective dynamics in neural progenitor cell cultures,” Nature, vol. 545, pp. 327–331, May 2017.
  * (22) T. Sarkar, V. Yashunsky, L. Brézin, C. B. Mercader, T. Aryaksama, M. Lacroix, T. Risler, J.-F. Joanny, and P. Silberzan, “Crisscross multilayering of cell sheets,” bioRxiv, preprint, 2021.
  * (23) K. Copenhagen, R. Alert, N. S. Wingreen, and J. W. Shaevitz, “Topological defects promote layer formation in Myxococcus xanthus colonies,” Nature Physics, vol. 17, pp. 211–215, Feb. 2021.
  * (24) P. de Gennes and J. Prost, The Physics of Liquid Crystals.  Oxford Science Publications, 1974.
  * (25) L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek, and M. Cristina Marchetti, “Defect dynamics in active nematics,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 372, p. 20130365, Nov. 2014.
  * (26) J. Rønning, C. M. Marchetti, M. J. Bowick, and L. Angheluta, “Flow around topological defects in active nematic films,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 478, no. 2257, p. 20210879, 2022.
  * (27) G. Toth, C. Denniston, and J. M. Yeomans, “Hydrodynamics of topological defects in nematic liquid crystals,” Physical Review Letters, vol. 88, p. 105504, Feb. 2002.
  * (28) K. Kruse, J. F. Joanny, F. Jülicher, J. Prost, and K. Sekimoto, “Generic theory of active polar gels: A paradigm for cytoskeletal dynamics,” The European Physical Journal E, vol. 16, pp. 5–16, Jan. 2005.
  * (29) J. F. Joanny, F. Jülicher, K. Kruse, and J. Prost, “Hydrodynamic theory for multi-component active polar gels,” New Journal of Physics, vol. 9, pp. 422–422, Nov. 2007.
  * (30) J.-F. Joanny and J. Prost, “Active gels as a description of the actin-myosin cytoskeleton,” HFSP Journal, vol. 3, pp. 94–104, Apr. 2009.
  * (31) G. Duclos, S. Garcia, H. G. Yevick, and P. Silberzan, “Perfect nematic order in confined monolayers of spindle-shaped cells,” Soft Matter, vol. 10, no. 14, pp. 2346–2353, 2014.
  * (32) G. Duclos, C. Blanch-Mercader, V. Yashunsky, G. Salbreux, J.-F. Joanny, J. Prost, and P. Silberzan, “Spontaneous shear flow in confined cellular nematics,” Nature Physics, vol. 14, pp. 728–732, July 2018.
  * (33) R. Alert and X. Trepat, “Physical Models of Collective Cell Migration,” Annual Review of Condensed Matter Physics, vol. 11, pp. 77–101, Mar. 2020.
  * (34) P. C. Martin, O. Parodi, and P. S. Pershan, “Unified Hydrodynamic Theory for Crystals, Liquid Crystals, and Normal Fluids,” Physical Review A, vol. 6, pp. 2401–2420, Dec. 1972.
  * (35) D. Forster, T. C. Lubensky, P. C. Martin, J. Swift, and P. S. Pershan, “Hydrodynamics of Liquid Crystals,” Physical Review Letters, vol. 26, pp. 1016–1019, Apr. 1971.
  * (36) P. M. Chaikin and T. C. Lubensky, Principles of Condensed Matter Physics.  Cambridge: Cambridge University Press, 1995.
  * (37) F. Jülicher, S. W. Grill, and G. Salbreux, “Hydrodynamic theory of active matter,” Reports on Progress in Physics, vol. 81, p. 076601, July 2018\.
  * (38) M. Abramowitz and I. Stegun, Handbook of Mathematical functions.  National Bureau of Standards, Applied Mathematics Series - 55, 1964\.
  * (39) S. Garcia, E. Hannezo, J. Elgeti, J.-F. Joanny, P. Silberzan, and N. S. Gov, “Physics of active jamming during collective cellular motion in a monolayer,” Proceedings of the National Academy of Sciences, vol. 112, pp. 15314–15319, Dec. 2015.
  * (40) M. Basan, T. Risler, J.-F. Joanny, X. Sastre-Garau, and J. Prost, “Homeostatic competition drives tumor growth and metastasis nucleation,” HFSP Journal, vol. 3, pp. 265–272, Aug. 2009.
  * (41) F. Montel, M. Delarue, J. Elgeti, L. Malaquin, M. Basan, T. Risler, B. Cabane, D. Vignjevic, J. Prost, G. Cappello, and J.-F. Joanny, “Stress Clamp Experiments on Multicellular Tumor Spheroids,” Physical Review Letters, vol. 107, p. 188102, Oct. 2011.
  * (42) M. Delarue, F. Montel, O. Caen, J. Elgeti, J.-M. Siaugue, D. Vignjevic, J. Prost, J.-F. Joanny, and G. Cappello, “Mechanical Control of Cell flow in Multicellular Spheroids,” Physical Review Letters, vol. 110, p. 138103, Mar. 2013.
  * (43) M. Delarue, F. Montel, D. Vignjevic, J. Prost, J.-F. Joanny, and G. Cappello, “Compressive Stress Inhibits Proliferation in Tumor Spheroids through a Volume Limitation,” Biophysical Journal, vol. 107, pp. 1821–1828, Oct. 2014.
  * (44) F. Bowman, Introduction to Bessel Functions.  Dover Publications Inc., 1958.
  * (45) R. Voituriez, J.-F. Joanny, and J. Prost, “Spontaneous flow transition in active polar gels,” Europhysics Letters (EPL), vol. 70, pp. 404–410, May 2005.
  * (46) S. A. Edwards and J. M. Yeomans, “Spontaneous flow states in active nematics: A unified picture,” EPL (Europhysics Letters), vol. 85, p. 18008, Jan. 2009.
  * (47) L. Giomi, “Geometry and Topology of Turbulence in Active Nematics,” Physical Review X, vol. 5, p. 031003, July 2015.
  * (48) A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, “Active nematics,” Nature Communications, vol. 9, p. 3246, Dec. 2018.
  * (49) S. Shankar, S. Ramaswamy, and M. C. Marchetti, “The low noise phase of a 2d active nematic,” Physical Review E, vol. 97, p. 012707, Jan. 2018.
  * (50) F. Vafa and L. Mahadevan, “Active nematic defects and epithelial morphogenesis,” arXiv, preprint, May 2021.
  * (51) L. A. Hoffmann, L. N. Carenza, J. Eckert, and L. Giomi, “Defect-mediated morphogenesis,” arXiv, preprint, May 2021.

##  Appendix A Modified Bessel functions

###  A.1 Definition

The modified Bessel functions $I_{\alpha},K_{\alpha}$ are the two general
solutions of the following equation:

$r^{2}f^{\prime\prime}(r)+rf^{\prime}(r)-(\alpha^{2}+r^{2})f(r)=0\ .$ (38)

A related equation to the modified Bessel equation (38) is given by Bowman
bowman1958 :

$r^{2}f^{\prime\prime}(r)+(1-2a)rf^{\prime}(r)-(b^{2}c^{2}r^{2c}-a^{2}+d^{2}c^{2})f(r)=0\
.$ (39)

Two independent solutions of eq. (39) are $r^{a}I_{d}(br^{c})$ and
$r^{a}K_{d}(br^{c})$.

###  A.2 Derivatives

The first-order derivatives of the modified Bessel functions are given by:

$\displaystyle I_{\alpha}^{\prime}(r)$
$\displaystyle=I_{\alpha-1}(r)-\frac{\alpha}{r}I_{\alpha}$ (40) $\displaystyle
K_{\alpha}^{\prime}(r)$
$\displaystyle=-\left(K_{\alpha-1}(r)+\frac{\alpha}{r}K_{\alpha}\right)$ (41)
$\displaystyle I_{0}^{\prime}(r)$ $\displaystyle=I_{1}(r)$ (42) $\displaystyle
K_{0}^{\prime}(r)$ $\displaystyle=-K_{1}(r)$ (43)

###  A.3 Asymptotic expansions

We also use the following asymptotic expansions:

$\displaystyle I_{\alpha}(r)$ $\displaystyle\underset{r\rightarrow
0}{\simeq}\frac{1}{\Gamma(\alpha+1)}\left(\frac{r}{2}\right)^{\alpha}$ (44)
$\displaystyle K_{\alpha}(r)$ $\displaystyle\underset{r\rightarrow
0}{\simeq}\begin{cases}\frac{\Gamma(\alpha)}{2}\left(\frac{2}{r}\right)^{\alpha}\quad\text{if
}\alpha>0\\\ -\ln(\frac{r}{2})-\gamma\quad\text{if }\alpha=0\end{cases}$ (45)

where $\gamma$ denotes here the Euler’s constant. At infinity, we use

$\displaystyle I_{\alpha}(r)$
$\displaystyle\underset{r\rightarrow+\infty}{\simeq}\frac{e^{r}}{\sqrt{2\pi
r}}$ (46) $\displaystyle K_{\alpha}(r)$
$\displaystyle\underset{r\rightarrow+\infty}{\simeq}\sqrt{\frac{\pi}{2r}}e^{-r}$
(47)

##  Appendix B Velocity field

###  B.1 Limit of a vanishing rotational viscosity

This appendix is dedicated to the full computation of the velocity starting
from eq. (9). As justified in the main text, the stream function reads
$\tilde{\psi}(\tilde{r},\theta)=\tilde{\psi}(\tilde{r})\sin\theta$.
Integrating one Laplace operator gives:

$\Delta\tilde{\psi}(\tilde{r},\theta)-\tilde{\psi}(\tilde{r},\theta)=\left[A\tilde{r}+\frac{B}{\tilde{r}}+s\right]\sin\theta\
,$ (48)

where $A$ and $B$ are integration constants. The radial dependence of the
stream function then satisfies:

$\frac{\mathop{}\\!\mathrm{d}^{2}\tilde{\psi}(r)}{\mathop{}\\!\mathrm{d}\tilde{r}^{2}}+\frac{1}{\tilde{r}}\frac{\mathop{}\\!\mathrm{d}\tilde{\psi}(\tilde{r})}{\mathop{}\\!\mathrm{d}\tilde{r}}-\left(1+\frac{1}{\tilde{r}^{2}}\right)\tilde{\psi}(r)=A\tilde{r}+\frac{B}{\tilde{r}}+s\
.$ (49)

The homogeneous solution $\tilde{\psi}^{0}$ to this equation reads:

$\tilde{\psi}^{0}(\tilde{r})=A_{0}I_{1}(\tilde{r})+B_{0}K_{1}(\tilde{r})\ ,$
(50)

where $A_{0}$ and $B_{0}$ are integration constants. The Wronskian associated
to eq. (49) reads

$I_{1}(\tilde{r})K^{\prime}_{1}(\tilde{r})-I^{\prime}_{1}(\tilde{r})K_{1}(\tilde{r})=-\frac{1}{\tilde{r}}\
,$ (51)

which leads to:

$\displaystyle\tilde{\psi}(\tilde{r})$
$\displaystyle=\left\\{I_{1}(\tilde{r})\left(A_{0}-\int_{\tilde{r}}^{+\infty}{K_{1}(u)(Au^{2}+su+B)\>\mathop{}\\!\mathrm{d}u}\right)\right.$
$\displaystyle+\left.K_{1}(\tilde{r})\left(B_{0}-\int_{0}^{\tilde{r}}{I_{1}(u)(Au^{2}+su+B)\>\mathop{}\\!\mathrm{d}u}\right)\right\\}\
.$ (52)

The velocity field
$\mathbf{\tilde{v}}=\tilde{v}_{r}(\tilde{r})\cos\theta\>\mathbf{e_{r}}+\tilde{v}_{\theta}(\tilde{r})\sin\theta\>\mathbf{e_{\theta}}$
is obtained from the derivatives of the stream function:
$\tilde{v}_{r}(\tilde{r})=\tilde{\psi}(\tilde{r})/\tilde{r}$ and
$\tilde{v}_{\theta}(\tilde{r})=-\mathop{}\\!\mathrm{d}\tilde{\psi}(\tilde{r})/\mathop{}\\!\mathrm{d}\tilde{r}$.
We obtain

$\displaystyle\tilde{v}_{r}(\tilde{r})$
$\displaystyle=\tilde{r}^{-1}\left[I_{1}(\tilde{r})\left(A_{0}-\int_{\tilde{r}}^{+\infty}{K_{1}(u)(Au^{2}+su+B)\>\mathop{}\\!\mathrm{d}u}\right)\right.$
$\displaystyle\qquad\qquad\quad+\left.K_{1}(\tilde{r})\left(B_{0}-\int_{0}^{\tilde{r}}{I_{1}(u)(Au^{2}+su+B)\>\mathop{}\\!\mathrm{d}u}\right)\right]$
(53) $\displaystyle\tilde{v}_{\theta}(\tilde{r})$
$\displaystyle=-\left[\left(I_{0}(\tilde{r})-\frac{I_{1}(\tilde{r})}{\tilde{r}}\right)\left(A_{0}-\int_{\tilde{r}}^{+\infty}{K_{1}(u)(Au^{2}+su+B)\>\mathop{}\\!\mathrm{d}u}\right)\right.$
$\displaystyle\qquad\quad-\left.\left(K_{0}(\tilde{r})+\frac{K_{1}(\tilde{r})}{\tilde{r}}\right)\left(B_{0}-\int_{0}^{\tilde{r}}{I_{1}(u)(Au^{2}+su+B)\>\mathop{}\\!\mathrm{d}u}\right)\right]\
.$ (54)

Imposing a finite velocity at the origin leads to $B_{0}=0$ and $B=0$, as a
finite velocity at infinity leads to $A_{0}=0$. We then obtain
$\tilde{v}_{r}(\tilde{r}=0)=-(A+s\pi/4)$ and
$\tilde{v}_{\theta}(\tilde{r}=0)=A+s\pi/4$ at the origin, and
$\tilde{v}_{r}(\tilde{r}=\infty)=-A$ and
$\tilde{v}_{\theta}(\tilde{r}=\infty)=A$ at infinity. Imposing
$\tilde{v}_{r}(0)=\tilde{v}_{\theta}(0)=0$ yields $A=-s\pi/4$, and imposing
$\tilde{v}_{r}(\tilde{r}=\infty)=-\tilde{v}_{0}$ and
$\tilde{v}_{\theta}(\tilde{r}=\infty)=\tilde{v}_{0}$ yields $A=\tilde{v}_{0}$.
Coming back to the physical units, these lead to the self-advection velocity
given by eq. (12).

###  B.2 Velocity field with a finite rotational viscosity

We start the computation from eq. (27) with the dependence
$\tilde{\psi}(\tilde{r},\theta)=\tilde{\psi}(\tilde{r})\sin\theta$.
Integrating one Laplace operator gives:

$\frac{\mathop{}\\!\mathrm{d}^{2}\tilde{\psi}(r)}{\mathop{}\\!\mathrm{d}\tilde{r}^{2}}+\frac{(1-\lambda)}{\tilde{r}}\frac{\mathop{}\\!\mathrm{d}\tilde{\psi}(r)}{\mathop{}\\!\mathrm{d}\tilde{r}}-(1+\frac{1}{\tilde{r}^{2}})\tilde{\psi}=A\>\tilde{r}+\frac{B}{\tilde{r}}+s\
.$ (55)

The left-hand side of this equation is of the form (39) with $a=\lambda/2$,
$b=1$, $c=1$, and $d=\sqrt{1+\lambda^{2}/4}$. Following A, two independent
homogeneous solutions to eq. (55) are
$\tilde{r}^{\lambda/2}I_{\alpha}(\tilde{r}),\tilde{r}^{\lambda/2}K_{\alpha}(\tilde{r})$,
with $\alpha=\sqrt{1+\lambda^{2}/4}$. Following a similar procedure as in B.1,
the dimensionless velocity field has the form:

$\displaystyle\tilde{v}_{r}=\tilde{r}^{\lambda/2-1}$
$\displaystyle\left\\{I_{\alpha}(\tilde{r})\left(A_{0}-\int_{\tilde{r}}^{+\infty}{u^{-\lambda/2}K_{\alpha}(u)(Au^{2}+su+B)\>\mathop{}\\!\mathrm{d}u}\right)\right.$
$\displaystyle+\left.K_{\alpha}(\tilde{r})\left(B_{0}-\int_{0}^{\tilde{r}}{u^{-\lambda/2}I_{\alpha}(u)(Au^{2}+su+B)\>\mathop{}\\!\mathrm{d}u}\right)\right\\}$
$\displaystyle\times\cos\theta$ (56)
$\displaystyle\tilde{v}_{\theta}=\tilde{r}^{\lambda/2}$
$\displaystyle\left\\{\left(\frac{\alpha-\lambda/2}{\tilde{r}}I_{\alpha}(\tilde{r})-I_{\alpha-1}(\tilde{r})\right)\times\right.$
$\displaystyle\left.\quad\left(A_{0}-\int_{\tilde{r}}^{+\infty}{u^{-\lambda/2}K_{\alpha}(u)(Au^{2}+su+B)\>\mathop{}\\!\mathrm{d}u}\right)\right.$
$\displaystyle+\left.\left(\frac{\alpha-\lambda/2}{\tilde{r}}K_{\alpha}(\tilde{r})+K_{\alpha-1}(\tilde{r})\right)\right.\times$
$\displaystyle\left.\quad\left(B_{0}-\int_{0}^{\tilde{r}}{u^{-\lambda/2}I_{\alpha}(u)(Au^{2}+su+B)\>\mathop{}\\!\mathrm{d}u}\right)\right\\}\sin\theta\
.$ (57)

The boundary conditions still set $A_{0}=0$, $B_{0}=0$, $B=0$, and
$A=\tilde{v}_{0}$, similarly as in B.1. The asymptotic expansion of the
velocity field close to the core reads:

$\displaystyle\tilde{v}_{r}(\tilde{r})$ $\displaystyle\underset{\tilde{r}\ll
1}{\simeq}-\frac{2^{-\alpha}}{\Gamma(\alpha+1)}(sC_{1}^{\lambda}+\tilde{v}_{0}C_{2}^{\lambda})\tilde{r}^{\lambda/2+\alpha-1}$
(58) $\displaystyle\tilde{v}_{\theta}(\tilde{r})$
$\displaystyle\underset{\tilde{r}\ll
1}{\simeq}\frac{2^{-\alpha}(\lambda/2+\alpha)}{\Gamma(\alpha+1)}(sC_{1}^{\lambda}+\tilde{v}_{0}C_{2}^{\lambda})\tilde{r}^{\lambda/2+\alpha-1}\
,$ (59)

with

$\displaystyle C_{1}^{\lambda}$
$\displaystyle=\int_{0}^{\infty}{\mathop{}\\!\mathrm{d}u\>u^{1-\lambda/2}K_{\alpha}(u)}$
$\displaystyle=2^{-\lambda/2}\Gamma\left[1-\frac{1}{2}\left(\alpha+\frac{\lambda}{2}\right)\right]\>\Gamma\left[1+\frac{1}{2}\left(\alpha-\frac{\lambda}{2}\right)\right]$
(60) $\displaystyle C_{2}^{\lambda}$
$\displaystyle=\int_{0}^{\infty}{\mathop{}\\!\mathrm{d}u\>u^{2-\lambda/2}K_{\alpha}(u)}$
$\displaystyle=2^{1-\lambda/2}\Gamma\left[\frac{1}{2}\left(3-\alpha-\frac{\lambda}{2}\right)\right]\>\Gamma\left[\frac{1}{2}\left(3+\alpha-\frac{\lambda}{2}\right)\right]\
.$ (61)

Contrary to the case of a vanishing rotational viscosity, the velocity field
in the defect reference frame vanishes close to the core for any finite value
of $\tilde{v}_{0}$, since the exponent $\lambda/2+\alpha-1$ is positive. To
set the value of $\tilde{v}_{0}$, we must consider the tangential stress
$\sigma_{r\theta}$. Using the angular dependence of the stream function and
force balance, its dimensionless version reads

$\displaystyle\tilde{\sigma}_{r\theta}(\tilde{r},\theta)=\left(2\frac{\lambda-1}{\tilde{r}}(\tilde{v}_{\theta}(\tilde{r})+\tilde{v}_{r}(\tilde{r}))-\tilde{r}\tilde{v}_{r}(\tilde{r})-\tilde{v}_{0}\>\tilde{r}\right)\sin\theta\
.$ (62)

Since $\lambda/2+\alpha-2$ is negative, $\tilde{v}_{0}$ must equal
$-sC_{1}^{\lambda}/C_{2}^{\lambda}$ for this tangential stress to remain
finite at the origin. This leads to eq. (28) in physical units.

##  Appendix C Stall force

###  C.1 Limit of vanishing rotational viscosity

In the limit of vanishing rotational viscosity, the stall force is determined
by the velocity and pressure solutions to eqs. (14) and (15), together with
the incompressibility condition $\nabla\cdot\mathbf{v}=0$. Quantities defined
inside the core $r<a$ bear the superscript ‘c’. Using the dimensionless units
of sec. 2, the stream function satisfies:

$\displaystyle\Delta\left[\Delta\tilde{\psi}-\tilde{\psi}\right]$
$\displaystyle=-s\frac{\sin\theta}{\tilde{r}^{2}}\qquad\tilde{r}>\tilde{a}$
(63)
$\displaystyle\Delta\left[\Delta\tilde{\psi}^{\mathrm{c}}-\tilde{\psi}^{\mathrm{c}}\right]$
$\displaystyle=-\tilde{f}\frac{\sin\theta}{\tilde{r}^{2}}\qquad\tilde{r}<\tilde{a}\
,$ (64)

with $\tilde{a}=a/L$ and $\tilde{f}=f/(\pi a\zeta\Delta\mu)$ the normalized
core radius and overall force applied to the core. The solution for the
velocity field outside the core is given by eqs. (53) and (54) with $A_{0}=0$
and $A=\tilde{v}_{0}$ as in B.1. The velocity field inside the core however
reads:

$\displaystyle\tilde{v}_{r}^{\mathrm{c}}(\tilde{r})$
$\displaystyle=\tilde{r}^{-1}\left[I_{1}(\tilde{r})\left(A_{0}^{\mathrm{c}}+\int_{0}^{\tilde{r}}{K_{1}(u)(A^{\mathrm{c}}u^{2}+\tilde{f}u+B^{\mathrm{c}})\>\mathop{}\\!\mathrm{d}u}\right)\right.$
$\displaystyle\qquad+\left.K_{1}(\tilde{r})\left(B^{\mathrm{c}}_{0}-\int_{0}^{\tilde{r}}{I_{1}(u)(A^{\mathrm{c}}u^{2}+\tilde{f}u+B^{\mathrm{c}})\>\mathop{}\\!\mathrm{d}u}\right)\right]$
(65) $\displaystyle\tilde{v}^{\mathrm{c}}_{\theta}(\tilde{r})$
$\displaystyle=-\left[\left(I_{0}(\tilde{r})-\frac{I_{1}(\tilde{r})}{\tilde{r}}\right)\times\right.$
$\displaystyle\qquad\quad\left.\left(A^{\mathrm{c}}_{0}+\int_{0}^{\tilde{r}}{K_{1}(u)(A^{\mathrm{c}}u^{2}+\tilde{f}u+B^{\mathrm{c}})\>\mathop{}\\!\mathrm{d}u}\right)\right.$
$\displaystyle\quad-\left(K_{0}(\tilde{r})+\frac{K_{1}(\tilde{r})}{\tilde{r}}\right)\times$
$\displaystyle\qquad\quad\left.\left(B^{\mathrm{c}}_{0}-\int_{0}^{\tilde{r}}{I_{1}(u)(A^{\mathrm{c}}u^{2}+\tilde{f}u+B^{\mathrm{c}})\>\mathop{}\\!\mathrm{d}u}\right)\right]\
.$ (66)

The boundary condition $\mathbf{v}^{\mathrm{c}}(\tilde{r}=0)=\mathbf{0}$
imposes that $B^{\mathrm{c}}$, $B^{\mathrm{c}}_{0}$, and $A_{0}^{\mathrm{c}}$
vanish. The other integration constants are set by imposing the continuity of
the velocity and stress fields at the boundary of the core $r=a$. Introducing
$\tilde{P}$, the pressure normalized by $|\zeta|\Delta\mu/2$ and taking the
divergence of the force-balance eqs. (14) and (15), we obtain:

$\displaystyle\Delta\tilde{P}$
$\displaystyle=s\frac{\cos\theta}{\tilde{r}^{2}}\qquad\tilde{r}>\tilde{a}$
(67) $\displaystyle\Delta\tilde{P}^{\mathrm{c}}$
$\displaystyle=\tilde{f}\frac{\cos\theta}{\tilde{r}^{2}}\qquad\tilde{r}<\tilde{a}\
.$ (68)

These equations are solved by:

$\displaystyle\tilde{P}(\tilde{r},\theta)$
$\displaystyle=\left[A_{*}\tilde{r}+\frac{B_{*}}{\tilde{r}}-s\right]\cos\theta$
(69) $\displaystyle\tilde{P}^{\mathrm{c}}(r,\theta)$
$\displaystyle=\left[A_{*}^{\mathrm{c}}\tilde{r}+\frac{B_{*}^{\mathrm{c}}}{\tilde{r}}-\tilde{f}\right]\cos\theta\
,$ (70)

where $A_{*}$, $B_{*}$, $A_{*}^{\mathrm{c}}$, and $B_{*}^{\mathrm{c}}$ are
integration constants. Pressure and velocity are linked by force balance,
imposing

$\displaystyle B_{*}$ $\displaystyle=B$ (71) $\displaystyle
B_{*}^{\mathrm{c}}$ $\displaystyle=B^{\mathrm{c}}$ (72) $\displaystyle A_{*}$
$\displaystyle=0$ (73) $\displaystyle A_{*}^{\mathrm{c}}$
$\displaystyle=A^{\mathrm{c}}-\tilde{v}_{0}\ .$ (74)

The normal and tangential components of the stress read
$\sigma_{rr}=2\eta\partial_{r}v_{r}-P$ and
$\sigma_{r\theta}=\eta[(\partial_{\theta}v_{r}-v_{\theta})/r+\partial_{r}v_{\theta}]$,
respectively. We impose the continuity of the velocity and stress fields at
the boundary of the core in $r=a$ in the limit of a small core size
$\tilde{a}\ll 1$, in which case we can make use of the following asymptotic
expressions:

$\displaystyle\tilde{v}_{r}^{\mathrm{c}}(\tilde{r})$
$\displaystyle\underset{\tilde{r}\ll
1}{\simeq}\frac{A^{\mathrm{c}}}{8}\tilde{r}^{2}+\frac{\tilde{f}}{3}\tilde{r}$
(75) $\displaystyle\tilde{v}_{r}(\tilde{r})$
$\displaystyle\underset{\tilde{r}\ll
1}{\simeq}-s\frac{C^{0}_{1}}{2}-\tilde{v}_{0}\frac{C^{0}_{2}}{2}+\frac{B_{0}}{\tilde{r}^{2}}+\frac{B}{2}\log\tilde{r}$
(76) $\displaystyle\tilde{v}_{\theta}^{\mathrm{c}}(\tilde{r})$
$\displaystyle\underset{\tilde{r}\ll
1}{\simeq}-\frac{3A^{\mathrm{c}}}{8}\tilde{r}^{2}-\frac{2\tilde{f}}{3}\tilde{r}$
(77) $\displaystyle\tilde{v}_{\theta}(\tilde{r})$
$\displaystyle\underset{\tilde{r}\ll
1}{\simeq}s\frac{C^{0}_{1}}{2}+\tilde{v}_{0}\frac{C^{0}_{2}}{2}+\frac{B_{0}}{\tilde{r}^{2}}-\frac{B}{2}\log\tilde{r}\
,$ (78)

with $C_{1}^{0}$ and $C_{2}^{0}$ given respectively by eqs. (60) and (61) with
$\lambda=0$. Finally, we obtain the stall force by imposing $\tilde{v}_{0}=0$.
The result is given by eq. (16), to leading order in
$\tilde{a}=a/L=a\sqrt{\xi/\eta}$.

###  C.2 Finite rotational viscosity

The procedure to determine the stall force with a finite rotational viscosity
resembles that presented in C.1. Since we assume no nematic order within the
core region, the fields inside the core are unchanged as compared with C.1.
The pressure field outside the core does not depend on the rotational
viscosity and is unchanged. The quantity that changes is the velocity field
outside the core region, now given by eqs. (58)–(61). Close to the core, for
$\bar{a}<\bar{r}\ll 1$, the velocity components have the following asymptotic
expressions:

$\displaystyle\bar{v}_{r}(\bar{r})\underset{\bar{r}\ll 1}{\simeq}$
$\displaystyle-\frac{2^{-\alpha}}{\Gamma(\alpha+1)}(s\,C^{\lambda}_{1}+\tilde{v}_{0}\,C^{\lambda}_{2})\bar{r}^{\lambda/2+\alpha-1}$
$\displaystyle+2^{\alpha-1}B_{0}\Gamma(\alpha)\bar{r}^{\lambda/2-\alpha-1}-\frac{B}{\lambda}$
(79) $\displaystyle\bar{v}_{\theta}(\bar{r})\underset{\bar{r}\ll 1}{\simeq}$
$\displaystyle\frac{2^{-\alpha}(\lambda/2+\alpha)}{\Gamma(\alpha+1)}(s\,C^{\lambda}_{1}+\tilde{v}_{0}\,C^{\lambda}_{2})\bar{r}^{\lambda/2+\alpha-1}$
$\displaystyle+2^{\alpha-1}(\alpha-\lambda/2)B_{0}\Gamma(\alpha)\bar{r}^{\lambda/2-\alpha-1}+\frac{B}{\lambda}\
,$ (80)

with the notations of B.2. Using the expressions eqs. (69), (70), (75), (77),
(79), and (80), continuity at the core boundary leads to the stalling force
given by eq. (30) for $\tilde{v}_{0}=0$111Note that special care must be taken
when matching the fields inside and outside the core, because normalizations
have been done in these two instances using two different characteristic
lengths, $L=\sqrt{\eta/\xi}$ and $\bar{L}=\sqrt{(\eta+\gamma/4)/\xi}$,
respectively..

##  Appendix D Cell division/extrusion

In section 2.3, we introduce the Helmholtz decomposition
$\mathbf{v}=\mathbf{\nabla}\times(\psi\ \mathbf{e}_{z})+\mathbf{\nabla}\phi$.
In dimensionless units, $\tilde{\phi}$ satisfies eq. (21), which is solved by:

$\displaystyle\tilde{\phi}(\tilde{r})$
$\displaystyle=\left\\{I_{1}\left(\frac{\tilde{r}}{\delta}\right)\left(A^{0}_{\phi}-\int_{\tilde{r}/\delta}^{+\infty}{K_{1}(u)(A_{\phi}u^{2}+su+B_{\phi})\>\mathop{}\\!\mathrm{d}u}\right)\right.$
$\displaystyle+\left.K_{1}\left(\frac{\tilde{r}}{\delta}\right)\left(B^{0}_{\phi}-\int_{0}^{\tilde{r}/\delta}{I_{1}(u)(A_{\phi}u^{2}+su+B_{\phi})\>\mathop{}\\!\mathrm{d}u}\right)\right\\}\,,$
(81)

where $\delta=\sqrt{(\eta+\kappa)/\eta}$. The velocity then reads:

$\displaystyle\tilde{v}_{r}(\tilde{r})$
$\displaystyle=\left(\frac{\tilde{\psi}(\tilde{r})}{\tilde{r}}+\frac{\mathop{}\\!\mathrm{d}\phi(\tilde{r})}{\mathop{}\\!\mathrm{d}\tilde{r}}\right)\cos\theta$
(82) $\displaystyle\tilde{v}_{\theta}(\tilde{r})$
$\displaystyle=-\left(\frac{\mathop{}\\!\mathrm{d}\tilde{\psi}(\tilde{r})}{\mathop{}\\!\mathrm{d}\tilde{r}}+\frac{\tilde{\phi}(\tilde{r})}{\tilde{r}}\right)\sin\theta\
,$ (83)

where $\tilde{\psi}$ is given by eq. (B.1) as in B.1. The force-balance
condition (8) then imposes:

$\displaystyle B$ $\displaystyle=\frac{B_{\phi}}{\delta}$ (84) $\displaystyle
A$ $\displaystyle+\frac{A_{\phi}}{\delta}=\tilde{v}_{0}\ .$ (85)

Imposing that the divergence of the velocity field—or equivalently the net
division rate $k$—does not diverge at the core nor at infinity, we find that
$A^{0}_{\phi}$, $A_{0}$, $B^{0}_{\phi}$, $B_{0}$, $B$, and $B_{\phi}$ must all
vanish. At the center of the defect $\tilde{r}=0$, the components of velocity
then read $\tilde{v}_{r}(\tilde{r}=0)=-(\tilde{v}_{0}+s(\pi/4)(1+1/\delta))$
and $\tilde{v}_{\theta}(\tilde{r}=0)=\tilde{v}_{0}+s(\pi/4)(1+1/\delta)$.
Imposing $\tilde{\mathbf{v}}(\tilde{r}=0)=\mathbf{0}$ in the reference frame
of the default yields the self-advection velocity given by eq. (22).

The divergence of the velocity field is then given by:

$\displaystyle\nabla\cdot\tilde{\mathbf{v}}=-\frac{s}{\delta^{2}}$
$\displaystyle\left[I_{1}\left(\frac{\tilde{r}}{\delta}\right)\int_{\tilde{r}/\delta}^{+\infty}{K_{1}(u)\,u\>\mathop{}\\!\mathrm{d}u}\right.$
$\displaystyle\left.+K_{1}\left(\frac{\tilde{r}}{\delta}\right)\int_{0}^{\tilde{r}/\delta}{I_{1}(u)\,u\>\mathop{}\\!\mathrm{d}u}-1\right]\cos\theta\
.$ (86)

It is represented in fig. 3.

##  Appendix E -1/2 defects

###  E.1 Velocity field

This appendix is dedicated to the computation of the velocity field around a
$-1/2$ defect, plotted on fig. 2. We add the superscript ‘-’ to denote the
quantities associated to a $-1/2$ defect. Force balance reads:

$\eta\Delta\mathbf{v}^{-}-\mathbf{\nabla}P^{-}-\xi\mathbf{v}^{-}-\frac{\zeta\Delta\mu}{2r}\left(-\cos(3\theta)\>\mathbf{e}_{r}+\sin(3\theta)\>\mathbf{e}_{\theta}\right)=0\
.$ (87)

Here the defect is immotile, such that the velocity field satisfies:

$\displaystyle\mathbf{v}^{-}(+\infty,\theta)$ $\displaystyle=0$ (88)
$\displaystyle\mathbf{v}^{-}(0,\theta)$ $\displaystyle=0\ .$ (89)

The curl of eq. (87) gives:

$\Delta\left[\eta\Delta\psi^{-}-\xi\psi^{-}\right]=\frac{3}{2}\zeta\Delta\mu\frac{\sin
3\theta}{r^{2}}\ .$ (90)

Using a normalization by a characteristic length $L=\sqrt{\eta/\xi}$ and a
characteristic time $\tau^{-}=2\eta/(3|\zeta|\Delta\mu)$, we introduce the
dimensionless stream function $\tilde{\psi}^{-}=(\tau^{-}/L^{2})\psi^{-}$ and
spatial variable $\tilde{r}=r/L$. We then have:

$\Delta\left[\Delta\tilde{\psi}^{-}-\tilde{\psi}^{-}\right]=s\frac{\sin
3\theta}{r^{2}}\ ,$ (91)

where $s=\textrm{sign}(\zeta)$ and $\psi^{-}$ is of the form
$\tilde{\psi}^{-}(\tilde{r},\theta)=\tilde{\psi}^{-}(\tilde{r})\sin 3\theta$.
Integrating one Laplace operator in this equation leads to:

$\frac{\mathop{}\\!\mathrm{d}^{2}\tilde{\psi}^{-}(\tilde{r})}{\mathop{}\\!\mathrm{d}\tilde{r}^{2}}+\frac{1}{\tilde{r}}\frac{\mathop{}\\!\mathrm{d}\tilde{\psi}^{-}(\tilde{r})}{\mathop{}\\!\mathrm{d}\tilde{r}}-\left(1+\frac{3}{\tilde{r}^{2}}\right)\tilde{\psi}^{-}(\tilde{r})=A^{-}\tilde{r}+\frac{B^{-}}{\tilde{r}}-s\
.$ (92)

Solving this equation leads to the following dimensionless velocity field:

$\displaystyle\tilde{v}_{r}^{-}=s$
$\displaystyle\left\\{I_{3}(\tilde{r})\int_{\tilde{r}}^{+\infty}{K_{3}(u)u\>\mathrm{d}u}\right.$
$\displaystyle+\left.K_{3}(\tilde{r})\int_{0}^{\tilde{r}}{I_{3}(u)u\>\mathrm{d}u}\right\\}\cos
3\theta$ (93) $\displaystyle\tilde{v}_{\theta}^{-}=s$
$\displaystyle\left\\{\left(\frac{3}{\tilde{r}}I_{3}(\tilde{r})-I_{2}(\tilde{r})\right)\int_{\tilde{r}}^{+\infty}{K_{3}(u)u\>\mathrm{d}u}\right.$
$\displaystyle+\left.\left(\frac{3}{\tilde{r}}K_{3}(\tilde{r})+K_{2}(\tilde{r})\right)\int_{0}^{\tilde{r}}{I_{3}(u)u\>\mathrm{d}u}\right\\}\sin
3\theta\ ,$ (94)

accounting for the boundary conditions eqs. (88) and (89). This velocity field
is plotted on fig. 2b,d.

The pressure field $P^{-}$ is obtained by taking the divergence of eq. (87):

$\Delta P^{-}=\frac{3\zeta\Delta\mu}{2r^{2}}\cos 3\theta\ ,$ (95)

which yields, with the boundary conditions (88)–(89):

$P^{-}(r,\theta)=-\frac{3\zeta\Delta\mu}{2}\cos 3\theta\ .$ (96)

###  E.2 Cell division/extrusion

We derive in this section the divergence of the velocity field represented on
fig. 3. The derivation resembles that for +1/2 defects as presented in section
2.3 with a pressure-dependent division rate given by Eq. (19). The velocity
field is decomposed into $\mathbf{v}^{-}=\mathbf{\nabla}\times(\psi^{-}\
\mathbf{e}_{z})+\mathbf{\nabla}\phi^{-}$. The divergence-free and curl-free
parts of the velocity respectively satisfy:

$\displaystyle\Delta\left[\Delta\tilde{\psi}^{-}-\tilde{\psi}^{-}\right]$
$\displaystyle=s\frac{\sin 3\theta}{\tilde{r}^{2}}$ (97)
$\displaystyle\Delta\left[\Delta\tilde{\phi}^{-}-\frac{1}{\delta^{2}}\tilde{\phi}^{-}\right]$
$\displaystyle=s\frac{\cos{3\theta}}{\delta^{2}\tilde{r}^{2}}\ .$ (98)

Solving eq. (98), we get the divergence profile of the velocity field
represented on fig. 3b,d as:

$\displaystyle\nabla\cdot\tilde{\mathbf{v}}^{-}=\frac{s}{\delta^{2}}$
$\displaystyle\left[I_{3}\left(\frac{\tilde{r}}{\delta}\right)\int_{\tilde{r}/\delta}^{+\infty}{K_{3}(u)\,u\>\mathop{}\\!\mathrm{d}u}\right.$
$\displaystyle\left.+K_{3}\left(\frac{\tilde{r}}{\delta}\right)\int_{0}^{\tilde{r}/\delta}{I_{3}(u)\,u\>\mathop{}\\!\mathrm{d}u}-1\right]\cos
3\theta\ .$ (99)