# QuantumSkynet: A High-Dimensional Quantum Computing Simulator

Andres Giraldo-Carvajal1, Daniel A. Duque-Ramirez1,
Jose A. Jaramillo-Villegas1,2 1 Universidad Tecnológica de Pereira, Pereira,
Risaralda, Colombia
2 Laboratory for Research in Complex Systems, Menlo Park, California, USA
<EMAIL_ADDRESS>
###### Abstract

The use of classical computers to simulate quantum computing has been
successful in aiding the study of quantum algorithms and circuits that are too
complex to examine analytically. Current implementations of quantum computing
simulators are limited to two-level quantum systems. Recent advances in high-
dimensional quantum computing systems have demonstrated the viability of
working with multi-level superposition and entanglement. These advances allow
an agile increase in the number of dimensions of the system while maintaining
quantum entanglement, achieving higher encoding of information and making
quantum algorithms less vulnerable to decoherence and computational errors. In
this paper, we introduce QuantumSkynet, a novel high-dimensional cloud-based
quantum computing simulator. This platform allows simulations of qudit-based
quantum algorithms. We also propose a unified generalization of high-
dimensional quantum gates, which are available for simulations in
QuantumSkynet. Finally, we report simulations and their results for qudit-
based versions of the Deutsch–Jozsa and quantum phase estimation algorithms
using QuantumSkynet.

  * Septiembre 2020

##  1 Introduction

Quantum computing has made significant progress in recent years. The race to
build a computer capable of solving problems that classical computers cannot
is reaching an intense level [1, 2].

Cryptography [3, 4], drug design and discovery [5, 6, 7], the study of complex
molecular structures [8, 9], and optimization and financial forecasting [10,
11] are some areas that hope to benefit from the construction of a quantum
computer. However, today’s quantum computers are highly susceptible to noise,
which can lead to errors in the calculations they perform [12].

An emerging alternative that allows increased resistance to noise is high-
dimensional quantum computing [13]. Working with quantum states that handle
multiple dimensions not only enables systems that are more robust to noise but
also more secure systems with greater information capacity [14].

Different physical platforms can be used to implement high-dimensional quantum
systems [15], such as photons [16, 17, 18, 13, 19, 20], trapped ions [21, 22],
superconducting systems [23, 24], and molecules [25, 26, 27]. The practical
application of these qudit-based systems is a very active research field with
great potential.

One tool that is instrumental in the progress in this area is quantum
computing simulators. With these simulators, it is possible to study quantum
circuits and measure their performance in the presence of noise, which helps
make decisions that can improve the design of quantum hardware [28].

Currently, different quantum computing simulators and frameworks, such as
qHiPSTER [28], Quantum Composer [29], Strawberry Fields [30], Rigetti
Computing [31], Libquantum [32], and Quirk [33], among others, allow the
execution of quantum algorithms and circuits to be efficiently simulated based
on qubits; however, there is currently no quantum computing tool that allows
simulation for high-dimensional quantum circuits.

In this paper, we propose QuantumSkynet, the first high-dimensional quantum
computing simulator, which will enable scientists to evaluate, analyze, and
tune high-dimensional quantum algorithms in a cloud-based environment.
Additionally, we propose a generalization for high-dimensional quantum gates,
where we show a relationship between Weyl’s adjoint operator and high-
dimensional Pauli gates. Furthermore, we propose a more generalized equation
to describe high-dimensional Fourier quantum transform.

##  2 High-dimensional quantum gates

During the research and development of QuantumSkynet, a compendium and
generalization of quantum gates were employed to build arbitrary high-
dimensional quantum algorithms, as described in the following.

###  2.1 Weyl operators

Weyl operators are a type of gate that apply to a single qudit. They can be
used to create states of maximum entanglement [34] [35]. These unitary
operators form an orthonormal basis of a $d$-dimensional Hilbert–Schmidt space
and are represented by $d^{2}$ matrices generated from the following
expression:

$W_{p,q}=\sum_{k=0}^{d-1}\omega^{k\hskip 0.5ptp}\Ket{k}\Bra{k\oplus q}$ (1)

where $k\oplus q$ is equal to $(k+q)\,mod\,d$, $\omega$ and is defined by
$e^{(\frac{2\pi}{d})i}$, and $p,q$ are subscripts that take values from the
standard base of the Hilbert space of dimension $d$.

###  2.2 Weyl’s adjoint operator

One of the analyses conducted in this study involved establishing a
correspondence between Weyl’s adjoint operator,

$W_{q,p}^{\dagger}=\sum_{k=0}^{d-1}\omega^{k\hskip 0.5ptq}\Ket{k\oplus
p}\Bra{k}$ (2)

and the $I,X,Y$, and $Z$ gates in high dimensions:

$\begin{gathered}I^{m}=W_{0,0}^{\dagger}\\\\[10.00002pt]
X^{m}=W_{0,m}^{\dagger}\\\\[10.00002pt]
Z^{m}=W_{m,0}^{\dagger}\\\\[10.00002pt]
Y^{m}=i^{(m\%d)}\;W_{m,m}^{\dagger}\\\\[10.00002pt] \end{gathered}$ (3)

In this manner, it is possible to find the $I,X,Y$, and $Z$ gates for any
dimension raised to any exponent. The identity operator for high dimensions
can be found using equation (2) by replacing (q,p) with (0,0):

$\begin{gathered}I^{m}=W_{0,0}^{\dagger}\\\\[10.00002pt]
=\sum_{k=0}^{d-1}\omega^{0}\Ket{k\oplus 0}\Bra{k}\\\\[10.00002pt]
=\sum_{k=0}^{d-1}\Ket{k}\Bra{k}\\\\[10.00002pt]
=\sum_{k=0}^{d-1}P_{k}\\\\[10.00002pt] \end{gathered}$ (4)

where $P_{k}$ is known as the projection operator, which is equal to
$\Ket{k}\Bra{k}$.

By replacing (q,p) with (0,m) in equation (2), it is possible to find the
expression that generalizes the $X$ Pauli matrix for high dimensions and
raising to any exponent:

$\begin{gathered}X^{m}=W_{0,m}^{\dagger}\\\\[10.00002pt]
=\sum_{k=0}^{d-1}\omega^{0}\Ket{k\oplus m}\Bra{k}\\\\[10.00002pt]
=\sum_{k=0}^{d-1}\Ket{k\oplus m}\Bra{k}\\\\[10.00002pt] \end{gathered}$ (5)

To find the generalized form of the $Z$ Pauli matrix, simply replace (q,p)
with (m,0) in equation (2):

$\begin{gathered}Z^{m}=W_{m,0}^{\dagger}=\sum_{k=0}^{d-1}\omega^{k\hskip
0.5ptm}P_{k}\\\\[10.00002pt] =\sum_{k=0}^{d-1}\omega^{k\hskip
0.5ptm}\Ket{k\oplus 0}\Bra{k}\\\\[10.00002pt] =\sum_{k=0}^{d-1}\omega^{k\hskip
0.5ptm}\Ket{k}\Bra{k}\\\\[10.00002pt] =\sum_{k=0}^{d-1}\omega^{k\hskip
0.5ptm}P_{k}\\\\[10.00002pt] \end{gathered}$ (6)

The $Y$ matrix can be viewed as a combination of the $X$ and $Z$ matrices. It
is possible to find its generalized equation for high dimensions and raising
to any exponent by using equation (2) again, replacing (q,p) with (m,m) and
multiplying by a cyclical global phase of $i^{(m\%d)}$:

$\begin{gathered}Y^{m}=i^{(m\%d)}\;W_{m,m}^{\dagger}\\\\[10.00002pt]
Y^{m}=i^{(m\%d)}\;\sum_{k=0}^{d-1}\omega^{k\hskip 0.5ptm}\Ket{k\oplus
m}\Bra{k}\end{gathered}$ (7)

In this manner, we find the $X^{m}$, $Y^{m}$, and $Z^{m}$ gates generalized to
$d$ dimensions of the qubit-based $X$, $Y$, and $Z$ gates, respectively,
raised to the $m$ power. Below is a summary of the gate generalizations
described thus far:

Figure 1: Single-qudit high-dimensional gates $X^{m}$, $Y^{m}$, $Z^{m}$, and
$W_{p,q}$.
$P_{k}$ is the projection operator equal to $\Ket{k}\Bra{k}$, and $\omega$ is
defined by $e^{(\frac{2\pi}{d})i}$

###  2.3 Controlled gate

The $CNOT$ gate, which acts on qubits, can be generalized to high dimensions
by the $SUM$ gate, which is defined by the following expression:

$SUM=\sum_{k=0}^{d-1}P_{k}\otimes X^{k}$ (8)

This gate can also be seen as a particular case of the arbitrary generalized
controlled gate:

$CU=\sum_{k=0}^{d-1}P_{k}\otimes U^{k}$ (9)

where $U$ can be $X,Y,Z$, or any high-dimensional operator.

In this manner, the high-dimensional controlled gates available in
QuantumSkynet are as follows:

Figure 2: Multi-qudit generalized controlled gate and $SUM$ gate.
U is an arbitrary unit matrix, and $P_{k}$ is the projection operator equal to
$\Ket{k}\Bra{k}$

###  2.4 Generalized quantum Fourier transform

The classic discrete Fourier transform maps a series of $N$ complex numbers
into another series of complex numbers as follows [36]:

$y_{k}=\frac{1}{\sqrt{N}}\sum_{j=0}^{N-1}\omega_{N}^{k\hskip
0.5ptj}x_{j}\;,\quad k=0,1,...,N-1$ (10)

where
$x=[x_{0},...,x_{N}-1]^{T}\in\mathbb{C}^{N}\;,\;y=[y_{0},...,y_{N}-1]^{T}\in\mathbb{C}^{N}\;,\;\omega_{N}=e^{(\frac{2\pi}{N})i}$.

Similarly, the quantum Fourier transform acts on a qudit with quantum state of
dimension $d$ (where $N=d$) and maps it to another $d$-dimensional quantum
state according to the following expression:

