# Dirac-vortex topological photonic crystal fibre

Hao Lin  Institute of Physics, Chinese Academy of Sciences/Beijing National
Laboratory for Condensed Matter Physics, Beijing 100190  School of Physical
Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
Ling Lu<EMAIL_ADDRESS>Institute of Physics, Chinese Academy of
Sciences/Beijing National Laboratory for Condensed Matter Physics, Beijing
100190  Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China

###### Abstract

The success of photonic crystal fibres relies largely on the endless variety
of two-dimensional photonic crystals in the cross-section. Here, we propose a
topological bandgap fibre whose bandgaps along in-plane directions are opened
by generalized $\rm{Kekul\acute{e}}$ modulation of a Dirac lattice with a
vortex phase. Then, the existence of mid-gap defect modes is guaranteed to
guide light at the core of this Dirac-vortex fibre, where the number of
guiding modes equals the winding number of the spatial vortex. The single-
vortex design provides a single-polarization single mode for a bandwidth as
large as one octave.

##  I Introduction

Topological photonics Lu  _et al._ (2014); Khanikaev and Shvets (2017); Ozawa
_et al._ (2019), initiated with the idea of robust waveguiding, is inspiring
novel fibre concepts, such as a one-way fibre inside a magnetic three-
dimensional photonic crystal Lu  _et al._ (2018) and a Bragg fibre with
nontrivial edge modes Pilozzi  _et al._ (2020). In this article, we introduce
the topological photonic crystal fibre (PCF) whose invariant cross-section
resembles the recent Dirac-vortex topological cavity Gao  _et al._ (2020) in
two-dimensional photonic crystals. Such a topological bound state can be
traced back to the Jackiw-Rossi zero mode in the 2D Dirac equation Jackiw and
Rossi (1981), and has been realized in honeycomb lattices Hou  _et al._
(2007); Iadecola  _et al._ (2016) in a couple of systems Menssen  _et al._
(2020); Noh  _et al._ (2020); Gao  _et al._ (2019); Chen  _et al._ (2019).
This Dirac-vortex silica fibre can support an _arbitrary_ number of nearly
degenerate guiding modes by varying the winding number ($w$) of the spatial
vortex. When $w=\pm 1$, the fibre can support a single-polarization single
mode (SPSM) with a large bandwidth.

The SPSM fibre supports truly one mode, while traditional single-mode fibres
and polarization-maintaining fibres both support two polarizations, either
degenerate or nondegenerate, respectively. Such fibre birefringence (dual
polarization) broadens the optical pulses being transmitted, known as
polarization-mode dispersion. To solve this limitation, SPSM fibres have been
designed to separate the degenerate cutoff frequencies by lowering the
symmetry of the fibre cross-section, which can be achieved by structural
asymmetry or nonuniform stress. This dominant asymmetric approach is applied
mostly to the lowest-frequency index-guided fundamental modes (polarization-
degenerate) Okoshi and Oyamada (1980); Eickhoff (1982); Simpson  _et al._
(1983); Kubota  _et al._ (2004); Folkenberg  _et al._ (2005); Lee  _et al._
(2008), but also works for the nonfundamental degenerate modes inside a
bandgap Ferrando and Miret (2001); Eguchi and Tsuji (2012); Szpulak  _et al._
(2007).

The SPSM bandwidth is very limited due to the amount of asymmetry that can be
applied; the best experimental value is approximately 30% (frequency span over
central frequency), which is achieved in a stressed PCF Folkenberg  _et al._
(2005). (In integrated waveguides, an SPSM bandwidth over one octave [66.67%]
has been demonstrated Chiles and Fathpour (2016).) For example, in Fig. 1(a),
the SPSM bandwidth is bounded from above by the cutoff frequency of the other
polarization (blue) and bounded from below by the confinement loss, even if
this guiding mode (red line) has no cutoff frequency. However, there exists a
symmetric approach to an SPSM by operating at a singly degenerate mode inside
a bandgap so that the fibre cross-section can remain highly symmetric. Shown
in Fig. 1(b), this symmetric approach has thus far only been proposed
theoretically in a hollow-core Bragg fibre Bassett and Argyros (2002); Argyros
_et al._ (2004), in which many other guided modes (blue line) can be made much
lossier than the TE01 mode (red line). The SPSM bandwidth of this design is
even more limited.

The Dirac-vortex fibre is an ideal design for an ultrabroadband SPSM by
ensuring singlet mid-gap dispersion inside the bandgap (the symmetric
approach), as illustrated in Fig. 1(c). Without the topological mechanism of
the Dirac vortex, it is generally difficult to stabilize a defect mode at the
middle of the bandgap for every wave vector.

Figure 1:  Asymmetric (a) and symmetric (b, c) approaches to an SPSM. (a)
Dominant approach to split the degenerate fundamental modes by structural
asymmetry. (b) Proposed Bragg fibre design to guide light by a singly
degenerate mode with lower loss than other modes. (c) The topological PCF
provides singlet mid-gap dispersion for a broadband SPSM.

In the rest of the article, we first show how to gap the nodal line without
leaving residual Weyl points. Such a supercell ${\rm{Kekul\acute{e}}}$
modulation Hou  _et al._ (2007); Gao  _et al._ (2020) is generalized to
continuous $2\pi$ phase angles that are used to construct a vortex gap around
the fibre core that can confine any number of fibre modes. To ease the
fabrication, a simplified design of only four capillary silica tubes is
introduced. Finally, we enlarge the vortex size to eliminate the index-guided
modes and achieve an octave-spanning SPSM.

Figure 2: Band diagram of silica ($\varepsilon=2.1$) photonic crystal lattices
uniform in the out-of-plane direction ($z$). (a) Projected band diagram of the
triangular photonic crystal, in which the nodal line degeneracy is
highlighted. (b) An extra air hole in the primitive cell breaks the inversion
symmetry, and the nodal line is lifted into Weyl points. Insets: cross-section
structures and in-plane band structures at $k_{z}a/2\pi=2.02$. Two different
primitive-cell choices are drawn in (a). Figure 3: Dirac-vortex fibre obtained
by continuous Kekul$\rm\acute{e}$ modulations. (a) Example of how an “atom”
can be shifted in any direction ($\arg[\boldsymbol{\delta}]=\phi$) with finite
amplitude ($|\boldsymbol{\delta}|$) by changing the widths of the three
struts. (b1), (c1), (d1) Supercell examples of three coordinated atoms (A1,
A2, A3) with $\left|\boldsymbol{\delta}\right|=\sqrt{3}a/80$,
$\phi=120^{\circ}$; $|\boldsymbol{\delta}|=0$; and
$\left|\boldsymbol{\delta}\right|=\sqrt{3}a/80$, $\phi=60^{\circ}$,
respectively. $\Delta t=2\sqrt{3}|\boldsymbol{\delta}|t_{0}/a=0.012a$. The
corresponding band structures are plotted in (b2), (c2), (d2), respectively.
(e) Structure of a continuous Dirac-vortex PCF, in which every strut is
coloured according to its width. (f) Band diagram of the fibre plotted in
reference to the frequency of the original nodal line (central dashed line).
The inset shows the intensity patterns ($\hat{\boldsymbol{z}}\cdot\rm
Re[\boldsymbol{E}^{*}\times\boldsymbol{H}]$) of the topological mode and one
local-defect mode. The single-polarization topological mode (red line) spans
over two octaves.  Figure 4: Continuous Dirac-vortex fibres with winding
number $w=+1,+2,+3$. (a1), (b1), and (c1) are the fibre structures. The colour
wheels show the corresponding phases of the generalized $\rm{Kekul\acute{e}}$
modulation. (a2), (b2), and (c2) are the corresponding band diagrams, and the
insets show the mode profiles ($\hat{\boldsymbol{z}}\cdot\rm
Re[\boldsymbol{E}^{*}\times\boldsymbol{H}]$) of modes at $k_{z}a/2\pi=4$.

##  II Results

###  II.1  Nodal lines and Weyl points in a PCF

We start with the most common PCF structure Knight (2003); Russell (2006), a
silica photonic crystal with a triangular lattice of air holes, shown in Fig.
2(a). There are two nodal lines of 2D Dirac points Xie  _et al._ (2015);
Biswas  _et al._ (2016), at the $\pm K$ points in the Brillouin zone. The
Dirac points are frequency isolated in the in-plane 2D band structure for
$k_{z}a/2\pi>1.6$. Each nodal line gaps into Weyl points Lu  _et al._ (2013),
if we break the inversion symmetry by adding an extra small air hole in the
primitive cell, as shown in Fig. 2(b). Since this PCF structure has mirror
symmetry along $z$, the Weyl surface state does not have handedness
propagating along the $z$ axis, unlike that observed in coupled spiral
waveguide arrays Rechtsman  _et al._ (2013); Noh  _et al._ (2017) containing
similar type-II Weyl points. To gap the entire nodal line, without leaving
Weyl degeneracies, we apply supercell modulations.

###  II.2 Generalized Kekul$\mathbf{\acute{e}}$ modulation

As shown in Fig. 3(a-d), we perturb the nodal-line lattice with the
generalized $\rm{Kekul\acute{e}}$ pattern discussed in Ref. Gao  _et al._
(2020). The idea is to couple the two nodal lines (of Dirac points) together
in an enlarged supercell and annihilate them into a bandgap through supercell
modulation. Since each supercell has three primitive cells, we label each
primitive cell [Fig. 2(a) inset] as an artificial “atom” consisting of three
struts [Fig. 3(a)]. In practice, we shift the three “atoms” in the supercell
with identical amplitude and constant phase difference (120∘) between each two
“atoms”. To move each “atom”, we can shift its centre of mass
$\boldsymbol{\delta}(\phi)$, in any direction $\phi$, by adjusting the
thickness of the three struts ($t_{1}$, $t_{2}$, $t_{3}$) without changing the
total mass of the “atom” ($\propto 3t_{0}$). Therefore, the strut thickness
for any modulation vector $\boldsymbol{\delta}(\phi)$ can be determined by the
following equation:

$\sum_{i=1}^{3}t_{i}\boldsymbol{R}_{i}=\boldsymbol{\delta}(\phi)\sum_{i=1}^{3}t_{i}=\boldsymbol{\delta}(\phi)3t_{0}$
(1)

, where $\boldsymbol{R}_{i}$ ($|\boldsymbol{R}_{i}|=\sqrt{3}/6a$) is the
position vector of the centre of mass of each strut. Examples of the lattices
before and after the $\rm{Kekul\acute{e}}$ modulations are drawn in Fig.
3(b1-d1), and their corresponding band structures are plotted in (b2-d2).

After the generalized modulation ($|\boldsymbol{\delta}(\phi)|\neq 0$), the
entire nodal lines gap out for arbitrary $\phi\in(0,2\pi)$, and the bandgaps
share a common frequency range, which is eventually determined by the band
edges of the two most symmetric supercells of $\phi$=60∘ and 120∘. The
lattices of the rest of the modulation phases can be understood as
interpolations between the two. Note that the modulation phase $\phi$ of each
supercell is labelled by the shift angle of “atom” A1 in Fig. 3(c1).

###  II.3 Continuous modulation

As shown in Fig. 3(e), we arrange the series of
$\rm{Kekul\acute{e}}$-modulated 2D lattices angularly around a chosen core
point in the fibre cross-section as a function of their modulation phase
$\phi$.As a result, a single defect mode appears in the middle of the bandgap
and is spatially localized at the fibre core, as shown in Fig. 3(f). This
result is expected from the Dirac-vortex cavity result Gao  _et al._ (2020)
because a PCF can be regarded as a 2D cavity for every propagation constant
($k_{z}$). Since the fibre is $C_{3v}$ symmetric, one-sixth of the structure
is sufficient to study this topological mode.

The large momentum [$k_{z}a/(2\pi)>5$], or short wavelength, behaviour of the
fibre band diagram in Fig. 3(f) is worth discussion. First, the common bandgap
frequency is lower than that of the nodal line (central dashed line in the
figure) at short wavelengths because the short-wavelength modes can localize
preferentially at struts that are thicker than $t_{0}$ (rather than at those
that are thinner) in the modulated lattice, resulting in lowering of the
overall band frequencies. Second, the topological mode is very sensitive to
the central three struts, where the modal intensity peaks. In this case, the
three central red struts are the thickest, and consequently, the fibre mode
cuts off, in the short-wavelength region, at a frequency lower than the
bandgap. Third, although we “continuously” vary the modulation angle $\phi$ in
the surrounding lattice, the discretization is still limited to a single strut
in our design. A consequence of the strut-discrete vortex is that a set of
extra modes is trapped very close to the spectral edges of the bandgap. We
call them local-defect modes, shaded in grey colour in the figure. These modes
tend to arise at short wavelengths and are more prominent for large modulation
amplitudes (large bandgaps), as shown in the example in Fig. 3(f). Since these
local-defect modes are located near the bulk band edges with an extended mode
profile, their confinement losses are at least one order of magnitude higher
than that of the guided topological mode with $w=+1$.

Figure 5: Dirac-vortex fibre obtained by discrete Kekul$\rm\acute{e}$
modulations. (a) Fibre preform with four capillary silica tubes of the same
outer diameter. (b) Fibre structure, having seven different strut widths,
resulting from the preform in (a). The black rectangle is the unit cell for
the line-defect modes. (c) Band diagram with the frequency referenced to the
original nodal line of the unmodulated lattice. The higher-order modes are
plotted with blue and green lines according their degeneracy. Insets show the
mode profiles ($\hat{\boldsymbol{z}}\cdot\rm
Re[\boldsymbol{E}^{*}\times\boldsymbol{H}]$). (d) Purple dispersion and mode
profiles of the polarization-degenerate index-guided fundamental mode. Other
index-guide modes are plotted with green dashes.

