# Coexistence of defect morphologies in three dimensional active nematics

Pasquale Digregorio<EMAIL_ADDRESS>Departament de Física de la
Matèria Condensada, Universitat de Barcelona, Carrer de Martí i Franqués 1,
08028 Barcelona, Spain  Universitat de Barcelona Institute of Complex Systems
(UBICS), Universitat de Barcelona, 08028 Barcelona, Spain  Cecilia Rorai
CECAM, Centre Européen de Calcul Atomique et Moléculaire, École Polytechnique
Fédérale de Lausanne (EPFL), Batochime, Avenue Forel 2, 1015 Lausanne,
Switzerland  Ignacio Pagonabarraga<EMAIL_ADDRESS>Departament de
Física de la Matèria Condensada, Universitat de Barcelona, Carrer de Martí i
Franqués 1, 08028 Barcelona, Spain  Universitat de Barcelona Institute of
Complex Systems (UBICS), Universitat de Barcelona, 08028 Barcelona, Spain
Federico Toschi  Department of Applied Physics and Science Education,
Eindhoven University of Technology, Den Dolech 2, 5600 MB Eindhoven,
Netherlands  CNR-IAC, I-00185 Rome, Italy

###### Abstract

We establish how active stress globally affects the morphology of disclination
lines of a three dimensional active nematic liquid crystal under chaotic flow.
Thanks to a defect detection algorithm based on the local nematic orientation,
we show that activity selects a crossover length scale in between the size of
small defect loops and the one of long and tangled defect lines of fractal
dimension $2$. This length scale crossover is consistent with the scaling of
the average separation between defects as a function of activity. Moreover, on
the basis of numerical simulation in a 3D periodic geometry, we show the
presence of a network of regular defect loops, contractible onto the
$3$-torus, always coexisting with wrapping defect lines. While the length of
regular defects scales linearly with the emerging active length scale, it
verifies an inverse quadratic dependence for wrapping defects. The shorter the
active length scale, the more the defect lines wrap around the periodic
boundaries, resulting in extremely long and buckled structures.

Introduction. Local active stresses in a highly viscous nematic liquid crystal
lead to the formation of sustained chaotic flow. Such striking behavior has
been observed in experiments of microtubules and kinesin molecular motors
powered by ATP Sanchez  _et al._ (2012); Guillamat  _et al._ (2017); Lemma
_et al._ (2019); Martínez-Prat  _et al._ (2019); Tan  _et al._ (2019), as well
as in bacterial suspensions Dombrowski  _et al._ (2004); Wensink  _et al._
(2012), epithelial cells Blanch-Mercader  _et al._ (2018), synthetic self-
propelled particles Nishiguchi and Sano (2015). As internal energy is
converted into a stirring at the scale of the single constituents, these
systems are referred to as active nematics Marchetti  _et al._ (2013); Gompper
et al (2020).

Activity induces out-of-equilibrium chaotic flow patterns at low Reynolds
number, and this spatiotemporal state shares qualitative properties with
inertial turbulence Rorai  _et al._ (2022). Many attempts have been made to
identify a general framework and the universality classes for these phenomena.
In the context of active systems with nematic order, it has been shown that
scaling regimes exist for the kinetic energy and enstrophy spectra,
characterized by universal exponents Alert  _et al._ (2022). However,
differently than classically driven fluids at high Reynolds number, where a
mechanism of energy cascade plays a key role in sustaining the turbulent
behaviour, no energy transfers across scales have been observed so far for
active chaotic flow Alert  _et al._ (2020); Urzay  _et al._ (2017); Carenza
_et al._ (2020a, b). The transition to a chaotic state in active nematics
proceeds as an instability at the large scale of the nematic orientation,
driven by the presence of local active stresses Edwards and Yeomans (2009). As
the magnitude of the stress increases, smaller wavelength modes become
unstable, generating spatiotemporal chaos.

Such transition generally leads to a proliferation of topological defects.
Contrary to passive systems where defects of opposite sign annihilate, in
active nematics they are also steadily generated because of the energy input
at small scales, and their non trivial dynamics is strongly correlated with
the chaotic flow Thampi and Yeomans (2016); Thampi  _et al._ (2014); Giomi
_et al._ (2013); Shankar  _et al._ (2018). The study of the statistics and
dynamics of topological defects provides fundamental understanding of the
phenomenology and the critical behaviour of systems at equilibrium Mermin
(1979), as well as active systems Digregorio  _et al._ (2022); Shankar and
Marchetti (2019). Regarding chaotic flow driven by activity, numerous studies
have focused on 2D active nematics, revealing that the key properties of the
chaotic flow can be inferred from the behaviour of the $\pm 1/2$
disclinations, as they are closely correlated with the fluid local vorticity
Giomi (2015). Recent results in 3D active nematics show that disclination
lines are spontaneously generated in a regime of well-developed chaotic flow.
The primary defects excited in bulk active nematics are closed loops that
undergo complex dynamics Duclos  _et al._ (2020), as they are locally, and
non-uniformly advected by the flow Kralj  _et al._ (2023); Houston and
Alexander (2022); Duclos  _et al._ (2020); Long  _et al._ (2021); Binysh  _et
al._ (2020); Čopar  _et al._ (2019). Alongside with isolated loops, a complex
network of disclination lines has been observed in experiments and simulations
Duclos  _et al._ (2020); Krajnik  _et al._ (2020). The volume density occupied
by defects increases with activity, and disclination lines continuously
recombine chaotically, providing a coarse-grained qualitative description of
the underlying turbulent flow. However the properties of these defects lines
remain poorly understood and a clear connection with the transition to chaotic
flow is still lacking.