$\Ket{j}\xrightarrow{\mathscr{F}}\frac{1}{\sqrt{d}}\sum_{k=0}^{d-1}\omega^{kj}\Ket{k}$
(11)

where $j$ is a quantum state belonging to the standard base of a
$d$-dimensional Hilbert space with $\omega=e^{\frac{2\pi i}{d}}$.

When $d$ is equal to 2, the quantum Fourier transform is the Hadamard gate
acting on 1 or $N$ qubits. $d$ greater than 2 gives a generalized Hadamard
gate that acts on 1 or $N$ qudits, which is known as $DFT$ in some papers [18,
19].

We grouped all different combinations in which the quantum Fourier transform
can be applied for any number of qudits and any number of dimensions, which we
call a Generalized $QFT$ ($GQFT$). For this, it is necessary to cover all
elements of the quantum state resulting from the interaction between these
qudits. That is, the size $N$ of the vector that will be affected by the
$GQFT$ will be $d^{n}$, which is the same size as the tensor product between
the qudits affected by the gate:

$\Ket{j}\xrightarrow{\mathscr{F}}\frac{1}{\sqrt{d^{n}}}\sum_{k=0}^{d^{n}-1}\omega^{kj}\Ket{k}$
(12)

where $n$ refers to the number of qudits of $d$ dimensions that will be
affected by the gate.

The following figure shows the proposed $GQFT$.

Figure 3: Proposed GQFT

As can be seen in the previous figure, the Hadamard gate for the case of
qubits is a special case of the $QFT$ gate, while for the case of qudits, the
$GH$ or Generalized Hadamard gate (commonly misnamed $DFT$) is a special case
of the $GQFT$ gate.

Therefore, this is another gate available for use with QuantumSkynet:

Figure 4: GQFT gate

##  3 Methods

QuantumSkynet was written using the C++ programming language. It is divided
into two core classes, Qudit and QRegister, as described below.

###  3.1 Qudit class

Consider the following definition of a qudit:

$\Ket{\Psi}=\sum_{k=0}^{d-1}\alpha_{k}\Ket{k}$ (13)

where $d$ is the number of dimensions and $\alpha_{k}$ are the complex
amplitudes that accompany the base states of the qudit. A class for the
construction of qudits was implemented in the simulator as follows:

⬇

#include <vector>

#include <complex>

class Qudit {

private:

vector<complex<double>> amplitudes;

public:

Qudit(vector<complex<double>> coefficients) {

amplitudes = coefficients;

};

vector<complex<double>> getAmplitudes() {

return amplitudes;

};

};

This class receives a vector of complex values that correspond to the
amplitudes $\alpha_{k}$ that make up the quantum state of the qudit. This
vector is assigned to the private attribute amplitudes, which can be accessed
through the getAmplitudes() function.

###  3.2 QRegister class

To store and process the input qudits and allow their interaction between
themselves and other functions in the quantum circuit, the QRegister class was
created to define the structure of the quantum register:

⬇

class QRegister {

private:

int quditsDimensions;

int quditsNumber;

vector<Qudit> initialQudits;

vector<complex<double>> statesVector;

public:

QRegister(

int dimensions,

vector<vector<complex<double>>> initialQuditsState

) {

quditsDimensions = dimensions;

quditsNumber = initialQuditsState.size();

if (quditsNumber == 1) {

Qudit qudit(initialQuditsState[0]);

initialQudits.push_back(qudit);

statesVector = initialQuditsState[0];

} else {

for (int i=0; i<quditsNumber; i++) {

Qudit qudit(initialQuditsState[i]);

initialQudits.push_back(qudit);

}

statesVector = tensorProduct(

dimensions,

initialQudits

);

}

statesVector = normalize(statesVector);

}

int getDimensions() {

return quditsDimensions;

};

int getQuditsNumber() {

return quditsNumber;

};

vector<Qudit> getInitialQudits() {

return initialQudits;

};

vector<complex<double>> getStatesVector() {

return statesVector;

};

};

The $statesVector$ variable defined within the QRegister class stores the
vector of amplitudes with the global quantum state of the circuit. This
variable is modified by the application of the high-dimensional quantum gates
supported by the simulator.

###  3.3 Cloud-based architecture

To allow access to the different functions of the simulator, cloud computing
services provided by Amazon Web Services (AWS) were used. The following
diagram shows the set of AWS-specific technologies that allowed the simulator
to be deployed on the web:

Figure 5: Cloud-based architecture

As shown in the previous diagram, this project was deployed in the form of an
application programming interface (API) to be used as a gateway to the data
and functionalities of the simulator.

A Cloud9 instance was used to debug the code and configure the simulator
deployment on the cloud. Cloud9 automatically creates an EC2 instance and
manages the lifecycle of this instance as required.

From the EC2 instance, a Lambda function was created in C++ with the logic to
process the requests received by the client and return a response from the
simulator. An API was then created using the Amazon API Gateway service and
integrated with the pre-built Lambda function.

To request information from the API, an HTTP request is sent from the client
with the quantum circuit to be simulated in JSON format in the body of the
request, and a response that also arrives in JSON format is expected.

##  4 Results and discussion

###  4.1 Quantum phase estimation algorithm

Quantum phase estimation is one of the most important subroutines in quantum
computation [37]. This algorithm is a central part of many other quantum
algorithms, including the Shor factoring algorithm.

The high-dimensional version of this algorithm, which was presented in [18],
was studied and implemented in QuantumSkynet. The quantum circuit associated
with the phase estimation algorithm (PEA) is shown below.

Figure 6: Phase estimation algorithm

The $DFT$ gates are qudit Hadamard gates, which can be represented by the
$GQFT$ gate for a qudit in QuantumSkynet, and the $U$ gate used to test this
algorithm was the same as that was used as a proof of concept in [18], the
generalized 3-dimensional $Z$ gate.

The eigen-phases that correspond to each of the eigen-states of the
generalized gate $Z$ are as follows:

EigenState | $\Ket{0}$ | $\Ket{1}$ | $\Ket{2}$
---|---|---|---
True Phase, $\phi$ | 0 | $\frac{2\pi}{3}$ | $\frac{4\pi}{3}$

The simulation of the above circuit in QuantumSkynet was conducted as follows:

Figure 7: PEA in QuantumSkynet

where $\Ket{\Psi}$ represents each of the eigen-states
$\Ket{0},\Ket{1}y\Ket{2}$ required to find the respective eigen-phases.

Below are the respective resulting states based on the vector of amplitudes
and value $\tau$ equal to the measured qudit value with index 0 based on the
vector of probabilities, calculated by QuantumSkynet for each eigen-state (the
qubit at the bottom of the circuit is the least significant):

EigenState | $\Ket{0}$ | $\Ket{1}$ | $\Ket{2}$
---|---|---|---
Result | $\Ket{00}$ | $\Ket{11}$ | $\Ket{22}$
$\tau$ | 0 | 1 | 2

The phase $\phi$ is equal to $2\pi\tau/3$. Therefore, the values of the phases
found by QuantumSkynet for the respective eigen-states of the 3-dimensional
$Z$ matrix were as follows:

EigenState | $\Ket{0}$ | $\Ket{1}$ | $\Ket{2}$
---|---|---|---
True Phase, $\phi$ | 0 | $\frac{2\pi}{3}$ | $\frac{4\pi}{3}$

which are effectively the values of the phases presented in [18].

###  4.2 Deutsch–Jozsa algorithm

The Deutsch–Jozsa algorithm is a remarkable example that demonstrates quantum
speed-up. This algorithm can determine whether a function $f(x)$ is constant
or balanced using only one function evaluation.

Figure 8: Deutsch–Jozsa algorithm

In the case of qudits, the Hadamard gate is replaced by the $GQFT$ gate, but
the analysis to identify a function as constant or balanced is similar to that
for the two-dimensional case; if each x-register qudit returns 0, it is a
constant function; otherwise, the function is balanced [15].

As an example, the following high-dimensional quantum circuits were tested on
QuantumSkynet to evaluate whether the functions $f(x)$ are constant or
balanced:

Figure 9: Deutsch–Jozsa algorithm on QuantumSkynet

For the previous circuit, the following probability vector was obtained for
each qudit:

| Chance result
---|---
| $\Ket{0}$ | $\Ket{1}$ | $\Ket{2}$
q[0] | 1 | 0 | 0
q[1] | 1 | 0 | 0
q[2] | 0.33 | 0.33 | 0.33

The probability of obtaining $\Ket{0}$ for the qudits at the top of the
circuit was $100\%$, which means that the function is constant. A similar
occurrence is seen with the following circuit, where all the probability
values for the two least significant qudits of the circuit are equal to 1:

Figure 10: Deutsch–Jozsa algorithm on QuantumSkynet

| Chance result
---|---
| $\Ket{0}$ | $\Ket{1}$ | $\Ket{2}$ | $\Ket{3}$ | $\Ket{4}$
q[0] | 1 | 0 | 0 | 0 | 0
q[1] | 1 | 0 | 0 | 0 | 0
q[2] | 0.2 | 0.2 | 0.2 | 0.2 | 0.2

Now we show two examples of quantum circuits in which the function $f(x)$ is
balanced:

Figure 11: Deutsch–Jozsa algorithm on QuantumSkynet

| Chance result
---|---
| $\Ket{0}$ | $\Ket{1}$ | $\Ket{2}$
q[0] | 0 | 1 | 0
q[1] | 0 | 1 | 0
q[2] | 0 | 1 | 0
q[3] | 0.33 | 0.33 | 0.33

Because the three least significant qudits have zero probability of being at
$\Ket{0}$, the function $f(x)$ is balanced.

Figure 12: Deutsch–Jozsa algorithm on QuantumSkynet

| Chance result
---|---
| $\Ket{0}$ | $\Ket{1}$ | $\Ket{2}$ | $\Ket{3}$
q[0] | 0 | 0 | 1 | 0
q[1] | 0.25 | 0.25 | 0.25 | 0.25

In this last example, the function $f(x)$ is also balanced because the least
significant qudit has no chance of being $\Ket{0}$.

##  5 Conclusion