###  II.4 Arbitrary degenerate modes

One main topological feature of the Dirac-vortex fibre is the ease of creating
multiple near-degenerate modes by simply increasing the winding number ($w$)
of the vortex Gao  _et al._ (2020). We demonstrate this in Fig. 4 for
$w=+1,+2,+3$, where $w$ is the topological invariant of the system and can be
an arbitrary integer. The sign of $w$ determines where the field peaks around
one of the two sublattices of a honeycomb lattice. The two sublattices can be
viewed as the two joints, each having three struts, in the primitive cell of
Fig. 2(a). An example of $w=-1$ is presented in Supplementary Part I.

In these multimode examples, we decrease the modulation amplitude (as well as
the bandgap size) to eliminate the local-defect modes in the momentum range
plotted. Additionally, we keep the core centre at the $C_{3v}$ centre of the
$w=+1$ vortex. If we choose a w-dependent centre, then any vortex can be
$C_{3v}$ symmetric, so some pairs of fibre dispersions will be rigorously
degenerate due to the doublet representation of the point group.

In principle, a continuous, single-mode or multimode, Dirac-vortex PCF can be
fabricated either from 3D-printed preforms or by the stack-and-draw method
with over a hundred tubes with different tube thicknesses. Neither of the two
solutions are convenient. Therefore, we present the discrete version of the
fibre design.

###  II.5 Discrete modulation with four tubes

We discretize the continuous vortex design [Fig.3(e)] into six discrete
modulation phase angles ($\phi$). Due to the $C_{3v}$ symmetry of the cross-
section, only two phases of the six have distinct structures. The three phases
that differ by 120∘ are actually identical in the unit-cell structure (but
with different origins and orientations). As a result, we only need four tubes
to stack-and-draw the Dirac-vortex PCF, which is very reasonable for
fabrication. Another discrete version is presented in Supplementary Part II
using four strut thicknesses.

As shown in Fig. 5(a), the four silica capillary tubes have the same outer
diameter $d_{\rm{outer}}$ to maintain the lattice but different inner
diameters $d_{i}$ for modulation. The four exact ratios of
$d_{i}/d_{\rm{outer}}$ can be determined from Eq. 2, satisfying the thickness
and area proportionality.

$\begin{cases}\frac{d_{\rm{outer}}-(d_{1,4}+d_{2,3})/2}{t_{2,4}}&=\frac{d_{\rm{outer}}-d_{2,3}}{t_{1,3}}\\\
\frac{\pi(3d_{\rm{outer}}^{2}-d_{1,4}^{2}-2d_{2,3}^{2})}{4d_{\rm{outer}}^{2}}&=\frac{3\sqrt{3}t_{0}(a-t_{0}/2)}{a^{2}}\end{cases}$
(2)

The resulting discrete Dirac-vortex fibre and its band structure are plotted
in Fig. 5(b) and (c). Compared to the continuous version in Fig. 3(e) and (f),
the structural nonuniformity now only exists at the six identical interfaces
(also determined by Eq. 2) between the two distinct lattices Wu and Hu (2015),
at which the extra defect modes can locate. In addition to these line-defect
modes inside the bandgap, there are also higher-order vortex modes plotted in
green and blue lines, according to their $C_{3v}$ point group representations.
The mode profiles of the higher-order modes are much larger than that of the
topological mode.

We have been focusing on the modes inside the bandgap. However, as shown in
Fig. 5(d) [also in Fig. 3(e, f)], there could also be index-guided modes in
the Dirac-vortex fibre whose frequency is the lowest for a common wave vector.
Index-guided modes exist wherever there is a sharp local maximum of the strut
thickness, equivalent to a local rise of the effective refractive index. The
loss of the index-guided modes is usually much lower than that of the in-gap
modes.

###  II.6 Nonzero vortex size with continuous modulation

To operate at the topological mid-gap dispersion, we have to remove the index-
guided modes in the design. We do this by enlarging the vortex size and
smoothing the sharp lattice (strut) interfaces at the fibre core. The smooth
modulation envelope function is $\Delta t(r)=\rm{tanh}(r/R)$, where $r$ is the
radial distance and $R$ is the vortex radius of choice, similar to the 2D
cavity design in Ref. Gao  _et al._ (2020). Note that in the previous examples
of this paper, the vortex size is $R=0$.

The fibre cross-section with vortex size $R=3a$ is shown in Fig. 6(a), and the
corresponding band structure is shown in Fig. 6(b) and (c). There are no
index-guided modes at the bottom of the bands. For completeness, we also plot
the blue mid-gap dispersion inside the second topological bandgap, originating
from a higher frequency Dirac point (nodal line) illustrated in Supplementary
Part III. We did not discuss it in the previous examples because it has a much
higher loss than the red mid-gap modes in the first topological gap, as shown
in Fig. 6(d).

Figure 6: Octave SPSM in a continuous Dirac-vortex fibre with a finite vortex
size. (a) Fibre structure with sixteen cladding periods in radius. The colour
wheel represents the phase and amplitude of the generalized
$\rm{Kekul\acute{e}}$ modulation. (b) Complete fibre dispersion in absolute
frequency. First topological mode (red line) and higher-order doublet modes
(green dotted line) in the first topological bandgap, as well as a second
topological mode (blue line) in the second topological bandgap at higher
frequency. The mode profiles ($\hat{\boldsymbol{z}}\cdot\rm
Re[\boldsymbol{E}^{*}\times\boldsymbol{H}]$) of the two topological modes are
shown in the insets, circled with different colours for clarity. (c) Fibre
dispersion in frequency relative to the original nodal-line frequency. (d)
Confinement losses of the guided modes. (e) Dispersion parameter and effective
area of the first topological mode. (f) Bending loss of the first topological
mode at $\lambda=1550$ $nm$.

###  II.7 Octave SPSM

The continuously modulated Dirac-vortex PCF with a finite vortex diameter
($2R=6a$) has an SPSM. We evaluate its potential performance in terms of the
confinement loss, dispersion parameter, effective area and bending loss in
Fig. 6(d,e,f). We choose the wave vector ($k_{z}a/2\pi$) range from 1.8 to 5.7
and the period $a=3$ $\mu m$, so that the corresponding wavelength ranges from
$1700$ $nm$ to $700$ $nm$ — the low-loss window of silica.

The modes with the lowest confinement losses are plotted in Fig. 6(d),
computed with absorbing boundary conditions. The loss of the topological mode
(red) is always the lowest for the whole wavelength range over one octave.
Most of its confinement loss is less than 1 $dB/km$ with sixteen layers of air
holes along the fibre radius [Fig. 6(a)].

The dispersion parameter $D$ of the topological mode can take large negative
values of -1500 $ps/nm/km$, and the zero-dispersion wavelength is
approximately 1.3 $\mu m$. The effective area ($A_{\rm{eff}}$ defined in Ref.
Joannopoulos  _et al._ (2008)) is 23 $\mu m^{2}$ at a 0.7 $\mu m$ wavelength
and increases to 243 $\mu m^{2}$ at a 1.7 $\mu m$ wavelength. The bending loss
is calculated using a conformal transformation analysis Heiblum and Harris
(1975) with an equivalent refractive index profile $n_{eq}(x,y)=n(x,y){\rm
exp}(x/X)$ Beravat  _et al._ (2016) for the bending radius $X$. The results in
Fig. 6(f) show that the fibre can be bent to a radius as small as 15 $mm$
without experiencing significant losses. These specifications of the Dirac-
vortex PCF are similar to those of previous PCFs Finazzi  _et al._ (2003);
Tsuchida  _et al._ (2005), with the key difference being that the topological
mode is a single-polarization mode.

##  III Discussion

We numerically investigate the Dirac-vortex topological PCF in terms of its
principle, construction, and potential performance. This fibre could be made
with the standard stack-and-draw process with silica glass tubes or 3D-printed
preforms. The design is tolerant of structural details such as the interface
curvatures, not considered in the main text but discussed in Supplementary
Part IV. It is also possible to move the topological dispersion above the
light line, as discussed in Supplementary Part V.

Similar to the case of a three-dimensional one-way fibre Lu  _et al._ (2018),
where the topological invariant is a four-dimensional second Chern number, the
topological invariant of the Dirac-vortex PCF can be written as a three-
dimensional integral (equation 4.2 in Ref. Teo and Kane (2010) under chiral
symmetry) where the coordinates are $k_{x},k_{y}$ and the polar angle (a form
of synthetic dimension Yuan  _et al._ (2018)).

The advantage of the Dirac-vortex PCF over previous fibres is the ability to
guide any number of near-degenerate modes at will. The single-mode design
provides an SPSM fibre with an octave bandwidth. The effective mode area is
easily tuneable by changing the vortex size ($R$). This work suggests PCFs as
a new platform for topological photonics.

##  IV Acknowledgements

We thank Wei Ding, YingYing Wang, Xiaomei Gao, Zhong Wang and Guoqing Chang
for helpful discussions. This work was supported by the National Key R&D
Program of China (2017YFA0303800, 2016YFA0302400), the Natural Science
Foundation of China (11721404), the Strategic Priority Research Program
(XDB33000000) and the international partnership program (112111KYSB20200024)
of the Chinese Academy of Sciences, and the Beijing Natural Science Foundation
(Z200008).

##  V Author contributions

Both authors contributed to the design of the fibre and the writing of the
manuscript. H.L. performed the numerical simulations.