In this work, we design an algorithm to detect nematic topological defects and
carefully reconstruct the disclination lines in 3D active nematics. This gives
us access to the entire spectrum of defect lengths, from small loops to
complex and tangled structures that span the whole system. Depending on their
topology, we distinguish a collection of regular defect loops, contractible
onto a $3$-torus, coexisting with a background of wrapping defects.
Remarkably, these two classes of defects show different behaviour: as activity
increases regular loops get smaller, while the length of the wrapping defects
increases. Furthermore, we study the morphology of regular defect loops and
recognize a characteristic length scale, controlled by activity, that
separates two different scaling regimes.

The model. The description of the nematic orientational order is given in
terms of the symmetric and traceless nematic tensor field
$Q_{ij}=q(n_{i}n_{j}-\delta_{ij}/3)$. $\underline{\bm{Q}}$ is invariant under
local inversion of the nematic director $\bm{n}$, which represents the local
orientation of the nematic field, $q$ being the magnitude of the local nematic
tensor. The nematic is embedded in an incompressible fluid of density $\rho$.
We make use of the well known Beris-Edwards model for the coupled dynamics of
the nematic field and the solvent de Gennes and Prost (1993); Beris and
Edwards (1994); Thampi and Yeomans (2016)

$\displaystyle(\partial_{t}+\bm{u}\cdot\bm{\nabla})\underline{\bm{Q}}-\underline{\bm{S}}$
$\displaystyle=\Gamma\underline{\bm{H}}\rm{,}$ (1)
$\displaystyle(\partial_{t}+\bm{u}\cdot\bm{\nabla})\bm{u}$
$\displaystyle=\frac{1}{\rho}\bm{\nabla}\cdot\underline{\bm{\Pi}}^{\rm{T}}\rm{.}$
(2)

The nematic field is coupled to the flow, and its relaxation is governed by
the Landau-de Gennes free energy with elastic coefficient $K$. Activity
induces local active stresses,
$\underline{\bm{\sigma}}^{\rm{act}}=-\alpha\underline{\bm{Q}}$, with magnitude
$\alpha$ (see Supplementary Material (SM) sm  for a detailed description of
the model). We integrate numerically Eq. (1) and (2) in a 3D periodic cubic
domain using a hybrid lattice Boltzmann (LB) finite-difference method
Marenduzzo  _et al._ (2007); Rorai  _et al._ (2021). The linear size of the
simulation box is set to $L=192$, and a range of box sizes, from $L=96$ up to
$L=768$, is used to perform finite-size analysis. All lengths are given in
lattice units.

The elastic and bulk free energy coefficients are chosen such that $\alpha$
swipes over two orders of magnitude, from $10^{-4}$ to $3\times 10^{-2}$,
still having the stationary state in the nematic phase. Accordingly, the
characteristic scale associated to the active stress,
$l_{\rm{a}}=\sqrt{K/\alpha}$, ranges from $l_{\rm{a}}\sim 10^{-1}$ to
$l_{\rm{a}}\sim 10$. We consider extensile active stresses and flow aligning
nematics, but quantitatively similar results hold for contractile activity and
flow tumbling nematics, as shown in the SM sm .

Defects’ tracking. In order to identify topological defects as singular points
in the nematic orientational order, we extend the general idea introduced in
Zapotocky  _et al._ (1995). For uniaxial nematics in 3D, defects are organised
into disclination lines. Winding around a closed contour that surrounds the
defect in its vicinity, concatenated with the defect line, the director field
defines a map from the path to the projective plane $\mathbb{RP}^{2}$, which
is the unit sphere with antipodal points identified. After continuous
deformation over the sphere, the maps associated to disclinations are the ones
that connect two antipodal points Mermin (1979); Alexander  _et al._ (2012).
Traversing a contour that does not concatenate any defects returns a map that
is contractible to a single point. On this basis, we build a defects’ tracking
algorithm that is parameter-free and returns a set of defect points in space.
Once all defect points have been identified, they are clustered according to a
prescribed criterion of proximity between points on the lattice sm . Within
the clusters, all defects and the corresponding neighbours network define a
connected undirected graph. Notably, each graph built according to this
protocol has an even degree. Under these circumstances, Euler’s theorem
ensures that a map exists between the cluster of defect points and a closed
path that visits all defects exactly once Biggs  _et al._ (1976). All clusters
can therefore be described as closed polygonal curves.

Figure 1: (a) Probability density function of the length of defect lines, for
different values of the active length $l_{\rm a}$ (reported in the legend)
with corresponding errors (shaded area). (b) Volume occupation for all defect
sizes, for the same values of activities. The vertical dashed lines indicate
the location of the average defect length, $\langle\l_{\rm{d}}\rangle$, as
computed from the distributions in panel (a). Panels (c) and (d) show that the
distributions in the first row, with the defect size, collapse onto a
universal curve, when rescaled by the active length $l_{\rm{a}}$.