Due to the advantages of qudit-based circuits over qubit-based ones, high-
dimensional quantum computing simulation will continue to be an area of high
interest.

The calculations derived from the analysis of quantum circuits based on qudits
can become complex and grow as the number of dimensions increases. In this
study, we explored a very useful tool to easily and intuitively simulate high-
dimensional quantum circuits.

QuantumSkynet opens to the possibility of helping extend qubit-based quantum
circuits to high-dimensional circuits as well as exploring and analyzing new
qudit-based quantum circuits.

As future work, we intend to add support for the distinction between pure and
mixed states and allow the simulation of different levels of quantum
decoherence.

## References

  * [1] Mario Piattini, Guido Peterssen and Ricardo Pérez-Castillo  “Quantum Computing: A New Software Engineering Golden Age”  In _ACM SIGSOFT Software Engineering Notes_ 45.3 ACM New York, NY, USA, 2020, pp. 12–14
  * [2] Gideon Lichfield  “Inside the race to build the best quantum computer on Earth”  URL: https://www.technologyreview.com/2020/02/26/916744/quantum-computer-race-ibm-google/
  * [3] Dan Boneh and Mark Zhandry  “Secure signatures and chosen ciphertext security in a quantum computing world”  In _Annual Cryptology Conference_ , 2013, pp. 361–379  Springer
  * [4] Kevin Marshall et al.  “Continuous-variable quantum computing on encrypted data”  In _Nature communications_ 7.1 Nature Publishing Group, 2016, pp. 1–7
  * [5] Yudong Cao, Jhonathan Romero and Alán Aspuru-Guzik  “Potential of quantum computing for drug discovery”  In _IBM Journal of Research and Development_ 62.6 IBM, 2018, pp. 6–1
  * [6] Saad M Darwish, Tamer A Shendi and Ahmed Younes  “Chemometrics approach for the prediction of chemical compounds’ toxicity degree based on quantum inspired optimization with applications in drug discovery”  In _Chemometrics and Intelligent Laboratory Systems_ 193 Elsevier, 2019, pp. 103826
  * [7] Ting Zhou, Danzhi Huang and Amedeo Caflisch  “Quantum mechanical methods for drug design”  In _Current topics in medicinal chemistry_ 10.1 Bentham Science Publishers, 2010, pp. 33–45
  * [8] Yudong Cao et al.  “Quantum chemistry in the age of quantum computing”  In _Chemical reviews_ 119.19 ACS Publications, 2019, pp. 10856–10915
  * [9] Rongxin Xia and Sabre Kais  “Qubit coupled cluster singles and doubles variational quantum eigensolver ansatz for electronic structure calculations.”  In _Quantum Science and Technology_ IOP Publishing, 2020
  * [10] Roman Orus, Samuel Mugel and Enrique Lizaso  “Quantum computing for finance: overview and prospects”  In _Reviews in Physics_ 4 Elsevier, 2019, pp. 100028
  * [11] Román Orús, Samuel Mugel and Enrique Lizaso  “Forecasting financial crashes with quantum computing”  In _Physical Review A_ 99.6 APS, 2019, pp. 060301
  * [12] Benjamin Nachman, Miroslav Urbanek, Wibe A Jong and Christian W Bauer  “Unfolding quantum computer readout noise”  In _npj Quantum Information_ 6.1 Nature Publishing Group, 2020, pp. 1–7
  * [13] Manuel Erhard, Robert Fickler, Mario Krenn and Anton Zeilinger  “Twisted photons: new quantum perspectives in high dimensions”  In _Light: Science & Applications_ 7.3 Nature Publishing Group, 2018, pp. 17146–17146
  * [14] Daniele Cozzolino, Beatrice Da Lio, Davide Bacco and Leif Katsuo Oxenløwe  “High-Dimensional Quantum Communication: Benefits, Progress, and Future Challenges”  In _Advanced Quantum Technologies_ 2.12 Wiley Online Library, 2019, pp. 1900038
  * [15] Yuchen Wang, Zixuan Hu, Barry C Sanders and Sabre Kais  “Qudits and high-dimensional quantum computing”  In _arXiv preprint arXiv:2008.00959_ , 2020
  * [16] Nurul T Islam et al.  “Provably secure and high-rate quantum key distribution with time-bin qudits”  In _Science advances_ 3.11 American Association for the Advancement of Science, 2017, pp. e1701491
  * [17] Poolad Imany et al.  “50-GHz-spaced comb of high-dimensional frequency-bin entangled photons from an on-chip silicon nitride microresonator”  In _Optics express_ 26.2 Optical Society of America, 2018, pp. 1825–1840
  * [18] Hsuan-Hao Lu et al.  “Quantum Phase Estimation with Time-Frequency Qudits in a Single Photon”  In _Advanced Quantum Technologies_ 3.2 Wiley Online Library, 2020, pp. 1900074
  * [19] Poolad Imany et al.  “High-dimensional optical quantum logic in large operational spaces”  In _npj Quantum Information_ 5.1 Nature Publishing Group, 2019, pp. 1–10
  * [20] Joseph M Lukens  “Quantum information processing with frequency-bin qubits: progress, status, and challenges”  In _CLEO: Applications and Technology_ , 2019, pp. JTu4A–3  Optical Society of America
  * [21] Pei Jiang Low et al.  “Practical trapped-ion protocols for universal qudit-based quantum computing”  In _Physical Review Research_ 2.3 APS, 2020, pp. 033128
  * [22] Brendan Bramman  “Measuring Trapped Ion Qudits”, 2019
  * [23] EO Kiktenko, AK Fedorov, AA Strakhov and VI Man’Ko  “Single qudit realization of the Deutsch algorithm using superconducting many-level quantum circuits”  In _Physics Letters A_ 379.22-23 Elsevier, 2015, pp. 1409–1413
  * [24] Xinsheng Tan et al.  “Experimental observation of tensor monopoles with a superconducting qudit”  In _Physical Review Letters_ 126.1 APS, 2021, pp. 017702
  * [25] Rahul Sawant et al.  “Ultracold polar molecules as qudits”  In _New Journal of Physics_ 22.1 IOP Publishing, 2020, pp. 013027
  * [26] Wolfgang Wernsdorfer and Mario Ruben  “Synthetic Hilbert space engineering of molecular qudits: Isotopologue chemistry”  In _Advanced Materials_ 31.26 Wiley Online Library, 2019, pp. 1806687
  * [27] Eufemio Moreno-Pineda et al.  “Molecular spin qudits for quantum algorithms”  In _Chemical Society Reviews_ 47.2 Royal Society of Chemistry, 2018, pp. 501–513
  * [28] Gian Giacomo Guerreschi, Justin Hogaboam, Fabio Baruffa and Nicolas PD Sawaya  “Intel Quantum Simulator: A cloud-ready high-performance simulator of quantum circuits”  In _Quantum Science and Technology_ 5.3 IOP Publishing, 2020, pp. 034007
  * [29] Shaeema Zaman Ahmed et al.  “Quantum composer: A programmable quantum visualization and simulation tool for education and research”  In _American Journal of Physics_ 89.3 American Association of Physics Teachers, 2021, pp. 307–316
  * [30] Nathan Killoran et al.  “Strawberry fields: A software platform for photonic quantum computing”  In _Quantum_ 3 Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften, 2019, pp. 129
  * [31] Eyob A Sete, William J Zeng and Chad T Rigetti  “A functional architecture for scalable quantum computing”  In _2016 IEEE International Conference on Rebooting Computing (ICRC)_ , 2016, pp. 1–6  IEEE
  * [32] B Butscher and H Weimer  “The libquantum library”
  * [33] “Quirk Quantum Simulator” Accessed: 2021-04-27, https://algassert.com/quirk
  * [34] Reinhold A Bertlmann and Philipp Krammer  “Bloch vectors for qudits”  In _Journal of Physics A: Mathematical and Theoretical_ 41.23 IOP Publishing, 2008, pp. 235303
  * [35] Bernhard Baumgartner, Beatrix Hiesmayr and Heide Narnhofer  “A special simplex in the state space for entangled qudits”  In _Journal of Physics A: Mathematical and Theoretical_ 40.28 IOP Publishing, 2007, pp. 7919
  * [36] Daan Camps, Roel Van Beeumen and Chao Yang  “Quantum Fourier Transform Revisited”  In _arXiv preprint arXiv:2003.03011_ , 2020
  * [37] Stefano Paesani et al.  “Experimental Bayesian quantum phase estimation on a silicon photonic chip”  In _Physical review letters_ 118.10 APS, 2017, pp. 100503


# Symmetry-restoring crossover from defect-free to defect-laden turbulence in
polar active matter

Benjamin H. Andersen  Niels Bohr Institute, University of Copenhagen,
Blegdamsvej 17, Copenhagen, Denmark  Julian Renaud  École Normale Supérieure,
PSL Research University, 45 rue d’Ulm, 75005 Paris, France  Jonas Rønning
Njord Centre, Department of Physics, University of Oslo, P. O. Box 1048, 0316
Oslo, Norway  Luiza Angheluta<EMAIL_ADDRESS>Njord Centre,
Department of Physics, University of Oslo, P. O. Box 1048, 0316 Oslo, Norway
Amin Doostmohammadi<EMAIL_ADDRESS>Niels Bohr Institute, University
of Copenhagen, Blegdamsvej 17, Copenhagen, Denmark

###### Abstract

Coherent flows of self-propelled particles are characterized by vortices and
jets that sustain chaotic flows, referred to as active turbulence. Here, we
reveal a crossover between defect-free active turbulence and active turbulence
laden with topological defects. Interestingly, we show that concurrent to the
crossover from defect-free to defect-laden active turbulence is the
restoration of the previously broken $\operatorname{\mathrm{SO}}(2)$-symmetry
signaled by the fast decay of the two-point correlations. By stability
analyses of the topological charge density field, we provide theoretical
insights on the criterion for the crossover to the defect-laden active
turbulent state. Despite the distinct symmetry features between these two
active turbulence regimes, the flow fluctuations exhibit universal statistical
scaling behaviors at large scales, while the spectrum of polarity fluctuations
decays exponentially at small length scales compared to the active energy
injection length. These findings reveal a new dynamical crossover between
distinct spatiotemporal organization patterns in polar active mater.