##  VI Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * Lu  _et al._ (2014) Ling Lu, John D Joannopoulos, and Marin Soljačić, “Topological photonics,” Nature Photonics 8, 821–829 (2014).
  * Khanikaev and Shvets (2017) Alexander B Khanikaev and Gennady Shvets, “Two-dimensional topological photonics,” Nature Photonics 11, 763 (2017).
  * Ozawa  _et al._ (2019) Tomoki Ozawa, Hannah M Price, Alberto Amo, Nathan Goldman, Mohammad Hafezi, Ling Lu, Mikael C Rechtsman, David Schuster, Jonathan Simon, Oded Zilberberg, _et al._ , “Topological photonics,” Reviews of Modern Physics 91, 015006 (2019).
  * Lu  _et al._ (2018) Ling Lu, Haozhe Gao, and Zhong Wang, “Topological one-way fiber of second chern number,” Nature communications 9, 1–7 (2018).
  * Pilozzi  _et al._ (2020) L Pilozzi, Daniel Leykam, Zhigang Chen, and Claudio Conti, “Topological photonic crystal fibers and ring resonators.” Optics Letters 45, 1415–1418 (2020).
  * Gao  _et al._ (2020) Xiaomei Gao, Lechen Yang, Hao Lin, Lang Zhang, Jiafang Li, Fang Bo, Zhong Wang, and Ling Lu, “Dirac-vortex topological cavities,” Nature Nanotechnology , 1–7 (2020).
  * Jackiw and Rossi (1981) R Jackiw and Paolo Rossi, “Zero modes of the vortex-fermion system,” Nuclear Physics B 190, 681–691 (1981).
  * Hou  _et al._ (2007) Chang-Yu Hou, Claudio Chamon, and Christopher Mudry, “Electron fractionalization in two-dimensional graphenelike structures,” Physical review letters 98, 186809 (2007).
  * Iadecola  _et al._ (2016) Thomas Iadecola, Thomas Schuster, and Claudio Chamon, “Non-abelian braiding of light,” Physical Review Letters 117, 073901 (2016).
  * Menssen  _et al._ (2020) Adrian J Menssen, Jun Guan, David Felce, Martin J Booth, and Ian A Walmsley, “Photonic topological mode bound to a vortex,” Physical Review Letters 125, 117401 (2020).
  * Noh  _et al._ (2020) Jiho Noh, Thomas Schuster, Thomas Iadecola, Sheng Huang, Mohan Wang, Kevin P Chen, Claudio Chamon, and Mikael C Rechtsman, “Braiding photonic topological zero modes,” Nature Physics 16, 989–993 (2020).
  * Gao  _et al._ (2019) Penglin Gao, Daniel Torrent, Francisco Cervera, Pablo San-Jose, José Sánchez-Dehesa, and Johan Christensen, “Majorana-like zero modes in kekulé distorted sonic lattices,” Physical review letters 123, 196601 (2019).
  * Chen  _et al._ (2019) Chun-Wei Chen, Natalia Lera, Rajesh Chaunsali, Daniel Torrent, Jose Vicente Alvarez, Jinkyu Yang, Pablo San-Jose, and Johan Christensen, “Mechanical analogue of a majorana bound state,” Advanced Materials 31, 1904386 (2019).
  * Okoshi and Oyamada (1980) T Okoshi and K Oyamada, “Single-polarisation single-mode optical fibre with refractive-index pits on both sides of core,” Electronics Letters 16, 712–713 (1980).
  * Eickhoff (1982) W Eickhoff, “Stress-induced single-polarization single-mode fiber,” Optics letters 7, 629–631 (1982).
  * Simpson  _et al._ (1983) J Simpson, R Stolen, F Sears, William Pleibel, J MacChesney, and Richard Howard, “A single-polarization fiber,” Journal of Lightwave Technology 1, 370–374 (1983).
  * Kubota  _et al._ (2004) Hiirokazu Kubota, Satoki Kawanishi, Shigeki Koyanagi, Masatoshi Tanaka, and Shyunichiro Yamaguchi, “Absolutely single polarization photonic crystal fiber,” IEEE Photonics Technology Letters 16, 182–184 (2004).
  * Folkenberg  _et al._ (2005) JR Folkenberg, MD Nielsen, and C Jakobsen, “Broadband single-polarization photonic crystal fiber,” Optics letters 30, 1446–1448 (2005).
  * Lee  _et al._ (2008) Karen KY Lee, Yehuda Avniel, and Steven G Johnson, “Design strategies and rigorous conditions for single-polarization single-mode waveguides,” Optics express 16, 15170–15184 (2008).
  * Ferrando and Miret (2001) Albert Ferrando and Juan José Miret, “Single-polarization single-mode intraband guidance in supersquare photonic crystals fibers,” Applied Physics Letters 78, 3184–3186 (2001).
  * Eguchi and Tsuji (2012) Masashi Eguchi and Yasuhide Tsuji, “Single-polarization elliptical-hole lattice core photonic-bandgap fiber,” Journal of lightwave technology 31, 177–182 (2012).
  * Szpulak  _et al._ (2007) M Szpulak, T Martynkien, J Olszewski, W Urbanczyk, T Nasilowski, Francis Berghmans, and H Thienpont, “Single-polarization single-mode photonic band gap fiber,” ACTA PHYSICA POLONICA SERIES A 111, 239 (2007).
  * Chiles and Fathpour (2016) Jeff Chiles and Sasan Fathpour, “Demonstration of ultra-broadband single-mode and single-polarization operation in t-guides,” Optics Letters 41, 3836–3839 (2016).
  * Bassett and Argyros (2002) Ian M Bassett and Alexander Argyros, “Elimination of polarization degeneracy in round waveguides,” Optics Express 10, 1342–1346 (2002).
  * Argyros  _et al._ (2004) A Argyros, N Issa, I Bassett, and MA Van Eijkelenborg, “Microstructured optical fiber for single-polarization air guidance,” Optics letters 29, 20–22 (2004).
  * Knight (2003) Jonathan C Knight, “Photonic crystal fibres,” nature 424, 847–851 (2003).
  * Russell (2006) Philip St J Russell, “Photonic-crystal fibers,” Journal of lightwave technology 24, 4729–4749 (2006).
  * Xie  _et al._ (2015) Kang Xie, Wei Zhang, Allan D Boardman, Haiming Jiang, Zhijia Hu, Yong Liu, Ming Xie, Qiuping Mao, Lei Hu, Qian Li, _et al._ , “Fiber guiding at the dirac frequency beyond photonic bandgaps,” Light: Science & Applications 4, e304–e304 (2015).
  * Biswas  _et al._ (2016) Tushar Biswas, Rik Chattopadhyay, and Shyamal K Bhadra, “Dirac-mode guidance in silica-based hollow-core photonic crystal fiber with high-index dielectric rings,” physica status solidi (b) 253, 1898–1906 (2016).
  * Lu  _et al._ (2013) Ling Lu, Liang Fu, John D Joannopoulos, and Marin Soljačić, “Weyl points and line nodes in gyroid photonic crystals,” Nature photonics 7, 294 (2013).
  * Rechtsman  _et al._ (2013) Mikael C Rechtsman, Julia M Zeuner, Yonatan Plotnik, Yaakov Lumer, Daniel Podolsky, Felix Dreisow, Stefan Nolte, Mordechai Segev, and Alexander Szameit, “Photonic floquet topological insulators,” Nature 496, 196–200 (2013).
  * Noh  _et al._ (2017) Jiho Noh, Sheng Huang, Daniel Leykam, Yi Dong Chong, Kevin P Chen, and Mikael C Rechtsman, “Experimental observation of optical weyl points and fermi arc-like surface states,” Nature Physics 13, 611–617 (2017).
  * Wu and Hu (2015) Long-Hua Wu and Xiao Hu, “Scheme for achieving a topological photonic crystal by using dielectric material,” Physical Review Letters 114, 223901 (2015).
  * Joannopoulos  _et al._ (2008) John D Joannopoulos, Steven G Johnson, Joshua N Winn, and Robert D Meade, “Molding the flow of light,” Princeton Univ. Press, Princeton, NJ [ua] , 167 (2008).
  * Heiblum and Harris (1975) Mordehai Heiblum and Jay Harris, “Analysis of curved optical waveguides by conformal transformation,” IEEE Journal of Quantum Electronics 11, 75–83 (1975).
  * Beravat  _et al._ (2016) Ramin Beravat, Gordon KL Wong, Michael H Frosz, Xiao Ming Xi, and Philip St J Russell, “Twist-induced guidance in coreless photonic crystal fiber: A helical channel for light,” Science Advances 2, e1601421 (2016).
  * Finazzi  _et al._ (2003) Vittoria Finazzi, Tanya M Monro, and David J Richardson, “Small-core silica holey fibers: nonlinearity and confinement loss trade-offs,” JOSA B 20, 1427–1436 (2003).
  * Tsuchida  _et al._ (2005) Yukihiro Tsuchida, Kunimasa Saitoh, and Masanori Koshiba, “Design and characterization of single-mode holey fibers with low bending losses,” Optics express 13, 4770–4779 (2005).
  * Teo and Kane (2010) Jeffrey CY Teo and Charles L Kane, “Topological defects and gapless modes in insulators and superconductors,” Physical Review B 82, 115120 (2010).
  * Yuan  _et al._ (2018) Luqi Yuan, Qian Lin, Meng Xiao, and Shanhui Fan, “Synthetic dimension in photonics,” Optica 5, 1396–1405 (2018).

See pages 1 of TopofibreAppend See pages 2 of TopofibreAppend See pages 3 of
TopofibreAppend See pages 4 of TopofibreAppend See pages 5 of TopofibreAppend
See pages 6 of TopofibreAppend See pages 7 of TopofibreAppend See pages 8 of
TopofibreAppend


# Grand Motzkin paths and $\\{0,1,2\\}$-trees – a simple bijection

Helmut Prodinger  Department of Mathematics, University of Stellenbosch 7602,
Stellenbosch, South Africa and NITheCS (National Institute for Theoretical and
Computational Sciences), South Africa<EMAIL_ADDRESS>
###### Abstract.

A well-known bijection between Motzkin paths and ordered trees with outdegree
always $\leq 2$, is lifted to Grand Motzkin paths (the nonnegativity is
dropped) and an ordered list of an odd number of such $\\{0,1,2\\}$ trees.
This offers an alternative to a recent paper by Rocha and Pereira Spreafico.

###### Key words and phrases:

Motzkin paths, unary-binary trees, bijection

###### 2010 Mathematics Subject Classification:

05A15

##  1\. Introduction

Motzkin paths appear first in [6]. In the encyclopedia [9] they are enumerated
by sequence A001006, with many references given. They consist of up-steps
$U=(1,1)$, down-steps $D=(1,-1)$ and horizontal (flat) steps $F=(1,0)$. They
start at the origin and must never go below the $x$-axis. Usually one requires
the path to end on the $x$-axis as well, but occasionally one uses the term
_Motzkin path_ also for paths that end on a different level. Figure 1 shows
all Motzkin paths of 4 steps (=length 4).

Figure 1. All 9 Motzkin of 4 steps (length 4).

The enumeration of Motzkin paths is done using the generating function
$M=M(z)$ and a decomposition according to the first return to the $x$-axis,
viz.

$M=1+zM+z^{2}M^{2},$

this can be found in many books, e.g. in [5]. Solving,

$M(z)=\frac{1-z-\sqrt{1-2z-3z^{2}}}{2z^{2}}.$

The other combinatorial structure that plays a role in this note are _ordered
trees_. They are enumerated by an equation for the generating function
(according to the number of nodes)

$P=z+zP+zP^{2}+zP^{3}+\cdots=\frac{z}{1-P}\quad\text{and therefore}\quad
P=P(z)=\frac{1-\sqrt{1-4z}}{2}.$

The subclass of $\\{0,1,2\\}$-trees of interest only allows outdegrees $0,1,2$
and so we get again a generating function

$Q=z+zQ+zQ^{2}\quad\text{and therefore}\quad
Q=Q(z)=\frac{1-z-\sqrt{1-2z-3z^{2}}}{2z}=zM(z).$

So there should be a bijection of $\\{0,1,2\\}$-trees with $n+1$ nodes (= $n$
edges) and Motzkin paths of length $n$; the simplest I know is from [3]:

One runs through the $\\{0,1,2\\}$-tree in _pre-order_ ; if one sees an edge
for the _first_ time, one translates a single edge (degree 1) into a flat
step, a left edge into an up-step and a right edge into a down-step. It is
easy to see that the process is reversible, which is the desired bijection.

##  2\. Grand Motzkin paths

As a first step, we need the generating function of Motzkin paths, ending on
level $k$, not just the usual case $0$. This is a standard argument, by
decomposing such a path according to the last time you visit level $0$, then
an up-step, and we wait until we visit level $1$ for the last time, and so
one.

$\cdot$$\cdot$$\cdot$$\cdot$$\cdot$$\cdot$$\cdot$$\cdot$$\cdot$Motzkin
pathMotzkin pathMotzkin path Figure 2. The decomposition of a Motzkin path
that ends at level $2$.

From this we find the generating function as $z^{k}M^{k+1}$. This is well-
known.

Now we come to Grand Motzkin paths, a notation I picked up from [4]. This more
general family has the same steps as Motzkin paths, returns to the $x$-axis at
the end, but the condition that the paths must stay (weakly) above the
$x$-axis is dropped.

Let $k$ be the unique negative number such that the Grand Motzkin paths
reaches the level $-k$ when going from left to right; for $k=0$ they are just
ordinary Motzkin paths. We consider the first such point $(a,-k)$ and the last
such point $(b,-k)$; it is possible that $a=b$. Then the paths decomposes
canonically into 3 parts:

$\bullet$$\bullet$$\bullet$$\bullet$ Figure 3. The decomposition of a Grand
Motzkin path that has $-3$ as its minimal level.

This decomposition produces the generating function of such paths as a product
of 3 terms: The first part corresponds to $z\cdot z^{k-1}M^{k}$ (the last step
must be a down-step), the second part to $M$, and the third part again to
$z\cdot z^{k-1}M^{k}$ (the first step must be an up-step). The product is then
$z^{2k}M^{2k+1}$. Summing over all possible values of $k$, the enumeration of
Grand Motzkin paths is done via the generating function

$M+\sum_{k\geq 1}{z^{2k}M^{2k+1}}=\sum_{k\geq
0}{z^{2k}M^{2k+1}}=\frac{1}{\sqrt{1-2z-3z^{2}}}.$

In the world of $\\{0,1,2\\}$-trees, we form a new super-root, with $2k+1$
successors, each of which is a $\\{0,1,2\\}$-tree. The generating function is
then

$z\sum_{k\geq 0}Q^{2k+1}=z^{2}\sum_{k\geq
0}z^{2k}M^{2k+1}=\frac{z^{2}}{\sqrt{1-2z-3z^{2}}}.$

The previous bijection takes over to the new situation: if a grand Motzkin
path consists of $2k+1$ ordinary Motzkin paths, then each of them corresponds
to a $\\{0,1,2\\}$-tree as before. The extra factor $z^{2}$ stems from the
fact that for $\\{0,1,2\\}$-trees, the edges correspond to the steps. And now,
there is the super-root, which should not be counted when considering the
corresponding Grand Motzkin path.

The paper [8] has a correspondence between Grand Motzkin paths and
$\\{0,1,2\\}$-trees with a super-root with an odd number of successors.
However the arguments are perhaps less direct than the present ones.

##  3\. Trinomial coefficients and enumeration

The trinomial coefficients (notation from Comtet [1]) are given by

$\binom{n,3}{k}=[z^{k}](1+z+z^{2})^{n}.$

These coefficients are intimately related to the Motzkin-world, as we will
discuss for the reader’s benefit. We use the substitution
$z=\dfrac{v}{1+v+v^{2}}$ as we first did in [7]. Using this substitution, all
generating functions become much easier, like

$Q(z)=v,\quad\frac{1}{\sqrt{1-2z-3z^{2}}}=\frac{1+v+v^{2}}{1-v^{2}}.$

Coefficients can be extracted via contour integration, which is a variant of
the Lagrange inversion formula. See [2] and [7] as illustrations of the
technique. Note that when $z$ runs around the origin in a small circle, $v$
runs around the origin once as well, in a deformed circle. We will show two
sample computations:

$\displaystyle[z^{n}]\frac{1}{\sqrt{1-2z-3z^{2}}}$
$\displaystyle=\frac{1}{2\pi
i}\oint\frac{dz}{z^{n+1}}\frac{1}{\sqrt{1-2z-3z^{2}}}$
$\displaystyle=\frac{1}{2\pi
i}\oint\frac{dv(1-v^{2})}{(1+v+v^{2})^{2}}\frac{(1+v+v^{2})^{n+1}}{v^{n+1}}\frac{1+v+v^{2}}{1-v^{2}}$
$\displaystyle=\frac{1}{2\pi
i}\oint\frac{dv}{v^{n+1}}(1+v+v^{2})^{n}=[v^{n}](1+v+v^{2})^{n}=\binom{n,3}{n};$

and