Results. Fig. 1 (a) displays the probability density function (pdf)
$P(l_{\rm{d}})$ of the defect lines length, measured as the sum of the length
of all edges in the corresponding polygonal curve, for different values of
$l_{\rm{a}}$. All cluster size distributions are peaked around a preferred
size, $\langle l_{\rm{d}}\rangle$, which decreases with activity. Hence,
activity selects the average length of the disclination lines. Fig. 1 (b)
shows the distribution $l_{\rm{d}}P(l_{\rm{d}})$, which is peaked at the same
value $\langle l_{\rm{d}}\rangle(l_{\rm{a}})$, meaning that the latter is also
the size of defects that statistically occupy most of the volume in the
system, while Figs. 1 (c)-(d) show that the defect lines length pdf scale with
the active length $l_{a}$.

Figure 2: (a) Illustrative examples of one regular defect (blue) and one
wrapping (red), respectively contractible and non contractible onto a
$2$-torus, corresponding to a 2D periodic system. (b,c) Two disclination lines
for active length $l_{\rm{a}}\simeq 3$. The thick black line corresponds to
the disclination line tracked, while coloured lines show the same defect,
unwrapped with respect to the periodic boundaries. In panel (a) the blue line
corresponds to a regular loop, while in panel (b) the red line corresponds to
a wrapping one. Note that only one defect line, out of the large number
present, is shown for each configuration.

In turn, Figs 1(a) and 1(b) show that the system is populated by a wide
distribution of defect lines, characterised by wide tail of the pdfs, whose
length largely exceeds the system size, and whose excitation is allowed due to
the system periodicity. Under periodic boundary conditions (PBC), a defect
line lives on a $3$-torus and can either close connecting to its initial point
or to any other periodic images of the initial point, defining two different
topological classes, as illustrated in Fig. 2 (a) for a $2$-torus. In the
former case, referred to as regular defect, the contour can be contracted
continuously to a point, while in the latter, referred to as wrapping, this is
not possible. Fig. 2 (b) and (c) show an example of a regular and wrapping
defect, respectively.

These two classes of defects contribute in a qualitatively different way to
the defect tangle. Fig. 1 (c) shows that the defect pdf scales as
$l_{\rm{a}}P(l_{\rm{d}}/l_{\rm{a}})$, except for the pfds’ tails. The pdf of
the occupied volume fulfils the same scaling, as shown in Fig. 1 (d).

Figure 3:  (a) Average defect length, corresponding to the first moment of the
length distribution of Fig. 1, as a function of the active length (grey
crosses). The same quantity is reported, separately, for wrapping defects (red
circles) and regular defects (blue triangles), with corresponding error bars.
The corresponding scaling behaviour is also shown with dashed lines. (b)
Average number of defect lines per configuration with the corresponding
scaling behaviour. (c) Probability density function of defects length at
$l_{\rm{a}}=1$, for different system sizes, indicated in the legend. The
continuous lines represent the distribution for all defects, while symbols
represent the distribution for regular defects only. (d) Volume fraction
occupied by wrapping and regular defects, as a function of the system size,
for $l_{\rm{a}}=1$.

Moreover, Figs. 3 (a,b) display the single defect average length, and average
number of defect lines, respectively, as a function of $l_{\rm{a}}$, for
wrapping defects (red), for regular defects (blue), and for all defects
together (grey). The average linear size of regular defects scales as $\langle
l_{\rm{d}}\rangle\sim l_{\rm{a}}$, and the average number of lines scale as
$\langle N_{\rm{d}}\rangle\sim l_{\rm{a}}^{-3}$ while wrapping defects’
average size grows as activity increases. In particular, we verify that
$\langle l_{\rm{d}}\rangle\sim l_{\rm{a}}^{-2}$, while the average number of
these defects is always of order $\langle N_{\rm{d}}\rangle\sim 1$,
independent on activity.

Note that, for $l_{\rm{a}}\gtrsim 3$, the average number of defects with the
two topologies is approximately the same, consistent with the fact that the
system is close to the transition from chaotic to laminar flow. Larger active
lengths are no longer able to excite defects and the system relaxes to a long-
range ordered state. This behaviour of $N_{\rm{d}}$ results in a crossover
from $\langle l_{\rm{d}}\rangle\sim l_{a}$ to $\langle l_{\rm{d}}\rangle\sim
l_{a}^{-2}$ for the average length of all defects in the regime of large
active length, as evident in the behaviour of the grey points in Fig. 3 (a).
Close to the transition to the chaotic state the average length for the two
classes of defects is comparable, meaning that the easiest defects to sustain
are system-size long line defects, either regular or wrapping. As activity
increases, the characteristic sizes of the two classes move away from each
other, such that smaller regular defects proliferate and wrapping defects get
longer.

Figs. 3 (c,d) display the impact of system size on both classes of defects.
Panel (c) displays $P(l_{\rm{d}})$ for the same active length $l_{\rm{a}}=1$
and different system sizes, ranging from $L=96$ to $L=768$. The symbols
represent the distribution of regular defects, while the continuous lines
represent the distribution for all defects. The shape of the pdf for regular
defects does not change with $L$ and the peak is located always at the same
value of $l_{\rm{d}}$. However, larger system sizes allow to explore a wider
range in the defects length, as is clear from the power law tail, and the
average size of the wrapping defects is displaced to larger lengths as the
system size increases. In fact, the total length of a single wrapping defect,
although it can exceed the system size with no a priori limits, is bounded
from below: a defect shorter that the system size cannot wrap. Even though the
probability of finding wrapping defects of a given length $l_{\rm{d}}$
decreases upon increasing the system size, it does not vanish, as $\langle
N_{\rm{d}}\rangle\sim 1$ up to the largest system sizes considered. This
behaviour is further clarified in Fig. 3 (d), which displays the overall
volume fraction occupied by regular and wrapping defects as a function of the
system size, defined as the ratio between the total length of the defects
$L_{\rm d}$ and the system volume, $\rho_{\rm d}\sim L_{\rm d}/L^{3}$. For all
sizes considered and $\l_{\rm a}\simeq 1$, the contribution to the volume
fraction from the wrapping defects is always larger than the one coming from
the regular defects, and the two are approximately constant, ruling out finite
size effects. Hence, in the regime where defects proliferate, wrapping defects
are always present and they can predominate in terms of total volume occupied.
The number of wrapping defects is small, suggesting that they behave as a
topologically-preserved condensate that coexists with the broad distribution
of regular defects. We stress that wrapping disclinations are present and
statistically relevant even though the topology of defects lines is not fixed
and can turn from non-contractible to contractible and vice versa through
local recombination of the tangle (see the movie in SM sm ).