##  I Introduction

Recent developments in active matter have shed light on a form of chaotic
flows sustained by living microswimmers at low Reynolds numbers, often known
as active turbulence [1]. A distinctive feature of active turbulence, compared
to classic inertial turbulence, is the continuous injection of energy at the
scale of individual constituents of the active material. Striking examples are
dense bacterial suspensions [2, 3], cellular monolayers [4, 5], and assemblies
of subcellular filaments [6, 7], that are all composed of active elements -
i.e., individual bacterium, single cell, or motor proteins walking on single
filaments - each capable of converting chemical energy to mechanical work [8,
9, 10]. Although, most of the work so far has been limited to visual
resemblance with inertial turbulence, recent works suggest the existence of
universal scaling features in active turbulence [11, 7] and even show [12]
that crossover to active turbulence belongs to the same directed percolation
universality class as in the inertial turbulence [13].

Nevertheless, a majority of studies of active turbulence so far, including the
works suggesting universal scaling laws and universality classes of turbulence
transition, have focused on a subclass of active materials known as active
nematics [14, 15, 16, 17, 18, 19], which models flows generated by dense
assemblies of elongated particles [20]. Within the active nematic framework,
the particles are essentially characterized as shakers [21]: they neither have
any polarity or ability to propel themselves, and instead generate head-tail
symmetric (nematic) active stresses in the fluid. The active stresses in turn
drive hydrodynamic instabilities in the flow and the orientation field of
particles, resulting in a chaotic flow state that is interleaved with
topological defects - singular points in the orientation field where the order
breaks down [8, 10, 22].

Previous studies of polar active matter have characterized spontaneous flows
[23, 24], and associated nonequilibrium steady-states [25, 26]. On the other
hand, studies of self-propulsion effects on active turbulence have mostly
neglected orientational couplings. Instead, the focus was put on generalized
Navier-Stokes equations, where activity is introduced through addition of
phenomenological higher-order terms in the momentum equation, corresponding to
a second order negative viscosity and a fourth-order hyper-viscosity, to give
characteristic vortex length to active turbulence [2, 27, 28, 29]. As such,
characteristics of active turbulence in models of polar active fluids in the
presence of topological defects were only marginally explored in the context
of polar active emulsions [30] and polar flocks with inertia [31].
Additionally, notwithstanding the interesting recent characterization of the
annihilation of topological charges in colloidal flocks [32] and despite the
emerging roles of topological defects in various biological fluids [33, 34,
35, 36], the majority of research so far has focused on nematic topological
defects [10, 35, 37] and studies of active turbulence in polar fluids in the
presence of polar topological defects are lacking.

Figure 1: Crossover from defect-free to defect-laden active turbulence upon
increasing strength of active stresses. Shown here are the averaged defect
density (blue dots) and the rms-velocity (orange dots) in units of
$V_{p}=\sqrt{AK}/\eta$, characterizing the passive relaxation velocity of the
polar particles. Active turbulence arises as soon as the activity is non-zero,
but topological defects form and persist in the polarity field only beyond a
finite threshold both in contractile ($\zeta<0$) and extensile systems
($\zeta>0$). Results were averaged over time and for 10 simulations from
different initial conditions, the error bars standing for the corresponding
standard deviation. The excess kurtosis $\kappa_{i}$ of both quantities are
also represented at the bottom: for large activity levels,
$\kappa(\rho_{d})\simeq 0$, which indicates the population of defects is
typically equilibrium-like and follows a normal distribution law. For small
activity levels, there are no defects at all and $\kappa$ is set by convention
to $-3$. In between, the crossover is characterized by very rare formation and
annihilation events, as evidenced by the peaks in $\kappa$ (blue highlighted
areas). There are no such peaks for the kurtosis of the rms-velocity, which
indicates there is no equivalent crossover regime. Insets illustrate polarity
field (left panels) and the corresponding velocity field (right panels)
characterized by streamlines (solid black lines) colored by the normalized
vorticity $\omega/\omega_{\text{max}}$. Only a small subset of the full domain
($1/16th$ zoomed-in region) is shown. Within the polarity field positive and
negative integer topological defects are marked by red asters and blue
squares, respectively.

Here we use numerical simulations of continuum polar active matter together
with analytical arguments to shed light on the turbulence characteristics of
active fluids laden with topological defects, accounting for both hydrodynamic
effects and polar ordering. We first show how increasing active stresses
results in a regime shift from defect-free active turbulence to active
topological turbulence, where chaotic flows are interleaved with chaotic
arrangements of full-integer topological defects. We provide theoretical
arguments to predict the activity threshold for the crossover to active
topological turbulence, based on the stability analyses of the topological
charge density field. We further show that ordering and symmetry
characteristics are different within these two dynamical regimes: while
$\operatorname{\mathrm{SO}}(2)$-symmetry is broken within the defect-free
active turbulence, the defect-laden turbulence restores the global
$\operatorname{\mathrm{SO}}(2)$-symmetry. Additionally, we demonstrate the
existence of universal scaling behavior in the power spectra of enstrophy and
polarity. This is observed for all activities and within both defect-free and
defect-laden active turbulence.

##  II Methods

###  II.1 Hydrodynamic model

We consider an incompressible polar active fluid in two spatial dimensions,
with the local orientational order described by a two-component order
parameter $\mathbf{p}$ pointing in the direction of self-propulsion. Let then
$\mathbf{u}$ and $\rho$ be the velocity and density of the polar fluid,
respectively. Coarse-grained hydrodynamic equations can be derived by
phenomenological considerations [8, 31, 38] and take the form

$\displaystyle\frac{\mathrm{D}u_{i}}{\mathrm{D}t}=\frac{1}{\rho}\partial_{j}\sigma_{ij},$
(1a) $\displaystyle\frac{\mathrm{D}p_{i}}{\mathrm{D}t}-\left(\lambda
E_{ij}+\Omega_{ij}\right)p_{j}=\frac{1}{\gamma}h_{i},$ (1b)

along with the incompressibility condition $\partial_{i}u_{i}=0$. Here
$\frac{\mathrm{D}}{\mathrm{D}t}=\partial_{t}+u_{j}\partial_{j}$ is the usual
advective derivative. The strain rate tensor
$E_{ij}=(\partial_{i}u_{j}+\partial_{j}u_{i})/2$ and the vorticity tensor
$\Omega_{ij}=(\partial_{i}u_{j}-\partial_{j}u_{i})/2$ are respectively the
symmetric and anti-symmetric parts of the velocity gradient tensor. In the
momentum balance equation (1a), the stress tensor $\sigma_{ij}$ is broken into
a sum of three parts; viscous $\sigma_{ij}^{\text{vis.}}=2\eta E_{ij}$,
passive $\sigma_{ij}^{\text{pas.}}=-P\delta_{ij}+C_{ijkl}p_{k}h_{l}$, and
active stresses
$\sigma_{ij}^{\text{act.}}=-\zeta\left(p_{i}p_{j}-\frac{p^{2}}{2}\delta_{ij}\right)$
[39, 40]. The first term in the passive stress is the usual hydrodynamic
pressure. The second term accounts for elastic stresses through the
anisotropic tensor
$C_{ijkl}=\frac{\lambda+1}{2}\delta_{ik}\delta_{jl}-\frac{\lambda}{2}\delta_{ij}\delta_{kl}+\frac{\lambda-1}{2}\delta_{il}\delta_{jk}$,
with $\lambda$ the flow alignment parameter and
$\mathbf{h}=-\delta\mathcal{F}/\delta\mathbf{p}$ [41, 42] the molecular field
defined here from the free energy

$\mathcal{F}=\int\mathrm{d}^{2}\mathbf{x}\left\\{A\left(-\frac{|\vec{p}|^{2}}{2}+\frac{|\vec{p}|^{4}}{4}\right)+\frac{K}{2}\partial_{i}p_{j}\partial_{i}p_{j}\right\\}.$
(2)

The free energy contains a local energy density with an energy scale $A$ that
controls the isotropic-polar transition favoring the emergence of finite
polarity at $|\mathbf{p}|=1$, and a non-local energy contribution with an
elastic constant $K$ that penalizes deformations in the polarity field [43].
It is important to note that we do not treat the polarity as a unit vector,
with fixed magnitude. Instead, in this formulation the polarity modulus is an
important dynamical variable.

###  II.2 Numerical method and the simulation parameters

We simulate equations (1a,1b) using a hybrid lattice-Boltzmann method,
combining finite-difference method for the evolution of polarity vector Eq.
(1b), and the lattice-Boltzmann method for solving the incompressible Navier-
Stokes equation Eq. (1a) with $\rho=40$ and $\eta=3.6$ in lattice Boltzmann
units, ensuring that the Reynolds number in the simulations is negligible
($\operatorname{\mathrm{Re}}\ll 1$) [44, 12]. The other relevant dimensionless
numbers describing the system are: (i) the dimensionless ratio of the
viscosities $\eta/\gamma$, (ii) ratio of micro to macro length scales
$(\sqrt{K/A})/$L, (iii) the flow alignment parameter $\lambda$, and (iv) the
dimensionless active stress $\bar{\zeta}=\zeta/A$. Unless otherwise stated, we
fix the viscosity ratio to $\eta/\gamma=3.6$, micro to macro length scale to
$(\sqrt{K/A})/L=2\times 10^{-4}$ (assuring that the coherence length
$l_{p}=\sqrt{K/A}$ is significantly smaller than the domain size $L$), and the
flow alignment parameter to $\lambda=0.1$

Simulations were initialized with quiescent velocity field and noisy polar
alignments close to the uniformly oriented state
$\mathbf{p}=\hat{\mathbf{e}}_{x}$ under periodic boundary conditions, on
quadratic domains of linear dimension $L=1024$, unless otherwise is stated.

##  III Results

###  III.1 Activity-induced crossover to defect-laden turbulence