$\displaystyle[z^{n}]Q^{j}$ $\displaystyle=\frac{1}{2\pi
i}\oint\frac{dz}{z^{n+1}}v^{j}=\frac{1}{2\pi
i}\oint\frac{dv(1-v^{2})}{(1+v+v^{2})^{2}}\frac{(1+v+v^{2})^{n+1}}{v^{n+1}}v^{j}$
$\displaystyle=\frac{1}{2\pi
i}\oint\frac{dv(1-v^{2})}{v^{n+1-j}}(1+v+v^{2})^{n-1}$
$\displaystyle=[v^{n-j}](1+v+v^{2})^{n-1}-[v^{n-j-2}](1+v+v^{2})^{n-1}=\binom{n-1,3}{n-j}-\binom{n-1,3}{n-j-2}.$

## References

  * [1] Louis Comtet.  Advanced Combinatorics.  D. Reidel Publishing Co., Dordrecht, enlarged edition, 1974.  The art of finite and infinite expansions.
  * [2] Nicolaas G. de Bruijn, Donald E. Knuth, and Stephen O. Rice.  The average height of planted plane trees.  In Graph Theory and Computing, pages 15–22. Academic Press, New York, 1972.
  * [3] Emeric Deutsch and Louis W. Shapiro.  A bijection between ordered trees and 2-Motzkin paths and its many consequences.  volume 256, pages 655–670. 2002.  LaCIM 2000 Conference on Combinatorics, Computer Science and Applications (Montreal, QC).
  * [4] Luca Ferrari and Emanuele Munarini.  Enumeration of edges in some lattices of paths.  Journal of Integer Sequences, 17:14.1.5 (22 pages), 2014.
  * [5] Philippe Flajolet and Robert Sedgewick.  Analytic Combinatorics.  Cambridge University Press, Cambridge, 2009.
  * [6] Theodore S. Motzkin.  Relations between hypersurface cross ratios, and a combinatorial formula for partitions of a polygon, for permanent preponderance, and for non-associative products.  Bulletin of the American Mathematical Society, 54 (4):352–360, 1948\.
  * [7] Helmut Prodinger.  The average height of a stack where three operations are allowed and some related problems.  J. Combin. Inform. System Sci., 5(4):287–304, 1980.
  * [8] Leandro Rocha and Elen V. Pereira Spreafico.  A combinatorial bijection between ordered trees and lattice paths.  Trends in Computational and Applied Mathematics, 24(3):427–436, 2023.
  * [9] Neil J. A. Sloane and The OEIS Foundation Inc.  The On-Line Encyclopedia of Integer Sequences, 2023.


|  Spontaneous flows and dynamics of full-integer topological defects in polar
active matter
---|---
|  Jonas Rønning,a Julian Renaud,b Amin Doostmohammadi,∗c and Luiza
Angheluta,∗∗a
|  Polar active matter of self-propelled particles sustain spontaneous flows
through the full-integer topological defects. We study theoretically the
effect of both polar and dipolar active forces on the flow profile around $\pm
1$ defects and their interaction in the presence of both viscosity and
frictional dissipation. The vorticity induced by the active stress is non-zero
at the $+1$ defect contributing to the active torque acting on the defect. A
near-core flow reversal is predicted in absence of hydrodynamic screening
(zero friction) as observed in numerical simulations. While $\pm 1$ defects
are sources of spontaneous flows due to active stresses, they become sinks of
flows induced by the polar active forces. We show analytically that the flow
velocity induced by polar active forces increases away from a $\pm 1$ defect
towards the uniform far-field, while its associated vorticity field decays as
$1/r$ in the far-field. In the friction-dominated regime, we demonstrate that
the flow induced by polar active forces enhances defect pair annihilation, and
depends only on the orientation between a pair of oppositely charged defects
relative to the orientation of the background polarization field.
Interestingly, we find that this annihilation dynamics through mutual defect-
defect interactions is distance independent, in contradiction with the effect
of dipolar active forces which decay inversely proportional to the defect
separation distance. As such, our analyses reveals a new, truly long-ranged
mechanism for the pairwise interaction of oppositely-charged topological
defects in polar active matter.

##

††footnotetext: a Njord Centre, Department of Physics, University of Oslo, P.
O. Box 1048, 0316 Oslo, Norway††footnotetext: b École Normale Supérieure, PSL
Research University, 45 rue d’Ulm, 75005 Paris, France††footnotetext: c Niels
Bohr Institute, University of Copenhagen, Blegdamsvej 17, Copenhagen,
Denmark††footnotetext: ${\dagger}$ Electronic Supplementary Information (ESI)
available: [https://github.com/jonasron/SM_Flowfield_Polar]. See DOI:
10.1039/cXsm00000x/††footnotetext: $\ast$ Email:
<EMAIL_ADDRESS>$\ast\ast$ Email<EMAIL_ADDRESS>
##  1 Introduction

Active matter refers to non-equilibrium systems of interacting, self-propelled
entities that consume energy from their surrounding in the form of persistent
motion and their collective interactions lead to emergent, dynamical patterns,
and self-sustained flows 1, 2. Models of active matter are largely inspired by
biological systems from bacterial suspensions 1, 2 and cell monolayers 3, 4
and down to subcellular active systems such as mixtures of cytoskeletal
filaments and motor proteins 5, 6, 7, 8. However, this also pertains to non-
living active systems such as layers of vibrated granular matter, microrobots,
or synthetic catalytic nanomotors 9, 2.

Several hydrodynamic models have been proposed to capture the macroscopic
dynamics and emergent phenomena corresponding to a collection of active (self-
propelled) particles with different symmetries and alignment interactions 10,
2, 11, 12. The prototypical models are based on the analogy to liquid crystals
formed by rod-like particles with polar or nematic symmetries in their
alignment interactions. Active rods with only orientational alignment act as
headless “shakers" and form nematic phases described by a slowly-varying
director field $\vec{n}$ which has head-tail (nematic) symmetry 13, 2, 1, 14.
By contrast, active rods that align their direction of motion tend to flock
into polar systems that are described instead by a slowly-varying polar vector
field $\vec{p}$ 11, 15, 10, 12. Active rod-like particles generate persistent
flows sustained by the active stress originating from extensile/contractile
dipolar forces. On hydrodynamic scales, this active stress is proportional to
the nematic $\mathbf{Q}$ tensor order parameter, i.e.
$\sigma^{a}=\alpha_{0}\mathbf{Q}$ with a proportionality coefficient
$\alpha_{0}$ as an effective activity parameter 2, 16. Topological defects,
innate to ordered systems with broken continuous symmetries, are also present
in active systems. The interplay between active stresses and nematic
distortions feeds into self-sustain flows and the proliferation of topological
defects to generate chaotic flows also known as active turbulence 17, 18,
chandragiri2020flow.

From the rotational symmetry modulo $\pi$ of the $\mathbf{Q}$ tensor, it
follows that the lowest energy orientational defects have a $\pm 1/2$
topological charge corresponding to $\pm\pi$ jump in the orientational phase
around them. The $+1/2$ defects acquire a self-propulsion due to the net
active flow passing through their cores 19, 20. The stable, low-energy defects
in polar active systems have instead $\pm 1$ topological charges corresponding
to a $2\pi$ phase jumps around the defects. This is analogous to vortices in
the XY-model of 2D ferromagnets 15, 21.

The bulk of recent studies have focused on the formation and characterization
of half-integer nematic defects (see recent reviews 22, 23). This is in part
due to the ubiquitous emergence of the nematic defects in a wide range of
biological systems from subcellular filaments 24, 8 to bacterial colonies 25,
26 and assemblies of eukaryotic cells 27, 28. This is despite the fact that
several biological active entities, such as bacteria or eukaryotic cells, are
endowed with a clear head-tail asymmetry and directional self-propulsion,
which characterizes a polar order for these systems. Such a polar order at the
scale of collective is apparent from flocking domains within bacterial
colonies 26 and eukaryotic cells 29. Nevertheless, because of the appearance
of half-integer defects, at the coarse-grained level, these systems are often
modelled as active nematics neglecting the polarity of the self-propelled
particles. There have been proposed models that couple the evolution of polar
and nematic order parameters 30, 31, 32 to allow for the coexistence of both
types of symmetries. Recently, in Ref. 33, a hydrodynamic model was proposed
for the coexistence of both nematic and polar alignment interactions in the
polarization field $\vec{p}$, and was used to study different active
turbulence regimes sustained by both half- and full- integer topological
defects. In addition to dipolar (nematic) active forces, the model also
includes a polar active force $\vec{f}=\alpha_{p}\vec{p}$, which describes the
self-propulsion direction and is shown to suppress defect-laden active
turbulence 34.

More recently, theoretical and experimental studies have revealed the
importance of full-integer topological defects in cell assemblies. In
particular, it is shown that positive full-integer defects formed due to
collision of two nematic half-integer defects in fast-moving bacterial
colonies can lead to the verticalization of bacteria and escape to the out-of-
plane directions 26, 35. Furthermore, positive full-integer defects induced by
confinement of myoblast cells in circular geometries are shown to activate
cell differentiation and formation of 3D helical structures 36. A
corresponding theoretical analyses, in the limit of compressible flows inside
the core region of defects in small confinements, have revealed the
corresponding flow fields and force patterns, and shows how confinement-
induced topological defects can be used to probe the material properties of
the cell layers 37.

In this paper, we carry out a theoretical analysis that reveals subtle cross-
talks between polar and dipolar active forces in generating spontaneous
incompressible flows around both $\pm 1$ defects. We theoretically predict
spontaneous vortical flows induced by dipolar active forces around $+1$
defects which correspond to isotropic active stress and pressure fields.
Interestingly, the competition between dipolar force and viscous force leads
to a flow reversal close to the core of the $+1$ defect. This effect is also
confirmed by numerical simulations of the full hydrodynamic model including
additional passive stresses. However, it turns out the hydrodynamic screening
induced by friction lifts up this flow reversal effect. By contrast, the $-1$
defects which have an innate $4$-fold symmetry in the polarity field lead to
$8$-fold symmetries of the active flows induced by dipolar active forces. The
same $8$-fold symmetry is present also in the profile of the pressure field.
We demonstrate that polar active force trigger a distinct active flow pattern
characterised by uniform flow in the far-field of its source, i.e. $\pm 1$
defects. This is very different that the active flows sustained by dipolar
forces which decay algebraically with distance. Polar active forces have a
drastic effect on the mutual interactions between defects by promoting fast
annihilation of defect pairs through non-local and non-reciprocal attraction
forces. We show that this distance-independent mechanism of defect pair
annihilation is responsible for the suppression of active turbulence by polar
forces as recently reported in Ref. 34.

The paper is organized as follows: We begin in Section 2 by introducing the
flow equations within a minimal polar hydrodynamic model and derive the
corresponding defect kinematic equations using Halperin-Mazenko formalism 38.
The main analytical results on the active flow velocity induced around
isolated $\pm 1$ are discussed in Section 3. We also compare the analytical
predictions with direct numerical simulations and find very good agreement of
the flow profiles around defects and the flow reversal pattern near the $+1$
defect cores. In Section 4, we discuss the polar active force and its effect
on the defect kinematics. In particular, we consider the motion of a pair of
oppositely charged full-integer defects under the polar active flows induced
by each defect and demonstrate that polar active forces enhance the defect
annihilation rate. In a recent study Ref. 34, it was numerically evidenced
that polar active forces suppress defect-laden active turbulence and tend to
restore polar order. Here, we demonstrate theoretically that polar active
forces have a net effect on the defect kinematics to promote defect pair
binding and subsequent annihilation of defects of opposite topological
charges. The final section provides a summary of the theoretical insights and
concluding remarks.

##  2 Hydrodynamic model of polar active matter

We consider the hydrodynamic model 33 that describes the collective dynamics
of self-propelled entities in terms of the evolution of the polar order
parameter $\vec{p}$ given by

$\displaystyle\partial_{t}\vec{p}+\vec{u}\cdot\nabla\vec{p}+\lambda\mathbf{E}\cdot\vec{p}+\mathbf{\Omega}\cdot\vec{p}$
$\displaystyle=-\frac{1}{\gamma}\frac{\delta\mathcal{F}}{\delta\vec{p}},$ (1)

where $\gamma$ is the rotational viscosity, $\lambda$ is the flow-aligning
parameter 39, and the free energy favoring the polar order is described as

$\displaystyle\mathcal{F}=\int
d\vec{r}\left\\{A\left(-\frac{|\vec{p}|^{2}}{2}+\frac{|\vec{p}|^{4}}{4}\right)+\frac{K_{p}}{2}|\nabla\vec{p}|^{2}\right\\}.$
(2)

Here, $K_{p}$ is the isotropic elastic constant for distortions in the
polarity field and $A$ is the height of local energy barrier. The polarity is
coupled with the flow field $\mathbf{u}$ which is described by the
incompressible Stokes equations

$\displaystyle(\Gamma-\eta\nabla^{2})\vec{u}$
$\displaystyle=\alpha_{0}\nabla\cdot\mathbf{Q}-\nabla P,$ (3)
$\displaystyle\nabla\cdot\vec{u}$ $\displaystyle=0,$ (4)

where the dipolar active forces are mediated by the nematic order through the
nematic tensor $Q_{ij}=(p_{i}p_{j}-\frac{p^{2}}{2}\delta_{ij})$ and the
activity parameter $\alpha_{0}$. The incompressibility constraint determines
the fluid pressure $P$. The viscosity is set by $\eta$ while the friction with
a substrate is introduced by the frictional drag $\Gamma$ 40. Since we focus
on theoretical derivations of the active flows induced by polar and dipolar
forces, we hereby neglect the additional passive stresses which depend on the
polarity and that are typically present in the Stokes equations simulated
numerically.

We consider the dimensionless forms of these equations following the
appropriate rescalings in units of length $\xi=\sqrt{K_{p}/A}$ and time
$\tau=\gamma/A$. The system is then controlled by three dimensionless
parameters: the scale number $\zeta=\ell_{d}/\xi$ is the ratio between the
hydrodynamic dissipation length $\ell_{d}=\sqrt{\eta/\Gamma}$ and the nematic
coherence length $\xi$, and the rescaled activity
$\tilde{\alpha}_{0}=\alpha_{0}\tau/(\xi^{2})=\alpha_{0}\gamma/(K_{p})$. Thus,
the dimensionless flow equations read as

$\displaystyle(1-\zeta^{2}\nabla^{2})\vec{u}=\vec{F}_{a}-\nabla P,$ (5)
$\displaystyle\nabla^{2}P=\nabla\cdot\vec{F}_{a},$ (6)

where the dipolar active force induced by the active stress is

$\vec{F}_{a}=\frac{\tilde{\alpha}_{0}}{\Gamma}\nabla\cdot\mathbf{Q}.$ (7)

For infinite systems, the flow velocity and pressure originated from dipolar
active forces can be calculated from convolution integrals of the source terms
and the corresponding Green’s functions as

$\displaystyle P(\vec{r})$ $\displaystyle=\frac{1}{2\pi}\int
d\vec{r^{\prime}}\ln(|\vec{r}-\vec{r^{\prime}}|)\left[\nabla^{\prime}\cdot\vec{F}_{a}\right],$
(8)

and

$\displaystyle\vec{u}(\vec{r})$
$\displaystyle=\int\frac{d\vec{r^{\prime}}}{2\pi\zeta^{2}}K_{0}\left(\frac{|\vec{r}-\vec{r^{\prime}}|}{\zeta}\right)\left(\vec{F}_{a}(\vec{r^{\prime}})-\nabla^{\prime}P(\vec{r^{\prime}})\right).$
(9)

To derive analytical expressions, we evaluate the source term due to the
active stress using the order parameter for an isolated point defect located
at the origin and given by

$\vec{p}(r,\theta)=\chi(r)[\cos(q\theta+\phi)\vec{e}_{x}+\sin(q\theta+\phi)\vec{e}_{y}],$

where $\theta=\arctan(y/x)$ is the polar angle in a coordinate system centered
at the defect position and $\chi(r)$ is the core function, where we assume
that the core size is much smaller than any other length scales and set
$\chi(r)=1$. The angle of the polar vector order parameter for an ideal defect
is given by $\Theta(\theta)=q\theta+\phi$ 21, where $q\theta$ is the singular
part giving the winding number $q$ when performing the integral $\oint
d\Theta=2\pi q$ on a contour surrounding the defect and $\phi$ is a constant.
For defects of charge $q\neq+1$ this constant sets the orientation of the
defect and can be ignored since we can always transform to a system with
$\phi=0$ by a change of basis. When $q=+1$ it is impossible to remove a non-
zero $\phi$ by changing the basis 41. The positive defects are rotationally
symmetric and the constant $\phi$ becomes important since it distinguishes
different types of defects: $\phi=0$ corresponds to an aster and $\phi=\pi/2$
gives rise to a vortex, and any value in-between corresponds to a spiral
defect. For a more intuitive depiction of this, we plot in Fig. 1 a) and c)
the $\vec{p}$ field on the circumference of a circle centered at $\pm 1$
defects to show that $\phi$ corresponds to the constant angle that $\vec{p}$
makes with the radial direction. For the negative defect one can define a
polarisation from the line where $\vec{p}$ is pointing radially outward from
the core, the angle of this polarisation is half of $\phi$. The different
patterns of the $\vec{p}$ field around a $\pm 1$ defect are also illustrated
in Fig. 1 b) and d).