Figure 4: (a) Squared radius of gyration, ${R_{\rm
g}^{2}=\sum_{i,j}{|\bm{r}_{i}-\bm{r}_{j}|^{2}}/2N^{2}}$, as a function of the
defect length, for different values of activity. (b) Average curvature of
defect lines as a function of $l_{\rm{d}}$. (c) Defect length and radius of
gyration at the crossover between the two scaling regimes of panel (a), as a
function of activity, with the corresponding scaling behaviour shown with
dashed lines. (d) Average value of the estimated defects’ separation, as a
function of $l_{a}$, represented by the grey crosses, with the corresponding
error. A power law $l_{\rm a}^{0.8}$ is shown with a dashed grey line, while
$\sim l_{\rm{a}}$ is shown with a dashed red line for comparison. The scaling
behaviour $\lambda_{\rm d}\sim l_{\rm a}^{0.8}$ is the best fit, as shown in
the inset.

Activity strongly affects the defects’ local curvature and their overall
morphology. Fig. 4 (a) displays the average squared radius of gyration,
$\langle R_{{\rm g}}^{2}\rangle$, of defect lines as a function of the their
length, $l_{{\rm d}}$, for different values of $l_{a}$. For all activities
considered, we clearly observe two separate length regimes, correlated to two
different exponents in the power law scaling
$R_{\rm{g}}^{2}\sim{l_{\rm{d}}}^{2/d}$, with $d$ the fractal dimension of the
defect structures. Small defects are characterised by $d=1$, which means that
they are arranged into closed loops with a persistence length of the same
order of the defect length. From direct observation we conclude that these
objects are rather planar, rigidly moving, small defect rings of regular
shape. Large defects show $d=2$, which, on the contrary, is associated to
buckled structures with no preferred spatial direction. In particular, the
value of the exponent establishes that, independent of activity, in the limit
of very large sizes, disclination lines are arranged as closed random walks.

Fig. 4 (b) displays the defects’ average local curvature,
$\langle\kappa\rangle$, as a function of $l_{{\rm d}}$. For very small defects
$\langle\kappa\rangle$ decreases, approaching the limiting power law
$\langle\kappa\rangle\sim\l_{\rm{d}}^{-1}$, characteristic of planar shapes.
At larger sizes $\langle\kappa\rangle$ reaches a plateau, corresponding to
constant curvature, tangled structures. The value of the plateau depends on
activity and, consistently with the results of $\langle R_{{\rm
g}}^{2}\rangle$, it increases with activity. For the largest active length
$l_{\rm{a}}=5.23$, we do not observe the first branch of the curve at smaller
defects length, compatibly with the fact that, close to the transition to the
chaotic regime, defects are long stretched lines of length comparable with the
system size. The crossover scale $l_{\rm{d}}^{*}$ between the two different
structures scales with activity, $l_{\rm{d}}^{*}\sim{l_{\rm{a}}}^{\gamma}$
with $\gamma\simeq 0.8$, as shown in Fig. 4 (c).

In light of the above characterisation of the defects’ morphology, we now look
at the way the whole tangle is spatially arranged. We measure the total length
$L_{\rm{d}}$ of the defects tangle and the corresponding defect density
$\rho_{\rm d}$, and we compute the average separation between defects as
$\lambda_{\rm d}\sim\rho_{\rm d}^{-d/2}$. Assuming that $\lambda_{\rm d}\sim
l_{\rm a}$, as proposed in previous works Hemingway  _et al._ (2016); Krajnik
_et al._ (2020), we expect $\rho_{\rm d}^{-1/2}\sim l_{\rm a}$ for small
linear defects of $d=1$ and $\rho_{\rm d}^{-1/2}\sim l_{\rm a}^{1/2}$ for
large buckled ones of fractal dimension $d=2$. Comparing the location of the
crossover length with the defects’ length pdfs (see Fig. 1), most defects
populate the region around the crossover length. No one of the two asymptotic
morphologies dominate the defects tangle and we find an intermediate value of
the exponent, $\rho_{\rm d}^{-1/2}\sim l_{\rm a}^{0.8}$.

Conclusions. As activity drives chaotic flows in nematic liquid crystals
disclination lines proliferate. We have thoroughly characterized the nature of
these defects and shown their morphological complexity. In a 3D periodic
system defects belong to two topologically different classes, regular and
wrapping disclination lines. While the number of regular defect lines
increases with activity, the number of wrapping defects is of order one and
can be regarded as a topologically preserved condensate.

