# Helical Flow States in Active Nematics

Ryan Keogh  School of Physics and Astronomy, The University of Edinburgh,
Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK.  Santhan Chandragiri
Department of Chemical Engineering, Indian Institute of Technology Madras,
Chennai 600036, India.  Benjamin Loewe  School of Physics and Astronomy, The
University of Edinburgh, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK.
Tapio Ala-Nissila  MSP Group, QTF Centre of Excellence, Department of Applied
Physics, Aalto University, P.O. Box 11000, FI-00076 Aalto, Espoo, Finland
Interdisciplinary Centre for Mathematical Modelling, Department of
Mathematical Sciences, Loughborough University, Loughborough LE11 3TU, United
Kingdom  Sumesh Thampi  Department of Chemical Engineering, Indian Institute
of Technology Madras, Chennai 600036, India.  Tyler N. Shendruk
<EMAIL_ADDRESS>School of Physics and Astronomy, The University of
Edinburgh, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK.

###### Abstract

We show that confining extensile nematics in 3D channels leads to the
emergence of two self-organized flow states with nonzero helicity. The first
is a pair of braided anti-parallel streams—this double helix occurs when the
activity is moderate, anchoring negligible and reduced temperature high. The
second consists of axially aligned counter-rotating vortices—this grinder
train arises between spontaneous axial streaming and the vortex lattice. These
two unanticipated helical flow states illustrate the potential of active
fluids to break symmetries and form complex but organized spatio-temporal
structures in 3D fluidic devices.

Figure 1: Active dissipative flow structures exhibiting rotational flow. (a)
Double helix: Two streams of anti-parallel flow braid around one another for
parameters
$\tilde{P}=\\{\tilde{K}^{-1/2},\tilde{W},\tilde{T}^{1/2}\\}=\\{25,1.5,9.47\\}$.
Coloured by axial velocity $u_{\parallel}$ (Movie 1). (b) Grinder train:
Procession of spontaneously drifting, counter-rotating vortices aligned
axially down the duct for $\tilde{P}=\\{20,250,0\\}$. Coloured by axial
vorticity $\omega_{\parallel}$ (Movie 2). (c) Ceilidh vortex: Vortex lattice
oriented in the transverse direction for $\tilde{P}=\\{22,250,0\\}$. Coloured
by transverse vorticity $\omega_{\perp}$ (Movie 3).

Because active materials are composed of microscopic constituents, which
locally transmute internal energy into mechanical forces that are dissipated
by the surrounding medium, they can possess broken symmetries and emergent
dynamic properties on macroscopic scales. These self-organized, active
dissipative structures have been observed to take many forms. The premier
example of activity-induced dissipative structures is flocking of polar
particles Vicsek  _et al._ (1995); Toner and Tu (1995, 1998), in which a
liquid-gas-like transition Solon  _et al._ (2015) from disorderly to
collective motion is associated with spontaneous spatial phase separation
Chaté  _et al._ (2008). In addition to coexistence, hexatic, smectic and solid
phases can originate from activity Digregorio  _et al._ (2018); Chattopadhyay
_et al._ (2021). Likewise, activity can generate orientational order Santhosh
_et al._ (2020); Giordano  _et al._ (2021), including active alignment of
motile rods Bär  _et al._ (2020); Nagel  _et al._ (2020) and swimming
bacilliforms Nishiguchi  _et al._ (2017).

Figure 2: Probability distributions of vorticity for active dissipative flow
structures at parameters
$\tilde{P}=\\{\tilde{K}^{-1/2},\tilde{W},\tilde{T}^{1/2}\\}$. The components
are the vorticity taken about the center line in the axial direction
$\omega_{\parallel}$ (red lines) and in the two transverse directions
$\omega_{\perp_{1,2}}$ (blue and cyan lines) averaged over ten random
initializations. (a) Double helix ($\tilde{P}=\\{25,1.5,9.47\\}$) possesses
non-zero vorticity in both transverse directions $\omega_{\perp_{1,2}}$. (b)
Grinder train ($\tilde{P}=\\{20,250,0\\}$). (c) Ceilidh vortex lattice
($\tilde{P}=\\{21,250,0\\}$), possesses non-zero vorticity only in one of the
two transverse directions. (d) Turbulent flow ($\tilde{P}=\\{25,250,0\\}$) is
nearly isotropic.

Activity can also induce self-organized flow states Voituriez  _et al._
(2006). In confining channels, an extensile active stress induced instability
to bend perturbations leads to self-sustained spontaneous flows Voituriez  _et
al._ (2005); Ramaswamy and Rao (2007); Edwards and Yeomans (2009); Fielding
_et al._ (2011); Ravnik and Yeomans (2013) and spontaneous circulations
Woodhouse and Goldstein (2012); Jiang  _et al._ (2017); Theillard  _et al._
(2017) when the characteristic activity length scale $\ell_{\text{act}}$ is
comparable to the confinement size $L$, as seen in numerous experimental
systems Wioland  _et al._ (2013); Lushi  _et al._ (2014); Beppu  _et al._
(2017); Norton  _et al._ (2018); Hamby  _et al._ (2018); Opathalage  _et al._
(2019); Hardoüin  _et al._ (2020); Thijssen  _et al._ (2021). Vortex lattices
arise when the characteristic length scale matches the confinement, which can
be due to an array of obstacles or cavities Thampi  _et al._ (2016); Wioland
_et al._ (2016a, b); Nishiguchi  _et al._ (2018); Reinken  _et al._ (2020);
James  _et al._ (2021), bound to spherical surfaces Khoromskaia and Alexander
(2017), or channel confinement Shendruk  _et al._ (2017); Doostmohammadi  _et
al._ (2017); Zhang  _et al._ (2020); Samui  _et al._ (2021). Local
circulations arise because $\ell_{\text{act}}$ represents the characteristic
vorticity length scale Hemingway  _et al._ (2016), associated with a peak in
the enstrophy distribution Wensink  _et al._ (2012); Bratanov  _et al._
(2015); Carenza  _et al._ (2020); Alert  _et al._ (2020), and a lattice of
vortices represents stress-free solutions that lie at the interface of the
stable and unstable modes in minimal continuum models of active fluids Słomka
and Dunkel (2017). Such 2D flow states represent the emergence of a broken
translational symmetry of the vorticity field.

However, 3D fluids possess the capacity for flows that are not allowed in 2D.
Helicity, for instance, is identically zero in 2D and has not previously been
observed to spontaneously break symmetry in response to activity. One reason
for this apparent absence is simply because much of the theoretical,
computational and experimental work on active nematics is on 2D films.
However, the recent advent of experimentally realizable 3D active nematics
Duclos  _et al._ (2020) has spurred immediate interest in active liquid
crystals beyond films Shendruk  _et al._ (2018); Čopar  _et al._ (2019);
Chandrakar  _et al._ (2020); Chandragiri  _et al._ (2020); Binysh  _et al._
(2020); Ruske and Yeomans (2021); Varghese  _et al._ (2020); Houston and
Alexander (2021). By removing the dimensional limitations of confinement in
2D, additional classes of self-organized flows states may be possible. We
report long-lived, non-zero helical active dissipative 3D structures that
spontaneously break chiral symmetry and have not to our knowledge been
previously been observed.

We numerically confine active nematics in square ducts of size $L$ defined by
four impermeable, no-slip walls. The square-duct geometry is utilized to
produce quasi-1D confinement without curved boundaries Ellis  _et al._ (2018);
Pearce  _et al._ (2019); Maroudas-Sacks  _et al._ (2021). We vary the
extensile activity, reduced temperature, and anchoring strength on the channel
walls. When the anchoring is weak and reduced temperature is high, we find a
double helix flow structure, in which two streams of anti-parallel flow braid
around one another spontaneously breaking chiral symmetry (Fig. 1a; Movie 1).
When the anchoring strength is strong, we find a grinder train flow structure,
in which a lattice of axially aligned, counter-rotating vortices drift down
the channel with non-zero helicity (Fig. 1b; Movie 2). Both flow states are
helical, span the entire system and are long-lived. The self-organization of
helicity from a system with nematic symmetry highlights the potential of
activity as a pathway for designing emergent, material dynamics. Spontaneous
chiral symmetry breaking is characteristic of helix formation, which extends
beyond active fluids.

In our model of an active nematic fluid, the non-dimensionalized velocity
$\vec{u}(\vec{r};t)$ and nematic order $\boldsymbol{\mathbf{Q}}(\vec{r};t)$
are coupled, and density is constant Marenduzzo  _et al._ (2007); Thampi  _et
al._ (2014); Rivas  _et al._ (2020). The dimensionless Cauchy equation
balances the material derivative of velocity against the divergence of the
stress, $\tilde{\rho}D_{t}\vec{u}=\vec{\nabla}\cdot\boldsymbol{\mathbf{\Pi}}$.
Since the active contribution to the stress $-\zeta\boldsymbol{\mathbf{Q}}$
drives spontaneous nematic distortions on length scales $\ell_{\text{act}}$
and flows of speed $U$, we non-dimensionalize each term by the activity
$\zeta$, resulting in the active inertia number $\tilde{\rho}=\rho
U^{2}/\zeta$ for density $\rho$. Two additional dimensionless numbers appear
from the stress divergence — dimensionless viscosity (Eq. SI(5)) and, more
relevantly, the activity number $\tilde{K}^{-1/2}=L/\ell_{\text{act}}$, where
$\ell_{\text{act}}\sim\sqrt{K/\zeta}$ is the characteristic activity length
due to the competition of activity $\zeta$ and nematic elasticity $K$.

The evolution equation for the nematic orientation
$D_{t}\boldsymbol{\mathbf{Q}}-\boldsymbol{\mathbf{\mathcal{S}}}=\tilde{\Gamma}\boldsymbol{\mathbf{\mathcal{H}}}$,
in which the covariant derivative (including a co-rotational advection term
$\boldsymbol{\mathbf{\mathcal{S}}}$), is balanced by the relaxation towards
equilibrium, characterized by an inverse nematic Péclet number
$\tilde{\Gamma}$. The molecular field $\boldsymbol{\mathbf{\mathcal{H}}}$ has
contributions due to the bulk (characterized by the lowest order Landau free
energy coefficient $A$), distortions (by $K$) and surface anchoring (by
anchoring strength $W$). Non-dimensionalizing these by the distortion free
energy scale produces two dimensionless numbers — the distance from the
passive-thermodynamic isotropic-nematic transition, which is a reduced
temperature $\tilde{T}^{1/2}=\sqrt{A/K}L$; and the strength of the degenerate
planar anchoring, which can be measured as $\tilde{W}=WL/K$, as in colloidal
liquid crystals Aryasova  _et al._ (2004). Further technical details of the
model are provided in SI A.

Non-dimensionalized parameter space is denoted by

$\displaystyle\tilde{P}$
$\displaystyle=\left\\{\tilde{K}^{-1/2},\tilde{W},\tilde{T}^{1/2}\right\\}\equiv\left\\{\sqrt{\frac{\zeta}{K}}L,\frac{W}{K}L,\sqrt{\frac{A}{K}}L\right\\},$
(1)

each component of which can be interpreted as a ratio of system size to a
characteristic length scale (SI B). By varying anchoring strength $\tilde{W}$
and reduced temperature $\tilde{T}^{1/2}$, the simulations explore regions of
parameters space that have previously not received adequate attention.

For sufficiently small extensile activity ($\tilde{K}^{-1/2}\lesssim 18$), we
find spontaneous streaming in the axial direction. These are unidirectional
and oscillating flows, which have been well documented in 2D channels Giomi
_et al._ (2008); Fielding  _et al._ (2011); Shendruk  _et al._ (2017). Here,
we do not delineate between these but refer to both as axial streaming. At
sufficiently large activity, mesoscale or active turbulence occurs. In active
turbulence, the components of vorticity
$\vec{\omega}=\vec{\nabla}\times\vec{u}$ in the axial and two transverse
directions (denoted by subscript $\parallel$ and $\perp_{1,2}$ respectively)
are isotropic—the probability distributions of the axial $\omega_{\parallel}$
and transverse $\omega_{\perp_{1,2}}$ components of the vorticity are
equivalent (Fig. 2d). The distributions are symmetric about zero and normal.