Fig. 1: Polarity field around the $\pm 1$ topological defects. a) and b) show
a $+1$ defect, while c) and d) illustrate a $-1$ defect. The drawings a) and
c) show the $\vec{p}$ plotted on a circle with the angle $\phi$ marked. This
illustrates why $\phi$ can be set to zero by a change of basis for the $-1$
defect and not for the $+1$ defect. b) and d) show the polarity fields around
a $+1$ defect with $\phi=\pi/3$ (spiral) and a $-1$ defect with $\phi=0$
respectively (Changing $\phi$ rotates the negative defect).

We use the Halperin-Mazenko formalism 42, 43 for tracking topological defects
as zeros of the polar order parameter to derive the corresponding equations of
motion of defects from the evolution of the $\vec{p}$ field similar the
approach from Ref. 38 to describe orientational defects in active nematic
films. The basic idea is that since the $\pm 1$ defects are associated
simultaneously with topological singularities in the orientation field and
zero magnitude of the vector order parameter, we can track their position by
Dirac delta functions centered at the zeros of the $\vec{p}$ field. Hence, a
configuration of well-separated defects punctuating the $\vec{p}$-field
corresponds to a defect charge density field, which can be written
equivalently either as a superposition of the delta functions associated with
the topological singularities located in the physical space or as the zeros in
the order parameter space,

$\rho(\vec{r},t)=\sum_{n}q^{(n)}\delta(\vec{r}-\vec{r}^{(n)})=D\delta(\vec{p}),$
(10)

where $D$ is the determinant of the polarity distortion tensor
$\nabla\vec{p}$, i.e. $D=\epsilon_{ij}\partial_{i}p_{x}\partial_{j}p_{y}$
which can be expressed equivalently in the complex representation
$\psi=p_{x}+ip_{y}$ as
$2iD=\epsilon_{ij}\partial_{i}\bar{\psi}\partial_{j}\psi$. The $D$-field is
zero in regions of uniform polar order and becomes nonzero where there are
distortions in the orientation field. For configurations of well-separated
defects punctuating a uniform polar order state, the $D$ field is zero
everywhere except at the defect positions labeled by the index $n$, where the
topological charge $q_{n}=\pm 1$ is determined by the sign of $D$. Thus, the
$D$ field represents a non-singular charge density field.

It can be shown that the $D$-field follows the conservation law 42

$\partial_{t}D+\partial_{i}\vec{J}_{i}=0,$ (11)

with the current density
$J=\epsilon_{ij}\Im(\partial_{t}\psi\partial_{j}\bar{\psi})$ determined by the
evolution of the polar order. Combining this with the conservation of
topological charge density $\rho$, we can derive a general expression for the
defect velocity in terms of the $D$ and its current,
$\vec{v}_{n}=\vec{J}/D|_{\vec{r}_{n}}$ 38. We parameterise the polarity field
from the viewpoint of the $n$-th defect as $\psi=\psi_{n}e^{i\phi_{n}}$, where
$\phi_{n}$ is the smooth phase at the defect position $\mathbf{r}_{n}$ and
$\psi_{n}=|\vec{r}-\vec{r}_{n}|e^{i\theta_{n}}$ using that the polarization
vanishes linearly while its phase $\theta_{n}$ is singular at the pointwise
defect position 44, 45. Within this approach, the general expression of the
defect current density can be reduced to a closed expression for the defect
velocity given by

$\displaystyle\vec{v}_{n}=\vec{u}(\vec{r}_{n})+2q_{n}\nabla^{\perp}\phi_{n}\bigg{|}_{\vec{r}=\vec{r}_{n}}.$
(12)

Using the stationary phase approximation, i.e. that the constant (equilibrium)
phase remains stationary when it is punctuated by moving defect singularities,
we can further simplify the defect velocity and express it in terms of the net
spontaneous flow and the forces induced by the other defects as 46, 38

$\vec{v}_{n}=\vec{u}(\vec{r}_{n})+2\sum_{k\neq
n}q_{n}q_{k}\frac{\vec{r}_{n}-\vec{r}_{k}}{|\vec{r}_{n}-\vec{r}_{k}|^{2}}.$
(13)

The topological defects interact through Coulomb-like forces, where like-
signed defects repel, and opposite-signed defects attract each other. However,
there are additional interactions through the flow field $\vec{u}$ which
depends on the dipolar active forces. In the next section, we derive analytic
expressions for this flow velocity and discuss its effect on the defect
motion.

##  3 Active flow fields around $\pm 1$ topological defects

From the parameterization of the $\vec{p}$-field for a pointwise defect, we
can evaluate the dipolar active force induced by an isolated defect with
charge $q=\pm 1$ in an otherwise homogeneous polarity field with constant
background orientation $\phi$ as

$\displaystyle\vec{F}_{a}^{+}=\frac{\tilde{\alpha}_{0}}{\Gamma
r^{2}}\left[\vec{r}\cos(2\phi)-\vec{r}^{\perp}\sin(2\phi)\right],$ (14)
$\displaystyle\vec{F}_{a}^{-}=-\frac{\tilde{\alpha}_{0}}{\Gamma
r}\mathbf{R}_{2\phi}\left[\cos(3\theta)\vec{e}_{x}-\sin(3\theta)\vec{e}_{y}\right],$
(15)

where $\vec{r}^{\perp}=(y,-x)$, and $\mathbf{R}_{2\phi}$ is a matrix that
rotates the vector by $2\phi$, and can be removed by a change of basis. Notice
that the first term in $\vec{F}_{a}^{+}$ is a source of gradient flow, which
however is removed by pressure through the incompressibility condition. The
second term proportional to $\sin(2\phi)$ is related to a rotated gradient and
induces a purely vortical flow.

We insert these forces into the integrals in Eq. (8) and (9), and solve them
as described in the Supplementary Material 47. The resulting expressions for
the active flow velocity and pressure for a $+1$ defect reduce to

$\displaystyle u_{a}^{+}(r,\theta)$
$\displaystyle=\frac{\tilde{\alpha}_{0}}{\zeta\Gamma\hat{r}}(1-\hat{r}K_{1}(\hat{r}))e^{i(\theta+\pi/2)}\sin(2\phi),$
(16) $\displaystyle P_{a}^{+}$
$\displaystyle=-\frac{\tilde{\alpha}_{0}}{\Gamma}\cos(2\phi)\ln{\frac{L}{r}}.$
(17)

Here $L$ is a cutoff scale set by the system size. We have written the
velocity field as a complex field defined as $u=u_{x}+iu_{y}$, where $r$ and
$\theta$ are the polar coordinates centered at the defects position. We define
the scaled radial coordinate as $\hat{r}=r/\zeta$, which is the same as
changing the length scale to the hydrodynamic dissipation length $\ell_{d}$
(see Fig. 2, with the cutoff scale $L=50$). Notice that the pressure is
isotropic and its gradient force cancels the radial component of the dipolar
active force, thus no net pressure flow. Furthermore, the corresponding
vorticity is also isotropic

$\omega^{+}_{a}=\frac{\tilde{\alpha}_{0}}{\Gamma\zeta^{2}}K_{0}(\hat{r})\sin(2\phi),$
(18)

and is non-zero at the centre of the defect as evidenced also in Fig. 2. This
has important consequences for the stability of vortex, spiral, or aster
shaped positive defects as discussed in detail in Ref. 15. The sign of the
global circulation is modulated by the character of the $+1$ defect through
$\phi$.

Fig. 2: Incompressible active flow streamlines (a) and corresponding pressure
field (b) around an isolated $+1$ vortex with $\sin(2\phi)>0$ in an extensile
($\tilde{\alpha}_{0}<0$) system. The colormap shows the magnitude of the
velocity and pressure normalized by the maximum. Note that the pressure in the
far-field diverges with the system size $L$. Since the pressure diverges at
the centre the pressure inside the core $r<1$ has been set to be constant. c)
shows the normalised pressure, for an extensile system, at a fixed radius as a
function of $\phi$. Notice that for $\phi=\pm\pi/4+\pi n$ the pressure
vanishes.

Similarly, we find analytical expressions for the active flow velocity,
pressure and vorticity related to the $-1$ defect, which are written in a
compact form as

$\displaystyle u^{-}_{a}(r,\theta)$
$\displaystyle=-\frac{\tilde{\alpha}_{0}}{2\Gamma\zeta}\left(f_{3}^{a}(\hat{r})e^{-3i\theta}+f_{5}^{a}(\hat{r})e^{5i\theta}\right),$
(19) $\displaystyle P_{a}^{-}$
$\displaystyle=-\frac{\tilde{\alpha}_{0}}{4\Gamma
r^{4}}(x^{4}-6x^{2}y^{2}+y^{4}),$ (20) $\displaystyle\omega_{a}^{-}(r,\theta)$
$\displaystyle=-\frac{\tilde{\alpha}_{0}}{\Gamma\zeta^{2}}\sin{4\theta}f_{\omega}^{a}(\hat{r}).$
(21)

The functions $f_{3}^{a}$, $f_{5}^{a}$, and $f_{\omega}^{a}$ giving the radial
dependence of the velocity and vorticity are listed in the Appendix.

Fig. 3: Incompressible active flow streamlines (a) and corresponding pressure
field (b) around an isolated $-1$ vortex with $\phi=0$ in an extensile
($\tilde{\alpha}_{0}<0$) system. The colormaps show the normalised magnitude
of velocity (a) and pressure (b).