The average length of the wrapping defects is proportional to the system size
and increases upon decreasing the characteristic length scale, $l_{\rm a}$,
while regular disclination loops get smaller and their length depends linearly
on $l_{\rm a}$. The morphology and spatial organization of the regular defects
strongly depends on their length. Small defect rings of planar shape and
finite self-propelling speed coexist with long buckled defects that show the
typical scaling exponent of random walks.

Acknowledging this richness in defect morphology allows us to better
rationalize the global properties of the defect tangle. Specifically, the fact
that the majority of defects belong to the crossover region between these two
asymptotic classes, impacts the observed scalings, such as the dependence of
the defect average separation length with activity. This opens new questions
about the dynamics of disclination lines in relation to the heterogeneity in
morphology observed and their topology in a periodic 3D system.

Acknowledgments P.D. acknowledges support from the Juan de la Cierva-Formación
program. C.R. acknowledges funding from the European Union’s Horizon 2020
research and innovation program under the Marie Sklodowska-Curie Grant
Agreement No. 754462. I.P. acknowledges support from Ministerio de Ciencia,
Innovación y Universidades MCIU/AEI/FEDER for financial support under grant
agreement PID2021-126570NB-100 AEI/FEDER-EU, from Generalitat de Catalunya
under Program Icrea Acadèmia and project 2021SGR-673. This work was supported
by grants from the Swiss National Supercomputing Centre (CSCS), project IDs
s1079 and s1194, and MareNostrum Supercomputer at Barcelona Supercomputing
Center (BSC).

## References

  * Sanchez  _et al._ (2012) T. Sanchez, D. Chen, and S. DeCamp et al., Nature 491, 431–434 (2012).
  * Guillamat  _et al._ (2017) P. Guillamat, J. Ignés-Mullol, and F. Sagués, Nat Commun 8 (2017), 10.1038/s41467-017-00617-1.
  * Lemma  _et al._ (2019) L. M. Lemma, S. J. DeCamp, Z. You, L. Giomi, and Z. Dogic, Soft Matter 15, 3264 (2019).
  * Martínez-Prat  _et al._ (2019) B. Martínez-Prat, J. Ignés-Mullol, and J. Casademunt et al, Nat. Phys. 15, 362–366 (2019).
  * Tan  _et al._ (2019) A. Tan, E. Roberts, and S. Smith et al, Nat. Phys. 15, 1033–1039 (2019).
  * Dombrowski  _et al._ (2004) C. Dombrowski, L. Cisneros, S. Chatkaew, R. E. Goldstein, and J. O. Kessler, Phys. Rev. Lett. 93, 098103 (2004).
  * Wensink  _et al._ (2012) H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, Proceedings of the National Academy of Sciences 109, 14308 (2012).
  * Blanch-Mercader  _et al._ (2018) C. Blanch-Mercader, V. Yashunsky, S. Garcia, G. Duclos, L. Giomi, and P. Silberzan, Phys. Rev. Lett. 120, 208101 (2018).
  * Nishiguchi and Sano (2015) D. Nishiguchi and M. Sano, Phys. Rev. E 92, 052309 (2015).
  * 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. 85, 1143 (2013).
  * Gompper et al (2020) G. Gompper et al, J. Phys.: Condens. Matter 32 (2020), 10.1088/1361-648X/ab6348.
  * Rorai  _et al._ (2022) C. Rorai, F. Toschi, and I. Pagonabararga, Physical Review Letters 129, 218001 (2022).
  * Alert  _et al._ (2022) R. Alert, J. Casademunt, and J.-F. Joanny, Annual Review of Condensed Matter Physics 13, 143 (2022).
  * Alert  _et al._ (2020) R. Alert, J. Joanny, and J. Casademunt, Nat. Phys. 16, 682–688 (2020).
  * Urzay  _et al._ (2017) J. Urzay, A. Doostmohammadi, and J. Yeomans, Journal of Fluid Mechanics 822, 762–773 (2017).
  * Carenza  _et al._ (2020a) L. Carenza, L. Biferale, and G. Gonnella, EPL 132, 44003 (2020a).
  * Carenza  _et al._ (2020b) L. N. Carenza, L. Biferale, and G. Gonnella, Phys. Rev. Fluids 5, 011302 (2020b).
  * Edwards and Yeomans (2009) S. A. Edwards and J. M. Yeomans, Europhysics Letters 85, 18008 (2009).
  * Thampi and Yeomans (2016) S. Thampi and J. Yeomans, Eur. Phys. J. Special Topics 225, 651 (2016).
  * Thampi  _et al._ (2014) S. Thampi, R. Golestanian, and Y. J.M., EPL 105, 18001 (2014).
  * Giomi  _et al._ (2013) L. Giomi, M. J. Bowick, X. Ma, and M. C. Marchetti, Phys. Rev. Lett. 110, 228101 (2013).
  * Shankar  _et al._ (2018) S. Shankar, S. Ramaswamy, M. C. Marchetti, and M. J. Bowick, Phys. Rev. Lett. 121, 108002 (2018).
  * Mermin (1979) N. D. Mermin, Rev. Mod. Phys. 51, 591 (1979).
  * Digregorio  _et al._ (2022) P. Digregorio, D. Levis, L. F. Cugliandolo, G. Gonnella, and I. Pagonabarraga, Soft Matter 18, 566 (2022).
  * Shankar and Marchetti (2019) S. Shankar and M. C. Marchetti, Phys. Rev. X 9, 041047 (2019).
  * Giomi (2015) L. Giomi, Phys. Rev. X 5, 031003 (2015).
  * Duclos  _et al._ (2020) G. Duclos, R. Adkins, D. Banerjee, M. Peterson, M. Varghese, I. Kolvin, A. Baskaran, R. Pelcovits, T. Powers, A. Baskaran, F. Toschi, M. Hagan, S. Streichan, V. Vitelli, D. Beller, and Z. Dogic, Science 367, 1120 (2020).
  * Kralj  _et al._ (2023) N. Kralj, M. Ravnik, and i. c. v. Kos, Phys. Rev. Lett. 130, 128101 (2023).
  * Houston and Alexander (2022) A. J. H. Houston and G. P. Alexander, Phys. Rev. E 105, L062601 (2022).
  * Long  _et al._ (2021) C. Long, X. Tang, R. L. B. Selinger, and J. V. Selinger, Soft Matter 17, 2265 (2021).
  * Binysh  _et al._ (2020) J. Binysh, i. c. v. Kos, S. Čopar, M. Ravnik, and G. P. Alexander, Phys. Rev. Lett. 124, 088001 (2020).
  * Čopar  _et al._ (2019) S. Čopar, J. Aplinc, i. c. v. Kos, S. Žumer, and M. Ravnik, Phys. Rev. X 9, 031051 (2019).
  * Krajnik  _et al._ (2020) Z. Krajnik, Z. Kos, and M. Ravnik, Soft Matter 16, 9059 (2020).
  * de Gennes and Prost (1993) P.-G. de Gennes and J. Prost,  _The Physics of Liquid Crystals_ (Oxford University Press, New York, 1993).
  * Beris and Edwards (1994) A. N. Beris and B. J. Edwards,  _Thermodynamics of Flowing Systems: with Internal Microstructure_ (Oxford University Press, 1994).
  * (36) See Supplemental Material at [url], which contains Refs. Rorai  _et al._ (2021); Hierholzer and Wiener (1873); Fleischner (1990), for additional details on the numerical methods used and the analysis of defects shown.
  * Marenduzzo  _et al._ (2007) D. Marenduzzo, E. Orlandini, M. E. Cates, and J. M. Yeomans, Phys. Rev. E 76, 031921 (2007).
  * Rorai  _et al._ (2021) C. Rorai, F. Toschi, and I. Pagonabarraga, Phys. Rev. Fluids 6, 113302 (2021).
  * Zapotocky  _et al._ (1995) M. Zapotocky, P. M. Goldbart, and N. Goldenfeld, Phys. Rev. E 51, 1216 (1995).
  * Alexander  _et al._ (2012) G. P. Alexander, B. G.-g. Chen, E. A. Matsumoto, and R. D. Kamien, Rev. Mod. Phys. 84, 497 (2012).
  * Biggs  _et al._ (1976) N. Biggs, E. Lloyd, and R. Wilson,  _Graph Theory, 1736–1936_ (Clarendon Press, 1976).
  * Hemingway  _et al._ (2016) E. J. Hemingway, P. Mishra, M. C. Marchetti, and S. M. Fielding, Soft Matter 12, 7943 (2016).
  * Hierholzer and Wiener (1873) C. Hierholzer and C. Wiener, Mathematische Annalen 6, 30 (1873).
  * Fleischner (1990) H. Fleischner,  _X.1 Algorithms for Eulerian Trails, in Eulerian Graphs and Related Topics: Part 1, Volume 2_ (Elsevier, 1990).
  * (45) Movie available at https://www.dropbox.com/scl/fi/seze0j6gqh001tygy2i57/regular-wrapping_defs_alpha0.001.mov?rlkey=0jp95zigxhtzhwlky7elx9uww&dl=0.