The active dissipative structures exist in the intermediate regime between
axial streaming and active turbulence, where the active length scale
$\ell_{\text{act}}$ is comparable to the confinement length $L$. The most
common of these dissipative structures is a non-helical lattice of counter-
rotating vortices oriented transverse to the channel, breaking translational
symmetry (Fig. 1c). This is the 3D equivalent of the 2D Ceilidh dynamic state
Shendruk  _et al._ (2017); Doostmohammadi  _et al._ (2017). The 3D Ceilidh
lattice exhibits dancing disclinations; however, these are now curved
disclination lines that span the channel (Movie 4). As a vortex lattice, the
vorticity distribution is strongly bimodal (Fig. 2c). The $\omega_{\perp_{2}}$
distribution is symmetric about zero and bimodal with prominent peaks,
representing the lattice of counter-rotating vortices with spontaneous
symmetry breaking between the two transverse directions. There is negligible
vorticity in both the other directions.

While the grinder train is similar to the Ceilidh lattice, the crucial
distinction is that the counter-rotating vortices are oriented axially, rather
than transversely (Fig. 1b; Movie 2). Like axial streaming states Voituriez
_et al._ (2006); Chandragiri  _et al._ (2020), the grinder train has a net
flow along the channel. This flow structure is a train of axially aligned
counter-rotating vortices drifting down the axis of the channel that exists
when anchoring is strong (Movie 2). Thus, the distributions of vorticity are
similar to the Ceilidh vortex, except the bimodal distribution is narrower and
in the axial direction $\omega_{\parallel}$, rather than a transverse
direction (Fig. 2b). Though both the Ceilidh lattice and the grinder train
manifest counter-rotating vortices and both break translational symmetry, the
Ceilidh lattice does not possess local helicity $H=\vec{u}\cdot\vec{\omega}$
in contrast to the grinder train (Fig. 3a). The instantaneous centreline
helicity forms a well-defined wave (Fig. 3c) and the drift exhibits temporal
oscillation (Fig. 3d; Movie 5), reflecting the steady motion of the helicity
train. The steady translation, but oscillating sign, of the helicity indicates
the grinder is a lattice of counter-rotating vortices.

Figure 3: Local helicity $H=\vec{u}\cdot\vec{\omega}$ for the double helix
($\tilde{P}=\\{25,0.0015,6.5\\}$) and grinder train
($\tilde{P}=\\{25,250,0\\}$ flow states). (a) Snapshots showing the helicity
field as a 3D isosurface through the center plane of the duct. Top: Double
helix with negative chirality (Movie 6). Bottom: Grinder train with
alternating chirality (Movie 5). (b) Distributions of helicity down the
channel centreline averaged over ten random initializations. (c) Spatial
variation of instantaneous cross-sectional averaged helicity. (d) Temporal
variation of cross-sectional averaged helicity at a single point.

Qualitatively, the grinder train possesses characteristics of both axial
streaming and vorticity translational symmetry breaking. Indeed, as the axial
streaming begins to oscillate on scales comparable to the confinement,
counter-rotating vortices can align axially, resulting in non-zero local
helicity—though the average is zero (Movie 6), the distribution is bimodal
(Fig. 3b). It also generates smaller vorticity (Fig. 2b-c) due to the fact
that it possesses less nematic distortions and is entirely free of
disclinations (Movie 4).