We begin by introducing global measures of the flow fluctuations and polar
order parameter as functions of the activity parameter $\bar{\zeta}$. The flow
is globally characterized by the dimensionless root mean squared (rms)
velocity $\overline{V}_{\text{rms}}=V_{\text{rms}}/V_{p}$ that is normalized
by the characteristic velocity of passive relaxation of polarity
$V_{p}=\sqrt{AK}/\eta$. The polar order parameter is associated with the SO(2)
symmetry and carries full-integer topological defects. Its global measure is
given by the average density of topological defects $\rho_{d}=\langle
N_{d}\rangle/L^{2}$. The averaging is done both over space and time in the
statistical steady-state regime. As evident from Fig. 1, increasing activity
beyond a certain threshold results in a continuous increase in the defect
density for both extensile ($\zeta>0$) and contractile ($\zeta<0$) systems. We
have carefully checked that the activity thresholds do not depend on the
system size $L$, by simulating different domain sizes of $512\times 512$,
$1024\times 1024$, and $4096\times 4096$ and finding precisely the same
quantitative dependence of defect density on the activity (see Fig. 2; square
and circle symbols). Interestingly, simultaneous measurement of the rms
velocity $\overline{V}_{\text{rms}}$ indicates that even below the critical
activity for defect nucleation, active stresses disturb the polarity field
such that spontaneous flows are generated within the system, for any non-zero
activity (Fig. 1C). A closer look at the velocity field of the system below
the defect nucleation threshold shows that chaotic flows characterized by flow
vortices and jets span the system (Fig. 1A; only a fraction of the entire
domain is shown), although no topological defects are present in the polarity
field. Qualitatively similar chaotic flows manifest at higher activities with
a smaller typical length scale $L_{\text{act.}}$ defined below, though with
the distinctive feature that the flows are now laden with the presence of
full-integer topological defects (unbound vortex and antivortex pairs) in the
polarity field (Fig. 1B). Together, these results establish that there is a
well-defined activity threshold for the crossover from defect-free and defect-
laden active turbulence in polar active matter.

Figure 2: Global symmetry restoration. (A) Averaged defect density $\rho_{d}$,
blue solid line, and averaged magnitude of the global polarization
$\mathbf{\Psi}$ (Eq. 3), orange solid line. (B) The correlation length
$\bar{\xi}$ (in units of $l_{p}=\sqrt{K/A}$) defined by fitting an exponential
to the pair-correlation function $\Gamma(r)$. The square and circle symbols in
(B) correspond to the simulations results for $4096\times 4096$ and
$1024\times 1024$ domains, respectively, showing no detectable difference.

###  III.2 Global symmetry is restored in the defect-laden active turbulence

We next investigate whether the emergence of topological defects leads to
alterations in the global ordering features of the polar active turbulence. To
this end, we first introduce the global polarization $\mathbf{\Psi}$ as the
spatiotemporal average value of the polarity field $\mathbf{p}$

$\mathbf{\Psi}=\int\frac{\mathrm{d}t}{T}\int\frac{\mathrm{d}^{2}\mathbf{x}}{L^{2}}\,\mathbf{p}(t,\mathbf{x}).$
(3)

The magnitude of this global polarization $|\mathbf{\Psi}|$ is a macroscopic
order parameter for polar order akin to total magnetization in spin systems:
$|\mathbf{\Psi}|=0$ is associated with disordered states (i.e. continuous
rotational symmetry) while $|\mathbf{\Psi}|>0$ is associated with ordered
states (discrete/broken rotational symmetry). We use this global measure to
probe the crossover between the two active turbulence regime as shown in Fig.
2A; orange solid line. Upon increasing the active stress from zero, the active
turbulence state immediately manifests a broken
$\operatorname{\mathrm{SO}}(2)$-symmetry evident from a non-zero value of
$|\mathbf{\Psi}|$ in the absence of topological defects. The global
polarisation decreases in magnitude with increasing activity and, remarkably,
drops to zero beyond a critical activity for the nucleation of topological
defects. This is indicative of the restoration of the
$\operatorname{\mathrm{SO}}(2)$-symmetry in the active defect-laden
turbulence.

To further characterize the signature of this crossover, we also measure the
equal-time, spatial pair-correlation function for the polarity field
$\Gamma(\mathbf{r})=\langle\delta p_{i}(\mathbf{x})\delta
p_{i}(\mathbf{x}+\mathbf{r})\rangle$, where $\delta
p_{i}(\mathbf{r})=p_{i}(\mathbf{r})-\langle p_{i}(\mathbf{r})\rangle$ is the
deviation from the local mean value and the averaging is performed over $10$
different realizations, with distinct initial conditions and over time, once
the system has reached a statistical steady state. Accordingly, the
correlation length $\xi$ is obtained by fitting an exponential to the the
pair-correlation function $\Gamma(r)$. Similar to the global polarization, the
correlation length $\bar{\xi}=\xi/l_{p}$, normalized with respect to the
coherence length $l_{p}=\sqrt{K/A}$, decreases with increasing activity,
showing a fast decay that is followed by a significantly slower decay of the
correlation length at higher activities (Fig. 2). The regime shift is marked
by the crossover from the fast to slow decay of the correlation length
$\bar{\xi}$ which matched with the onset of defect nucleation and vanishing of
$|\mathbf{\Psi}|=0$ as evidenced in Fig. (2A, B). We have carefully checked
that correlation measurements are independent of the domain sizes for very
large systems up to $4096\times 4096$. We note that the behavior of the
correlation length shares interesting similarity to that in frustrated two-
dimensional Heisenberg magnets, where a crossover between a topological
defects dominated regime at high temperatures and a spin-wave regime at low
temperatures is observed [45]. Together, these measurements of the global
polarization and the correlation length indicate that the crossover from
defect-free to defect-laden active turbulence is marked by significant
alterations in both the global symmetry of the collective organization and in
the local correlations between the polarities of active constituents.

It is important to note that the observed restoration of the
$\operatorname{\mathrm{SO}}(2)$-symmetry is solely activity-driven and is in
the noise-less limit of the equations of motion. Furthermore, the result
presented in Fig. 2 demonstrate that, point-by-point, for any activity, (i)
the correlation length, (ii) polarisation, and (iii) defect density for domain
sizes $1024\times 1024$ and $4096\times 4096$ fall exactly on top of each
other. This strongly suggests a system-size independent behavior, at least for
the numerically accessible system sizes. The question about the existence of a
well-defined thermodynamic limit of this symmetry-breaking/restoring
transition remains open for future studies. Additionally, it is constructive
to note that we do not observe any difference between longitudinal and
transverse correlations in the polarity. The correlation length does not
diverge with the system size within the ordered state, which we conjecture is
due to the non-equilibrium nature of the system and the breaking of detailed
balance [46] through activity-induced flows that couple to the dynamics of the
order parameter.

Figure 3: Defect pair nucleation process. (A,B) Temporal evolution of the
average bulk free-energy
$\bar{F}_{\mathrm{bulk}}=(F_{\mathrm{bulk}}+A/4)/AL_{\mathrm{act.}}$ (top)
measured in units of $AL_{\mathrm{act.}}$ relative to the background field
$-A/4$; along with (C,D) simulation snapshots depicting the polarity field and
diffuse charge density (heat plot) before, at, and after the nucleation of one
pair of oppositely charged topological defects (bottom). Column (A,C)
contractile and (B,D) extensile. The blue regions in A, B mark the period of
existence of the defect pair considered. Time $\bar{t}$ is measured in units
of the passive relaxation time of polarity over the simulation domain
$t_{p}=\frac{\gamma L}{\sqrt{AK}}$ (see also Supplementary Movies [47]).

###  III.3 Mechanism of defect pair nucleation in active turbulence of polar
particles

Having established the significant impact of the topological defects on the
collective organization of polar active matter and its spontaneous flows, we
next explore the mechanism of how topological defects are formed within the
active turbulence state as the activity approaches the threshold value. Figure
3 illustrates zoomed-in snapshots of the polarity field at the onset of one
defect pair nucleation. Three snapshots are shown at simulation times
corresponding to: before, at, and after a single pair nucleation for both
extensile and contractile active stresses. To identify and track topological
defects, we make use of the diffuse charge density, defined in the next
section as
$D=\frac{1}{2\pi}\epsilon_{ij}\epsilon_{kl}\partial_{i}p_{k}\partial_{j}p_{l}$:
it carries the $\pm 1$ charge of the integer topological defect at the
singularities in the polarity field, and is zero elsewhere. Shortly before the
onset of defect pair nucleation, the polarity field is characterized by
locally ordered domains separated by lines of kink walls (Fig. 3 and
Supplementary Movies [47]). At the onset of defect pair nucleation the
polarity within one of the kink walls flips locally, leading to spontaneous
emergence of a pair of full-integer defects (Fig. 3). After the nucleation
event, the defect pair gets separated by an ordered region of polarity
alignment (Fig. 3). While the negatively charged antivortices ($-1$
topological defects) have similar structures in both extensile and contractile
systems, the positively charged $+1$ topological defects take the form of
asters and vortices in contractile and extensile systems, respectively (Fig.
3C,D), in agreement with earlier analytical predictions on the stability of
defect structures in polar active matter [48].

In order to gain more insight into the defect pair nucleation process, we
measure the averaged bulk free energy of the system from Eq. (2) during the
time leading to a single defect pair nucleation event, as depicted in (Fig.
3A-B). For both extensile and contractile systems the onset of defect
nucleation coincides with a burst in the bulk free energy. This is consistent
with the spontaneous flipping of the polarity within the kink walls that
occurs when activity-induced fluctuations in the polarity overcome the energy
barrier of local flipping, set by the bulk free energy. Once the defect pair
is formed, the bulk free energy drops as the high energy stored in the kink
walls is released. In the following section, we look more closely at this
defect nucleation process and derive an analytical description of the activity
threshold that is set by the competition between bulk and elastic energies of
the system, and activity-induced flows.

###  III.4 Onset of active defect-laden turbulence