These fields are also plotted in Fig. 3 showing that while the positive defect
has closed streamlines also in an infinite system, the negative defect has an
8-fold symmetry of the vorticity, underlying a 4-fold symmetry of the polarity
field around the defect. The streamlines do not close in an infinite system,
but in confinement they might close due to boundary conditions or influence
from other defects similarly to the vortices formed around $\pm 1/2$ defects
in an active nematic 48, 20. It is important to note the distinctions of the
analytical descriptions provided herein with a recent calculations of
compressible flow fields inside circular confinements 37. In the latter, the
analyses is restricted to circular domains with radius smaller than the
coherence length $R^{2}\ll K_{p}/A$, i.e the limit where
relaxation/penetration length in the free energy is larger than the system
size. Those results are therefore in the opposite limit of the calculation
performed here.

In addition to the characteristic flow fields, our closed form analytic
descriptions provide an insight into the isotropic stress patterns, i.e., half
of the trace of the stress tensor, around full-integer topological defects.
Such isotropic stresses have been shown to be a determining factor for the
biological functionality of nematic defects, where concentration of
compressive stress around $+1/2$ defects was shown to lead to cell death and
extrusion 27, while the tensile stresses around $-1/2$ defects has been shown
to lead to spontaneous gap opening in epithelial cell layers 49. Similarly, we
find distinct isotropic stress patterns around positive and negative full-
integer defects: while around a negative ($-1$) defect alternating regions of
tension and compression appear in a 4-fold symmetric pattern (Fig. 3b), a
strong augmentation of compressive stress is observed at the core of positive
$+1$ defects (Fig. 2b). We conjecture that such a concentrated compressive
stress at $+1$ defects could contribute to the activation of mechanosensitive
signals in cell layer and potentially be linked to the recent observation of
the cellular differentiation at +1 defect cores in cartilage cells 50.

###  3.1 Vortical flow reversal around $+1$ defect

The results presented above provide closed-form analytical expressions for
incompressible flows around full-integer topological defects in polar active
matter. Interestingly, in the viscous-dominated regime ($\Gamma=0$), we find
that the activity-induced flow around the $+1$ defect forms two counter-
rotating vortices around the defect core, and thus the rotational flow
reverses sign on a length scale comparable with the coherence length outside
of the defect core. In this limit, the flow velocity becomes instead

$\displaystyle
u_{a}^{+}=\frac{i\tilde{\alpha}_{0}}{4\tilde{\eta}}r(1-2\ln(r))e^{i\theta}\sin(2\phi),$
(22)

and its corresponding vorticity is thus

$\omega^{+}_{a}=-\frac{\tilde{\alpha}_{0}}{\tilde{\eta}}\sin(2\phi)\ln(r),$
(23)

where $\tilde{\eta}=\eta/\xi^{2}$. We notice that the vorticity changes
handedness for $r=\sqrt{e}\approx 1.7$ giving rise to two counter-rotating
regions as shown in Fig. 4. Notice that the flow diverges at
$r\rightarrow\infty$ since there is no screening lengthscale. The mathematical
reason for this flow reversal is that the logarithmic Greens function changes
sign when $r<1$, unlike the Bessel function $K_{0}$ kernel which appears when
both friction and viscosity are present. Note that in the absence of friction,
the flow equations have no natural lengthscale, thus the rescaling must be
done with an appropriate lengthscale coming from the free energy, such as the
coherence length which is consistent with the numerical observations of the
flow reversal discussed next. Similar flow reversal is also expected in
systems confined to circular domains that are small enough such that the
energy scale can be treated as linearly dependent on the radial distance from
the center of singularity $\chi(r)\sim r$ in the entire domain 37. For such
systems the flow reversal happens inside the core region.

Fig. 4: (a) Streamlines of the theoretical active flow around an isolated $+1$
vortex for $\sin(2\phi)>0$ for an extensive system and $\Gamma=0$. The
colormap shows the flow magnitude. The change in color of the streamlines from
white to black highlight the flow reversal at $r\approx 1.7$. This is more
apparent in (b) which is a cross section of the azimuthal velocity component
along the positive $x$ axis. The velocity changes sign when it crosses the
zero line (doted black).

###  3.2 Comparison with numerical simulations

(a) $\tilde{\alpha}_{0}=+0.4$

(b) $\tilde{\alpha}_{0}=+1$

(c) $\tilde{\alpha}_{o}=-0.4$

(d) $\tilde{\alpha}_{0}=-1$

Fig. 5: Numerical average polarity and velocity fields around $(+1)-$defects,
for various activities. $\vec{p}$ (resp. $\vec{u}$) field is plotted on the
right (resp. left) of each panel. Averages were computed on $\simeq 100$
defects for $\tilde{\alpha}_{0}=\pm 0.4$, and on $\simeq 1000$ defects for
$\tilde{\alpha}_{0}=\pm 1$. The polarity field consists each time of a
statistical combination of asters, spirals and vortices.

(a) $\tilde{\alpha}_{0}=+0.4$

(b) $\tilde{\alpha}_{0}=+1$

(c) $\tilde{\alpha}_{0}=-0.4$

(d) $\tilde{\alpha}_{0}=-1$

Fig. 6: Numerical average polarity and velocity fields around $(-1)-$defects,
for various activities. $\vec{p}$ (resp. $\vec{u}$) field is plotted on the
right (resp. left) of each panel. Averages were computed on $\simeq 100$
defects for $\tilde{\alpha}_{0}=\pm 0.4$, and on $\simeq 1000$ defects for
$\tilde{\alpha}_{0}=\pm 1$. Abnormal points (max. 8, on the edge of the
averaging window) were set down to 0 to retain only significant flow patterns.

Here we discuss the comparison of the analytical predictions of active flow
induced by $\pm 1$ defects with those obtained from direct numerical
simulations of the polarization evolution Eq. (1-2) coupled with the Navier-
Stokes equations

$\rho(\partial_{t}\vec{u}+\vec{u}\cdot\nabla\vec{u})=\alpha_{0}\nabla\cdot\mathbf{Q}-\nabla
P+\nabla\cdot\mathbf{\sigma}_{p},\quad\nabla\cdot\vec{u}=0$ (24)

with the additional passive stresses $\mathbf{\sigma}_{p}$ including the
viscous stress as defined in Ref. 18. In numerical simulations, we use a
hybrid lattice-Boltzmann method, combining finite-difference for the evolution
of polarity and the lattice-Boltzmann method for that of velocity. Using the
same prescription as in Ref. 34 for density and viscosity, i.e. $\rho=40$ and
$\eta=3.6$, we ensure that the Reynolds number in the simulations remains
negligible ($Re\ll 1$) 51, 18 so that the dynamics of velocity virtually
reduces to the incompressible Stokes equations Eq. (3-4) considered here. We
fix the viscosity ratio to $\eta/\gamma=3.6$, micro to macro length scale to
$(\sqrt{K_{p}/A})/L=2\times 10^{-3}$ (assuring that the coherence length
$\sqrt{K_{p}/A}$ is significantly smaller than the domain size $L$), and the
flow alignment parameter to $\lambda=0.1$. Dimensionless parameters in the
simulations are defined similar to the theoretical parameters defined in
Section 2.

Simulations were initialized with quiescent velocity field and noisy polar
alignments close to the uniformly oriented state $\vec{p}=\vec{e}_{x}$ under
periodic boundary conditions, on square domains of linear dimension $L=256$.

We ran simulations with sufficiently large values of the activity parameter
$\tilde{\alpha}_{0}$ for $\pm 1$ topological defects to form spontaneously 34.
We compile the flow profiles in the vicinity of each defect for a large sample
of defects (ranging from $\simeq 10^{2}-10^{3}$) for each simulation. From
these individual flow profiles, we compute the statistically-averaged flows
generated by $\pm 1$ topological defects. However, because every defect had
its own orientation, the flows had to be carefully reoriented upstream of the
averaging procedure, as we detail below.

The orientation of $-1$ defects was defined from the surrounding polarity
field. Given the 2-fold symmetry of the $-1$ defect (see fig. 1d)), we search
for the two principal axes along which the polarity vector points inwards and
outwards, respectively, and align the local fields accordingly.

Around $+1$ defects, both the polarity field and the velocity field are
rotationally-symmetric (see Fig. 1a,b and Fig. 4). However, the vortical
velocity around $+1$ defects may have different chirality: negative
("clockwise") when $u_{\theta}<0$, or positive ("counterclockwise") when
$u_{\theta}>0$ as in Fig. 4. When simply summed up, these two chirality would
cancel out, even though they have equivalent dynamic properties. To avoid
this, we identify the flow chirality of each defect before averaging the
fields. The flow chirality of each $+1$ defect was probed systematically by
computing the sign of the circulation of $\vec{u}$ along small circular
contours enclosing the defect center within the core of the defect, i.e. with
the contour radius smaller than the length scale of the flow reversal
$\sqrt{e}$. Defects with clockwise chirality ($\oint\ u_{\theta}<0$) are
reversed ($p_{\theta}\leftarrow-p_{\theta},u_{\theta}\leftarrow-u_{\theta}$),
while those with counterclockwise chirality are kept with the same
orientation, so that all fields are consistently averaged. Note that no
further distinction is made with regard to the polarity field, which therefore
aggregates asters, spirals and vortices, since they lead to the same flow
chirality.

The numerical results of the average polarity and flow profiles around $\pm 1$
defects are shown in Fig. 5 and Fig. 6. First, for both types of defects, the
profile of the polarity field (right panels in Fig. 5 and Fig. 6) is
consistent with the ideal vortex ansatz made by neglecting the backflow and
the finite size of the core. The symmetries – of the average – are notably the
same as in the passive case (see Fig. 1): positive defects are rotationally
invariant – seemingly with a mixture of asters, spirals and vortices –, while
negative defects feature a 2-fold symmetry. Nonetheless, since our model
allows the polarity vector to shrink and vanish, the local magnitude of
polarity constitutes an additional degree of freedom that yields qualitatively
robust patterns around $(-1)-$defects, yet with a shift by a $\pi/4-$rotation
between contractile ($\bar{\alpha}_{0}>0$) and extensile
($\bar{\alpha}_{0}<0$) systems, as can be seen from the colored fonts in fig.
6( 6(b), 6(b) and 6(d), 6(d) ). These patterns become smaller as the net
activity increases since the active length-scale is inversely proportional to
the activity $l^{\text{active}}\sim\sqrt{K/\alpha_{0}}$ (see fig. 6 6(b), 6(d)
and 6(b), 6(d) ).

We now discuss the flow fields generated by negative defects. As predicted
analytically above, the negative full-integer defects feature 8-fold
rotational symmetry, with regions where the velocity points alternatively
inwards and outwards, interspersed with counter-rotating swirls. This
structure is relatively robust, and strengthens as the level of activity
increases. The flow structure is also very similar for contractile and
extensile systems, up to a rotation by the same angle $\pi/4-$ as above. In
other words, equivalently, their profiles match, provided that one reverses
the direction of the polarity vector according to the sign of the activity,
which is consistent with the factor $\bar{\alpha}_{0}$ in the expression of
the field $\vec{u}$ we previously derived.

Next, we discuss the velocity field generated by the positive defects. As
expected from our theoretical considerations (fig. 2, fig. 4), the flow field
exhibits a full-rotational symmetry. Since there was no friction, the dynamics
was plainly viscous-dominated, and we have indeed striking evidence of the
flow reversal predicted by our calculations in this regime, for any level of
activity: as in the previous fig. 4, one can identify $2$ inner and outer
regions rotating contrariwise in fig. 5 5(b)-5(d), that we will refer to as
‘core’ and ‘shell’ respectively. (In our captions, the core is rotating
counterclockwise due to our averaging selection process). Taken together, the
averaged velocity fields found from numerical simulation of polar active
matter corroborate analytical predictions of the flow structures around full-
integer defects and existence of counter-rotating vortices around positively
charged defects. Furthermore, the simulation results were obtained from
averaging the flow fields around full-integer defects within the active
turbulence state and their close agreement with the analytical predictions
around isolated defects give credit to the assumption of neglecting defect-
defect interactions in the analytical calculations.

###  3.3 Defect pair interaction modulated by dipolar active forces

As seen in the above section the dipolar active force does not lead to any
defect self-propulsion for isolated $\pm 1$ defects in an infinite domain. We
now consider a pair of $\pm 1$ defects and show that the dipolar active force
leads to an interaction between them through the active flow. For analytic
tractability, we consider the friction-dominated regime $\Gamma\gg\eta$ and
set the viscosity to zero, $\ell_{d}\rightarrow 0$. The defect velocity is
calculated using the methods from 38 (see 47 for details) and is given by

$\displaystyle u^{+}_{a}=\frac{\tilde{\alpha}_{0}}{\Gamma
R}e^{i(3\varphi-2\phi)}\left(\ln\left(\frac{R}{\Lambda}\right)-\frac{3}{2}\right),$
(25) $\displaystyle u^{-}_{a}=\frac{\tilde{\alpha}_{0}}{\Gamma
R}\left(\frac{1}{6}e^{i(3\varphi-2\phi)}+\frac{1}{4}e^{i(2\phi-\varphi)}\right),$
(26)