Figure 4: Active dissipative structure phase diagrams. (a) Spatial velocity-
velocity correlations along the channel centre line for
$\tilde{P}=\\{\tilde{K}^{-1/2},250,0\\}$ for various $\tilde{K}^{-1/2}$.
Example curves for axial streaming (AS; $\tilde{K}^{-1/2}=15$), grinder train
(GT; $\tilde{K}^{-1/2}=20$), vortex lattice (VL; $\tilde{K}^{-1/2}=15$) and
mesoscale turbulence (T; $\tilde{K}^{-1/2}=15$). Open symbols denote non-zero
long-range correlations, while square symbols denote long-range oscillations
with non-zero amplitude. (b) Characteristic length scale $\ell_{u}$ determined
from the autocorrelation functions (SI C) for
$\tilde{P}=\\{\tilde{K}^{-1/2},250,0\\}$. $\boxplus$ marks large values of
peak resolution (Eq. SI(27)). (c) The length scales $\ell_{u}$ as a function
of dimensionless extensile activity $\tilde{K}^{-1/2}$ and planar anchoring
$\tilde{W}$ for reduced temperature $\tilde{T}^{1/2}=0$. (d) Length scales as
a function of $\tilde{W}$ and $\tilde{T}^{1/2}$ for $\tilde{K}^{-1/2}=25$.
Markers denote for average system helicity (SI E): $\boxtimes$ denotes large
average helicity $H_{\left|\langle\rangle\right|}$ and
$\operatorname{\leavevmode\hbox to7.15pt{\vbox
to7.15pt{\pgfpicture\makeatletter\hbox{\hskip 0.2pt\lower-0.2pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}
{}{{}}{}{}{}{}{{}}{}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{6.75pt}\pgfsys@lineto{6.75pt}{6.75pt}\pgfsys@lineto{6.75pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{6.75pt}{6.75pt}\pgfsys@stroke\pgfsys@invoke{
} {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{6.75pt}{6.75pt}\pgfsys@stroke\pgfsys@invoke{
} } \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}$ indicates large
values of $H_{\langle\left|\right|\rangle}$, which respectively identify
regions for the double helix (DH) and intermittent double helix (IDH).

On the other hand, the double helix does possess large nematic distortions: a
clear cork-screw disclination line spirals through the centre of the channel
(Movie 4). This is because the grinder train exists at strong anchoring, while
the double helix exists at weak anchoring and high reduced temperature, which
leads to more disorderly helicity fields (Fig. 3d; Movie 6). The helicity is
largest far from the no-slip boundaries. We find double helices most
frequently in regions of high reduced temperature $\tilde{T}^{1/2}$,
suggesting that the double helix emerges when the thermodynamic persistence
length is comparable to the confinement size. Unlike the Ceilidh or grinder
lattices, the double helix is not composed of isolated vortices and so the
centreline helicity does not exhibit clear spatio-temporal structure (Fig.
3c-d). Additionally, the vorticity of the double helix state is smallest in
the axial direction and not strongly bimodal in the two transverse directions
(Fig. 2a). Unlike the grinder train, the helicity distribution is primarily
unimodal (Fig. 3b). Since it is composed of a pair of helices with opposite
handedness but also opposite streaming directions, the double helix
spontaneously breaks chiral symmetry.

To map the existence of helical flow in parameter space, we characterize the
flow states by their associated length scales and helicity. The non-helical
Ceilidh vortex lattice is long-lived and system-spanning, with spatial
persistence clear from the velocity-velocity autocorrelation function $C_{uu}$
(SI C) down the length of the channel (Fig. 4a; blue squares), exhibiting
long-range oscillating correlations. Dissimilarly, the correlation function of
the turbulent state rapidly decorrelates (Fig. 4a; blue circles). While the
correlation functions of these two intercept zero (Fig. 4; filled symbols),
the axial streaming state initially decays to a long-range constant value
(Fig. 4a; open pink circles) and the grinder train possesses long-range
oscillating correlations about a non-zero value (Fig. 4a; open purple
squares). The correlation function of the grinder train has a longer period
but smaller amplitude than the Ceilidh lattice and does not decorrelate (Fig.
4a) since the grinder train state appears between axial streaming and the
vortex lattice (Fig. 4b).

The active dissipative structures are identified by their characteristic
length scales $\ell_{u}$, measured from the velocity correlations (SI C). The
length scales reveal distinct transitions between each flow state (Fig. 4b).
The Ceilidh lattice occupies the largest region of the activity-anchoring
phase plane (Fig. 4c). The phase diagram demonstrates that axial streaming
exists for all anchoring numbers at sufficiently small activity. The activity
at which the flow transitions from axial streaming is only weakly dependent on
anchoring strength across many orders of magnitude (Fig. 4c). For moderate to
weak anchoring, the Ceilidh lattice is the only dissipative structure between
axial streaming and active turbulence.

The grinder train arises because the anchoring suppresses the formation of the
dancing disclination lines that are present in the vortex lattice. This
reduction in nematic distortion lowers the vorticity magnitude and makes the
flow deterministic, such that the transverse directions coincide (Fig. 2b).
However, it does not alter the match between the characteristic vortex size
$\ell_{\text{act}}$ and the confinement $L$. Thus, the vortices preferentially
lie axially along the channel.

When the reduced temperature $\tilde{T}^{1/2}$ is increased above the
isotropic-nematic transition, the fluid can still be dynamically oriented by
the active flows Santhosh  _et al._ (2020); Giordano  _et al._ (2021). Indeed,
vortex lattices exist for a wide range of $\tilde{T}^{1/2}$ and $\tilde{W}$ at
$\tilde{K}^{-1/2}=25$ (Fig. 4d). However, in the low anchoring limit and at
moderate reduced temperatures, the double helix structure exists as a second
region for which there is non-zero helicity. While the grinder train and
double helix both exist for intermediate activity, the grinder train requires
strong anchoring at $\tilde{T}^{1/2}=0$ (Fig. 4c) and the double helix occurs
for weak anchoring at high reduced temperatures (Fig. 4d).

In addition to the orderly double helix and grinder train, we identify a less
structured, noisy regime at high reduced temperatures and strong anchoring
(Fig. 4d; slashes). In this regime, we find intermittent double helices
interspersed with more chaotic behaviour (Movie 7). Naturally, the chaos in
these systems results in inconsistent helicity but higher fluid speeds,
resulting in lower average helicity.

To understand the origin of active helical structures, we derive the transport
equation for helicity

$\displaystyle\tilde{\rho}D_{t}H+\vec{\nabla}\cdot\vec{J}$
$\displaystyle=\Sigma$ (2)

with helicity flux $\vec{J}$ and source terms $\Sigma$ (see SI D). By
calculating the dimensionless numbers associated with each source (Eqs.
SI(20)-(23)), we find that the nematic and viscous sources are negligible
compared to the activity-induced helicity source
$\Sigma^{\text{a}}=-\vec{\omega}\cdot[\vec{\nabla}\cdot\boldsymbol{\mathbf{Q}}]-\vec{u}\cdot(\vec{\nabla}\times[\vec{\nabla}\cdot\boldsymbol{\mathbf{Q}}])$.
The first term emerges from the projection of active force on to vorticity
field, while the second term represents projection of the curl of the active
force onto the velocity field. In both helical flow states, the two active
terms are of the same magnitude (Fig. SI10). Immediately following
initialization, the transient grinder train structure begins as oscillating
discoidal regions of $\Sigma^{\text{a}}$ (Movie 8). Once the drifting grinder
train forms, $\Sigma^{\text{a}}$ takes the form of a travelling wave that
produces the drifting train of alternating helicity (Fig. 3). Similarly for
the double helix, the active source forms local regions of production and
elimination (Movie 9), which are associated with the cork-screw disclination.

Recent work on 3D active nematics Shendruk  _et al._ (2018); Čopar  _et al._
(2019); Chandrakar  _et al._ (2020); Chandragiri  _et al._ (2020); Binysh  _et
al._ (2020); Ruske and Yeomans (2021); Varghese  _et al._ (2020); Houston and
Alexander (2021) has shown them to be complicated by their tangle of 3D
disclination lines. However, they also possess the potential for an exciting
range of possible dissipative structures. We have shown a simple case with the
emergence of two previously unobserved spontaneous, long-lived, non-zero
helicity structures in simple confinement for extensile activity, planar
anchoring and non-negative reduced temperature. In driven systems, double
helical flow states can exhibit complex non-linear dynamics Delbende  _et al._
(2021) and control the dispersivity of particles eluting in streams Song and
Brenner (2009). In fact, surface-activity driven helical flows have been
identified and analyzed in the context of cytoplasmic streaming Goldstein  _et
al._ (2008); van de Meent  _et al._ (2008); Verchot-Lubicz and Goldstein
(2010); Woodhouse and Goldstein (2013). In these cases of long plant cells,
helical flow enables significant transport, mixing and enhanced rates of
nutrient exchange with the surrounding membrane Goldstein  _et al._ (2008).
Indeed, recent work has considered defect dynamics Pearce (2020) and helical
flows Napoli and Turzi (2020) on cylindrical surfaces. Spontaneous chiral
symmetry breaking is a characteristic aspect of helix formation in such
material dynamics, which extends beyond the sphere of active fluids. For
example, chirality has been shown to be relevant in morpohological processes,
such as the gastrulation of Drosophila embryos, in which gut formation
acquires spontaneous left-right asymmetry that leads to twisting Hozumi  _et
al._ (2006). Likewise, the mitotic spindle, which has been modelled as an
active nematic Oriola  _et al._ (2020), exhibits chirality Novak  _et al._
(2018). Furthermore, spontaneous chiral symmetry breaking has been observed in
the director field of passive achiral nematics subject to strong confinements
Tortora and Lavrentovich (2011); Koning  _et al._ (2014); Nayani  _et al._
(2015); Jeong  _et al._ (2015).

The experimental realization of 3D active nematics Duclos  _et al._ (2020)
promises many opportunities to explore active dissipative structures that were
not possible in 2D and helicity may play a key in many such structures, as
suggested by these results. Helicity is involved in energy cascades in active
turbulence Słomka and Dunkel (2017) and generally plays an important role in
characterizing the topological nature of 3D flows Moffatt (2021). This letter
reports spontaneous chiral symmetry breaking in active nematics leading to a
steady state helical flow. Indeed, in stark contrast to the defect dancing
observed in the Ceilidh vortex lattice state, the grinder train does not
possess any topological singularities, while the double helix possesses a
single disclination line that winds through the two helical flows.

## Acknowledgements

We thank Kevin Stratford and Kristian Thijssen for invaluable help. This
research has received funding (TNS) from the European Research Council under
the European Union’s Horizon 2020 research and innovation programme (Grant
agreement No. 851196). SPT acknowledges the support by the Department of
Science and Technology, India, via the research grant CRG/2018/000644. TA-N
has been supported in part by the Academy of Finland PolyDyna (no. 307806) and
QFT Center of Excellence Program grant (no. 312298).

##  Appendix A Numerical Methods

Simulations numerically solve the incompressible active nematohydrodynamics
equations of motion Marenduzzo  _et al._ (2007); Thampi  _et al._ (2014) for
the local non-dimensionalized velocity $\vec{u}(\vec{r};t)$ and nematic order
$\boldsymbol{\mathbf{Q}}(\vec{r};t)$, for constant density,
$\vec{\nabla}\cdot\vec{u}=0$. The non-dimensionalized Cauchy equation balances
the material derivative $D_{t}=\partial_{t}+\vec{u}\cdot\vec{\nabla}$ of
velocity against the divergence of the stress $\boldsymbol{\mathbf{\Pi}}$

$\displaystyle\tilde{\rho}D_{t}\vec{u}$
$\displaystyle=\vec{\nabla}\cdot\boldsymbol{\mathbf{\Pi}}.$ (3)

Since activity drives the fluid flow, the stress is normalized by the
magnitude of the activity $\zeta$, resulting in a dimensionless active inertia
number

$\displaystyle\tilde{\rho}$ $\displaystyle=\frac{\rho U^{2}}{\zeta},$ (4)

which compares inertial to active forces in an active fluid of density $\rho$
and spontaneous characteristic speed $U$. The non-dimensionalized stresses in
Eq. (3) have four contributions Denniston  _et al._ (2004):

(i)

active stress $\boldsymbol{\mathbf{\Pi}}^{\text{a}}=-\boldsymbol{\mathbf{Q}}$;

(ii)

viscous stress
$\boldsymbol{\mathbf{\Pi}}^{\text{v}}=2\tilde{\eta}\boldsymbol{\mathbf{E}}$,
where
$2\boldsymbol{\mathbf{E}}=\vec{\nabla}\vec{u}+\left(\vec{\nabla}\vec{u}\right)^{\intercal}$
is the strain rate and

$\displaystyle\tilde{\eta}$ $\displaystyle=\frac{\eta U}{\zeta L}$ (5)

is the dimensionless viscosity number;

(iii)

nematic contribution
$\boldsymbol{\mathbf{\Pi}}^{\text{n}}=\tilde{K}\boldsymbol{\mathbf{\pi}}^{\text{n}}$,
where
$\boldsymbol{\mathbf{\pi}}^{\text{n}}=2\lambda(\boldsymbol{\mathbf{Q}}:\boldsymbol{\mathbf{\mathcal{H}}})\boldsymbol{\mathbf{\mathcal{Q}}}-\lambda\left\\{\boldsymbol{\mathbf{\mathcal{H}}},\boldsymbol{\mathbf{\mathcal{Q}}}\right\\}_{+}+\left\\{\boldsymbol{\mathbf{Q}},\boldsymbol{\mathbf{\mathcal{H}}}\right\\}_{-}-\vec{\nabla}\left(\boldsymbol{\mathbf{Q}}:\tfrac{\delta
F}{\delta\vec{\nabla}\boldsymbol{\mathbf{Q}}}\right)$ in which
$\boldsymbol{\mathbf{\mathcal{Q}}}=\boldsymbol{\mathbf{Q}}+\boldsymbol{\mathbf{\delta}}/3$,
$\boldsymbol{\mathbf{\mathcal{H}}}=-(\frac{\delta
F}{\delta\boldsymbol{\mathbf{Q}}}-\frac{1}{3}\boldsymbol{\mathbf{I}}\;\text{Tr}\frac{\delta
F}{\delta\boldsymbol{\mathbf{Q}}})$ is the molecular field, $\lambda$ is an
alignment parameter, and
$\left\\{\boldsymbol{\mathbf{A}},\boldsymbol{\mathbf{B}}\right\\}_{\pm}=\boldsymbol{\mathbf{A}}\cdot\boldsymbol{\mathbf{B}}\pm\boldsymbol{\mathbf{B}}\cdot\boldsymbol{\mathbf{A}}$
is a (anti)commutator; and

(iv)

hydrodstatic pressure
$\boldsymbol{\mathbf{\Pi}}^{\text{p}}=-P\boldsymbol{\mathbf{\delta}}$, where
pressure is non-dimensionalized by activity.

The dimensionless number

$\displaystyle\tilde{K}$ $\displaystyle=\frac{K}{\zeta L^{2}}$ (6)

assumes a single elastic Frank coefficient $K$ and is the squared ratio of the
characteristic activity length in bulk $\ell_{\text{act}}\sim\sqrt{K/\zeta}$
and the scale of confinement $L$. A convenient dimensionless activity number
is then

$\displaystyle\tilde{K}^{-1/2}$ $\displaystyle=\frac{L}{\ell_{\text{act}}},$
(7)

which is employed in the main text as the dimensionless measure of activity.

The transport equation for the nematic order tensor balances the relaxation
towards equilibrium against the generalized material derivative

$\displaystyle D_{t}\boldsymbol{\mathbf{Q}}-\boldsymbol{\mathbf{\mathcal{S}}}$
$\displaystyle=\tilde{\Gamma}\boldsymbol{\mathbf{\mathcal{H}}},$ (8)

including co-rotational advection,

$\displaystyle\boldsymbol{\mathbf{\mathcal{S}}}$
$\displaystyle=\left\\{\lambda\boldsymbol{\mathbf{E}}+\boldsymbol{\mathbf{\mathcal{W}}},\boldsymbol{\mathbf{\mathcal{Q}}}\right\\}_{+}-2\lambda\left(\boldsymbol{\mathbf{Q}}:\vec{\nabla}\vec{u}\right)\boldsymbol{\mathbf{\mathcal{Q}}},$
(9)

where
$2\boldsymbol{\mathbf{\mathcal{W}}}=\left(\vec{\nabla}\vec{u}\right)^{\intercal}-\vec{\nabla}\vec{u}$
is the vorticity tensor. An inverse nematic Péclet number

$\displaystyle\tilde{\Gamma}$ $\displaystyle=\frac{\Gamma K}{LU}$ (10)

relates elastic rotational relaxation rate $\Gamma K/L^{2}$ to advection rate
$L/U$.

Figure 5:  Average root mean squared velocities $U$ in LB simulation units for
reduced temperature $\tilde{T}^{1/2}=0$.

The non-dimensionalized molecular field
$\boldsymbol{\mathbf{\mathcal{H}}}=-\tfrac{\delta
F}{\delta\boldsymbol{\mathbf{Q}}}+\tfrac{\boldsymbol{\mathbf{I}}}{3}\text{Tr}\left[\tfrac{\delta
F}{\delta\boldsymbol{\mathbf{Q}}}\right]$ appears in both Eq. (3) and (8) and
belongs to the free energy $F$, which is comprised of Landau-de Gennes and
deformation contributions $F=\tilde{A}_{0}\int d^{3}r\
\tfrac{\left(1-\gamma/3\right)}{2}\boldsymbol{\mathbf{Q}}:\boldsymbol{\mathbf{Q}}-\tfrac{\gamma}{3}\boldsymbol{\mathbf{Q}}:(\boldsymbol{\mathbf{Q}}\cdot\boldsymbol{\mathbf{Q}})+\tfrac{\gamma}{4}\left(\boldsymbol{\mathbf{Q}}:\boldsymbol{\mathbf{Q}}\right)^{2}+\left(\vec{\nabla}\boldsymbol{\mathbf{Q}}\right)^{2}$.
The bulk free energy density has been normalized by the elasticity
$\tilde{A}_{0}=A_{0}L^{2}/K$ and includes the dimensionless parameter
$\gamma$, which determines the magnitude of the nematic order with
$\gamma=2.7$ the coexistence point between isotropic and nematic phases
Marenduzzo  _et al._ (2007). Within the nematic phase, the characteristic
defect core size can be approximated by the thermodynamic nematic persistence
length $\ell_{Q}\sim\sqrt{K/A}$ for $A=A_{0}\left(1-\gamma/3\right)$ Hemingway
_et al._ (2016), such that we can define a reduced temperature in terms of a
squared ratio of length scales

$\displaystyle\tilde{T}$
$\displaystyle=\left(\frac{L}{\ell_{Q}}\right)^{2}=\tilde{A}_{0}\left(\frac{3-\gamma}{3}\right).$
(11)

This is convenient, as it allows us to directly compare the effects of
thermodynamic nematic ordering to confinement. We utilize $\tilde{T}^{1/2}$ as
a measure of reduced temperature in the main text.

In addition to the bulk, the walls impose degenerate planar anchoring Fournier
and Galatola (2005) via a surface free energy density

$\displaystyle f_{\text{wall}}$
$\displaystyle=W_{1}\left(\boldsymbol{\mathbf{Q}}-\boldsymbol{\mathbf{Q}}^{\perp}\right):\left(\boldsymbol{\mathbf{Q}}-\boldsymbol{\mathbf{Q}}^{\perp}\right)$
$\displaystyle\qquad+W_{2}\left(\boldsymbol{\mathbf{Q}}:\boldsymbol{\mathbf{Q}}-S_{\text{eq}}^{2}\right)^{2},$
(12)

where $6S_{\text{eq}}=1+3\sqrt{1-8/(3\gamma)}$ is the equilibrium value of the
scalar order parameter. The projection of the orientation tensor onto the
surface is
$\boldsymbol{\mathbf{Q}}^{\perp}=\boldsymbol{\mathbf{P}}\cdot\tilde{\boldsymbol{\mathbf{Q}}}\cdot\boldsymbol{\mathbf{P}}$,
where
$\tilde{\boldsymbol{\mathbf{Q}}}=\boldsymbol{\mathbf{Q}}+S_{\text{eq}}\boldsymbol{\mathbf{\delta}}/3$,
the projection operator is
$\boldsymbol{\mathbf{P}}=\boldsymbol{\mathbf{\delta}}-\vec{n}_{\perp}\vec{n}_{\perp}$
and the surface normal is $\vec{n}_{\perp}$. We set both the surface anchoring
constants for the orientation and ordering to the same value $W=W_{1}=W_{2}$
and non-dimensionalize the anchoring

$\displaystyle\tilde{W}$ $\displaystyle=\frac{WL}{K}.$ (13)

As with activity number $\tilde{K}^{-1/2}$ and the reduced temperature
$\tilde{T}^{1/2}$, $\tilde{W}$ can be interpreted as the ratio of confinement
size $L$ to the de Gennes-Kleman extrapolation length $\ell_{\text{dGK}}\sim
K/W$ Suh  _et al._ (2019). In the main text, $\tilde{W}$ serves as the measure
of anchoring strength.

Figure 6:  Non-dimensional viscosity number, $\tilde{\eta}=(\eta U)(\zeta
L)\sim 10^{-2}$, for reduced temperature $\tilde{T}^{1/2}=0$.

The equations of motion are solved using a hybrid scheme Marenduzzo  _et al._
(2007), in which the Cauchy equation (Eq. (3)) is solved using the lattice
Boltzmann algorithm, while Eq. (8) is solved using a finite difference
predictor-corrector algorithm. Simulations are performed in 3D square channels
of size $L=25$ defined by impermeable, no-slip walls and periodic boundary
conditions separated by a channel length of $\mathcal{L}\in\\{128,200\\}$.
Simulations are run for a minimum of $2\times 10^{5}$ time steps to a maximum
of $5\times 10^{5}$. The density $\rho=1$, rotational diffusivity
$\Gamma=0.3375$ and alignment parameter $\lambda=1$ are all held constant. On
the other hand, Landau-de Gennes energy density scale
$A_{0}\in\left\\{0.1,1\right\\}$, distance from the isotropic-nematic
transition $\gamma\in[2.668,3]$, dynamic viscosity $\eta=4/3$, Frank
elasticity $K\in\left\\{0.05,0.05555\right\\}$, anchoring energy
$W\in\left[3\times 10^{-6},0.5\right]$ and activity
$\zeta\in\left[0.018,0.122\right]$ are varied throughout.

##  Appendix B Magnitude of non-dimensional numbers

Not all the terms in the active nematohydrodynamic equations of motion are
significant. The parameters produce flows with average root mean squared
velocities $U\sim 10^{-2}$ (Fig. 5). The lowest values ($U\approx 3.4\times
10^{-3}$) occur in the vortex lattice phase (intermediate $\tilde{K}^{-1/2}$
for low anchoring $\tilde{W}$), as well as axial streaming for intermediate
anchoring $\tilde{W}$. The largest ($U\approx 5.5\times 10^{-2}$) occur for
high activities and low anchoring — in which the fluid streams at high
velocities with non-zero long-range correlations. Throughout the main text, we
take $U$ to be the characteristic velocity scale and employ it to normalize
the quantities of interest.

Figure 7:  Non-dimensional inertia number, $\tilde{\rho}=\rho U^{2}/\zeta\sim
10^{-2}$, for reduced temperature $\tilde{T}^{1/2}=0$.

The small values of $U$ cause many of the non-dimensionless numbers from SI A
to be small. In particular, the dimensionless viscosity number
$\tilde{\eta}=(\eta U)/(\zeta L)$ (Fig. 6) can be interpreted as a ratio of
the measured spontaneous flow speed in the channel in units of the
characteristic active speed $\zeta L/\eta$. Because of this, the map of
$\tilde{\eta}$ is identical to the map for $U$ (Fig. 5 and Fig. 6). The fact
that $\tilde{\eta}\sim 10^{-2}$ highlights the difference between confined
active nematics and unconfined, bulk active turbulence: In bulk, one expects
that the spontaneous flow directly results from the balance viscous forces to
active force and so $\tilde{\eta}\sim 1$. However, the no-slip boundaries
modify this expectation and cause $\tilde{\eta}$ to be substantially smaller.

The dimensionless inertia number $\tilde{\rho}=\rho U^{2}/\zeta$ (Fig. 7) is
qualitatively similar to the maps of $\tilde{\eta}$ and $U$. Indeed, the order
of magnitude of $\tilde{\rho}\sim 10^{-2}$ is the same as the dimensionless
viscosity $\tilde{\eta}$. Because $\tilde{\rho}\sim U^{2}$, the dimensionless
inertia number is much larger in the region with the fastest flows (high
activities and low anchoring) than elsewhere. The small scale of
$\tilde{\rho}$ indicates that hydrodynamic flows that result from Eq. (3) is
primarily a balance of stresses
$\vec{\nabla}\cdot\boldsymbol{\mathbf{\Pi}}\approx 0$.

The inverse Peclet number $\tilde{\Gamma}=(\Gamma K)/(LU)$ is also relatively
small, with values $\tilde{\Gamma}\sim 10^{-1}$ (Fig. 8), despite
$\tilde{\Gamma}\sim U^{-1}$. This is because $\gamma$ and $K$ are both small
and $L$ is large. Because $\tilde{\Gamma}<1$, the generalized material
derivative in Eq. (8) is more significant than the thermodynamic relaxation
towards equilibrium, which indicates
$D_{t}\boldsymbol{\mathbf{Q}}\approx\boldsymbol{\mathbf{\mathcal{S}}}$. Thus,
the nematic orientational order field are primarily governed by the active
flows.

Thus, the other three dimensionless numbers from SI A act as our dimensionless
parameter space $\tilde{P}=\\{\tilde{K}^{-1/2},\tilde{W},\tilde{T}^{1/2}\\}$.
Each of these numbers can be interpreted as a ratio of length scales. We
consider dimensionless activity numbers $\tilde{K}^{-1/2}=L/\ell_{\text{act}}$
in the range $10$ to $40$, which is the regime in which flows transition from
axial streaming to mesoscale turbulence. The anchoring number
$\tilde{W}=L/(K/W)$ varies over five orders of magnitude, from $10^{-3}$ to
$250$. The reduced temperature $\tilde{T}^{1/2}=L/\ell_{Q}$ is varied from $0$
(representing the thermodynamic isotropic-nematic transition point) to $10$
(well within the isotropic phase).

While the Reynolds and Ericksen numbers do not explicitly appear, they are
both intermediate in magnitude (Fig. 9). The Ericksen number
$\text{Er}=\tilde{\eta}/\tilde{K}$ rises as high as $35$. However, this is
restricted to the high activity/low anchoring region of active streaming. In
the regions of interest (vorticity lattices, grinder train and double helix),
$\text{Er}\sim 3$. This suggests that active flows are moderately more
significant than elastic response. Likewise, the Reynolds number
$\text{Re}=\tilde{\rho}/\tilde{\eta}$ can approach unity in the high activity
region but is $\text{Re}\sim 0.1$ in the regions of interest. This indicates
that viscous forces are somewhat more prominent than inertial ones. However,
as described above, $\tilde{\rho}$ is the ratio that expresses the relative
importance of inertia to the activity that generates spontaneous flows; rather
than to viscosity. Thus, it is the smallness of $\tilde{\rho}$ that indicates
the relevance of inertial non-linearities.

Figure 8:  Non-dimensional inverse Péclet number, $\tilde{\Gamma}=(\eta
U)/(\zeta L)\sim 10^{-2}$, for reduced temperature $\tilde{T}^{1/2}=0$.

##  Appendix C Characteristic velocity lengths

Figure 9:  Ericksen numbers $\text{Er}=\tilde{\eta}/\tilde{K}$ (left) and
Reynolds numbers $\text{Re}=\tilde{\rho}/\tilde{\eta}$ (right) for reduced
temperature $\tilde{T}^{1/2}=0$.

To characterize the various flow states, we examine the spatial velocity-
velocity correlation function

$\displaystyle C_{uu}(x)$
$\displaystyle=\frac{\left\langle\vec{u}\left(\vec{r}_{0};t\right)\cdot\vec{u}\left(\vec{r}_{0}+x\hat{\vec{e}}_{\parallel};t\right)\right\rangle}{\left\langle\vec{u}\left(\vec{r}_{0};t\right)\cdot\vec{u}\left(\vec{r}_{0};t\right)\right\rangle}$
(14)

down the length of the channel for positions $\vec{r}_{0}$ within $\delta
r=\pm 2$ of the centre line of the channel. The autocorrelations $C_{uu}$
reveal four discernible families of behaviours: (i) Those with an initial
exponential decay but a constant, positive, long-range correlation; (ii) those
that pass through $C_{uu}=0$; (iii) those with an initial decay but long-range
oscillating correlations; and (iv) those without long-range oscillating
correlations.

To describe all four behaviours, the autocorrelation functions are fit by the
empirical form

$\displaystyle C_{uu}(x)$
$\displaystyle=\left(c_{1}e^{-x/\ell_{1}}+c_{2}\right)\cos\left(\frac{2\pi
x}{\ell_{2}}+\phi\right)+c_{0}.$ (15)

The autocorrelation function is a decaying cosine, which so possesses two
length scales: a decay length $\ell_{1}$ and a wavelength $\ell_{2}$. The
amplitude of the oscillations decays from $c_{1}+c_{2}$ to $c_{2}$, and the
autocorrelation oscillates about $c_{0}$. When the exponential decay does not
approach zero, the flow structure has long-range persistence, which is
identified in cases for which $c_{0}$ is significantly greater than zero. This
absence of any anti-correlation is indicative of a flow state streaming along
the channel, which is the case for the axial streaming and grinder states. In
Fig. 4, cases for which $c_{0}>1/3$ are represented by open symbols. For high
activity flow states (Ceilidh lattice and turbulence), $C_{uu}$ crosses zero
and $c_{0}<1/3$, which is denoted by filled symbols.

In the main text (Fig. 4), the spontaneous flow structures are quantified
according to the characteristic length scale $\ell_{u}$, which is determined
by assessing the fit of the correlation functions to Eq. (15). Whether the
characteristic length scale $\ell_{u}$ corresponds to the decay length
$\ell_{1}$ or the wavelength $\ell_{2}$, is depends on the value of $c_{2}$.
We found that all velocity autocorrelations had a short-range drop, which is
successfully fit using a decaying cosine of amplitude $c_{1}$ and decay length
$\ell_{1}$. In cases of non-periodic flow (axial flow and turbulence), there
is no long-range amplitude and so $c_{2}\approx 0$. In these cases, the
wavelength is unimportant; thus, the reported length scale is the decay length
$\ell_{u}=\ell_{1}$ and a circle symbol is used. However, both the Ceilidh
lattice and grinder train exhibited long-range oscillations in $C_{uu}$ of
long-range amplitude $c_{2}$. When $c_{2}>0$ the flow states are spatially
ordered with translational symmetry breaking (the Ceilidh lattice and grinder
train). These cases are shown by square symbols in Fig. 4. When $c_{2}>0.05$,
the spatial ordering with periodicity characterized by the wavelength
$\ell_{2}$. In these cases, the value characteristic lengthscale of the flow
is taken to be $\ell_{u}=\ell_{2}$.

##  Appendix D Helicity transport equations

In order to understand the source of helicity in double helix and grinder
train helical structures, we derive the transport equation for helicity for
active nematics. The local helicity $H=\vec{u}\cdot\vec{\omega}$ is a
pseudoscalar field and the integral of $H$ over the entire system is the total
helicity, such that $H$ can be referred to as the helicity density field
Moffatt and Tsinober (1992). As a first step towards the helicity transport
equation, we derive the transport equation for the vorticity. Starting from
Eq. (3), we note that the advection term can be rewritten as
$(\vec{u}\cdot\vec{\nabla})\vec{u}=\frac{1}{2}\vec{\nabla}u^{2}-\vec{u}\times\vec{\omega}$,
and take the curl of both sides of the equation. We next employ the fact that
$\vec{\nabla}\times(\vec{u}\times\vec{\omega})=\vec{\omega}\cdot\vec{\nabla}\vec{u}-\vec{u}\cdot\vec{\nabla}\vec{\omega}$
and that the curl of a gradient vanishes, we obtain

$\displaystyle\tilde{\rho}D_{t}\vec{\omega}$
$\displaystyle=\tilde{\eta}\vec{\nabla}^{2}\vec{\omega}+\tilde{\rho}(\vec{\omega}\cdot\vec{\nabla})\vec{u}$
$\displaystyle\qquad+\vec{\nabla}\times\left[\vec{\nabla}\cdot(\tilde{K}\boldsymbol{\mathbf{\pi}}^{n}-\boldsymbol{\mathbf{Q}})\right].$
(16)

This is the non-dimensionalized transport equation for vorticity in an active
nematic fluid. On the right hand side, the first two terms exist for passive
isotropic fluids fluids, whereas the last term is the curl of the divergence
of the passive nematic stresses Chono  _et al._ (1998) and activity.

We can employ the vorticity transport equation to derive the transport
equation for the helicity $H=\vec{u}\cdot\vec{\omega}$. We consider the rate
of change of helicity

$\partial_{t}H=(\vec{\omega}\cdot\partial_{t}\vec{u}+\vec{u}\cdot\partial_{t}\vec{\omega}),$
(17)

and use Eqs. (3) and (D) in Eq. (17). We notice four facts:

(i)

$\vec{\omega}\cdot(\vec{u}\cdot\vec{\nabla})\vec{u}+\vec{u}\cdot(\vec{u}\cdot\vec{\nabla})\vec{\omega}=(\vec{u}\cdot\vec{\nabla})H$;

(ii)

$\vec{\omega}\cdot\vec{\nabla}^{2}\vec{u}+\vec{u}\cdot\vec{\nabla}^{2}\vec{\omega}=\vec{\nabla}^{2}H-2[\vec{\nabla}\vec{\omega}:\vec{\nabla}\vec{u}]$;

(iii)

$\vec{\nabla}\cdot\vec{\omega}=0$; and

(iv)

$\vec{u}\cdot(\vec{\omega}\cdot\vec{\nabla})\vec{u}=\vec{\nabla}\cdot\left(\frac{1}{2}u^{2}\vec{\omega}\right)$.

Employing these facts, we obtain the dimensionless transport equation

$\displaystyle\tilde{\rho}D_{t}H+\vec{\nabla}\cdot\vec{J}$
$\displaystyle=\Sigma,$ (18)

where we have explicitly written the transport equation as resulting from the
divergence of flux $\vec{J}$ and source terms $\Sigma$. The flux is composed
of two terms

$\displaystyle\vec{J}$
$\displaystyle=-\tilde{\eta}\vec{\nabla}H+P^{\prime}\vec{\omega},$ (19)

where $P^{\prime}=P-\frac{1}{2}\tilde{\rho}\,u^{2}$ is an effective pressure.
There are three sources of helicity

$\displaystyle\Sigma$
$\displaystyle=\tilde{\eta}\Sigma^{\text{v}}+\tilde{K}\Sigma^{\text{n}}+\Sigma^{\text{a}}.$
(20)

The first contribution is a viscous source

$\displaystyle\Sigma^{\text{v}}$
$\displaystyle=-2\vec{\nabla}\vec{u}:\vec{\nabla}\vec{\omega},$ (21)

which is present in passive isotropic fluids. However, the second passive
contribution to the helicity is nematohydrodynamic in nature, with the form

$\displaystyle\Sigma^{\text{n}}$
$\displaystyle=\vec{u}\cdot\left(\vec{\nabla}\times\left[\vec{\nabla}\cdot\boldsymbol{\mathbf{\pi}}^{\text{n}}\right]\right)+\vec{\omega}\cdot\left[\vec{\nabla}\cdot\boldsymbol{\mathbf{\pi}}^{\text{n}}\right].$
(22)

The Beris-Edwards model of passive nematohydrodynamics is not expected to
conserve helicity Gay-Balmaz and Tronci (2011). Because $\tilde{K}\sim
10^{-3}$ and $\tilde{\eta}\sim 10^{-2}$ (Fig. 6), both the nematic and viscous
sources are insignificant. Most relevantly, the active contribution to the
helicity takes the form

$\displaystyle\Sigma^{\text{a}}$
$\displaystyle=-\vec{u}\cdot\left(\vec{\nabla}\times\left[\vec{\nabla}\cdot\boldsymbol{\mathbf{Q}}\right]\right)-\vec{\omega}\cdot\left[\vec{\nabla}\cdot\boldsymbol{\mathbf{Q}}\right]$
(23) $\displaystyle=\Sigma^{\text{a}}_{1}+\Sigma^{\text{a}}_{2}.$

The first of these two terms
$\Sigma^{\text{a}}_{1}=-\vec{u}\cdot(\vec{\nabla}\times[\vec{\nabla}\cdot\boldsymbol{\mathbf{Q}}])$
is the projection of the curl of the active force density
($-\zeta\vec{\nabla}\cdot\boldsymbol{\mathbf{Q}}$) on the velocity. On the
other hand, the second
$\Sigma^{\text{a}}_{2}=-\vec{\omega}\cdot[\vec{\nabla}\cdot\boldsymbol{\mathbf{Q}}]$
is the curl of velocity (i.e. the vorticity
$\vec{\omega}=\vec{\nabla}\times\vec{u}$) projected on the force density.

To test the expectation that the non-active contributions to $\Sigma$ are
negligible, we measure the viscous source field for the grinder train and
compare it to the active terms (Fig. 10a-c). We find that the viscous
contribution is two orders of magnitude smaller than the active contribution,
which corresponds directly to the scale of $\tilde{\eta}$. The viscous source
$\tilde{\eta}\Sigma^{\text{n}}\sim 10^{0}$ follows the form of the grinder
train helicity quite closely (Fig. 10a). The two active sources are seen to be
of the same order of magnitude, but the
$\Sigma^{\text{a}}_{1}=-\vec{u}\cdot(\vec{\nabla}\times[\vec{\nabla}\cdot\boldsymbol{\mathbf{Q}}])$
contribution is roughly twice as strong as the
$\Sigma^{\text{a}}_{2}=-\vec{\omega}\cdot[\vec{\nabla}\cdot\boldsymbol{\mathbf{Q}}]$
contribution (Fig. 10b-c). This is found to be true of both the grinder train
and double helix.

When the active sources are summed for the grinder train, they are seen to
produce a travelling wave of helicity production and elimination. However, the
net active source of helicity is small and oscillates about zero, explaining
why the net helicity in the system is zero. On the other hand for the double
helix, a net value is found when the active source term is averaged. The sign
of the average source is found to correspond to the sign of the helicity of
each double helix (Fig. 10d)

Figure 10:  Most significant non-dimensionalized sources of helicity $\Sigma$
for $\tilde{P}=\\{20,250,0\\}$. (a) Viscous source given by Eq. (21). (b) The
$\Sigma^{\text{a}}_{1}=-\vec{u}\cdot(\vec{\nabla}\times[\vec{\nabla}\cdot\boldsymbol{\mathbf{Q}}])$
contribution to activity-induced sources (first term from Eq. (23)). (c) The
$\Sigma^{\text{a}}_{2}=-\vec{\omega}\cdot[\vec{\nabla}\cdot\boldsymbol{\mathbf{Q}}]$
contribution to activity-induced sources (second term from Eq. (23)), the
magnitude of which is roughly half the first term. (d) System averaged active
helicity source $\left\langle\Sigma^{\text{a}}\right\rangle_{\mathbf{r}}(t)$
and averaged helicity $\left\langle HL/U^{2}\right\rangle_{\mathbf{r}}(t)$ as
a function of time for three different double helix flow states. The red
system occurs at $\tilde{P}=\\{25,1.5,9.47\\}$, blue at
$\tilde{P}=\\{25,0.15,4.7\\}$, and cyan at $\tilde{P}=\\{25,0.15,6.5\\}$.

##  Appendix E Helicity averages

In order to differentiate between helical flow states, multiple statistical
measures of the helicity field
$H\left(\vec{r};t\right)=\vec{u}\cdot\vec{\omega}$ are considered. The
instantaneous spatial average $\langle H\rangle_{\mathbf{r}}(t)$ is computed
by taking the mean of helicity over space at time $t$. To aid comparison in
the active dissipative structure phase diagram (Fig. 4c-d), spatially averaged
helicty is normalized by $\langle U\rangle^{2}_{\mathbf{r}}/L$, where $\langle
U\rangle_{\mathbf{r}}(t)$ is the instantaneous mean fluid speed at time $t$
and $L$ is the confinement length scale. Following this normalisation, the
absolute value of the non-dimensional temporal average is taken

$\displaystyle\langle\tilde{H}\rangle_{|\langle\rangle|}$
$\displaystyle=\left\langle\frac{\left|\langle
H\rangle_{\mathbf{r}}\right|}{\langle
U\rangle^{2}_{\mathbf{r}}/L}\right\rangle_{t},$ (24)

where $\langle\ldots\rangle_{t}$ denotes the temporal average. In Eq. (24),
the spatial average $\langle H\rangle_{\mathbf{r}}(t)$ is taken for the entire
channel and the temporal average of the absolute value of the spatial average
is evaluated. This is appropriate for the double helix flow states, which
spontaneously has positive or negative helicity. However, in the case of the
grinder train, the helicity oscillates symmetrically about zero in space (see
Fig. 3c) and so $\langle\tilde{H}\rangle_{|\langle\rangle|}\simeq 0$. Thus,
$\langle\tilde{H}\rangle_{|\langle\rangle|}>0$ highlights the double helix in
Fig. 11. As a cutoff value to obtain the $\otimes$ markers in Fig. 4d, we use
$\langle\tilde{H}\rangle_{|\langle\rangle|}>4$.

Because it is the temporal average of the absolute value of the spatial
average, $\langle\tilde{H}\rangle_{|\langle\rangle|}=0$ for the grinder state.
Thus, we also consider a separate measure of helicity by taking the modulus
before spatially averaging

$\displaystyle\langle\tilde{H}\rangle_{\langle||\rangle}$
$\displaystyle=\bigg{\langle}\frac{\langle\left|H\right|\rangle_{\mathbf{r}}L^{2}}{\langle
U\rangle_{\mathbf{r}}}\bigg{\rangle}_{t},$ (25)

where $\langle H\rangle_{\mathbf{r}}(t)$ has been taken over the channel
centreline. Both the double helicies and grinder trains exhibit non-zero
$\langle\tilde{H}\rangle_{\langle||\rangle}$.

Figure 11:  Non-dimensionalized, temporal average of the absolute value of the
spatially averaged helicity $\langle\tilde{H}\rangle_{|\langle\rangle|}$ given
by Eq. (24), for activity number $\tilde{K}^{-1/2}=25$. All values are close
to zero, except instances of the double helix dissipative structure.

The grinder train points can be seen in Fig. 12 as the points at high
anchoring and moderate activity $\tilde{K}^{-1/2}\approx 20$ with
$\langle\tilde{H}\rangle_{|\langle\rangle|}>4$. The average helicity
$\langle\tilde{H}\rangle_{\langle||\rangle}$ is small in the surrounding
parameter space. However, in addition to the grinder train, there exists an
extended region of large $\langle\tilde{H}\rangle_{\langle||\rangle}$ at high
activity and high anchoring number — this is the region of active turbulence.
The turbulent regime exhibits large values of
$\langle\tilde{H}\rangle_{\langle||\rangle}$ owing to high fluid speeds. To
differentiate between the grinder train and turbulent flow states, we would
thus need to combine the helicity measure with the periodic nature of the
velocity autocorrelation functions (Eq. (14)). However, we find that peak
resolution (SI F; Eq. (27)) is a clearer metric.

For the double helix both $\langle\tilde{H}\rangle_{|\langle\rangle|}$ (Fig.
11) and $\langle\tilde{H}\rangle_{\langle||\rangle}$ (Fig. 13) are large.
However, the intermittent double helix, which exists for high reduced
temperatures and anchoring numbers, has near zero
$\langle\tilde{H}\rangle_{|\langle\rangle|}$. As cutoffs for the intermittent
double helix, marked by $\oslash$ in Fig. 4d, we check that
$\langle\tilde{H}\rangle_{|\langle\rangle|}<1.0$ but
$\langle\tilde{H}\rangle_{\langle||\rangle}>4.5$, which is a value that is
even greater than in the orderly double helix.

Figure 12:  Non-dimensionalized, average of the absolute value of helicity
$\langle\tilde{H}\rangle_{\langle||\rangle}$ given by Eq. (25), for reduced
temperature $\tilde{T}^{1/2}=0$. In addition to non-zero values in regions of
active turbulence, the grinder train regime exists as well defined region of
large $\langle\tilde{H}\rangle_{\langle||\rangle}$ surrounded by small values
at strong anchoring and moderate activity.

##  Appendix F Helicity Peak Resolution

The distributions of helicity across the examined parameter space are
generally Gaussians, symmetric about zero. However, this is not true of the
helical regimes presented in Fig. 3b of the main text. The double helix does
have a Gaussian distribution of helicity but with non-zero mean. Furthermore,
the grinder train possesses a strongly peaked, symmetric, bimodal distribution
of helicity, which can be fit to the sum of two skewed normal distributions.
To each of the peaks in the distribution we find a fit of the form

$\displaystyle\mathcal{G}_{i}$
$\displaystyle=A_{i}\exp\left(-\frac{(x-\mu_{i})^{2}}{2\sigma_{i}^{2}}\right)\left[1+\textrm{erf}\left(\frac{S^{k}_{i}(x-\mu_{i})}{\sigma_{i}\sqrt{2}}\right)\right].$
(26)

The peaks are defined by their amplitude $A_{i}$, mean $\mu_{i}$, standard
deviation $\sigma_{i}$, and skew $S^{k}_{i}$, $i=1,2$. To quantify if the two
peaks are separated, and thus the distribution bimodal, we define the peak
resolution as

$\displaystyle\Delta$
$\displaystyle=\frac{\left|\mu_{1}-\mu_{2}\right|}{2\sqrt{\sigma_{1}\sigma_{2}}}.$
(27)

It is found that only the grinder train flow states possess strong bimodal
helicity (Fig. 14). We mark the grinder train with $\oplus$ in Fig. 4c when
$\Delta>1.2$.

Figure 13:  Non-dimensionalized, average of the absolute value of helicity
$\langle\tilde{H}\rangle_{\langle||\rangle}$ given by Eq. (25), for
$\tilde{K}^{-1/2}=25$. Non-zero values exist for the double helix, but also
active turbulent regions of parameter space.  Figure 14:  Peak resolution
$\Delta$ of the helicity distributions given by Eq. (27), for reduced
temperature $\tilde{T}^{1/2}=0$. Only the grinder train regime possesses
significant bimodality to produce large peak resolution.

##  Appendix G Movie Captions

  1. 1.

Double helix flow state for parameters
$\tilde{P}=\\{\tilde{K}^{-1/2},\tilde{W},\tilde{T}^{1/2}\\}=\\{25,1.5,9.47\\}$.
Two streams of anti-parallel flow braid around one another. Coloured by axial
velocity $u_{\parallel}$. Initialized at rest and isotropic.

  2. 2.

Grinder train for
$\tilde{P}=\\{\tilde{K}^{-1/2},\tilde{W},\tilde{T}^{1/2}\\}=\\{20,250,0\\}$. A
Procession of spontaneously drifting, counter-rotating vortices aligned
axially down the duct. Coloured by axial vorticity $\omega_{\parallel}$.
Initialized at rest and isotropic.

  3. 3.

Ceilidh vortex for
$\tilde{P}=\\{\tilde{K}^{-1/2},\tilde{W},\tilde{T}^{1/2}\\}=\\{22,250,0\\}$.
Vortex lattice oriented in the transverse direction do not exhibit net axial
flow. Coloured by transverse vorticity $\omega_{\perp}$. Initialized at rest
and isotropic.

  4. 4.

Disclinations for the flow states shown in Fig. 1 and Movie 1-3). Colourmaps
show the same quantities as in Fig. 1 in the central plane. Top: The double
helix flow state is accompanied by a single helical disclination line that
turns like a corkscrew. Middle: The grinder train is entirely devoid of
defects. Bottom: The vortex lattice possesses Ceilidh dancing defects which in
3D are arced disclination lines, as well as stationary disclination lines at
the walls. Time is stated in units of $10^{2}$ LB time steps.

  5. 5.

Local helicty density field for the grinder train for
$\tilde{P}=\\{\tilde{K}^{-1/2},\tilde{W},\tilde{T}^{1/2}\\}=\\{20,250,0\\}$.
Initialized at rest and isotropic.

  6. 6.

Local helicty density field for the double helix flow state for parameters
$\tilde{P}=\\{\tilde{K}^{-1/2},\tilde{W},\tilde{T}^{1/2}\\}=\\{25,0.0015,6.5\\}$.
Initialized at rest and isotropic.

  7. 7.

Intermittent double helix flow state for parameters
$\tilde{P}=\\{\tilde{K}^{-1/2},\tilde{W},\tilde{T}^{1/2}\\}=\\{25,250,9.47\\}$.
More chaotic behaviours, in which one helix sporadically dominates or the
helical structure temporarily disperses. Coloured by axial velocity
$u_{\parallel}$. Initialized at rest and isotropic. Time intervals in units of
$10^{3}$ LB time steps between frames.

  8. 8.

Active source of helicity $\Sigma^{\text{a}}$ for the grinder train
(parameters
$\tilde{P}=\\{\tilde{K}^{-1/2},\tilde{W},\tilde{T}^{1/2}\\}=\\{20,250,0\\}$).
Red and blue contours show $\Sigma^{\text{a}}=\pm 175$. Time is stated in
units of $10^{2}$ LB time steps. Initialized at rest and isotropic.

  9. 9.

Active source of helicity $\Sigma^{\text{a}}$ for the double helix (parameters
$\tilde{P}=\\{\tilde{K}^{-1/2},\tilde{W},\tilde{T}^{1/2}\\}=\\{25,1.5,9.47\\}$.
Red and blue contours show $\Sigma^{\text{a}}=\pm 80$. The disclination lines
are highlighted as yellow contours of $S=0.2$. Time is stated in units of
$10^{2}$ LB time steps. Initialized at rest and isotropic.

## References

  * Vicsek  _et al._ (1995) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Phys. Rev. Lett. 75, 1226 (1995).
  * Toner and Tu (1995) J. Toner and Y. Tu, Phys. Rev. Lett. 75, 4326 (1995).
  * Toner and Tu (1998) J. Toner and Y. Tu, Phys. Rev. E 58, 4828 (1998).
  * Solon  _et al._ (2015) A. P. Solon, H. Chaté, and J. Tailleur, Phys. Rev. Lett. 114, 068101 (2015).
  * Chaté  _et al._ (2008) H. Chaté, F. Ginelli, G. Grégoire, and F. Raynaud, Phys. Rev. E 77, 046113 (2008).
  * Digregorio  _et al._ (2018) P. Digregorio, D. Levis, A. Suma, L. F. Cugliandolo, G. Gonnella, and I. Pagonabarraga, Phys. Rev. Lett. 121, 098003 (2018).
  * Chattopadhyay  _et al._ (2021) J. Chattopadhyay, S. Pannir-Sivajothi, K. Varma, S. Ramaswamy, C. Dasgupta, and P. K. Maiti, Phys. Rev. E 104, 054610 (2021).
  * Santhosh  _et al._ (2020) S. Santhosh, M. R. Nejad, A. Doostmohammadi, J. M. Yeomans, and S. P. Thampi, J. Stat. Phys. , 1 (2020).
  * Giordano  _et al._ (2021) M. G. Giordano, F. Bonelli, L. N. Carenza, G. Gonnella, and G. Negro, EPL (Europhysics Letters) 133, 58004 (2021).
  * Bär  _et al._ (2020) M. Bär, R. Großmann, S. Heidenreich, and F. Peruani, Annual Review of Condensed Matter Physics 11, 441 (2020).
  * Nagel  _et al._ (2020) A. M. Nagel, M. Greenberg, T. N. Shendruk, and H. W. de Haan, Scientific Reports 10, 1 (2020).
  * Nishiguchi  _et al._ (2017) D. Nishiguchi, K. H. Nagai, H. Chaté, and M. Sano, Phys. Rev. E 95, 020601(R) (2017).
  * Voituriez  _et al._ (2006) R. Voituriez, J. F. Joanny, and J. Prost, Phys. Rev. Lett. 96, 028102 (2006).
  * Voituriez  _et al._ (2005) R. Voituriez, J. F. Joanny, and J. Prost, EPL (Europhysics Letters) 70, 404 (2005).
  * Ramaswamy and Rao (2007) S. Ramaswamy and M. Rao, New Journal of Physics 9, 423 (2007).
  * Edwards and Yeomans (2009) S. A. Edwards and J. M. Yeomans, EPL (Europhysics Letters) 85, 18008 (2009).
  * Fielding  _et al._ (2011) S. M. Fielding, D. Marenduzzo, and M. E. Cates, Phys. Rev. E 83, 041910 (2011).
  * Ravnik and Yeomans (2013) M. Ravnik and J. M. Yeomans, Phys. Rev. Lett. 110, 026001 (2013).
  * Woodhouse and Goldstein (2012) F. G. Woodhouse and R. E. Goldstein, Phys. Rev. Lett. 109, 168105 (2012).
  * Jiang  _et al._ (2017) H. Jiang, H. Ding, M. Pu, and Z. Hou, Soft Matter 13, 836 (2017).
  * Theillard  _et al._ (2017) M. Theillard, R. Alonso-Matilla, and D. Saintillan, Soft Matter 13, 363 (2017).
  * Wioland  _et al._ (2013) H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler, and R. E. Goldstein, Phys. Rev. Lett. 110, 268102 (2013).
  * Lushi  _et al._ (2014) E. Lushi, H. Wioland, and R. E. Goldstein, PNAS 111, 9733 (2014).
  * Beppu  _et al._ (2017) K. Beppu, Z. Izri, J. Gohya, K. Eto, M. Ichikawa, and Y. T. Maeda, Soft Matter 13, 5038 (2017).
  * Norton  _et al._ (2018) M. M. Norton, A. Baskaran, A. Opathalage, B. Langeslay, S. Fraden, A. Baskaran, and M. F. Hagan, Phys. Rev. E 97, 012702 (2018).
  * Hamby  _et al._ (2018) A. E. Hamby, D. K. Vig, S. Safonova, and C. W. Wolgemuth, Science Advances 4 (2018).
  * Opathalage  _et al._ (2019) A. Opathalage, M. M. Norton, M. P. N. Juniper, B. Langeslay, S. A. Aghvami, S. Fraden, and Z. Dogic, PNAS 116, 4788 (2019).
  * Hardoüin  _et al._ (2020) J. Hardoüin, J. Laurent, T. Lopez-Leon, J. Ignés-Mullol, and F. Sagués, Soft Matter 16, 9230 (2020).
  * Thijssen  _et al._ (2021) K. Thijssen, D. A. Khaladj, S. A. Aghvami, M. A. Gharbi, S. Fraden, J. M. Yeomans, L. S. Hirst, and T. N. Shendruk, PNAS 118 (2021).
  * Thampi  _et al._ (2016) S. P. Thampi, A. Doostmohammadi, T. N. Shendruk, R. Golestanian, and J. M. Yeomans, Science Advances 2 (2016).
  * Wioland  _et al._ (2016a) H. Wioland, F. G. Woodhouse, J. Dunkel, and R. E. Goldstein, Nat. Phys. 12, 341 (2016a).
  * Wioland  _et al._ (2016b) H. Wioland, E. Lushi, and R. E. Goldstein, New Journal of Physics 18, 075002 (2016b).
  * Nishiguchi  _et al._ (2018) D. Nishiguchi, I. S. Aranson, A. Snezhko, and A. Sokolov, Nat. Comm. 9, 1 (2018).
  * Reinken  _et al._ (2020) H. Reinken, D. Nishiguchi, S. Heidenreich, A. Sokolov, M. Bär, S. H. Klapp, and I. S. Aranson, Communications Physics 3, 1 (2020).
  * James  _et al._ (2021) M. James, D. A. Suchla, J. Dunkel, and M. Wilczek, Nat. Comm. 12, 1 (2021).
  * Khoromskaia and Alexander (2017) D. Khoromskaia and G. P. Alexander, New Journal of Physics 19, 103043 (2017).
  * Shendruk  _et al._ (2017) T. N. Shendruk, A. Doostmohammadi, K. Thijssen, and J. M. Yeomans, Soft Matter 13, 3853 (2017).
  * Doostmohammadi  _et al._ (2017) A. Doostmohammadi, T. N. Shendruk, K. Thijssen, and J. M. Yeomans, Nat. Comm. 8 (2017).
  * Zhang  _et al._ (2020) B. Zhang, B. Hilton, C. Short, A. Souslov, and A. Snezhko, Phys. Rev. Research 2, 043225 (2020).
  * Samui  _et al._ (2021) A. Samui, J. M. Yeomans, and S. P. Thampi, Soft Matter (2021).
  * Hemingway  _et al._ (2016) E. J. Hemingway, P. Mishra, M. C. Marchetti, and S. M. Fielding, Soft Matter 12, 7943 (2016).
  * Wensink  _et al._ (2012) H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, PNAS 109, 14308 (2012).
  * Bratanov  _et al._ (2015) V. Bratanov, F. Jenko, and E. Frey, PNAS 112, 15048 (2015).
  * Carenza  _et al._ (2020) L. N. Carenza, L. Biferale, and G. Gonnella, EPL (Europhysics Letters) 132, 44003 (2020).
  * Alert  _et al._ (2020) R. Alert, J.-F. Joanny, and J. Casademunt, Nat. Phys. 16, 682 (2020).
  * Słomka and Dunkel (2017) J. Słomka and J. Dunkel, Phys. Rev. Fluids 2, 043102 (2017).
  * Duclos  _et al._ (2020) G. Duclos, R. Adkins, D. Banerjee, M. S. Peterson, M. Varghese, I. Kolvin, A. Baskaran, R. A. Pelcovits, T. R. Powers, A. Baskaran, F. Toschi, M. F. Hagan, S. J. Streichan, V. Vitelli, D. A. Beller, and Z. Dogic, Science 367, 1120 (2020).
  * Shendruk  _et al._ (2018) T. N. Shendruk, K. Thijssen, J. M. Yeomans, and A. Doostmohammadi, Phys. Rev. E 98, 010601(R) (2018).
  * Čopar  _et al._ (2019) S. Čopar, J. Aplinc, i. c. v. Kos, S. Žumer, and M. Ravnik, Phys. Rev. X 9, 031051 (2019).
  * Chandrakar  _et al._ (2020) P. Chandrakar, M. Varghese, S. Aghvami, A. Baskaran, Z. Dogic, and G. Duclos, Phys. Rev. Lett. 125, 257801 (2020).
  * Chandragiri  _et al._ (2020) S. Chandragiri, A. Doostmohammadi, J. M. Yeomans, and S. P. Thampi, Phys. Rev. Lett. 125, 148002 (2020).
  * Binysh  _et al._ (2020) J. Binysh, i. c. v. Kos, S. Čopar, M. Ravnik, and G. P. Alexander, Phys. Rev. Lett. 124, 088001 (2020).
  * Ruske and Yeomans (2021) L. J. Ruske and J. M. Yeomans, Phys. Rev. X 11, 021001 (2021).
  * Varghese  _et al._ (2020) M. Varghese, A. Baskaran, M. F. Hagan, and A. Baskaran, Phys. Rev. Lett. 125, 268003 (2020).
  * Houston and Alexander (2021) A. J. Houston and G. P. Alexander, arXiv preprint arXiv:2106.15424 (2021).
  * Ellis  _et al._ (2018) P. W. Ellis, D. J. Pearce, Y.-W. Chang, G. Goldsztein, L. Giomi, and A. Fernandez-Nieves, Nat. Phys. 14, 85 (2018).
  * Pearce  _et al._ (2019) D. J. G. Pearce, P. W. Ellis, A. Fernandez-Nieves, and L. Giomi, Phys. Rev. Lett. 122, 168002 (2019).
  * Maroudas-Sacks  _et al._ (2021) Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun, and K. Keren, Nature Physics 17, 251 (2021).
  * Marenduzzo  _et al._ (2007) D. Marenduzzo, E. Orlandini, M. E. Cates, and J. M. Yeomans, Phys. Rev. E 76, 031921 (2007).
  * Thampi  _et al._ (2014) S. P. Thampi, R. Golestanian, and J. M. Yeomans, Philos. Trans. Royal Soc. A 372, 20130366 (2014).
  * Rivas  _et al._ (2020) D. P. Rivas, T. N. Shendruk, R. R. Henry, D. H. Reich, and R. L. Leheny, Soft Matter 16, 9331 (2020).
  * Aryasova  _et al._ (2004) N. Aryasova, Y. Reznikov, and V. Reshetnyak, Molecular Crystals and Liquid Crystals 412, 351 (2004).
  * Giomi  _et al._ (2008) L. Giomi, M. C. Marchetti, and T. B. Liverpool, Phys. Rev. Lett. 101, 198101 (2008).
  * Delbende  _et al._ (2021) I. Delbende, C. Selçuk, and M. Rossi, Phys. Rev. Fluids 6, 084701 (2021).
  * Song and Brenner (2009) Y. S. Song and H. Brenner, J. Chem. Phys. 131, 224108 (2009).
  * Goldstein  _et al._ (2008) R. E. Goldstein, I. Tuval, and J.-W. van de Meent, PNAS 105, 3663 (2008).
  * van de Meent  _et al._ (2008) J.-W. van de Meent, I. Tuval, and R. E. Goldstein, Phys. Rev. Lett. 101, 178102 (2008).
  * Verchot-Lubicz and Goldstein (2010) J. Verchot-Lubicz and R. E. Goldstein, Protoplasma 240, 99 (2010).
  * Woodhouse and Goldstein (2013) F. G. Woodhouse and R. E. Goldstein, PNAS 110, 14132 (2013).
  * Pearce (2020) D. J. G. Pearce, New Journal of Physics 22, 063051 (2020).
  * Napoli and Turzi (2020) G. Napoli and S. Turzi, Phys. Rev. E 101, 022701 (2020).
  * Hozumi  _et al._ (2006) S. Hozumi, R. Maeda, K. Taniguchi, M. Kanai, S. Shirakabe, T. Sasamura, P. Spéder, S. Noselli, T. Aigaki, R. Murakami, and K. Matsuno, Nature 440, 798 (2006).
  * Oriola  _et al._ (2020) D. Oriola, F. Jülicher, and J. Brugués, PNAS 117, 16154 (2020).
  * Novak  _et al._ (2018) M. Novak, B. Polak, J. Simunić, Z. Boban, B. Kuzmić, A. W. Thomae, I. M. Tolić, and N. Pavin, Nat. Comm. 9, 1 (2018).
  * Tortora and Lavrentovich (2011) L. Tortora and O. D. Lavrentovich, Proceedings of the National Academy of Sciences 108, 5163 (2011).
  * Koning  _et al._ (2014) V. Koning, B. C. van Zuiden, R. D. Kamien, and V. Vitelli, Soft Matter 10, 4192 (2014).
  * Nayani  _et al._ (2015) K. Nayani, R. Chang, J. Fu, P. W. Ellis, A. Fernandez-Nieves, J. O. Park, and M. Srinivasarao, Nature Communications 6, 1 (2015).
  * Jeong  _et al._ (2015) J. Jeong, L. Kang, Z. S. Davidson, P. J. Collings, T. C. Lubensky, and A. G. Yodh, Proceedings of the National Academy of Sciences 112, E1837 (2015).
  * Słomka and Dunkel (2017) J. Słomka and J. Dunkel, PNAS 114, 2119 (2017).
  * Moffatt (2021) H. Moffatt, J. Fluid Mech. 914 (2021).
  * Denniston  _et al._ (2004) C. Denniston, D. Marenduzzo, E. Orlandini, and J. Yeomans, Philos. Trans. Royal Soc. A 362, 1745 (2004).
  * Fournier and Galatola (2005) J.-B. Fournier and P. Galatola, EPL (Europhysics Letters) 72, 403 (2005).
  * Suh  _et al._ (2019) A. Suh, M.-J. Gim, D. Beller, and D. K. Yoon, Soft Matter 15, 5835 (2019).
  * Moffatt and Tsinober (1992) H. K. Moffatt and A. Tsinober, Annual Review of Fluid Mechanics 24, 281 (1992).
  * Chono  _et al._ (1998) S. Chono, T. Tsuji, and M. M. Denn, Journal of Non-Newtonian Fluid Mechanics 79, 515 (1998).
  * Gay-Balmaz and Tronci (2011) F. Gay-Balmaz and C. Tronci, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 467, 1197 (2011).


# Evaluating $BDP\\_FRAME$ extension for QUIC

Nicolas Kuhn  CNES Francklin Simo  VIVERIS TECHNOLOGIES David Pradas  VIVERIS
TECHNOLOGIES Emile Stephan  ORANGE LABS

###### Abstract

The first version of QUIC has recently been standardized by the IETF. The
framework of QUIC enables the proposition, negociation and exploitation of
extensions to adapt some of its mechanisms. As one example, the DATAGRAM
extension enables the unreliable transmission of data.

The $BDP\\_FRAME$ extension is a method that can improve traffic delivery by
allowing a QUIC connection to remember the knowledge of path characteristics
and exploit them when resuming a session.

This technical report presents the rationale behind fast convergence in SATCOM
systems and evaluate the $BDP\\_FRAME$ extension in emulated and live
environments.

###### Index Terms:

VPN, SATCOM, PEP, TCP

##  I Introduction

The protocols deployed in the extremities can hardly be relevant for each of
the links available on the Internet, due to their diversity ranging from “very
high speed low latency” from data centers to “high throughput high latency”
from satellite systems.

Systems exploiting satellites in geostationary orbit see their throughput
increased to provide offers comparable to those of terrestrial systems. This
increase in throughput, combined with the intrinsic latency of these systems,
impacts congestion controls such as TCP. In order to make full use of the
available capacity, these systems split the TCP connections into sub-segments
to use suitable congestion control on the satellite segment [1].

These solutions can be applied at different levels of the protocol stack and
are generally at the level of the transport layer in SATellite COMmunication
(SATCOM), and in in particular by adapting the Transmission Control Protocol
(TCP). With a TCP PEP, the packet losses are distributed over three sub-
segments and the control of congestion can be adapted on the satellite link.
This can result in a reduction by half the loading time of a web page.

The end-to-end deployment of QUIC challenges these adaptations. Indeed, with
QUIC (RFC8999 [2], RFC9000 [3] and RFC9001 [4]), the functionalities
previously distributed between Hypertext Transfer Protocol (HTTP)1/1.1/2,
Transport Layer Security (TLS) and TCP are shared between HTTP3, QUIC and UDP.
As the UDP protocol is not in connected mode, it is not possible to cut out
end-to-end communication as with TCP. In addition, QUIC was measured as taking
a not insignificant part of the volume transmitted in a broadband access, in
particular because of the actors implementing it, which are Google or
Facebook. It is therefore necessary understand the performance of the QUIC
protocol in a deployment context that has not been, a priori, considered in
its design.

The authors of [5] identify that the main challenges ahead of QUIC in SATCOM
systems are: (1) tolerance to packet loss, (2) adapted buffer sizes at the end
points, (3) adequate congestion control and (4) quickly exploiting the
available capacity. As the server does not have a priori knowledge of the
characteristics of the underlying system, the convergence of congestion
control can introduce a significant delay before the flow of the communication
is operated. An extension do QUIC currently under discussion at the IETF[6] is
to record the RTT and congestion window parameters during a first connection
and send them to the client in a $BDP\\_FRAME$. When the client wishes to
reconnect to the server, the client returns this $BDP\\_FRAME$ and the
parameters of the previous connection can be used.

This technical report presents the rationale behind fast convergence in SATCOM
systems and evaluate the 0rtt-bdp extension in emulated and live environments.

##  II Configurations

This section presents the different configurations that are considered in this
technical report.

The generic SATCOM system for a broadband access is shown in Figure 1. When
QUIC is exploited at the end points, the TCP-PEP component disapears as shown
in Figure 2.

Figure 1: End-to-end SATCOM broadband access - with TCP-PEP Figure 2: End-to-
end SATCOM broadband access - with QUIC

The configuration of the considered QUIC implementations are shown in Figure
3, the configuration of the end points TCP in Figure 4 and the configuration
of the PEP in Figure 5. The PEP exploits PEPSal.

Figure 3: QUIC implementations and configurations Figure 4: TCP end points
configurations Figure 5: PEP configurations

Unless specify otherwise, the bottleneck bandwidth on the forward link of the
SATCOM system is set to $50$ Mbps and the return link to $10$ Mbps.

##  III Rationale behind congestion control convergence

This section provides some emulation results to illustrate the rationale
behind better congestion control convergence and the relation between achieved
throughput and transmitted file size.

The tests presented in this section exploits a two peers network. More details
on the set up and the exploited platform can be found in GITHUB
111https://github.com/NicoKos/openbach-example-simple.

QUIC implementation in picoquic using BBR congestion control is exploited to
transfer 500 KB, 1 MB, 10 MB and 100 MB files over RTT $\in$ [100 ms; 500 ms]
and bottlenecks of 1 Mbps (forward) / 100 kbps (return), 10 Mbps (forward) / 2
Mbps (return), 50 Mbps (forward) / 25 Mbps (return) and 200 Mbps (forward) /
100 Mbps (return). Figures 6, 7, 8 and 9 represent the used bandwidth for
various forward link bottleneck rate as a function of the RTT. The used
bandwidth is computed as the ratio between the experienced goodput (ratio
between the file size and the time it required to download) and the bottleneck
data rate.

Figure 6: 500KB file size - used bandwidth as a function of the RTT for
various bottleneck data rates Figure 7: 1MB file size - used bandwidth as a
function of the RTT for various bottleneck data rates Figure 8: 10MB file size
- used bandwidth as a function of the RTT for various bottleneck data rates
Figure 9: 100MB file size - used bandwidth as a function of the RTT for
various bottleneck data rates

The main conclusions are the following:

  * •

with a 10 MB file and a data rate of 1 Mbps, the bottleneck is used for all
RTT;

  * •

for shorter files (e.g. 1 MB), increasing the RTT severely impacts link
utilization;

  * •

when the data rate is high (e.g. 250 Mbps), even a 100 MB transfer does not
utilize the available capacity;

  * •

increasing the file size increases the link utilization.

Moreover, the results presented in this section show that there are many cases
where the available capacity is not exploited. This is the case for large RTT
but also the case for small RTT. The convergence of the congestion control has
an impact on the transfer time and in general, the simple estimation of the
transfer time using the file size and the data rate is wrong. Indeed, this
estimation considers instant convergence of the congestion control.

##  IV Performance improvements of $BDP\\_FRAME$

###  IV-A $BDP\\_FRAME$ extension in a nutshell

[6] presents three methods to exploit the transport parameters of a previous
session when resuming a session: (1) local storage, where the server stores
the parameters without negociations with the client, (2) NEW_TOKEN, where a
token that the client can not read is exploited and (3) $BDP\\_FRAME$
extension, where a token that the client can read is exploited. The core idea
of the $BDP\\_FRAME$ extension is illustrated in Figure 10 and resides in the
following :

  * •

during a previous session, the RTT, the congestion window and the IP of the
peers are registered by the server and stored in a $BDP\\_FRAME$;

  * •

during this same previous session or at the end of it, the $BDP\\_FRAME$ is
sent to the client;

  * •

when resuming a session to the same server, the client sends the $BDP\\_FRAME$
to the server;

  * •

the server can exploit the parameters contained in the $BDP\\_FRAME$ to adapt
the congestion control parameters.

Figure 10: Illustration of the $BDP\\_FRAME$ exploitation

More details on the $BDP\\_FRAME$ and the “0RTTBDP” activity in QUIC can be
found in [6]. The local storage option, where the server stores the parameters
when resuming a session is implemented in picoquic, along with a safety check
that the parameters have not change since previous session. The $BDP\\_FRAME$
option is also available in picoquic.

###  IV-B Emulated performance of $BDP\\_FRAME$

Figure 11 presents the received data rate as a function of time for QUIC, QUIC
with 0RTT, QUIC with the $BDP\\_FRAME$ extension, TCP with a PEP and TCP
without a PEP. This illustrates that the 0RTT enables a faster connection
establishment than QUIC without 0RTT but also illustrates how the
$BDP\\_FRAME$ furthers helps in quickly exploiting the available capacity.

Figure 11: Received data as a function of time with QUIC or TCP, with or
without 0RTT

Figure 12 shows the transfer time of $500$ KB with three different
implementations of QUIC with and without 0RTT, with TCP and with TCP-PEP.

Figure 12: 500KB transfer time in SATCOM context with QUIC or TCP, with or
without 0RTT

A good news for SATCOM systems is that the three QUIC implementations provide
results close to those of TCP. There are even cases (quicly without 0RTT or
quicly with 0RTT or ngtcp2) where QUIC performs as good as TCP-PEP
configuration.

The download time is reduced by 13 % with the usage of the 0RTT with ngtcp2
and the reduction is neglectible with quicly. Despite the important RTT of
SATCOM systems, using the 0RTT does not seem to contribute much to the
reduction of the transfer time of short files.

With picoquic implementation of QUIC, the exploitation of the 0RTT reduces the
transfer time by 14 % and by 58 % with the $BDP\\_FRAME$ extension.

##  V 0RTTBDP and variable network conditions

Section IV measured that the $BDP\\_FRAME$ extension can result in very
important reduction of file transfer time. However, the approach may be
considered as aggressive and it is necessary to assess its performance in
presence of competitive flows.

###  V-A Performance of $BDP\\_FRAME$ in congested environments

The scenario that is considered in this section is the following. At $t=0$
seconds, a first QUIC connection transmits 100 MB between two peers (purple
curve). Then, at $t=60$ seconds, the QUIC connection resumes the connection
(blue curve). However, in the meantime, at $t=45$ seconds, a TCP flow has
grabbed the available capacity (green curve). This scenario mesures the
received bit rate as a function of time in Figure 13 (without 0RTT), Figure 14
(with 0RTT) and Figure 15 (with $BDP\\_FRAME$). As a reminder, the 0RTT
variant also exploits the characteristics of the previous session and both
0RTT and $BDP\\_FRAME$ variants introduce safety checks as those proposed in
[6].

Figure 13: Resume session without 0RTT Figure 14: Resume session with 0RTT
Figure 15: Resume session with $BDP\\_FRAME$

Contrary to what could have been expected, because the exploitation of
previous session parameters comes with safety nets, the resumed session is
less agressiv than a standard slow-start approach.

###  V-B Performance of $BDP\\_FRAME$ in real life

To further assess the performance of the $BDP\\_FRAME$ and its deployability,
we have tested it using a real satellite broadband access. The exploitated
architecture is shown in Figure 16.

Figure 16: Real satellite acess : architecture

The picoquic server is hosted in private owned servers on the Internet. We
exploited the KaSat satellite and upload 500 KB or 1 MB files, 50 times each,
and use different picoquic clients.

The results of the experiments are shown in Table I. Using the 0RTT as opposed
to the 1RTT can result in up to 33 % reduction in transfer time for 500 KB.
With the exploitation of the $BDP\\_FRAME$, the transfer time is reduced by up
to 67 %.

TABLE I: Real satellite acess : Download Time (DT) of 500 kB and 1 MB | Without 0RTT | With 0RTT | With $BDP\\_FRAME$
---|---|---|---
File size (MB) | 0.5 | 1 | 0.5 | 1 | 0.5 | 1
Min DT (s) | 3.12 | 3.87 | 2.43 | 3.19 | 1.87 | 2.78
Avg DT (s) | 11.34 | 17.15 | 7.59 | 10.24 | 4.24 | 6.47
Max DT (s) | 47.82 | 61.43 | 33.69 | 33.92 | 15.55 | 23.88

##  VI Conclusion

This technical report illustrates the need for considering congestion control
convergence when determining the transfer time of short objects. Just dividing
the file size by the bottleneck data rate may provide wrong assumptions on the
actual transfer time.

To improve the convergence of the congestion control in SATCOM environments,
this technical report measures the performance of the $BDP\\_FRAME$ extension
for QUIC that exploits the previously measured path characteristics when
resuming a session. With picoquic implementation of QUIC, the exploitation of
the 0RTT reduces the transfer time by 14 % and by 58 % with the $BDP\\_FRAME$
extension.

The approach may be considered as aggressive and it is necessary to assess its
performance in presence of competitive flows. The evaluations illustrate that
exploiting the previous parameters with a safety check can be less agressive
than a standard slow-start mechanism.

## Acknowlegments

The authors would like to thank Olivier Martinez for his help in setting up
the architecture for the real satellite experiments. The authors thank
Christian Huitema for his fruitfull contribution to the 0RTTBDP and the
implementation of the safety guidelines in picoquic. The authors would like to
thank Gorry Fairhurst and Tom Jones for their significant discussion on
0RTTBDP.

## References

  * [1] J. Border and M. Kojo and J. Griner and G. Montenegro and Z. Shelby, “Performance Enhancing Proxies Intended to Mitigate Link-Related Degradations,” Internet Requests for Comments, RFC Editor, RFC 3135, June 2001\.
  * [2] M. Thomson, “Version-independent properties of quic,” Internet Requests for Comments, RFC Editor, RFC 8999, May 2021.
  * [3] J. Iyengar and M. Thomson, “Quic: A udp-based multiplexed and secure transport,” Internet Requests for Comments, RFC Editor, RFC 9000, May 2021.
  * [4] M. Thomson and S. Turner, “Using tls to secure quic,” Internet Requests for Comments, RFC Editor, RFC 9001, May 2021.
  * [5] N. Kuhn, F. Michel, L. Thomas, E. Dubois, and E. Lochin, “Quic: Opportunities and threats in satcom,” in _2020 10th Advanced Satellite Multimedia Systems Conference and the 16th Signal Processing for Space Communications Workshop (ASMS/SPSC)_ , 2020, pp. 1–7.
  * [6] N. Kuhn, E. Stephan, G. Fairhurst, T. Jones, and C. Huitema, “Transport parameters for 0-rtt connections,” Working Draft, IETF Secretariat, Internet-Draft draft-kuhn-quic-0rtt-bdp-09, June 2021, https://www.ietf.org/archive/id/draft-kuhn-quic-0rtt-bdp-09.txt. [Online]. Available: https://www.ietf.org/archive/id/draft-kuhn-quic-0rtt-bdp-09.txt