Supplemental Material for
“Coexistence of defect morphologies in three dimensional active nematics”

##  S1 Model and numerical simulations

Figure S1: (a) A $3\times 3$ vertices portion of the system with the
corresponding director field. Arrows are double-headed to represent the
nematic symmetry. The face highlighted in blue colour is in the ordered phase,
while the one highlighted in red corresponds to a defect, identified by the a
red asterisk. Black faces also correspond to topological defect points, shown
as black asterisks. The black dashed line is a disclination line joining the
four point defects identified in the picture. (b) Mapping onto a unit sphere
of the nematic field around the faces highlighted in blue and red in panel
(a). The black line follows the path of the red arrows onto the spherical
surface.

We write the equation of motion used for the evolution of a nematic liquid
crystal embedded in a fluid in the presence of local active stress, and
provide additional details about the model and the set of parameters adopted.

$\displaystyle(\partial_{t}+\bm{u}\cdot\bm{\nabla})\underline{\bm{Q}}-\underline{\bm{S}}$
$\displaystyle=\Gamma\underline{\bm{H}}\rm{,}$ (S1)
$\displaystyle(\partial_{t}+\bm{u}\cdot\bm{\nabla})\bm{u}$
$\displaystyle=\frac{1}{\rho}\bm{\nabla}\cdot\underline{\bm{\Pi}}^{\rm{T}}\rm{.}$
(S2)

Here, the nematic field $\underline{\bm{Q}}$ is advected by the flow and
coupled to the flow gradients through the co-rotational term