To gain more theoretical insights into the onset of defect pair nucleation, we
perform a stability analysis of the topological charge density field. A
stability analysis of the polarity field to probe the nucleation of
topological defects is unfeasible since it is intrinsically nonlinear.
However, we take advantage of the property of the topological charge density
field being zero for low-energy excitations of the polarity field (aka kinks,
domain walls) to study the onset of defect nucleation as an instability to a
nonzero and localised charge density field. To this end, we employ the
Halperin-Mazenko method [49, 50, 51] of topological defects as zeros in the
polarity vector field $\mathbf{p}$, i.e. localised regions where $\mathbf{p}$
vanishes in magnitude while its orientation is undefined (multi-valued). The
associated topological charge is a quanta of the phase jump measured around an
arbitrary contour $C_{\alpha}$ enclosing the defect

$2\pi
q_{\alpha}=\oint\limits_{C_{\alpha}}\mathrm{d}\mathbf{l}\cdot\nabla\theta,$
(4)

where $\theta=\arg(\mathbf{p})$. By Stokes’ theorem, we see that topological
defects correspond to phase singularities where the phase ceases to be
irrotational,

$2\pi
q_{\alpha}\int\limits_{A_{\alpha}}\mathrm{d}^{2}\mathbf{x}\,\delta(\mathbf{x}-\mathbf{x}_{\alpha})=\int\limits_{A_{\alpha}}\mathrm{d}^{2}\mathbf{x}\,\nabla\times\nabla\theta.$
(5)

As shown in Ref. [51], the topological content inside an area region
$A_{\alpha}$ can be also obtained from the corresponding integration in the
order parameter space $(p_{1},p_{2})$

$2\pi
q_{\alpha}\int\limits_{A_{\alpha}}\mathrm{d}^{2}\mathbf{x}\,\delta(\mathbf{x}-\mathbf{x}_{\alpha})=\int\mathrm{d}^{2}\mathbf{p}\mathrm{det}(\nabla\mathbf{p})\delta(\mathbf{p}).$
(6)

We relate the topological defect density $D$ with the Jacobi determinant of
the transformation to the order parameter space, which is the same as the
determinant of the distortion tensor $\nabla\mathbf{p}$ to define

$D=\frac{1}{2}\epsilon_{ij}\epsilon_{kl}\partial_{i}p_{k}\partial_{j}p_{l},$
(7)

which is a smooth scalar field that vanishes where $\mathbf{p}$ is smooth and
becomes non-zero and localized at the core topological defect picking up the
phase singularity. We have validated the defects detected by D-filed against
our routine approach of using winding angles and have obtained identical
results, in agreement with previous studies on active nematics, where the
$D$-field was defined as the diffuse charge density [52, 53, 54]. The
$D$-field has the physical interpretation of a non-singular topological charge
density field which follows the conservation law [51]

$\partial_{t}D+\partial_{i}J_{i}=0,$ (8)

with the corresponding topological current density

$J_{i}=-\epsilon_{ij}\epsilon_{kl}\partial_{t}p_{k}\partial_{j}p_{l},$ (9)

which is fully determined by the evolution of the polarity $\mathbf{p}$ and
the flow field $\mathbf{u}$ through the main Eqs. (1b-1a).

We perform a stability analysis of Eq. (8) around the ground state of zero $D$
associated with slowly-varying polar order to estimate the critical activity
for the onset of defect nucleation, identified as an instability where $D$
becomes non-zero because polar order vanishes locally. In general, the
evolution of $D$, resulting from inserting the evolution of $\mathbf{p}$ from
Eq. (1b) into Eq. (8), cannot be closed only in terms of the $D$-field.
However, in the limit of retaining only the shear flow alignment contribution
and the molecular term, the evolution of $\mathbf{p}$ from Eq. (8) can be
reduced to

$\partial_{t}D=\epsilon_{ij}\epsilon_{kl}\left[\lambda\partial_{i}(E_{kn}p_{n})\partial_{j}p_{l}+\frac{1}{\gamma}\partial_{i}h_{k}\partial_{j}p_{l}\right].$
(10)

Now, by balancing viscous stress with the active stress in Eq. (1a), we relate
the traceless strain rate $E_{kn}$ directly to the nematic order parameter,
namely

$E_{kn}=\frac{\zeta}{2\eta}\left(p_{n}p_{k}-\frac{p^{2}}{2}\delta_{nk}\right),$
(11)

such that the flow alignment contribution dominated by orientational
distortion has a closed form given by

$\displaystyle\frac{\zeta\lambda}{4\eta}\epsilon_{ij}\epsilon_{kl}\partial_{i}(p^{2}p_{k})\partial_{j}p_{l}=$
$\displaystyle\frac{\zeta\lambda}{4\eta}\left(p^{2}\epsilon_{ij}\epsilon_{kl}\partial_{i}p_{k}\partial_{j}p_{l}+2\epsilon_{ij}\epsilon_{kl}p_{k}p_{m}\partial_{i}p_{m}\partial_{j}p_{l}\right)=$
$\displaystyle\frac{\zeta\lambda}{2\eta}p^{2}D+\frac{\zeta\lambda}{2\eta}p^{2}(\partial_{1}p_{1}\partial_{2}p_{2}-\partial_{2}p_{1}\partial_{1}p_{2})=\frac{\zeta\lambda}{\eta}p^{2}D.$
(12)

The molecular field is due to the elastic distortions and the ordering
potential

$h_{k}=K\nabla^{2}p_{k}+A(1-p^{2})p_{k},$

and leads to the contribution

$\epsilon_{ij}\epsilon_{kl}\partial_{i}h_{k}\partial_{j}p_{l}=K\nabla^{2}D+2A(1-2p^{2})D.$

Thus, within this limiting case, the evolution of the $D$-field reduces to a
diffusion-reaction equation

$\partial_{t}D\approx\frac{K}{\gamma}\nabla^{2}D+\left(\frac{2A}{\gamma}(1-2p^{2})+\frac{\zeta\lambda}{\eta}p^{2}\right)D.$
(13)

First, we notice that the polarity stiffness $K$ enters as en effective
diffusivity coefficient and sets the scale of the diffusive core. The kinetic
rate coefficients arise from the stabilizing passive contribution through the
energy scale parameter $A$ of the local order potential and its relaxation
timescale $\gamma$, as well as from the destabilizing contribution through the
flow alignment $\lambda$, activity $\zeta$ and viscosity $\eta$. At the first
sight, Eq. (13) looks non-conservative due to the source/sink term. This is,
however, not the case since this term is derived from a conservative current
and is related to the fact that the $D$-field itself is the divergence of a
vector field,
$D=\partial_{i}(\epsilon_{ij}\epsilon_{kl}p_{k}\partial_{j}p_{l})/2$. We
estimate the critical activity from the linear stability analysis of Eq. (13)
about the ordered state $p=1+\delta p$ and $D=\delta D$. By balancing out the
two contribution to the reaction term we find that the equation becomes
unstable at

$\zeta_{c}\approx\frac{2\eta}{\gamma\lambda}A.$ (14)

Thus, in the limit where we neglect the convective and rotational terms to the
evolution of the $D$-field (i.e. $\mathbf{u}\cdot\nabla\mathbf{p}$ and
$[\mathbf{\Omega},\mathbf{p}]$), the topological charge density becomes
unstable to defect nucleation only for extensile systems where $\zeta>0$. This
instability is due to the flow alignment.

Advection and rotation of polarity $\mathbf{p}$, due to flow velocity and
vorticity, could also trigger instabilities in the $D$-field. Alas, their
contributions are non-locally through viscosity and incompressibility
constraint and cannot be reduced to linear operators acting on $D$. However,
it turns out that the contribution of vorticity to the $D$-evolution is
symmetric with respect to the sign of $\zeta$, and, from this, we can infer
that an instability induced by vorticity would occur at a critical activity
both for extensile/contractile systems, also symmetric with respect to the
sign of $\zeta$.

It is important to note that in numerical simulations all these destabilizing
forces orchestrate and compete with the relaxation to uniform order.
Therefore, the value of the critical activity is different than the estimated
one corresponding to an isolated triggering factor. However, this stability
analysis is informative and predicts generic properties of $\zeta_{c}$: i) it
depend on the energy scale $A$ controlling polar order, and ii) that the
threshold activity has an asymmetry with respect to extensile/contractile
active stresses in the presence of flow alignment.

To test the theoretical predictions on the activity threshold for defect
nucleation, we next explore the dependence of defect density on the flow-
aligning parameter $\lambda$ and the energy scale $A$ for the bulk free
energy. Figure 4A compares the variation of defect density $\rho_{d}$ with
activity for three values of the flow-aligning parameter $\lambda$.
Interestingly, and consistently with the theoretical prediction, for
$\lambda=0$ the activity threshold for defect nucleation is symmetric with
respect to the sign of the activity. Additionally, as predicted from the
theory, a positive (negative) flow-aligning parameter shifts the defect
density curve such that the activity threshold becomes smaller - in absolute
value - for extensile (contractile) activities compared to their contractile
(extensile) counterparts, leading to a more asymmetric profile of the defect
density with respect to the sign of the activity. Moreover, as predicted from
the theory, decreasing the energy scale in the bulk free energy $A$ lowers the
activity threshold for defect nucleation (Fig. 4B). This is also consistent
with the mechanism identified for the nucleation of pairs of defects, relying
on the spontaneous flipping of the polarity within the kink walls (Fig. 3C,D
and Supplementary Movies [47]): since the parameter $A$ sets the depth of a
double-well potential in the bulk free energy, decreasing its value amounts to
lowering the energy barrier for flipping the polarity, and as such leads to a
reduced activity threshold for defect pair nucleation.

Figure 4: Dependence of the steady-state density of defects on (A) the flow-
aligning parameters $\lambda$ and (B) the bulk free-energy strength A. Inset
in (B) shows that normalizing the active stress by $\sqrt{A}$ results in a
data collapse near the crossover point.

###  III.5 Universal scaling in defect-free and defect-laden active turbulence