for the $+1$ defect and $-1$ defect, respectively. Here $\Lambda$ is a
regularising lower cut-off corresponding to the finite defect core due to the
divergent pressure at the $+1$ defect. $R$ is the distance between the
defects, $\varphi$ is the angle of $\vec{R}=\vec{r}_{-}-\vec{r}_{+}$. $\phi$
is the uniform background orientation field. This effective interaction
induced by dipolar active forces decays inversely proportional with the
distance between defects similar to the Coulomb-like force induced by the
phase gradients alone. However, these interactions are anisotropic since they
depend on the orientation of the background polarization as well as the
orientation of the defect pair.

##  4 The effect of polar active forces

Next, we also include the polar active force in the main flow equations Eq.
(5) and study its contribution to the flow velocity $\mathbf{u}_{p}$ and
pressure $P_{p}$. In dimensionless units, the polar active force is

$\vec{F}_{p}=\tilde{\alpha}_{p}\vec{p}.$ (27)

where $\tilde{\alpha}_{p}=\alpha_{p}\tau/(\xi\Gamma)$ is a rescaled parameter
that measures the strength of polar active forces relative the frictional
drag. We consider $\alpha_{p}>0$ corresponding to polar particles moving in
the direction of their head. The corresponding flow velocity and pressure
fields can now be written compactly as (see details in SM 47)

$\displaystyle u_{p}^{+}(\hat{r},\theta)$
$\displaystyle=i\tilde{\alpha}_{p}e^{i\theta}\sin(\phi)f^{p}_{1}(\hat{r}),$
(28) $\displaystyle P_{p}^{+}(r,\theta)$
$\displaystyle=\tilde{\alpha}_{p}r\cos(\phi),$ (29) $\displaystyle
u_{p}^{-}(\hat{r},\theta)$
$\displaystyle=\frac{\tilde{\alpha}_{p}}{2}(f^{p}_{1}(\hat{r})e^{-i\theta}+f^{p}_{3}(\hat{r})e^{3i\theta}),$
(30) $\displaystyle P_{p}^{-}(r,\theta)$
$\displaystyle=\frac{\tilde{\alpha}_{p}}{3}\frac{(x^{2}-y^{2})}{r}.$ (31)

for the positive and the negative defect respectively. The velocity profiles
are plotted in the Fig. 7 while the functions $f_{1}$ and $f_{3}$ are plotted
in Fig. 8 together with their asymptotic values. It is straightforward to show
that both $f_{1}$ and $f_{3}$ tend to zero at the defect origin, thus the
polar active forces do not contribute to defect self-induced motility.

Fig. 7: Velocity streamlines due to polar active forces for a) the positive
defect (vortex flow) and b) the negative defect (saddle-point flow). Colormap
represents the magnitude of the velocity field normalized by its maximum
value.

Fig. 8: Plot of the functions (blue, whole lines) $f^{p}_{1}$ (a) and
$f^{p}_{3}$ (b) against their asymptotic limits (black, doted lines) as given
in eq. (32) and (33). Note that the functions asymptotically approach to non-
zero values in the far-field.

It is important to note that the flow induced by polar active forces around
the $+1$ defect is proportional to $\sin(\phi)$, instead of the $\sin(2\phi)$
dependence that was found for the flow induced by dipolar active forces (Eq.
16). This means that while the polar active forces do not produce any flow
around an ideal aster ($\phi=0$), there will be a polar activity-induced flow
around an ideal vortex ($\phi=\pi/2$). Furthermore, Eq. 29 shows that the
pressure is proportional to the radial distance $\sim r$, indicating that, in
the absence of any far-field screenings from other defects, the pressure can
become quite large. This is related to the flow induced by polar active forces
being constant in magnitude in the fare-field, as illustrated in Fig. 8, or by
explicitly writing the asymptotic limits

$\displaystyle u_{p}^{-}(\hat{r},\theta)$
$\displaystyle=\frac{\tilde{\alpha}_{p}}{2}\Bigg{[}\left(e^{-i\theta}+\frac{1}{3}e^{3i\theta}\right)-\frac{1}{\hat{r}^{2}}\left[e^{-i\theta}+3e^{3i\theta}\right]\Bigg{]},$
(32) $\displaystyle u_{p}^{+}(\hat{r},\theta)$
$\displaystyle=i\tilde{\alpha}_{p}e^{i\theta}\sin{\phi}\left(1-\frac{1}{\hat{r}^{2}}\right),$
(33)

for the negative and positive defect respectively. As we will show explicitly
for a pair of defects in the next section, Eq. 32-33, together with Eq. (13),
indicate that there are strong interactions between defects regardless of the
distance between them.

Having found closed form formulas for the velocity field around the defects,
the far-field vorticity induced by the positive defect can be calculated as

$\omega^{+}_{p}=\frac{\tilde{\alpha}_{p}}{\zeta}\sin{\phi}\left(\frac{1}{\hat{r}}+\frac{1}{\hat{r}^{3}}\right),$
(34)

and, for the negative defect, this is

$\omega_{p}^{-}=\frac{2\tilde{\alpha}_{p}}{\zeta}\sin{2\theta}\left(\frac{1}{\hat{r}}-\frac{3}{\hat{r}^{3}}\right).$
(35)

Interestingly the vorticity of the positive defect is independent on the polar
angle $\theta$, which is also true in the near-field, because we can write the
velocity field in the form $\vec{u}=\vec{r}^{\perp}f(r)$. As such, close to
the defect center the vorticity of the negative defect vanishes at the center,
while the positive defects has a finite vorticity

$\omega_{p}^{+}(r=0)=\frac{\tilde{\alpha}_{p}\pi}{2\zeta}\sin(\phi),$ (36)

indicating that the positive full-integer defect will be endowed with a spin.

###  4.1 Defect-pair interaction in friction-dominated system

Taking again the analytically tractable limit of zero viscosity, we can find
the velocity induced by polar active forces for a pair of oppositely-charged
defects following the approach from Sec. 3.3. The flow velocities in the
centre of the $\pm 1$ defects reduce to

$\displaystyle\vec{u}^{+}_{p}$
$\displaystyle=\frac{\tilde{\alpha}_{p}}{2|\vec{R}|^{2}}\vec{R}^{\perp}(\vec{R}^{\perp}\cdot\hat{p}_{0}),$
(37) $\displaystyle\vec{u}_{p}^{-}$
$\displaystyle=\frac{\tilde{\alpha}_{p}}{6|\vec{R}|^{2}}\vec{R}^{\perp}(\vec{R}^{\perp}\cdot\vec{p}_{0})+\frac{\alpha_{p}}{3|\vec{R}|^{2}}\vec{R}(\vec{R}\cdot\vec{p}_{0}),$
(38)

where $\vec{R}=\vec{r}^{-}-\vec{r}^{+}$ is the separation vector between the
negative defect and the positive one, and
$\vec{p}_{0}=\cos(\phi)\vec{e}_{x}+\sin(\phi)\vec{e}_{-}$ is the uniform
background polarity. Interestingly, the polar active forces induce non-
reciprocal and non-local mutual interactions that depend only on the
orientation of $\vec{R}$ relative to $\phi$, independent of the separation
distance $R$ between the defects. This suggests a truly long-ranged
interaction between oppositely-charged defects in the presence of polar active
forces.

Furthermore, the $-1$ defect tends to move towards/away from $+1$ depending on
the orientation angle of $\vec{R}$ with respect to $\vec{p}$. Both defects
move perpendicular to $\vec{R}$ but at different rates inducing the pair to
rotate. We can see this behavior more clearly by looking at the defect pair
velocity under polar active forces alone,
$\dot{\vec{R}}=\vec{u}_{p}^{-}-\vec{u}_{p}^{+}$ given by

$\frac{d}{dt}\vec{R}=\frac{\tilde{\alpha}_{p}}{3R^{2}}(\vec{R}\cdot\vec{p}_{0})\vec{R}-\frac{\tilde{\alpha}_{p}}{3R^{2}}(\vec{R}^{\perp}\cdot\vec{p}_{0})\vec{R}^{\perp},$
(39)

where the first term ($\parallel\vec{R}$) is an attraction/repulsion force
between the defects, while the second term ($\perp\vec{R}$) rotates the defect
pair. The rotation is zero, when the defect pair aligns with the background
polarization in either directions ($\vec{R}\parallel\vec{p}_{0}$). We show in
SM. 47 that the defect pair rotates until annihilation, unless $\vec{R}$ is
initially parallel to $\vec{p}_{0}$. This shows that the polar active force
leads to pair annihilation of oppositely-charged defects to promote large-
scale polar order.

To check this interesting dynamics numerically, we compare the pair trajectory
determined by Eq. (13) with that predicted by the full hydrodynamic model from
integrating Eqs. (1) (3) and (4) with the same initial configuration. The
initial uniform polarisation $\vec{p}=\vec{e}_{y}$ is seeded with a defect
pair in the $x$-direction, i.e. $\vec{r}_{-}=-20\vec{e}_{y}$ and
$\vec{r}_{+}=20\vec{e}_{x}$, such that $\vec{R}\perp\vec{p}_{0}$ at $t=0$. The
model parameters are set to $\tilde{\alpha}_{p}=0.5$, $\tilde{\alpha}_{0}=0$,
$\lambda=0$ and $\zeta^{2}=0.01$. The hydrodynamic equations are solved with
periodic boundary conditions in a domain $256\times 256$ with a spatial
discretisation $\Delta=0.5$, using spectral methods and an exponential time
differentiation scheme 52. Under the polar active force, the separation vector
$\vec{R}$ rotates relative to $\vec{p}_{0}$, and this changes the shape of the
$+1$ defect from an initial vortex to an aster at the annihilation time. This
effect is lost in the absence of polar active forces, i.e.at
$\tilde{\alpha}_{p}=0$ (see animations of the defect pair annihilation in 47).

The two trajectories are shown in Fig. 9a) and agree very well for large
$R=|\vec{R}|$ in consistency with the pointwise approximation. It also shows
that polar active forces increase the annihilation rate. This is further
evidenced in Fig. 9 b) where we compare the evolution of the separation
distance $R(t)$ when the two defects interact through the Coulomb-like forces
with or without the presence of polar active forces. At
$\tilde{\alpha}_{p}=0$, $R(t)$ decreases in the far-field as $\sqrt{t_{0}-t}$
($t_{0}$ being the annihilation time) due to the $1/R$ interaction forces.
However, for $\tilde{\alpha}_{p}\neq 0$, the annihilation timescale is greatly
reduced and the deviations from the $\sqrt{t_{0}-t}$ behavior arise due the
non-local and non-reciprocal attraction force. The rotation rate is also
contributing to aligning the defects to increase the attraction between the
defects.

Fig. 9: a) Phase portraits of the defect pair trajectory obtained from (blue)
the kinematic law for point defects (PD) from Eq. (13), and (red) by direct
simulation of the full hydrodynamic model (HM) starting from the same initial
condition. The positive defect is marked by $+$, the negative by $-$.
$\tilde{\alpha}_{p}=0.5$, $\tilde{\alpha}_{0}=0$, $\lambda=0$ and
$\zeta^{2}=0.01$. b) Distance $R$ between a defect pair as a function of time
for Coulomb-like interactions only (blue curve) and in the presence of polar
active forces (orange curve). The black doted lines are fitting curves with
the $R(t)=\sqrt{At+B}$ found by linear regression on $R^{2}$.

##  5 Discussion/conclusion

In summary, we present theoretical derivations of incompressible flow fields
and dynamics of $\pm 1$ defects under dipolar and polar active forces. These
active forces do not endow the full-integer defects with any self-propulsion,
as expected from the symmetry of the defects. However, the $+1$ defect
acquires a non-zero active torque due to both polar and dipolar active forces.
We show that the strength of this spin is dependent on whether the defect is a
vortex, aster or spiral. In the absence of hydrodynamic screening due to
friction, the vortical active flow around a $+1$ vortex changes sign on a
length scale set by the coherence length as observed in numerical simulations
and predicted analytically. For both defects, the dipolar active force
contribution to the flow field vanishes with $1/r$ in the defect far-field,
whereas the polar active force contribution approach a constant value.
Remarkably, polar active forces mediate mutual interaction between oppositely-
charged defect pairs in a manner that renders defect-defect interactions
independent of the distance between the defect pair. We have shown that, under
these polar active forces, a pair of oppositely-charged defects rotates to
align with the background polarization field, while the negative defect chases
the positive one until annihilation. The rate of annihilation is greatly
enhanced by polar active forces, and this is the main underlying mechanism for
the suppression of defect-laden active turbulence in polar active matter as
reported numerically in Ref. 34.

## Author Contributions

J.Rø contributed with analytical derivations, numerical simulations,
visualisation and data analysis. J.Re. contributed with lattice-Boltzmann
simulations, visualisation and data analysis. A.D. and L.A. formulated the
project and contributed with supervision. A.D. provided computational
resources. L.A. checked the analytical calculations. All authors contributed
with drafting and writing the manuscript.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

J.Rø. and L.A. acknowledge support from the Research Council of Norway through
the Center of Excellence funding scheme, Project No. 262644 (PoreLab). A. D.
acknowledges funding from the Novo Nordisk Foundation (grant No.
NNF18SA0035142 and NERD grant No. NNF21OC0068687), Villum Fonden Grant no.
29476, and the European Union via the ERC-Starting Grant PhysCoMeT. Views and
opinions expressed are however those of the authors only and do not
necessarily reflect those of the European Union or the European Research
Council. Neither the European Union nor the granting authority can be held
responsible for them.

## Appendix:

Here we write down the radially-dependent functions associated with the flow
fields induced by polar and dipolar active forces. Derivation details are
found in the SM. 47. The radial functions for the dipolar active flow velocity
are:

$\displaystyle f_{3}^{a}(\hat{r})$
$\displaystyle=\frac{1}{\hat{r}}[1-\hat{r}K_{1}(\hat{r})]-\frac{2}{\hat{r}^{3}}[4-\hat{r}^{3}K_{1}(\hat{r})-2\hat{r}^{2}K_{2}(\hat{r})],$
(40) $\displaystyle f_{5}^{a}(\hat{r})$
$\displaystyle=\frac{1}{\hat{r}}[1-\hat{r}K_{1}(\hat{r})]-\frac{6}{\hat{r}^{3}}[4-\hat{r}^{3}K_{1}(\hat{r})-2\hat{r}^{2}K_{2}(\hat{r})]+\frac{6}{\hat{r}^{5}}\left[64-4\hat{r}^{2}(8+\hat{r}^{2})K_{0}(\hat{r})-\hat{r}(8+\hat{r}^{2})^{2}K_{1}(\hat{r})\right].$
(41)

The radial function for the corresponding vorticity field

$f_{\omega}^{a}(\hat{r})=\frac{1}{\hat{r}^{4}}\left(-48+4\hat{r}^{2}+\hat{r}^{2}(24+\hat{r}^{2})K_{0}(\hat{r})+8\hat{r}(6+\hat{r}^{2})K_{1}(\hat{r})\right).$
(42)

The radial functions for the polar active flow velocity are:

$f^{p}_{1}(\hat{r})=1-\hat{r}K_{1}(\hat{r})+\frac{\pi}{2}I_{1}(\hat{r})+\sum_{k,n=0}^{\infty}\kappa_{1}(n,k)\frac{\hat{r}^{2k+2}}{((2k)!!)^{2}}+\sum_{k,n=0}^{\infty}\kappa_{2}(n,k)\left[\ln\left(\frac{\hat{r}}{2}\right)-\psi^{(0)}(k+1)\right]\frac{\hat{r}^{2k+2}}{((2k)!!)^{2}},$
(43)

and

$f^{p}_{3}(\hat{r})=-\frac{\pi}{2}I_{3}(\hat{r})+\sum_{k,n=0}^{\infty}\kappa_{3}(n,k)\frac{\hat{r}^{2k+2}}{((2k)!!)^{2}}+\sum_{k,n=0}^{\infty}\kappa_{4}(n,k)\frac{\hat{r}^{2k+2}}{((2k)!!)^{2}}\left[\ln\left(\frac{\hat{r}}{2}\right)-\psi^{(0)}(k+1)\right].$
(44)

Where $\psi^{(0)}$ is the digamma function. The coefficients of the infinite
power series

$\displaystyle\kappa_{1}(n,k)$
$\displaystyle=-\left(\frac{(2n-1)!!}{(2n)!!}\right)^{2}\frac{1}{2n+2}\left(\frac{1}{(2n-1-2k)^{2}}+\frac{1}{4}\frac{2n+1}{(2n+2)(n+2+k)^{2}}\right)$
(45) $\displaystyle\kappa_{2}(n,k)$
$\displaystyle=\left(\frac{(2n-1)!!}{(2n)!!}\right)^{2}\frac{1}{2n+2}\left(\frac{2n+1}{2(2n+2)(n+2+k)}-\frac{1}{2n-1-2k}\right),$
(46) $\displaystyle\kappa_{3}(n,k)$
$\displaystyle=\left(\frac{(2n-1)!!}{(2n)!!}\right)^{2}\left(\frac{2n+1}{4(3-2n)(1-2n)(n+k+1)^{2}}+\frac{2n}{(2n+4)(2n+2)(2n-1-2k)^{2}}\right),$
(47) $\displaystyle\kappa_{4}(n,k)$
$\displaystyle=\left(\frac{(2n-1)!!}{(2n)!!}\right)^{2}\left(\frac{2n}{(2n+4)(2n+2)(2n-1-2k)}-\frac{2n+1}{2(3-2n)(1-2n)(n+k+1)}\right).$
(48)

Note that when summing over $n$ the coefficients $\kappa_{4}(n,k)$ tends to
zero for any $k$. The infinite series are convergent, but slowly.

## Notes and references

  * Doostmohammadi _et al._ 2018 A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans and F. Sagués, _Nat. Commun._ , 2018, 9, 1–13.
  * Marchetti _et al._ 2013 M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, _Rev. Mod. Phys._ , 2013, 85, 1143.
  * Blanch-Mercader _et al._ 2018 C. Blanch-Mercader, V. Yashunsky, S. Garcia, G. Duclos, L. Giomi and P. Silberzan, _Physical review letters_ , 2018, 120, 208101\.
  * Lin _et al._ 2021 S.-Z. Lin, W.-Y. Zhang, D. Bi, B. Li and X.-Q. Feng, _Communications Physics_ , 2021, 4, 1–9.
  * Sanchez _et al._ 2012 T. Sanchez, D. T. Chen, S. J. DeCamp, M. Heymann and Z. Dogic, _Nature_ , 2012, 491, 431–434.
  * Guillamat _et al._ 2017 P. Guillamat, J. Ignés-Mullol and F. Sagués, _Nature communications_ , 2017, 8, 1–8.
  * Needleman and Dogic 2017 D. Needleman and Z. Dogic, _Nat. Rev. Mat._ , 2017, 2, 17048.
  * Kumar _et al._ 2018 N. Kumar, R. Zhang, J. J. de Pablo and M. L. Gardel, _Sci. Adv._ , 2018, 4, eaat7779.
  * Kudrolli _et al._ 2008 A. Kudrolli, G. Lumay, D. Volfson and L. S. Tsimring, _Physical review letters_ , 2008, 100, 058001.
  * Kruse _et al._ 2005 K. Kruse, J.-F. Joanny, F. Jülicher, J. Prost and K. Sekimoto, _The European Physical Journal E_ , 2005, 16, 5–16.
  * Toner and Tu 1995 J. Toner and Y. Tu, _Phys. Rev. Lett._ , 1995, 75, 4326.
  * Juelicher _et al._ 2007 F. Juelicher, K. Kruse, J. Prost and J.-F. Joanny, _Physics reports_ , 2007, 449, 3–28.
  * Giomi _et al._ 2013 L. Giomi, M. J. Bowick, X. Ma and M. C. Marchetti, _Phys. Rev. Lett._ , 2013, 110, 228101.
  * Ramaswamy _et al._ 2003 S. Ramaswamy, R. A. Simha and J. Toner, _Europhy. Lett. (EPL)_ , 2003, 62, 196.
  * Kruse _et al._ 2004 K. Kruse, J.-F. Joanny, F. Jülicher, J. Prost and K. Sekimoto, _Physical review letters_ , 2004, 92, 078101.
  * Simha and Ramaswamy 2002 R. A. Simha and S. Ramaswamy, _Physical review letters_ , 2002, 89, 058101.
  * Thampi _et al._ 2014 S. P. Thampi, R. Golestanian and J. M. Yeomans, _Philos. Trans. Royal Soc. A_ , 2014, 372, 20130366.
  * Doostmohammadi _et al._ 2017 A. Doostmohammadi, T. N. Shendruk, K. Thijssen and J. M. Yeomans, _Nature communications_ , 2017, 8, 1–7.
  * Giomi _et al._ 2013 L. Giomi, M. J. Bowick, X. Ma and M. C. Marchetti, _Physical review letters_ , 2013, 110, 228101.
  * Rønning _et al._ 2022 J. Rønning, C. M. Marchetti, M. J. Bowick and L. Angheluta, _Proceedings of the Royal Society A_ , 2022, 478, 20210879.
  * Chaikin and Lubensky 1995 P. M. Chaikin and T. C. Lubensky, in _Frontmatter_ , Cambridge University Press, 1995, pp. i–vi.
  * Doostmohammadi and Ladoux 2021 A. Doostmohammadi and B. Ladoux, _Trends in Cell Biology_ , 2021.
  * Shankar _et al._ 2022 S. Shankar, A. Souslov, M. J. Bowick, M. C. Marchetti and V. Vitelli, _Nature Reviews Physics_ , 2022, 4, 380–398.
  * DeCamp _et al._ 2015 S. J. DeCamp, G. S. Redner, A. Baskaran, M. F. Hagan and Z. Dogic, _Nat. Mater._ , 2015, 14, 1110–1115.
  * Duclos _et al._ 2018 G. Duclos, C. Blanch-Mercader, V. Yashunsky, G. Salbreux, J.-F. Joanny, J. Prost and P. Silberzan, _Nat. Phys._ , 2018, 14, 728–732.
  * Meacock _et al._ 2021 O. J. Meacock, A. Doostmohammadi, K. R. Foster, J. M. Yeomans and W. M. Durham, _Nature Phys._ , 2021, 17, 205–210.
  * Saw _et al._ 2017 T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, _Nature_ , 2017, 544, 212–216.
  * Duclos _et al._ 2017 G. Duclos, C. Erlenkämper, J.-F. Joanny and P. Silberzan, _Nat. Phys._ , 2017, 13, 58–62.
  * Malinverno _et al._ 2017 C. Malinverno, S. Corallino, F. Giavazzi, M. Bergert, Q. Li, M. Leoni, A. Disanza, E. Frittoli, A. Oldani, E. Martini _et al._ , _Nature materials_ , 2017, 16, 587–596.
  * Baskaran and Marchetti 2008 A. Baskaran and M. C. Marchetti, _Physical Review E_ , 2008, 77, 011920\.
  * Baskaran and Marchetti 2010 A. Baskaran and M. C. Marchetti, _Journal of Statistical Mechanics: Theory and Experiment_ , 2010, 2010, P04019.
  * Baskaran and Marchetti 2012 A. Baskaran and M. C. Marchetti, _The European Physical Journal E_ , 2012, 35, 1–8.
  * Amiri _et al._ 2022 A. Amiri, R. Mueller and A. Doostmohammadi, _Journal of Physics A: Mathematical and Theoretical_ , 2022, 55, 094002.
  * Andersen _et al._ 2022 B. H. Andersen, J. Renaud, J. Rønning, L. Angheluta and A. Doostmohammadi, _arXiv preprint arXiv:2209.10916_ , 2022.
  * Ardaševa and Doostmohammadi 2022 A. Ardaševa and A. Doostmohammadi, _Nature Reviews Physics_ , 2022, 4, 354–356.
  * Guillamat _et al._ 2022 P. Guillamat, C. Blanch-Mercader, G. Pernollet, K. Kruse and A. Roux, _Nature materials_ , 2022, 21, 588–597.
  * Blanch-Mercader _et al._ 2021 C. Blanch-Mercader, P. Guillamat, A. Roux and K. Kruse, _Physical Review E_ , 2021, 103, 012405.
  * Angheluta _et al._ 2021 L. Angheluta, Z. Chen, M. C. Marchetti and M. J. Bowick, _New Journal of Physics_ , 2021, 23, 033009.
  * Chandragiri _et al._ 2019 S. Chandragiri, A. Doostmohammadi, J. M. Yeomans and S. P. Thampi, _Soft matter_ , 2019, 15, 1597–1604.
  * Doostmohammadi _et al._ 2016 A. Doostmohammadi, M. F. Adamer, S. P. Thampi and J. M. Yeomans, _Nat. Commun._ , 2016, 7, 1–9.
  * Tang and Selinger 2017 X. Tang and J. V. Selinger, _Soft matter_ , 2017, 13, 5481–5490.
  * Mazenko 1997 G. F. Mazenko, _Physical review letters_ , 1997, 78, 401.
  * Halperin 1981 B. Halperin, _published in Physics of Defects, proceedings of Les Houches, Session XXXV 1980 NATO ASI, editors Balian, Kléman and Poirier_ , 1981\.
  * Mazenko 2001 G. F. Mazenko, _Physical Review E_ , 2001, 64, 016110.
  * Pismen _et al._ 1999 L. M. Pismen _et al._ , _Vortices in nonlinear fields: From liquid crystals to superfluids, from non-equilibrium patterns to cosmic strings_ , Oxford University Press, 1999, vol. 100.
  * Vafa _et al._ 2020 F. Vafa, M. J. Bowick, M. C. Marchetti and B. I. Shraiman, _arXiv:2007.02947_ , 2020.
  * 47 _Supplementary Materials_ , https://github.com/jonasron/SM_Flowfield_Polar.
  * Giomi _et al._ 2014 L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek and M. Cristina Marchetti, _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 2014, 372, 20130365.
  * Sonam _et al._ 2022 S. Sonam, L. Balasubramaniam, S.-Z. Lin, Y. M. Y. Ivan, I. Pi-Jaumà, C. Jebane, M. Karnat, Y. Toyama, P. Marcq, J. Prost _et al._ , _Nature Physics_ , 2022, 1–10.
  * Makhija _et al._ 2022 E. Makhija, Y. Zheng, J. Wang, H. R. Leong, R. B. Othman, E. X. Ng, E. H. Lee, L. Tucker-Kellogg, Y. H. Lee, H. Yu _et al._ , _bioRxiv_ , 2022.
  * Thampi _et al._ 2014 S. P. Thampi, R. Golestanian and J. M. Yeomans, _EPL (Europhysics Letters)_ , 2014, 105, 18001.
  * Cox and Matthews 2002 S. M. Cox and P. C. Matthews, _Journal of Computational Physics_ , 2002, 176, 430–455.