$\displaystyle\underline{\bm{S}}$
$\displaystyle=(\xi\underline{\bm{E}}+\underline{\bm{\Omega}})\cdot(\underline{\bm{Q}}+\underline{\bm{I}}/3)+(\underline{\bm{Q}}+\underline{\bm{I}}/3)\cdot(\xi\underline{\bm{E}}-\underline{\bm{\Omega}})$
(S3) $\displaystyle-2\xi(\underline{\bm{Q}}+\underline{\bm{I}}/3)\
{\rm{tr}}(\underline{\bm{Q}}\cdot\underline{\bm{\nabla}\bm{u}})\rm{.}$

The character of this coupling is set by the flow aligning parameter $\xi$,
while $\underline{\bm{E}}$ and $\underline{\bm{\Omega}}$ are the strain rate
and the vorticity tensor, respectively. All the results in the main text are
given for the choice $\xi=0.7$, which corresponds to a regime of flow
aligning, where the director field tends to align to an underlying shear flow.
For values $\xi<0.5$ the nematic director field tends to rotate with respect
to a shear flow. This regime is called flow tumbling, and we report later in
this Supplementary Material results for this different coupling regime. The
nematic relaxation is governed by the molecular field
$\underline{\bm{H}}=-\delta\mathcal{F}/\delta\underline{\bm{Q}}+\underline{\bm{I}}\
{\rm{tr}}(\delta\mathcal{F}/\delta\underline{\bm{Q}})$, with Landau-de Gennes
free energy

$\mathcal{F}=\int{\rm{dV}}\biggl{[}\frac{K}{2}(\bm{\nabla}\underline{\bm{Q}})^{2}+\frac{A}{2}{\rm{tr}}(\underline{\bm{Q}}^{2})+\frac{B}{3}{\rm{tr}}(\underline{\bm{Q}}^{3})+\frac{C}{4}({\rm{tr}}(\underline{\bm{Q}}^{2}))^{2}\biggr{]}\rm{.}$
(S4)

The stress tensor is given by
$\underline{\bm{\Pi}}=-P\underline{\bm{I}}+2\eta\underline{\bm{E}}+\underline{\bm{\sigma}}+\underline{\bm{\sigma}}^{\rm{act}}$,
with

$\displaystyle\underline{\bm{\sigma}}$
$\displaystyle=2\xi(\underline{\bm{Q}}+\underline{\bm{I}}/3)\
{\rm{tr}}(\underline{\bm{Q}}\cdot\underline{\bm{H}})-\xi\underline{\bm{H}}\cdot(\underline{\bm{Q}}+\underline{\bm{I}}/3)$
(S5)
$\displaystyle-\xi(\underline{\bm{Q}}+\underline{\bm{I}}/3)\cdot\underline{\bm{H}}-\bm{\nabla}\underline{\bm{Q}}\colon\frac{\delta\mathcal{F}}{\delta\bm{\nabla}\underline{\bm{Q}}}$
$\displaystyle+\underline{\bm{Q}}\cdot\underline{\bm{H}}-\underline{\bm{H}}\cdot\underline{\bm{Q}}\
\rm{,}$

and $\underline{\bm{\sigma}}^{\rm{act}}=-\alpha\underline{\bm{Q}}$. The
density and viscosity of the fluid are set to $\rho=1.0$ and $\eta=0.167$,
respectively. For the nematic liquid crystal, we set $K=0.002738$, $A=-0.2K$,
$C=20.0K$, $B=-C$, $\Gamma=8.28$. All parameters values are given in lattice
units. Times are expressed in terms of the characteristic relaxation time
$\tau=L^{2}/(K\Gamma)$, with $L$ the linear size of the system. This choice of
the parameters returns $q_{0}\simeq 0.5$ in the ordered phase, and it allows
us to keep a constant characteristic defect core size $l_{\rm{c}}\simeq 1$
Rorai  _et al._ (2021).

All measures are performed over $5$ independent simulations for each set of
parameters used. Each simulation runs over $10^{6}$ steps and configurations
are sampled every $10^{4}$ steps.

##  S2 The defect tracking algorithm

Topological defects in 3D nematics are named disclinations, as they correspond
to singularity in the orientational nematic field. In order to identify
singular points of the nematic orientation field, we examine all square
contours corresponding to the faces of each 3D voxel in the lattice and
establish if the face encloses a disclination, based on the construction
depicted in Fig. S1. While circuiting in an ordered way around the four
vertices of the face, we construct a map of the local director field onto a
unit sphere, as shown in Fig. S1 (b). As we deal with nematic symmetry, all
opposite points on the sphere are identified. Therefore, once we move to a new
vertex, we choose, between the two equivalent (opposite) orientations of the
director, the one that is closer to the previous one. Using this procedure,
the last step from the forth vertex back to the initial one returns two
possible outcomes: either we go back to the initial point on the sphere (see
the set of blue arrows in Fig. S1 (b)), or we reach its antipodal (red arrows
and the black path in Fig. S1 (b)). The latter case corresponds to the
presence of a topological defect, which we place by convention in the center
of the square, as indicated by the red asterisk in Fig. S1 (a).

Figure S2: Defect lines, as constructed with the method described in text for
active length (a) $l_{\rm a}=1.65$, (b) $l_{\rm a}=0.955$, (c) $l_{\rm
a}=0.523$. Regular defects are shown in blue, while wrapping defects are shown
in red.