Figure 5: Universal scaling in defect-free and defect-laden active turbulence.
Power spectra normalized with its corresponding maximum value for (A) polarity
$\hat{E}_{p}(k)$ and (B) kinetic energy $\hat{E}_{v}(k)$ on semilog plots,
demarcating the exponential decay at small scales (large wave numbers) that
imply the dynamics is dominated by viscous dissipation. (C,D) Power spectra
normalized with its corresponding maximum value for (C) the enstrophy
$\hat{E}_{w}(k)$ and (D) kinetic energy again, on log-log plots, showing the
power-law scaling at scales larger than the active length scale
$L_{\mathrm{act.}}$. Insets show the non-scaled power spectra and the colorbar
indicates the dimensionless active stress strength $\bar{\zeta}$.

Having established the mechanism of defect nucleation and the difference
between defect-free and defect-laden active turbulence, we now investigate the
common features between these two active turbulent states by characterizing
the energetic features of the flow and order parameter (Fig. 5). To this end,
we first measure the power spectrum of the polarity order parameter
$\hat{E}_{p}(k)=\tfrac{1}{2}\langle\hat{p}_{i}(k)\hat{p}_{i}(-k)\rangle$, with
$k$ the wave number averaged over azimuthal direction, for a range of
activities spanning both defect-free and defect-laden turbulent states.
Interestingly, when the wave number is non-dimensionalized by the active
length scale $L_{\text{act.}}=\sqrt{K/|\zeta|}$, which characterizes the scale
of energy injection into the system by active stresses, the power spectrum
$\hat{E}_{p}$ for all activity values shows an exponential decay at wave
numbers larger than the inverse active length scale, as evident from a semilog
plot of the order parameter spectrum versus the normalized wave number (Fig.
5A). Notice that the exponential cutoff scale still depends on activity and
approaches a constant value at sufficiently high activity. This is related to
the initial fast decay of the correlation function $\Gamma(r)$, with a
correlation length that is activity dependent as discussed earlier. Notice
that the exponential dependence of $\hat{E}_{p}$ on $k$ breaks down for small
wavenumbers. On the other hand, robust exponential tails at large wave numbers
are characteristic for the spectral properties of the spontaneous flow,
measured by the kinetic energy spectrum
$\hat{E}_{v}(k)=\tfrac{1}{2}\langle\hat{v}_{i}(k)\hat{v}_{i}(-k)\rangle$, with
$\hat{v}_{i}$ the Fourier transforms of flow velocity (Fig. 5B). However, for
wave numbers below the inverse active length scale, the kinetic energy
spectrum $\hat{E}_{v}(k)$ and enstrophy spectrum
$\hat{E}_{\omega}(k)=\tfrac{1}{2}\langle\hat{\omega}(k)\hat{\omega}(-k)\rangle\sim
k^{2}\hat{E}_{v}(p)$, both exhibit power-law behavior (Fig. 5C,D). This is
important, because these activity levels span both the defect-free and defect-
laden active turbulence states, indicating that for a polar active matter,
despite distinct symmetry features and correlation lengths, there exist
universal scaling of flow and order parameter spectrum in active turbulence
state with and without topological defects. In particular, the enstrophy power
spectrum that relates to the flow field shows a universal power-law scaling at
length scales larger than the active length $L_{\text{act.}}$, while the
polarity power spectrum that relates to the order parameter field shows a
universal exponential decay at length scales smaller than the active length
scale. These two distinct scaling relate defect-free and defect-laden
turbulence in polar active matter. It is important to note that the scaling
observed here is different from the power-law scaling suggested for active
nematic turbulence without topological defects, where the orientational field
is treated as a unit vector field with fixed modulus [11]. The heuristic
argument of the scaling behavior reported in Ref. [11] is based on a spectral
analysis of the equation for the vorticity source to argue for the power-law
scaling of $\hat{E}_{p}(k)$ and $\hat{E}_{v}(k)$ at small $k$. The scaling
argument hinges on two important assumptions: 1) vorticity source is dominated
by the active stress determined only by phase gradients; 2) at small $k$ (long
wavelength) phase fluctuations become uncorrelated. Both of these assumptions
fail to apply when the polarity magnitude is allowed to vary, since the
vorticity fluctuations couple with fluctuations in both magnitude and phase of
the polarity. Thus, the power spectrum of vorticity source is not trivially
related to polarity power spectrum by extending the scaling analysis from Ref.
[11]. Therefore, we do not expect to observe the same scaling behavior at
long-wavelength as in Ref. [11].

##  IV Discussion

The results presented herein reveal distinct flow fields and patterns of
collective self-organization of polar active particles that are primarily
controlled by active stresses. We show that increasing active stress in a
polar active fluid leads to a dynamical crossovers from a defect-free to a
defect-laden active turbulence state. Above a critical activity threshold, the
turbulent flow is seeded with topological defects nucleated in the polar order
parameter. Importantly, we find that the proliferation of defects screens the
global polar order. Thus, in the defect-laden turbulent states the
$\operatorname{\mathrm{SO}}(2)$-symmetry is restored, whereas the defect-free
active turbulence is endowed with broken
$\operatorname{\mathrm{SO}}(2)$-symmetry.

Interestingly, we show that in spite of their distinct topological content,
the two active turbulence regimes share similar statistical properties with
robust exponential tails for power spectra of polar order, flow velocity and
vorticity on length scales smaller than the active length scale. On the other
hand, over the length scales that are larger than the active length scale,
enstrophy and energy spectra of flow fluctuations exhibit self-similarity with
power-law scaling exponents. These statistical properties are also universal
in the sense that they do not depend on the presence of topological defects,
which may indicate that non-topological dissipative structures (e.g. kink
walls) are the dominant driving force to fluctuations in order and flow of
polar active fluids.

It is noteworthy that despite the ubiquitous presence of polarity in living
and synthetic active materials [8, 9], majority of the current understanding
of topological defects and flow features within the active turbulence is based
on studies of active nematics. This is in part due to the growing number of
biological systems that are identified with half-integer, nematic, topological
defects [33, 35, 22] and discovery of potential biological functions for such
defects [55, 56, 57]. It is important to note, however, that full-integer
defects have also been identified as possible organization centers for mitotic
spindles in microtubule-motor protein assemblies [58] and more recently as
potential sites for mechanotransduction [59], and for cell differentiation in
mouse myoblasts [34] and cartilage cells [60]. Our results provide a
quantitative characterisation of the dynamical crossovers, as well as the
statistical imprints of flow and topological defects in polar active matter.
Future work should focus on details of the dynamic distinctions in topological
and flow features between the active turbulence in polar and nematic active
materials, as well as generalized frameworks of the coupling between polar and
nematic symmetries that have already been suggested theoretically [61, 62, 38]
and observed experimentally [58, 60].

## Acknowledgement and Funding

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. J.R. and
L.A. acknowledge support from the Research Council of Norway through the
Center of Excellence funding scheme, Project No. 262644 (PoreLab).