Once defect points have been found and a neighbourhood relationship has been
assigned, we can group them into connected clusters. In particular, we
consider as neighbours two defect points if they belong to two adjacent faces
or two opposite faces of the same voxel. As described in the main text, we
note that all clusters have even degree, which means that all nodes are
associated to an even number of edges. We can thus look for a closed path that
visits all defect points in a cluster, travelling along each edge only once.
We identify such path by making use of the so-called Hierholzer’s algorithm
Hierholzer and Wiener (1873); Fleischner (1990). As shown in Fig. S2, this
procedure returns a complete set of defect lines, that perfectly matches the
defect pattern identified with regions of small scalar order parameter $q$.
See in the main text a detailed description of regular and wrapping defect
lines, represented respectively in light blue and red in Fig. S2.

##  S3 Flow tumbling regime and contractile activity

Figure S3: (a,d) Probability density function of disclinations’ length, as a
function of $l_{\rm d}/l_{\rm a}$. (b,e) Average length of the disclinations
and (c,f) their average number, as a function of the active length. The first
row shows the results with contractile activity $\alpha<0$ and $\xi=0.7$,
while $\alpha>0$ and $\xi=0.3$ in the second row.

All the results reported in the main text refer to a regime of positive active
parameter $\alpha$ and flow aligning coupling between the nematic field and
underlying shear flow $\xi>0.5$. There are however other choices available,
and we will provide in this Section a set of results, corresponding separately
to negative $\alpha$ and $\xi<0.5$. It is well known that two main categories
of swimmers exist, which are named pushers and pullers, depending whether they
push or pull the flow along their swimming direction, respectively. The sign
of the parameter $\alpha$ in our model captures this different behavior. In
particular, positive values of $\alpha$ generate a puller-type local flow
field and in this regime the active stress $-\alpha\underline{\bm{Q}}$ is
called extensile. On the other hand, negative values of $\alpha$ correspond to
puller swimmers and the activity is named contractile. We provide in Fig. S3
(a-c) the scaled probability density function of the defect length, the
average defect length and the average number of disclination lines, in the
regime of contractile activity, for different values of active length $l_{\rm
a}=\sqrt{K/|\alpha|}$. Fig. S3 (d-f) show the same quantities, in a regime of
extensile activity, for the same set of values of $\l_{\rm a}$, but for
$\xi=0.3$, in the flow tumbling regime.

Both the results with contractile activity and in the flow tumbling regime do
not show notable differences from the ones with extensile activity and flow
aligning, which are the focus of the analysis in the main text. The maximum in
the probability density function of the defect length scales with active
length, identifying a clear dependence of the length of the regular defects on
activity, is not modified either by the contractile nature of the activity or
the flow-tumbling character of the coupling to the flow. We also observe the
presence of system-spanning wrapping defects, at all values of activity
explored. As shown in Fig S3 (b,c) for contractile activity and (e,f) for flow
tumbling, the average length of these percolating network of defects increases
as $\alpha$ grows in magnitude, with a power law dependence on active length
$\langle l_{\rm d}\rangle\sim l_{\rm a}^{-2}$. The average number of
percolating defects is always of order $\langle N_{\rm d}\rangle\sim 1$,
independently on activity.

##  S4 Scaling crossover

Figure S4: (a) Disclination loops of length $l_{\rm d}<l_{d}^{*}$ and (b)
$l_{\rm d}<l_{d}^{*}$. (c) Average radius of gyration as a function of the
defect length, for different values of activity, given in the legend. The same
curves as in Fig. 4 are shown, but in separate panels and with the
corresponding statistical uncertainty. (d) Average local curvature, as a
function of the defect length, for the same values of activity, as shown in
Fig. 4, with corresponding errors.

We provide supplemental pictures concerning the crossover length between the
two scaling regimes of defect lines, reported in Fig. 4 of the main text and
therein commented. Small defects are closed loops of approximately planar
shape, and are coherently described bu a scaling exponent $R_{\rm g}^{2}\sim
l_{\rm d}^{2}$. See the typical shape of these defects in Fig. S4 (a). On the
other hand, larger defects show a scaling exponent $R_{\rm g}^{2}\sim l_{\rm
d}$. Therefore, their morphology is the one of a random walk in 3D, with a
given curvature that depends on activity. Examples of these defects are given
in Fig. S4 (b). The crossover length between these two groups of defects
strongly depends on activity, according $l_{\rm d}^{*}\sim l_{\rm a}^{0.8}$
(see Fig. 4 in the main text). We locate the crossing point with an iterative
protocol. We make a first guess of the point $l_{\rm{d}}^{*}$ and fit the two
branches restricted to the range $l_{\rm{d}}<l_{\rm{d}}^{*}$ and
$l_{\rm{d}}>l_{\rm{d}}^{*}$, respectively. This returns a new estimate of the
crossing point $l_{\rm{d}}^{*}$. The same procedure is iterated until
convergence. Fig. S4 (c,d) show the average radius of gyration and the average
local curvature, respectively, as a function of the defect length, with
corresponding uncertainty on the average.

##  S5 Supplemental movie

The supplemental movie “regular-wrapping_defs_alpha0.001.mov” mov  show the
typical dynamics of the disclination lines, as tracked with our algorithm. The
evolution starts from a nematically ordered state and active length $l_{\rm
a}=1.65$. Very quickly the systems starts to develop topological defects, that
grow in number and length, until they reach a stationary global density with
regular and wrapping lines. The panel on the left shows the whole defect
tangle with regular defects coloured in blue and wrapping defects coloured in
red. In panel on the right the same evolution is shown, with only regular
defects visible. Defect lines evolve in time and continuously recombine
between each other, in a way that the length and shape of single lines are not
conserved quantities. The topology of defect is not protected and they can
turn from regular to wrapping and vice versa with no a priori constraints.