## References

  * [1] Ricard Alert, Jaume Casademunt, and Jean-François Joanny.  Active turbulence.  Annual Review of Condensed Matter Physics, 13(1):143–170, 2022\.
  * [2] Henricus Wensink, Jörn Dunkel, Sebastian Heidenreich, Knut Drescher, Raymond E Goldstein, Hartmut Löwen, and Julia M Yeomans.  Meso-scale turbulence in living fluids.  PNAS, 109(36):14308–14313, 2012.
  * [3] O. J Meacock, A. Doostmohammadi, K. R Foster, J. M Yeomans, and W. M Durham.  Bacteria solve the problem of crowding by moving slowly.  Nature Phys., 17(2):205–210, 2021.
  * [4] C Blanch-Mercader, V Yashunsky, S Garcia, G Duclos, L Giomi, and Pascal Silberzan.  Turbulent dynamics of epithelial cell cultures.  Physical review letters, 120(20):208101, 2018.
  * [5] Shao-Zhen Lin, Wu-Yang Zhang, Dapeng Bi, Bo Li, and Xi-Qiao Feng.  Energetics of mesoscale cell turbulence in two-dimensional monolayers.  Communications Physics, 4(1):1–9, 2021.
  * [6] T. Sanchez, D. T. N. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic.  Spontaneous motion in hierarchically assembled active matter.  Nature, 491:431, 2012.
  * [7] Berta Martínez-Prat, Ricard Alert, Fanlong Meng, Jordi Ignés-Mullol, Jean-François Joanny, Jaume Casademunt, Ramin Golestanian, and Francesc Sagués.  Scaling regimes of active turbulence with external dissipation.  arXiv preprint arXiv:2101.11570, 2021.
  * [8] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, Madan Rao, and R. Aditi Simha.  Hydrodynamics of soft active matter.  Rev. Mod. Phys., 85:1143–1189, 2013.
  * [9] Clemens Bechinger, Roberto Di Leonardo, Hartmut Löwen, Charles Reichhardt, Giorgio Volpe, and Giovanni Volpe.  Active particles in complex and crowded environments.  Rev. Mod. Phys., 88:045006, 2016.
  * [10] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués.  Active nematics.  Nat. Commun., 9(1):3246, 2018.
  * [11] Ricard Alert, Jean-François Joanny, and Jaume Casademunt.  Universal scaling of active nematic turbulence.  Nat. Phys., pages 1–7, 2020.
  * [12] Amin Doostmohammadi, Tyler N Shendruk, Kristian Thijssen, and Julia M Yeomans.  Onset of meso-scale turbulence in active nematics.  Nature communications, 8(1):1–7, 2017.
  * [13] Grégoire Lemoult, Liang Shi, Kerstin Avila, Shreyas V Jalikop, Marc Avila, and Björn Hof.  Directed percolation phase transition to sustained turbulence in couette flow.  Nature Physics, 12(3):254–258, 2016.
  * [14] Sandrine Ngo, Anton Peshkov, Igor S Aranson, Eric Bertin, Francesco Ginelli, and Hugues Chaté.  Large-scale chaos and fluctuations in active nematics.  Physical review letters, 113(3):038302, 2014.
  * [15] Javier Urzay, Amin Doostmohammadi, and Julia M Yeomans.  Multi-scale statistics of turbulence motorized by active matter.  Journal of Fluid Mechanics, 822:762–773, 2017.
  * [16] Žiga Krajnik, Žiga Kos, and Miha Ravnik.  Spectral energy analysis of bulk three-dimensional active nematic turbulence.  Soft Matter, 16(39):9059–9068, 2020.
  * [17] Livio Nicola Carenza, Luca Biferale, and Giuseppe Gonnella.  Cascade or not cascade? energy transfer and elastic effects in active nematics.  Europhysics Letters, 132(4):44003, 2020.
  * [18] Rodrigo CV Coelho, Nuno AM Araújo, and Margarida M Telo da Gama.  Propagation of active nematic–isotropic interfaces on substrates.  Soft Matter, 16(17):4256–4266, 2020.
  * [19] Santhan Chandragiri, Amin Doostmohammadi, Julia M Yeomans, and Sumesh P Thampi.  Flow states and transitions of an active nematic in a three-dimensional channel.  Physical Review Letters, 125(14):148002, 2020.
  * [20] Amin Doostmohammadi, Jordi Ignés-Mullol, Julia M Yeomans, and Francesc Sagués.  Active nematics.  Nat. Commun., 9(1):1–13, 2018.
  * [21] Sriram Ramaswamy.  The mechanics and statistics of active matter.  Annu. Rev. Condens. Matter Phys., 1(1):323–345, 2010.
  * [22] Aleksandra Ardaševa and Amin Doostmohammadi.  Topological defects in biological matter.  Nature Reviews Physics, 4(6):354–356, 2022.
  * [23] R Voituriez, JF Joanny, and J Prost.  Generic phase diagram of active polar films.  Physical review letters, 96(2):028102, 2006.
  * [24] Luca Giomi, M Cristina Marchetti, and Tanniemola B Liverpool.  Complex spontaneous flows and concentration banding in active polar films.  Physical review letters, 101(19):198101, 2008.
  * [25] Elsen Tjhung, Michael E Cates, and Davide Marenduzzo.  Nonequilibrium steady states in polar active fluids.  Soft Matter, 7(16):7453–7464, 2011.
  * [26] Luca Giomi and M Cristina Marchetti.  Polar patterns in active fluids.  Soft Matter, 8(1):129–139, 2012.
  * [27] Vasil Bratanov, Frank Jenko, and Erwin Frey.  New class of turbulence in active fluids.  Proceedings of the National Academy of Sciences, 112(49):15048–15053, 2015.
  * [28] Moritz Linkmann, Guido Boffetta, M Cristina Marchetti, and Bruno Eckhardt.  Phase transition to large scale coherent structures in two-dimensional active matter turbulence.  Physical review letters, 122(21):214503, 2019.
  * [29] Moritz Linkmann, M Cristina Marchetti, Guido Boffetta, and Bruno Eckhardt.  Condensate formation and multiscale dynamics in two-dimensional active suspensions.  Physical Review E, 101(2):022609, 2020.
  * [30] Livio Nicola Carenza, Luca Biferale, and Giuseppe Gonnella.  Multiscale control of active emulsion dynamics.  Physical Review Fluids, 5(1):011302, 2020.
  * [31] Rayan Chatterjee, Navdeep Rana, R Aditi Simha, Prasad Perlekar, and Sriram Ramaswamy.  Inertia drives a flocking phase transition in viscous active fluids.  Physical Review X, 11(3):031063, 2021.
  * [32] Amélie Chardac, Ludwig A Hoffmann, Yoann Poupart, Luca Giomi, and Denis Bartolo.  Topology-driven ordering of flocking matter.  Physical Review X, 11(3):031069, 2021.
  * [33] Amin Doostmohammadi and Benoit Ladoux.  Physics of liquid crystals in cell biology.  Trends in Cell Biology, 2021.
  * [34] Pau Guillamat, Carles Blanch-Mercader, Guillaume Pernollet, Karsten Kruse, and Aurélien Roux.  Integer topological defects organize stresses driving tissue morphogenesis.  Nature materials, 21(5):588–597, 2022.
  * [35] Suraj Shankar, Anton Souslov, Mark J Bowick, M Cristina Marchetti, and Vincenzo Vitelli.  Topological active matter.  Nature Reviews Physics, 4(6):380–398, 2022.
  * [36] Kirsten D Endresen, MinSu Kim, Matthew Pittman, Yun Chen, and Francesca Serra.  Topological defects of integer charge in cell monolayers.  Soft matter, 17(24):5878–5887, 2021.
  * [37] Farzan Vafa, Mark J Bowick, M Cristina Marchetti, and Boris I Shraiman.  Multi-defect dynamics in active nematics.  arXiv:2007.02947, 2020.
  * [38] Aboutaleb Amiri, Romain Mueller, and Amin Doostmohammadi.  Unifying polar and nematic active matter: emergence and co-existence of half-integer and full-integer topological defects.  Journal of Physics A: Mathematical and Theoretical, 55(9):094002, feb 2022.
  * [39] Jacques Prost, Frank Jülicher, and Jean-François Joanny.  Active gel physics.  Nature physics, 11(2):111–117, 2015.
  * [40] Frank Jülicher, Stephan W Grill, and Guillaume Salbreux.  Hydrodynamic theory of active matter.  Rep. Prog. Phys., 81(7):076601, 2018.
  * [41] Pierre-Gilles De Gennes and Jacques Prost.  The physics of liquid crystals, volume 83.  Oxford university press, 1993.
  * [42] Santhan Chandragiri, Amin Doostmohammadi, Julia M Yeomans, and Sumesh P Thampi.  Active transport in a channel: stabilisation by flow or thermodynamics.  Soft matter, 15(7):1597–1604, 2019.
  * [43] Frederick C Frank.  I. liquid crystals. on the theory of liquid crystals.  Discussions of the Faraday Society, 25:19–28, 1958.
  * [44] Sumesh P Thampi, Ramin Golestanian, and Julia M Yeomans.  Instabilities and topological defects in active nematics.  EPL (Europhysics Letters), 105(1):18001, 2014.
  * [45] N Hasselmann and Andreas Sinner.  Interplay of topology and geometry in frustrated two-dimensional heisenberg magnets.  Physical Review B, 90(9):094404, 2014.
  * [46] Hal Tasaki.  Hohenberg-mermin-wagner-type theorems for equilibrium models of flocking.  Physical Review Letters, 125(22):220601, 2020.
  * [47] Supplementary Movies.  https://www.dropbox.com/sh/9e9p1jbkxdoeec3/AACeQ5dVtds0S0x9SleOA9Waa?dl=0.
  * [48] Karsten Kruse, Jean-François Joanny, Frank Jülicher, Jacques Prost, and Ken Sekimoto.  Asters, vortices, and rotating spirals in active gels of polar filaments.  Physical review letters, 92(7):078101, 2004.
  * [49] BI Halperin.  published in physics of defects, proceedings of les houches, session xxxv 1980 nato asi, editors balian, kléman and poirier, 1981.
  * [50] Gene F Mazenko.  Vortex velocities in the o (n) symmetric time-dependent ginzburg-landau model.  Physical review letters, 78(3):401, 1997.
  * [51] Vidar Skogvoll, Jonas Rønning, Marco Salvalaglio, and Luiza Angheluta.  A unified field theory of topological defects and non-linear local excitations.  arXiv preprint arXiv:2302.03035, 2023.
  * [52] Matthew L Blow, Sumesh P Thampi, and Julia M Yeomans.  Biphasic, lyotropic, active nematics.  Physical review letters, 113(24):248303, 2014.
  * [53] Amin Doostmohammadi, Sumesh P Thampi, and Julia M Yeomans.  Defect-mediated morphologies in growing cell colonies.  Phys. Rev. Lett., 117(4):048102, 2016.
  * [54] Thuan Beng Saw, Amin Doostmohammadi, Vincent Nier, Leyla Kocgozlu, Sumesh Thampi, Yusuke Toyama, Philippe Marcq, Chwee Teck Lim, Julia M Yeomans, and Benoit Ladoux.  Topological defects in epithelia govern cell death and extrusion.  Nature, 544(7649):212–216, 2017.
  * [55] T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T Lim, J. M. Yeomans, and B. Ladoux.  Topological defects in epithelia govern cell death and extrusion.  Nature, 544:212–216, 2017.
  * [56] K. Kawaguchi, R. Kageyama, and M. Sano.  Topological defects control collective dynamics in neural progenitor cell cultures.  Nature, 545:327–331, 2017.
  * [57] Katherine Copenhagen, Ricard Alert, Ned S Wingreen, and Joshua W Shaevitz.  Topological defects promote layer formation in myxococcus xanthus colonies.  Nature Physics, pages 1–5, 2020.
  * [58] Johanna Roostalu, Jamie Rickman, Claire Thomas, François Nédélec, and Thomas Surrey.  Determinants of polar versus nematic organization in networks of dynamic microtubules and mitotic motors.  Cell, 175(3):796–808, 2018.
  * [59] Kirsten D Endresen, MinSu Kim, Matthew Pittman, Yun Chen, and Francesca Serra.  Topological defects of integer charge in cell monolayers.  Soft matter, 17(24):5878–5887, 2021.
  * [60] Ekta Makhija, Yang Zheng, Jiahao Wang, Han Ren Leong, Rashidah Binte Othman, Ee Xien Ng, Eng Hin Lee, Lisa Tucker-Kellogg, Yie Hou Lee, Hanry Yu, et al.  Topological defects govern mesenchymal condensations, offering a morphology-based tool to predict cartilage differentiation.  bioRxiv, 2022.
  * [61] Aparna Baskaran and M Cristina Marchetti.  Nonequilibrium statistical mechanics of self-propelled hard rods.  Journal of Statistical Mechanics: Theory and Experiment, 2010(04):P04019, 2010.
  * [62] Aurelio Patelli, Ilyas Djafer-Cherif, Igor S Aranson, Eric Bertin, and Hugues Chaté.  Understanding dense active nematics from microscopic models.  Phys. Rev. Lett., 123(25):258001, 2019.